# Bottling the Champagne: Dynamics and Radiation Trapping of Wind-Driven
Bubbles around Massive Stars

Sam Geen1 , Alex de Koter1,2
1 Anton Pannekoek Institute for Astronomy, Universiteit van Amsterdam, Science
Park 904, 1098 XH Amsterdam, The Netherlands
2 Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven,
Belgium
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this paper we make predictions for the behaviour of wind bubbles around
young massive stars using analytic theory. We do this in order to determine
why there is a discrepancy between theoretical models that predict that winds
should play a secondary role to photoionisation in the dynamics of Hii
regions, and observations of young Hii regions that seem to suggest a driving
role for winds. In particular, regions such as M42 in Orion have neutral
hydrogen shells, suggesting that the ionising radiation is trapped closer to
the star. We first derive formulae for wind bubble evolution in non-uniform
density fields, focusing on singular isothermal sphere density fields with a
power law index of -2. We find that a classical “Weaver”-like expansion
velocity becomes constant in such a density distribution. We then calculate
the structure of the photoionised shell around such wind bubbles, and
determine at what point the mass in the shell cannot absorb all of the
ionising photons emitted by the star, causing an “overflow” of ionising
radiation. We also estimate perturbations from cooling, gravity, magnetic
fields and instabilities, all of which we argue are secondary effects for the
conditions studied here. Our wind-driven model provides a consistent
explanation for the behaviour of M42 to within the errors given by
observational studies. We find that in relatively denser molecular cloud
environments around single young stellar sources, champagne flows are unlikely
until the wind shell breaks up due to turbulence or clumping in the cloud.

###### keywords:

stars: massive – stars: formation $<$ Stars – ISM: H ii regions – ISM: clouds
$<$ Interstellar Medium (ISM), Nebulae – methods: analytical $<$ Astronomical
instrumentation, methods, and techniques

††pubyear: 2021††pagerange: Bottling the Champagne: Dynamics and Radiation
Trapping of Wind-Driven Bubbles around Massive Stars–C

##  1 Introduction

Ionising radiation from stars plays an important role in the interstellar
medium. Lyman continuum photons leave the star and photoionise the gas around
them in a region bounded by an ionisation front, which separates this
photoionised gas and the neutral gas outside. These volumes of photoionised
gas are called Hii regions.

This photoionised gas recombines to neutral hydrogen over time, and hence a
certain number of ionising photons are required to keep the gas ionised. The
flux of photons reaching the ionisation front is also reduced by a geometric
dilution factor, which describes the fact photons leaving the surface of a
star in radial directions are spread across increasingly large spherical
shells. The ionisation front travels outwards if there are any remaining
ionising photons reaching the front. Since the photoionised gas is warmer
($\sim 10^{4}~{}$K) than the neutral gas outside (10-1000 K), this also causes
thermal expansion of the Hii region. This expansion reduces the density of the
photoionised gas, lowering its recombination rate and hence allowing ionising
photons to reach larger radii before being absorbed, and hence the ionisation
front moves outwards as the Hii region expands. The thermal expansion produces
a leading shock wave that accelerates and condenses the surrounding gas into a
dense shell of partially or wholly neutral gas. Models including these
processes have been produced by, e.g., Kahn (1954), Spitzer (1978) and Dyson &
Williams (1980).

In a uniform medium with densities similar to those found in star-forming
molecular clouds, the behaviour of these solutions is typically for the
ionisation front to move outwards rapidly over the recombination time in the
gas, until none of the ionising photons arrive at the ionisation front, and so
the Hii region reaches photoionisation equilibrium, where the total
recombination rate balances the ionising photon emission rate from the star.
The recombination time of dense gas inside molecular clouds is typically very
short compared to the lifetime of massive stars that produce significant
quantities of ionising radiation. After photoionisation equilibrium is
achieved, the Hii region grows primarily due to thermal expansion at speeds
slower than the sound speed in the photoionised gas.

However, if the ionisation front moves into regions with lower gas densities,
such as towards the outskirts of a molecular cloud, the outward movement of
the ionisation front can accelerate since the more diffuse gas is less
efficient at absorbing the ionising photons than the denser cloud material. In
some cases, the Hii region leaves photoionisation equilibrium and the
ionisation front travels outwards supersonically. This mode of Hii region
evolution has been referred to as a “champagne” flow by Tenorio-Tagle (1979),
who invoke a step function in density to model a flow leaving a dense cloud
and entering the diffuse medium outside.

Comeron (1997) produces analytic descriptions and 2D numerical simulations of
such a case for a pair of massive stars embedded in the denser part of the
step function, including both photoionisation and stellar winds. They find
that winds produce complex structures as the wind bubble enters the diffuse
medium, but that their presence does not significantly affect the evolution of
a rapid champagne flow driven by photoionisation into the diffuse medium.

Alternatively, Franco et al. (1990) show that a champagne flow occurs if the
density gradient in the neutral cloud has a power law index steeper than
$-3/2$, (i.e. one in which the density halves when the radius increases by a
factor of roughly 1.6). In this mode, the Hii region never reaches
photoionisation equilibrium. Simulations of star formation predict that stars
are born in density peaks in clouds with a power law index around $-2$ (e.g.
Bate, 2012; Lee & Hennebelle, 2018), and so are likely to form champagne flows
as ionising radiation escapes the protostellar environment around the young
star.

Simulations of feedback in molecular clouds with self-consistent star
formation frequently reproduce such champagne flows. For example, Dale et al.
(2012), Ali et al. (2018) and Zamora-Avilés et al. (2019) include ionising
radiation in their simulations and report champagne flows resulting from such
conditions. Ali et al. (2018) target the Orion nebula, and find some
similarities with the observed region, including far ultra-violet (FUV)
emission and expansion of gas above 10 km/s. Simulations including winds and
photoionising radiation from massive stars, such as Dale et al. (2014) and
Geen et al. (2021) also report similar flows, with winds embedded inside them.
Geen et al. (2015) find that magnetic fields can constrain the fragmentation
of the shells around Hii regions and hence the escape of ionising radiation to
some extent, although even with magnetic fields, ionising radiation and warm
($\sim 10^{4}~{}$K) gas drives rapid expansion of champagne-like flows into
the external medium.

However, observers such as Guedel et al. (2008) and Pabst et al. (2019) find
evidence that nearby Hii regions such as the Orion Nebula are filled with
x-ray emitting gas surrounded by a rapidly-expanding neutral shell. This is
characteristic for stellar winds, which are emitted from the star at thousands
of km/s (Lamers et al., 1993; Vink et al., 2011) and subsequently shock-heats
the gas around the star to millions of degrees (e.g. Weaver et al., 1977;
Dunne et al., 2003). Since the shell around this wind bubble is neutral, it
implies that any ionising radiation from the star is trapped by the stellar
wind bubble. The hypothesis is that the stars in such regions are not
producing champagne flows driven by ionising radiation despite having the
conditions for such flows according to the previously described theoretical
models.

The purpose of this work is to use algebraic prescriptions to identify how
such systems should behave and provide a stepping stone between observed
behaviour of feedback structures around young massive stars, and more detailed
numerical models. We apply analytic theory to power law density fields where
champagne flows can occur (e.g. Franco et al., 1990). We first give an
overview of the physical systems modelled. We then derive a classical “Weaver-
like” solution (Weaver et al., 1977) for an expanding wind bubble in power law
density fields. In the next Section, we introduce semi-analytic and algebraic
solutions for the evolution of the photoionised shell around the wind bubble
and the point at which ionising photons can overflow the neutral shell. We
then discuss possible perturbations to this model, and compare our results to
M42 in Orion.

The key argument of this paper is that there exists a set of solutions for a
wind-driven bubble expanding into a power law density field where the ionising
radiation remains trapped, preventing a “champagne” flow as described by
Franco et al. (1990). Mathematically, there is a point where this trapping
should end. For typical conditions in molecular clouds around individual
massive stars, our model predicts that any possible champagne flow is likely
to remain “bottled” by the shell around the wind bubble for the duration of
the model’s applicability, unless the cloud is sufficiently diffuse. Other 3D
environmental effects are more likely to allow ionising radiation to break out
from the cloud, such as nearby wind bubbles merging (e.g. Calderón et al.,
2020), the density profile in the cloud flattening out as the wind expands
into the wider cloud environment, or a step function to a much lower density
outside the cloud as described by Tenorio-Tagle (1979).

For the rest of this Section we review the literature regarding the modelling
of photoionised Hii regions with embedded wind bubbles, and lay out the
structure of the rest of this paper.

###  1.1 Theory of Stellar Wind Feedback

Early theoretical work by Avedisova (1972), Castor et al. (1975) and Weaver et
al. (1977) established a basis for understanding adiabatic wind bubbles, as
well as initial work on the behaviour of the bubbles in response to radiative
cooling and ionising radiation.

More recent work by Capriotti & Kozminski (2001) looked in more detail at the
dynamical evolution of wind bubbles inside photoionised Hii regions, and
argued that wind bubbles play a secondary role to photoionisation. Haid et al.
(2018) found similar results in controlled hydrodynamic simulations in a
uniform medium. This was also the conclusion of Geen et al. (2020) in a power
law density field where an efficiently cooling wind bubble was invoked inside
an existing photoionised Hii region.

However, in conditions where wind bubbles retain most of their energy, they
can efficiently drive expansion of hot ($>~{}10^{6}~{}$K) bubbles around
massive stars. Whether or not wind bubbles retain most or all of the energy
deposited by the star is an open question. Fierlinger et al. (2016) argue that
properly resolving the contact discontinuity between the hot, diffuse wind
bubble and the denser shell around it is critical in setting the energetics of
the wind bubble. Gentry et al. (2017) argue that the same is true for
supernovae, and that numerical diffusion in grid hydrodynamic codes can cause
artificial energy losses. Fielding et al. (2020) and Tan et al. (2021) argue
that turbulent mixing can increase the cooling rate of the wind bubble,
provided the turbulence is strong enough. Increased wind pressure also
increases the density of the shell around the wind bubble, more effectively
trapping the ionising radiation.

The role of radiation pressure in Hii regions has been explored in models by
Mathews (1967), Krumholz & Matzner (2009), Murray et al. (2010), Draine (2011)
and Kim et al. (2016), both with and without an embedded wind bubble. In the
model of Pellegrini et al. (2007), which includes an extended wind bubble,
radiation pressure further reduces the efficiency of photoionisation feedback
by pushing the ionised gas into a thinner shell, where it recombines faster.
This analysis has been further applied to observed feedback structures around
massive stars and clusters with an embedded wind bubble in Pellegrini et al.
(2009), Yeh & Matzner (2012), Verdolini et al. (2013), Rahner et al. (2017)
and Pellegrini et al. (2020).

How winds and radiation interact is a timely problem due to the recent
inclusion of stellar winds in self-consistent
(radiative-)(magneto-)hydrodynamic simulations of star formation on a cloud
scale. Rogers & Pittard (2013) simulated stellar winds from a small cluster of
O stars in an inhomogeneous cloud. Meanwhile, Dale et al. (2014) included
winds in simulations of stars forming in turbulent clouds with
photoionisation. Since then, Wall et al. (2019), Wall et al. (2020), Decataldo
et al. (2020), Geen et al. (2021), and Grudić et al. (2020) have all included
stellar winds in star formation simulations, as well as numerous other works
studying winds from single sources or in idealised settings (e.g. Gallegos-
Garcia et al., 2020) or on larger scales (e.g. Agertz et al., 2013; Gatto et
al., 2017). To facilitate further development of such studies, more careful
work is needed to determine whether simulations properly resolve the
interaction between winds, radiation and star-forming clouds.

Simulations of stellar winds are significantly more expensive than simulations
with just photoionisation due to the characteristic speeds and temperatures,
and the subsequent effect on the simulation timestep. The sound speed in
photoionised gas is typically 10 km/s (e.g. Oort & Spitzer, 1955), whereas
stellar winds have a terminal velocity of up to around 1% of the speed of
light (Groenewegen et al., 1989; Lamers et al., 1993; Howarth et al., 1997;
Massey et al., 2005). Given this, careful work is required to explore the
parameter space of the influence of winds, and ensure that simulations
accurately capture the interaction between each of the processes that occur
inside a star-forming cloud. For example, in Geen et al. (2021), it was noted
that the behaviour of the wind bubble had some differences with the Orion
nebula described in Pabst et al. (2019), although there were similarities
between the conditions in the two systems. Establishing a firm basis for
feedback physics requires a convergence between complex and costly numerical
simulations and observations, for which analytic theory is a valuable stepping
stone.

##  2 Overview

Figure 1: Diagram showing a schematic of the model used in this paper. 1) A
hot ($>10^{6}~{}$K) wind bubble (pink) expands from a stellar source at $r=0$
into a conic section subtended by a solid angle $\Omega$ out to radius
$r_{\mathrm{w}}$. If $\Omega<4\pi$, a wall of dense gas blocks its expansion
in all other directions (blue, left part of the diagram). For the most part we
ignore the “‘free-streaming” phase described in Weaver et al. (1977) since it
does not affect our result. The wind bubble expands into a neutral external
density field with a power law distribution
($\rho=\rho_{\mathrm{0}}(r/1\mathrm{pc})^{-2}$, light blue gradient) and
temperature $\sim 10-100~{}$K. A dense shell containing neutral hydrogen forms
between $r_{\mathrm{i}}$ and $r_{\mathrm{s}}$ from the material swept up by
the wind. The inner part of the shell (mauve) is photoionised by radiation
from the star to $\sim 10^{4}~{}$K between $r_{\mathrm{w}}$ and
$r_{\mathrm{i}}$, while the outer part (blue) between $r_{\mathrm{i}}$ and
$r_{\mathrm{s}}$ remains neutral at a similar temperature to the external
medium. 2) As the wind bubble expands, more of the shell becomes photoionised
until a point where the whole shell is photoionised. 3) The ionising radiation
overflows the shell and creates a rapidly-expanding “champagne” flow into the
external medium. In this paper we argue that the transition to phase 3 is
unlikely in dense molecular cloud conditions until the wind shell is disrupted
by some other effect, such as encountering another wind bubble or being broken
up by asphericities in the cloud.

