The Energy and Dynamics of Trapped Radiative Feedback with Stellar Winds

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 Msun 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'.


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. 2022), shaping the interstellar medium of galaxies (McKee & Ostriker 1977;Rathjen et al. 2021;Bieri et al. 2023) and even affecting the composition of the circumgalactic medium (Tumlinson et al. 2017) and driving cosmic reionisation (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 photoionising 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 radiation field of the star drives mass loss from the star as stellar winds, which can reach outflow ★ E-mail: s.t.geen@uva.nlvelocities over 1000 km/s for massive stars (Puls et al. 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, magnetised, self-gravitating clouds.Simulations of photoionisation in a cloud context have been carried out by, e.g., Dale et al. (2005), Gritschneder et al. (2009), Peters et al. (2010), among many others.Photoionisation-only feedback simulations with self-consistent star formation, e.g.Dale et al. (2012), typically display rapid photoionisation of large portions of the cloud material followed by dispersal of the cloud due to thermal expansion.Ali et al. (2018) also reproduce this for a system similar to that in M42, and find some similarities between their simulations and M42 in, e.g.ionised gas velocities.
A key feature of photoionisation-driven cloud dispersal is the "champagne" flow (Tenorio-Tagle 1979), in which a density gradient in the cloud triggers a rapid breakout of ionised 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 et al. 1990).Comeron (1997) also study this effect with winds included, which produces an embedded wind plume within the photoionised 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 sizes of feedback-driven shells around massive stars ( > 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 to photoionisation.This is also found in simulations of larger spiral arm regions by Ali et al. (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 onto protostars.Supernovae generally occur too late to have a principal role in shaping star formation for molecular clouds with a global freefall time ≲ 7 Myr.

Feedback Efficiency and Radiative Cooling
The efficiency of stellar feedback depends not only on the total quantity of energy deposited into the interstellar medium by stars, but also the ability for the interstellar medium to retain this energy as kinetic flows that impact larger scales.For example, Walch et al. (2012) argues that typically less than 0.1% of the energy from ionising radiation is converted to kinetic energy in the gas via photoionisation.This is because while this process is important in thermalising 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 ionising radiation can have an additional impact via direct radiation pressure, this is typically weak at Hii region radii more than a few pc (e.g.Olivier et al. 2021).
Stellar winds rely on the overpressurisation 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% for an adiabatic wind bubble in a uniform medium (Weaver et al. 1977) or 33% in a singular isothermal sphere density profile ( ∝  −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.Firstly, 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).Kruĳssen et al. (2019) argues that such cooling is inefficient on the timescales for feedback in molecular clouds.Thermal evaporation can be enhanced by thermal conduction via electrons (Spitzer 1962), which Fierlinger et al. (2016) finds reduces the energy in a typical wind bubble by 10% when running 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 photoionised 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 generalised environments (Fielding et al. 2020;Tan et al. 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 photoionisation (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 et al. 2021).The inclusion of winds creates an additional need to resolve the short timesteps demanded by including hot and/or fast wind-driven flows, with characteristic speeds exceeding 1000 km/s.
Additionally, certain numerical features of simulation codes can cause problems in accurately resolving the energetics of hot structures in the interstellar medium.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 to 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 the OB1 cluster's massive stars and supernovae (Brown et al. 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  1 Ori C, with secondary contributions from the O9.5V star  2 Ori A (O'Dell et al. 2017).Simón-Díaz et al. (2006) analyse this cluster and find an age of 2.5 ± 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) analyses 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 starforming 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 HII 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 multiwavelength 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 Hii 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.

Paper Overview
We run a series of simulations of wind bubbles around a selfconsistently formed star, similar to the star  1 Ori C that powers the Orion Nebula M42, in a cloud similar to Orion, which has been the focus of multiple recent multi-wavelength 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 photoionised Hii region and the neutral cloud, and the role of randomness in setting the initial conditions.Finally, in we present our conclusions in Section 4.

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 magnetohydrodynamic (RMHD) Eulerian Adaptive Mesh Refinement (AMR) code ramses (Teyssier 2002;Fromang et al. 2006;Rosdahl et al. 2013).Details of the full data reproduction package are available in Section 5.
The simulations are similar to the setup described in Geen et al. (2021), in which we simulated a set of massive stars between 30 and 120 M ⊙ 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 Hii 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.

Initial Conditions
Our initial conditions use the same setup as the diffuse cloud in Geen et al. (2021), since it matches closely the gas density distribution found in Orion (Geen et al. 2017).We impose an initially spherically symmetric density field where  0 = 1.80 × 10 −21 g cm −3 (equivalent to a hydrogen number density  0 = 1078 cm −3 ),  c = 2.533 pc and  ISM = 2.24 × 10 −24 g cm −3 (equivalent to a hydrogen number density  ISM = 1.34 cm −3 ).The initial temperature is set to 10 K inside 6  c and 8000 K outside.The simulation volume is a cubic box of length 24  c = 60.8 pc.The total mass inside 6  c from the centre of the box is 10 4 M ⊙ .We set up an initial magnetic field, oriented in the  direction.The magnetic field strength along each line of sight in this direction is given by where  max,ini = 3.76 G is the peak initial magnetic field strength, Σ x is the gas column density in the  direction and Σ max,ini is the peak column density.This is a nominal initial magnetic field strength corresponding to an initial Alfvén crossing time to freefall 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  = 0, with no further turbulent forcing applied.The turbulent velocity field's root-mean sphere velocity,  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.

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 3  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 et al. 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 12  c .
Our maximum refinement level is 11, giving a minimum cell length of Δ = 0.03 pc.This can be reached under three conditions.Firstly, we refine any gas that is denser than ten times the Jeans density,  J ≡ ( s /Δ) 2 /, where  s is the sound speed in the cell.Secondly, we refine any cell more massive than 0.25 M ⊙ .Thirdly, we refine any cell whose neighbour has a 99% difference in pressure (i.e.2| left −  right |/( left +  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.

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 photoionisation, dust absorption and direct radiation pressure.The handling of sources of stellar radiation are discussed in Section 2.6.
We track photons in groups representing far ultraviolet (FUV, between 11.2 and 13.6 eV) and 3 extreme ultraviolet (EUV) groups capable of ionising Hi to Hii (between 13.6 eV and 24.59 eV), Hei to Heii (24.59 eV to 54.4 eV) and Heii to Heiii (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 ionisation front speeds.The radiation field can evolve multiple times in one (magneto)hydrodynamic timestep (termed "subcycling") if the radiation timestep is shorter than the MHD timestep, each governed by the typical crossing time of radiation and gas flows in a cell.
Each gas cell tracks the ionisation state of hydrogen and helium.Gas becomes ionised through photoionisation and collisional ionisation, 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 ( abs ) for all EUV groups, and scattering ( sc ) for all groups, is given by 10 3 / ref cm 2 /g, 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  ref = 0.02.We do not track the creation and destruction of dust as a separate fluid, as in, e.g., Lebreuilly et al. (2019).Instead, we approximate the effect of dust by including an additional absorption and scattering opacity, which are given by 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. R = ( R /) 1/4 is the local radiation temperature, where  R is the energy of the photons in each group in the cell and  ≡ 7.566 × 10 −15 erg cm −3 K −4 is the radiation constant,  HII is the local hydrogen ionisation fraction and  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.

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 timestep.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 ionised 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 ionised 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  = 0.014, i.e.Solar metallicity.We ignore metal enrichment by stellar mass loss as we expect this to be small compared to 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 preshock.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  = 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".

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% 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 timestep, this sink accretes 90% 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 ⊙ has accreted in total onto 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 a 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 ⊙ , similar to  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 timestep, 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).

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 summarise the module here.
We do not consider the pre-main sequence phase of the 35 M ⊙ 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 (1 FUV, 3 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% 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 onto 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 onto 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 Hii region around the star, which in turn affects how radiation from the star is allowed to propagate into the surrounding medium at  ∼ 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 thermalises, 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 ac-counts 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 timestep 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 nonetheless 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 Hii 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.

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.

Combining Winds and Radiative Feedback
We now compare the influence of stellar winds and UV radiation feedback in shaping Hii 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.
Figure 1 shows a sequence of images of each simulation highlighting the presence of cold, neutral gas (purple), warm, photoionised gas (orange) and hot, wind-shocked gas (light yellow), with contours outlining the wind bubble (cyan) and photoionised gas (red) separating the phases of the Hii 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 ionising radiation.This is also seen in the mean spherical-equivalent radial expansion of the region (Figure 2, where the radius is calculated as (3 II /4) 1/3 , where  II is the Hii region volume including both photoionised and hot, wind-shocked gas above 10 5 K) and momentum in outflows as a result of feedback (Figure 3, including only components of gas flows moving radially away from the star).In the UV Only run, there is a very rapid expansion of the Each image shows a projection of emission from collisional excitation from cold gas ( 2 for cells where  < 1000 K), hot gas ( 2 for cells where  > 106 K), and recombination emission ( 2  HII , ignoring cells where  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 ionisation fraction  HII > 0.1, i.e. outlining where the extent of the Hii region lies.Some gas with a lower ionisation fraction is visible outside this region, as the ionisation front is not a sharp discontinuity in runs that exhibit champagne flows.The columns from left to right show runs with winds and UV photoionisation, winds only, UV photoionisation 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 normalisation effects and does not reflect a drop in photoionisation.
Hii 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 ionising radiation any more, causing the ionisation front to move outwards supersonically without significant gas flows.This can be seen in the UV Only run in Figure 2, where the mean radius of the Hii 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 Figure 3 advances more slowly with no real sign of a slowing at 0.1 Myr as the photoionised gas begins to respond hydrodynamically (see Franco et al. 1990).
The trapping of ionising 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 Hii 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 Figure 2), the ionisation 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 Figure 4, where we plot slices through the position of the star showing gas density.In the UV Only run, the gas is rapidly photoionised as described in Franco et al. (1990), and as such there is no shell around parts of the Hii region that form a Champagne flow.Rather, the photoionised cloud expands thermally.By contrast, although the champagne phase has begun in the Wind & UV run in this figure, there is still a notice- able dense shell around the diffuse wind bubble that acts to absorb significant quantities of ionising 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 Figure 1, which use cut-off thresholds of  HII < 10 −1 and  HII < 10 −4 respectively.Despite this, nearly all of the ionising radiation is trapped, compared to the rapid champagne flow seen in the UV Only runs.

The Physics at the Hii Region Interface
A key component in Hii regions is the interface between different phases, in particular the contact discontinuity between the wind bub-ble and the gas around it, which is comprised of either a dense, partially-or a fully-photoionised shell of material.
To this end, we study two recipes designed to improve our simulations' ability to capture the physics of this region compared to 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).

Wind-only Comparison
We first compare results in a test setup without ionising radiation from the star, in order to simplify the comparison between runs with and without the wind contact discontinuity cooling mask.
In Figure 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  w above 10 5 K is a sphere, i.e. (3 w /4) 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 photoionisation.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 Figure 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 Figure 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 summarise, 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 sub-grid mixing, which is suppressed by the mask.The downside is that the mask suppresses both physical, unresolved cooling (Rosen et   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 panel shows the mean spherical-equivalent radius of the wind bubble  w (where the wind bubble volume  w = 4/3   3 w , and  w is the volume of gas in cells above 10 5 K).The right 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.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) and feedback (right) sets.Including UV photoionisation 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.
although this is unavoidable without fully resolving the wind bubble contact discontinuity, which in turn becomes computationally expensive.

Comparison including Photoionisation and Refinement
We now consider a comparison of the effect of the cooling mask including photoionisation feedback and refinement on pressure gradients, i.e. at the edge of the Hii region.In Figure 7 we plot the energy retained in the wind bubble with and without the cooling mask and refinement on pressure gradients (left) as well as with and without photoionisation (right).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 photoionisation also increases the energy retained in the wind bubble.Photoionisation 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, photoionised Hii region and surrounding cold cloud in Figure 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 Hii 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 photoionised gas fills in the volume between the star and the edge of the Hii region.In the "Mask, Refine" simulation, the photoionised 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 Hii region in more detail in Figure 9.We plot here the maximum radius, i.e. the distance of the furthest cell in the wind bubble or Hii region from the star.Measuring the radius by using the Hii region volume  HII , where  HII = (3 HII /4) 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 Hii region than for runs without the mask.Unlike the Wind Only simulations, UV radiation photoionises any substructure trapped by the wind bubble.As the Hii region expands further, the wind bubble also expands into the underdensity created by the Hii 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 Hii 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.

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 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, we show gas density, bulk gas speed, the inverse of gas vorticity (i.e. the eddy turnover timescale at a scale of 0.12 pc (4 times the smallest cell size)), and the ratio of characteristic eddy speed (ℓ ∇ × v) to bulk gas speed.By comparing the eddy speed to 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 Hii region.The circle around the star is the free-streaming radius of the wind.The white dots in the bottom right figure are where the lower grid resolution creates points where a cell is compared to itself, resulting in zero eddy speed.
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 (Kruĳssen 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 idealised 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 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 ( > 10 5 K) or photoionised gas ( HII > 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, photoionised gas back-fills the solid angles subtended by the plume at smaller radii, since the ionising radiation travels directly outwards rather than seeking pressure gradients as the wind bubble does.to resolve.Nonetheless, 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 characterised by fluid eddies at multiple scales which diffuse energy from large to small scales.There are a number of ways to characterise this in simulations.The most simple 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 (∇ × v), which gives the turnover frequency of fluid eddies at a given scale.
In Figure 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 ℓ of 4 simulation cells at the highest level of refinement (ℓ = 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 timescale (1/∇ × v) and the ratio between the characteristic eddy speed ℓ∇ × v and the bulk gas speed in each cell.Some grid artefacts can be seen due to the lower cell resolution outside the wind bubble where we do not fully refine the grid.
The eddy turnover timescale is longest (i.e.darker green in the bottom left panel of Figure 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 photoionised 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 remains a tough limitation to overcome, particularly given the already difficult computational requirements in capturing the timesteps needed to trace flows travelling at 1000s of km/s.

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 Figure 7, we find a factor of a few difference 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 % of the energy injected by the star in the gas (including both thermal and kinetic components).This is compared to 1-2 % when no mask is implied.This is significantly lower energy retention than the 33% 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% 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 Figure 11, showing the solid angle subtended by the wind bubble (in red) and the total Hii 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  of the full 4 of the solid angle.This is an average across all directions from the star including the smaller "back-blast" seen in the bottom left of each panel in Figure 10 and earlier figures showing the same simulation.
Additionally, by examining both Figures 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 photoionised 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 Hii regions appearing to have large wind bubble volume filling fractions may in fact exhibit a large plume near the edge of the Hii 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.

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 initialise the turbulent velocity field in the cloud, which in turn sets the distribution of density peaks that form in the cloud, causes feedback to behave in different ways.
In Figure 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 Hii 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 Hii regions purely by varying the initial seed.In some cases the photoionised 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 Hii region for each seed is visible in Figure 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'Enghien 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 Figure 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 ionising 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 Figure 14 we plot the energy emitted by the star as stellar winds retained in the wind bubble.An adiabatic wind should retain 45% of its energy in the hot wind bubble for a uniform cloud (Weaver et al. 1977) or 33% in a  ∝  −2 power law density field (Geen & de Koter 2022).Three of the seeds retain 10-20% 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.Nonetheless, despite the differences in the wind bubble evolution, 3 out of 5 of the simulations display evolution in their energetics after the first 0.1 Myr.

"Hot" Champagne
Seeds Ampere, Bellecour and d'Enghien 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 onto 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 Hii regions entering a breakout champagne flow phase.Figure 15 displays the progression of the Hii region at the point a champagne flow begins.The depressurisation of the photoionised region allows the wind bubble to expand more quickly and mix with the photoionsed gas, causing a slightly hotter Hii 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 6 K) wind-shocked and warm (∼ 10 4 K) photoionised 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 Figure 15.This is combined with other effects such as photoionisation in the champagne flow being temporarily disrupted.As the phenomenon is linked to out-of-dynamical-equilibrium behaviour in the Hii region following its encountering a discontinuity, this effect is also temporary.
If found in physical Hii 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 photoionised region to the spectrum of the star, assuming that residual energy from the wind bubble remains in the photoionised gas for an appreciable amount of time.A Hii 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.

Shell Structure
Comparisons with full synthetic observations are left to future work, given the complexity of matching the full velocity structure of line emission from tracers such as C + (e.g.Schneider et al. 2020).However, we approximate the observable features of the neutral shell around the Hii regions by measuring the column density of gas in cells with a total speed faster than 5 km/s, in order to isolate 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 Strong wind bubbles retain up to 10% of the energy from stellar winds, while weak bubbles retain around 1%, with some results moving between the two limits.(For reference, an adiabatic wind bubble evolving in a singular isothermal sphere retains 33% 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.typically moving at lower velocities.We plot this quantity in Figure 16.
In the top image (x-projection, used in the other figures in this paper), large-scale photoionised 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 (zprojection), 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 primarly 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 Figure 15.

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  ini,max = √︁  wind /(4 amb ), where  wind is the momentum of winds from the star and  amb is the ambient pressure.The reason for this is that at larger injection radii, the initial phase of overpressurisation 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 × 10 27 g cm/s.In our Fiducial model, the typical ambient pressure is 10 −9 erg/cm −3 .This gives  ini,max = 0.22 pc.By comparison, our injection radius of 5 cells at the highest refinement level gives an injection radius  inj of 0.15 pc.This is sufficient to produce a wind, however Pittard et al. (2021) argue that a ratio  ini / ini,max < 0.1 is desirable to reproduce the correct wind bubble dynamics.
By comparison, in Seed Duhamel, the ambient pressure is roughly 10× higher,  ini,max = 0.069 pc, or 2-3 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 Hii 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 highresolution simulations that include stellar winds, an outstanding problem remains in identifying methods for producing stellar wind bubbles accurately in complex, high-pressure environments.

Additional Effects
It is worth pointing out that even with multiple realisations 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, Verliat et al. 2022and Guszejnov et al. 2022).Similarly, we do not include cosmic rays (see review by Girichidis et al. 2020), although Morlino et al. (2021) argues 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.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.Photoionised gas fills out the volume up to the plume and ionises part of the neutral shell driven by the nebula.This shell traps the ionising radiation, slowing the emergence of a rapid "champagne" flow of photoionised gas.

CONCLUSIONS
The goal of this paper is to understand how winds shape Hii 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 Hii regions evolve, and thus to inform how more complex systems should be tackled.We self-consistently form a star similar to the ionising source of M42 ( 1 Ori C) inside a turbulent cloud.We vary the feedback processes included, physical recipes for the treatment of interfaces in the Hii 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 ionising 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 ionisation 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-photoionised) 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 Hii 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 Figure 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 Hii 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 Hii region.
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 do show some convergence in 3 of the 5 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 Hii region is highly non-linear, restricting the expansion of the photoionised region and displaying turbulent mixing.This adds to the overall cost of simulations that wish to reproduce stellar feedback in a molecular cloud environment.Further work remains to explore a wider parameter space relevant to star formation in galaxies.

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.

Figure 1 .
Figure1.Sequence of normalised emission maps showing the evolution of the nebula around the star in the feedback set, where we vary whether winds or UV photoionisation are included.Each image shows a projection of emission from collisional excitation from cold gas ( 2 for cells where  < 1000 K), hot gas ( 2 for cells where  > 106 K), and recombination emission ( 2  HII , ignoring cells where  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 ionisation fraction  HII > 0.1, i.e. outlining where the extent of the Hii region lies.Some gas with a lower ionisation fraction is visible outside this region, as the ionisation front is not a sharp discontinuity in runs that exhibit champagne flows.The columns from left to right show runs with winds and UV photoionisation, winds only, UV photoionisation 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 normalisation effects and does not reflect a drop in photoionisation.

Figure 2 .Figure 3 .
Figure 2. Mean spherical-equivalent Hii region radius  HII (including wind and photoionisation-heated gas, where the region volume  HII = 4 3  3 HII ) as a function of time in each simulation in the feedback set.Winds act to constrain runaway photoionisation of the whole cloud by trapping the radiation in the shell around the wind bubble.

Figure 4 .
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 photoionisation are included.The columns show the same simulations as Figure 1.Runs with winds contain dense shells around the wind bubble, which absorbs UV radiation from the star.

Figure 6 .
Figure6.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 photoionisation).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.

Figure 8 .Figure 9 .
Figure8.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 photoionised and colder gas is neutral cloud material.A translucent blue contour outlines the masked cells in the runs that include a cooling mask.

Figure 11 .
Figure11.Solid angle subtended by the wind bubble (red) and the overall Hii region (black, includes both photoionised 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 ( > 10 5 K) or photoionised gas ( HII > 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, photoionised gas back-fills the solid angles subtended by the plume at smaller radii, since the ionising radiation travels directly outwards rather than seeking pressure gradients as the wind bubble does.

Figure 12 .Figure 13 .Figure 14 .
Figure 12.As the Wind & UV simulation shown in Figure 1 but with the seeds set showing the effect of different random seeds in the initial turbulent field on the resulting nebula.

Figure 15 .
Figure15.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 photoionised fraction (ignoring collisional ionisation, which occurs at higher temperatures in the wind bubble).Time evolves from top to bottom, showing snapshots as photoionised gas breaks out a champagne flow, where the wind bubble mixes efficiently with the photoionised gas to the point where no hot (> 10 6 K) gas remains visible.The phenomenon is transient and lasts for less than 0.1 Myr.

Figure 16 .
Figure 16.Projections 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.Photoionised champagne flows are evident where dense bubble-like structures have been 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.

Figure 17 .
Figure17.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.Photoionised gas fills out the volume up to the plume and ionises part of the neutral shell driven by the nebula.This shell traps the ionising radiation, slowing the emergence of a rapid "champagne" flow of photoionised gas.

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 setup of the simulations.UV means that UV photoionisation is included.Direct radiation pressure is also included in the UV runs except in the UV Only runs -see