Simulating Jet-Feedback on Kpc Scales: A Review

1. Introduction

1.1. A Brief Overview of Classical AGN Feedback

Feedback from supermassive blackholes (SMBH) in large early type galaxies has been strongly established as a major influencer of galaxy evolution [1,2]. However, the exact mechanism of how the active galactic nuclei (AGN) affect the galaxy and its environment, and the different implications of this on the galaxy’s properties, is still not settled. From a historical perspective, since the advent of X-ray observations of galaxy clusters, cooling flows of gas cooled via thermal Bremsstrahlung from the cluster environment [3] have been both postulated from theoretical modeling [4] and observationally confirmed [5,6]. However, the fate of the gas cooled below X-ray emitting temperatures ($\lesssim 1-2$ keV) hwas left uncertain, due to lack of distinct observational signatures [7,8]. This prompted considerations of re-heating of the gas by some mechanism, with feedback from the AGN being a viable source[1,9]. The early concept of AGN feedback had been primarily proposed to investigate two major implications: a) explain the the well known $M-\sigma$ relation due to the co-evolution of the SMBH and the galaxy-core’s [10] b) a heating mechanism to offset over-cooling of cluster cores [9,11]. These two different tracks, eventually led to the evolution of the concept of dual-mode feedback by AGN, viz. a) Quasar or Establishment mode related to the local impact of AGN driven outflows and co-evolution of the SMBH and galaxy mass, and b) Radio or Maintainence mode, catering to the large scale heating of gas reservoirs external to the galaxy and preventing cooling flows. In this dual mode scenario, the role of relativistic jets have been largely confined to their impact on extra-galactic gas, for the Radio/Establishment mode, whereas non-relativistic winds in high Eddington ratio systems have been considered to be the primary driver of Quasar/Establishment mode feedback. However, in the recent decades, a large body of studies have demonstrated both from theory and observations that jets can have a significant impact on the ISM of the host galaxy. This makes the earlier dual mode distinction ambiguous in some cases, requiring re-thinking of the traditional defintions [see [2] for a discusison].

1.2. Scope of the Current Review

Over the last few decades there have been some excellent reviews on different aspects of the topic of relativistic jets and their feedback, by various authors. However, their scope and focus have been different and often non-overlapping. Some have highlighted the advances on the physics of relativistic jets themselves [12] or their simulations [13]. Others have given excellent overviews on the astrophysical implications of jets and outflows in general [14,15], and on the topic of AGN feedback in particular [1,2,16,17,18,19,20]. Some related reviews have focused on on other diverse issues related to the topic, such as AGN sources and their radio-loud counterparts [21,22,23,24], gas in and around AGN host galaxies leading to feeding and feedback [25,26,27,28], episodic nature of AGN outbursts [17] etc. The above works provide broad overviews of the various complex astrophysical processes interlinked with the topics of feedback and galaxy evolution. However, there have been only a few detailed reviews that have discussed the complex issues of the interactions of such outflows, specifically jets, with the host galaxy itself [29,30,31,32]. In this review we primarily focus on the development of simulation techniques involved in studying jets in general, with a focus on developments over the past few decades dedicated to the study interaction of jets with their environment. The review is not meant to be a comprehensive summary of all accumulated results till date. Rather, it aims to provide broad summary of the major achievements in this field and their historical developments, to place them in the context of the general question of AGN feedback and galaxy evolution.

2. Modelling Jet Driven Feedback at Galactic Scales

2.1. Jets in Homogeneous Medium

