# The energy and dynamics of trapped radiative feedback with stellar winds

Sam Geen, Rebekka Bieri, Alex de Koter, Taysun Kimm, Joakim A Rosdahl

► To cite this version:

Sam Geen, Rebekka Bieri, Alex de Koter, Taysun Kimm, Joakim A Rosdahl. The energy and dynamics of trapped radiative feedback with stellar winds. Monthly Notices of the Royal Astronomical Society, 2023, 526 (2), pp.1832-1849. 10.1093/mnras/stad2667. hal-04305960

HAL Id: hal-04305960

https://cnrs.hal.science/hal-04305960v1

Submitted on 29 Dec 2023

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 CC BY 4.0 - Attribution - International License
# The energy and dynamics of trapped radiative feedback with stellar winds

Sam Geen<sup>1,2,★</sup> Rebekka Bieri<sup>1,3</sup> Alex de Koter,<sup>1,4</sup> Taysun Kimm<sup>1,5</sup> and Joakim Rosdahl<sup>1,6</sup>

Accepted 2023 August 31. Received 2023 August 31; in original form 2023 January 31

## ABSTRACT

In this paper, we explore the significant, non-linear impact that stellar winds have on H II regions. We perform a parameter study using three-dimensional radiative magnetohydrodynamic simulations of wind and ultraviolet radiation feedback from a 35 M☉ star formed self-consistently in a turbulent, self-gravitating cloud, similar to the Orion Nebula (M42) and its main ionizing source $\theta^1$ Ori C. Stellar winds suppress early radiative feedback by trapping ionizing radiation in the shell around the wind bubble. Rapid breakouts of warm photoionized gas ('champagne flows') still occur if the star forms close to the edge of the cloud. The impact of wind bubbles can be enhanced if we detect and remove numerical overcooling caused by shocks crossing grid cells. However, the majority of the energy in the wind bubble is still lost to turbulent mixing between the wind bubble and the gas around it. These results begin to converge if the spatial resolution at the wind bubble interface is increased by refining the grid on pressure gradients. Wind bubbles form a thin chimney close to the star, which then expands outwards as an extended plume once the wind bubble breaks out of the dense core the star formed in, allowing them to expand faster than a spherical wind bubble. We also find wind bubbles mixing completely with the photoionized gas when the H II region breaks out of the cloud as a champagne flow, a process we term 'hot champagne'.

**Key words:** methods: numerical – stars: formation – stars: massive – stars: winds, outflows – ISM: clouds – H II regions.

## 1 INTRODUCTION

Massive stars strongly impact their environment via a set of processes called 'feedback', terminating their own formation through a combination of protostellar feedback processes (e.g. Kuiper & Hosokawa 2018; Bate 2019), clearing out their immediate formation environment (Chevance et al. 2023), shaping the interstellar medium (ISM) of galaxies (McKee & Ostriker 1977; Rathjen et al. 2021; Bieri et al. 2023), and even affecting the composition of the circumgalactic medium (Tumlinson, Peeples & Werk 2017) and driving cosmic reionization (e.g. Rosdahl et al. 2018).

The problem of feedback in molecular clouds stems from both the complexity of the environment and the interaction between multiple feedback mechanisms. High-energy radiation from stars plays a strong role in molecular clouds, both photoionizing the gas around the stars, causing thermal expansion and dispersal of the cloud (Kahn 1954) and pushing the gas away via direct radiation pressure (Mathews 1967).

Moreover, the strong ultraviolet (UV) radiation field of the star drives mass-loss from the star as stellar winds, which can reach out-flow velocities over 1000 km s⁻<fnref n="1" /> for massive stars (Puls, Sundqvist & Markova 2015). This shock-heats the gas around the star, creating a hot, X-ray emitting wind bubble that expands into the surrounding medium (Weaver et al. 1977).

Numerical simulations are a powerful tool for studying feedback because they self-consistently capture many of the complex features that emerge non-linearly from the interaction of different mechanisms in turbulent, magnetized, self-gravitating clouds. Simulations of photoionization in a cloud context have been carried out by, e.g. Dale et al. (2005), Gritschneder et al. (2009), and Peters et al. (2010), among many others. Photoionization-only feedback simulations with self-consistent star formation, e.g. Dale, Ercolano & Bonnell (2012), typically display rapid photoionization of large portions of the cloud material followed by dispersal of the cloud due to thermal expansion. Ali, Harries & Douglas (2018) also reproduce this for a system similar to that in M42, and find some similarities between their simulations and M42 in, e.g. ionized gas velocities.

A key feature of photoionization-driven cloud dispersal is the 'champagne' flow (Tenorio-Tagle 1979), in which a density gradient in the cloud triggers a rapid breakout of ionized gas. This can occur due to a sharp step in density as the cloud edge is reached (Whitworth 1979) or due to a steep density gradient away from the star in the cloud itself (Franco, Tenorio-Tagle & Bodenheimer 1990). Comeron (1997) also studied this effect with winds included, which produces an embedded wind plume within the photoionized outflow.

Similarly, simulations of winds in molecular clouds have been carried out by, e.g. Rogers & Pittard (2013), Rey-Raposo et al. (2017), Geen et al. (2021), and Lancaster et al. (2021c). Rosen et al. (2021) argue using a suite of simulations that winds alone can explain the

© 2023 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Anton Pannekoek Institute for Astronomy, Universiteit van Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Institute for Computational Science, University of Zurich, 8057 Zurich, Switzerland</footnote>

<footnote n="4"><sup>4</sup>Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven, Belgium</footnote>

<footnote n="5"><sup>5</sup>Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="6"><sup>6</sup>Centre de Recherche Astrophysique de Lyon UMR 5574, CNRS, Université de Lyon 1, ENS de Lyon, Université de Lyon, F-69230 Saint-Genis-Laval, France</footnote>

<footnote>* E-mail: s.t.geen@uva.nl</footnote>
sizes of feedback-driven shells around massive stars ($M > 10$ M$_{⊙}$), but not around lower mass stars.

The combination of winds and radiation was explored by Dale et al. (2014), who found that winds are relatively ineffective in contributing to stellar feedback in molecular clouds compared with photoionization. This is also found in simulations of larger spiral-arm regions by Ali, Bending & Dobbs (2022). Simulations by Grudić et al. (2022) and Guszejnov et al. (2022) establish that winds and radiation in concert effectively regulate star formation, though they also name protostellar jets as having a key role in both quenching star formation and in removing and disrupting matter accreting on to protostars. Supernovae generally occur too late to have a principal role in shaping star formation for molecular clouds with a global free-fall time $\le 7$ Myr.

## 1.1 Feedback efficiency and radiative cooling

The efficiency of stellar feedback depends not only on the total quantity of energy deposited into the ISM by stars, but also the ability for the ISM to retain this energy as kinetic flows that impact larger scales. For example, Walch et al. (2012) argue that typically $<0.1$ per cent of the energy from ionizing radiation is converted to kinetic energy in the gas via photoionization. This is because while this process is important in thermalizing the gas and driving outflows, most of the energy in photons is lost to radiative cooling of the gas at lower photon energies. The cooling rate becomes drastically higher in denser gas, particularly in very dense molecular clouds such as those in the central molecular zone (Barnes et al. 2020) or in dense shells swept up by feedback (e.g. Rahner et al. 2017; Geen & de Koter 2022). While ionizing radiation can have an additional impact via direct radiation pressure, this is typically weak at H II region radii more than a few parsecs (e.g. Olivier et al. 2021).

Stellar winds rely on the overpressurization of the hot bubble shocked by the stellar wind to drive expansion of the feedback-driven structures. The energy retained inside the wind bubble, which is mostly thermal energy, is roughly 45 per cent for an adiabatic wind bubble in a uniform medium (Weaver et al. 1977) or 33 per cent in a singular isothermal sphere density profile ($\rho \propto r^{-2}$, see Geen & de Koter 2022). The rest is used to do thermodynamic work on the gas swept up by the wind bubble.

There are various channels via which the pressure in the wind bubble can be lost. First, thermal evaporation of material from the dense wind-driven shell can cause mixing between hot, diffuse gas in the wind bubble and cold, dense gas, which enhances radiative cooling (Mac Low & McCray 1988). Kruijssen et al. (2019) argues that such cooling is inefficient on the time-scales for feedback in molecular clouds. Thermal evaporation can be enhanced by thermal conduction via electrons (Spitzer 1962), which Fierlinger et al. (2016) find reduces the energy in a typical wind bubble by 10 per cent when running one-dimensional (1D) simulations that include it. Finally, turbulent mixing between the hot, wind-shocked gas and denser gas outside (whether a cold, dense shell or a warm, moderately dense photoionized region) can be highly effective at enhancing cooling in the wind bubble. Early work by Kahn (1980) suggested that Kelvin Helmholz instabilities in the wind bubble cannot grow fast enough to encourage strong cooling. However, more recent observational analysis (Rosen et al. 2014) and analytic work (Lancaster et al. 2021a) confirmed using simulations (Lancaster et al. 2021b), as well as work on shocks in more generalized environments (Fielding et al. 2020; Tan, Oh & Gronke 2021) find that turbulent mixing, particularly in the context of structured or turbulent cloud environments, can effectively remove nearly all of

the energy injected by stellar winds from the wind bubble, causing it to expand in a momentum-driven rather than pressure-driven mode. This is also found in simulations that include both stellar winds and photoionization (Dwarkadas 2022).

However, numerical simulations have significant limitations in reproducing all of these mechanisms. Spatial resolution, which is at a premium for 3D simulations, limits how effectively each of these processes can be captured given extant limits on computational power (Pittard, Wareing & Kupilas 2021). The inclusion of winds creates an additional need to resolve the short time-steps demanded by including hot and/or fast wind-driven flows, with characteristic speeds exceeding 1000 km s$^{-1}$.