The goal of this paper is to describe the behaviour of feedback structures
around young massive stars. Initially, we invoke a wind bubble around the star
expanding in a power law density field. Since the medium is dense, it traps
ionising radiation initially, which forms a thin photoionised layer inside the
dense neutral shell swept up around the wind bubble. This is reasonable
because as authors such as Kuiper & Hosokawa (2018) find, the initial
protostellar outflows around stars are not driven by photoionisation, so it
follows that a pre-existing shell trapping ionising radiation can form.

As the dense shell expands, the fraction of the shell that is photoionised
becomes larger, until the whole shell is photoionised. After this point, the
ionising radiation overtakes the shell around the wind bubble and causes an
“overflow” of ionising radiation into the surrounding medium, beginning a
“champagne flow” that rapidly photoionises the unperturbed neutral gas
outside.

We aim to describe (i) the early evolution of the wind bubble in a power law
density field, (ii) the point at which ionising radiation overflow occurs, and
(iii) identify where other effects not included in the main model lead to the
requirement for more complex semi-analytic or numerical models.

In Figure 1 we show a schematic for the structure of the shell around the wind
bubble prior to overflow. The inside of the wind bubble is shock-heated to a
high temperature by the energy from winds injected at high velocity by the
star. Weaver et al. (1977) describes the internal structure of this wind
bubble, which includes a “free-streaming” region in the centre of the wind
bubble where the material from the wind travels outwards at supersonic speeds
before shocking against the gas inside the bubble and creating temperatures >
$10^{6}~{}$K. We return to the internal structure of the bubble later in the
paper where we discuss radiative cooling, but for most of this paper it is
only relevant that the bubble is a hot adiabatic volume of gas.

The interior of the wind bubble is collisionally ionised, and hence the flux
of ionising radiation from the star is not noticeably depleted until it
encounters the inner edge of the dense shell around the wind bubble at
$r_{\mathrm{w}}$ and begins photoionising the gas there. Between
$r_{\mathrm{w}}$ and $r_{\mathrm{i}}$, there is a thin photoionised shell.
This gas has a temperature of roughly $10^{4}~{}$K, determined by the
equilibrium of cooling and heating of the photionised gas.

In the absence of radiation pressure, this region has a uniform density and is
supported by thermal pressure alone. With the inclusion of radiation pressure,
there is a pressure differential inside the region due to the larger surface
area at larger radii and the presence of dust.

Finally, between $r_{\mathrm{i}}$ and $r_{\mathrm{s}}$, there is a cold, dense
neutral shell of material swept up by the wind bubble. At the point of
overflow, the quantity of the shell that remains neutral tends to zero, at
which point ionising radiation can enter the neutral power law density field
outside.

We note that this picture applies for steeper power law density fields as
found around young massive stars. In uniform density fields, Comeron (1997)
and Silich & Tenorio-Tagle (2013) note that overflow becomes less possible the
larger the wind bubble becomes. They argue that as a shell around the wind
bubble grows, the Hii region reaches a “trapping” point where ionising
radiation emitted by the star becomes trapped by the shell around the wind
bubble. For much of this paper we thus adopt a fiducial “singular isothermal
sphere” density field with a power law index of $-2$, which is typical for the
conditions around very young massive stars sitting in star-forming cores that
have just ended the protostellar phase and is adopted by other authors such as
Shu et al. (2002).

The models in this paper also allow for a wind bubble that is completely
constrained by dense gas along certain lines of sight, as in a blister region
or a case where outflow occurs only in a small solid angle around the star. In
this case we invoke a solid angle subtended by the wind bubble $\Omega$, which
tends towards $4\pi$ as the wind bubble expands in all directions into a
spherically symmetric density field.

###  2.1 Initial Environment around the star

Due to the large parameter space of possible environments into which feedback
from young stars can evolve, we confine ourselves to a set of physically
motivated quasi-1D models for the conditions for young massive stellar objects
in star-forming cores. A similar analysis can be applied to other conditions,
such as a moving source or a density step function for flows breaking out of a
denser cloud environment.

We define the young massive star as sitting in a cold neutral distribution of
gas with a spherically symmetric power law density distribution given by

$n(r)=n_{\mathrm{0}}(r/r_{\mathrm{0}})^{-\omega}$ (1)

where $n(r)$ is the hydrogen number density at a radius $r$, $n_{\mathrm{0}}$
is the density at $r_{\mathrm{0}}$ and $\omega$ is the power law index.
$\omega=0$ gives a uniform density field, and $\omega=2$ is described as a
singular isothermal sphere.

We consider this distribution over an open angle $\Omega$, with a maximum
value of $4\pi$ if the density distribution extends across all lines of sight.
Smaller values can be found if much denser gas constrains the bubble along
certain lines of sight. The mass enclosed in radius $r$ in the initial
conditions is thus

$M(<r)=\int_{\mathrm{0}}^{r}\Omega r^{\prime
2}n(r^{\prime})\frac{m_{\mathrm{H}}}{X}.dr^{\prime}=\frac{\Omega}{3-\omega}n_{\mathrm{0}}\,r_{\mathrm{0}}^{\omega}\,r^{3-\omega}\frac{m_{\mathrm{H}}}{X}$
(2)

where $m_{\mathrm{H}}/X$ is conversion factor from hydrogen number density to
mass density, and $X$ is the hydrogen mass fraction in the gas. We also define
a characteristic mass density $\rho_{\mathrm{0}}\equiv
n_{\mathrm{0}}m_{\mathrm{H}}/X$.

Franco et al. (1990) use this distribution above a threshold “core” radius of
$r_{\mathrm{c}}$, which defines the edge of the protostellar core. For
$r<r_{\mathrm{c}}$ they assume that the density is constant, approximating a
protostellar core. For the sake of simplicity, we neglect this core radius in
further calculations in this paper. We justify this by noting that simulations
of protostellar core formation by Lee & Hennebelle (2018) find radii where the
density profile begins to flatten at radii of around $20~{}\mathrm{AU}\simeq
10^{-4}~{}\mathrm{pc}$. This is much smaller than the scales studied here (on
the order of parsecs).

However, the presence of a small core at the point the star is formed remains
an important “initial condition” in our assumptions due to the fact that in a
pure power law density field, ionising radiation would quickly photoionise the
whole cloud as in Franco et al. (1990). We also note that in reality, main
sequence stellar feedback is preceded by feedback from protostellar outflows,
where photoionisation is a secondary process (Kuiper & Hosokawa, 2018).
Therefore, in a realistic physical context, winds and ionising photons from
the star will be emitted into a medium initially preprocessed by a kinetic
outflow rather than a pristine neutral gas field.

###  2.2 Stellar Models

Throughout this work, we use the same stellar evolution tables as Geen et al.
(2021). These are based on Ekström et al. (2012) for the stellar parameters,
using Leitherer et al. (2014) to give the stellar spectra and a correction
from Vink et al. (2011) to convert the escape velocity of the star into a wind
terminal velocity. Relevant physical values for each star are listed in
Appendix C and Tables 2 and 3.

##  3 Stellar Wind Bubbles

In this Section we discuss the evolution of a wind bubble around a young
massive star, and the role that environment plays in shaping its evolution.

###  3.1 Revisiting Adiabatic Wind Bubbles

Weaver et al. (1977) give a solution for the expansion of an adiabatic wind
bubble in a uniform medium surrounded by a shell of swept-up interstellar gas.
We now re-derive this solution for an adiabatic wind bubble expanding into a
power law density distribution, and make observations about how such systems
evolve differently to the Weaver et al. (1977) solution.

The dynamical equations governing the wind bubble’s expansion are

$E_{\mathrm{b}}=\frac{3}{2}\frac{\Omega}{3}r_{\mathrm{w}}^{3}P_{\mathrm{w}},$
(3)
$\frac{\mathrm{d}}{\mathrm{d}t}\left(M(<r_{\mathrm{w}})\frac{\mathrm{d}r_{\mathrm{w}}}{\mathrm{d}t}\right)=\Omega
r_{\mathrm{w}}^{2}P_{\mathrm{w}},$ (4)

and

$\frac{\mathrm{dE_{\mathrm{b}}}}{\mathrm{d}t}=L_{\mathrm{w}}-\Omega
r_{\mathrm{w}}^{2}P_{\mathrm{w}}\frac{\mathrm{dr_{\mathrm{w}}}}{\mathrm{d}t},$
(5)

where $E_{\mathrm{b}}$ is the thermal energy in the shocked gas inside the
wind bubble and $P_{\mathrm{w}}$ is the pressure from the wind bubble acting
on the inside of the shell. Equation 4 describes the conservation of momentum,
and Equation 5 describes the conservation of energy. The first term on the
right hand side of this latter equation, $L_{\mathrm{w}}$, is the wind
luminosity of the central source, i.e. the mechanical power of the material
ejected by the star, assumed constant here for simplicity. The second right
hand side term is the work done by the gas pressure $P_{\mathrm{w}}$. The
temperature and density of the wind bubble are not immediately important in
deriving this solution. We give estimates for these quantities in Section 5.3.

To solve these equations, we require the energy stored inside the wind bubble,
which is the fraction of energy emitted by the star not used to do work on the
shell around the wind bubble. We derive a quantity for this in Appendix A,
which is given by

$E_{\mathrm{b}}=\left(\frac{5-\omega}{11-\omega}\right)L_{\mathrm{w}}t$ (6)

where $t$ is the time since the onset of the wind.

To produce a tractable solution, we further assume that
$r_{\mathrm{w}}(t)\propto t^{B}$, where $B$ is some power law index. Solving
for these equations plus Equation 2, we find

$r_{\mathrm{w}}(\omega,t)=\left(A_{\mathrm{w}}(\omega,\Omega)L_{\mathrm{w}}\rho_{\mathrm{0}}^{-1}r_{\mathrm{0}}^{-\omega}t^{3}\right)^{1/(5-\omega)}$
(7)

where

$A_{\mathrm{w}}(\omega,\Omega)=\frac{4\pi}{\Omega}\frac{(1-\omega/3)(1-\omega/5)^{3}}{(1-2\omega/7)(1-\omega/11)}\frac{250}{308\pi}$
(8)

This becomes the “Weaver” solution (equation 21 of Weaver et al., 1977) if
$\Omega=4\pi$ and $\omega=0$. We verify Equations 6 and 7 using simple 1D
hydrodynamic simulations in Appendix B, and find good agreement.

For $\omega=2$ and $\Omega=4\pi$, Equation 7 becomes

$r_{\mathrm{w,2}}(t)=13.9~{}\mathrm{pc}\left(\frac{L_{\mathrm{w}}}{10^{36}~{}\mathrm{erg/s}}\right)^{1/3}\left(\frac{n_{\mathrm{0}}}{1000~{}\mathrm{cm}^{-3}}\right)^{-1/3}\left(\frac{t}{1~{}\mathrm{Myr}}\right),$
(9)

which has a constant expansion velocity

$v_{\mathrm{w,2}}=13.5~{}\mathrm{km/s}~{}\left(\frac{L_{\mathrm{w}}}{10^{36}~{}\mathrm{erg/s}}\right)^{1/3}\left(\frac{n_{\mathrm{0}}}{1000~{}\mathrm{cm}^{-3}}\right)^{-1/3}.$
(10)

This is a useful result for the regions immediately around young massive
stars, since stars typically form in cores where $\omega\simeq 2$ as matter
accretes onto the site of protostar formation. (e.g. Lee & Hennebelle, 2018).
At larger radii where this steep profile merges into the cloud background, we
expect the solution to become more Weaver-like again, or break up due to
clumping and turbulence in the cloud. Where this happens depends on the
specific environment around the star.

###  3.2 Wind Bubble Expansion Rate

In this Section we discuss solutions to Equation 7.

####  3.2.1 Dependence on Environment

Figure 2: Expansion rate of the wind bubble with time around an example star
to demonstrate the role of environment. The left panel shows the density field
power law index $\omega$, which describes how steep the density around the
star is. The right panel shows the solid angle that the wind bubble expands
into, $\Omega$, where $4\pi$ describes a full sphere with the star at the
centre. In the left panel, we fix $\Omega=4\pi$, and in the right panel we fix
$\omega=2$. For both panels we use a star of mass 35 M⊙,
$n_{\mathrm{0}}=4000~{}$cm-3 and $r_{\mathrm{0}}=1~{}$pc.

We begin by studying the effect of the immediate environment around the star
on the behaviour of solutions to Equation 7. The first parameter is the
density field power law index $\omega$, which describes how the density around
the star drops with radius. The second is the solid angle subtended by the
wind bubble, $\Omega$. We discuss characteristic cloud density
$n_{\mathrm{0}}$ below.

In Figure 2 we plot the evolution of an example wind bubble around a rotating
35 M⊙ star at Solar metallicity with
$\mathrm{log}(L_{\mathrm{w}}/\mathrm{erg/s})=35.9$ (see Table 2). We use
values of $n_{\mathrm{0}}=4000~{}$cm-3 and $r_{\mathrm{0}}=1~{}$pc, varying
$\omega$ and $\Omega$. These values are chosen to approximate the conditions
discussed subsequently in the observational comparison in Section 5.5.

