<figure><img src="image_15.png" /><figcaption>**Figure 14.** Mass-weighted 2D histogram of velocity versus density at 0.85$t_{cc}$ for self-gravity simulations with different wind power. The physical time for each simulation is indicated in the legends. The contours in each panel represent the initial distribution. The vertical dashed line indicates the gas density of 10 cm$^{-3}$, while a horizontal line marks $v_x = 100$ km s$^{-1}$.</figcaption></figure>

### 3.3.2 Multiphase structure

To analyse the gas phases of the outflowing material, we have constructed phase diagrams representing gas temperature versus number density in Fig. 15, where the colourbar in the top and bottom rows show the total mass and the mass-weighted mean velocity of the cloud material, respectively. From the top row, it is evident that most of the cloud material is in the cold phase ($T \lesssim 10^3$ K) for the winds with lower powers, whereas there exists a notable amount of moderately dense gas ($10 \lesssim 100$ cm$^{-3}$) in warm ($10^3 \lesssim T \lesssim 10^4$ K) and hot ($10^4 \lesssim T \lesssim 10^7$ K) phases in the higher power cases. Cloud materials in these phases experience acceleration up to 400 km s$^{-1}$. Interestingly, comparatively diffused gas with densities (0.1–10 cm$^{-3}$), which notably contributes to the ionized outflows, exists mostly in the hot phase, with velocities up to several 1000 km s$^{-1}$. On the other hand, gas in the dense and cold phase exhibits very low velocities, less than 100 km s$^{-1}$, even in the highest power case. This suggests that accelerating dense, cold interstellar gas to high velocities ($\sim 1000$ km s$^{-1}$) via AGN-driven winds is exceedingly challenging.

In the simulations featuring lower power winds (GC42 and GC43) shown in Fig. 14, the presence of high-velocity gas is notably limited, due to the low ram pressure of the wind, which primarily accelerates the cloud.

### 3.3.3 Outflow rate scaling with wind power

In order to quantify the outflow driven by the winds and its dependence on the wind power, we calculate the total mass-outflow rates ($\dot{M}_{OF}$) of the cloud material as a function of time for self-gravitating simulations with different wind power, including the simulation with the kinetic wind (GC45_kinetic), which is illustrated in the left panel of Fig. 16. The outflow rates are calculated across the y-z plane at $x = 75$ pc. In the global picture, this corresponds to the outflow rate through an area of $A = 100 \times 100$ pc$^2$ from a single cloud situated at a distance $R_{OF} \sim 1.075$ kpc from the AGN as per our configuration. Notably, the mass outflow initiates after the compression phase when the ablation and stripping of the cloud start to dominate the evolution and progressively approaches a near-steady state value over time. The amount of outflow increases with increasing wind power, as expected. Furthermore, during the steady-state phase, the outflow rate induced by a kinetic wind (dashed-dotted magenta) surpasses that of a thermal wind (solid magenta) with equivalent power due to the higher momentum transfer efficiency. It is important to note

that, the mass-outflow rates, calculated at $x = 75$ pc, are entirely contributed by the stripped and entrained materials from the cloud.

To establish the relationship between the mass outflow rate and wind power, we determine the average mass outflow rate ($\dot{M}_{OF,avg}$) in each simulation after $\dot{M}_{OF}$ reaches a near-steady state, as indicated by the dotted lines (red for thermal and magenta for kinetic wind) in the left panel, which is plotted against the wind power in the right panel of Fig. 16. The open squares correspond to $\dot{M}_{OF,avg}$ values for self-gravitating simulations with different powers, as indicated in the legend. Additionally, the grey cross marks represent values obtained from non-self-gravitating simulations with corresponding wind powers.

From the right panel of Fig. 16, it is evident that the winds with power in the range $10^{42}-10^{45}$ erg s$^{-1}$ cause mass-outflow rates of $\sim 10^{-4}-10^{-2}$ M$_{\odot}$ yr$^{-1}$ from a single cloud through an area of $A = 100 \times 100$ pc. If we consider multiple clouds are distributed around the AGN globally, the global mass-outflow rates ($\dot{M}_{OF,glob}$) at a distance of $R_{OF}$ from the AGN can be calculated as

$$ \dot{M}_{OF,glob} \approx f_V \dot{M}_{OF,avg} \times \frac{4\pi R_{OF}^2}{A} \left(\frac{\Omega}{4\pi}\right), \quad (30) $$

where $f_V$ is the volume filling factor of the dense cloud, typically of the order of 0.01–0.1 (de Avillez & Breitschwerdt 2004; Mukherjee et al. 2018) and $\Omega$ is the solid angle covered by the outflowing region. For a spherical outflowing region with a maximum outflow rate $\Omega = 4\pi$, and $f_V = 1$. Thus, equation (30) becomes $\dot{M}_{OF,glob} \approx \dot{M}_{OF,avg} \times 4\pi R_{OF}^2/A$. Therefore, assuming a simple spherical arrangement of the clouds, the global mass outflow rates at a distance of $R_{OF} \sim 1.1$ kpc from the AGN are of the order of 0.1–10 M$_{\odot}$ yr$^{-1}$ for the considered wind powers in this study. These values are similar to what has been found in observations (see Fiore et al. 2017, and references therein).

A noteworthy feature in Fig. 16 is the tight correlation between the mass-outflow rate and the wind power. As the power of the wind increases, there is a corresponding increase in the mass-outflow rate, a trend that is in agreement with our intuitive expectations. Fitting a simple PL ($\dot{M}_{OF,avg} \propto P^\kappa$) to the results from the self-gravitating simulations (indicated by the red dashed line) yields a PL exponent of $\kappa \sim 0.52$. However, from the left panel, it is important to note the outflow rates for low-power winds (GC42 and GC43) have not reached a steady state by the end of the simulation and are expected to rise further before becoming constant. Therefore, the calculated average mass outflow rates for the lower power winds presented in the right panel serve as lower bounds. Thus, the value of $\kappa \sim 0.52$
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Mass-weighted 2D histogram of velocity versus density at 0.85$t_{cc}$ for self-gravity simulations with different wind power. The physical time for each simulation is indicated in the legends. The contours in each panel represent the initial distribution. The vertical dashed line indicates the gas density of 10 cm$^{-3}$, while a horizontal line marks $v_x = 100$ km s$^{-1}$.</figcaption></figure>

### 3.3.2 Multiphase structure

To analyse the gas phases of the outflowing material, we have constructed phase diagrams representing gas temperature versus number density in Fig. 15, where the colourbar in the top and bottom rows show the total mass and the mass-weighted mean velocity of the cloud material, respectively. From the top row, it is evident that most of the cloud material is in the cold phase ($T \lesssim 10^3$ K) for the winds with lower powers, whereas there exists a notable amount of moderately dense gas ($10 \lesssim 100$ cm$^{-3}$) in warm ($10^3 \lesssim T \lesssim 10^4$ K) and hot ($10^4 \lesssim T \lesssim 10^7$ K) phases in the higher power cases. Cloud materials in these phases experience acceleration up to 400 km s$^{-1}$. Interestingly, comparatively diffused gas with densities (0.1–10 cm$^{-3}$), which notably contributes to the ionized outflows, exists mostly in the hot phase, with velocities up to several 1000 km s$^{-1}$. On the other hand, gas in the dense and cold phase exhibits very low velocities, less than 100 km s$^{-1}$, even in the highest power case. This suggests that accelerating dense, cold interstellar gas to high velocities (~1000 km s$^{-1}$) via AGN-driven winds is exceedingly challenging.

In the simulations featuring lower power winds (GC42 and GC43) shown in Fig. 14, the presence of high-velocity gas is notably limited, due to the low ram pressure of the wind, which primarily accelerates the cloud.

### 3.3.3 Outflow rate scaling with wind power

In order to quantify the outflow driven by the winds and its dependence on the wind power, we calculate the total mass-outflow rates ($\dot{M}_{OF}$) of the cloud material as a function of time for self-gravitating simulations with different wind power, including the simulation with the kinetic wind (GC45_kinetic), which is illustrated in the left panel of Fig. 16. The outflow rates are calculated across the y-z plane at $x = 75$ pc. In the global picture, this corresponds to the outflow rate through an area of $A = 100 \times 100$ pc$^2$ from a single cloud situated at a distance $R_{OF} \sim 1.075$ kpc from the AGN as per our configuration. Notably, the mass outflow initiates after the compression phase when the ablation and stripping of the cloud start to dominate the evolution and progressively approaches a near-steady state value over time. The amount of outflow increases with increasing wind power, as expected. Furthermore, during the steady-state phase, the outflow rate induced by a kinetic wind (dashed-dotted magenta) surpasses that of a thermal wind (solid magenta) with equivalent power due to the higher momentum transfer efficiency. It is important to note

that, the mass-outflow rates, calculated at $x = 75$ pc, are entirely contributed by the stripped and entrained materials from the cloud.

To establish the relationship between the mass outflow rate and wind power, we determine the average mass outflow rate ($\dot{M}_{OF,avg}$) in each simulation after $\dot{M}_{OF}$ reaches a near-steady state, as indicated by the dotted lines (red for thermal and magenta for kinetic wind) in the left panel, which is plotted against the wind power in the right panel of Fig. 16. The open squares correspond to $\dot{M}_{OF,avg}$ values for self-gravitating simulations with different powers, as indicated in the legend. Additionally, the grey cross marks represent values obtained from non-self-gravitating simulations with corresponding wind powers.

From the right panel of Fig. 16, it is evident that the winds with power in the range $10^{42}-10^{45}$ erg s$^{-1}$ cause mass-outflow rates of ~$10^{-4}-10^{-2}$ M$_{\odot}$ yr$^{-1}$ from a single cloud through an area of $A = 100 \times 100$ pc. If we consider multiple clouds are distributed around the AGN globally, the global mass-outflow rates ($\dot{M}_{OF,glob}$) at a distance of $R_{OF}$ from the AGN can be calculated as

$$ \dot{M}_{OF,glob} \approx f_V \dot{M}_{OF,avg} \times \frac{4\pi R_{OF}^2}{A} \left(\frac{\Omega}{4\pi}\right), \quad (30) $$

where $f_V$ is the volume filling factor of the dense cloud, typically of the order of 0.01–0.1 (de Avillez & Breitschwerdt 2004; Mukherjee et al. 2018) and $\Omega$ is the solid angle covered by the outflowing region. For a spherical outflowing region with a maximum outflow rate $\Omega = 4\pi$, and $f_V = 1$. Thus, equation (30) becomes $\dot{M}_{OF,glob} \approx \dot{M}_{OF,avg} \times 4\pi R_{OF}^2/A$. Therefore, assuming a simple spherical arrangement of the clouds, the global mass outflow rates at a distance of $R_{OF} \sim 1.1$ kpc from the AGN are of the order of 0.1–10 M$_{\odot}$ yr$^{-1}$ for the considered wind powers in this study. These values are similar to what has been found in observations (see Fiore et al. 2017, and references therein).

A noteworthy feature in Fig. 16 is the tight correlation between the mass-outflow rate and the wind power. As the power of the wind increases, there is a corresponding increase in the mass-outflow rate, a trend that is in agreement with our intuitive expectations. Fitting a simple PL ($\dot{M}_{OF,avg} \propto P^\kappa$) to the results from the self-gravitating simulations (indicated by the red dashed line) yields a PL exponent of $\kappa \sim 0.52$. However, from the left panel, it is important to note the outflow rates for low-power winds (GC42 and GC43) have not reached a steady state by the end of the simulation and are expected to rise further before becoming constant. Therefore, the calculated average mass outflow rates for the lower power winds presented in the right panel serve as lower bounds. Thus, the value of $\kappa \sim 0.52$
<figure><img src="image_16.png" /><figcaption>**Figure 15.** 2D histograms of temperature versus density of the cloud material for the self-gravitating simulation of different powers at 0.85$t_{cc}$, which corresponds to different physical times for different wind power as indicated in the legends. The colour bars in the top and bottom rows represent the mass and mass-weighted mean velocity of the cloud material along the direction of the wind, respectively. The vertical dashed line indicates the gas density of 10 cm$^{-3}$, while a horizontal line marks $T = 1000$ K.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Left panel: the time evolution of mass outflow rates ($\dot{M}_{OF}$) through $x = 75$ pc surface for self-gravitating simulations with different powers as indicated in the legend. Right panel: the average mass outflow rate ($\dot{M}_{OF,avg}$) as a function of wind power. These $\dot{M}_{OF,avg}$ values are determined after reaching a near-steady state, marked by vertical dotted lines in the left panel (red for thermal winds and magenta for kinetic winds). The grey cross marks denote values for non-self-gravitating simulations at corresponding wind powers. Additionally, the red-dashed line represents a PL fit to $\dot{M}_{OF,avg}$ obtained from the self-gravitating runs.</figcaption></figure>

represents the upper bound of the PL exponent estimated in this study.

We present a simple complementary picture, where the mass outflow, driven by the AGN, is due to the stripped-off material from the embedded clouds inside the flow, rather than the swept-up shell model (e.g. Faucher-Giguère & Quataert 2012; King & Pounds 2015). We assume that the mass-loss due to the ablation of the cloud by the hot wind is driven by the pressure gradient at the interface (see Hartquist et al. 1986). For a purely hydrostatic ablation, the ablation rate ($\dot{M}_{abl}$) of one cloud embedded in a hot wind of Mach number

$M_w$ is given as (Hartquist et al. 1986; Wagner et al. 2012),

$$ \dot{M}_{abl} = \alpha \min(1, M_w^{4/3}) (M_c c_{s,c})^{2/3} (\rho_w v_w)^{1/3}, \quad (31) $$

where $M_c$ and $c_{s,c}$ are the mass and internal isothermal sound speed of the cloud and $\alpha$ is a constant factor of the order of unity for a spherical cloud. Keeping the cloud configuration identical, the mass ablation rate depends on the wind velocity as

$$ \dot{M}_{abl} \propto v_w^{5/3}. \quad (32) $$
<figure><img src="image_16.png" /><figcaption>**Figure 15.** 2D histograms of temperature versus density of the cloud material for the self-gravitating simulation of different powers at 0.85$t_{cc}$, which corresponds to different physical times for different wind power as indicated in the legends. The colour bars in the top and bottom rows represent the mass and mass-weighted mean velocity of the cloud material along the direction of the wind, respectively. The vertical dashed line indicates the gas density of 10 cm$^{-3}$, while a horizontal line marks $T = 1000$ K.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Left panel: the time evolution of mass outflow rates ($\dot{M}_{OF}$) through $x = 75$ pc surface for self-gravitating simulations with different powers as indicated in the legend. Right panel: the average mass outflow rate ($\dot{M}_{OF,avg}$) as a function of wind power. These $\dot{M}_{OF,avg}$ values are determined after reaching a near-steady state, marked by vertical dotted lines in the left panel (red for thermal winds and magenta for kinetic winds). The grey cross marks denote values for non-self-gravitating simulations at corresponding wind powers. Additionally, the red-dashed line represents a PL fit to $\dot{M}_{OF,avg}$ obtained from the self-gravitating runs.</figcaption></figure>

represents the upper bound of the PL exponent estimated in this study.

We present a simple complementary picture, where the mass outflow, driven by the AGN, is due to the stripped-off material from the embedded clouds inside the flow, rather than the swept-up shell model (e.g. Faucher-Giguère & Quataert 2012; King & Pounds 2015). We assume that the mass-loss due to the ablation of the cloud by the hot wind is driven by the pressure gradient at the interface (see Hartquist et al. 1986). For a purely hydrostatic ablation, the ablation rate ($\dot{M}_{abl}$) of one cloud embedded in a hot wind of Mach number

$M_w$ is given as (Hartquist et al. 1986; Wagner et al. 2012),

$$ \dot{M}_{abl} = \alpha \min(1, M_w^{4/3}) (M_c c_{s,c})^{2/3} (\rho_w v_w)^{1/3}, \quad (31) $$

where $M_c$ and $c_{s,c}$ are the mass and internal isothermal sound speed of the cloud and $\alpha$ is a constant factor of the order of unity for a spherical cloud. Keeping the cloud configuration identical, the mass ablation rate depends on the wind velocity as

$$ \dot{M}_{abl} \propto v_w^{5/3}. \quad (32) $$
Thus, by combining equations (32) and (12), we can show that

$$ \dot{M}_{abl} \propto P_w^{5/9}, \qquad (33) $$

where $P_w$ is the power of the wind. Therefore, if we assume the ablated mass is the source of the majority of the mass-loading in the outflow, then the outflow rate depends on the wind power as $\dot{M}_{OF} \propto P_w^{5/9} \propto P_w^{0.55}$, which is close to the value we are estimating from the simulations.

It is worth noting that similar correlations between the outflow rates and the bolometric luminosity of the AGN have been consistently reported in numerous observational investigations, spanning both the ionized and cold molecular phases. For instance, studies focusing on ionized outflows have documented a PL exponent of about $\kappa \sim 1.29$ (Fiore et al. 2017), whereas the exponent for cold molecular outflows varies between $\kappa \sim 0.68-0.76$ in different instances (Cicone et al. 2014; Fiore et al. 2017; Fluetsch et al. 2019).

While our analysis does not differentiate between different phases of the outflow, the scaling relation we estimate for the overall mass-outflow rate, encompassing both the hot, warm, and cold phases, as a function of wind power, is fairly close to the scaling relations to the molecular outflow. However, it's crucial to emphasize that despite the observed strong correlation between the average mass outflow rate ($\dot{M}_{OF,avg}$) and wind power as found in this study, there exist several sources of scatter contributing to this correlation. For instance, a cloud impacted by a kinetic wind tends to show a higher mass outflow rate than in the case of a thermal wind with similar power, as illustrated by the magenta squares in the right panel (filled for kinetic and open for thermal). Additionally, factors such as the morphology (i.e. whether the cloud is porous or compact) of the cloud and whether self-gravity plays a dominant role in the cloud's evolution at the impact stage can introduce considerable scatter to the correlation. Moreover, on a global scale, the geometric distribution (e.g. spherical, disc-like, etc.) of clouds around the driving source, as represented by the $f_V$ factor in equation (30), can also cause significant scatter in the global mass outflow rate. Furthermore, comparing outflow rates for different systems with the same power but at different evolutionary stages can also introduce scatter as the temporal evolution of $\dot{M}_{OF}$ can vary significantly (left panel of Fig. 16). Therefore, acknowledging these factors while interpreting the correlation between the mass outflow rate and wind power is crucial. However, it is noteworthy that our simplified estimates do present qualitative similarities to the observed results, reinforcing the idea that an ablation-based model of mass-loss can likely explain the observed correlations.