Additionally, certain numerical features of simulation codes can cause problems in accurately resolving the energetics of hot structures in the ISM. One consequence specific to Eulerian hydrodynamic codes (i.e. codes that use static computational elements to simulate moving flows) is that shocks must cross cells. In doing so, the contact discontinuity between the hot, shocked gas and cold, unshocked gas must exist inside a cell or set of cells. Gentry et al. (2017) show that, in the case of supernova-driven superwind bubbles, the artificial mixing of hot and cold gas that occurs in these cells causes considerable numerical overcooling compared with simulations using Lagrangian hydrodynamics (i.e. codes that use computational elements that move with flows, hence reducing while not completely eliminating the effect of this artificial mixing). It is thus crucial to understand in detail how the energy from feedback transfers to gas on larger scales if we are to have an accurate picture for how stars shape the wider universe.

The Orion region provides an excellent observational case study for feedback in a nearby star-forming region. Star formation has been ongoing for up to 10 Myr in this region (Alves & Bouy 2012; Rio et al. 2012; Drass et al. 2016), with evidence of ongoing massive stellar feedback in the Orion–Eridanus superbubble, powered by the OB1 cluster’s massive stars and supernovae (Brown, de Geus & de Zeeuw 1994). In this paper, we focus on the case study of the young M42 nebula in this region, powered by the Trapezium cluster, and in particular by the O7V star $\theta^1$ Ori C, with secondary contributions from the O9.5V star $\theta^2$ Ori A (O’Dell, Kollatschny & Ferland 2017). Simón-Díaz et al. (2006) analyse this cluster and find an age of $2.5 \pm 0.5$ Myr for the cluster as a whole, although they urge caution given the uncertainties in measuring the age of young stars where initial rotation velocities are poorly understood. By contrast, Pabst et al. (2019) analyse the velocity structure of the M42 nebula itself using 3D spectroscopic observations of $C^+$, finding an age of 0.3 Myr based on an analysis encompassing the expansion rate and spatial size of the nebula. Interpreting the age of young star-forming regions precisely remains a challenge, and the input from theoretical models of both stars and nebulae remains important in piecing together the history of our Galactic neighbourhood.

The goal of this paper is to assess how well numerical simulations capture the evolution and properties of H II regions powered by massive stars by focussing on the example of a single star in conditions similar to the solar neighbourhood and in particular M42 in Orion, where we have a wealth of high-resolution multi-wavelength and spectroscopic observations. Rather than exploring a parameter space of stellar properties as we do in Geen et al. (2021), we explore a set of numerical choices that affect the evolution of the wind bubble in the simulation, and whether we can recover a consistent set of behaviours for the wind bubble. These are grid refinement on pressure gradients (ensuring that interfaces in the H II region are well resolved), masking the numerical cooling described above, the inclusion or otherwise of radiative
feedback, and the random seed used to generate the cloud initial conditions.

## 1.2 Paper overview

We run a series of simulations of wind bubbles around a self-consistently formed star, similar to the star $\theta^1$ Ori C that powers the Orion Nebula M42, in a cloud similar to Orion, which has been the focus of multiple recent multiwavelength observational studies. We run these simulations with varying parameters, described above, in order to assess whether we can reproduce the broad observed features of the Orion nebula, and hence motivate further more detailed studies.

We begin by laying out the methods used in Section 2. We then discuss the results in Section 3, focussing on the combination of wind and radiative feedback, the physics at the interfaces between the wind bubble, the photoionized H II region and the neutral cloud, and the role of randomness in setting the initial conditions. Finally, we present our conclusions in Section 4.

## 2 METHODS

We simulate a set of isolated molecular clouds with an initial magnetic and turbulent velocity field, self-gravity, and stellar feedback. The simulations are performed using the radiative magnetohydrodynamical (MHD) Eulerian Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002; Fromang, Hennebelle & Teyssier 2006; Rosdahl et al. 2013). Details of the full data reproduction package are available in Section 5.

The simulations are similar to the set-up described in Geen et al. (2021), in which we simulated a set of massive stars between 30 and 120 M$_{\odot}$ using one set of initial conditions. The main differences in the physics implemented are the full refinement of cells across steep pressure gradients, and a module that identifies and masks out cooling in cells across the wind bubble contact discontinuity. The goal of this is to explore in depth whether these numerical choices affect the resulting H II region properties, particularly when used with different combinations of feedback processes from the star. Full details are given below for completeness.

We list a summary of simulations used in this project in Table 1. We divide our suite of simulations into different labelled sets, used to identify a specific parameter space exploration. The simulations in this study form part of the PRALINES suite (Kimm et al. 2022), with the wider goal of understanding how radiative and mechanical feedback propagate to larger, (extra)galactic scales.

In the rest of this section, we discuss the numerical methods used to construct these simulations and the physical models employed.

### 2.1 Initial conditions

Our initial conditions use the same set-up as the DIFFUSE cloud in Geen et al. (2021), since it matches closely the gas density distribution found in Orion (Geen, Soler & Hennebelle 2017). We impose an initially spherically symmetric density field:

$$ 
\begin{align*}
\rho(r) &= \rho_0/(1 + (r/r_c)^2) & r < 3r_c \\
\rho(r) &= 0.01\rho_0 & 3r_c < r < 6r_c \\
\rho(r) &= \rho_{\text{ISM}} & r > 6r_c
\end{align*}
\quad (1) $$

where $\rho_0 = 1.80 \times 10^{-21} \text{ g cm}^{-3}$ (equivalent to a hydrogen number density $n_0 = 1078 \text{ cm}^{-3}$), $r_c = 2.533 \text{ pc}$, and $\rho_{\text{ISM}} = 2.24 \times 10^{-24} \text{ g cm}^{-3}$ (equivalent to a hydrogen number density $n_{\text{ISM}} = 1.34 \text{ cm}^{-3}$). The initial temperature is set to 10 K inside $6r_c$ and

<table><caption>**Table 1.** List of simulations included in this paper according to the set they are included in. See Section 2 for more detail concerning the set-up of the simulations. UV means that UV photoionization is included. Direct radiation pressure is also included in the UV runs except in the UV Only runs—see Section 2.6 for a discussion of this choice. WIND means that stellar winds are included. MASK means that the cooling mask at the contact discontinuity is included. REFINE means that refinement on pressure gradients is included (refinement on spatial coordinates, the Jeans criterion, and around the stellar source are always included). SEED refers to the initial random seed used to initialize the turbulent velocity field, described in Section 2.1. As these seeds are not ordered sequentially in a meaningful way, we give each of them an arbitrary name to reference them throughout this work.</caption><thead><tr><th>Simulation name</th><th>UV</th><th>WIND</th><th>MASK</th><th>REFINE</th><th>SEED</th></tr></thead><tbody><tr><td>Fiducial Run<sup>a</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>BELLECOUR</td></tr><tr><td colspan="6">Set WIND ONLY</td></tr><tr><td>Wind Only, Mask On</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td><td>BELLECOUR</td></tr><tr><td>Wind Only, Mask Off</td><td>✗</td><td>✓</td><td>✗</td><td>✓</td><td>BELLECOUR</td></tr><tr><td colspan="6">Set FEEDBACK</td></tr><tr><td>No Feedback</td><td>✗</td><td>✗</td><td>✗</td><td>✗</td><td>BELLECOUR</td></tr><tr><td>Wind Only</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td><td>BELLECOUR</td></tr><tr><td>UV Only</td><td>✓</td><td>✗</td><td>✗</td><td>✓</td><td>BELLECOUR</td></tr><tr><td>Wind & UV<sup>a</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>BELLECOUR</td></tr><tr><td colspan="6">Set PHYSICS</td></tr><tr><td>No Mask, No Refine</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td><td>BELLECOUR</td></tr><tr><td>No Mask, Refine</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td><td>BELLECOUR</td></tr><tr><td>Mask, No Refine</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>BELLECOUR</td></tr><tr><td>Mask, Refine<sup>a</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>BELLECOUR</td></tr><tr><td colspan="6">Set SEEDS</td></tr><tr><td>Seed AMPERE<sup>b</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>AMPERE</td></tr><tr><td>Seed BELLECOUR<sup>a</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>BELLECOUR</td></tr><tr><td>Seed CARNOT<sup>c</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>CARNOT</td></tr><tr><td>Seed DUHAMEL</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>DUHAMEL</td></tr><tr><td>Seed D'ENGHIEN</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>D'ENGHIEN</td></tr></tbody></table>

Note.

<sup>a</sup>These simulations are identical, and labelled differently in each set to clarify how they differ within the set.

<sup>b</sup>This seed was used in Geen et al. (2021), without the cooling mask or pressure gradient refinement turned on.

<sup>c</sup>This simulation is used in the 'hot champagne' analysis.

8000 K outside. The simulation volume is a cubic box of length $24r_c = 60.8$ pc. The total mass inside $6r_c$ from the centre of the box is $10^4 M_{\odot}$.

We set up an initial magnetic field, oriented in the $x$ direction. The magnetic field strength along each line of sight in this direction is given by

$$ B_x = B_{\text{max, ini}}(\Sigma_x / \Sigma_{\text{max, ini}}) \quad (2) $$

where $B_{\text{max, ini}} = 3.76 \mu\text{G}$ is the peak initial magnetic field strength, $\Sigma_x$ is the gas column density in the $x$ direction, and $\Sigma_{\text{max, ini}}$ is the peak column density. This is a nominal initial magnetic field strength corresponding to an initial Alfvén crossing time to free-fall time ratio of 0.2. The magnetic field strength evolves over time as density structures form in the cloud.

We set an initial turbulent velocity field in the cloud following the prescription given in Geen et al. (2021) at $t = 0$, with no further turbulent forcing applied. The turbulent velocity field's root-mean-sphere velocity, $v_{\text{RMS}}$ is nominally set to give a crossing time of half the free-fall time. The turbulent velocity field is created by convolving a white noise field with a Kolmogorov power law. The white noise field is initialised with a given random seed. A list of the seeds used
to generate specific white noise fields, and hence repeatable sets of initial conditions, is given in Table 1. All of the initial conditions generated by these seeds are created identically except for the value of the seed.