The wind velocity is intially higher for shallower power law density fields
(lower $\omega$), although this drops over time as the wind bubble expands and
the mass of a spherical shell of fixed thickness in the external density field
increases. As given in Equation 10, the expansion velocity of a wind bubble is
constant where $\omega=2$. We note that the mass of a spherical shell of fixed
thickness is constant, and thus both the energy in the bubble and the mass
swept up by the wind bubble increase linearly with time. This is a useful
result since it removes any time dependence when comparing with observations
in young cores where $\omega\simeq 2$, while the radius becomes a simple
linear function of the main sequence age of the star.

We also plot the expansion velocity of the wind bubble as a function of the
solid angle subtended by the wind bubble. A solid angle of $4\pi$ means that
the wind bubble is a sphere centred on the star. However, some wind bubbles
expand preferentially in certain directions, being confined by much denser gas
in other directions. As the solid angle decreases, the volume of the wind
bubble also decreases and hence the pressure inside a bubble of a given radius
rises, causing a larger expansion rate. Equation 7 demonstrates that the shell
radius and expansion rate vary with $\Omega^{-1/(5-\omega)}$, so for a power
law density field where $\omega=2$, decreasing $\Omega$ by a factor of 4 (i.e.
confining the wind bubble’s outflow) increases the wind velocity by around
60%.

####  3.2.2 Dependence on Stellar Mass, Metallicity and Cloud Density

Figure 3: Expansion rate of the wind bubble as a function of stellar mass and
characteristic background density $n_{\mathrm{0}}$ at 1 pc where $\omega=2$
and $\Omega=4\pi$. Solid lines show results at a Solar metallicity Z=0.014,
and dashed lines the results at a sub-Solar metallicity Z=0.002. This
expansion rate is constant with time such that
$r_{\mathrm{w,2}}(t)=v_{\mathrm{w,2}}t$. Winds are weaker in the sub-Solar
case and so the bubble expands more slowly.

In Figure 3 we plot the dependence of the wind bubble expansion rate
$v_{\mathrm{w,2}}$ in Equation 10 (i.e. for a power law density profile where
$\omega=2$). Wind luminosity increases with stellar mass and metallicity,
giving faster expansion rates. Meanwhile, the expansion rate drops with
increasing characteristic density. We include this figure as a quick visual
reference for physically relevant conditions.

##  4 Expansion and Overflow of the Photoionised Shell

The dense shell swept up by the wind bubble interacts with the ionising
radiation emitted by the star, and contains an ionised and a neutral component
(see Figure 1). In this Section we use the results of Section 3 to calculate
the structure and evolution of the photoionised component of the shell,
determine at which point the dense shell cannot absorb all of the ionising
radiation from the star, and discuss what happens after this point.

###  4.1 Structure of the Ionised Shell

The dense, wind-swept shell is typically dense enough that the recombination
time of the gas is short, and hence we assume photoionisation equilibrium,
i.e. the photons emitted by the star in the solid angle of the shell are all
absorbed by the shell. We can write the balance of photon emission and
absorption as

$\frac{\Omega}{4\pi}Q_{\mathrm{H}}=\int_{\mathrm{r_{\mathrm{w}}}}^{r_{\mathrm{i}}}\Omega
r^{2}n_{\mathrm{i}}^{2}(r)\alpha_{\mathrm{B}}.dr$ (11)

where $Q_{\mathrm{H}}$ is the ionising photon emission rate,
$n_{\mathrm{i}}(r)$ is the hydrogen number density in the photoionised gas at
a radius $r$ and $\alpha_{\mathrm{B}}$ is the Baker & Menzel (1938) case B
recombination rate of the gas for Hii, which is 2 to $3\times 10^{-13}$ cm3
s-1 at solar metallicity. The exact formula used to calculate
$\alpha_{\mathrm{B}}$ depends on the temperature of the photoionised gas and
is given in Appendix E2 of Rosdahl et al. (2013).

Equation 11 has an algebraic solution if either (i) the shell is sufficiently
thin that $r_{\mathrm{i}}\rightarrow r_{\mathrm{w}}$, or (ii) the density of
the photoionised shell is uniform, i.e. $n_{\mathrm{i}}$ is independent of
$r$. Case (i) becomes valid where photoionisation is negligible compared to
radiation pressure and winds and case (ii) becomes valid where radiation
pressure is negligible. A more general solution, however, requires a numerical
approach.

Numerical solutions for photoionsied shell structures have been used widely in
the literature to calculate the dynamical and observable properties of
photoionised shells, e.g. Pellegrini et al. (2007), Yeh & Matzner (2012) and
Kim et al. (2016). A form of these equations is given in Draine (2011) for
dusty Hii regions, which was used by Martínez-González et al. (2014) to
include a central wind bubble. Draine (2011) gives the conditions for
hydrostatic pressure balance, ionising photon balance and dust absorption
optical depth respectively at radius $r$ inside the photoionised shell as

$n_{\mathrm{i}}\sigma_{\mathrm{d}}\frac{(L_{\mathrm{n}}e^{-\tau}+L_{\mathrm{i}}\phi)}{4\pi
r^{2}c}+\alpha_{\mathrm{B}}n_{\mathrm{i}}^{2}\frac{\langle
h\nu\rangle_{\mathrm{i}}}{c}-\frac{\mathrm{d}}{\mathrm{d}r}\left(2n_{\mathrm{i}}(r)k_{\mathrm{B}}T_{\mathrm{i}}\right)=0,$
(12)
$\frac{\mathrm{d}\phi}{\mathrm{d}r}=-\frac{1}{Q_{\mathrm{H}}}\alpha_{\mathrm{B}}n_{\mathrm{i}}^{2}4\pi
r^{2}-n_{\mathrm{i}}\sigma_{\mathrm{d}}\phi,$ (13)

and

$\frac{\mathrm{d}\tau}{\mathrm{d}r}=n_{\mathrm{i}}\sigma_{\mathrm{d}},$ (14)

where $\sigma_{\mathrm{d}}$ is the dust absorption cross section per hydrogen
nucleon (we use 10${}^{-21}~{}$cm-2 $Z/Z_{\mathrm{solar}}$ as in Draine, 2011,
where $Z$ is the gas metallicity and $Z_{\mathrm{solar}}=0.014$),
$L_{\mathrm{n}}$ and $L_{\mathrm{i}}$ are the luminosities of the star for
non-hydrogen-ionising and hydrogen-ionising photon energies respectively,
$\phi$ is the fraction of ionising photons reaching radius $r$, $c$ is the
speed of light, $\langle h\nu\rangle_{\mathrm{i}}$ is the average energy of
ionising photons (where $h\nu$ > 13.6 eV), $T_{\mathrm{i}}$ is the temperature
of the photoionised gas (where
$2k_{\mathrm{B}}T_{\mathrm{i}}=c_{\mathrm{i}}^{2}m_{\mathrm{H}}/X$) $\simeq
10^{4}~{}$K, and $\tau$ is the dust absorption optical depth.

We solve these equations beginning at $r_{\mathrm{w}}$, where the density
$n_{\mathrm{i}}$ is set by the pressure balance with the wind bubble,

$P_{\mathrm{w}}=P_{\mathrm{i}}\equiv\frac{m_{\mathrm{H}}}{X}n_{\mathrm{i}}c_{\mathrm{i}}^{2}$
(15)

(e.g. Dyson & Williams, 1980), where $P_{\mathrm{w}}$ is given in Equation 5.
We integrate for increasing $r$ until $\phi$ reaches zero, i.e. all the
ionising photons have been absorbed, which we define as the edge of the
photoionised shell $r_{\mathrm{i}}$.

For each value of $r_{\mathrm{w}}$ we use, we calculate the mass of the
photoionised shell

$M_{\mathrm{i}}=\int_{\mathrm{r_{\mathrm{w}}}}^{r_{\mathrm{i}}}\Omega
r^{2}n_{\mathrm{i}}(r)\frac{m_{\mathrm{H}}}{X}.dr,$ (16)

and compare it to the total mass of gas swept up by the Hii region out to
$r_{\mathrm{i}}$, $M(<r_{\mathrm{i}})$ using Equation 2. If the value found
for Equation 16 exceeds $M(<r_{\mathrm{i}})$, the solution breaks down because
there is no more mass in the neutral shell to absorb the ionising photons. We
refer to this stage of the solution as “overflowed”, since ionising photons
are able to overflow the shell. We discuss the consequences of this later in
this Section.

###  4.2 Numerical Solutions for Shell Structure and Overflow

Figure 4: Point at which ionising photon overflow occurs as a function of
initial stellar mass, with varying metallicity and characteristic background
hydrogen number density $n_{\mathrm{0}}$, where
$n(r)=n_{\mathrm{0}}(r~{}/~{}1~{}\mathrm{pc})^{-2}$. The left plot shows the
outer radius of the ionised shell $r_{\mathrm{i}}$ at which overflow occurs,
the right plot shows the time of overflow. The denser the environment, the
denser the shell around the wind bubble and hence the longer it takes for
ionising photons to overflow the shell. Due to weaker winds, overflow occurs
sooner and at smaller radii at lower metallicity. Solid lines show results at
Solar metallicity ($Z=0.014$) and dashed lines show results at sub-Solar
metallicity ($Z=0.002$). Line colour shows characteristic cloud density
$n_{\mathrm{0}}$.

Figure 5: Relative thickness of the photoionised shell
$r_{\mathrm{i}}-r_{\mathrm{w}}$ against the wind bubble radius
$r_{\mathrm{w}}$, plotted as a function of $r_{\mathrm{w}}$ with varying
initial stellar mass and characteristic background hydrogen number density
$n_{\mathrm{0}}$, where $n(r)=n_{\mathrm{0}}(r~{}/~{}1~{}\mathrm{pc})^{-2}$.
The left plot shows results at $Z=0.014$, and the right plot shows results at
$Z=0.002$. Dashed lines show results for a 120 M⊙ star, solid lines show
results for a 30 M⊙ star. Line colour shows characteristic cloud density
$n_{\mathrm{0}}$. The relative thickness of the ionised shell is typically
small until near the point of overflow, which is shown as a circle at the end
of each line. We omit results after 250 pc since this is larger than most
molecular clouds in our Galaxy and nearby galaxies. We also note that
substructure in complex star-forming environments is likely to cause the power
law density distribution to break down at smaller radii than this. In this
scenario, photon breakout will be caused by 3D effects before the Hii region
reaches the overflow point.

In this section we plot a set of solutions to the equations in Section 4.1. In
order to provide a tractable subset of the parameter space of possible inputs,
we focus on cases where $\omega=2$ for a set of physically-motivated
conditions. As in Section 3.2, we use values from Table 2 for the stellar
properties. For these solutions, we require $Q_{\mathrm{H}}$,
$L_{\mathrm{n}}$, $L_{\mathrm{i}}\equiv Q_{\mathrm{H}}\langle
h\nu\rangle_{\mathrm{i}}$ , $T_{\mathrm{i}}$ and $L_{\mathrm{w}}$. The latter
quantity is used to calculate $r_{\mathrm{w}}$, which is the inner boundary of
the photoionised shell. We use fiducial densities
$n_{\mathrm{0}}=10^{2},10^{3}~{}$and $10^{4}~{}$cm-3, where
$n(r)=n_{\mathrm{0}}(r/1~{}$pc$)^{2}$ for the purposes of this analysis.

Figure 4 shows the time and ionised shell radius at which the ionising
radiation overflows the shell (i.e. there is not enough mass in the shell to
absorb all the ionising radiation). There is a small dependence on
metallicity, and a much larger dependence on cloud density. Note that the
lifetime of the stars described here is on the order of 10 Myr, or less in the
case of very massive stars. We thus expect the star to reach the end of its
main sequence before overflow conditions are reached in very dense clouds. By
contrast, for more diffuse clouds, overflow may occur when the star is
younger. We discuss other scenarios outside the scope of the 1D model in
Section 5. Nonetheless, this model suggests that, provided certain plausible
initial conditions are achieved, there is a range of physically relevant
conditions for which the ionising radiation from massive stars is trapped for
some or all of a star’s life, preventing a champagne flow and reducing the
influence of ionising radiation in molecular clouds.

We plot the relative thickness of the ionised shell compared to the radius of
the wind bubble, $(r_{\mathrm{i}}-r_{\mathrm{w}})/r_{\mathrm{w}}$ in Figure 5.
This is initially small compared to the overall size of the wind bubble,
growing towards the point of overflow to be similar in radius, or even larger
for sub-Solar metallicities. The relative thickness is smaller for more
massive stars due to the stronger winds. The relative expansion velocity of
the ionised shell compared to the wind velocity
$(\dot{r}_{\mathrm{i}}-v_{\mathrm{w,2}})/v_{\mathrm{w,2}}$ follows a very
similar relation, having values between 1.5 and 2 times
$(r_{\mathrm{i}}-r_{\mathrm{w}})/r_{\mathrm{w}}$. Therefore, the expansion
rate of the edge of the ionised shell, where C+ emission is located, begins at
around $v_{\mathrm{w,2}}$, but can be up to $2v_{\mathrm{w,2}}$ at the point
of overflow.

The effect of radiation pressure on the structure of the ionised shell is
typically small compared to the influence of winds. We can measure this by
measuring the density gradient in a photoionised region around a star (Draine,
2011; Martínez-González et al., 2014). If radiation pressure is non-
negligible, the density should be constant with radius. This is because the
pressure balance in the photoionised gas can be written as

$P_{\mathrm{i}}(r)=\rho_{\mathrm{i}}(r)c_{\mathrm{i}}^{2}+P_{\mathrm{rad}}(r),$
(17)

