# Probing the role of self-gravity in clouds impacted by AGN-driven winds

Ankush Mandal, Dipanjan Mukherjee, Christoph Federrath, Geoffrey V Bicknell, Nicole P.H Nesvadba, Andrea Mignone

► To cite this version:

Ankush Mandal, Dipanjan Mukherjee, Christoph Federrath, Geoffrey V Bicknell, Nicole P.H Nesvadba, et al. Probing the role of self-gravity in clouds impacted by AGN-driven winds. Mon.Not.Roy.Astron.Soc., 2024, 531 (1), pp.2079-2110. 10.1093/mnras/stae1295. hal-04600041

HAL Id: hal-04600041

https://hal.science/hal-04600041

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing the role of self-gravity in clouds impacted by AGN-driven winds

Ankush Mandal, Dipanjan Mukherjee, Christoph Federrath, Geoffrey V Bicknell, Nicole P.H Nesvadba, Andrea Mignone

► To cite this version:

Ankush Mandal, Dipanjan Mukherjee, Christoph Federrath, Geoffrey V Bicknell, Nicole P.H Nesvadba, et al.. Probing the role of self-gravity in clouds impacted by AGN-driven winds. Monthly Notices of the Royal Astronomical Society, 2024, 531 (1), pp.2079-2110. 10.1093/mnras/stae1295. hal-04600041

HAL Id: hal-04600041

https://hal.science/hal-04600041v1

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing the role of self-gravity in clouds impacted by AGN-driven winds

Ankush Mandal<sup>1,★</sup> Dipanjan Mukherjee<sup>1,★</sup> Christoph Federrath,<sup>2,3</sup> Geoffrey V. Bicknell<sup>1,2</sup> Nicole P. H. Nesvadba<sup>4</sup> and Andrea Mignone<sup>5</sup>

Accepted 2024 May 16. Received 2024 May 16; in original form 2024 January 18

## ABSTRACT

The impact of winds and jet-inflated bubbles driven by active galactic nuclei (AGN) are believed to significantly affect the host galaxy’s interstellar medium (ISM) and regulate star formation. To explore this scenario, we perform a suite of hydrodynamic simulations to model the interaction between turbulent star-forming clouds and highly pressurized AGN-driven outflows, focusing on the effects of self-gravity. Our results demonstrate that the cloudlets fragmented by the wind can become gravitationally bound, significantly increasing their survival time. While external pressurization leads to a global collapse of the clouds in cases of weaker winds ($10^{42}-10^{43}$ erg s⁻<fnref n="1" />), higher power winds ($10^{44}-10^{45}$ erg s⁻<fnref n="1" />) disperse the gas and cause localized collapse of the cloudlets. We also demonstrate that a kinetic energy-dominated wind is more efficient in accelerating and dispersing the gas than a thermal wind with the same power. The interaction can give rise to multiphase outflows with velocities ranging from a few 100 to several 1000 km s⁻<fnref n="1" />. The mass outflow rates are tightly correlated with the wind power, which we explain by an ablation-based mass-loss model. Moreover, the velocity dispersion and the virial parameter of the cloud material can increase by up to one order of magnitude through the effect of the wind. Even though the wind can suppress or quench star formation for about 1 Myr during the initial interaction, a substantial number of gravitationally bound dense cloudlets manage to shield themselves from the wind’s influence and subsequently undergo rapid gravitational collapse, leading to an enhanced star formation rate.

**Key words:** hydrodynamics – ISM: clouds – galaxies: active – galaxies: star formation – methods: numerical.

## 1 INTRODUCTION

The feedback from active galactic nuclei (AGNs) on the overall evolution of their host galaxies is thought to be a dominant mechanism in galaxy evolution theory (Silk & Rees 1998; Fabian 2012). It is postulated that the large-scale outflow in the form of ‘jets’ from the AGN heats up the intracluster medium and stops the cooling flow towards the centre of the cluster, therefore regulating the star-forming fuel (McNamara & Nulsen 2007, 2012). Indeed, in recent cosmological simulations, it is necessary to include various models of feedback from the AGN by injecting thermal or kinetic energy (Springel, Di Matteo & Hernquist 2005; Schaye et al. 2015; Weinberger et al. 2017; Davé et al. 2019; Schaye et al. 2023), in order to regulate star formation in massive galaxies and reproduce various observed scaling relations, including the luminosity functions and the $M_{BH}-\sigma$ relation (see the recent review Vogelsberger et al. 2020, and references therein).

While modern cosmological simulations successfully replicate the statistical characteristics and the redshift evolution of galaxies, they lack the ability to predict the feedback’s impact on individual host galaxies due to challenges in accurately modelling the multiphase

interstellar medium (ISM) and associated small-scale physics, i.e. these processes are included as sub-grid recipes in cosmological simulations. There is increasing observational evidence that wind and young radio jets originating from the central AGN significantly affect the host galaxy’s ISM by driving multiphase outflows, which expel gas from the central region, driving turbulence, and potentially diminishing the star-forming fuel (Nesvadba et al. 2010; García-Burillo et al. 2014; Harrison et al. 2014; Fluetsch et al. 2019; Girdhar et al. 2022; Ramos Almeida et al. 2022; Leftley et al. 2024). This is also demonstrated by dedicated hydrodynamic (HD) simulations (Sutherland & Bicknell 2007; Wagner, Bicknell & Umemura 2012; Mukherjee et al. 2018; Mandal et al. 2021; Tanner & Weaver 2022; Meenakshi et al. 2022b). These phenomena have a direct impact on the star formation activity inside the host as demonstrated by several observational studies where it has been found that some galaxies hosting radio-loud AGN show a lower star formation rate (SFR) compared to main-sequence galaxies, which follow the standard Kennicutt–Schmidt relation (Schmidt 1959, 1963; Kennicutt 1998a,b) between gas–mass and SFR surface density (Ogle et al. 2007, 2010; Nesvadba et al. 2010, 2011, 2021; Alatalo et al. 2014, 2015; Lanz et al. 2016).

Conversely, the overpressurized winds/jet can cause significant compression of the ISM and may trigger collapse to rapidly form stars (Silk 2005; Gaibler et al. 2012; Zubovas & King 2014;

<footnote n="1"><sup>1</sup>Inter-University Centre for Astronomy and Astrophysics, Pune 411007, India</footnote>

<footnote n="2"><sup>2</sup>Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia</footnote>

<footnote n="3"><sup>3</sup>ARC Centre of Excellence for Astronomy in Three Dimensions (ASTRO-3D), Canberra, ACT 2601, Australia</footnote>

<footnote n="4"><sup>4</sup>Université de la Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France</footnote>

<footnote n="5"><sup>5</sup>Dipartimento di Fisica, Università di Torino via Pietro Giuria, 1 (I-10125) Torino, Italy</footnote>

<footnote>*E-mail: ankush@iucaa.in (AM); dipanjan@iucaa.in (DM)</footnote>
# Probing the role of self-gravity in clouds impacted by AGN-driven winds

Ankush Mandal<sup>1,★</sup> Dipanjan Mukherjee<sup>1,★</sup> Christoph Federrath,<sup>2,3</sup> Geoffrey V. Bicknell<sup>1,2</sup> Nicole P. H. Nesvadba<sup>4</sup> and Andrea Mignone<sup>5</sup>

Accepted 2024 May 16. Received 2024 May 16; in original form 2024 January 18

## ABSTRACT

The impact of winds and jet-inflated bubbles driven by active galactic nuclei (AGN) are believed to significantly affect the host galaxy’s interstellar medium (ISM) and regulate star formation. To explore this scenario, we perform a suite of hydrodynamic simulations to model the interaction between turbulent star-forming clouds and highly pressurized AGN-driven outflows, focusing on the effects of self-gravity. Our results demonstrate that the cloudlets fragmented by the wind can become gravitationally bound, significantly increasing their survival time. While external pressurization leads to a global collapse of the clouds in cases of weaker winds ($10^{42}$–$10^{43}$ erg s⁻<fnref n="1" />), higher power winds ($10^{44}$–$10^{45}$ erg s⁻<fnref n="1" />) disperse the gas and cause localized collapse of the cloudlets. We also demonstrate that a kinetic energy-dominated wind is more efficient in accelerating and dispersing the gas than a thermal wind with the same power. The interaction can give rise to multiphase outflows with velocities ranging from a few 100 to several 1000 km s⁻<fnref n="1" />. The mass outflow rates are tightly correlated with the wind power, which we explain by an ablation-based mass-loss model. Moreover, the velocity dispersion and the virial parameter of the cloud material can increase by up to one order of magnitude through the effect of the wind. Even though the wind can suppress or quench star formation for about 1 Myr during the initial interaction, a substantial number of gravitationally bound dense cloudlets manage to shield themselves from the wind’s influence and subsequently undergo rapid gravitational collapse, leading to an enhanced star formation rate.

**Key words:** hydrodynamics – ISM: clouds – galaxies: active – galaxies: star formation – methods: numerical.

## 1 INTRODUCTION

The feedback from active galactic nuclei (AGNs) on the overall evolution of their host galaxies is thought to be a dominant mechanism in galaxy evolution theory (Silk & Rees 1998; Fabian 2012). It is postulated that the large-scale outflow in the form of ‘jets’ from the AGN heats up the intracluster medium and stops the cooling flow towards the centre of the cluster, therefore regulating the star-forming fuel (McNamara & Nulsen 2007, 2012). Indeed, in recent cosmological simulations, it is necessary to include various models of feedback from the AGN by injecting thermal or kinetic energy (Springel, Di Matteo & Hernquist 2005; Schaye et al. 2015; Weinberger et al. 2017; Davé et al. 2019; Schaye et al. 2023), in order to regulate star formation in massive galaxies and reproduce various observed scaling relations, including the luminosity functions and the $M_{BH}-\sigma$ relation (see the recent review Vogelsberger et al. 2020, and references therein).

While modern cosmological simulations successfully replicate the statistical characteristics and the redshift evolution of galaxies, they lack the ability to predict the feedback’s impact on individual host galaxies due to challenges in accurately modelling the multiphase

interstellar medium (ISM) and associated small-scale physics, i.e. these processes are included as sub-grid recipes in cosmological simulations. There is increasing observational evidence that wind and young radio jets originating from the central AGN significantly affect the host galaxy’s ISM by driving multiphase outflows, which expel gas from the central region, driving turbulence, and potentially diminishing the star-forming fuel (Nesvadba et al. 2010; García-Burillo et al. 2014; Harrison et al. 2014; Fluetsch et al. 2019; Girdhar et al. 2022; Ramos Almeida et al. 2022; Leftley et al. 2024). This is also demonstrated by dedicated hydrodynamic (HD) simulations (Sutherland & Bicknell 2007; Wagner, Bicknell & Umemura 2012; Mukherjee et al. 2018; Mandal et al. 2021; Tanner & Weaver 2022; Meenakshi et al. 2022b). These phenomena have a direct impact on the star formation activity inside the host as demonstrated by several observational studies where it has been found that some galaxies hosting radio-loud AGN show a lower star formation rate (SFR) compared to main-sequence galaxies, which follow the standard Kennicutt–Schmidt relation (Schmidt 1959, 1963; Kennicutt 1998a,b) between gas–mass and SFR surface density (Ogle et al. 2007, 2010; Nesvadba et al. 2010, 2011, 2021; Alatalo et al. 2014, 2015; Lanz et al. 2016).

Conversely, the overpressurized winds/jet can cause significant compression of the ISM and may trigger collapse to rapidly form stars (Silk 2005; Gaibler et al. 2012; Zubovas & King 2014;

<footnote n="1"><sup>1</sup>Inter-University Centre for Astronomy and Astrophysics, Pune 411007, India</footnote>

<footnote n="2"><sup>2</sup>Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia</footnote>

<footnote n="3"><sup>3</sup>ARC Centre of Excellence for Astronomy in Three Dimensions (ASTRO-3D), Canberra, ACT 2601, Australia</footnote>

<footnote n="4"><sup>4</sup>Université de la Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice cedex 4, France</footnote>

<footnote n="5"><sup>5</sup>Dipartimento di Fisica, Università di Torino via Pietro Giuria, 1 (I-10125) Torino, Italy</footnote>

<footnote>*E-mail: ankush@iucaa.in (AM); dipanjan@iucaa.in (DM)</footnote>
Dugan, Gaibler & Silk 2017b; Mukherjee et al. 2018). Observational evidence also supports this hypothesis, where compact radio jets or quasar winds are found to enhance star formation activity (Bicknell et al. 2000; Inskip et al. 2008; Zinn et al. 2013; Kalfountzou et al. 2014; Salomé, Salomé & Combes 2015; Lacy et al. 2017; Salomé et al. 2017).

Advancements in observational techniques and improved modelling of star formation physics within HD simulations are beginning to shed light on the distinction between 'negative' and 'positive' feedback from AGN. Recent observations indicate the coexistence of both types of feedback within a single system (Cresci & Maiolino 2018; Shin et al. 2019). Indeed, recent hydrodynamical simulations of jet-ISM interactions have revealed that while jet-inflated bubbles globally reduce star formation by enhancing turbulence, they can cause local regions of enhanced SFR due to the compression near the nuclear region (Mandal et al. 2021; Mercedes-Feliz et al. 2023). Thus, how AGN-driven outflows affect star formation is a complex competition between various phenomena on different scales.

However, a complete understanding of star formation as well as the survivability of the dense gas subjected to powerful AGN outflows remains elusive without the effect of the self-gravity of the gas. With typical densities around $100 \text{ cm}^{-3}$ (Miville-Deschênes, Murray & Lee 2017), star-forming giant molecular clouds (GMCs) have a freefall time-scale of a few Myr, which is comparable to or shorter than the typical duration of AGN episodes, lasting between 10 and 100 Myr (e.g. Marconi et al. 2004). Moreover, the presence of self-gravity can increase/prolong the survival time of the clouds, when faced with strong outflows from AGN, by making them dense and compact, effectively shielding them from erosion caused by the outflow and/or AGN radiation. Conversely, the fragmentation induced by self-gravity can give rise to numerous smaller cloudlets that may be susceptible to evaporation or entrainment by the hot wind/jet cocoon, leading to the formation of multiphase outflows, which may regulate the available fuel to form stars.

Therefore, the significance of self-gravity at the cloud level can influence how AGN-driven winds/jet cocoons impact the host galaxy on larger scales. Thus, well-resolved simulations modelling the interaction between AGN-driven winds and individual star-forming clouds may offer a supplementary perspective to both observations and global-scale simulations. After all, the ultimate fate of the clouds depends on small-scale processes. Additionally, the results from these small-scale studies are important for building better sub-resolution prescriptions of different mechanisms in global (galaxy and cosmological) simulations.

In this study, we revisit the classical 'cloud-crushing' problem (Klein, McKee & Colella 1994) with the help of a suite of three-dimensional (3D), self-gravitational HDs simulations in the context of the interaction between AGN-driven winds/jet cocoons and star-forming clouds. There have been extensive studies of the effects of external shocks or winds on clouds in various different contexts, with a primary emphasis on supersonic winds/shocks from galactic winds (e.g. Klein et al. 1994; Nakamura et al. 2006; Pittard et al. 2009; Pittard, Hartquist & Falle 2010; Scannapieco & Brüggen 2015; Banda-Barragán et al. 2016, 2018, 2019, 2020, 2021; Cottle, Scannapieco & Brüggen 2018; Gronke & Oh 2018; Cottle et al. 2020). However, wind-driven bubbles or jet-cocoons are known to be highly pressurized during the energy-driven phase (e.g. Begelman & Cioffi 1989; Wagner & Bicknell 2011), and thus can be subsonic depending on the density of the wind, while also exhibiting extreme velocities of up to tens of thousands of km s$^{-1}$. None the less, only a limited number of studies have taken into account the parameters of shocks/winds (e.g. density, velocity, and pressure) which can reach

extremes comparable to those generated by AGN-jet cocoons or quasar winds (Mellema, Kurk & Röttgering 2002; Fragile et al. 2004; Cooper et al. 2009, 2018; Dugan et al. 2017a).

Mellema et al. (2002) examined the influence of a radio jet cocoon on uniform spherical and elliptical clouds using 2D simulations. They identified three significant phases in the evolution (Klein et al. 1994) of these interactions: (i) the initial impact between the blast wave and the cloud, (ii) compression induced by the thermal and ram pressure of the wind, and (iii) fragmentation of the cloud. Additionally, their work revealed that in the presence of radiative cooling, the growth rate of the Kelvin-Helmholtz instability was highly suppressed and the mixed gas fraction was considerably reduced (less than 1 per cent of the original cloud mass), implying a slow evaporation process and a prolonged lifetime of the cloudlets. Cooper et al. (2009) reached a similar conclusion and further demonstrated that when dealing with a fractal cloud, the fragmentation induced by the wind is more pronounced compared to a uniform cloud structure. The fragmented cloudlets are compressed to high densities, which can then cool very efficiently and survive for a much longer period. While they discussed the potential role of self-gravity in this scenario, self-gravity was not included explicitly in the simulations.

Considering the effect of external pressurization on a spherical cloud that may result from AGN winds or jet cocoons, Zubovas & King (2014) demonstrated that the pressure confinement triggers the collapse of the cloud, leading to an enhanced SFR. Using adaptive mesh refinement (AMR), self-gravitating simulations of the interaction between more realistic AGN-driven winds and Bonnor-Ebert spheres (resembling star-forming cores), Dugan et al. (2017a) also arrived at a similar conclusion. Additionally, they identified a threshold ram pressure of the wind, above which the cloud will be destroyed before significant amounts of star formation can take place. Another recent study by Li et al. (2020) concluded that, if the cloud is initially Jeans' unstable, the interaction will eventually enhance the collapse of the cloud in the presence of self-gravity. While these studies have individually delved into various significant aspects and distinct physical processes, there is a lack of a comprehensive global perspective that incorporates all the crucial parameters and physics.

This study seeks to extend previous research on the interaction between AGN-driven outflows with more realistic fractal star-forming interstellar clouds, including radiative cooling and self-gravity. We explore a wide parameter space, systematically varying parameters including the wind power, the average cloud density, the internal fractal density distribution within the cloud, and whether the wind is primarily dominated by kinetic or thermal energy. This approach enables us to investigate diverse facets of the cloud's evolution. The paper is organized as follows. In Section 2, we describe the simulation method and the choice of initial conditions. In Section 3, we present the main results of this study. We discuss the implications of this work in Section 4. Finally, in Section 5, we summarize and conclude.

## 2 METHOD

### 2.1 Simulation code

The numerical simulations presented in this study are performed using the grid-based code PLUTO v4.4 (Mignone et al. 2007, 2012) in 3D (x, y, z) Cartesian geometry. We use the HLLC Riemann solver (Gurski 2004; Li 2005) along with a piecewise parabolic reconstruction scheme (PPM; Colella & Woodward 1984) for solving the self-gravitating HD equations:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (1) $$
Dugan, Gaibler & Silk 2017b; Mukherjee et al. 2018). Observational evidence also supports this hypothesis, where compact radio jets or quasar winds are found to enhance star formation activity (Bicknell et al. 2000; Inskip et al. 2008; Zinn et al. 2013; Kalfountzou et al. 2014; Salomé, Salomé & Combes 2015; Lacy et al. 2017; Salomé et al. 2017).

Advancements in observational techniques and improved modelling of star formation physics within HD simulations are beginning to shed light on the distinction between ‘negative’ and ‘positive’ feedback from AGN. Recent observations indicate the coexistence of both types of feedback within a single system (Cresci & Maiolino 2018; Shin et al. 2019). Indeed, recent hydrodynamical simulations of jet-ISM interactions have revealed that while jet-inflated bubbles globally reduce star formation by enhancing turbulence, they can cause local regions of enhanced SFR due to the compression near the nuclear region (Mandal et al. 2021; Mercedes-Feliz et al. 2023). Thus, how AGN-driven outflows affect star formation is a complex competition between various phenomena on different scales.

However, a complete understanding of star formation as well as the survivability of the dense gas subjected to powerful AGN outflows remains elusive without the effect of the self-gravity of the gas. With typical densities around $100 \text{ cm}^{-3}$ (Miville-Deschênes, Murray & Lee 2017), star-forming giant molecular clouds (GMCs) have a freefall time-scale of a few Myr, which is comparable to or shorter than the typical duration of AGN episodes, lasting between 10 and 100 Myr (e.g. Marconi et al. 2004). Moreover, the presence of self-gravity can increase/prolong the survival time of the clouds, when faced with strong outflows from AGN, by making them dense and compact, effectively shielding them from erosion caused by the outflow and/or AGN radiation. Conversely, the fragmentation induced by self-gravity can give rise to numerous smaller cloudlets that may be susceptible to evaporation or entrainment by the hot wind/jet cocoon, leading to the formation of multiphase outflows, which may regulate the available fuel to form stars.

Therefore, the significance of self-gravity at the cloud level can influence how AGN-driven winds/jet cocoons impact the host galaxy on larger scales. Thus, well-resolved simulations modelling the interaction between AGN-driven winds and individual star-forming clouds may offer a supplementary perspective to both observations and global-scale simulations. After all, the ultimate fate of the clouds depends on small-scale processes. Additionally, the results from these small-scale studies are important for building better sub-resolution prescriptions of different mechanisms in global (galaxy and cosmological) simulations.

In this study, we revisit the classical ‘cloud-crushing’ problem (Klein, McKee & Colella 1994) with the help of a suite of three-dimensional (3D), self-gravitational HDs simulations in the context of the interaction between AGN-driven winds/jet cocoons and star-forming clouds. There have been extensive studies of the effects of external shocks or winds on clouds in various different contexts, with a primary emphasis on supersonic winds/shocks from galactic winds (e.g. Klein et al. 1994; Nakamura et al. 2006; Pittard et al. 2009; Pittard, Hartquist & Falle 2010; Scannapieco & Brüggen 2015; Banda-Barragán et al. 2016, 2018, 2019, 2020, 2021; Cottle, Scannapieco & Brüggen 2018; Gronke & Oh 2018; Cottle et al. 2020). However, wind-driven bubbles or jet-cocoons are known to be highly pressurized during the energy-driven phase (e.g. Begelman & Cioffi 1989; Wagner & Bicknell 2011), and thus can be subsonic depending on the density of the wind, while also exhibiting extreme velocities of up to tens of thousands of km s$^{-1}$. None the less, only a limited number of studies have taken into account the parameters of shocks/winds (e.g. density, velocity, and pressure) which can reach