## 2.2 Initial evolution and refinement

We evolve the simulation using the HLLD solver for MHD (Miyoshi & Kusano 2005) as described in Fromang et al. (2006).

The average free-fall time of the gas inside $3r_c$ is 4.22 Myr. We allow the cloud to evolve without self-gravity for half of this time (2.1 Myr) in order to allow the density and turbulent velocity field to mix (see Klessen, Heitsch & Mac Low 2000). After this time we turn on self-gravity.

We use a cubic octree mesh that refines adaptively under defined conditions. Each refinement level involves a cell subdividing $2^3$ times. We set a coarse cubic grid of $2^7 = 128$ cells per side across the simulation volume. We further fully refine the grid up to level 9 inside a sphere of diameter $12r_c$.

Our maximum refinement level is 11, giving a minimum cell length of $\Delta x = 0.03$ pc. This can be reached under three conditions. First, we refine any gas that is denser than 10 times the Jeans density, $\rho_J \equiv (c_s/\Delta x)^2/G$, where $c_s$ is the sound speed in the cell. Secondly, we refine any cell more massive than 0.25 M$_\odot$. Thirdly, we refine any cell whose neighbour has a 99 per cent difference in pressure [i.e. $2|P_{\text{left}} - P_{\text{right}}|/(P_{\text{left}} + P_{\text{right}}) > 0.99$]. Simulations without this third refinement criterion are labelled ‘No Refine’, however all simulations including those labelled as such include the first two types of refinement.

## 2.3 Radiative transfer

We use fully coupled multigroup radiative transfer with the M1 method. Full details of the method as described here are given in Rosdahl et al. (2013). Photons interact with the gas via photoionization, dust absorption, and direct radiation pressure. The handling of sources of stellar radiation is discussed in Section 2.6.

We track photons in groups representing far-ultraviolet (FUV, between 11.2 and 13.6 eV) and three extreme-ultraviolet (EUV) groups capable of ionizing H I to H II (between 13.6 and 24.59 eV), He I to He II (24.59 to 54.4 eV), and He II to He III (above 54.4 eV). All photons in each group are assumed to have the same energy and cross-sections.

We store the photon density and flux for each group in each cell. We use a reduced speed of light of 0.01 c, designed to capture high-ionization front speeds. The radiation field can evolve multiple times in one (magneto)hydrodynamic time-step (termed ‘subcycling’) if the radiation time-step is shorter than the MHD time-step, each governed by the typical crossing time of radiation and gas flows in a cell.

Each gas cell tracks the ionization state of hydrogen and helium. Gas becomes ionized through photoionization and collisional ionization, and becomes neutral through recombination.

We implement radiation pressure on gas and dust as described in Rosdahl & Teyssier (2015). The local gas opacity to absorption ($\kappa_{\text{abs}}$) for all EUV groups, and scattering ($\kappa_{\text{sc}}$) for all groups, is given by $10^3 Z/Z_{\text{ref}} \text{ cm}^2 \text{ g}^{-1}$, where Z is the metallicity of the gas (=0.014 in all runs, matching the Solar metallicity used in Ekström et al. 2012) and $Z_{\text{ref}} = 0.02$. We do not track the creation and destruction of dust as a separate fluid, as in, e.g. Lebreuilly, Commerçon & Laibe (2019). Instead, we approximate the effect of dust by including an additional

absorption and scattering opacity, which are given by

$$ \{\kappa_{\text{abs}}, \kappa_{\text{sc}}\}_{\text{dust}} = \{\kappa_{\text{abs}}, \kappa_{\text{sc}}\}_{\text{gas}} \left(\frac{T_R}{10 \text{ K}}\right)^2 \times \exp\left(\frac{-T_R}{10^3 \text{ K}}\right) (1 - x_{\text{HII}})\quad(3) $$

where ‘abs’ refers to absorption and ‘sc’ refers to scattering. The term $10^3$ K is designed to mimic dust destruction by sputtering above this temperature. $T_R = (E_R/a)^{1/4}$ is the local radiation temperature, where $E_R$ is the energy of the photons in each group in the cell and $a \equiv 7.566 \times 10^{-15} \text{ erg cm}^{-3} \text{ K}^{-4}$ is the radiation constant, $x_{\text{HII}}$ is the local hydrogen ionization fraction, and $\rho$ is the local gas density. A full dust model coupling to radiation pressure is left to future work, though based on our previous models in Geen et al. (2021) and Geen & de Koter (2022), we do not expect this to change our results significantly for the systems studied here.

## 2.4 Radiative cooling

In addition to radiative transfer of high-energy photons, a certain fraction of the gas thermal energy is considered lost to low-energy radiation in each time-step. Hydrogen and helium follow the cooling and heating functions described in Rosdahl et al. (2013). Metal cooling is accounted for with separate functions following the neutral and ionized hydrogen fractions in each cell. Metal cooling for the neutral fraction of a cell is given by the model of Audit & Hennebelle (2005) which tracks cooling of various elements in the ISM. Metal cooling in the ionized fraction uses a piecewise fit to the cooling curve given in Ferland (2003), with a fit to Sutherland & Dopita (1993) above $10^4$ K. For all of our simulations, the metal fraction is $Z = 0.014$, i.e. Solar metallicity. We ignore metal enrichment by stellar mass-loss as we expect this to be small compared with the initial cloud metallicity.

Shocks moving on Eulerian grid codes are well captured when adiabatic, but when cooling is included, artificial numerical mixing occurs as the shock crosses a cell and a contact discontinuity exists between the hot, diffuse gas post-shock and the cold, dense gas pre-shock. This creates an artificially warm, dense cell that leads to spurious overcooling. This is discussed in some detail in Fierlinger et al. (2016) and Gentry et al. (2017), the latter for the case of supernova remnants rather than main-sequence stellar winds.

In order to mitigate this effect, we invoke a very simple fix that masks cells in which neighbouring cells (identified in a 3x3x3 cube centred on a given cell) contain both cells above $10^6$ K and cells below $10^6$ K, and turns off cooling in this mask. We use $T = 10^6$ K as a reliable heuristic for whether gas is in the wind-shocked bubble or outside it. We discuss the implications of our simplified choice later in the paper. Our simulations are considered to have this fix by default, or are labelled ‘Mask’. Simulations without this contact discontinuity mask are labelled ‘No Mask’.

## 2.5 Sinks and star formation

In our simulations, we follow star formation using sink particles that accrete from the surrounding gas. We first identify if cells reach 10 per cent of the Jeans density at the highest level of refinement. We then identify peaks amongst these cells using the watershed method described in Bleuler et al. (2014). If the peak is denser than the Jeans density at the highest level of refinement, we form a sink particle. Every time-step, this sink accretes 90 per cent of the mass
in the clump above the Jeans density. For a complete description, see Bleuler & Teyssier (2014).

We use these sink particles to track material that collapses to scales smaller than those we can track with the AMR grid, and that become the sites of star formation. We consider a small cluster to have formed once $120 M_\odot$ has accreted in total on to sink particles. At this point, we create a stellar object and assign it to the most massive sink at the time of formation. This object is not in itself an *N*-body particle and rather sits on the position of this sink at all times.

We track the age of the object from the moment of its creation and consider it to be a star of $35 M_\odot$, similar to $\theta^1$ Ori C, the source of the Orion Nebula, M42 (Kraus et al. 2009; Balega et al. 2014). This mass is used by the stellar evolution track to extract accurate wind and radiation emission rates for a given star in each time-step, and does not have any influence on the gravity in the simulation aside from the mass of the sink it is a part of. We do not allow any other stellar objects to form that produce feedback, in order to isolate the feedback from one star, although we allow sinks to continue accreting material, modelling the ongoing formation of lower mass stars in the cloud. This is representative of the situation in M42. Future work will include feedback from multiple self-consistently forming stars, as in Geen et al. (2018).

## 2.6 Stellar feedback

The age and mass of the stellar object is used to track the radiation and winds emitted by the star using a stellar feedback module. A full description of the feedback module with figures showing the tracks is given in Geen et al. (2021). We summarize the module here.

We do not consider the pre-main-sequence phase of the $35 M_\odot$ star since we do not track the scales relevant to protostar formation. The star is considered to immediately arrive at the zero-age main sequence once the stellar object is created.

Our star emits energy as radiation in the four modelled groups (one FUV and three EUV) and experiences mass, momentum, and energy loss from stellar winds. Our simulations end before the star should go supernova (at an age of 6.8 Myr in our model).

We base our stellar evolution tracks on the Geneva Model (Ekström et al. 2012; Georgy et al. 2012). We use the rotating tracks, i.e. stars rotating at 40 per cent of their critical velocity. Stellar spectra are extracted from STARBURST99 (Leitherer et al. 2014) with the Geneva models as inputs.

Our wind model is similar to that in Gatto et al. (2017). Mass-loss rates are taken from the Geneva model, with a corrected wind terminal velocity calculated from the escape velocity as in Vink et al. (2011). A full description of the model used is given in Geen et al. (2021).

Radiation is deposited on to the grid in a single cell where the host sink of the star is located. This radiation then propagates outwards following the M1 method.

Winds are deposited on to the grid as mass, momentum, and energy in a thick shell between 3.5 and 5 cells in radius at the maximum level of refinement around the host sink of the star. We do this to ensure that the initial state of the feedback region captures swept-up material by the unresolved compact H II region around the star, which in turn affects how radiation from the star is allowed to propagate into the surrounding medium at $t \sim 0$. We inject all of the energy from the wind as kinetic energy.

Initially, the background density of the cells that the wind is injected into is high. This means that the flow thermalizes, since the momentum added by the wind bubble is negligible. As the wind sweeps out a low-density bubble, the mass injected by the wind