where $P_{\mathrm{i}}$ is the total pressure at $r$ in the photoionised gas
and $P_{\mathrm{rad}}$ is the contribution from radiation pressure at $r$. For
a hydrostatic shell, $P_{\mathrm{i}}$ is constant with $r$, as is
$c_{\mathrm{i}}$ due to the cooling equilibrium in photoionised gas.
$P_{\mathrm{rad}}$ drops as the flux of radiation decreases over the radius of
the shell due to geometric dilution and absorption by atoms and dust.
Therefore, in order to maintain a constant $P_{\mathrm{i}}$,
$\rho_{\mathrm{i}}$ must increase with $r$.

Draine (2011) measures the density contrast in a photoionised region by
comparing the root-mean-square (RMS) density, the mean density and the maximum
density, found at the outer edge of the region at $r_{\mathrm{i}}$. The ratio
between the RMS density and the mean density is negligible in all cases in our
model. The maximum density is initially 1-5% higher than the RMS density at
Solar metallicities, dropping towards smaller values for larger wind bubbles
where the radiation pressure is spread across a larger surface area. For
$Z=0.002$, this difference is initially 2-13%, again dropping towards
negligible values at larger radii. Radiation pressure can thus have a second-
order effect on the structure of compact Hii regions, but this importance
drops at larger Hii region radii towards the point of overflow.

###  4.3 Analytic Description

In order to interpret the results in Section 4.2, we produce an approximate
and illustrative analytic solution for the expansion of the photoionised shell
that allows us to more clearly explain the behaviour of the numerical
solutions to Equations 12 to 14 presented in Section 4.2.

####  4.3.1 Ionised Shell Structure

Our first assumption is that radiation pressure causes only a small change in
$n_{\mathrm{i}}$, and thus we assume a constant $n_{\mathrm{i}}$, giving a
solution to Equation 11 that reads

$Q_{\mathrm{H}}=\frac{4\pi}{3}n_{\mathrm{i}}^{2}(r_{\mathrm{i}}^{3}-r_{\mathrm{w}}^{3})\alpha_{\mathrm{B}},$
(18)

noting the dependence on $4\pi$ and not $\Omega$ since the radiation travels
in all directions, not just the solid angle subtended by the wind bubble.

Using Equation 15 and substituting for $P_{\mathrm{w}}$ using Equations 3 and
6, we can write the number density of the photoionised gas as

$n_{\mathrm{i}}=\frac{2}{\Omega
r_{\mathrm{w}}^{3}}\left(\frac{5-\omega}{11-\omega}\right)\frac{X}{m_{\mathrm{H}}}\frac{1}{c_{\mathrm{i}}^{2}}L_{\mathrm{w}}t.$
(19)

Substituting for $n_{\mathrm{i}}$ in Equation 18, we can write

$r_{\mathrm{i}}^{3}=r_{\mathrm{w}}^{3}+\frac{3}{4\pi}\frac{Q_{\mathrm{H}}}{\alpha_{\mathrm{B}}}\left(\frac{\Omega}{2}\left(\frac{11-\omega}{5-\omega}\right)\frac{m_{\mathrm{H}}}{X}\frac{c_{\mathrm{i}}^{2}}{L_{\mathrm{w}}t}\right)^{2}r_{\mathrm{w}}^{6}$
(20)

Substituting for $t$ in Equation 7, we can write this as

$r_{\mathrm{i}}=r_{\mathrm{w}}\left(1+3B_{\mathrm{w}}~{}r_{\mathrm{w}}^{(2\omega-1)/3}\right)^{\frac{1}{3}}$
(21)

where

$B_{\mathrm{w}}=\pi\left(\frac{\Omega}{4\pi}\right)^{2}\left(\frac{11-\omega}{5-\omega}\right)^{2}\left(\frac{m_{\mathrm{H}}}{X}\right)^{2}\frac{Q_{\mathrm{H}}}{\alpha_{\mathrm{B}}}c_{\mathrm{i}}^{4}\left(\frac{A_{\mathrm{w}}(\omega,\Omega)}{L_{\mathrm{w}}^{2}\rho_{\mathrm{0}}r_{\mathrm{0}}^{\omega}}\right)^{2/3}$
(22)

In the thin shell limit (where $r_{\mathrm{i}}\rightarrow r_{\mathrm{w}}$),
this equation can be written using the binomial approximation as

$\Delta r\equiv r_{\mathrm{i}}-r_{\mathrm{w}}\simeq
B_{\mathrm{w}}~{}r_{\mathrm{w}}^{(2\omega-1)/3}$ (23)

Some general comments may be made about the form of this equation. Firstly,
the thickness of the ionised shell increases with $r_{\mathrm{w}}$ for
$\omega>1/2$, and decreases for $\omega<1/2$. Secondly, the thickness
increases with photoionisation-related quantities $Q_{\mathrm{H}}$ and
$c_{\mathrm{i}}$ and decreases with the wind luminosity $L_{\mathrm{w}}$ and
characteristic background density $\rho_{\mathrm{0}}$.

It is also possible to calculate the thickness of the neutral shell from the
remaining swept-up mass that is not photoionised. This also requires a cooling
function to calculate its density from the hydrostatic pressure in the shell,
which in turn allows an estimate of the thickness, as well as including the
effects of FUV radiation pressure on dust in the shell. However, we do not
discuss this here since it does not affect the results of this work.

####  4.3.2 Condition for Overflow

In order to calculate the overflow radius, we make a further simplifying
assumption that the mass swept up by the shell $M(<r_{\mathrm{i}})\simeq
M(<r_{\mathrm{w}})$. If we do not apply this assumption, a cubic equation in
$r_{\mathrm{w}}$ can be found by using Equation 21 that gives a similar
result, although this becomes more difficult to interpret.

Since the density inside the wind bubble is typically very low, the mass in
the dense shell around the wind bubble is assumed to be equal to the swept-up
mass of unshocked gas at $t=0$ in the background medium, i.e.
$M(<r_{\mathrm{w}})$ using Equation 2. In the limit where this mass equals the
mass of the photoionised shell, all of the shell is photoionised. Beyond this
point, the shell cannot absorb all of the ionising photons, which begin
leaking into the surrounding neutral gas. We can write a condition for such
photon “overflow” as

$\frac{\Omega}{3-\omega}n_{\mathrm{0}}r_{\mathrm{0}}^{\omega}r_{\mathrm{w}}^{3-\omega}\frac{m_{\mathrm{H}}}{X}\leq\frac{\Omega}{3}\left(r_{\mathrm{i}}^{3}-r_{\mathrm{w}}^{3}\right)n_{\mathrm{i}}\frac{m_{\mathrm{H}}}{X}$
(24)

Making the same substitutions for $n_{\mathrm{i}}$ and $t$ as in Section 4.1,
we can produce an expression for the wind bubble radius at which overflow is
possible, where

$r_{\mathrm{w}}^{\frac{5-4\omega}{3}}\leq\frac{\Omega}{4\pi}\left(\frac{3-\omega}{2}\right)\left(\frac{11-\omega}{5-\omega}\right)\frac{Q_{\mathrm{H}}c_{\mathrm{i}}^{2}}{\alpha_{\mathrm{B}}}\left(\frac{m_{\mathrm{H}}}{X}\frac{A_{\mathrm{w}}^{1/2}(\omega,\Omega)}{L_{\mathrm{w}}n_{\mathrm{0}}^{2}r_{\mathrm{0}}^{2\omega}}\right)^{2/3}$
(25)

We define an overflow radius $r_{\mathrm{o}}$ equal to the radius
$r_{\mathrm{w}}$ at the limit where overflow is just possible, i.e. the mass
of the photoionised shell is exactly equal to the mass of the swept-up
material.

There are again some major points to note concerning the behaviour of this
equation. The right hand side increases for larger values of $Q_{\mathrm{H}}$
and $c_{\mathrm{i}}$, and decreases for larger values of $L_{\mathrm{w}}$ and
$n_{\mathrm{0}}$. However, the left hand side behaves differently depending on
the value of $\omega$. For $\omega<5/4$, as in an initially uniform density
field where $\omega=0$, overflow becomes less likely with increasing radius,
assuming a dense shell can accumulate at $r_{\mathrm{w}}$ rather than being
dispersed by hydrodynamic mixing. This mode is described in Comeron (1997) and
Silich & Tenorio-Tagle (2013), who argue for a “trapping” radius, where
ionising radiation initially forms an extended photoionised region, but over
time a wind shell grows in the warm photoionised medium at $r_{\mathrm{w}}$
until its density and mass are sufficient to trap the radiation.

However, for a power law density field above this value, such as where
$\omega=2$, the left hand side decreases with radius, and so overflow becomes
more likely as the wind bubble grows.

This is important because for a wind bubble around a young massive star
expanding into a steep power law density field there exists a family of
solutions for which the ionising radiation cannot escape until some overflow
radius. This requires that the initial conditions are such that an instant
champagne flow as described by Franco et al. (1990) is prevented, e.g. by a
small shallow core in the density profile around the star or an initial shell
set up by protostellar outflows that can trap ionising radiation during the
main sequence stage of the star’s evolution. If these conditions are met, this
offers a plausible explanation for why we observe some young Hii regions that
do not exhibit photoionised champagne flows, such as M42 in Orion (Pabst et
al., 2019), and rather appear to be filled with hot, wind-shocked gas (Guedel
et al., 2008) and are surrounded by neutral shells (Pabst et al., 2020).

In the following Section, we discuss perturbations to the model from other
physical effects, and cases where photon escape may be possible before the
overflow radius is reached.

###  4.4 Post-Overflow Conditions

Reaching the overflow radius does not guarantee an immediate and strong
champagne flow. Firstly, just after the point of overflow, while the neutral
shell and any observational tracers from it disappear, only a small number of
ionising photons will leak into the surrounding medium, causing only a weak
decoupling of the thermally-expanding Hii region and ionisation front.
Initially, the shell does not respond to the change in conditions outside.
However, as the rarefaction wave that moves backwards from $r_{\mathrm{i}}$
approaches $r_{\mathrm{w}}$, pressure gradients will cause the shell around
the wind bubble to begin to be dispersed into the surrounding medium. Under
these conditions, the recombination rate in the shell drops and more ionising
photons can leak, causing the shell to disintegrate and a strong champagne
flow to begin. This will likely be enhanced if discontinuities in the cloud
density are encountered, e.g. as the shell reaches a sharp step in density
between the molecular cloud and the diffuse interstellar medium around it, as
in Tenorio-Tagle (1979), or if clumping in the cloud breaks up the wind-driven
shell. Under these conditions, the wind bubble will likely rapidly break out
of the cloud in an unstable plume as found in Comeron (1997) and Geen et al.
(2021).

For very massive stars in more diffuse molecular clouds, the wind bubble can
expand much faster than the speed of sound in the photoionised gas, so under
these conditions the wind bubble can expand more rapidly than the gas can
expand thermally after it is photoionised, preventing a strong champagne flow.
The precise behaviour of the Hii region after it reaches the overflow radius
is mostly beyond the scope of this paper, since we primarily wish to explain
why certain Hii regions appear to retain neutral shells that suggest they
precede the emergence of a champagne flow. We thus leave a more detailed
description of the post-overflow dynamics for future work.

##  5 Discussion

In this Section we discuss effects neglected so far in our model, and the
impact they may have on the solutions.

###  5.1 Magnetic Fields

In this paper so far we have neglected magnetic fields. We now discuss their
likely effect on the evolution of the shell around the wind bubble.

We can modify Equation 15 to include a term to describe the magnetic pressure
in the shell. We then get

$P_{\mathrm{w}}=\frac{m_{\mathrm{H}}}{X}n_{\mathrm{i}}(r_{\mathrm{w}})c_{\mathrm{i}}^{2}+\frac{B^{2}}{8\pi},$
(26)

where the first term of the right hand side is the thermal pressure of the
inner edge of the shell $P_{\mathrm{therm}}$ and the second is the magnetic
pressure $P_{\mathrm{mag}}$.

To get an estimate of the likely magnetic field in the shell, we use an
approximation to Crutcher (2012) and assume a critical B-field, where

$\left(\frac{B}{10~{}\mu\mathrm{G}}\right)^{2}=\frac{n_{\mathrm{i}}(r_{\mathrm{w}})}{300~{}\mathrm{cm}^{-3}}$
(27)

in the regime studied here. The magnetic field saturates at
$10~{}\mu\mathrm{G}$ below 300 cm-3. Equation 26 then becomes:

$P_{\mathrm{w}}=\frac{m_{\mathrm{H}}}{X}n_{\mathrm{i}}(r_{\mathrm{w}})c_{\mathrm{i}}^{2}+(10~{}\mu\mathrm{G})^{2}\frac{n_{\mathrm{i}}(r_{\mathrm{w}})}{300~{}\mathrm{cm}^{-3}}.$
(28)

The plasma beta ($\beta\equiv P_{\mathrm{therm}}/P_{\mathrm{mag}}$) of the
ionised shell assuming $c_{\mathrm{i}}\simeq$10 km/s is 169. The pressure from
magnetic fields is thus under 1% of the thermal energy, and so we justify
neglecting it for the rest of the calculations in the paper. The importance of
magnetic fields is likely instead to be to prevent the fragmentation of the
shell and other gaseous structures, as found in Hennebelle (2013) and Geen et
al. (2015). Models taking into account geometry effects on the magnetic field
may also lead to a non-negligible role in shaping the photoionised shell, as
discussed in Pellegrini et al. (2007).

###  5.2 Gravity

Figure 6: Ratio of pressure from gravity on the shell compared to the pressure
from the wind bubble as a function of stellar mass and characteristic
background density $n_{\mathrm{0}}$ at 1 pc where $\omega=2$. Gravity is more
important in denser environments, but typically does not overcome the pressure
inside the wind bubble for this range of densities. A horizontal grey line
shows the point where $P_{\mathrm{grav}}=P_{\mathrm{w}}$.

