Role of AGN and star formation feedback in the evolution of galaxy outflows

ABSTRACT

We conducted 3D-magnetohydrodynamic simulations to investigate the feedback processes in the central 1-kpc scale of galaxies hosting both active star formation (SF) and an active galactic nucleus (AGN) wind. Our simulations naturally generated a turbulent and clumpy interstellar medium driven by SF evolution. We found that the AGN-wind duty cycle plays a crucial role in shaping the evolution of the outflows. A single duty cycle (which can repeat several times over the galaxy lifetime) consists of an active, a remnant and an inactive phase, lasting up to 1.5 Myr in our simulations. The duration of the cycle increases with larger star formation rate (SFR) and smaller AGN-wind power (tested for luminosities 10⁴²–10⁴⁴ erg s⁻¹ and SFR = 1–1000 M_⊙ yr⁻¹). The feedback on SF, whether positive or negative, depends on various factors, including the AGN outflow opening angle, power, and phase of activity, as well as the initial SFR. The passage of the AGN wind enhances SF in a ring around it, resembling the structures observed in ULIRGs (Ultraluminous Infrared Galaxies) and LINERS (low-ionization nuclear emission-line region), and is stronger for larger AGN power or SFR. Also, a higher SFR enhances the mixing of interstellar matter with the AGN wind, resulting in a greater number of colder, denser structures with volume filling factors ∼0.02 to 0.12 and velocities comparable to those observed in Seyferts and LINERs, but smaller than those observed in ULIRGs. The efficiency of the AGN wind in transporting mass to kiloparsec distances diminishes with increasing SFR. The mass-loss rates range from 50 to 250 M_⊙ yr⁻¹ within the initial 2 Myr of evolution, which aligns with observed rates in nearby Seyferts and ULIRGs.

1 INTRODUCTION

The impact of supermassive black holes (SMBHs) on the evolution of their host galaxies is known as feedback. The (momentum-driven) jets from active galactic nucleus (AGN) galaxies are generally believed to be the main agents to prevent the cooling of the gas at the large (dozen-kpc) scales where these jets impinge into the warm and hot intergalactic medium (IGM, e.g. Silk & Rees 1998; Netzer 2015; Falceta-Gonçalves et al. 2010). This process, often denominated ‘maintenance mode’, is complemented at smaller (sub-kpc) scales, by the so-called quasar mode where (radiation-driven) outflows are believed to be responsible for sweeping the gas from the galaxy inner radii (e.g. McNamara & Nulsen 2012; Morganti et al. 2021b; Morganti, Oosterloo & Tadhunter 2021a) and even quench star formation (SF, Nesvadba et al. 2010, Wylezalek & Zakamska 2016, Mukherjee et al. 2018b). On the other hand, with the recent increase of higher resolution observations revealing multiple phase structures of gas in different classes of AGN, the aforementioned picture of AGN feedback has become more complex (e.g. Harrison et al. 2020; Scholtz et al. 2020; Morganti et al. 2021a; Aalto et al. 2020; Pereira-Santaella et al. 2018). The contribution in various classes of AGN to the feedback and the scales at which they operate are not yet completely understood (e.g. Murthy et al. 2022). For instance, it is not clear yet whether galactic winds driven by SF (also denominated supernova (SN)-driven winds) play some role in the feedback process too, blurring the effects of the AGN winds (e.g. Melioli & de Gouveia Dal Pino 2015; de Gouveia Dal Pino, Clavijo-Bohórquez & Melioli 2018; Silk & Rees 1998; Netzer 2015; Harrison et al. 2018; Gowardhan et al. 2018; Wylezalek & Morganti 2018; Pereira-Santaella et al. 2021). At the same time, cosmological simulations, specially those focused on the formation of SMBHs (e.g. Barai & de Gouveia Dal Pino 2019) or the production of cosmic rays and diffuse emission (e.g. Hussain et al. 2023), have demanded more realistic constraints on the processes of energy transfer from the AGN and SB winds to the interstellar medium (ISM) and IGM.

The gaseous outflows observed in a high variety of AGN (e.g. Tadhunter 2008; Morganti et al. 2015; Harrison et al. 2018; Veilleux et al. 2020), are often limited to the central kiloparsec or sub-kiloparsec regions of the galaxy and carry relatively small amounts of gas (see e.g. Holt, Tadhunter & Morganti 2008; Rose et al. 2017, 2018; Oosterloo et al. 2017, 2019; Baron & Netzer 2019; Bischetti et al. 2019; Santoro et al. 2018, 2020; Scholtz et al. 2020). Their impact effects may change during the various stages in the life of the AGN (e.g. Morganti et al. 2021a). Previous numerical studies have indicated that the compression by the jet-driven shocks can enhance SF in the disc, but the enhanced turbulent dispersion in the disc can also lead to quenching of SF, as indicated by observations. Whether positive or negative feedback prevails depends on jet power, ISM density, jet orientation with respect to the disc, and the time-scale under consideration (Sutherland & Bicknell 2007; Wagner, Bicknell & Umemura 2012; Melioli & de Gouveia Dal Pino 2015; Mukherjee et al. 2016, 2018a,b; Tanner & Weaver2022).

Another crucial aspect that may largely influence the feedback process is the observed variability in the AGN outflow activity. In optical observations of nearby AGN sources, the analysis of short time-scales reveals that the duration of activity cycles can range from 10⁵ to 10⁶ yr (e.g. Schawinski et al. 2015). Conversely, when observing in the radio band, the complete cycle can span from 10⁷ to 10⁹ yr, leading to the categorization of distinct object groups with differing time-scales. Nevertheless, a harmonious explanation can be found by considering that the shorter-term activity cycles are modulated by the longer-term ones. For instance, the presence of extensive outflows on a scale of approximately 10 kpc or more implies prolonged episodes of activity needed to create and sustain such expansive regions. However, it appears improbable that an SMBH can consistently emit the required luminosity over such extended periods. To reconcile this, it is postulated that multiple short-term activity cycles intervene to drive these massive outflows over such vast distances. This process referred to as the ‘duty cycle’ mechanism (Schawinski et al. 2015; Zubovas & King 2016; Morganti 2017a,b; Zubovas 2018), determines a sort of AGN-galaxy co-evolution, and though it is not yet fully understood, it is widely accepted that the gas content plays a crucial role in triggering both AGN and SF activity.

In this work, in order to improve our understanding on the formation of these outflows and on how they impact on the galaxy evolution at the inner regions in different stages of the AGN activity, we present three-dimensional (3D) MHD (magnetohydrodynamic) simulations of the sub-kpc-scale core region of active galaxies with intense SF. Distinctly from previous works that have simulated the interaction of an AGN jet with a clumpy ISM with large filling factor fractal distribution of high-density clouds (e.g. Wagner & Bicknell 2011; Wagner et al. 2012; Mukherjee et al. 2016, 2018a,b; Murthy et al. 2022; Tanner & Weaver 2022), we here consider fiducial values of SFR and SN rates (SNR) in an initially smooth multiphase galactic disc in order to allow for natural formation of a turbulent and structured ISM where the AGN jet is then injected. A galactic wind driven by SF is also naturally produced in this environment and we explore the interactions and relative role of both the AGN- and the SF-driven outflows on the feedback processes, including the temporary suppression of gas infall due to the AGN outflow. We also explore the effects of initial conditions, by launching the AGN outflow in different stages of the galactic ISM evolution.

This work extends upon previous one by Melioli & de Gouveia Dal Pino (2015), where the combined effects of SF and AGN outflows were explored in the context of Seyfert galaxies and the main conclusion was that the SF is the main responsible for the mass loading in the outflows. The AGN jet alone would be unable to drive a massive gas outflow, at least in this class of systems, though it may have power enough to drag and accelerate clumps to very high velocities.¹ We quantify, in particular, observed quantities like the mass outflow rate, the filling factors of different gas components, and the duty cycle of the AGN activity, and also compare our results with previous simulations and several observed classes of active galaxies, including young radio galaxies with molecular and ionized outflows (e.g. Morganti et al. 2013, 2015; Dasyra et al. 2016; Oosterloo et al. 2017, 2019; Maccagni et al. 2018; Schulz et al. 2018, 2021; Murthy et al. 2019, 2022; Morganti et al. 2021a), Ultraluminous Infrared Galaxies (ULIRGs, Sturm et al. 2011; Dasyra & Combes 2012; Veilleux et al. 2017; Cazzoli 2017; Pereira-Santaella et al. 2018, 2021; Chen et al. 2020; Herrera-Camus et al. 2020; Ramos Almeida et al. 2022; Fluetsch et al. 2021), and Seyfert galaxies (Fluetsch et al. 2019; Tombesi et al. 2012, 2015, 2017).

This paper is organized as follows: in Section 2, we describe our numerical model. Section 3 is devoted to describe and analyse the results of our simulations. In Section 4, we discuss and compare our results with previous studies and with observations, and finally in Section 5, we draw our conclusions.

2 NUMERICAL MODEL

We integrate numerically the MHD equations given by,

where ρ is the density, |$\boldsymbol v$| is the velocity, B is the magnetic field, p* is a diagonal effective pressure tensor with components p* = p + B²/2, and E is the total energy density given by,