## 3.4 Star formation rate

Outflows from AGNs are believed to be one of the major drivers behind the star formation activity in galaxies. The local input of energy and momentum from these outflows into the star-forming regions can induce turbulence, which may have a dual role. On the one hand, the induced turbulence increases the stability of the clouds against gravitational forces preventing global collapse, which can reduce the star formation inside such clouds. On the other hand, it can promote overdensities via shock compression (e.g. see Mac Low & Klessen 2004), which may result in an enhanced star formation. Therefore, the effect of these AGN-driven winds on the star formation activity is determined by a complex interplay between different physical processes acting on the cloud scale.

We estimate the SFR in the simulation using a semi-analytical model of turbulence-regulated star formation (see Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund

2011; Federrath & Klessen 2012; Burkhart & Mocz 2019), which take into account the local physical quantities such as the velocity dispersion, virial parameter, local sound speed, turbulence driving mode, and magnetic field (if present) to estimate the SFR. In this framework, the SFR is calculated by integrating the density PDF from a critical density. The SFR of a cloud of mass $M_c$ is calculated as

$$ SFR = SFR_{ff} \frac{M_c}{t_{ff}(\rho_0)}, \qquad (34) $$

where $t_{ff}(\rho_0)$ is the freefall time at the mean density ($\rho_0$) of the cloud. The quantity $SFR_{ff}$, defined as the SFR per freefall time, is given by (Federrath & Klessen 2012),

$$ SFR_{ff} = \frac{\epsilon}{\phi_t} \int_{s_{crit}}^{\infty} \frac{t_{ff}(\rho_0)}{t_{ff}(\rho)} \frac{\rho}{\rho_0} p(s) ds, \qquad (35) $$

where $s = \ln(\rho/\rho_0)$ is the logarithmic density contrast, and $p(s)$ is the density PDF expressed in terms of $s$. On the GMC scale, the parameter $\epsilon$ is interpreted as the fraction of the global mass of the whole cloud that eventually turns into stars, and is typically 1–2 per cent (see section 3.1.1 of Mandal et al. 2021, for a discussion). Therefore, we set $\epsilon = 0.015$ in the calculation of $SFR_{ff}$. The $\phi_t$ parameter is a numerical factor of the order of unity to account for uncertainties in the integral and is set to $\phi_t = 2.04$ as calibrated in Federrath & Klessen (2012).

The critical logarithmic density ($s_{crit}$) of collapse, used as the lower limit of the integral of equation (35), is estimated by comparing the Jeans length to the sonic length (where the turbulent velocity dispersion is of the order of the local sound speed) and is given by (Krumholz & McKee 2005; Federrath & Klessen 2012)

$$ s_{crit} = \ln \left[ \left( \frac{\pi^2}{5} \right) \phi_x \alpha_{vir} M^2 \right], \qquad (36) $$

where $\phi_x$ is a numerical factor of the order of unity to account for slight differences in the exact equality between the Jeans length and the sonic scale (Krumholz & McKee 2005; Federrath & Klessen 2012) and set to 0.19 in this study (calibrated in Federrath & Klessen 2012). For the detailed theoretical background of the turbulence-regulated star formation model, we refer interested readers to Krumholz & McKee (2005), Federrath & Klessen (2012), Mandal et al. (2021), and references therein.

In order to estimate the SFR in our simulations, we use equations (24) and (26) for calculating the Mach number ($M$) and virial parameter ($\alpha_{vir}$), which gives us the value of $s_{crit}$ via equation (36). We then construct the density PDF $p(s)$ from the simulation data, which is subsequently integrated according to equation (35) to obtain $SFR_{ff}$. Finally, we calculate the SFR from equation (34), where equation (23) is used to estimate the cloud mass.

### 3.4.1 Dependence on wind power

Fig. 17 shows the time evolution of $SFR_{ff}$ (left panel) and SFR (right panel) for different simulations, as denoted in the legend. During the initial stages, as the wind begins to interact with the cloud, it injects a substantial amount of kinetic and thermal energy into it, leading to an increase in the velocity dispersion and virial parameter, as depicted in Fig. 10. Thus, the turbulent motions of the gas impede the gravitational collapse of many regions, resulting in a decreased SFR compared to the simulation without a wind. This behaviour is consistent among all the wind simulations; however, the degree of initial suppression in $SFR_{ff}$ increases with rising wind power. In fact, for the most powerful wind (magenta), the star formation
Thus, by combining equations (32) and (12), we can show that

$$ \dot{M}_{abl} \propto P_w^{5/9}, \quad (33) $$

where $P_w$ is the power of the wind. Therefore, if we assume the ablated mass is the source of the majority of the mass-loading in the outflow, then the outflow rate depends on the wind power as $\dot{M}_{OF} \propto P_w^{5/9} \propto P_w^{0.55}$, which is close to the value we are estimating from the simulations.

It is worth noting that similar correlations between the outflow rates and the bolometric luminosity of the AGN have been consistently reported in numerous observational investigations, spanning both the ionized and cold molecular phases. For instance, studies focusing on ionized outflows have documented a PL exponent of about $\kappa \sim 1.29$ (Fiore et al. 2017), whereas the exponent for cold molecular outflows varies between $\kappa \sim 0.68-0.76$ in different instances (Cicone et al. 2014; Fiore et al. 2017; Fluetsch et al. 2019).

While our analysis does not differentiate between different phases of the outflow, the scaling relation we estimate for the overall mass-outflow rate, encompassing both the hot, warm, and cold phases, as a function of wind power, is fairly close to the scaling relations to the molecular outflow. However, it's crucial to emphasize that despite the observed strong correlation between the average mass outflow rate ($\dot{M}_{OF,avg}$) and wind power as found in this study, there exist several sources of scatter contributing to this correlation. For instance, a cloud impacted by a kinetic wind tends to show a higher mass outflow rate than in the case of a thermal wind with similar power, as illustrated by the magenta squares in the right panel (filled for kinetic and open for thermal). Additionally, factors such as the morphology (i.e. whether the cloud is porous or compact) of the cloud and whether self-gravity plays a dominant role in the cloud's evolution at the impact stage can introduce considerable scatter to the correlation. Moreover, on a global scale, the geometric distribution (e.g. spherical, disc-like, etc.) of clouds around the driving source, as represented by the $f_V$ factor in equation (30), can also cause significant scatter in the global mass outflow rate. Furthermore, comparing outflow rates for different systems with the same power but at different evolutionary stages can also introduce scatter as the temporal evolution of $\dot{M}_{OF}$ can vary significantly (left panel of Fig. 16). Therefore, acknowledging these factors while interpreting the correlation between the mass outflow rate and wind power is crucial. However, it is noteworthy that our simplified estimates do present qualitative similarities to the observed results, reinforcing the idea that an ablation-based model of mass-loss can likely explain the observed correlations.

## 3.4 Star formation rate

Outflows from AGNs are believed to be one of the major drivers behind the star formation activity in galaxies. The local input of energy and momentum from these outflows into the star-forming regions can induce turbulence, which may have a dual role. On the one hand, the induced turbulence increases the stability of the clouds against gravitational forces preventing global collapse, which can reduce the star formation inside such clouds. On the other hand, it can promote overdensities via shock compression (e.g. see Mac Low & Klessen 2004), which may result in an enhanced star formation. Therefore, the effect of these AGN-driven winds on the star formation activity is determined by a complex interplay between different physical processes acting on the cloud scale.

We estimate the SFR in the simulation using a semi-analytical model of turbulence-regulated star formation (see Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund

2011; Federrath & Klessen 2012; Burkhart & Mocz 2019), which take into account the local physical quantities such as the velocity dispersion, virial parameter, local sound speed, turbulence driving mode, and magnetic field (if present) to estimate the SFR. In this framework, the SFR is calculated by integrating the density PDF from a critical density. The SFR of a cloud of mass $M_c$ is calculated as

$$ \text{SFR} = \text{SFR}_{\text{ff}} \frac{M_c}{t_{\text{ff}}(\rho_0)}, \quad (34) $$

where $t_{ff}(\rho_0)$ is the freefall time at the mean density ($\rho_0$) of the cloud. The quantity SFR$_{ff}$, defined as the SFR per freefall time, is given by (Federrath & Klessen 2012),

$$ \text{SFR}_{\text{ff}} = \frac{\epsilon}{\phi_t} \int_{s_{\text{crit}}}^{\infty} \frac{t_{\text{ff}}(\rho_0) \rho}{t_{\text{ff}}(\rho) \rho_0} p(s) ds, \quad (35) $$

where $s = \ln(\rho/\rho_0)$ is the logarithmic density contrast, and $p(s)$ is the density PDF expressed in terms of $s$. On the GMC scale, the parameter $\epsilon$ is interpreted as the fraction of the global mass of the whole cloud that eventually turns into stars, and is typically 1–2 per cent (see section 3.1.1 of Mandal et al. 2021, for a discussion). Therefore, we set $\epsilon = 0.015$ in the calculation of SFR$_{ff}$. The $\phi_t$ parameter is a numerical factor of the order of unity to account for uncertainties in the integral and is set to $\phi_t = 2.04$ as calibrated in Federrath & Klessen (2012).

The critical logarithmic density ($s_{crit}$) of collapse, used as the lower limit of the integral of equation (35), is estimated by comparing the Jeans length to the sonic length (where the turbulent velocity dispersion is of the order of the local sound speed) and is given by (Krumholz & McKee 2005; Federrath & Klessen 2012)

$$ s_{\text{crit}} = \ln \left[ \left( \frac{\pi^2}{5} \right) \phi_x \alpha_{\text{vir}} M^2 \right], \quad (36) $$

where $\phi_x$ is a numerical factor of the order of unity to account for slight differences in the exact equality between the Jeans length and the sonic scale (Krumholz & McKee 2005; Federrath & Klessen 2012) and set to 0.19 in this study (calibrated in Federrath & Klessen 2012). For the detailed theoretical background of the turbulence-regulated star formation model, we refer interested readers to Krumholz & McKee (2005), Federrath & Klessen (2012), Mandal et al. (2021), and references therein.

In order to estimate the SFR in our simulations, we use equations (24) and (26) for calculating the Mach number ($M$) and virial parameter ($\alpha_{vir}$), which gives us the value of $s_{crit}$ via equation (36). We then construct the density PDF $p(s)$ from the simulation data, which is subsequently integrated according to equation (35) to obtain SFR$_{ff}$. Finally, we calculate the SFR from equation (34), where equation (23) is used to estimate the cloud mass.

### 3.4.1 Dependence on wind power

Fig. 17 shows the time evolution of SFR$_{ff}$ (left panel) and SFR (right panel) for different simulations, as denoted in the legend. During the initial stages, as the wind begins to interact with the cloud, it injects a substantial amount of kinetic and thermal energy into it, leading to an increase in the velocity dispersion and virial parameter, as depicted in Fig. 10. Thus, the turbulent motions of the gas impede the gravitational collapse of many regions, resulting in a decreased SFR compared to the simulation without a wind. This behaviour is consistent among all the wind simulations; however, the degree of initial suppression in SFR$_{ff}$ increases with rising wind power. In fact, for the most powerful wind (magenta), the star formation
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Time evolution of the SFR per freefall time (left panel) and the absolute SFR (right panel), for different simulations as indicated in the legend. For comparison, the results for simulations without a wind are also shown (red lines). The grey-out parts of the plots in self-gravitating simulations correspond to the regime where the Jeans length of the highest density cell in the simulations is not resolved by at least four computational cells.</figcaption></figure>

is completely quenched. This is because higher wind power leads to a greater transfer of energy into the cloud, resulting in more pronounced turbulent motion within the gas.

Subsequent to this phase, as the wind fragments the cloud, it becomes easier for the wind to pass through the intercloudlet channels due to the smaller column depth. Hence, the strength of the interaction between the wind and cloud reduces. Additionally, the dense cloudlets embedded within the hot wind are more resilient to propagation of the transmitted shocks into the cores, thus, weakening the effect of the wind within the dense cores. Furthermore, the cloudlets that have been shock-compressed achieve significantly higher densities, resulting in a deeper gravitational potential. Therefore, the gravitational energy becomes dominant over the kinetic support, leading to an increased SFR<sub>ff</sub> after the initial suppression phase. For simulations without self-gravity, except for the C.45 case, the rebound reaches a level comparable to that in simulations without the wind, as depicted by the dashed lines. For the C.45 case, the cloud is destroyed and entrained by the wind before reaching that value of SFR<sub>ff</sub>.

However, in the presence of self-gravity, the shock compression causes the cloudlets to attain a much higher density. Additionally, the cloudlets merge due to the attractive gravitational force, forming a massive and highly dense central core, which undergoes runaway collapse. As a consequence, the SFR<sub>ff</sub> experiences a rapid increase following the initial decline (solid lines), to the extent that it surpasses the SFR<sub>ff</sub> value in the simulation without wind (red solid line). It is important to note that the stronger the wind, the earlier and the more rapid this increase in SFR<sub>ff</sub> occurs, as the higher power winds quickly compress the cloud, effectively reducing the free-fall time. Therefore, in all the simulations with self-gravity, the wind triggers the collapse of the cloud after the initial suppression. While the morphological features in Fig. 4 show that a high-power wind significantly disrupts the cloud, we note that a substantial portion of cloud material still undergoes compression and experiences rapid collapse.

It is crucial to emphasize that our simulations do not include self-consistent modelling of star formation using sink particles, which we aim to explore in a forthcoming study. Therefore, the reported

SFR values in this context should be interpreted as a qualitative trend between various wind parameters.

### 3.4.2 Dependence on cloud properties

The results discussed in the earlier sections focused on the impact of winds with different powers on the same initial cloud structures and density. However, clouds with different densities and morphology can react differently to the wind. Therefore, to investigate this, in this section, we present the results for the wind-cloud interaction for a fixed value of wind power ($10^{43} \text{ erg s}^{-1}$) but varying cloud properties, viz. different cloud wavenumber ($\text{GC43\_k}\{1,3,6,10\}$) and mean density ($\text{GC43\_k3\_low}$).

The top left panel of Fig. 18 shows the evolution of the cold gas fraction with time. Notably, the evolution appears similar for high-density clouds. Initially, there is a decrease in the cold gas fraction attributed to energy injection from the wind, followed by a subsequent increase due to the efficient cooling of the compressed gas. However, for the low-density cloud ($\text{GC43\_k3\_low}$), the initial cold gas fraction is lower compared to the other cases. This is because, initially, the clouds in our simulations are in pressure equilibrium with the ambient medium, as described in Section 2.3. Due to the lower mean density, the equilibrium condition leads to a higher mean initial temperature (~2700 K) of the low-density cloud, resulting in a lower fraction of cold gas mass compared to the high-density clouds. However, as the evolution progresses, the gas rapidly cools, and the cold gas fraction increases. None the less, as the wind impacts the cloud, the fraction decreases, and the subsequent evolution closely resembles that of the high-density cases. The major difference lies in the fact that the cold gas fraction remains low and decreases rapidly due to higher degrees of ablation and mixing. However, until the end of the simulation (~3 Myr), a significant amount of cold gas survives.

The top right and bottom left panels of Fig. 18 illustrate the time evolution of velocity dispersion and virial parameter, respectively. All the values are normalized to the initial values in order to examine the effect of the wind on the evolution of the parameters with respect to the initial values. The qualitative evolution of both parameters
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Time evolution of the SFR per freefall time (left panel) and the absolute SFR (right panel), for different simulations as indicated in the legend. For comparison, the results for simulations without a wind are also shown (red lines). The grey-out parts of the plots in self-gravitating simulations correspond to the regime where the Jeans length of the highest density cell in the simulations is not resolved by at least four computational cells.</figcaption></figure>

is completely quenched. This is because higher wind power leads to a greater transfer of energy into the cloud, resulting in more pronounced turbulent motion within the gas.

Subsequent to this phase, as the wind fragments the cloud, it becomes easier for the wind to pass through the intercloudlet channels due to the smaller column depth. Hence, the strength of the interaction between the wind and cloud reduces. Additionally, the dense cloudlets embedded within the hot wind are more resilient to propagation of the transmitted shocks into the cores, thus, weakening the effect of the wind within the dense cores. Furthermore, the cloudlets that have been shock-compressed achieve significantly higher densities, resulting in a deeper gravitational potential. Therefore, the gravitational energy becomes dominant over the kinetic support, leading to an increased SFR<sub>ff</sub> after the initial suppression phase. For simulations without self-gravity, except for the C.45 case, the rebound reaches a level comparable to that in simulations without the wind, as depicted by the dashed lines. For the C.45 case, the cloud is destroyed and entrained by the wind before reaching that value of SFR<sub>ff</sub>.

However, in the presence of self-gravity, the shock compression causes the cloudlets to attain a much higher density. Additionally, the cloudlets merge due to the attractive gravitational force, forming a massive and highly dense central core, which undergoes runaway collapse. As a consequence, the SFR<sub>ff</sub> experiences a rapid increase following the initial decline (solid lines), to the extent that it surpasses the SFR<sub>ff</sub> value in the simulation without wind (red solid line). It is important to note that the stronger the wind, the earlier and the more rapid this increase in SFR<sub>ff</sub> occurs, as the higher power winds quickly compress the cloud, effectively reducing the free-fall time. Therefore, in all the simulations with self-gravity, the wind triggers the collapse of the cloud after the initial suppression. While the morphological features in Fig. 4 show that a high-power wind significantly disrupts the cloud, we note that a substantial portion of cloud material still undergoes compression and experiences rapid collapse.