We have until now ignored gravity in our model. The pressure from gravity on a
spherical shell of inside $r_{\mathrm{i}}$ is given by

$P_{\mathrm{grav}}=\frac{GM(<r_{\mathrm{w}})^{2}}{4\pi r_{\mathrm{i}}^{4}}$
(29)

We calculate $P_{\mathrm{grav}}/P_{\mathrm{w}}$ at the point of overflow for
the same conditions as Figure 4 and plot it for various stellar masses and
densitites in Figure 6. Gravity has only a small effect at low densities, high
metallicities and large stellar masses. Only for the highest densities and
lowest stellar masses does gravity become comparable to the pressure from
winds. We therefore justify ignoring it for this work, although in certain
regimes, particularly if wind cooling becomes efficient and hence
$P_{\mathrm{w}}$ drops, it may act to compress the shell or even stall the
expansion of the Hii region.

For regimes such as older Hii regions around massive clusters, gravity has
been found to be more important where the pressure inside the bubble has begun
to drop due to expansion and radiative cooling (e.g. Rahner et al., 2017).

###  5.3 Radiative Cooling

Figure 7: Ratio between cooling rate of the wind bubble
$\mathrm{d}E_{\mathrm{cool}}/\mathrm{d}t$ and wind luminosity $L_{\mathrm{w}}$
as a function of time, stellar mass and external gas density profile, as
discussed in Section 5.3. The cooling rate predicted here is far lower than
the wind luminosity, particularly at low metallicity, suggesting that the wind
bubble should be largely adiabatic under such conditions, provided increased
mixing with the background does not occur. The effect of increasing the
stellar mass from 30 to 120 M⊙ is roughly the same as decreasing
$n_{\mathrm{0}}$ by a factor of 10.

In this section we discuss the role of radiative cooling on the wind bubble.
To do this, we adopt a similar approach to Mac Low & McCray (1988). The total
energy loss rate from cooling inside the wind bubble can be written as

$\frac{\mathrm{dE_{\mathrm{cool}}}}{\mathrm{d}t}=\int_{\mathrm{0}}^{r_{\mathrm{lim}}}n_{\mathrm{b}}^{2}(r)\Lambda(T_{\mathrm{b}}(r),Z).\Omega
r^{2}dr$ (30)

where the cooling function in the range being considered is given by

$\Lambda(T_{\mathrm{b}}(r,t),Z)=10^{-22}~{}\mathrm{erg~{}cm}^{3}~{}\mathrm{s}^{-1}~{}\left(\frac{T_{\mathrm{b}}(r,t)}{10^{6}~{}\mathrm{K}}\right)^{-7/10}\left(\frac{Z}{Z_{\mathrm{solar}}}\right),$
(31)

and $T_{\mathrm{b}}(r)$ and $n_{\mathrm{b}}(r)$ are the temperature and
hydrogen number density inside the wind bubble as a function of $r$ (distinct
from the external density profile $n(r)$), and $r_{\mathrm{lim}}\rightarrow
r_{\mathrm{w}}$ is the radius inside which the cooling is considered. Weaver
et al. (1977) gives the temperature inside the bubble as

$T_{\mathrm{b}}(r,t)=\left(\frac{P_{\mathrm{w}}r^{2}_{\mathrm{w}}}{tC}\right)^{2/7}(1-r/r_{\mathrm{w}})^{2/5}=T_{\mathrm{c}}(t)(1-r/r_{\mathrm{w}})^{2/5},$
(32)

where $T_{\mathrm{c}}$ is the temperature at the centre of the bubble and
$C=6\times 10^{-7}$ erg s-1 cm-1 K-7/2. Using the ideal gas equation
$P_{\mathrm{w}}=n_{\mathrm{b}}(r)k_{\mathrm{B}}T_{\mathrm{b}}(r)$ together
with Equations 3 and 6, we can find a similar expression for
$n_{\mathrm{b}}(r)$,

$n_{\mathrm{b}}(r,t)=\frac{P_{\mathrm{w}}^{5/7}}{k_{\mathrm{B}}}\left(\frac{tC}{r^{2}_{\mathrm{w}}}\right)^{2/7}(1-r/r_{\mathrm{w}})^{-2/5}=n_{\mathrm{c}}(t)(1-r/r_{\mathrm{w}})^{-2/5},$
(33)

where $n_{\mathrm{c}}$ is the particle density at the centre of the bubble.

Since these equations create a singularity in Equation 30 at $r_{\mathrm{w}}$,
Mac Low & McCray (1988) pick
$r_{\mathrm{lim}}=r_{\mathrm{w}}(1-(T_{\mathrm{cut}}/T_{\mathrm{c}})^{5/2})$
where $T_{\mathrm{cut}}=10^{5}~{}$K, arguing that below this temperature the
cooling function $\Lambda$ changes shape. In addition, if the wind bubble
cooled down to $\sim~{}10^{4}~{}$K, it would reach temperature equilibrium
with the photoionised shell outside.

Since $r_{\mathrm{lim}}$ is a function of $T_{\mathrm{c}}$, we do not give a
fully algebraic solution to Equation 30, but instead plot solutions for it as
a function of time in the case where $\omega=2$ in Figure 7.

We find that for all cases studied, the cooling rate from the wind bubble is
1% of the wind luminosity or lower, dropping over time as the bubble expands
into a lower-density environment, and as such this cooling channel should not
affect these calculations.

There are a few other cooling channels that can become significant. Fierlinger
et al. (2016) invoke thermal conduction across the wind bubble interface, as
do Gentry et al. (2017) and El-Badry et al. (2019) for supernova remnants.
Correct thermal transfer across the contact discontinuity between the wind
bubble and the (photoionised) shell is crucial in determining the energy
budget of the wind bubble. Further, Fierlinger et al. (2016) argue that
cooling in the outer part of the shell as it shock heats against the dense gas
should be correctly treated to account for artificial numerical mixing.
Turbulent mixing can also be important under the conditions described in Tan
et al. (2021). Under these conditions, fractal effects on the structure of the
wind shell can lead to enhanced cooling (Lancaster et al., 2021a, see also
Fielding et al., 2020), causing a momentum-driven wind bubble with negligible
thermal pressure (Silich & Tenorio-Tagle, 2013). Rosen et al. (2014) argue
that cooling is effective in analyses of observed Hii regions, as do Lancaster
et al. (2021b) in simulations of wind bubbles in a turbulent density field.

One other cooling mechanism is described in Arthur (2012), who performed 1D
spherically-symmetric hydrodynamic simulations using a uniform external
density distribution. They find hotter gas in their simulations than is
observed in regions such as the Orion Nebula. They argue that the evaporation
of dense embedded proplyds (protoplanetary disks affected by feedback) causes
mixing with the hot wind bubble, cooling it and slowing its expansion.

We caution that the cooling rate of the wind bubble is very sensitive to the
ability for the hot bubble to mix with the cooler, denser gas outside the
bubble. More careful hydrodynamic simulations of this phenomenon are needed to
constrain the role of cooling further.

###  5.4 Gravitational Shell Instability

Ionising photons can escape if instabilities in the shell cause it to break up
sufficiently for low-density channels to form. Gravitational instabilities
provide one such channel. For a wind shell moving through a uniform medium
($w=0$), Elmegreen (1994) give the following condition for gravitational
instabilities:

$\frac{\pi
G\rho_{\mathrm{0}}}{3c_{\mathrm{s}}}>\frac{8^{1/2}v_{\mathrm{0}}}{r_{\mathrm{w}}^{2}}$
(34)

where $c_{\mathrm{s}}$ is the sound speed in the neutral shell and
$v_{\mathrm{0}}$ is the speed of the shell in the uniform medium. In this
case, the column density of the shell
$\sigma_{\mathrm{0}}=\rho_{\mathrm{0}}r_{\mathrm{w}}/3$. For a $\omega=2$
profile, we can use Equation 2 to give a column density of the shell
$\sigma_{\mathrm{2}}=\rho_{\mathrm{0}}r_{\mathrm{0}}^{2}/r_{\mathrm{w}}$.
Modifying this equation to use the new column density, we arrive at a limiting
velocity for instabilities of

$v_{\mathrm{inst}}=\frac{\pi
G\rho_{\mathrm{0}}r_{\mathrm{0}}^{2}}{8^{1/2}c_{\mathrm{s}}}$ (35)

where $v_{\mathrm{w,2}}<v_{\mathrm{inst}}$ for instability to occur. This
gives $v_{\mathrm{inst}}\simeq(n_{\mathrm{0}}/1000~{}\mathrm{cm}^{-3})~{}$km/s
for $c_{\mathrm{s}}=0.2~{}$km/s. This is much lower than typical values for
$v_{\mathrm{w,2}}$ except for the lower mass stars studied here, as in Figure
3, where it is still at least an order of magnitude lower.

Thus we do not expect the wind bubble to experience gravitational
instabilities before other forms of photon breakout occur. We further discount
Rayleigh-Taylor instabilities since the shell is not accelerating, although
even if it were, the growth rate would typically be small compared to the
timescales studied in Section 4.

However, this does not discount other forms of instability or interactions
with clumps and other aspherical structures in the surrounding medium
triggering some form of breakup. There remains a crucial role for detailed 3D
hydrodynamic simulations to explore such behaviour.

###  5.5 Application to Orion

Recent observations of nebula M42 in Orion, henceforth referred to as the
Orion Nebula, suggest that it is a wind-driven bubble with a neutral shell
travelling at 13 km/s with a derived age of 0.2 Myr (Pabst et al., 2019,
2020). The fact that the shell contains neutral gas implies that the shell
traps ionising radiation from the source star, since otherwise the shell would
be photoionised. In this section we apply our model to the region to determine
whether we predict this behaviour.

The nebula is powered by the O7V star $\theta^{1}$ Ori C. This star has an
evolutionary mass of $33\pm 5$ M⊙ (Balega et al., 2014; Simón-Díaz et al.,
2006). The star is situated close to a dense filament and the nebula is
expanding away from the star out to a distance of 4 pc (Guedel et al., 2008)
in a roughly circular shape on the sky, with the source star near the edge of
the nebula. We assume for the purposes of this work that the wind bubble is
expanding rapidly into a $\omega=2$ density field, since this is similar to
what we find in simulations of young massive stars in clouds similar to this
region (Geen et al., 2021).

We adopt two models for the initial geometry of the region for the purposes of
calculating solutions to the equations in Sections 3 and 4, illustrated with a
schematic in Figure 8. In the first model, the initial gas mass distribution
is centred in the middle of the nebula, with a radius of 2 pc and solid angle
subtending the whole sphere around the centre, $\Omega=4\pi$. This model is
used by authors such as Arthur (2012), who invoke an Hii region expanding into
a uniform medium.

In addition to this model, we adopt a second model. This is because the
position of the star is not at the centre of the nebula, but at one edge
bounded by denser gas that constrains the nebula in that direction (see Figure
8). To better capture this view of the region, the second model places the
star at the corner of a conic section of opening angle $\pi$, reaching out to
4 pc at the edge of the nebula.

Both models adopt a $\omega=2$ density profile. We consider this to be
relatively unlikely in the first model since the remaining neutral dense gas
and star are to one edge of the nebula, although we include it for the
purposes of comparison with uniform background models in the literature that
invoke a spherically symmetric solution with radius 2 pc.

Since the star is rotating only at a few percent of critical velocity (Simón-
Díaz et al., 2006), we use the non-rotating Geneva tracks to model the region
(Ekström et al., 2012), with stellar spectra obtained from Starburst99
(Leitherer et al., 2014) (see Table 3). We give solutions from the 30 and 35
M⊙ tracks to provide estimates of the uncertainty induced by variations in the
stellar evolution model for the star.

Pabst et al. (2020) derive from observations masses for the shell around the
Orion Nebula between 600 and 2600 M⊙ (or up to 3400 M⊙ with the upper bound on
errors given in Pabst et al., 2019), with an expected value of 1500 M⊙. The
uncertainty depends on the method used to calculate the shell mass and which
parts of the shell are considered. In this mass range,
$n_{\mathrm{0}}=1921^{+1200}_{\mathrm{-1200}}~{}$cm-3 if $\Omega=4\pi$ and
$n_{\mathrm{0}}=3842^{+2400}_{\mathrm{-2400}}~{}$cm-3 if $\Omega=\pi$.

From this mass range, we can use Equation 2 to find $\Omega
n_{\mathrm{0}}r_{\mathrm{0}}^{2}$. We use these values as inputs to our semi-
analytic model in Section 4 for the shell expansion. Since the C+ emission
measured by Pabst et al. (2020) comes from the part of the shell after
$r_{\mathrm{i}}$, we use $r_{\mathrm{i}}$ and
$\mathrm{d}r_{\mathrm{i}}/\mathrm{d}t$ as the shell radius and expansion
velocity, and iterate solutions to the model to find $r_{\mathrm{w,2}}$ and
$v_{\mathrm{w,2}}$. However, we find that $r_{\mathrm{i}}$ is roughly equal to
$r_{\mathrm{w,2}}$ to within a few percent in all cases listed, so the results
should be similar to if we used $r_{\mathrm{w,2}}$ and $v_{\mathrm{w,2}}$.

