TY - JOUR
AU - de Koter, Alex
AB - ABSTRACT Winds from young massive stars contribute a large amount of energy to their host molecular clouds. This has consequences for the dynamics and observable structure of star-forming clouds. In this paper, we present radiative magnetohydrodynamic simulations of turbulent molecular clouds that form individual stars of 30, 60, and 120 solar masses emitting winds and ultraviolet radiation following realistic stellar evolution tracks. We find that winds contribute to the total radial momentum carried by the expanding nebula around the star at 10 per cent of the level of photoionization feedback, and have only a small effect on the radial expansion of the nebula. Radiation pressure is largely negligible in the systems studied here. The 3D geometry and evolution of wind bubbles is highly aspherical and chaotic, characterized by fast-moving ‘chimneys’ and thermally driven ‘plumes’. These plumes can sometimes become disconnected from the stellar source due to dense gas flows in the cloud. Our results compare favourably with the findings of relevant simulations, analytic models and observations in the literature while demonstrating the need for full 3D simulations including stellar winds. However, more targeted simulations are needed to better understand results from observational studies. methods: numerical, stars: formation, stars: massive, stars: winds, outflows, H ii regions 1 INTRODUCTION Stars form from gas that collapses under gravity. The more massive stars eject large quantities of energy and mass over their lifetimes. This process is often called ‘feedback’, because it has the capability to regulate future star formation by driving away or evaporating dense gas. These feedback processes include radiation at multiple wavelengths, jets, and other protostellar outflows, winds, and supernovae. Feedback is a ‘multiscale’ process, that is, it affects a large range of spatial scales. On the smallest scales, stars regulate their own formation. This has been modelled by Kuiper & Hosokawa (2018) for massive stars above 50 M⊙, and by Bate (2019) for clusters of less massive stars (up to 10 M⊙), amongst others. Stars also regulate gas flows on cloud scales of 1–100 pc (see review by Dale 2015, and later references in this section), collectively cause the thermal phases of the interstellar medium (ISM, e.g. McKee & Ostriker 1977; Gatto et al. 2017) and even drive flows out of galaxies (see review on galactic winds by Veilleux, Cecil & Bland-Hawthorn 2005) and shape the ionization state of the Universe (e.g. Rosdahl et al. 2018). In this paper, we focus on the interaction between two categories of feedback processes on cloud scales, namely high-energy radiation and winds, in the first Myr of the main sequence of massive stars. The first, photoionization feedback, is driven by the ionization of interstellar material by photons above the ionization energy of hydrogen (13.6 eV). This heats the gas to approximately 104 K, which creates a pressure difference between the ionized gas and neutral material outside, causing the ionized bubble to expand (Oort & Spitzer 1955). Photons can also affect the gas via direct radiation pressure when it interacts with gas and dust in the ambient medium (Mathews 1967). The second, stellar wind feedback, is the ejection of material from the surface of the star through radiation pressure exerted on the gas in the star’s atmosphere. Around young massive stars, this material can be accelerated to thousands of km s−1, according to models of, for example, Castor, Weaver & McCray (1975), Abbott & Lucy (1985), Krtička & Kubát (2004), and Muijres et al. (2012). This has been confirmed observationally by, for example, Groenewegen & Lamers (1989), Prinja, Barlow & Howarth (1990), and Crowther et al. (2016). However, it is not immediately obvious how radiation and wind feedback interact, and which process has the biggest impact on its surroundings. In this paper, we analyse 3D feedback on scales of 1–100 pc from single massive stars of up to 120 M⊙ in the first million years of their life using magnetohydrodynamic (MHD) simulations, disentangling and quantifying the effects of photoionization, radiation pressure, and stellar winds. The interplay of these modes of feedback, together with initial conditions, create complex 3D geometries of cold neutral, hot ionized, and ultra-hot x-ray emitting gas characterized by kinetic energy-driven ‘chimneys’ and thermally driven ‘plumes’ that may be disconnected from the stellar source. The wind bubble and photoionized gas combined are referred to as an H ii region, since both contain ionized hydrogen. We will use this definition throughout the paper. 1.1 Analytic models for winds and photoionization Early analytic work by Kahn (1954), Spitzer (1978), Whitworth (1979), and others confirms that photoionization feedback is capable of driving gas flows into the ISM. The same is true for analytic work focused on adiabatic stellar wind bubbles (bubbles of hot gas driven by stellar winds), by Avedisova (1972), Castor et al. (1975), and Weaver et al. (1977). More recently, authors have studied the interaction between these two processes in more detail. Analytic calculations by Capriotti & Kozminski (2001) suggest that winds are not likely to have a significant dynamical effect on the expansion of H ii regions compared to photoionization. Krumholz & Matzner (2009) argue that leakage from stellar wind bubbles due to fragmentation or inhomogeneities in the shell further reduces the dynamical input from winds. Haid et al. (2018) confirm this with analytic models and simulations using a uniform dense, neutral background, but argue that once the wind has entered the already ionized ISM, photoionization cannot drive further outflows and winds are required. This is because the warm ionized medium outside the cloud has a similar ionization state and temperature to the H ii region, and so there is no pressure difference across the H ii region radius. Geen et al. (2020) use analytic models of winds and radiation expanding into a power-law density field, as is expected for recently formed massive stars (Lee & Hennebelle 2018). In this model, winds only become more important than photoionization close to massive stars (< 0.1 pc). This is because the energy or momentum from winds is spread across a smaller surface area. Rahner et al. (2017) argue that winds from massive stellar clusters do contribute a significant fraction of the force acting on the H ii region, peaking at around 3 Myr when massive stars begin to emit stronger winds during the Wolf–Rayet (WR) phase. However, Silich & Tenorio-Tagle (2017) argue that this depends on the ability for the wind bubbles around individual stars to merge, with isolated wind bubbles being less effective. Fierlinger et al. (2016) argue that winds deposit around 2–3 times the energy from supernovae into the surrounding material, and further that careful modelling of the mixing of hot and cold gas at the bubble interface is crucial for determining how much of the energy deposited by winds is lost to radiative cooling. 1.2 Structure and observability Harper-Clark & Murray (2009), Yeh & Matzner (2012), and Yeh et al. (2013) construct quasi-static 1D models of H ii regions including photoionization, radiation pressure, and winds. They argue that nearby observed H ii regions are consistent with models in which winds are not dynamically significant. Pellegrini et al. (2007, 2012) and Pellegrini, Baldwin & Ferland (2011) argue that winds are required to explain the observed structure of these regions. This is dynamically significant since winds shape the density of the photoionized region, which in turn affects the thermal pressure of the photoionized gas. They also argue that there are certain regions, such as the Orion Veil nebula, for which pressure equilibrium has not been reached and which are not well described by such quasi-static models. Guedel et al. (2007) find extended X-ray emission inside the Orion Veil nebula, arguing that winds fill the H ii region. Pabst et al. (2019) find that the shell around the Orion Veil nebula travels faster than the sound speed in ionized gas, while the cooling from X-ray emission is low, and thus the expansion of the region is best explained by adiabatic wind models versus a photoionization-driven model. Kruijssen et al. (2019) also argue that the dispersal of molecular clouds by adiabatic stellar winds and photoionization should happen at the same rate. 1.3 Numerical simulations These (semi-)analytic analyses do not always agree, since by their nature they rely on simplified assumptions, both to reduce computational cost and to create solutions that are easy to follow from first principles. Their geometry is often simplified to enable this, neglecting the complex 3D distribution of structures in star-forming regions in galaxies or using only simple clumping factors to account for their omission. There is value in this genre of models for their ability to populate a large parameter space quickly, but they also have innate limitations that must be overcome in order to create more realistic theoretical models for feedback in star-forming regions. There is thus a need for more comprehensive, if costly, self-consistent 3D radiation-hydrodynamic simulations to study this phenomenon. When stellar winds are included, the cost of these simulations is increased dramatically. Due to the high temperature of wind bubbles (106–108 K, or even higher), satisfying the Courant condition forces the time-step of the hydrodynamic simulations to be much lower than for simulations with just photoionization (with characteristic temperatures of ∼104 K). As smaller spatial scales are resolved, this time-step drops further. The problem of combining stellar winds and photoionization has thus to this point still not been explored fully. In this paper, we focus on molecular cloud scales. On protostellar scales, protostellar jets and outflows are the main feedback sources, with winds and ionizing radiation expanding into the cavities created by these processes (Kuiper & Hosokawa 2018), while at larger scales (Agertz et al. 2013; Gatto et al. 2017) winds add to the thermal pressure in hot gas in the galaxy and reduce the Galactic star formation efficiencies (SFE). Grudić et al. (2020) include winds and radiation in simulations of very massive (>106 M⊙) clouds, considering star particles to be well-sampled populations for the purpose of feedback. Our work should thus be seen as a bridge between these scales, tracing the flow of wind-driven structures from subparsec to ∼100 pc scales. Simulations of wind outflows on cloud scales by Rogers & Pittard (2013) and Rey-Raposo et al. (2017) demonstrate that winds escape preferentially through low-density channels, reducing their effectiveness at dispersing clouds. Dale et al. (2014), who for the first time include both photoionization and wind feedback with self-consistent star formation on a molecular cloud scale, find that the dynamical role of winds is small compared to photoionization. Mackey, Langer & Gvaramadze (2013) and Mackey et al. (2015) argue, using simulations of stars moving at varying speeds with respect to the background, that winds lose most of their energy to evaporation and mixing, with photoionization being the principal driver of H ii regions around the star. None the less, emission from the interface between the wind bubble and the gas around it is an important observational tracer (Green et al. 2019). Magnetic fields have also often been omitted from simulations with radiative and wind feedback on a cloud scale due to the additional computational cost. However, as we showed in Geen et al. (2015), magnetic fields are important for the structure of H ii regions since they limit the breakup of filaments and shells (see also Hennebelle 2013). Recent work by Wall et al. (2020) simulates self-consistent photoionizaton and winds with MHD, although since the paper focuses on resolving stellar multiplicity, in their highest resolution model they do not form stars larger than 10 M⊙, which have weaker winds and ionizing photon emission rates than more massive stars. 1.4 Outline of this work In this work, we present radiative MHD simulations of photoionization, radiation pressure, and wind feedback in turbulent molecular clouds. We follow the formation of massive stars self-consistently through sink particle accretion. However, in order to isolate the effects of stellar winds in controlled conditions, we allow only one massive star to form of a pre-selected mass of either 30, 60, and 120 M⊙. The advantage of this approach is that the source of early stellar feedback is in a realistic position within the cloud, that is, on a gas density peak. From this star, we track feedback according to a full single-star evolution model (see Section 2). Sink particle accretion, representing the formation of lower mass stars, continues. We then follow the evolution of the wind bubble and photoionized region. Our goal is to study the complex geometry of the wind bubble, the feedback efficiencies and interactions of stellar winds and radiation. In Section 2, we discuss the methods used to set up and run our simulations. In Section 3, we present the results of these simulations, focusing on the evolution of the wind bubble. In Section 4, we compare our results to analytic models and observations, and discuss the results of our simulations in the context of the wider literature. Finally, in Section 5, we summarize our conclusions. 2 NUMERICAL SIMULATIONS this section, we describe the numerical setup of the simulations used in this paper (see Table 1 for a full list). Each of the simulations describes an isolated molecular cloud with an initial turbulent velocity field, magnetic field, self-gravity, and stellar feedback. All of the simulations are performed with the radiative MHD Eulerian adaptive mesh refinement (AMR) code ramses (Teyssier 2002; Fromang, Hennebelle & Teyssier 2006; Rosdahl et al. 2013). The total CPU time used in these simulations was approximately 500 kh. More specific details about the setup of the code can be found in the Data Availability statement in Section 1.1. Table 1. List of simulations included in this paper. Cloud refers to the cloud setup used (see Table 2). M* refers to the mass of the star used as a source of winds and radiation. NOFB indicates that no feedback is included. UV indicates that UV photoionization is included. WIND indicates that stellar winds are included. PRESS indicates that radiation pressure is included. See Section 2 for a discussion of how these effects are implemented. Simulation name . Cloud . M* / M⊙ . UV . WIND . PRESS . NOFB diffuse – ✗ ✗ ✗ UV_30 diffuse 30 ✓ ✗ ✗ UVWIND_30 diffuse 30 ✓ ✓ ✗ UVWINDPRESS_30 diffuse 30 ✓ ✓ ✓ UV_60 diffuse 60 ✓ ✗ ✗ UVWIND_60 diffuse 60 ✓ ✓ ✗ UVWINDPRESS_60 diffuse 60 ✓ ✓ ✓ UV_120 diffuse 120 ✓ ✗ ✗ UVWIND_120 diffuse 120 ✓ ✓ ✗ UVWINDPRESS_120 diffuse 120 ✓ ✓ ✓ NOFB_DENSE dense – ✗ ✗ ✗ UV_120_DENSE dense 120 ✓ ✗ ✗ UVWIND_120_DENSE dense 120 ✓ ✓ ✗ UVWINDPRESS_120_DENSE dense 120 ✓ ✓ ✓ Simulation name . Cloud . M* / M⊙ . UV . WIND . PRESS . NOFB diffuse – ✗ ✗ ✗ UV_30 diffuse 30 ✓ ✗ ✗ UVWIND_30 diffuse 30 ✓ ✓ ✗ UVWINDPRESS_30 diffuse 30 ✓ ✓ ✓ UV_60 diffuse 60 ✓ ✗ ✗ UVWIND_60 diffuse 60 ✓ ✓ ✗ UVWINDPRESS_60 diffuse 60 ✓ ✓ ✓ UV_120 diffuse 120 ✓ ✗ ✗ UVWIND_120 diffuse 120 ✓ ✓ ✗ UVWINDPRESS_120 diffuse 120 ✓ ✓ ✓ NOFB_DENSE dense – ✗ ✗ ✗ UV_120_DENSE dense 120 ✓ ✗ ✗ UVWIND_120_DENSE dense 120 ✓ ✓ ✗ UVWINDPRESS_120_DENSE dense 120 ✓ ✓ ✓ Open in new tab Table 1. List of simulations included in this paper. Cloud refers to the cloud setup used (see Table 2). M* refers to the mass of the star used as a source of winds and radiation. NOFB indicates that no feedback is included. UV indicates that UV photoionization is included. WIND indicates that stellar winds are included. PRESS indicates that radiation pressure is included. See Section 2 for a discussion of how these effects are implemented. Simulation name . Cloud . M* / M⊙ . UV . WIND . PRESS . NOFB diffuse – ✗ ✗ ✗ UV_30 diffuse 30 ✓ ✗ ✗ UVWIND_30 diffuse 30 ✓ ✓ ✗ UVWINDPRESS_30 diffuse 30 ✓ ✓ ✓ UV_60 diffuse 60 ✓ ✗ ✗ UVWIND_60 diffuse 60 ✓ ✓ ✗ UVWINDPRESS_60 diffuse 60 ✓ ✓ ✓ UV_120 diffuse 120 ✓ ✗ ✗ UVWIND_120 diffuse 120 ✓ ✓ ✗ UVWINDPRESS_120 diffuse 120 ✓ ✓ ✓ NOFB_DENSE dense – ✗ ✗ ✗ UV_120_DENSE dense 120 ✓ ✗ ✗ UVWIND_120_DENSE dense 120 ✓ ✓ ✗ UVWINDPRESS_120_DENSE dense 120 ✓ ✓ ✓ Simulation name . Cloud . M* / M⊙ . UV . WIND . PRESS . NOFB diffuse – ✗ ✗ ✗ UV_30 diffuse 30 ✓ ✗ ✗ UVWIND_30 diffuse 30 ✓ ✓ ✗ UVWINDPRESS_30 diffuse 30 ✓ ✓ ✓ UV_60 diffuse 60 ✓ ✗ ✗ UVWIND_60 diffuse 60 ✓ ✓ ✗ UVWINDPRESS_60 diffuse 60 ✓ ✓ ✓ UV_120 diffuse 120 ✓ ✗ ✗ UVWIND_120 diffuse 120 ✓ ✓ ✗ UVWINDPRESS_120 diffuse 120 ✓ ✓ ✓ NOFB_DENSE dense – ✗ ✗ ✗ UV_120_DENSE dense 120 ✓ ✗ ✗ UVWIND_120_DENSE dense 120 ✓ ✓ ✗ UVWINDPRESS_120_DENSE dense 120 ✓ ✓ ✓ Open in new tab 2.1 Initial conditions and refinement criteria We use two sets of initial conditions, one of a diffuse cloud similar to the nearby Gould Belt clouds, and one denser (Geen, Soler & Hennebelle 2017). In both of them, we define a cloud with an initially spherically symmetric density profile n(r) defined by $$\begin{eqnarray*} n(r) = n_0 / (1 + (r/r_c)^2) \end{eqnarray*}$$(1) where n0 and rc are the central hydrogen number density and characteristic radius, respectively. This profile is imposed out to a radius rini = 3rc, where n(rini) = 0.1n0. Outside this, a uniform sphere is imposed up to 2 rini, with a density 0.1 times that just inside rini, or 0.01 n0, to provide a reservoir of material to accrete on to the cloud. The initial temperature inside 2 rini is set to 10 K. Outside this radius, the hydrogen number density is set to 1 cm−3 and the temperature to 8000 K. The total length of the cubic volume simulated Lbox = 16 rini. Note that the cloud evolves significantly between the start of the simulation and the time the first star forms. There are two clouds used in this study, both with an initial mass of 104 M⊙. One is a cloud similar to the nearby Gould belt as established in Geen et al. (2017) by comparing the column density distributions from our simulated clouds and the observed clouds. The other is denser, to test the effects of feedback in different environments. These are, respectively, the ‘L’ and ‘S’ clouds in Geen et al. (2017). We list the properties of both of these clouds in Table 2. Table 2. List of cloud setups included in this paper, as described in Section 2. Mc is the initial cloud mass in M⊙. tff is the initial free-fall time of the cloud as a whole. tsound is the sound crossing time. tA is the Alfvén crossing time. trms is the crossing time of the root mean square velocity of the initial turbulence of the cloud. Lbox is the box length. Δx is the minimum cell size. rc is the characteristic radius of the central isothermal part of the cloud at t = 0 (see equation 1). n0 is the central hydrogen number density. Bmax, ini is the maximum magnetic field strength in the initial seed magnetic field. Cloud name . log(Mc / M⊙) . tff / Myr . tff/tsound . tff/tA . tff/trms . Lbox / pc . Δxmin / pc . rc / pc . n0 / cm−3 . Bmax, ini / |$\mu \mathrm{G}$| . diffuse 4 4.22 0.15 0.2 2.0 122 0.03 2.533 823.4 3.76 dense 4 0.527 0.075 0.2 2.0 30.4 0.03 0.6335 52700 60.1 Cloud name . log(Mc / M⊙) . tff / Myr . tff/tsound . tff/tA . tff/trms . Lbox / pc . Δxmin / pc . rc / pc . n0 / cm−3 . Bmax, ini / |$\mu \mathrm{G}$| . diffuse 4 4.22 0.15 0.2 2.0 122 0.03 2.533 823.4 3.76 dense 4 0.527 0.075 0.2 2.0 30.4 0.03 0.6335 52700 60.1 Open in new tab Table 2. List of cloud setups included in this paper, as described in Section 2. Mc is the initial cloud mass in M⊙. tff is the initial free-fall time of the cloud as a whole. tsound is the sound crossing time. tA is the Alfvén crossing time. trms is the crossing time of the root mean square velocity of the initial turbulence of the cloud. Lbox is the box length. Δx is the minimum cell size. rc is the characteristic radius of the central isothermal part of the cloud at t = 0 (see equation 1). n0 is the central hydrogen number density. Bmax, ini is the maximum magnetic field strength in the initial seed magnetic field. Cloud name . log(Mc / M⊙) . tff / Myr . tff/tsound . tff/tA . tff/trms . Lbox / pc . Δxmin / pc . rc / pc . n0 / cm−3 . Bmax, ini / |$\mu \mathrm{G}$| . diffuse 4 4.22 0.15 0.2 2.0 122 0.03 2.533 823.4 3.76 dense 4 0.527 0.075 0.2 2.0 30.4 0.03 0.6335 52700 60.1 Cloud name . log(Mc / M⊙) . tff / Myr . tff/tsound . tff/tA . tff/trms . Lbox / pc . Δxmin / pc . rc / pc . n0 / cm−3 . Bmax, ini / |$\mu \mathrm{G}$| . diffuse 4 4.22 0.15 0.2 2.0 122 0.03 2.533 823.4 3.76 dense 4 0.527 0.075 0.2 2.0 30.4 0.03 0.6335 52700 60.1 Open in new tab In the initial conditions, we impose a supersonic turbulent velocity field over the cloud. We do not apply further turbulent forcing to the cloud. Each cloud has a global free-fall time |$t_{\mathrm{ ff}} \equiv \sqrt{3 \pi / 32 G \rho _{av}}$|⁠, defined by the average density of the isothermal sphere ρav inside rini. The radius of the cloud is set via the sound crossing time tsound for a fiducial neutral gas sound speed of 0.28 km s−1. The balance of turbulence and gravity is set via the turbulent root mean square (rms) velocity Vrms, which has a crossing time trms. The magnetic field strength is set via the Alfvén wave crossing time tA. The magnetic field is initially oriented along the x-direction. We calculate a maximum initial magnetic field strength Bmax, ini using the value of the Alfvén wave crossing time |$t_{A} \equiv r_c \sqrt{\rho _0} / B_{\mathrm{ max},\mathrm{ ini}}$|⁠, where ρ0 = n0mH/X for a hydrogen mass of mH and hydrogen mass fraction X. We assign Bmax, ini to the density peak at the centre of the cloud. We calculate the gas column density Σx along each line of sight in the x-direction, and calculate the magnetic field strength Bx of each cell along a given line of sight as $$\begin{eqnarray*} B_x = B_{\mathrm{ max},\mathrm{ ini}}(\Sigma _x / \Sigma _{\mathrm{ max},\mathrm{ ini}}) \end{eqnarray*}$$(2) where Σmax, ini is the initial maximum column density in the x-direction, which intersects the density peak of the cloud. At t > 0, the magnetic field then evolves with time according to the Harten-Lax-van Leer with contact and Alfven mode (HLLD) scheme described in Fromang et al. (2006). The maximum magnetic field strength grows considerably larger than Bmax, ini over time. We ‘relax’ the clouds by running the simulations without self-gravity for 0.5 tff, in order to mix the turbulent velocity and density fields, since the density field is initially spherically symmetric (see Klessen, Heitsch & Mac Low 2000; Lee & Hennebelle 2016, amongst others). After 0.5tff, we apply self-gravity to the cloud, which allows the gas to collapse to form sink particles as described in Section 2.3. We trace the gas dynamics on an octree mesh that refines adaptively when certain conditions are met. Every time a cell at level l fulfills certain criteria, it subdivides itself into 8 child cells at level l + 1. The cell size is given by Δx = Lbox/2l. We minimally refine everywhere up to level 7, giving a cube with 27 = 128 cells on a side. Everywhere inside a sphere of diameter 8rini we fully refine up to level 9, that is, two further levels. Finally, any gas cell that is 10 times denser than the Jeans density1 or has a mass above 0.25 M⊙ anywhere in the simulation volume is refined, down to a minimum cell size Δxmin of 0.03 pc. This corresponds to level 12, or level 10 in the ‘dense’ clouds. 2.2 Cooling and radiative transfer We track the propagation of radiation across the full AMR grid using the M1 method (Rosdahl et al. 2013). The radiation is coupled to the gas via photoionization, dust absorption, and direct pressure from the transfer of photon momentum to the gas. In runs labelled ‘UV’ (see Table 1), we track extreme ultraviolet (EUV) photons above the ionization energy of hydrogen, with photons binned into three groups bounded by the ionization energies of H i, He i, and He ii. In these runs, we do not automatically include radiation pressure. Runs with radiation pressure included are labelled ‘PRESS’. In these runs, we include an additional far-ultraviolet (FUV) group between 5.7 and 13.6 eV, which interacts only via radiation pressure on dust (see below). In each grid cell, the code stores the photon density and flux for each group, and couples the photons to the gas at every time-step via the cooling function. Radiation travels at a reduced speed of light of 0.01 c, in order to reduce the cost of the radiation transport steps. This value is chosen to be similar to the maximum speed of stellar winds, and to ensure that the code can capture the speed of ionization fronts in the simulation. We subcycle the radiation step, so that the hydrodynamic time-step is not limited by the (reduced) speed of light. More details on how this is done are given in Rosdahl et al. (2013). In Bisbas et al. (2015), we use this value for the reduced speed of light in our code to compute a known test problem and find excellent agreement with other radiative transfer codes that use different techniques, with differences on the order of 1 per cent, or one grid cell. Each grid cell tracks the ionization state of hydrogen and helium. Ionization fractions change via photoionization, recombination, and collisional ionization as calculated in the radiative transfer module described in Rosdahl et al. (2013). Each group uses a ‘grey’ approximation, that is, all photons in the group are considered to have the same energy, energy- and number-weighted cross-sections, using representative values from a starburst99 (Leitherer et al. 2014) stellar population as in Geen et al. (2017). Radiation pressure is calculated in runs labelled ‘PRESS’ according to Rosdahl & Teyssier (2015). Direct radiation pressure is applied for each photon absorption event in the gas. The local gas opacity to the radiation in all ionizing photon groups is given by κabs = 103 Z/Zref cm2 g−1, where Z is the metallicity of the gas and Zref = 0.02. We use the reduced flux approximation described in (Rosdahl & Teyssier 2015, appendix B) to ensure that the correct radiation momentum pressure is applied. Each absorption event transfers momentum from the photons to the gas. In addition to the EUV bins interacting with atoms, all radiation bins interact with dust. We do not track the formation, destruction, and advection of dust self-consistently, and instead assume a fixed coupling to the gas proportional to the gas metallicity and neutral hydrogen fraction. As the EUV photons are absorbed by dust grains, this reduces the budget of photons that can ionize hydrogen and helium. Cooling rates are applied on a per-cell basis. The temperature of the photoionized hydrogen and helium evolves at each time-step following the recombination cooling and photoionization heating functions described in Rosdahl et al. (2013). To account for metals, we add a cooling rate Λmetal = Λmetal, ionx + Λmetal, neutral(1 − x) where x is the hydrogen ionization fraction. Λmetal, neutral follows the cooling module of Audit & Hennebelle (2005) that uses fits to various coolants in an ISM environment with a heating term from a radiation background below the ionization energy of hydrogen. Cooling for photoionized metals, Λmetal, ion, is described by a fit to Ferland (2003). Above 104 K, we use a fit to Sutherland & Dopita (1993) where the cooling rate is higher than our simple photoionized model. The gas is set to solar metallicity, which we take to be Z = 0.014 as in Ekström et al. (2012). We consider metal enrichment from a single star’s winds to be negligible at this metallicity. 2.3 Sinks and star formation If a gas cell is above 10 per cent of the Jeans density at the highest refinement level, it is assigned to a ‘clump’, that is, a patch of dense gas. Clump peaks are identified using the ‘watershed’ method, in which contours from high to low density are drawn, with clumps merged by identifying saddle points in the density field (Bleuler et al. 2014). If a clump is denser than the Jeans density at the highest refinement level, a sink particle is formed, and every time-step, 90 per cent of the mass in the clump above the Jeans density is accreted on to the sink particle (Bleuler & Teyssier 2014). This 90 per cent is a heuristic quantity to prevent zero or negative densities in accreting cells. Once the total mass of all sink particles in the simulation exceeds 120 M⊙, we create a ‘stellar object’ representing a massive star of mass M* and attach it to the most massive sink, as in Geen et al. (2018). The rest of the mass is considered to be stars below 8 M⊙, which produces negligible high-energy radiation or winds. None of the mass in the sink is removed, and this stellar object is simply a book-keeping tool to track the age and number of massive stars in the simulation. This age is used by the stellar evolution and feedback model to calculate the radiative and mechanical feedback from the star. This radiation, mass, momentum, and energy is deposited at the position of the host sink every time-step. We use this method, rather than assigning a fixed position to the stellar object, because as Geen et al. (2020) show, the density profile around the star has a dramatic influence on the dynamics of the radiation and wind feedback. All sinks continue to accrete to a larger mass if they continue to fulfill the criteria described above. In this paper, we run simulations where M* is either 30, 60, or 120 M⊙. We also run simulations with no stellar object (i.e. M* = 0). We pick these masses to sample a range of masses of stars where winds are expected to have a greater effect (Geen et al. 2020). Based on the estimates of Weidner, Kroupa & Bonnell (2009), these stellar masses are possible given the reservoir of gas in our simulation, although a 120 M⊙ star is significantly more likely to be formed in a more massive cloud. We do not include more stellar masses in this study due to limits in the masses included in the stellar tracks we use (see Section 2.4), although these stars have been reported (e.g. Crowther et al. 2010; Bestenlehner et al. 2011). In each simulation we only allow one massive star to form, in order to study in detail the response of the gas to winds and radiation from a single source. 2.4 Stellar evolution and feedback We implement feedback from the massive star as the emission of UV radiation and winds. The star is considered to start on the main sequence from the moment of formation, since we do not have the resolution to properly track the protostellar phase, which is typically ∼105 yr (e.g. Hosokawa & Omukai 2009). The radiation and winds are emitted from the position of the sink that the star is attached to. We follow the evolution of massive star models at solar metallicity (Z = 0.014) computed using the Geneva model (Ekström et al. 2012), assuming the stars are rotating at 0.4 of the critical velocity. For completeness, we show the resulting photon emission rates and wind properties in Appendix A. Though in the current paper we focus on the early stage of stellar evolution up to 1 Myr, we include a description of our stellar evolution tracks including older and less massive stars than the ones featured in this paper. We stop after 1 Myr due to computational cost, although by this time in all of the simulations, ionizing radiation has begun to escape the simulation volume. At each time-step we deposit radiation and winds on to the grid. The number of photons emitted per unit time in each radiation group is calculated using individual stellar spectra extracted from starburst99 (Leitherer et al. 2014), using the Geneva model as inputs for the stellar atmosphere calculations. To calculate the number of photons emitted between time t and t + Δt, we interpolate linearly between pre-computed tables for each photon group at intervals of 5 M⊙. We inject winds every time-step in the same fashion. Stellar mass-loss rates |$\dot{m}_w$| and escape velocities vesc are taken from Ekström et al. (2012). Note that mass-loss rates are uncertain by a factor of 2–3 (e.g. Muijres et al. 2012; Smith 2014; Puls, Sundqvist & Markova 2015; Keszthelyi, Puls & Wade 2017), and these models should be considered in this light. We convert the escape velocity vesc at the stellar surface to terminal wind velocities vw using the corrections given in Gatto et al. (2017). We list these here for clarity. We first calculate an ‘effective’ escape velocity veff, $$\begin{eqnarray*} v_{\mathrm{ eff}}^2 = (1 - \Gamma _e) v_{\mathrm{ esc}}^2, \end{eqnarray*}$$(3) where Γe is the Eddington factor, described in Lamers et al. (1993) as a correction factor to the Newtonian gravity set by radiation pressure on free electrons from the star. Γe is given by $$\begin{eqnarray*} \Gamma _e = \frac{\sigma _e \sigma _{\mathrm{ SB}} T_e^4}{g c}, \end{eqnarray*}$$(4) where σSB is the Stefan–Boltzmann constant, g is the surface gravity, c is the speed of light, and σe is the cross-section for electron–photon scattering per unit mass given by $$\begin{eqnarray*} \sigma _{e}= 0.4 (1 + I_{\mathrm{ He}}Y_{\mathrm{ He}}) / (1 + 4 Y_{\mathrm{ He}})~\mathrm{cm^2 \, g^{-1}} \end{eqnarray*}$$(5) where IHe is the number of free electrons per He atom or ion, and YHe is the Helium abundance by number (approx 0.1). IHe is zero below an effective surface temperature Te = 104 K, 2 above Te = 2.5 × 104 K, and 1 otherwise. All of these parameters are calculated for each time-step of each of the different Geneva stellar tracks used in this study. Γe is typically between 0.4 and 0.95 (Vink et al. 2011). We divide massive stars into different classifications as in Crowther (2007) and Georgy et al. (2012). Stars with Te > 104 K and a surface hydrogen mass fraction of less than 0.3 are WR stars. Stars below TRSG = 5000 K are red supergiants (RSG). Stars above TBSG = 8700 K (up to 104 K) are blue supergiants (BSG) and stars between 5000 and 8700 K are yellow supergiants (YSG). Stars that do not fall into these categories are OB stars. We subdivide WR stars into categories WNL and WNE, and WC and WO, depending on the surface abundances of H, He, C, N, and O. For these stars, vw is given using a clamped linear interpolation $$\begin{eqnarray*} v_w &=& v_0 \mathrm{~if~}~T_e \lt T_0 \nonumber\\ v_w &=& v_1 \mathrm{~if~}~T_e \gt T_1 \nonumber\\ v_w &=& v_0 + (v_1 - v_0) \times (T_e - T_0) / (T_1 - T_0)~\mathrm{~otherwise} \end{eqnarray*}$$(6) where v0 and v1 are reference wind velocities in km s−1. T0 and T1 are surface temperatures in units of 104 K. For OB stars we use (v0, v1, T0, T1) = (1.3 veff, 2.45 veff, 1.8, 2.3). For WO and WC stars, we use (v0, v1, T0, T1) = (700, 2800, 2.0, 8.0), and for WNL and WNE stars (v0, v1, T0, T1) = (700, 2100, 2.0, 5.0). For RSG stars, we use vw = 10 km s−1 × (L/Lref)1/4, where Lref ≡ 3 × 104 L⊙. For YSG stars, we use $$\begin{eqnarray*} v_w = 10 \mathrm{~km\,s^{-1}~} \times 10^{[(\mathrm{log}(T_e) - \mathrm{log}(T_{\mathrm{ RSG}})) / (\mathrm{log}(T_{\mathrm{ BSG}}) - \mathrm{log}(T_{\mathrm{ RSG}}))]}, \end{eqnarray*}$$(7) in order to fit the argument of Gatto et al. (2017) that the geometric mean velocity in this range is 50 km s−1, while RSG winds are typically somewhere above 10 km s−1 and BSG winds are 100 km s−1. Mass-loss rates from RSG and YSG stars are comparable to other massive stars, but the vw is lower. Thus the momentum and energy deposition rates from stellar winds are typically much weaker for these stars than OB or WR stars. We further note that this field is subject to ongoing study, and as such this model does not represent the final word in winds from massive stars. We force the cells around the sink particle with the massive star to be at the highest refinement level. We inject winds into a five cell radius around the sink. Mass, momentum, and energy are injected evenly in all cells inside this radius. The injected momentum and energy are calculated as |$\dot{m}_w v_w~\mathrm{d}t$| and |$\frac{1}{2} \dot{m}_w v_w^2~\mathrm{d}t$| respectively, where dt is the time-step. When injected into a cell with a low density, the injected mass and momentum dominate, and the solution becomes free-streaming. If the cell has a high density, the injected mass and momentum have less of an impact on the final velocity of the cell, and the injected energy effectively becomes thermalized. This means that the free streaming phase of the wind described in Weaver et al. (1977) appears as the wind bubble grows and the density of gas around the star is largely from the wind itself and not the pre-existing circumstellar medium. 3 RESULTS In this section, we present the results of the simulations and explore the effect that winds have on the photoionization region. We begin with an overview of the simulations. We then discuss the influence that winds have on bulk properties of the system. Finally, we look in detail at the evolution of the wind bubble itself, and on the role that radiative cooling plays in the wind bubble’s evolution. 3.1 Overview In Fig. 1, we plot column density maps for each simulation containing photoionization, radiation pressure, and stellar winds at 0.1 Myr intervals after the formation of the star. We overplot contours showing the position of the wind bubble and the ionization front. The wind bubble, outlined in cyan, includes all cells above 2 × 104 K or travelling faster than 100 km s−1. These thresholds are selected to exclude photoionized gas and include cooler parts of the wind bubble, although a larger temperature threshold is possible with no change in the results. The edge of the ionization front is shown in red. Figure 1. Open in new tabDownload slide Maps of column density projected in the y-axis in each simulation including photoionization, stellar winds, and radiation pressure. From top to bottom are snapshots at 0.1, 0.2, 0.3, and 0.4 Myr after the star is formed. From left to right, we plot results for the 30, 60, and 120 M⊙ star, and the 120 M⊙ star in the dense cloud. The red contours show the extent of the ionized gas around the star. The cyan contours show the extent of the wind bubble. See Section 3.1 for a description of how these are defined. The ionization front expands rapidly outwards in a champagne flow. The wind bubble expands inside this region, and has a highly aspherical geometry. Figure 1. Open in new tabDownload slide Maps of column density projected in the y-axis in each simulation including photoionization, stellar winds, and radiation pressure. From top to bottom are snapshots at 0.1, 0.2, 0.3, and 0.4 Myr after the star is formed. From left to right, we plot results for the 30, 60, and 120 M⊙ star, and the 120 M⊙ star in the dense cloud. The red contours show the extent of the ionized gas around the star. The cyan contours show the extent of the wind bubble. See Section 3.1 for a description of how these are defined. The ionization front expands rapidly outwards in a champagne flow. The wind bubble expands inside this region, and has a highly aspherical geometry. The photoionized H ii region expands faster than the wind bubble in all cases, rapidly escaping the cloud in a ‘champagne flow’, or a directional escape of photoionized gas from the cloud described by Tenorio-Tagle (1979). Franco, Tenorio-Tagle & Bodenheimer (1990) argue that supersonic ionization front expansion occurs for clouds with a power-law density profile with index w > 3/2. The wind bubble expands against the photoionized gas, which is thermalized to around 104 K. Gas pushed by the wind into the surrounding medium is already photoionized. This gas mixes with the rest of the photoionized gas rather than forming a dense neutral shell. This picture is consistent with the 1D hydrostatic spherically symmetric model we present in Geen et al. (2020). However, there are divergences between this simple analytic picture and the full 3D simulations in this work, which we now discuss. The wind bubble size increases with stellar mass and time. The shape of the bubble is highly aspherical, and follows channels of low density in the gas. At early stages, the wind bubble is confined, but as it expands it develops a fast-flowing chimney structure that reaches beyond the cloud. Once it does so, it develops Rayleigh–Taylor-like plumes of hot gas that extend into the ISM. The chimney expands in size as the wind bubble expands and the neutral cloud is dispersed. This chimney-and-plume structure is especially evident in the dense cloud. In addition, part of the plume in the dense cloud at 0.2 Myr is cut-off from the rest of the wind bubble, and cools rapidly. This is caused by dense flows in the cloud temporarily cutting off part of the wind bubble from the star. This is somewhat similar to the flickering effect seen by Peters et al. (2010) in their simulations of ultracompact H ii regions. This effect also occurs in the UVWIND_120 simulation. The shape of the wind bubble is not stable, but changes chaotically and appears to be influenced by turbulent motions in the rest of the cloud. The wind bubble’s characteristic speed (either its sound speed or wind terminal velocity, which are similar in magnitude) is around 1000 km s−1 or more, so it can rapidly respond to changes in the density and pressure in the surrounding gas. The growth of the total volume of the wind bubble is more stable, however. We do not run simulations with varying magnetic fields due to the already considerable cost of these simulations. At the point the star forms, the magnetic energy in the cloud is similar to the thermal energy, although smaller than the kinetic energy in turbulence. As the H ii region expands, the magnetic energy grows, although it lags behind the thermal and kinetic energy from feedback. We will perform a more detailed parameter study of the role of magnetic fields in future work. 3.2 Radial evolution of feedback structures The 3D shape of the wind bubble is highly aspherical and evolves chaotically. In order to trace the bulk evolution of the wind bubble using a simple, stable diagnostic, we plot the ‘sphericized’ radius of the H ii region and wind bubble in Fig. 2 for each simulation. The sphericized radius of the ionization front ri, s is defined via the total volume Vi of the H ii region or wind bubble |$V_i = \frac{4}{3} \pi r_{i,s}^3$|⁠. A similar radius rw, s is calculated for the wind bubble with volume Vw. The photoionized and wind bubble cells counted inside Vi and Vw, respectively, are calculated as in Section 3.1. Figure 2. Open in new tabDownload slide ‘Sphericized’ radius of the H ii region and wind bubble as a function of time in each simulation. Sphericized radius ri, s is defined via the total volume Vi of the H ii region |$V_i = \frac{4}{3} \pi r_{i,s}^3$|⁠. A similar radius rw, s is calculated for the wind bubble with volume Vw. Vi is the total volume of all cells above a hydrogen ionization fraction of 0.1. Vw is the total volume of all cells with a bulk velocity above 100 km s−1 or a temperature above 2 × 104 K. Dashed–dotted lines show the H ii region radius in the UV photoionization-only simulations, dotted lines show the result of simulations including UV photoionization and stellar winds, and solid lines show the results of simulations including UV photoionization, radiation pressure, and stellar winds. Thick outlined lines show the radius of the ionization front, and thin lines show the wind bubble radius in simulations including winds. On the left are the results from the diffuse cloud, and on the right are results from the dense cloud. The wind bubble radius is always smaller than the ionization front radius. Figure 2. Open in new tabDownload slide ‘Sphericized’ radius of the H ii region and wind bubble as a function of time in each simulation. Sphericized radius ri, s is defined via the total volume Vi of the H ii region |$V_i = \frac{4}{3} \pi r_{i,s}^3$|⁠. A similar radius rw, s is calculated for the wind bubble with volume Vw. Vi is the total volume of all cells above a hydrogen ionization fraction of 0.1. Vw is the total volume of all cells with a bulk velocity above 100 km s−1 or a temperature above 2 × 104 K. Dashed–dotted lines show the H ii region radius in the UV photoionization-only simulations, dotted lines show the result of simulations including UV photoionization and stellar winds, and solid lines show the results of simulations including UV photoionization, radiation pressure, and stellar winds. Thick outlined lines show the radius of the ionization front, and thin lines show the wind bubble radius in simulations including winds. On the left are the results from the diffuse cloud, and on the right are results from the dense cloud. The wind bubble radius is always smaller than the ionization front radius. The H ii regions grow slowly at first in simulations with the lower mass stars. However, in all simulations, the regions eventually accelerate outwards from the star into lower density neutral cloud material, before reaching a plateau phase. According to the expansion model of Franco et al. (1990), this phase is unstable and the ionization front can expand supersonically. The expansion eventually plateaus as the whole volume becomes ionized. The radius of the H ii region in simulations with winds is only slightly larger than the radius without. The wind bubble itself is only a small fraction of the H ii region’s radius in all cases, although it is a larger fraction in the dense cloud. We discuss the ratio between these two radii later in the paper. The flickering of the wind bubble in the dense cloud is seen at around 300 kyr. Radiation pressure has almost no effect on the results. The main difference is in the volume of the wind bubble in the dense cloud. This is because in the UVWIND_120_DENSE simulation, the whole wind bubble plume is cut-off from the star and cools, whereas in the UVWINDPRESS_120_DENSE simulation, only part of the plume is cut-off. We attribute this to the chaotic nature of the wind bubble geometry rather than any systematic effect. 3.3 Radial momentum added to the cloud The momentum in the cloud as a whole in flows directed away from the star in each simulation is shown in Fig. 3. This includes outflows from the star, flows driven outwards in the H ii region by feedback and the component of turbulent flows in the neutral gas in the outward radial direction. For each star, simulations containing different stellar feedback physics are shown in different line styles to demonstrate the contribution of each effect. Figure 3. Open in new tabDownload slide Radial outwards momentum from the position of the star as a function of time in each simulation. Dashed lines show the results of the NOFB simulations, dashed–dotted lines show the UV photoionization-only simulations, dotted lines show simulations with UV photoionization and stellar winds, and solid lines show the result of simulations including UV photoionization, radiation pressure, and stellar winds. On the left are the results from the diffuse cloud, and on the right are results from the dense cloud. The momentum in the NOFB run evolves due to the interplay between turbulence (some of which is oriented in a radial direction) and gravity. The contribution from winds is at most 10 per cent in any simulation set, while the contribution from radiation pressure is mostly negligible. The results only show momentum in directions away from the star, and ignores flows moving towards the star. Figure 3. Open in new tabDownload slide Radial outwards momentum from the position of the star as a function of time in each simulation. Dashed lines show the results of the NOFB simulations, dashed–dotted lines show the UV photoionization-only simulations, dotted lines show simulations with UV photoionization and stellar winds, and solid lines show the result of simulations including UV photoionization, radiation pressure, and stellar winds. On the left are the results from the diffuse cloud, and on the right are results from the dense cloud. The momentum in the NOFB run evolves due to the interplay between turbulence (some of which is oriented in a radial direction) and gravity. The contribution from winds is at most 10 per cent in any simulation set, while the contribution from radiation pressure is mostly negligible. The results only show momentum in directions away from the star, and ignores flows moving towards the star. The boost in momentum from winds or radiation pressure in the 30 M⊙ case is negligible. The boost from winds becomes larger as the stellar mass is increased, but is no more than 10 per cent in the 120 M⊙ star’s case. Winds never add more momentum to the flows around the star than photoionization. Radiation pressure provides a negligible contribution to the outflowing momentum, adding an imperceptible amount of momentum on top of the momentum added by photoionization and winds. Our dust model is simplified, and a more complex model may give a different result. However, we do not expect a different model to change our results significantly given that radiation pressure plays a negligible role in our current simulations. Similar studies by Kim, Kim & Ostriker (2018) and Fukushima et al. (2020) find that radiation pressure only affects SFEs in much denser regions than the ones studied here. The biggest effect would likely be to modify the budget of ionizing photons for the photoionization process as dust grains in the H ii region absorb ionizing photons that would otherwise ionize hydrogen or helium (Krumholz & Matzner (2009) estimate a reduction of up to 27 per cent for a H ii region with typical Milky Way dust fractions). The momentum includes all flows in the cloud in the direction away from the star, and thus contains a component of the turbulent flow, visible in the simulation without feedback (NOFB). It ignores flows moving towards or perpendicular to the star. These flows can cause the expansion of feedback structures around stars to stall if the feedback source is weak or the cloud is dense, as we found in similar simulations in Geen et al. (2015, 2016). This turbulent neutral gas flow momentum also accounts for the ‘floor’ in momentum in the dense cloud and for the 30 M⊙ star. 3.4 The evolution of the embedded wind bubble 3.4.1 Structure To illustrate the typical internal structure of the wind bubble, in Fig. 4 we plot slices through simulation UVWINDPRESS_120 at 0.4 Myr after the star is formed. The slice is taken through the position of the star, from which the winds and radiation are emitted. Figure 4. Open in new tabDownload slide Slices through the simulation volume in the UVWINDPRESS_120 simulations at 0.4 Myr after the star is formed, showing various properties of the gas. The images are shown in the y-axis of the Cartesian volume, with the slice taken through the y-position of the star. Each image is 61 pc on-a-side, that is, half of the total box length. From top left to bottom right, we plot mass density, temperature, thermal pressure, total radiative cooling rate (as luminosity per unit volume), hydrogen ionization fraction, and ratio of kinetic energy Ekin to thermal energy Etherm. The radiative cooling rate is shown to illustrate where thermal energy is lost. The positions of sink particles are shown as white dots, and the star as a red star-shaped icon with a white outline. The low pressure region to the right is the shadow behind the remaining neutral gas in the cloud. In the bottom right plot, we overplot contours around the wind bubble in cyan, the photoionized gas in red and gas moving faster than 1000 km s−1 in magenta. See Section 3.1 for a description of how these are defined. At 0.4 Myr in this simulation, only a small section of the neutral cloud remains. Figure 4. Open in new tabDownload slide Slices through the simulation volume in the UVWINDPRESS_120 simulations at 0.4 Myr after the star is formed, showing various properties of the gas. The images are shown in the y-axis of the Cartesian volume, with the slice taken through the y-position of the star. Each image is 61 pc on-a-side, that is, half of the total box length. From top left to bottom right, we plot mass density, temperature, thermal pressure, total radiative cooling rate (as luminosity per unit volume), hydrogen ionization fraction, and ratio of kinetic energy Ekin to thermal energy Etherm. The radiative cooling rate is shown to illustrate where thermal energy is lost. The positions of sink particles are shown as white dots, and the star as a red star-shaped icon with a white outline. The low pressure region to the right is the shadow behind the remaining neutral gas in the cloud. In the bottom right plot, we overplot contours around the wind bubble in cyan, the photoionized gas in red and gas moving faster than 1000 km s−1 in magenta. See Section 3.1 for a description of how these are defined. At 0.4 Myr in this simulation, only a small section of the neutral cloud remains. The structure of the bubble agrees in qualitative terms with the schematic presented in the theoretical work of Weaver et al. (1977). In this paper, the authors describe different structures inside the wind bubble. Around the star, a free-streaming wind leaves the star as material travels outwards at the terminal velocity of the wind. This region has a lower thermal pressure than its surroundings, as seen in Fig. 4, since the flow is mostly kinetic. At some radius, the wind shocks against the ambient medium, creating a hot diffuse bubble. Since the hot wind bubble is collisionally ionized, it is almost completely transparent to ionizing radiation. A warm (104 K) photoionized nebula is found outside the wind bubble. Despite this qualitative agreement, the 3D geometry of the wind is very different from the purely spherical models of Weaver et al. (1977). This is due to the complex 3D geometry and motions in the cloud. The wind bubble moves preferentially in certain directions of lower density out of the cloud. ‘Chimneys’ of gas moving at above 1000 km s−1 are visible in the slice pointing along the direction the wind bubble escapes. These are shown as purple contours in the bottom right image in Fig. 4. These chimneys are not free-streaming, but are also not hydrostatic, and show evidence of bulk outwards flow. Closer to the edge of the bubble, a mostly thermalized region of shocked gas can be seen. This part of the wind bubble is in rough thermal pressure equilibrium with the surroundings. The surrounding gas is made of ionized cloud material that has not yet expanded into the surrounding medium after being overtaken by the photoionized champagne flow. 3.4.2 Energetics In this section, we discuss how the ability of the wind bubble to retain energy from the star determines its dynamics. We focus on two aspects of this issue. First, we calculate how much energy from the stellar wind is retained in the gas, and where this energy is lost. Secondly, we look at how energy from the winds influences the dynamics of the system. This allows us to disentangle the role of winds from that of other processes in the cloud. Throughout this section, we consider two extremes for the expansion of the wind bubble, from strongest to weakest influence on the cloud, as discussed in Silich & Tenorio-Tagle (2013). In the strongest case, the wind bubble is adiabatic and thus retains all of the energy deposited by the star as winds. It expands due to its internal overpressure as described in Avedisova (1972) and Weaver et al. (1977). In the weakest case, the bubble cools very efficiently and thus loses all of the deposited wind energy very quickly. In this mode, the bubble is accelerated only at a rate equivalent to the momentum output rate of the star. With moderate cooling rates, the bubble evolves between one of these two extremes. It is crucial to understand which best describes a given wind bubble, since both models produce drastically different outcomes for the expansion of the wind bubble, with the adiabatic mode expanding considerably faster than the efficiently cooled mode. In Fig. 5, we plot the energy lost to radiative cooling as a function of time in the hot wind bubble (i.e. ignoring radiative losses in the warm or cold gas) and compare it to the wind energy input from the star as a function of time. We identify cells as being in the wind bubble using the same method as in Section 3.1. The cooling rate inside the wind bubble is at most ∼10  per cent of the emitted wind luminosity. We further plot the cooling from cells above 106  K, which will predominantly emit in X-rays, and find an even lower fraction of the wind bubble’s energy is lost in this gas phase. This can also be extracted from the panel displaying Lcool in Fig. 4, where it is clearly visible that the cooling rate inside the bubble is indeed low. One might thus expect the wind bubble to be close to adiabatic. However, when we measure the total energy stored in the wind bubble, it is typically only 1  per cent of the total wind energy injected. This is explained by a region of enhanced cooling along the interface between the wind bubble and the denser photoionized gas outside. In an adiabatic wind bubble solution, some of the energy in the wind bubble is stored in the dense shell around the wind bubble. In our simulations, we do not find a dense shell, and it appears that this energy is instead lost in mixing with the rapidly cooling photoionized gas. Figure 5. Open in new tabDownload slide Wind luminosity of the star at each simulation output (solid line) versus rate of energy loss to radiative cooling in the wind bubble (dashed line) and the amount of that cooling in gas above 106 K (dotted line) as a function of time in each simulation containing photoionization, winds, and radiation pressure. The wind bubble is identified as described in Section 3.1. Emitted luminosity values are sampled from the stellar evolution tables (see Section 2) for the star’s age in each simulation output. Cooling in gas cells above 106 K is given as an upper bound on the amount of X-rays that can be emitted by the wind bubble. The gap from 0 to 0.1 Myr is due to limited frequency of simulation outputs. There is a further delay in cooling in UVWIND_30 due to the time taken to heat the cells around the star to above 2 × 104 K. The wind bubble itself does not radiate away the majority of the energy emitted by the star as winds. Figure 5. Open in new tabDownload slide Wind luminosity of the star at each simulation output (solid line) versus rate of energy loss to radiative cooling in the wind bubble (dashed line) and the amount of that cooling in gas above 106 K (dotted line) as a function of time in each simulation containing photoionization, winds, and radiation pressure. The wind bubble is identified as described in Section 3.1. Emitted luminosity values are sampled from the stellar evolution tables (see Section 2) for the star’s age in each simulation output. Cooling in gas cells above 106 K is given as an upper bound on the amount of X-rays that can be emitted by the wind bubble. The gap from 0 to 0.1 Myr is due to limited frequency of simulation outputs. There is a further delay in cooling in UVWIND_30 due to the time taken to heat the cells around the star to above 2 × 104 K. The wind bubble itself does not radiate away the majority of the energy emitted by the star as winds. The efficient cooling of wind bubbles through mixing with the interface between the hot wind bubble and the cooler, denser gas outside the bubble is modelled analytically in Mac Low & McCray (1988). They give a cooling time tcool over which the wind bubble loses the majority of its energy, $$\begin{eqnarray*} t_{\mathrm{ cool}} = 2.3 \times 10^4~n_0^{-0.71} L_{38}^{0.29}~\mathrm{yr} \end{eqnarray*}$$(8) where n0 is the density of the ambient medium around the wind bubble and L38 is the wind luminosity in units of 1038 erg s−1. We measure the average density just outside the wind bubble to find n0, and obtain L38 from our stellar evolution tables. We thus obtain a value of tcool for our simulations, which we plot in Fig. 6. tcool is typically between 0.1 per cent and 1 per cent of the star’s age. This short cooling time is consistent with our finding that the wind bubble typically retains only a small fraction of the energy injected by stellar winds (⁠|$\sim 1~{{\ \rm per\ cent}}$|⁠), and that the wind energy of the star is mainly lost in the interface between the wind bubble and the photoionized gas outside it. Figure 6. Open in new tabDownload slide Estimated cooling time tcool of the wind bubble as a function of time in each simulation containing winds. The tcool is calculated using equation (8). The external density n0 is taken to be the median density in cells immediately neighbouring the wind bubble, and the wind luminosity is taken from our stellar evolution tables. See Section 3.4.2 for a discussion of this calculation and some caveats. Lines begin where a wind bubble is first identified in an output, which depends on the output times of the simulation. tcool is typically between 0.1 and 1 per cent of the age of the star, and thus we expect the wind bubble to lose around 90–99 per cent of its energy. Figure 6. Open in new tabDownload slide Estimated cooling time tcool of the wind bubble as a function of time in each simulation containing winds. The tcool is calculated using equation (8). The external density n0 is taken to be the median density in cells immediately neighbouring the wind bubble, and the wind luminosity is taken from our stellar evolution tables. See Section 3.4.2 for a discussion of this calculation and some caveats. Lines begin where a wind bubble is first identified in an output, which depends on the output times of the simulation. tcool is typically between 0.1 and 1 per cent of the age of the star, and thus we expect the wind bubble to lose around 90–99 per cent of its energy. It should be noted that despite the difference in cloud densities between the dense and diffuse clouds, the wind bubbles cool at a similar rate in both clouds. This can be explained by noting the geometry of the wind bubble in Fig. 5 – by 0.1 Myr, the wind bubble chimney has punched through the dense cloud and created a plume in the more diffuse gas around the cloud. This means that the wind bubble is less affected by mixing with the denser cloud material. We therefore expect our wind bubble to behave most similarly to the efficiently cooled model in Silich & Tenorio-Tagle (2013) and now compare directly to their model to determine whether this is the case. There are two factors contributing to the pressure driving the wind bubble. One is the direct momentum injection rate from the star, |$\dot{p}_w$|⁠, against the surface of the bubble. The second is the stored energy, Eb, acting to overpressure the bubble. The total pressure inside the bubble from both factors can be written as $$\begin{eqnarray*} P_b = \dot{p}_w / A + E_b / V , \end{eqnarray*}$$(9) where A is the surface area of the wind bubble and V is its volume. From classical mechanics, the force on the bubble’s surface Fb is equal to its rate of change of momentum |$\dot{p}_b = P_b A$|⁠. We thus arrive at $$\begin{eqnarray*} F_b = \dot{p}_b = \dot{p}_w + E_b (A / V) . \end{eqnarray*}$$(10) If Fb is similar to |$\dot{p}_w$|⁠, the wind bubble contains negligible stored energy, and it is driven purely by direct momentum deposition. If Fb is much larger than |$\dot{p}_w$|⁠, the wind bubble retains a large quantity of energy input by the star, that is, is partially or completely adiabatic, and its expansion is driven by overpressure inside the wind bubble. In Fig. 7, we plot this comparison as a function of time. We calculate the difference in momentum between the UVWIND and UV simulations to give the contribution to the momentum from winds, assuming this is a linear perturbation to the total momentum. We then calculate the rate of change in this difference, and compare it to the momentum output rate |$\dot{p}_w$| in winds from the star. From the time, the wind bubble is established up to a stellar age of 0.3 Myr, winds from the 60 and 120 M⊙ stars in the diffuse cloud add an order of magnitude more momentum to the wind bubble than |$\dot{p}_w$|⁠, suggesting a stored quantity of energy. However, this energy is quickly lost and |$\dot{p}_b$| drops to a similar value to |$\dot{p}_w$|⁠. Winds in the dense cloud only ever expand on the level of |$\dot{p}_w$|⁠. The 30 M⊙ star does expand faster than if it were driven purely by |$\dot{p}_w$|⁠, but does not expand smoothly, suggesting some influence of pressure from external forces such as turbulence. Figure 7. Open in new tabDownload slide The figure demonstrating the influence of stellar winds on the momentum of the H ii region, and the role of stored energy in the bubble. The output wind momentum deposition rate of the star |$\dot{p}_{w} \equiv v_w \dot{m}_w$| is plotted as a dotted line. To directly compare to this quantity, we calculate the difference in momentum Δp between the UVWIND and UV simulations in Fig. 3. We then calculate the differential of this with time, dΔp(UVWIND − UV)/dt. A solid line means that the UVWIND simulation has a larger increase in momentum per unit time than the UV simulation for the same star (i.e. adding winds increases momentum). A dashed line means that the UV simulation has a larger increase in momentum per unit time (i.e. a negative version of the solid line). If the solid line lies above the dotted line, it implies that the wind bubble is expanding faster than the efficiently cooled wind bubble model described in Section 3.4.2, and is thus partially driven by pressure induced by stored energy in the wind bubble. We find that in most cases, our wind bubble follows the efficiently cooled model. Figure 7. Open in new tabDownload slide The figure demonstrating the influence of stellar winds on the momentum of the H ii region, and the role of stored energy in the bubble. The output wind momentum deposition rate of the star |$\dot{p}_{w} \equiv v_w \dot{m}_w$| is plotted as a dotted line. To directly compare to this quantity, we calculate the difference in momentum Δp between the UVWIND and UV simulations in Fig. 3. We then calculate the differential of this with time, dΔp(UVWIND − UV)/dt. A solid line means that the UVWIND simulation has a larger increase in momentum per unit time than the UV simulation for the same star (i.e. adding winds increases momentum). A dashed line means that the UV simulation has a larger increase in momentum per unit time (i.e. a negative version of the solid line). If the solid line lies above the dotted line, it implies that the wind bubble is expanding faster than the efficiently cooled wind bubble model described in Section 3.4.2, and is thus partially driven by pressure induced by stored energy in the wind bubble. We find that in most cases, our wind bubble follows the efficiently cooled model. In our simulations, the wind bubble evolves according to an efficiently cooled model, driven only by the momentum from the current stellar wind output, with negligible stored energy. The high temperature inside the wind bubble is counterbalanced by its low density, and the stored energy is only around 1 per cent of the input wind energy. 4 DISCUSSION We have so far established a picture of a wind bubble that broadly follows the classical picture in Weaver et al. (1977), where the wind bubble expands spherically with a free-streaming volume embedded within a hot, shocked volume. However, in our simulations, the wind bubble cools rapidly and has a highly aspherical geometry that responds chaotically to structures in the photoionized cloud in which it is embedded. In this section, we compare this picture from simulations with a sample of previous analytic models and observations, and discuss how our results in Section 3 should be viewed in this context. 4.1 Comparison with analytic models Geen et al. (2020) provide algebraic expressions to describe the evolution of an efficiently cooled wind bubble embedded inside a photoionized H ii region using a hydrostatic 1D model to describe its internal structure. They characterize the behaviour of the wind bubble with a coefficient Cw. Cw > 1 if winds contribute more to the expansion of the H ii region than photoionization, and Cw < 1 if photoionization contributes more. Using equation (26) of Geen et al. (2020), we can write $$\begin{eqnarray*} C_w &=& 0.0119 \left(\frac{\dot{p}_w}{10^{29}~\mathrm{g\,cm\,s^{-2}}} \right)^{3/2} \left(\frac{Q_H}{10^{49}~\mathrm{s}^{-1}} \right)^{-3/4}\nonumber\\ &&\times \left(\frac{r_i}{1~\mathrm{pc}} \right)^{-3/4} \left(\frac{c_i}{10~\mathrm{km\,s^{-1}}} \right)^{-3}, \end{eqnarray*}$$(11) where |$\dot{p}_w$| is the momentum deposition rate from the wind, QH is the ionizing photon emission rate, ri is the radius of the ionization front, and ci is the isothermal sound speed in the photoionized gas. We set ci to 11.1 km s−1, which is the sound speed in the ionized gas in our simulations. In Geen et al. (2020), we find that Cw > 1 for stars above 60 M⊙below 0.01 pc, and above 120 M⊙ below 0.1 pc. Using values for our stellar sources at 1 pc, we find Cw = 0.005 for the 30 M⊙ star, 0.065 for the 60 M⊙ star, and 0.28 for the 120 M⊙ star. At larger radii, this value drops. This is broadly consistent with the simple hydrostatic model in Geen et al. (2020), where the contribution from winds from a 30 M⊙star is negligible, while as stellar mass increases, the contribution from winds becomes more apparent but is never the primary source of momentum from outflows on cloud scales. For more massive clusters, this may change, since the equation scales as |$\dot{p}_w^2 / Q_H$|⁠. As we add more stars to the source of the H ii region, winds will play a larger role. Geen et al. (2020) give an equation for the ratio of the ionization front radius ri and the wind bubble radius rw $$\begin{eqnarray*} \left(\frac{r_i}{r_w} \right)^{4} = 2^{1/3} C_w^{-4/3} + \frac{r_i}{r_w} . \end{eqnarray*}$$(12) In Fig. 8, we compare this expression to the relative radii of the wind bubble and ionized nebula over time in each simulation with winds. We again take the sphericized radii ri, s and rw, s shown in Fig. 2 as representative radii, where rw, s/ri, s ≡ (Vw/Vi)1/3. Vw is the volume of the wind bubble and Vi is the volume of ionized gas. Figure 8. Open in new tabDownload slide Ratio between radius of wind bubble and radius of ionized region. The thick outlined solid lines show results for our simulations, that is, (Vw/Vi)1/3, where Vw is the volume of the wind bubble and Vi is the volume of ionized gas. The thin solid lines show the results of equation (12). Figure 8. Open in new tabDownload slide Ratio between radius of wind bubble and radius of ionized region. The thick outlined solid lines show results for our simulations, that is, (Vw/Vi)1/3, where Vw is the volume of the wind bubble and Vi is the volume of ionized gas. The thin solid lines show the results of equation (12). In general, the analytic model in equation (12) somewhat overestimates of the radius of the wind bubble at early times. At late times, the analytic equation matches the simulations more closely. Note that once the H ii region reaches the edge of the box, the simulation is unable to track its whole evolution. One difference between the analytic model and our simulations is the geometry of the wind bubble. An elongated wind bubble has a higher surface area to volume ratio than a sphere. This means that the highly aspherical wind bubbles in our simulations will, for a given internal pressure, typically occupy a smaller volume than they would if they were purely spherical. The smaller, trapped wind bubble around the 30 M⊙ star provides a better match to the analytic model. A second difference is that the spherical 1D analytic model misses 3D dense gas clumps and filaments in the cloud that follow the turbulent gas flows in the cloud. The effect of this is seen in the dense cloud, where the wind radius drops temporarily as the chimney structure is cut-off and the hot plume cools. However, in this simulation the analytic model is otherwise a reasonable match to the simulation. Our results agree qualitatively with other analytic models of the interaction between winds and photoionization. Capriotti & Kozminski (2001) create a set of analytic models taking into account the cooling rate of the wind bubble, while Haid et al. (2018) produce a model that assumes zero cooling to compare to their simulations. In both cases, photoionization remains the main driver of the H ii region except in cases where the medium is already ionized, such as the diffuse ISM. Dale et al. (2014) produce a model that assumes the wind bubble is efficiently cooled, as in the Geen et al. (2020) model. However, they use a Spitzer (1978) model for the expansion of the H ii region, which requires a uniform background gas density. Dale et al. (2014) find that their model underpredicts rw/ri, which they argue is due to leakage of the H ii region into the gas outside the cloud. The Geen et al. (2020) model uses a power-law density field ρ∝r−2, which captures some of this behaviour and produces a closer fit between the simulation results and analytic model. Note that the Dale et al. (2014) simulations include multiple wind and radiation sources, which adds extra complications to the comparison. The picture thus far is one in which winds are a secondary effect in the expansion of H ii regions and the destruction of molecular clouds. They produce complex, chaotic structures, but these follow the structures shaped by other processes rather than setting the conditions in the cloud themselves. 4.2 Comparison with observations In this section, we confront our simulations with observations of the Orion nebula by Pabst et al. (2019). A subset of our simulations are to some extent representative of the observed Orion nebula. The diffuse cloud in this paper was chosen to match the global density distribution of the nearby Gould belt clouds, including Orion (see Geen et al. 2017). These cloud conditions are typical for the Milky Way. Density, metallicity, and ISM pressure can vary for extragalactic environments (e.g. Gurvich et al. 2020). The Orion Veil nebula is driven largely by a single young massive star, θ1 Ori C, which has a mass of 33 M⊙ (Balega et al. 2014), with error bars of 5 M⊙ that cover the 30 M⊙ star we simulate. The conclusions of Pabst et al. (2019) are that the Orion Veil nebula is a wind-driven bubble that cools inefficiently with no sign of influence from photoionization feedback. This is at odds with our conclusions, and the conclusions of other theoretical works. We thus confront our simulations with the Pabst et al. (2019) observations of Orion to try to determine possible reasons for this discrepancy. 4.2.1 Bulk wind bubble properties We first compare the radius, velocity, and age of the wind bubble to the values in Pabst et al. (2019). We plot these values in Fig. 9. We use the maximum radius of the wind bubble, since the Orion Veil nebula expands in only one direction away from the star, constrained by dense gas in the opposite direction. See Pellegrini et al. (2007) for a discussion of the constrained part of the nebula. The maximum radius is a lot less stable than the spherically averaged radius shown in Fig. 2, since the geometry of the wind bubble is chaotic. The velocity evolution is particularly chaotic, since wind bubbles can grow and collapse rapidly in certain directions. Figure 9. Open in new tabDownload slide Evolution of the wind bubble’s maximum radius from the star, rw, max and the expansion velocity of the wind bubble vw, max in UVWINDPRESS30, the simulation most closely matching θ1 Ori C (Balega et al. 2014) and its host cloud (Geen et al. 2017). vw, max is defined as the rate of change in rw, max. Left: vw, max versus rw, max. Right: rw, max versus the age of the star. The observationally derived results of Pabst et al. (2019) are overlaid as a point. The maximum radius of the wind bubble varies non-linearly significantly with time, due to the complex geometry and behaviour of the wind bubble. Our results overlap the velocity and radius measurement of Pabst et al. (2019). Figure 9. Open in new tabDownload slide Evolution of the wind bubble’s maximum radius from the star, rw, max and the expansion velocity of the wind bubble vw, max in UVWINDPRESS30, the simulation most closely matching θ1 Ori C (Balega et al. 2014) and its host cloud (Geen et al. 2017). vw, max is defined as the rate of change in rw, max. Left: vw, max versus rw, max. Right: rw, max versus the age of the star. The observationally derived results of Pabst et al. (2019) are overlaid as a point. The maximum radius of the wind bubble varies non-linearly significantly with time, due to the complex geometry and behaviour of the wind bubble. Our results overlap the velocity and radius measurement of Pabst et al. (2019). The results of Pabst et al. (2019) sit on top of our results in velocity–radius space. Our bubble reaches this radius at a later time than the predictions of Pabst et al. (2019), approx 0.3–0.4 Myr versus 0.2 Myr. However, since the age in Pabst et al. (2019) is estimated using a simple Weaver et al. (1977) analytic calculation, and our initial conditions may vary appreciably from those in Orion, it is not unreasonable to expect different ages for the wind bubble. They also argue for a negligible role from radiation pressure, as do we. Pabst et al. (2019) estimate a shell mass of 2600 M⊙, with a lower bound at 900 M⊙ and upper bound at 3400 M⊙ depending on assumptions in how it is calculated. For comparison, our dense cloud contains roughly 1000 M⊙in the 2 pc radius around the star, which is close to the lower bound of the observed estimate. Assuming the observed mass estimate is correct, this suggests that the volume immediately around θ1 Ori C is relatively dense. In Fig. 10, we show slices through the region. The geometry of the wind bubble is roughly similar to the Orion Veil nebula, that is, a hemisphere bounded on one side by denser gas. The gas on the other side of the nebula is already photoionized, unlike the Orion Veil nebula. However, we do not expect a perfect match between the simulations and observations in terms of geometry, since the system is turbulent and chaotic. We do not run several realizations of the cloud with different initial seeds to test the effect of chaos on the system as we do in Geen et al. (2018) for reasons of computational cost. Figure 10. Open in new tabDownload slide Slices through the position of the star in the simulation containing a 30 M⊙ star with UV photoionization, radiation pressure, and winds showing, from left to right, the gas density, temperature, and ionization state. There is no sign of a dense shell around the wind bubble, and ionizing radiation escapes from the cloud. This is in contrast to the Orion Veil nebula, which has a dense shell containing neutral hydrogen. Figure 10. Open in new tabDownload slide Slices through the position of the star in the simulation containing a 30 M⊙ star with UV photoionization, radiation pressure, and winds showing, from left to right, the gas density, temperature, and ionization state. There is no sign of a dense shell around the wind bubble, and ionizing radiation escapes from the cloud. This is in contrast to the Orion Veil nebula, which has a dense shell containing neutral hydrogen. Pabst et al. (2019) find that the X-ray emission rate of the bubble is approximately 4 × 1031 erg s−1. This is roughly an order of magnitude lower than the cooling rate of hot gas in Fig. 5, although they claim that many X-rays in Orion are absorbed by intervening gas in the denser lines of sight. Guedel et al. (2007) also argue that Orion appears to be mostly filled with wind-shocked gas using X-ray observations. They argue that regions of the Orion Veil nebula that do not appear to emit X-rays are surrounded by denser gas that absorbs these photons. This is plausible, although as we show, the geometry of wind bubbles can be highly aspherical and it is possible for the wind bubble to reach large radii while not filling volumes closer to the star. In Section 3.4.2, we find the energy injected by stellar winds is mostly lost to radiative cooling, and so the wind bubble follows a model in which its expansion is purely driven by the momentum deposition from the star. By integrating the momentum over the star’s lifetime, we can determine whether the wind bubble follows this path or retains a significant amount of energy. The total momentum deposited by the 30 M⊙ star in winds after 0.2 Myr is around 3 × 1040 g cm s−1. Using the observed shell mass of Ms = 2600 M⊙ at 0.2 Myr travelling at 13 km s−1, we get a shell momentum of 7 × 1042 g cm s−1. This is much higher than the direct momentum injection from the star. Using the lower shell mass estimate Ms = 900 M⊙also gives a much higher momentum. Assuming that the Orion Veil nebula is driven principally by winds as Pabst et al. (2019) suggest, and using equation (10), this implies that the Orion Veil nebula is driven by a store of thermal energy in the wind bubble, rather than being completely cooled. 4.2.2 Trapping of ionizing photons The observed Orion Veil nebula appears to have a large, neutral shell around the wind bubble (van der Werf, Goss & O’Dell 2013). By comparison, the wind bubbles in our simulations sit inside a larger photoionized H ii region, with little or no sign of a denser shell around the wind bubble (see Fig. 10). We discuss briefly why we do not find this phenomenon in our simulation results. We take the observational results for the shell, Ms = 2600 M⊙ at rs = 2 pc, moving at vs = 13 km s−1 into a medium with hydrogen number density n0 = 1400 cm−3. Using equation (67) in Weaver et al. (1977) and making the simplifying assumption that this is a spherical shell, we can estimate the hydrogen number density of the shell $$\begin{eqnarray*} n_s = n_0 (v_s / c_0)^2 \end{eqnarray*}$$(13) where c0 is the sound speed in the neutral gas outside, and which we approximate to be 1 km s−1. Using the values given above, we find ns ≃ 2.4 × 105 cm−3. The volume of the shell Vs is given by |$4 \pi r_s^2 \Delta r$|⁠, where Δr is the shell thickness. Approximating Ms/Vs = nsmH/X, we find Δr ≃ 0.007 pc, which is 4–5 times smaller than our finest grid cell size. The total recombination rate of such a shell is 1053 photons s−1, or 104 times the photon emission rate of the star. In other words, the shell should easily absorb all ionizing photons coming from the star. If we assume that the shell around the wind bubble can only reach a thickness of 0.1 pc, or ∼3 cells at our highest level of refinement, the total photon absorption rate is still ∼1052 photons s−1, which is higher than the photon emission rate from the star. In principle, with sufficient adaptive refinement, our simulations should be able to capture the trapping of ionizing photons by the shell. As we noted earlier, a shell mass of Ms = 2600 M⊙implies that the background density around θ1 Ori C is significantly higher than we find in our simulations. This may explain why such a dense shell is produced, and hence why the ionizing photons are trapped. It also raises the question of why the wind bubble does not cool more efficiently, however, as modelled in equation (8). More targeted simulations that make direct comparisons to specific observations are needed to understand why the case of Orion appears to differ from theoretical predictions to date about how photoionization and winds should interact. 4.3 Further considerations Full 3D hydrodynamic simulations are still relatively expensive due to the cost of simulating fast, hot flows such as stellar winds. This presents a problem for exploring a larger parameter space. By reducing the problem to 1D and making certain simplifications, some of these limitations can be overcome. Recent 1D analytic models by, for example, Rahner et al. (2017) and Pellegrini et al. (2020) are able to match certain observed properties of nearby H ii regions around massive clusters. There is still a need for 3D simulations to capture the full behaviour of molecular clouds with embedded stellar wind bubbles. We have shown that the behaviour of wind bubbles around single stellar sources is already complex. Various other authors have already begun to explore the interactions between wind bubbles, which adds a further layer of complication. Rogers & Pittard (2013) and Dale et al. (2014) perform simulations of multiple massive stellar wind sources, finding that the ablation of dense clumps as multiple wind bubbles merge around them provides an additional cooling channel. Clumping and shell fragmentation also occurs in the interaction between free-streaming winds on small scales, such as around close binaries (Calderón et al. 2020a) and the Galactic centre (Calderón et al. 2020b). A further channel for clumping is the time variability of wind velocities in the WR phase (see review by Wade 2012). Cooling rates from wind bubbles are also a matter of debate. Our simulations match the analytic model of Mac Low & McCray (1988), in which material evaporated from the wind bubble shell mixes with the hot gas and causes efficient cooling. However, a more detailed understanding of the microphysics of the shell is needed to determine whether this happens in all cases. As we resolve smaller scales, thermal conduction and thermal instabilities become important (Koyama & Inutsuka 2004). Authors such as Gentry et al. (2017) argue for lower cooling rates in hot superwind bubbles formed by multiple supernovae. Cosmic rays can also retain some dynamically important energy in hot interstellar bubbles, as they can interact with gas at larger radii (e.g. Wadepuhl & Springel 2011; Dashyan & Dubois 2020). The microphysics of gas cooling is thus important in understanding the dynamics of hot bubbles such as those driven by stellar winds. The stellar evolution framework described in this work can be extended to cover longer time-scales and larger spatial scales. In this regime, the late stage behaviour of massive stars becomes more important. As wind bubbles expand and merge into superbubbles, winds are boosted at late times relative to ionizing radiation emission, as found in Rahner et al. (2017). However, this depends on the stellar evolution model used. Sana et al. (2012) observed that the majority of massive stars are in interacting binaries, strongly affecting their evolution. One of the results of this is that the extended envelopes of these stars can be removed by binary interactions, allowing for higher ionizing emission rates at late times as in Götberg et al. (2019). 5 CONCLUSIONS We simulate a set of molecular clouds containing a single massive star of either 30, 60, or 120 M⊙ formed self-consistently using sink particles in clouds of two different densities. We track the expansion of H ii regions due to over pressure caused by photoionization and radiation pressure from photons and stellar winds produced by the star. We find that winds contribute at most 10 per cent to the outflow momentum in the first Myr of the lifetime of the star, and have only a small impact on the radial expansion of the H ii region. The contribution from winds in our simulations shows limited or no evidence for large quantities of stored energy driving expansion, as expected in models of adiabatic wind bubbles. Radiation pressure has a negligible effect on the evolution of the systems modelled in this paper. While the volume and momentum of the simulated wind bubbles evolves smoothly, the geometry of wind bubbles is highly aspherical and chaotic. The high characteristic velocity of wind bubbles means that they can rapidly evolve to fill pressure gradients in clouds and H ii regions. Outside of the classical free-streaming radius, the structure of these wind bubbles is characterized by kinetic-energy-driven ‘chimneys’ and thermally driven ‘plumes’. These plumes can be cut-off by changes in the denser gas flows in the cloud and H ii region, in some cases leading to hot bubbles not connected to a stellar source. Our simulations provide good agreement to previous simulations and analytic models in key aspects while demonstrating the need for full 3D simulations to capture the complex behaviour of stellar winds. Comparison to the Orion Veil nebula match certain bulk properties, but differences between the two systems suggest that new simulations designed to match the specific environment of observed regions are needed to close the gap between observations and theory. ACKNOWLEDGEMENTS The authors would like to thank Jo Puls, Xander Tielens, Cornelia Pabst, Patrick Hennebelle, Eric Pellegrini, Daniel Rahner, Ralf Klessen, Lex Kaper, and Zsolt Keszthelyi for useful discussions during the writing of this paper. The authors would like to further thank the anonymous referee for their careful and insightful comments that helped improve the clarity of the manuscript. The simulations in this paper were performed on the Dutch National Supercomputing cluster Cartesius at SURFsara. The authors gratefully acknowledge 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. Postprocessing was performed on the ISMSIM computing resources hosted by ITA, ZAH, University of Heidelberg. The project also made use of the SURFsara ResearchDrive facility for remote data sharing. SG acknowledges support from a NOVA grant for the theory of massive star formation. DATA AVAILABILITY This paper has been prepared according to the Research Data Management plan of the Anton Pannekoek Institute for Astronomy at the University of Amsterdam. Details of the data products and scripts used to generate the figures in this paper are found via DOI reference 10.5281/zenodo.3696806. Footnotes 1 Jeans density ρJ ≡ (cs/Δx)2/G, where cs is the sound speed in the cell and Δx is the length of the cell. Note that this does not include support from magnetic pressure. REFERENCES Abbott D. , Lucy L., 1985 , ApJ , 288 , 679 10.1086/162834 Crossref Search ADS Agertz O. , Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013 , ApJ , 770 , 25 10.1088/0004-637X/770/1/25 Crossref Search ADS Audit E. , Hennebelle P., 2005 , A&A , 433 , 1 10.1051/0004-6361:20041474 Crossref Search ADS Avedisova V. S. , 1972 , SvA , 15 , 708 Balega Y. Y. , Chentsov E., Leushin V., Rzaev A. K., Weigelt G., 2014 , Astrophys. Bull. , 69 , 46 10.1134/S1990341314010052 Crossref Search ADS Bate M. R. , 2019 , MNRAS , 484 , 2341 10.1093/mnras/stz103 Crossref Search ADS Bestenlehner J. et al. , 2011 , A&A , 530 , L14 10.1051/0004-6361/201117043 Crossref Search ADS Bisbas T. G. et al. , 2015 , MNRAS , 453 , 1324 10.1093/mnras/stv1659 Crossref Search ADS Bleuler A. , Teyssier R., 2014 , MNRAS , 445 , 4015 10.1093/mnras/stu2005 Crossref Search ADS Bleuler A. , Teyssier R., Carassou S., Martizzi D., 2014 , Comput. Astrophys. Cosmol. , 2 , 16 10.1186/s40668-015-0009-7 Crossref Search ADS Calderón D. , Cuadra J., Schartmann M., Burkert A., Prieto J., Russell C. M. P., 2020a , MNRAS , 493 , 447 10.1093/mnras/staa090 Crossref Search ADS Calderón D. , Cuadra J., Schartmann M., Burkert A., Russell C. M. P., 2020b , ApJ , 888 , 5pp 10.3847/2041-8213/ab5e81 Crossref Search ADS Capriotti E. , Kozminski J., 2001 , PASP , 113 , 677 10.1086/320809 Crossref Search ADS Castor J. , Weaver R., McCray R., 1975 , ApJ , 200 , L107 10.1086/181908 Crossref Search ADS Crowther P. A. , 2007 , ARA&A , 45 , 177 10.1146/annurev.astro.45.051806.110615 Crossref Search ADS Crowther P. A. , Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P., Kassim H. A., 2010 , MNRAS , 408 , 731 10.1111/j.1365-2966.2010.17167.x Crossref Search ADS Crowther P. A. et al. , 2016 , MNRAS , 458 , 624 10.1093/mnras/stw273 Crossref Search ADS Dale J. E. , 2015 , Phys. Rev. B , 68 , 1 10.1103/PhysRevB.92.184504 Crossref Search ADS Dale J. E. , Ngoumou J., Ercolano B., Bonnell I. A., 2014 , MNRAS , 442 , 694 10.1093/mnras/stu816 Crossref Search ADS Dashyan G. , Dubois Y., 2020 , A&A , 638 , A123 10.1051/0004-6361/201936339 Crossref Search ADS Ekström S. et al. , 2012 , A&A , 537 , A146 10.1051/0004-6361/201117751 Crossref Search ADS Ferland G. J. , 2003 , ARA&A , 41 , 517 10.1146/annurev.astro.41.011802.094836 Crossref Search ADS Fierlinger K. M. , Burkert A., Ntormousi E., Fierlinger P., Schartmann M., Ballone A., Krause M. G. H., Diehl R., 2016 , MNRAS , 456 , 710 10.1093/mnras/stv2699 Crossref Search ADS Franco J. , Tenorio-Tagle G., Bodenheimer P., 1990 , ApJ , 349 , 126 10.1086/168300 Crossref Search ADS Fromang S. , Hennebelle P., Teyssier R., 2006 , A&A , 457 , 371 10.1051/0004-6361:20065371 Crossref Search ADS Fukushima H. , Yajima H., Sugimura K., Hosokawa T., Omukai K., Matsumoto T., 2020 , MNRAS , 497 , 3830 10.1093/mnras/staa2062 Crossref Search ADS Gatto A. et al. , 2017 , MNRAS , 466 , 1903 10.1093/mnras/stw3209 Crossref Search ADS Geen S. , Hennebelle P., Tremblin P., Rosdahl J., 2015 , MNRAS , 454 , 4484 10.1093/mnras/stv2272 Crossref Search ADS Geen S. , Hennebelle P., Tremblin P., Rosdahl J., 2016 , MNRAS , 463 , 3129 10.1093/mnras/stw2235 Crossref Search ADS Geen S. , Soler J. D., Hennebelle P., 2017 , MNRAS , 471 , 4844 10.1093/mnras/stx1765 Crossref Search ADS Geen S. , Watson S. K., Rosdahl J., Bieri R., Klessen R. S., Hennebelle P., 2018 , MNRAS , 481 , 2548 10.1093/mnras/sty2439 Crossref Search ADS Geen S. , Pellegrini E., Bieri R., Klessen R., 2020 , MNRAS , 492 , 915 10.1093/mnras/stz3491 Crossref Search ADS Gentry E. S. , Krumholz M. R., Dekel A., Madau P., 2017 , MNRAS , 465 , 2471 10.1093/mnras/stw2746 Crossref Search ADS Georgy C. , Ekström S., Meynet G., Massey P., Levesque E. M., Hirschi R., Eggenberger P., Maeder A., 2012 , A&A , 542 , 19 10.1051/0004-6361/201118340 Crossref Search ADS Götberg Y. , de Mink S., Groh J., Leitherer C., Norman C., 2019 , A&A , 629 , A134 10.1051/0004-6361/201834525 Crossref Search ADS Green S. , Mackey J., Haworth T. J., Gvaramadze V. V., Duffy P., 2019 , A&A , 625 , A4 10.1051/0004-6361/201834832 Crossref Search ADS Groenewegen M. , Lamers H., 1989 , A&AS , 79 , 359 Grudić M. Y. , Kruijssen J. D., Faucher-Giguère C.-A., Hopkins P. F., Ma X., Quataert E., Boylan-Kolchin M., 2020 , preprint (arXiv:2008.04453) Guedel M. , Briggs K. R., Montmerle T., Audard M., Rebull L., Skinner S. L., 2007 , Science , 319 , 309 10.1126/science.1149926 Crossref Search ADS PubMed Gurvich A. B. et al. , 2020 , MNRAS , 498 , 3664 10.1093/mnras/staa2578 Crossref Search ADS Haid S. , Walch S., Seifried D., Wünsch R., Dinnbier F., Naab T., 2018 , MNRAS , 478 , 4799 10.1093/mnras/sty1315 Crossref Search ADS Harper-Clark E. , Murray N., 2009 , ApJ , 693 , 1696 10.1088/0004-637X/693/2/1696 Crossref Search ADS Hennebelle P. , 2013 , A&A , 556 , A153 10.1051/0004-6361/201321292 Crossref Search ADS Hosokawa T. , Omukai K., 2009 , ApJ , 703 , 1810 10.1088/0004-637X/703/2/1810 Crossref Search ADS Kahn F. D. , 1954 , Bull. Astron. Inst. Netherlands , 12 , 187 Keszthelyi Z. , Puls J., Wade G., 2017 , A&A , 598 , 20 10.1051/0004-6361/201629468 Crossref Search ADS Kim J.-G. , Kim W.-T., Ostriker E. C., 2018 , ApJ , 859 , 68 10.3847/1538-4357/aabe27 Crossref Search ADS Klessen R. S. , Heitsch F., Mac Low M.-M., 2000 , ApJ , 535 , 887 10.1086/308891 Crossref Search ADS Koyama H. , Inutsuka S.-I., 2004 , ApJ , 602 , L25 10.1086/382478 Crossref Search ADS Krtička J. , Kubát J., 2004 , A&A , 417 , 1003 10.1051/0004-6361:20034030 Crossref Search ADS Kruijssen J. M. D. et al. , 2019 , Nature , 569 , 519 10.1038/s41586-019-1194-3 Crossref Search ADS PubMed Krumholz M. R. , Matzner C. D., 2009 , ApJ , 703 , 1352 10.1088/0004-637X/703/2/1352 Crossref Search ADS Kuiper R. , Hosokawa T., 2018 , A&A , 616 , A101 10.1051/0004-6361/201832638 Crossref Search ADS Lamers H. J. G. L. M. , Leitherer C., Lamers H. J. G. L. M., Leitherer C., 1993 , ApJ , 412 , 771 10.1086/172960 Crossref Search ADS Lee Y.-N. , Hennebelle P., 2016 , A&A , 30 , 1 10.1051/0004-6361/201527981 Crossref Search ADS Lee Y.-N. , Hennebelle P., 2018 , A&A , 611 , A89 10.1051/0004-6361/201731523 Crossref Search ADS Leitherer C. , Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014 , ApJS , 212 , 14 10.1088/0067-0049/212/1/14 Crossref Search ADS Mac Low M.-M. , McCray R., 1988 , ApJ , 324 , 776 10.1086/165936 Crossref Search ADS Mackey J. , Langer N., Gvaramadze V. V., 2013 , MNRAS , 436 , 859 10.1093/mnras/stt1621 Crossref Search ADS Mackey J. , Gvaramadze V. V., Mohamed S., Langer N., 2015 , A&A , 573 , 14 10.1051/0004-6361/201424716 Crossref Search ADS Mathews W. G. , 1967 , ApJ , 147 , 965 10.1086/149087 Crossref Search ADS McKee C. F. , Ostriker J. P., 1977 , ApJ , 218 , 148 Crossref Search ADS Muijres L. , Vink J. S., de Koter A., Müller P., Langer N., 2012 , A&A , 537 , A37 10.1051/0004-6361/201015818 Crossref Search ADS Oort J. H. , Spitzer L., 1955 , ApJ , 121 , 6 10.1086/145958 Crossref Search ADS Pabst C. et al. , 2019 , Nature , 565 , 618 10.1038/s41586-018-0844-1 Crossref Search ADS PubMed Pellegrini E. W. et al. , 2007 , ApJ , 658 , 1119 10.1086/511258 Crossref Search ADS Pellegrini E. W. , Baldwin J. A., Ferland G. J., 2011 , ApJ , 738 , 20 10.1088/0004-637X/738/1/34 Crossref Search ADS Pellegrini E. W. , Oey M. S., Winkler P. F., Points S. D., Smith R. C., Jaskot A. E., Zastrow J., 2012 , ApJ , 755 10.1088/0004-637X/755/1/40 Crossref Pellegrini E. W. , Rahner D., Reissl S., Glover S. C. O., Klessen R. S., Rousseau-Nepton L., Herrera-Camus R., 2020 , MNRAS , 496 , 339 10.1093/mnras/staa1473 Crossref Search ADS Peters T. , Banerjee R., Klessen R. S., Low M.-M. M., Galván-Madrid R., Keto E. R., 2010 , ApJ , 711 , 1017 10.1088/0004-637X/711/2/1017 Crossref Search ADS Prinja R. K. , Barlow M., Howarth I. D., 1990 , ApJ , 361 , 607 10.1086/169224 Crossref Search ADS Puls J. , Sundqvist J. O., Markova N., 2015 , in Meynet G., Georgy C., Groh J. H., Stee Ph., eds, Proc. IAU Symp 307, New Windows on Massive Stars: Asteroseismology, Interferometry, and Spectropolarimetry . Cambridge Univ. Press , Cambridge , p. 25 10.1017/S174392131400622X Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Rahner D. , Pellegrini E. W., Glover S. C. O., Klessen R. S., 2017 , MNRAS , 470 , 4453 10.1093/mnras/stx1532 Crossref Search ADS Rey-Raposo R. , Dobbs C., Agertz O., Alig C., 2017 , MNRAS , 464 , 3536 10.1093/mnras/stw2607 Crossref Search ADS Rogers H. , Pittard J. M., 2013 , MNRAS , 431 , 1337 10.1093/mnras/stt255 Crossref Search ADS Rosdahl J. , Teyssier R., 2015 , MNRAS , 449 , 4380 10.1093/mnras/stv567 Crossref Search ADS Rosdahl J. , Blaizot J., Aubert D., Stranex T., Teyssier R., 2013 , MNRAS , 436 , 2188 10.1093/mnras/stt1722 Crossref Search ADS Rosdahl J. et al. , 2018 , MNRAS , 479 , 994 10.1093/mnras/sty1655 Crossref Search ADS Sana H. et al. , 2012 , Science , 337 , 444 10.1126/science.1223344 Crossref Search ADS PubMed Silich S. , Tenorio-Tagle G., 2013 , ApJ , 765 , 43 10.1088/0004-637X/765/1/43 Crossref Search ADS Silich S. , Tenorio-Tagle G., 2017 , MNRAS , 465 , 1375 10.1093/mnras/stw2879 Crossref Search ADS Smith N. , 2014 , ARA&A , 52 , 487 10.1146/annurev-astro-081913-040025 Crossref Search ADS Spitzer L. , 1978 , Physical Processes in the Interstellar Medium . Wiley-Interscience , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sutherland R. S. , Dopita M. A., 1993 , ApJS , 88 , 253 10.1086/191823 Crossref Search ADS Tenorio-Tagle G. , 1979 , A&A , 71 , 59 Teyssier R. , 2002 , A&A , 385 , 337 10.1051/0004-6361:20011817 Crossref Search ADS van der Werf P. P. , Goss W., O’Dell C., 2013 , ApJ , 762 , 101 10.1088/0004-637X/762/2/101 Crossref Search ADS Veilleux S. , Cecil G., Bland-Hawthorn J., 2005 , ARA&A , 43 , 769 10.1146/annurev.astro.43.072103.150610 Crossref Search ADS Vink J. S. , Muijres L. E., Anthonisse B., de Koter A., Graefener G., Langer N., 2011 , A&A , 531 , 132 10.1051/0004-6361/201116614 Crossref Search ADS Wade G. A. , 2012 , in Drissen L., Robert C., St-Louis N., Moffat A. F. J., eds, ASP Conf. Ser. Vol. 465, Four Decades of Research on Massive Stars . Astron. Soc. Pac , San Francisco , p. 33 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wadepuhl M. , Springel V., 2011 , MNRAS , 410 , 1975 10.1111/j.1365-2966.2010.17576.x Crossref Search ADS Wall J. E. , McMillan S. L. W., Mac Low M.-M., Klessen R. S., Zwart S. P., 2020 , ApJ , 904 , 192 10.3847/1538-4357/abc011 Crossref Search ADS Weaver R. , McCray R., Castor J., Shapiro P., Moore R., 1977 , ApJ , 218 , 377 10.1086/155692 Crossref Search ADS Weidner C. , Kroupa P., Bonnell I., 2009 , MNRAS , 401 , 275 10.1111/j.1365-2966.2009.15633.x Crossref Search ADS Whitworth A. , 1979 , MNRAS , 186 , 59 10.1093/mnras/186.1.59 Crossref Search ADS Yeh S. C. C. , Matzner C. D., 2012 , ApJ , 757 , 18 10.1088/0004-637X/757/2/108 Crossref Search ADS Yeh S. C. C. , Verdolini S., Krumholz M. R., Matzner C. D., Tielens A. G. G. M., 2013 , ApJ , 769 , 7 10.1088/0004-637X/769/1/11 Crossref Search ADS APPENDIX A: STELLAR TRACKS In this section, we plot radiation and wind properties for each of the stellar tracks featured in this paper (see Section 2). In Fig. A1, we plot cumulative wind energy output, wind luminosity, wind mass loss, and photon emission rates binned in photon energy bands corresponding to the ionizing continua of H i(> 13.6 eV), He i(> 24.6 eV), and He ii(> 54.2). Figure A1. Open in new tabDownload slide Wind and radiation outputs for each of the stellar tracks featured in this paper (30, 60, and 120 M⊙ stars). The left-hand column shows wind properties: from top to bottom, cumulative energy output from winds, wind luminosity (⁠|$1/2 \dot{m}_w v_w^2$|⁠), and wind mass-loss rate. The right-hand column shows photon emission rate in groups bounded by the ionization energy of, from top to bottom, H ii, He ii, and He iii. See Section 2 for more details. Figure A1. Open in new tabDownload slide Wind and radiation outputs for each of the stellar tracks featured in this paper (30, 60, and 120 M⊙ stars). The left-hand column shows wind properties: from top to bottom, cumulative energy output from winds, wind luminosity (⁠|$1/2 \dot{m}_w v_w^2$|⁠), and wind mass-loss rate. The right-hand column shows photon emission rate in groups bounded by the ionization energy of, from top to bottom, H ii, He ii, and He iii. See Section 2 for more details. © 2020 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - The geometry and dynamical role of stellar wind bubbles in photoionized H ii regions
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/staa3705
DA - 2021-01-05
UR - https://www.deepdyve.com/lp/oxford-university-press/the-geometry-and-dynamical-role-of-stellar-wind-bubbles-in-zyIkplTBFa
SP - 1352
EP - 1369
VL - 501
IS - 1
DP - DeepDyve
ER -