© The Authors 2026

1. Introduction

It has been known for more than two decades that the masses of supermassive black holes (SMBHs) correlate with the properties of their host galaxies (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000). Of the numerous relations proposed so far, the M − σ relation, which links SMBH mass to the velocity dispersion of the host-galaxy spheroid (Ferrarese & Merritt 2000; Gültekin et al. 2009; McConnell & Ma 2013; Saglia et al. 2016; de Nicola et al. 2019), is the most fundamental (Shankar et al. 2013; Marsden et al. 2020). Its existence is readily explained as the result of active galactic nucleus (AGN) feedback. The energy released via luminous accretion onto the SMBH during the AGN phase couples to the interstellar medium (ISM) and regulates the growth of both the SMBH and the host galaxy, establishing various correlations between the SMBH mass and galaxy properties. Feedback manifests in the form of radiation, narrow jets, and/or wide-angle outflows (Fabian 2012; King & Pounds 2015; Harrison & Ramos Almeida 2024). Of these, the wide-angle wind-driven outflow model (see King & Pounds 2015; Zubovas & King 2019 for recent reviews) is the most successful in explaining galaxy-scale effects, such as the properties of observed outflows, their influence on galaxy properties, and the M − σ relation.

Within the framework of AGN wind-driven feedback, AGN luminosity (L_AGN) is communicated to the surrounding gas via a quasi-relativistic wind emanating from the accretion disc (King 2010a). The wind moves with velocities (v_w) of ∼0.1c and carries a kinetic power (Ė_w) of ∼0.05L_AGN. Upon encountering the relatively static ISM, the wind develops a strong shock with a post-shock temperature (T_sh) of ∼10¹⁰ K. The pressure in the shocked wind bubble is significantly higher than that of the ISM, so the bubble begins to expand. At the same time, the inverse-Compton process–the only efficient cooling process at this temperature–cools the shocked gas. Depending on which of the two processes–expansion or cooling–operates on a shorter timescale, the resulting outflow can be either momentum- or energy-driven. An energy-driven outflow is approximately adiabatic and carries most of the wind energy (Ė_kin,out), i.e. ∼0.02L_AGN (the rest of the wind energy is used up to do work against gravity and pdV work). This results in a large-scale (R > 1 kpc), fast (v_out ∼ 1000 km s⁻¹), and massive (Ṁ_out ∼ 1000 M_⊙ yr⁻¹) outflow. These properties agree well with those of observed outflows in AGN host galaxies (Zubovas & King 2012; Cicone et al. 2014; Fluetsch et al. 2019; Lutz et al. 2020). The momentum-driven outflow carries only the direct momentum of the wind, ṗ_out ∼ L_AGN/c, and is almost two orders of magnitude less powerful. The condition for this type of outflow to overcome the gravitational potential of the host galaxy leads to the M − σ relation (King 2010a).

How do we reconcile the fact that the M − σ relation requires transferring only the wind momentum to the gas, while large-scale outflows appear to be driven by adiabatic expansion of the same shocked wind bubble? At least two solutions exist. The first suggests that the wind cools efficiently close to the SMBH, where the AGN radiation field is stronger (King 2003, 2010a), and transitions to an adiabatic state beyond the cooling radius (R_cool) of ∼0.5 kpc. Whether this happens is a matter of some debate: although at least one galaxy, NGC 4051, shows a spectral feature that resembles the expected signature of the cooling wind (Pounds & Vaughan 2011), in general it is not found in AGNs (Bourne & Nayakshin 2013).

The other possibility is that the outflow type depends on the density of impacted gas. If the ions and electrons in the shocked wind form a two-temperature plasma, the cooling rate diminishes significantly and the cooling radius reduces to < 1 pc (Faucher-Giguère & Quataert 2012). This is consistent with most observed AGN spectra (Bourne & Nayakshin 2013). However, the outflow expands primarily through low-density channels in the clumpy ISM, leaving cool dense gas behind (Zubovas & Nayakshin 2014; Ward et al. 2024). The dense gas is then affected primarily by the wind momentum, and the condition that this gas must be pushed away establishes the M − σ relation.

This description is qualitative because quantitative analytical calculations are impossible when one considers a multi-phase, turbulent, and non-spherical ISM. In a recent paper (Zubovas et al. 2024, hereafter Paper I), we investigated the effects of turbulence and cooling on energy-driven outflows and found that the inclusion of cooling leads to the formation of multi-phase outflows, where the dense gas may be affected mostly by the wind momentum rather than its energy. Here, we extend the analysis of dense gas motion in outflows driven by AGNs with different luminosities. We show that AGNs with luminosities well below the Eddington limit are unable to prevent the significant accretion of cold dense gas, while those with L_AGN > L_Edd efficiently remove cold gas and stifle further SMBH growth, as predicted by the momentum-driven outflow formalism. We also show that dense gas clusters mostly experience feedback consistent with pure momentum driving.

We briefly review the analytical derivation of outflow parameters in Sect. 2, with an emphasis on the M − σ relation and conditions for the removal of cold gas. We present the numerical setup in Sect. 3 and the simulation results in Sect. 4. We discuss the implications of our results in Sect. 5, focusing on establishing and maintaining the M − σ relation across cosmic time, the expected properties of multi-phase outflows, and the importance of gas self-gravity and magnetic fields. We summarise the main results and conclude in Sect. 6.

2. The M − σ relation in different kinds of outflows

In the momentum-conserving regime, the equation of motion of a spherical outflow shell at distance r is (following King 2010a; King & Pounds 2015)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{p}_{\mathrm{out}}}{\mathrm{d}t}& \equiv \frac{\mathrm{d}}{\mathrm{d}t}\left[M_{\mathrm{g}}\left(r\right)\dot{r}\right]\hfill \\ \hfill & =\frac{L_{\mathrm{AGN}}}{c}-\frac{{\mathit{GM}}_{\mathrm{g}}\left(r\right)[M_{\mathrm{BH}}+M_{\mathrm{b}}\left(r\right)+M_{\mathrm{g}}\left(r\right)/2]}{r^2},\hfill \end{array}$$(1)

where M_g(r) is the mass of the swept-up gas within the outflow radius r, M_BH is the SMBH mass and M_b and M_g represent the non-gaseous and gaseous components of the distributed mass of the host, all measured within the outflow radius. Assuming that the AGN luminosity is equal to the Eddington luminosity L_Edd = 4πGM_BHc/κ, where κ ≈ 0.4 cm² g⁻¹ is the electron scattering opacity, and that the gas and total matter are distributed isothermally, with M(r) = 2σ²r/G, M_g(r) = f_gM(r) and M_b(r) = (1−f_g)M(r), we can show that the outflow can only expand to an arbitrarily large radius if the SMBH mass exceeds

$$\begin{array}{c}\hfill M_{\sigma }=\frac{f_{\mathrm{c}}\kappa }{\pi G^2}{\sigma }^4\approx 3.8\times {10}^8{\sigma }_{200}^4{\mathrm{M}}_{\odot };\end{array}$$(2)

here, f_c = 0.16 is the cosmological baryon fraction and we scale the velocity dispersion to σ₂₀₀ ≡ σ/(200 km s⁻¹). At large radii, the outflow velocity can be expressed as

$$\begin{array}{c}\hfill {\dot{r}}^2\approx 2{\sigma }^2(\frac{\mathit{lm}}{f}-1),\end{array}$$(3)

