Numerical Simulations of Jets from Active Galactic Nuclei

Numerical simulations have been playing a crucial role in the understanding of jets from active galactic nuclei (AGN) since the advent of the first theoretical models for the inflation of giant double radio galaxies by continuous injection in the late 1970s. In the almost four decades of numerical jet research, the complexity and physical detail of simulations, based mainly on a hydrodynamical/magneto-hydrodynamical description of the jet plasma, have been increasing with the pace of the advance in theoretical models, computational tools and numerical methods. The present review summarizes the status of the numerical simulations of jets from AGNs, from the formation region in the neighborhood of the supermassive central black hole up to the impact point well beyond the galactic scales. Special attention is paid to discuss the achievements of present simulations in interpreting the phenomenology of jets as well as their current limitations and challenges.


Introduction
Theoretical models of jets from active galactic nuclei (AGN) motivated by observations have been the subject of thorough testing by numerical simulations for almost forty years, since the simulations of supersonic jets performed in the early 1980s by Norman and collaborators [1] confirmed the viability of the beam model proposed by Scheuer [2], and Blandford and Rees [3] one decade before to explain the powering of lobes of extended radio sources and quasars (classic doubles or symmetric doubles, as these sources were known at those dates).
According to our present understanding [4], jets are produced on scales of a few gravitational radii of the central black hole (BH) powering the AGN (R BH ≈ 10 −4 (M BH /10 9 M ) pc, where M BH is the BH mass) but extend to hundreds of kiloparsecs. This disparity in scales is exacerbated when considering those microphysical processes involved in the formation of the jet or the particle acceleration. The viscosity of the accretion disc from which the jet originates is thought to be generated by the magneto-rotational instability [4] affecting the magnetic field in the disc on scales orders of magnitude smaller than R BH . Besides that, the acceleration of particles in the jet is ultimately governed by processes occurring on scales of the gyroradius of mildly relativistic electrons/protons, again several orders of magnitude smaller than the radius of the central BH, for typical jet non-thermal particle energies and magnetic fields. Interestingly, however, the fact that the gyroradius (that establishes the effective collisional mean free paths by suppressing particle diffusion perpendicular to the magnetic fields) and the Debye length (governing the scales of effective charge neutrality) are much smaller than the jet width, R j , supports a continuum (i.e., magneto-hydrodynamical) approximation of the otherwise collisionless plasma forming the jet [5,6].
Incorporating all the relevant microphysics (shock acceleration, magnetic reconnection, radiative processes, etc.) into a global (general relativistic) magneto-hydrodynamical simulation describing the formation of the jet from an accreting BH and its propagation across the interstellar and intergalactic media along hundreds of kiloparsecs represents a daunting task for present-day computational tools and numerical methods. However, numerical simulations have been accumulating important (although partial) triumphs along their decades of development. This review is aimed to summarize them as well as their current limitations and challenges.
The book edited by Böttcher, Harris and Krawczynski [7] is a basic reference covering many theoretical, observational and numerical aspects of the physics of AGN jets of relevance for the present review. A major break-through in the field of Computational Relativistic Astrophysics (and, in particular, in the simulation of jets from AGN) was accomplished in the 1990s when the so-called high-resolution shock-capturing (HRSC) methods were applied to integrate the equations of relativistic hydrodynamics (RHD), and few years later, those of relativistic magneto-hydrodynamics (RMHD) overcoming the traditional difficulties of standard finite-differencing methods in simulating flows with high Lorentz factors. Two recent reviews [8,9] summarize the implementation of these methods in special-relativistic and general-relativistic magnetohydrodynamics (GRMHD) codes and a description of relevant numerical applications in different areas of Relativistic Astrophysics.
The review is organized as follows. In Section 2, we give an overview of the status of the observations and theoretical models of jets that establish the framework for the numerical simulations.
Simulations have traditionally divided the study of the jet phenomenon into separate problems. Section 3 goes into the domain of the simulations of jet formation. Section 4 focuses on simulations of jets at parsec scales up to the smallest scales accessible to observations with present-day instruments. Section 5 analyzes the present status of the simulations of jets at the largest spatial and temporal scales including the connection with the formation and evolution of galaxies and clusters of galaxies. The contribution ends with a summary of the achievements, limitations and challenges of contemporary numerical simulations in Section 6.

Models of Jet Formation
Models proposed to explain the origin of the relativistic jets found in several astrophysical scenarios involve accretion in the form of a disc onto a compact central object. In the case of the jets emanating from AGNs, the central object is a rotating supermassive BH fed by interstellar gas and gas from tidally disrupted stars. There is a general agreement that MHD processes are responsible for the formation, collimation and acceleration up to relativistic speeds of the outflows. In the models of magnetically driven outflows [10][11][12][13], poloidal magnetic fields anchored at the basis of the accretion disc generate a toroidal field component and consequently a poloidal electromagnetic flux of energy (Poynting flux) that accelerates the magnetospheric plasma and plasma from the disc along the poloidal magnetic field lines, converting the Poynting flux into kinetic energy of bulk motion reaching relativistic speeds at extended (≈pc) scales [14][15][16][17]. Energy can also be extracted from rotating BHs via the Blandford-Znajek (BZ) mechanism [18,19]. Several parameters are potentially important for powering the jets: the BH's mass and spin, the accretion rate, the type of accretion disc, the properties of the magnetic field, and the environment of the source [20].
Since in any version of a central engine powering an AGN the accretion disc is magnetized, we expect a wind to be driven from its surface (close to the BH) shrouding the BZ jet, producing transversely stratified outflows in both composition (outer electron-proton wind and inner electron-positron jet) and speed. Furthermore, the disc wind can provide the initial confinement of the jets.
Even at the highest angular resolutions achievable today (tens of micro-arcseconds with space very long baseline interferometry (VLBI) imaging using RadioAstron [21]; 1 µas ≈ 1 (D/20 Mpc)(10 9 M /M BH )R BH , where D is the distance to the source), the scales of jet launching are still inaccessible to observations (although they are the target of the Event Horizon Telescope (EHT) [22], for the M 87 jet −3.9 µas BH). In these conditions, numerical simulations boosted by the development of specific HRSC techniques to solve the GRMHD equations in the last two decades still provide the main guide for further theoretical advances.