In the mid and late 1970s, there were several seminal developments of theoretical models to explain the dynamics and emission from extra-galactic relativistic jets [such as the `twin-exhaust’ & beam models, the B-Z jet-launch mechanism, diffusive shock acceleration models [34,35,36], etc.], which defined the studies of jets and non-thermal emission in future. Attempts at simulating such jet beams have been made even at such early stages as well, although with limited resolution [37]. The first detailed 2D simulations of propagation of hypersonic jet beams and their structures were presented by Yokosawa et al. [38] and Norman et al. [39], published nearly at the similar times in 1982, the latter paper being more widely recognized in literature. The Yokosawa et al. [38] paper showed that nature of the jet beam (ballistic vs turbulent) and formation of well defined backflows depends on the relative density contrast of the jet and ambient media, and jet velocity. Norman et al. [39] presented more detailed exploration of jetted beams which confirmed the structure of the jet-beam, working surface and backflow; as proposed in Blandford and Rees [34] (see Figure 1). These works spawned several numerical works to probe different aspects of dynamics of supersonic jet beams such as beamed synchrotron emission [40], 3D generalisation [41,42,43,44], stability of slab jets [45], MHD simulations [46,47] etc. Future works have built on the early success of such numerical simulations with larger domain, grid sizes and resolution, although true convergence of the cocoon and beam structures remain elusive [48,49] due to the small scale structures generated with increase in resolution.

Relativistic nature of jet flows had been inferred very early [50,51,52,53,54] from observations of super-luminal motion of jet knots in radio studies. Numerical simulations of steady state jets with relativistic solvers had been presented as early as 1987 [55]. Full fledged dynamic simulations with relativistic solvers were presented a decade later [56,57,58,59,60,61]. A key focus of these simulations were firstly to re-confirm the proposed model of the jet structure in the relativistic limit, and explore a larger parameter space of some key properties of the jet. These initial studies established that internal structures of relativistic jets show significant dependence on the Mach number of the jet beam, jet Lorentz factor and internal pressure. Highly relativistic flows are more stable due to longer growth time scales (as later demonstrated in linear stability analyses [62]). On the other hand, lower mach number (progressively hotter jets) show a difference in their behavior; varying from more internal structures in warm jets ($\mathcal{M}\sim 2$[61]) to more stable cocoons for even hotter jets ($\mathcal{M}\lesssim 1.6$ [61]) as the small scale perturbations are either ill-resolved by the numerical grid or KH instabilities do not couple well with the jet flow. Attempts at simulating magnetised relativistic jets have been made been in tandem as well [59,63,64,65]. However, such simulations of jets received more momentum with the development of efficient high resolution shock capturing schemes in the late 90s and early 2000s [66,67,68,69,70,71,72,73].

In later years, improvement of computational facilities and high order numerical schemes has led to a wide range of AGN jet simulations probing diverse topics such as impact of fluid instabilities and on the jet and cocoon structure [e.g., [49,74,75,76,77,78,79,80,81], etc.], the origin of turbulent structures and dynamics of lower power FR-I jets [e.g., [82,83,84,85,86,87,88,89]] and the larger more powerful FR-II counterparts [such as [90,91,92]] etc. Another focus of such works have been to compare the jet dynamics with predictions of semi-analytical models of jet evolution [e.g., [93,94]]. Models based on self-similar expansion has been proposed for analytical simplicity [90,94,95], although such evolution remain doubtful for all of the life-time of the jet [96,97], especially for the early phase of jets [74]. In many cases, the general scaling laws predicted by Begelman and Cioffi [93], duly modified for a power-law background atmosphere, show good fit with simulated results [74,90].

Simulations of large scale jets have been further driven by efforts to understand impact of jets on cluster scale environment and follow the evolution of non-thermal emission [e.g., [98,99,100,101,102,103,104,105,106] etc.]. Such works have given detailed results of the dynamics of jets in large scale medium, energy transfer to the environment and evolution of synchrotron surface brightness and polarization characteristics as a function of jet length. These have further motivated numerical model based scaling laws predicting synchrotron power for a given mechanical power of the jet [107]. In recent years, large scale simulations of jets have also been utilized to address other science goals such as production of ultra-high energy cosmic rays from large scale shocks driven by jets [108,109,110], impact of multi-species fluid on jet dynamics and emission [e.g., [111,112,113]], jet precession and production of X-shaped structures [114,115,116,117,118,119,120], impact of in-situ particle acceleration of non-thermal electrons [121,122,123,124,125,126,127,128] etc., demonstrating the diverse areas of interest on this topic.

2.2. Jets in Inhomogeneous Medium

2.2.1. Non-Relativistic Simulations

The early phase: The earliest of such suggestions of jets interacting with intervening gas clouds was proposed to explain the observed knots in the jet of Cen A, as early as 1979 [54,129]1. Focus of a jet’s impact on a dense environment was renewed in the 90s by studies of Gigahertz Peaked Spectrum (GPS), Compact Steep Spectrum (CSS) or Compact Symmetric object (CSO), as evolutionary phases of radio galaxies [22,134,135,136]. One of the proposed absorption methods to explain the turnover in the spectrum is due to free-free absorption by intervening ionised gas; either as swept up matter in the forward shock or pre-existing clouds engulfed by the evolving bubble of a radio jet [135,137]. This motivated several theoretical simulations to probe the evolutionary stages of a jet through the host’s ISM, as outlined below.

Jet-single cloud interactions:

Some of the earliest 2D simulation of jets drilling through an inhomogeneous ISM were by DeYoung [139], Steffen et al. [140], who considered a random distribution of spherical (or point-like) dense structures to mimic an inhomogeneous ISM. However, the several aspects of the physics of jet-ISM interaction were first elucidated by more simpler configurations of jets piercing an oblique density discontinuity [141,142,143]. These papers highlighted how the jet’s mach disc tilts and regular structure of the jet-cocoon is broken into turbulent vortices, that also promote instabilities. Such results, gave the early hint of disruption of the jet-beam demonstrated by more complex simulations later. Following on, several simulations explored the impact of jets on individual clouds, usually modeled as a sphere [144,145,146,147,148,149], and some with more updated version of the earlierDeYoung [139] papers evolved 3D jets with multiple randomly distributed clouds [e.g., [150]]. Such simulations2 were a first step to understand the impact of jets on an inhomogeneous ISM, besides probing other science goals such as explaining sources with bent radio jet (wide angle tailed) or those with asymmetric hybrid morphologies. Later simulations included more involved physics such as atomic and molecular cooling [160,161], idealised set ups of shear-layers and mixing [162], self-gravity and starformation [163,164] and realistic morphological models of inhomogeneous molecular clouds including all other previously described physics based models[138,165].

Resolved simulations of jet-cloud3 are insightful in providing the details of how the jets/outflows are affected by the presence of a cloud [143,145,166,167] and more importantly the various evolutionary stages of the clouds themselves (e.g., see Figure 2). Different fluid instabilities can be well captured in such resolved simulations[138,160,165], which otherwise become difficult to follow on global scales. Recent works have also included upgraded models of starformation [138,165] to quantify the positive feedback expected due to the compression by radiative shocks from the AGN outflows. However, a drawback of such individual simulations of jet-cloud interaction is that they do not probe the global impact on the ISM at larger scales, and the evolutionary stages of the jet through the ISM, as it interacts with individual clouds.

The first studies of Jets in fractal inhomogeneous medium: A separate line of simulations probed jets in large scale inhomogeneous ISM; a more realistic depiction than the early works of DeYoung [139], Steffen et al. [140]. Although such simulations have moderate resolutions ($\sim 10-20$ cell across a cloud diameter [168,169]) than the previously mentioned single jet-cloud works, they probe the global impact of the outflow on the turbulent structures in the central few kpc of the galaxy. Hence, such simulations are a bridge between the highly resolved jet-ISM interaction of single clouds, to large scale cosmological simulations with much poorer resolutions [$\sim 100$ pc ] which cannot capture the internal structures of molecular clouds in any detail. A key new introduction of such simulations were the use of a fractal density distribution as a realistic model of the ISM. 2D simulations of jets through a fractal ISM were first introduced in the early 2000s [170,171]. The first detailed 3D simulations were presented in Sutherland and Bicknell [33, hereafter SB07], where a non-relativistic jet was injected through a two-phase ISM. This was a pioneering paper in many aspects, as it laid both the technical foundation for several future publications, as well as elucidating the basic evolutionary stages of a jet breaking out through an inhomogeneous ISM. These works were later improved upon by relativistic simulations of jet-ISM interaction, outlined in the next section (Section 2.2.2).

Other non-relativistic jet-ISM simulations with improved physics: Beyond the SB07 works, there has been several other papers, that have explored different facets of the process of jet-ISM interaction with slightly different initialization. One such series viz. Dugan et al. [163], Gaibler et al. [172,173], Dugan et al. [174] focused on simulating larger gas disks ($\sim 30\mathrm{kpc}$), as opposed to a few kpc disks considered in other papers [33,169,175,176]. However, the larger scales necessitated a more modest resolution $\sim 67$ pc, lower by nearly 10 times that of the other resolved studies mentioned earlier. Nonetheless, these simulations truly probed the global impact of an AGN jet on the large scale disk.

These works outlined several important results such as asymmetric jet morphologies due to local inhomogeneities [172], large scale compression driven positive feedback over the galactic scales [173], a ring-like morphology of enhanced SFR surrounding the central cavity and appearance of hyper-velocity stars with strong non-circular velocities [174], comparison of jet vs wind driven feedback and impact of jet orientation [163] which was a pre-cursor to the later study by Mukherjee et al. [169]. These provide definitive predictions to test the applicability of such models in observations, along with identifying past AGN activity from the disturbed stellar kinematics. While the above papers considered an existing fractal density distribution, other recent works have also explored generation of a self-consistent inhomogeneous ISM by stellar feedback, before jet launch [177]. Most of the above works have only explored hydrodynamic simulations, without magnetic fields. Only a handful of papers [177,178] have explored the impact of magnetic fields on jet-cloud interaction, which remains an area to be explored in future.

In recent years, a more detailed self-consistent evolution of jets and their environment have been carried out in another series of publications viz. Fiacconi et al. [179], Talbot et al. [180,181,182]. The authors have a introduced a novel sub-grid blackhole accretion-ejection model based on a thin accretion disk model, duly accounting for mass build up of the accretion disk while following the angular momentum exhange between the in-falling gas, blackhole and the disk. The jet power is determined from a Blandford-Znajeck mechanism, with the efficiency parameterised from GRMHD results [179,180]. The innovative model has been employed to study the mutual evolution of jets and a circum-nuclear disk of height (radius$\sim 70$ pc, height $\sim 9$ pc) [180,181]. Jets in larger kpc scale gas disk were explored in a subsequent paper [182]. The simulations probed outflows from lower-mass blackholes ($\sim {10}^6{M}_{\odot }$) and hence had typically low kinetic power of jets ($\sim {10}^{42}{ \mathrm{ergs}}^{-1}$), although higher power jets were probed for simulations exploring AGN triggering from mergers [182] which also had larger kpc scale gas reservoirs.

Although these simulations explore a different range of parameters, more suited to Seyferts, than typical massive radio-loud AGN, the qualitative evolution of the jet and its interaction with the multi-phase gas is similar to other works (see Section 3 for a general summary). However, one of the key outcomes of Talbot et al. [181] works is the self-consistent evolution of the jet angular momentum and the Bardeen-Peterson effect driven re-orientation of the jets for inclined jets. The simulations also predicted significant cold ($T<{10}^4$ K) outflows from the circum-nuclear disk, with an increase in rates for inclined jets, which conforms well with observations [183]. Although the authors did not find significant blackhole spin evolution during a single outburst, the general applicability of the method makes it very suitable to be implemented in large scale cosmological simulations as demonstrated in Talbot et al. [182], which can potentially track the blackhole growth over cosmic time.

2.2.2. Relativistic Simulations

Why relativistic MHD? A drawback of non-relativistic simulations probing dynamics of jets that are inherently relativistic in their bulk flows, is the difference in the momentum exchange with the external environment. The momentum conservation equation in relativistic hydrodynamics is

$${\displaystyle \frac{\partial }{\partial t}}\left({\gamma }^2\rho h\mathbf{v}\right)+\nabla ·\left({\gamma }^2\rho h\mathbf{v}\mathbf{v}+\mathbf{I}p\right)=0,$$

(1)

where $\rho h=\rho c^2+\rho ϵ+p$ is the relativistic enthalpy, $\rho ϵ$ the internal energy; and $\mathbf{I}$ an identity tensor. The Lorentz factor of the flow velocity is $\gamma$. A jet with a given rest-frame density, pressure and bulk velocity, a relativistic formulation implies higher momentum imparted by the jet beam, at least by a factor of ${\gamma }^2$. For example, the difference (${\gamma }^2-1$) becomes $\sim 10\%$ even for mildly relativistic flows of $\beta \sim 0.3c$ ($\gamma \sim 1.05$). Thus non-relativistic formulation of the fluid equations will under evaluate the momentum advantage imparted by the jet. Of course the total energy flux of a relativistic and non-relativistic jet, with identical fluid parameters, are thus not the same. Hence, non-relativistic simulations of jets have often employed an equivalent jet beam by constraining the jet pressure, velocity and injection radius to be same, but deriving the density to match the total energy flux. The relations between the densities for such an equivalent non-relativistic (${\rho }_{\mathrm{nr}}$) and its relativistic counterpart (${\rho }_r$) were derived by Komissarov and Falle [184] as

$${\rho }_{\mathrm{nr}}=2{\rho }_r{\gamma }^2\left({\displaystyle \frac{\gamma }{\gamma +1}}+{\displaystyle \frac{1}{\chi }}\right),\chi ={\displaystyle \frac{\rho c^2}{\rho ϵ+p}}={\displaystyle \frac{(\Gamma -1)}{\Gamma }}{\displaystyle \frac{\rho c^2}{p}},$$

(2)

where $\chi$ is the ratio of the rest mass energy and the non-relativistic part of the rest-frame enthalpy ($\rho h-\rho c^2$). Equation 2 shows that a flux matched non-relativistic jet has a higher density and hence a heavier jet. This results in narrower jet cocoons, faster jet propagation, higher mach numbers [61,184] and lower cavity pressures [185] than that of a relativistic jets. Thus irrespective of the choice of initial jet parameters, whether done by strictly ignoring relativistic effects, or by deriving effective jet parameters by matching fluxes, the momentum balance is strongly affected by the neglect of relativistic solvers while evolving AGN jets[13,61,184]. However, accuracy of the momentum exchange is crucial for the physics of jet-ISM interaction and the implications for local scale AGN feedback effects by the jets. This necessitates the usage of relativistic solvers in simulations of jet-feedback.

Jets in static fractal ISM: The first such simulations were presented in Wagner and Bicknell [168] and later expanded with a larger set of simulations [186] probing different volume filling factors of the dense gas. These simulations modeled a relativistic jet ploughing through a static fractal ISM [33] immersed in a constant density background halo. The simulations did not have an external gravitational field of the galaxy. The dense ISM was assumed to be distributed spherically, unlike SB07 [33] who considered a disk. The suite of simulations probed several different parameters of the simulations relevant for studying jet-ISM interaction, such as i) jet power: ${10}^{43}-{10}^{46}{ \mathrm{ergs}}^{-1}$, ii) mean cloud density: ${10}^2-{10}^3{ \mathrm{cm}}^{-3}$, iii) volume filling factor4: $f_V\sim 0.027-0.4$ and iv) cloud sizes: $10-50$ pc (See Table 2 of Wagner et al. [186]). All the simulations were run for a jet of $\gamma =10$ and $\chi =1.6$ (see Equation (2)).

These simulations, along with SB07 [33], were the first to explicitly identify the various stages of the jet’s evolution as it channels through an inhomogeneous ISM (discussed in more detail in Section 3.1). They discussed the impact of the jet on the cloud ablation and acceleration, dynamics of the jet driven bubble, their impact on ISM energetics, and placed such cases in the context of AGN feedback processes in galaxies, showing definitively that as jets can significantly impact the host’s ISM. Future simulations have expanded and confirmed these results. Besides the above, these papers highlighted two other primary impacts of jets, not often highlighted in other works:

• ISM with small size clouds ($\lambda \lesssim 10-20$ pc) could easily be cleared by the jet and result in mean radial velocities higher than the stellar velocity dispersion with moderate values of Eddington ratios ($\eta =P_{\mathrm{jet}}/L_{\mathrm{Edd}}\lesssim {10}^{-4}-{10}^{-3}$). However, sufficiently accelerating larger clouds, typical of galactic GMC[187,188,189] ($\gtrsim 50$ pc), would require higher Eddington ratio ($\eta \gtrsim 0.01-0.1$) for jets of power ($P_{\mathrm{jet}}\gtrsim {10}^{45}{ \mathrm{ergs}}^{-1}$), and hence more efficient outburst from larger mass SMBH. This implies stronger confinement of jets in ISM with larger clouds, which would be more difficult to ablate, as also confirmed in later works [175].
• The jets were found to provide a strong mechanical advantage, higher than unity (defined as the ratio of total outward momentum of clouds to the net momentum imparted by the jet). This also correlated with kinetic energy transfer to the ISM $\sim 20-30\%$, later refined to slightly lower values by future works [175,190].

Jet-ISM interaction in dynamic environments: The above works were taken forward by four subsequent papers: Mukherjee et al. [175], Bicknell et al. [191],Mukherjee et al. [192] and Mukherjee et al. [169], as a direct continuation of the earlier efforts. These papers had a few new additions: i) An external gravitational potential (double isothermal) with a hydrostatic atmosphere, which enabled more correct calculations of long term gas kinematics (see Section 3.3). ii) A dynamic turbulent velocity dispersion in the dense fractal ISM. This added further realism to the set ups, which successfully enabled ready comparison with several observed sources (see Section 4). The results from these simulations have formed the primary benchmark for studies of jet-ISM interactions in recent years.

The above works primarily focused on high jet powers $P_{\mathrm{jet}}\sim {10}^{44}-{10}^{46}{ \mathrm{ergs}}^{-1}$. An update on these for a broader range of powers were carried out by Tanner and Weaver [176], where a key focus was extending the jet-ISM simulations to much lower power ($P_{\mathrm{jet}}\lesssim {10}^{42}{ \mathrm{ergs}}^{-1}$) and compare the dynamics with higher power cases in the same numerical framework. While these simulations recovered the earlier results of [169], one important distinction was that the confinement of the jet was seen to less strongly depend on the morphology of the ISM for higher power jets $P_{\mathrm{jet}}\gtrsim {10}^{44}{ \mathrm{ergs}}^{-1}$, which eventually drill through the ISM. Lower power jets were prone to disruption and break up, likely resulting from the lower jet density and pressure resulting in lower momentum flux; which would in turn result in longer ablation time scales of the clouds. Similar results have also been supported in more recent non-relativistic simulations [193]. The nature of the jet-ISM interaction for a given jet power was in general however, found to proceed along similar lines of earlier works [169].

3. Summary of Key Results

Although specific individual results often vary between different simulations, due to different choices of parameters of the jet, the ambient medium or the physics involved; a common set of general outcomes can be ascertained. In the following the sections we outline the broad summary of some general conclusions and key results.

3.1. Evolutionary Stages of the Jet Through an Inhomogeneous Medium

A common set of evolutionary phases [first identified by [33] of a jet moving through the ISM of a gas rich host can be identified, as illustrated in Figure 3.

• The Confined phase: The jet remains confined within the clumpy ISM ($\sim 0.5-1$ kpc), resulting in the formation of a flood-channel scenario (see right panel of Figure 1). The jet-plasma is diverted to low density channels through the clouds, percolating into the ISM. The jet beam’s forward progress is halted, leading to a temporary stalling of the jet-head. However, the jet’s energy is dispersed over a quasi-spherical volume as the forward shock sweeps through the ambient medium in the form a energy driven bubble. Simulations find the time of the confined phase can last from a few hundred kilo-years to a $\lesssim 2$ Myr, depending on the power of the jet, the density of the ambient medium and extent of the dense gas. An approximate analytical analysis of the duration of confinement is presented in Appendix A. Since these conditions can vary over wide range between different galaxies, the impact of jet and the efficiency of coupling with the ISM can have a wide variation as well.
• Jet breakout phase: The jet breaks free from the confinement of the dense ISM and the hemi-spherical bubble, to proceed onwards. During this phase the jet driven bubble can still indirectly impact the dense ISM. The bubble remains over-pressurised and eventually engulfs the ISM. This is accentuated by the drafts of backflow from the tip of the jet. This drives shocks into the clouds away from the jet axis, which can raise turbulence [169] and impact starformation over the inner few kpc of the galaxy [169,194,195]. During this phase the jet and its ensuing bubble still has significant impact on the ISM. However, as the decrease in pressure with the expansion of the bubble weakens the impact on the dynamical evolution of the dense gas.
• The classical phase: Beyond the break-out phase, the jet carves a clear path through the ISM. Subsequent energy flows have less impact on the ISM. The jet proceeds into the low density stratified homogeneous halo gas. Beyond this, the dynamics of the jet are similar to the conventional models of jet propagation into a static homogeneous medium. The dynamics of the ISM and perturbed velocity dispersion of the clouds start to decay back to the pre-jet levels [195].

Of the above stages, the confined phase is of primary interest in the context of AGN feedback. During this phase, there is significant coupling of the jet with the dense gas, which results in transfer of $\sim 10-20\%$ of the jet’s energy flux into the ISM in the form of kinetic energy, creating local outflows [175,186,190] and also additionally heating the gas to create radiative shocks. The efficiency of such an interaction depends primarily on the following criteria:

• The volume filling factor of the dense gas ($f_V$).
• Jet’s orientation with respect the ISM morphology, with jets inclined to a gas disk being more productive (${\theta }_j$).
• Jet power ($P_j$).
• Mean density of the clouds in the ISM ($n_c$).

Thus the impact of jet driven feedback on kpc scales depends on a four-dimensional parameter space. The maximal impact of course is in a case with jet directly oriented into the dense ISM (e.g., jets pointed into the disk [169,192,196]), higher jet power and higher cloud densities, leading to longer confinement. More detailed discussion on the impact of parameters 1, 3 and 4 on jet confinement are discussed later in Appendix A. However, the interactions can also be gentle if one of the parameters are weak, even though others are prominent; such as in the case of NGC 3100 [197] where in spite of dense gas disk observed in CO 1-0 with a moderately powerful jet ($P_j\sim {10}^{44}{ \mathrm{ergs}}^{-1}$), the impact of the jet on the disk’s kinematics is minimal. This is likely due to weak coupling arising from the relative orientation of the jet away from the disk’s plane. However, on the other hand, several detailed spatially resolved observations have uncovered more telltale smoking gun signatures of strong impact of the jet on the confining ISM, as discussed later in Section 4.

3.2. Global Impact on the ISM

Simulations of jets through the inhomogeneous ISM have strongly supported that jets can cause a large scale effect on the central few kpcs of the galaxy, contrary to earlier beliefs that such thin collimated structures are less important in the global context [198]. There are again three distinct types of impact of the jet, which has varied effects on the ISM.

• Direct impact of jet-beam ($\lesssim 1$ kpc): Clouds directly along the path of the jet are strongly impacted by the flow and eventually ablated. Such an interaction affects both the clouds as well as the jet. For large clouds (e.g., a GMC of size $\gtrsim 50$ pc) directly along the jet-beam, which may nearly cover the jet’s width, the jet is strongly decelerated till the cloud moves away from the jet’s path or is completely disintegrated. The region of such impact is usually confined to $\lesssim 1$ kpc, where the jet-beam and its ensuing backflow directly interacts with the ISM. This region experiences much higher turbulent velocity dispersion and density enhancement [195] due to stronger ram pressure driven shocks. In addition, the stronger interaction in the central region also results in mass removal and formation of a cavity [169,173,195]. However, simulations that better resolve the cloud structures show that such cavities are not completely devoid of dense gas [192,196]. Strands of dense cloud cores, with a radiative shock enhanced high density outer shell remain embedded inside such cavities, that are slowly ablated by the jet driven flows [169,192].
• Indirect impact by energy bubble ($\gtrsim 1$ kpc): As mentioned earlier in Section 3.1, the confined jet’s energy spreads out in the form of an energy bubble sweeping through the ISM. The indirect interaction operates differently depending on the evolutionary phase of the jet (see Section 3.1). During the jet confinement phase, the forward shock sweeps through the ISM. The embedded clouds face a steady outward radial flow of the jet plasma being re-directed in lateral directions from the jet axis, through the flood-channel mechanism. This results in outward radial flows inside the ISM, away from the jet axis. In the jet breakout phase and beyond, the jet expands beyond the immediate confines along its path and the over-pressured cocoon engulfs the ISM. This is more prominent for gas disks, as shown in the right panel of Figure 3. Such indirect interactions are responsible for more large scale impact of the jet, beyond the central 1 kpc range. This raises the velocity dispersion of the gas in general all through the ISM and also shocks a larger volume of the ISM [199,200]. Inclined jets strongly confined within the ISM [163,169,176,192] are also able to process a large volume as the jet-plasma spreads beyond the decelerated jet-head and drive radial inflows through the ISM.

The fact that jets can in principle affect a larger volume of the ISM than their apparent width near their launch axis, has strong implications for AGN feedback. This demonstrates that jets can creates strong global outflows and impact gas kinematics and starformation rate, as discussed in further below.

3.3. Impact on ISM Kinematics

As discussed above, the jets can affect a significant fraction of the ISM during the confined and break-out phase. This results in launch of fast multiphase outflows [33,163,175,176,181,186,190,192]. The nature of the outflows depend on the four key factors listed in Section 3.1. We list below some of the broad major inferences that can be drawn from these theoretical simulations.

3.3.1. Multi-Phase Outflow:

In a realistic system, one would expect a wide range of gas phases to co-exist: dilute hot gas in the halo of galaxies, dense gas ionised by shocks (collisional) or photons from a central source like AGN [199] or shock-precursors [201], warm molecular gas likely representing cooling fronts or shocks/turbulence that have penetrated dense clouds [202,203], cold dense molecular gas [204,205] and neutral HI gas [31]. However, most simulations consider a single fluid system and do not explicitly track the chemical evolution of the constituents of the fluid due to the computational complexities, although some recent works have started some preliminary investigations to model differences in composition [206,207,208].

Nonetheless, simulations that can resolve the density sub-structures inside an inhomogeneous ISM can track the variation of density, temperature and velocity structures of the fluid, which act as a proxy for the different phases [192,199,209]. Some theoretical papers have attempted to understand the multi-phase nature of the ISM and disentangle the relative contributions of different gas properties in the outflows by evaluating 2D histograms [138,175,190] or analysing the simulations based on temperature thresholds [181]. A better representation of the multi-dimensional nature of the phase space is shown in Figure 4, where the mass distribution is represented in terms of the three primary variables of interest, viz. density (n), temperature (T) and outward radial velocity ($v_r$) corresponding to outflowing gas. The corresponding 2D distributions, obtained by summing along each axis of the 3D distribution, are plotted on the right. However, the summed 2D distributions often fail to capture variation in the phase space visible in the 3D image. The plot is for the last panel of Figure 14 of simulation D from Mukherjee et al. [169], which corresponds to a jet of power $P_j={10}^{45}{ \mathrm{ergs}}^{-1}$ launched at ${45}^{\circ }$ to the axis of the disk.

Several distinct identifiable regions have been highlighted in the 3D figure (Figure 4).

• Cloud cores: There is collection of mass at $T\sim 1000$ K, with high densities near the left face of the 3D figure. This corresponds to the cores of the clouds, with a temperature near the cooling floor of the simulation ($T={10}^3$K). The clouds have some positive radial velocity ($v_r\lesssim 100{ \mathrm{kms}}^{-1}$) which likely correspond to the turbulent bulk velocity of the clouds injected at initialisation and also mild acceleration after jet-ISM interaction.
• Dense warm outflow: There is a distinct collection of mass in Figure 4, shifted from the cloud cores, that is centred around $T\sim {10}^4$ K in temperature and extending from $v_r\sim 100-1000{ \mathrm{kms}}^{-1}$ in velocity and density $n\gtrsim 100{ \mathrm{cm}}^{-3}$. This corresponds to dense shock heated gas that has cooled and accelerated to high velocities. This phase has the highest mass among all of the outflowing gas, and hence accounts for the dominant contributor to the kinetic energy budget of the outflows. In observational studies, this phase would correspond to the warm molecular gas [202,210,211,212] and may also proxy the cold gas outflows [204,209,213] as modelled in Mukherjee et al. [192].
• Shocked cloud layers: Beyond the dense warm phase, there is another distinct, but small, collection of mass peaking between $T\sim {10}^4-{10}^5$ K, and at a lower density ($n\sim 10-100{ \mathrm{cm}}^{-3}$) than the dense warm phase. The temperature range corresponds to the peak of the cooling curve. This phase belongs to either the outskirts of the clouds being shocked by the enveloping pressure bubble or shocked dense cloud-lets ablated from large clouds [169,200]. This phase accounts for the majority of the observed emission in optical lines used as diagnostics of shock ionisation such as O[II], O[III], S[II] etc. [192,200]. It should be noted that the mass represented in this phase is small compared to the dense phase. Hence, inferences of ionised gas mass from shocked gas are often lower limits to the gas mass actually contained in the ISM, but often missed due to lack of multi-wavelength coverage of the system.
• Hot tenuous outflow: The jet driven outflows pushes out the ablated gas in a tenuous hot form ($n\lesssim 10{ \mathrm{cm}}^{-3}$, $T>{10}^6$K). This gas forms a tail of the distribution extending to very low densities and high velocities. Such a hot tenuous gas is predicted to be seen in X-rays wavebands [33]. Detecting the soft X-rays from such thermal gas with sufficient spatial resolution to distinguish it from the central nucleus is challenging, owning to the contribution from the AGN, but has been tenatively confirmed in some sources [214,215].

3.3.2. Galactic Fountain:

Although the jets can launch strong local outflows, there is no significant large scale blow-out of the dense ISM, as often required by semi-analytical models of galaxy evolution in the context of AGN feedback to quench starformation by mass removal. Instead, the total mass weighted mean velocities are in fact negative [175,190]. This indicates, that although there are strong localised outflows, a significant fraction of the ISM stays within the gravitational potential of the host galaxy. This is best demonstrated by the escape fraction plot in Figure 20 of Mukherjee et al. [175]. It shows that only $\lesssim 10\%$ of the ISM moves beyond the central few kpc. Thus outflows, without escape, result in a localised galactic fountain scenario, where gas will likely be re-cycled within the galaxy’s own confines [175].

3.3.3. Turbulent Velocity Dispersion:

A key focus of local AGN feedback studies in both observational and theoretical domains has been to understand the influence of AGN driven outflows on the turbulence of the ISM. Resolved simulations of jet-ISM interactions have shown that as the jet driven bubble sweeps over the central few kpc of the ISM, it can significantly raise the velocity dispersion by up to an order of magnitude from its initial value. However, this seems to depend on the phase of the gas. For example, Mukherjee et al. [169] shows that hotter shock ionised gas ($T>{10}^4$ K) will have higher velocity dispersion ($\sim 400-600{ \mathrm{kms}}^{-1}$) than the colder component ($\lesssim 100{ \mathrm{kms}}^{-1}$). This is because the shocks progress very slowly within the dense cores. The primary impact of the jet driven bubble is on the ablated cloud-lets stripped from the larger clouds. Random bulk motions of such clouds adds to the velocity dispersion. However, the dense gas is not completely undisturbed. Detailed comparison of kinematics of dense gas in galaxies such as IC 5063 [204] and B2-0258 [209,216] with simulations [169,192,200] have shown excellent correspondence of the spread in the observed gas kinematics. An interesting new feature identified in both simulations [192,200] and observations [217,218,219,220,221] is the the appearance of high velocity dispersion in directions perpendicular to the jet. Such features are conjectured to arise from deceleration of a jet strongly inclined in to the gas disk, resulting in outflows of plasma both along the minor axis following the path of least resistance [192,200], as well as into the plane of the gas disk [222]. A combined effect of both types of motions are predicted to result in such apparent enhanced widths perpendicular to the jet [200,222].

3.4. Impact on Starformation Rate

One of the primary focus of the studies of AGN feedback has been to understand the impact of AGN driven outflows on the starformation rate (SFR) of the galaxy, both instantaneous and long term. In simulations starformation rates have been primarily estimated as

$$\mathrm{SFR}=ϵ{\displaystyle \frac{\rho }{t_{\mathrm{ff}}}};t_{\mathrm{ff}}={\left({\displaystyle \frac{3\pi }{32G\rho }}\right)}^{1/2},$$

(3)

where $ϵ$ is often assumed to be a constant efficiency factor ($\lesssim 0.01$) and $t_{\mathrm{ff}}$ is the local free-fall time. Since gas density was the only fluid parameter affecting the above, reduction in SFR was modeled to occur primarily due to evacuation of gas by the outflows and enhancement due to compression by shocks. Spatially resolved simulations of jet-ISM interaction have predicted a central cavity (hence negative feedback) surrounded by a ring of SFR enhanced region due to compression in the immediate rim of the cavity [163,169,173,177]. Strong compression from the ensuing pressure bubble have been predicted to strongly promote starformation in extended regions of the galaxy [173,194]. Such predictions of positive feedback do indeed align with some observed sources, especially with the aligment of starforming streams with the outflows [223,224,225,226]. However, the simplistic density threshold based starformation model as described above often lead to very high SFR ($\gtrsim 300M_{\odot }{yr}^{-1}$), which would imply very high addition of new stellar material ($\sim {10}^9{M}_{\odot }$ [173]) for each outburst. Such high rates of starformation are unlikely for Radio loud AGN and likely point to the deficiencies in the quantitative sub-grid SFR models. Predictions of global negative feedback have been limited in such studies, except in some simulations with low gas mass/density that suffers very strong ablation from the jet [e.g., simulation E of [169]]. However, several radio-loud galaxies have been demonstrated to have significant reduction of their SFR [210,227].

In a recent study, Mandal et al. [195] proposed a new turbulence regulated framework to estimate the SFR in simulations where turbulent gas structures can be resolved. A key difference in this approach is to duly account for variation in the local free-fall time as a function of density, the virial parameter (${\alpha }_{\mathrm{vir}}=2E_{\mathrm{kin}}/E_{\mathrm{grav}}$) and the mach number of the gas ($\mathcal{M}$). The work extends the theoretical framework of turbulence regulated starformation in molecular clouds [228,229]. This method assumes that only densities beyond a certain threshold are gravitationally unstable to form stars. This occurs for densities with Jean’s length lower than the sonic scale (${\lambda }_J={\left(\pi c_s^2/\left(G\rho \right)\right)}^{1/2}\lesssim {\lambda }_s$), where the sonic scale corresponds to the length scales at which the turbulent velocity dispersion is lower than the sound speed. At scales higher than the sonic scale, turbulent pressure offsets the gravitational collapse. Using the above criteria, a turbulence regulated estimate of SFR has been carried out in postprocess of jet-ISM interaction simulations by Mandal et al. [195]. In a departure from either pure positive or negative feedback, the physics driven approach reveals some new aspects, such as i) a mild global reduction in the SFR during the onset of the jet-ISM interaction, ii) in-efficient positive feedback occurs in inner regions directly impacted by the jet that suffer both jet driven compression as well as enhanced turbulence, iii) the ISM going through a sequence of evolutionary phases in the Kennicutt-Schmidt (KS) plot, till the jet-break out phase, beyond which turbulent velocities return to pre-jet levels. Although in its infancy, the inclusion of turbulence regualated starformation in large scale simulations shows promise in addressing several existing issues related to impact of AGN feedback processes on SFR.

4. Observational Implications

As outlined earlier in Section 1.2, several reviews have summarized the observational evidence and implications of AGN feedback. We refer the reader to Harrison and Ramos Almeida [2], a more recent addition to the series. However, in the following sections we briefly summarize the observational results pertaining in particular to jet-ISM interaction, with a narrower focus than the discussion in the reviews mentioned earlier.

4.1. Observations of Jet-ISM Interactions

Observations of jet-ISM interactions have been reported from the very early start of studies of relativistic jets in galaxies. One of the earliest suggestions of expulsion ionised gas co-incident with radio emitting ridges in NGC 4258 by a “nuclear explosion” was reported as early as 1972 [230]. However, the authors at that time discounted that relativistic electrons were streamed directly from the nucleus, assuming propagation at Alfvén speeds, as the concept of non-thermal plasma moving at bulk relativistic speeds was still not well-developed. Later works of this well studied galaxy have established it to harbor radio jets co-planar with the gas disk [231,232,233]. Later in the early 80s, several studies of individual galaxies with co-located radio and ionised gas emission, often with enhanced kinematics and evidence of de-polarisation of radio emission were presented, suggestive of jet-ISM interaction e.g., 3C 66B, 3C 31, NGC 315, 3C 449 [234], 3C 277.3 [235], 3C 305 [236], 3C 171 [237], 3C 293 [238], 3C 310 with a proposed radio bubble from a trapped jet [239] etc. These works were later updated with multi-object studies of radio-optical alignment [240,241,242,243], which provided strong evidence of wide prevalence of impact of jets on ionised gas and potentially positive feedback in starformation. More on the alignment effect is presented in the reviews by McCarthy [244] and O’Dea and Saikia [22].

In recent years, with better access to spatially resolved and more sensitive observations, there has been a plethora of multi-wavelength multi-phase observations of jets interacting with their ambient medium (see Figure 5 for an example). A collection of sources have been presented in Table A1 of Appendix B, along with comments on the nature of the interaction where applicable. The sources listed show signs of interaction of a jet with the ambient gas by demonstration of spatially resolved outflows and/or enhanced velocity dispersion of the gas. The list is not meant to be a complete collection of all observed cases, but a representative sample showing the diverse nature of observations of jet-ISM interaction. In many cases however, presence of strong central AGN also permits driving of outflows from an energetic point of view. However, several factors, especially resolve spatial morphology of the gas kinematics and jet, lends credence to the jet-driven feedback scenario. Besides such spatially resolved studies of individual galaxies, several other papers have explored the tentative observations of jet driven outflows in a broader sample of radio loud galaxies [e.g., [245,246,247,248], etc.], confirming that jet driven outflows strongly influence the gas in their host.

Such observations provide strong support to the results of the spatially resolved simulations outlined earlier in this review. These results have been further augmented by direct predictions of observable signatures of jet-ISM interaction from simulations, such as predictions for ionised gas kinematics by Meenakshi et al. [200]. Similar efforts of direct comparison of observed sources with simulations have been very fruitful in recent years e.g., IC 5063 [192], 4C 31.04 [203], B2 0258+35 [209,214,216,249], 2MASSX J23453269-044925 [227,250], Tea Cup galaxy [222] etc., highlighting the growing synergy between spatially resolved simulations and observations.

4.2. Implications for Compact and Peaked Sources (CSS/GPS/CSO)

Studies of compact jets are directly related to and motivated by the confined young radio jets discussed in this review . The origin and nature of compact peaked spectrum sources, on whether they are young newly evolving sources or trapped by the host’s ISM are still debatable [22,136]. Synchrotron based age estimates of such sources with linear size $\sim 1$ kpc, have been found to support the “young” source scenario [251]. However, the short confinement times can also be due to the powerful ($P_{1.4}GHz\gtrsim {10}^{25}{WHz}^{-1}$) nature of these sources [252,253]. There are many ancillary evidence to suggest that such sources are gas rich, such as high rotation measures [see [22,136,254,255,256], and references therein], high column depths in X-ray observations [257,258], enhanced IR emission as indication of reprocessing from dust [253,259]. There are also direct detection of atomic and molecular gas in many such systems, as described in reviews by O’Dea [136], Fanti [260], O’Dea and Saikia [22]. Recent observations [261] have also uncovered strong ionized outflows in a large number of such compact systems, indicating strong coupling of the jet with the gas, as predicted by simulations outlined in Section 3.3.1. Hence, given the abundance of possibilities of jet-ISM interaction in these systems, such galaxies may indeed follow the evolutionary sequence outlined in Section 3.1 and remain confined for some duration (see Appendix A for a derivation of the confinement period). Besides the traditional peaked spectrum sources, recent surveys have uncovered a large sample of compact FR0 sources [24,262] where the jet remains compact. Compact jet-like systems have also been uncovered in resolved radio imaging of traditional radio-quiet systems as well [263,264]. Longer confinement time ($\tau \gtrsim 5$ Myr) of low power jets ($P_j\lesssim {10}^{41}{ \mathrm{ergs}}^{-1}$), even for modest mean densities of the ISM clouds ($n_c\sim {10}^2-{10}^3$) can explain compactness of such sources. Indeed, significant detection of molecular gas has been detected in several such galaxies with compact radio emission [265,266]. This makes such compact radio sources ideal test beds for the physical processes of jet-ISM interaction, as discussed in this review.

5. Concluding Perspectives

The review outlines the development of numerical simulations of relativistic jets propagating through their environment, with a particular focus on jet-ISM interaction. In addition, the observational implications of such processes have also been outlined. The last few decades has seen a very prominent growth of studies in this domain. This is in contrast to earlier general skepticism (e.g., see the arguments in Section 1 of Ostriker et al. [198]) regarding the impact of jets in the Establishment phase of AGN feedback. However, as discussed in this review, such perceptions are starting to change. Some of the major points that have emerged in this context are summarized below.

• Are RLAGN gas rich? One of the major concerns of the role of jets on their host galaxy was whether radio loud galaxies have enough gas in the first place to be affected by jets. The traditional view has been that in the nearby universe, powerful radio jets are usually found in early type galaxies (ETG), which were considered to be gas poor. However, systematic surveys of such systems have uncovered a significant fraction ($\sim 25\%$[267]) to host dense gas, with higher fractions for radio loud AGN ($\gtrsim 34\%$ [268]). A summary of the various surveys can be found in Table 4 of Tadhunter et al. [268]; also see recent the recent review Ruffa and Davis [269] for more details on molecular gas in local ETG. Interestingly, fractions of radio loud galaxies with molecular gas and the estimated $H_2$ masses, have also been found to increase in with redshift [270], with masses spanning ${10}^7-{10}^{10}{M}_{H2}$. Thus, dense gas is present in a significant fraction of radio galaxies, opening up the potential for jet-driven local feedback, should jets have significant coupling with the ISM.
• Radio detected fraction of AGN? Earlier studies of AGN populations had demonstrated that the radio-detected fraction of AGN reaches up to $\sim 30\%$ for high mass galaxies [271,272]. This, though small, is non-negligible. More recent sensitive radio surveys [273] have extended these to lower radio luminosities and have substantially increased the fractions, to point towards more ubiquitous distribution of nuclear radio activity. Although the higher fractions correspond to radio powers an order of magnitude or more lower than considered earlier, implying weaker nature of the AGN, they nonetheless, provide credence to wide-spread presence of radio activity. Furthermore, it is important to distinguish between the traditional definition of “radio-loudness” inferred from correlations of [OIII] and 1.4 GHz radio luminosity may not often imply radio “silent”. As demonstrated in recent surveys [263,264] a significant fraction of traditional “radio-quiet” sources may harbor nuclear radio emission driven by an AGN, or a jet. They may also demonstrate jet-ISM interaction, as shown in Figure 5 [219].
• Extent of the jet’s impact: Although the apparent beam of radio jets are often found to be thin, collimated structures, the fact that they can have a wider influence, has been well demonstrated by simulations and observations, as outlined in this review. The observed size in radio wavelengths, may often under-represent the large scale impact (e.g., in 4C 31.04 [203], and several other sources in Section B) as the jets are broken into low density streams in the flood-channel phase. However, the even though such impacts may extend beyond the immediate confines of the jet, in most cases, the impact has been seen to be in the central few kpc. Though non-negligible, there needs to be better proof for more wider scale impact, to confirm/discard the predictions from simulations.
• Nature of jet’s impact: outflows, turbulence, starformation: Observed studies have strongly established the presence of jet driven outflows in the central few kpc of several sources, in line with predictions from simulations as outlined in Section 4. The outflows are coincident with broad line widths indicating turbulent motions, demonstrating the jet’s ability to affect the local gas. However, the broader long term implications for such actions in the context of galaxy evolution, especially with regards to starformation, is an open question. Although, several prominent radio loud sources are known to show deficiency in star formation rate [210,225], the ubiquitousness of such cases of jet feedback have been questioned in other recent studies that find less impact on wider scale molecular gas [266]. Even from a theoretical point of view, the theory of turbulence regulated starformation applied to large scale simulations is in its infancy. It should also be noted that any impact from a given jet/AGN feedback episode may not have an instantaneous impact on the SFR, but will have accumulated effect, as stressed in recent reviews [2]. Hence, more systematic studies of long term impact of jets in particular, and AGN in general are needed in future to answer these questions.
• Need for theoretical improvements: As outlined in the review, recent simulations efforts have reached high levels of sophistication and realism in modeling jet-ISM interaction and its impact on galaxy evolution. There however do remain several lacunae that needs improvement. A primary drawback of kpc scale simulations of jet-ISM interaction is the inability to resolve cooling length scales at the outer surface of dense clouds. For example, as outlined in Appendix A of Meenakshi et al. [200], the typical cooling length5 in multi-phase simulations of Mukherjee et al. [169] ranges from $\sim 0.014-1$ pc, below the resolution of the simulations. Achieving such resolutions will require an order of magnitude increase in current resources, which remains a challenging task. Such resolutions are further required also to better understand the shock-cloud interaction, as demonstrated in Figure 2. Such intricate substructures of the cloudlets are not resolved in current simulations.
  Besides the need for better resolutions, most of the simulations in this domain have been carried out for only a few Myr, due to limitations of computational time requirements. However, this explores only a very short phase of the jet and galaxy’s lifetime. Larger scale feedback studies exploring the heating-cooling cycles of jet driven large scale feedback [20,274,275,276] have explored longer run times up to a Gyr. However, they do not resolve the multi-phase gas structures internal to the ISM. Future efforts have to explore at least few tens of Myr of run time, with self-consistent injection of AGN power to account for at least one duty cycle of the AGN. All of these would require larger computational resources, which is expected to become available in near future.
  In addition to the above, new physics based modules need to be incorporated to augment current capabilities of simulations. One such of primary importance is the need for inclusion of the chemistry of ionized and molecular gas phases and other species such as dust. Most numerical codes follow a single fluid prescription, with cooling of matter primary driven by pre-computed tables based on gas densities and temperatures. Very works [208,277] have included more sophisticated treatments of individual fluid elements. In addition to this, impact of photoionizing radiation from the central AGN have been largely unexplored in large scale simulations of AGN feedback, barring a few works [196,199,278,279]. Although well explored for studying cloud dynamics in broad line regions or close to wind lauch zones [e.g., see [280,281,282], and references therein], their effect on larger kpc scale simulations are yet to fully explored.
  Another ill-explored parameter is the effect of magnetic field on shock-cloud dynamics and starformation. Very simulations have included the evolution of magnetic fields [177,178] in simulations of jet-ISM interaction. Magnetic fields can potentially change the nature of shock-cloud interaction by affecting Kelvin-Helmholtz growth rates and also affect estimates of turbulence regulated starformation rates [229], and should be explored in more detail.
  Lastly, another key ingredient overlooked in the current literature is the effect of cosmic rays on the fluid dynamics of jet-ISM interaction in particular, and AGN feedback in general. Active interaction of the jet with dense clouds are expected to be strong sites of production of cosmic rays, due to diffusive shock acceleration at jet-cloud interfaces. Cosmic rays are expected to provide additional momentum and pressure to the fluid, which would in turn affect the local dynamics of the gas. This has been tentatively explored in some cases, e.g., for IC 5063 [283]. Inclusion of cosmic ray diffusion and heating in MHD simulations of galaxy formation is being actively explored by several groups [e.g., [284,285,286]]. However, their impact is yet to be explored in the context of multi-phase AGN feedback.