where p is the thermal pressure and γ is the ratio of the specific heats. We assumed an equation of state for an ideal gas, p = (γ−1)e, where e is the internal energy density and γ = 5/3. The source term |$\dot{m_S}$| represents the rate of mass injection by supernovae, SNe Ia and SNe II, f_G represents the external gravitational forces due to dark mater, the bulge, and the disc of the galaxy, and E_S = E_SNIa + E_SNII + E_AGN + E_DM + E_disc + E_bulge−Λ(ρ, T) represents the sources of energy density coming from SN type Ia and SN type II explosions, the AGN wind, and the external potentials due to the dark matter halo, the disc, and the stellar bulge as well as the radiative cooling losses, respectively. In order to integrate these equations numerically, we used the Godunov-based code amun (Kowal, Lazarian & Beresnyak 2007),² which was modified in order to introduce the effects of radiative cooling over a wide range of temperatures, 50 K < T < 10⁹ K (Raga 2000; Melioli & de Gouveia Dal Pino 2004; Schure et al. 2009; Townsend 2009; Melioli & de Gouveia Dal Pino 2015). We considered an optically thin gas in ionization equilibrium with a metal abundance in the mid-plane given by z/z_⊙ = 1.0 and decreasing linearly with the height of the disc up to a minimum value of z/z_⊙ = 0.1, and also included molecular hydrogen H₂. For integration in space we used the Harten–Lax–van-Leer C (HLLC) Riemann solver for the HD simulations, and the HLL D (HLLD) for the MHD simulations (see Toro 2009). A second-order Runge–Kutta scheme (Liu & Osher 1998) was employed for time integration. The radiative cooling function Λ(ρ, T) is calculated implicitly in each time-step and implies errors smaller than 10 per cent for the gas cooling down to a minimum (floor) temperature. In all our simulations, we have considered a minimum temperature T = 50 K (see Section 3).

2.1 Boundary conditions

The computational domain has physical dimensions 1 × 1 × 1 kpc³ in the x-, y-, and z-directions for most of the models with a collimated AGN wind, and 1 × 1 × 2 kpc³ dimensions for the models with a spherical (360°) AGN wind as well as for a few models with 10° opening angle (see more details in Section 2.2.2). We have run most of the models with 256³ resolution (corresponding to a cell size ∼3.9 pc). Simulations performed by Melioli & de Gouveia Dal Pino (2015) for a similar domain considering resolutions between 128³ and 512³ have indicated a good degree of convergence in the results for 256³ resolution. For this reason, we have adopted this resolution, that is, a cell size ∼3.9 pc, in this work for most of the runs, but have also run a few models with resolution twice as large, with cell size ∼1.95 pc. Our simulations have outflow boundaries conditions, which ensure that the fluid reaching the boundaries escapes from the system without coming back.

2.2 Initial setup

The initial setup is also similar to Melioli & de Gouveia Dal Pino (2015), except for the introduction of the magnetic field into the system.

The initial mass for the galaxy includes the contribution of a bulge and a stellar disc. The rotating ISM gas in the disc is set in equilibrium with the galaxy gravitational potential, Φ(r, z), given by the sum of the dark matter halo, the bulge, and the disc contributions (see equations 8–10 of Melioli & de Gouveia Dal Pino 2015). We determine the disc rotation velocity by solving the equation v(r, z) = (rd(r, z)/dr)^0.5 (also referred to in Melioli et al. 2008, 2009; Melioli, de Gouveia Dal Pino & Geraissate 2013; Melioli & de Gouveia Dal Pino 2015). The velocities obtained through this approach closely resemble the rotation curve typically observed in a spiral galaxy. The galaxy parameters adopted in this study are presented in Table 1.

We consider an initially stratified multiphase IS gas distribution, with each i-phase defined by a typical temperature T_i, which is initially in hydrostatic equilibrium in the gravitational potential as described above. The initial gas density distribution for this smooth setup is shown on the left panel of Fig. 1. Also, we consider initially four typical temperature phases: very cold gas with a temperature of T₁ = 50 K, cold gas with a temperature of T₂ = 200 K, hot gas with a temperature of T₃ = 10⁴ K, and ionized gas (extragalactic halo) with a temperature of T₄ = 10⁶ K.

In order to account for the presence of magnetic fields in our simulations, we assume the fluid to be initially embedded in a horizontal magnetic field either constant or stratified assuming an initial constant ratio of the thermal to magnetic pressure β = P_th/Pm. With this prescription, equation (11) of Melioli & de Gouveia Dal Pino (2015) is modified as follows:

where |$c^2_{s,i}$| is the isothermal sound speed of the i-phase of the gas and e_i quantifies the fraction of rotational support of the ISM (see Table 1). In our model, the i-phase density is replaced by the (i+1)-phase density wherever the (i+1)-phase total (thermal+magnetic) pressure is larger than the (coldest and densest) pressure of the i-phase.

For models where we assumed initially constant magnetic field, we adopted only one value B₀ = 0.76 µG. In the runs with stratified magnetic fields, we considered three initial values of β = p/(B²/8π) = ∞, 300, and 30, corresponding to initial magnetic fields at the plane of the disc, B₀ = 0.0 (HD models), 7.6, and 24 µG, respectively (where B₀ is the magnetic field strength at the galactic plane).

These values are compatible with observed values in nearby Seyfert galaxies and ULIRGs which point to average values of ∼1–100 µG, with the latter value occurring only in starburst (SB) systems with very large SFRs (Thompson et al. 2006; Drzazga et al. 2011). They are also comparable with those adopted in previous 3D-MHD studies (Rodenbeck & Schleicher 2016; Ntormousi 2018).

2.2.1 SF wind parameters

SF drives turbulence and the development of a clumpy ISM. In order to account for this and also for the possibility of occurrence of SF-driven outflows in the central regions of the galaxy, we have considered the feedback from massive stars in the late stages of their lifetime when they explode as type II SNe (SNe II), and also from low-mass white dwarfs in binary systems which explode in type Ia SNe, particularly in the bulge.

The rate of SN Ia explosions is proportional to the bulge mass, |$\text{SNIa} \sim 0.01({M}_{{\mathrm{ bulge}}}/10^{10}\text{M}_{\odot })\text{yr}^{-1}$| (Pain et al. 1996) which in turn is proportional to the mass of the central SMBH (⁠|${M}_{\text {BH}}$|⁠) given by the relation |${M}_{\text {BH}}/ {M}_{\text {bulge}} \sim 10^{-3}$| (Häring & Rix 2004), so that we can obtain this rate in terms of the SMBH mass (Melioli & de Gouveia Dal Pino 2015):

To introduce SNe II, we consider an SB in the central region of the host galaxy distributed in the disc. Assuming the Salpeter initial mass function, the number N of stars with mass larger than 8 M_⊙ (i.e. the SNe II number) is N ∼ 0.01(M_b/M_⊙), where M_b is the total mass of the stars in the SB. The SNe II activity lasts up to t_b ∼ 3 × 10⁷ yr, which is the lifetime of an 8 M_⊙ star. An effective average SN rate is thus given by SNR ≃ N/t_b, which is in good agreement with more accurate evaluations given by Leitherer et al. (1999). In order to inject the SNe in space and time, each ith-type 1 and type 2 SN explosion is randomly associated with a position P_i = (r_i, h_i) in the bulge (in the case of SNIa) and in the disc (in the case of SNII) at a random time t_i, within the interval 0 < t_i < 30 Myr. This methodology generates explosions with a spatial frequency that aligns with the stellar density of the bulge (following a Plummer profile) and the disc (encompassing a radius of approximately 300 pc and a typical half-height, above and below, between 80 and 160 pc, depending on the specific model under consideration). We have assigned SNe events with an average luminosity in accordance with the explosion rates specified by each model. More details can be found in Melioli et al. (2009) and Melioli & de Gouveia Dal Pino (2015).

We have constrained the initial amount of gas in the disc by the star formation rate (SFR) using the Kennicutt–Schmidt relation, which matches observations for a wide range of galaxies (Kennicutt 1998):

where SFRD = SFR/A is the star formation rate density, A is the disc area, and Σ_gas = N_H × m_H is the gas surface density. Therefore, for a given disc area A, both the SFR and the SNR depend only on the disc column density, that is,

where N_H,23 is the column density in units of 10²³ cm⁻². Consequently, using the relationships presented above, each model is characterized by a specific initial column density to which it is associated a given SFR and SNR. We explore here models with column densities between 10²³ and 10²⁵ cm⁻³ to which we have matched SFR between 1 and 1000 M_⊙ yr⁻¹.

2.2.2 AGN-wind setup

We assume an AGN outflow driven by the SMBH activity in the centre of the galaxy. We consider three different geometries for this injection: (i) a collimated outflow with 0° opening angle, (ii) a conical outflow with an opening angle of 10°, and (iii) a spherical outflow injection. In the cases (i) and (ii), a total luminosity (or mechanical power) L1 = 3.5 × 10⁴² erg s⁻¹ or L2 = 2.35 × 10⁴³ erg s⁻¹ is injected symmetrically in a volume of ∼(7.8 pc)³, which corresponds to eight cells in the centre of the computational domain (four cells above and four cells below the mid-plane of the galaxy) characterized by an injection velocity v_AGN perpendicular to the plane of the disc and a total rate of injected matter |$\dot{M}_{\text{AGN}}$|⁠, wich are given in Table 2. In the case (iii), we also consider two possible values for the total injected power L2 = 2.35 × 10⁴³ erg s⁻¹ or L3 = 2.35 × 10⁴⁴ erg s⁻¹, introduced in the system as thermal power, which represents a good physical approximation to an adiabatic (energy-driven) propagation of the AGN wind from inner parsec scale (Zubovas & King 2012). In both cases, the wind velocity is given by|$v_{\mathrm{ AGN}}=(2 L/\dot{M}_{\text{AGN}})^{1/2}$|⁠.