It is crucial to emphasize that our simulations do not include self-consistent modelling of star formation using sink particles, which we aim to explore in a forthcoming study. Therefore, the reported

SFR values in this context should be interpreted as a qualitative trend between various wind parameters.

### 3.4.2 Dependence on cloud properties

The results discussed in the earlier sections focused on the impact of winds with different powers on the same initial cloud structures and density. However, clouds with different densities and morphology can react differently to the wind. Therefore, to investigate this, in this section, we present the results for the wind-cloud interaction for a fixed value of wind power ($10^{43} \text{ erg s}^{-1}$) but varying cloud properties, viz. different cloud wavenumber ($\text{GC43\_k}\{1,3,6,10\}$) and mean density ($\text{GC43\_k3\_low}$).

The top left panel of Fig. 18 shows the evolution of the cold gas fraction with time. Notably, the evolution appears similar for high-density clouds. Initially, there is a decrease in the cold gas fraction attributed to energy injection from the wind, followed by a subsequent increase due to the efficient cooling of the compressed gas. However, for the low-density cloud ($\text{GC43\_k3\_low}$), the initial cold gas fraction is lower compared to the other cases. This is because, initially, the clouds in our simulations are in pressure equilibrium with the ambient medium, as described in Section 2.3. Due to the lower mean density, the equilibrium condition leads to a higher mean initial temperature (~2700 K) of the low-density cloud, resulting in a lower fraction of cold gas mass compared to the high-density clouds. However, as the evolution progresses, the gas rapidly cools, and the cold gas fraction increases. None the less, as the wind impacts the cloud, the fraction decreases, and the subsequent evolution closely resembles that of the high-density cases. The major difference lies in the fact that the cold gas fraction remains low and decreases rapidly due to higher degrees of ablation and mixing. However, until the end of the simulation (~3 Myr), a significant amount of cold gas survives.

The top right and bottom left panels of Fig. 18 illustrate the time evolution of velocity dispersion and virial parameter, respectively. All the values are normalized to the initial values in order to examine the effect of the wind on the evolution of the parameters with respect to the initial values. The qualitative evolution of both parameters
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Time evolution of the cold gas fraction (top left panel), velocity dispersion (top right panel), virial parameter (bottom left panel) of the cloud material and the SFR per freefall time (bottom right panel) in the simulations with same wind power of $10^{43}$ erg s$^{-1}$ but different cloud wavenumber (GC43.k{1,3,6,10})$ and lower mean density (GC43.k3_low).</figcaption></figure>

is similar in all the simulations to what we observe in Fig. 10. Interestingly, the evolution for $k_{min} = 1, 6$, and 10 is quite similar apart from some initial deviation due to different initialization. As discussed in Section 3.1.4, while in $k_{min} = 1$ case the single large clump experiences rapid global collapse, for $k_{min} = 6, 10$ cases, the narrow intercloudlets channels prevent the wind from penetrating and dispersing the cloud significantly, eventually leading to global collapse. Despite different reasons, the ultimate fates of the clouds in these cases are similar. In contrast, for the $k_{min} = 3$ case, the intercloudlets channels are wide enough for the wind to penetrate into the cloud, which directly compresses the cloudlets locally but prevents the global collapse for a longer period.

For the cloud with a lower mean density (GC43.k3_low), the velocity dispersion increases by more than an order of magnitude as the wind induces stronger turbulent motion due to a lower density. The virial parameter also increases by two orders of magnitude in this case by the impact of wind. However, here we observe a second increase of the virial parameter at ~0.25$t_{ff}$. This is primarily a result of the increase in turbulent velocity together with the decrease in mean density due to the disintegration of the cloud, which contributes to the decrease in gravitational binding energy.

The bottom right panel of Fig. 18 shows the evolution of the normalized SFR per freefall time (SFR$_{ff}$) for all these simulations. Similar to Fig. 17, the initial increase of the virial parameter reduces the SFR$_{ff}$ for all simulations. However, the consequent rebound of SFR$_{ff}$ depends on the cloud properties. For the case with $k_{min} = 1$, one big clump (see top panel of Fig. 5) rapidly collapses via compression, leading to an early growth of SFR$_{ff}$. For $k_{min} = 6$ and 10 cases, the initial small-scale velocity field of the cloud, characterized by the same values of $k_{min}$, causes an enhanced amount of intercloud

motion. This motion persists for a longer period before decaying, keeping the SFR relatively low during this time. However, as the wind compresses the cloud globally, the majority of the cloudlets merge, effectively forming a large, compact clump, which undergoes rapid collapse globally. As a result, the SFR increases rapidly after the initial suppression.

In the case of the cloud with a lower mean density, we observe a complete cessation of star formation for approximately 0.1 $t_{ff}$ (~1.3 Myr) due to the initial rise of the virial parameter by two orders of magnitude. However, as some of the dense clump gets compressed by the wind, we observe some intermittent star formation activity after the initial quenching. None the less, the SFR values during these intervals remain more than one order of magnitude lower than the initial value until the end of the simulation.

# 4 DISCUSSION

## 4.1 Comparison to previous studies

There have been a few studies concerning the effect of highly pressurized winds possibly originating from AGN on the interstellar clouds (Mellema et al. 2002; Zubovas, Sabulis & Naujalis 2014; Dugan et al. 2017a). The results obtained in the present work add to the previous studies of AGN-driven wind-cloud interaction by including additional physics such as self-gravity and more realistic cloud morphologies (fractal structures) and carrying out simulations over a wide range of parameters for both winds and clouds. Some of the key findings from our study align with the results from previous works in the context of shock/wind-cloud interaction. Our results reveal the emergence of dense filaments due to ram pressure stripping
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Time evolution of the cold gas fraction (top left panel), velocity dispersion (top right panel), virial parameter (bottom left panel) of the cloud material and the SFR per freefall time (bottom right panel) in the simulations with same wind power of $10^{43}$ erg s$^{-1}$ but different cloud wavenumber (GC43.k{1,3,6,10})$ and lower mean density (GC43.k3_low).</figcaption></figure>

is similar in all the simulations to what we observe in Fig. 10. Interestingly, the evolution for $k_{min} = 1, 6$, and 10 is quite similar apart from some initial deviation due to different initialization. As discussed in Section 3.1.4, while in $k_{min} = 1$ case the single large clump experiences rapid global collapse, for $k_{min} = 6, 10$ cases, the narrow intercloudlets channels prevent the wind from penetrating and dispersing the cloud significantly, eventually leading to global collapse. Despite different reasons, the ultimate fates of the clouds in these cases are similar. In contrast, for the $k_{min} = 3$ case, the intercloudlets channels are wide enough for the wind to penetrate into the cloud, which directly compresses the cloudlets locally but prevents the global collapse for a longer period.

For the cloud with a lower mean density (GC43.k3_low), the velocity dispersion increases by more than an order of magnitude as the wind induces stronger turbulent motion due to a lower density. The virial parameter also increases by two orders of magnitude in this case by the impact of wind. However, here we observe a second increase of the virial parameter at $\sim 0.25 t_{ff}$. This is primarily a result of the increase in turbulent velocity together with the decrease in mean density due to the disintegration of the cloud, which contributes to the decrease in gravitational binding energy.

The bottom right panel of Fig. 18 shows the evolution of the normalized SFR per freefall time (SFR$_{ff}$) for all these simulations. Similar to Fig. 17, the initial increase of the virial parameter reduces the SFR$_{ff}$ for all simulations. However, the consequent rebound of SFR$_{ff}$ depends on the cloud properties. For the case with $k_{min} = 1$, one big clump (see top panel of Fig. 5) rapidly collapses via compression, leading to an early growth of SFR$_{ff}$. For $k_{min} = 6$ and 10 cases, the initial small-scale velocity field of the cloud, characterized by the same values of $k_{min}$, causes an enhanced amount of intercloud

motion. This motion persists for a longer period before decaying, keeping the SFR relatively low during this time. However, as the wind compresses the cloud globally, the majority of the cloudlets merge, effectively forming a large, compact clump, which undergoes rapid collapse globally. As a result, the SFR increases rapidly after the initial suppression.

In the case of the cloud with a lower mean density, we observe a complete cessation of star formation for approximately $0.1 t_{ff}$ ($\sim 1.3$ Myr) due to the initial rise of the virial parameter by two orders of magnitude. However, as some of the dense clump gets compressed by the wind, we observe some intermittent star formation activity after the initial quenching. None the less, the SFR values during these intervals remain more than one order of magnitude lower than the initial value until the end of the simulation.

# 4 DISCUSSION

## 4.1 Comparison to previous studies

There have been a few studies concerning the effect of highly pressurized winds possibly originating from AGN on the interstellar clouds (Mellema et al. 2002; Zubovas, Sabulis & Naujalis 2014; Dugan et al. 2017a). The results obtained in the present work add to the previous studies of AGN-driven wind-cloud interaction by including additional physics such as self-gravity and more realistic cloud morphologies (fractal structures) and carrying out simulations over a wide range of parameters for both winds and clouds. Some of the key findings from our study align with the results from previous works in the context of shock/wind-cloud interaction. Our results reveal the emergence of dense filaments due to ram pressure stripping
and KH instabilities acting upon the interface between the cloud and the wind, similar to the results from previous studies (Mellema et al. 2002; Orlando et al. 2005; Cooper et al. 2009; Pittard et al. 2010; Dugan et al. 2017a; Banda-Barragán et al. 2019, 2020; Li et al. 2020; Banda-Barragán et al. 2021).

However, it is worth noting that in our simulations, these filaments display a more clumpy structure as opposed to the extended, continuous filaments observed in earlier studies. This difference arises from our incorporation of a fractal medium characterized by the presence of dense cores and a diffuse intercore medium. This configuration facilitates the infiltration of the wind into the cloud by clearing out low-density channels and enveloping the cores, effectively breaking the cloud into numerous cloudlets. Moreover, as the KH instability grows faster at smaller wavelengths, the fragmentation of the cloud is much higher compared to a uniform or smooth cloud (Cooper et al. 2009). Therefore, these fragmented cloudlets undergo acceleration to shape the clumpy filaments observed in our simulations. In this way, the morphology of the filaments in our study resembles more the findings by Cooper et al. (2009) and Banda-Barragán et al. (2020, 2021), where a similar fractal cloud configuration is considered.

Another interesting phenomenon to be compared with previous studies is the lifetime of the dense clumps that are generated in the process of fragmentation. Cooper et al. (2009) concluded that radiative cooling plays a significant role in prolonging the survival of dense clumps entrained within a hot wind. This cooling process leads to the creation of a dense protective layer around the clump, enhancing the cloud's lifetime. In such scenarios, these clouds can survive for more than 1 Myr, which aligns with our results. In all simulations performed in this study, the majority of the fragmented clumps survive until the end of the simulations (1.2–5.5 Myr, depending on the wind power). Additionally, as demonstrated previously in Sections 3.1.2 and 3.2.1, self-gravity further increases the density of the shock-compressed cloudlets, and since the growth rate of instabilities is inversely proportional to the density contrast (see equation 21), the growth of instabilities is diminished, therefore, further prolonging the lifetime of the cloudlets. Additionally, the increased density contrast due to self-gravity reduces the velocity of the transmitted shock ($v_{ts} \sim v_w/\chi^{1/2}$) into the cloudlets, consequently reducing the overall heating of the cores.

The effect of self-gravity has been considered in a few previous studies (Zubovas et al. 2014; Falle, Vaidya & Hartquist 2017; Dugan et al. 2017a; Li et al. 2020; Girichidis et al. 2021; Kupilas et al. 2022) in the context of the cloud wind/shock interaction. All of these studies have demonstrated that when self-gravity is taken into account, the compression resulting from the wind contributes to an increased rate of fragmentation and eventual gravitational collapse (Zubovas et al. 2014), which agrees with our findings. Furthermore, Dugan et al. (2017a) have shown that under specific conditions, a sufficiently high ram pressure in the wind can lead to the complete destruction of the cloud. However, our study diverges from this result as we do not observe any instance of a cloud being entirely destroyed by the wind, even when subjected to the most powerful winds in our simulations. Moreover, despite the substantial disruption of the cloud by the wind in the simulation with an initial lower mean density (GC43.k3_low), pockets of dense material persist, fostering intermittent star formation activity (see bottom panel of Fig. 18). This discrepancy can be attributed to two major factors. First, in order to obtain that large value of ram pressure of the wind, Dugan et al. (2017a) consider quite higher values of the wind density. Therefore, the momentum imparted by the wind into the cloud is orders of magnitude higher than in our highest wind-power case. As a result, the cloud in their study undergoes rapid acceleration, and the combination of

Rayleigh-Taylor instabilities and ram pressure stripping leads to the cloud's disintegration. Secondly, we consider a fractal cloud in our simulations, which has many dense gravitationally bound cores at the initiation of the simulations. These cores are further compressed by the wind, rendering them self-shielded and resistant to ablation, which contributes to the cloud's survival, despite the wind's powerful effects.

## 4.2 Implication for AGN feedback

The impact of AGN feedback on star formation activity remains a complex and elusive process in current astrophysical research. While numerous studies support the idea of negative feedback, attributed to the turbulence and thermal energy enhancement induced by AGN winds or the jet-inflated cocoon, there are also proponents of positive feedback. In this scenario, the overpressurized wind associated with AGN activity has the effect of compressing and fragmenting star-forming clouds. This compression leads to an increased SFR by triggering the collapse of gas clouds into new stars.

While there is little consensus on negative feedback from the theoretical point of view, a few studies have shown that AGN activity can indeed suppress the SFR globally inside the host galaxy by overall induction of turbulence and thermal energy (Mandal et al. 2021; Mercedes-Feliz et al. 2023). Furthermore, Wagner et al. (2012) have argued that the gas within AGN-driven outflows can become unbound, effectively escaping the host galaxy's potential and thereby reducing the available star-forming fuel. This long-term process may contribute to a negative feedback mechanism.

Numerous theoretical studies have presented scenarios in which the triggering of star formation by AGN activity can be a viable mechanism (e.g. Wagner & Bicknell 2011; Gaibler et al. 2012; Nayakshin & Zubovas 2012; Dugan et al. 2017a, 2017b; Mukherjee et al. 2018; Mandal et al. 2021). For instance, Gaibler et al. (2012) and Mukherjee et al. (2018) have shown that the compression resulting from the high-pressure bubble inflated by an AGN jet can enhance the SFR in a disc galaxy by a factor of 2–3. Similarly, Zubovas et al. (2014) arrived at a similar conclusion, demonstrating that external pressurization within the ISM can confine and compress star-forming regions and accelerate the onset of gravitational collapse and subsequent star formation. Moreover, the study by Dugan et al. (2017a) revealed that the star formation efficiency depends on the ram pressure of the wind. They identified a critical threshold value of the ram pressure, below which the clouds experience rapid collapse, leading to an enhancement in SFR. However, above this threshold, ram pressure becomes strong enough that the clouds are ablated before gravitational forces can significantly influence their evolution.

The result we present in this study also leans toward a positive feedback scenario. While there exists a short period of suppressed/quenched star formation during the initial interaction between the wind and the cloud, the compression due to the overpressurized wind dominates in the long term. The radiative shocks compress many massive cloudlets to significantly higher densities, and the presence of self-gravity intensifies this process. As a result, the shock propagating into these cores decelerates swiftly as the gas densities progressively increase, and the cloudlets effectively become self-shielded from the wind, such that the cloudlets can survive for a long time in the hot wind. Without the support from internal turbulence inside these cores, which is crucial for their stability, they eventually collapse to form stars, which will stop once all the gas in the core is consumed.

Therefore, from our present results, positive feedback by the AGN is inevitable even for high-power winds. Indeed, some observational
and KH instabilities acting upon the interface between the cloud and the wind, similar to the results from previous studies (Mellema et al. 2002; Orlando et al. 2005; Cooper et al. 2009; Pittard et al. 2010; Dugan et al. 2017a; Banda-Barragán et al. 2019, 2020; Li et al. 2020; Banda-Barragán et al. 2021).

However, it is worth noting that in our simulations, these filaments display a more clumpy structure as opposed to the extended, continuous filaments observed in earlier studies. This difference arises from our incorporation of a fractal medium characterized by the presence of dense cores and a diffuse intercore medium. This configuration facilitates the infiltration of the wind into the cloud by clearing out low-density channels and enveloping the cores, effectively breaking the cloud into numerous cloudlets. Moreover, as the KH instability grows faster at smaller wavelengths, the fragmentation of the cloud is much higher compared to a uniform or smooth cloud (Cooper et al. 2009). Therefore, these fragmented cloudlets undergo acceleration to shape the clumpy filaments observed in our simulations. In this way, the morphology of the filaments in our study resembles more the findings by Cooper et al. (2009) and Banda-Barragán et al. (2020, 2021), where a similar fractal cloud configuration is considered.

Another interesting phenomenon to be compared with previous studies is the lifetime of the dense clumps that are generated in the process of fragmentation. Cooper et al. (2009) concluded that radiative cooling plays a significant role in prolonging the survival of dense clumps entrained within a hot wind. This cooling process leads to the creation of a dense protective layer around the clump, enhancing the cloud's lifetime. In such scenarios, these clouds can survive for more than 1 Myr, which aligns with our results. In all simulations performed in this study, the majority of the fragmented clumps survive until the end of the simulations (1.2–5.5 Myr, depending on the wind power). Additionally, as demonstrated previously in Sections 3.1.2 and 3.2.1, self-gravity further increases the density of the shock-compressed cloudlets, and since the growth rate of instabilities is inversely proportional to the density contrast (see equation 21), the growth of instabilities is diminished, therefore, further prolonging the lifetime of the cloudlets. Additionally, the increased density contrast due to self-gravity reduces the velocity of the transmitted shock ($v_{ts} \sim v_w/\chi^{1/2}$) into the cloudlets, consequently reducing the overall heating of the cores.