extremes comparable to those generated by AGN-jet cocoons or quasar winds (Mellema, Kurk & Röttgering 2002; Fragile et al. 2004; Cooper et al. 2009, 2018; Dugan et al. 2017a).

Mellema et al. (2002) examined the influence of a radio jet cocoon on uniform spherical and elliptical clouds using 2D simulations. They identified three significant phases in the evolution (Klein et al. 1994) of these interactions: (i) the initial impact between the blast wave and the cloud, (ii) compression induced by the thermal and ram pressure of the wind, and (iii) fragmentation of the cloud. Additionally, their work revealed that in the presence of radiative cooling, the growth rate of the Kelvin-Helmholtz instability was highly suppressed and the mixed gas fraction was considerably reduced (less than 1 per cent of the original cloud mass), implying a slow evaporation process and a prolonged lifetime of the cloudlets. Cooper et al. (2009) reached a similar conclusion and further demonstrated that when dealing with a fractal cloud, the fragmentation induced by the wind is more pronounced compared to a uniform cloud structure. The fragmented cloudlets are compressed to high densities, which can then cool very efficiently and survive for a much longer period. While they discussed the potential role of self-gravity in this scenario, self-gravity was not included explicitly in the simulations.

Considering the effect of external pressurization on a spherical cloud that may result from AGN winds or jet cocoons, Zubovas & King (2014) demonstrated that the pressure confinement triggers the collapse of the cloud, leading to an enhanced SFR. Using adaptive mesh refinement (AMR), self-gravitating simulations of the interaction between more realistic AGN-driven winds and Bonnor-Ebert spheres (resembling star-forming cores), Dugan et al. (2017a) also arrived at a similar conclusion. Additionally, they identified a threshold ram pressure of the wind, above which the cloud will be destroyed before significant amounts of star formation can take place. Another recent study by Li et al. (2020) concluded that, if the cloud is initially Jeans’ unstable, the interaction will eventually enhance the collapse of the cloud in the presence of self-gravity. While these studies have individually delved into various significant aspects and distinct physical processes, there is a lack of a comprehensive global perspective that incorporates all the crucial parameters and physics.

This study seeks to extend previous research on the interaction between AGN-driven outflows with more realistic fractal star-forming interstellar clouds, including radiative cooling and self-gravity. We explore a wide parameter space, systematically varying parameters including the wind power, the average cloud density, the internal fractal density distribution within the cloud, and whether the wind is primarily dominated by kinetic or thermal energy. This approach enables us to investigate diverse facets of the cloud’s evolution. The paper is organized as follows. In Section 2, we describe the simulation method and the choice of initial conditions. In Section 3, we present the main results of this study. We discuss the implications of this work in Section 4. Finally, in Section 5, we summarize and conclude.

# 2 METHOD

## 2.1 Simulation code

The numerical simulations presented in this study are performed using the grid-based code PLUTO v4.4 (Mignone et al. 2007, 2012) in 3D (x, y, z) Cartesian geometry. We use the HLLC Riemann solver (Gurski 2004; Li 2005) along with a piecewise parabolic reconstruction scheme (PPM; Colella & Woodward 1984) for solving the self-gravitating HD equations:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (1) $$
$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot [\rho\mathbf{v}\mathbf{v} + p\mathbf{I}] = \rho\mathbf{g}, \qquad (2) $$

$$ \frac{\partial E_t}{\partial t} + \nabla \cdot [(E_t + p)\mathbf{v}] = \rho\mathbf{v} \cdot \mathbf{g} - \left(\frac{\rho}{\mu m_H}\right)^2 \Lambda(T), \qquad (3) $$

$$ \frac{\partial(\rho C)}{\partial t} + \nabla \cdot \rho C \mathbf{v} = 0, \qquad (4) $$

$$ \nabla^2 \Phi = 4\pi G \rho, \qquad (5) $$

where $\rho$ is the mass density, $v$ is the velocity, $p$ is the thermal pressure, $C$ is a Lagrangian scalar used to track gas in different components (i.e. cloud and wind), $g = -\nabla\Phi$ (where $\Phi$ is the gravitational potential) is the acceleration due to gravity, and $E_t$ is the total energy density given by

$$ E_t = \rho e + \frac{\rho \mathbf{v}^2}{2}. \qquad (6) $$

The above equations are closed by an ideal gas equation of state (EOS):

$$ p = (\gamma - 1)\rho e, \qquad (7) $$

where we consider $\gamma = 5/3$ throughout. The energy conservation equation (equation 3) also includes the cooling term ($\Lambda$) in order to account for the radiative losses, discussed in Section 2.5. Time evolution of equations (1)–(4) is performed using a third-order Runge–Kutta time-stepping scheme. For simulations involving gravity, we incorporate the self-gravity module<fnref n="1" /> developed by Mandal, Mukherjee & Mignone (2023), which employs a Runge–Kutta–Legendre-based Poisson solver coupled to a V-cycle multigrid algorithm to solve the Poisson equation for the gravitational potential. Details of the numerical implementation of the self-gravity module are presented in Mandal et al. (2023). We also solve for the potential in the non-self-gravitating runs but do not couple the gravitational acceleration terms to the HDs (equations 2 and 3). In this way, we calculate the gravitational potential-related quantities for these simulations and compare them with the corresponding self-gravitating runs.

## 2.2 Computational domain

In our simulations, we employ a uniform Cartesian domain ($x-y-z$) with a physical range of $-50 \le x \le 150$ pc, $-50 \le y \le 50$ pc, and $-50 \le z \le 50$ pc. The domain is discretized into a grid with $1024 \times 512 \times 512$ cells, resulting in a computational cell size (resolution) of 0.195 pc. In order to study the effect of the numerical resolution, we also perform two lower resolution simulations with a grid size of $512 \times 256 \times 256$ and $256 \times 128 \times 128$, respectively, and present the results in Appendix D. The fractal clouds (see Section 2.3) are initially positioned at the origin (0, 0, 0) of the domain. The clouds possess a core radius of 25 pc and an envelope of width 5 pc. Therefore, the core radius of the cloud is resolved with ~128 cells, this being adequate to capture the overall evolution (Klein et al. 1994; Fujita et al. 2009; Banda-Barragán et al. 2018). The wind is launched from the $y-z$ boundary at $x = -50$ pc in the positive $x$-direction. The details of the cloud setup are described in Section 2.3 and the wind parameters and their injection in Sections 2.4 and 2.6, respectively.

## 2.3 Cloud initialization

Several theoretical and observational studies have shown that the density structures in the turbulent molecular cloud are well described by a lognormal distribution function (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Federrath et al. 2010a; Price, Federrath & Brunt 2011; Federrath & Klessen 2012; Federrath & Banerjee 2015; Kritsuk, Ustyugov & Norman 2017; Mandal, Federrath & Körtgen 2020). Hence, we model the cloud in our simulation such that the PDF in terms of the logarithmic density $s = \ln(\rho/\rho_0)$ (where $\rho_0$ is the mean density of the cloud) is given by a lognormal distribution:

$$ P(s) = \frac{1}{\sqrt{2\pi\sigma_s^2}} \exp\left[-\frac{(s-s_0)^2}{2\sigma_s^2}\right]. \qquad (8) $$

Here, $s_0$ and $\sigma_s$ are the mean and dispersion of the logarithmic density fluctuation, respectively, which from normalization constraints ($\int e^s P(s) ds = 1$) are related by $s_0 = -\sigma_s^2/2$ (Ostriker, Stone & Gammie 2001; Li, Klessen & Mac Low 2003; Federrath & Klessen 2012).

The lognormal density field of the cloud is constructed using the pyFC library,<fnref n="2" /> which generates a periodic random scalar field in 3D Cartesian space from a given PDF and power-law (PL) spectrum, $D(k) \propto k^{-\beta}$, in Fourier space, where $k$ is the dimensionless wavenumber. The two-point fractal distribution is characterized by the slope of the PL ($\beta$), the Nyquist limit $k_{max}$ and a lower cutoff wavenumber, $k_{min}$, which corresponds to the largest spatially correlated scale ($\lambda_{max} \approx L/k_{min}$) for a positive value of $\beta$. For a given value of $k_{min}$, the largest size of the perturbations or ‘cloudlets’ is $r_{cloudlets}, k_{min} \approx L/(2k_{min})$, where $L$ is the size of the periodic box (Lewis & Austin 2002; Sutherland & Bicknell 2007; Wagner et al. 2012). In this study, we have set the value of $\beta$ to be 1.66 for all the clouds. This falls in the range of cloud density spectral indices for supersonic turbulence (e.g. see Federrath, Klessen & Schmidt 2009; Federrath & Klessen 2013). The $k_{min}$ value is primarily set to 3 ($\lambda_{max} \approx 20$ pc) for most of the clouds while varying the wind properties to investigate their effects. Nevertheless, we also explore different values of $k_{min} = 1, 3, 6$, and 10, while keeping the wind parameters identical. This allows us to examine the impact of variations in the density distribution within the cloud.

In addition to the fractal density distribution of the cloud, we also initialize a Gaussian random field for each component of the velocity with zero mean and 1D velocity dispersion of $\sigma_v$. For a particular cloud setup, the values of $k_{min}$ and $\beta$ for the velocity field are the same as the lognormal density field. We choose the value of $\sigma_v$ for most of the simulations (except for the cloud with lower mean density) such that the 3D velocity dispersion ($\sigma_v = \sqrt{3}\sigma_v$) of the cloud is ~8 km s⁻<fnref n="1" />, which is typical of the observed velocity dispersion of GMCs in the Milky Way and nearby galaxies on this scale (e.g. Hughes et al. 2010, 2013; Miville-Deschênes et al. 2017; Faesi, Lada & Forbrich 2018).

We set the mean number density of the fractal cloud to 200 cm⁻³, which gives us an initial virial parameter ($\alpha_{vir} = 2E_{kin}/E_{grav}$) of ~0.9, again typical for star-forming clouds (Miville-Deschênes et al. 2017; Faesi et al. 2018).

The standard deviation of the density PDF in equation (8) of our clouds is calculated using the well-established relation:

$$ \sigma_s \approx [\ln(1 + b^2 M^2)]^{1/2} \qquad (9) $$

<footnote n="1">¹The self-gravity patch for the PLUTO code is publicly available at https://bitbucket.org/mankush/pluto-4.4-self-gravity-patch</footnote>

<footnote n="2">²https://www2.ccs.tsukuba.ac.jp/Astro/Members/ayw/code/pyFC</footnote>
$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot [\rho\mathbf{v}\mathbf{v} + p\mathbf{I}] = \rho\mathbf{g}, \qquad (2) $$

$$ \frac{\partial E_t}{\partial t} + \nabla \cdot [(E_t + p)\mathbf{v}] = \rho\mathbf{v} \cdot \mathbf{g} - \left(\frac{\rho}{\mu m_H}\right)^2 \Lambda(T), \qquad (3) $$

$$ \frac{\partial(\rho C)}{\partial t} + \nabla \cdot \rho C \mathbf{v} = 0, \qquad (4) $$

$$ \nabla^2 \Phi = 4\pi G \rho, \qquad (5) $$

where $\rho$ is the mass density, $v$ is the velocity, $p$ is the thermal pressure, $C$ is a Lagrangian scalar used to track gas in different components (i.e. cloud and wind), $\mathbf{g} = -\nabla\Phi$ (where $\Phi$ is the gravitational potential) is the acceleration due to gravity, and $E_t$ is the total energy density given by

$$ E_t = \rho e + \frac{\rho \mathbf{v}^2}{2}. \qquad (6) $$

The above equations are closed by an ideal gas equation of state (EOS):

$$ p = (\gamma - 1)\rho e, \qquad (7) $$

where we consider $\gamma = 5/3$ throughout. The energy conservation equation (equation 3) also includes the cooling term ($\Lambda$) in order to account for the radiative losses, discussed in Section 2.5. Time evolution of equations (1)–(4) is performed using a third-order Runge-Kutta time-stepping scheme. For simulations involving gravity, we incorporate the self-gravity module<fnref n="1" /> developed by Mandal, Mukherjee & Mignone (2023), which employs a Runge-Kutta-Legendre-based Poisson solver coupled to a V-cycle multigrid algorithm to solve the Poisson equation for the gravitational potential. Details of the numerical implementation of the self-gravity module are presented in Mandal et al. (2023). We also solve for the potential in the non-self-gravitating runs but do not couple the gravitational acceleration terms to the HDs (equations 2 and 3). In this way, we calculate the gravitational potential-related quantities for these simulations and compare them with the corresponding self-gravitating runs.

## 2.2 Computational domain

In our simulations, we employ a uniform Cartesian domain ($x-y-z$) with a physical range of $-50 \le x \le 150$ pc, $-50 \le y \le 50$ pc, and $-50 \le z \le 50$ pc. The domain is discretized into a grid with $1024 \times 512 \times 512$ cells, resulting in a computational cell size (resolution) of 0.195 pc. In order to study the effect of the numerical resolution, we also perform two lower resolution simulations with a grid size of $512 \times 256 \times 256$ and $256 \times 128 \times 128$, respectively, and present the results in Appendix D. The fractal clouds (see Section 2.3) are initially positioned at the origin (0, 0, 0) of the domain. The clouds possess a core radius of 25 pc and an envelope of width 5 pc. Therefore, the core radius of the cloud is resolved with $\sim 128$ cells, this being adequate to capture the overall evolution (Klein et al. 1994; Fujita et al. 2009; Banda-Barragán et al. 2018). The wind is launched from the $y-z$ boundary at $x = -50$ pc in the positive $x$-direction. The details of the cloud setup are described in Section 2.3 and the wind parameters and their injection in Sections 2.4 and 2.6, respectively.

## 2.3 Cloud initialization

Several theoretical and observational studies have shown that the density structures in the turbulent molecular cloud are well described by a lognormal distribution function (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Federrath et al. 2010a; Price, Federrath & Brunt 2011; Federrath & Klessen 2012; Federrath & Banerjee 2015; Kritsuk, Ustyugov & Norman 2017; Mandal, Federrath & Körtgen 2020). Hence, we model the cloud in our simulation such that the PDF in terms of the logarithmic density $s = \ln(\rho/\rho_0)$ (where $\rho_0$ is the mean density of the cloud) is given by a lognormal distribution:

$$ P(s) = \frac{1}{\sqrt{2\pi\sigma_s^2}} \exp\left[-\frac{(s - s_0)^2}{2\sigma_s^2}\right]. \qquad (8) $$

Here, $s_0$ and $\sigma_s$ are the mean and dispersion of the logarithmic density fluctuation, respectively, which from normalization constraints ($\int e^s P(s) ds = 1$) are related by $s_0 = -\sigma_s^2/2$ (Ostriker, Stone & Gammie 2001; Li, Klessen & Mac Low 2003; Federrath & Klessen 2012).

The lognormal density field of the cloud is constructed using the pyFC library,<fnref n="2" /> which generates a periodic random scalar field in 3D Cartesian space from a given PDF and power-law (PL) spectrum, $D(k) \propto k^{-\beta}$, in Fourier space, where $k$ is the dimensionless wavenumber. The two-point fractal distribution is characterized by the slope of the PL ($\beta$), the Nyquist limit $k_{max}$ and a lower cutoff wavenumber, $k_{min}$, which corresponds to the largest spatially correlated scale ($\lambda_{max} \approx L/k_{min}$) for a positive value of $\beta$. For a given value of $k_{min}$, the largest size of the perturbations or ‘cloudlets’ is $r_{cloudlets}, k_{min} \approx L/(2k_{min})$, where $L$ is the size of the periodic box (Lewis & Austin 2002; Sutherland & Bicknell 2007; Wagner et al. 2012). In this study, we have set the value of $\beta$ to be 1.66 for all the clouds. This falls in the range of cloud density spectral indices for supersonic turbulence (e.g. see Federrath, Klessen & Schmidt 2009; Federrath & Klessen 2013). The $k_{min}$ value is primarily set to 3 ($\lambda_{max} \approx 20$ pc) for most of the clouds while varying the wind properties to investigate their effects. Nevertheless, we also explore different values of $k_{min} = 1, 3, 6$, and 10, while keeping the wind parameters identical. This allows us to examine the impact of variations in the density distribution within the cloud.

In addition to the fractal density distribution of the cloud, we also initialize a Gaussian random field for each component of the velocity with zero mean and 1D velocity dispersion of $\sigma_v$. For a particular cloud setup, the values of $k_{min}$ and $\beta$ for the velocity field are the same as the lognormal density field. We choose the value of $\sigma_v$ for most of the simulations (except for the cloud with lower mean density) such that the 3D velocity dispersion ($\sigma_v = \sqrt{3}\sigma_v$) of the cloud is $\sim 8$ km s$^{-1}$, which is typical of the observed velocity dispersion of GMCs in the Milky Way and nearby galaxies on this scale (e.g. Hughes et al. 2010, 2013; Miville-Deschênes et al. 2017; Faesi, Lada & Forbrich 2018).

We set the mean number density of the fractal cloud to $200$ cm$^{-3}$, which gives us an initial virial parameter ($\alpha_{vir} = 2E_{kin}/E_{grav}$) of $\sim 0.9$, again typical for star-forming clouds (Miville-Deschênes et al. 2017; Faesi et al. 2018).

The standard deviation of the density PDF in equation (8) of our clouds is calculated using the well-established relation:

$$ \sigma_s \approx [\ln(1 + b^2 M^2)]^{1/2} \qquad (9) $$

<footnote n="1">¹The self-gravity patch for the PLUTO code is publicly available at https://bitbucket.org/mankush/pluto-4.4-self-gravity-patch</footnote>

<footnote n="2">²https://www2.ccs.tsukuba.ac.jp/Astro/Members/ayw/code/pyFC</footnote>
which connects the standard deviation of the log-density ($\sigma_s$) and the turbulent Mach number ($M = \sigma_v/c_{s,rms}$, where $c_{s,rms}$ is the root mean square sound speed) (Federrath, Klessen & Schmidt 2008; Federrath et al. 2010a; Federrath & Klessen 2012). The parameter $b$ is the driving parameter of turbulence, which represents the ratios of energies in solenoidal and compressive modes turbulence and varies between 1/3 (purely solenoidal) and 1 (purely compressive), respectively (Federrath et al. 2008). Here, we set $b = 0.4$ for a mixed mode of turbulence driving (Federrath et al. 2010a), as often observed in different environments (Federrath et al. 2016; Menon et al. 2021; Dhawalikar et al. 2022; Sharda et al. 2022; Gerrard et al. 2023). The resulting value of $\sigma_s$ then serves as the input parameter in the pyFC routine which generates the fractal density field.

For the cloud density, we create a 310³-sized lognormal data cube with a mean of 1, which corresponds to a physical domain of $(x, y, z) \in [-30 \text{ pc}, 30 \text{ pc}]^3$ following the method described above. The density cube is then multiplied by the desired mean density of the cloud and a spherical volume of a radius of 30 pc is extracted from the cube. The sphere is then tapered with a radially decreasing function in order to ensure a smooth transition between the cloud's edge and the ambient medium:

$$ \rho_c(r) = \rho_a + \frac{\rho_{\text{cube}}(r)}{\cosh\left[\left(\frac{r}{r_{\text{core}}}\right)^8\right]}, \quad (10) $$

where $\rho_a$, $\rho_{\text{cube}}$, and $\rho_c$ are the ambient density, the original density values in the fractal data cube, and the final density for the cloud material that is used for the simulations, respectively. Here, $r_{\text{core}} = 25 \text{ pc}$ is the core radius of the cloud where the density remains the same as the density cube, along with an envelope with 5 pc radially decreasing density. Finally, the turbulent cloud is placed at the origin (0, 0, 0) of the computational domain using a trilinear interpolation scheme.

We also perform the same procedure for each velocity component where a 310³-sized Gaussian random data cube is generated for each component of the velocity field and mapped into the computational grid in a similar way to the density initialization. The remaining portion of the computational domain is initialized with a static ambient medium having a density of $0.1 \text{ cm}^{-3}$ and a temperature of $10^6 \text{ K}$, while the cloud is set to be in pressure equilibrium with this ambient medium initially.

## 2.4 Wind parameters

In this study, our main objective is to examine how AGN-driven 'fast' winds or jet-inflated cocoons affect star-forming complexes within host galaxies. Therefore, we focus on parameters typical of such scenarios, particularly considering the pressure and velocity ranges associated with these AGN-driven processes. Numerous observational studies have demonstrated that the velocities of ionized winds (Rupke & Veilleux 2013; Genzel et al. 2014; Harrison et al. 2014; Brusa et al. 2015; Carniani et al. 2015; Bischetti et al. 2017) and broad absorption line (BAL) winds (Korista et al. 2008; Moe et al. 2009; Borguet et al. 2013) driven by the central AGN on kiloparsec scales range from a few hundred to several thousand km s⁻¹ (see Fiore et al. 2017, for a review). In this work, we consider the cloud to be located at ~1 kpc away from the AGN and the velocity of

the injected winds at ($v_w$) this distance are in the range of 400 to 4000 km s⁻¹, with a total (kinetic + thermal) wind power ($P_w$) ranging from $10^{42}$ to $10^{45}$ erg s⁻¹. From theoretical and numerical investigations, these types of winds are found to be very hot and highly pressurized during the energy-conserving phase with pressure ranging from $10^{-10}$ to $10^{-7}$ dyne cm⁻² depending on the wind power (e.g. Begelman & Cioffi 1989; Sutherland & Bicknell 2007; Wagner & Bicknell 2011; Mukherjee et al. 2016, 2018; Richings & Faucher-Giguère 2018; Costa, Pakmor & Springel 2020). Within these parameter ranges the thermal energy of the winds dominates the total energy budget, making them subsonic. Interestingly, the self-similar solution describing an expanding bubble propelled by a central source, as proposed by Weaver et al. (1977), aligns with these pressure values. For an AGN-driven bubble with a given injected power ($P_w$), expanding in a spatially homogeneous ambient medium of density $\rho_a$, the pressure ($p_w$) at a distance $R_w$ from the central source is expressed as