accounts for the majority of the mass inside the injection radius, and so the flow becomes kinetic. We thus begin to capture the free-streaming radius described in Weaver et al. (1977).

We perform runs with both winds and radiation, winds only and radiation only (see Table 1). Runs with radiation only ('UV Only') omit direct radiation pressure on gas and dust. This is because of a physical effect where radiation pressure creates a central empty cavity around the star (see e.g. Draine 2011), leading to an artificially short time-step in this region as the simulation attempts to balance the pressure in the region. This causes the simulation to slow down considerably. This effect does not occur when stellar winds of any kind are included, as the ambient volume around the star becomes filled with free-streaming wind material, and so in runs with winds, radiation pressure is also included. In this way we consider the effect to be unphysical, but include the UV Only run none the less for the purposes of comparison with the runs including winds. In any case, we note that radiation pressure has a negligible effect on our results, in agreement with our previous simulation work (Geen et al. 2021) and analysis of observed H II regions of a similar size (Olivier et al. 2021).

Additionally, for the UV Only run, we disable refinement on pressure gradients, since the simulation attempts to refine over a large volume once the 'champagne flow' phase begins. This causes the simulation to request excessive memory and CPU resources.

## 3 RESULTS

In this section, we describe the results of combining wind and radiative feedback when using the numerical recipes implemented in this work. We first present the results of combining winds and radiative feedback using a fiducial set of these recipes, finding that winds are capable of trapping radiative feedback as described in Geen & de Koter (2022). We then study the effect of correcting for numerical overcooling at the wind bubble contact discontinuity, studying in particular the role of turbulent mixing when this correction is applied. We then study how the initial random seed used to generate the initial conditions plays a strong role in the evolution of the wind bubble.

### 3.1 Combining winds and radiative feedback

We now compare the influence of stellar winds and UV radiation feedback in shaping H II regions. To do this, we analyse the FEEDBACK set of simulations, which explores the effect of individual feedback processes included in the simulations. Specifically, we review the combined impact of UV radiation and winds on the cloud.

Fig. 1 shows a sequence of images of each simulation highlighting the presence of cold, neutral gas (purple); warm, photoionized gas (orange); and hot, wind-shocked gas (light yellow), with contours outlining the wind bubble (cyan) and photoionized gas (red) separating the phases of the H II region around the star (red star icon), and with time running from top to bottom with snapshots at 0.1, 0.2, and 0.3 Myr.

The notable result of this comparison is that combining winds and radiation from the star gives it a much weaker impact than simply including ionizing radiation. This is also seen in the mean spherical-equivalent radial expansion of the region (Fig. 2, where the radius is calculated as $(3V_{\text{II}}/4\pi)^{1/3}$, where $V_{\text{II}}$ is the H II region volume including both photoionized and hot, wind-shocked gas above $10^5$ K) and momentum in outflows as a result of feedback (Fig. 3, including only components of gas flows moving radially away from the star). In the UV Only run, there is a very rapid
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Sequence of normalized emission maps showing the evolution of the nebula around the star in the FEEDBACK set, where we vary whether winds or UV photoionization are included. Each image shows a projection of emission from collisional excitation from cold gas ($\rho^2$ for cells where $T < 1000$ K), hot gas ($\rho^2$ for cells where $T > 106$ K), and recombination emission ($\rho^2 x_{HII}$, ignoring cells where $x_{HII} < 10^{-4}$). A cyan contour encloses pixels with gas above $10^6$ K in them. A red contour encloses pixels containing cells with a hydrogen ionization fraction $x_{HII} > 0.1$, i.e. outlining where the extent of the HII region lies. Some gas with a lower ionization fraction is visible outside this region, as the ionization front is not a sharp discontinuity in runs that exhibit champagne flows. The columns from left to right show runs with winds and UV photoionization, winds only, UV photoionization only and no feedback. The rows from top to bottom show the outputs in each simulation at a stellar age of 0.1, 0.2, and 0.3 Myr. The apparent drop in recombination line emission at 0.3 Myr in the third column is due to normalization effects and does not reflect a drop in photoionization.</figcaption></figure>

expansion of the H II region at a stellar age of 40 kyr, corresponding to a 'champagne' flow (Tenorio-Tagle 1979; Whitworth 1979), in which the dense gas around the star cannot absorb all of the ionizing radiation any more, causing the ionization front to move outwards supersonically without significant gas flows. This can be seen in the UV Only run in Fig. 2, where the mean radius of the H II region expands rapidly after ~0.04 Myr, before slowing at 0.1 Myr where parts of the region reach the edge of the box. By comparison, the momentum of the gas in the UV Only run in Fig. 3 advances more slowly with no real sign of a slowing at 0.1 Myr as the photoionized gas begins to respond hydrodynamically (see Franco et al. 1990).

The trapping of ionizing radiation by a wind-blown shell in a power-law density field is described analytically in Geen & de Koter (2022), where the presence of a hot wind bubble acts to increase the pressure inside the H II region, creating a denser shell around the region and hence trapping a large quantity of radiation. It is

worth noting that this is created by swept-up material rather than the mass directly from the star, which is typically much smaller than the mass of the circumstellar medium in a dense cloud. Even when the champagne flow phase occurs in simulations where winds are included (around 0.17 Myr in the Wind & UV simulation, seen in the increase in expansion rate in Fig. 2), the ionization front grows more slowly when winds are included. This suggests that a significantly reduced amount of radiation is able to escape even once the champagne phase begins.

The presence of a dense shell can be seen by examining Fig. 4, where we plot slices through the position of the star showing gas density. In the UV Only run, the gas is rapidly photoionized as described in Franco et al. (1990), and as such there is no shell around parts of the H II region that form a Champagne flow. Rather, the photoionized cloud expands thermally. By contrast, although the champagne phase has begun in the Wind & UV run in this figure, there is still a noticeable dense shell around the diffuse
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Mean spherical-equivalent H II region radius $R_{HII}$ (including wind and photoionization-heated gas, where the region volume $V_{HII} = \frac{1}{2}\pi R_{HII}^3$) as a function of time in each simulation in the FEEDBACK set. Winds act to constrain runaway photoionization of the whole cloud by trapping the radiation in the shell around the wind bubble.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Momentum in flows directed away from the position of the star as a function of time in each simulation in the FEEDBACK set. As with the radial evolution in Fig. 2, winds drastically reduces the impact of radiation feedback and hence the total effect of feedback on the cloud.</figcaption></figure>

wind bubble that acts to absorb significant quantities of ionizing radiation.

There is some very limited leakage through parts of the shell in the Wind & UV simulation. This can be seen when comparing the position of the red contour and the edge of the orange region in Fig. 1, which use cutoff thresholds of $x_{HII} < 10^{-1}$ and $x_{HII} < 10^{-4}$, respectively. Despite this, nearly all of the ionizing radiation is trapped, compared with the rapid champagne flow seen in the UV Only runs.

## 3.2 The physics at the H II region interface

A key component in H II regions is the interface between different phases, in particular the contact discontinuity between the wind bubble and the gas around it, which is comprised of either a dense, partially, or a fully photoionized shell of material.

To this end, we study two recipes designed to improve our simulations' ability to capture the physics of this region compared with Geen et al. (2021), namely our masking out of cooling around the contact discontinuity between the wind bubble and the gas outside it (simulations labelled 'Mask' when included, otherwise 'No Mask'), and the use of grid refinement to better resolve strong pressure gradients in our simulations (labelled 'Refine' when pressure gradient refinement is included, or 'No Refine' otherwise).

### 3.2.1 Wind-only comparison

We first compare results in a test set-up without ionizing radiation from the star, in order to simplify the comparison between runs with and without the wind contact discontinuity cooling mask.

In Fig. 5, we compare the radial extent of the wind bubble as a function of time and the energy from stellar winds retained in the wind bubble with and without the cooling mask. The spherical-equivalent radius (calculated assuming the volume of hot gas $V_w$ above $10^5$ K is a sphere, i.e. $(3V_w/4\pi)^{1/3}$) at 300 kyr is typically half as large when the cooling mask is applied versus without, although the stellar wind energy retained is higher when the cooling mask is applied. This slower radial expansion rate is somewhat counterintuitive, as one might expect more retained energy in the wind bubble to equate to a faster expansion rate. Indeed, this effect is reversed when we include photoionization. We discuss the reasons for this below. Similarly, analysing the simulations further, the momentum added to the cloud by the wind is roughly doubled when the cooling mask is applied, although the pressure in both regions is very similar.

Another major difference between the behaviour of the wind bubbles with and without the cooling mask is seen by analysing a slice through the bubble in Fig. 6. Though certain parts of the wind bubble have extended further when the cooling mask is applied, there is considerably more substructure inside the bubble itself, seen as the yellow ridge inside the red bubble. This in turn reduces the spherical-equivalent mean radius, which does not include any dense, neutral substructures inside the wind bubble, only cells heated to above $10^5$ K.

As seen in Fig. 6, runs where we suppress cooling at the contact discontinuity exhibit increased macroscopic mixing with the cloud material around the wind bubble. This is likely because energy losses via the turbulent cascade into thermal energy are suppressed. However, environment around the star is turbulent with large density gradients, complicating a simple linear analysis, and it is equally possible that a small change in the simulation leads to diverging behaviour as the turbulent system evolves, and differences over time become the result of non-linearity inherent in the fluid equations. Overall, applying the cooling mask thus does not prevent the wind bubble cooling entirely—rather, it sets a lower bound on unresolved 'subgrid' cooling, where not applying the mask gives an upper bound that assumes that any unphysical numerical cooling is equivalent to that from any unresolved mixing or conduction processes.