The effect of self-gravity has been considered in a few previous studies (Zubovas et al. 2014; Falle, Vaidya & Hartquist 2017; Dugan et al. 2017a; Li et al. 2020; Girichidis et al. 2021; Kupilas et al. 2022) in the context of the cloud wind/shock interaction. All of these studies have demonstrated that when self-gravity is taken into account, the compression resulting from the wind contributes to an increased rate of fragmentation and eventual gravitational collapse (Zubovas et al. 2014), which agrees with our findings. Furthermore, Dugan et al. (2017a) have shown that under specific conditions, a sufficiently high ram pressure in the wind can lead to the complete destruction of the cloud. However, our study diverges from this result as we do not observe any instance of a cloud being entirely destroyed by the wind, even when subjected to the most powerful winds in our simulations. Moreover, despite the substantial disruption of the cloud by the wind in the simulation with an initial lower mean density (GC43.k3_low), pockets of dense material persist, fostering intermittent star formation activity (see bottom panel of Fig. 18). This discrepancy can be attributed to two major factors. First, in order to obtain that large value of ram pressure of the wind, Dugan et al. (2017a) consider quite higher values of the wind density. Therefore, the momentum imparted by the wind into the cloud is orders of magnitude higher than in our highest wind-power case. As a result, the cloud in their study undergoes rapid acceleration, and the combination of

Rayleigh-Taylor instabilities and ram pressure stripping leads to the cloud's disintegration. Secondly, we consider a fractal cloud in our simulations, which has many dense gravitationally bound cores at the initiation of the simulations. These cores are further compressed by the wind, rendering them self-shielded and resistant to ablation, which contributes to the cloud's survival, despite the wind's powerful effects.

## 4.2 Implication for AGN feedback

The impact of AGN feedback on star formation activity remains a complex and elusive process in current astrophysical research. While numerous studies support the idea of negative feedback, attributed to the turbulence and thermal energy enhancement induced by AGN winds or the jet-inflated cocoon, there are also proponents of positive feedback. In this scenario, the overpressurized wind associated with AGN activity has the effect of compressing and fragmenting star-forming clouds. This compression leads to an increased SFR by triggering the collapse of gas clouds into new stars.

While there is little consensus on negative feedback from the theoretical point of view, a few studies have shown that AGN activity can indeed suppress the SFR globally inside the host galaxy by overall induction of turbulence and thermal energy (Mandal et al. 2021; Mercedes-Feliz et al. 2023). Furthermore, Wagner et al. (2012) have argued that the gas within AGN-driven outflows can become unbound, effectively escaping the host galaxy's potential and thereby reducing the available star-forming fuel. This long-term process may contribute to a negative feedback mechanism.

Numerous theoretical studies have presented scenarios in which the triggering of star formation by AGN activity can be a viable mechanism (e.g. Wagner & Bicknell 2011; Gaibler et al. 2012; Nayakshin & Zubovas 2012; Dugan et al. 2017a, 2017b; Mukherjee et al. 2018; Mandal et al. 2021). For instance, Gaibler et al. (2012) and Mukherjee et al. (2018) have shown that the compression resulting from the high-pressure bubble inflated by an AGN jet can enhance the SFR in a disc galaxy by a factor of 2–3. Similarly, Zubovas et al. (2014) arrived at a similar conclusion, demonstrating that external pressurization within the ISM can confine and compress star-forming regions and accelerate the onset of gravitational collapse and subsequent star formation. Moreover, the study by Dugan et al. (2017a) revealed that the star formation efficiency depends on the ram pressure of the wind. They identified a critical threshold value of the ram pressure, below which the clouds experience rapid collapse, leading to an enhancement in SFR. However, above this threshold, ram pressure becomes strong enough that the clouds are ablated before gravitational forces can significantly influence their evolution.

The result we present in this study also leans toward a positive feedback scenario. While there exists a short period of suppressed/quenched star formation during the initial interaction between the wind and the cloud, the compression due to the overpressurized wind dominates in the long term. The radiative shocks compress many massive cloudlets to significantly higher densities, and the presence of self-gravity intensifies this process. As a result, the shock propagating into these cores decelerates swiftly as the gas densities progressively increase, and the cloudlets effectively become self-shielded from the wind, such that the cloudlets can survive for a long time in the hot wind. Without the support from internal turbulence inside these cores, which is crucial for their stability, they eventually collapse to form stars, which will stop once all the gas in the core is consumed.

Therefore, from our present results, positive feedback by the AGN is inevitable even for high-power winds. Indeed, some observational
studies support this scenario. For instance, Maiolino et al. (2017) and Gallagher et al. (2019) detected a substantial amount of star formation (~15 M☉ yr⁻¹) inside a galactic outflow, with a significant young (~10 Myr) stellar population and higher stellar velocities (≤ 100 km s⁻¹). This indicates that the stars have been formed inside the outflows, which have been triggered by the compression induced by the out-flowing gas. A recent investigation by Duggal et al. (2023) identified a young (1–10 Myr) stellar population inside young compact steep spectrum radio galaxies. Notably, the dynamical age of the radio sources in these galaxies corresponds well to the stellar ages, suggesting that star formation was triggered by jet activity. There are various other observational studies where enhanced star formation by AGN activity (jets or winds) has been reported (e.g. Bicknell et al. 2000; Zinn et al. 2013; Salomé et al. 2015; Bernhard et al. 2016; Salomé et al. 2017), in agreement with our findings.

Moreover, as outlined in Section 2.4, the wind velocity and pressure experienced by the cloud at a distance of 1 kpc for different powers can be mapped to various distances from the AGN for a specific wind power. For example, the wind parameters corresponding to powers of 10⁴², 10⁴³, 10⁴⁴, and 10⁴⁵ erg s⁻¹ at 1kpc (refer to Table 1) can be replicated by considering a wind with a power of 10⁴³ erg s⁻¹ at distances of 3.2, 1, 0.3, and 0.1 kpc from the AGN. Therefore, our findings hold relevance for comprehending the impact of AGN-driven winds on a galactic-scale environment. Importantly, even though a lower power wind (e.g. 10⁴², 10⁴³ erg s⁻¹) may not appear to significantly affect the cloud at a distance of 1 kpc according to our study, the closer environment of the AGN on a scale of hundreds of parsecs will likely be strongly influenced by the wind, akin to the outcomes observed in higher power simulations (e.g. 10⁴⁴, 10⁴⁵ erg s⁻¹). Similarly, for a higher power wind, the outskirts of the galaxy will be mildly affected, mirroring the outcomes of simulations with lower power winds presented in this study. In any case, the presence of self-gravity is likely to trigger star formation inside the clouds on a wide length scale while impacted by AGN-driven winds, albeit at different time-scales, as observed in simulations with different powers (Fig. 17).

The reported suppression of star formation by AGN activity, as observed in several studies (e.g. Ogle et al. 2007, 2010; Nesvadba et al. 2010; Alatalo et al. 2014, 2015; Lanz et al. 2016), continues to await theoretical confirmation. Our study, which incorporates comprehensive modelling of the interaction between a star-forming cloud and an AGN-driven wind, does not reveal any significant long-term suppression of star formation activity over the duration of an AGN lifetime (typically a few Myr). Based on our results, it appears that the only plausible means of reducing the star formation would be to employ mechanisms capable of disrupting the dense cores from the inside. Stellar feedback in the form of winds, jets, radiation pressure, and photoionization appears to be a viable mechanism (Menon, Federrath & Krumholz 2023), at least on the scale of a few parsecs, which corresponds to the typical size of the dense cloudlets we have identified in our study. Indeed, previous studies on star-cluster formation in clouds of size ~1–10 pc have demonstrated that stellar feedback can slow down the SFR (Federrath 2015; Grudić et al. 2018). Therefore, the underlying concept is that following the initial burst of star formation induced by wind-driven compression, the feedback from young and massive stars acts to disrupt the cloud, which in turn makes wind material entrainment more favourable. Therefore, in tandem with stellar feedback, the AGN wind could potentially contribute to the destruction of dense cores, breaking them into even smaller structures that are no longer prone to gravitational collapse due to their reduced size. This presents an intriguing scenario that we plan to investigate in future studies.

Additionally, the wind parameters examined in this study represent fast and hot AGN-driven outflows. However, the impact of slow and relatively cold winds, which are typically found at greater distances from the AGN, is not considered here and should be investigated in future.

## 5 CONCLUSION

In this study, we have performed a series of three-dimensional hydrodynamical simulations of the interaction between AGN-driven winds and star-forming interstellar clouds, including radiative cooling and realistic cloud morphology such as fractal geometry. We have conducted two sets of simulations with and without self-gravity to examine the effect of self-gravity on the evolution of the clouds. We consider a large range of parameter space for investigating various aspects of the evolution process, including the power of the wind, the mean density of the cloud, the fractal density distribution of the cloud, and whether the wind is dominated by kinetic or thermal energy. In the following, we summarize the main results of this study as follows:

(i) *Interaction of the wind with the fractal cloud:* When the wind interacts with the fractal cloud consisting of dense cores separated by low-density channels, it rapidly erases the low-density areas, resulting in the formation of numerous dense cores. Subsequently, these dense cores undergo compression due to radiative shocks.

(ii) *Effect of self-gravity:* While the cloudlets get compressed by the wind irrespective of the presence of self-gravity with the same wind power, the cloudlets formed in the self-gravitating simulations attain much higher densities and become gravitationally bound, compared to the cloudlets in the simulations without self-gravity. In the absence of self-gravity, after attaining the maximum possible compression, the clouds start to disintegrate and expand as a result of the momentum transfer from the wind to the cloud material.

(iii) *Dependence on wind power:* The amount of cloud material that is retained and accreted by the gravitationally bound cloudlets depends on the strength of the wind. For lower power winds, the momentum transfer is reduced, and as a result, a significant portion of the cloud material does not gain enough acceleration to overcome the gravitational potential of the cloud. Consequently, this material falls back into the potential well, forming a central, massive clump that undergoes rapid gravitational collapse. Conversely, for the higher power winds, although some initially high-density clumps collapse under self-gravity, a significant number of relatively low-density clumps are accelerated and dispersed by the wind prior to self-gravity becoming a significant factor.

(iv) *Effect of cloud morphology:* The size of the cloudlets (determined by the minimum wave number $k_{min}$ used to generate the fractal density distribution) inside the cloud significantly affects the evolution when impacted by a wind with the same power. For the $k_{min} = 1$ case, the whole cloud roughly contains a single large dense clump and is therefore large enough to prevent the shock from penetrating into its cores before it undergoes global gravitational collapse due to external pressurization. In contrast, for $k_{min} = 6$ and 10, the cloud is characterized by numerous dense cores, separated by narrow channels. However, when the initial shock sweeps across the cloud, it effectively blocks these channels with swept-up material, preventing the wind material from infiltrating the cloud. As a result, in the absence of significant dispersal of the cloudlets, the cloudlets accumulate near the centre by the gravitational force, ultimately experiencing global collapse. Thus, there exists a narrow range of $k_{min}$, for which the interclumps channels are wide enough for the
studies support this scenario. For instance, Maiolino et al. (2017) and Gallagher et al. (2019) detected a substantial amount of star formation (~15 M☉ yr⁻¹) inside a galactic outflow, with a significant young (~10 Myr) stellar population and higher stellar velocities (≤ 100 km s⁻¹). This indicates that the stars have been formed inside the outflows, which have been triggered by the compression induced by the out-flowing gas. A recent investigation by Duggal et al. (2023) identified a young (1–10 Myr) stellar population inside young compact steep spectrum radio galaxies. Notably, the dynamical age of the radio sources in these galaxies corresponds well to the stellar ages, suggesting that star formation was triggered by jet activity. There are various other observational studies where enhanced star formation by AGN activity (jets or winds) has been reported (e.g. Bicknell et al. 2000; Zinn et al. 2013; Salomé et al. 2015; Bernhard et al. 2016; Salomé et al. 2017), in agreement with our findings.

Moreover, as outlined in Section 2.4, the wind velocity and pressure experienced by the cloud at a distance of 1 kpc for different powers can be mapped to various distances from the AGN for a specific wind power. For example, the wind parameters corresponding to powers of 10⁴², 10⁴³, 10⁴⁴, and 10⁴⁵ erg s⁻¹ at 1kpc (refer to Table 1) can be replicated by considering a wind with a power of 10⁴³ erg s⁻¹ at distances of 3.2, 1, 0.3, and 0.1 kpc from the AGN. Therefore, our findings hold relevance for comprehending the impact of AGN-driven winds on a galactic-scale environment. Importantly, even though a lower power wind (e.g. 10⁴², 10⁴³ erg s⁻¹) may not appear to significantly affect the cloud at a distance of 1 kpc according to our study, the closer environment of the AGN on a scale of hundreds of parsecs will likely be strongly influenced by the wind, akin to the outcomes observed in higher power simulations (e.g. 10⁴⁴, 10⁴⁵ erg s⁻¹). Similarly, for a higher power wind, the outskirts of the galaxy will be mildly affected, mirroring the outcomes of simulations with lower power winds presented in this study. In any case, the presence of self-gravity is likely to trigger star formation inside the clouds on a wide length scale while impacted by AGN-driven winds, albeit at different time-scales, as observed in simulations with different powers (Fig. 17).

The reported suppression of star formation by AGN activity, as observed in several studies (e.g. Ogle et al. 2007, 2010; Nesvadba et al. 2010; Alatalo et al. 2014, 2015; Lanz et al. 2016), continues to await theoretical confirmation. Our study, which incorporates comprehensive modelling of the interaction between a star-forming cloud and an AGN-driven wind, does not reveal any significant long-term suppression of star formation activity over the duration of an AGN lifetime (typically a few Myr). Based on our results, it appears that the only plausible means of reducing the star formation would be to employ mechanisms capable of disrupting the dense cores from the inside. Stellar feedback in the form of winds, jets, radiation pressure, and photoionization appears to be a viable mechanism (Menon, Federrath & Krumholz 2023), at least on the scale of a few parsecs, which corresponds to the typical size of the dense cloudlets we have identified in our study. Indeed, previous studies on star-cluster formation in clouds of size ~1–10 pc have demonstrated that stellar feedback can slow down the SFR (Federrath 2015; Grudić et al. 2018). Therefore, the underlying concept is that following the initial burst of star formation induced by wind-driven compression, the feedback from young and massive stars acts to disrupt the cloud, which in turn makes wind material entrainment more favourable. Therefore, in tandem with stellar feedback, the AGN wind could potentially contribute to the destruction of dense cores, breaking them into even smaller structures that are no longer prone to gravitational collapse due to their reduced size. This presents an intriguing scenario that we plan to investigate in future studies.

Additionally, the wind parameters examined in this study represent fast and hot AGN-driven outflows. However, the impact of slow and relatively cold winds, which are typically found at greater distances from the AGN, is not considered here and should be investigated in future.

## 5 CONCLUSION

In this study, we have performed a series of three-dimensional hydrodynamical simulations of the interaction between AGN-driven winds and star-forming interstellar clouds, including radiative cooling and realistic cloud morphology such as fractal geometry. We have conducted two sets of simulations with and without self-gravity to examine the effect of self-gravity on the evolution of the clouds. We consider a large range of parameter space for investigating various aspects of the evolution process, including the power of the wind, the mean density of the cloud, the fractal density distribution of the cloud, and whether the wind is dominated by kinetic or thermal energy. In the following, we summarize the main results of this study as follows:

(i) *Interaction of the wind with the fractal cloud:* When the wind interacts with the fractal cloud consisting of dense cores separated by low-density channels, it rapidly erases the low-density areas, resulting in the formation of numerous dense cores. Subsequently, these dense cores undergo compression due to radiative shocks.

(ii) *Effect of self-gravity:* While the cloudlets get compressed by the wind irrespective of the presence of self-gravity with the same wind power, the cloudlets formed in the self-gravitating simulations attain much higher densities and become gravitationally bound, compared to the cloudlets in the simulations without self-gravity. In the absence of self-gravity, after attaining the maximum possible compression, the clouds start to disintegrate and expand as a result of the momentum transfer from the wind to the cloud material.

(iii) *Dependence on wind power:* The amount of cloud material that is retained and accreted by the gravitationally bound cloudlets depends on the strength of the wind. For lower power winds, the momentum transfer is reduced, and as a result, a significant portion of the cloud material does not gain enough acceleration to overcome the gravitational potential of the cloud. Consequently, this material falls back into the potential well, forming a central, massive clump that undergoes rapid gravitational collapse. Conversely, for the higher power winds, although some initially high-density clumps collapse under self-gravity, a significant number of relatively low-density clumps are accelerated and dispersed by the wind prior to self-gravity becoming a significant factor.

(iv) *Effect of cloud morphology:* The size of the cloudlets (determined by the minimum wave number $k_{min}$ used to generate the fractal density distribution) inside the cloud significantly affects the evolution when impacted by a wind with the same power. For the $k_{min} = 1$ case, the whole cloud roughly contains a single large dense clump and is therefore large enough to prevent the shock from penetrating into its cores before it undergoes global gravitational collapse due to external pressurization. In contrast, for $k_{min} = 6$ and 10, the cloud is characterized by numerous dense cores, separated by narrow channels. However, when the initial shock sweeps across the cloud, it effectively blocks these channels with swept-up material, preventing the wind material from infiltrating the cloud. As a result, in the absence of significant dispersal of the cloudlets, the cloudlets accumulate near the centre by the gravitational force, ultimately experiencing global collapse. Thus, there exists a narrow range of $k_{min}$, for which the interclumps channels are wide enough for the
wind to penetrate into the cloud and provide stability against global collapse by transferring energy and momentum, which is the case for $k_{min} = 3$.