$$ \begin{aligned} p_w = & \frac{7}{25} \left(\frac{125}{154\pi}\right)^{2/3} \rho_a^{1/3} P_w^{2/3} R_w^{-4/3}, \\ & \approx 5.5 \times 10^{-10} \left(\frac{n_a}{0.1 \text{ cm}^{-3}}\right)^{1/3} \left(\frac{P_w}{10^{43} \text{ erg s}^{-1}}\right)^{2/3} \\ & \times \left(\frac{R_w}{\text{kpc}}\right)^{-4/3} \text{dyne cm}^{-2}. \end{aligned} \quad (11) $$

In this study, we set the pressure of the wind from the above equation (equation 11) for $R_w \sim 1 \text{ kpc}$, which yields values in range $p_w \sim 10^{-10}–10^{-8}$ dyne cm⁻². The pressure is kept constant in time at the injection region, unlike the true Weaver et al. (1977) wind solution. We discuss the implications of this choice in later paragraphs. Moreover, the velocity of the bubble's forward shock at $R_w = 1 \text{ kpc}$, as given by Weaver et al. (1977),

$$ \begin{aligned} v_w = & \left(\frac{243}{3850\pi}\right)^{1/5} \rho_a^{-1/5} P_w^{1/5} t^{-2/5}, \\ & \approx 850 \left(\frac{n_a}{0.1 \text{ cm}^{-3}}\right)^{-1/3} \left(\frac{P_w}{10^{43} \text{ erg s}^{-1}}\right)^{1/3} \left(\frac{R_w}{\text{kpc}}\right)^{-2/3} \text{km s}^{-1}, \end{aligned} \quad (12) $$

also predicts values within the considered velocity ranges ~400–4000 km s⁻¹. In Table 1, we list the wind parameters for each simulation.

We also consider a wind of power $10^{45}$ erg s⁻¹ whose velocity is set to 15 000 km s⁻¹ and whose pressure is such that it is supersonic. The total energy, in this case, is dominated by the kinetic energy (see Table 1). This kind of extreme outflow parameter can appear in the scenario where a cloud directly lies along the path of an AGN jet. Additionally, this choice allows us to explore the effect of thermal versus kinetic (subsonic versus supersonic) winds of the same power in an otherwise identical cloud setup.

In order to completely specify the wind state, the wind density is calculated from energy conservation, i.e. the total energy flux of the thermal and kinetic components are equal to the injected power of the wind,

$$ \left(\frac{1}{2} \rho_w v_w^2 + \frac{P_w}{\gamma - 1}\right) 4\pi R_w^2 v_w = P_w. \quad (13) $$

<footnote n="3">³For this choice, the scale largest correlated density structure for a given value of $k_{\text{min}}$ is $\lambda_{\text{max},k_{\text{min}}} \approx (60/k_{\text{min}})$ pc.</footnote>
which connects the standard deviation of the log-density ($\sigma_s$) and the turbulent Mach number ($M = \sigma_v/c_{s,rms}$, where $c_{s,rms}$ is the root mean square sound speed) (Federrath, Klessen & Schmidt 2008; Federrath et al. 2010a; Federrath & Klessen 2012). The parameter $b$ is the driving parameter of turbulence, which represents the ratios of energies in solenoidal and compressive modes turbulence and varies between 1/3 (purely solenoidal) and 1 (purely compressive), respectively (Federrath et al. 2008). Here, we set $b = 0.4$ for a mixed mode of turbulence driving (Federrath et al. 2010a), as often observed in different environments (Federrath et al. 2016; Menon et al. 2021; Dhawalikar et al. 2022; Sharda et al. 2022; Gerrard et al. 2023). The resulting value of $\sigma_s$ then serves as the input parameter in the pyFC routine which generates the fractal density field.

For the cloud density, we create a 310³-sized lognormal data cube with a mean of 1, which corresponds to a physical domain of $(x, y, z) \in [-30 \text{ pc}, 30 \text{ pc}]^3$ following the method described above. The density cube is then multiplied by the desired mean density of the cloud and a spherical volume of a radius of 30 pc is extracted from the cube. The sphere is then tapered with a radially decreasing function in order to ensure a smooth transition between the cloud's edge and the ambient medium:

$$ \rho_c(r) = \rho_a + \frac{\rho_{\text{cube}}(r)}{\cosh\left[\left(\frac{r}{r_{\text{core}}}\right)^8\right]}, \quad (10) $$

where $\rho_a$, $\rho_{\text{cube}}$, and $\rho_c$ are the ambient density, the original density values in the fractal data cube, and the final density for the cloud material that is used for the simulations, respectively. Here, $r_{\text{core}} = 25 \text{ pc}$ is the core radius of the cloud where the density remains the same as the density cube, along with an envelope with 5 pc radially decreasing density. Finally, the turbulent cloud is placed at the origin (0, 0, 0) of the computational domain using a trilinear interpolation scheme.

We also perform the same procedure for each velocity component where a 310³-sized Gaussian random data cube is generated for each component of the velocity field and mapped into the computational grid in a similar way to the density initialization. The remaining portion of the computational domain is initialized with a static ambient medium having a density of $0.1 \text{ cm}^{-3}$ and a temperature of $10^6 \text{ K}$, while the cloud is set to be in pressure equilibrium with this ambient medium initially.

## 2.4 Wind parameters

In this study, our main objective is to examine how AGN-driven 'fast' winds or jet-inflated cocoons affect star-forming complexes within host galaxies. Therefore, we focus on parameters typical of such scenarios, particularly considering the pressure and velocity ranges associated with these AGN-driven processes. Numerous observational studies have demonstrated that the velocities of ionized winds (Rupke & Veilleux 2013; Genzel et al. 2014; Harrison et al. 2014; Brusa et al. 2015; Carniani et al. 2015; Bischetti et al. 2017) and broad absorption line (BAL) winds (Korista et al. 2008; Moe et al. 2009; Borguet et al. 2013) driven by the central AGN on kiloparsec scales range from a few hundred to several thousand km s⁻¹ (see Fiore et al. 2017, for a review). In this work, we consider the cloud to be located at ~1 kpc away from the AGN and the velocity of

the injected winds at ($v_w$) this distance are in the range of 400 to 4000 km s⁻¹, with a total (kinetic + thermal) wind power ($P_w$) ranging from $10^{42}$ to $10^{45}$ erg s⁻¹. From theoretical and numerical investigations, these types of winds are found to be very hot and highly pressurized during the energy-conserving phase with pressure ranging from $10^{-10}$ to $10^{-7}$ dyne cm⁻² depending on the wind power (e.g. Begelman & Cioffi 1989; Sutherland & Bicknell 2007; Wagner & Bicknell 2011; Mukherjee et al. 2016, 2018; Richings & Faucher-Giguère 2018; Costa, Pakmor & Springel 2020). Within these parameter ranges the thermal energy of the winds dominates the total energy budget, making them subsonic. Interestingly, the self-similar solution describing an expanding bubble propelled by a central source, as proposed by Weaver et al. (1977), aligns with these pressure values. For an AGN-driven bubble with a given injected power ($P_w$), expanding in a spatially homogeneous ambient medium of density $\rho_a$, the pressure ($p_w$) at a distance $R_w$ from the central source is expressed as

$$ \begin{aligned} p_w = & \frac{7}{25} \left( \frac{125}{154\pi} \right)^{2/3} \rho_a^{1/3} P_w^{2/3} R_w^{-4/3}, \\ & \approx 5.5 \times 10^{-10} \left( \frac{n_a}{0.1 \text{ cm}^{-3}} \right)^{1/3} \left( \frac{P_w}{10^{43} \text{ erg s}^{-1}} \right)^{2/3} \\ & \times \left( \frac{R_w}{\text{kpc}} \right)^{-4/3} \text{dyne cm}^{-2}. \end{aligned} \quad (11) $$

In this study, we set the pressure of the wind from the above equation (equation 11) for $R_w \sim 1 \text{ kpc}$, which yields values in range $p_w \sim 10^{-10}–10^{-8}$ dyne cm⁻². The pressure is kept constant in time at the injection region, unlike the true Weaver et al. (1977) wind solution. We discuss the implications of this choice in later paragraphs. Moreover, the velocity of the bubble's forward shock at $R_w = 1 \text{ kpc}$, as given by Weaver et al. (1977),

$$ \begin{aligned} v_w = & \left( \frac{243}{3850\pi} \right)^{1/5} \rho_a^{-1/5} P_w^{1/5} t^{-2/5}, \\ & \approx 850 \left( \frac{n_a}{0.1 \text{ cm}^{-3}} \right)^{-1/3} \left( \frac{P_w}{10^{43} \text{ erg s}^{-1}} \right)^{1/3} \left( \frac{R_w}{\text{kpc}} \right)^{-2/3} \text{km s}^{-1}, \end{aligned} \quad (12) $$

also predicts values within the considered velocity ranges ~400–4000 km s⁻¹. In Table 1, we list the wind parameters for each simulation.

We also consider a wind of power $10^{45}$ erg s⁻¹ whose velocity is set to 15 000 km s⁻¹ and whose pressure is such that it is supersonic. The total energy, in this case, is dominated by the kinetic energy (see Table 1). This kind of extreme outflow parameter can appear in the scenario where a cloud directly lies along the path of an AGN jet. Additionally, this choice allows us to explore the effect of thermal versus kinetic (subsonic versus supersonic) winds of the same power in an otherwise identical cloud setup.

In order to completely specify the wind state, the wind density is calculated from energy conservation, i.e. the total energy flux of the thermal and kinetic components are equal to the injected power of the wind,

$$ \left( \frac{1}{2} \rho_w v_w^2 + \frac{P_w}{\gamma - 1} \right) 4\pi R_w^2 v_w = P_w. \quad (13) $$

<footnote n="3">³For this choice, the scale largest correlated density structure for a given value of $k_{\text{min}}$ is $\lambda_{\text{max},k_{\text{min}}} \approx (60/k_{\text{min}})$ pc.</footnote>
<table><caption>**Table 1.** Initial conditions of all the simulations in this study.</caption><thead><tr><th>Name</th><th>P<sub>w</sub><sup>(1)</sup><br>(erg s<sup>-1</sup>)</th><th>v<sub>w</sub><sup>(2)</sup><br>(km s<sup>-1</sup>)</th><th>p<sub>w</sub><sup>(3)</sup><br>(dyne cm<sup>-2</sup>)</th><th>M<sub>w</sub><sup>(4)</sup></th><th>n&#773;<sub>c</sub><sup>(5)</sup><br>(cm<sup>-3</sup>)</th><th>Self-gravity?</th><th>k<sub>min</sub><sup>(6)</sup></th><th>t<sub>cc</sub><sup>(7)</sup><br>(Myr)</th><th>t<sub>fl</sub><sup>(8)</sup><br>(Myr)</th><th>t<sub>KH</sub><sup>(9)</sup><br>(Myr)</th></tr></thead><tbody><tr><td>C_no wind.k3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>200</td><td>No</td><td>3</td><td>—</td><td>3.65</td><td>—</td></tr><tr><td>GC_no wind.k3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>200</td><td>Yes</td><td>3</td><td>—</td><td>3.65</td><td>—</td></tr><tr><td>C42.k3</td><td>10<sup>42</sup></td><td>400</td><td>1.32 &times; 10<sup>-10</sup></td><td>0.286</td><td>200</td><td>No</td><td>3</td><td>8.64</td><td>3.65</td><td>4.65</td></tr><tr><td>GC42.k3</td><td>10<sup>42</sup></td><td>400</td><td>1.32 &times; 10<sup>-10</sup></td><td>0.286</td><td>200</td><td>Yes</td><td>3</td><td>8.64</td><td>3.65</td><td>4.65</td></tr><tr><td>C43.k3</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>No</td><td>3</td><td>3.46</td><td>3.65</td><td>1.86</td></tr><tr><td>GC43.k3</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>3</td><td>3.46</td><td>3.65</td><td>1.86</td></tr><tr><td>C44.k3</td><td>10<sup>44</sup></td><td>1800</td><td>2.95 &times; 10<sup>-9</sup></td><td>0.277</td><td>200</td><td>No</td><td>3</td><td>1.92</td><td>3.65</td><td>1.03</td></tr><tr><td>GC44.k3</td><td>10<sup>44</sup></td><td>1800</td><td>2.95 &times; 10<sup>-9</sup></td><td>0.277</td><td>200</td><td>Yes</td><td>3</td><td>1.92</td><td>3.65</td><td>1.03</td></tr><tr><td>C45.k3</td><td>10<sup>45</sup></td><td>4000</td><td>1.32 &times; 10<sup>-8</sup></td><td>0.285</td><td>200</td><td>No</td><td>3</td><td>0.86</td><td>3.65</td><td>0.46</td></tr><tr><td>GC45.k3</td><td>10<sup>45</sup></td><td>4000</td><td>1.32 &times; 10<sup>-8</sup></td><td>0.285</td><td>200</td><td>Yes</td><td>3</td><td>0.86</td><td>3.65</td><td>0.46</td></tr><tr><td>GC43.k1</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>1</td><td>3.46</td><td>3.65</td><td>5.58</td></tr><tr><td>GC43.k6</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>6</td><td>3.46</td><td>3.65</td><td>0.93</td></tr><tr><td>GC43.k10</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>10</td><td>3.46</td><td>3.65</td><td>0.56</td></tr><tr><td>GC43_uniform</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>—</td><td>3.46</td><td>3.65</td><td>4.61</td></tr><tr><td>GC43.k3_low</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>20</td><td>Yes</td><td>3</td><td>1.09</td><td>11.55</td><td>0.59</td></tr><tr><td>GC45.k3_kinetic</td><td>10<sup>45</sup></td><td>15 000</td><td>1.00 &times; 10<sup>-10</sup></td><td>12.16</td><td>200</td><td>Yes</td><td>3</td><td>0.23</td><td>3.65</td><td>0.13</td></tr></tbody></table>

Notes. (1) Power of the expanding bubble. (2) Speed of the wind material. (3) Pressure of the wind. (4) Mach number of the wind defined as $M_w = v_w/c_{s,w}$, where $c_{s,w}$ is the sound speed inside the wind. (5) Mean number density of the cloud. (6) The minimum normalized wavenumber ($k_{min} \equiv k$) of the density distribution of the cloud. (7) Cloud-crushing time-scale. (8) Freefall time-scale. (9) Growth time-scale for the Kelvin-Helmholtz instabilities.

Therefore, for the given values of wind power ($P_w$), velocity ($v_w$), and pressure ($p_w$), the density ($n_w$) of the wind is given by

$$n_w \approx 3 \times 10^{-2} \left( \frac{v_w}{10^3 \text{ km s}^{-1}} \right)^{-2} \left[ 5.57 \left( \frac{P_w}{10^{43} \text{ erg s}^{-1}} \right) \times \left( \frac{v_w}{10^3 \text{ km s}^{-1}} \right)^{-1} \left( \frac{R_w}{\text{kpc}} \right)^{-2} - \left( \frac{p_w}{10^{-10} \text{ dyne cm}^{-2}} \right) \right] \text{cm}^{-3}, \quad (14)$$

which is typically $n_w \approx 0.01 \text{ cm}^{-3}$ for the wind powers with corresponding velocity and pressure considered in this study. Therefore, we set the density of the wind to $n_w = 0.01 \text{ cm}^{-3}$ for all the simulations. With these wind parameters, the mass outflow rate over a surface of $4\pi R_w^2$ at $R_w = 1$ kpc lies within the range of ~0.7–7 M$_\odot$ yr$^{-1}$, which is typical of the mass outflow rates of ionized winds found in observations (e.g. see Fiore et al. 2017, and references therein).

Although we set the pressure and velocity of the winds with different powers assuming that the cloud is located at $R_w = 1$ kpc from the AGN, from equations (11) and (12), it becomes apparent that similar values of $p_w$ and $v_w$ can be achieved for a wind of different power when considering an alternative $R_w$ value. For instance, $p_w$ and $v_w$ are very similar for winds with $P_w = 10^{44} \text{ erg s}^{-1}$ at $R_w = 1$ kpc and $P_w = 10^{42} \text{ erg s}^{-1}$, but at $R_w = 0.1$ kpc. Therefore, the wind parameters considered in this study not only reflect the strength of the wind at a particular location (~1 kpc) for varying wind powers, but can also be mapped to different locations relative to the central AGN for a particular wind power.

It is essential to emphasize that Weaver et al. (1977) solution is only used to set the reference values of wind parameters, which are kept constant throughout the simulation. In reality, the radius of the bubble ($R_w$) increases with time as it expands. Therefore, the wind solutions are time-dependent. Ideally, one should consider a self-consistent evolution of the bubble with time for the wind injection. Additionally, the bubble solution of an AGN-driven wind solution consists of several distinct internal structures, i.e. the forward shock, shocked ambient medium, shocked wind and wind material (e.g. Weaver et al. 1977; King 2003; Faucher-Giguère & Quataert 2012;

Zubovas & King 2012; Costa, Sijacki & Haehnelt 2014), which are not generally easy to implement in local box simulations like ours. Moreover, as our main focus is to investigate the effect of a steady wind whose parameters are similar to AGN-driven outflows on kpc scales on star-forming complexes, we adopt the simplistic wind injection model. Importantly, the main purpose of this study is to investigate the impact of the wind on clouds with and without self-gravity. Thus, while details of the wind modelling are simplified, we can still make a meaningful comparison, as clouds in the simulation with and without self-gravity face the exact same wind.

## 2.5 Cooling

In order to account for the energy losses due to radiative cooling (see equation 3), we use the non-equilibrium cooling function calculated using MAPPINGS V code (Sutherland & Dopita 2017; Sutherland et al. 2018). This code utilizes a comprehensive database of atomic data to self-consistently compute the optically thin cooling rate for various gas phases, including cold neutral, warm neutral, partially ionized, and fully ionized gas in the temperature range $10^2-10^9$ K. For temperatures exceeding ~$10^9$ K, the cooling function is extended by assuming bremsstrahlung emission (e.g. Krause & Alexander 2007; Mukherjee et al. 2018). We also impose a lower temperature threshold of $T_{\text{floor}} = 100$ K, below which the cooling is turned off and the temperature of any cell falling below $T_{\text{floor}}$ is enforced to stay at $T_{\text{floor}}$. No additional heating terms, such as the galactic ultraviolet (UV) background or ionizing photons (UV, soft, and hard X-ray) from the AGN, are included. In our simulations, the gas that can cool below $10^4$ K is sufficiently dense ($n_c > 500 \text{ cm}^{-3}$). In this density range, the high-density cores can be safely assumed to be self-shielded from the external UV and soft X-ray photons (Rahmati et al. 2013; Meenakshi et al. 2022a). Conversely, due to the very small photoabsorption cross-section of hard X-ray photons ($E \gtrsim 20$ keV, primarily due to K-shell ionization of Fe and Ni ions) (Band et al. 1990; Verner et al. 1993), the cloudlets remain optically thin, resulting in negligible heating. Through explicit radiative-transfer calculations using CLOUDY (Ferland et al. 2017) with the AGN spectral energy distribution (SED) as used in Meenakshi et al.
<table><caption>**Table 1.** Initial conditions of all the simulations in this study.</caption><thead><tr><th>Name</th><th>P<sub>w</sub><sup>(1)</sup><br>(erg s<sup>-1</sup>)</th><th>v<sub>w</sub><sup>(2)</sup><br>(km s<sup>-1</sup>)</th><th>p<sub>w</sub><sup>(3)</sup><br>(dyne cm<sup>-2</sup>)</th><th>M<sub>w</sub><sup>(4)</sup></th><th>n&#773;<sub>c</sub><sup>(5)</sup><br>(cm<sup>-3</sup>)</th><th>Self-gravity?</th><th>k<sub>min</sub><sup>(6)</sup></th><th>t<sub>cc</sub><sup>(7)</sup><br>(Myr)</th><th>t<sub>fl</sub><sup>(8)</sup><br>(Myr)</th><th>t<sub>KH</sub><sup>(9)</sup><br>(Myr)</th></tr></thead><tbody><tr><td>C_no wind.k3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>200</td><td>No</td><td>3</td><td>—</td><td>3.65</td><td>—</td></tr><tr><td>GC_no wind.k3</td><td>—</td><td>—</td><td>—</td><td>—</td><td>200</td><td>Yes</td><td>3</td><td>—</td><td>3.65</td><td>—</td></tr><tr><td>C42.k3</td><td>10<sup>42</sup></td><td>400</td><td>1.32 &times; 10<sup>-10</sup></td><td>0.286</td><td>200</td><td>No</td><td>3</td><td>8.64</td><td>3.65</td><td>4.65</td></tr><tr><td>GC42.k3</td><td>10<sup>42</sup></td><td>400</td><td>1.32 &times; 10<sup>-10</sup></td><td>0.286</td><td>200</td><td>Yes</td><td>3</td><td>8.64</td><td>3.65</td><td>4.65</td></tr><tr><td>C43.k3</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>No</td><td>3</td><td>3.46</td><td>3.65</td><td>1.86</td></tr><tr><td>GC43.k3</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>3</td><td>3.46</td><td>3.65</td><td>1.86</td></tr><tr><td>C44.k3</td><td>10<sup>44</sup></td><td>1800</td><td>2.95 &times; 10<sup>-9</sup></td><td>0.277</td><td>200</td><td>No</td><td>3</td><td>1.92</td><td>3.65</td><td>1.03</td></tr><tr><td>GC44.k3</td><td>10<sup>44</sup></td><td>1800</td><td>2.95 &times; 10<sup>-9</sup></td><td>0.277</td><td>200</td><td>Yes</td><td>3</td><td>1.92</td><td>3.65</td><td>1.03</td></tr><tr><td>C45.k3</td><td>10<sup>45</sup></td><td>4000</td><td>1.32 &times; 10<sup>-8</sup></td><td>0.285</td><td>200</td><td>No</td><td>3</td><td>0.86</td><td>3.65</td><td>0.46</td></tr><tr><td>GC45.k3</td><td>10<sup>45</sup></td><td>4000</td><td>1.32 &times; 10<sup>-8</sup></td><td>0.285</td><td>200</td><td>Yes</td><td>3</td><td>0.86</td><td>3.65</td><td>0.46</td></tr><tr><td>GC43.k1</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>1</td><td>3.46</td><td>3.65</td><td>5.58</td></tr><tr><td>GC43.k6</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>6</td><td>3.46</td><td>3.65</td><td>0.93</td></tr><tr><td>GC43.k10</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>10</td><td>3.46</td><td>3.65</td><td>0.56</td></tr><tr><td>GC43_uniform</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>200</td><td>Yes</td><td>—</td><td>3.46</td><td>3.65</td><td>4.61</td></tr><tr><td>GC43.k3_low</td><td>10<sup>43</sup></td><td>1000</td><td>5.25 &times; 10<sup>-10</sup></td><td>0.329</td><td>20</td><td>Yes</td><td>3</td><td>1.09</td><td>11.55</td><td>0.59</td></tr><tr><td>GC45.k3_kinetic</td><td>10<sup>45</sup></td><td>15 000</td><td>1.00 &times; 10<sup>-10</sup></td><td>12.16</td><td>200</td><td>Yes</td><td>3</td><td>0.23</td><td>3.65</td><td>0.13</td></tr></tbody></table>