Figure 8: Schematic showing the model used to describe the Orion Nebula (M42). The image on the left shows the location of the principal components of the nebula, including the source of the feedback, star $\theta^{1}$ Ori C in the Trapezium cluster, at the left edge of the nebula. The nearby nebula M43 is not treated in this work but is included for visual reference. The first model, a sphere centred on the nebula but not the star, is shown in the centre panel. The second model, a cone with its point located on the star with flows moving in a radial direction rightwards in a solid angle of $\pi$, is shown in the right panel. We discuss both models in Section 5.5. Geometry | Shell velocity | Age | Overflow radius
---|---|---|---
$r_{\mathrm{i}}$ / pc | $\Omega$ |  $\mathrm{d}r_{\mathrm{i}}/\mathrm{d}t$ / km/s | t / Myr |  $r_{\mathrm{o,2}}$ / pc
Stellar mass = 30 M⊙
2.0 | $4\pi$ | $9.6^{+4.8}_{\mathrm{-1.7}}$ | $0.215^{+0.042}_{\mathrm{-0.067}}$ | $44.4^{+40.5}_{\mathrm{-32.4}}$
4.0 | $1\pi$ | $11.1^{+4.7}_{\mathrm{-1.7}}$ | $0.359^{+0.065}_{\mathrm{-0.103}}$ | $282.2^{+257.0}_{\mathrm{-205.9}}$
Stellar mass = 35 M⊙
2.0 | $4\pi$ | $11.7^{+5.6}_{\mathrm{-2.0}}$ | $0.176^{+0.033}_{\mathrm{-0.053}}$ | $37.7^{+34.3}_{\mathrm{-27.5}}$
4.0 | $1\pi$ | $13.7^{+5.7}_{\mathrm{-2.1}}$ | $0.290^{+0.052}_{\mathrm{-0.083}}$ | $239.2^{+217.8}_{\mathrm{-174.5}}$
Table 1: Table showing the results of our model for the conditions in the
Orion Nebula. We follow two stellar tracks, one at 30 M⊙ and one at 35 M⊙ in
order to capture the uncertainty in the stellar mass measurement for
$\theta^{1}$ Ori C. We further adopt two models for the nebula structure,
visualised in Figure 8. In the first model, the initial gas mass distribution
and the stellar wind and radiation source are centred in the middle of the
nebula, with a radius of 2 pc and solid angle subtending the whole sphere
around the centre, $\Omega=4\pi$. In the second model, the nebula expands in
only $\Omega=\pi$ around the star, the rest of the nebula being constrained by
denser gas, with the star at the point of the conical volume and the nebula
reaching out to 4 pc. We consider the second model to be more representative
of the real Orion nebula, where the source of the winds and radiation is to
one side of the region. For each model setup, the shell velocity
$v_{\mathrm{w,2}}$, age of the nebula $t$ and ionising radiation overflow
radius $r_{\mathrm{o,2}}$ are calculated. The error in the last three
quantities comes from the uncertainty in shell mass used - see Section 5.5 for
more information. The models predict that the overflow radius of the Orion
Nebula is much larger than its current radius, reproducing the observed result
that the ionising radiation should be trapped by the shell around the wind
bubble.

Table 1 lists the model setups used and their outcomes. All of the models give
velocities and ages that agree with the Pabst et al. (2020) results to within
the error given by uncertainties in the shell mass used. Our preferred model,
with the star situated at one end of a Hii region with solid angle
$\Omega=\pi$, predicts a slightly older region of around 0.3 Myr, expanding
around 1-2 km/s faster, than the model where $\Omega=4\pi$ due to the extra
time taken to travel 4 pc instead of 2 pc.

A few additional processes may affect our result. Internal cooling by proplyds
as described in Arthur (2012) may reduce the overall expansion rate of the Hii
region. We also note that there are uncertainties in stellar atmosphere and
evolution models that may affect our result, as well as uncertainties in the
shell masses used that limit a more accurate velocity estimate. Nonetheless,
to within the given observational errors, our model is consistent with the
observed shell velocity.

The other main result of this comparison is to note that our model predicts
that in all cases, the Orion Nebula has an overflow radius well above its
current radius, suggesting that the ionising radiation should be easily
trapped by the shell and cause the Orion Nebula to be a wind-driven structure.
The overflow radius is significantly larger in the case where $\Omega=\pi$
since the wind is compressed into a smaller volume, increasing its relative
pressure and hence reducing the likelihood of overflow.

We note that the lack of 3D information in the model limits our ability to
predict the exact behaviour of the shell when turbulence and instabilities
occur. A better constrained model would require closer comparison to the
observations, perhaps using 3D hydrodynamic simulations to account for the
role of instabilities, turbulence, and other 3D effects in the system. Rather,
the purpose of this work is to establish the feasibility for such behaviour of
wind bubbles in trapping ionising radiation and preventing photoionised
champagne flows as observed in systems like the Orion Nebula.

There are a number of possible reasons why simulations to date do not
reproduce systems behaving in this way, such as the Orion Nebula, in addition
to the cooling arguments made in Section 5.3. One is that they lack the
spatial resolution to properly capture the dense shell, and hence the ionising
photon trapping. Another is the role of magnetic fields, which act to prevent
the breakup of extended structures in molecular clouds (Hennebelle, 2013), and
hence the fragmentation of feedback-driven shells (Geen et al., 2015).
Finally, the precise conditions of feedback in the cloud may simply determine
whether the wind bubble is surrounded by sufficient dense cloud material to
prevent a breakout into the diffuse external medium. Further detailed
comparison is needed to determine whether the discrepancies between
conclusions from observations and numerical theory can be resolved.

##  6 Conclusions

In this paper we study the expansion of wind bubbles in power law density
fields around young massive stars, with a particular focus on “singular
isothermal sphere” density fields $(\rho\propto r^{-2})$ as predicted to exist
around young stars.

The goal of this paper is to a) derive analytic solutions for wind bubble
expansion in non-uniform density fields, b) calculate at what point ionising
radiation can escape such a system and ionise the surrounding medium and c)
determine at what point effects not included in this model grow to require
more detailed numerical models. The reason for undertaking this study is to
reconcile certain observations, which suggest feedback around young massive
stars is driven principally by winds and not ionising radiation, and theory,
which has often argued that winds are not dynamically important.

We find that adiabatic wind bubbles expanding into density fields where
$\rho\propto r^{-2}$ do so at a constant rate $v_{\mathrm{w,2}}$, which
depends on the wind luminosity $L_{\mathrm{w}}$ and characteristic cloud
density $\rho_{\mathrm{0}}$. Radiation pressure contributes only a small
amount to the expansion of the Hii region. In denser cloud environments,
ionising radiation cannot overflow the shell around the wind bubble and escape
into the surrounding medium is prevented until the shell reaches a size of
several parsecs and an age of $\gtrsim 1~{}$Myr. A more likely scenario is
that the conditions in the external medium change to allow shell fragmentation
and ionising photon escape.

Radiative cooling, magnetic fields, gravity and instabilities do not appear to
have a strong influence on the solution, although gravity can become a larger
perturbation around less massive O stars in dense clouds.

Our model provides a framework for understanding the behaviour of the Orion
Nebula. Our preferred model is one in which the 30-35 M⊙ star $\theta^{1}$ Ori
C sits at one edge of the nebula, as is observed, and the nebula expands into
a singular isothermal sphere over a solid angle of $\sim\pi$, constrained by
denser gas in all other directions. Using this model, we recover an expansion
velocity consistent with the observed shell velocity of 13 km/s to within the
errors given by the observational studies. We also explain that the wind-
driven shell should easily trap ionising radiation under the conditions found
in Orion, as observed.

The solution found in this work is likely to break down as the wind bubble
enters the wider structured cloud environment. This highlights the need for
detailed 3D radiative magnetohydrodynamic simulations. At the same time,
sufficient resolution and care with sub-grid feedback injection recipes is
crucial in correctly capturing the interaction between stellar winds,
radiation and the star-forming environment. To reproduce this phenomenon,
simulations should correctly capture the initial wind shell formation before a
champagne flow occurs.

In summary, stellar winds around young massive stars should not be discounted
from a theoretical perspective. Careful interaction between analytic theory,
numerical simulations and observations are needed to determine precisely how
massive stars shape their environment.

## Acknowledgements

The authors would like to thank the anonymous referee for comments that helped
greatly improve the manuscript. The authors would like to thank Xander
Tielens, Cornelia Pabst, and Eric Pellegrini for useful discussions during the
writing of this paper. SG acknowledges support from a NOVA grant for the
theory of massive star formation.

## Data Availability

The code and data products used in this work are stored according to the Anton
Pannekoek Institute for Astronomy Research Data Management plan, and can be
found via DOI 10.5281/zenodo.5646533.