(v) *Relative effect of thermal versus kinetic wind*: The total energy partition of the wind in the thermal and kinetic components has a major effect on the evolution of the cloud. A thermal energy-dominated wind primarily affects the cloud material by increasing its internal energy, which in turn is converted into kinetic energy of the gas. However, in the presence of strong radiative cooling, the internal energy of the gas quickly dissipates, therefore the effective momentum transfer from the wind to the cloud is reduced for a thermal wind. A kinetic energy-dominated wind directly transfers momentum to the cloud material, resulting in a higher level of expansion and acceleration.

(vi) *Evolution of the velocity dispersion*: The velocity dispersion inside the cloud in all the wind simulations increases due to energy transfer from the wind to the cloud material. This effect is more pronounced with stronger winds. In cases without self-gravity, the velocity dispersion stabilizes at a roughly constant value as the cloud disperses. With self-gravity, the velocity dispersion starts to increase again when gravity becomes dominant as it pulls fragmented cloudlets towards the core, inducing additional motions between the cloudlets.

(vii) *Effect on the virial parameter*: In all the simulations, an initial increase in the virial parameter occurs as the wind imparts kinetic energy, initially surpassing the gravitational energy. As the cloud compresses and the gas density rises, the gravitational potential deepens, reducing $\alpha_{vir}$. Simulations without self-gravity reach a steady state with balanced kinetic and gravitational energy. On the other hand, with self-gravity, the cloud becomes gravitationally bound and undergoes collapse, and $\alpha_{vir}$ drops significantly, especially in cases with low-power winds, due to the dominant gravitational binding energy.

(viii) *Generation of multiphase outflow*: The ablation of the cloud material by the wind can give rise to multiphase outflows with velocities from a few 100 to several 1000 km s$^{-1}$ over a huge range of temperatures ($10^2-10^7$ K), consisting of cold, warm, and hot gas. The calculated mass-outflow rates correlate tightly with the wind power ($M_{OF} \propto P^\kappa$). We find a PL exponent of $\kappa \approx 0.52$.

(ix) *Impact on the SFR*: In the presence of self-gravity, which is very important within the environment we are interested in, our results favour a positive feedback scenario triggered by the AGN-driven winds, at least within the parameter space we consider in this study. Even though the wind can suppress or quench the star formation for about 1 Myr during the initial interaction, a substantial number of shock-compressed, dense cloudlets manage to shield themselves from the wind's influence and subsequently undergo rapid gravitational collapse. This process ultimately leads to an increased SFR.

## ACKNOWLEDGEMENTS

We thank the referee, Dr Tiago Costa, for his thorough and insightful comments, which have improved the quality of the paper significantly. AM would like to thank Moun Meenakshi and Kavita Kumari for useful discussions about the effect of the AGN radiation field on interstellar gas. AM further thanks Moun Meenakshi for the help with the radiative transfer calculation using CLOUDY. CF acknowledges funding provided by the Australian Research Council (Discovery Project DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). We gratefully acknowledge high-performance computing resources provided by

the Australian National Computational Infrastructure (grant n72 and ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme, the Pegasus<fnref n="5" /> high-performance computing facilities of IUCAA, and the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391). DM and NN acknowledge support from the IFCPAR/CEFIPRA (project no. 6504-2) for collaborative research. This work has received funding from the European High Performance Computing Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany, Greece, Italy, Norway, and Spain under grant agreement No 101093441. This work is supported by 'Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC)', project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2. Spoke 3, Astrophysics and Cosmos Observations.

## DATA AVAILABILITY

Data related to this work will be shared upon reasonable request to the corresponding author.

## REFERENCES

Alatalo K. et al., 2014, ApJ, 780, 186

Alatalo K. et al., 2015, ApJ, 798, 31

Appel S. M., Burkhart B., Semenov V. A., Federrath C., Rosen A. L., 2022, ApJ, 927, 75

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Band I. M., Trzhaskovskaia M. B., Verner D. A., Iakovlev D. G., 1990, A&A, 237, 267

Banda-Barragán W. E., Parkin E. R., Federrath C., Crocker R. M., Bicknell G. V., 2016, MNRAS, 455, 1309

Banda-Barragán W. E., Federrath C., Crocker R. M., Bicknell G. V., 2018, MNRAS, 473, 3454

Banda-Barragán W. E., Zertuche F. J., Federrath C., García Del Valle J., Brüggen M., Wagner A. Y., 2019, MNRAS, 486, 4526

Banda-Barragán W. E., Brüggen M., Federrath C., Wagner A. Y., Scannapieco E., Cottle J., 2020, MNRAS, 499, 2173

Banda-Barragán W. E., Brüggen M., Heesen V., Scannapieco E., Cottle J., Federrath C., Wagner A. Y., 2021, MNRAS, 506, 5658

Bayliss A., Turkel E., 1980, Commun. Pure Appl. Math., 33, 707

Begelman M. C., Cioffi D. F., 1989, ApJ, 345, L21

Bernhard E., Mullaney J. R., Daddi E., Ciesla L., Schreiber C., 2016, MNRAS, 460, 902

Bertoldi F., McKee C. F., 1992, ApJ, 395, 140

Bicknell G. V., Sutherland R. S., van Breugel W. J. M., Dopita M. A., Dey A., Miley G. K., 2000, ApJ, 540, 678

Bischetti M. et al., 2017, A&A, 598, A122

Borguet B. C. J., Arav N., Edmonds D., Chamberlain C., Benn C., 2013, ApJ, 762, 49

Brusa M. et al., 2015, MNRAS, 446, 2394

Burkhart B., Mocz P., 2019, ApJ, 879, 129

Carniani S. et al., 2015, A&A, 580, A102

Cazzoli S., Arribas S., Maiolino R., Colina L., 2016, A&A, 590, A125

Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability. International Series of Monographs on Physics, Oxford University Press, Amen House, London E.C.

Cicone C. et al., 2014, A&A, 562, A21

Colella P., Woodward P. R., 1984, J. Comput. Phys., 54, 174

Collins D. C., Padoan P., Norman M. L., Xu H., 2011, ApJ, 731, 59

<footnote n="5">5 http://hpc.iucaa.in</footnote>
wind to penetrate into the cloud and provide stability against global collapse by transferring energy and momentum, which is the case for $k_{min} = 3$.

(v) *Relative effect of thermal versus kinetic wind*: The total energy partition of the wind in the thermal and kinetic components has a major effect on the evolution of the cloud. A thermal energy-dominated wind primarily affects the cloud material by increasing its internal energy, which in turn is converted into kinetic energy of the gas. However, in the presence of strong radiative cooling, the internal energy of the gas quickly dissipates, therefore the effective momentum transfer from the wind to the cloud is reduced for a thermal wind. A kinetic energy-dominated wind directly transfers momentum to the cloud material, resulting in a higher level of expansion and acceleration.

(vi) *Evolution of the velocity dispersion*: The velocity dispersion inside the cloud in all the wind simulations increases due to energy transfer from the wind to the cloud material. This effect is more pronounced with stronger winds. In cases without self-gravity, the velocity dispersion stabilizes at a roughly constant value as the cloud disperses. With self-gravity, the velocity dispersion starts to increase again when gravity becomes dominant as it pulls fragmented cloudlets towards the core, inducing additional motions between the cloudlets.

(vii) *Effect on the virial parameter*: In all the simulations, an initial increase in the virial parameter occurs as the wind imparts kinetic energy, initially surpassing the gravitational energy. As the cloud compresses and the gas density rises, the gravitational potential deepens, reducing $\alpha_{vir}$. Simulations without self-gravity reach a steady state with balanced kinetic and gravitational energy. On the other hand, with self-gravity, the cloud becomes gravitationally bound and undergoes collapse, and $\alpha_{vir}$ drops significantly, especially in cases with low-power winds, due to the dominant gravitational binding energy.

(viii) *Generation of multiphase outflow*: The ablation of the cloud material by the wind can give rise to multiphase outflows with velocities from a few 100 to several 1000 km s$^{-1}$ over a huge range of temperatures ($10^2-10^7$ K), consisting of cold, warm, and hot gas. The calculated mass-outflow rates correlate tightly with the wind power ($M_{OF} \propto P^\kappa$). We find a PL exponent of $\kappa \approx 0.52$.

(ix) *Impact on the SFR*: In the presence of self-gravity, which is very important within the environment we are interested in, our results favour a positive feedback scenario triggered by the AGN-driven winds, at least within the parameter space we consider in this study. Even though the wind can suppress or quench the star formation for about 1 Myr during the initial interaction, a substantial number of shock-compressed, dense cloudlets manage to shield themselves from the wind's influence and subsequently undergo rapid gravitational collapse. This process ultimately leads to an increased SFR.

## ACKNOWLEDGEMENTS

We thank the referee, Dr Tiago Costa, for his thorough and insightful comments, which have improved the quality of the paper significantly. AM would like to thank Moun Meenakshi and Kavita Kumari for useful discussions about the effect of the AGN radiation field on interstellar gas. AM further thanks Moun Meenakshi for the help with the radiative transfer calculation using CLOUDY. CF acknowledges funding provided by the Australian Research Council (Discovery Project DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). We gratefully acknowledge high-performance computing resources provided by

the Australian National Computational Infrastructure (grant n72 and ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme, the Pegasus<fnref n="5" /> high-performance computing facilities of IUCAA, and the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391). DM and NN acknowledge support from the IFCPAR/CEFIPRA (project no. 6504-2) for collaborative research. This work has received funding from the European High Performance Computing Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany, Greece, Italy, Norway, and Spain under grant agreement No 101093441. This work is supported by 'Italian Research Center on High Performance Computing Big Data and Quantum Computing (ICSC)', project funded by European Union - NextGenerationEU - and National Recovery and Resilience Plan (NRRP) - Mission 4 Component 2. Spoke 3, Astrophysics and Cosmos Observations.

## DATA AVAILABILITY

Data related to this work will be shared upon reasonable request to the corresponding author.

## REFERENCES

Alatalo K. et al., 2014, ApJ, 780, 186

Alatalo K. et al., 2015, ApJ, 798, 31

Appel S. M., Burkhart B., Semenov V. A., Federrath C., Rosen A. L., 2022, ApJ, 927, 75

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Band I. M., Trzhaskovskaia M. B., Verner D. A., Iakovlev D. G., 1990, A&A, 237, 267

Banda-Barragán W. E., Parkin E. R., Federrath C., Crocker R. M., Bicknell G. V., 2016, MNRAS, 455, 1309

Banda-Barragán W. E., Federrath C., Crocker R. M., Bicknell G. V., 2018, MNRAS, 473, 3454

Banda-Barragán W. E., Zertuche F. J., Federrath C., García Del Valle J., Brüggen M., Wagner A. Y., 2019, MNRAS, 486, 4526

Banda-Barragán W. E., Brüggen M., Federrath C., Wagner A. Y., Scannapieco E., Cottle J., 2020, MNRAS, 499, 2173

Banda-Barragán W. E., Brüggen M., Heesen V., Scannapieco E., Cottle J., Federrath C., Wagner A. Y., 2021, MNRAS, 506, 5658

Bayliss A., Turkel E., 1980, Commun. Pure Appl. Math., 33, 707

Begelman M. C., Cioffi D. F., 1989, ApJ, 345, L21

Bernhard E., Mullaney J. R., Daddi E., Ciesla L., Schreiber C., 2016, MNRAS, 460, 902

Bertoldi F., McKee C. F., 1992, ApJ, 395, 140

Bicknell G. V., Sutherland R. S., van Breugel W. J. M., Dopita M. A., Dey A., Miley G. K., 2000, ApJ, 540, 678

Bischetti M. et al., 2017, A&A, 598, A122

Borguet B. C. J., Arav N., Edmonds D., Chamberlain C., Benn C., 2013, ApJ, 762, 49

Brusa M. et al., 2015, MNRAS, 446, 2394

Burkhart B., Mocz P., 2019, ApJ, 879, 129

Carniani S. et al., 2015, A&A, 580, A102

Cazzoli S., Arribas S., Maiolino R., Colina L., 2016, A&A, 590, A125

Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability. International Series of Monographs on Physics, Oxford University Press, Amen House, London E.C.

Cicone C. et al., 2014, A&A, 562, A21

Colella P., Woodward P. R., 1984, J. Comput. Phys., 54, 174

Collins D. C., Padoan P., Norman M. L., Xu H., 2011, ApJ, 731, 59

<footnote n="5">5 http://hpc.iucaa.in</footnote>
Collins D. C., Kritsuk A. G., Padoan P., Li H., Xu H., Ustyugov S. D., Norman M. L., 2012, ApJ, 750, 13

Cooper J. L., Bicknell G. V., Sutherland R. S., Bland-Hawthorn J., 2009, ApJ, 703, 330

Costa T., Sijacki D., Haehnelt M. G., 2014, MNRAS, 444, 2355

Costa T., Pakmor R., Springel V., 2020, MNRAS, 497, 5229

Cottle J. N., Scannapieco E., Brüggen M., 2018, ApJ, 864, 96

Cottle J., Scannapieco E., Brüggen M., Banda-Barragán W., Federrath C., 2020, ApJ, 892, 59

Cresci G., Maiolino R., 2018, Nat. Astron., 2, 179

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

de Avillez M. A., Breitschwerdt D., 2004, A&A, 425, 899

Dhawalikar S., Federrath C., Davidovits S., Teyssier R., Nagel S. R., Remington B. A., Collins D. C., 2022, MNRAS, 514, 1782

Dugan Z., Gaibler V., Bieri R., Silk J., Rahman M., 2017a, ApJ, 839, 103

Dugan Z., Gaibler V., Silk J., 2017b, ApJ, 844, 37

Duggal C. et al., 2023, ApJ, 965, 40

Engquist B., Majda A., 1977, Math. Comput., 31, 629

Engquist B., Majda A., 1979, Commun. Pure Appl. Math., 32, 313

Fabian A. C., 2012, ARA&A, 50, 455

Faesi C. M., Lada C. J., Forbrich J., 2018, ApJ, 857, 19

Falle S. A. E. G., Vaidya B., Hartquist T. W., 2017, MNRAS, 465, 260

Faucher-Giguère C.-A., Quataert E., 2012, MNRAS, 425, 605

Federrath C., 2015, MNRAS, 450, 4035

Federrath C., 2018, Phys. Today, 71, 38

Federrath C., Banerjee S., 2015, MNRAS, 448, 3297

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Federrath C., Klessen R. S., 2013, ApJ, 763, 51

Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688, L79

Federrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692, 364

Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010a, A&A, 512, A81

Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010b, ApJ, 713, 269

Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, ApJ, 731, 62

Federrath C. et al., 2016, ApJ, 832, 143

Ferland G. J. et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385

Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel H., Lamastra A., Fiore F., 2010, A&A, 518, L155

Fiore F. et al., 2017, A&A, 601, A143

Fluetsch A. et al., 2019, MNRAS, 483, 4586

Fragile P. C., Murray S. D., Anninos P., van Breugel W., 2004, ApJ, 604, 74

Fujita A., Martin C. L., Mac Low M.-M., New K. C. B., Weaver R., 2009, ApJ, 698, 693

Gaibler V., Khochfar S., Krause M., Silk J., 2012, MNRAS, 425, 438

Gallagher R., Maiolino R., Belfiore F., Drory N., Riffel R., Riffel R. A., 2019, MNRAS, 485, 3409

García-Burillo S. et al., 2014, A&A, 567, A125

Genzel R. et al., 2014, ApJ, 796, 7

Gerrard I. A. et al., 2023, MNRAS, 526, 982

Girdhar A. et al., 2022, MNRAS, 512, 1608

Girichidis P., Konstandin L., Whitworth A. P., Klessen R. S., 2014, ApJ, 781, 91

Girichidis P., Naab T., Walch S., Berlok T., 2021, MNRAS, 505, 1083

Glover S. C. O., Federrath C., Mac Low M. M., Klessen R. S., 2010, MNRAS, 404, 2

Goldsmith P. F., Langer W. D., 1978, ApJ, 222, 881

Gronke M., Oh S. P., 2018, MNRAS, 480, L111

Gronke M., Oh S. P., 2020, MNRAS, 492, 1970

Grudić M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018, MNRAS, 475, 3511

Guillard P. et al., 2012, ApJ, 747, 95

Gurski K. F., 2004, SIAM J. Sci. Comput., 25, 2165

Harrison C. M., Alexander D. M., Mullaney J. R., Swinbank A. M., 2014, MNRAS, 441, 3306

Hartquist T. W., Dyson J. E., Pettini M., Smith L. J., 1986, MNRAS, 221, 715

Hedstrom G., 1979, J. Comput. Phys., 30, 222

Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

Holt J., Tadhunter C. N., Morganti R., 2008, MNRAS, 387, 639

Hughes A. et al., 2010, MNRAS, 406, 2065

Hughes A. et al., 2013, ApJ, 779, 46

Inskip K. J., Villar-Martín M., Tadhunter C. N., Morganti R., Holt J., Dicken D., 2008, MNRAS, 386, 1797

Jaupart E., Chabrier G., 2020, ApJ, 903, L2

Kalfountzou E. et al., 2014, MNRAS, 442, 1181

Kennicutt R. C., Jr, 1998a, ARA&A, 36, 189

Kennicutt R. C., Jr, 1998b, ApJ, 498, 541

Khullar S., Federrath C., Krumholz M. R., Matzner C. D., 2021, MNRAS, 507, 4335

King A., 2003, ApJ, 596, L27

King A., Pounds K., 2015, ARA&A, 53, 115

King A. R., Zubovas K., Power C., 2011, MNRAS, 415, L6

Klein R. I., McKee C. F., Colella P., 1994, ApJ, 420, 213

Klessen R. S., 2000, ApJ, 535, 869

Korista K. T., Bautista M. A., Arav N., Moe M., Costantini E., Benn C., 2008, ApJ, 688, 108