2.2.3 AGN-wind duty cycle setup

In our simulations, we have accounted for the fact that the propagation of the AGN wind can eventually prevent the inflow of interstellar (IS) gas to the core, hence interrupting temporarily the AGN activity and the outflow itself. As stressed in Introduction, AGN flux variability is observed at different time-scales (e.g. Schawinski et al. 2015). This is most probably related to variations in the nuclear activity due to the accretion flow (e.g. Ho 2008). In our simulations, we have looked at the time at which the total mass in the central region (within a radius ∼ 40 pc) decreases to 1 per cent of its initial mass, due to the passage of the AGN wind. Then, we have turned-off the AGN wind at this time in order to follow the evolution of the system over a single cycle of the outflow under this condition. As we will show in Section 3.2.2, this ‘duty cycle’, which can repeat for several times over the system lifetime, consists of the active phase (before the outflow is turned-off), followed by a remnant and an inactive phases.

We should note that this turning-off condition is influenced by the resolution of our calculations. Given cell sizes ranging from 7 to 14 pc, we are unable to model the interruption of the inflow to the nuclear region very accurately. The inflow of very small dense structures into the core, as the ones detected, for instance, in Cielo et al. (2018) studies, might enable the persistence of extremely narrow jets in the central region. Nevertheless, we contend that these mechanisms would not suffice to maintain the broader, less collimated AGN winds discussed in this work.

2.2.4 Initially smooth and clumpy environments

The AGN wind is injected in the system considering two different initial setups. In a family of models denominated smooth (S), the wind is injected when the system is in an initially stratified state, in hydrostatic equilibrium (see left panel of Fig. 1). In this case, the AGN wind propagates at the same time that the outflow driven by SF with SN explosions starts to develop. In another family of models denominated clumpy (C), the AGN wind is launched when the system has already developed a clumpy state due to SF and the SNe-induced outflow (see right panel of Fig. 1).

The initial conditions and parameters of the models analysed are provided in Table 2. The model labels use S or C to indicate initially smooth or clumpy environments, respectively; HD represents pure hydrodynamical models, MHD models with an initial constant magnetic field (B₀ = 0.76 µG), and MHD β_number denotes MHD models with initial constant β = 30 and 300. The number SF stands for the SFR in units of M_⊙ yr⁻¹. The subsequent label in the model specifies the AGN-wind opening angle (0°, 10°, or 360°); and the last label indicates the AGN-wind luminosity: L1, L2, or L3, corresponding to L_AGN = 3.5, 23.5, or 235 × 10⁴² erg s⁻¹, respectively. Finally, the only model presented in the table that was run at higher resolution is labelled as high-res. We have grouped the models into two samples. One includes the models with SFR 1 M_⊙ yr⁻¹, named Low SFR, most of which have collimated AGN winds, that is, with 0° or 10°, characteristic of Seyfert-like systems. This sample has been subdivided in two categories smooth and clumpy models. The second big sample (labelled High SFR) includes all models with higher SFR and spherical AGN-wind injection, which are typical characteristics of ULIRG-like systems. All the models in this group are initially clumpy.

3 RESULTS

In this section, we discuss the main results or our study, focusing on the feedback of the stellar activity and AGN outflows into the 1 kpc³ central region of the host galaxy over the first 2.0 Myr of evolution. In particular, we will analyse the influence of the initial conditions in the evolution of the galaxy gas.

3.1 Convergence of our simulations

Let us start discussing briefly the convergence of our results. Fig. 2 compares the density, temperature, and velocity maps for model C.MHD.1000SF.360°.L2 of Table 2 (i.e. with an initial SFR = 1000 M_⊙ yr⁻¹, an initial constant horizontal magnetic field with value 0.76 µG, a spherical (360°) AGN wind and luminosity L2) for two distinct resolutions 512² × 1024, corresponding to a cell size Δl = 1.95 pc in each direction, and 256² × 512, corresponding to a cell size Δl = 3.9 pc in each direction. Fig. 3 compares the evolution of volume-averaged quantities for these two models.

We note that the general features of the evolving system are comparable in both resolutions, as well as the volume averaged quantities. We will discuss all these properties in more detail in the next sections, but the results in these figures reveal good convergence of the results already for a resolution of 256² × 512, and for this reason, most of the models considered in this work were run with this resolution. The higher resolution maps in Fig. 2 evidence the natural presence of more features, particularly in the temperature and velocity maps in the halo regions, but the macrodynamical properties are essentially preserved, as we see also in Fig. 3. We only note a slight delay in the formation and the strength of the peak of the average velocity, which is slightlty higher and is attained earlier in the higher resolution run. This is due to the AGN-wind injection which encounters more resistance to break into the smoother (less porous) lower resolution environment, but this will not affect the overall evolution of the system.

3.2 Effect of ISM initial conditions on AGN-wind propagation

3.2.1 Smooth versus clumpy environments

In this section, we show how the porosity of the ISM induced by SF affects the feeding and the dynamics of the AGN wind as it evolves. Porosity is one of the most important characteristics of this study where the multiphase ISM is not induced by an imposed initial density distribution function, but is instead, a natural consequence of the stellar feedback. It is important to emphasize that the porosity of the ISM does not depend only on the stellar feedback, but more generally on gas turbulence (see e.g. Barreto-Mota et al. 2021) and also on gravitational instabilities and fragmentation (Toomre-like instabilities) in the disc. Nevertheless, stellar activity and in particular SN explosions play a key role in the generation of the turbulence, and gravitational disc instabilities act at larger gas mass structures and on longer time-scales than those investigated in this study (see e.g. Elmegreen 2009). For this reason, we believe that in this study the stellar feedback is the most important parameter in regulating the porosity of the medium.

The initial gas distribution determines the evolution of the system. The main differences between an initially smooth and an SF-induced clumpy environment are shown in Fig. 4. Here, we compare density, temperature, and velocity distributions, at t ≃ 300 kyr, of two systems where a cone-shaped AGN wind with 10° opening angle is injected in an initially smooth system in hydrostatic equilibrium and in an evolved clumpy ISM, both characterized by an initial horizontal magnetic field B = 0.76 µG and an SFR = 1 M_⊙ yr⁻¹. We clearly see that in the smooth environment, the AGN outflow expands faster, easily disrupting the disc, reaching the halo, and forming a low density, hot regular bipolar-shaped wind. On the other hand, in the case of the clumpy environment, an asymmetric, wiggling outflow develops on both hemispheres, propagating at slower velocity through the irregular structures that form in the thicker disc due to SF activity. In this case, the AGN outflow barely reaches the halo. The SF-driven galactic wind is in turn much broader, but much slower. This can be particularly distinguished in the density and temperature diagrams associated to the expanding SN bubbles arising from the disc in the clumpy system. In the smooth system, this SF-induced wind is much less developed at the time depicted because the SNe have just started to explode in this case. Also striking is the fact that when the disc gas is smooth and not dominated by stellar feedback (left panel), the AGN wind expands through the galaxy disc without modifying it significantly and no clouds (or cold gas) are transported above the disc. On the contrary, when the AGN wind impinges directly into the inhomogeneous, clumpy ISM, formed as a consequence of the stellar feedback (right panel), the AGN wind has more difficulties in expanding in the disc due to the action of external ISM pressure gradients. These include all the outside ISM material which is interacting with the AGN and also the SB wind material. The gas evolution is thus dominated by the stellar activity, and high-density, low-temperature structures, like filaments and clouds, are more easily pushed into the halo, both by the SB and by the AGN wind.

SFR (and the related gas column density of the environment) is therefore, a fundamental parameter in these systems because it determines the efficiency of the stellar feedback and the consequent multiphase structure of the gas in the disc. A galaxy with low density and thus negligible stellar activity will produce no significant amount of clouds or filaments in the disc and, consequently will hardly produce, for instance, cold (molecular) gas outflows. For this reason, in the analysis that follows we consider, in general, reference models from Table 2 with a galactic disc setup in which the stellar feedback has already created a clumpy environment and a multiphase ISM.

3.2.2 Duty cycle of the AGN activity

Both, the stellar feedback and the AGN wind inject into the ISM of the disc a large amount of energy which accelerates a fraction of the gas to the halo of the galaxy, and part of it returns, but to the external regions of the disc. As stressed in Section 2.2.3, after a given time, which depends on the initial conditions of the system and the amount of energy injected, the core region of the galaxy in our simulations will lose almost completely its gas mass, forcing temporarily an interruption of the active phase of the AGN. What is evident from our simulations is that the active phase of the AGN cannot be continuous, but must follow alternating cycles of outflow activity and inactivity.

Indeed, as the AGN wind sweeps through, the IS gas is completely expelled from the core (as we can see in Fig. 6). However, after some time of the inactive phase of the AGN wind, the IS matter is able to recommence its inflow towards the central region (Fig. 6). This sets the stage for a new active phase of the AGN.