Notes. (1) Power of the expanding bubble. (2) Speed of the wind material. (3) Pressure of the wind. (4) Mach number of the wind defined as $M_w = v_w/c_s w$, where $c_s w$ is the sound speed inside the wind. (5) Mean number density of the cloud. (6) The minimum normalized wavenumber ($k_{min} \equiv k$) of the density distribution of the cloud. (7) Cloud-crushing time-scale. (8) Freefall time-scale. (9) Growth time-scale for the Kelvin-Helmholtz instabilities.

Therefore, for the given values of wind power ($P_w$), velocity ($v_w$), and pressure ($p_w$), the density ($n_w$) of the wind is given by

$$n_w \approx 3 \times 10^{-2} \left( \frac{v_w}{10^3 \text{ km s}^{-1}} \right)^{-2} \left[ 5.57 \left( \frac{P_w}{10^{43} \text{ erg s}^{-1}} \right) \times \left( \frac{v_w}{10^3 \text{ km s}^{-1}} \right)^{-1} \left( \frac{R_w}{\text{kpc}} \right)^{-2} - \left( \frac{P_w}{10^{-10} \text{ dyne cm}^{-2}} \right) \right] \text{cm}^{-3}, \quad (14)$$

which is typically $n_w \approx 0.01 \text{ cm}^{-3}$ for the wind powers with corresponding velocity and pressure considered in this study. Therefore, we set the density of the wind to $n_w = 0.01 \text{ cm}^{-3}$ for all the simulations. With these wind parameters, the mass outflow rate over a surface of $4\pi R_w^2$ at $R_w = 1$ kpc lies within the range of ~0.7–7 M$_\odot$ yr$^{-1}$, which is typical of the mass outflow rates of ionized winds found in observations (e.g. see Fiore et al. 2017, and references therein).

Although we set the pressure and velocity of the winds with different powers assuming that the cloud is located at $R_w = 1$ kpc from the AGN, from equations (11) and (12), it becomes apparent that similar values of $p_w$ and $v_w$ can be achieved for a wind of different power when considering an alternative $R_w$ value. For instance, $p_w$ and $v_w$ are very similar for winds with $P_w = 10^{44} \text{ erg s}^{-1}$ at $R_w = 1$ kpc and $P_w = 10^{42} \text{ erg s}^{-1}$, but at $R_w = 0.1$ kpc. Therefore, the wind parameters considered in this study not only reflect the strength of the wind at a particular location (~1 kpc) for varying wind powers, but can also be mapped to different locations relative to the central AGN for a particular wind power.

It is essential to emphasize that Weaver et al. (1977) solution is only used to set the reference values of wind parameters, which are kept constant throughout the simulation. In reality, the radius of the bubble ($R_w$) increases with time as it expands. Therefore, the wind solutions are time-dependent. Ideally, one should consider a self-consistent evolution of the bubble with time for the wind injection. Additionally, the bubble solution of an AGN-driven wind solution consists of several distinct internal structures, i.e. the forward shock, shocked ambient medium, shocked wind and wind material (e.g. Weaver et al. 1977; King 2003; Faucher-Giguère & Quataert 2012;

Zubovas & King 2012; Costa, Sijacki & Haehnelt 2014), which are not generally easy to implement in local box simulations like ours. Moreover, as our main focus is to investigate the effect of a steady wind whose parameters are similar to AGN-driven outflows on kpc scales on star-forming complexes, we adopt the simplistic wind injection model. Importantly, the main purpose of this study is to investigate the impact of the wind on clouds with and without self-gravity. Thus, while details of the wind modelling are simplified, we can still make a meaningful comparison, as clouds in the simulation with and without self-gravity face the exact same wind.

## 2.5 Cooling

In order to account for the energy losses due to radiative cooling (see equation 3), we use the non-equilibrium cooling function calculated using MAPPINGS V code (Sutherland & Dopita 2017; Sutherland et al. 2018). This code utilizes a comprehensive database of atomic data to self-consistently compute the optically thin cooling rate for various gas phases, including cold neutral, warm neutral, partially ionized, and fully ionized gas in the temperature range $10^2-10^9$ K. For temperatures exceeding ~$10^9$ K, the cooling function is extended by assuming bremsstrahlung emission (e.g. Krause & Alexander 2007; Mukherjee et al. 2018). We also impose a lower temperature threshold of $T_{\text{floor}} = 100$ K, below which the cooling is turned off and the temperature of any cell falling below $T_{\text{floor}}$ is enforced to stay at $T_{\text{floor}}$. No additional heating terms, such as the galactic ultraviolet (UV) background or ionizing photons (UV, soft, and hard X-ray) from the AGN, are included. In our simulations, the gas that can cool below $10^4$ K is sufficiently dense ($n_c > 500 \text{ cm}^{-3}$). In this density range, the high-density cores can be safely assumed to be self-shielded from the external UV and soft X-ray photons (Rahmati et al. 2013; Meenakshi et al. 2022a). Conversely, due to the very small photoabsorption cross-section of hard X-ray photons ($E \gtrsim 20$ keV, primarily due to K-shell ionization of Fe and Ni ions) (Band et al. 1990; Verner et al. 1993), the cloudlets remain optically thin, resulting in negligible heating. Through explicit radiative-transfer calculations using CLOUDY (Ferland et al. 2017) with the AGN spectral energy distribution (SED) as used in Meenakshi et al.
(2022a), we confirm that ionization within the dense region ($n_c > 500 \text{ cm}^{-3}$) by an external radiation field (AGN and UV background) is insignificant, even for an AGN with a bolometric luminosity of $L_{bol} = 10^{45} \text{ erg s}^{-1}$, where radiation can only penetrate to a depth of $\lesssim 0.5 \text{ pc}$ from the illuminated surface. This depth is even smaller for lower luminosities or higher densities. Therefore, the impact of ionizing radiation can be safely disregarded.

We create a table of the cooling rate $\Lambda$ as a function of temperature in the range of $10^2-10^{10} \text{ K}$ assuming solar metallicities (Asplund et al. 2009). The values are then interpolated to every computational cell at runtime and are treated as a source term in equation (3). The equation is solved using a fractional step formalism, where the HD evolution and source step are solved separately through operator splitting. The energy losses from radiative cooling are computed by integrating the internal energy equation

$$ \frac{\partial(\rho e)}{\partial t} = - \left( \frac{\rho}{\mu m_{\text{H}}} \right)^2 \Lambda(T), \quad (15) $$

using an adaptively chosen, explicit or semi-implicit embedded Runge-Kutta method, depending on the 'stiffness' of the equation (see e.g. Teşileanu, Mignone & Massaglia 2008).

We note that we do not account for the explicit density dependence of the cooling function ($\Lambda$), which becomes significant in gas with $T \lesssim 10^4 \text{ K}$ and $n \gtrsim 10^4 \text{ cm}^{-3}$ in the case for non-equilibrium cooling. In our simulations, however, the initial temperature of gas with $n > 10^3 \text{ cm}^{-3}$ is below the floor temperature of $T < 100 \text{ K}$, where radiative cooling is turned off. Additionally, the shock-heated gas in our simulations rarely reaches the density and temperature ranges where the density dependence of the cooling curve would be significant. Therefore, the simplification of the cooling curve as a function of temperature only does not affect our results. Moreover, our simulations do not include molecular cooling, which might be a dominant mechanism below $T = 1000 \text{ K}$, and completely dominates the cooling process below 100 K (e.g. Goldsmith & Langer 1978). Therefore, for more accurate modelling of gas cooling and star formation, one must include the low-temperature processes of the cold molecular phase (Koyama & Inutsuka 2000; Vázquez-Semadeni et al. 2007; Glover et al. 2010). None the less, as we impose a temperature floor of 100 K, the temperature does not fall below this range. Hence, the absence of molecular cooling does not affect the result significantly.

## 2.6 Wind injection and boundary conditions

We approximate the AGN-driven spherically symmetric wind at a distance of 1 kpc from the central AGN as a planar wind on the scale of the cloud, propagating in the positive x-direction. The initial position of the forward shock is set at $x = -40 \text{ pc}$ and the domain within $-50 \le x \le -40 \text{ pc}$ is initialized with the wind state. The wind is constantly injected from the left boundary along the x-direction using an inflow boundary condition, where we populate the fluid quantities (pressure, velocity, and density) in the ghost cells with the wind properties as listed in Table 1.

However, one complication of subsonic inflows is that only two out of three characteristic waves enter the domain, while the third leaves. This implies that only two out of three primitive variables (ideally density-pressure or density-velocity pairs for well-poised conditions, e.g. see section 19.3 of Laney 1998) can be specified physically at the boundary, with the remaining variable set numerically from the interior solution (Thompson 1990; LeVeque 2002). Yet, a difficulty arises due to the finite size of the computational domain: the outgoing wave, which should freely exit the domain, experiences numerical

reflection at the boundary, thereby reducing the accuracy of the interior solution (Majda & Osher 1975). Therefore, the value of the third primitive variable should be chosen to allow the wave to exit the domain with minimal reflection, a condition commonly referred to as the non-reflecting boundary condition (NRBC), which itself is a broad area of research (Engquist & Majda 1977; Engquist & Majda 1979; Hedstrom 1979; Bayliss & Turkel 1980).

However, in our simulations, we do not include such a treatment to minimize boundary reflections. As a result, the subsonic winds injected at the ghost-cell layers, do not emerge with the same parameters (density, velocity, pressure, Mach, and power) as tabulated in Table 1. Fig. 1 depicts the time evolution of the wind density (top left panel), velocity (top right panel), pressure (bottom left panel), and power (bottom right panel) for the subsonic cases, as indicated in the legend. The values are calculated by taking the average over a y-z slice at $i = 1$, i.e. the first interior cell from the x-left boundary. We observe that all wind parameters experience an initial transient phase when the wind starts to progress through the stationary ambient medium, sweeping up the material. However, after $t \gtrsim 0.4 t_{cc}$, the parameters saturate to constant values, which are different from the injected values at the ghost zone. In Table 2, we tabulate the injected values of the wind parameters and their saturated values. None the less, the deviation of the parameters from the intended values is almost similar for all the powers. The wind velocities exhibit the most significant deviation, reaching $\lesssim 75 \text{ per cent}$ higher than the intended values. Conversely, pressure values decrease by $\lesssim 40 \text{ per cent}$, therefore increasing the Mach number of the winds (see Table 2), which are still subsonic. However, as the winds are dominated by thermal energy, the deviations of the wind power from the intended values are less, increased by $\lesssim 30 \text{ per cent}$.

Nevertheless, since our main objective in this study is to investigate the interaction between a cloud and winds possessing parameters akin to AGN-driven winds at specific powers, the qualitative outcomes among winds of varying powers, as well as the impact of self-gravity, remain unaffected. Moreover, the saturated values of the wind parameters, i.e. density, velocity, pressure, and power (see Table 2) fall well within the desired AGN wind parameter space, as laid out in Section 2.4. Furthermore, implementing non-reflecting boundary conditions at the inflow boundary does not ensure that the wind will emerge with the intended power, as one of the three primitive variables remains unconstrained and is set numerically. Indeed, as demonstrated in Appendix A, even though the physically specified primitive variables (density and velocity) remain constant, the pressure and power of the wind deviate significantly ($\sim 90 \text{ per cent}$) from the intended value, exceeding that of the simulation using the wind injection method employed in this study.

Except for the inflow boundary, all other boundaries of the computational domain are set to diode boundary conditions – the modification of the outflow boundary condition that prevents inflow into the domain so that gas can only leave the computational domain. We employ isolated boundary conditions for the gravitational potential ($\Phi$) on all sides of the domain, which is calculated using a multipole expansion of the density distribution up to order $l = 4$ (see appendix C of Mandal et al. 2021, for details).

## 2.7 Simulations and naming convention

In this study, we conduct a total of 13 three-dimensional (3D) simulations (including 10 simulations with self-gravity) covering a large parameter space of both the wind and cloud. An additional two simulations (with and without self-gravity) are performed without wind in order to ascertain the effects of winds. For all these
(2022a), we confirm that ionization within the dense region ($n_c > 500 \text{ cm}^{-3}$) by an external radiation field (AGN and UV background) is insignificant, even for an AGN with a bolometric luminosity of $L_{bol} = 10^{45} \text{ erg s}^{-1}$, where radiation can only penetrate to a depth of $\lesssim 0.5 \text{ pc}$ from the illuminated surface. This depth is even smaller for lower luminosities or higher densities. Therefore, the impact of ionizing radiation can be safely disregarded.

We create a table of the cooling rate $\Lambda$ as a function of temperature in the range of $10^2-10^{10} \text{ K}$ assuming solar metallicities (Asplund et al. 2009). The values are then interpolated to every computational cell at runtime and are treated as a source term in equation (3). The equation is solved using a fractional step formalism, where the HD evolution and source step are solved separately through operator splitting. The energy losses from radiative cooling are computed by integrating the internal energy equation

$$ \frac{\partial(\rho e)}{\partial t} = - \left( \frac{\rho}{\mu m_{\text{H}}} \right)^2 \Lambda(T), \quad (15) $$

using an adaptively chosen, explicit or semi-implicit embedded Runge-Kutta method, depending on the 'stiffness' of the equation (see e.g. Teşileanu, Mignone & Massaglia 2008).

We note that we do not account for the explicit density dependence of the cooling function ($\Lambda$), which becomes significant in gas with $T \lesssim 10^4 \text{ K}$ and $n \gtrsim 10^4 \text{ cm}^{-3}$ in the case for non-equilibrium cooling. In our simulations, however, the initial temperature of gas with $n > 10^3 \text{ cm}^{-3}$ is below the floor temperature of $T < 100 \text{ K}$, where radiative cooling is turned off. Additionally, the shock-heated gas in our simulations rarely reaches the density and temperature ranges where the density dependence of the cooling curve would be significant. Therefore, the simplification of the cooling curve as a function of temperature only does not affect our results. Moreover, our simulations do not include molecular cooling, which might be a dominant mechanism below $T = 1000 \text{ K}$, and completely dominates the cooling process below 100 K (e.g. Goldsmith & Langer 1978). Therefore, for more accurate modelling of gas cooling and star formation, one must include the low-temperature processes of the cold molecular phase (Koyama & Inutsuka 2000; Vázquez-Semadeni et al. 2007; Glover et al. 2010). None the less, as we impose a temperature floor of 100 K, the temperature does not fall below this range. Hence, the absence of molecular cooling does not affect the result significantly.

## 2.6 Wind injection and boundary conditions

We approximate the AGN-driven spherically symmetric wind at a distance of 1 kpc from the central AGN as a planar wind on the scale of the cloud, propagating in the positive x-direction. The initial position of the forward shock is set at $x = -40 \text{ pc}$ and the domain within $-50 \le x \le -40 \text{ pc}$ is initialized with the wind state. The wind is constantly injected from the left boundary along the x-direction using an inflow boundary condition, where we populate the fluid quantities (pressure, velocity, and density) in the ghost cells with the wind properties as listed in Table 1.

However, one complication of subsonic inflows is that only two out of three characteristic waves enter the domain, while the third leaves. This implies that only two out of three primitive variables (ideally density-pressure or density-velocity pairs for well-poised conditions, e.g. see section 19.3 of Laney 1998) can be specified physically at the boundary, with the remaining variable set numerically from the interior solution (Thompson 1990; LeVeque 2002). Yet, a difficulty arises due to the finite size of the computational domain: the outgoing wave, which should freely exit the domain, experiences numerical

reflection at the boundary, thereby reducing the accuracy of the interior solution (Majda & Osher 1975). Therefore, the value of the third primitive variable should be chosen to allow the wave to exit the domain with minimal reflection, a condition commonly referred to as the non-reflecting boundary condition (NRBC), which itself is a broad area of research (Engquist & Majda 1977; Engquist & Majda 1979; Hedstrom 1979; Bayliss & Turkel 1980).

However, in our simulations, we do not include such a treatment to minimize boundary reflections. As a result, the subsonic winds injected at the ghost-cell layers, do not emerge with the same parameters (density, velocity, pressure, Mach, and power) as tabulated in Table 1. Fig. 1 depicts the time evolution of the wind density (top left panel), velocity (top right panel), pressure (bottom left panel), and power (bottom right panel) for the subsonic cases, as indicated in the legend. The values are calculated by taking the average over a $y-z$ slice at $i = 1$, i.e. the first interior cell from the x-left boundary. We observe that all wind parameters experience an initial transient phase when the wind starts to progress through the stationary ambient medium, sweeping up the material. However, after $t \gtrsim 0.4 t_{cc}$, the parameters saturate to constant values, which are different from the injected values at the ghost zone. In Table 2, we tabulate the injected values of the wind parameters and their saturated values. None the less, the deviation of the parameters from the intended values is almost similar for all the powers. The wind velocities exhibit the most significant deviation, reaching $\lesssim 75 \text{ per cent}$ higher than the intended values. Conversely, pressure values decrease by $\lesssim 40 \text{ per cent}$, therefore increasing the Mach number of the winds (see Table 2), which are still subsonic. However, as the winds are dominated by thermal energy, the deviations of the wind power from the intended values are less, increased by $\lesssim 30 \text{ per cent}$.

Nevertheless, since our main objective in this study is to investigate the interaction between a cloud and winds possessing parameters akin to AGN-driven winds at specific powers, the qualitative outcomes among winds of varying powers, as well as the impact of self-gravity, remain unaffected. Moreover, the saturated values of the wind parameters, i.e. density, velocity, pressure, and power (see Table 2) fall well within the desired AGN wind parameter space, as laid out in Section 2.4. Furthermore, implementing non-reflecting boundary conditions at the inflow boundary does not ensure that the wind will emerge with the intended power, as one of the three primitive variables remains unconstrained and is set numerically. Indeed, as demonstrated in Appendix A, even though the physically specified primitive variables (density and velocity) remain constant, the pressure and power of the wind deviate significantly ($\sim 90 \text{ per cent}$) from the intended value, exceeding that of the simulation using the wind injection method employed in this study.

Except for the inflow boundary, all other boundaries of the computational domain are set to diode boundary conditions – the modification of the outflow boundary condition that prevents inflow into the domain so that gas can only leave the computational domain. We employ isolated boundary conditions for the gravitational potential ($\Phi$) on all sides of the domain, which is calculated using a multipole expansion of the density distribution up to order $l = 4$ (see appendix C of Mandal et al. 2021, for details).

## 2.7 Simulations and naming convention

In this study, we conduct a total of 13 three-dimensional (3D) simulations (including 10 simulations with self-gravity) covering a large parameter space of both the wind and cloud. An additional two simulations (with and without self-gravity) are performed without wind in order to ascertain the effects of winds. For all these
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Time evolution of the wind parameters, i.e. density (top left panel), velocity (top right panel), pressure (bottom left panel) and power (bottom right panel) in the simulations with different power as indicated in the legend. The bottom section of each panel shows the fractional deviation of the parameter from the initial values.</figcaption></figure>

<table><caption>**Table 2.** The saturated values of the wind parameters for the subsonic cases.</caption><thead><tr><th>P<sub>w, 0</sub><br>(erg s<sup>-1</sup>)</th><th>P<sub>w, sat</sub><br>(erg s<sup>-1</sup>)</th><th>n<sub>w, 0</sub><br>(cm<sup>-3</sup>)</th><th>n<sub>w, sat</sub><br>(cm<sup>-3</sup>)</th><th>v<sub>w, 0</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>w, sat</sub><br>(km s<sup>-1</sup>)</th><th>p<sub>w, 0</sub><br>(dyne cm<sup>-2</sup>)</th><th>p<sub>w, sat</sub><br>(dyne cm<sup>-2</sup>)</th><th>M<sub>w, 0</sub></th><th>M<sub>w, sat</sub></th><th>t<sub>cc, sat</sub><br>(Myr)</th></tr></thead><tbody><tr><td>10<sup>42</sup></td><td>1.09 × 10<sup>42</sup></td><td>10<sup>-2</sup></td><td>8.01 × 10<sup>-3</sup></td><td>400</td><td>623</td><td>1.32 × 10<sup>-10</sup></td><td>8.04 × 10<sup>-11</sup></td><td>0.286</td><td>0.621</td><td>6.20</td></tr><tr><td>10<sup>43</sup></td><td>1.28 × 10<sup>43</sup></td><td>10<sup>-2</sup></td><td>8.06 × 10<sup>-3</sup></td><td>1000</td><td>1481</td><td>5.25 × 10<sup>-10</sup></td><td>3.84 × 10<sup>-10</sup></td><td>0.329</td><td>0.678</td><td>2.60</td></tr><tr><td>10<sup>44</sup></td><td>1.10 × 10<sup>44</sup></td><td>10<sup>-2</sup></td><td>7.62 × 10<sup>-3</sup></td><td>1800</td><td>3075</td><td>2.95 × 10<sup>-9</sup></td><td>1.60 × 10<sup>-9</sup></td><td>0.277</td><td>0.671</td><td>1.29</td></tr><tr><td>10<sup>45</sup></td><td>1.09 × 10<sup>45</sup></td><td>10<sup>-2</sup></td><td>7.98 × 10<sup>-3</sup></td><td>4000</td><td>6281</td><td>1.32 × 10<sup>-8</sup></td><td>8.02 × 10<sup>-9</sup></td><td>0.285</td><td>0.625</td><td>0.62</td></tr></tbody></table>