To summarize, the cooling mask does increase the amount of energy retained in the wind bubble. However, due to interactions with the turbulent cloud, this does not necessarily lead to a faster radial expansion of the wind bubble, and the mask instead enhances the role of macroscopic turbulent mixing on the scale of a few cells as opposed to subgrid mixing, which is suppressed by the mask. The downside is that the mask suppresses both physical, unresolved cooling (Rosen et al. 2014; Lancaster et al. 2021a) and non-physical numerical overcooling (Fierlinger et al. 2016; Gentry et al.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Sequence of slices through the position of the star showing gas density in the FEEDBACK set, where we vary whether winds or UV photoionizations are included. The columns show the same simulations as Fig. 1. Runs with winds contain dense shells around the wind bubble, which absorbs UV radiation from the star.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of runs in the WIND ONLY set where the cooling mask around the wind bubble contact discontinuity is either turned on or off. The left-hand panel shows the mean spherical-equivalent radius of the wind bubble $R_w$ (where the wind bubble volume $V_w = 4/3\pi R_w^3$, and $V_w$ is the volume of gas in cells above $10^5$ K). The right-hand panel shows the fraction of energy from stellar winds retained in the wind bubble (kinetic + thermal). Despite retaining less energy, the wind bubble in the run without the cooling mask has a larger volume. We provide an explanation for this difference in Section 3.2.1.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Slices in temperature through the position of the star in the WIND ONLY set at a stellar age of 0.3 Myr, comparing the effects of the cooling mask around the wind bubble contact discontinuity with only stellar winds included (no UV photoionization). The left image with the cooling mask applied shows more internal substructure than the run shown on the right without the cooling mask. A translucent blue contour outlines the masked cells in the left image. The wind bubble is considered to be all gas cells above $10^5$ K, which does not include any substructure inside the wind bubble at a lower temperature.</figcaption></figure>

2017), although this is unavoidable without fully resolving the wind bubble contact discontinuity, which in turn becomes computationally expensive.

### 3.2.2 Comparison including photoionization and refinement

We now consider a comparison of the effect of the cooling mask including photoionization feedback and refinement on pressure gradients, i.e. at the edge of the H II region. In Fig. 7, we plot the energy retained in the wind bubble with and without the cooling mask and refinement on pressure gradients (left-hand panel) as well as with and without photoionization (right-hand panel). We define the wind bubble as including all cells with a temperature above $10^5$ K, where energy is the total thermal and kinetic energy.

Unsurprisingly, the cooling mask causes the wind bubble to retain more energy, although the difference shrinks if refinement on pressure gradients is included, suggesting that the results should converge if increased resolution across the contact discontinuity is achieved (with the corresponding increase in computational requirements).

The inclusion of photoionization also increases the energy retained in the wind bubble. Photoionization heats the gas around the wind bubble to ~$10^4$ K, which both reduces the temperature difference between this gas and the hot wind bubble and causes it to thermally expand to a lower density. This means the wind bubble mixes with a relatively warmer, more diffuse medium than the cold, dense neutral shell, leading to reduced cooling in the wind bubble.

We plot slices through the wind bubble, photoionized H II region and surrounding cold cloud in Fig. 8. The wind bubble is highly aspherical, forming a chimney of hot, fast-flowing gas that pushes through low-density channels in the cloud. This chimney becomes a wider plume of hot gas as it encounters the interface between the H II

region and the unshocked gas outside. As refinement is turned on, this plume becomes more structured with fluid instabilities, in some cases breaking up into isolated hot bubbles. Warm photoionized gas fills in the volume between the star and the edge of the H II region. In the ‘Mask, Refine’ simulation, the photoionized gas does not extend further than the wind bubble, although in all cases a full champagne flow is prevented.

We study the radial evolution of the wind bubble and H II region in more detail in Fig. 9. We plot here the maximum radius, i.e. the distance of the furthest cell in the wind bubble or H II region from the star. Measuring the radius by using the H II region volume $V_{HII}$, where $R_{HII} = (3V_{HII}/4\pi)^{1/3}$, gives a similar result. However, in certain simulations, e.g. ‘No Mask, No Refine’, the wind bubble is less extended but has a larger volume. We use this comparison here since the wind bubble is highly aspherical and its most observable characteristic is its extent.

Masking cooling at the contact discontinuity allows for a more extended H II region than for runs without the mask. Unlike the Wind Only simulations, UV radiation photoionizes any substructure trapped by the wind bubble. As the H II region expands further, the wind bubble also expands into the underdensity created by the H II region, expanding further in runs with the cooling mask applied. Increased refinement either gives similar or a larger radius. However, in the case of the H II region radius, the ‘Mask, No Refine’ run reaches the edge of the cloud and begins to break out in one direction, despite having a similar wind bubble extent to the run with refinement.

The complex, non-linear gas dynamics displayed makes it hard to demonstrate a consistent linear response to the different physical models used. However, all models display a similar behaviour, with increased refinement leading to the wind bubble displaying more turbulent substructures. Including a cooling mask also leads to a generally wider plume as well as a longer chimney.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The energy in the wind bubble as a fraction of the total energy injected as stellar winds, as a function of time in each simulation in the PHYSICS (left-hand panel) and FEEDBACK (right-hand panel) sets. Including UV photoionization reduces wind cooling, as does including a cooling mask around the wind bubble contact discontinuity, though including refinement on pressure gradients causes some modest convergence between the results with and without the cooling mask.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Slices through the position of the stellar source showing gas temperature in the PHYSICS set at a stellar age of 0.3 Myr. Gas above $10^5$ K is wind-shocked gas, gas at ~$10^4$ K is photoionized and colder gas is neutral cloud material. A translucent blue contour outlines the masked cells in the runs that include a cooling mask.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Maximum extent away from the star of the H II region (left-hand panel) and the wind bubble $R_w$ (right-hand panel) as a function of time in each simulation in the PHYSICS set. Winds act to constrain runaway photoionization of the whole cloud by trapping the radiation in the shell around the wind bubble.</figcaption></figure>
### 3.2.3 The role of turbulence

In the classical Weaver et al. (1977) model of wind bubbles, the majority of cooling takes place via evaporation of material from the dense shell into the hot wind bubble (Mac Low & McCray 1988), which is slow even if the wind bubble is embedded in a dense molecular cloud environment (Kruijssen et al. 2019; Geen & de Koter 2022). However, wind bubbles evolve into complex structured environments, and develop fluid instabilities that are impossible to simulate directly in simple 1D models, and which can enhance cooling via turbulent mixing (Rosen et al. 2014). Lancaster et al. (2021a) present a model of the evolution of wind bubbles with a fractal surface that assumes efficient cooling via turbulent mixing with the medium around the wind bubble. The authors describe a number of diagnostics to track where this efficient cooling should occur, assuming a globally uniform cloud density field. This compares well to simulations of wind bubbles representing stellar populations in molecular clouds with idealized fixed sources (Lancaster et al. 2021b) and with sources forming on sink particles (Lancaster et al. 2021c).

In our work, we focus on a single stellar source (versus multiple sources in Lancaster et al. 2021c) but also include a cooling mask designed to remove spurious numerical cooling. However, while this mask does remove artificial cooling, it also prevents all cooling at the wind bubble’s contact discontinuity on the scale of 1–2 grid cells, including possible cooling from unresolved turbulent mixing. While Tan et al. (2021) argue that the largest eddy scale dominates turbulent mixing, Lancaster et al. (2021b) find that in simulations of wind bubbles in molecular clouds, the relevant cooling scale is difficult to resolve. None the less, they argue that it is possible to determine if efficient turbulent mixing can occur if the turbulent velocity of gas at the wind bubble interface allows faster diffusion of energy than can be injected by the bulk velocity of material from the star entering this interface.

Turbulence is a chaotic phenomenon characterized by fluid eddies at multiple scales which diffuse energy from large to small scales. There are a number of ways to characterize this in simulations. The simplest method for cases where there is one source is to assume that the source will generally produce radially aligned outflows, and thus measuring non-radial flows gives an estimate of the levels of turbulence in the outflows as opposed to laminar radial flows. However, this is not possible in our case since even the broadly laminar chimneys flowing from the star become non-radial due to following decreasing density contours in the cloud. Another is to measure the velocity dispersion, which is a statistical measure of how different fluid velocities are in surrounding computational elements. Finally, one can measure the vorticity, i.e. curl of the velocity field ($\nabla \times \mathbf{v}$), which gives the turnover frequency of fluid eddies at a given scale.

In Fig. 10, we show slices through the wind bubble in the Fiducial Run to analyse the impact of turbulence on the wind bubble using vorticity as the measure of turbulence. We sample across a nominal length scale $\ell$ of four simulation cells at the highest level of refinement ($\ell = 0.12$ pc) in order to smooth out noise on a cell-by-cell level. In this figure we plot the gas density, the bulk gas speed, the eddy turnover time-scale ($1/\nabla \times \mathbf{v}$), and the ratio between the characteristic eddy speed $\ell \nabla \times \mathbf{v}$ and the bulk gas speed in each cell. Some grid artifacts can be seen due to the lower cell resolution outside the wind bubble where we do not fully refine the grid.

The eddy turnover time-scale is longest (i.e. darker green in the bottom-left panel of Fig. 10) in the chimney of the wind bubble, as well as the extended plume at the edge of the bubble, while the

bulk gas velocity is highest in the chimney. There are regions with strong eddies on the shearing interface between the chimney and the photoionized gas outside, as well as on the complex surface where the plume interacts with the surrounding material. There are also strong eddies in the dense shell itself, although this can also be a measure of the velocity difference between the expanding wind bubble and the undisturbed material outside it.

The bulk speed of the gas in the wind bubble, and the expansion rate of the wind bubble itself, is significantly higher than the speed of the gas outside the wind bubble. This suggests that much of the turbulence at the wind bubble interface is shear from fast-moving chimneys of gas expanding into regions of lower density.