Fig. 5 illustrates this duty cycle for the higher resolution model C.MHD.1000SF.360°.L2.high-res. The figure depicts maps of density, temperature, and velocity for three different snapshots of the AGN activity. At t = 300 kyr (active phase), the AGN has disrupted the disc, pushing IS gas from the core, and reached the halo forming a spectacular hot bubble that would be probably observable in X-rays. After a ∼350 kyr active period, the AGN wind is turned-off because the inflow of ISM fuel to the nuclear region has been interrupted by the AGN outflow itself. The panel at t = 450 kyr shows the remnant phase of the wind, which is followed by a fully inactive phase (starting after t ∼ 525 kyr and ending around t ∼ 1125 kyr) where there is no more trace of the AGN-wind remnant, only a broader outflow in the halo which is mainly driven by the SN bursts. At this stage, the inflow of fuel to the core from the active SF of the ISM is resuming and a new AGN-wind phase should start.

Now, let us examine in more detail the evolution of the gas mass within the central region of the galaxy, in order to determine a more accurate estimate of the duration of the active phase of the AGN in our models.

Fig. 6 shows the time evolution of the average mass fraction (M_core(t)/M₀) in the central region of the galaxy, for the smooth model S.MHD.SF.10°.L1, where M_core(t) gives the gas mass within a radius of 40 pc, and M₀ = M_core(0) is its initial value (see Table 2). The model depicted by the green line illustrates an interrupted AGN wind after a duration of ∼2 × 10⁵ yr, whereas the black line represents a continuous injection of the AGN wind. In the case of continuous AGN wind, the gas within the central region is entirely expelled after approximately 2 × 10⁵ yr. The high pressure exerted by the AGN wind prevents the inward flow of IS gas into the central region, thus interrupting the fuel supply to the SMBH and subsequently halting the AGN outflow. This suggests that a continuous AGN wind injection beyond approximately 2 × 10⁵ yr is no longer realistic. Consequently, we conducted the same model with a constrained active phase duration of the AGN wind at around 2 × 10⁵ yr (green line).

In the absence of AGN-wind energy injection, the gas stellar feedback generates an inflow of gas that gradually fills the central region of the galaxy. Eventually, the replenished gas reaches a maximum value, which may trigger a new phase of AGN-wind activity. This pattern exemplifies the duty cycle of the AGN wind, which, depending on the initial conditions of each model, may exhibit slight variations in duration. However, our simulations consistently indicate a range of a few 10⁵ yr for both smooth and clumpy environments (see Table 2).

We note that these values are consistent with observed variability phenomena in such systems (see e.g. Schawinski et al. 2015 and references therein).

In view of the results above, we focus in this work on models in which the injection of the AGN wind is turned-off when the fuel provided by the IS gas inflow is interrupted by the wind itself. This allows us to follow at least one duty cycle of the AGN feedback into the galaxy. An illustration of this is shown in Fig. 7 for several initially clumpy models of Table 2 with different SFR. The core mass evolution is followed until the end of an AGN-wind cycle, that is, right before a new cycle of gas inflow to the nuclear region and a new activity of AGN wind should resume. In the figure, we also highlight (in the vertical blue, grey, yellow, and brown lines) the characteristic time (t_active) at which the AGN wind is turned-off due to the gas emptying in the central region. These values are also listed in Table 2.

It is evident that the duration of AGN-wind activity exhibits a clear correlation with the SFR, reflecting the amount of gas present in the disc or the column density, denoted as N_H. Additionally, we observe that the active phase of the AGN decreases, as expected, with the power of the AGN wind. For instance, when considering an SFR of 1000 M_⊙yr⁻¹, we note that the increase of the AGN-wind luminosity by a factor of 10 (from L2 to L3) leads to a 49 per cent decrease in the AGN activity time, reducing it from 375 to 190 kyr. In contrast, for the same luminosity L2, decreasing the SFR by a factor of 10 (from 1000 to 100 M_⊙yr⁻¹) results in a reduction of the activity time from 375 to 225 kyr, amounting to a 40 per cent decrease. Fig. 7 also shows another important result. The continuous lines highlight the gas evolution in the central region for models hosting stellar feedback only, with no AGN wind. Although stellar feedback plays a key role in the creation of a multiphase medium in the disc of the galaxy (as shown in Section 3.2.1) and is, therefore, essential in promoting the formation of dense, filamentary structures (see discussion in the next section), it is evident that its action alone results to be much less efficient in removing gas from the central region of the galaxy, taking much longer time for that, of at least ∼3 Myr. For this reason, the duty cycle in our models is clearly more correlated with the active phase of the AGN wind.³

3.2.3 AGN-wind opening angle

We are also interested in understanding the role that the AGN wind plays on the evolution of the host galaxy. For this reason, we have considered three different setups for the AGN wind: (i) a collimated outflow perpendicular to the disc of the galaxy with 0° opening angle; (ii) a collimated outflow perpendicular to the disc of the galaxy with 10° opening angle; and (iii) and a spherical wind corresponding to an opening angle of 360° (which resembles observed ultra fast outflows – UFOs).

The simulations with 0° and 10° opening angle revealed no substantial differences, a result which is also sustained in higher resolution simulation (performed for model C.MHD.1SF.10°.L1) with resolutions 256² × 512 and 512² × 1024, compared with its 0° counterpart at resolution 256³).

However, it is essential to conduct a more meticulous analysis of the differences between collimated and spherical wind models. To achieve this, we examine two specific models from Table 2: C.MHD.SF1.10°.L2 and C.MHD.SF1.360°.L2, both with identical initial conditions, except for the wind opening angle, which is 10° and 360°, respectively. Fig. 8 shows the evolution of the density, temperature, and velocity, for these two models. During the first 775 kyr, the energy deposited by the AGN wind heats the ISM, which in turn starts to expand doing work on the circumnuclear gas. For the 360° AGN wind, the working surface is larger and thus it drags more material during the expansion. At the end of the activity phase (at t = 0.225 Myr), the remnant of the 10° AGN wind propagates rapidly through the disc and is deposited in the halo region, mixing with the hot gas. The same occurs for the spherical AGN wind, but at smaller speed, and it mixes with the hotter gas in the halo region only after about 1 Myr. The reason why this happens is that, unlike the 10° AGN wind, for which the energy is deposited mostly in the halo region, part of the energy of the spherical AGN wind is used to push the dense gas in all directions of the disc plane.

3.2.4 Star formation rate

In this section, we explore the evolution of the gas of the galaxy as a function of the SFR, in order to understand how it influences the formation of high-density structures, such as clouds and filaments, in the outflow. The relation between the amount of gas column density available in the disc to form stars and the star formation rate is well known: SFR |$\propto N_\mathrm{ H}^{1.4}$| (Kennicutt 1998; see section 2.2.1). As remarked, we have performed numerical simulations considering SFR ranging from 1to 1000 M_⊙yr⁻¹ (Table 2). In order to ensure the correspondence with the Kennicutt–Schmidt law (Kennicutt 1998), the central gas density in these models varies between 10²–10⁴ cm⁻³, respectively.

In Fig. 9, we show the gas density distribution evolution for three models with spherical AGN-wind ejection with different initial SFR = 1, 100, and 1000 M_⊙yr⁻¹. We note three main differences. First of all, the higher density of the disc in the model with higher SFR naturally slows down the expansion of the AGN wind. This effect is also enhanced by the higher pressure of the ISM induced by the higher number of SN explosions, directly connected with the SFR. The second difference is the larger number of dense structures formed in the disc in this case. This characterizes how the porosity (or ‘clumpiness’) of the disc, increases with the SFR (and N_H) and will be studied in detail in Section 3.4, where we evaluate the filling factors of the multiphase components of the gas, but in Fig. 10, we can already distinguish the larger number of colder, denser structures in the model with larger SFR. Lastly, the third distinction pertains to the evacuation process of the central region of the galaxy. The larger density and SFR facilitate the inward movement of IS gas towards the nuclear region, consequently delaying the emptying process and promoting a more effective mixing of the ISM with the AGN wind. This phenomenon also accounts for the longer active phase of the AGN wind observed in models with higher SFR (see Section 3.2.2 and Fig. 7; see also Table 2).

Fig. 11 shows the column density maps integrated along the normal direction to the disc for the same models of Fig. 9. The column density, which is a quantity that can be directly observed, obviously blurs the structuring features we see in the density maps. None the less, we can still appreciate the intrinsic differences due to the increase of the SFR, particularly in the disc region.

SFR (and the associated column density) is, therefore, one of the most sensitive parameters of our models, and for this reason in Section 3.4 we will focus mainly on this parameter for the study of the evolution of the gas of the galaxy and in the comparison with observations in Section 4.1.

3.2.5 HD versus MHD models

In this section, we compare HD and MHD models to better understand the role that the presence of magnetic fields can play on the evolution of the system, specially in SF-induced clumpy environments.

In Fig. 12, we compare the evolution of volume-averaged quantities of the system, such as velocity (⁠|$\bar{v}$|⁠), density (⁠|$\bar{n}$|⁠), temperature (⁠|$\bar{T}$|⁠), magnetic field intensity (⁠|$\bar{B}$|⁠), and thermal pressure (⁠|$\bar{p}$|⁠).

