Backflows by AGN jets: Global properties and influence on SMBH accretion

TL;DR

This study investigates how AGN jet-induced backflows influence the gas dynamics and accretion processes near supermassive black holes, revealing their intermittent nature and potential to significantly modulate jet power.

Contribution

It provides the first detailed numerical and analytical analysis of AGN jet backflows, linking their properties to jet parameters and their impact on SMBH accretion and jet variability.

Findings

Backflows are spatially coherent and intermittent regardless of jet power.

Backflows can enhance SMBH accretion rates by up to ten times during episodes.

Backflows influence the density profile of the circumnuclear disk, creating a flat core.

Abstract

Jets from Active Galactic Nuclei (AGN) inflate large cavities in the hot gas environment around galaxies and galaxy clusters. The large-scale gas circulation promoted within such cavities by the jet itself gives rise to backflows that propagate back to the center of the jet-cocoon system, spanning all the physical scales relevant for the AGN. Using an Adaptive Mesh Refinement code, we study these backflows through a series of numerical experiments, aiming at understanding how their global properties depend on jet parameters. We are able to characterize their mass flux down to a scale of a few kiloparsecs to about $0.5\,\mathrm{M_\odot/y}$for as long as $15$ or $20$ Myr, depending on jet power. We find that backflows are both spatially coherent and temporally textbf{intermittent}, independently of jet power in the range $10^{43-45}$ erg/s. Using the mass flux thus measured, we model analytically the effect of backflows on the central accretion region, where a Magnetically Arrested Disk lies at the center of a thin circumnuclear disk. Backflow accretion onto the disk modifies its density profile, producing a flat core and tail. We use this analytic model to predict how accretion beyond the BH magnetopause is modified, and thus how the jet power is temporally modulated. Under the assumption that the magnetic flux stays frozen in the accreting matter, and that the jets are always launched via the Blandford-Znajek (1977) mechanism, we find that backflows are capable of boosting the jet power up to tenfold during relatively short time episodes (a few Myr).

Figures

Tables

Table 1. Table 1: Parameters, timings and bubble characteristics. All simulation parameters: run specifications (name, smallest cell side, simulation time), halo parameters (mass, central cooling time), jet parameters (power of each of the two jets, Mach number, lifetime, mass and momentum injection fluxes) and total backflow gas mass at the given time (i.e. the total mass in the backflow region isolated in Figures LABEL:fig:large-scaleEC to LABEL:fig:RC ).

Simulation			Halo		Jet					Backflowing mass
										(at given time)
Name	Resolution	$t_{\max }$	$M_{200}$	$t_{cool,0}$	$P_{\mathrm{jet}}$	${ℳ}_{\mathrm{jet}}$	$\Delta t_{\mathrm{jet}}$	$\dot{m}$	$\dot{p}$	$Total$	$central$
	[pc]	[Myr]	[$M_{\odot }$]	[yr]	[erg/s]		[Myr]	[$M_{\odot }$ /yr]	[$M_{\odot }$ /yr $km/s$]	[$M_{\odot }$]	[$M_{\odot }$]
Elongated Cavity series
EC42	78.125	473	$1.7\times {10}^{[}12]$	$6\times {10}^{[}8]$	${10}^{42}$	5	79	0.0088	167.22	$4.84\times {10}^{[}5]$	$1.28\times {10}^{[}4]$
										20 Myr	20 Myr
EC43	78.125	140	$1.7\times {10}^{[}12]$	$6\times {10}^{[}8]$	${10}^{43}$	5	42	0.0190	776.15	$4.69\times {10}^{[}5]$	$7.11\times {10}^{[}3]$
										10 Myr	10 Myr
EC44	78.125	115	$1.7\times {10}^{[}12]$	$6\times {10}^{[}8]$	${10}^{44}$	5	21	0.0409	3603	$9.92\times {10}^{[}5]$	$1.9\times {10}^{[}4]$
										7 Myr	7 Myr
Round Cavity series
RC44	78.125	23.1	$2.6\times {10}^{[}12]$	$4\times {10}^{[}8]$	$1.12\times {10}^{[}44]$	5	5.37	0.0237	2900	$6.90\times {10}^{[}4]$	$1.04\times {10}^{[}5]$
										10 Myr	10 Myr
RC45	78.125	22.2	$2.6\times {10}^{[}12]$	$4\times {10}^{[}8]$	$1.12\times {10}^{[}45]$	5	5.37	0.0510	13461	$4.84\times {10}^{[}4]$	$6.80\times {10}^{[}4]$
										8 Myr	8 Myr
RC46	78.125	22.2	$2.6\times {10}^{[}12]$	$4\times {10}^{[}8]$	$1.12\times {10}^{[}46]$	5	5.37	0.1098	62548	$2.71\times {10}^{[}5]$	$4.34\times {10}^{[}4]$
										7 Myr	7 Myr

Table 2. Table 2: Parameters of the backflow mass profile model: the linear combination is shown as the template (black continuous line) in Figure 6 .

i	${\dot{\sigma }}_0$	q	${\tau }_0$	${\tau }_1$	ts
	[$M_{\odot }/p{c}^{2}y$]		[Myr]	[Myr]	[Myr]
1	$0.14$	$5$	$0.25$	$0.2$	$5$
2	$3\times {10}^{-4}$	$3$	$3.5\times {10}^{-2}$	$0.4$	$6.5$
3	$2.5\times {10}^{-5}$	$3$	$1.5\times {10}^{-2}$	$0.4$	$8$
4	$3.5\times {10}^{-5}$	$2$	$1.5\times {10}^{-2}$	$2$	$10$

Equations

$$\frac{\partial \mathbf{v}}{\partial t}-\mathbf{v}\times \nabla \times \mathbf{v}=-\nabla h+T\nabla S$$

$${\rho }_g(r)={\rho }_0\exp \left(-\mu m_p\frac{\left[\Phi (r)-\Phi (0)\right]}{k_B{T}_0}\right)$$

$$r_c=\frac{L^2}{GM\left(r_c\right)}$$

$$\delta r$$

$$l$$

$$t_{ex}\approx 6.34\times 10^5{L}_{44}^{-1}{R}_s^{3/2}{M}_8^{1/2}\left(15/\ln \left(\Lambda \right)\right)T_9^{3/2}\mathrm{yr}$$

$$t_d$$

$$p_{gas}$$

$$p_{kin}\equiv \rho v^2$$

$$\frac{\partial \Sigma }{\partial t}+v_r\frac{\partial \Sigma }{\partial r}+\frac{\Sigma }{r}\frac{\partial }{\partial r}\left(r{v}_r\right)={\left[-2{\rho }_{bf}{v}_{bf\mid z}\right]}_{\mid z=H/2}\equiv A(t)$$