where m ≡ M_BH/M_σ, f ≡ f_g/f_c is the current gas fraction scaled to the cosmological value, and l ≡ L_AGN/L_Edd. This equation clearly has no solution unless lm/f > 1, i.e. the outflow can only expand to large radii if the black hole mass is large enough, the AGN is super-Eddington or the gas density is low.

The kinetic power of the outflow is

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{kin},\mathrm{m}}=\frac{{\dot{M}}_{\mathrm{out}}{\dot{r}}^2}{2}\lesssim \frac{L_{\mathrm{AGN}}}{2}\frac{\dot{r}}{c},\end{array}$$(4)

where we used the upper limit Ṁ_out ṙ ≲ L_AGN/c obtained from Eq. (1) by neglecting gravity. Using Eq. (3), the fraction of AGN luminosity transferred to the outflow is then

$$\begin{array}{c}\hfill \frac{{\dot{E}}_{\mathrm{kin},\mathrm{m}}}{L_{\mathrm{AGN}}}\approx \frac{\sigma }{\sqrt{2}c}\sqrt{\frac{\mathit{lm}}{f}-1}\approx 5\times {10}^{-4}{\sigma }_{200}\sqrt{\frac{\mathit{lm}}{f}-1}.\end{array}$$(5)

In the energy-driven regime, the equation of motion is similar, but the driving term is replaced with 4πr²P, where P is the pressure of the shocked wind bubble. It is calculated using the energy equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\mathit{PV}}{\gamma -1}\right)=\frac{\eta }{2}{L}_{\mathrm{AGN}}-P\frac{\mathrm{d}V}{\mathrm{d}t}-\frac{G{M}_{\mathrm{g}}(M_{\mathrm{b}}+M_{\mathrm{g}}/2)}{r^2}\dot{r},\end{array}$$(6)

where γ is the specific heat ratio of the gas. We have dropped the (r) notation for brevity. The derivation of the full equation of motion is algebraically involved, but straightforward; we refer the reader to Appendix A of Zubovas et al. (2022) for details. In its general form, the equation of motion has no analytical solutions. If we make the same simplifying assumptions regarding the gas distribution as before, a constant-velocity solution emerges (King & Pounds 2015, cf.):

$$\begin{array}{c}\hfill v_{\mathrm{e}}\equiv \dot{r}\approx {\left(\frac{2\eta l{\sigma }^{2}c}{3}\frac{f_{\mathrm{c}}}{f_{\mathrm{g}}}\frac{M_{\mathrm{BH}}}{M_{\sigma }}\right)}^{1/3}\approx 925{(lm/f)}^{1/3}{\sigma }_{200}^{2/3}{\mathrm{km}\mathrm{s}}^{-1}.\end{array}$$(7)

The kinetic power of the outflow is ∼1/3 of the total energy injected into the ISM, i.e.

$$\begin{array}{c}\hfill \frac{{\dot{E}}_{\mathrm{kin},\mathrm{e}}}{L_{\mathrm{AGN}}}\approx 0.016,\end{array}$$(8)

which is ∼30/σ₂₀₀ times higher than in the momentum-driven case. It is clear that an outflow driven by all of the AGN wind energy can escape from the galaxy much more easily, and so the SMBH growth would be quenched at a much lower mass. This unrealistic outcome can be prevented in several ways. The AGN can always be significantly sub-Eddington, although this appears unlikely (King 2010b; Mountrichas & Georgantopoulos 2024). The gas can be very dense, with f ≫ 1, but this leads to rapid fragmentation into stars. Finally, the ISM can be ‘porous’ (Silk & Rees 1998), allowing the energy of the outflow to leak out through low-density channels (Zubovas & Nayakshin 2014; Ward et al. 2024). Cold gas clumps left behind the fast outflow then experience both a confining pressure of being embedded in a shocked wind bubble and a combined drag and directional pressure force from the wind. These forces are of order ρ_wv_w² per unit cross-sectional area of the cold gas clump (King & Pounds 2015), which is precisely the momentum rate of the AGN wind. So it is possible that dense gas is affected primarily, or even solely, by the AGN wind momentum, while the energy of the shocked wind is carried by the outflow of diffuse gas.

3. Numerical simulations

Our numerical setup is very similar to that of turbulent cooling simulations of Paper I. To recap, we used a modified version of GADGET-3 (Springel 2005) with the SPHS formulation of the main hydrodynamical equations (Read et al. 2010; Read & Hayfield 2012) and a Wendland kernel (Wendland 1995; Dehnen & Aly 2012) with particle neighbour number N_ngb = 100. The initial conditions consist of an SMBH and a turbulent gas shell embedded in a static background gravitational potential. The SMBH particle has a mass of M_BH = 10⁸ M_⊙ and is fixed at the centre of the simulation. Using Eq. (2), we derived the velocity dispersion σ_b = 142 km s⁻¹ and used it to generate the background potential. The potential corresponds to an enclosed mass M_b(< 1 kpc) = 9.4 × 10⁹ M_⊙ but extends indefinitely far. The gas is initially distributed in a shell between r_in = 0.1 kpc and r_out = 1 kpc. The total mass is M_g(< 1 kpc) = 9.4 × 10⁸ M_⊙ = 0.1 M_b(< 1 kpc), i.e. our modelled system has a gas fraction f_g = 0.1. A real gas-rich bulge in a galaxy with a 10⁸ M_⊙ SMBH would extend farther and be more massive, but we are only interested in gas dynamics within the central several hundred parsecs, so truncating the gas distribution closer in does not affect our results. The gas is tracked with N ≈ 1.08 × 10⁶ smoothed-particle hydrodynamics (SPH) particles (the exact number differs in simulations with different initial conditions, as explained below), giving a particle mass m_SPH ≈ 875 M_⊙. Globally, the azimuthally averaged gas density falls approximately as ρ ∝ r⁻², but the full density distribution is produced by giving the gas a turbulent velocity field (following the method outlined in Dubinski et al. 1995 and Hobbs et al. 2011) and allowing the gas shell to relax for 1 Myr. The initial characteristic turbulent velocity is σ_t = 149 km s⁻¹.

As the simulation proceeds, the gas is affected by external gravity (we neglected self-gravity in the simulations; the importance of this assumption is discussed in Sect. 5.2) and AGN wind feedback. The wind was not modelled hydrodynamically, but instead propagated on a static spherical grid, as described in Sect. 3.3 of Paper I. This prescription ensures a spherically symmetric injection of AGN wind energy and momentum and conserves these quantities, as we show in Appendix A. We employed two sub-resolution prescriptions for gas cooling. At temperatures T_g > 10⁴ K, we used the prescription of Sazonov et al. (2005), which models the heating and cooling of optically thin gas due to a typical AGN radiation field. It includes photoionisation heating, the Compton effect and metal line cooling. We modified this function by neglecting the effect of Compton cooling, a change appropriate for the expected two-temperature plasma nature of the hottest outflowing gas (see Faucher-Giguère & Quataert 2012; Bourne & Nayakshin 2013). Overall, the heating effect was negligible, especially for dense gas, as expected (cf. Sazonov et al. 2005, who showed that direct thermal feedback is only relevant for gas with densities 1–2 orders of magnitude lower than the average density in our simulations). Below T_g = 10⁴ K, the gas cooled according to the prescription of Mashchenko et al. (2008), which approximates the cooling effect of metastable C, N, O, Fe, S, and Si lines in ionisation equilibrium maintained by locally produced cosmic rays at approximately Solar metallicity. This function allowed the gas to cool to 20 K. All gas particles that fell closer than r_acc = 0.01 kpc from the black hole particle were removed from the simulation. We used this to track the expected accretion onto the SMBH, but did not adjust the AGN luminosity because the time for gas to fall from ∼10 pc through the accretion disc to the SMBH is comparable to the duration of our simulations. A summary of the main simulation parameters and salient results is given in Table 1.