For the set of initial horizontal stratified magnetic fields considered here, with β = 30, 300, and ∞ (the latter corresponding to the hydrodynamical system), or an initial constant field with B₀ = 0.76 µG, the comparison of the average quantities indicates no substantial differences between the models. In fact, the system with an initial constant B₀ = 0.76 µG, which is the value adopted in most of the models of this work (see Table 2) has essentially the same average density, temperature, pressure, and vertical velocity of the HD counterpart. The same applies to the stratified model with β = 300, which has a maximum initial magnetic field 10 times larger at the disc. Some minor differences on the average velocity, density, temperature, and pressure become more noticeable in the model with β = 30, which corresponds to a maximum initial magnetic field B₀ = 24 µG at the disc. The average total velocity decreases from 5000 km s⁻¹ in the B₀ = 0.76 µG and HD models, to less than 2000 km s⁻¹ in the β = 30 model, caused by the deceleration of the AGN wind due to the impact into the magnetized gas which absorbs part of the kinetic energy of the terminal shock. The same applies to the average temperature and thermal pressure which decrease with increasing magnetic field (decreasing β), due to the absorption of part of the internal energy density in the compression of the magnetic field. Additionally, we have found that when we intensify the initial magnetic field to B₀ = 76 µG in the disc (β = 3), the aforementioned distinctions become significantly more pronounced, especially after the impact of the AGN wind, reflecting tangible dynamical influence of the magnetic fields in the system’s evolution. While we have not presented the results of these stronger magnetized models in this study, we intend to delve into them in a forthcoming work.

Figs 13–15 confirm the trends above revealed by the average quantities. Figs 13 and 14 compare central 2D slices of the magnetic field intensity and density distribution, respectively, for the models with initial stratified magnetic fields β = 300 and 30, and the HD model (β = ∞). Fig. 15, compares 2D histograms of the number density versus vertical velocity, number density versus temperature, and temperature versus vertical velocity of the entire gas distribution for the same models. Clearly, no substantial differences are seen among these models, except for the occasional presence of small filaments and clumps with larger-than-average magnetic field strength accumulated there by compression, specially for the model with β = 30. In this case, there is a slightly larger quantity of high velocity (up to few ∼100 km s⁻¹), colder and denser structures (see Fig. 15). In summary, the results above indicate that magnetic fields with average values up to a few tens µG, as observed in most of the sources investigated here, have no major effect on the evolution or on the formation of high-velocity large dense, cold structures in these systems.

3.3 Gas mass outflow

As discussed previously, the expulsion of gas from the galactic disc towards the halo is driven by both stellar feedback (the SB-driven wind) and the AGN wind. This evolutionary progression is highly influenced by the initial conditions of the system (refer to Section 3.2) and the radiative cooling of the gas. Conversely, as emphasized in Introduction, observations reveal the existence of dense and cold (molecular) gas structures with high velocities outside the discs of numerous active galaxies, indicating that the AGN wind predominantly contributes to the presence of such features at elevated altitudes. Consequently, our focus lies in quantifying the mass-loss rate within our models to gain a deeper understanding of which combination of initial conditions and physical processes can effectively describe and substantiate these observational outcomes.

In Fig. 16, we show the mass-loss rate from the disc region (|z| < 200 pc) for the models labelled as low SFR in Table 2.

First and foremost, these plots provide compelling evidence for a bimodal evolution of mass outflow, contingent upon whether the galactic disc exhibits an initial smooth or clumpy gas distribution. It is clearly evident that the gas mass outflow is considerably smaller in the smooth models. Conversely, no significant differences are observed between models with AGN-wind opening angles of 0° or 10°, as well as their counterparts without AGN wind (SB curves). However, as expected (refer to Section 3.2.3), a higher mass-loss rate is observed in the case of a 360° AGN wind.

Models with stellar feedback+AGN wind characterized by both larger SFR (and thus larger gas column density) and larger opening angle, present the most influential effects on the gas outflow and structure formation. This is shown in Fig. 17 for the set of models labelled High SFR in Table 2, which are also compared with their counterpart with SFR = 1 M_⊙yr⁻¹. For these models, the total gas mass lost by the central region of the disc (within R = 200 pc) varies between 10  per cent (for SFR = 100–1000 M_⊙ yr⁻¹) and more than 40  per cent for SFR ≤ 10 M_⊙ yr⁻¹, of the initial total mass of that region. The smaller percentage for larger SFR is due to the comparatively larger total amount of gas in the system for increasing SFR. A more complex, nonlinear behaviour is denoted in the corresponding mass-loss rates shown in the bottom diagram. Models with same AGN-wind luminosity, but different initial SFR reveal no clear trend, particularly the model with SFR = 100 M_⊙ yr⁻¹, which has a comparatively smaller rate than its counterpart models. This can be attributed to the fact that for small SFR, there are less SN explosions and structure formation to prevent the outflow expansion to the halo, while for a large enough SFR, as in the case of SFR = 100 M_⊙ yr⁻¹, these obstacles become larger causing the decrease of the mass-loss rate, but if one further increases the SFR (to SFR = 1000 M_⊙ yr⁻¹), then the expansion of the structures due to much larger number of SN explosions will instead help the outflow propagation thus increasing again the mass-loss rate (as seen in the bottom panel), evidencing a clear cooperation of both the SN-driven and the AGN winds. This is what we observe from an inspection of the maps of outflow evolution in these models. Moreover, the mass-loss rate is generally smaller in models where only the SB wind is present, compared to their counterparts where the AGN wind is also considered. The increase of the AGN-wind power also affects the mass-loss rate from the disc. This can be observed by comparing the two models in Fig. 17 with SFR = 100 M_⊙ yr⁻¹ and AGN-wind powers L2 (dashed orange curve) and L3 (dotted orange curve). Notably, the model with larger power has a mass-loss rate twice as large.

The same applies to the models in the Fig. 17 with SFR = 1000 M_⊙ yr⁻¹. Overall, the models indicate mass-loss rates from the disc between 50 and 250 M_⊙ yr⁻¹.

We also notice that for the magnetic field strengths investigated here, there is no significant influence on the mass outflow.

The results above are consistent with those of the previous sections, that is, they confirm that the models which are able to transport a large amount of gas to outside the disc and allow for the formation of dense structures in the halo are those wherein there is a combined action of AGN and stellar feedback, being sensitive to both the SFR, and the power and opening angle of the AGN wind, though the dependence on the SFR reveals a more complex trend, that we will try to better disentangle in the next section.

3.4 Filling factor

The analysis above of the gas mass outflow is not enough to understand the evolution of the gas above the disc, and more specifically, its multiphase composition and velocity. Therefore, in order to trace the formation of structures in the outflow, we compute the volume filling factor f_V given by the ratio of the volume V′ occupied by a given set of N′ structures with a density above a given threshold, to the total volume V of the system (V′ ≤ V). The total volume is given by the box of our computational domain and considering that the volume of the cell is v_cell, we have

For this study, we have considered three different gas number density (n) ranges, namely: (low) f_l = f_V(n < 1 cm⁻³); (medium) f_m = f_V(1 cm⁻³ < n < 100 cm⁻³), and (high) f_h = f_V(n > 100 cm⁻³). The volume filling factor for the models with spherical AGN wind of Table 2 which differ mainly in their SFR, are depicted in Fig. 18.

The multiphase composition of the gas influences the evolution of the outflow and the appearance of dense clouds outside the disc. In previous studies (e.g. Wagner & Bicknell 2011; Asahina, Nomura & Ohsuga 2017; Cielo et al. 2018; Mukherjee et al. 2018a, b; Tanner & Weaver 2022), the distribution of the multiphase medium was chosen as a free parameter in the initial conditions, whereas in this study this distribution is a direct consequence of stellar feedback. As shown in Fig. 18, this distribution evolves with time and changes completely depending on the initial conditions of the system as well as the interactions between stellar and AGN wind. For models with SFR higher than 1, the AGN wind reduces the fraction of low-density gas, n < 1 cm⁻³ (Fig. 18, top) and together with the continuous action of the stellar feedback, it increases the fraction of medium, 1 < n < 100 cm⁻³ (Fig. 18, middle) and high-density gas, n > 100 cm⁻³ (Fig. 18, bottom).

It is interesting to compare our results in Fig. 18 with those of Wagner & Bicknell (2011). In their models, they imposed two fixed values for the initial filling factors, namely 0.13 and 0.42, for clumps with densities in the range 30–1000 cm⁻³, distributed at heights of the order of 500 kpc above the disc. In our models, where the filling factor is naturally built by SF, we obtain initial values that are at most ∼ 0.12, for the same range of density structures, distributed up to a ∼ 200 kpc height, at the initial states of the clumpy models. As time evolves, this quantity decreases to values of the order of 0.01–0.02.

While Fig. 17 lacks information on how the amount of high altitude dense gas in the galaxy relates to the SFR, Fig. 18 provides more detailed information. As expected, the concentration of dense material will be much more significant for SFR = 100 and 1000 M_⊙yr⁻¹, than for SFR = 1 or 10 M_⊙yr⁻¹, as clearly indicated in the bottom diagram of the Fig. 18. In fact, after about 2.5 Myr, the filling factor of the gas with number density ≥ 100 cm⁻³ above the disc reaches values between 11 per cent and 2 per cent for SFR = 1000 and SFR = 1 M_⊙yr⁻¹, respectively. The intensity of the stellar feedback is therefore, one of the main ingredients to load the galactic outflow with cold (molecular) gas. We also note that the β = 30 model improves slightly the filling factor of the high-density component by at most a factor ∼ 1.7 with regard to the β =300 counterpart.