*Notes.* The subscripts with 0 refer to the intended values of the wind parameters as tabulated in Table 1, while ‘sat’ subscripts are the saturated values of the parameters after $\ge 0.4t_{cc}$.

simulations, we initialize the cloud with an initial virial parameter, $\alpha_{vir, 0} = 0.9$, and using this value, we determine the initial mean cloud density and velocity dispersion.

In our fiducial simulations, the mean density of the cloud ($\bar{n}_c$) is set to $200 \text{ cm}^{-3}$ with $k_{\min} = 3$ and an initial velocity dispersion of $8 \text{ km s}^{-1}$ to achieve $\alpha_{vir, 0} = 0.9$. Therefore, the mass of the standard cloud in our simulation is $M_c = 2.873 \times 10^5 \text{ M}_\odot$, which is typical of the average mass of GMCs in the Milky Way and nearby galaxies (Hughes et al. 2013; Miville-Deschênes et al. 2017). To investigate the effect of a similar power wind on a cloud with lower density, we consider one simulation with a lower mean density of $20 \text{ cm}^{-3}$ and velocity dispersion of $2.5 \text{ km s}^{-1}$. The fiducial wind in our simulation

is thermal energy dominated. The wind parameters and the minimum wavelength of the cloud are varied as discussed in Sections 2.3 and 2.4. The initial conditions and the parameters along with the name for all the simulations are listed in Table 1

It is useful to clarify the naming convention of the simulations. For example, in GC43_k3, ‘G’ implies self-gravity is present and ‘C’ stands for cooling, which is common in all the simulations. The number ‘43’ represents the total power of the wind of $10^{43} \text{ erg s}^{-1}$. the label ‘k3’ represents the minimum wavenumber of the cloud; in this case $k_{\min} = 3$. We explicitly added the label ‘_low’ to ‘GC43_k3’ for the case of the lower mean density cloud with wind power of $10^{43} \text{ erg s}^{-1}$. Similarly, ‘_kinetic’ is added to ‘GC45_k3’ to indicate
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Time evolution of the wind parameters, i.e. density (top left panel), velocity (top right panel), pressure (bottom left panel) and power (bottom right panel) in the simulations with different power as indicated in the legend. The bottom section of each panel shows the fractional deviation of the parameter from the initial values.</figcaption></figure>

<table><caption>**Table 2.** The saturated values of the wind parameters for the subsonic cases.</caption><thead><tr><th>P<sub>w, 0</sub><br>(erg s<sup>-1</sup>)</th><th>P<sub>w, sat</sub><br>(erg s<sup>-1</sup>)</th><th>n<sub>w, 0</sub><br>(cm<sup>-3</sup>)</th><th>n<sub>w, sat</sub><br>(cm<sup>-3</sup>)</th><th>v<sub>w, 0</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>w, sat</sub><br>(km s<sup>-1</sup>)</th><th>p<sub>w, 0</sub><br>(dyne cm<sup>-2</sup>)</th><th>p<sub>w, sat</sub><br>(dyne cm<sup>-2</sup>)</th><th>M<sub>w, 0</sub></th><th>M<sub>w, sat</sub></th><th>t<sub>cc, sat</sub><br>(Myr)</th></tr></thead><tbody><tr><td>10<sup>42</sup></td><td>1.09 × 10<sup>42</sup></td><td>10<sup>-2</sup></td><td>8.01 × 10<sup>-3</sup></td><td>400</td><td>623</td><td>1.32 × 10<sup>-10</sup></td><td>8.04 × 10<sup>-11</sup></td><td>0.286</td><td>0.621</td><td>6.20</td></tr><tr><td>10<sup>43</sup></td><td>1.28 × 10<sup>43</sup></td><td>10<sup>-2</sup></td><td>8.06 × 10<sup>-3</sup></td><td>1000</td><td>1481</td><td>5.25 × 10<sup>-10</sup></td><td>3.84 × 10<sup>-10</sup></td><td>0.329</td><td>0.678</td><td>2.60</td></tr><tr><td>10<sup>44</sup></td><td>1.10 × 10<sup>44</sup></td><td>10<sup>-2</sup></td><td>7.62 × 10<sup>-3</sup></td><td>1800</td><td>3075</td><td>2.95 × 10<sup>-9</sup></td><td>1.60 × 10<sup>-9</sup></td><td>0.277</td><td>0.671</td><td>1.29</td></tr><tr><td>10<sup>45</sup></td><td>1.09 × 10<sup>45</sup></td><td>10<sup>-2</sup></td><td>7.98 × 10<sup>-3</sup></td><td>4000</td><td>6281</td><td>1.32 × 10<sup>-8</sup></td><td>8.02 × 10<sup>-9</sup></td><td>0.285</td><td>0.625</td><td>0.62</td></tr></tbody></table>

*Notes.* The subscripts with 0 refer to the intended values of the wind parameters as tabulated in Table 1, while ‘sat’ subscripts are the saturated values of the parameters after $\ge 0.4t_{cc}$.

simulations, we initialize the cloud with an initial virial parameter, $\alpha_{vir, 0} = 0.9$, and using this value, we determine the initial mean cloud density and velocity dispersion.

In our fiducial simulations, the mean density of the cloud ($\bar{n}_c$) is set to $200 \text{ cm}^{-3}$ with $k_{\min} = 3$ and an initial velocity dispersion of $8 \text{ km s}^{-1}$ to achieve $\alpha_{vir, 0} = 0.9$. Therefore, the mass of the standard cloud in our simulation is $M_c = 2.873 \times 10^5 \text{ M}_\odot$, which is typical of the average mass of GMCs in the Milky Way and nearby galaxies (Hughes et al. 2013; Miville-Deschênes et al. 2017). To investigate the effect of a similar power wind on a cloud with lower density, we consider one simulation with a lower mean density of $20 \text{ cm}^{-3}$ and velocity dispersion of $2.5 \text{ km s}^{-1}$. The fiducial wind in our simulation

is thermal energy dominated. The wind parameters and the minimum wavelength of the cloud are varied as discussed in Sections 2.3 and 2.4. The initial conditions and the parameters along with the name for all the simulations are listed in Table 1

It is useful to clarify the naming convention of the simulations. For example, in GC43_k3, ‘G’ implies self-gravity is present and ‘C’ stands for cooling, which is common in all the simulations. The number ‘43’ represents the total power of the wind of $10^{43} \text{ erg s}^{-1}$. the label ‘k3’ represents the minimum wavenumber of the cloud; in this case $k_{\min} = 3$. We explicitly added the label ‘_low’ to ‘GC43_k3’ for the case of the lower mean density cloud with wind power of $10^{43} \text{ erg s}^{-1}$. Similarly, ‘_kinetic’ is added to ‘GC45_k3’ to indicate
that the wind of power $10^{45}$ erg s$^{-1}$ is kinetic energy dominated. If not specified, the default minimum wavenumber of the cloud is $k_{min} = 3$.

## 2.8 Relevant time-scales

There are various dynamical time-scales involved in the problem that we consider in this study and the absolute values for different simulations are tabulated in Table 1.

(i) The shock-passing time ($t_{sp}$), i.e. the approximate time the initial shock takes to sweep over the whole cloud (Klein et al. 1994; Fragile et al. 2004):

$$t_{sp} \approx \frac{2R_c}{v_w}, \qquad (16)$$

where $R_c$ is the radius of the cloud and $v_w$ is the velocity of the wind.

(ii) The cloud-crushing time ($t_{cc}$), the typical time in which the shock will compress the cloud (Fragile et al. 2004):

$$t_{cc} \approx \frac{R_c}{v_{st}} \approx \chi^{1/2} \frac{R_c}{v_w}, \qquad (17)$$

where $v_{st}$ is the velocity of the transmitted shock into the cloud and $\chi = \rho_c/\rho_w$ is the density contrast between the cloud and wind.

(iii) The freefall time-scale:

$$t_{ff} = \sqrt{\frac{3\pi}{32G\rho_c}}, \qquad (18)$$

(iv) The cooling time-scale (Klein et al. 1994; Mellema et al. 2002; Fragile et al. 2004):

$$t_{\text{cool}} = \frac{3k_B n_c \langle T_c \rangle}{2n_c^2 \Lambda(T)}, \qquad (19)$$

where $k_B$ is the Boltzmann constant, $n_c$ and $\langle T_c \rangle$, the average initial number density and temperature of the cloud. For the fiducial simulations with an initial mean number density of $n_c = 200 \text{ cm}^{-3}$ and mean temperature of $\langle T_c \rangle \approx 3 \times 10^3 \text{ K}$, as estimated from the pressure equilibrium condition with the ambient medium, the typical cooling time-scale is of the order of $3 \times 10^2 \text{ yr}$, which is significantly shorter than relevant other time-scales (see Table 1). However, it is essential to recognize that this estimate relies on average temperature and density values, whereas the fractal nature of the clouds in the simulations introduces a wide range of temperatures due to the inhomogeneous density distribution. Therefore, the cooling time-scale can exhibit significant variations across different regions within the cloud.

(v) The drag time-scale ($t_{drag}$), by which the cloud is accelerated to a similar velocity as that of the wind (Klein et al. 1994; Fragile et al. 2004):

$$t_{\text{drag}} \approx \chi^{1/2} t_{cc} \qquad (20)$$

(vi) The time ($t_{KH}$) for the Kelvin–Helmholtz (KH) instability to grow for $\chi \gg 1$ (Chandrasekhar 1961):

$$t_{KH} \sim \frac{\chi^{1/2}}{k_{KH} v_{rel}} \approx \frac{4 r_{\text{cloudlet, kmin}} \chi^{1/2}}{3 v_w}, \qquad (21)$$

where $k_{KH}$ and $v_{rel}$ are the wavenumber of the KH perturbations and the relative velocity between the post-shock background and cloud, respectively. For $\chi \gg 1$, the relative velocity is approximately equal to the post-shock velocity of the background, i.e. $v_{rel} \approx 3v_w/4$. We set $k_{KH} \sim 1/r_{\text{cloudlet, kmin}}$ because even though the instabilities at the smallest wavelengths grow more rapidly, the most detrimental

<figure><img src="image_3.png" /><figcaption>Figure 2. Number density slice through the $z=0$ plane for simulations with a fractal with $k_{min} = 3$ (GC43_k3, top panel) and uniform (GC43_uniform, bottom panel) cloud at 0.63 Myr. The wind power in both the simulation is $10^{43}$ erg s$^{-1}$. The wind readily penetrates the fractal cloud through low-density intercloudlet channels (top panel), leading to the formation of numerous small and dense cloudlets. On the other hand, in the simulation featuring a uniform spherical cloud (bottom panel), the initial compression forms a highly dense shell, which cools efficiently. This shell slows down the transmitted shock into the cloud, significantly diminishing the wind's impact on the overall cloud.</figcaption></figure>

wavelengths are those approximately equal to the cloudlet radius (e.g. see Klein et al. 1994; Poludnenko, Frank & Blackman 2002).

# 3 RESULTS

## 3.1 Morphological evolution of the cloud

### 3.1.1 General features

The evolution of a cloud impacted by a highly pressurized wind can be divided into three phases. Initially, as the wind comes into contact with the cloud, it initiates an internal shock propagating from the cloud's surface toward its centre. The force exerted by the wind's ram pressure causes compression within the cloud, leading to a rise in the average density of the cloud material. The lognormal density distribution of the clouds results in many high-density fragments (cloudlets) and low-density channels inside the fractal clouds. This aids the wind to propagate deeper into the clouds, resulting in a stronger interaction with the cloud and deeper penetration when compared to cases with a uniform spherical cloud, as demonstrated in Fig. 2 (also see Cooper et al. 2009; Schneider & Robertson 2017; Banda-Barragán et al. 2018).

Subsequently, as the wind continues to flow downstream wrapping around the cloud, a shear layer forms at the wind–cloud interface, leading to the onset of Kelvin–Helmholtz (KH) instability. As a result, the outer layer of the cloud gets stripped, mixed with the wind, and funnelled downstream where it condensates (middle left panel of Fig. 3), giving rise to a series of cold and dense cloudlets that form a trailing tail behind the main cloud. During this phase, if no cooling mechanism is present, the energy transport from the wind to the cloud in the form of thermal energy should increase the temperature of the cloud material, which would result in an expansion
that the wind of power $10^{45}$ erg s$^{-1}$ is kinetic energy dominated. If not specified, the default minimum wavenumber of the cloud is $k_{min} = 3$.

## 2.8 Relevant time-scales

There are various dynamical time-scales involved in the problem that we consider in this study and the absolute values for different simulations are tabulated in Table 1.

(i) The shock-passing time ($t_{sp}$), i.e. the approximate time the initial shock takes to sweep over the whole cloud (Klein et al. 1994; Fragile et al. 2004):

$$t_{sp} \approx \frac{2R_c}{v_w}, \qquad (16)$$

where $R_c$ is the radius of the cloud and $v_w$ is the velocity of the wind.

(ii) The cloud-crushing time ($t_{cc}$), the typical time in which the shock will compress the cloud (Fragile et al. 2004):

$$t_{cc} \approx \frac{R_c}{v_{st}} \chi^{1/2} \frac{R_c}{v_w}, \qquad (17)$$

where $v_{st}$ is the velocity of the transmitted shock into the cloud and $\chi = \rho_c/\rho_w$ is the density contrast between the cloud and wind.

(iii) The freefall time-scale:

$$t_{ff} = \sqrt{\frac{3\pi}{32G\rho_c}}, \qquad (18)$$

(iv) The cooling time-scale (Klein et al. 1994; Mellema et al. 2002; Fragile et al. 2004):

$$t_{cool} = \frac{3k_B n_c \langle T_c \rangle}{2n_c^2 \Lambda(T)}, \qquad (19)$$

where $k_B$ is the Boltzmann constant, $n_c$ and $\langle T_c \rangle$, the average initial number density and temperature of the cloud. For the fiducial simulations with an initial mean number density of $n_c = 200$ cm$^{-3}$ and mean temperature of $\langle T_c \rangle \approx 3 \times 10^3$ K, as estimated from the pressure equilibrium condition with the ambient medium, the typical cooling time-scale is of the order of $3 \times 10^2$ yr, which is significantly shorter than relevant other time-scales (see Table 1). However, it is essential to recognize that this estimate relies on average temperature and density values, whereas the fractal nature of the clouds in the simulations introduces a wide range of temperatures due to the inhomogeneous density distribution. Therefore, the cooling time-scale can exhibit significant variations across different regions within the cloud.

(v) The drag time-scale ($t_{drag}$), by which the cloud is accelerated to a similar velocity as that of the wind (Klein et al. 1994; Fragile et al. 2004):

$$t_{drag} \approx \chi^{1/2} t_{cc} \qquad (20)$$

(vi) The time ($t_{KH}$) for the Kelvin–Helmholtz (KH) instability to grow for $\chi \gg 1$ (Chandrasekhar 1961):

$$t_{KH} \sim \frac{\chi^{1/2}}{k_{KH} v_{rel}} \approx \frac{4 r_{cloudlet, k_{min}} \chi^{1/2}}{3 v_w}, \qquad (21)$$

where $k_{KH}$ and $v_{rel}$ are the wavenumber of the KH perturbations and the relative velocity between the post-shock background and cloud, respectively. For $\chi \gg 1$, the relative velocity is approximately equal to the post-shock velocity of the background, i.e. $v_{rel} \approx 3v_w/4$. We set $k_{KH} \sim 1/r_{cloudlet, k_{min}}$ because even though the instabilities at the smallest wavelengths grow more rapidly, the most detrimental

<figure><img src="image_3.png" /><figcaption>Figure 2. Number density slice through the $z=0$ plane for simulations with a fractal with $k_{min} = 3$ (GC43_k3, top panel) and uniform (GC43_uniform, bottom panel) cloud at 0.63 Myr. The wind power in both the simulation is $10^{43}$ erg s$^{-1}$. The wind readily penetrates the fractal cloud through low-density intercloudlet channels (top panel), leading to the formation of numerous small and dense cloudlets. On the other hand, in the simulation featuring a uniform spherical cloud (bottom panel), the initial compression forms a highly dense shell, which cools efficiently. This shell slows down the transmitted shock into the cloud, significantly diminishing the wind’s impact on the overall cloud.</figcaption></figure>

wavelengths are those approximately equal to the cloudlet radius (e.g. see Klein et al. 1994; Poludnenko, Frank & Blackman 2002).

# 3 RESULTS

## 3.1 Morphological evolution of the cloud

### 3.1.1 General features

The evolution of a cloud impacted by a highly pressurized wind can be divided into three phases. Initially, as the wind comes into contact with the cloud, it initiates an internal shock propagating from the cloud’s surface toward its centre. The force exerted by the wind’s ram pressure causes compression within the cloud, leading to a rise in the average density of the cloud material. The lognormal density distribution of the clouds results in many high-density fragments (cloudlets) and low-density channels inside the fractal clouds. This aids the wind to propagate deeper into the clouds, resulting in a stronger interaction with the cloud and deeper penetration when compared to cases with a uniform spherical cloud, as demonstrated in Fig. 2 (also see Cooper et al. 2009; Schneider & Robertson 2017; Banda-Barragán et al. 2018).

Subsequently, as the wind continues to flow downstream wrapping around the cloud, a shear layer forms at the wind–cloud interface, leading to the onset of Kelvin–Helmholtz (KH) instability. As a result, the outer layer of the cloud gets stripped, mixed with the wind, and funnelled downstream where it condensates (middle left panel of Fig. 3), giving rise to a series of cold and dense cloudlets that form a trailing tail behind the main cloud. During this phase, if no cooling mechanism is present, the energy transport from the wind to the cloud in the form of thermal energy should increase the temperature of the cloud material, which would result in an expansion
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Time evolution (row-wise) of the projected number density along the z-axis of the simulations without self-gravity (left panel) and with self-gravity (right panel). The power of the wind for both simulations is $10^{42}$ erg s$^{-1}$. It is evident that in the absence of self-gravity (left panel), the cloud undergoes expansion after the initial compression, whereas in the self-gravity case (right panel) global collapse takes place.</figcaption></figure>

of the cloud, making the cloud more prone to destruction (Orlando et al. 2005; Cooper et al. 2009). However, with the presence of radiative cooling, the compressed cloudlets cool efficiently, leading to a cessation of internal pressure support, contraction, and the formation of denser cloudlets. This contraction enhances the density contrast ($\chi$) between the cloudlets and the surrounding hot shear flow, consequently reducing the growth rate of KH instabilities (equation 21). Therefore, in the presence of radiative cooling, the cloudlets are relatively protected from instabilities induced by the shear flow compared to a non-radiative scenario (Mellema et al. 2002; Cooper et al. 2009; Scannapieco & Brüggen 2015; Banda-Barragán et al. 2021).

While the presence of radiative cooling reduces the growth rate of instabilities, it does not provide complete protection. Therefore, the wind-induced shear flow along the cloudlet surface still causes the entrainment of cloud material, albeit to a lesser extent compared to a non-radiative scenario (Cooper et al. 2009; Banda-Barragán et al. 2021). Additionally, the direct momentum transfer from the wind accelerates the cloudlets. Combined with the shredded cloudlets

from the main cloud due to the KH instability, these fragments persistently contribute to the elongation of the cloud, eventually leaving the computational domain (bottom left panel of Fig. 3).

Although the general evolutionary stages of a cloud impacted by a highly pressurized wind are similar for all cases, the corresponding time-scales and the overall evolution can depend on several factors which are discussed in subsequent sections.

### 3.1.2 Gravity versus no-gravity

Here, we examine the effect of self-gravity on the morphological evolution of the cloud by comparing it with a simulation that shares identical initial conditions and wind strength but lacks self-gravity. In Fig. 3, we show the column density map in the x-y plane of the simulations both without (C42.k3; left panels) and with (GC42.k3; right panels) self-gravity at three different times. The power of the wind is $10^{42}$ erg s$^{-1}$.

During the initial compression phase, the mean density of the cloud in both cases (with and without self-gravity) increases. After
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Time evolution (row-wise) of the projected number density along the z-axis of the simulations without self-gravity (left panel) and with self-gravity (right panel). The power of the wind for both simulations is $10^{42}$ erg s$^{-1}$. It is evident that in the absence of self-gravity (left panel), the cloud undergoes expansion after the initial compression, whereas in the self-gravity case (right panel) global collapse takes place.</figcaption></figure>

of the cloud, making the cloud more prone to destruction (Orlando et al. 2005; Cooper et al. 2009). However, with the presence of radiative cooling, the compressed cloudlets cool efficiently, leading to a cessation of internal pressure support, contraction, and the formation of denser cloudlets. This contraction enhances the density contrast ($\chi$) between the cloudlets and the surrounding hot shear flow, consequently reducing the growth rate of KH instabilities (equation 21). Therefore, in the presence of radiative cooling, the cloudlets are relatively protected from instabilities induced by the shear flow compared to a non-radiative scenario (Mellema et al. 2002; Cooper et al. 2009; Scannapieco & Brüggen 2015; Banda-Barragán et al. 2021).

While the presence of radiative cooling reduces the growth rate of instabilities, it does not provide complete protection. Therefore, the wind-induced shear flow along the cloudlet surface still causes the entrainment of cloud material, albeit to a lesser extent compared to a non-radiative scenario (Cooper et al. 2009; Banda-Barragán et al. 2021). Additionally, the direct momentum transfer from the wind accelerates the cloudlets. Combined with the shredded cloudlets

from the main cloud due to the KH instability, these fragments persistently contribute to the elongation of the cloud, eventually leaving the computational domain (bottom left panel of Fig. 3).

Although the general evolutionary stages of a cloud impacted by a highly pressurized wind are similar for all cases, the corresponding time-scales and the overall evolution can depend on several factors which are discussed in subsequent sections.

### 3.1.2 Gravity versus no-gravity

Here, we examine the effect of self-gravity on the morphological evolution of the cloud by comparing it with a simulation that shares identical initial conditions and wind strength but lacks self-gravity. In Fig. 3, we show the column density map in the x-y plane of the simulations both without (C42.k3; left panels) and with (GC42.k3; right panels) self-gravity at three different times. The power of the wind is $10^{42}$ erg s$^{-1}$.