## References

  * Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, The Astrophysical Journal, 770, 25
  * Ali et al. (2018) Ali A., Harries T. J., Douglas T. A., 2018, Monthly Notices of the Royal Astronomical Society, 477, 5422
  * Arthur (2012) Arthur S. J., 2012, Monthly Notices of the Royal Astronomical Society, 421, 1283
  * Avedisova (1972) Avedisova V. S., 1972, Soviet Astronomy, 15
  * Baker & Menzel (1938) Baker J. G., Menzel D. H., 1938, The Astrophysical Journal, 88, 52
  * Balega et al. (2014) Balega Y. Y., Chentsov E., Leushin V., Rzaev A. K., Weigelt G., 2014, Astrophysical Bulletin, 69, 46
  * Bate (2012) Bate M. R., 2012, Monthly Notices of the Royal Astronomical Society, 419, 3115
  * Calderón et al. (2020) Calderón D., Cuadra J., Schartmann M., Burkert A., Prieto J., Russell C. M. P., 2020, Monthly Notices of the Royal Astronomical Society, Advance Access, 493, 447
  * Capriotti & Kozminski (2001) Capriotti E., Kozminski J., 2001, Publications of the Astronomical Society of the Pacific, 113, 677
  * Castor et al. (1975) Castor J., Weaver R., McCray R., 1975, The Astrophysical Journal, 200, L107
  * Comeron (1997) Comeron F., 1997, Astronomy & Astrophysics, 326, 1195
  * Crutcher (2012) Crutcher R. M., 2012, Annual Review of Astronomy and Astrophysics, 50, 29
  * Dale et al. (2012) Dale J. E., Ercolano B., Bonnell I. A., 2012, Monthly Notices of the Royal Astronomical Society, 424, 377
  * Dale et al. (2014) Dale J. E., Ngoumou J., Ercolano B., Bonnell I. A., 2014, Monthly Notices of the Royal Astronomical Society, 442, 694
  * Decataldo et al. (2020) Decataldo D., Lupi A., Ferrara A., Pallottini A., Fumagalli M., 2020, Monthly Notices of the Royal Astronomical Society, 497, 4718
  * Draine (2011) Draine B. T., 2011, The Astrophysical Journal, 732, 100
  * Dunne et al. (2003) Dunne B. C., Chu Y. H., Chen C. H. R., Lowry J. D., Townsley L., Gruendl R. A., Guerrero M. A., Rosado M., 2003, The Astrophysical Journal, 590, 306
  * Dyson & Williams (1980) Dyson J. E., Williams D. A., 1980, Physics of the interstellar medium. New York, Halsted Press
  * Ekström et al. (2012) Ekström S., et al., 2012, Astronomy & Astrophysics, 537, A146
  * El-Badry et al. (2019) El-Badry K., Ostriker E. C., Kim C.-G., Quataert E., Weisz D. R., 2019, Monthly Notices of the Royal Astronomical Society, 490, 1961
  * Elmegreen (1994) Elmegreen B. G., 1994, The Astrophysical Journal, 427, 384
  * Fielding et al. (2020) Fielding D. B., Ostriker E. C., Bryan G. L., Jermyn A. S., 2020, Astrophysical Journal Letters, 894, L24
  * Fierlinger et al. (2016) Fierlinger K. M., Burkert A., Ntormousi E., Fierlinger P., Schartmann M., Ballone A., Krause M. G. H., Diehl R., 2016, Monthly Notices of the Royal Astronomical Society, 456, 710
  * Franco et al. (1990) Franco J., Tenorio-Tagle G., Bodenheimer P., 1990, The Astrophysical Journal, 349, 126
  * Gallegos-Garcia et al. (2020) Gallegos-Garcia M., Burkhart B., Rosen A. L., Naiman J. P., Ramirez-Ruiz E., 2020, Astrophysical Journal Letters, 899, L30
  * Gatto et al. (2017) Gatto A., et al., 2017, Monthly Notices of the Royal Astronomical Society, 466, 1903
  * Geen et al. (2015) Geen S., Hennebelle P., Tremblin P., Rosdahl J., 2015, Monthly Notices of the Royal Astronomical Society, 454, 4484
  * Geen et al. (2020) Geen S., Pellegrini E., Bieri R., Klessen R., 2020, Monthly Notices of the Royal Astronomical Society, 492, 915
  * Geen et al. (2021) Geen S., Bieri R., Rosdahl J., de Koter A., 2021, Monthly Notices of the Royal Astronomical Society, 501, 1352
  * Gentry et al. (2017) Gentry E. S., Krumholz M. R., Dekel A., Madau P., 2017, Monthly Notices of the Royal Astronomical Society, 465, 2471
  * Groenewegen et al. (1989) Groenewegen M., Lamers H., Pauldrach A., 1989, Astronomy & Astrophysics, 221, 78
  * Grudić et al. (2020) Grudić M. Y., Kruijssen J. D., Faucher-Giguère C.-A., Hopkins P. F., Ma X., Quataert E., Boylan-Kolchin M., 2020, arXiv e-prints, p. arXiv:2008.04453
  * Guedel et al. (2008) Guedel M., Briggs K. R., Montmerle T., Audard M., Rebull L., Skinner S. L., 2008, Science, 319, 309
  * Haid et al. (2018) Haid S., Walch S., Seifried D., Wünsch R., Dinnbier F., Naab T., 2018, Monthly Notices of the Royal Astronomical Society, 478, 4799
  * Hennebelle (2013) Hennebelle P., 2013, Astronomy & Astrophysics, 556, A153
  * Howarth et al. (1997) Howarth I. D., Siebert K. W., Hussain G. A., Prinja R. K., 1997, Monthly Notices of the Royal Astronomical Society, 284, 265
  * Kahn (1954) Kahn F. D., 1954, Bulletin of the Astronomical Institutes of the Netherlands, 12, 187
  * Kim et al. (2016) Kim J.-G., Kim W.-T., Ostriker E. C., 2016, The Astrophysical Journal, 819, 23
  * Krumholz & Matzner (2009) Krumholz M. R., Matzner C. D., 2009, The Astrophysical Journal, 703, 1352
  * Kuiper & Hosokawa (2018) Kuiper R., Hosokawa T., 2018, Astronomy & Astrophysics, 616, 22
  * Lamers et al. (1993) Lamers H. J. G. L. M., Leitherer C., Lamers H. J. G. L. M., Leitherer C., 1993, The Astrophysical Journal, 412, 771
  * Lancaster et al. (2021a) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021a, The Astrophysical Journal, 914, 89
  * Lancaster et al. (2021b) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021b, The Astrophysical Journal, 914, 90
  * Lee & Hennebelle (2018) Lee Y.-N., Hennebelle P., 2018, Astronomy & Astrophysics, 611, A88
  * Leitherer et al. (2014) Leitherer C., Ekström S., Meynet G., Schaerer D., Agienko K. B., Levesque E. M., 2014, The Astrophysical Journal Supplement Series, 212, 14
  * Mac Low & McCray (1988) Mac Low M.-M., McCray R., 1988, The Astrophysical Journal, 324, 776
  * Martínez-González et al. (2014) Martínez-González S., Silich S., Tenorio-Tagle G., 2014, Astrophysical Journal, 785, 164
  * Massey et al. (2005) Massey P., Puls J., Pauldrach A., Bresolin F., Kudritzki R. P., Simon T., 2005, Astrophysical Journal, 627, 477
  * Mathews (1967) Mathews W. G., 1967, The Astrophysical Journal, 147, 965
  * Murray et al. (2010) Murray N., Quataert E., Thompson T. A., 2010, The Astrophysical Journal, 709, 191
  * Oort & Spitzer (1955) Oort J. H., Spitzer L., 1955, The Astrophysical Journal, 121, 6
  * Pabst et al. (2019) Pabst C., et al., 2019, Nature, 565, 618
  * Pabst et al. (2020) Pabst C. H. M., et al., 2020, Astronomy & Astrophysics, 639, A2
  * Pellegrini et al. (2007) Pellegrini E. W., et al., 2007, The Astrophysical Journal, 658, 1119
  * Pellegrini et al. (2009) Pellegrini E. W., Baldwin J. A., Ferland G. J., Shaw G., Heathcote S., 2009, The Astrophysical Journal, 693, 285
  * Pellegrini et al. (2020) Pellegrini E. W., Rahner D., Reissl S., Glover S. C. O., Klessen R. S., Rousseau-Nepton L., Herrera-Camus R., 2020, Monthly Notices of the Royal Astronomical Society, 496, pp339
  * Rahner et al. (2017) Rahner D., Pellegrini E. W., Glover S. C. O., Klessen R. S., 2017, Monthly Notices of the Royal Astronomical Society, 470, 4453
  * Rogers & Pittard (2013) Rogers H., Pittard J. M., 2013, Monthly Notices of the Royal Astronomical Society, 431, 1337
  * Rosdahl et al. (2013) Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, Monthly Notices of the Royal Astronomical Society, 436, 2188
  * Rosen et al. (2014) Rosen A. L., Lopez L. A., Krumholz M. R., Ramirez-Ruiz E., 2014, Monthly Notices of the Royal Astronomical Society, 442, 2701
  * Shu et al. (2002) Shu F., Lizano S., Galli D., Canto J., Laughlin G., 2002, The Astrophysical Journal, 580, 969
  * Silich & Tenorio-Tagle (2013) Silich S., Tenorio-Tagle G., 2013, The Astrophysical Journal, 765
  * Simón-Díaz et al. (2006) Simón-Díaz S., Herrero A., Esteban C., Najarro F., 2006, Astronomy & Astrophysics, 448, 351
  * Spitzer (1978) Spitzer L., 1978, Physical processes in the interstellar medium. New York Wiley-Interscience
  * Tan et al. (2021) Tan B., Oh S. P., Gronke M., 2021, Monthly Notices of the Royal Astronomical Society, 502, 3179
  * Tenorio-Tagle (1979) Tenorio-Tagle G., 1979, Astronomy and Astrophysics, 71, 59
  * Verdolini et al. (2013) Verdolini S., et al., 2013, The Astrophysical Journal, 769, 12
  * Vink et al. (2011) Vink J. S., Muijres L. E., Anthonisse B., de Koter A., Graefener G., Langer N., 2011, Astronomy & Astrophysics, 531, 132
  * Wall et al. (2019) Wall J. E., McMillan S. L. W., Mac Low M.-M., Klessen R. S., Zwart S. P., 2019, Astrophysical Journal, 887
  * Wall et al. (2020) Wall J. E., Mac Low M.-M., McMillan S. L., Klessen R. S., Portegies Zwart S., Pellegrino A., 2020, Astrophysical Journal, 904, 192
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, The Astrophysical Journal, 218, 377
  * Yeh & Matzner (2012) Yeh S. C. C., Matzner C. D., 2012, The Astrophysical Journal, 757, 18
  * Zamora-Avilés et al. (2019) Zamora-Avilés M., et al., 2019, Monthly Notices of the Royal Astronomical Society, 487, 2200

##  Appendix A Calculation of the Energy Partition in the Wind Bubble

In this Appendix we derive the fraction of energy from the wind that remains
in the wind bubble, as opposed to being used to do work on the dense shell.
This is a generalised re-derivation of Section II of Weaver et al. (1977).

We begin with a dimensional analysis which establishes the power law
relationships between the radius of the wind bubble $r_{\mathrm{w}}$, the wind
luminosity $L_{\mathrm{w}}$, the background density at $r_{\mathrm{w}}$,
$\rho(r_{\mathrm{w}})=\rho_{\mathrm{0}}(r_{\mathrm{w}}/r_{\mathrm{0}})^{-\omega}$
and time $t$:

$r_{\mathrm{w}}^{5}(t)\rho(r_{\mathrm{w}})L_{\mathrm{w}}^{-1}t^{-3}\propto 1$
(36)

and thus

$r_{\mathrm{w}}(t)\propto\left(L_{\mathrm{w}}t^{3}r_{\mathrm{0}}^{-\omega}\rho_{\mathrm{0}}^{-1}\right)^{1/(5-\omega)}.$
(37)

We further assume a constant partition of the energy added by the wind to the
bubble

$E_{\mathrm{b}}=f_{\mathrm{b}}L_{\mathrm{w}}t$ (38)

where $f_{\mathrm{b}}$ is a fraction to be derived.

Weaver et al. (1977) divide the wind bubble into three phases. The first, from
$0\leq r\leq R_{\mathrm{1}}$, is the free-streaming stellar wind. At
$R_{\mathrm{1}}$, the wind shocks against the surrounding medium, and from
$R_{\mathrm{1}}\leq r\leq R_{\mathrm{c}}$ is the shocked stellar wind.
Finally, from $R_{\mathrm{c}}\leq r\leq R_{\mathrm{2}}$ is the shocked
interstellar material, which becomes a dense shell once radiative cooling
becomes efficient and the internal pressure reduces. Due to the lower density
in the region $R_{\mathrm{1}}\leq r\leq R_{\mathrm{c}}$, this volume remains
roughly adiabatic until later times.

We adopt the hydrostatic approximation as in Weaver et al. (1977) for the
region $R_{\mathrm{1}}\leq r\leq R_{\mathrm{c}}$ due to the typically high
sound speeds and short sound-crossing times. Thus the internal pressure
$P(r,t)=P(t)=P(R_{\mathrm{c}})$.

The continuity equation for mass is given by

$\frac{D\rho}{Dt}+\frac{\rho}{r^{2}}\frac{\partial(r^{2}v)}{\partial r}=0.$
(39)

This is the standard equation for mass conservation in spherical coordinates,
assuming all flows occur in radial directions, where the material derivative
$D/Dt$ is defined as

$\frac{D}{Dt}\equiv\frac{\partial}{\partial t}+v\frac{\partial}{\partial r}.$
(40)

We further invoke the energy conservation equation for the adiabatic flow

$\frac{D(P\rho^{-\gamma})}{Dt}=0$ (41)

where $\gamma=5/3$, i.e. the gas is monatomic.

Expanding Equation 41 gives us

$\frac{D\rho}{Dt}=\frac{\rho}{P\gamma}\frac{DP}{Dt},$ (42)

and hence we can use Equation 39 to write

$\frac{1}{r^{2}}\frac{\partial(r^{2}v)}{\partial
r}=\frac{1}{P\gamma}\frac{DP}{Dt}.$ (43)

To solve the right hand side of this equation, we note that the pressure in
the bubble is the energy divided by the volume, and thus $P\propto
L_{\mathrm{w}}t/r^{3}$. Using Equation 37, we can write

$P\propto t^{-(4+\omega)/(5-\omega)}$ (44)

and so

$\frac{1}{P\gamma}\frac{DP}{Dt}=\frac{1}{P\gamma}\frac{\partial P}{\partial
t}=-\frac{3}{5}.\frac{4+\omega}{5-\omega}\frac{1}{t}.$ (45)

Substituting this into Equation 43 and integrating, we can write

$v(r,t)=\frac{f(t)}{r^{2}}+\frac{4+\omega}{25-5\omega}\frac{r}{t}.$ (46)

By differentiating Equation 37 in $t$, we have
$v(R_{\mathrm{c}},t)={3R_{\mathrm{c}}/(5-\omega)t}$ at $R_{\mathrm{c}}$, and
so we can solve for $f(t)$ to give111Note that Equation 11 of Weaver et al.
(1977) contains a typo, and should read $R_{\mathrm{c}}^{3}$, not
$R_{\mathrm{c}}^{2}$

$v(r,t)=\frac{11-\omega}{25-5\omega}\frac{R_{\mathrm{c}}^{3}}{r^{2}t}+\frac{4+\omega}{25-5\omega}\frac{r}{t}.$
(47)

We adopt the solutions for an adiabatic wind in Weaver et al. (1977) in the
region $R_{\mathrm{1}}\ll r\ll R_{\mathrm{c}}$

$v\sim\frac{v_{\mathrm{w}}}{4}\left(\frac{15}{16}\right)^{3/2}\frac{R_{\mathrm{1}}^{2}}{r^{2}}$
(48)

and

$P\sim\frac{3}{4}\left(\frac{16}{15}\right)^{5/2}\frac{v_{\mathrm{w}}\dot{M}_{\mathrm{w}}}{\Omega
R_{\mathrm{1}}^{2}}$ (49)

where $v_{\mathrm{w}}$ is the wind terminal velocity, and
$\dot{M}_{\mathrm{w}}$ is the wind mass loss rate. The wind luminosity is
given by
$L_{\mathrm{w}}\equiv\frac{1}{2}\dot{M}_{\mathrm{w}}v_{\mathrm{w}}^{2}$. We
also adopt a solid angle for a spherical segment $\Omega$, which equals $4\pi$
for a full spherical wind bubble as in Weaver et al. (1977).

Noting the condition that $r\ll R_{\mathrm{c}}$ as in Weaver et al. (1977), we
neglect the second term in Equation 47 and combine this with Equation 48 to
give

$R_{\mathrm{1}}=\left(\frac{44-4\omega}{25-5\omega}\right)^{1/2}\left(\frac{16}{15}\right)^{3/4}\left(\frac{R_{\mathrm{c}}^{3}}{v_{\mathrm{w}}t}\right)^{1/2}$
(50)

Neglecting the kinetic energy in the wind bubble (since the solution is
hydrostatic), the energy in the region $R_{\mathrm{1}}<r<R_{\mathrm{c}}$ is
given by

$E_{\mathrm{b}}=\frac{3}{2}\frac{\Omega}{3}\left(R_{\mathrm{c}}^{3}-R_{\mathrm{1}}^{3}\right)P(t)\simeq\frac{\Omega}{2}R_{\mathrm{c}}^{3}P(t)$
(51)

Substituting Equations 49 and 50 into this, we arrive at

$E_{\mathrm{b}}=\left(\frac{5-\omega}{11-\omega}\right)L_{\mathrm{w}}t$ (52)

which reduces to the Weaver et al. (1977) result for the case $\omega=0$.

##  Appendix B Comparison with Numerical Simulations

In Appendix A we re-derive the fraction of energy from the wind stored in the
hot wind bubble, versus the amount used to do work on the dense shell around
it, as in Weaver et al. (1977) but for a power law density field. In this
Section we verify the calculation using a numerical simulation.

The simulations were performed using
vh-1222http://wonka.physics.ncsu.edu/pub/VH-1/index.php, a freely available
hydrodynamics code. We use a 1D spherically symmetric grid with radius 4 pc,
subdivided into 2048 cells, numbered 0 to 2047, with identical radial
thickness $\Delta x=403~{}$AU, and employ reflective boundary conditions. The
initial gas state is a power law density field as in Equation 1, with
$n_{\mathrm{0}}=1000~{}$cm-3, $r_{\mathrm{0}}=1~{}$pc, and three simulations
performed with $\omega=0,1,2$ respectively. We set the density in cell 0 to
the density in cell 1 to prevent a singularity at $r=0$. The temperature of
the gas is set initially to 10 K everywhere, and $\gamma=5/3$. No gravity,
radiative transfer or radiative cooling is included.