4. Results

Qualitatively, the simulations can be divided into three groups. When L_AGN/L_Edd ≤ 0.5, the AGN is too weak to significantly affect gas infall. When 0.7 ≤ L_AGN/L_Edd ≤ 1.7, the effect is moderate and allows some infall, although most of the gas is pushed away. AGNs with L_AGN/L_Edd ≥ 1.9 essentially completely shut down gas infall. We therefore mostly focus on presenting the results of three simulations–L0.5, L1.0 and L2.0–in detail, including the control simulation where appropriate, and limit the depiction of other simulations’ results to those plots where we show direct parameter dependence on AGN luminosity.

We start by showing density maps to describe the qualitative evolution of the gas distribution. We then analyse the properties of infalling gas, using kinematic phase maps, radial mass profiles and time evolution of infalling gas mass to highlight the differences among the three simulations. Then we consider the feedback effect on individual dense gas clusters. Finally, we show the differences in SMBH particle growth and determine the AGN luminosity that would produce a self-consistent growth rate.

4.1. Gas morphology

Figure 1 shows the evolution of gas density in runs L0.5, L1.0, and L2.0 over 0.25–0.75 Myr. In L0.5, the AGN disturbs the gas and inflates bubbles that gradually expand to a radius ∼0.3 kpc in all directions by 0.75 Myr, although gas disturbances are seen in some directions as far as ∼0.5 kpc. This corresponds to an average velocity ∼400–670 km s⁻¹, somewhat lower than the analytical expectation v_e ≈ 680 km s⁻¹ (Eq. (7)). However, the bubbles are faint and only affect the diffuse gas. The region around the SMBH is not cleared out; dense gas continuously falls inwards and feeds the SMBH particle, as we show below.

Fig. 1. Gas density integrated through a 20° wedge around the mid-plane (z = 0), defined by |z| < 0.18(x2+y2), in simulations L0.3 (top row), L1.0 (middle), and L1.7 (bottom). Columns show t = 0.5, 0.75, and 1.0 Myr from left to right. Brighter colours indicate higher gas densities, in the range ∼10−25 g cm−3 < ρ < 10−19 g cm−3.

In simulation L1.0 (middle row), the outflow is significantly faster (v_r ∼ 670–1000 km s⁻¹) and almost entirely clears the central ∼100 pc by 1 Myr. SMBH accretion is essentially shut down during the activity phase. However, dense gas ‘fingers’ remain and even approach the nucleus, so once the AGN switches off, they may provide significant fuel to restart the activity and continue growing the SMBH.

Finally, in simulation L2.0, the outflow breaks out of the initial gas shell by t = 0.75 Myr. The effect on dense gas is particularly noticeable: the outflow clears a ∼200 pc-wide region by 0.75 Myr and pushes out both dense and more diffuse gas. Additionally, dense gas is confined to isolated elongated clouds rather than ‘fingers’ as in L1.0.

4.2. Mass infall

To get a better sense of the properties of infalling gas, we show in Fig. 2 kinematic phase diagrams across the three luminosity runs together with our adopted criterion for selecting ‘cold’ (T < 3 × 10⁴ K) and ‘rapidly infalling’ (v_rad ≤ −σ_b) gas. The latter criterion identifies gas with inward motion that is not the result of turbulence. All distributions are qualitatively similar, with a very high phase density region corresponding to the outer parts of the initial distribution that have not been affected by the outflow yet. Gas radial velocities range from ∼−1000 km s⁻¹ to > 4000 km s⁻¹, although most of the gas has velocities much closer to zero. Some gas has more negative radial velocities than the lowest values in the control simulation. Its low density and high temperature indicate hot gas expanding into the region evacuated by the outflow. This effect may be a numerical artefact caused by the lack of hydrodynamical particles representing the AGN wind (see Zubovas et al. 2024). The gas density distribution does not change very significantly compared with the control model, except for the maximum gas density decreasing by a factor of ∼1.5–2 in all simulations. This suggests that AGN feedback prevents the formation of the densest clumps. Nevertheless, a lot of dense gas with ρ > 10⁻²¹ g cm⁻³(n ≳ 10³ cm⁻³) exists in all simulations. The temperature distribution is also qualitatively similar in all simulations with AGNs, but differ markedly from the control model, because even a low-luminosity AGN can rapidly heat and accelerate diffuse gas. The multiple peaks seen in the temperature histogram are imposed by the adopted cooling function, which has minima at T ∼ 10⁴ K and T ∼ 5 × 10⁴ K. The highest temperatures result from shock heating of the outflowing gas. We used the approximate midpoint between the two peaks, T = 3 × 10⁴ K, to divide the gas into ‘cold’ and ‘warm-hot’ phases. Subsequently, we focused on the infall of cold gas, which can be gravitationally bound to the SMBH.

Fig. 2. Kinematic phase diagrams at t = 0.5 Myr in simulations L0.5 (left), L1.0 (middle), and L2.0 (right). In each panel, the top 2D histogram shows radial velocity against density, and the bottom the radial velocity against temperature. 1D histograms at the top and right show distributions of individual gas properties. Grey contours show corresponding distributions in the control simulation. Vertical dotted lines show vr = 0 and vr = −σb; horizontal dotted line shows T = 3 × 104 K.

For the purposes of our investigation, the most significant result is the change of dense gas (ρ > 10⁻²¹ g cm⁻³) velocities with L_AGN. This gas is exclusively cold and lacks thermal energy to escape the SMBH gravitational potential. As a result, it can feed the SMBH efficiently when inflowing. In the control simulation, most of the dense gas is infalling by t = 0.5 Myr, as expected, as the gravitational potential overcomes decaying turbulent velocities. In L0.5, the balance shifts towards higher velocities with essentially no gas having v_rad < −280 km s⁻¹. However, most dense gas retains negative radial velocities. Higher AGN luminosities gradually change this distribution: in L1.0, the lowest velocities of infalling dense gas at t = 0.5 Myr increase to v_rad ∼ −180 km s⁻¹, while in L2.0, they further increase to v_rad ∼ −145 km s⁻¹. Additionally, in L1.0, the total mass of infalling dense gas becomes comparable to the outflowing mass, while in L2.0, the outflow clearly dominates.

The radial gas distribution (Fig. 3) also reveals the evolution of outflowing and inflowing material. In L0.5, the outflow is initially able to clear a ∼30 pc region (also seen in Fig. 1, top left). However, the infalling gas quickly returns and fills the central region, reaching (0.7–1) × 10⁶ M_⊙ per radial bin, approximately half the value in the control simulation. In L1.0, the situation is markedly different. Initially, the outflow clears the central ∼40 pc. Beyond that, the gas mass decreases, but persistent dense clumps remain at r ≲ 250 pc after the outflow has passed. Some of them gradually approach the SMBH, reaching the sink radius around t = 0.5 Myr. Outside this region, cold gas is gradually evacuated and does not return until the end of the simulation. The total infalling mass remaining in the central region is only a few per cent of the control simulation value. In L2.0, evacuation is more efficient: the bulk of the material is cleared within t ≈ 0.5 Myr, and the outflow maintains a clear zone of r ∼ 50 pc until the end of the simulation. The total mass of cold rapidly infalling gas in each bin never exceeds 1% of that in the control simulation.