During the initial compression phase, the mean density of the cloud in both cases (with and without self-gravity) increases. After
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Projected number density along the z-direction at $t = 0.85 t_{cc}$ of simulations with self-gravity and different wind power: $10^{42}$ (top left panel), $10^{43}$ (top right panel), $10^{44}$ (bottom left panel) and $10^{45}$ erg s$^{-1}$ (bottom right panel). In the unit of the Kelvin–Helmholtz time-scale, this time corresponds to $t = 1.51 t_{KH}$. The corresponding physical time for each simulation is shown in the bottom left corner of each panel. Interesting to note that, as wind power increases, the cloud becomes more disrupted. While in low-power cases, the cloud experiences global collapse, the higher power winds cause localized collapses of the dense cloudlets, simultaneously still disrupting the cloud.</figcaption></figure>

the initial shock-induced fragmentation phase, in the absence of self-gravity, the pressure-confined fragments undergo the highest attainable compression due to the transmitted shock into the cloudlets (middle left panel of Fig. 3). Once this compression limit is reached, the density of the clumps cannot increase any further. Hence, the cloud as a whole commences disintegration, resulting in cloud expansion (bottom left panel).

We point out that the expansion is not a result of adiabatic heating in the absence of radiative cooling, as outlined in various studies (Banda-Barragán et al. 2018, 2019). Instead, the expansion occurs because the low-density channels within the fractal cloud facilitate the infiltration of high-velocity wind material, which induces turbulence and vorticity within the cloud (Klein et al. 1994), leading to the transfer of momentum from the wind to the cloud material. This momentum transfer is what ultimately causes the cloud to expand.

In contrast, when self-gravity is present, the shock-compressed cloudlets become gravitationally bound and start to accrete material from the surroundings (middle right panel of Fig. 3), thereby attaining significantly higher density and mass compared to the case without self-gravity. Additionally, when acting in an inhomogeneous medium, self-gravity causes additional local fragmentation of the cloud, which creates many low-density channels for the wind to propagate inside the cloud, resulting in a stronger interaction. The increased average cloud density produced by compression deepens the gravitational potential, resulting in a more compact and tightly bound cloud structure. As a result, the interplay between the gravitational pull and the strength of the wind–cloud interaction becomes

the pivotal factor in deciding whether the cloud will experience a runaway collapse or will be disintegrated by the force of the wind.

For a low-power wind, for example, $10^{42}$ erg s$^{-1}$ as shown in Fig. 3, the momentum transfer from the wind to the cloud is significantly lower and therefore less gas is pushed out of the potential well of the cloud (which becomes deeper with time as more gas is accreted) compared to a high-power wind. Therefore, the wind triggers a runaway collapse of the cloud in the presence of self-gravity as can be seen from the bottom right panel of Fig. 3. In contrast, in the simulation without self-gravity, the cloud expands (bottom left panel of Fig. 3) and eventually will be elongated and destroyed by the wind.

### 3.1.3 Effect of varying wind power

Here, we show the morphological differences of the cloud in simulations in the presence of self-gravity with different wind velocities and pressure, which are the proxy for the power of the wind driven by AGN. Fig. 4 illustrates the column density map on the x–y plane of simulations with wind power $10^{42}$ (GC42_k3), $10^{43}$ (GC43_k3), $10^{44}$ (GC44_k3), and $10^{45}$ erg s$^{-1}$ (GC45_k3) at $t = 0.85 t_{cc}$. The corresponding physical time is also shown in each panel.

For the lower power cases (GC42_k3 and GC43_k3), the initial freefall time is shorter than the growth time-scale for the KH instability due to the lower wind velocity ($t_{KH} \propto 1/v_w$), which can further be reduced by the dramatic increase of density through compression. In this scenario, the cloud as a whole undergoes
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Projected number density along the z-direction at $t = 0.85 t_{cc}$ of simulations with self-gravity and different wind power: $10^{42}$ (top left panel), $10^{43}$ (top right panel), $10^{44}$ (bottom left panel) and $10^{45}$ erg s$^{-1}$ (bottom right panel). In the unit of the Kelvin–Helmholtz time-scale, this time corresponds to $t = 1.51 t_{KH}$. The corresponding physical time for each simulation is shown in the bottom left corner of each panel. Interesting to note that, as wind power increases, the cloud becomes more disrupted. While in low-power cases, the cloud experiences global collapse, the higher power winds cause localized collapses of the dense cloudlets, simultaneously still disrupting the cloud.</figcaption></figure>

the initial shock-induced fragmentation phase, in the absence of self-gravity, the pressure-confined fragments undergo the highest attainable compression due to the transmitted shock into the cloudlets (middle left panel of Fig. 3). Once this compression limit is reached, the density of the clumps cannot increase any further. Hence, the cloud as a whole commences disintegration, resulting in cloud expansion (bottom left panel).

We point out that the expansion is not a result of adiabatic heating in the absence of radiative cooling, as outlined in various studies (Banda-Barragán et al. 2018, 2019). Instead, the expansion occurs because the low-density channels within the fractal cloud facilitate the infiltration of high-velocity wind material, which induces turbulence and vorticity within the cloud (Klein et al. 1994), leading to the transfer of momentum from the wind to the cloud material. This momentum transfer is what ultimately causes the cloud to expand.

In contrast, when self-gravity is present, the shock-compressed cloudlets become gravitationally bound and start to accrete material from the surroundings (middle right panel of Fig. 3), thereby attaining significantly higher density and mass compared to the case without self-gravity. Additionally, when acting in an inhomogeneous medium, self-gravity causes additional local fragmentation of the cloud, which creates many low-density channels for the wind to propagate inside the cloud, resulting in a stronger interaction. The increased average cloud density produced by compression deepens the gravitational potential, resulting in a more compact and tightly bound cloud structure. As a result, the interplay between the gravitational pull and the strength of the wind–cloud interaction becomes

the pivotal factor in deciding whether the cloud will experience a runaway collapse or will be disintegrated by the force of the wind.

For a low-power wind, for example, $10^{42}$ erg s$^{-1}$ as shown in Fig. 3, the momentum transfer from the wind to the cloud is significantly lower and therefore less gas is pushed out of the potential well of the cloud (which becomes deeper with time as more gas is accreted) compared to a high-power wind. Therefore, the wind triggers a runaway collapse of the cloud in the presence of self-gravity as can be seen from the bottom right panel of Fig. 3. In contrast, in the simulation without self-gravity, the cloud expands (bottom left panel of Fig. 3) and eventually will be elongated and destroyed by the wind.

### 3.1.3 Effect of varying wind power

Here, we show the morphological differences of the cloud in simulations in the presence of self-gravity with different wind velocities and pressure, which are the proxy for the power of the wind driven by AGN. Fig. 4 illustrates the column density map on the x–y plane of simulations with wind power $10^{42}$ (GC42_k3), $10^{43}$ (GC43_k3), $10^{44}$ (GC44_k3), and $10^{45}$ erg s$^{-1}$ (GC45_k3) at $t = 0.85 t_{cc}$. The corresponding physical time is also shown in each panel.

For the lower power cases (GC42_k3 and GC43_k3), the initial freefall time is shorter than the growth time-scale for the KH instability due to the lower wind velocity ($t_{KH} \propto 1/v_w$), which can further be reduced by the dramatic increase of density through compression. In this scenario, the cloud as a whole undergoes
runaway gravitational collapse before any instability has a chance to act. This is clearly visible in the top left panel, where the cloud in GC42.k3 has collapsed by this time, becoming very compact and highly dense. None the less, due to the continuous ablation, the loosely bound outer layer of the cloud has been gradually stripped away by the wind, giving rise to a tail in the direction of the wind. In GC43.k3 (top right panel), the cloud is undergoing collapse but at a somewhat lower rate than in the $10^{42}$ erg s⁻¹ case due to the stronger wind and comparatively quicker growth of shear instabilities (such as KH instability). Thus, even though the cloud's central region is experiencing collapse, the wind's influence extends to a larger section of the cloud, resulting in a more extended structure compared to the $10^{42}$ erg s⁻¹ case.

The cloud disruption time-scales in GC44.k3 (bottom left panel) and GC45.k3 (bottom right panel) are much shorter than the free-fall time due to increased velocity and ram pressure (see Table 1). Therefore, a high-power wind induces much stronger turbulence and percolation of gas inside the inhomogeneous cloud, preventing the runaway collapse globally, compared to a low-power wind. Thus, in this scenario, the cloud is ablated and disrupted significantly before the influence of gravity becomes important. However, even in the case of a strong wind, the small cloudlets formed by the fragmentation become gravitationally bound and shielded from the external wind, by forming a high-density post-shock outer layer both due to compression and radiative cooling, preventing the wind material from infiltrating these cloudlets. Therefore, in the absence of any opposing mechanism, very high-density cloudlets collapse locally, while the wind ablates comparatively low-density material, shaping the cloud into an elongated and extended structure of cloudlets.

However, it is important to note that the density contrast between the cloud and wind material in these simulations is very high ($\chi \sim 2 \times 10^4$), which results in a very high value of the drag time-scale ($t_{drag} \sim \sqrt{\chi} t_{cc}$), which is the time by which the cloudlets are expected to attain a similar velocity as the wind. Therefore, even in scenarios of high power, where the cloud's overall collapse is impeded by the wind, a considerable portion of the fragmented cloudlets persists and remains gravitationally bound to the central potential well due to the inability to attain escape velocity, owing to the large drag time-scale.

### 3.1.4 Dependence on cloud fractal wavenumber

The maximum size ($\lambda_{max} \approx L/2\pi k_{min}$) of individual cloudlets inside the cloud is parameterized by the minimum wavenumber $k_{min}$. Therefore, a higher value of $k_{min}$ results in a larger number of small cloudlets as well as low-density intercloudlet channels. As discussed earlier, the interaction between the wind and a fractal cloud is influenced by the extent to which the wind material can permeate the medium between individual clumps separated by low-density cloud material, in addition to its impact on the fractal surfaces of the cloud. Thus, one can anticipate that in the case of a higher value of $k_{min}$ with a greater abundance of low-density channels, it is easier for the wind to percolate through these channels into the cloudlets and mix with cloud material, thereby transferring energy and momentum. Additionally, due to the reduced size of the cloudlets in this case, KH instabilities grow faster ($t_{KH} \propto 1/k$) in a cloud with lower $k_{min}$. Consequently, it is expected that a cloud with a higher $k_{min}$ value would undergo more pronounced disruption.

Hence, in order to investigate the impact of the cloud fractal wavenumber, we simulate wind-cloud interaction for clouds with

varying values of $k_{min}$ (=1, 3, 6, 10), while keeping the mass and mean density of the clouds fixed. The simulations are initiated with self-gravity and a wind of power $10^{43}$ erg s⁻¹. In Fig. 5, we show the number density slice<fnref n="4" /> through the $z=0$ plane for these four cases (row-wise) at 0 Myr (left panel) and 1.27 Myr (right panel).

If we first consider the $k_{min} = 1$ case ('k1', first row), the cloud contains approximately one big cloudlet with increasing density towards its centre, surrounded by low-density outer layers, which are rapidly eroded by the wind through ablation. Nevertheless, owing to its larger size and the added shielding from compression and cooling, the shock is unable to penetrate deep into the clump, at least not at early times compared to the other $k_{min}$ cases. Furthermore, external overpressurization triggers the gravitational collapse of the clump to form a compact, gravitationally bound structure.

In the $k_{min} = 3$ case ('k3', second row), the cloud is initially composed of many small, dense clumps separated by wide interclump channels. These channels help the high-velocity wind material to percolate into the cloud, transferring momentum. This results in more fragmentation of the cloud, preventing the global collapse, unlike in the $k_{min} = 1$ case.

Clouds with smaller scale density perturbations as in the $k_{min} = 6$ ('k6') and 10 ('k10') cases (third and fourth row), the numbers of cloudlets are much higher, and the intercloud channels are also consequently narrower. During the initial interaction, the swept-up material by the wind creates a dense layer near the interaction area, blocking the entrance of these narrow channels. This prevents the wind from dispersing the cloudlets efficiently, resulting in the accumulation near the centre. Thus, external overpressurization triggers the global collapse of the cloud in these cases similar to the 'k1' case. However, a striking difference between the 'k1' case and the $k_{min} = 6, 10$ cases is evident. In the former, cloudlets form through fragmentation and stripping of the large clump. In contrast, in the 'k6' and 'k10' cases, cloudlets were initially seeded by the fractal generator itself at the initialization of the simulation and they remain organized in a more spherically symmetric manner at later times.

### 3.1.5 Dependence on the mean cloud density

In order to understand the influence of the density contrast between the cloud and wind, we consider two simulations including self-gravity with identical wind power of $10^{43}$ erg s⁻¹, but two different values of the initial mean cloud density, namely, $\bar{n}_c = 200 \text{ cm}^{-3}$ ($\chi = 2 \times 10^4$, GC43.k3) and $20 \text{ cm}^{-3}$ ($\chi = 2 \times 10^3$, GC43.k3_low). Fig. 6 shows the column density map on the x-y plane of the simulation with $\bar{n}_c = 200 \text{ cm}^{-3}$ (left column) and $20 \text{ cm}^{-3}$ (right column) at different times (row-wise). We observed that the cloud in GC43.k3 undergoes gravitational collapse due to the compression from the wind. However, in the case of $\bar{n}_c = 20$, the value of $\chi$ is one order of magnitude lower, resulting in a much smaller cloud-crushing time ($\sim 1$ Myr) and drag time ($\sim 44$ Myr). Hence, the wind rapidly disintegrates the cloud in the $\bar{n}_c = 20 \text{ cm}^{-3}$ case, giving rise to a filamentary structure before the gravitational force has a chance to significantly impact the evolution. Additionally, due to the relatively shorter drag time, the fragmented cloudlets attain a sufficiently high velocity to effectively overcome the gravitational potential and get dispersed.

<footnote n="4">⁴Here, we show the density slices instead of the column depth, as presented earlier, to better demonstrate the percolation of the wind within the clouds for different values of the $k_{min}$ parameter.</footnote>
runaway gravitational collapse before any instability has a chance to act. This is clearly visible in the top left panel, where the cloud in GC42.k3 has collapsed by this time, becoming very compact and highly dense. None the less, due to the continuous ablation, the loosely bound outer layer of the cloud has been gradually stripped away by the wind, giving rise to a tail in the direction of the wind. In GC43.k3 (top right panel), the cloud is undergoing collapse but at a somewhat lower rate than in the $10^{42}$ erg s⁻¹ case due to the stronger wind and comparatively quicker growth of shear instabilities (such as KH instability). Thus, even though the cloud’s central region is experiencing collapse, the wind’s influence extends to a larger section of the cloud, resulting in a more extended structure compared to the $10^{42}$ erg s⁻¹ case.

The cloud disruption time-scales in GC44.k3 (bottom left panel) and GC45.k3 (bottom right panel) are much shorter than the free-fall time due to increased velocity and ram pressure (see Table 1). Therefore, a high-power wind induces much stronger turbulence and percolation of gas inside the inhomogeneous cloud, preventing the runaway collapse globally, compared to a low-power wind. Thus, in this scenario, the cloud is ablated and disrupted significantly before the influence of gravity becomes important. However, even in the case of a strong wind, the small cloudlets formed by the fragmentation become gravitationally bound and shielded from the external wind, by forming a high-density post-shock outer layer both due to compression and radiative cooling, preventing the wind material from infiltrating these cloudlets. Therefore, in the absence of any opposing mechanism, very high-density cloudlets collapse locally, while the wind ablates comparatively low-density material, shaping the cloud into an elongated and extended structure of cloudlets.

However, it is important to note that the density contrast between the cloud and wind material in these simulations is very high ($\chi \sim 2 \times 10^4$), which results in a very high value of the drag time-scale ($t_{drag} \sim \sqrt{\chi} t_{cc}$), which is the time by which the cloudlets are expected to attain a similar velocity as the wind. Therefore, even in scenarios of high power, where the cloud’s overall collapse is impeded by the wind, a considerable portion of the fragmented cloudlets persists and remains gravitationally bound to the central potential well due to the inability to attain escape velocity, owing to the large drag time-scale.

### 3.1.4 Dependence on cloud fractal wavenumber

The maximum size ($\lambda_{max} \approx L/2\pi k_{min}$) of individual cloudlets inside the cloud is parameterized by the minimum wavenumber $k_{min}$. Therefore, a higher value of $k_{min}$ results in a larger number of small cloudlets as well as low-density intercloudlet channels. As discussed earlier, the interaction between the wind and a fractal cloud is influenced by the extent to which the wind material can permeate the medium between individual clumps separated by low-density cloud material, in addition to its impact on the fractal surfaces of the cloud. Thus, one can anticipate that in the case of a higher value of $k_{min}$ with a greater abundance of low-density channels, it is easier for the wind to percolate through these channels into the cloudlets and mix with cloud material, thereby transferring energy and momentum. Additionally, due to the reduced size of the cloudlets in this case, KH instabilities grow faster ($t_{KH} \propto 1/k$) in a cloud with lower $k_{min}$. Consequently, it is expected that a cloud with a higher $k_{min}$ value would undergo more pronounced disruption.

Hence, in order to investigate the impact of the cloud fractal wavenumber, we simulate wind-cloud interaction for clouds with

varying values of $k_{min}$ (=1, 3, 6, 10), while keeping the mass and mean density of the clouds fixed. The simulations are initiated with self-gravity and a wind of power $10^{43}$ erg s⁻¹. In Fig. 5, we show the number density slice<fnref n="4" /> through the $z=0$ plane for these four cases (row-wise) at 0 Myr (left panel) and 1.27 Myr (right panel).

If we first consider the $k_{min} = 1$ case ('k1', first row), the cloud contains approximately one big cloudlet with increasing density towards its centre, surrounded by low-density outer layers, which are rapidly eroded by the wind through ablation. Nevertheless, owing to its larger size and the added shielding from compression and cooling, the shock is unable to penetrate deep into the clump, at least not at early times compared to the other $k_{min}$ cases. Furthermore, external overpressurization triggers the gravitational collapse of the clump to form a compact, gravitationally bound structure.

In the $k_{min} = 3$ case ('k3', second row), the cloud is initially composed of many small, dense clumps separated by wide interclump channels. These channels help the high-velocity wind material to percolate into the cloud, transferring momentum. This results in more fragmentation of the cloud, preventing the global collapse, unlike in the $k_{min} = 1$ case.

Clouds with smaller scale density perturbations as in the $k_{min} = 6$ ('k6') and 10 ('k10') cases (third and fourth row), the numbers of cloudlets are much higher, and the intercloud channels are also consequently narrower. During the initial interaction, the swept-up material by the wind creates a dense layer near the interaction area, blocking the entrance of these narrow channels. This prevents the wind from dispersing the cloudlets efficiently, resulting in the accumulation near the centre. Thus, external overpressurization triggers the global collapse of the cloud in these cases similar to the 'k1' case. However, a striking difference between the 'k1' case and the $k_{min} = 6, 10$ cases is evident. In the former, cloudlets form through fragmentation and stripping of the large clump. In contrast, in the 'k6' and 'k10' cases, cloudlets were initially seeded by the fractal generator itself at the initialization of the simulation and they remain organized in a more spherically symmetric manner at later times.

### 3.1.5 Dependence on the mean cloud density

In order to understand the influence of the density contrast between the cloud and wind, we consider two simulations including self-gravity with identical wind power of $10^{43}$ erg s⁻¹, but two different values of the initial mean cloud density, namely, $\bar{n}_c = 200 \text{ cm}^{-3}$ ($\chi = 2 \times 10^4$, GC43.k3) and $20 \text{ cm}^{-3}$ ($\chi = 2 \times 10^3$, GC43.k3_low). Fig. 6 shows the column density map on the x-y plane of the simulation with $\bar{n}_c = 200 \text{ cm}^{-3}$ (left column) and $20 \text{ cm}^{-3}$ (right column) at different times (row-wise). We observed that the cloud in GC43.k3 undergoes gravitational collapse due to the compression from the wind. However, in the case of $\bar{n}_c = 20$, the value of $\chi$ is one order of magnitude lower, resulting in a much smaller cloud-crushing time ($\sim 1$ Myr) and drag time ($\sim 44$ Myr). Hence, the wind rapidly disintegrates the cloud in the $\bar{n}_c = 20 \text{ cm}^{-3}$ case, giving rise to a filamentary structure before the gravitational force has a chance to significantly impact the evolution. Additionally, due to the relatively shorter drag time, the fragmented cloudlets attain a sufficiently high velocity to effectively overcome the gravitational potential and get dispersed.

<footnote n="4">⁴Here, we show the density slices instead of the column depth, as presented earlier, to better demonstrate the percolation of the wind within the clouds for different values of the $k_{min}$ parameter.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Number density slice through the $z=0$ plane of self-gravity simulations with $P_w = 10^{43}$ erg s$^{-1}$ and different cloud porosity: $k_{\min} = 1$ (first row), 3 (second row), 6 (third row), and 10 (fourth row). The left column shows the initial number density slice of each simulation. The right column represents the distribution at 1.27 Myr, which corresponds to 0.38 $t_{cc}$ and 0.33 $t_{ff}$.</figcaption></figure>

### 3.1.6 Thermal wind versus kinetic wind

Until now, the energy budget of the wind – specifically, whether the energy carried by the wind is primarily thermally dominated (subsonic) or kinetically dominated (supersonic) – has received limited attention in the context of the traditional ‘cloud-crushing’ problem. Most of the previous studies focus on the impact of a supersonic wind on the cloud, which generally holds true for galactic or starburst-driven winds. However, as indicated by various