Throughout the simulation, we inject a constant source of mass
$\dot{M}_{\mathrm{w}}$, momentum $\dot{M}_{\mathrm{w}}v_{\mathrm{w}}$ and
kinetic energy
$L_{\mathrm{w}}=\frac{1}{2}\dot{M}_{\mathrm{w}}v_{\mathrm{w}}^{2}$ into cell 0
directed outwards. We use $L_{\mathrm{w}}=10^{36}~{}$erg/s and
$\dot{M}_{\mathrm{w}}=2.5\times 10^{19}~{}$g/s, equivalent to a $\sim 30~{}$M⊙
star. In addition to the inbuilt Courant condition in vh-1, we restrict the
simulation timestep to be smaller than $\Delta x/(100~{}$km/s), in order to
ensure that the first few timesteps before the wind bubble is established are
not too long. We end the simulations at 200 kyr, before the edge of the swept-
up material around the wind bubble encounters the edge of the simulation
volume in any simulation.

The energy in the bubble $E_{\mathrm{b}}$ is defined to be the total energy in
all cells above $10^{5}~{}$K, a threshold we select to distinguish between the
hot wind bubble and the denser shell around it. We plot this quantity as a
fraction of $L_{\mathrm{w}}t$ against time in Figure 9 for each simulation
where $\omega=0,1,2$, and overplot the value of Equation 52 in each case. Once
the wind bubble is established, the results converge towards the analytic
solution. The results converge faster where $\omega$ is smaller. We also note
that the solution in Equation 52, as with Weaver et al. (1977), makes various
approximations and is not an exact solution. A full analysis of the
convergence to the analytic solution, including the role of sound waves (i.e.
non-hydrostatic effects), is beyond the scope of this paper.

We additionally plot the radial expansion of the wind bubble in these
simulations against the analytic solution given in Equation 7, with $\omega=0$
reducing to the canonical Weaver et al. (1977) solution. For the simulations,
we calculate the radius of the wind bubble as being the largest radius of any
cell where the temperature is above $10^{5}~{}$K. There is good agreement in
all cases. The radii converge around 100 kyr, although the expansion rates are
still somewhat different.

Figure 9: Comparison between Equation 52 and 1D numerical hydrodynamic
simulations for the fraction of energy from stellar winds retained in the wind
bubble, as opposed to the energy used to do work on the shell around the wind
bubble. We show three simulations using initial power law density fields whose
index is $\omega=0,1,2$ as solid lines. Details of the simulations are given
in Section B. Dashed lines show the value of Equation 52 in each case. This
figure shows that the energy fraction in the wind bubble converges to the
analytic solution. Figure 10: Comparison between Equation 7 and 1D numerical
hydrodynamic simulations for the radius of the wind bubble, using the same
simulations and line styles as in Figure 9.

##  Appendix C Stellar Properties

log$(M_{\mathrm{ini}}$ / M${}_{\mathrm{\odot}})$ | log$(L_{\mathrm{w}}$ / erg/s$)$ | log$(\dot{M}_{\mathrm{w}}$ / M⊙ / yr$)$ | log$(L_{\mathrm{n}}$ / erg/s$)$ | log$(L_{\mathrm{i}}$ / erg/s$)$ | log$(Q_{\mathrm{H}}/$s${}^{-1})$ | log$(T_{\mathrm{i}}/$K$)$
---|---|---|---|---|---|---
120 | 37.51 | (36.8) | -5.01 | (-5.01) | 39.41 | (39.17) | 39.57 | (39.69) | 50.07 | (50.1) | 3.885 | (4.184)
115 | 37.47 | (36.75) | -5.058 | (-5.058) | 39.38 | (39.15) | 39.55 | (39.67) | 50.04 | (50.07) | 3.884 | (4.182)
110 | 37.42 | (36.71) | -5.107 | (-5.107) | 39.35 | (39.12) | 39.52 | (39.64) | 50.01 | (50.04) | 3.882 | (4.181)
105 | 37.38 | (36.67) | -5.156 | (-5.156) | 39.32 | (39.08) | 39.48 | (39.6) | 49.98 | (50.01) | 3.88 | (4.179)
100 | 37.33 | (36.62) | -5.205 | (-5.205) | 39.29 | (39.05) | 39.45 | (39.57) | 49.94 | (49.98) | 3.879 | (4.178)
95 | 37.28 | (36.57) | -5.253 | (-5.253) | 39.25 | (39.01) | 39.41 | (39.53) | 49.9 | (49.94) | 3.877 | (4.176)
90 | 37.23 | (36.52) | -5.301 | (-5.301) | 39.2 | (38.97) | 39.37 | (39.49) | 49.86 | (49.89) | 3.876 | (4.175)
85 | 37.17 | (36.46) | -5.35 | (-5.35) | 39.16 | (38.92) | 39.32 | (39.44) | 49.82 | (49.85) | 3.874 | (4.173)
80 | 37.1 | (36.39) | -5.428 | (-5.428) | 39.11 | (39.1) | 39.27 | (39.29) | 49.77 | (49.74) | 3.871 | (4.17)
75 | 37.02 | (36.32) | -5.507 | (-5.507) | 39.06 | (39.05) | 39.22 | (39.24) | 49.71 | (49.68) | 3.867 | (4.167)
70 | 36.94 | (36.24) | -5.586 | (-5.586) | 39.01 | (39.0) | 39.18 | (39.2) | 49.67 | (49.64) | 3.864 | (4.164)
65 | 36.85 | (36.16) | -5.664 | (-5.664) | 38.95 | (38.94) | 39.11 | (39.13) | 49.61 | (49.58) | 3.86 | (4.161)
60 | 36.77 | (36.07) | -5.743 | (-5.743) | 38.94 | (38.87) | 38.99 | (39.06) | 49.51 | (49.51) | 3.857 | (4.158)
55 | 36.63 | (35.94) | -5.882 | (-5.882) | 38.86 | (38.79) | 38.91 | (38.98) | 49.43 | (49.43) | 3.849 | (4.153)
50 | 36.48 | (35.81) | -6.021 | (-6.021) | 38.78 | (38.7) | 38.83 | (38.9) | 49.34 | (49.34) | 3.842 | (4.147)
45 | 36.33 | (35.66) | -6.161 | (-6.161) | 38.7 | (38.62) | 38.75 | (38.82) | 49.26 | (49.26) | 3.834 | (4.142)
40 | 36.18 | (35.51) | -6.301 | (-6.301) | 38.58 | (38.51) | 38.63 | (38.7) | 49.14 | (49.15) | 3.827 | (4.136)
35 | 35.94 | (35.29) | -6.525 | (-6.525) | 38.55 | (38.47) | 38.4 | (38.51) | 48.94 | (48.97) | 3.818 | (4.127)
30 | 35.66 | (35.02) | -6.793 | (-6.793) | 38.39 | (38.38) | 38.24 | (38.29) | 48.78 | (48.76) | 3.808 | (4.117)
25 | 35.3 | (34.68) | -7.127 | (-7.127) | 38.2 | (38.19) | 38.05 | (38.1) | 48.59 | (48.57) | 3.799 | (4.108)
20 | 34.79 | (34.2) | -7.598 | (-7.598) | 38.03 | (38.07) | 37.65 | (37.63) | 48.22 | (48.14) | 3.792 | (4.077)
Table 2: Stellar properties used in this paper, for rotating stars (see Section 2.2). Figures outside brackets are for stars at Solar metallicity ($Z=0.014$), with figures in brackets for stars at sub-Solar metallicity ($Z=0.02$). $M_{\mathrm{ini}}$ is the initial mass of the star, $L_{\mathrm{w}}$ is luminosity (i.e. mechanical power injection $\equiv 1/2\dot{M}_{\mathrm{w}}V_{\mathrm{w}}^{2}$) of the wind, $\dot{M}_{\mathrm{w}}$ is the mass loss rate, $L_{\mathrm{n}}$ is the non-ionising radiative luminosity of the star (photons below 13.6 eV), $L_{\mathrm{i}}$ is the total ionising luminosity, $Q_{\mathrm{H}}$ is the hydrogen-ionising photon emission rate and $T_{\mathrm{i}}$ is the temperature of the photionised gas around the star (assumed to be density-independent to good approximation). log$(M_{\mathrm{ini}}$ / M${}_{\mathrm{\odot}})$ | log$(L_{\mathrm{w}}$ / erg/s$)$ | log$(\dot{M}_{\mathrm{w}}$ / M⊙ / yr$)$ | log$(L_{\mathrm{n}}$ / erg/s$)$ | log$(L_{\mathrm{i}}$ / erg/s$)$ | log$(Q_{\mathrm{H}}/$s${}^{-1})$ | log$(T_{\mathrm{i}}/$K$)$
---|---|---|---|---|---|---
120 | 37.5 | (36.77) | -5.035 | (-5.035) | 39.42 | (39.17) | 39.58 | (39.69) | 50.07 | (50.1) | 3.893 | (4.191)
115 | 37.46 | (36.73) | -5.083 | (-5.083) | 39.39 | (39.15) | 39.55 | (39.67) | 50.05 | (50.07) | 3.892 | (4.189)
110 | 37.41 | (36.69) | -5.13 | (-5.13) | 39.36 | (39.12) | 39.52 | (39.64) | 50.02 | (50.04) | 3.89 | (4.187)
105 | 37.37 | (36.65) | -5.177 | (-5.177) | 39.33 | (39.08) | 39.49 | (39.6) | 49.99 | (50.01) | 3.888 | (4.186)
100 | 37.32 | (36.6) | -5.224 | (-5.224) | 39.29 | (39.05) | 39.46 | (39.57) | 49.95 | (49.98) | 3.887 | (4.184)
95 | 37.27 | (36.55) | -5.271 | (-5.271) | 39.26 | (39.01) | 39.42 | (39.53) | 49.91 | (49.94) | 3.885 | (4.183)
90 | 37.22 | (36.51) | -5.318 | (-5.318) | 39.21 | (38.97) | 39.37 | (39.49) | 49.87 | (49.89) | 3.883 | (4.181)
85 | 37.17 | (36.45) | -5.366 | (-5.366) | 39.17 | (38.92) | 39.33 | (39.44) | 49.82 | (49.85) | 3.882 | (4.179)
80 | 37.09 | (36.38) | -5.443 | (-5.443) | 39.12 | (39.1) | 39.28 | (39.29) | 49.78 | (49.74) | 3.878 | (4.176)
75 | 37.02 | (36.31) | -5.519 | (-5.519) | 39.07 | (39.05) | 39.23 | (39.24) | 49.72 | (49.68) | 3.875 | (4.173)
70 | 36.94 | (36.24) | -5.596 | (-5.596) | 39.03 | (39.0) | 39.19 | (39.2) | 49.68 | (49.64) | 3.871 | (4.17)
65 | 36.86 | (36.16) | -5.673 | (-5.673) | 38.96 | (38.94) | 39.12 | (39.13) | 49.62 | (49.58) | 3.867 | (4.167)
60 | 36.77 | (36.07) | -5.75 | (-5.75) | 38.89 | (38.87) | 39.05 | (39.06) | 49.55 | (49.51) | 3.864 | (4.164)
55 | 36.64 | (35.95) | -5.886 | (-5.886) | 38.88 | (38.79) | 38.92 | (38.98) | 49.44 | (49.43) | 3.856 | (4.158)
50 | 36.5 | (35.81) | -6.022 | (-6.022) | 38.79 | (38.7) | 38.84 | (38.9) | 49.35 | (49.34) | 3.849 | (4.153)
45 | 36.35 | (35.67) | -6.157 | (-6.157) | 38.71 | (38.62) | 38.76 | (38.82) | 49.27 | (49.26) | 3.842 | (4.147)
40 | 36.2 | (35.53) | -6.293 | (-6.293) | 38.59 | (38.51) | 38.64 | (38.7) | 49.16 | (49.15) | 3.835 | (4.142)
35 | 35.96 | (35.3) | -6.517 | (-6.517) | 38.47 | (38.47) | 38.52 | (38.51) | 49.03 | (48.97) | 3.824 | (4.133)
30 | 35.68 | (35.04) | -6.783 | (-6.783) | 38.4 | (38.38) | 38.25 | (38.29) | 48.79 | (48.76) | 3.813 | (4.123)
25 | 35.33 | (34.71) | -7.111 | (-7.111) | 38.21 | (38.19) | 38.06 | (38.1) | 48.6 | (48.57) | 3.802 | (4.112)
20 | 34.82 | (34.24) | -7.593 | (-7.593) | 38.1 | (38.07) | 37.49 | (37.63) | 48.03 | (48.14) | 3.797 | (4.088)
Table 3: As in Table 2 but for non-rotating stars.

We give the stellar properties used in this paper in Table 2 for rotating
stars, and 3 for non-rotating stars. For simplicity, in this paper the values
are taken to be constant throughout the star’s lifetime, using the values at a
stellar age of $10^{5}~{}$years. Stellar properties are taken from the Geneva
stellar evolution models (Ekström et al., 2012), with atmospheres calculated
using Starburst99 (Leitherer et al., 2014). Wind luminosities are calculated
by converting the surface escape velocity of the star into terminal velocity
$V_{\mathrm{w}}$ by using the prescription described in Gatto et al. (2017).
Further detailed descriptions of the techniques used to calculate the stellar
models and plots of the evolutions of these properties at late times can be
found in Geen et al. (2021). Based on Figure A1 of that paper, we justify our
choice to use constant wind luminosities and photon emission rates by noting
that these values do not change significantly for the first 1-2 Myr of the
stars’ lives.