Ultimately, the ‘correct’ modelling of the effect of turbulence on the scale of a few cells around the wind bubble interface remains a difficult task, and the true fraction of energy lost to cooling will lie somewhere between the models with and without the cooling mask. The ‘true’ cooling rates and hence efficiency of wind bubbles as feedback sources thus remains a subject of difficulty, despite promising work in recent years. Some of these difficulties can be mitigated by using Lagrangian (i.e. moving mesh) methods versus static Eulerian grids. However, spatial resolution requirements remain a tough limitation to overcome, particularly given the already difficult computational requirements in capturing the time-steps needed to trace flows travelling at thousands of kilometre per second.

### 3.2.4 How do physical wind bubbles expand rapidly?

Analysis of the Orion Nebula M42 by Pabst et al. (2019) finds that the wind bubble expands at a speed consistent with an adiabatic wind bubble, or ~13 km s⁻¹. Conversely, following our simulations, we do not expect the wind bubble interface to be adiabatic as the wind bubble interface exhibits large quantities of turbulence on all scales, which we expect to lead to enhanced cooling. This in turn should lead the wind bubble to expand more slowly, as it has a lower amount of energy available to drive its expansion.

A cooling mask does help to mitigate this by removing sources of overcooling from artificial numerical mixing. However, it does not completely remove sources of cooling to an adiabatic level. When comparing simulations with and without a cooling mask in Fig. 7, we find a factor of a few differences in cooling when applying the cooling mask versus not in cases where we also include refinement on pressure gradients. Applying the cooling mask in these cases allows the wind bubble to retain ~5–7 per cent of the energy injected by the star in the gas (including both thermal and kinetic components). This is compared with 1–2 per cent when no mask is implied. This is significantly lower energy retention than the 33 per cent of wind energy retained (after work done on the surrounding medium is included) found in the adiabatic model of Geen & de Koter (2022) assuming a singular isothermal sphere power law density field, or 45 per cent retention when using the Weaver et al. (1977) model in a uniform medium.

One explanation is that the opening angle of the wind chimney is significantly smaller than the total opening angle of the wind bubble. This creates a wind bubble that has a smaller volume than a spherical wind bubble with the same extent away from the star, and may explain why the wind bubble is able to reach larger radii. We plot this in Fig. 11, showing the solid angle subtended by the wind bubble (in red) and the total H II region (in black) at each radius away from the star for a set of stellar ages. As the wind bubble evolves, the region around the star remains fully wind-driven (partly caused by the free-streaming phase of the wind), before forming chimneys that take up under $\pi$ of the full $4\pi$ of the solid angle. This is an average across
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Slices through the Fiducial Run at a stellar age of 0.3 Myr, depicting the turbulent properties of the gas. From top-left to bottom-right panels, we show gas density, bulk gas speed, the inverse of gas vorticity [i.e. the eddy turnover time-scale at a scale of 0.12 pc (4 times the smallest cell size)], and the ratio of characteristic eddy speed ($\ell\nabla \times \mathbf{v}$) to bulk gas speed. By comparing the eddy speed with the bulk gas speed we aim to identify where the gas flows are mostly laminar at a 0.12 pc scale, and which flows are highly turbulent. The interior of the wind bubble is indeed highly turbulent, with mixing evident between the wind bubble and the rest of the H II region. The circle around the star is the free-streaming radius of the wind. The white dots in the bottom-right panel are where the lower grid resolution creates points where a cell is compared with itself, resulting in zero eddy speed.</figcaption></figure>

all directions from the star including the smaller ‘back-blast’ seen in the bottom-left of each panel in Fig. 10 and earlier figures showing the same simulation.

Additionally, by examining both Figs 8 and 10, we see that the picture in our 3D turbulent cloud is very different to that in spherically symmetric 1D models where a spherical wind bubble sits inside a photoionized shell. While the hot wind bubble follows regions of low density, radiation from the star must travel in straight lines, and so fills a larger solid angle around the star than the wind bubble, which travels in a narrow chimney through the cloud before expanding into a wide plume. This implies that some H II regions appearing to have large wind bubble volume filling fractions may in fact exhibit a

large plume near the edge of the H II region with a narrow embedded chimney. It also suggests that the wind bubble may expand faster than in a purely spherical model if it is allowed to expand over narrow solid angles.

## 3.3 The role of chaos and initial conditions in shaping wind bubbles

Geen et al. (2018) demonstrate how the role of chaos plays a significant role in the outcomes of feedback in star-forming clouds. The use of different random seeds to initialize the turbulent velocity field in the cloud, which in turn sets the distribution of density
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Solid angle subtended by the wind bubble (red) and the overall H II region (black, includes both photoionized and wind-shocked gas) as a function of radius in the Fiducial Run. Results taken at stellar ages of 0.1 Myr are shown as a solid line, 0.2 Myr as a dashed line and at 0.3 Myr as a dotted line. We project 1000 rays away from the star, and measure the fraction of the full solid angle around the star (4π) that contains either wind-shocked gas (T > 10⁵ K) or photoionized gas (xHII > 0.1) at each radius from the star. Over time, a wide plume develops that has a larger solid angle than the chimney that feeds it. In many cases, photoionized gas back-fills the solid angles subtended by the plume at smaller radii, since the ionizing radiation travels directly outwards rather than seeking pressure gradients as the wind bubble does.</figcaption></figure>

peaks that form in the cloud, causes feedback to behave in different ways.

In Fig. 12, we plot the results of the SEEDS simulation set where we vary the random seed used to generate the initial conditions, as described in Section 2.1, and explore what effect this has on the resulting wind bubble and H II region. We form the stellar source self-consistently following sink particles, and use the same stellar tracks (a 35 M☉ star) in each simulation.

There is a large diversity of behaviours of the H II regions purely by varying the initial seed. In some cases the photoionized region is effectively contained near the wind bubble, even after a few hundred kyr, but in other cases there is still a rapid champagne flow out of the cloud. Similarly, some simulations have a prominent wind bubble, while others have a small wind bubble. The wind bubble in Seed CARNOT even grows and disappears, before appearing again 0.1 Myr later, which we discuss below.

The emergence of a champagne flow breakout or otherwise in the H II region for each seed is visible in Fig. 13, which is displayed as a sudden increase in radial expansion rate much faster than the gas sound speed of ~10 km s⁻¹ ≃ 1 pc/0.1 Myr. Seeds BELLECOUR and D'ENGIHEN never display a strong champagne breakout. By comparison, in Seed AMPERE, the seed used in Geen et al. (2021), the breakout begins almost immediately (see also the top row of Fig. 12).

A key aspect of whether a champagne flow forms is whether the star forms on the edge of the cloud or not. The trapping of ionizing radiation by the wind bubble described in Geen & de Koter (2022) requires a supply of swept-up cloud material to maintain the trapping. If the wind bubble encounters a density drop at the edge of the cloud, the shell can 'burst', leading to a champagne flow as described in Tenorio-Tagle (1979) and Whitworth (1979) (see also Comeron 1997, for models including winds).

In Fig. 14, we plot the energy emitted by the star as stellar winds retained in the wind bubble. An adiabatic wind should retain 45 per cent of its energy in the hot wind bubble for a uniform cloud (Weaver et al. 1977) or 33 per cent in a $\rho \propto r^{-2}$ power-law density field (Geen & de Koter 2022). Three of the seeds retain 10–20 per cent of their energy at some point in their evolution, which approaches adiabatic assuming the power-law density field around the star, with the retained energy slowly decreasing over time. We discuss the cases of Seeds CARNOT and DUHAMEL in the following two sections, and discuss their behaviour. Seed CARNOT undergoes a hydrodynamically unstable process, which drastically increases mixing of the wind bubble. In the case of Seed DUHAMEL, the low-energy retention is likely to be due to resolution limits rather than a physical lack of energy retention. None the less, despite the differences in the wind bubble evolution, three out of five of the simulations display evolution in their energetics after the first 0.1 Myr.

### 3.3.1 ‘Hot’ champagne

Seeds AMPERE, BELLECOUR, and D'ENGIHEN display large-developed wind bubbles. By contrast, Seed CARNOT loses its wind bubble around 0.3 Myr, while Seed DUHAMEL never develops a strong wind bubble at all. This is not a function of strong accretion on to the sink in these cases—although Seed DUHAMEL's star does continue to be embedded inside a large filament up to 0.3 Myr, accretion rates are similar in all runs, except for Seed BELLECOUR which is higher despite having a well-developed wind bubble. Similarly, the ram pressure of flows in the neutral gas around the sinks at the time of formation is similar.

In the case of Seed CARNOT, the disappearance of the wind bubble is linked to the H II regions entering a breakout champagne flow phase. Fig. 15 displays the progression of the H II region at the point a champagne flow begins. The depressurization of the photoionized region allows the wind bubble to expand more quickly and mix with the photoionized gas, causing a slightly hotter H II region with the same temperature as the now cooled wind bubble. As the wind is injected on top of the background gas, turbulent mixing near the star can cause the wind bubble to disappear temporarily even if the wind energy is still being injected.

We refer to this as a ‘hot champagne’ flow, since it involves the efficient mixing of hot (> 10⁶ K) wind-shocked and warm (~10⁴ K) photoionized gas. The phenomenon is transient and lasts for less than 0.1 Myr. Once the initial rapid expansion phase has occurred, the wind bubble gradually reappears (see the lower panel of Fig. 15). This is combined with other effects such as photoionization in the champagne flow being temporarily disrupted. As the phenomenon is linked to out-of-dynamical-equilibrium behaviour in the H II region following its encountering a discontinuity, this effect is also temporary.