$$v_r=-\frac{\beta A(t)}{r}$$

$${\Sigma }_0\left(r\right)={\sigma }_0\frac{r_0}{r}.$$

$$\Sigma (r,t)={\int }_0^{t}d\tau A\left(\tau \right)+{\Sigma }_0\left(\sqrt[]{r^2+2\beta {\int }_0^{t}d\tau A\left(\tau \right)}\right)$$

$$A_i(t)={\dot{\sigma }}_0{\left(\frac{t-t_s}{{\tau }_0}\right)}^q\exp \left(-\frac{t-t_s}{{\tau }_1}\right)\vartheta \left(t_s\right)$$

$$\Delta =\frac{\Delta T_9}{\Sigma }\frac{\partial \Sigma }{\partial T_9}$$

$$\frac{B_m^2}{8\pi }=\frac{G{M}_{BH}{\Sigma }_m}{r_m^2}$$

$$B(r)=B_m{\left(\frac{r_m}{r}\right)}^{1/2}$$

$${\Phi }_m={\int }_{r_g}^{r_m}dr2\pi rB(r)\backsimeq \frac{4\pi }{3}{B}_m{r}_m^2$$

$${\Phi }_m={\left(\frac{128{\pi }^{2}G}{9}\right)}^{1/2}{M}_{BH}^{1/2}{r}_m{\Sigma }_m^{1/2}$$

$$r_m={\left(\frac{\eta }{2\pi }\right)}^2{G}^{-1}\frac{M_{BH}{\dot{m}}^2}{{ϵ}^2{\Sigma }_m^2}$$

$${\Phi }_m=\kappa G^{-1/2}{M}_{BH}^{3/2}{\dot{m}}^2{ϵ}^{-2}{\Sigma }_m^{-3/2}$$

$$\dot{m}=\frac{\dot{M}}{{\dot{M}}_{Edd}}=\frac{4\pi r_m^2\rho (r_m)v_r}{\eta M_{BH}}$$

$$\Phi \propto M_{BH}^{-1/2}{ϵ}^{-2}{\Sigma }_m^{1/2}{A}^2(t)$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0$$

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla \cdot \left(\rho \mathbf{vv}\right)+\nabla P=\nabla \left(\rho \mathbf{g}\right)$$

$$\frac{\partial \rho E}{\partial t}+\nabla \cdot \left[\left(\rho E+P\right)\mathbf{v}\right]=\rho \mathbf{v}\cdot \mathbf{g}-\Lambda \left(\rho ,T\right)$$

Full text

Using an Adaptive Mesh Refinement code, we study these backflows through a series of numerical experiments, aiming at understanding how their global properties depend on jet parameters. We are able to characterize their mass flux down to a scale of a few kiloparsecs to about $0.5$ ${\rm M}_{\odot}\,{\rm yr}^{-1}$ for as long as $15$ or $20$ Myr, depending on jet power. We find that backflows are both spatially coherent and temporally intermittent, independently of jet power in the range $10^{43-45}$ erg/s.

Using the mass flux thus measured, we model analytically the effect of backflows on the central accretion region, where a Magnetically Arrested Disc lies at the centre of a thin circumnuclear disc. Backflow accretion onto the disc modifies its density profile, producing a flat core and tail.

We use this analytic model to predict how accretion beyond the BH magnetopause is modified, and thus how the jet power is temporally modulated. Under the assumption that the magnetic flux stays frozen in the accreting matter, and that the jets are always launched via the Blandford-Znajek (1977) mechanism, we find that backflows are capable of boosting the jet power up to tenfold during relatively short time episodes (a few Myr).

††pagerange: Backflows by AGN jets: Global properties and influence on SMBH accretion–A††pubyear: 2016

1 Introduction: backflow morphology and AGN jet self-regulation

The propagation of AGN jets inflates large, hot, turbulent cavities in the interstellar medium of their host galaxies. Circulation of gas in such cavities gives rise to pronounced streams of hot gas flowing back from the hot spot (if present, as in FRII radio galaxies), along the cavity boundaries to the central plane.

Such backflows are driven by the thermodynamics of the gas, and — once in the central plane — consist of very low angular momentum gas, which potentially reaches down to very small scales, contributing to the mass and energy supply in the accretion region around the SMBH. Backflows carry very hot, high pressure gas; they can thus heavily affect circumnuclear star formation and the properties of the accretion disc, as a self-regulating feedback mechanism. Backflows as a feature of jet-cocoon systems were already noticed in the first simulations of the propagation of relativistic jets into homogeneous atmospheres (Norman et al., 1982), and confirmed in more recent simulations (Rossi et al., 2008; Perucho & Martí, 2007). Mizuta et al. (2010) distinguish between two types of backflows, according to the different geometries of the flow itself: a straight backflow, with flow lines extending from the tip of the hotspot back to the origin, and a bent backflow, where the flow lines acquire curvature. In early 2D simulations, precursors of what we present in this work, Antonuccio-Delogu & Silk (2010) also noticed the formation of this feature, and that the backflow evolved from a bent to a straight geometry. In that work, backflow was described as large-scale vorticity created by sharp gradients in the thermodynamic state of the gas at the hot spot and cavity boundaries, precisely as stated by a fundamental theorem of fluid dynamics, known as Crocco’s theorem (Crocco, 1937). This can be understood from the Euler momentum equation:

[TABLE]

Here S is the entropy and $h=U+p/\rho+v^{2}/2\,$ is the stagnation enthalpy (Cap, 2001; Shu, 1992). Even for a stationary flow, Crocco’s theorem states that vorticity can only be created by finite gradients of enthalpy $h$ and/or entropy $S$.

Antonuccio-Delogu & Silk (2010) pointed to the connection between backflows with large-scale vorticity in the cavity. The flow begins near the hot spot (HS), where a curved shock front induces a jump in entropy and a gradient in the Bernoulli constant transverse to the shock. The downstream gas thus gains a vorticity (Shu, 1992), and its flow is then confined between the dense and hot bow shock from the outer side, and the hot turbulent cavity gas from the inner side.

This goes on until the gas falls back to the central plane and follows the cavity edge (or collides with mirror backflows in a bipolar jet) falling down towards the jet origin with very low impact parameter (and thus angular momentum), although its inflow velocities reach up to several hundreds km s*-1*.

In three dimensions, however, this mechanism loses some effectiveness, as with the additional degree of freedom the velocity could be directed (in absence of other constraints) anywhere in the $z=0$ plane. Also, the flow is subject to more effective hydrodynamic instability, which can slow and disrupt it.

In previous 3D simulations, Cielo et al. (2014) showed that despite the unarguably reduced efficiency, substantial backflows (always around $1$ ${\rm M}_{\odot}\,{\rm yr}^{-1}$ ) reach the central few hundred parsecs. The duration of such backflows varies with jet power (higher powers move cavities away from the centre at earlier times, killing backflows), but always encompasses a few Myr. Furthermore, the backflow gas was found to be stable against hydrodynamics instabilities, although the simulations covered just the first few million years.