theoretical and numerical investigations (King 2003; King, Zubovas & Power 2011; Faucher-Giguère & Quataert 2012; Zubovas & King 2012; Costa et al. 2014), the progression of AGN-driven winds can encompass various evolutionary stages – such as pressure-dominated and kinetic energy-dominated phases – contingent upon diverse parameters, including the black hole mass, the launch velocity of the wind at the accretion scale, ambient density profile, and various cooling mechanisms, among others (e.g. Faucher-Giguère
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Number density slice through the $z = 0$ plane of self-gravity simulations with $P_w = 10^{43}$ erg s$^{-1}$ and different cloud porosity: $k_{\min} = 1$ (first row), 3 (second row), 6 (third row), and 10 (fourth row). The left column shows the initial number density slice of each simulation. The right column represents the distribution at 1.27 Myr, which corresponds to 0.38 $t_{cc}$ and 0.33 $t_{ff}$.</figcaption></figure>

### 3.1.6 Thermal wind versus kinetic wind

Until now, the energy budget of the wind – specifically, whether the energy carried by the wind is primarily thermally dominated (subsonic) or kinetically dominated (supersonic) – has received limited attention in the context of the traditional ‘cloud-crushing’ problem. Most of the previous studies focus on the impact of a supersonic wind on the cloud, which generally holds true for galactic or starburst-driven winds. However, as indicated by various

theoretical and numerical investigations (King 2003; King, Zubovas & Power 2011; Faucher-Giguère & Quataert 2012; Zubovas & King 2012; Costa et al. 2014), the progression of AGN-driven winds can encompass various evolutionary stages – such as pressure-dominated and kinetic energy-dominated phases – contingent upon diverse parameters, including the black hole mass, the launch velocity of the wind at the accretion scale, ambient density profile, and various cooling mechanisms, among others (e.g. Faucher-Giguère
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Time evolution of projected number density along the z-direction of self-gravity simulations with mean cloud density of 200 cm⁻³ (left panel) and 20 cm⁻³ (right panel). The power of the outflow for both simulations is 10⁴³ erg s⁻¹. Here, the lower density cloud gets quickly disrupted by the wind before self-gravity becomes important compared to the higher density cloud where we observe a significantly higher fraction of high-density gas.</figcaption></figure>

& Quataert 2012). Thus, it is useful to investigate how different kinds of wind, despite having the same power, affect the evolution of the cloud, as pressure and momentum couple differently to the HDs.

Thus, to examine the effect of thermal versus kinetic wind on the evolution of the cloud, we consider two different simulations with the same initial cloud configuration and wind power of 10⁴⁵ erg s⁻¹, but varying wind properties. One simulation involves a thermal wind (GC45_k3), characterized by a Mach number ($M_w$) of 0.28 (as used in the previous simulation comparisons in this study so far), while the other initializes a kinetic wind (GC45_k3_kinetic, see Table 1), with a higher velocity (15 000 km s⁻¹) but lower pressure (10⁻¹⁰ dyne cm⁻²) such that total power is 10⁴⁵ erg s⁻¹, resulting in $M_w$ = 12.16. Fig. 7 displays the column density map in the x-y plane of the simulation with the thermal wind (left column) and the kinetic wind (right column) at different times (row-wise). Evidently, the cloud impacted by the kinetic wind undergoes a higher fraction of ablation across all stages of its evolution, compared to the thermal wind's effect. We

observe a significantly larger amount of gas in the low-density tail in the kinetic wind case, which has been stripped from the original cloud by the wind. Although possessing the same power, the differences in the effect of the thermal and kinetic wind are due to the fact that cloud material is primarily entrained and accelerated by the direct momentum transfer from the wind compared to the PdV work (Wagner et al. 2012; Wagner, Umemura & Bicknell 2013; Nayakshin 2014). With the thermal wind containing less kinetic energy and consequently, lower ram pressure compared to the kinetic energy-dominated wind, the mixing and acceleration of the cloudlets are significantly lower for the thermal wind case. In contrast, due to higher ram pressure, the direct momentum transfer is much higher for the case of the kinetic wind, leading to increased turbulence and vorticity, and causing the cloud to expand kinetically. As a result, the kinetic wind, despite possessing the same power, is expected to exhibit more destructive behaviour than a thermal wind with equivalent power.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Time evolution of projected number density along the z-direction of self-gravity simulations with mean cloud density of 200 cm⁻³ (left panel) and 20 cm⁻³ (right panel). The power of the outflow for both simulations is 10⁴³ erg s⁻¹. Here, the lower density cloud gets quickly disrupted by the wind before self-gravity becomes important compared to the higher density cloud where we observe a significantly higher fraction of high-density gas.</figcaption></figure>

& Quataert 2012). Thus, it is useful to investigate how different kinds of wind, despite having the same power, affect the evolution of the cloud, as pressure and momentum couple differently to the HDs.

Thus, to examine the effect of thermal versus kinetic wind on the evolution of the cloud, we consider two different simulations with the same initial cloud configuration and wind power of 10⁴⁵ erg s⁻¹, but varying wind properties. One simulation involves a thermal wind (GC45_k3), characterized by a Mach number ($M_w$) of 0.28 (as used in the previous simulation comparisons in this study so far), while the other initializes a kinetic wind (GC45_k3_kinetic, see Table 1), with a higher velocity (15 000 km s⁻¹) but lower pressure (10⁻¹⁰ dyne cm⁻²) such that total power is 10⁴⁵ erg s⁻¹, resulting in $M_w$ = 12.16. Fig. 7 displays the column density map in the x-y plane of the simulation with the thermal wind (left column) and the kinetic wind (right column) at different times (row-wise). Evidently, the cloud impacted by the kinetic wind undergoes a higher fraction of ablation across all stages of its evolution, compared to the thermal wind's effect. We

observe a significantly larger amount of gas in the low-density tail in the kinetic wind case, which has been stripped from the original cloud by the wind. Although possessing the same power, the differences in the effect of the thermal and kinetic wind are due to the fact that cloud material is primarily entrained and accelerated by the direct momentum transfer from the wind compared to the PdV work (Wagner et al. 2012; Wagner, Umemura & Bicknell 2013; Nayakshin 2014). With the thermal wind containing less kinetic energy and consequently, lower ram pressure compared to the kinetic energy-dominated wind, the mixing and acceleration of the cloudlets are significantly lower for the thermal wind case. In contrast, due to higher ram pressure, the direct momentum transfer is much higher for the case of the kinetic wind, leading to increased turbulence and vorticity, and causing the cloud to expand kinetically. As a result, the kinetic wind, despite possessing the same power, is expected to exhibit more destructive behaviour than a thermal wind with equivalent power.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time evolution (row-wise) of the projected number density along the z-direction of self-gravity simulations with thermal (left panel) and kinetic (right panel) winds of the same power ($10^{45}$ erg s⁻¹). The cloud impacted by the kinetic wind is more extended and elongated. Also, there exists a higher fraction of low-density mixed gas (blue colour) in the kinetic wind case.</figcaption></figure>

## 3.2 Cloud dynamics

In this section, our main emphasis lies in examining how various initial conditions influence the dynamical evolution of the cloud, i.e. mass-loss, gas turbulence, cloud elongation, acceleration, etc. In order to extract a particular quantity ($\Psi$) for the cloud material from the whole simulation domain, we use the definition of the mass-weighted volume average of $\Psi$ by (Banda-Barragán et al. 2018)

$$ \langle \Psi \rangle = \frac{\int \Psi \rho C dV}{\int \rho C dV}, \qquad (22) $$

where $\rho$ is the density and $C$ is the cloud tracer (a passive scalar that traces cloud material, i.e. if a cell contains only cloud material then $C = 1$; if the cell contains no cloud material $C = 0$, such that any mixture of cloud and non-cloud material in a cell can be represented by $C$). To exclude significant fractions of wind and hot mixed material from the cloud, we impose a threshold value on the temperature, i.e. only cells with $T < 10^4$ K are used in equation (22). This gives us a

reliable estimate of a particular quantity of cold gas embedded in a hot wind.

One common convention in all the line plots presented in this section is that the greyed-out parts of the lines (if they exist) correspond to the regime where the Jeans length of the highest density cell is not resolved by four resolution elements (e.g. see Fig. 9). This occurs exclusively in the self-gravitating simulations, where few cells can reach very large densities, whose Jeans length can not be resolved by at least four computational cells, which is necessary to avoid artificial fragmentation (Truelove et al. 1997).

### 3.2.1 Density PDF

In this section, we investigate the influence of self-gravity and wind power on the distribution of cloud density. The left panel of Fig. 8 illustrates the number density distribution at $t = 0.9 t_{ff} = 3.48$ Myr in simulations with self-gravity (solid line) and without (dashed line), considering a wind power of $10^{42}$ erg s⁻¹. The grey dashed-dotted line represents the initial probability distribution function (PDF) for
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time evolution (row-wise) of the projected number density along the z-direction of self-gravity simulations with thermal (left panel) and kinetic (right panel) winds of the same power ($10^{45}$ erg s⁻¹). The cloud impacted by the kinetic wind is more extended and elongated. Also, there exists a higher fraction of low-density mixed gas (blue colour) in the kinetic wind case.</figcaption></figure>

## 3.2 Cloud dynamics

In this section, our main emphasis lies in examining how various initial conditions influence the dynamical evolution of the cloud, i.e. mass-loss, gas turbulence, cloud elongation, acceleration, etc. In order to extract a particular quantity ($\Psi$) for the cloud material from the whole simulation domain, we use the definition of the mass-weighted volume average of $\Psi$ by (Banda-Barragán et al. 2018)

$$ \langle \Psi \rangle = \frac{\int \Psi \rho C dV}{\int \rho C dV}, \qquad (22) $$

where $\rho$ is the density and $C$ is the cloud tracer (a passive scalar that traces cloud material, i.e. if a cell contains only cloud material then $C = 1$; if the cell contains no cloud material $C = 0$, such that any mixture of cloud and non-cloud material in a cell can be represented by $C$). To exclude significant fractions of wind and hot mixed material from the cloud, we impose a threshold value on the temperature, i.e. only cells with $T < 10^4$ K are used in equation (22). This gives us a

reliable estimate of a particular quantity of cold gas embedded in a hot wind.

One common convention in all the line plots presented in this section is that the greyed-out parts of the lines (if they exist) correspond to the regime where the Jeans length of the highest density cell is not resolved by four resolution elements (e.g. see Fig. 9). This occurs exclusively in the self-gravitating simulations, where few cells can reach very large densities, whose Jeans length can not be resolved by at least four computational cells, which is necessary to avoid artificial fragmentation (Truelove et al. 1997).

### 3.2.1 Density PDF

In this section, we investigate the influence of self-gravity and wind power on the distribution of cloud density. The left panel of Fig. 8 illustrates the number density distribution at $t = 0.9 t_{ff} = 3.48$ Myr in simulations with self-gravity (solid line) and without (dashed line), considering a wind power of $10^{42}$ erg s⁻¹. The grey dashed-dotted line represents the initial probability distribution function (PDF) for
<figure><img src="image_9.png" /><figcaption>**Figure 8.** PDFs of log-density of the cloud material for different simulations. Left panel: PDFs in the simulations with (solid) and without (dashed) self-gravity for the wind power of $10^{42}$ erg s$^{-1}$ at $0.9 t_{ff} = 3.48$ Myr. Right panel: the PDFs of the cloud material in self-gravitating simulations with different wind power at $t = 0.26 t_{ff} \approx 1$ Myr. The initial density PDF of the cloud material is shown in both panels by the grey dashed-dotted lines. The vertical black dotted lines in both panels mark the density value beyond which the Jeans length cannot be resolved by at least four grid cells with the current computational setup.</figcaption></figure>

the cloud number density. Notably, we observe an increase in the low-density tail of the PDFs due to the stripping and mixing of the cloud material with the wind. This behaviour remains consistent regardless of the presence of self-gravity. On the other hand, the wind significantly compresses a substantial portion of intermediate-density gas (~10–100 cm$^{-3}$) to higher densities. Interestingly, the highest density regions, which are surrounded by comparatively lower density gas, exhibit minimal influence from the wind, as the high-density tail of the PDF in the simulation without self-gravity at 3.48 Myr (blue dashed line) coincides with the initial PDF. However, in the presence of self-gravity, the densities of these self-gravitating cores increase by accreting lower density gas from the surroundings, forming an extended, high-density PL tail (solid line), which is a prominent feature in turbulent clouds when self-gravity becomes important (e.g. Klessen 2000; Collins et al. 2011, 2012; Kritsuk, Norman & Wagner 2011; Federrath & Klessen 2013; Girichidis et al. 2014; Jaupart & Chabrier 2020; Khullar et al. 2021; Appel et al. 2022).

The right panel of Fig. 8 depicts the density PDF of the self-gravitating cloud at $t = 0.26 t_{ff} = 1$ Myr, for simulations with different wind power. This time approximately corresponds to the stopping point of the highest power simulation (GC45), where due to increased pressure, the simulation necessitates substantially lower time steps, rendering computations beyond this time prohibitively resource-intensive. None the less, due to the shorter cloud-crushing time-scales (see Table 1), all the major evolutionary stages of the clouds are well captured even for the relatively short duration of the simulations.

Notably, as the wind power increases, the compression time-scale by the wind decreases, resulting in more rapid growth of the high-density tail of the PDF. However, at this early stage of the simulations, the influence of self-gravity on the cloud’s evolution has not yet given rise to the emergence of the PL tail in the density PDF (see Fig. B1 where we show the density PDF at the same cloud-crushing time, which corresponds to different absolute times for different simulations). None the less, the mean of the PDF shifts to higher values as the wind power increases due to elevated compression ratios. Interestingly, for the higher power winds (GC44 and GC45), a secondary peak in the low-density regime of the PDF can be observed

(at $n \sim 10^{-2}$ and 1 cm$^{-3}$ for GC.44 and GC.45, respectively), owing to the presence of a significant portion of hot, diffuse, and mixed cloud material which has been stripped off from the cloud by the wind.

### 3.2.2 Evolution of the cold gas mass

A commonly employed method for understanding cloud survival is to calculate the evolution of the cloud mass with $n \ge \bar{n}_c/3$, where $\bar{n}_c$ is the mean number density of the cloud (e.g. Scannapieco & Brüggen 2015; Gronke & Oh 2018). However, as we are interested in the evolution of the cloud gas fraction (which is essential for potential star formation of the cloud), this particular definition might encompass certain portions of mixed and shock-heated dense gas. Thus, we define the cold gas mass as

$$M_{\text{cold}} = \int_{T < 10^4 \text{ K}} \rho C dV. \qquad (23)$$

In Fig. 9, we show the evolution of the cold gas mass for simulations with different wind power. The dashed (solid) lines correspond to simulations without (with) self-gravity. For comparison, we also present the results of simulations without a wind (red lines).

A general pattern for all the simulations is an initial rapid decline followed by subsequent growth over a brief time span and finally a gradual decline. During the initial interaction, the wind transfers a significant fraction of thermal energy into the cloud, resulting in an increase in the overall temperature of the cloud, which reduces the amount of cold gas. The fall-off becomes more pronounced with increasing wind power, as a larger amount of energy is transported into the cloud. It is important to note that, this decrease in cold gas mass is not due to cloud ablation, which is not significant during this phase.

After the initial interaction, when the shock-compressed gas cools efficiently, the excess internal energy is removed via radiative losses, resulting in the decrease of the mean temperature and the cold gas mass increases to aid the growth phase. None the less, the value of the cold cloud mass to which a particular simulation reaches after the initial impact with the cloud decreases with increasing wind power.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** PDFs of log-density of the cloud material for different simulations. Left panel: PDFs in the simulations with (solid) and without (dashed) self-gravity for the wind power of $10^{42}$ erg s$^{-1}$ at $0.9 t_{ff} = 3.48$ Myr. Right panel: the PDFs of the cloud material in self-gravitating simulations with different wind power at $t = 0.26 t_{ff} \approx 1$ Myr. The initial density PDF of the cloud material is shown in both panels by the grey dashed-dotted lines. The vertical black dotted lines in both panels mark the density value beyond which the Jeans length cannot be resolved by at least four grid cells with the current computational setup.</figcaption></figure>

the cloud number density. Notably, we observe an increase in the low-density tail of the PDFs due to the stripping and mixing of the cloud material with the wind. This behaviour remains consistent regardless of the presence of self-gravity. On the other hand, the wind significantly compresses a substantial portion of intermediate-density gas (~10–100 cm$^{-3}$) to higher densities. Interestingly, the highest density regions, which are surrounded by comparatively lower density gas, exhibit minimal influence from the wind, as the high-density tail of the PDF in the simulation without self-gravity at 3.48 Myr (blue dashed line) coincides with the initial PDF. However, in the presence of self-gravity, the densities of these self-gravitating cores increase by accreting lower density gas from the surroundings, forming an extended, high-density PL tail (solid line), which is a prominent feature in turbulent clouds when self-gravity becomes important (e.g. Klessen 2000; Collins et al. 2011, 2012; Kritsuk, Norman & Wagner 2011; Federrath & Klessen 2013; Girichidis et al. 2014; Jaupart & Chabrier 2020; Khullar et al. 2021; Appel et al. 2022).

The right panel of Fig. 8 depicts the density PDF of the self-gravitating cloud at $t = 0.26 t_{ff} = 1$ Myr, for simulations with different wind power. This time approximately corresponds to the stopping point of the highest power simulation (GC45), where due to increased pressure, the simulation necessitates substantially lower time steps, rendering computations beyond this time prohibitively resource-intensive. None the less, due to the shorter cloud-crushing time-scales (see Table 1), all the major evolutionary stages of the clouds are well captured even for the relatively short duration of the simulations.

Notably, as the wind power increases, the compression time-scale by the wind decreases, resulting in more rapid growth of the high-density tail of the PDF. However, at this early stage of the simulations, the influence of self-gravity on the cloud’s evolution has not yet given rise to the emergence of the PL tail in the density PDF (see Fig. B1 where we show the density PDF at the same cloud-crushing time, which corresponds to different absolute times for different simulations). None the less, the mean of the PDF shifts to higher values as the wind power increases due to elevated compression ratios. Interestingly, for the higher power winds (GC44 and GC45), a secondary peak in the low-density regime of the PDF can be observed

(at $n \sim 10^{-2}$ and 1 cm$^{-3}$ for GC.44 and GC.45, respectively), owing to the presence of a significant portion of hot, diffuse, and mixed cloud material which has been stripped off from the cloud by the wind.

### 3.2.2 Evolution of the cold gas mass

A commonly employed method for understanding cloud survival is to calculate the evolution of the cloud mass with $n \ge \bar{n}_c/3$, where $\bar{n}_c$ is the mean number density of the cloud (e.g. Scannapieco & Brüggen 2015; Gronke & Oh 2018). However, as we are interested in the evolution of the cloud gas fraction (which is essential for potential star formation of the cloud), this particular definition might encompass certain portions of mixed and shock-heated dense gas. Thus, we define the cold gas mass as

$$M_{\text{cold}} = \int_{T<10^4 \text{ K}} \rho C dV. \qquad (23)$$

In Fig. 9, we show the evolution of the cold gas mass for simulations with different wind power. The dashed (solid) lines correspond to simulations without (with) self-gravity. For comparison, we also present the results of simulations without a wind (red lines).

A general pattern for all the simulations is an initial rapid decline followed by subsequent growth over a brief time span and finally a gradual decline. During the initial interaction, the wind transfers a significant fraction of thermal energy into the cloud, resulting in an increase in the overall temperature of the cloud, which reduces the amount of cold gas. The fall-off becomes more pronounced with increasing wind power, as a larger amount of energy is transported into the cloud. It is important to note that, this decrease in cold gas mass is not due to cloud ablation, which is not significant during this phase.

After the initial interaction, when the shock-compressed gas cools efficiently, the excess internal energy is removed via radiative losses, resulting in the decrease of the mean temperature and the cold gas mass increases to aid the growth phase. None the less, the value of the cold cloud mass to which a particular simulation reaches after the initial impact with the cloud decreases with increasing wind power.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution of the cold gas mass defined in equation (23) as a function of time for simulations with different wind power, with and without self-gravity. The solid lines correspond to simulations with self-gravity, while the dashed lines with the same colour represent the same wind-power case, but without self-gravity. The red lines are the result of a simulation without any wind. The grey-out portions of the lines indicate the regime in which the Jeans length of the highest density cell in the simulations is not resolved by at least four computational cells.</figcaption></figure>

Subsequent to the growth phase, the cold gas mass in all the simulations decreases, primarily due to ablation, mixing and removal of the gas from the computational domain by the wind. As the wind continues to strip the cloud material due to various instabilities and cooling-driven pressure gradients (Gronke & Oh 2020), the cloud material gets mixed and heated by the wind. The effectiveness of mixing escalates in tandem with increasing wind power, leading to a quicker decline in the amount of surviving cold gas under the influence of the hot wind. However, in the simulations without a wind (red lines), the late time decrease in the cold gas mass is due to the escape of a fraction of mass through the computational boundary, which is stirred by the initial random velocity field. This is more prominent in the case without self-gravity (red-dashed line) due to the absence of attractive gravitational forces.

Interestingly, the simulations with self-gravity (solid lines) retain more cold gas compared to the cases without self-gravity (dashed lines) at the same wind power. The presence of self-gravity renders the shock-compressed clouds gravitationally bound and more compact, thus, reducing wind's ability to ablate the cloud and induce subsequent mixing. Furthermore, the compact cloud provides less surface area for the wind to interact with, compared to that of an expanded cloud in the absence of self-gravity. This combination of factors contributes to a higher fraction of cold gas in simulations with self-gravity. However, the difference between the ensuing evolution in simulations with and without self-gravity becomes less pronounced for high-power winds ($P_w \gtrsim 10^{44}$ erg s⁻¹). due to the shorter cloud-crushing time-scales compared to the freefall time of the cloud.

### 3.2.3 Gas turbulence

Quantifying the generation of turbulence inside a cloud impacted by the wind is one of the important factors because this can regulate star formation (e.g. Federrath 2018). As the wind progresses

through the fractal cloud, various shear layers are created inside the intercloudlet medium, depending on the local density, which leads to the generation of vorticity ($\omega = \nabla \times \mathbf{v}$). Additionally, when acting on an inhomogeneous medium, the self-gravitational force causes internal motion between the clumps, creating vorticity inside the cloud (Federrath et al. 2011), which can further be amplified by the wind. In order to quantify the turbulence inside the cloud, we define the 3D mass-weighted velocity dispersion ($\sigma_V$) as

$$ \sigma_V = \sqrt{\sum_{i=1}^{3} \sigma_{v_i}^2}, \qquad (24) $$

with

$$ \sigma_{v_i} = \sqrt{\langle v_i^2 \rangle - \langle v_i \rangle^2}, \qquad (25) $$

where $\sigma_{v_i}$ is the 1D velocity dispersion along each axis.