If found in physical H II regions that have just begun to undergo the champagne phase, this behaviour is likely to pose a problem for observational studies that map the temperature of the photoionized region to the spectrum of the star, assuming that residual energy from the wind bubble remains in the photoionized gas for an appreciable amount of time. An H II region without a visible wind bubble may be an indication that such a phenomenon has occurred. More simulation work is needed to establish how common this phenomenon is, as well as whether it can be found in other cosmic conditions such as at lower metallicities.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** As for the Wind & UV simulation shown in Fig. 1 but with the SEEDS set showing the effect of different random seeds in the initial turbulent field on the resulting nebula.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Maximum extent of the H II region away from the star (including wind and photoionization-heated gas) as a function of time in each simulation in the SEEDS set. Winds act to constrain runaway photoionization of the whole cloud by trapping the vast majority of the radiation in the shell around the low-density wind bubble, as described in Section 3.1. Fast radial expansion in the AMPERE, CARNOT, and DUHAMEL runs is caused by the H II region encountering the edge of the cloud and expanding into the low-density medium outside as a 'champagne' flow.</figcaption></figure>

### 3.3.2 Shell structure

Comparisons with full synthetic observations are left to future work, given the complexity of matching the full velocity structure of line

<figure><img src="image_15.png" /><figcaption>**Figure 14.** The energy in the wind bubble as a fraction of the total energy injected as stellar winds, as a function of time in each simulation in the SEEDS set. Strong wind bubbles retain up to 10 per cent of the energy from stellar winds, while weak bubbles retain around 1 per cent, with some results moving between the two limits. (For reference, an adiabatic wind bubble evolving in a singular isothermal sphere retains 33 per cent of its energy inside the bubble, see Geen & de Koter 2022). The loss of the wind bubble in the Seed CARNOT run due to the emergence of a 'hot champagne' flow is visible as the sharp decline around 200 kyr.</figcaption></figure>

emission from tracers such as C⁺ (e.g. Schneider et al. 2020). However, we approximate the observable features of the neutral shell around the H II regions by measuring the column density of gas in cells with a total speed faster than 5 km s⁻¹, in order to isolate
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Slices through the star in the Seed CARNOT run, depicting the ‘hot champagne’ outflow and disappearance of the wind bubble. Fields from left to right show gas density, temperature, and photoionized fraction (ignoring collisional ionization, which occurs at higher temperatures in the wind bubble). Time evolves from top to bottom, showing snapshots as photoionized gas breaks out into a champagne flow, where the wind bubble mixes efficiently with the photoionized gas to the point where no hot (>10⁶ K) gas remains visible. The phenomenon is transient and lasts for less than 0.1 Myr.</figcaption></figure>
material in the dense, expanding shell. This is also where tracers such as C⁺ would be found, while also removing ambient gas in the cloud that is typically moving at lower velocities. We plot this quantity in Fig. 16.

In the top image (x-projection, used in the other figures in this paper), large-scale photoionized champagne flows, evidenced by the lack of a shell over a large solid angle, are visible in Seed AMPERE, but less so in other images. By contrast, in the bottom image (z-projection), the champagne flow is more visible in seeds DUHAMEL and D'ENGHIEN. Analysis of the full velocity cube may be necessary to find gaps in the shell.

While we do not perform a direct comparison here, the multiple bubbles observed by Beuther et al. (2022) in NGC7538, powered primarily by an O3 star, are recovered. This is caused by preferential expansion in multiple directions by the wind bubble creating chimneys following density minima in the cloud, followed by spreading out in a plume outside the minima. This is visible in a 2D slice in Fig. 15.

### 3.3.3 The role of resolution

A further point that will likely affect our results is the requirement for a minimum injection radius for stellar winds. Pittard et al. (2021) argue that in order for stellar winds to be produced, the wind must be injected within a radius of $r_{ini, max} = \sqrt{p_{wind}/(4\pi P_{amb})}$, where $p_{wind}$ is the momentum of winds from the star and $P_{amb}$ is the ambient pressure. The reason for this is that at larger injection radii, the initial phase of overpressurization in the wind bubble is not resolved, and thus the resulting bubble risks stalling or expanding too slowly. Our star has a momentum injection rate of $5.7 \times 10^{27} \text{ g cm s}^{-1}$. In our Fiducial model, the typical ambient pressure is $10^{-9} \text{ erg cm}^{-3}$. This gives $r_{ini, max} = 0.22 \text{ pc}$. By comparison, our injection radius of five cells at the highest refinement level gives an injection radius $r_{inj}$ of 0.15 pc. This is sufficient to produce a wind, however Pittard et al. (2021) argue that a ratio $r_{ini}/r_{ini, max} < 0.1$ is desirable to reproduce the correct wind bubble dynamics.

By comparison, in Seed DUHAMEL, the ambient pressure is roughly $10 \times$ higher, $r_{ini, max} = 0.069 \text{ pc}$, or two-three cells at our highest refinement level. This offers a plausible explanation for why a wind bubble struggles to form in this environment.

The conditions during the 'hot champagne' event in Seed CARNOT should not in principle violate the pressure-resolution requirement given above, with the temporary disappearance of the wind bubble being due to a large growth in turbulent mixing in the H II region, including the location of the source star. However, we caution that resolution limits in general may affect results in surprising ways, and this phenomenon may be affected by both numerical and physical conditions around the star.

Noting the already strenuous requirements on producing high-resolution simulations that include stellar winds, an outstanding problem remains in identifying methods for producing stellar wind bubbles accurately in complex, high-pressure environments.

### 3.3.4 Additional effects

It is worth pointing out that even with multiple realizations of the simulated system, our work misses some effects that likely shape regions such as Orion. For example, protostellar jets have been argued to shape similar star-forming nebulae, both in observations (Kavak et al. 2022) and simulations (Grudić et al. 2022; Guszejnov et al.

2022; Verliat et al. 2022). Similarly, we do not include cosmic rays (see review by Girichidis et al. 2020), although Morlino et al. (2021) argue that wind bubbles may not be significant sources of cosmic rays. We also do not simulate a wide range of molecular cloud densities, masses, and metallicities as found across galaxies, nor do we allow multiple massive stars to form as feedback sources. Expanding this work thusly by expanding the parameter space of the problem, though likely costly in computational resources, will likely be necessary to expand our understanding of how massive stars transfer energy to the wider universe.

## 4 CONCLUSIONS

The goal of this paper is to understand how winds shape H II regions in detail by simulating an object similar to the Orion nebula (M42) and exploring a set of physical parameters that influence such systems. The goal is to establish how robust our results are to changes in numerical and physical recipes, to identify whether there are significant gaps in our basic understanding of how H II regions evolve, and thus to inform how more complex systems should be tackled. We self-consistently form a star similar to the ionizing source of M42 ($\theta^1$ Ori C) inside a turbulent cloud. We vary the feedback processes included, physical recipes for the treatment of interfaces in the H II region, and the random seed used to generate the initial conditions.

The surprising outcome of this effect is that feedback becomes less effective when winds are included. This is because the ionizing photon budget is reduced by recombination, which occurs at a rate roughly proportional to density squared. The recombination rate in a dense, wind-swept shell, is thus higher than a typical density profile initially around a young star that decreases with radius. The dense shell is thus found to delay the rapid expansion of ionization fronts around young stars, in agreement with the analytic solutions found in Geen & de Koter (2022).

We discuss whether we reproduce observational metrics for M42, including a dense shell containing neutral (non-photoionized) hydrogen (Pabst et al. 2019), an embedded wind bubble (Guedel et al. 2008) and a similar radial expansion (a maximum extent of 4 pc with an age of roughly 200 kyr; see Pabst et al. 2020). We are to first order broadly successful with these criteria. A key aspect missing from previous work is the presence of a neutral shell rather than a rapidly escaping champagne flow. We argue that champagne flows can still occur, but mainly when the H II region encounters a sharp density discontinuity as in Tenorio-Tagle (1979) and Whitworth (1979). Champagne flows stemming from steep density gradients (Franco et al. 1990) can effectively be trapped by the shell around the wind bubble (see Geen & de Koter 2022). We plot a schematic depicting the general structure of the resulting nebula in Fig. 17.

We find that even at our small scales, the wind bubble is highly turbulent, with strong eddy mixing leading to thermal pressure losses in the hot wind-shocked gas. To correct for numerical overcooling of expanding shock-fronts in static grids (Fierlinger et al. 2016; Gentry et al. 2017), we run simulations with a mask over the contact discontinuity between the hot wind bubble and the gas outside and turn off cooling inside the mask. We also add a criterion that increases the resolution on pressure gradients, i.e. in the interface around the H II region. The cooling mask does reduce cooling, with the energy in the wind bubble increased by a factor of a few, but not to completely adiabatic levels, as argued from observational estimates by, e.g. Pabst et al. (2019). Instead, we argue the expansion rate of a wind bubble away from the star can be increased if the wind bubble is not spherical, but rather composed of a narrow chimney with a wider plume of hot wind-shocked gas at the edge of the H II region.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Projections of column density of gas cells with a speed greater than 5 km s⁻¹ in the SEEDS set. The top row is projected along the x-axis, and the bottom row along the z-axis. The goal is to track mass swept up by feedback and identify where dense shell-like structures remain. Photoionized champagne flows are evident where dense bubble-like structures have been burst, e.g. in the top-left panel, where the expanding gas on the left side of the nebula merges into the background without a clear transition indicating a dense shell.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Diagram showing the schematic behaviour of an Orion (M42)-like nebula as modelled in our simulations. The winds stream away from the source star in a collimated chimney, constrained in other angles by denser gas around the star, before expanding into a turbulent plume. Photoionized gas fills out the volume up to the plume and ionizes part of the neutral shell driven by the nebula. This shell traps the ionizing radiation, slowing the emergence of a rapid 'champagne' flow of photoionized gas.</figcaption></figure>

Finally, some of our simulations with different initial random seeds exhibit wind bubbles that shrink or disappear entirely. We track a simulation that does this, and find that a form of catastrophic mixing occurs once a champagne flow breaks out of the cloud, forming a 'hot champagne' flow, where the wind bubble as evidenced by gas above $10^6$ K disappears entirely. In one of our simulations a wind bubble struggles to form, likely because the environment it forms in is dense enough for the simulation to fall below resolution requirements set by Pittard et al. (2021). This raises a difficult prospect for resolving stellar winds with consistent accuracy in molecular cloud simulations. Despite this, the energetics of the wind bubbles does show some convergence in three of the five seeds studied.