The study of the volume filling factor of the gas is also particularly important when we analyse the temperature distribution of the gas above the disc. To do this, we have classified the gas in the following ranges: cold: (50–200) K, warm: (200–1500) K, hot: (1500, 10⁴) K, ionized: (10⁴–10⁶) K, bremsstrahlung: (10⁶–10⁹) K, and res: (residual) >10⁹ K. In Figs 19 and 20, we explore the mass fraction distribution above the disc region (|z| > 200 pc) for the two models that, according to Fig. 18 are able to transport cold dense (molecular) gas to outside the disc more efficiently, namely, C.MHD.100SF.360° and C.MHD.1000SF.360°.

More specifically, in Fig. 19, we compare the multiphase gas distribution both, in the absence and presence of the AGN wind (left and right panels, respectively), calculated for the model with SFR = 100 M_⊙yr⁻¹. The figure shows not only the volume filling factor (top), but also the mass filling factor (bottom panels). At t = 0 Myr, the gas distribution above the disc is the same in both cases and the most massive component is the hot gas (1.500–10.000 K) followed by the ionized component, although their volume filling factor is |$\sim 5~{{\ \rm per\ cent}}$|⁠. During the system’s evolution, it becomes evident that the SB wind alone, as well as the combined influence of the SB and AGN wind, effectively generate and carry a substantial amount of cold gas beyond the disc’s boundaries. Nevertheless, these two winds fulfill distinct roles in the gas evolution. The SB wind primarily contributes to filament formation, whereas the AGN wind predominantly facilitates the transportation of denser structures located above the disc. Consequently, in the absence of the AGN wind in the model, the proportion of cold gas transported above the disc after 2 Myr amounts to approximately 30 per cent. Conversely, when the AGN wind is included in the model, this percentage increases to around 70 per cent.

Similarly, Fig. 20 depicts the multiphase gas distribution for the model with SFR = 1000 M_⊙ yr⁻¹, but considering two different injection luminosities for the AGN wind. Clearly, the model with 10 times larger luminosity is able to remove a grater amount of cold gas, approximately 2.7 times larger, corresponding to a total mass of ∼ 1.75 × 10⁷ M_⊙.

Moreover, we have seen in Section 3.2.5 that the AGN wind is able to accelerate cold gas to high velocities. Cold structures with a few 100 km s⁻¹ are identified for AGN-wind luminosities L2 (see Fig. 15) and should increase further with increasing luminosity. On the other hand, when comparing the model C.MHD.1000SF.360.L2 in Fig. 20 with its counterpart in Fig. 19 differing only by the SFR = 100 M_⊙ yr⁻¹, we observe that the cold gas above the disc after 2 Myr increases from 70 per cent to about 90 per cent from the lower to the larger SFR (and thus denser) model.

In Fig. 21, we show the evolution of cold gas above the disc in models with different magnetic field strength: β = 300 and 30, depicted in the left and right panels, respectively. It is worth noting that a stronger magnetic field presence enhances the formation of cold gas. At t =0.75 Myr, the mass filling factor is approximately 15 per cent in contrast to the quasi-hydrodynamic model, which exhibits a value of around 5 per cent. However, by the end of the simulation, the magnetic field’s increased intensity does not yield a significantly larger quantity of cold gas. Consequently, we can infer that within the range of magnetic field values explored in this study, where thermal pressure surpasses magnetic pressure throughout the entire domain, the distribution of cold gas remains unaffected.

Fig. 22 shows the filling factors for two models with AGN-wind opening angle 10° (left) and 360° (right), in different times. As time goes by, the fraction of cold gas at high altitudes is clearly larger in the 360° model, as expected. The net amount of mass transported from the disc for the cold gas component in the 360° case is |$\Delta {M}^{360^{\circ }}_{\mathrm{ cold}} \simeq 1.55\times 10^7\text{M}_{\odot }$|⁠, while it decreases with time for the 10° case (by |$\Delta {M}^{10^{\circ }}_{\mathrm{ cold}} \simeq -2.04\times 10^6\text{M}_{\odot }$|⁠).

Figs 23 and 24 summarize well the results presented above, for the final stages of the systems. Fig. 23 exhibits face-on column density maps obtained by excluding the disc within different ranges: ±200, ±500, and ±750 pc, from the top to the bottom panels, respectively. The maps provide a comparative analysis of the impact of two systems: one featuring only the SB wind (model C.MHD.100SF, depicted on the left), and the other incorporating both the SB and AGN winds (model C.MHD.100SF.360°.L2, shown on the right). Essentially, these maps are like a ‘tomography’ of the region above and below the disc, showing the quantity and morphological characteristics of the denser (molecular) structures transported to varying heights above the disc. It is evident that both models, with either the SB wind alone or the combination of SB and AGN winds, are capable of transporting dense gas to heights of approximately 1 kpc. However, the presence of the AGN wind noticeably enhances both the quantity and the maximum height of these structures.

Fig. 24 presents similar column density maps, but for two models that differ only on the AGN luminosity. Clearly, the more luminous AGN wind can account for the transport of more massive structures to the highest altitudes, outside the disc.

We note that investigating the evolution of the filling factor distinguishes our approach from models that assume fixed initial conditions and makes direct comparisons of probability density functions and power spectra, such as those aforementioned works (e.g. Mukherjee et al. 2016, 2018a, b), or studies devoted to ISM modelling (e.g. Barreto-Mota et al. 2021). In fact, we believe that tracking the filling factor, specifically exploring how stellar feedback and AGN wind influence the formation of cold clouds, is crucial for understanding the transport and evolution of these clouds above the galactic disc. Consequently, in this study, we have preferred to maintain a more pragmatic approach, associating distinct cloud distributions with varying SFRs (or column densities), as depicted in the diagrams illustrating the filling factor’s evolution over time.

4 DISCUSSION

The primary aim of this study was to determine the conditions under which an AGN wind can produce gas outflows that agree with observed data. To achieve this, the study sought to identify the essential mechanisms responsible for feedback and for generating the observed cold (molecular) gas above the disc, as well as the initial conditions that must be considered in theoretical models to provide a justification for the observations. Key observables, such as the mass of the gas outflow, the mass of the extraplanar cold (molecular) gas, the SFR of the galaxy, and the luminosity of the AGN wind, were measured in a large sample of recent observations, providing a macroscopic yardstick for selecting the best and most realistic models from those under consideration here.

It is important to remark that this study has been built upon the previous work by Melioli & de Gouveia Dal Pino (2015), which examined the combined effects of SF and AGN outflows within the context of Seyfert galaxies. The primary conclusion of that study was that SF plays the dominant role in mass loading the outflows. The AGN jet alone appeared insufficient to drive significant gas outflows, at least in the case of this class of systems. However, it did demonstrate the capability to entrain and accelerate clumps to very high velocities.

In the present investigation, we have expanded the scope considerably. We have considered a broader range of AGN luminosity, spanning from 10⁴² to 10⁴⁴ erg s⁻¹, a wider range of disc gas mass, ranging from 10⁸ to 5 × 10⁹ M_⊙, core gas mass varying from 10⁵ to 10⁷ M_⊙, and surface number density between 10²³ and 10²⁵ cm⁻². Consequently, the SFR spans from 1 to 1000 M_⊙ yr⁻¹. Additionally, in this study, we have examined both collimated jets and spherical winds. The latter exhibit characteristics akin to UFOs, which are regarded as some of the most potent disc-driven winds in AGN. Furthermore, unlike Melioli & de Gouveia Dal Pino (2015), all the simulations conducted in this work are MHD, encompassing magnetic fields ranging from 0 to 24 μG.

4.1 Comparison with observations

In the upper part of Table 3, we summarize the main properties of the gas outflow obtained from our models with a special focus on the cold (molecular) gas component (50–200 K) and its kinematics. In the lower part of the table, for comparison, we depict the same quantities for some selected nearby sources which show active star-forming regions and AGN winds and which have, in general, similar values to those computed from our models (Pereira-Santaella et al. 2018; Fluetsch et al. 2019), see also Pereira-Santaella et al. (2020, 2021), Veilleux et al. (2020), and Fluetsch et al. (2021).

In Fig. 25, we plot these quantities, namely, the mass-loss rate |$\dot{M}_{\mathrm{ out}}=M_{\mathrm{ out}}/t_{\mathrm{ dyn}}=M_{\mathrm{ out}}/(R_{\mathrm{ out}}/v_{\mathrm{ out}})$| (left), the momentum rate |$\dot{P}_{\mathrm{ out}}=\dot{M}_{\mathrm{ out}}v_{\mathrm{ out}}$| (middle), and the luminosity (mechanical power) |$L_{\mathrm{ out}}=\dot{M}_{\mathrm{ out}}v^2_{\mathrm{ out}}/2$| (right), where R_out is defined as the outflow distance from the centre, and we compare them with a larger sample of observed sources. To calculate these quantities from our simulations, we compute the outflow gas mass, M_out, with temperatures between T = 50 and 200 K, and vertical velocities greater than 100 km s⁻¹. Then, the outflow velocity v_out is calculated from the average of the vertical velocity of the cold structures. For most of the sources used for comparison here, the outflow distance R_out ranges from ∼ 100 to ∼ 400 pc, the velocity v_out from 100 to 1000 km s⁻¹ and the corresponding outflow dynamical time-scale, t_dyn = R_out/v_out ranges from 10⁵ to 10⁶ yr.