The left panel of Fig. 10 shows the evolution of $\sigma_V$ for different simulations with and without self-gravity as a function of $t_{ff}$ (see Fig. C1 for the evolutionary trends in terms of the cloud-crushing time in Appendix C). Irrespective of the wind power, the velocity dispersion in all the wind simulations increases due to energy transfer from the wind to the cloud material. Furthermore, the magnitude of the enhancement increases with increasing wind power. In all cases, there is an initial enhancement of the velocity dispersion during the compression phase. Consequently, in the simulations without self-gravity, when the wind disperses the cloud, it becomes comparatively easier for the wind to traverse through the intercloudlet channels. Therefore, the velocity dispersion settles to an approximately constant value, depending on the strength of the wind. In the cases, where self-gravity is included, the evolution deviates when gravity starts to dominate. In these instances, the fragmented cloudlets are pulled to the central gravitational potential well, towards the core, which induces additional motions between the cloudlets, resulting in a higher velocity dispersion compared to the cases without self-gravity. It is worth mentioning that, as the wind power increases, the disparity in velocity dispersion evolution between simulations with and without self-gravity diminishes, as a stronger wind disperses the cloud before gravity can significantly influence its evolution.

To thoroughly gauge the wind's influence on star formation, it's inadequate to solely consider velocity dispersion. The crucial factor lies in the balance between turbulent kinetic energy ($E_{kin}$) and gravitational binding energy ($E_{grav}$), which determines whether a gas cloud experiences runaway collapse or maintains stability through turbulent support. Thus, we consider the virial parameter ($\alpha_{vir}$) of the cloud, defined as the ratio between the kinetic and gravitational energy (Bertoldi & McKee 1992; Federrath & Klessen 2012),

$$ \alpha_{vir} = \frac{2E_{kin}}{|E_{grav}|}, \qquad (26) $$

where $E_{kin}$ and $E_{grav}$ are defined as

$$ E_{kin} = \frac{1}{2} M_{\text{cold}} \sigma_V^2, \qquad (27) $$

$$ E_{grav} = \int_{T < 10^4 \text{ K}} \Phi \rho C dV, \qquad (28) $$

where $\Phi$ is the gravitational potential. In the simulations without self-gravity, we also solve for the gravitational potential at runtime, but do not couple it to the HDs.

The right-hand panel of Fig. 10 shows the evolution of the virial parameter with time. In most of the simulations, except for the cases with wind power $10^{42}$ erg s⁻¹, the virial parameter initially increases beyond unity because of the infusion of kinetic energy from the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution of the cold gas mass defined in equation (23) as a function of time for simulations with different wind power, with and without self-gravity. The solid lines correspond to simulations with self-gravity, while the dashed lines with the same colour represent the same wind-power case, but without self-gravity. The red lines are the result of a simulation without any wind. The grey-out portions of the lines indicate the regime in which the Jeans length of the highest density cell in the simulations is not resolved by at least four computational cells.</figcaption></figure>

Subsequent to the growth phase, the cold gas mass in all the simulations decreases, primarily due to ablation, mixing and removal of the gas from the computational domain by the wind. As the wind continues to strip the cloud material due to various instabilities and cooling-driven pressure gradients (Gronke & Oh 2020), the cloud material gets mixed and heated by the wind. The effectiveness of mixing escalates in tandem with increasing wind power, leading to a quicker decline in the amount of surviving cold gas under the influence of the hot wind. However, in the simulations without a wind (red lines), the late time decrease in the cold gas mass is due to the escape of a fraction of mass through the computational boundary, which is stirred by the initial random velocity field. This is more prominent in the case without self-gravity (red-dashed line) due to the absence of attractive gravitational forces.

Interestingly, the simulations with self-gravity (solid lines) retain more cold gas compared to the cases without self-gravity (dashed lines) at the same wind power. The presence of self-gravity renders the shock-compressed clouds gravitationally bound and more compact, thus, reducing wind's ability to ablate the cloud and induce subsequent mixing. Furthermore, the compact cloud provides less surface area for the wind to interact with, compared to that of an expanded cloud in the absence of self-gravity. This combination of factors contributes to a higher fraction of cold gas in simulations with self-gravity. However, the difference between the ensuing evolution in simulations with and without self-gravity becomes less pronounced for high-power winds ($P_w \gtrsim 10^{44}$ erg s⁻¹). due to the shorter cloud-crushing time-scales compared to the freefall time of the cloud.

### 3.2.3 Gas turbulence

Quantifying the generation of turbulence inside a cloud impacted by the wind is one of the important factors because this can regulate star formation (e.g. Federrath 2018). As the wind progresses

through the fractal cloud, various shear layers are created inside the intercloudlet medium, depending on the local density, which leads to the generation of vorticity ($\omega = \nabla \times \mathbf{v}$). Additionally, when acting on an inhomogeneous medium, the self-gravitational force causes internal motion between the clumps, creating vorticity inside the cloud (Federrath et al. 2011), which can further be amplified by the wind. In order to quantify the turbulence inside the cloud, we define the 3D mass-weighted velocity dispersion ($\sigma_V$) as

$$ \sigma_V = \sqrt{\sum_{i=1}^{3} \sigma_{v_i}^2}, \qquad (24) $$

with

$$ \sigma_{v_i} = \sqrt{\langle v_i^2 \rangle - \langle v_i \rangle^2}, \qquad (25) $$

where $\sigma_{v_i}$ is the 1D velocity dispersion along each axis.

The left panel of Fig. 10 shows the evolution of $\sigma_V$ for different simulations with and without self-gravity as a function of $t_{ff}$ (see Fig. C1 for the evolutionary trends in terms of the cloud-crushing time in Appendix C). Irrespective of the wind power, the velocity dispersion in all the wind simulations increases due to energy transfer from the wind to the cloud material. Furthermore, the magnitude of the enhancement increases with increasing wind power. In all cases, there is an initial enhancement of the velocity dispersion during the compression phase. Consequently, in the simulations without self-gravity, when the wind disperses the cloud, it becomes comparatively easier for the wind to traverse through the intercloudlet channels. Therefore, the velocity dispersion settles to an approximately constant value, depending on the strength of the wind. In the cases, where self-gravity is included, the evolution deviates when gravity starts to dominate. In these instances, the fragmented cloudlets are pulled to the central gravitational potential well, towards the core, which induces additional motions between the cloudlets, resulting in a higher velocity dispersion compared to the cases without self-gravity. It is worth mentioning that, as the wind power increases, the disparity in velocity dispersion evolution between simulations with and without self-gravity diminishes, as a stronger wind disperses the cloud before gravity can significantly influence its evolution.

To thoroughly gauge the wind's influence on star formation, it's inadequate to solely consider velocity dispersion. The crucial factor lies in the balance between turbulent kinetic energy ($E_{kin}$) and gravitational binding energy ($E_{grav}$), which determines whether a gas cloud experiences runaway collapse or maintains stability through turbulent support. Thus, we consider the virial parameter ($\alpha_{vir}$) of the cloud, defined as the ratio between the kinetic and gravitational energy (Bertoldi & McKee 1992; Federrath & Klessen 2012),

$$ \alpha_{vir} = \frac{2E_{kin}}{|E_{grav}|}, \qquad (26) $$

where $E_{kin}$ and $E_{grav}$ are defined as

$$ E_{kin} = \frac{1}{2} M_{\text{cold}} \sigma_V^2, \qquad (27) $$

$$ E_{grav} = \int_{T < 10^4 \text{ K}} \Phi \rho C dV, \qquad (28) $$

where $\Phi$ is the gravitational potential. In the simulations without self-gravity, we also solve for the gravitational potential at runtime, but do not couple it to the HDs.

The right-hand panel of Fig. 10 shows the evolution of the virial parameter with time. In most of the simulations, except for the cases with wind power $10^{42}$ erg s⁻¹, the virial parameter initially increases beyond unity because of the infusion of kinetic energy from the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Time evolution of the velocity dispersion (left panel) and virial parameter (right panel) as for different simulations. The grey-out parts of the plots 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>

wind, which initially surpasses the overall gravitational energy. For the highest power simulation (GC45), this increase in $\alpha_{vir}$ is more than one order of magnitude (~20) higher than the initial value of ~0.9. Nevertheless, following the compression phase, there is a rise in the average cloud density, intensifying the gravitational potential well, ultimately leading to a reduction in the value of $\alpha_{vir}$. In the simulations without self-gravity, the subsequent evolution becomes nearly steady-state. This suggests that the kinetic and gravitational energy maintain a delicate balance, remaining relatively constant over time. In the presence of self-gravity, the highly dense collapsing cores give rise to a very deep gravitational potential, leading to the domination of gravitational binding energy over the kinetic energy. This results in a very low value of $\alpha_{vir}$ below 1, which is particularly prominent in low-power wind cases.

It is important to acknowledge that due to the absence of a sink particles algorithm (Federrath et al. 2010b) in our simulations, we cannot take advantage of removing gas from a cell whose local Jeans length becomes unresolved by the computational grid (Truelove et al. 1997). Therefore, in scenarios where a region experiences runaway collapse, the density within the specific cell at the bottom of the gravitational potential well becomes exceptionally high as it accumulates matter. This situation results in a violation of the Truelove et al. (1997) criterion, which suggests that the Jeans length of any region should be resolved by at least four grid cells to prevent artificial fragmentation. None the less, even if we do not resolve the small-scale structures, the qualitative behaviour in the virial parameter will not change significantly as the potential well created by this very high-density region will still be very deep, even in higher resolution simulations.

### 3.2.4 Cloud elongation

Here, we discuss the role of self-gravity on the compression and subsequent elongation of the cloud, impacted by the wind, from a quantitative point of view. We calculate the effective length ($\mathcal{X}$) of the cloud along the x-direction (i.e. along the direction of the wind) as follows (see Klein et al. 1994; Banda-Barragán et al. 2018):

$$ \mathcal{X} \propto \sqrt{\langle (x^2) - \langle x \rangle^2 \rangle}, \qquad (29) $$

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Time evolution of the longitudinal extent (equation 29, normalized to the initial value) of the cloud for different simulations. We see that, in the presence of self-gravity, the clouds are more compact.</figcaption></figure>

where equation (22) is used to calculate the average quantities. Fig. 11 shows the evolution of the cloud length (normalized to the initial value), parallel to the wind (x-direction), for different simulations. The solid and dashed lines correspond to simulations with and without self-gravity, respectively. We see that in all scenarios, the cloud width initially decreases as a result of the compression. Consequently, the cloud undergoes elongation due to the combined effects of shock- and turbulence-induced expansion. In the absence of self-gravity, this elongation is much more prominent for the lower power winds, as depicted by the blue dashed line (C.42) in Fig. 11. In contrast, due to the global collapse of the clouds in the presence of self-gravity in cases with lower wind power, the clouds are more compact leading to a smaller longitudinal width as can be seen from the solid blue line (GC.42). Nevertheless, as the wind power increases, the distinction between these two scenarios becomes less pronounced, as the wind
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Time evolution of the velocity dispersion (left panel) and virial parameter (right panel) as for different simulations. The grey-out parts of the plots 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>

wind, which initially surpasses the overall gravitational energy. For the highest power simulation (GC45), this increase in $\alpha_{vir}$ is more than one order of magnitude (~20) higher than the initial value of ~0.9. Nevertheless, following the compression phase, there is a rise in the average cloud density, intensifying the gravitational potential well, ultimately leading to a reduction in the value of $\alpha_{vir}$. In the simulations without self-gravity, the subsequent evolution becomes nearly steady-state. This suggests that the kinetic and gravitational energy maintain a delicate balance, remaining relatively constant over time. In the presence of self-gravity, the highly dense collapsing cores give rise to a very deep gravitational potential, leading to the domination of gravitational binding energy over the kinetic energy. This results in a very low value of $\alpha_{vir}$ below 1, which is particularly prominent in low-power wind cases.

It is important to acknowledge that due to the absence of a sink particles algorithm (Federrath et al. 2010b) in our simulations, we cannot take advantage of removing gas from a cell whose local Jeans length becomes unresolved by the computational grid (Truelove et al. 1997). Therefore, in scenarios where a region experiences runaway collapse, the density within the specific cell at the bottom of the gravitational potential well becomes exceptionally high as it accumulates matter. This situation results in a violation of the Truelove et al. (1997) criterion, which suggests that the Jeans length of any region should be resolved by at least four grid cells to prevent artificial fragmentation. None the less, even if we do not resolve the small-scale structures, the qualitative behaviour in the virial parameter will not change significantly as the potential well created by this very high-density region will still be very deep, even in higher resolution simulations.

### 3.2.4 Cloud elongation

Here, we discuss the role of self-gravity on the compression and subsequent elongation of the cloud, impacted by the wind, from a quantitative point of view. We calculate the effective length ($\mathcal{X}$) of the cloud along the x-direction (i.e. along the direction of the wind) as follows (see Klein et al. 1994; Banda-Barragán et al. 2018):

$$ \mathcal{X} \propto \sqrt{(\langle x^2 \rangle - \langle x \rangle^2)}, \qquad (29) $$

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Time evolution of the longitudinal extent (equation 29, normalized to the initial value) of the cloud for different simulations. We see that, in the presence of self-gravity, the clouds are more compact.</figcaption></figure>

where equation (22) is used to calculate the average quantities. Fig. 11 shows the evolution of the cloud length (normalized to the initial value), parallel to the wind (x-direction), for different simulations. The solid and dashed lines correspond to simulations with and without self-gravity, respectively. We see that in all scenarios, the cloud width initially decreases as a result of the compression. Consequently, the cloud undergoes elongation due to the combined effects of shock- and turbulence-induced expansion. In the absence of self-gravity, this elongation is much more prominent for the lower power winds, as depicted by the blue dashed line (C.42) in Fig. 11. In contrast, due to the global collapse of the clouds in the presence of self-gravity in cases with lower wind power, the clouds are more compact leading to a smaller longitudinal width as can be seen from the solid blue line (GC.42). Nevertheless, as the wind power increases, the distinction between these two scenarios becomes less pronounced, as the wind
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Time evolution of the *x*-component of the centre of mass velocity of the cloud in different simulations as indicated in the legend.</figcaption></figure>

initiates the stripping of the cloud before gravity has a significant influence.

### 3.2.5 Cloud acceleration

Fig. 12 illustrates the evolution of the *x*-component of the cloud’s centre of mass velocity across various simulations. It’s evident that higher powered winds quickly accelerate the clouds due to a greater net momentum transfer from the wind to the cloud. Notably, despite a similar initial evolution, clouds in simulations incorporating self-gravity at the same wind power begin to decelerate (solid lines) as gravitational forces toward the central potential well counteract the wind-induced acceleration. It is essential to note that, within the chosen simulation parameters, acceleration time-scales significantly exceed the duration of the simulations themselves. As a result, the net velocity achieved by the cloud’s centre of mass remains significantly lower than the flow velocity. Nevertheless, the discernible impact of self-gravity on the cloud’s acceleration is clearly evident.

In Fig. 13, we depict the cloud’s net centre of mass velocity for GC45_k3_thermal (solid line) and GC45_k3_kinetic (dashed-dotted line) to analyse the impact of wind energy composition on cloud acceleration. As anticipated, the cloud exposed to the kinetic wind demonstrates significantly higher velocities compared to the thermal wind of equal power. This difference arises from the more efficient direct momentum transfer in the kinetic wind, contrasting with the thermal wind where injected internal energy dissipates rapidly due to strong radiative loss, diminishing momentum transfer efficiency.

## 3.3 Multiphase outflow

While the findings presented in Section 3.2.5 provide insights into the acceleration of the ‘cloud as a whole’ by the wind, it becomes increasingly evident that the definition of the cloud as a single object diminishes as the wind fragments the initial cloud within a spherical zone into numerous smaller pieces. This effect is particularly prominent for a fractal cloud, which is the case in this study. The fragmented and stripped cloud material with different densities and temperatures are entrained in the hot wind, resulting

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Same as Fig. 12 but for a simulation with a thermal wind (GC45_k3_thermal; solid line) compared to a kinetic wind (GC45_k3_kinetic; dashed-dotted line), with the same wind power ($10^{45}$ erg s⁻¹).</figcaption></figure>

in a multiphase outflow that can extend up to a few kpc from the centre of a galaxy, with velocities from a few 100 to several 1000 km s⁻¹. These kinds of ionized (Nesvadba et al. 2006, 2008; Holt, Tadhunter & Morganti 2008; Westmoquette et al. 2012; Harrison et al. 2014; Rupke, Gültekin & Veilleux 2017), neutral atomic (Morganti, Tadhunter & Oosterloo 2005; Nesvadba et al. 2010; Guillard et al. 2012; Cazzoli et al. 2016; Morganti et al. 2016), and molecular (Feruglio et al. 2010; Cicone et al. 2014; Fiore et al. 2017; Fluetsch et al. 2019) outflows in galaxies, resulting from AGN, have been observed and characterized by a plethora of observational studies (see Morganti 2017, for a review). Moreover, numerical simulations of the interaction between relativistic jets and interstellar clouds have demonstrated that atomic gas can form *in situ* inside the post-shock material which cools rapidly, at least for low-power jets (e.g. Perucho 2024). Thus, from a theoretical point of view, it is important to investigate and quantify how AGN-driven winds with different powers lead to such multiphase outflows.

### 3.3.1 Velocity distribution

Fig. 14 presents the mass-weighted 2D histogram of velocity versus number density for the cloud material, for four self-gravitating simulations with different wind power at $0.85t_{cc}$, which corresponds to different physical times for different power as indicated in the legends. The contours in each panel represent the initial distribution of the cloud material in the density–velocity plane. The impact of the wind on the cloud can be readily seen from the diagram, particularly in the high-power cases (GC44, GC45), where a significant portion of the cloud material has been compressed to higher densities and accelerated to higher velocities. The velocities of the dense gas with number densities of ~10–100 cm⁻³ (right side of the vertical line) is enhanced up to ~400 km s⁻¹, which constitutes a considerable fraction of the cloud’s mass. A notable amount of diffuse gas with densities of ~0.01–1 cm⁻³ is accelerated to velocities of several 1000 km s⁻¹. However, the velocities of the very dense gas (n ~ $10^3$ cm⁻³) are only mildly affected by the wind, as the gas with larger column densities is more resistant to acceleration through direct momentum transfer.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Time evolution of the *x*-component of the centre of mass velocity of the cloud in different simulations as indicated in the legend.</figcaption></figure>

initiates the stripping of the cloud before gravity has a significant influence.

### 3.2.5 Cloud acceleration

Fig. 12 illustrates the evolution of the *x*-component of the cloud’s centre of mass velocity across various simulations. It’s evident that higher powered winds quickly accelerate the clouds due to a greater net momentum transfer from the wind to the cloud. Notably, despite a similar initial evolution, clouds in simulations incorporating self-gravity at the same wind power begin to decelerate (solid lines) as gravitational forces toward the central potential well counteract the wind-induced acceleration. It is essential to note that, within the chosen simulation parameters, acceleration time-scales significantly exceed the duration of the simulations themselves. As a result, the net velocity achieved by the cloud’s centre of mass remains significantly lower than the flow velocity. Nevertheless, the discernible impact of self-gravity on the cloud’s acceleration is clearly evident.

In Fig. 13, we depict the cloud’s net centre of mass velocity for GC45_k3_thermal (solid line) and GC45_k3_kinetic (dashed-dotted line) to analyse the impact of wind energy composition on cloud acceleration. As anticipated, the cloud exposed to the kinetic wind demonstrates significantly higher velocities compared to the thermal wind of equal power. This difference arises from the more efficient direct momentum transfer in the kinetic wind, contrasting with the thermal wind where injected internal energy dissipates rapidly due to strong radiative loss, diminishing momentum transfer efficiency.

## 3.3 Multiphase outflow

While the findings presented in Section 3.2.5 provide insights into the acceleration of the ‘cloud as a whole’ by the wind, it becomes increasingly evident that the definition of the cloud as a single object diminishes as the wind fragments the initial cloud within a spherical zone into numerous smaller pieces. This effect is particularly prominent for a fractal cloud, which is the case in this study. The fragmented and stripped cloud material with different densities and temperatures are entrained in the hot wind, resulting

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Same as Fig. 12 but for a simulation with a thermal wind (GC45_k3_thermal; solid line) compared to a kinetic wind (GC45_k3_kinetic; dashed-dotted line), with the same wind power ($10^{45}$ erg s⁻¹).</figcaption></figure>

in a multiphase outflow that can extend up to a few kpc from the centre of a galaxy, with velocities from a few 100 to several 1000 km s⁻¹. These kinds of ionized (Nesvadba et al. 2006, 2008; Holt, Tadhunter & Morganti 2008; Westmoquette et al. 2012; Harrison et al. 2014; Rupke, Gültekin & Veilleux 2017), neutral atomic (Morganti, Tadhunter & Oosterloo 2005; Nesvadba et al. 2010; Guillard et al. 2012; Cazzoli et al. 2016; Morganti et al. 2016), and molecular (Feruglio et al. 2010; Cicone et al. 2014; Fiore et al. 2017; Fluetsch et al. 2019) outflows in galaxies, resulting from AGN, have been observed and characterized by a plethora of observational studies (see Morganti 2017, for a review). Moreover, numerical simulations of the interaction between relativistic jets and interstellar clouds have demonstrated that atomic gas can form *in situ* inside the post-shock material which cools rapidly, at least for low-power jets (e.g. Perucho 2024). Thus, from a theoretical point of view, it is important to investigate and quantify how AGN-driven winds with different powers lead to such multiphase outflows.

### 3.3.1 Velocity distribution

Fig. 14 presents the mass-weighted 2D histogram of velocity versus number density for the cloud material, for four self-gravitating simulations with different wind power at $0.85t_{cc}$, which corresponds to different physical times for different power as indicated in the legends. The contours in each panel represent the initial distribution of the cloud material in the density–velocity plane. The impact of the wind on the cloud can be readily seen from the diagram, particularly in the high-power cases (GC44, GC45), where a significant portion of the cloud material has been compressed to higher densities and accelerated to higher velocities. The velocities of the dense gas with number densities of ~10–100 cm⁻³ (right side of the vertical line) is enhanced up to ≲400 km s⁻¹, which constitutes a considerable fraction of the cloud’s mass. A notable amount of diffuse gas with densities of ~0.01–1 cm⁻³ is accelerated to velocities of several 1000 km s⁻¹. However, the velocities of the very dense gas (n ≳ $10^3$ cm⁻³) are only mildly affected by the wind, as the gas with larger column densities is more resistant to acceleration through direct momentum transfer.