Observations of backflows have been quite challenging for a long time, as the gas is hot but very sparse, and only mildly relativistic, so easily out shined by the jets. However, observational characterization of backflows has recently been emerging; in particular Laing & Bridle (2012) observed backflows in two low-luminosity jetted radio-galaxies. In particular, a mildly backflowing component around the lobes is needed in order to fit the emission and polarization radio maps.

AGN are multi-scale systems, and as the backflows get to closer to the central BH, they experience all the relevant physics at different scales (Antonuccio-Delogu & Silk, 2007): after the thermodynamics-dominated circulation in the cavities, they will collide with a central structure (likely, a circumnuclear-nuclear disc extended up to about $100$ pc), where they may generate further inflows.

After this stage, this secondary inflow may enter the magnetosphere of the BH, where the dominating energy source is the magnetic field, ultimately responsible of launching the jets.

Throughout this work, we aim to follow backflows from start to end. This is important when we consider their effect as a source of hot gas accretion onto the central BH, capable of triggering further jet activity or increasing the power of an already running jet in a self-regulating context.

We start from the largest scales (kpc or larger) by making use of the hydrodynamic simulations, described in Section 2 and interpret the results on the basis of Crocco’s theorem. For this purpose, we will provide visual snapshots of the density field, as well as histograms of the gas spatial distribution along the z-axis (Section 3).

Next, we proceed by investigating, with similar methods, the flow of gas which has already reached the $z=0$ plane. In order to quantify the impact parameter of the gas at this stage, we plot the evolution of its circularization radius and calculate the mass flux within $2$ kpc (Section 4).

We then explore the effects of this infall onto the circumnuclear region, not resolved in our simulations, assuming that the innermost magnetized structure is a Magnetically Arrested Disc (Section 5).

Finally, in Section 6 we develop some quantitative considerations of how the processes described above influence the rate of mass accretion onto the central SMBH and the production of the jet itself (via the Blandford-Znajek mechanism from Blandford & Znajek, 1977).

2 Setup and simulation description

We run our simulation using the hydrodynamic, Adaptive Mesh Refinement (AMR) code FLASH v. 4.2 (Fryxell et al., 2000). In our computational setup, we solve the non-relativistic Euler equations for an ideal gas with specific heat ratio $\gamma=5/3$ (see Appendix A for more details). In order to properly model the energetic of the system, we include a static, spherical gravitational potential of the host dark matter halo (following a NFW profile), as well as the self-gravity of the gas. This set-up is an evolution of that used in Cielo et al. (2014) — henceforth C14.

We choose a cubic 3D computational domain and Cartesian coordinates. The side of the box is in all runs fixed to $L=640$ kpc, much larger than the maximum extent of the jets, in order to encompass most of the halo, whose typical size $R_{200}$ is $\sim 250$ kpc. The resolution, defined as the size of the smallest computational element, is in all cases fixed to $\Delta r=78.125\,{\rm pc}$ (see Appendix A). All the boundary conditions are set to the FLASH outflow (i.e. zero-gradient) default value. We present two families of runs, mainly differing in jet power $P_{jet}$. In both series, the dark matter halo density follows a spherical NFW profile with $R_{\mathrm{200}}=0.25$ Mpc, $M_{\mathrm{200}}=1.71\times 10^{12}\,M_{\odot}$ and a concentration parameter $c_{\mathrm{200}}=7.8$ (Dutton & Macciò, 2014). We then fill this potential well with hot coronal gas, having a uniform temperature $T_{0}$ and metallicity ([Fe/H] = -1.0) in hydrostatic equilibrium within the dark matter potential. We achieve the latter condition by adopting the following constant temperature profile, introduced by Capelo et al. (2010):

[TABLE]

where $\mu$ is the mean molecular weight of the gas, $\Phi(r)$ the gravitational potential and $k_{B}$ the Boltzmann constant. For the normalization, we follow McCarthy et al. (2008) in setting the ratio of gas to dark matter mass within the halo $R_{\mathrm{500}}$ to 0.85 times the cosmic baryon fraction (taken in turn from Komatsu et al. 2011). Once the normalization of the gas profile is fixed, the gas properties are thus completely specified by the chosen $T_{0}$, or alternatively by the gas central cooling time, which we report in Table 1 (column 5).

The first family is the EC series (for Elongated Cavities), where the DM halo has a mass fixed to 1.7\times 10^{[}12]$${\rm M}_{\odot} and a central cooling time of $6\times 10^{[}8]$ yrs. For this series, we launch three runs, differing only in the jet mechanical power $P_{jet}$, which assumes values of $10^{42}$, $10^{43}$ and $10^{44}$ erg/s (run EC42, EC43 and EC44, respectively). The second series (denoted the RC series, for Round Cavities), differs from the EC series mostly in that it features a higher halo mass, a slightly shorter cooling time (a consequence of a denser hot gas phase) and higher values of $P_{jet}$ (runs RC44, RC45 and RC46; see Table 1 for a complete list of all parameters).

The dependence of the cavity’s shape on the jet/halo physical parameters was previously investigated in C14. Briefly, the cavity shape is linked to these parameters via the jet’s thermal pressure: a jet having higher thermal pressure (or propagating in a colder environment) will originate rounder cavities as the over-pressure determines isothermal expansion; on the contrary a cold jet with a small ratio between internal and kinetic pressure will inflate elongated cavities. This being said, in our case the RC jets create rounder cavities than their EC counterparts for three reasons:

• •

they are intrinsically hotter, as they are more powerful (i.e. faster) but at injection they have all the same internal Mach number $\mathcal{M}_{jet}=5$ (in other words, the Mach number relative to the environment changes).

• •

the core of the halo in which they propagate is slightly colder;

• •

they thermalize a larger part of their total energy, as more powerful jets have shorter lifetimes, and thermalization is more efficient at early stages.

The jets are modelled by injecting gas in the central region for a prescribed lifetime. The injection power, mass and momentum flux are also reported in Table 1 for each run (columns 6, 9 and 10, respectively); a given kinetic power also corresponds to a specific lifetime $\Delta t_{jet}$ (column 8 of Table 1). In the EC series, backflows on the central plane significantly fade before the end of the jet’s active phase; in the RC series, however, there are residual backflows even after the jets are switched off. In some cases the injection velocities resulting from these parameters choices are slightly superluminal. The ram pressure of the local environment where the jet emerges, however, brings the jet to non-relativistic velocities in the first few cells111We had also tried increasing the jet’s cross section and lowered the velocity to keep the injected power constant, but the results are indistinguishable after sufficiently long times.. The overall jet/Hot Spot advancement velocity is however much slower than that (maximum about $1.5\times{10}^{4}$ km/s, Figure 1), so this high nominal injection speed does not affect the evolution of the cavities.