Koyama H., Inutsuka S.-I., 2000, ApJ, 532, 980

Krause M., Alexander P., 2007, MNRAS, 376, 465

Kritsuk A. G., Norman M. L., Wagner R., 2011, ApJ, 727, L20

Kritsuk A. G., Ustyugov S. D., Norman M. L., 2017, New J. Phys., 19, 065003

Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250

Kupilas M. M., Pittard J. M., Wareing C. J., Falle S. A. E. G., 2022, MNRAS, 513, 3345

Lacy M., Croft S., Fragile C., Wood S., Nyland K., 2017, ApJ, 838, 146

Laney C. B., 1998, Computational Gasdynamics. Cambridge Univ. Press, Cambridge

Lanz L., Ogle P. M., Alatalo K., Appleton P. N., 2016, ApJ, 826, 29

LeVeque R. J., 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, Cambridge Univ. Press, Cambridge

Leftley J. H., Nesvadba N. P. H., Bicknell G., Janssen R. M. J., Mukherjee D., Petrov R., Shende M. B., Zovaro H. R. M., 2024, preprint (arXiv:2404.04341)

Lewis G. M., Austin P. H., 2002, An Iterative Method for Generating Scaling Log-Normal Simulations. in 11th Conference on Atmospheric Radiation, American Meteorological Association, Columbia University, New York, https://ams.confex.com/ams/11AR11ICP/techprogram/paper.42772.htm

Li S., 2005, J. Comput. Phys., 203, 344

Li Y., Klessen R. S., Mac Low M.-M., 2003, ApJ, 592, 975

Li Z., Hopkins P. F., Squire J., Hummels C., 2020, MNRAS, 492, 1841

McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117

McNamara B. R., Nulsen P. E. J., 2012, New J. Phys., 14, 055023

Mac Low M.-M., Klessen R. S., 2004, Rev. Mod. Phys., 76, 125

Maiolino R. et al., 2017, Nature, 544, 202

Majda A., Osher S., 1975, Commun. Pure Appl. Math., 28, 479

Mandal A., Federrath C., Körtgen B., 2020, MNRAS, 493, 3098

Mandal A., Mukherjee D., Federrath C., Nesvadba N. P. H., Bicknell G. V., Wagner A. Y., Meenakshi M., 2021, MNRAS, 508, 4738

Mandal A., Mukherjee D., Mignone A., 2023, ApJS, 268, 40

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, MNRAS, 351, 169

Meenakshi M., Mukherjee D., Wagner A. Y., Nesvadba N. P. H., Morganti R., Janssen R. M. J., Bicknell G. V., 2022a, MNRAS, 511, 1622

Meenakshi M. et al., 2022b, MNRAS, 516, 766

Mellema G., Kurk J. D., Röttgering H. J. A., 2002, A&A, 395, L13

Menon S. H., Federrath C., Klaassen P., Kuiper R., Reiter M., 2021, MNRAS, 500, 1721

Menon S. H., Federrath C., Krumholz M. R., 2023, MNRAS, 521, 5160

Mercedes-Feliz J. et al., 2023, MNRAS, 524, 3446

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228

Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, ApJS, 198, 7

Miville-Deschênes M.-A., Murray N., Lee E. J., 2017, ApJ, 834, 57
Collins D. C., Kritsuk A. G., Padoan P., Li H., Xu H., Ustyugov S. D., Norman M. L., 2012, *ApJ*, 750, 13

Cooper J. L., Bicknell G. V., Sutherland R. S., Bland-Hawthorn J., 2009, *ApJ*, 703, 330

Costa T., Sijacki D., Haehnelt M. G., 2014, *MNRAS*, 444, 2355

Costa T., Pakmor R., Springel V., 2020, *MNRAS*, 497, 5229

Cottle J. N., Scannapieco E., Brüggen M., 2018, *ApJ*, 864, 96

Cottle J., Scannapieco E., Brüggen M., Banda-Barragán W., Federrath C., 2020, *ApJ*, 892, 59

Cresci G., Maiolino R., 2018, *Nat. Astron.*, 2, 179

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, *MNRAS*, 486, 2827

de Avillez M. A., Breitschwerdt D., 2004, *A&A*, 425, 899

Dhawalikar S., Federrath C., Davidovits S., Teyssier R., Nagel S. R., Remington B. A., Collins D. C., 2022, *MNRAS*, 514, 1782

Dugan Z., Gaibler V., Bieri R., Silk J., Rahman M., 2017a, *ApJ*, 839, 103

Dugan Z., Gaibler V., Silk J., 2017b, *ApJ*, 844, 37

Duggal C. et al., 2023, *ApJ*, 965, 40

Engquist B., Majda A., 1977, *Math. Comput.*, 31, 629

Engquist B., Majda A., 1979, *Commun. Pure Appl. Math.*, 32, 313

Fabian A. C., 2012, *ARA&A*, 50, 455

Faesi C. M., Lada C. J., Forbrich J., 2018, *ApJ*, 857, 19

Falle S. A. E. G., Vaidya B., Hartquist T. W., 2017, *MNRAS*, 465, 260

Faucher-Giguère C.-A., Quataert E., 2012, *MNRAS*, 425, 605

Federrath C., 2015, *MNRAS*, 450, 4035

Federrath C., 2018, *Phys. Today*, 71, 38

Federrath C., Banerjee S., 2015, *MNRAS*, 448, 3297

Federrath C., Klessen R. S., 2012, *ApJ*, 761, 156

Federrath C., Klessen R. S., 2013, *ApJ*, 763, 51

Federrath C., Klessen R. S., Schmidt W., 2008, *ApJ*, 688, L79

Federrath C., Klessen R. S., Schmidt W., 2009, *ApJ*, 692, 364

Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2010a, *A&A*, 512, A81

Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010b, *ApJ*, 713, 269

Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011, *ApJ*, 731, 62

Federrath C. et al., 2016, *ApJ*, 832, 143

Ferland G. J. et al., 2017, *Rev. Mex. Astron. Astrofis.*, 53, 385

Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel H., Lamastra A., Fiore F., 2010, *A&A*, 518, L155

Fiore F. et al., 2017, *A&A*, 601, A143

Fluetsch A. et al., 2019, *MNRAS*, 483, 4586

Fragile P. C., Murray S. D., Anninos P., van Breugel W., 2004, *ApJ*, 604, 74

Fujita A., Martin C. L., Mac Low M.-M., New K. C. B., Weaver R., 2009, *ApJ*, 698, 693

Gaibler V., Khochfar S., Krause M., Silk J., 2012, *MNRAS*, 425, 438

Gallagher R., Maiolino R., Belfiore F., Drory N., Riffel R., Riffel R. A., 2019, *MNRAS*, 485, 3409

García-Burillo S. et al., 2014, *A&A*, 567, A125

Genzel R. et al., 2014, *ApJ*, 796, 7

Gerrard I. A. et al., 2023, *MNRAS*, 526, 982

Girdhar A. et al., 2022, *MNRAS*, 512, 1608

Girichidis P., Konstandin L., Whitworth A. P., Klessen R. S., 2014, *ApJ*, 781, 91

Girichidis P., Naab T., Walch S., Berlok T., 2021, *MNRAS*, 505, 1083

Glover S. C. O., Federrath C., Mac Low M. M., Klessen R. S., 2010, *MNRAS*, 404, 2

Goldsmith P. F., Langer W. D., 1978, *ApJ*, 222, 881

Gronke M., Oh S. P., 2018, *MNRAS*, 480, L111

Gronke M., Oh S. P., 2020, *MNRAS*, 492, 1970

Grudić M. Y., Hopkins P. F., Faucher-Giguère C.-A., Quataert E., Murray N., Kereš D., 2018, *MNRAS*, 475, 3511

Guillard P. et al., 2012, *ApJ*, 747, 95

Gurski K. F., 2004, *SIAM J. Sci. Comput.*, 25, 2165

Harrison C. M., Alexander D. M., Mullaney J. R., Swinbank A. M., 2014, *MNRAS*, 441, 3306

Hartquist T. W., Dyson J. E., Pettini M., Smith L. J., 1986, *MNRAS*, 221, 715

Hedstrom G., 1979, *J. Comput. Phys.*, 30, 222

Hennebelle P., Chabrier G., 2011, *ApJ*, 743, L29

Holt J., Tadhunter C. N., Morganti R., 2008, *MNRAS*, 387, 639

Hughes A. et al., 2010, *MNRAS*, 406, 2065

Hughes A. et al., 2013, *ApJ*, 779, 46

Inskip K. J., Villar-Martín M., Tadhunter C. N., Morganti R., Holt J., Dicken D., 2008, *MNRAS*, 386, 1797

Jaupart E., Chabrier G., 2020, *ApJ*, 903, L2

Kalfountzou E. et al., 2014, *MNRAS*, 442, 1181

Kennicutt R. C., Jr, 1998a, *ARA&A*, 36, 189

Kennicutt R. C., Jr, 1998b, *ApJ*, 498, 541

Khullar S., Federrath C., Krumholz M. R., Matzner C. D., 2021, *MNRAS*, 507, 4335

King A., 2003, *ApJ*, 596, L27

King A., Pounds K., 2015, *ARA&A*, 53, 115

King A. R., Zubovas K., Power C., 2011, *MNRAS*, 415, L6

Klein R. I., McKee C. F., Colella P., 1994, *ApJ*, 420, 213

Klessen R. S., 2000, *ApJ*, 535, 869

Korista K. T., Bautista M. A., Arav N., Moe M., Costantini E., Benn C., 2008, *ApJ*, 688, 108

Koyama H., Inutsuka S.-I., 2000, *ApJ*, 532, 980

Krause M., Alexander P., 2007, *MNRAS*, 376, 465

Kritsuk A. G., Norman M. L., Wagner R., 2011, *ApJ*, 727, L20

Kritsuk A. G., Ustyugov S. D., Norman M. L., 2017, *New J. Phys.*, 19, 065003

Krumholz M. R., McKee C. F., 2005, *ApJ*, 630, 250

Kupilas M. M., Pittard J. M., Wareing C. J., Falle S. A. E. G., 2022, *MNRAS*, 513, 3345

Lacy M., Croft S., Fragile C., Wood S., Nyland K., 2017, *ApJ*, 838, 146

Laney C. B., 1998, Computational Gasdynamics. Cambridge Univ. Press, Cambridge

Lanz L., Ogle P. M., Alatalo K., Appleton P. N., 2016, *ApJ*, 826, 29

LeVeque R. J., 2002, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, Cambridge Univ. Press, Cambridge

Leftley J. H., Nesvadba N. P. H., Bicknell G., Janssen R. M. J., Mukherjee D., Petrov R., Shende M. B., Zovaro H. R. M., 2024, preprint (arXiv:2404.04341)

Lewis G. M., Austin P. H., 2002, An Iterative Method for Generating Scaling Log-Normal Simulations. in 11th Conference on Atmospheric Radiation, American Meteorological Association, Columbia University, New York, https://ams.confex.com/ams/11AR11ICP/techprogram/paper.42772.htm

Li S., 2005, *J. Comput. Phys.*, 203, 344

Li Y., Klessen R. S., Mac Low M.-M., 2003, *ApJ*, 592, 975

Li Z., Hopkins P. F., Squire J., Hummels C., 2020, *MNRAS*, 492, 1841

McNamara B. R., Nulsen P. E. J., 2007, *ARA&A*, 45, 117

McNamara B. R., Nulsen P. E. J., 2012, *New J. Phys.*, 14, 055023

Mac Low M.-M., Klessen R. S., 2004, *Rev. Mod. Phys.*, 76, 125

Maiolino R. et al., 2017, *Nature*, 544, 202

Majda A., Osher S., 1975, *Commun. Pure Appl. Math.*, 28, 479

Mandal A., Federrath C., Körtgen B., 2020, *MNRAS*, 493, 3098

Mandal A., Mukherjee D., Federrath C., Nesvadba N. P. H., Bicknell G. V., Wagner A. Y., Meenakshi M., 2021, *MNRAS*, 508, 4738

Mandal A., Mukherjee D., Mignone A., 2023, *ApJS*, 268, 40

Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., Salvati M., 2004, *MNRAS*, 351, 169

Meenakshi M., Mukherjee D., Wagner A. Y., Nesvadba N. P. H., Morganti R., Janssen R. M. J., Bicknell G. V., 2022a, *MNRAS*, 511, 1622

Meenakshi M. et al., 2022b, *MNRAS*, 516, 766

Mellema G., Kurk J. D., Röttgering H. J. A., 2002, *A&A*, 395, L13

Menon S. H., Federrath C., Klaassen P., Kuiper R., Reiter M., 2021, *MNRAS*, 500, 1721

Menon S. H., Federrath C., Krumholz M. R., 2023, *MNRAS*, 521, 5160

Mercedes-Feliz J. et al., 2023, *MNRAS*, 524, 3446

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, *ApJS*, 170, 228

Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, *ApJS*, 198, 7

Miville-Deschênes M.-A., Murray N., Lee E. J., 2017, *ApJ*, 834, 57
Moe M., Arav N., Bautista M. A., Korista K. T., 2009, ApJ, 706, 525

Morganti R., 2017, Front. Astron. Space Sci., 4, 42

Morganti R., Tadhunter C. N., Oosterloo T. A., 2005, A&A, 444, L9

Morganti R., Veilleux S., Oosterloo T., Teng S. H., Rupke D., 2016, A&A,
593, A30

Mukherjee D., Bicknell G. V., Sutherland R., Wagner A., 2016, MNRAS,
461, 967

Mukherjee D., Bicknell G. V., Wagner A. Y., Sutherland R. S., Silk J., 2018,
MNRAS, 479, 5544

Nakamura F., McKee C. F., Klein R. I., Fisher R. T., 2006, ApJS, 164, 477

Nayakshin S., 2014, MNRAS, 437, 2404

Nayakshin S., Zubovas K., 2012, MNRAS, 427, 372

Nesvadba N. P. H., Lehnert M. D., Eisenhauer F., Gilbert A., Tecza M., Abuter
R., 2006, ApJ, 650, 693

Nesvadba N. P. H., Lehnert M. D., De Breuck C., Gilbert A. M., van Breugel
W., 2008, A&A, 491, 407

Nesvadba N. P. H. et al., 2010, A&A, 521, A65

Nesvadba N. P. H., Boulanger F., Lehnert M. D., Guillard P., Salome P., 2011, A&A, 536, L5

Nesvadba N. P. H. et al., 2021, A&A, 654, A8

Ogle P., Antonucci R., Appleton P. N., Whysong D., 2007, ApJ, 668, 699

Ogle P., Boulanger F., Guillard P., Evans D. A., Antonucci R., Appleton P.
N., Nesvadba N., Leipski C., 2010, ApJ, 724, 1193

Orlando S., Peres G., Reale F., Bocchino F., Rosner R., Plewa T., Siegel A., 2005, A&A, 444, 505

Ostriker E. C., Stone J. M., Gammie C. F., 2001, ApJ, 546, 980

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Passot T., Vázquez-Semadeni E., 1998, Phys. Rev. E, 58, 4501

Perucho M., 2024, A&A, 684, 7

Pittard J. M., Falle S. A. E. G., Hartquist T. W., Dyson J. E., 2009, MNRAS,
394, 1351

Pittard J. M., Hartquist T. W., Falle S. A. E. G., 2010, MNRAS, 405, 821

Poludnenko A. Y., Frank A., Blackman E. G., 2002, ApJ, 576, 832

Price D. J., Federrath C., Brunt C. M., 2011, ApJ, 727, L21

Rahmati A., Pawlik A. H., Raičević M., Schaye J., 2013, MNRAS, 430,
2427

Ramos Almeida C. et al., 2022, A&A, 658, A155

Richings A. J., Faucher-Giguère C.-A., 2018, MNRAS, 474, 3673

Rupke D. S. N., Veilleux S., 2013, ApJ, 768, 75

Rupke D. S. N., Gültekin K., Veilleux S., 2017, ApJ, 850, 40

Salomé Q., Salomé P., Combes F., 2015, A&A, 574, A34

Salomé Q., Salomé P., Miville-Deschênes M. A., Combes F., Hamer S., 2017, A&A, 608, A98

Scannapieco E., Brüggen M., 2015, ApJ, 805, 158

Schaye J. et al., 2015, MNRAS, 446, 521

Schaye J. et al., 2023, MNRAS, 526, 4978

Schmidt M., 1959, ApJ, 129, 243

Schmidt M., 1963, ApJ, 137, 758

Schneider E. E., Robertson B. E., 2017, ApJ, 834, 144

Sharda P. et al., 2022, MNRAS, 509, 2180

Shin J., Woo J.-H., Chung A., Baek J., Cho K., Kang D., Bae H.-J., 2019, ApJ, 881, 147

Silk J., 2005, MNRAS, 364, 1337

Silk J., Rees M. J., 1998, A&A, 331, L1

Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776

Sutherland R. S., Bicknell G. V., 2007, ApJS, 173, 37

Sutherland R. S., Dopita M. A., 2017, ApJS, 229, 34

Sutherland R., Dopita M., Binette L., Groves B., 2018, Astrophysics Source
Code Library, record ascl:1807.005

Tanner R., Weaver K. A., 2022, AJ, 163, 134

Teşileanu O., Mignone A., Massaglia S., 2008, A&A, 488, 429

Thompson K. W., 1990, J. Comput. Phys., 89, 439

Truelove J. K., Klein R. I., McKee C. F., Holliman J. H. II, Howell L. H.,
Greenough J. A., 1997, ApJ, 489, L179

Vazquez-Semadeni E., 1994, ApJ, 423, 681

Vázquez-Semadeni E., Gómez G. C., Jappsen A. K., Ballesteros-Paredes J.,
González R. F., Klessen R. S., 2007, ApJ, 657, 870

Verner D. A., Yakovlev D. G., Band I. M., Trzhaskovskaya M. B., 1993,
Atom. Data Nucl. Data Tables, 55, 233

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys.,
2, 42