Fig. 3. Mass of cold rapidly infalling gas in radial 3.3-pc-wide bins at 0.25, 0.5, and 0.75 Myr (solid red, blue, and green lines and shading, respectively) in simulations L0.5, L1.0, and L2.0 (top, middle, and bottom panels, respectively). Solid lines show the mean of the four stochastically different simulations; shading encompasses their full range. Dotted lines in each panel show equivalent results from the control simulation.

The top panel of Fig. 4 shows the time evolution of rapidly infalling cold gas mass within the central 500 pc in the four simulations, including control. We interpret these values as the maximum gas mass that could feed the SMBH once the AGN episode has finished (neglecting possible star formation, discussed in Sect. 5.2). The total includes the gas mass accreted by the SMBH particle, though its contribution remains below 10% except in the control simulation. In the control simulation, the infalling gas mass increases continuously as gas from outer regions approaches the SMBH. In L0.5, the outflow delays this process, but cannot prevent cold dense clumps from falling inwards, maintaining a large total infalling mass. L1.0 shows a gradual decrease from an initial value just above 10⁸ M_⊙ to a minimum of ∼7.8 × 10⁶ M_⊙ around t = 1 Myr. Later, when dense clumps overtaken by the outflow fall inwards, the infalling gas mass slowly increases to ∼2 × 10⁷ M_⊙. Note that the minimum mass is reached at significantly different times in stochastically different simulations: two reach minima already at t ∼ 0.75 Myr; the minimum mass varies by a factor of ∼3. In L2.0, the decrease is faster, down to ∼3.1 × 10⁵ M_⊙ by t = 0.6 Myr, but later the infalling mass increases back up to ∼4 × 10⁶ M_⊙ for the same reasons as L1.0. The time of minimum value shows a similar scatter among stochastically different simulations as L1.0, but the minimum mass ranges over more than two orders of magnitude.

Fig. 4. Top: Total mass of cold rapidly infalling gas within the central 500 pc as a function of time in the control (red), L0.5 (blue), L1.0 (green), and L2.0 (magenta) simulations. We added the mass accreted by the SMBH particle to this total. Solid lines show the mean of the four stochastically different simulations; shading encompasses their full range. Squares mark the times when the mean mass reaches its minimum value. Bottom: Minimum value of cold rapidly infalling gas mass in all simulations (circles) as a function of luminosity, with colours denoting the time when the minimum is reached. The solid red line is the best linear fit between log M and LAGN of all simulations except the control and L0.1; the dashed lines indicate 95% confidence limits on the line parameters, estimated via bootstrapping.

The bottom panel shows the minimum rapidly infalling cold gas mass within the central 500 pc in each simulation versus AGN luminosity. In the control and L0.1 simulations, the infalling mass only increases, so the minimum occurs at t = 0 and is equal to the initial cold rapidly infalling gas mass M_{cold, in, init}. With increasing luminosity, the minimum mass follows an exponential trend: M_{cold, in, min} ≈ 4.3 × 10⁸exp(−1.67L/L_Edd) M_⊙. Note that M_{cold, in, min} < 0.1 × M_{cold, in, init} in all simulations with L ≥ L_Edd and lower than M_{cold, in, min} < 0.01 × M_{cold, in, init} when L ≥ 1.9L_Edd. We conclude that AGNs with luminosities L ≳ L_Edd(M_σ) can significantly disrupt and remove dense gas reservoirs that would otherwise feed the SMBH.

4.3. Feedback on dense gas clumps

Each gas particle experiences feedback differently not only because of its properties, such as density or temperature, but also due to environment. For example, diffuse gas can be shielded by dense clumps, while particles can enter and leave such clumps. To quantify the effect of AGN feedback on cold dense gas, we identified clusters of particles at t = 0.75 Myr and tracked their evolution between two snapshots, i.e. over Δt = 0.05 Myr. Using DBSCAN from sklearn (Ester et al. 1996; Pedregosa et al. 2011), we searched for clusters defined as groups of at least 100 particles with density ρ > 10⁻²² g cm⁻³ and T < 3 × 10⁴ K, whose particles have nearest neighbours closer than δd < 10 pc and relative velocities δv < 0.5σ. The top panel of Fig. 5 shows detected clusters in simulation L1.0 superimposed on the grey scale density map. Not all dense regions are assigned to clusters, usually because they contain fewer than 100 particles.

Fig. 5. Top: Clusters (coloured by gas density) detected at t = 0.75 Myr in simulation L1.0, plotted over the gas density map (grey scale). Middle panel: Momentum loading factor (fp) distributions (Eq. (9); box and whisker plots) in simulations with different AGN luminosities. Circles highlight outliers beyond 1.5 times the interquartile range covered by the whiskers. The horizontal dashed line highlights fp = 1, which corresponds to pure momentum driving, for ease of comparison. Bottom panels: Momentum loading factor (fp) as a function of dcl (left) and mean gas density (right) in simulations L0.5 (red), L1.0 (green), and L2.0 (blue).

We quantified the AGN wind effect on each cluster using Eq. (1) to calculate the momentum loading factor:

$$\begin{array}{c}\hfill f_{\mathrm{p}}=\frac{\mathrm{\Delta }\left(M_{\mathrm{cl}}{v}_{\mathrm{rad},\mathrm{cl}}\right)/\mathrm{\Delta }t+G{M}_{\mathrm{cl}}{M}_{\mathrm{pot}}(<d_{\mathrm{cl}})/d_{\mathrm{cl}}^2}{L_{\mathrm{AGN}}/c\times {\mathrm{\Omega }}_{\mathrm{cl}}/4\pi },\end{array}$$(9)

where M_cl is the cluster mass, v_rad, cl is its mean radial velocity, d_cl is the distance between the origin and its centre of mass and Ω_cl is the solid angle subtended by the cluster when viewed from the origin. The numerator contains the change in cluster radial momentum over time and the gravitational force from the background potential; the denominator is the fraction of the AGN wind momentum L_AGN/c directly impinging on the cluster.

In the middle panel of Fig. 5, we show the distributions of f_p in simulations with different AGN luminosities at t = 0.75 Myr using box-and-whisker plots. A spherically symmetric energy-driven outflow has f_p ∼ v_w/v_out ≳ 30, higher than almost all clusters. Despite the large scatter, most values fall in the range 0.1 < f_p < 1.4, with a mean of f_p ≈ 0.33 and a slight positive correlation with luminosity. So most clusters experience an effective force lower than or comparable to that of the AGN wind momentum, even though the shocked wind does not cool. This results from self-shielding: gas on the inner edge of the cluster is exposed to the wind and absorbs its momentum, while gas farther away is shielded even from this modest effect. Values of f_p > 1 occur when clusters are pushed by the shocked wind plasma before it escapes through low-density channels.

In the bottom two panels of Fig. 5, we plot f_p versus d_cl and mean cluster gas density. Each simulation shows a slight negative correlation between f_p and d_cl. This results from the shocked AGN wind confinement: the hot plasma attempts to escape through low-density channels, but cannot do so indefinitely fast and pushes against dense gas too. Near the SMBH, clusters subtend relatively large solid angles, so the shocked wind plasma takes longer to escape and pushes them more. There is a positive correlation between f_p and cluster density. In denser clusters, self-shielding is less important because the material on the inner edge of the cluster efficiently pushes against that further out. We find no significant correlations with other cluster parameters, such as their sizes, temperatures, masses or radial velocities.