Parsec-Scale Jets and Superluminal Radio Sources
At parsec scales, AGN jets, imaged via their synchrotron emission at radio frequencies with VLBI networks, appear to be highly collimated with a bright spot (the core) at one end of the jet and a series of components which separate from the core, sometimes at superluminal speeds [23]. In the standard model of Blandford and Königl [24], these speeds originally predicted by Rees [25] are the consequence of relativistic bulk motion in jets propagating at small angles to the line of sight with Lorentz factors up to 20 or more. Moving components in these jets, usually appearing after outbursts in emission at radio wavelengths, are interpreted in terms of traveling shock waves [24]. The evolution of the continuum spectrum of the ejected components fits nicely with this interpretation [26].
An ongoing, important debate is concerned with the nature of the radio core. Whereas in the standard Blandford and Königl's conical jet model the core corresponds to the location near the BH where the jet becomes optically thin, recent multi-wavelength observations of several sources (e.g., 3C 120 [27], BL Lac [28], 3C 111 [29], and 1803+784 [30]) suggest that the radio core can be a physical feature in the jet (as, e.g., a recollimation shock in  model) placed probably parsecs (i.e., tens of thousands of gravitational radii of the central BH) away from the central engine.
In the acceleration region the width of the jet broadens rather slowly with distance (with a parabolic-like surface [14,16,17]) and the jet collimates. The collimation and acceleration continues until the flow ceases to be Poynting dominated. At this point, the flow has reached its terminal speed and becomes inertial. If the jet is not confined by the external pressure, then the flow becomes conical. Signatures of the acceleration/collimation region can be investigated in those jets where the BH position is matched with the radio core. Parabolic/conical transitions in the jet shape have been observed in M 87 at scales of 10 5 -10 6 R BH [32] and Cyg A at scales ≈ 2.5 × 10 4 R BH [33].
Observations of jets at subparsec scales also show signs of transversal stratification which could be revealing the duality of the jet formation mechanism. The transversal structure could also be the observational counterpart of the growth of different kinds of instabilities triggered at the jet base and/or at the jet/wind or jet/ambient medium interfaces. In 3C 84, the ridge-to-limb brightening change in a scale of years [34] is interpreted as a change in the viewing angle of a jet with a highly relativistic spine and a mildly relativistic sheath. A fast central spine with a slow sheath is also inferred in Cyg A [33] and M 87 [35]. In this last case, the limb brightened structure points towards a subrelativistic layer associated either with an instability pattern speed or an outer wind, and a fast, accelerating stream. Moreover, the systematic difference of the apparent speeds in the northern and southern limbs of the jet provides evidence for jet rotation about the jet axis. The angular velocity of the magnetic field line associated with this rotation suggests that the jet in M 87 is launched in the inner part of the disc, at a distance ≈5 R BH of the central engine. In 3C 273, the double-rail structure has been consistently interpreted as created by helical Kelvin-Helmholtz instabilities [36].
The composition of jets is still a major undecided issue. The presence of protons in the jet can be inferred from a weak component of circular polarization (CP) in the synchrotron emission (which would be exactly zero in the case of a pure pair plasma jet due to the cancellation of the electron and positron separated contributions). However, the CP can be contaminated from Faraday conversion of linear polarization. Recently, VLBI imaging of both circular and linear polarization have been carried out for a few blazars on sub-parsec scales [37,38]. Gabuzda and collaborators [37] detected CP in eight AGN, the most likely origin being Faraday conversion in helical magnetic fields. Homan and collaborators [38] obtained the full polarization spectra of 3C 279 and modeled them by radiative transfer simulations to constrain the magnetic field and particle properties of the pc-scale jet and concluded that the jet is kinetically dominated by electron-proton plasma with a non-negligible contribution of pair plasma to the jet radiation. However, this conclusion is based on an estimate of the low-energy cutoff of the power-law electron distribution, a parameter prone to large uncertainties.
Another way to determine the jet composition is fitting its multiwavelength spectrum. Potter and Cotter [39] developed a multizone inhomogeneous fluid jet emission model to fit with unprecedented accuracy the entire multiwavelength spectrum of a large sample of quiescent blazar spectra with a pure electron-positron plasma trhough their synchrotron, self-Compton and external Compton emission. The model can be viewed as providing an estimate of the power contained in the magnetic and non-thermal electrons in the jet hence establishing a lower limit of the jet power.
Finally, hadronic emission models [40] currently seem to require implausibly large proton kinetic luminosities to match the data (although the recent observation of a highly energetic neutrino from the blazar TXS 0506+056 is best explained by a hadronic emission process [41]).
The quest for the signatures of helical magnetic fields (as those invoked by the jet formation mechanisms) in the inner regions of extragalactic jets represents another observational challenge. Observational signatures for the existence of such helical magnetic fields can be obtained by looking for Faraday rotation measure (RM) gradients across the jet width, produced by the systematic change in the net line-of-sight magnetic field [42] and have been reported for 3C 273 [43,44] and more recently for 3C 454.3 [45] and BL Lac [21]. The paper by Gabuzda in this volume reviews the observational evidences for helical magnetic fields associated with AGN jets.

Kiloparsec-Scale Jets
At kiloparsec scales, jets divide into two main classes following the classification established by Fanaroff and Riley [46] in the 1970s according to their large scale morphology. Fanaroff-Riley type I (FR I) jets are decollimated at kiloparsec scales and display extended lobes of diffuse emission (e.g., 3C 31 [47]). On the contrary, jets from Fanaroff-Riley type II (FR II) objects (e.g., Cyg A [48]) remain highly collimated until the hot-spots, the regions where the jet flow impacts with the ambient medium. Whereas current models [49] interpret FR I morphologies as the result of a smooth deceleration from relativistic to non-relativistic, transonic speeds on kpc scales, flux asymmetries between jets and counter-jets in FR II radio galaxies and quasars indicate that relativistic motion extends up to kpc scales in these sources [50].
Beyond the original divide in power between FR I and FR II sources [46,51], and the differences in their hosts, the widely accepted explanation for the FR I jet deceleration is the entrainment of cold and dense gas through a mixing-layer [52][53][54] or from stars within its volume [55]. Parameterized models assuming intrinsically symmetrical, axisymmetric, relativistic, stationary jets (see [49] and references therein) have been successful in interpreting the brightness distributions of FR I jets characterized by an extended flaring region. However, the ultimate validation/rejection of the different models must rely on long-term (M)HD simulations to account for the non-linear processes leading to the flow deceleration and decollimation.