Wagner A. Y., Bicknell G. V., 2011, ApJ, 728, 29

Wagner A. Y., Bicknell G. V., Umemura M., 2012, ApJ, 757, 136

Wagner A. Y., Umemura M., Bicknell G. V., 2013, ApJ, 763, L18

Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377

Weinberger R. et al., 2017, MNRAS, 465, 3291

Westmoquette M. S., Clements D. L., Bendo G. J., Khan S. A., 2012, MNRAS,
424, 416

Yirak K., Frank A., Cunningham A. J., 2010, ApJ, 722, 412

Zinn P. C., Middelberg E., Norris R. P., Dettmar R. J., 2013, ApJ, 774, 66

Zubovas K., King A., 2012, ApJ, 745, L34

Zubovas K., King A. R., 2014, MNRAS, 439, 400

Zubovas K., Sabulis K., Naujalis R., 2014, MNRAS, 442, 2837

APPENDIX A: NON-REFLECTING BOUNDARY
CONDITION FOR SUBSONIC INFLOW

In this section, we present a simplified approximation of the
non-reflecting boundary conditions for subsonic inflow for one-
dimensional flow.

A1 Characteristic form of the conservation equations

The one-dimensional (1D) conservation equations in the primitive
variable are given as

$$
\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial x} + \rho \frac{\partial u}{\partial x} = 0, \qquad (A1)
$$

$$
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{1}{\rho} \frac{\partial p}{\partial x} = 0, \quad (A2)
$$

$$
\frac{\partial p}{\partial t} + \rho a^2 \frac{\partial u}{\partial x} + u \frac{\partial p}{\partial x} = 0, \quad (A3)
$$

where *a* is the sound speed and defined as $a = \sqrt{\gamma p / \rho}$. Equations (A1)–(A3) can be written in terms of the primitive variable vector Q as

$$
\frac{\partial Q}{\partial t} + A(Q) \frac{\partial Q}{\partial x} = 0, \qquad (A4)
$$

where $Q$ and $A(Q)$ are defined as

$$
Q = \begin{bmatrix} \rho \\ u \\ p \end{bmatrix}, \quad A(Q) = \begin{bmatrix} u & \rho & 0 \\ 0 & u & \frac{1}{\rho} \\ 0 & \rho a^2 & u \end{bmatrix}. \qquad (A5)
$$

Applying a similarity transformation to A, one obtains

$$
A = S \Lambda S^{-1}, \qquad (A6)
$$

where S = [K¹, K², K³] is the matrix consisting of the right
eigenvectors of A such that AK⁽¹⁾ = λᵢK⁽¹⁾, where λᵢ are the eigenvalues
corresponding to K⁽¹⁾. Λ is a diagonal matrix with its diagonal
elements being λᵢ. Replacing equation (A6), in equation (A4), we
obtain

$$
\frac{\partial Q}{\partial t} + S \Lambda S^{-1} \frac{\partial Q}{\partial x} = 0. \qquad (A7)
$$

Multiplying equation (A7) with S⁻¹ from the left, we obtain

$$
S^{-1} \frac{\partial Q}{\partial t} + \Lambda S^{-1} \frac{\partial Q}{\partial x} = 0. \qquad (A8)
$$

If we approximate A(Q) to be locally constant, so is S⁻¹. Therefore,
defining

$$
W = S^{-1} Q, \tag{A9}
$$
Moe M., Arav N., Bautista M. A., Korista K. T., 2009, ApJ, 706, 525

Morganti R., 2017, Front. Astron. Space Sci., 4, 42

Morganti R., Tadhunter C. N., Oosterloo T. A., 2005, A&A, 444, L9

Morganti R., Veilleux S., Oosterloo T., Teng S. H., Rupke D., 2016, A&A,
593, A30

Mukherjee D., Bicknell G. V., Sutherland R., Wagner A., 2016, MNRAS,
461, 967

Mukherjee D., Bicknell G. V., Wagner A. Y., Sutherland R. S., Silk J., 2018,
MNRAS, 479, 5544

Nakamura F., McKee C. F., Klein R. I., Fisher R. T., 2006, ApJS, 164, 477

Nayakshin S., 2014, MNRAS, 437, 2404

Nayakshin S., Zubovas K., 2012, MNRAS, 427, 372

Nesvadba N. P. H., Lehnert M. D., Eisenhauer F., Gilbert A., Tecza M., Abuter R., 2006, ApJ, 650, 693

Nesvadba N. P. H., Lehnert M. D., De Breuck C., Gilbert A. M., van Breugel W., 2008, A&A, 491, 407

Nesvadba N. P. H. et al., 2010, A&A, 521, A65

Nesvadba N. P. H., Boulanger F., Lehnert M. D., Guillard P., Salome P., 2011, A&A, 536, L5

Nesvadba N. P. H. et al., 2021, A&A, 654, A8

Ogle P., Antonucci R., Appleton P. N., Whysong D., 2007, ApJ, 668, 699

Ogle P., Boulanger F., Guillard P., Evans D. A., Antonucci R., Appleton P. N., Nesvadba N., Leipski C., 2010, ApJ, 724, 1193

Orlando S., Peres G., Reale F., Bocchino F., Rosner R., Plewa T., Siegel A., 2005, A&A, 444, 505

Ostriker E. C., Stone J. M., Gammie C. F., 2001, ApJ, 546, 980

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Passot T., Vázquez-Semadeni E., 1998, Phys. Rev. E, 58, 4501

Perucho M., 2024, A&A, 684, 7

Pittard J. M., Falle S. A. E. G., Hartquist T. W., Dyson J. E., 2009, MNRAS,
394, 1351

Pittard J. M., Hartquist T. W., Falle S. A. E. G., 2010, MNRAS, 405, 821

Poludnenko A. Y., Frank A., Blackman E. G., 2002, ApJ, 576, 832

Price D. J., Federrath C., Brunt C. M., 2011, ApJ, 727, L21

Rahmati A., Pawlik A. H., Raičević M., Schaye J., 2013, MNRAS, 430, 2427

Ramos Almeida C. et al., 2022, A&A, 658, A155

Richings A. J., Faucher-Giguère C.-A., 2018, MNRAS, 474, 3673

Rupke D. S. N., Veilleux S., 2013, ApJ, 768, 75

Rupke D. S. N., Gültekin K., Veilleux S., 2017, ApJ, 850, 40

Salomé Q., Salomé P., Combes F., 2015, A&A, 574, A34

Salomé Q., Salomé P., Miville-Deschênes M. A., Combes F., Hamer S., 2017, A&A, 608, A98

Scannapieco E., Brüggen M., 2015, ApJ, 805, 158

Schaye J. et al., 2015, MNRAS, 446, 521

Schaye J. et al., 2023, MNRAS, 526, 4978

Schmidt M., 1959, ApJ, 129, 243

Schmidt M., 1963, ApJ, 137, 758

Schneider E. E., Robertson B. E., 2017, ApJ, 834, 144

Sharda P. et al., 2022, MNRAS, 509, 2180

Shin J., Woo J.-H., Chung A., Baek J., Cho K., Kang D., Bae H.-J., 2019, ApJ, 881, 147

Silk J., 2005, MNRAS, 364, 1337

Silk J., Rees M. J., 1998, A&A, 331, L1

Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776

Sutherland R. S., Bicknell G. V., 2007, ApJS, 173, 37

Sutherland R. S., Dopita M. A., 2017, ApJS, 229, 34

Sutherland R., Dopita M., Binette L., Groves B., 2018, Astrophysics Source
Code Library, record ascl:1807.005

Tanner R., Weaver K. A., 2022, AJ, 163, 134

Teşileanu O., Mignone A., Massaglia S., 2008, A&A, 488, 429

Thompson K. W., 1990, J. Comput. Phys., 89, 439

Truelove J. K., Klein R. I., McKee C. F., Holliman J. H. II, Howell L. H., Greenough J. A., 1997, ApJ, 489, L179

Vazquez-Semadeni E., 1994, ApJ, 423, 681

Vázquez-Semadeni E., Gómez G. C., Jappsen A. K., Ballesteros-Paredes J., González R. F., Klessen R. S., 2007, ApJ, 657, 870

Verner D. A., Yakovlev D. G., Band I. M., Trzhaskovskaya M. B., 1993, Atom. Data Nucl. Data Tables, 55, 233

Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat. Rev. Phys., 2, 42

Wagner A. Y., Bicknell G. V., 2011, ApJ, 728, 29

Wagner A. Y., Bicknell G. V., Umemura M., 2012, ApJ, 757, 136

Wagner A. Y., Umemura M., Bicknell G. V., 2013, ApJ, 763, L18

Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377

Weinberger R. et al., 2017, MNRAS, 465, 3291

Westmoquette M. S., Clements D. L., Bendo G. J., Khan S. A., 2012, MNRAS, 424, 416

Yirak K., Frank A., Cunningham A. J., 2010, ApJ, 722, 412

Zinn P. C., Middelberg E., Norris R. P., Dettmar R. J., 2013, ApJ, 774, 66

Zubovas K., King A., 2012, ApJ, 745, L34

Zubovas K., King A. R., 2014, MNRAS, 439, 400

Zubovas K., Sabulis K., Naujalis R., 2014, MNRAS, 442, 2837

APPENDIX A: NON-REFLECTING BOUNDARY CONDITION FOR SUBSONIC INFLOW

In this section, we present a simplified approximation of the non-reflecting boundary conditions for subsonic inflow for one-dimensional flow.

A1 Characteristic form of the conservation equations

The one-dimensional (1D) conservation equations in the primitive variable are given as

$$
\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial x} + \rho \frac{\partial u}{\partial x} = 0, \qquad (A1)
$$

$$
\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{1}{\rho} \frac{\partial p}{\partial x} = 0, \qquad (A2)
$$

$$
\frac{\partial p}{\partial t} + \rho a^2 \frac{\partial u}{\partial x} + u \frac{\partial p}{\partial x} = 0, \quad (A3)
$$

where *a* is the sound speed and defined as $a = \sqrt{\gamma p / \rho}$. Equations (A1)–(A3) can be written in terms of the primitive variable vector Q as

$$
\frac{\partial Q}{\partial t} + A(Q) \frac{\partial Q}{\partial x} = 0, \qquad (A4)
$$

where $Q$ and $A(Q)$ are defined as

$$
Q = \begin{bmatrix} \rho \\ u \\ p \end{bmatrix}, \quad A(Q) = \begin{bmatrix} u & \rho & 0 \\ 0 & u & \frac{1}{\rho} \\ 0 & \rho a^2 & u \end{bmatrix}. \qquad (A5)
$$

Applying a similarity transformation to A, one obtains

$$
A = S \Lambda S^{-1}, \qquad (A6)
$$

where S = [K¹, K², K³] is the matrix consisting of the right
eigenvectors of A such that AK⁽¹⁾ = λᵢK⁽¹⁾, where λᵢ are the eigenvalues
corresponding to K⁽¹⁾. Λ is a diagonal matrix with its diagonal
elements being λᵢ. Replacing equation (A6), in equation (A4), we
obtain

$$
\frac{\partial Q}{\partial t} + S \Lambda S^{-1} \frac{\partial Q}{\partial x} = 0. \qquad (A7)
$$

Multiplying equation (A7) with S⁻¹ from the left, we obtain

$$
S^{-1} \frac{\partial Q}{\partial t} + \Lambda S^{-1} \frac{\partial Q}{\partial x} = 0. \qquad (A8)
$$

If we approximate $A(Q)$ to be locally constant, so is $S^{-1}$. Therefore,
defining

$$
W = S^{-1} Q, \tag{A9}
$$
equation (A8) reduces to

$$ \frac{\partial W}{\partial t} + \Lambda \frac{\partial W}{\partial x} = 0, \qquad (A10) $$

which is the characteristic form of the conservation equations, and $W = [w^1, w^2, w^3]^T$ is the characteristic variable vector. In component form, equation (A10) can be written as

$$ \frac{\partial w^i}{\partial t} + \lambda_i \frac{\partial w^i}{\partial x} = 0, \qquad (A11) $$

which is a set of wave equations with characteristic velocities $\lambda_i$.

In order to compute the characteristic variables, first we find the eigenvalues ($\lambda_i$) and right eigenvectors ($K^i$) of A, which are given by

$$ \lambda_1 = u - a, \qquad (A12) $$

$$ \lambda_2 = u, \qquad (A13) $$

$$ \lambda_3 = u + a, \qquad (A14) $$

and the corresponding eigenvectors are given as

$$ K^1 = \alpha_1 \begin{bmatrix} 1 \\ -\frac{a}{\rho} \\ \frac{a}{a^2} \end{bmatrix}, \quad K^2 = \alpha_2 \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \quad K^3 = \alpha_3 \begin{bmatrix} 1 \\ \frac{a}{\rho} \\ \frac{a^2}{a^2} \end{bmatrix}, \qquad (A15) $$

where $\alpha_i$ are the scale factors. Assuming $\alpha_1 = \rho/2a$, $\alpha_2 = 1$, and $\alpha_3 = \rho/2a$, the S matrix and its inverse is calculated as

$$ S = \begin{bmatrix} \frac{\rho}{2a} & 1 & \frac{\rho}{2a} \\ -\frac{1}{2} & 0 & \frac{1}{2} \\ \frac{\rho}{2} & \frac{\rho}{a} & 0 \end{bmatrix}, \quad S^{-1} = \begin{bmatrix} 0 & -1 & \frac{1}{\rho a} \\ 1 & 0 & -\frac{1}{a^2} \\ 0 & 1 & \frac{1}{\rho a} \end{bmatrix}. \qquad (A16) $$

Therefore, from equations (A5), (A9), and (A16), the form of characteristic variables ($W = [w^1, w^2, w^3]^T$) for the 1D Euler equations are given as

$$ w^1 = \frac{p}{\rho a} - u, \quad \text{corresponding to } \lambda_1 = u - c, \qquad (A17) $$

$$ w^2 = \rho - \frac{p}{a^2}, \quad \text{corresponding to } \lambda_2 = u, \qquad (A18) $$

$$ w^3 = \frac{p}{\rho a} + u, \quad \text{corresponding to } \lambda_3 = u + c. \qquad (A19) $$

## A2 Boundary conditions

The primary idea behind the non-reflecting boundary condition is that any outgoing wave ($w^{\text{out}}$) that is leaving the computational domain through the boundary should leave the domain without getting reflected back into the interior region, which means the wave amplitude associated with the outgoing wave is constant at the boundary, i.e.

$$ w^{\text{out}} = \text{const}, \quad \text{or} \quad \Delta w^{\text{out}} = 0, \text{ at the boundary,} \qquad (A20) $$

and the boundary condition associated with the outgoing characteristic variables should be extrapolated from the interior solution. The other waves that are incoming into the computational domain require physical boundary conditions.

For a subsonic inflow boundary (let us assume x-beg), one wave ($w^1$) is leaving the domain ($\lambda_1 = u - a < 0$), and the other two ($w^2, w^3$) are incoming, which means only two out of three primitive variables (ideally ($\rho, u$) or ($\rho, p$) for well-poised conditions, e.g. see section 19.3 of Laney 1998) can be specified physically at the

boundary. The remaining variable has to be set numerically from the interior solution in such a way that $\Delta w^1 = 0$ at the boundary. Let's say, we want to specify the value of $\rho_B$ and $u_B$ physically, then the pressure value ($p_B$) at the boundary is calculated by setting

$$ \Delta w^1 = 0, \qquad (A21) $$

$$ w_B^1 - w_I^1 = 0, \qquad (A22) $$

$$ p_B = \gamma \rho_B (w_I^1 + u_B)^2, \qquad (A23) $$

where $w_I^1$ is the value of $w^1$ at the first active zone ($i = 1$) of the computational domain, which is to be calculated from the interior solution. Therefore, for target density and velocity values of $\rho_t$ and $u_t$, the ghost cell values of the primitive variables are given by

$$ \rho_B = \rho_t, \qquad (A24) $$

$$ u_B = u_t, \qquad (A25) $$

$$ p_B = \gamma \rho_t (w_I^1 + u_t)^2. \qquad (A26) $$

Equations (A24)–(A26) ensure that $\Delta w^1 = 0$ at the boundary, and the wave-reflection is minimized.

## A3 Test and comparison

In order to demonstrate how the wind emerges in the interior solution with and without NRBC, we perform one simulation with NRBC implemented at the subsonic inflow boundary with the same setup as described in Section 2. The wind and cloud initializations are the same as in GC45.k3 in Table 1, with half the resolution of the fiducial runs. We compare the evolution of the wind parameters with a simulation that uses the same initial condition and resolution, but using the method of wind injection adopted in this paper (Section 2.6). In the simulation with NRBC, we specify the density and velocity values at the ghost zone, while the pressure of the wind is calculated using equation (A26).

In Fig. A1, we present the time evolution of the wind density (top left panel), velocity (top right panel), pressure (bottom left panel), and power (bottom right panel), at the first active cell in simulations with (red lines) and without (blue lines) NRBC. The lower section of each panel illustrates the fractional deviation of each parameter from its intended value. As observed, the evolution of the physically specified wind parameters – namely density and velocity – remains consistent throughout the simulations with NRBC. However, the pressure of the injected wind shows a significant deviation compared to the simulation without NRBC. Consequently, the deviation of wind power from its intended value is substantial (~90 per cent). On the other hand, although all wind parameters in the simulation without NRBC deviate from their intended values (blue lines), the cumulative effect results in more stable wind power (with a variation of ~30 per cent) compared to the scenario with NRBC. Moreover, we have performed additional simulations using the saturated values of the density and velocity for the GC45 case (see Table 2) with the NRBC. We confirmed that, although the quantitative results differ from the result presented in this study during the transient phase, they closely resemble the results (e.g. the diagnostics presented in Figs 9, 10, and 17) at the saturation phase ($t \gtrsim 0.4t_{cc}$). Therefore, the qualitative results presented in this study are robust.
equation (A8) reduces to