Our work succeeds in reproducing some of the observed features of nebulae such as M42 in Orion, including an extended wind bubble surrounded by a dense shell not overrun by a champagne flow. However, the influence of winds in the H II region is highly non-linear, restricting the expansion of the photoionized region and displaying turbulent mixing. This adds to the overall cost of simulations that wish to reproduce stellar feed-

back in a molecular cloud environment. Further work remains to explore a wider parameter space relevant to star formation in galaxies.

## ACKNOWLEDGEMENTS

The authors would like to thank Xander Tielens, Lex Kaper, Ben Keller, and Eric Pellegrini for useful discussions during the writing of this paper. The authors would additionally like to thank the anonymous referee for their careful and detailed comments that helped improve the clarity of the manuscript. The simulations in this paper were performed on the Dutch National Supercomputing cluster Snellius at SURFsara. The authors gratefully acknowledge the SURFsara ResearchDrive facility for remote data sharing. We gratefully acknowledge the Common Computing Facility (CCF) of the LABEX Lyon Institute of Origins (ANR-10 LABX-0066). Copies of the data were also stored on the data storage service SDS@hd supported by the Ministry of Science, Research and the Arts Baden-Württemberg (MWK) and the German Research Foundation (DFG) through grant INST 35/1314-1 FUGG. SG acknowledges support from a NOVA grant for the theory of massive star formation and a Spinoza award of the Dutch Science Organization (NWO) for research on the physics and chemistry of the interstellar medium. TK was supported by the National Research Foundation of Korea (NRF-2020R1C1C1007079 and 2022R1A6A1A03053472).

## 5 DATA AVAILABILITY

The data products associated with this paper, including simulation code, initial conditions, analysis code, intermediate data products, and (given adequate storage) the full simulation outputs, can be made available by the authors upon reasonable request.

## REFERENCES

* Ali A., Harries T. J., Douglas T. A., 2018, MNRAS, 477, 5422

* Ali A. A., Bending T. J. R., Dobbs C. L., 2022, MNRAS, 510, 5592
Alves J., Bouy H., 2012, A&A, 547, A97

Audit E., Hennebelle P., 2005, A&A, 433, 1

Balega Y. Y., Chentsov E. L., Leushin V. V., Rzaev A. K., Weigelt G., 2014, Astrophys. Bull., 69, 46

Barnes A. T., Longmore S. N., Dale J. E., Krumholz M. R., Kruijssen J. M. D., Bigiel F., 2020, MNRAS, 498, 4906

Bate M. R., 2019, MNRAS, 484, 2341

Beuther H. et al., 2022, A&A, 659, A77

Bieri R., Naab T., Geen S., Coles J. P., Pakmor R., Walch S., 2023, MNRAS, 523, 6336

Bleuler A., Teyssier R., 2014, MNRAS, 445, 4015

Bleuler A., Teyssier R., Carassou S., Martizzi D., 2014, Comput. Astrophys. Cosmology, 2, 16

Brown A. G. A., de Geus E. J., de Zeeuw P. T., 1994, A&A, 289, 101

Chevance M., Krumholz M. R., McLeod A. F., Ostriker E. C., Rosolowsky E. W., Sternberg A., 2023, in Inutsuka S.-i., Aikawa Y., Muto T., Tomida K., Tamura M., eds, ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 1

Comeron F., 1997, A&A, 326, 1195

Dale J. E., Bonnell I. A., Clarke C. J., Bate M. R., 2005, MNRAS, 358, 291

Dale J. E., Ercolano B., Bonnell I. A., 2012, MNRAS, 424, 377

Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, MNRAS, 442, 694

Draine B. T., 2011, ApJ, 732, 100

Drass H., Haas M., Chini R., Bayo A., Hackstein M., Hoffmeister V., Godoy N., Vogt N., 2016, MNRAS, 461, 1734

Dwarkadas V. V., 2022, Galaxies, 10, 37

Ekström S. et al., 2012, A&A, 537, A146

Ferland G. J., 2003, ARA&A, 41, 517

Fielding D. B., Ostriker E. C., Bryan G. L., Jermyn A. S., 2020, ApJ, 894, L24

Fierlinger K. M., Burkert A., Ntormousi E., Fierlinger P., Schartmann M., Ballone A., Krause M. G. H., Diehl R., 2016, MNRAS, 456, 710

Franco J., Tenorio-Tagle G., Bodenheimer P., 1990, ApJ, 349, 126

Fromang S., Hennebelle P., Teyssier R., 2006, A&A, 457, 371

Gatto A. et al., 2017, MNRAS, 466, 1903

Geen S., de Koter A., 2022, MNRAS, 509, pp.4498

Geen S., Soler J. D., Hennebelle P., 2017, MNRAS, 471, 4844

Geen S., Watson S. K., Rosdahl J., Bieri R., Klessen R. S., Hennebelle P., 2018, MNRAS, 481, 2548

Geen S., Bieri R., Rosdahl J., de Koter A., 2021, MNRAS, 501, 1352

Gentry E. S., Krumholz M. R., Dekel A., Madau P., 2017, MNRAS, 465, 2471

Georgey C., Ekström S., Meynet G., Massey P., Levesque E. M., Hirschi R., Eggenberger P., Maeder A., 2012, A&A, 542, A29

Girichidis P. et al., 2020, Space Sci. Rev., 216, 68

Gritschneder M., Naab T., Walch S., Burkert A., Heitsch F., 2009, ApJ, 694, L26

Grudić M. Y., Guszejnov D., Offner S. S. R., Rosen A. L., Raju A. N., Faucher-Giguère C.-A., Hopkins P. F., 2022, MNRAS, 512, 216

Guedel M., Briggs K. R., Montmerle T., Audard M., Rebull L., Skinner S. L., 2008, Science, 319, 309

Guszejnov D., Grudić M. Y., Offner S. S., Faucher-Giguère C. A., Hopkins P. F., Rosen A. L., 2022, MNRAS, 515, 4929

Kahn F. D., 1954, Bull. Astron. Inst. Neth., 12, 187

Kahn F. D., 1980, A&A, 83, 303

Kavak U., Goicoechea J. R., Pabst C. H. M., Bally J., van der Tak F. F. S., Tielens A. G. G. M., 2022, A&A, 660, A109

Kimm T., Bieri R., Geen S., Rosdahl J., Blaizot J., Michel-Dansac L., Garel T., 2022, ApJS, 259, 21

Klessen R. S., Heitsch F., Mac Low M.-M., 2000, ApJ, 535, 887

Kraus S. et al., 2009, A&A, 497, 195

Kruijssen J. M. D. et al., 2019, Nature, 569, 519

Kuiper R., Hosokawa T., 2018, A&A, 616, 22

Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, ApJ, 914, 89

Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, ApJ, 914, 90

Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021c, ApJ, 922, L3

Lebreuilly U., Commerçon B., Laibe G., 2019, A&A, 626, A96

Leitherer C., Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014, ApJS, 212, 14

Mac Low M.-M., McCray R., 1988, ApJ, 324, 776

McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148

Mathews W. G., 1967, ApJ, 147, 965

Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315

Morlino G., Blasi P., Peretti E., Cristofari P., 2021, MNRAS, 504, 6096

O'Dell C. R., Kollatschny W., Ferland G. J., 2017, ApJ, 837, 151

Olivier G. M., Lopez L. A., Rosen A. L., Nayak O., Reiter M., Krumholz M. R., Bolatto A. D., 2021, ApJ, 908, 68

Pabst C. et al., 2019, Nature, 565, 618

Pabst C. H. M. et al., 2020, A&A, 639, A2

Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., Galván-Madrid R., Keto E. R., 2010, ApJ, 711, 1017

Pittard J. M., Wareing C. J., Kupilas M. M., 2021, MNRAS, 508, 1768

Puls J., Sundqvist J. O., Markova N., 2015. Proc. IAU Symp. Vol. 307, New Windows on Massive Stars: Asteroseismology, Interferometry, and Spectropolarimetry. Cambridge Univ. Press, Cambridge, p. 25

Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2017, MNRAS, 470, 4453

Rathjen T.-E. et al., 2021, MNRAS, 504, 1039

Rey-Raposo R., Dobbs C., Agertz O., Alig C., 2017, MNRAS, 464, 3536

Rio N. D., Robberto M., Hillenbrand L. A., Henning T., Stassun K. G., 2012, ApJ, 748, 14

Rogers H., Pittard J. M., 2013, MNRAS, 431, 1337

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Rosdahl J. et al., 2018, MNRAS, 479, 994

Rosen A. L., Lopez L. A., Krumholz M. R., Ramirez-Ruiz E., 2014, MNRAS, 442, 2701

Rosen A. L., Offner S. S. R., Foley M. M., Lopez L. A., 2021, preprint (arXiv:2107.12397)

Schneider N. et al., 2020, PASP, 132, 104301

Simón-Díaz S., Herrero A., Esteban C., Najarro F., 2006, A&A, 448, 351

Spitzer L., 1962, Physics of Fully Ionized Gases, 2nd edn. Interscience Publishers, New York.

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Tan B., Oh S. P., Gronke M., 2021, MNRAS, 502, 3179

Tenorio-Tagle G., 1979, A&A, 71, 59

Teyssier R., 2002, A&A, 385, 337

Tumlinson J., Peebles M. S., Werk J. K., 2017, ARA&A, 55, 389

Verliat A., Hennebelle P., González M., Lee Y.-N., Geen S., 2022, A&A, 663, A6

Vink J. S., Muijres L. E., Anthonisse B., de Koter A., Graefener G., Langer N., 2011, A&A, 531, 132

Walch S. K., Whitworth A. P., Bisbas T., Wünsch R., Hubber D., 2012, MNRAS, 427, 625

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

Whitworth A., 1979, MNRAS, 186, 59

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