These results show that cold dense gas clusters experience almost exclusively momentum feedback from the AGN wind, even while the shocked wind plasma remains hot and drives an energy-conserving outflow of diffuse gas.

4.4. SMBH growth

Despite AGN outflow feedback, some gas falls through the sink radius and is accreted by the SMBH particle. We show this Fig. 6. In the top panel, we plot the evolution of total mass accreted by the sink particle. In all simulations except control, the SMBH rapidly consumes a few particles (< 10⁴ M_⊙) heated by the AGN, followed by a quiescent period of zero growth, when the outflow has cleared the inner region and dense gas clumps have not returned. In L0.5, this period lasts around 0.4 Myr; later, the SMBH particle begins growing at an accelerating rate, consuming ∼5% of its initial mass within < 1 Myr, equivalent to an average accretion rate ∼2.5 times higher than Eddington. However, the SMBH itself would not grow at this rate, because the SMBH particle growth represents infall of gravitationally bound matter through the sink radius at r_sink = 10 pc. In L1.0, the quiescent period lasts 0.75 Myr, followed by two short yet significant bursts of accretion, swallowing ∼6 × 10⁴ M_⊙ in ∼0.1 Myr and another ∼10⁶ M_⊙ in 0.2 Myr about 0.25 Myr later. These bursts correspond to accretion rates ∼0.3–2.5 times Eddington. Nevertheless, the SMBH particle growth is significantly limited by the available mass supply, which in this simulation is only ∼7.8 × 10⁶ M_⊙ ∼ 0.08M_BH. The SMBH cannot grow by more than a few per cent of its initial mass because the AGN episode has expelled most of the surrounding gas. In L2.0, the quiescent period lasts ∼1.2 Myr, followed by only one short accretion burst of a few times 10⁵ M_⊙. Based on the evolution of radial gas profiles (see the previous subsection), we expect at most a few more such bursts powered by the remaining available gas.

Fig. 6. Top: Same as the top panel of Fig. 4 but showing the change in SMBH particle mass. The horizontal dashed line shows M = 10mSPH. Bottom: Change in SMBH particle mass between the start of the simulation and t = 0.25, 0.5, 0.75, and 1 Myr (red, blue, green, and magenta circles, respectively), in simulations with different AGN luminosities. Error bars represent the range of values in stochastically different simulations; they are slightly offset horizontally for clarity. Solid lines show the SMBH mass change expected due to accretion powering an AGN of the given Eddington fraction over the same time intervals, assuming a radiative efficiency (η) of 0.1; the shaded region represents the effect of varying radiative efficiency between 0.057 and 0.42. The horizontal dashed line is the same as in top panel.

The behaviour described above suggests that an AGN shining at L = 0.5L_Edd is unable to prevent Eddington-level accretion onto the SMBH, while a luminosity exceeding L = L_Edd is sufficient for this. To quantify this trend, we show, in the bottom panel of Fig. 6, the SMBH particle mass change against AGN luminosity and the mass change that would result from accretion producing that luminosity. It is obvious that in simulations with L ≤ 0.5L_Edd, the SMBH particle accretes more gas that would be needed to maintain the luminosity, i.e. feedback is unable to regulate the accretion rate. Conversely, when L_AGN ≳ 0.7L_Edd(M_σ), the AGN severely limits the accretion rate on to the SMBH, i.e. feedback is effective at preventing further SMBH growth. After an AGN episode with L_AGN ≳ 0.7L_Edd(M_σ), the SMBH would stay quiescent for a prolonged period, potentially indefinitely, unless external factors (e.g. instabilities, galaxy mergers, tidal disruption events) reignite the AGN later. So the critical AGN luminosity required for feedback to become efficient in regulating SMBH growth is very close to the Eddington luminosity for a black hole with mass given by M_σ (i.e. Eq. (2)). This fact strongly supports the hypothesis that momentum feedback is responsible for establishing the M − σ relation, even when the shocked AGN wind transfers all of its kinetic power to the surrounding gas.

5. Discussion

5.1. Establishing and maintaining the M − σ relation

We have shown that an energy-driven outflow expanding through a turbulent gas distribution affects cold dense gas primarily through momentum injection. This, in principle, recovers the critical luminosity condition, based on the momentum-driven outflow paradigm, for shutting down further SMBH growth (Murray et al. 2005), which leads to the establishment of the M − σ relation (King 2010a). Below, we outline several additional steps between our results and the establishment of the relation that our simulations are not designed to track.

Once the dense clumps and filaments are removed, further accretion of material on to the SMBH must be prevented or, at least, its average rate reduced significantly. In simulations with high AGN luminosities, L > 0.7L_Edd, the amount of dense infalling gas is reduced to < 0.1M_BH. However, once the AGN switches off, some of the outflowing gas should inevitably decelerate and fall back on to the SMBH, even if the outflow itself persists for many times longer than the AGN episode (King et al. 2011; Zubovas & Maskeliūnas 2023). This can restart the AGN after as little as t_fall ∼ R/σ ∼ 3.5 Myr, assuming gas infall from outside the 500 pc radius. The duty cycle implied by this timescale, compared with the duration of our simulations (t_AGN ∼ 1–1.5 Myr), is very large, δ_AGN ∼ t_AGN/(t_AGN+t_fall) ∼ 0.2.

In between episodes, some cold clumps disperse and their gas is removed more efficiently during later AGN phases. This leads to a decrease in the duty cycle, until eventually the bulge is essentially empty of cold gas, terminating SMBH growth. We can expect the SMBH mass to grow by a factor of a few over these multiple episodes. Our simulation results imply self-regulation of SMBH growth at around L = 0.7L_Edd (Fig. 6, bottom panel), which is equivalent to an Eddington luminosity for a black hole with M_BH = 0.7M_σ. Increasing this by a factor of a few leads to a final mass close to, or slightly above, the value given by the M − σ relation. The fact that the M − σ relation is established over multiple AGN episodes is consistent with the conclusion that SMBH mass is a better predictor of star formation quenching than any instantaneous AGN property, as indicated by both simulations and observations (Harrison 2017; Martín-Navarro et al. 2018, 2021; Piotrowska et al. 2022; Bessiere et al. 2024).

The M − σ relation arose early in cosmic history; it was already present at z > 6 (Maiolino et al. 2024; Juodžbalis et al. 2025). Since then, subsequent galaxy and SMBH growth maintained it. The kinematic phase plots (Fig. 2) show an anti-correlation between density and velocity in the dense gas. This is qualitatively expected because denser gas has a higher weight and needs more force to be pushed away. It follows that in gas-rich galaxies, the luminosity required to shut down accretion is higher than in gas-poor ones. Given that galaxies at earlier cosmic epochs had higher gas masses (Carilli & Walter 2013), we expect high-redshift galaxies to reach higher SMBH masses before long-term quenching can occur. At later times, when gas is consumed and ejected, the luminosity required to eject new infalling gas clouds decreases. As a result, the SMBH no longer grows significantly, even over gigayear timescales, except during major mergers, which also increase the velocity dispersion. Furthermore, since self-regulation of SMBH feeding occurs at lower luminosities as the Universe evolves, the typical Eddington ratios should decrease with time, as observed (Trump et al. 2009; Willott et al. 2010; Mazzucchelli et al. 2023).