$$ \frac{\partial W}{\partial t} + \Lambda \frac{\partial W}{\partial x} = 0, \qquad (A10) $$

which is the characteristic form of the conservation equations, and $W = [w^1, w^2, w^3]^T$ is the characteristic variable vector. In component form, equation (A10) can be written as

$$ \frac{\partial w^i}{\partial t} + \lambda_i \frac{\partial w^i}{\partial x} = 0, \qquad (A11) $$

which is a set of wave equations with characteristic velocities $\lambda_i$.

In order to compute the characteristic variables, first we find the eigenvalues ($\lambda_i$) and right eigenvectors ($K^i$) of A, which are given by

$$ \lambda_1 = u - a, \qquad (A12) $$

$$ \lambda_2 = u, \qquad (A13) $$

$$ \lambda_3 = u + a, \qquad (A14) $$

and the corresponding eigenvectors are given as

$$ K^1 = \alpha_1 \begin{bmatrix} 1 \\ -\frac{a}{\rho} \\ \frac{a}{a^2} \end{bmatrix}, \quad K^2 = \alpha_2 \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \quad K^3 = \alpha_3 \begin{bmatrix} 1 \\ \frac{a}{\rho} \\ \frac{a^2}{a^2} \end{bmatrix}, \qquad (A15) $$

where $\alpha_i$ are the scale factors. Assuming $\alpha_1 = \rho/2a$, $\alpha_2 = 1$, and $\alpha_3 = \rho/2a$, the S matrix and its inverse is calculated as

$$ S = \begin{bmatrix} \frac{\rho}{2a} & 1 & \frac{\rho}{2a} \\ -\frac{1}{2} & 0 & \frac{1}{2} \\ \frac{\rho \bar{a}}{2} & 0 & \frac{\rho a}{2} \end{bmatrix}, \quad S^{-1} = \begin{bmatrix} 0 & -1 & \frac{1}{\rho a} \\ 1 & 0 & -\frac{1}{a^2} \\ 0 & 1 & \frac{1}{\rho a} \end{bmatrix}. \qquad (A16) $$

Therefore, from equations (A5), (A9), and (A16), the form of characteristic variables ($W = [w^1, w^2, w^3]^T$) for the 1D Euler equations are given as

$$ w^1 = \frac{p}{\rho a} - u, \quad \text{corresponding to } \lambda_1 = u - c, \qquad (A17) $$

$$ w^2 = \rho - \frac{p}{a^2}, \quad \text{corresponding to } \lambda_2 = u, \qquad (A18) $$

$$ w^3 = \frac{p}{\rho a} + u, \quad \text{corresponding to } \lambda_3 = u + c. \qquad (A19) $$

## A2 Boundary conditions

The primary idea behind the non-reflecting boundary condition is that any outgoing wave ($w^{\text{out}}$) that is leaving the computational domain through the boundary should leave the domain without getting reflected back into the interior region, which means the wave amplitude associated with the outgoing wave is constant at the boundary, i.e.

$$ w^{\text{out}} = \text{const}, \quad \text{or} \quad \Delta w^{\text{out}} = 0, \quad \text{at the boundary}, \qquad (A20) $$

and the boundary condition associated with the outgoing characteristic variables should be extrapolated from the interior solution. The other waves that are incoming into the computational domain require physical boundary conditions.

For a subsonic inflow boundary (let us assume x-beg), one wave ($w^1$) is leaving the domain ($\lambda_1 = u - a < 0$), and the other two ($w^2, w^3$) are incoming, which means only two out of three primitive variables (ideally ($\rho, u$) or ($\rho, p$) for well-poised conditions, e.g. see section 19.3 of Laney 1998) can be specified physically at the

boundary. The remaining variable has to be set numerically from the interior solution in such a way that $\Delta w^1 = 0$ at the boundary. Let's say, we want to specify the value of $\rho_B$ and $u_B$ physically, then the pressure value ($p_B$) at the boundary is calculated by setting

$$ \Delta w^1 = 0, \qquad (A21) $$

$$ w_B^1 - w_I^1 = 0, \qquad (A22) $$

$$ p_B = \gamma \rho_B (w_I^1 + u_B)^2, \qquad (A23) $$

where $w_I^1$ is the value of $w^1$ at the first active zone ($i=1$) of the computational domain, which is to be calculated from the interior solution. Therefore, for target density and velocity values of $\rho_t$ and $u_t$, the ghost cell values of the primitive variables are given by

$$ \rho_B = \rho_t, \qquad (A24) $$

$$ u_B = u_t, \qquad (A25) $$

$$ p_B = \gamma \rho_t (w_I^1 + u_t)^2. \qquad (A26) $$

Equations (A24)–(A26) ensure that $\Delta w^1 = 0$ at the boundary, and the wave-reflection is minimized.

## A3 Test and comparison

In order to demonstrate how the wind emerges in the interior solution with and without NRBC, we perform one simulation with NRBC implemented at the subsonic inflow boundary with the same setup as described in Section 2. The wind and cloud initializations are the same as in GC45.k3 in Table 1, with half the resolution of the fiducial runs. We compare the evolution of the wind parameters with a simulation that uses the same initial condition and resolution, but using the method of wind injection adopted in this paper (Section 2.6). In the simulation with NRBC, we specify the density and velocity values at the ghost zone, while the pressure of the wind is calculated using equation (A26).

In Fig. A1, we present the time evolution of the wind density (top left panel), velocity (top right panel), pressure (bottom left panel), and power (bottom right panel), at the first active cell in simulations with (red lines) and without (blue lines) NRBC. The lower section of each panel illustrates the fractional deviation of each parameter from its intended value. As observed, the evolution of the physically specified wind parameters – namely density and velocity – remains consistent throughout the simulations with NRBC. However, the pressure of the injected wind shows a significant deviation compared to the simulation without NRBC. Consequently, the deviation of wind power from its intended value is substantial (~90 per cent). On the other hand, although all wind parameters in the simulation without NRBC deviate from their intended values (blue lines), the cumulative effect results in more stable wind power (with a variation of ~30 per cent) compared to the scenario with NRBC. Moreover, we have performed additional simulations using the saturated values of the density and velocity for the GC45 case (see Table 2) with the NRBC. We confirmed that, although the quantitative results differ from the result presented in this study during the transient phase, they closely resemble the results (e.g. the diagnostics presented in Figs 9, 10, and 17) at the saturation phase ($t \gtrsim 0.4t_{cc}$). Therefore, the qualitative results presented in this study are robust.
<figure><img src="image_20.png" /><figcaption>**Figure A1.** Time evolution of the wind parameters for a 3D model problem of wind–cloud interaction without (blue) and with (red) an NRBC at the inflow boundary. The bottom section of each panel shows the fractional deviation from the injected values at the ghost zone.</figcaption></figure>

## APPENDIX B: DENSITY PDF OF THE CLOUD MATERIAL AT THE SAME CLOUD-CRUSHING TIME FOR DIFFERENT POWER

Fig. B1 depicts the density PDF of the self-gravitating cloud at $t = 0.4t_{cc}$ for simulations with different wind power. This time corresponds to an absolute time of 3.46, 1.38, 0.78, and 0.34 Myr for GC42, GC43, GC44, and GC45, respectively.

We notice a distinctly contrasting trend in the high-density tail of the PDF compared to what is depicted in Fig. 8. In this case, the cloud in the lowest power simulation (GC42, blue) has the highest fraction of cloud material in the PL tail. As wind power increases, this fraction as well as the maximum cloud density decreases due to an increasingly smaller absolute time for self-gravity to influence the evolution. Notably, in the highest power simulation (GC45), there is no evident signature of the PL tail at $0.4t_{cc}$. Corresponding to an approximate time of $0.4t_{cc} = 0.34\text{Myr} \approx 0.1t_{\text{ff}}$ in this instance, the evolution is at a very early stage in terms of the impact of self-gravity.

<figure><img src="image_21.png" /><figcaption>**Figure B1.** The PDFs of the cloud material in self-gravitating simulations with different wind power as indicated in the legend at $t = 0.4t_{cc}$. The vertical black-dotted lines in both panels mark the density value beyond which Jeans length cannot be resolved by at least four grid cells with the current computational setup.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure A1.** Time evolution of the wind parameters for a 3D model problem of wind–cloud interaction without (blue) and with (red) an NRBC at the inflow boundary. The bottom section of each panel shows the fractional deviation from the injected values at the ghost zone.</figcaption></figure>

## APPENDIX B: DENSITY PDF OF THE CLOUD MATERIAL AT THE SAME CLOUD-CRUSHING TIME FOR DIFFERENT POWER

Fig. B1 depicts the density PDF of the self-gravitating cloud at $t = 0.4t_{cc}$ for simulations with different wind power. This time corresponds to an absolute time of 3.46, 1.38, 0.78, and 0.34 Myr for GC42, GC43, GC44, and GC45, respectively.

We notice a distinctly contrasting trend in the high-density tail of the PDF compared to what is depicted in Fig. 8. In this case, the cloud in the lowest power simulation (GC42, blue) has the highest fraction of cloud material in the PL tail. As wind power increases, this fraction as well as the maximum cloud density decreases due to an increasingly smaller absolute time for self-gravity to influence the evolution. Notably, in the highest power simulation (GC45), there is no evident signature of the PL tail at $0.4t_{cc}$. Corresponding to an approximate time of $0.4t_{cc} = 0.34\text{Myr} \approx 0.1t_{\text{ff}}$ in this instance, the evolution is at a very early stage in terms of the impact of self-gravity.

<figure><img src="image_21.png" /><figcaption>**Figure B1.** The PDFs of the cloud material in self-gravitating simulations with different wind power as indicated in the legend at $t = 0.4t_{cc}$. The vertical black-dotted lines in both panels mark the density value beyond which Jeans length cannot be resolved by at least four grid cells with the current computational setup.</figcaption></figure>
## APPENDIX C: TURBULENT PROPERTY EVOLUTION IN TERMS OF CLOUD-CRUSHING TIME

Here, we examine how the turbulent properties of the clouds in the simulations with different wind power evolve with time in unit of the cloud-crushing time-scale.

Fig. C1, shows the evolution of velocity dispersion (left panel) and the virial parameter (right panel) of the cloud material as a function of cloud-crushing time for the fiducial cloud setup with different wind power. It is evident that in the absence of self-gravity (dashed line), the evolutionary trends are similar for all the powers. The velocity

<figure><img src="image_22.png" /><figcaption>Figure C1. Same as Fig. 10 but expressed in terms of $t_{cc}$.</figcaption></figure>

dispersion as well as the virial parameter peaks at ~0.2$t_{cc}$, although the magnitudes are different. However, in the self-gravity runs, the freefall time ($t_{ff}$) is a more relevant time-scale to compare to. As the freefall time of the cloud in the lower power case (GC.42) is shorter than $t_{cc}$, the cloud collapses early at around ~0.5$t_{cc}$. With increasing power, $t_{cc}$ of the cloud becomes comparable or shorter than $t_{ff}$. Therefore, the collapse of the clouds gets delayed in terms of $t_{cc}$ with increasing wind power.

## APPENDIX D: EFFECT OF NUMERICAL RESOLUTION

In order to quantify the numerical resolution dependence of the results presented in this study, we perform two simulations with the same wind and cloud initialization as in GC43_k3, but with decreasing resolution. The widths of a cell in the considered simulations are $\Delta x = 0.19, 0.39$, and $0.78$ pc and the initial cloud radius ($r_c$) is resolved by 128 (high), 64 (medium), and 32 (low) cells, respectively.

Fig. D1 illustrates the time evolution of the cold gas mass (top left panel), velocity dispersion (top right panel), virial parameter (bottom left panel), and SFR per freefall time (bottom right panel), for simulations with $r_c/\Delta x = 128$ (solid), 64 (dashed-dotted), and 32 (dotted). Across all panels, there is a clear dependence on resolution regarding the evolution of these quantities. While the evolution in the high (solid) and medium (dashed-dotted) resolution simulations closely align, the low (dotted) resolution one exhibits more significant deviations from the medium resolution case. This systematic dependence on resolution in the scenario of shock-cloud interaction directly stems from how well the instabilities acting on the cloud surface are resolved (Klein et al. 1994; Cooper et al. 2009; Yirak, Frank & Cunningham 2010; Banda-Barragán et al. 2016, 2018, 2020). As resolution increases, we can resolve perturbations with smaller wavelengths ($\lambda$), i.e. higher wavenumber ($k$). Since the growth rates for the Kelvin-Helmholtz and Rayleigh-Taylor instabilities are directly proportional to the wavenumber of the perturbation, the clouds in higher resolution simulations are more susceptible to instabilities and get ablated faster and mixed with the wind. This systematic variation affects all derived quantities. For instance, the amount of cold gas in the top left panel of Fig. D1 in the higher resolution simulation is consistently lower compared to the lower resolution case, due to the increased amount of ablation and mixing. This systematic variation is reflected in all the other panels as well. However, we find that the differences of the results between $r_c/\Delta x = 128$ and 64 are much less compared to the $r_c/\Delta x = 64$ and 32 pair, implying that the resolution of the fiducial simulations ($r_c/\Delta x = 128$) presented in this paper is adequate, as found by various previous studies (e.g. Klein et al. 1994; Fujita et al. 2009; Banda-Barragán et al. 2018).
## APPENDIX C: TURBULENT PROPERTY EVOLUTION IN TERMS OF CLOUD-CRUSHING TIME

Here, we examine how the turbulent properties of the clouds in the simulations with different wind power evolve with time in unit of the cloud-crushing time-scale.

Fig. C1, shows the evolution of velocity dispersion (left panel) and the virial parameter (right panel) of the cloud material as a function of cloud-crushing time for the fiducial cloud setup with different wind power. It is evident that in the absence of self-gravity (dashed line), the evolutionary trends are similar for all the powers. The velocity

<figure><img src="image_22.png" /><figcaption>Figure C1. Same as Fig. 10 but expressed in terms of $t_{cc}$.</figcaption></figure>

dispersion as well as the virial parameter peaks at ~0.2$t_{cc}$, although the magnitudes are different. However, in the self-gravity runs, the freefall time ($t_{ff}$) is a more relevant time-scale to compare to. As the freefall time of the cloud in the lower power case (GC.42) is shorter than $t_{cc}$, the cloud collapses early at around ~0.5$t_{cc}$. With increasing power, $t_{cc}$ of the cloud becomes comparable or shorter than $t_{ff}$. Therefore, the collapse of the clouds gets delayed in terms of $t_{cc}$ with increasing wind power.

## APPENDIX D: EFFECT OF NUMERICAL RESOLUTION

In order to quantify the numerical resolution dependence of the results presented in this study, we perform two simulations with the same wind and cloud initialization as in GC43_k3, but with decreasing resolution. The widths of a cell in the considered simulations are $\Delta x = 0.19$, 0.39, and 0.78 pc and the initial cloud radius ($r_c$) is resolved by 128 (high), 64 (medium), and 32 (low) cells, respectively.

Fig. D1 illustrates the time evolution of the cold gas mass (top left panel), velocity dispersion (top right panel), virial parameter (bottom left panel), and SFR per freefall time (bottom right panel), for simulations with $r_c/\Delta x = 128$ (solid), 64 (dashed-dotted), and 32 (dotted). Across all panels, there is a clear dependence on resolution regarding the evolution of these quantities. While the evolution in the high (solid) and medium (dashed-dotted) resolution simulations closely align, the low (dotted) resolution one exhibits more significant deviations from the medium resolution case. This systematic dependence on resolution in the scenario of shock-cloud interaction directly stems from how well the instabilities acting on the cloud surface are resolved (Klein et al. 1994; Cooper et al. 2009; Yirak, Frank & Cunningham 2010; Banda-Barragán et al. 2016, 2018, 2020). As resolution increases, we can resolve perturbations with smaller wavelengths ($\lambda$), i.e. higher wavenumber ($k$). Since the growth rates for the Kelvin-Helmholtz and Rayleigh-Taylor instabilities are directly proportional to the wavenumber of the perturbation, the clouds in higher resolution simulations are more susceptible to instabilities and get ablated faster and mixed with the wind. This systematic variation affects all derived quantities. For instance, the amount of cold gas in the top left panel of Fig. D1 in the higher resolution simulation is consistently lower compared to the lower resolution case, due to the increased amount of ablation and mixing. This systematic variation is reflected in all the other panels as well. However, we find that the differences of the results between $r_c/\Delta x = 128$ and 64 are much less compared to the $r_c/\Delta x = 64$ and 32 pair, implying that the resolution of the fiducial simulations ($r_c/\Delta x = 128$) presented in this paper is adequate, as found by various previous studies (e.g. Klein et al. 1994; Fujita et al. 2009; Banda-Barragán et al. 2018).
<figure><img src="image_23.png" /><figcaption>**Figure D1.** Time evolution of the cold gas mass (top left panel), velocity dispersion (top right panel), virial parameter (bottom left panel), and SFR per freefall time (bottom right panel), for simulations with wind power of $10^{43}$ erg s$^{-1}$, but with different numerical grid resolutions. The widths of the cells in these three simulations are $\Delta x = 0.19$ pc (solid), 0.39 pc (dashed-dotted), and 0.78 pc (dotted), and the initial cloud radius ($r_c$) is resolved by 128, 64, and 32 cells, respectively.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.
<figure><img src="image_23.png" /><figcaption>**Figure D1.** Time evolution of the cold gas mass (top left panel), velocity dispersion (top right panel), virial parameter (bottom left panel), and SFR per freefall time (bottom right panel), for simulations with wind power of $10^{43}$ erg s$^{-1}$, but with different numerical grid resolutions. The widths of the cells in these three simulations are $\Delta x = 0.19$ pc (solid), 0.39 pc (dashed-dotted), and 0.78 pc (dotted), and the initial cloud radius ($r_c$) is resolved by 128, 64, and 32 cells, respectively.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.