AGN Jets in the Cosmological Context
Observational evidence is growing that the baryonic part of the low-redshift Universe has been shaped by the energy and momentum output of BHs, through AGN feedback, with profound implications for our understanding of galaxy, group, and cluster evolution [56]. In the so-called kinetic or radio-jet mode, the powerful jets emerging from the AGN push the galactic halo (as inferred from the observed anti-correlation between the radio lobes formed by the jets and the X-ray emission from the cluster gas) and are responsible for inhibiting radiative cooling, which would otherwise lead to unrealistically high star formation rates. The proof of this connection between radio jet heating and radiative cooling suppression is provided by the fact that most cool core clusters (i.e., those with the shortest central cooling times) contain radio lobes [57] and by the direct correlation found between the estimated energy content of the lobes and the cooling rate/luminosity of the intra-cluster medium (ICM) [58].
The morphology, particularly in the closest clusters, of the X-ray holes has led to the interpretation of these cavities as bubbles of relativistic gas [57] blown by the AGN and rising against the cluster's potential well by buoyancy. Whereas this description does apply to the so-called ghost cavities, recent observations of shocks with low Mach number surrounding the lobes in powerful radio sources (e.g., Her A [59], Hydra A [60], MS0735.6+7421 [61], and HCG 62 [62]), tell us that these cavities have not yet reached the buoyancy stage. Hence, although the gross energetics of the AGN/ICM feedback process is roughly understood, the details are not.
On the other hand, the modeling of the feedback processes induced by AGNs is a crucial ingredient to understand the mass distribution of galaxies and their morphologies, the star formation rates and some spatial distributions of the stellar populations (such as the radial gradients of age, metallicity and velocity dispersion).

Jet Formation Mechanisms on the Test Bench
With the advances in the numerical methods in RMHD incorporated into general relativistic codes in the late 1990s, the possibility to explore for the first time the formation mechanism of relativistic jets opened. The first papers considered the problem of jet formation from Schwarzschild and Kerr BHs surrounded by accretion discs. In the case of Schwarzschild black holes [63][64][65], jets were formed via Blandford-Payne's mechanism [12] with a two-layered concentric structure with an inner fast gas pressure-driven jet and an outer slow magnetically-driven outflow both being collimated by the global poloidal magnetic field that penetrates the disc. For Kerr BHs [66], qualitative differences arise between corotating and counter-rotating disc cases. Whereas corotating discs around Kerr BHs produce almost the same kind of outflows as those from Schwarzschild BHs, counter-rotating discs lead to powerful (although still subrelativistic) magnetically-driven jets accelerated by the magnetic field anchored in the ergospheric disc inside the gas pressure-driven jets. Although there is an agreement about the ultimate source of energy for the outflow (the rotational energy of the BH), whether this extraction was the result of a kind of Penrose process [19] that uses the magnetic field to extract rotational energy of the BH at the cost of swallowing plasma with negative energy at infinity, or a purely electromagnetic mechanism, was in dispute [67,68]. The debate now seems to have been settled on the basis of recent theoretical and numerical work [69] which establishes the absorption of negative energy and angular momentum as the necessary and sufficient conditions for arbitrary fields or matter (in the case under study, the magnetic field piercing the BH ergosphere) to extract the BH rotational energy (BZ/Penrose mechanism).
Despite the ground-breaking character of the simulations discussed in the previous paragraph, none of them was able to generate sustained, relativistic outflows. Simulations lasted typically for less than a couple of rotation cycles at the accretion disc's innermost stable circular orbit (ISCO). Besides that, in these simulations, both the accretion discs and the magnetic fields were prescribed in the initial conditions. As a next step, a couple of studies [70,71] focused on the influence of the initial magnetic field configuration around the rotating BH on the outflow properties and considered monopole magnetospheres as in the original BZ mechanism. Koide [70] obtained outflows with Lorentz factors of ∼2 but again the simulation was extremely short (a fraction of a rotation at the ISCO). In a longer simulation [71], the numerical solution evolved towards a stable steady-state solution very close to the force-free solution found by Blandford and Znajek. For the first time, numerical solutions showed the development of a Poynting-flux dominated ultrarelativistic particle wind (Lorentz factor ∼15). The wind was mostly radial and had the largest Lorentz factors on the equatorial plane. Finally, direct numerical simulations of the BZ mechanism were performed [72,73] by solving the time-dependent equations of force-free electrodynamics in a Kerr BH magnetosphere.