The jets advance through the Interstellar Medium (ISM), initially producing a hot, localized shock in the Hot Spot (HS) and a bow-shock region, accordingly to the previous findings of C14. In Figure 1 we show the location $z_{max}$ of the bow-shock region, simply defined as the geometrical extent along the jet axis of the hot spot; this gives an idea of the actual advance speed of the jets and of their effect on the hot gas.

The jets expand up to a maximum z-distance ranging from a few tens up to about $100$ kpc in the EC runs (with a clear turnover when the jets are switched off), while they reach about $20$ kpc in the RC ones, mainly due to the shorter simulation time. In any case, in this work ,we focus on the first few tens of Myr, as backflows arise at these times.

3 Galactic-scale backflows around cavities

The prediction of Crocco’s theorem that backflows originate from the HS and then flow along the lobe/bubble boundary is verified in the simulations run, as can be seen by looking at Figures LABEL:fig:large-scaleEC and LABEL:fig:large-scaleRC222Some animations of the simulations are available from https://blackerc.wordpress.com/people/salvatore-cielo/, where we show visual slices of the density field along the $y=0$ plane (left columns). The backflows follow a different colour scheme (colour legend on the right), in order to highlight them within their context in the cocoon; the backflow region is also recognizable as it is the only one where the velocity field is superimposed333See the colour legend at the bottom of each panel in the Figure for the magnitudes of the arrows..

We define the backflow regions to include all grid cells whose velocities point towards the jet origin within a $\pm 45$ degrees cone.

This (rather conservative) selection is necessary to view a cleaner flow, as backflows are otherwise contaminated by gas patches that ‘bounce” on the cavity walls. While this is contemplated in our Crocco theorem description, the resulting velocity field shows many spurious features. The constraints we can put on the mass flows with this analysis are for this reason only lower limits.

We set an additional threshold on the radial (w.r.t. the origin) velocity: $v_{r}\leq-225$ km/s. This is to exclude coherent cooling flows (or occasional highly-turbulent spots) at late times; this selection does not significantly change the backflow mass at the epochs considered in this study.

Cooling flows however can develop after about $50$ or $100$ Myr, i.e. later than the time intervals during which the jets are on. A discussion of the cooling flows is outside the scope of this paper; as we see from the slices shown, backflows reach high velocities (often $\leqslant-2500$ km/s) so this selection does not exclude any significant part of the backflows even at late times.

Finally, we cut out the innermost $\pm 2$ kpc of gas from the cavity selection (Figure LABEL:fig:large-scaleEC and LABEL:fig:large-scaleRC), which will be the object of Section 4 (as one can see from Figure 4 and 5).

The backflows initially appear as thin layers contained within the cavity/dense-shell interface (left panels in figures LABEL:fig:large-scaleEC and LABEL:fig:large-scaleRC); thus in 3D these flow layers wrap the entire inner cavity. Since the cavities at this stage reproduce the lobes observed in radio-galaxies (see C14), backflows can appear around most $10-20$ kpc galactic radio-lobes. After $10-20$ Myr, the cocoon develops an internal structure: the lobes detach from the central plane, and leave a gap filled by denser gas. Following C14, we call this the lobe phase. During this phase, turbulence develops as a consequence of shocks and shearing between the different gas layers, which creates turbulent eddies through Kelvin-Helmholtz (KH) instabilities.

The large-scale backflows are affected by this structure, and can converge back to the jet axis following the bubble boundaries (right panels in Figure LABEL:fig:large-scaleEC and LABEL:fig:large-scaleRC). In their path, they also take part in the cocoon’s turbulent motion, both near the HS and along the shearing cavity boundary; they also contribute to generate the shear, as they initially consist of laminar flows in relative motion with respect to both the inner cavities and the outer bow-shock.

These aspects were analysed by C14 and the flows were found to be stable against KH instability444In C14 the resolution was slightly better than in the present work, however the simulations lasted only $\sim 6$ Myr or less.. Regardless of how much they contribute to the generation of turbulent motions, the backflows are perturbed and fragmented by it : one can clearly see patches of coherent inward radial velocity in Figure LABEL:fig:large-scaleRC, more prominent at later times (panels on the right) and within the innermost $\sim 5$ kpc, where the cavities start to detach from the centre.

The backflows can gain vorticity and momentum at two different sites: first at the HS, where they start their journey around the lobes, and later near the $z=0$ plane, since after the lobe phase backflows must bend again, this time following the jet beam chimney that connects the lobe to the jet origin.

Near the central region, the backflow follows instead a rather straight pattern. The mass of the gas involved in the backflows (Table 1, column 12), as well as the time it takes to get back to the central plane depend on the cavity shape and size, and on the hotspot pressure (enthalpy) which gives the initial kick.

After a sufficiently long time, the HS are usually too far away from the centre, and the backflows stop halfway. A few Myr after the jet has been switched off ($t\geq\Delta t_{jet}$) , the lobes turn into roughly spherical bubbles and detach completely from the centre. We refer to this stage as the bubble phase. At this time, cooling flows can occur near the $z=0$ plane. Any residual backflows will then cease; however analogous circulation patterns will persist in the bubble for all its lifetime (bubbles from light, supersonic jets such as these create vortex-ring-shaped cavities; see e.g. Guo, 2016).

In the EC runs, large-scale backflows are more extended and carry more mass; although the jet power is on average lesser, it can drive the gas around the cavities efficiently (also because jet lifetimes are correspondingly longer).

In the RC runs instead, the presence of more spherical cavities forces streamlines to gain more curvature since the start, until they reach the $z=0$ plane. Also the flow appears more fragmented in the RC case, as more disconnected patches are clearly visible in the slices. This is probably due to the increased turbulence in the cocoon environment generated by hotter and more powerful jets. Such patches linger for up to about $10$ Myr near or within the central plane.

In order to estimate the backflow location and mass transport, in Figure LABEL:fig:large-scaleEC and LABEL:fig:large-scaleRC, we add density histograms of the backflowing gas distribution along the z-axis, for the same snapshots shown in the slices; the total mass in the region for the first snapshot is also reported in Table 1. In many EC runs, one sees that initially the distribution extends back to the first central few kpc and then moves farther away at later times. This is less true for the RC runs, although overall the mass involved is smaller by a factor $5$ or $10$.

As for the total gas mass accumulated at the centre, it seems to be approximately constant (about $5\times 10^{5}\ \mathrm{M}_{\odot}$ in the EC case, $5\times 10^{4}\ \mathrm{M}_{\odot}$ in the RC case), except for the most powerful (and shortest-lived) jet events, in which it increases by a large factor (about $2$ and $5$, respectively).