At the highest redshifts (z ≳ 6), dense gas streams falling from intergalactic space (Khandai et al. 2012; Heintz et al. 2024) can be completely resilient to AGN feedback, as any outflows generated expand around them, while the AGN wind momentum is insufficient to stop gas infall and push it away (also see Nayakshin & Power 2010). This can lead to very efficient SMBH growth, with duty cycles essentially equal to unity and frequent super-Eddington episodes, again as observed (Wu et al. 2022; Lupi et al. 2024; Suh et al. 2025).

5.2. Self-gravity and star formation

Our main simulations neglected gas self-gravity and star formation to conserve computational resources. To test the importance of these effects, we set up four additional simulations–L0-sg, L0.5-sg, L1.0-sg, and L2.0-sg. They included gas self-gravity and a crude star formation prescription: gas particles with Jeans mass lower than the resolved mass m_res = 100m_SPH ≈ 10⁵ M_⊙ were converted into star particles (see e.g. Tartėnas & Zubovas 2022, for details).

Figure 7 shows the gas morphology in simulation L2.0-sg at t = 0.75 Myr. This is directly comparable to the bottom right plot of Fig. 1. Overall, the outflow is remarkably similar, although slightly smaller due to self-gravity effectively increasing the gravitational potential by a factor of 1.1. This leads to a corresponding decrease in outflow velocity. The gas clumps are slightly smaller and denser, making them more resilient to feedback. As a result, the mass of infalling cold gas is slightly higher, and some clumps fall into the SMBH earlier than in the fiducial simulations. Nevertheless, the simulations evolve qualitatively similarly, leaving our main conclusions unaffected.

Fig. 7. Gas density distribution in simulation L2.0-sg at t = 0.75 Myr, directly comparable to the bottom-right panel of Fig. 1.

Some of the dense self-gravitating gas turns into star particles at approximately constant rates¹ throughout the simulations. By t = 1 Myr, the total mass of star particles is M_*/(10⁶ M_⊙) = 1.4, 5.9, 12, and 17 in simulations L0-sg, L0.5-sg, L1.0-sg, and L2.0-sg, respectively. The positive correlation between star particle mass and AGN luminosity qualitatively agrees with many other simulations (e.g. Nayakshin & Zubovas 2012; Gaibler et al. 2012; Silk 2013; Zubovas et al. 2013; Zubovas & King 2014; Schaye et al. 2015; Zubovas & Bourne 2017; Laužikas & Zubovas 2024) and observations (e.g. Cresci et al. 2015; Maiolino et al. 2017; Gallagher et al. 2019; Shin et al. 2019).

The results reveal an interesting relationship between SMBH feeding and star formation. In low-luminosity simulations, star particle masses are insignificant compared to infalling cold gas. However, in the two high-luminosity simulations, total star particle mass exceeds that of infalling cold gas. This suggests that star formation in outflow-compressed dense gas can reduce the rate of SMBH feeding and lead to a smaller final SMBH mass. Such a tradeoff between SMBH and stellar mass growth relates to competitive feedback and accretion in galaxy formation (Nayakshin et al. 2009) and leads to a somewhat counter-intuitive conclusion (see also Silk 2013): galaxies experiencing extreme, perhaps super-Eddington, AGN episodes at early cosmic times may end up with under-massive SMBHs compared to those that grow at more modest rates. Such a connection would be very important for constraining the models of first SMBH growth at redshifts z > 6. It should be checked with dedicated simulations, which we plan to perform in the future.

5.3. Coexistence of hot and cold flows

The hot diffuse gas expands much faster than the cold. When such a system is observed, producing a static snapshot, we expect to see ionised gas at larger radii than molecular gas. This is generally the case when we consider large observed outflow samples, i.e. ionised outflows typically have larger extents than neutral or molecular ones (see the compilations in e.g. Fiore et al. 2017, Lutz et al. 2020, and Fluetsch et al. 2021 and Fig. 2 in Laužikas & Zubovas 2024). The hottest plasma component has recently been probed using both thermal (Hall et al. 2019) and kinetic Sunyaev-Zeldovich effects (Hadzhiyska et al. 2024); it was found to extend farther than the dark matter component of galaxies, which is larger than almost all molecular and ionised outflow radii. In individual objects, the data are scarce, but ionised outflows are wider than molecular in PDS 456 (Travascio et al. 2024), several local ultra-luminous infrared galaxies (Fluetsch et al. 2021), and z ∼ 2 quasars (Vayner et al. 2021, although the trend is not universal, with at least one system showing a much larger molecular outflow extent than ionised).

Marconcini et al. (2025) recently detected evidence of ionised outflow acceleration beyond a few kiloparsecs in ten local AGNs. This is most likely the result of hot gas outflows escaping the gaseous bulge (Zubovas & Tartėnas 2025). Once this happens, the hot gas pressure decreases significantly, leading to less efficient cold gas clump confinement. This can result in filament expansion, which makes them more susceptible to AGN feedback, so SMBH feeding decreases as well (but see Sect. 5.4). Hot gas escape requires that the AGN episode lasts at least t_ep, min ∼ 0.5R_bulge/v_out, hot ∼ 0.5R_kpcv₈⁻¹ Myr, where R_bulge ≡ 1R_kpc kpc is the bulge radius and v_out, hot ≡ 1000v₈ km s⁻¹ is the hot outflow velocity. The factor 0.5 accounts for fossil outflow expansion after the episode (King et al. 2011; Zubovas & Maskeliūnas 2023). This timescale is only slightly longer than the expected typical AGN episode durations (King & Nixon 2015; Schawinski et al. 2015). So the changing influence of the hot gas on the cold filaments can be important for regulating SMBH feeding. In future, we plan to enhance our simulations by coupling AGN luminosity to the feeding rate, investigating the self-regulation aspect.

5.4. Impact of magnetic fields

We neglected the influence of magnetic fields in our simulations. It is known that starburst-driven outflows can lift magnetic fields away from galactic discs (Lopez-Rodriguez et al. 2021), and the same may be true for AGN-driven outflows and magnetic fields in galactic centres. The hot gas bubbles create loops of magnetic field that can compress the dense filaments in between by magnetic pressure. Additionally, magnetic field enhances the survival of cold gas (Shin et al. 2008; Sparre et al. 2020), especially when combined with radiative cooling (Hidalgo-Pineda et al. 2024). Thus, magnetic fields can facilitate SMBH feeding by maintaining narrower and denser filaments for longer. As a result, we expect that in a sample of galaxies with a similar velocity dispersion, those with stronger magnetic fields should have more massive SMBHs. It is known that galaxies with stronger magnetic fields have higher star formation rates (Chen et al. 2025) and higher dynamical masses (Tabatabaei et al. 2016), although correlations with SMBH mass and/or velocity dispersion have not been investigated, to the best of our knowledge. As the number of measurements of galactic magnetic fields grows, this prediction can provide a useful test of our model.

5.5. Comparison with other work

Several previous papers investigated the coupling between multi-phase gas and AGN outflows, generally making similar conclusions to ours. We compare our work with them below.

Zubovas & Nayakshin (2014) showed that dense gas in an outflow receives energy input consistent with momentum driving, while diffuse gas is driven by the whole shocked wind energy. The physical model of feedback used in that paper is equivalent to ours; however, the initial conditions were much simpler, with an initially smooth gas distribution and only an azimuthal density gradient used to create different outflow phases. Furthermore, that paper did not investigate the feeding of the SMBH or gas removal from the nuclear regions.