The grey symbols in Fig. 25, represent observed sources, while the colored ones, the results from our models. Diamonds stand for systems with AGN winds, while stars, for systems with SF wind only. The blue symbols stand for models with L_AGN = L2 = 2.35 × 10⁴³ erg s⁻¹, while the red and green diamonds represent models with higher AGN-wind luminosity, L_AGN = L3 = 2.35 × 10⁴⁴ erg s⁻¹ (and SFR = 100 and 1000 M_⊙ yr⁻¹, respectively).

The analysis of Fig. 25 points to four main results which are discussed below.

The impact of the AGN wind helps to release the gas heated by SF in the disc. The dashed lines in the left panel of Fig. 25, correspond to the mass loading factor, |$\eta = \dot{M}/\text{SFR}$|⁠. The leftmost blue diamonds in this panel, which represent systems with AGN wind with luminosity L_AGN = 2.35 × 10⁴³ erg s⁻¹, and low SFR (SFR = 1), have a mass loading factor η ∼ 10, while similar systems with SF wind only (blue stars), have η < 1. This means that the power generated by the SB wind alone is not enough to produce such a mass transport rate, as expected from our previous results.

Furthermore, we also see that maximum loading factors are obtained for the systems with AGN wind, and SFR = 1 M_⊙ yr⁻¹. For an AGN wind with same luminosity (L2 = 2.35 × 10⁴³ erg s⁻¹), but with larger SFR, the mass loading factor decreases from η ∼ 10 to ∼ 7, for SFR = 1 to 10 M_⊙ yr⁻¹, respectively, and for systems with SFR = 100 and 1000 M_⊙ yr⁻¹, η < 1. In other words, the AGN, looses efficiency in carrying out mass to kpc distances in systems with higher SFR and column densities, suggesting that for more massive systems, the mechanical power of the AGN should be higher, as indicated by the observations. In fact, the mass loading factor can be as high as η ∼ 4 for a larger luminosity AGN, as indicated by the red diamond in the left panel, which has a wind luminosity L3 = 2.35 × 10⁴⁴ erg s⁻¹. The green diamond model, with same luminosity, but larger SFR = 1000 M_⊙ yr⁻¹ also has a higher mass loading than its lower luminosity counterpart (blue diamond), but still smaller than 1.

The dashed lines in the right panel of Fig. 25 give the percentage a of the SFR luminositiy (L_SNe) that is converted into outflow luminosity, denoted as L_{out, SN} = a × L_SNe. We observe that L_out is approximately 10 per cent L_SNe, when injected by low SFR (low density) systems, and less than 1 per cent for higher SFR (denser) systems. Similarly to the mass loading factor, L_out also increases with the AGN-wind luminosity (as we see when comparing the lower luminosity blue diamonds with the higher luminosity red and green ones).

Finally, in Fig. 25 (left panel), we can see that there are several sources with the same SFR, but different mass-loss rates. This could be due to different properties of their environments, characterizing different stages specially of the AGN evolution and activity. In fact, although we have plotted in the figure, only the maximum values of the mass-loss rates obtained for the simulated models, we have seen in Fig. 17 that, as the AGN wind evolves with time in its duty cycle of activity, |$\dot{M}$| also increases, suggesting a vertical motion upwards in the |$\dot{M}$| diagram of Fig. 25.

Table 3 also gives values of the average velocity reached by clouds and filaments above the disc (cold structures with temperature between 50 and 200 K with a vertical velocity of at least 100 km s⁻¹) obtained from our models. These velocities are between about 100 and 270 km s⁻¹, which are values similar to those observed in some sources such as NGC 1266, NGC 4418, IRAS 15115 + 0208, M51, and Circinus Galaxy, which are highlighted in the bottom of the table and are classified as Seyferts or LINERS (low-ionization nuclear emission-line region, Fluetsch et al. 2019). Although in our simulations we have identified few cold structures with maximum velocity of several 100 km s⁻¹ (see e.g. Fig. 15), the average values obtained in the table are lower than those observed for the sources classified as ULIRG, as reported in Pereira-Santaella et al. (2018). These average values also differ somewhat from those obtained by, for example, Mukherjee et al. (2018a, b), where they detected dense clouds dragged along the jet with outflow velocities of ∼ 500 km s⁻¹.⁴ Overall, our results suggest that only a very few dense structures can accelerate to the observed velocities in ULIRGs, of about 500 km s⁻¹. However, we should remind that in our simulations the formation of cold structures in the disc is not imposed by the initial conditions, but instead, is naturally driven by the combination of the SF and radiative cooling in the ISM of the nuclear region of the galaxy, tested for different SFRs. So, what should be missing in our models? First of all, we have not considered an important ingredient, namely, the presence of dust. The radiative cooling down to temperatures of dust formation may help to improve and prevent the evaporation of several cold structures. Moreover, clouds may be further accelerated by the radiation pressure on dust relative to a hot, diffuse background, as shown in a series of numerical simulations (Sharma & Nath 2012; Thompson et al. 2015; Yu et al. 2020). Also, larger magnetic fields than those investigated in this work can help to preserve a larger amount of cold structures from evaporation (see e.g. Farber & Gronke 2022). Finally, higher resolution simulations could also allow for the formation of thinner shock front shells that could, therefore, accelerate cold filaments (that survive to evaporation) to high velocities, in situ at the high altitudes (Viegas & de Gouveia dal Pino 1992). All these effects will be investigated in forthcoming work.

Our simulations also indicate that the mechanism responsible for the acceleration of cold clouds is not only due to the mechanical interaction between the jet and the ISM. Of course, high SFR (or column density) is able to produce a large fraction of cold (molecular) structures, and the action of the AGN wind is essential to remove these clumps from the disc. Therefore, SF and AGN wind have to coexist to generate the observed molecular outflows. However, the velocities reached by the clouds are in general lower than the escape velocity (of about 700 km s⁻¹), hence, eventually most of these structures will fall back onto the disc, as we see in the simulations.

4.2 Comparison with other simulated models

The intense stellar feedback considered in our simulations is the main source of turbulence in the ISM of the disc, but further turbulence may also be driven transversely by the passage of the AGN wind. Both may contribute either to enhance SF or to quench it (e.g. Schawinski et al. 2006; Mukherjee et al. 2018b; Barreto-Mota et al. 2021), depending on features like the initial SFR (and thus the column density) and the AGN power, implying positive or negative feedback, respectively.

In order to compare our results with previous simulations (e.g. Mukherjee et al. 2018b), in Fig. 26 we show the star formation rate surface density (SFRD) distribution in the disc relative to the initial value for the models C.MHD.100SF.360° for two AGN-wind luminosities: L₂ = 2.35 × 10⁴³ erg s⁻¹ (top) and L₃ = 2.35 × 10⁴⁴ erg s⁻¹ (middle), for different times of evolution from the active to inactive phases of the AGN duty cycle. In these diagrams, SFRD is normalized to SFRD0, the star formation rate density at t = 0. SFR(t) is calculated from the Kenniccut–Schmidt law (equation 8), using the same normalization as in equation (3) of Melioli & de Gouveia Dal Pino (2015), SFRD ∼ 10⁻¹(n₁₀h_SB,200)^1.4 M_⊙yr⁻¹kpc⁻², where n₁₀ is the density of the disc at z = 0 in units of 10 cm⁻³ and h_SB,200 is the height of the SB region in units of 200 pc. The cavity in the middle of the panels is due to the passage of the AGN wind and it is larger for the AGN wind with larger luminosity (middle panel). The bottom diagram shows how the same system behaves when there is no AGN outflow. We see that in all diagrams, there is a decrease in the SFRD with time caused by the passage of both outflows. While the SF wind acts more smoothly over all radii, the AGN wind causes the formation of a ring of enhanced SF around the cavity during the active period (t ∼ 150 kyr), which disappears later on when the AGN wind disappears completely in the inactive phase (t ∼ 600 kyr). This result can be compared with that found by Mukherjee et al. (2018b) for a similar SFR ∼ 80 M_⊙yr⁻¹, but with collimated jet (see their model H). In their case, they detect the formation of a stronger ring of enhanced SF around the AGN cavity, and decreasing SF in the outer regions. Remarkably, a recent observation by the JWST–NIRSpec (Near-Infrared Spectrographer) IFS (integral field spectroscopy) of the luminous infrared galaxy NGC 7469, also classified as a Sey 1.5 nucleus, has revealed an enhanced star-forming ring surrounding its bright AGN (Bianchin et al. 2024). It has inner and outer radii of 330 pc and 616 pc, which resembles the rings in Fig. 26, though localized slightly further away from the nucleus. This difference could be due to intrinsic differences in the source properties.

In Fig. 27, we have plotted the same as in Fig. 26, but for the set of models with SFR 10 times larger. In this case, the SFR density surface is clearly larger over the entire system, as expected. The ring of enhanced SF is also much stronger and persists even when the AGN has completely disappeared in the inactive phase, though less intense.