4 Kiloparsec-scale backflows on the central disc

We now turn our attention to the central region. In Figures 4 and 5 ,we plot density slices along the $z=0$ plane, again with backflowing gas highlighted, and with velocity arrows superimposed.

The backflow in the central plane is more regular during the first few Myr, and more patchy afterwards, participating in the turbulence of the cocoon gas that affects the entire cavity. In this case, the mass transport is significant even for the low power jets. There is more mass involved in central disc backflows in the RC runs than in the EC series (see Table 1, last column) notably different from what is seen for the cavity-wide backflows described in Section 3.

Due to axial symmetry, we expect that the backflowing gas in the central disc should have on average little to zero angular momentum, and thus flow directly towards the BH accretion disc.

In order to test this prediction, we need to trace the angular momentum of the gas. In Figures 4 and 5 (right panels), we plot the evolution of the circularization radius — a proxy for angular momentum — within a small central cylindrical selection.

Let $L$ be the modulus of the specific angular momentum vector of a gas parcel in a computational cell. We define the circularization radius $r_{c}$ as the radius the parcel would have if it were on a circular orbit within its host Dark Matter halo:

[TABLE]

Here $M(r)$ is the mass of dark matter555We neglect self-gravity of the gas, although it is accounted for in the simulations within $r$. The evaluation of equation 3 is particularly simple as our dark matter halo is spherically symmetric, thus it is straightforward to find the implicit solution of equation 3. In using the full modulus of the 3D angular momentum of the gas in each cell $L$, rather than its z-component only ($L_{z}\leq L$), we are making a conservative estimate.

We evaluate $r_{c}$ for all the gas cells of the backflowing gas only, selected with the same velocity threshold as in Section 3. This time though we select a cylinder centred on the jet origin and having $2$ kpc radius in the xy plane, and a thickness of about $300$ pc (four simulation cells in total) along the z axis. The mass-weighted average value of $r_{c}$ is plotted for each snapshot time (thick lines in the plots of Figures 4 and 5) together with the minimum and maximum values in the selected region (shaded area around the lines) and a $2$ kpc line (in black). In general, gas having $r_{c}<r$ will always tend to migrate to smaller distances, thus falling towards the central BH. In this case, we can conclude that all gas parcels having $r_{c}<2$ kpc will always stay within the selection.

In almost all runs, the average of $r_{c}$ stays above the $2$ kpc line: we interpret this circumstance as arising from the fact that backflows have quite high characteristic velocities that do not necessarily have a negligible impact parameter, thus resulting in some floor values for $L$ and $r_{c}$. A noteworthy exception is run EC42, in which backflows reach the central region after $10$ Myr, while later the average $r_{c}$ stays always well below $2$ kpc. On the contrary, in both EC43 and EC44 they grow smoothly up to a relatively early peak at $5$ kpc around $7$ Myr then decline back to $2$ or $3$ kpc.

From the figures we can draw up two general conclusions:

• •

There is always backflow with $r_{c}<2$ kpc, as the shaded area always extends down to almost zero.

• •

Statistically, there is always a significant mass fraction able to migrate to smaller radii at all times.

About the second point, although flow masses are not indicated in the figures, we see from Table 1 that the total values are around $10^{4}$ $\mathrm{M_{\odot}}$ or a few times that. Note also that angular momentum is not necessarily conserved in the small-scale flows. It could be dissipated through viscosity, or again through the thermodynamic action described by Crocco’s theorem, as patches or streams of gas from the the opposite sides of the bipolar jet collide in the $z=0$ plane, as in the 2D-simulations by Antonuccio-Delogu & Silk, 2010; in support of this, the average $r_{c}$ clearly decreases with time. Thus the values in Figures 4 and 5 are just upper limits on RC.

In Figure 6 we show the mass flux666In our notation, a positive flux means inflow through the same cylindrical region used to estimate $r_{c}$.

We consider only the net flux through the side surface of the cylinder, as in most cases the flux through the bases is dominated by the jet outflow. This selection may miss some of the bent backflows at the lobe base, thus also in this respect we are just putting lower limits to the mass flux.

The flows in Figure 6 start with negative values (i.e. there is net outflow, not plotted) due to the predominance of the jet outflow, which still pushes the initial dense environment gas outwards for the first $5-10$ Myr.

All RC runs present several flux peaks between $5$ and $15$ Myr, reaching about $1\,\mathrm{M_{\odot}/y}$; these peaks originated from the patchy nature of the backflows, but in the RC44 and RC46 cases, we can see a background flux of about $0.5\,\mathrm{M_{\odot}/y}$ until $15$ or $20$ Myr. As noted in C14, similar values could indeed provide substantial central gas accretion, which could contribute to establish a jet self-regulation mechanism; this is the subject of Section 5.