Bieri et al. (2017) investigated the evolution of outflows in clumpy media and found that the momentum loading factor (in their terms, mechanical advantage) decreases as the outflow escapes from the initial cloud and only the dense clumps remain exposed to the AGN feedback. In this case, feedback was mediated by AGN radiation pressure on dusty gas and the authors did not investigate the effect on feeding the SMBH.

More recently, Ward et al. (2024) ran simulations of AGN wind driving an outflow through a clumpy disc. The underlying physical model is the same as ours. Again, they found that the dense gas is driven essentially by AGN wind momentum, while the diffuse gas forms an approximately energy-conserving outflow. They did not investigate the integrated effect on gas clusters or the effect on SMBH feeding directly. Their results agree very well with ours, despite them using a different wind velocity (10⁴ km s⁻¹ instead of 3 × 10⁴ km s⁻¹), cooling prescription and numerical scheme (moving mesh instead of SPH).

6. Summary

We simulated AGN wind-driven outflows in a turbulent multi-phase medium to analyse the effect of feedback on dense gas and SMBH feeding. Our simulations covered constant AGN luminosities in the range L_AGN = (1.3−32) × 10⁴⁵ erg s⁻¹, corresponding to 0.1–2.5 times the Eddington luminosity for a black hole on the M − σ relation given our galaxy setup, and stochastic variations in initial conditions. We show that dense gas is pushed almost exclusively by the AGN wind momentum, so SMBH feeding is suppressed only when L_AGN ≳ 0.7L_Edd. Since dense gas must be removed far enough to prevent re-accretion after the AGN episode ends, the SMBH mass can grow to or slightly exceed the M − σ relation prediction. Additionally, even at very high luminosities, some dense gas filaments lag behind the outflow, although the total mass of dense infalling gas decreases exponentially with AGN luminosity.

Our results align with numerous recent simulations that reveal a complex picture of multi-phase AGN outflows. A single feedback mechanism–a radiatively driven quasi-relativistic wind, with kinetic power equal to ∼0.05 of the AGN luminosity–produces outflows spanning wide ranges of radii, gas densities, temperatures, and velocities. Outflow properties can also vary substantially due to differences in host galaxy characteristics, particularly ISM clumpiness and turbulence. We intend to explore these multifaceted AGN–outflow connections in future work, building a consistent framework of AGN feedback applicable to galaxies of different masses, morphologies, gas densities, redshifts, and other relevant parameters.

Acknowledgments

MT and KZ are funded by the Research Council Lithuania grant no. S-MIP-24-100. Some simulations were performed on Galax, the computing cluster of the Centre for Physical Sciences and Technology in Vilnius, Lithuania.

References

Appendix A: AGN episode and feedback injection

The AGN affects the gas in two ways: by heating (described in Sect. 3) and by producing feedback in the form of a fast wind, which we tracked using a novel grid-based scheme, gridWind. The gridWind method works by propagating the AGN wind effect on a static grid. This allows the wind to be rapidly coupled with nearby SPH particles through a spatial hash-like method and a sorted distance matrix.

A.1. Feedback injection with gridWind

We used a simplified version of SREAG grid (Malkin 2019) with rectangular cells. To construct the grid, we subdivided the sphere’s surface into latitudinal strips of a constant Δθ = π/N_strip. Each of these strips was then subdivided into a number of rectangular cells, N_Φ(θ):

$$\begin{array}{c}\hfill N_{\mathrm{\Phi }}(\theta )=\lfloor \frac{2\pi cos(b(\theta ))}{\mathrm{\Delta }\theta }\rceil ,\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill b(\theta )=\lfloor \frac{\theta }{\mathrm{\Delta }\theta }\rfloor \mathrm{\Delta }\theta +\mathrm{\Delta }\theta /2-\pi /2.\end{array}$$(A.2)

Here ⌊x⌋ is the floor function and ⌊x⌉ denotes rounding to the nearest integer; θ is the latitude of the midpoint of a given strip. Now the longitudinal size of each rectangular cell in a given strip is Δϕ(θ) = 2π/N_Φ(θ). This results in a deviation of cell area of ±1 % from the mean except for the two polar strips, where the deviation is about 5 %; this is taken into account when distributing feedback instead of modifying the grid itself.

We can quickly determine the indices (n_r, n_θ, n_ϕ) of the nearest grid cell for a given SPH particle at spherical coordinates r_p, θ_p, ϕ_p:

$$\begin{array}{cc}\hfill n_r& =\lfloor r_{\mathrm{p}}/\mathrm{\Delta }r\rfloor ,\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill n_{\theta }& =\lfloor {\theta }_{\mathrm{p}}/\mathrm{\Delta }\theta \rfloor ,\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill n_{\varphi }& =\lfloor {\varphi }_{\mathrm{p}}/\mathrm{\Delta }\varphi (\theta )\rfloor ,\hfill \end{array}$$(A.5)

where Δr, Δθ, and Δϕ are the grid step sizes in the three directions. In practice, we derived a more practical set of indices for use in flattened arrays for each shell (n_shell) and the whole spherical volume (n_sphere):

$$\begin{array}{cc}\hfill n_{\mathrm{shell}}& =n_{\varphi }+{\displaystyle \sum _{k=0}^{n_{\theta }}}{N}_{\varphi }(\theta ),\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill n_{\mathrm{sphere}}& =n_{\mathrm{r}}{N}_{\mathrm{shell}}+n_{\mathrm{shell}},\hfill \end{array}$$(A.7)

where N_shell is the number of cells in each shell. In this work we used N_strip = 64, which corresponds to N_shell = 5216.

We used a simple discrete-step approach for wind propagation. We first injected the wind into the zeroth shell in proportion to the total energy produced by the AGN with luminosity L_AGN over the SMBH timestep; the amount injected into each cell is weighted by its area. We assumed that the wind travels radially outwards at a constant velocity, v_wind = 0.1c (King 2003, 2010a). We equalised the SMBH and wind timesteps by applying an additional constraint, Δt < CΔr/v_wind, where C = 0.4 is a Courant-type factor.

Feedback is distributed to SPH particles in proportion to those particles’ contribution to the overall density field at the centre of the grid cell. As each SPH particle has spatial extent expressed as the smoothing length h, we used a pre-calculated sorted distance matrix to quickly iterate over the nearest neighbouring cells. To prevent unnecessary steps, we checked that the distance from each cell centre is not greater than h. We found that, in practice, the radial extent of each particle over different shells can be safely neglected, greatly improving performance. This is safe as long as we only use the ratio between the contribution of the different particles to the overall density field. Each wind variable - in our case, energy ΔE = η/2 × E_cell and momentum Δp = E_cell/c - is transferred independently. When injecting momentum, the direction of injection is radial if the particle position is inside the cell; otherwise, the direction is parallel to the radial direction of the cell centre.

We tracked the energy injection over the entire run at each injection step in a log file; an example of the output from one L1.0 simulation is shown in Fig. A.1. Comparing the black and dotted red lines, we see an excellent agreement, with relative error approximately constant at 6 × 10⁻⁶. Momentum injection is equally precise because both quantities are proportional to the energy contained within a grid cell, which is the quantity tracked in the plot.