In Fig. 28, we have plotted the evolution of the SFR surface density for models with a much smaller initial SFR (1 M_⊙yr⁻¹) and different AGN-wind opening angle, namely, the conical C.MHD.1SF.10° and the spherical C.MHD.1SF.360° models, both with luminosity L2, which are compared with the counterpart model with no AGN wind. In this case, we also clearly notice the formation of the strong ring around the cavity with enhanced SF as time evolves, reflecting a positive feedback of the AGN on SF. Its counterpart model with opening angle 360° (bottom panel), shows a much stronger effect with the AGN-wind passage. The cavity is comparatively much bigger and the ring of enhanced SF around the wind is stronger.

Therefore, we conclude that the positive or negative feedback are highly sensitive to a non-linear combination of all the effects, that include the SF wind, the AGN outflow opening angle and power, the initial SFR (and column density), and the AGN-wind phase of activity. The passage of the SF wind causes a decrease in the SFRD with time over the entire disc, and this effect is more pronounced for larger SFR. The increase of the AGN power tends to enhance SF in a ring around the AGN wind that tends to disappear once the AGN wind faints, and this effect is also stronger for larger SFR. The wind opening angle also affects the SF feedback. A spherical injection creates a much larger cavity, but also a more enhanced ring of SF around the cavity.

In addition to earlier pioneering studies (e.g. Sutherland & Bicknell 2007), numerous recent investigations have explored the impact of AGN-wind feedback on the evolution of the ISM and SFR. Some of these studies, conducted on a larger scale than the current work, have produced results that generally align with our findings and offer complementary insights. For instance, Cielo et al. (2018) demonstrated that inclined jets can decrease star-forming gas. They observed that AGN-wind action removes the densest 20 per cent of the cold ISM by heating clump cores, generating substantial amounts of warm gas with densities ranging from 10 to 200 cm⁻¹, potentially reducing global SF. In contrast, hydrodynamical simulations by Bieri et al. (2016), revealed that the pressure enhancement from AGN action triggers instabilities leading to increased fragmentation compared to simulations without AGN influence. In such cases, enhanced fragmentation results in the formation of more clumps and higher SFR values. It is important to note that these simulations, conducted on a time-scale of 400 Myr and a size of 30 kpc, analyse galaxy evolutionary histories on scales significantly different from those explored in our study. Additionally, the influence of AGN wind on SFR evolution may vary depending on the distance from the galaxy’s centre. For instance, large-scale numerical simulations by Gaibler et al. (2012) indicated that the onset of jets increases SFR on the disc by a factor of 2 or more, simultaneously reducing SF in the central region. Generally, SFR at intermediate radii does not further increase and may slowly decrease as denser clumps are converted into stars. At larger radii, SFR continues to rise and is expected to remain high until the cocoon pressure of the jet approaches the original ambient pressure, as also observed by Gaibler, Krause & Camenzind (2009). In contrast, in galaxy mergers, the simulations conducted by Talbot, Sijacki & Bourne (2024) suggest a moderate suppression of SF throughout the entire merger process due to the presence of the AGN wind. The peak of the SF, achieved during the final coalescence of the galaxies, is reduced by nearly a factor of 2.

5 CONCLUSIONS

In this study, we have presented 3D MHD simulations of the core region of active galaxies at sub-kpc scales. Fiducial values of SF and SNR rates were considered to allow for the natural formation of a turbulent and structured ISM where the AGN jet was injected. In this environment, a galactic wind driven by SF is also naturally produced, and the interactions and relative roles of both the AGN and the SF-driven outflows on the feedback processes were explored. The study also investigated the effects of initial conditions, including the intensity of mild magnetic fields, the intermittency of the AGN activity, the AGN outflow opening angle, the SFR in the disc, and the porosity generated in the ISM by stellar feedback. The underlying idea of this study was that the AGN wind alone would be unable to drive a massive gas outflow in the class of systems investigated, but it may have enough power to drag and accelerate clumps to very high velocities. The results of this study are summarized below.

• The initial gas distribution determines the evolution of the system. In a smooth environment, the AGN outflow expands faster, easily disrupting the disc, reaching the halo, and forming a low density, hot regular bipolar-shaped wind. On the other hand, in the case of the clumpy environment, an asymmetric, wiggling outflow develops on both hemispheres, propagating at slower velocity through the irregular structures that form in the thicker disc due to SF activity.

• After being continuously injected, the AGN wind completely removes the gas from the central region in approximately a few 10⁵ yr, in all models investigated (e.g. Fig. 7). Moreover, the high pressure of the AGN wind against the IS gas interrupts the gas fuelling into the core which in turn, is responsible for the AGN outflow. This implies that the continuous injection of the AGN wind beyond about a few 10⁵yr is no longer realistic. Without the pressure and energy injected by the AGN wind, gas stellar feedback is able to resume gas inflow that again refills the central region of the galaxy. Over time, the replenished gas reaches a maximum value that could eventually trigger a new phase of AGN-wind activity. This evolution characterizes the duty cycle of the AGN. The entire cycle that includes the active, the remnant, and the inactive phases lasts at most ∼1.5 Myr in our models, being longer for larger SFR (larger column density) and decreasing AGN-wind power (Fig. 7). This is consistent with observed short term variabilities in active galaxies (Schawinski et al. 2015; Zubovas & King 2016; Morganti 2017a, b; Zubovas 2018). It is important to stress that this duty cycle pertains to a single phase within the entire AGN outflow activity. We expect a reactivation of this cycle once the remnants of the wind no longer impede the gas re-inflow into the core, thereby initiating a new duty cycle. It is worth noting that the formation of very large-scale jets is expected to result from multiple cycles of activity like the ones we describe.

• The impact of a collimated AGN wind (with an opening angle of approximately 0°–10°) on the evolution of the galaxy disc in the nuclear region is significantly less pronounced compared to a spherical injection. In the case of a spherical injection, a larger working surface results in a greater amount of material being dragged to the halo, albeit at a slower expansion rate (see e.g. Fig. 8).

• In the magnetized models, there is some deceleration of the AGN wind due to the impact into the magnetized gas which absorbs part of the kinetic energy at the terminal shock. However, for the magnetic fields we considered, with average values up to a few tens µG, as observed in most of the sources investigated here, we find no substantial effects in the macroscopic properties of the flow, or on the evolution and formation of high-velocity dense, cold structures.

• The AGN wind looses efficiency in carrying out mass to kpc distances in systems with increasing SFR (and thus column densities). For instance, for an AGN-wind luminostiy ∼10⁴³erg s⁻¹, the total gas mass lost by the disc (|z| < 200 pc) varies from 10 per cent for SFR = 100–1000 M_⊙ yr⁻¹, to over 40 per cent of the initial total mass for SFR ≤ 10 M_⊙ yr⁻¹, but the increase in the AGN-wind power by a factor 10 improves the mass-loss rate by about a factor 2. On the other hand, models with SB wind only, generally have smaller mass-loss rates than their counterparts with AGN wind present (e.g. Fig. 17). We find mass-loss rates of the order of 50–250 M_⊙ yr⁻¹ in the first 2 Myr of evolution, for SFR in the range 1–1000 M_⊙ yr⁻¹ and luminosities in the range 10^42–44 erg s⁻¹, which are comparable to observed values in nearby Seyferts and ULIRGs.

• Overall, we have found that the positive or negative feedback on SF is highly sensitive to a non-linear combination of various effects, including the AGN outflow opening angle, power, and phase of activity, as well as the initial SFR (and column density). The passage of the SF wind results in a decrease in the SFR surface density over the entire disc with time. Increasing the AGN power enhances SF in a ring around the AGN wind, which disappears once the AGN wind weakens, and this effect is more pronounced for larger SFR. The wind opening angle also affects the SF feedback, with a spherical injection creating a larger cavity and an enhanced ring of SF around it than a collimated wind (see Figs 26–28).

• Finally, our models indicate that a higher SFR (column density) in the disc favours a better mixing of the IS matter with the AGN wind increasing the porosity of the disc and providing larger number of colder, denser structures (Figs 9–11, and 19 and 20). The cold clouds and filaments, with temperatures between 50 and 200 K, transported to distances above the disc between 200 pc and 1 kpc, reach average velocities of approximately 100 to 270 km s⁻¹. These velocities are comparable to those observed in several Seyferts and LINERS (e.g. reported in Fluetsch et al. 2019). However, we have identified only very few cold structures with maximum velocities high enough to be compared to those observed in ULIRG systems, which are several hundred km s⁻¹ (e.g. Pereira-Santaella et al. 2018). Our average values differ from those obtained in other studies such as Mukherjee et al. (2018a, b). It is important to note that in our simulations, the formation of cold structures in the disc is naturally driven by the combination of SF and radiative cooling in the ISM of the nuclear region of the galaxy for different SFRs, instead of being imposed by initial conditions. In future work, we will include important missing ingredients such as the presence of dust and higher magnetic fields that may help to improve the formation and survival of high speed cold structures at high altitudes.

ACKNOWLEDGEMENTS

The authors acknowledge support from the Brazilian Funding Agency FAPESP (grants 13/10559-5 and 2021/02120-0). WECB also acknowledges support from the Brazilian Agency CNPq, and EMdGDP from CNPq grant (308643/2017-8). The simulations presented in this work were performed in the cluster of the Group of Plasmas and High-Energy Astrophysics (GAPAE), acquired with support from FAPESP (grant 2013/10559-5), and the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was also made possible by FAPESP (grant 2009/54006-4) and the INCT-A.

DATA AVAILABILITY

The simulated data generated during this study are available upon request to the authors.

Footnotes

REFERENCES