In comparison, backflows in the EC case involve masses smaller by a factor of a few (around $0.25\,\mathrm{M_{\odot}/y}$. They also tend to peak at later times ($10$ to $30$ Myr), possibly because of their different morphology which makes the backflow gas traverse a larger distance before approaching the central region.

5 Parsec-scale backflows and accretion disc kinematics

In order to estimate the impact of the backflow on the accretion region, in this section we present a model of the backflow-accretion disc interaction, which extends the analysis to scales too small to be reproduced in our numerical experiments.

In this model, the central accretion region is assumed to host a magnetically arrested disc (hereafter MAD , see Figure 7). The MAD occupies a small region around the centre of a circumnuclear-nuclear disc in the $z=0$ plane, around which the backflows accumulate.

However, some further analysis is required to prove that backflows are capable of interacting with the MAD , as the backflows we observe in the numerical experiments presented above have typically much smaller densities and much larger velocities than in a typical MAD . Our hot and sparse backflows have thermodynamic properties similar to those of winds, thus one may ask whether they will effectively accrete onto a disc rather than simply flowing past the $z=0$ plane.

In order to answer this question we will adopt an up-to-date model of the MAD modelled after recent GRMHD simulations, where a slim disc (Abramowicz et al., 1988) is threaded by an external magnetic field (Bisnovatyi-Kogan & Ruzmaikin, 1974, 1976; Narayan et al., 2003) (see Tchekhovskoy, 2015, for a recent review about MAD ).

The MAD is characterized by a magnetopause, defined as the region where magnetic and disc thermal pressures are comparable. In our case the magnetopause can be modelled as a sphere (Narayan et al., 2003) of radius $r_{mp}\simeq$ ${\rm r}_{acc}$ , where ${\rm r}_{acc}$ $\simeq 2GM_{BH}/\left(v_{a}^{2}+c_{s}^{2}\right)$, is the black hole’s accretion radius and $v_{a}$ and $c_{s}$ are the Alfvén and sound speed, respectively.

The backflow represents an additional flux component within the MAD . Arons & Lea (1976) have shown that a wind reaching the magnetopause is likely to be affected by the exchange instability, which modifies the morphology of the flow by creating knot-like structures aligned along the magnetic field lines. The transverse and longitudinal sizes of these knots are given by:

[TABLE]

where $M_{bh},R,$ and $m_{0}$ are the Black Hole mass, an independent distance variable, and the mode number of the exchange instability777For our estimate, we can safely take $m_{0}=1$, respectively. The characteristic time scale of this instability is instead given by:

[TABLE]

Here $L_{44}$ is the BH luminosity in units of $10^{44}$ erg s*-1*, $R_{s}$ is the standing shock’s distance888In the Arons and Lea model the infalling gas is hypersonic and forms a standing shock at a distance $R_{s}$. in pc, $M_{8}$ is the BH mass (in units of $10^{8}\rm{M}_{\sun}$), $\Lambda$ is the logarithmic Coulomb factor and $T_{9}$ $\sim[0.01,\,1]$ is the backflow temperature in units of $10^{9}\,$ K. We can easily check that $t_{ex}\ll t_{ff}$, thus at the magnetopause the exchange instability is effective on the backflow: blobs of typical sizes $\delta r,l$ are created and will fall towards the $z=0$ plane, flowing in between the poloidal magnetic field lines. The typical mass of such blobs will be: $m_{ex}\approx\rho_{b}\pi l^{2}\delta r$.

Near the central plane the MAD model predicts an almost azimuthal magnetic field, and the backflow blobs, being fully ionized, will become diamagnetic, a feature which will shield them from the further action of the magnetic fields and ease their survival along their path after having crossed the magnetosphere. Moreover, they will also feel a drag force with a characteristic time scale: (King, 1993; Vietri & Stella, 1998):

[TABLE]

Here $v_{A}$ and $\mu\simeq 0.62$ are the Alfvén velocity and mean molecular weight of a fully ionized hydrogen plasma, respectively, and: $n_{-2}=n/10^{-2}\,cm^{-3}$ . $B$ is measured in gauss and $l$ in parsec. In deriving Equation (5), we have assumed that the blobs can be modelled as cylinders of aspect ratio $\delta r/l=0.1$. We thus see that the backflow will be accreted onto the meridional disc on a very short time-scale. As is shown in Fig. 1 of Vietri & Stella (1998) this drag force is directed in the opposite direction to the blob magnetic field, and thus it will tend to further drive the blob away from the magnetic field lines.

We can now proceed to analyse the properties of the backflow-accreting disc. At variance with a standard thin disc, the accreting gas from the backflow is very hot ($T\gtrsim 10^{8-9}$ K), thus it will be fully ionized and will carry a frozen-in magnetic field. Taking typical values for the backflows from the numerical experiments presented in the previous sections, the internal and kinetic pressures of gas accreting on the disc are:

[TABLE]

(in cgs units), where $n_{-2},v_{3}\,$ and $T_{9}$ are the gas density, velocity and temperature in units of $10^{-2}\,{\rm cm}^{-3}\,$, $10^{3}\,{\rm km/sec}\,$ and $10^{9}\,$ K, respectively. Thus the internal and kinetic pressures are comparable, and the accreted disc will be hot even under the very conservative hypothesis that no fraction of the kinetic pressure would be dissipated and converted into internal energy.

The continuity equation for the surface density $\Sigma(r)=\int dz\rho$ in cylindrical coordinates can be written as:

[TABLE]

where the source term in the right-hand side denotes the mass flux contributed by the backflow. We now assume that the radial inflow velocity is given by:

[TABLE]

This velocity profile is predicted by magnetized Keplerian disc models (Kaburaki, 1986), which also predict an equilibrium surface density profile having the same radial dependence $\Sigma\propto r^{-1}$. Note that this assumption amounts to reducing the time dependence of the radial velocity to that of $A(t)$, the mass flux rate, implying then an instantaneous rearrangement of the accretion flow in the whole disc. We now assume that an initial low–mass, magnetized and Keplerian thin disc exists before backflow infall (i.e. for t$\leq 5\,$ Myr in our experiments), with an initial density profile as given by Kaburaki:

[TABLE]

Under these conditions, it is easy to show (Antonuccio-Delogu et al, in preparation) that eq. 6 admits an exact solution:

[TABLE]

The first term on the r.h.s. is always increasing with time: it describes homogeneous physical accretion, independent of position. On the other hand, the second term, for a monotonically decreasing initial density profile such as that considered here, is decreasing at any $r$. Physically, this is a consequence of the presence of an inflow (finite $\beta$) which tends to subtract mass, driving it towards the BH. The temporal evolution of the density profile will thus locally depend on the combination of these two opposing terms.

In Figure 6, we present the mass variation within a cylindrical surface of radius $r=2\,$ kpc and height 0.1 times the radius. We restrict ourselves only to positive values, i.e. inflow towards the central BH contributed from the backflow. We approximate the total backflow flux as the composition of few accretion episodes. The black curve in this figure presents a template model obtained by combining four profiles of the form:

[TABLE]

where $\vartheta$ is the Heaviside step function. The parameters of this fit are given in Table 2.

To determine the normalization constant $\beta$ we assume that the accretion radial inflow in the $z=0$ plane is everywhere subsonic: $v_{r}\leq c_{s}$, and that it reaches the sound speed precisely at the Innermost Stable Circular Orbit.

Using Eq. (7) we get: $\beta\leq rc_{s}A^{-1}_{max}$, where: $c_{s}=\left(kT/\mu m_{p}\right)^{1/2}$ is the sound speed and $A_{max}=max_{t\leq t_{max}}\left[A(t)\right]$, where we have set $t_{max}=20\,Myr$ (typical accretion duration from Figure 7). Assuming a solar composition plasma ($\mu^{-1}=1.613$) we finally get: $\beta=2.925\times 10^{-3}\left(T_{9}/\mu\right)^{1/2}$ pc4 M${}_{\odot}^{-1}$, where $T_{9}=T/10^{9}$ K is the only free parameter characterizing the backflow, within this model.

In Figure 8 we show the evolution of the density profile for two different values of $T_{9}$. In general, backflow accretion breaks the self-similarity of the density profile by creating a flattened tail and a central core. This latter feature arises as a result of the competition between accretion from the backflow and an increase of the radial inflow towards the BH (eq. 7). At large distances, the flattened density profile is a consequence of the convergence of the integral of $A(t)$ for $t\rightarrow\infty$ to a value which dominates the argument of the second term. This erases any memory of the initial density profile.

How robust is the solution given above, and in particular its asymptotic state? We have performed a sensitivity analysis w.r.t. $T_{9}$. In the insert in Figure 8 (right) we plot the total sensitivity (Dickinson & Gelinas, 1976, eqs. 12 and 17) which, for a single parameter, reduces to the parameter:

[TABLE]

The total sensitivity gives a measure of the relative variation of $\Sigma$ for a given relative variation of $T_{9}$. As we see from Figure 8, the most sensible relative variations are limited to the innermost few $10^{-2}$ pc, i.e. the region occupied by (or closest to) the central MAD : for larger radii, the final density profile is almost insensitive to variations in the temperature of the backflow, which is directly related to the accretion rate $\dot{m}$. These properties are physically consistent with a further, interesting feature of the exact solution given above (eq. 6): the source term $A(t)$ appears in the solution (9) only under integral, i.e. as the total accreted mass from the backflow. Thus the final state of the disc ($t\rightarrow\infty$) will be independent of the time history of backflow.

6 BH accretion and implications for jet power

In the magnetically arrested disc (MAD ) model there is a tight relationship between accretion and jet production, the latter being promoted through the Blandford-Znajek mechanism (hereafter BZ Blandford & Znajek, 1977). The details of this process involve physical mechanisms such as magnetic reconnection which are not quantitatively understood (see Tchekhovskoy, 2015; Punsly, 2015 for recent reviews). There is however general agreement that in the BZ mechanism the jet mechanical power $\rm{P}_{jet}$ scales as $\rm{P}_{jet}\propto\Phi_{m}^{2}\omega_{d}^{2}$, where $\Phi_{m}$ is the magnetic flux within the magnetosphere and $\omega_{d}$ is a function only of the black hole’s spin $a$ (see Tchekhovskoy, 2015, sect. 3.4).

In this section ,we adopt the mass inflow from the circumnuclear disc models described in Section 4 to estimate the resulting magnetic flux $\Phi_{m}$ seen by the BH magnetosphere. This calculation comes with a caveat: while the physical scale of interest for a MAD is of order $\simeq 100$ Schwarzschild radii, our circumnuclear disc model does not include MHD physics relevant on these scales, as it is instead smoothed down to $r=0$. Thus on intermediate scales the physical processes we do not consider here may introduce time delays in the accretion (arising from the MHD properties of the disc), or deviate part of the backflow mass into other forms of outflows. We will investigate the effects of some of these features in future, more systematical models of circumnuclear disc/MAD accretion (Antonuccio-Delogu et al. 2016, in prep.).

This first simple, exactly solvable model will give us a first-order estimate about how backflow accretion modifies $\rm{P}_{jet}$. In this section we will simply assume that backflow accretion does not modify the BH mass and spin, and only affects the magnetic flux within the BH magnetosphere. This is a reasonable assumption, as accretion of hot, sparse gas should have little effect on the growth of BH properties.

The magnetosphere of the BH is defined as the region where magnetic pressure and gravitational attraction are comparable:

[TABLE]

where $\Sigma_{m}\equiv\Sigma(r_{m})$. The size of the magnetic field threading the BH-disc system is not uniquely determined by the physics of accretion: however we know that the MAD region ($r\leq r_{m}$) will also be threaded by the external magnetic field which threads the disc, and it will have a density different from that of the disc itself. We assume that at the magnetopause discontinuity, the Alfvén velocity $v_{A}=\left(B^{2}/4\pi\rho\right)^{1/2}$ will be continuous, which is equivalent to assuming that the magnetosphere will be in an equilibrium thermal state, and will not be heated by Alfvén waves diffracted at $r_{m}$. Thus we have: $B\propto\rho^{1/2}$, and taking into account the disc density profile given before ($\rho\propto r^{-1}$), we eventually obtain:

[TABLE]

Thus, the magnetic flux within the magnetosphere will be given by:

[TABLE]

($r_{m}\gg r_{g}$) and, substituting $B_{m}$ from eq 12 above, we eventually get:

[TABLE]

The accretion rate onto the BH is defined as: $\dot{M}=2\pi r_{m}\epsilon v_{ff}=2\pi\epsilon\left(GM_{BH}r_{m}\right)^{1/2}$, where we have used the definition of free-fall velocity ($v_{ff}=\left(GM_{BH}/r_{m}\right)^{1/2}$).

Here $\epsilon\backsimeq 0.1\leq 1$ is a fudge factor which quantifies the average efficiency of accretion. Defining as usual the accretion rate in units of the Eddington accretion rate: $\dot{m}\equiv\dot{M}/\left(\eta M_{BH}\right)$, we finally arrive at an expression for $r_{m}$ as a function of the accretion rate:

[TABLE]

Inserting this expression into eq. 14 we get:

[TABLE]

where $\kappa$ is a purely numerical constant. To proceed further, we note that the backflow will modify the accretion rate at $r_{m}$ as given by eq. 7 above, thus we will have:

[TABLE]

and, using: $\rho(r_{m})=\Sigma_{m}/\left(r_{m}/2\right)$ and the expression for $v_{r}$ we eventually find that the reduced accretion rate during the backflow scales as: $\dot{m}\propto M_{BH}^{-1}\Sigma_{m}A(t)$. After having substituted in eq. 16, we eventually obtain:

[TABLE]

We plot in Figure 9 the temporal evolution of the jet’s mechanical power, in units of the (constant) power without backflow. The similarity with the evolution of the azimuthal flux in Figure 6 is a consequence of the hypothesis that the mass flux variation (the source term in eq. 6) generates variations of $\Sigma$ without any delay and/or modifications: the latter could arise for instance by finite conductivity effects which can threaten the ideal coupling between streamlines and magnetic fields, or by the enhancement of magnetic reconnection (see e.g. Punsly, 2015) when the backflow plasma carries a magnetic field. Thus, the results we have presented in this section represent only the simplest possible scenario, and we will consider more realistic ones in future papers (Antonuccio-Delogu et al. 2016, in prep.).

7 Conclusion

Backflows are a large-scale generic feature of jet propagation within the ISM of their host galaxies. The high-resolution numerical 3D experiments we have presented here confirm our previous findings concerning their physical origin (Antonuccio-Delogu & Silk, 2010), and in particular concerning their thermodynamic origin: it is the discontinuity in enthalpy near the hotspot, together with the finite curvature radius of this discontinuity, that generates a finite vorticity in the initially laminar jet flow, as predicted by Crocco’s theorem (Crocco, 1937).

There is also mounting observational evidence for backflows from analysis of FR-I radiogalaxies (Laing & Bridle, 2012), and from more recent studies concerning Cen-A: Neumayer et al. (2007) observed inflow in the kinematics of highly ionized gas in the nucleus of the source. They associate this feature with backflow of gas that was accelerated by the jet of Cen-A (hence the high ionization rate), confirming backflow presence down to parsec scales, as later confirmed by further studies (Hamer et al., 2015; Bicknell et al., 2013).

Until now, however, the most complete evidence for backflows has emerged from previous numerical simulations similar to those presented here (Laing & Bridle, 2014; Perucho & Martí, 2007), and also from simulations including a very inhomogeneous medium (Wagner et al., 2012; Wagner & Bicknell, 2011): in the latter, the backflow mostly originates from jet gas flows that have run into cold clouds and thus have an intermittent character. The high resolution, adaptive numerical 3D experiments we have presented here demonstrate instead that backflows are spatially coherent and can have a significant temporal extension, lasting for few tens of million years : both are structural features of jet propagation within their host galaxies.

The backflows will be fed as long as there are the physical conditions to generate them and support a finite curvature enthalpy discontinuity near the hotspot, either in the form of a shock or of a contact discontinuity. These features depend only slightly on $\rm{P}_{jet}$ , as in most cases we observe a highly intermittent ($\dot{M}\simeq 10^{-1}-1\,\rm{M}_{\sun}\rm{yr}^{-1}$) inflow towards the central BH region lasting about 15-20 Myrs. The morphology of the cavities may nonetheless change these numbers by a factor of order two (as seen from the differences between our EC and RC series, e.g. comparing Figures 4 and  5). Due to its axial symmetry, the backflow has a low azimuthal angular momentum w.r.t. the central accretion region. The combination of this circumstance and of the thermodynamics of its propagation conspire to drive a fraction of the backflow towards the very central regions.

We have explored the impact of backflows on the central unresolved accretion region near the supermassive BH, and in particular on the jet power, assuming backflow can be modelled as a hot gas accretion onto a MAD , this latter being a paradigm supported by many GRMHD simulations (Punsly, 2011; Ghisellini et al., 2010). The analysis we presented in Sect. 5 is not free of some assumptions. For instance, we have disregarded the dissipative processes in the innermost regions of the accretion disc, which will heavily affect the magnetic field strength and topology and dominate over the action of the time-dependent accretion flow. In our model the backflow will accrete into a thin, hot, high-$\beta$ disc, and we have shown that the jet power $\rm{P}_{jet}$ can be boosted by a factor of 10 or 20 for as long as 5 or 10 Myr. The temporal evolution of the disc density and jet power in this model is dependent on the profile and normalization of the initial disc: however, the final density profile will be independent of the details of the temporal evolution of backflow accretion. As is clear from Figure 6, all of the mass accretion episodes start about 3-5 Myr. after jet launching and decay after 20-30 Myr, independently of $\rm{P}_{jet}$ . However, as is clear from a visual inspection of Figure 8, despite large individual variations between different simulations, our ”template” model reproduces the generic features of accretion episodes.

From our analysis, it is evident that the backflow is the result of the interaction between the jet and the local host galaxy’s interstellar environment, and its contribution to the MAD demonstrates that a connection between galaxy-scale feedback and central accretion inevitably develops on time-scales of the order of $\sim 10^{6}$ years, or about $1/10$  of the AGN duty cycle. This *backflow accretion * time -scale however only refers to the typical time for the backflow to feed the MAD . The backflow phenomenon points to a deep connection between AGN feedback and SMBH accretion, as previously hinted at (Narayan & McClintock, 2008) 999Although magnetic phenomena may introduce further time delays, as briefly mentioned in Sec. 6.

Finally, we would like to point out that the backflow is a significant global dynamical feature of an AGN that is capable of ”bridging” the very large (kpc) scales, where jets propagate, with the accretion (subparsec) scales. The ”feedback” from the large scales is capable of modulating global properties of the jet, such as the mechanical power Pjet, which in turn can affect the thermodynamic properties near the hotspot, from where the backflow originates. This cycle creates a self-regulation mechanism which determines the duty cycle and other properties of the AGN, as we will show in a separate paper (Antonuccio-Delogu et al., in preparation).

Acknowledgements

We would like to acknowledge the anonymous referee for the careful reading and useful comments which improved the quality of the paper.

The work of S.C. has been supported by the ERC Project No. 267117 (DARK) hosted by Université Pierre et Marie Curie (UPMC) - Paris 6 and the ERC Project No. 614199 (BLACK) at Centre National De La Recherche Scientifique (CNRS). SC thanks Marta Volonteri for the profitable discussion and the precious advice.

The work of V.A.-D. …has partially been supported by the joint CNRS-INAF Project PICS 2013-2016 ”Modelling and Simulation of mechanical AGN Feedback”. V.A.-D. gratefully acknowledges the hospitality of IAP, Paris, during the completion of this work, particularly very useful conversations with G. Mamon and M. Volonteri. The work of JS was supported by ERC Project No. 267117 (DARK) hosted by Université Pierre et Marie Curie (UPMC) - Paris 6, PI J. Silk. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Centre at the University of Chicago.

Appendix A Numerical algorithms and code

The AMR code we have chosen to perform the numerical experiments described in this work (FLASH v. 4.2) solves the Euler equations system:

[TABLE]

[TABLE]

[TABLE]

where: $E=\epsilon+v^{2}/2$ and $\epsilon$ is the specific internal energy, $P$ is the (thermal) pressure and $\Lambda$ the cooling function.

Among the different numerical algorithms made available in FLASH we have adopted the shock-capturing PPM scheme (Woodward & Colella, 1984) which is particularly suited to model the shocks formed during the interaction of jets with the cocoon.

We use FLASH’s default ideal gamma equation of state for an ideal gas, $P=(\gamma-1)\rho\epsilon$, where $\gamma$ is the specific heat ratio ($5/3$ in our case). As mentioned in the text, we include a radiative cooling function extended to very high $T\geq 5\times 10^{9}$ K extended to very high temperatures to account for radiative losses due to $e^{+}-e^{-}$ annihilations (see the Appendix of Antonuccio-Delogu & Silk, 2008).

We adopt FLASH’s Multigrid Poisson Gravity Solver, with a custom modification for adding a static dark matter gravitational potential for the host NFW halo.

FLASH adopts the PARAMESH block-structured AMR decomposition of the spatial computational domain (MacNeice et al., 2000). In PARAMESH the initial grid is recursively refined dichotomically up to a maximum refinement level $l_{max}$. Moreover, each block is further divided into $n_{c}$ computational cells along each spatial direction. Thus, the minimum resolved block of cells has a linear size: $\Delta r=L/2^{l_{max}}n_{c}$. In all our runs we have: $\Delta r=\frac{640\ \mathrm{kpc}}{2^{10}\times 8}=78.125$ pc.

Finally, we allow refinements on the basis of a density and temperature gradient criterion, adopting FLASH default refinement strategy based on Loehner’s error estimator (Löhner, 1987, see FLASH user manual) set to $0.8$ for refinement and $0.6$ for de-refinement.