Long-Term Simulations of Jet Formation
In parallel, the first simulations of self-consistent jet production (i.e., without assuming a large-scale magnetic field right from the beginning) from accretion discs orbiting Kerr BHs in 2D (axisymmetric) [74] and 3D [75][76][77] were performed. In all the cases, the outflows (formed at the edge of a funnel about 0.5 rad wide around the BH's rotation axis) were sub-relativistic. However, in the axisymmetric case, McKinney [78] succeeded in generating a collimated, long-lived (∼1000 orbital periods at the BH horizon), super-fast magnetosonic, relativistic Poynting-flux dominated jet by tuning the floor model used to refill the evacuated funnel. Interestingly, the flow accelerates along paraboloidal field lines in agreement with analytic models [14,16,17] reaching up to a Lorentz factor ∼10 beyond 10 3 R BH (although still only a factor ∼10 −2 to 10 −3 of the estimated asymptotic value). Komissarov and collaborators [79] revisited McKinney's work with a simplified setup where the low-speed outflow wrapping around the relativistic Poynting-flux dominated jet was replaced by a rigid boundary of a prescribed shape in order to reduce the effects of the numerical dissipation, and fixed injection conditions at the jet base. In all cases, the outcome was a steady state characterized by a spatially extended acceleration region and an efficient transfer of Poynting flux into kinetic energy (much more efficient than that of McKinney's work).
The extended-acceleration scenario raised the question of the stability of the jets against the development of kink instabilities at large distances where the toroidal magnetic field dominates over the poloidal one [78]. Starting with a realistic setup, McKinney and Blandford [80] simulated the generation and propagation of a relativistic highly magnetized jet in 3D and explored both the stability of the jet against the development of the highly-disruptive non-axisymmetric helical kink mode (angular frequency of the radial perturbation, m = 1), and the stability of the jet formation process itself during accretion of dipolar and quadrupolar fields. In their dipolar model, despite strong non-axisymmetric disc turbulence, the jet reaches Lorentz factors of ∼10 with an opening half-angle ∼5 • at 10 3 R BH without significant disruption (see Figure 1).

Figure 1.
Inner ±100 M BH cubical region at t = 4000 M BH of a 3D global GRMHD simulation of the formation of a jet starting with an equilibrium magnetized matter torus, whose angular momentum is aligned with the BH spin. The figure shows the BH, accretion disc (pressure, yellow isosurface), outer disc and wind (log rest-mass density, low green, high orange, and volume rendering), relativistic jet (Lorentz factor Γ ≤ 4, low blue, high red, volume rendering) and magnetic field lines (green) threading BH. (Figure 1 of [80].) A series of papers [81][82][83] concentrate on the process of accretion onto rapidly rotating BHs and the role of magnetic fields in regulating both the accretion process and the jet production. The chosen numerical setup led to magnetically arrested discs [84,85] where the excess in magnetic flux near the BH impedes the accretion of gas hence maximizing the magnetic flux threading the BH per unit of accreted mass. As a result, very powerful outflows are produced, which in the case of rapidly rotating BHs reach efficiencies (power flowing out of the BH over rest-mass energy flux flowing into the BH) greater than 100%, demonstrating the net energy extraction from spinning BHs via the BZ/Penrose mechanism (see Figure 2). Recent numerical simulations along this line [86,87] have examined the interplay among the properties of the magnetized disc flow, the spin of the central BH, and the jet in producing warped discs and precessing jets as well as in controlling the alignment/misalignment properties of the jet with the disc/BH rotation axes. Let us note that the simulations presented in [87] (which qualify as the simulations of BH accretion discs at the highest resolutions to date) were carried out with a 3D-GRMHD code accelerated by GPUs performing 10 times faster than on last-generation multi-core CPUs.

Current-Driven Kink Instability and Magnetic Energy Dissipation
The noticeable large-scale stability of the newly-born fully three-dimensional jets emanating from accreting, rapidly rotating BHs in the first long-term simulations [78,80] prompted a series of studies addressed to re-evaluate the development of 3D instabilities in relativistic magnetized jets under controlled conditions. We concentrate on these studies here and leave for Section 4.2 a more general description of the instabilities arising in relativistic, magnetized flows and their role in the interpretation of the phenomenology in AGN jets at larger scales. Mizuno, Hardee and Nishikawa [88] focused on the stability of magnetized relativistic precessing spine-sheath jets with purely poloidal magnetic fields proving the stabilization effect of weakly relativistic sheaths against Kelvin-Helmholtz instabilities (KHI). In a series of papers, Mizuno and collaborators studied the development of the current-driven kink instability in relativistic force-free jets under idealized conditions (temporal analysis within the periodic box): comoving plasma columns [89]; sheared jets [90]; and differentially rotating jets [91]. The last papers in this series [92,93] have concentrated on the spatial growth of the current-driven instability (CDI) in relativistic jets. On its own, Porth [94] studied the stability of jets from rotating magnetospheres along the initial acceleration and collimation region by means of high-resolution adaptive mesh refinement (AMR) simulations in 3D. His analysis showed that the m = 1-5 modes saturate at a height of ∼20 inner disc radii.
The models of jet acceleration discussed thus far have considered a gradual acceleration of the jet by magnetic forces. However, dissipation via reconnection of alternating magnetic fields [95,96] has been suggested as an alternative energy conversion mechanism. On the other hand, the large-scale magnetic field may dissipate if the regular magnetic structure is destroyed as a result of a global MHD instability, the CD kink instability being the most plausible candidate. O'Neill and collaborators [97] considered the development of CDI in comoving plasma columns in initial radial force balance through various combinations of magnetic, pressure and rotational forces, and examined the resulting flow morphologies and energetics. Models in which the initial magnetic field is force-free deform, but do not become disrupted. Systems that achieve initial equilibrium by balancing pressure gradients and/or rotation against magnetic forces, however, tend to shred, mix and develop turbulence. Consistently with this result, CDI-driven kinetic energy amplification is slower and saturates at a lower value in force-free models. Singh and collaborators [93] also explored the possible link between CD kink instability and magnetic reconnection in their study of the spatial grow of CDI in relativistic jets. It is interesting to note however that, in these papers, based on (ideal) RMHD simulations, none of the solutions in the non-linear regime can be regarded as converged. For convergence, a physical dissipation scale provided by either, e.g., a Navier-Stokes viscosity or Ohmic resistivity (in most astrophysical scenarios, many orders of magnitude smaller than the numerical dissipation scales attainable by magneto-fluid dynamical codes with present-day computational resources) would have to be included in the problem. Additionally, magnetic reconnection can play a major role in regulating the accretion process and hence the production of jets [98][99][100].

General Relativistic Radiative Transfer Simulations
Impelled by the forthcoming EHT observations of Sgr A * and the radio core of M 87, spectra and images of the accretion disc/BH/jet system are currently being performed by post-processing GRMHD simulations using relativistic radiative transfer (for radiative processes such as synchrotron emission, self-absorption, and inverse-Compton processes) and ray-tracing including general relativistic effects ( [101] and references therein; see also Figure 3) after modeling the location-dependent distribution functions of the emitting particles. However, it is important to note that, despite the huge advances experienced by the numerical simulation of the jet formation process, they still suffer from important theoretical uncertainties. A remarkable example is the mass-loading of the jets, caused by the known failure of up-to-date GRMHD codes inside the highly magnetized funnel [102], with important implications in the jet composition issue (see Section 2.2) and the interpretation of the future observations from the EHT [103]. O'Riordan and collaborators [104] investigated the observational signatures of mass loading in the funnel by performing general-relativistic radiative transfer calculations on a range of 3D-GRMHD simulations of accreting BHs by removing the contribution to the spectrum of the artificially supplied floor material, i.e., restricting the analysis to the case where the funnel material is highly magnetized. Conversely, Broderick and Tchekhovskoy [105] considered the creation and acceleration of pairs at the stagnation surface of Poynting-dominated jets and the subsequent inverse Compton cascade as the mechanism for filling the jets with non-thermal particles near the horizon to explain the compact radio emission in M 87.

(Magneto-)Hydrodynamical Simulations
The first simulations of compact jets studied the structure of shocks of steady, relativistic jets propagating through pressure decreasing atmospheres [107,108] and the evolution of relativistically moving perturbations travelling down the underlying steady flow [109][110][111][112]. To account for the relativistic effects (relativistic Doppler boosting, light aberration, and time delays) dominating the emission of parsec-scale jets and compare with observations, calculations of the synchrotron radiation transfer were performed by post-processing the hydrodynamical models. These calculations led to the first synthetic radio images of compact jets reproducing the basic phenomenology of these sources (with steady quasi-periodic as well as superluminal radio components), and the first multifrequency (radio) light curves. In particular, Gómez and collaborators [110] simulated the interaction of standing shocks and relativistic perturbations mimicking the ejection of superluminal components from the VLBI core as in the Daly-Marscher's [31] model, and studied the dependence of the structure of radio components with frequency [113]. Next simulations [114][115][116][117] aimed to explain the complex behavior observed in many sources as, e.g., the dragging of steady components in 3C 279 [118], the presence of trailing components in 3C 120 [113,115] and 3C 111 [119], and the tangled evolution of components in 3C 111 [117]. In a very recent work, Fromm and collaborators [120] analyzed the influence of the thermal absorption of the obscuring torus in the appearance of radio jets.
The preceding simulations were based upon pure (relativistic) hydrodynamical models and the synchrotron emissivity was computed after adding an ad-hoc magnetic field with ordered and random components. Hence, as a natural improvement, subsequent simulations incorporated consistent magnetic fields from RMHD models of jets. Broderick and McKinney [121] used jets formed from rapidly rotating, accreting BHs through GRMHD simulations [80] extrapolated to parsec scales, and analyzed the RM generated in the disc wind wrapping the jet, containing ordered toroidally dominated magnetic fields. Porth and collaborators [122] studied the synchrotron emission from RMHD jets with large-scale helical magnetic fields along the acceleration and collimation region. Both simulations [121,122] reproduce the expected gradients in RM across the jet width due to the toroidal component of the helical magnetic field [42]. In addition, Porth et al. [122] provided detailed insights regarding the polarization structure throughout the jet, which depends strongly on the helical magnetic field pitch angle, jet viewing angle, Lorentz factor, and opacity.
Several papers have concentrated on the observational signatures of recollimation shocks in relativistic, magnetized, overpressured jets, commonly associated with the stationary components often seen in parsec-scale VLBI observations of AGN jets [21,[123][124][125]. Mizuno and collaborators [126] considered kinematically-dominated jets with axial, toroidal, and helical force-free magnetic fields to investigate the effects of different magnetic field topologies and strengths on the recollimation structures. Martí, Perucho and Gómez [127] extended the previous study to jet models with (non-force-free) helical magnetic fields and different energy compositions (internal, kinetic, and magnetic) including models with high magnetizations. The total-pressure mismatch between the jet and the ambient medium combines with the radial Lorentz force to generate jets with a complex superposition of periodical recollimation shocks and gentle expansions and compressions along the jet, and steep transversal gradients. In [128], the authors modeled the optically thin total and linearly polarized synchrotron emission from a selection of models of [127] (see Figure 4), with the ultimate goal to connect the properties of the magneto-hydrodynamical jets with the structures observed in extragalactic jets at parsec scales. In particular, these simulations show a top-down emission asymmetry along the jet produced by the helical magnetic field and a noticeable spine brightening for highly magnetized jets. The bright stationary components associated with the recollimation shocks (whose separation and strength depend on the magneto-hydrodynamical parameters) present relative intensities modulated by the Doppler boosting ratio between the pre-shock and post-shock states. Models at small viewing angles display a roughly bimodal distribution in the polarization angle, due to the helical structure of the magnetic field.
Complementing the above studies on overpressured jets propagating through homogeneous ambient media, the large-scale structure of initially free-expanding jets through pressure-decreasing atmospheres, p a ∝ z −κ , has been tackled using various approximate analytical approaches and numerical simulations for both purely hydrodynamical and magnetized relativistic jets (see, e.g., [129,130] and references therein). Values of κ 1 − 2 produce reconfinement shocks reaching the jet axis, whereas values of κ ≥ 4 keep the jets expanding freely. Most importantly, κ governs the global stability of the jets: beyond κ = 2, the jet suffers such a rapid lateral expansion that the causal communication across the jet is completely lost and hence global instabilities of any type become totally suppressed [131]. Fromm and collaborators [132] have computed synthetic radio maps of hot, mildly magnetized jets propagating down a pressure-decreasing atmosphere with κ = 1. Most interestingly, these authors incorporated the characteristics of the observing array (VLBA in this case) and of the observing experiment, and the imaging algorithm (Difmap's [133] CLEAN deconvolution [134]) in the computation of the radio maps.

Instabilities
Magneto-hydrodynamic flows can exhibit KH as well as CD instabilities. Whereas KHI grow in the contact discontinuity or shear layer of flows with relative speeds (and can develop in both magnetized and unmagnetized flows), the growth of CDI modes is driven by the distribution of the electric density current flowing along the jet in magnetized models. Strictly speaking, it is not possible to unambiguously separate their effects in a magnetized jet, although, in general, CDI prevail in magnetically dominated jets, whereas kinetically dominated jets are subject to KHI. In the case of CDI, the stability conditions can be predicted theoretically in the force-free approximation. In kinetically dominated jets, the KHI conditions follow from the linearized (R)MHD equations. Hardee [135] provided a condensed review of both KHI and CDI including references to papers on linear stability analyses (see also [8]). The analysis of the stability against these types of perturbations allows one to constrain the space of flow parameters leading to a successful jet launch and propagation. On the other hand, instability modes can grow to finite amplitudes without destroying the jet and leave observational imprints (in the form of radio knots, transverse structure, bends, etc.) whose properties can be used to pinpoint the jet flow parameters. This analysis has been successfully applied to probe the physical conditions in the jets of several sources (e.g., M 87 [136][137][138], S5 0836+710 [139][140][141][142], 3C 273 [36,143], 3C 120 [144][145][146], and Bl Lac [147]).
Beyond the linear regime, the analysis requires numerical (hydrodynamic or magnetohydrodynamic) simulations. Here, the main purpose is to assess the stability, collimation, and mass entrainment properties of jets at large (temporal and spatial) scales. We restrict ourselves to the relativistic case. In a series of papers, Perucho and collaborators studied the effects of relativistic dynamics and thermodynamics on the development of KH instabilities in relativistic slab jets-both in the vortex-sheet approximation [148,149] and for sheared flows [150]-covering the linear, saturation, and non-linear phase of the evolution by means of hydrodynamic simulations. In [150,151], the authors considered the effects of very high-order reflection modes on sheared relativistic (both kinematic and thermodynamic) slab jets, and in full 3D for initially cylindrical jets, in [152]. Other authors [153][154][155] have focused on the properties of the fully-developed turbulence flow driven by the KHI and check its consistency (velocity power spectrum, flow intermittency) with the classical Kolmogorov's model in the inertial range in which the turbulent eddies interact without external forcing or dissipation.
In the context of the magnetic KHI, Mizuno, Hardee and Nishikawa [88] focused on the stability of magnetized relativistic precessing spine-sheath jets with purely poloidal magnetic fields (see Section 3.3), whereas other papers [156,157] have been devoted to elucidate the basic properties (dynamo amplification of magnetic fields, power spectra) of the KHI-driven RMHD turbulence. In particular, Beckwith and Stone [157] presented preliminary results on the turbulent amplification of the magnetic field in the cocoon formed along the propagation of a kpc-scale jet. Let us note that, instead of starting from a shearing flow, Zrake and MacFadyen [158] studied the development of RMHD turbulence from initial spatially-uniform conditions.
The non-linear development of the CDI have concentrated on the most disruptive mode, the kink mode. The state of the art of these numerical studies as well as of the CDI-driven RMHD turbulence are summarized in Section 3.3 in the context of the magnetically-dominated flows arising at the jet formation region.
The KHI and CDI are by far the best studied instabilities in the context of (relativistic) jets. However, in the last years, other fluid instabilities have deserved the attention of theorists and numerical simulators in the context of reconfining jets. These instabilities are the Rayleigh-Taylor (RTI) and the centrifugal instabilities (CI). The radial force arising from the pressure mismatch between the jet and the surrounding ambient medium in the expanding/compressing jet flow is the driving force of the RTI [159][160][161]. In the case of the CI considered in [162], the driving force is the centrifugal force experienced by the fluid parcels moving along the poloidally curved trajectories in the course of the reconfinement process. The transverse structure of the jet is remarkably deformed by the non-linear growth of these radial oscillation-induced instabilities, which can evolve into turbulence leading to the complete jet disruption. These instabilities could be behind the division of AGN jets (although see Section 5.2) into two morphological types in the Fanaroff-Riley classification (see Figure 5). It is also worth noting that the CI is also present in the rotating flows forming the jets near the central engine [163][164][165]. Further work is needed to untangle the competition between the KHI and CI in the context of the FR division [162], and to analyze the effect of including strong magnetic fields which may inhibit the growth of CI and KHI modes and promote CDI [166]. The results correspond to simulation C1 in [162]. The jet to ambient density ratio is 1.6 × 10 −2 and the initial jet Lorentz factor is 5. The jet kinetic power is 2 × 10 44 ergs/s. (Figure adapted from [162] kindly provided by the authors.)

Microphysics
The simulations discussed previously rely on a pure macroscopic (magneto-)fluid dynamical modelization of the jet that accounts only for the evolution of the thermal plasma. However, the emission in these sources is produced by non-thermal particles (NTPs) accelerated at shocks [167] or in magnetic reconnection events [168] that require a microscopic description whose implementation in macroscopic fluid dynamical models is a challenge for present-day scientific computing.
Current models of blazar jet emission (e.g., [39] and references therein) focus on the microscopic processes (acceleration of NTPs, emission processes, and adiabatic losses) and a simplified jet structure model (in this case, a 1D time-independent relativistic fluid flow with a variable shape and bulk Lorentz factor) with the goal of fitting simultaneous multifrequency observations. Marscher [169] presented a parameterized model for variability of the flux and polarization of blazars in which turbulent plasma flowing at a relativistic speed down a jet crosses a standing conical shock representing the mm-wavelength core of these sources.
On the opposite side, RMHD simulations produce self-consistent models for the underlying jet and implement simplified recipes to describe the NTPs. Several papers [109][110][111]122,128,170] set the energy density of the NTPs proportional to the pressure of the thermal plasma. Since the pressure is a good tracer of shocks, this approach mimics the emission if the process of shock acceleration were dominant. Alternatively, other works [121,122,128] consider the energy density of the NTPs as proportional to the magnetic energy density more suitable to account for the particle acceleration in the case of magnetic reconnection. Anantua and collaborators [171] suggested some other similar recipes where the energy density of the NTPs is assumed proportional to a power of the magnetic energy density.
An alternative approach in the numerical modeling of non-thermal emission from astrophysical jets treats the population of NTPs as a separate population advected by the fluid. Transport of NTP populations in classical jets in radio galaxies were carried out in [172][173][174]. Mimica and collaborators [175] improved the approach of Gómez et al. [109,110] by incorporating a transport algorithm for the advection of the NTPs including synchrotron losses to study the spatial and temporal variations of the spectral index of the radio components. The same numerical treatment of the NTP transport and the synchrotron emission is used in [176] to interpret the light-curve and the VLBI observations of the 2006 flare of CTA 102 within the shock-shock interaction model devised in [110]. Vaidya and collaborators [177] incorporated a particle module in an RMHD code that includes particle acceleration at shocks and synchrotron and (external) inverse Compton emission.
Particle-in-cell (PIC) methods appear as the most realistic way to simulate the kinetic dynamics of the plasma forming the jet. PIC codes model plasmas as a collection of charged macroparticles (representing many physical particles) that are moved by integration of the Lorentz force. Currents associated with the macroparticles are deposited on a grid on which Maxwell's equations are discretized. Electromagnetic fields are then advanced via Maxwell's equations, with particle currents as the source term. Finally, the updated fields are extrapolated to the particle locations and used for the computation of the Lorentz force, so the loop is closed self-consistently. To ensure that kinetic effects are resolved in the simulation, it is necessary that the grid spacing be much smaller than the plasma skin depth, c/ω p , and that the timestep be much smaller than the corresponding timescale where ω p is the plasma frequency (ω p ≈ 5.3 × 10 5 / √ n e cm; n e is the electron density). PIC methods are described in this volume by Nishikawa. This approach is capable of treating all effects present in collisionless plasmas, including particle acceleration at shocks [178] and through magnetic reconnection [179] and incorporate them into global simulations of relativistic jets [180,181]. Finally, there are efforts to combine the dynamical evolution of the thermal and non-thermal components via MHD-PIC equations [182][183][184] (see Figure 6). In this approach, it is not necessary to resolve the plasma skin depth and it is enough to resolve the much larger Larmor (gyration) radius.

Morphology and Dynamics of kpc-Jets
The first numerical simulations of kiloparsec-scale jets aimed at probing the standard model of FR II radio sources [5]. Among them, the numerical simulations of supersonic (non-relativistic) jets propagating through a homogeneous atmosphere performed by Norman and collaborators [1,185] concentrated on the identification of the basic structural elements (supersonic beam, cocoon, terminal shock, bow shock) and the connection with the morphological elements in powerful radio sources (radio jet, radio lobes, hot spots). These early simulations also served to reveal the mechanism of the jet's working surface governing the jet propagation. The review paper by Burns, Norman and Clarke [186] summarizes the status of the numerical simulations of extragalactic radio sources up to 1990, comprising the first simulations of magnetized jets, the first synthetic images of kiloparsec-scale jets based on synchrotron emission, and the first 3D simulations.
The advent of relativistic codes based on HRSC techniques triggered the study of the effects of relativistic flow speeds and/or relativistic internal energies in the morphology and dynamics of jets in 2D [187][188][189][190][191][192][193] and 3D [194][195][196] simulations. These exploratory relativistic simulations soon incorporated the effects of magnetic fields. The first simulations focused on the propagation of relativistic jets with pure poloidal magnetic fields through ambient media with aligned [197,198] and oblique [199,200] magnetic fields to study how the fields affect the bending properties of relativistic jets. Later studies of large-scale RMHD jets explored the dependence of morphological and dynamic properties of jets on the magnetic field configuration ( [201][202][203]: toroidal fields; [204]: toroidal and poloidal fields), and on the ratios of magnetic energy density and thermal pressure, and magnetic energy density and rest-mass energy density.

Large-Scale Simulations of kpc-Jets and the FR I/FR II Dichotomy
The simulations discussed in the previous section concentrated on the characterization of the morphology and dynamics of jets in terms of the jet-to-ambient density ratio, jet Mach number, relativistic effects, magnetic field configuration and strength, or injection opening angle. However, none of those simulations have considered the evolution of powerful radio sources as a whole.
The long-term propagation of a powerful jet (FR II jet) through the ambient medium generates a characteristic and well understood structure, formed by: (i) a terminal or reverse shock at the head of the jet where the flow decelerates and heats (the hot-spot); (ii) a hot and light region shrouding the jet (the cocoon), inflated by the shocked jet particles; and (iii) a dense shell of shocked ambient medium. Begelman and Cioffi [205] described the expansion of this structure assuming a constant jet propagation speed, a complete and instantaneous conversion of the injected jet energy into internal energy in the cocoon and a sideways expansion mediated by a strong shock. Guided by observations showing that cocoons of powerful sources with very different physical size have similar axial ratios, several authors [206][207][208] developed self-similar expansion models by tuning parameters such as the initial jet opening angle and the ambient density gradient, ignored in the Begelman-Cioffi model.
The first (non-relativistic) long-term simulations of jets aimed to decipher the extended emission of powerful radio sources. Several authors [209][210][211] used postprocessed bremsstrahlung emission maps to explain the X-ray observations showing evidences of interaction between the radio emitting plasma and the X-ray emitting gas in the ambient medium (see Section 2.4). Other papers [172,174,212] concentrated on the interpretation of the emission at radio wavelengths from MHD simulations of radio galaxies including a detailed treatment of the microphysics: relativistic electron transport, diffusive acceleration at shocks and (synchrotron, inverse Compton from CMB photons) radiative and adiabatic cooling. More recent simulations [213][214][215] focused on the structure of the magnetic field in the radio lobes and the properties of the synchrotron linearly-polarized emission.
Among the first relativistic simulations of powerful radio sources, are those of Scheck and collaborators [216] aimed at looking for clues on the jet composition (e ± , e/p plasmas) in the large-scale morphology of (2D, axisymmetric, purely hydrodynamical) jets, and the first high-resolution (20 zones per jet radius) 3D simulations of relativistic magnetized jets presented by Mignone and collaborators [217].
Many papers based on simulations of powerful radio sources have focused on the consequences of the interaction of the jet with the interstellar and intergalactic media in an on-going task to comprehend the AGN feedback phenomenon at galactic and supra-galactic scales. These papers are described briefly in the next section. First, let us examine the contribution of numerical simulations of jets to the solution of the FR I/FR II dichotomy issue. As stated in Section 2.3, the ultimate process leading to the morphological division between FR I and FR II sources is the deceleration of the flow to transonic speeds in the low power sources at kpc scales. Numerical simulations have explored different mechanisms responsible of this deceleration. One of these mechanisms is the jet mass loading by stellar winds [55]. However, axisymmetric numerical simulations [218] show that only weak (L j ≈ 10 41 -10 42 erg s −1 ) FR I jets can be decelerated by mass entrainment rates consistent with present models of stellar mass loss in elliptical galaxies on scales of 1 kpc. Effective deceleration mechanisms for more powerful jets (L j ≥ 10 43 erg s −1 ), involve the development of instabilities and/or strong recollimation shocks. Rossi and collaborators [219] discussed the deceleration of FR I jets resulting from the growth of KH helical instabilities developing at the jet surface, whereas Perucho and collaborators [152] studied the properties of the KHI-driven entrainment in the case of sheared jets. However, none of these papers deal with the origin of the instabilities, which were imposed as initial/boundary conditions. Perucho and Martí [220] connected the jet deceleration with the entrainment of colder and denser ambient gas within the jet due to the growth of KHI pinching modes downstream of a strong recollimation shock. Krause and collaborators [221] studied the properties of this recollimation shock (and the resulting flows) as a function of the initial jet opening angle. A recent work [162] (see Section 4.2) proposes the development of the CI [165] in the outer layers of the reconfining flow passing a recollimation shock as the trigger of the turbulent entrainment of the jet and its subsequent deceleration. Relying on 3D, purely hydrodynamical, non-relativistic simulations of supersonic jets, Massaglia and collaborators [222] studied the morphologies of FR I jets. Interestingly, for low enough kinetic powers (L j ≤ 10 43 erg s −1 ), the energy of jets propagating down an inhomogeneous ambient medium, instead of being deposited at the terminal shock, is gradually dissipated by turbulence producing the plumes characteristic of FR I objects.
The simulations discussed in the previous paragraph ignore the effects of magnetic fields in the FR division. In the model of Tchekhovskoy and Bromberg [223], jets with powers below some critical value depending on the galaxy core mass and radius and the magnetic field strength can become kink-unstable within the core, stall, and inflate cavities filled with relativistically hot plasma leading to FR I morphologies (see Figure 7). Finally, let us note that a recent work [224] has considered the impact of the environment (galaxy cluster-like, poor group-like) in the jet morphology. Although focused on FR II sources, these and future simulations of this kind could shed light on how the observable properties of the FR I and FR II jets change due to jet-environment interaction.

AGN Feedback
The impact of AGN jets on their galactic scale environment was realized when the first X-ray observations of galaxies and clusters of galaxies showed a clear anticorrelation of the radio-and X-ray emitting plasmas [225]. At the same time, it became apparent that the energy deposition of the jets in the galactic environment could be an important agent in solving the cooling flow problem which manifests itself on a variety of scales, from isolated elliptical galaxies to large clusters of galaxies [226].
As mentioned in the previous section, some works [209][210][211] aimed at reproducing the X-ray cavities in powerful radio sources. Besides these active cavities associated to ongoing jet injection, we also find the so-called ghost cavities, which appear disconnected from the active jet-fed cocoon. This second type of cavities are associated with past periods of jet activity and/or low power (FR I) jets. The (axisymmetric) simulations of Churazov and collaborators [227] described the dynamics of buoyant bubbles inflated by an earlier phase of nuclear activity of the galaxy to explain the complex morphology of the X-ray and radio maps in the central 50 kpc region around the galaxy M 87. Quilis and collaborators [228] studied the effect of the energy injection on the cooling rate of a cluster core by means of 3D simulations of the evolution (rise, fall, and mixing) of a hot bubble inflated close to the core center. Their results confirmed this mechanism as a very efficient one for regulating cooling flows. Unlike these simulations, where the lobe inflation was modeled artificially by adding a low-density hot bubble [227] or injecting internal energy [228] within a region close to the gravitational center, Reynolds and collaborators [229] followed the evolution of cavities of powerful jets (L j ≥ 10 45 erg s −1 ) across the initial supersonic phase up to their fate as buoyantly rising plumes long after the jet activity has ceased. Let us note by passing that methods similar to that of Quilis and collaborators [228] (i.e., the injection of thermal energy into bubbles placed around the central black hole) have been subsequently implemented into full cosmological simulations [230] and used for example in the Illustris project [231], a series of large-scale hydrodynamical simulations of galaxy formation. Such feedback is a necessary ingredient in cosmological simulations to reproduce the properties of galaxy populations realistically. A more sophisticated AGN jet feedback method [232,233] which deposits mass, momentum and energy has been developed recently.
The early simulations discussed in the previous paragraph focused on the (hydro)dynamics of the cavity/bubble expansion in the cluster medium. However, the results by Bîrzan and collaborators [234] on a systematic study of radio-induced X-ray cavities in 18 systems show that the energy associated with the cavity inflation is about 10-50% the energy necessary to quench the cooling. The conclusion is that a major part of the AGN energy must be transferred to the ICM in other ways not related with the simple hydrodynamic expansion of the bubbles. Options include dissipation of weak shocks and sound waves, turbulence, thermal conduction, shock-heating, heating by cosmic rays, or magnetic fields [235]. Moreover, among the open problems remaining is also to find a robust mechanism by which the (intrinsically directional) AGN jets succeed in heating the ICM isotropically. In the following lines, we concentrate on a few representative papers facing this problem. The 3D-AMR simulations performed by Yang and Reynolds [236] set the state of the art of purely hydrodynamical, non-relativistic simulations including radiative cooling (no viscosity, no heat conduction). The (time-dependent) BH accretion/bipolar jet ejection is modeled via a simplified subgrid model. According to these authors, the isotropization of the heating is the result of a sustained circulation of ambient gas originally displaced by the jets towards the cluster center. Magneto-hydrodynamical simulations performed by the same authors [237] show that anisotropic (along the magnetic field lines) conductive heating is likely significant only for the most massive clusters. Martizzi and collaborators [238] performed (non-magnetized) simulations of jet heating similar to those of Yang and Reynolds and analyzed the dependence of the results on the diffusion of the numerical algorithms, specially the Riemann solver. The work reflects the difficulty to achieve a proper numerical convergence for the heating/cooling problem derived from the disparity in scales of the processes involved. Finally, recent papers (e.g., [239,240]; this last reference using a moving-mesh code) have simulated the heating by cosmic rays as a mechanism for the AGN energy thermalization in clusters. Supplied by the AGN jets, cosmic rays (likely protons) disperse in the magnetized ICM via streaming, and interact with the ICM via hadronic, Coulomb, and streaming instability heating. Moreover, jets energetically dominated by cosmic rays lose momentum more quickly causing the jet to expand laterally and to displace the ICM close to the center of the clsuter more isotropically.
The papers cited above consider static, idealized cluster cores. Heinz and collaborators [241,242] were the first to take into account the dynamic nature of the cluster gas and detailed cluster physics in their simulations of the interaction of AGN jets with galaxy clusters. The conclusion is that in a cosmologically evolved cluster, the motions of cluster gas effectively distribute the effects of the AGN over a wide angle. Using a more relaxed cluster, Mendygral and collaborators [243] performed MHD simulations confirming the distortion of the morphology of jets and lobes due to the ICM weather. Finally, the recent work by Bourne and Sijacki [244] (using a moving-mesh code) shows that the large-scale turbulence generated by orbiting substructures result in line-of-sight velocities and velocity dispersions consistent with the Hitomi observations of the Perseus cluster [245].
Following a different line of work, Perucho and collaborators [246][247][248] examined the properties of the cluster core heating by relativistic jets. The ultimate motivation behind this study is that modeling the energy input of a powerful jet within a non-relativistic framework leads to inconsistent sets of injection parameters. On one hand, the modeled jets are unrealistically slow (or unphysical). On the other hand, without the contributions of the flow Lorentz factor and the relativistic enthalpy to the inertia, jets are unrealistically massive [249]. The results from 2D axisymmetric simulations [246,247] show that the cavities of powerful relativistic jets are more symmetrical (less elongated) than their non-relativistic counterparts and produce more powerful and long-living bow shocks, and hence are less prone to break into buoyant bubbles, alleviating the problem of the heat isotropization (see Figure 8). Additionally, since the jet particles store much less energy in the form of rest-mass, a larger fraction of energy can be transferred (and, in fact, it is) to the ICM. These basic conclusions seems to survive in 3D simulations [248]. Three-dimensional simulations of lobe inflation in cluster environments by mildly-relativistic magnetized jets below equipartition have been recently performed by English and collaborators [250]. As important as unraveling the role of the jet-driven AGN feedback in the solution of the cooling flow process is to determine its influence in the evolution of the host galaxy, in particular through the galaxy's star formation history [251,252]. Gaibler and collaborators [253] performed 3D non-relativistic AMR simulations of the interaction of a powerful AGN jet with the massive gaseous disc of a high-redshift galaxy to asses the impact on the star formation. In [254], special attention is paid to the description of the galaxy ISM described as a two-phase fractal medium with different maximum cloud sizes and volume filling factors. The feedback efficiency, measured by the amount of cloud dispersal generated by the jet-ISM interactions, is investigated through 3D relativistic AMR simulations (see Figure 9).

Summary and Conclusions
The present review offers a perspective of the contribution of numerical simulations to the understanding of the AGN jet phenomenon along the last four decades of research in the field. Given the disparate scales involved, present simulations of jets concentrate on partial aspects of the problem. These can be broadly divided into the jet formation process (BH accretion, plasma acceleration and collimation), the understanding of the phenomenology of jets at parsec and sub-parsec scales, and the interplay between the jet and the ISM/ICM with important implications in the jet morphology, and the galaxy and cluster evolution (see Section 2). Since these three scenarios also require different physical ingredients, we have reviewed them separately, focusing in their achievements.
Over time, the complexity and physical detail of the simulations have been increasing with the pace of the advance in theoretical models, computational tools and numerical methods. The development of suitable numerical techniques for the (G)RMHD equations (mainly the HRSC methods; Section 1) brought about a revolution in the description of the thermal component of jets (Sections 3.1-3.3, 4.1, 4.2 and 5). Nevertheless, present simulations now demand incorporating microphysical processes (shock acceleration, magnetic reconnection, and radiative processes) in a realistic (i.e., consistent) way. In this regard, the recent attempts to combine (R)MHD and PIC methods into a single hybrid code (Section 4.3) seem very promising.
The AMR finite-volume/finite-difference methods and the more recent moving-mesh techniques (with great impact in the modeling of AGN feedback; Section 5.3) represent elegant and effective ways to improve the numerical resolution in critical regions of the flow, in combination with high-order finite differencing methods. Their use must extend to other areas of the jet modeling.
Finally, as important as the development of suitable numerical techniques is the increase in computing power. In this respect, the use of GPU-accelerated codes (Section 3.2) may represent a significant jump in computing speed.