Fig. A.1. Top: Energy injection over time in simulation L1.0. The solid black line shows the expected injection, Etot, exp = LAGN × t. The solid green line shows Ew, the amount of energy contained within the wind that is yet to be injected. The dotted red line shows Einj, the actual amount of energy injected into the particles divided by η/2. Bottom: Relative error, Δinj = |Etot, exp−(Einj+Ew)|/Etot, exp.

A.2. Momentum-only wind test

We illustrate the viability of our approach by performing the momentum-driven wind test as outlined in Nayakshin et al. (2009). In this setup, a luminous source is placed at the centre of a unit periodic box containing isothermal gas of constant density. Feedback is injected spherically into the surrounding medium in the form of outward momentum resulting in an expanding central cavity. Equating the rate of change of momentum of the expanding cavity to the momentum input from the luminous source gives

$$\begin{array}{c}\hfill \frac{d}{\mathit{dt}}\left[\left(\frac{4}{3}\pi R^3{\rho }_0\right)\dot{R}\right]=\frac{L}{c},\end{array}$$(A.8)

which is solved for the radius of the resultant cavity, R(t):

$$\begin{array}{c}\hfill R(t)={\left(\frac{3L}{2\pi c{\rho }_0}\right)}^{1/4}{t}^{1/2},\end{array}$$(A.9)

where L is the luminosity of the central source, c is the speed of light and ρ₀ is the initial density of the surrounding medium. A more complex estimate takes into account the finite velocity of the wind v and the restoring force due to external pressure:

$$\begin{array}{c}\hfill \left[\left(\frac{4}{3}\pi R^3{\rho }_0\right)\dot{R}\right]=\frac{L}{c}(t-\frac{R}{v})-4\pi R^2\rho c_{\mathrm{s}}^{2}t,\end{array}$$(A.10)

where c_s is the speed of sound. This expression is solved for R(t) numerically.

We performed two tests with Gadget, using N_sph = 1 × 10⁶ particles: one with virtual particles (Nayakshin et al. 2009) and one with gridWind. For gridWind we used N_strip = 64, which corresponds to N_shell = 5216 (the same as the main set of simulations) and the number of radial shells N_r = 200, with Δr = 0.0025, which is of the same order as the minimum smoothing length h = 0.001. For virtual particles we used reasonably good resolution with p_γ = 0.5m_sphc_s. The resulting expansion of the central cavity is shown in the top panel of Fig. A.2. After the initial few steps, both approaches converge on the expected values given by Eq. (A.10). In the bottom panel the agreement is highlighted by showing the relative difference from the analytical solution tracing both the density peak and its width. Our method produces a slightly finer peak, although this is more sensitive to the overall SPH simulation resolution. A density map of the spheres in both simulations (Fig. A.3) confirms that gridWind results in a slightly smoother density distribution with a more defined peak at the correct radius.

Fig. A.2. Top: Evolution of the radius of a sphere expanding into uniform isothermal gas due to momentum injection. Bottom: Deviation from the analytical solution. Shaded areas correspond to the width at half maximum from the density peak.

Fig. A.3. Density maps of the resultant spheres. The red dashed line shows the expected radius.

All Tables

All Figures

	Fig. 1. Gas density integrated through a 20° wedge around the mid-plane (z = 0), defined by |z| < 0.18(x2+y2), in simulations L0.3 (top row), L1.0 (middle), and L1.7 (bottom). Columns show t = 0.5, 0.75, and 1.0 Myr from left to right. Brighter colours indicate higher gas densities, in the range ∼10−25 g cm−3 < ρ < 10−19 g cm−3.
In the text

Fig. 2. Kinematic phase diagrams at t = 0.5 Myr in simulations L0.5 (left), L1.0 (middle), and L2.0 (right). In each panel, the top 2D histogram shows radial velocity against density, and the bottom the radial velocity against temperature. 1D histograms at the top and right show distributions of individual gas properties. Grey contours show corresponding distributions in the control simulation. Vertical dotted lines show vr = 0 and vr = −σb; horizontal dotted line shows T = 3 × 104 K.

In the text

Fig. 3. Mass of cold rapidly infalling gas in radial 3.3-pc-wide bins at 0.25, 0.5, and 0.75 Myr (solid red, blue, and green lines and shading, respectively) in simulations L0.5, L1.0, and L2.0 (top, middle, and bottom panels, respectively). Solid lines show the mean of the four stochastically different simulations; shading encompasses their full range. Dotted lines in each panel show equivalent results from the control simulation.

In the text

Fig. 4. Top: Total mass of cold rapidly infalling gas within the central 500 pc as a function of time in the control (red), L0.5 (blue), L1.0 (green), and L2.0 (magenta) simulations. We added the mass accreted by the SMBH particle to this total. Solid lines show the mean of the four stochastically different simulations; shading encompasses their full range. Squares mark the times when the mean mass reaches its minimum value. Bottom: Minimum value of cold rapidly infalling gas mass in all simulations (circles) as a function of luminosity, with colours denoting the time when the minimum is reached. The solid red line is the best linear fit between log M and LAGN of all simulations except the control and L0.1; the dashed lines indicate 95% confidence limits on the line parameters, estimated via bootstrapping.

In the text

Fig. 5. Top: Clusters (coloured by gas density) detected at t = 0.75 Myr in simulation L1.0, plotted over the gas density map (grey scale). Middle panel: Momentum loading factor (fp) distributions (Eq. (9); box and whisker plots) in simulations with different AGN luminosities. Circles highlight outliers beyond 1.5 times the interquartile range covered by the whiskers. The horizontal dashed line highlights fp = 1, which corresponds to pure momentum driving, for ease of comparison. Bottom panels: Momentum loading factor (fp) as a function of dcl (left) and mean gas density (right) in simulations L0.5 (red), L1.0 (green), and L2.0 (blue).

In the text

Fig. 6. Top: Same as the top panel of Fig. 4 but showing the change in SMBH particle mass. The horizontal dashed line shows M = 10mSPH. Bottom: Change in SMBH particle mass between the start of the simulation and t = 0.25, 0.5, 0.75, and 1 Myr (red, blue, green, and magenta circles, respectively), in simulations with different AGN luminosities. Error bars represent the range of values in stochastically different simulations; they are slightly offset horizontally for clarity. Solid lines show the SMBH mass change expected due to accretion powering an AGN of the given Eddington fraction over the same time intervals, assuming a radiative efficiency (η) of 0.1; the shaded region represents the effect of varying radiative efficiency between 0.057 and 0.42. The horizontal dashed line is the same as in top panel.

In the text

	Fig. 7. Gas density distribution in simulation L2.0-sg at t = 0.75 Myr, directly comparable to the bottom-right panel of Fig. 1.
In the text

Fig. A.1. Top: Energy injection over time in simulation L1.0. The solid black line shows the expected injection, Etot, exp = LAGN × t. The solid green line shows Ew, the amount of energy contained within the wind that is yet to be injected. The dotted red line shows Einj, the actual amount of energy injected into the particles divided by η/2. Bottom: Relative error, Δinj = |Etot, exp−(Einj+Ew)|/Etot, exp.

In the text

	Fig. A.2. Top: Evolution of the radius of a sphere expanding into uniform isothermal gas due to momentum injection. Bottom: Deviation from the analytical solution. Shaded areas correspond to the width at half maximum from the density peak.
In the text

	Fig. A.3. Density maps of the resultant spheres. The red dashed line shows the expected radius.
In the text