# Efficiently Cooled Stellar Wind Bubbles in Turbulent Clouds II.
Validation of Theory with Hydrodynamic Simulations

Lachlan Lancaster Department of Astrophysical Sciences, Princeton University,
4 Ivy Lane, Princeton, NJ 08544, USA  Eve C. Ostriker Department of
Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544,
USA  Jeong-Gyu Kim Department of Astrophysical Sciences, Princeton University,
4 Ivy Lane, Princeton, NJ 08544, USA  Chang-Goo Kim Department of
Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544,
USA

(Received XXX; Revised XXX; Accepted XXX)

###### Abstract

In a companion paper, we develop a theory for the evolution of stellar wind
driven bubbles in dense, turbulent clouds. This theory proposes that turbulent
mixing at a fractal bubble-shell interface leads to highly efficient cooling,
in which the vast majority of the input wind energy is radiated away. This
energy loss renders the majority of the bubble evolution momentum-driven
rather than energy-driven, with expansion velocities and pressures orders of
magnitude lower than in the classical Weaver et al. (1977) solution. In this
paper, we validate our theory with three-dimensional, hydrodynamic
simulations. We show that extreme cooling is not only possible, but is generic
to star formation in turbulent clouds over more than three orders of magnitude
in density. We quantify the few free parameters in our theory, and show that
the momentum exceeds the wind input rate by only a factor $\alpha_{p}\sim
1.2-4$. We verify that the bubble/cloud interface is a fractal with dimension
$\sim 2.5-2.7$. The measured turbulent amplitude ($v_{t}\sim
200-400\mathrm{\,km\,s^{-1}}$) in the hot gas near the interface is shown to
be consistent with theoretical requirements for turbulent diffusion to
efficiently mix and radiate away most of the wind energy. The fraction of
energy remaining after cooling is only $1-\Theta\sim 0.1-0.01$, decreasing
with time, explaining observations that indicate low hot-gas content and weak
dynamical effects of stellar winds.

ISM, Stellar Winds, Star forming regions

††journal: ApJ††software:  Athena (Stone et al., 2008; Stone & Gardiner,
2009), astropy (Astropy Collaboration et al., 2013, 2018), scipy (Virtanen et
al., 2020), numpy (Harris et al., 2020), IPython (Perez & Granger, 2007),
matplotlib (Hunter, 2007), xarray (Hoyer et al., 2017), pandas (Reback et al.,
2020), adstex (https://github.com/yymao/adstex)

##  1 Introduction

Feedback from massive stars is thought to be the dominant process setting the
lifetime efficiency of star formation on the scale of individual molecular
clouds (McKee & Ostriker, 2007; Lada & Lada, 2003). Luminous hot stars use
their strong radiation to disperse the surrounding dense gas in several ways:
imparting primarily-outward photon momentum to gas and dust where stellar
UV/optical radiation is first absorbed (direct radiation pressure); imparting
preferentially-outward photon momentum to dust where diffuse infrared is
absorbed (reprocessed radiation pressure); ionizing and heating gas, which
leads to over-pressured expansion of ionized gas and creates a rocket effect
on neutral structures where photoevaporation occurs (pressure of photo-ionized
gas); or directly depositing momentum in the star’s own atmosphere, resulting
in high-velocity stellar winds which shock and transfer momentum to the
ambient gas. It is much debated in the literature which of these mechanisms
dominates in which situations within the star-forming interstellar medium
(ISM) (Krumholz et al., 2019). While several recent efforts have used
numerical radiation hydrodynamic simulations for in-depth studies of the
effects of radiation feedback on star-forming clouds (Dale et al., 2012, 2013;
Walch et al., 2012; Raskutti et al., 2016, 2017; Howard et al., 2017; Kim et
al., 2018, 2019, 2020; Grudić et al., 2020; Geen et al., 2020b, a; Fukushima
et al., 2020), there have been fewer detailed numerical studies of stellar
wind feedback (reviewed below). In this paper, we focus on the thermal and
dynamical effects of stellar winds driven by star clusters in their natal
molecular clouds.

In an accompanying paper (Lancaster et al., 2021, hereafter, Paper I) we have
theoretically investigated the dynamics of expanding stellar wind bubbles in
the case that most of the wind’s energy is radiated away. There, we argue that
the interface between hot bubbles and the surrounding gas is subject to strong
turbulent mixing, which subsequently leads to efficient radiative cooling in
intermediate-temperature gas. We also argue that the surface of the hot bubble
is a fractal, and the large associated area enhances cooling. We hypothesize
that efficient cooling leads to the dominant phase of bubble evolution being
momentum-driven rather than energy-driven. In this work, we support the
theoretical model of Paper I with a suite of three-dimensional (3D)
hydrodynamic simulations, which we analyze to provide quantitative
corroboration of the arguments made there. Each simulation follows the
expansion of a hot bubble driven into the surrounding dense, turbulent ISM by
a constant luminosity point source.

Several numerical studies have previously explored the dynamical effects of
main sequence stellar winds on their environment. These have included one-
dimensional (1D) numerical models (Garcia-Segura et al., 1996a, b; Silich &
Tenorio-Tagle, 2013; Krause et al., 2016; Fierlinger et al., 2016; Rahner et
al., 2017) some of which have incorporated phenomenological accounts of three-
dimensional processes such as energy leakage (Harper-Clark & Murray, 2009) and
turbulent diffusivity that leads to cooling (El-Badry et al., 2019). These
problems have also been explored in two-dimensional (2D) studies, which have
sought to examine filamentation at reduced computational expense (Wünsch et
al., 2008; Ntormousi et al., 2011; Dwarkadas & Rosenberg, 2013). Finally,
numerical investigations of stellar wind feedback effects in 3D have ranged
from targeted studies of winds (Rogers & Pittard, 2013; Wareing et al., 2017;
Krause et al., 2013; Krause & Diehl, 2014; Dale & Bonnell, 2008; Dale et al.,
2013), to comparisons of winds and photo-ionization heating over a range of
parameter space (Dale et al., 2014; Geen et al., 2015; Haid et al., 2018; Geen
et al., 2020a). These simulations have contributed to assessments of the
relative contribution of winds in a number of scenarios and over a range of
densities and cloud masses. Some have additionally tried to constrain the
relative contribution of leakage and turbulent mixing.

In this paper, we explore a set of questions related to the detailed structure
and evolution of bubbles driven by stellar winds in turbulent clouds. Our
goals are to characterize the dependence of bubble evolution on input wind
power and ambient cloud properties, and to explain the physical mechanisms
controlling the evolution. Using numerical simulations, we quantify (1) the
temporal evolution of wind bubble sizes and the energy and momentum that
bubbles/shells contain; (2) the properties of the turbulence that drives
mixing at the interface between a bubble’s hot/diffuse interior and cool/dense
shell; (3) the fractal structure and foldedness of the cooling interface; and
(4) the total energy losses to cooling due to turbulent mixing. Confirming the
theory laid out in Paper I, we shall show that for the conditions within star-
forming clouds, the vast majority of wind energy is expected to be lost due to
turbulent mixing followed by radiative cooling, and that the large area
associated with the fractal geometry of the bubble/cloud interface is crucial
to enhancing losses. Our physical picture of turbulent, cooling mixing layers
is much informed by insights from recent numerical and analytic investigations
of the Kelvin-Helmholtz instability in the context of multiphase galactic
winds and the circumgalactic medium (Gronke & Oh, 2018; Fielding et al., 2020;
Tan et al., 2020).

The numerical model we adopt is intentionally idealized. It is designed to
provide a testing ground for the theory presented in Paper I that represents
key real-world features (to the extent that is computationally practical),
while still having that theory be applicable. To that end, we (i) have ignored
the effects of star formation that is extended in both space and time, (ii)
approximated the winds as having constant mechanical luminosity, (iii) ignored
the effects of magnetic fields and stellar radiation, and (iv) adopted the
assumption of collisional and ionization equilibrium for computing cooling
losses in ionized gas, and a simple analytic fit for cooling in warm/cold
neutral gas. Thus, our simulations are designed to test our theoretical
predictions, rather than faithfully represent the real world. By tackling an
idealized numerical problem, however, we are able to validate our theoretical
framework and quantitative predictions that can be extrapolated beyond the
necessarily limited scope of any single simulation suite.

The structure of this paper is as follows. In Section 2 we briefly review our
theory for evolution of wind-driven bubbles. The reader is referred to Paper I
for a complete description. In Section 3 we describe our numerical methods for
implementing winds in the Athena code, as well as our simulation set-up and
the parameters of our model suite. We present the analysis of our results in
Section 4 and conclude with a summary of our findings in Section 5.

##  2 Review of Paper I Theory: Wind Bubble Dynamics with Efficient Cooling

Our theory is developed in full in Paper I; for convenience we briefly
summarize key features here. Our global dynamical evolution model is based on
the premise that energy losses in the bubble/shell interface are so great that
the bubble is effectively momentum-driven rather than effectively energy-
driven111 We use the term “momentum-driven” (“energy-driven”) to denote a
solution in which momentum (energy) increases linearly in time. The Castor et
al. (1975)/Weaver et al. (1977) solution is also internally energy-conserving
in the sense that there are no radiative energy losses from the hot bubble or
its interface with the shell; however, the leading shock is assumed to be
fully radiative so that a fraction $27/77=0.35$ of the input wind energy is
radiated away as the shocked and accelerated ambient gas cools to join the
exterior of the shell.. We term this global solution the “Efficiently Cooled
(EC)” stellar wind-driven bubble.

In the EC model, there is a point source of constant mechanical luminosity
$\mathcal{L}_{w}$ and mass loss rate $\dot{M}_{w}$, which injects into the
surrounding cloud a wind with velocity

${\cal{V}}_{w}\equiv\left({\frac{2\mathcal{L}_{w}}{\dot{M}_{w}}}\right)^{1/2}\,$
(1)

and momentum input rate

$\dot{p}_{w}\equiv{\cal{V}}_{w}\dot{M}_{w}=2\mathcal{L}_{w}/{\cal{V}}_{w}=\left({2\mathcal{L}_{w}\dot{M}_{w}}\right)^{1/2}\,.$
(2)

The momentum of the system is mostly contained in the bubble shell, and
increases linearly in time as

$p_{r}=\alpha_{p}p_{\rm
EC}\equiv\alpha_{p}\dot{p}_{w}t=2\alpha_{p}\frac{\mathcal{L}_{w}}{{\cal{V}}_{w}}t\,.$
(3)

This allows for amplification by a factor $\alpha_{p}$ (assumed order-unity in
the EC theory) above the value $p_{\rm EC}$ that corresponds to the momentum
originally injected in the wind, which would apply if there is no thermal
energy buildup within the bubble (maximally efficient mixing/cooling).

We define $\mathcal{R}_{b}$ as the radius of the sphere that has the same
volume as the bubble:

$\mathcal{R}_{b}\equiv\left(\frac{3V_{b}}{4\pi}\right)^{1/3}\,.$ (4)

Under the assumptions that the background density variations are statistically
isotropic and homogeneous, and that the bubble’s momentum is given by Equation
3, the bubble’s effective radius evolves in time as

$\mathcal{R}_{b}(t)=\left(\frac{3\alpha_{R}\alpha_{p}}{2\pi}\frac{\dot{p}_{w}t^{2}}{\bar{\rho}}\right)^{1/4}\equiv(\alpha_{R}\alpha_{p})^{1/4}R_{\rm
EC}\,.$ (5)

Here, $\alpha_{R}$ is another order-unity parameter that accounts for
geometry; $R_{\rm EC}$ denotes the solution for an exactly spherical bubble
expanding in a uniform ambient medium with momentum increasing at a rate
$\dot{p}_{w}$. The scaling $\mathcal{R}_{b}\propto t^{1/2}$ is shallower than
$\mathcal{R}_{b}\propto t^{3/5}$ for the classical energy-driven solution.

The total energy contained within the bubble interior is

$E_{b}=\frac{1}{2}\dot{p}_{w}\mathcal{S}\mathcal{R}_{b}\,,$ (6)

where $\mathcal{S}\geq 1$ quantifies the enhancement in energy above the case
where the bubble is fully occupied by the free wind. For a spherically
expanding bubble, we can write $\mathcal{S}$ in terms of $\alpha_{p}$ (see
Equations A13 and A16 of Paper I), with $\mathcal{S}\approx\alpha_{p}$ within
6%.

The predicted pressure in the post-shock gas of the wind is

$\displaystyle P_{b}$ $\displaystyle=$
$\displaystyle\frac{3}{16\pi}\frac{\dot{p}_{w}}{\mathcal{R}_{\rm f}^{2}}$ (7)
$\displaystyle=$
$\displaystyle\frac{3}{16\pi}\frac{\dot{p}_{w}}{\mathcal{R}_{b}^{2}}\left\\{\frac{2}{3}\alpha_{p}+\left[\left(\frac{2}{3}\alpha_{p}\right)^{2}-\frac{1}{3}\right]^{1/2}\right\\}$

where $\mathcal{R}_{\rm f}$ is the radius of the free wind region, defined
analogously to Equation 4. From Equation A14 of Paper I, the term in curly
braces can be approximated as $(\mathcal{R}_{b}/\mathcal{R}_{\rm
f})^{2}\approx(3\alpha_{p}-1)/2$ within 4% for $1<\alpha_{p}<4$.

The radial kinetic energy in the bubble’s shell is given by

$E_{\rm r,sh}=\frac{\alpha_{p}}{4\alpha_{R}}\dot{p}_{w}\mathcal{R}_{b}.$ (8)

The shell will also have turbulent motion, such that $E_{\rm sh}=E_{\rm
r,sh}+E_{\rm turb,sh}$ is its total kinetic energy. We describe the level of
turbulent energy relative to the radial kinetic energy as

$f_{\rm turb}\equiv\frac{E_{\rm turb,sh}}{E_{\rm r,sh}}\,.$ (9)

Paper I shows that the above formulae can be combined to obtain a prediction
for the fraction of input energy retained after cooling,

$1-\Theta=\left(\frac{1}{2}(1+f_{\rm
turb})\frac{\alpha_{p}}{\alpha_{R}}+\mathcal{S}\right)\frac{\dot{\mathcal{R}}_{b}}{{\cal{V}}_{w}}\,,$
(10)

where $\Theta\equiv\dot{E}_{\rm cool}/\mathcal{L}_{w}$ is the complementary
fraction lost to cooling. To the extent that the terms in parentheses in
Equation 10 are relatively constant, we would expect $1-\Theta\propto
t^{-1/2}$.

In Paper I, we hypothesize that the bubble-shell interface is a fractal. Under
this assumption, the bubble’s total area when measured on scale $\ell$ becomes

$A_{b}(\mathcal{R}_{b};\ell)\equiv
4\pi\alpha_{A}\mathcal{R}_{b}^{2}\left(\frac{\mathcal{R}_{b}}{\ell}\right)^{d}\,,$
(11)

where $d$ is known as the “excess fractal dimension” of the surface and
$\alpha_{A}$ is an order-unity parameter meant to account for any minor
inconsistencies with this model.

Paper I argued that instabilities at the interface between the shocked wind
and the surroundings drive turbulence that feeds off the kinetic energy in the
shocked wind. The turbulence in the hot gas is assumed to follow a power-law
form, such that $v_{t}(\ell)=v_{t}(L)(\ell/L)^{p}$. Accounting for both the
turbulent velocity and fractal area of the interface, the effective enthalpy
flux at scale $\ell$ is then $\Phi_{\rm cool}\sim(5/2)Pv_{\rm equiv}$, where

$v_{\rm equiv}(\ell)=v_{t}(\ell)(\mathcal{R}_{b}/\ell)^{d}\,$ (12)

defines an equivalent radial velocity through the turbulent mixing/cooling
boundary layer for a bubble of size $\mathcal{R}_{b}$. This enthalpy flux
represents the capacity for mixing and radiative cooling, increasing to
smaller scale until reaching the scale where $v_{t}(\ell_{\rm cool})t_{\rm
cool}=\ell_{\rm cool}$, where

$t_{\rm
cool}\equiv\frac{P}{n^{2}\Lambda(T_{\mathrm{pk}})}=\frac{(k_{B}T_{\mathrm{pk}})^{2}}{P\Lambda(T_{\mathrm{pk}})}\,,$
(13)

for $T_{\mathrm{pk}}$ the temperature of peak cooling. With pressure in the
shocked gas given by Equation 7, $t_{\rm cool}\propto\mathcal{R}_{b}^{2}$.

Finally, we note that Paper I outlined the conditions under which the EC
regime should apply. One estimate, from comparing to the solution from El-
Badry et al. (2019), gives an upper limit on the energy retention fraction for
the EC solution to be valid:

$1-\Theta<4\left(\frac{5}{6}\right)^{1/4}\frac{\dot{R}_{\rm
EC}}{{\cal{V}}_{w}}\,.$ (14)

Another limit is that the velocity of shocked hot gas flowing into the
boundary layer,

$v_{\rm hot}(\mathcal{R}_{b})=v_{\rm ps}\left(\frac{\mathcal{R}_{\rm
f}}{\mathcal{R}_{b}}\right)^{2}\sim\frac{{\cal{V}}_{w}}{6\alpha_{p}-2},$ (15)

cannot exceed the equivalent velocity that turbulent diffusion can
accommodate, given in Equation 12.

##  3 Numerical Methods and Models

To test the theory presented in Paper I, we use the Athena code (Stone et al.,
2008) in running a series of three-dimensional (3D) hydrodynamic simulations
of constant luminosity stellar winds injected into turbulent clouds. We do not
include any magnetic fields in the present simulations, but we do include
cooling in the gas that is implemented as part of the Athena-TIGRESS code base
for the star-forming multiphase ISM (Kim & Ostriker, 2017, 2018). This
implements cooling in an operator-split manner following Koyama & Inutsuka
(2002) at $T<10^{4.2}{\rm K}$ and Sutherland & Dopita (1993) at
$T>10^{4.2}{\rm K}$.

There is also uniform background heating of $2\times 10^{-26}\,{\rm erg}\,{\rm
s}^{-1}\,{\rm H}^{-1}$ which decreases in hot, ionized gas; more details can
be found in Section 2.3.1 of Kim & Ostriker (2017). All of our simulations are
run using the linearized Roe Riemann solver (Roe, 1981), second-order spatial
reconstruction, and the unsplit van Leer integrator (Stone & Gardiner, 2009).

###  3.1 Wind Injection

Here we describe our treatment of stellar winds, which we implement in Athena
in an operator-split fashion. We provide a quantitative test of this
implementation in Appendix A. Each simulation contains a single wind source,
which we refer to as a “star particle,” although in practice it represents a
stellar cluster. This source particle is put in “by hand” and does not exert
any gravitational force on the surrounding gas. Our approach employs a hybrid
thermal/kinetic energy injection scheme which interpolates between pure
thermal injection close to a source and pure kinetic injection towards the
edge of the feedback region. There are multiple benefits to this hybrid energy
injection approach. From a physical perspective, it more directly represents
the reality of a bubble driven by a cluster of massive stars, wherein the
majority of the energy close to the center of the bubble is thermal (due to
shock thermalization through colliding winds of individual stars) while it is
mostly kinetic towards the edges (where this thermal energy has managed to
drive expansion). From a numerical perspective, for a purely kinetic feedback
implementation the vector momentum field cannot be properly resolved
immediately adjacent to the source. By transitioning to thermal energy
injection near the source, it is only necessary to resolve a scalar field. An
advantage of the hybrid approach over purely thermal energy injection is that
the latter can require a larger spatial scale for pressure gradients to
accelerate the flow modeling the primary wind.

Given a star particle at position $\mathbf{x}_{\rm star}$ within the domain of
the simulation, we deposit energy in the grid cells surrounding the star
particle using a subcell method based on the implementation of Ressler et al.
(2020). We specify $r_{\rm fb}$, the radius of the spherical feedback region,
and $N_{\rm sc}$ is the number of subcells per grid cell along a grid-aligned
direction. When initializing the simulation we loop through a cube of $N_{\rm
sc,tot}^{3}$ subcells where

$N_{\rm sc,tot}=2N_{\rm sc}\left\lceil r_{\rm fb}/\Delta x\right\rceil$ (16)

is the total number of subcells that could possibly lie within the feedback
sphere along a grid-aligned direction and $\lceil\cdot\rceil$ denotes the
ceiling function. The subcells within this cube are equally spaced along each
grid-aligned direction between $-\Delta x\lceil r_{\rm fb}/\Delta x\rceil$ and
$+\Delta x\lceil r_{\rm fb}/\Delta x\rceil$.

For each subcell in this cube with position $\mathbf{x}_{{\rm sc},i}$ relative
to the source, we determine if $|\mathbf{x}_{{\rm sc},i}|<r_{\rm fb}+\Delta
x/2N_{\rm sc}$ (i.e. if the subcell is within half a subcell spacing of the
feedback radius or fully within the feedback radius). If it is, we then record
the position of the subcell $\mathbf{x}_{{\rm sc},i}$ for use within the full
simulation. In order to account for the subcell volume only partially
overlapping the actual feedback region, we additionally assign each subcell a
weight $w_{{\rm sc},i}$ according to

$w_{{\rm sc},i}=\left\\{\begin{array}[]{ll}1&\quad|\mathbf{x}_{{\rm
sc},i}|\leq r_{\rm fb}-\frac{\Delta x}{2N_{\rm sc}}\\\
\frac{1}{2}\left(\frac{2N_{\rm sc}(r_{\rm fb}-|\mathbf{x}_{{\rm
sc},i}|)}{\Delta x}+1\right)&\quad|\mathbf{x}_{{\rm sc},i}|>r_{\rm
fb}-\frac{\Delta x}{2N_{\rm sc}}\end{array}{}\right.\,.$ (17)

This formula is simply a linear interpolation between 1 and 0 over the radial
range of $\pm$ half a subcell spacing from the feedback radius, which
approximates the true volume fraction (Jones & Williams, 2017). We
additionally store the total effective volume that wind energy is injected
into, defined as

$V_{\rm inj}=\left(\frac{\Delta x}{N_{\rm sc}}\right)^{3}\sum_{i}w_{{\rm
sc},i}\,.$ (18)

Once the positions and weights of the subcell template are initialized, they
are used for the hybrid thermal/kinetic energy injection.

Given a mechanical luminosity $\mathcal{L}_{w}$ and mass loss rate
$\dot{M}_{w}$ for the wind, the mean mass density $\rho_{w}$ and energy
density $\epsilon_{w}$ to be injected in a time $\Delta t$ are

$\rho_{w}=\frac{\dot{M}_{w}\Delta t}{V_{\rm
inj}}\,\,,\,\,\epsilon_{w}=\frac{\mathcal{L}_{w}\Delta t}{V_{\rm inj}}\,.$
(19)

For a given subcell in our template, indexed by $i$, we determine the grid
cell in which it resides at location $\mathbf{x}_{\rm star}+\mathbf{x}_{{\rm
sc},i}$, and increment the mass density by $\delta\rho_{i}=w_{{\rm
sc},i}\rho_{w}/N_{\rm sc}^{3}$. For the purposes of the current paper,
$\mathbf{x}_{\rm star}$ is always at the origin.

The total energy density in the target cell is similarly incremented by
$\delta\epsilon_{i}=w_{{\rm sc},i}\epsilon_{w}/N_{\rm sc}^{3}$. To compute the
injected momentum, we specify the variable $f_{\rm kin}$ at a given subcell
position $\mathbf{x}_{{\rm sc},i}$ according to

$f_{\rm kin}(|\mathbf{x}_{{\rm sc},i}|)={\rm max}\left(1,|\mathbf{x}_{{\rm
sc},i}|/r_{\rm fb}\right)\,,$ (20)

Once $f_{\rm kin}$ is determined we define the injected momentum density as

$\delta\mathbf{q}_{i}=\left(2f_{\rm kin}(|\mathbf{x}_{{\rm
sc},i}|)\delta\epsilon_{i}\delta\rho_{i}\right)^{1/2}\frac{\mathbf{x}_{{\rm
sc},i}}{|\mathbf{x}_{{\rm sc},i}|}\,.$ (21)

As was noted in Wall et al. (2020), the presence of momentum and mass within
the grid cell to which we are adding these quantities means that the kinetic
energy added by the above change in momentum density is not $f_{\rm
kin}\delta\epsilon_{i}$, but actually somewhat less depending on what the
initial density and momentum of the cell are. The discrepancy between the
actual injected kinetic energy density and $f_{\rm kin}\delta\epsilon_{i}$ is
especially important at the early times of wind onset, when the surrounding
density is high, but becomes negligible once the wind has evacuated the
feedback region. In the case where there is negligible pre-existing energy and
momentum, our prescription deposits 75% of the total energy as kinetic energy.

In order to make sure that our results are not affected by the details of our
feedback mechanism, we run tests with $f_{\rm kin}=0$ (purely thermal
feedback) in Appendix B.1. This scenario might be closer to the truth of a
cluster of stars whose colliding winds shock and turn initially kinetic energy
into bulk thermal energy. We found that this pure thermal feedback scenario
simply leads to a minor delay in the wind’s evolution, but broadly the results
remain unchanged.

When we inject the wind mass into the feedback region we add an equal amount
of “mass” to a passive scalar variable, which is simply passively advected
with the gas. We use this scalar to track the mixing of the wind material with
the surrounding gas as well as a means of separating mass that has been swept
up by the wind from the ambient medium, as we will explain in Section 4.

Table 1: Parameters of simulation suite.

Cloud Mass | Cloud Radius | $\overline{n}_{\rm H}$ | $v_{t}$ |  $\Delta x$aaOnly highest resolution is provided. | ResolutionaaOnly highest resolution is provided.
---|---|---|---|---|---
$[M_{\odot}]$ | $[{\rm pc}]$ | $[{\rm cm}^{-3}]$ | $[{\rm km}\,{\rm s}^{-1}]$ | $[{\rm pc}]$ |
$5\times 10^{4}$ | 20 | 43.1 | 3.59 | $0.15$ | $256^{3}$
$5\times 10^{4}$ | 10 | 345. | 5.08 | $0.08$ | $256^{3}$
$5\times 10^{4}$ | 5 | 2760 | 7.18 | $0.04$ | $256^{3}$
$5\times 10^{4}$ | 2.5 | 22800 | 10.2 | $0.02$ | $256^{3}$
$10^{5}$ | 20 | 86.3 | 5.08 | $0.08$ | $512^{3}$
$10^{5}$ | 10 | 690. | 7.18 | $0.04$ | $512^{3}$
$10^{5}$ | 5 | 5520 | 10.2 | $0.02$ | $512^{3}$
$10^{5}$ | 2.5 | 44200 | 14.4 | $0.01$ | $512^{3}$
$5\times 10^{5}$ | 20 | 431 | 11.4 | $0.15$ | $256^{3}$
$5\times 10^{5}$ | 10 | 3450 | 16.1 | $0.08$ | $256^{3}$
$5\times 10^{5}$ | 5 | 27600 | 22.7 | $0.04$ | $256^{3}$
$5\times 10^{5}$ | 2.5 | 228000 | 32.1 | $0.02$ | $256^{3}$

Note. — Each model is run with three different values for the mass of the wind
source particle, $M_{*}/M_{\rm cloud}\equiv\varepsilon_{*}=0.01,0.1,$ and $1$.

###  3.2 Model Parameters and Set-Up

Our goal is to span a large range of the potential parameter space for stellar
winds originating from young star clusters within massive molecular clouds,
which is where most star formation takes place. To that end, we consider three
cloud masses $M_{\rm cloud}$ of $5\times 10^{4}M_{\odot}$, $10^{5}M_{\odot}$,
and $5\times 10^{5}M_{\odot}$ and four different cloud radii $R_{\rm cloud}$
of $20\,{\rm pc}$, $10\,{\rm pc}$, $5\,{\rm pc}$, and $2.5\,{\rm pc}$. The
largest radii represent “normal” star-forming GMCs, cases with intermediate
radii may represent infrared dark clouds (IRDCs), and cases with the smallest
radii may represent extremely dense cluster-forming clumps within larger GMCs
or the natal clouds in which super star clusters are born in starbursting
galactic centers.

Each choice of $M_{\rm cloud}$ and $R_{\rm cloud}$ corresponds to a choice of
mean mass density $\bar{\rho}$ as

$\bar{\rho}=\frac{3M_{\rm cloud}}{4\pi R_{\rm cloud}^{3}}=\mu
m_{p}\overline{n}_{\rm H}\;,$ (22)

where $\mu_{\rm H}=1.4271$ is the mean molecular weight of the gas, $m_{p}$ is
the mass of a proton, and $\overline{n}_{\rm H}$ is the mean number density of
the Hydrogen nuclei. Rather than creating spherical clouds within our cubic
domain, we apply uniform density $\bar{\rho}$ everywhere within a domain of
side length $L_{\rm box}=2R_{\rm cloud}$. The “global cloud” parameters are
just used for reference in setting the ambient conditions with which the wind
interacts.

The grid has uniform spatial resolution in each dimension with $L_{\rm
box}/\Delta x=128,\,256,\,512$ for the simulations with $M_{\rm
cloud}=10^{5}\,M_{\odot}$ and only $L_{\rm box}/\Delta x=128,\,256$ for the
simulations at higher and lower cloud mass. The parameters describing the
simulation initial density conditions and numerical resolution are outlined in
Table 1. Only the highest resolution is listed for each model.

We initialize all the simulations at a uniform pressure of $P/k_{B}=3\times
10^{4}\,{\rm cm}^{-3}\,{\rm K}$. Since we evolve the simulations long enough
for thermal relaxation to occur before initiating the wind, this particular
choice of initial pressure is not important to our results. We also ran test
simulations at a initial pressure of $P/k_{B}=3\times 10^{3}\,{\rm
cm}^{-3}\,{\rm K}$ and found no differences.

Each simulation is initialized with a turbulent velocity field which is a
Gaussian random field with power spectrum $|{v}_{k}|^{2}\propto k^{-4}$ for
$2\leq kL_{\rm box}/(2\pi)\leq 64$. For higher resolution runs, the
realization of the power spectrum uses an identical set of modes to the lower-
resolution case. The amplitude of the turbulence is chosen so that the initial
kinetic energy per unit mass is equivalent to twice that of the gravitational
potential for a sphere of radius $R_{\rm cloud}$ and mass $M_{\rm cloud}$,
i.e.

$\tilde{E}_{\rm kin,i}=2|\tilde{W}_{\rm sphere}|=\frac{6GM_{\rm
cloud}}{5R_{\rm cloud}}\,,$ (23)

where the tildes are used to indicate energy per unit mass. We note, however,
that self-gravity is not included in the simulations.

For each simulation, we allow the turbulent velocity field to evolve and
decay, without any additional driving, until the average kinetic energy per
unit mass has decayed to $\tilde{W}_{\rm sphere}$ (i.e. it has been reduced by
half). At the end of this initial evolution each simulation has a turbulent
velocity scale $v_{t}=(\tilde{E}_{\rm kin,i})^{1/2}$. For each model, $v_{t}$
is listed in Table 1.

At the same time as turbulence is decaying, the Reynolds stresses create
inhomogeneous density structure throughout the simulation domain. The
turbulent decay period is generally one tenth of the initial flow crossing
time ($L_{\rm box}/v_{t,i}$) of the simulation domain, amounting to 0.03-0.7
Myr. Over this period, the thermal energy also relaxes, with temperature
approaching thermal equilibrium; in the cold gas this is $\sim 30-200\,{\rm
K}$.

After the initial period of turbulent decay, we create a star particle at the
center of the simulation domain of mass $M_{*}=\varepsilon_{*}M_{\rm cloud}$.
The parameter $\varepsilon_{*}$ can be thought of as akin to a total star
formation efficiency. For each set of initial conditions laid out in Table 1
we run simulations at $\varepsilon_{*}=0.01,0.1,$ and $1$. The star particles
inject energy and mass proportional to their own mass in accordance with the
method described in Section 3.1, with a feedback radius $r_{\rm fb}=L_{\rm
box}/40$ at all resolutions, meaning that the feedback radius is resolved by a
minimum of 3.2 resolution elements at the lowest resolution. We test this
choice of feedback radius in Appendix B.3 and find it does not affect our
results.

Figure 1: Ratio of the wind parameters derived using the population synthesis
code STARBURST99 (SB99) to the constant values for each parameter that we
adopt in our simulations.The SB99 mass-loss rate and momentum input rate are
respectively a factor $\sim 5$ and $\sim 2$ lower than the ones we adopt. We
validate that using the SB99 values does not change our conclusions (see
Appendix B.2).

We adopt a constant wind luminosity per unit mass of
$\mathcal{L}_{w}/M_{*}=10^{34}\,{\rm erg}\,{\rm s}^{-1}\,{M_{\odot}^{-1}}$ and
mass loss rate per unit mass $\dot{M}_{w}/M_{*}=10^{-2}\,{\rm Myr}^{-1}$. This
mass loss rate is about a factor of five larger than that expected for a
standard solar metallicity population, with a Kroupa IMF (Kroupa, 2001) as
calculated by STARBURST99 (SB99) (Leitherer et al., 1999). A comparison
between the constant values that we adopt and those determined from SB99 is
given in Figure 1. We use the larger mass loss rate because this reduces the
temperature of the shocked wind, yielding a less stringent time step.

Our adopted wind parameters correspond to a wind velocity of
${\cal{V}}_{w}=(2\mathcal{L}_{w}/\dot{M}_{w})^{1/2}=1780\,{\rm km}\,{\rm
s}^{-1}$ and a momentum injection rate per unit mass of
$\dot{p}_{w}/M_{*}={\cal{V}}_{w}\dot{M}_{w}/M_{*}=17.8\,{\rm km}\,{\rm
s}^{-1}\,{\rm Myr}^{-1}$. At the lower mass loss rates consistent with the
SB99 calculations, the wind velocity (momentum injection rate) is larger
(smaller) by a factor $\sim\sqrt{5}\approx 2.2$, i.e. ${\cal{V}}_{w}\sim
3750\mathrm{\,km\,s^{-1}}$ and $\dot{p}_{w}/M_{*}\sim
8.6\mathrm{\,km\,s^{-1}}\,\mathrm{\,Myr}^{-1}$. This only corresponds to a
difference in radial evolution (according to our theory) of a factor of
$5^{1/8}\approx 1.2$, so this should not effect the interpretation of our
results. We have additionally run simulations with a standard mass loss rate
and have obtained consistent results, as described in Appendix B.2.

Figure 2: An example snapshot of the wind bubble structure and phase
properties at time $t=2.23\,{\rm Myr}$ after the initiation of the wind for
the model with $M_{\rm cloud}=10^{5}\,M_{\odot}$, $R_{\rm cloud}=20\,{\rm
pc}$, $\varepsilon_{*}=1\%$, and $L_{\rm box}/\Delta x=512$ simulation. The
left and center columns show slices through the $z=0$ plane, while the right
column displays distributions. From left to right and top to bottom, panels
show: gas number density, pressure, a cooling-weighted phase diagram of
cooling time vs. pressure, radial momentum density, temperature, a volume-
weighted temperature-density phase diagram, the wind mass fraction ($f_{\rm
wind}$), the cooling rate per unit volume, and a cooling-weighted histogram of
the gas temperature distribution (which depicts the total cooling occurring in
each temperature bin). In the top-right panel, the highest concentration of
cooling is in gas that has $t_{\rm cool}\propto P^{-1}$; this slope is
indicated by a dashed red diagonal line for the coefficient in Equation 13
corresponding to $T_{\rm pk}=9650\,{\rm K}$. The last panel uses vertical
lines to delineate the phases described in Table 2, where the subscripts from
the table are indicated. The label “wind” includes both the free (f) and post-
shock (ps) wind phases. In this panel, $T_{\rm pk}$ is marked as a vertical
dashed line within the “w” phase.  Figure 3: Same as Figure 2 but for the
model with $M_{\rm cloud}=10^{5}\,M_{\odot}$, $R_{\rm cloud}=2.5\,{\rm pc}$,
$\varepsilon_{*}=100\%$, and $L_{\rm box}/\Delta x=512$ simulation at
$t=0.125\,{\rm Myr}$. All quantities are as specified in Figure 2, though the
range of many quantities are quite different.

##  4 Numerical Results

Here we lay out the results of the simulations described in Section 3.2. In
our presentation, we shall focus on cases with $M_{\rm
cloud}=10^{5}\,M_{\odot}$ and use the remaining simulations for validation of
the generality of our results, as summarized in Appendix C.

###  4.1 Cloud and Bubble Structure and Thermal Distributions

Example snapshots from our simulations are provided in Figure 2 ($R=20\,{\rm
pc}$ and $\varepsilon_{*}=1\%$ model) and Figure 3 ($R=2.5\,{\rm pc}$ and
$\varepsilon_{*}=100\%$ model) in order to provide a more concrete reference.
These two models are deliberately chosen to span the range of cloud size and
star formation efficiency parameters.

In both figures, the “wind-blown bubble” is the low-density, high-temperature
region filling the central region of the simulation domain, with fingers
extending outward. Outside of the bubble, the inhomogeneous density structure
generated by the background turbulence in the cold cloud is evident; higher
density portions of the ambient gas are also left behind in the bubble
interior, although they are ablated over time by KH instabilities as the high-
velocity wind flows past them.

As expected, the cooling rate (bottom center panel) is greatest at the
(fractal) interface between the hot bubble and cool shell, where mixing is
driven by turbulence. The outward momentum originally carried by the wind is
deposited at the bubble boundary by interface mixing, where it builds up in
the expanding cool shell (middle left panel). The degree of mixing is also
evident in the fraction of wind gas in each cell (bottom left panel).

Figure 2 and Figure 3 also include (right column) information on the
statistical distribution functions of gas in density, temperature, pressure,
and cooling time. The top right panel shows that the highest concentration of
cooling is in gas with $t_{\rm cool}\propto P^{-1}$, as predicted in Paper I
(see definition in Equation 13) . Quantitatively, we find that

$t_{\rm cool}=135\mathrm{\,yr}\left(\frac{P/k_{B}}{10^{6}\,{\rm K}{\rm
cm}^{-3}}\right)^{-1}$ (24)

follows well the prominent linear feature at short cooling time for both
cases, as well as our other models. The corresponding $T_{\rm pk}$ would be
$9650\,{\rm K}$.

Table 2: Definitions of gas phases.

Phase | Temperature Condition | Velocity Condition | Subscript | Schematic Region
---|---|---|---|---
Free Hyper-Sonic Wind | $T>10^{5}\,{\rm K}$ | $v_{r}>{\cal{V}}_{w}/2$ | ${\rm f}$ | ${\rm(i)}$
Shocked Stellar Wind | $T>10^{6}\,{\rm K}$ | $v_{r}<{\cal{V}}_{w}/2$ | ${\rm ps}$ | ${\rm(ii)}$
Ionized Gas | $2\times 10^{4}\,{\rm K}<T<10^{6}\,{\rm K}$ | $v_{r}<{\cal{V}}_{w}/2$ | ${\rm i}$ | ${\rm(iii)}$
Warm Neutral Gas | $5050\,{\rm K}<T<2\times 10^{4}\,{\rm K}$ |  | ${\rm w}$ | ${\rm(iii)/(iv)}$
Thermally Unstable Gas | $184\,{\rm K}<T<5050\,{\rm K}$ |  | ${\rm u}$ | ${\rm(iii)/(iv)}$
Cold Neutral Gas | $T<184\,{\rm K}$ |  | ${\rm c}$ | ${\rm(iv)}$

Note. — Names assigned to gas phases based on temperature and velocity
conditions, with (second-from-right column) subscript used to denote each
phase. The final column indicates the rough correspondence between gas phases
and locations in the schematic Figure 1b of Paper I.

###  4.2 Gas Phase Definitions

We define several different gas phases as an aid in quantifying the structure
and evolution of the wind-driven bubbles. These definitions can be thought of
as an expansion on the phases laid out for the classical stellar wind bubble
in Paper I, in order to better account for the cooling of the gas. These
definitions (see Table 2) are based on the temperature $T$ and radial velocity
$v_{r}$ of the gas. The table also lists subscripts used to denote quantities
associated with each phase (such as volume, momentum, energy, etc.). These
subscripts label the thermal phases of the wind in the bottom right panels of
Figure 2 and Figure 3.

The first two phases listed in Table 2 are analogous to the two wind phases
described in Weaver et al. (1977). The Ionized Gas is contained within the
cooling layer, and cooling in this phase is dominated by collisional
excitation of transitions of H, He, C, N, and O (primarily). This phase is
produced via mixing and subsequent cooling of shocked wind and shell gas, and
would not exist as part of the undisturbed/ambient ISM in either the uniform
or turbulent cases.

The Warm Neutral Gas has cooling dominated by collisionally-excited Ly$\alpha$
emission and recombination on dust grains. The Thermally Unstable Gas is in
the range between the stable equilibrium warm and cold phases in the static
case, for our adopted cooling function. This phase is continually populated
because of the turbulence in the system. Finally, at the lowest temperature
there is the Cold Neutral Gas, which in reality would include both atomic and
molecular phases, but here we do not follow the detailed chemistry.

The final three phases would all usually exist to varying degrees as part of
the background in a turbulent, dense cloud. We wish to separate the portions
that are ISM gas that has been shocked and then cooled after being swept into
the expanding bubble shell from the portions that are undisturbed ambient gas.
We do this using cuts based on the fraction of the mass in a given cell that
originates in the wind, $f_{\rm wind}$. The wind mass is tracked using a
scalar that is injected with the wind (see Section 3.1) and then passively
advected.

For each of the Ionized, Warm, Thermally Unstable, and Cold phases of the gas
we track all quantities of interest in gas with $f_{\rm
wind}>10^{-3},10^{-4},$ and $10^{-5}$ (each subsequent selection being a
superset of the previous selection). Unless otherwise stated, we use $f_{\rm
wind}>10^{-4}$ to separate the swept-up and ambient parts of these phases, as
this gives the best agreement between different methods of measuring the total
cooling, as described in Section 4.7. The tracking of gas properties at
different $f_{\rm wind}$ values also allows us to quantify the mixing of the
wind with the turbulent gas.

As is evident in the bottom right panels of Figure 2 and Figure 3, the
majority of the cooling is occurring in the Warm and Ionized phases. These are
the phases that primarily occupy the boundary region between the bubble and
surrounding cloud, as is clear from the temperature and cooling slices in
these figures.

###  4.3 Bubble Evolution Comparisons

Here we present results from our simulations of the temporal evolution of the
radial momentum of the gas ($p_{r}$), the wind bubble’s effective radius
($\mathcal{R}_{b}$), and the interior bubble energy ($E_{b}$). We compare
these results with the theory developed in Paper I and reviewed in Section 2.
When measuring quantities in the shell of the wind bubble, we include the
measurement of this quantity in all gas with a wind mass fraction greater than
$10^{-4}$. We found this was the best way of differentiating between the
swept-up shell and the background gas.

In this section we show results for the $M_{\rm cloud}=10^{5}\,M_{\odot}$
cases, while results for other cloud masses are presented in Appendix C
(Figure 24 and Figure 25).

Figure 4: Evolution of the total radial momentum of the bubble (interior plus
shell) in each of our simulations for cases with cloud mass $M_{\rm
cloud}=10^{5}M_{\odot}$. Model cloud radii (from least dense, $R_{\rm
cloud}=20\,{\rm pc}$, to most dense, $R_{\rm cloud}=2.5\,{\rm pc}$) are given
within the individual panels (a) to (d); see Table 1 for parameters. For each
case we show results for star cluster particles of mass
$M_{*}=10^{3}M_{\odot},\,10^{4}M_{\odot},\,$and$\,10^{5}M_{\odot}$
(corresponding to star formation efficiencies of
$\varepsilon_{*}=1\%,\,10\%,\,$and$\,100\%$) in yellow, orange, and red
respectively. At each value of $R_{\rm cloud}$ and $\varepsilon_{*}$ we show
three separate lines (solid, dashed, dotted) corresponding to varying
resolution, as shown in the key. For each simulation with $512^{3}$ cells, we
mark the time at which the first wind-polluted gas exits the simulation domain
$t_{\rm esc}$ with a square, and only show results from each simulation up
until $4t_{\rm esc}$ (though not all of the highest resolution simulations
were run for this long). For each value of $R_{\rm cloud}$ and
$\varepsilon_{*}$ we show $\dot{p}_{w}t$ (which should apply when cooling is
maximally efficient, per Equation 3) in black.

####  4.3.1 Shell Momentum

As is evident in Figure 2 and Figure 3, most of the radial momentum is carried
by the dense shell of swept-up gas, thus we define the momentum as:

$p_{r}\equiv p_{r,{\rm i}}+p_{r,{\rm w}}+p_{r,{\rm u}}+p_{r,{\rm c}}.$ (25)

As mentioned above, only cells with a wind mass fraction greater than
$10^{-4}$ are included in this definition. We found no significant changes in
our measurement when considering $f_{\rm wind}>10^{-5}$ gas or even $f_{\rm
wind}>0$ gas, which we took to indicate that there was no significant momentum
carried in “completely unpolluted” gas. In Figure 4 (for the $M_{\rm
cloud}=10^{5}\,M_{\odot}$ models), we compare this measurement of the momentum
in the simulations with $p_{r}=\dot{p}_{w}t$, as predicted by the EC theory
with $\alpha_{p}=1$ (Equation 3).

The results at different resolutions show that the radial momentum carried by
the bubble is extremely well converged in our simulations. Overall, results
are within a factor of 1.2-4 of the EC prediction. The results are
quantitatively closest to $\dot{p}_{w}t$ for more luminous winds
(corresponding to higher star formation efficiency).

The momentum in excess of $\dot{p}_{w}t$ can be attributed to nonzero buildup
of thermal energy in the shocked wind that aids in driving expansion (see
Appendix of Paper I). This is parameterized in our theory by $\alpha_{p}$ in
Equation 3.

We note that once the bubble has reached the edge of the simulation domain
(identified by nonzero outflow of wind-contaminated gas), we cannot expect the
numerical solution to continue to follow the EC theoretical prediction. We
mark this first “blowout” time in Figure 5 with large squares (for the
highest-resolution models). Some time after blowout, the rate of momentum
increase falls below $\dot{p}_{w}$ because a fraction of the wind exits the
domain without interacting with cloud gas.

Figure 5: Evolution of the effective radius $\mathcal{R}_{b}$ in each of our
simulations for cases with cloud mass $M_{\rm cloud}=10^{5}M_{\odot}$. Line
styles are as in Figure 4. For each value of $R_{\rm cloud}$ and
$\varepsilon_{*}$ we show (in black) the evolution prediction from the EC
theory (Equation 5, with $\alpha_{R}=\alpha_{p}=1$). The grey horizontal lines
indicate $R_{\rm cloud}$, while squares indicate the time when wind-
contaminated gas first leaves the domain ($t_{\rm esc}$). As in Figure 4, we
only show simulation results up until $4t_{\rm esc}$.

####  4.3.2 Effective Radius

The effective bubble radius $\mathcal{R}_{b}$ is defined from Equation 4 based
on the bubble volume, $V_{b}$. For this volume, we include just the hot phases
of the gas, the Free Wind (volume $V_{\rm f}$) and the Shocked Wind (volume
$V_{\rm ps}$):

$V_{b}\equiv V_{\rm f}+V_{\rm ps}\,.$ (26)

Comparisons between our simulations with $M_{\rm cloud}=10^{5}\,M_{\odot}$ and
the prediction given by Equation 5 (with $\alpha_{R}=\alpha_{p}=1$) are given
in Figure 5. It is clear that the EC theory very well explains the salient
features of the radial expansion. The only exceptions are the smallest two
clouds for $\epsilon_{*}=0.01$, which is the least realistic parameter regime
since high density regions are observed to have very high SFE (e.g. Leroy et
al., 2018). These deviations are mainly caused by the Reynolds stress in the
surrounding turbulent gas being non-negligible compared to the wind pressure,
given the low $\varepsilon_{*}$.

We note in particular that the scaling $\mathcal{R}_{b}\propto t^{1/2}$
represents the numerical results better than the steeper scaling $R\propto
t^{3/5}$ of Weaver et al. (1977) or El-Badry et al. (2019) for pressure-driven
expansion. For the majority of cases that follow the $\mathcal{R}_{b}\propto
t^{1/2}$ scaling well, the numerical result is within $20\%$ of Equation 5.
Quantitatively, the EC theory works especially well at higher densities (where
cooling is most efficient) and at higher $\mathcal{L}_{w}$ (i.e. higher
$\varepsilon_{*}$, where $\dot{p}_{w}$ and $\dot{{\cal R}}_{\rm b}$ are
larger, and therefore the constraint on $\Theta$ as given by Equation 14 is
not as stringent).

At the time of first blowout, the mean bubble expansion rates range from $\sim
10-100\mathrm{\,km\,s^{-1}}$ in the least dense cloud to $\sim
3-10\mathrm{\,km\,s^{-1}}$ in the most dense cloud, with larger velocities
applying in the cases with higher wind power
$\mathcal{L}_{w}\propto\varepsilon_{*}$. Note that the first blowout (shown
with colored squares) occurs much earlier than the time when
$\mathcal{R}_{b}(t)=R_{\rm cloud}$ (shown with a horizontal grey line) due to
the fractal nature of the bubble interface, where parts of the bubble surface
are at much larger radii than others.

Figure 6: Evolution of the total energy in the interior of the bubble in each
of our simulations for cases with cloud mass $M_{\rm cloud}=10^{5}M_{\odot}$.
Line styles are as in Figure 4. For each value of $R_{\rm cloud}$ and
$\varepsilon_{*}$ we show in black the energy given by Equation 6 (taking
$\mathcal{S}=\alpha_{p}=\alpha_{R}=1$ and using Equation 5).

####  4.3.3 Bubble Interior Energy

In keeping with the definition of the bubble interior given by Equation 26, we
measure the bubble’s interior’s energy using

$E_{b}\equiv E_{{\rm f},{\rm kin}}+E_{{\rm f},{\rm thm}}+E_{{\rm ps},{\rm
kin}}+E_{{\rm ps},{\rm thm}}$ (27)

where the “${\rm kin}$” and “${\rm thm}$” subscripts refer to kinetic and
thermal energy respectively. We note that the thermal energy within the free
wind ($E_{{\rm f},{\rm thm}}$) is negligible compared to its kinetic energy,
but we include it here for completeness.

Figure 6 compares the energy measurement from the simulations with the
prediction given by Equation 6, employing Equation 5 and taking
$\mathcal{S}=\alpha_{p}=\alpha_{R}=1$. When comparing the theory and
simulations, we see similar trends as those observed in the evolution of
bubble radii and momenta: the EC theory is most accurate at higher
$\varepsilon_{*}$ and higher density $\bar{\rho}$.

We also note that given the agreement with theoretical predictions for the
shell momentum (as evidenced by Figure 4) and the bubble internal energy (as
evidenced by Figure 6), we expect the shell’s radial kinetic energy to be half
that of the bubble’s interior energy, as predicted in Paper I (see Equation 23
there). This is indeed the case.

Figure 7: Evolution of the pressure in the shocked wind gas in cases with
cloud mass $M_{\rm cloud}=10^{5}M_{\odot}$. All panels show the evolution of
the volume averaged mean pressure (solid colored lines) along with shading the
$16^{\rm th}$ to $84^{\rm th}$ percentiles of the distribution of pressures in
post-shock gas (colored shaded regions). We show the theoretical prediction
(Equation 7, with $\alpha_{p}=1$) in black. The time at which wind first
escapes the domain is indicated by a colored square.

####  4.3.4 Hot Gas Pressure

Finally, we compare the prediction for the pressure in the shocked wind gas
given by Equation 7 and that measured in our simulations. We measure the
pressure by selecting the shocked wind gas (conditions specified in Table 2)
in simulation snapshots and calculating the volume-averaged mean thermal
pressure. We use the $L_{\rm box}/\Delta x=256$ simulations for these
comparisons as the pressure is computed directly from snapshots, which are
taken at higher cadence in the lower resolution simulations. We also compute
the $16^{\rm th}$, $50^{\rm th}$, and $84^{\rm th}$ percentiles of the
distribution of pressures in this gas.

In Figure 7 we compare the evolution of the pressure with the prediction of
Equation 7. For this comparison, we use $\alpha_{p}=\alpha_{R}=1$, and we see
that the theoretical prediction closely matches the evolution in the
simulations. As we discuss below, the simulations actually have
$\alpha_{p}>1$. From Equation 7 this would tend to increase the theoretical
prediction (black curves in Figure 7), bringing the prediction closer to the
‘observed’ simulation value for the low $\varepsilon_{*}$ cases (which have
the largest $\alpha_{p}$), but further away from agreement for high
$\varepsilon_{*}$ cases. The missing component here is the obliquity of the
shock surface: the more oblique the shock the less thermalised the shocked
wind becomes, hence lowering the pressure (discussed at the end of Appendix A
in Paper I). As is clear from Figure 2 and Figure 3, the high
$\varepsilon_{*}$ winds have more oblique shocks, explaining this discrepancy.

Figure 8: Evolution of the momentum ($\alpha_{p}$) and energy ($\mathcal{S}$)
enhancement factors shown as solid and dashed lines respectively. The
evolution of $\alpha_{p}$ is only shown up to the point of breakout, which is
indicated with colored squares on the curves for $\mathcal{S}$. We indicate
the ‘no enhancement’ case with the horizontal black line. We show the highest
resolution version of all simulations with $M_{\rm cloud}=10^{5}\,M_{\odot}$.
Figure 9: Evolution of the dimensionless parameter $\alpha_{R}$. Line colors
are associated with simulations of different $\varepsilon_{*}$. We use the
highest resolution version of all simulations with $M_{\rm
cloud}=10^{5}\,M_{\odot}$ and only show the evolution up to the point of
breakout (square markers in Figure 4 and others). Figure 10: Evolution of the
ratio of turbulent kinetic energy to kinetic energy in radial motion, for
shell gas. We show the highest resolution version of all simulations with
$M_{\rm cloud}=10^{5}\,M_{\odot}$. The square symbols mark the point at which
wind-polluted mass begins to leave the domain.

###  4.4 Dimensionless Parameters

In the above comparisons we explicitly set the dimensionless, order-unity
parameters of our theory equal to unity. However, the values and time
evolution of these parameters provides interesting insight into the validity
of the assumptions in the EC theory for different regimes (see Section 4.3).
As is explained in Appendix A of Paper I and summarized in Section 2, some of
these parameters are interdependent.

The first dimensionless parameter is $\alpha_{p}$, which is defined as the
rate of momentum input to the surrounding medium divided by the rate of
momentum injection by the wind. Specifically, we measure

$\alpha_{p}=\frac{p_{\rm r,b}}{\dot{p}_{w}t}$ (28)

where $p_{r}$ is as measured in Equation 25. This measurement of $\alpha_{p}$
is shown as solid lines in Figure 8.

Another quantity in the EC theory is the energy enhancement factor
$\mathcal{S}$, defined in Equation 6. As explored in depth in Appendix A of
Paper I, $\mathcal{S}$ and $\alpha_{p}$ are expected to track each other (with
$\mathcal{S}\approx\alpha_{p}$) because both are associated with buildup of
energy in the bubble interior. For the same reason, we expect $\mathcal{S}$ to
be higher when the shocked wind makes up a larger portion of the bubble
volume, i.e. larger $\mathcal{R}_{b}/\mathcal{R}_{\rm f}$. We also expect
larger $\mathcal{S}$ when shock surfaces are less oblique. In order to
directly compare $\mathcal{S}$ with $\alpha_{p}$, we show in Figure 8 a
measured value of $\mathcal{S}$ as dashed lines.

The measured value of $\mathcal{S}$ is taken as the ratio of the bubble energy
in the simulations, measured according to Equation 27, to
$\dot{p}_{w}\mathcal{R}_{b}/2$, where $\dot{p}_{w}$ is taken as the fixed wind
momentum input rate and $\mathcal{R}_{b}$ is computed from the simulations as
detailed in Section 4.3.2.

The next dimensionless parameter we introduce is $\alpha_{R}$, which encodes
geometric factors. Using Equation 5,

$\alpha_{R}=\frac{\mathcal{R}_{b}^{4}}{\alpha_{p}}\left(\frac{3}{2\pi}\frac{\dot{p}_{w}t^{2}}{\bar{\rho}}\right)^{-1}\,.$
(29)

We evaluate this using the measured $\mathcal{R}_{b}$ (see Section 4.3.2) and
$\alpha_{p}$ (Equation 28). We note that we have no explicit theoretical
prediction for $\alpha_{R}$ other than expecting it to be near unity. Figure 9
shows the measured values of $\alpha_{R}$ over time. We see that $\alpha_{R}$
remains quite close to unity, remaining in the range $0.5-1.5$ for the
majority of the evolution. There is also some indication that $\alpha_{R}$
takes on lower values in the higher density clouds.

The radial kinetic energy of the shell in the EC theory is determined by
momentum input from the wind, but there is no prediction for the kinetic
energy in non-radial, turbulent motion. Given that the theory relies on
efficient cooling facilitated through a turbulent interface, we expect a
significant fraction of the kinetic energy in the shell (and around the shell-
bubble interface) to be in turbulent motion. The simplest way to quantify this
is in terms of the $f_{\rm turb}$ parameter defined in Equation 9.

To measure this parameter we write the total radial kinetic energy in the
shell as

$E_{r,{\rm sh}}\equiv\sum_{j\in{\rm i,w,u,c}}\frac{p_{r,j}^{2}}{2M_{j}}$ (30)

where the sum is performed over the indicated phases, and as in Equation 25 we
sum over zones where $f_{\rm wind}>10^{-4}$. The total kinetic energy is
measured in an analogous way.

The ratio $f_{\rm turb}=E_{\rm tot,sh}/E_{\rm r,sh}-1$ is displayed in Figure
10. Surprisingly, for much of the bubble evolution and most of the parameter
space, the energy in turbulent motion is at least as large as the energy in
radial motion, with $f_{\rm turb}\approx 1$ (equipartition between radial and
turbulent motion) over much of the parameter space. This fact emphasizes how
important the turbulence is for the evolution of the bubble, especially the
turbulence driven by the wind. We observe that the more powerful winds instill
a larger fraction of their kinetic energy in radial motion.

In all simulations the fraction of turbulent kinetic energy in the shell
increases in time until seemingly reaching a set value before decreasing
again. As we will see upon further inspection of the turbulent motion below,
this reflects a saturation of turbulence as more and more radial energy is
provided to the wind.

Figure 11: Magnitude of the tangential velocity in snapshots of the $R_{\rm
cloud}=20\,{\rm pc}$ models at the time in each simulation when
$\mathcal{R}_{b}=R_{\rm cloud}/3$. The separate columns (left to right) are
the $\varepsilon_{*}=1\%,\,10\%,\,$and $100\%$ cases. The hot gas (Shocked
Stellar Wind and Ionized Gas) is shown with the purple-to-yellow color scale
while the cold gas (Warm, Thermally Unstable, and Cold Gas) is shown with the
green color scale. The top row of panels is a zoom-in of the bottom row, which
itself zooms in on just the bubble region, to highlight the structure of
turbulence in the bubble/shell interface. Figure 12: Temporal evolution of the
structure function of turbulent velocities in the hot gas (Shocked Stellar
Wind and Ionized Gas) in the first $1\,{\rm Myr}$ of the simulation for the
$R_{\rm cloud}=20\,{\rm pc}$, $M_{\rm cloud}=10^{5}M_{\odot}$ simulations. The
separate panels (left to right) are the $\varepsilon_{*}=1\%,\,10\%,\,$and
$100\%$ cases. Temporal evolution is indicated by the color of the lines in
each panel, moving from early times (dark blue) to the late times (bright
green) in increments of $10^{4}\,{\rm yrs}$. The characteristic velocity scale
of the turbulence decreases in time and varies by at most a factor of 2 over
the course of the simulation. The characteristic spatial scale of the
structure function increases over time. Turbulent velocities are overall
higher for higher wind power ($\mathcal{L}_{w}\propto\varepsilon_{*}$) models.
Figure 13: The structure function of turbulent velocities in the hot gas
(Shocked Stellar Wind and Ionized Gas) when $\mathcal{R}_{b}=R_{\rm cloud}/3$.
These are measured in the highest resolution versions of all simulations with
$M_{\rm cloud}=10^{5}M_{\odot}$ and we show the simulations with
$\varepsilon_{*}=1\%,10\%,100\%$ in yellow, orange, and red respectively. For
each panel we show scales from $4\Delta x$ to $R_{\rm cloud}$ (solid lines) as
well as the 25th to 75th percentiles of structure function measurements over
the ensembles of 200 samples from the hot gas. In each panel we show in black
the structure function of the background turbulent gas scaled up by a factor
$20(R_{\rm cloud}/20\mathrm{\,pc})^{1/2}$. The velocity scale of the
turbulence is set by the strength of the wind (represented here by
$\varepsilon_{*}$) while the characteristic length scale of the turbulence is
likely set by the bubble size and turbulent density structure of the
background.

###  4.5 Turbulent Structure Function

Both the amplitude and spectral shape of turbulence are important to the
mixing and cooling that dictate the bubble evolution. These are quantified via
the turbulent structure function.

Following the theory outlined in Paper I, we particularly wish to measure the
turbulence in the hot gas near the interface. To this end, we select gas with
$T>2\times 10^{4}\,{\rm K}$ and $v_{r}<{\cal{V}}_{w}/2$ (the Ionized Gas and
the Shocked Stellar Wind). In order to isolate the turbulent motion from the
bulk radial outflow in this gas, we apply our structure function analysis only
to the non-radial components of the velocity field. We will refer to this
proxy for the turbulent velocity field (which is really the tangential
velocity field) as $\mathbf{u}_{t}$ defined as

$\mathbf{u}_{t}\equiv\mathbf{v}-\left(\mathbf{v}\cdot\mathbf{\hat{r}}\right)\mathbf{\hat{r}}\,.$
(31)

where $\mathbf{v}$ is the full velocity field. Of course, there are also
turbulent contributions to the radial motion, but directly quantifying this is
problematic due to contamination by the strong background radial flow. Figure
11 shows example snapshots of the magnitude of $\mathbf{u}_{t}$ from the
$R_{\rm cloud}=20\mathrm{\,pc}$ model (in the purple to yellow color scheme).

To measure the turbulent structure function, we randomly select a cell $i$
from the region under consideration. We then compute two histograms. The first
histogram is the number of cells at position $j$ with a spatial separation
$\ell=|\mathbf{x}_{j}-\mathbf{x}_{i}|$ from cell $i$, binned by width $\Delta
x$. The second histogram is as above, except now the histogram is weighted by
the square of the difference in the tangential velocity field between
positions $i$ and $j$:

$|\delta\mathbf{u}_{t}|^{2}\equiv\left|\mathbf{u}_{t}(\mathbf{x}_{i})-\mathbf{u}_{t}(\mathbf{x}_{j})\right|^{2}\,.$
(32)

Dividing the second histogram by the first and taking a square root gives us a
single sampling of the root-mean-square velocity offset as a function of
separation scale $\ell$. In order to account for the (unmeasured) radial
component of the true turbulent velocity field, we additionally multiply by
$\sqrt{3/2}$ (assuming that the turbulence is isotropic). We repeat the above
process for 200 cells randomly selected from the Ionized Gas and Shocked
Stellar wind.

With this factor of $\sqrt{3/2}$ we refer to the measured structure function
as $\left\langle\mathbf{v}_{t}^{2}(\ell)\right\rangle^{1/2}$. With these 200
samples of $\left\langle\mathbf{v}_{t}^{2}(\ell)\right\rangle^{1/2}$ we take
the median value in each radial bin (over the ensemble of samples) to
represent the structure function. We also calculate the 25th to 75th
percentiles over these 200 samples to quantify the width of the distribution
of velocity deviations at each scale.

As the calculation described above is quite computationally expensive, we only
calculate the full time evolution of the structure function for a subset of
our simulations ($R_{\rm cloud}=20\,{\rm pc}$, $M_{\rm
cloud}=10^{5}\,M_{\odot}$, and $L_{\rm box}/\Delta x=512$). The results of
these calculations are illustrated in Figure 12, where the turbulent structure
function is shown in time increments of $10^{4}\,{\rm yrs}$ varying from early
times (shown in dark blue) to late times (shown in bright green). The
beginning time is $0.01\,{\rm Myr}$ and the ending times are $2.23,$ $1.21,$
and $0.88\,{\rm Myr}$ for the $\varepsilon_{*}=1\%,$ 10%, and 100% cases
respectively. At all star formation efficiencies the characteristic velocity
scale of the turbulence decreases in time but only by at most a factor of two.
In comparison, the shell velocity decreases by a factor of 10 over the same
time interval.

For the rest of our simulations with $M_{\rm cloud}=10^{5}\,M_{\odot}$ and
$L_{\rm box}/\Delta x=512$, we measure the structure function of the
turbulence at the individual time when $\mathcal{R}_{b}=R_{\rm cloud}/3$. We
show these structure functions, along with the 25th to 75th percentiles of
velocity deviations, in Figure 13. To provide a reference for the background
turbulence, we calculate in the same way the structure function of the
turbulent background in the full simulation volume just before the initiating
of the wind. This background turbulence has a much lower velocity than the
turbulence in the hot gas, so we scale the background measurements by a factor
$20(R_{\rm cloud}/20\mathrm{\,pc})^{1/2}$; this scale choice is purely for
plotting purposes. This scaling also adjusts for the differences in initial
turbulent velocity scales $v_{t}\propto R_{\rm cloud}^{-1/2}$ amongst the
different size clouds (as laid out in Table 1) so that all background
structure functions appear essentially identical.

It is clear from Figure 13 that the velocity scale of the turbulence is
determined primarily by the strength of the wind, which is represented here by
$\varepsilon_{*}$. However, the characteristic spatial scale of the
turbulence, which is where the structure function flattens out, is
consistently $\sim R_{\rm cloud}/10$, suggesting that this scale may be set by
the size of the bubble and/or the spatial scale of the background turbulence.
This variation in turbulent velocity magnitude with $\varepsilon_{*}$ could be
explained by the more oblique shocks in the higher $\varepsilon_{*}$ cases:
when a shock is more oblique the post-shock flow will have a larger speed.

Along those lines, we note that for the $\varepsilon_{*}=100\%$ case in the
$R_{\rm cloud}=2.5\;{\rm pc}$ cloud, the amplitude of the turbulent velocity
is lower than expected. This is because the $v_{r}<{\cal{V}}_{w}/2$
restriction cuts out much of the shocked wind due to the strongly-oblique
shocks in this case (see Figure 3).

Figure 14:  The excess fractal dimension of the bubble surface as a function
of scale using the two measurement techniques described in the text (see
Section 4.6.1 and Section 4.6.2). Left panel: The excess fractal dimension of
the bubble surface measured using a box-counting technique; boxes at different
scales mark where cooling is taking place. Right Panel: The excess fractal
dimension measured based on the logarithmic scaling of the area of an iso-
temperature surface ($T=10^{5}\,{\rm K}$ for all curves shown here) when
measured on different scales.

###  4.6 Fractal Structure of Interface

In this section we explore the fractal nature of the shell-bubble interface.
These results also inform the discussion of turbulence-induced mixing and
cooling in Paper I.

####  4.6.1 Fractal Dimension from Cooling

We measure the fractal dimension of the bubble surface in two different ways.
Since the details of the fractal structure of the bubble/shell interface
affect the cooling, it makes sense to measure the fractal dimension using the
cooling itself. To that end we determine the fractal dimension using a box-
counting or Minkowski-Bouligand algorithm (Schroeder, 1991) where membership
in or out of the fractal is determined based on cooling.222The measurement
technique we adopt was suggested to the authors by Drummond Fielding. To
maximize dynamic range we use simulations with $L_{\rm box}/\Delta x=512$. To
be consistent with the turbulent structure function measurements of Section
4.5, we analyze the fractal structure at the time when $\mathcal{R}_{b}=R_{\rm
cloud}/3$.

We first choose a fraction $f_{\rm cool}\in(0,1)$. For a given time snapshot
we measure the total amount of cooling, $\dot{E}_{\rm cool}$. We denote
$\dot{e}_{{\rm cool},i}$ as the cooling rate in the $i^{\rm th}$
hydrodynamical cell, for a sorted list starting with the maximal cooling rate.
We then compute the cumulative sum up to $I$ such that
$\sum_{i=0}^{I}\dot{e}_{{\rm cool},i}\approx f_{\rm cool}\dot{E}_{\rm cool}$.
The count $I$ of cells contributing to the specified cooling fraction at the
minimum scale $\ell=\Delta x$ (the resolution of the simulation) is hereafter
denoted by $I(\Delta x)$.

We next reduce the size of the simulation grid by a factor of two in each
direction, so that merged cells have length $\ell=2\Delta x$, and the cooling
of adjacent cells is combined. We repeat the above sort-and-accumulate process
to measure the number of cells $I(2\Delta x)$ needed for
$\sum_{i=0}^{I(2\Delta x)}\dot{e}_{{\rm cool},i}\approx f_{\rm
cool}\dot{E}_{\rm cool}$. We iteratively repeat this process for $\ell/\Delta
x=1,\,2,\,4,\,8,\,16,\,$ and $32$ and use these measurements to calculate the
fractal dimension as

$D=d+2=-\frac{d\log I(\ell)}{d\log\ell}\,.$ (33)

This quantity will vary with scale, $\ell$. However, we expect $d$ to take on
a roughly constant value on large scales, dropping below this at sufficiently
small scale due either to limited numerical resolution or to some physical
dissipation scale such as the cooling length, $\ell_{\rm cool}$. In practice
the former is more relevant as the cooling length is usually below our
resolution limit.

At sufficiently large $f_{\rm cool}$, $d$ is close to zero on large scales and
becomes very large (based on our fitting method) on small scales. This is
because $I(f_{\rm cool},\ell)$ becomes a very steep function of $f_{\rm cool}$
when transitioning beyond the concentrated cooling in the interface boundary
layer. Above a certain $f_{\rm cool}$, essentially the whole domain is
included, implying that the true $d\rightarrow 1$ at small $\ell$.

We believe that the most physical $d$ to associate with the fractal structure
of the cooling interface is the value at large scales for the largest value of
$f_{\rm cool}$ that does not exhibit the discontinuous behaviour explained
above. On the left hand side of Figure 14, we therefore show the dependence of
$d$ on scale for the largest choice of $f_{\rm cool}$ that does not result in
this discontinuous behavior. This choice of $f_{\rm cool}$ varies between
simulations.

The scale dependent nature of the “excess” fractal dimension $d$ is evident in
Figure 14, where $d$ decreases at small scales. We see that $d$ tends to be
larger at smaller values of $\varepsilon_{*}$ (this is clearest in the
larger/lower density clouds). In part, this may be because turbulence at the
bubble-shell interface, which is responsible for creating the fractal
structure, makes up a much larger fraction of the kinetic energy at low
$\varepsilon_{*}$, as was shown in Figure 10.

Overall, this measurement method suggests a fractal dimension $d\approx 1/2$
in most simulations, with larger values (up to $\sim 0.6-0.7$) in models with
very low $\varepsilon_{*}$ and large $R_{\rm cloud}$.

####  4.6.2 Fractal Dimension from Iso-Temperature Surfaces

Our second method of measuring the fractal dimension, paralleling Fielding et
al. (2020), uses the area of iso-temperature surfaces for temperatures that
should be characteristic of the interface between the wind bubble and the
shell. Specifically, we measure the area of iso-temperature surfaces for
$T=10^{4}$, $10^{4.5}$, $10^{5}$, $10^{5.5},$ and $10^{6}\,{\rm K}$. We use
the marching_cubes algorithm from the scikit-image package to measure these
surface areas, using step sizes $\ell/\Delta x=1,\,2,\,4,\,8,\,16,\,$and $32$.

The excess dimension $d$ is then given by the logarithmic derivative of the
area with respect to the scale on which the area is measured

$d=-\frac{d\log A_{b}(\ell)}{d\log\ell}\,.$ (34)

As with the previous method, this value should be a constant over a large
range of scales, decreasing at some small scale due to numerical or physical
dissipation (usually the former). Our measurements are made at the time when
$\mathcal{R}_{b}=R_{\rm cloud}/3$. The scale-dependent results are shown at
the right side of Figure 14. For the sake of brevity we only show the results
for the iso-temperature surfaces measured at $T=10^{5}\,{\rm K}$ for each
simulation.

Figure 14 clearly shows that this measurement technique consistently finds
$d\rightarrow 0$ at small scale, as expected. We also found (not shown) that
the value of $d$ on large scales depends on temperature: the higher
temperature surfaces have more fractal (higher $d$) structures than the lower
temperature surfaces, with the effect most extreme in the
smallest-$\varepsilon_{*}$, largest-$R_{\rm cloud}$ models. On large scales,
the range of $d$ from this technique falls roughly within $d\sim 1/2-3/4$.

It is also interesting to note that at all cloud sizes, the bubble surface
tends to become less fractal at higher wind luminosity (or $\varepsilon_{*}$),
as seen in the previous measurement technique, but now to a higher degree.
Physically, the more-fractal structure in the lower-$\varepsilon_{*}$ models
at large scales may potentially be explained by their greater buildup of hot
gas (higher $\mathcal{S}$), which can more effectively create “fingers” of
high-temperature gas within the lower-density parts of the turbulent cloud at
large scales.

Figure 15: The ratio of bubble surface area $A_{b}(T,\ell)$ to fractal area
$4\pi\mathcal{R}_{b}^{2}(\mathcal{R}_{b}/\ell)^{1/2}$ as $\mathcal{R}_{b}(t)$
evolves in time. We only show the evolution up to the point of wind breakout.
We set $\ell=32\Delta x$ and measure the area for several iso-temperature
surfaces as marked in key. Results for $\alpha_{A}$ from taking the area ratio
(see Equation 35) are shown for $M_{\rm cloud}=10^{5}M_{\odot}$, $R_{\rm
cloud}=10\,{\rm pc}$ and $\varepsilon_{*}=1\%,\,10\%,\,$and$\,100\%$. The
near-constancy of $\alpha_{A}$ in time, especially at $\varepsilon_{*}=10\%$
and $100\%$,is an excellent indication that a fractal with $d=1/2$ describes
the geometry of the bubble well.

####  4.6.3 Area Measurements

Using the results of Section 4.6.2, we can now assess the relationship between
the bubble’s surface area and its effective radius, $\mathcal{R}_{b}$, as
defined by Equation 11. In particular, this allows us to test whether a single
fractal dimension can describe the surface, and to evaluate the free parameter
$\alpha_{A}$ in that equation. For this analysis we consider the models with
$R_{\rm cloud}=10\,{\rm pc}$ and $M_{\rm cloud}=10^{5}\,M_{\odot}$ at all
values of $\varepsilon_{*}$.

We set

$\alpha_{A}(\mathcal{R}_{b},\ell)=\frac{A_{b}(T,\ell)}{4\pi\mathcal{R}_{b}^{2}(\mathcal{R}_{b}/\ell)^{1/2}}$
(35)

where $A_{b}(T,\ell)$ denotes the scale dependent area measurement of the iso-
temperature surface (at a range of $T$), and $\ell$ is the measurement scale.
We shall use $\ell=32\Delta x$, the scale on which the fractal structure
begins to saturate (as is evident in Figure 14). Equation 35 compares the
actual area with the prediction for a fractal with excess dimension $d=1/2$ as
this value is broadly consistent across the range of parameters investigated
and measurement techniques used.

We present the time evolution of Equation 35 in Figure 15. Especially for the
larger-$\varepsilon_{*}$ cases, we conclude that the simple $d=1/2$ fractal
agrees extremely well with bubble’s true surface area.

Figure 16: The average value of the dot product of the outward normal vector
to the bubble surface with the radial vector, which we term the “foldedness”
of the bubble surface. The degree to which this is less than one quantifies
the amount that the surface is “inward facing” or “crumpled.” The same models
are used as for Figure 15, along with the same calculations of the bubble iso-
temperature surfaces. We also show $4\pi\mathcal{R}_{b}^{2}/A_{b}(10^{6}\,{\rm
K},32\Delta x)$ in black. The degree to which the black curve matches the
other curves reflects how close the effective area, $A_{\rm eff}$ is to the
equivalent spherical surface $4\pi\mathcal{R}_{b}^{2}$.

####  4.6.4 Bubble Geometry

Finally, we consider an additional characteristic of the bubble geometry,
which we term its “foldedness.” We define this as the value of the dot product
of the unit normal to the bubble surface ($\mathbf{\hat{n}}$) and the unit
radial vector ($\mathbf{\hat{r}}$). There are physical reasons that this
quantity matters. In the spherically symmetric, pressure-driven bubble of
Weaver et al. (1977), El-Badry et al. (2019), and others, the bubble expands
due to a pressure gradient force normal to the bubble surface. Since the
bubble is spherical, the normal vector to the surface is always directed
exactly radially outwards, so that the pressure force and expansion velocity
are exactly aligned. In the case of a bubble in a turbulent medium, the
surface will have a complicated geometry, as Figure 2 and Figure 3 illustrate.
For any pressure force acting normal to the bubble surface, only a fraction of
the force will act radially outward. The radial contribution from a normal
force can be then be quantified by this foldedness quantity
$\left\langle\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}\right\rangle$.

The quantity $\left\langle\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}\right\rangle$
is our measure of the foldedness of the surface in the case of our wind-driven
bubbles, and this is important in setting the radial force that drives the
overall expansion. Specifically, if there is a uniform internal pressure $P$,
the effective outward radial force on the bubble surface can be written as

$F_{r}\approx PA_{\rm eff}(T,\ell)\,,$ (36)

for an appropriate effective area, $A_{\rm eff}(T,\ell)$. We define $A_{\rm
eff}(T,\ell)$ using the area of iso-temperature surfaces $A_{b}(T,\ell)$ from
Section 4.6.2 and the foldedness
$\left\langle\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}\right\rangle$ as

$A_{\rm
eff}(T,\ell)=A_{b}(T,\ell)\left\langle\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}\right\rangle\,.$
(37)

Again, using the marching-cubes algorithm, we evaluate
$\left\langle\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}\right\rangle$ by taking the
dot product of the unit normal to each triangular face output by the algorithm
with the unit radial vector pointing to the center of that triangular face. We
then average these values over the set of faces. The results are displayed in
Figure 16 for the $M_{\rm cloud}=10^{5}M_{\odot}$, $R_{\rm cloud}=10\,{\rm
pc}$ model at all values of $\varepsilon_{*}$. In Figure 16, we also show that
the scaling of
$\left\langle\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}\right\rangle$ in time is
very similar to the temporal behaviour of $4\pi\mathcal{R}_{b}^{2}/A_{b}$.
Both of these quantities decrease in time $\propto
t^{-1/4}\propto\mathcal{R}_{b}^{-1/2}$. But more important than the specific
time dependence is the fact that the scaling of the excess fractal dimension
is compensated by the foldedness of the surface, such that $A_{\rm
eff}(T,\ell)\propto\mathcal{R}_{b}^{2}$.

Our findings on scalings imply that the effective area for radial pressure
forces is the same as for a spherical surface. We posit that this property is
generally true. That is, for any bubble, the quantity $A_{\rm eff}$ defined in
Equation 37 scales as the square of its associated linear scale defined
through the cube root of volume. Though the divergence theorem can be used to
show that $A_{\rm eff}\equiv\int\hat{r}\cdot\hat{n}dA=4\pi\langle R_{\rm
max}^{2}\rangle_{\theta,\phi}$ for $R_{\rm max}$ the bubble radius at a given
spherical polar angle, we have not found a proof to connect this to Equation
4; nevertheless, it seems intuitively reasonable.

The above proposition that $A_{\rm eff}\propto\mathcal{R}_{b}^{2}$ points to
the self-consistency of the assumption in the EC solution that the angle-
averaged radial force in the momentum equation does not explicitly depend on
detailed geometry, even for a fractal bubble.

It is the combination $\rho v_{r}^{2}+P$ that drives outward expansion of the
bubble. The Reynolds stress term is radially directed and therefore always
contributes to the shell momentum in the $\hat{r}$ direction, but the pressure
term would produce a normal force where it acts on the shell, and therefore
would contribute to the radial momentum as in Equation 36. Our conclusion that
$A_{\rm eff}\propto\mathcal{R}_{b}^{2}$ is then what allows us to treat the
total momentum input rate via an equivalent spherical calculation, with the
breakdown of the two contributing terms provided in the Appendix of Paper I
(see Equation A11 there).

Figure 17: The fraction $1-\Theta$ of the input wind energy that is retained
in the bubble and shell, including both kinetic and thermal terms, as a
function of time. We show results for the $R_{\rm cloud}=20\,{\rm pc}$ in the
top row and $R_{\rm cloud}=10\,{\rm pc}$ in the bottom row (both with $M_{\rm
cloud}=10^{5}M_{\odot}$), for $\varepsilon_{*}=1\%,\,10\%,\,$and$\,100\%$. The
time when the wind breaks out of the simulation box is denoted in each panel
by a colored vertical line. The colored solid lines in each panel denote
$1-\Theta$ measured via the energy conservation method, while the dotted lines
denote the same quantity measured using cooling in wind-polluted ($f_{\rm
wind}>10^{-4}$) gas. The shaded areas in each panel delimit the range expected
from EC wind theory (see text). The black line in each panel shows the
theoretical prediction from Equation 10 with values of $f_{\rm turb}$,
$\mathcal{S}$, $\alpha_{p}$, and $\alpha_{R}$ determined as described in
Section 4.4.

###  4.7 Cooling and Energetics

####  4.7.1 Measured Cooling and Retained Energy

Energy inputs from the wind are split between the (thermal and kinetic) energy
in the interior of the bubble, $E_{b}$, the radial kinetic energy of the
shell, $E_{\rm r,sh}$, the turbulent kinetic energy in the shell, $E_{\rm
turb}$, and energy that is lost to cooling, $\dot{E}_{\rm cool}$. The EC
theory presented in Paper I and reviewed in Section 2 provides predictions for
$E_{b}$ and $E_{\rm r,sh}$ and thus (through the total input energy
$\mathcal{L}_{w}t$) for the sum of the energy lost to cooling and the energy
in turbulent motion. The exact split between turbulent energy dissipation and
cooling is not predicted by our theory, but we can account for it in
simulations using the quantity $f_{\rm turb}$, introduced in Section 2 and
shown in Figure 10. In this section we will endeavor to measure the cooling
that is occurring in our simulations directly and check that this is
consistent with the prediction given by Equation 10.

First, we must measure the cooling directly from the simulations. Since,
however, there is cooling present in the background gas, mixing between the
background gas and the wind, and cooling occurring in the background gas that
is shock-accelerated by the wind, this is not a trivial task. The measurement
is further complicated by the fact that the cooling will be a large fraction
of $\mathcal{L}_{w}$, so that an inclusion of only a small amount of
background cooling in our measurement can make our inference of
$\Theta\equiv\dot{E}_{\rm cool}/\mathcal{L}_{w}$ greater than 1. To address
these issues we proceed by measuring the cooling in two separate ways so that
we may check for consistency between our methods.

The first way follows our method for measuring momentum, as laid out in
Section 4.3.1. Specifically, we sum the cooling in all gas with $f_{\rm
wind}>10^{-4}$ in our simulations. The justification of this otherwise
arbitrary cut in $f_{\rm wind}$ is based on the excellent agreement among
measurements of the radial momentum at different resolutions when using this
criterion, as evidenced in Figure 4 and also its agreement with our second
measurement method, detailed below.

We measure the cooling in another way by using conservation of energy. The
time derivative of the total energy in the simulation domain can be written as

$\dot{E}=\mathcal{L}_{w}-\dot{E}_{\rm cool}-\dot{E}_{\rm out}-\dot{E}_{\rm
bkgnd}\,,$ (38)

where $\dot{E}_{\rm cool}$ is the sum of cooling that occurs in the swept up
ambient gas and the bubble/shell interface, $\dot{E}_{\rm out}$ is the rate at
which energy moves out of the simulation domain through the boundaries, and
$\dot{E}_{\rm bkgnd}$ is the net cooling (= cooling - heating) in the ambient
gas that is not due to the wind. While it is straightforward to measure the
total energy in the simulation domain, and hence its time derivative
$\dot{E}$, as well as the rate at which energy leaves the box, it is not so
easy to define what the cooling due to the background ambient gas is.

We address this by running the simulations that we initialized before turning
on a wind (as described in Section 3.2) forward in time without a wind and
measuring the net cooling that occurs throughout the volume. We will denote
the cooling that occurs in these simple turbulent evolution simulations as
$\dot{E}_{\rm tev}$. To account for the fact that part of the simulation
volume that would normally be cooling as background gas will have been
displaced by the wind bubble in our simulations, we measure $\dot{E}_{\rm
bkgnd}$ as

$\dot{E}_{\rm bkgnd}=\left(1-\frac{V_{b}}{L_{\rm box}^{3}}\right)\dot{E}_{\rm
tev}\,.$ (39)

With this measurement of the background cooling we can use Equation 38 to
determine $\dot{E}_{\rm cool}$.

The results of the two methods for calculating $\Theta=\dot{E}_{\rm
cool}/\mathcal{L}_{w}$ are displayed in Figure 17 as solid and dotted lines
for the $M_{\rm cloud}=10^{5}M_{\odot}$, $R_{\rm cloud}=20$ and $10\,{\rm pc}$
cases, for all $\varepsilon_{*}$ values. We only show these lower density
cases as, in the higher density simulations, the background cooling becomes
too strong and difficult to separate from the cooling associated with the
wind, especially at low values of $\varepsilon_{*}$. In fact, this effect is
already evident in the $R_{\rm cloud}=10\,{\rm pc}$, $\varepsilon_{*}=1\%$
case shown in the bottom left panel of Figure 17.

There is excellent agreement between our two measurement methods, with only
moderate deviations at low $\varepsilon_{*}$ values. The dramatic increase in
$1-\Theta$ (decrease in cooling) towards the end of each simulation can be
attributed to breakout and venting of the wind outside of the simulation
domain, the onset of which is indicated by the colored vertical line.

In Figure 17 we also show, as indicated by the shaded region in each panel,
where we expect the solutions for $1-\Theta$ to lie. This is bounded above by
the locus at which the EC condition applies as given by Equation 14, and
bounded below the locus of maximum cooling as given by Equation 10 with
$f_{\rm turb}=0$ and $\mathcal{S}=\alpha_{R}=\alpha_{p}=1$, so that the right-
hand side becomes $(3/2)\dot{R}_{\rm EC}/{\cal{V}}_{w}$. Both of these limits
have retained energy fraction $1-\Theta\propto\dot{R}_{\rm EC}/{\cal{V}}_{w}$
with different coefficients (3.8 and 1.5, respectively), decreasing in time as
$1-\Theta\propto t^{-1/2}$. Finally, as black lines, we show the prediction of
Equation 10 in its entirety, with time-varying values for $f_{\rm turb}$,
$\mathcal{S}$, $\alpha_{p}$, and $\alpha_{R}$ measured as described in Section
4.4.

Overall, Figure 17 shows that the different numerical measurements of
$1-\Theta$ are in good agreement with each other and with theoretical
estimates. At early times ($\sim 10^{4}\mathrm{\,yr}$), the energy retained in
the bubble amounts to $\sim 10\%$ of the input, slightly increasing for lower
density clouds and for more powerful winds (larger $\varepsilon_{*}$). This
decreases $\propto t^{-1/2}$ until the time breakout occurs, $t\sim
10^{5}-10^{6}\mathrm{\,yr}$ (decreasing at higher wind power and for denser
clouds).

Figure 18: Evolution of quantities related to the flow into and through the
turbulent mixing/cooling boundary layer between the bubble and the background.
The expected radial velocity (Equation 15) of shocked wind gas advected into
the fractal interface is shown as dashed lines. The “equivalent velocity” of
gas transiting the mixing/cooling layer (Equation 12 with $d=1/2$) is shown as
solid lines, including a factor of $0.25$ chosen to match the dashed lines.
These curves are shown for the model with $M_{\rm cloud}=10^{5}M_{\odot}$,
$R_{\rm cloud}=20$, and all $\varepsilon_{*}$. We use $\ell_{\rm min}=8\Delta
x$ for the measurement scale of the turbulent structure function.

####  4.7.2 Total Possible Cooling

Finally, we connect our results to the expected requirements for the EC
solution based on the theory of cooling and mixing at turbulent, fractal
interfaces, as discussed in Paper I. To this end, we measure two quantities.
The first is the expected radial velocity of shocked wind gas being advected
to the boundary mixing layer ($v_{\rm hot}(\mathcal{R}_{b})$, given by
Equation 15) The second is the “equivalent velocity” of gas flowing through
the mixing/cooling layer, as given by Equation 12; this takes into account the
fractal structure of the bubble surface, but is only really defined up to a
multiplicative factor from our theoretical analysis.

For this exercise, we use the time dependent measurements for the evolution of
the bubble’s effective radius (Figure 5), the shock radius ($\mathcal{R}_{\rm
f}$), and the turbulent structure function (Figure 12) in models with $M_{\rm
cloud}=10^{5}M_{\odot}$, $R_{\rm cloud}=20$, and all $\varepsilon_{*}$. We
choose $\ell_{\rm min}=8\Delta x$, and calculate $v_{\rm
equiv}=v_{t}(\ell_{\rm min})({\mathcal{R}_{b}/\ell_{\rm min}})^{1/2}$. Note
that we have implicitly adopted a fractal dimension of $d=1/2$, which is
broadly consistent with our results in Section 4.6.1 and Section 4.6.2. As
noted above, the “equivalent velocity” in Equation 12 is only predicted up to
an order-unity coefficient. We therefore multiply $v_{\rm equiv}$ by a factor
$0.25$, and show the results of this calculation in Figure 18.

Figure 18 shows that there is excellent agreement between the expected
velocity at which energy and mass arrives at the surface of the bubble (an
advection speed), and the effective flow velocity through the mixing/cooling
boundary layer (a diffusion speed). This directly demonstrates a key feature
of our model: all thermal energy that is delivered to the turbulent interface
is efficiently mixed in and radiated away as rapidly as it arrives.

##  5 Summary and Conclusion

In Paper I we described a theory for the expansion of a stellar wind-driven
bubble into the dense, turbulent ISM, characterized by strong cooling losses
due to turbulent mixing of the hot gas with denser gas at the bubble surface.
We posited that this cooling was large enough to cause the dominant phase of
the bubble’s evolution to be momentum driven.

A solution in which momentum input (rather than energy input) controls bubble
evolution has been discussed in the past by several authors in various
contexts (Steigman et al., 1975; Ostriker & McKee, 1988; Koo & McKee, 1992a,
b; Kim et al., 2017; Silich & Tenorio-Tagle, 2013), and others have suggested
that such solutions could be associated with efficient mixing and cooling at
boundary layers (Garcia-Segura et al., 1996a; Dale & Bonnell, 2008; Mackey et
al., 2015; Fierlinger et al., 2016), but none have previously demonstrated
that this regime generally applies within star-forming molecular clouds. We do
this by conducting a large suite of 3D hydrodynamic simulations with winds
injected into dense, turbulent ISM material. Analysis of our simulations shows
that the predictions of our theory very accurately describe the evolution of
the wind bubble’s volume, the momentum that it carries, and its energetics.
Our simulations demonstrate that the limit of maximally efficient cooling in
our theory is most appropriate for the strongest stellar winds (parameterized
in our model by high values of the star formation efficiency
$\varepsilon_{*}$).

Our theory and simulations explore in detail where energy is stored and
explain physically how most of it is radiated away, via processes analogous to
those that have been investigated in recent simulations of Kelvin-Helmholtz
unstable mixing layers with fractal geometries (e.g. Fielding et al., 2020;
Tan et al., 2020). The fractal theory of the bubble interface also accurately
predicts essential quantities of the bubble such as the area of its interface
with the ambient ISM and properties of the geometry of this interface. We
additionally demonstrate that it is the fractal nature of the bubble that
allows for such efficient cooling by showing that the cooling capacity under
our fractal theory is greater than the energy being provided by the stellar
wind.

The main conclusions from our simulations can be summarized as follows:

  1. 1.

The effective or mean bubble radius expands in a momentum-driven manner as
$\mathcal{R}_{b}\propto t^{1/2}$ (Equation 5). This is shallower than the
original Weaver et al. (1977) solution or the modified El-Badry et al. (2019)
solution ($\mathcal{R}_{b}\propto t^{3/5}$).

  2. 2.

The total momentum carried by the bubble increases approximately linearly in
time (Equation 3) and is typically only slightly larger (at most a factor 4)
than the total momentum injected by the wind, with the smallest enhancement
for the most luminous winds (largest $\varepsilon_{*}$), as shown in Figure 4.
The momentum enhancement factor $\alpha_{p}$ reflects buildup of hot shocked
gas within the bubble; higher turbulence levels in models with more powerful
winds drive stronger interface mixing and limit this buildup, keeping
$\alpha_{p}$ very close to unity.

  3. 3.

From analysis of the shells in our simulations, we find (see Figure 10) that
the energy carried in tangential motion is typically greater than that carried
in radial motion (turbulence dominated). These results show that shells
typically become more turbulence-dominated in time up to the point of bubble
‘breakout’ from the simulation domain. We also find that cases with less
luminous winds and denser environments are more turbulence dominated.

  4. 4.

We measure the turbulence in the hot gas directly, and find typical velocities
of $150-350\,{\rm km/s}$ (amplitude increasing with input momentum rate; see
Figure 13), which are about $10-15\%$ of the wind velocity ${\cal{V}}_{w}$.
This turbulence has characteristic energy-containing or outer scale that is
larger in larger clouds (Figure 13), suggesting that the background
inhomogeneity induced by turbulence helps to define it. The outer scale is
typically $\sim 10\%$ of the bubble size, growing with the bubble (Figure 12).

  5. 5.

We quantify the excess fractal dimension $d$ of the shell-bubble interface in
two different ways (see Figure 14), and find that this generally falls in the
range of $d=0.4-0.7$. Lower luminosity winds generally have higher fractal
dimensions, which is consistent with the shells driven by these winds being
more turbulence dominated, as it is the turbulence that induces the fractal
structure. We also find that the fractal characterization of the shell
interface is in excellent agreement with the measured area (Figure 15) and
geometry (Figure 16).

  6. 6.

We use $1-\Theta$ to denote the fraction of the energy input rate from the
wind that remains in the bubble as either thermal or kinetic energy at any
given time. We find that $1-\Theta\sim 0.1-0.01$, decreasing in time $\propto
t^{-1/2}$ (Equation 10) for all models (Figure 17). The very small value of
$1-\Theta$ is consistent with observational constraints (see Paper I
discussion).

We have demonstrated that our simulations are well converged in all of the
quantities that we investigate here. Higher resolution should only increase
cooling, so we believe our results for the efficient-cooling bubble evolution
solution are robust. However, our present simulations do not resolve the scale
$\ell_{\rm cool}$ at the wind bubble interface where mixing and cooling
timescales would be comparable in the real ISM. Full validation of the small
scale processes discussed here will therefore require higher resolution
simulations, an important direction for future work.

Finally, we reiterate that the simulations presented here include only a
limited set of the relevant physical processes. In particular, we have not
included radiation or magnetic fields. EUV radiation would photoevaporate gas
from the shell surfaces that “face” the cluster, while FUV radiation absorbed
by dust (where it is not destroyed) directly deposits photon momentum; both
effects drive expansion of the bubble surrounding the cluster. The free wind
and shocked wind regions would then in general interact with photoionized gas
that has been evaporated from the shell, rather than directly with the denser
shell gas (e.g. Dwarkadas & Rosenberg, 2013; Geen et al., 2020a). With a lower
density contrast, Kelvin-Helmholz instabilities at hot/warm interfaces would
have higher growth rates than for hot/cold interfaces at the same pressure.
The turbulent mixing process would likely occupy a larger volume of
photoionized gas compared to the situation here, since the turbulence levels
in the interaction region depend on the density contrast (Fielding et al.,
2020). Magnetic fields presumably also affect the mixing/cooling process.
Given the extremely high velocities of the shocked wind, field strengths in
the cloud would have to be very (unrealistically) high to prevent primary
instabilities at the interface, but magnetic tension could still limit the
turbulent cascade at small scales (and render it anisotropic). Numerical
magnetohydrodynamic studies including both winds and radiation will be needed
to assess how the turbulent mixing/cooling process and overall dynamical
evolution are quantitatively affected.

We thank Drummond Fielding and Erin Kado-Fong for useful discussions. We thank
the referee for their many insightful comments that improved the quality of
the manuscript. This work was partly supported by the National Science
Foundation (AARG award AST-1713949) and NASA (ATP grant No. NNX17AG26G).
J.-G.K. acknowledges support from the Lyman Spitzer, Jr. Postdoctoral
Fellowship at Princeton University. Computational resources were provided by
the Princeton Institute for Computational Science and Engineering (PICSciE)
and the Office of Information Technology’s High Performance Computing Center
at Princeton University.

##  Appendix A Validation of Wind Implementation

We test our implementation of stellar winds by comparing to analytic solutions
given by Weaver et al. (1977) for a spherical wind expanding into a uniform,
static background medium. In our test the background medium has a number
density of $n_{\rm H}=2\,{\rm cm}^{-3}$. The wind has a mechanical luminosity
of $\mathcal{L}_{w}=10^{37}\,{\rm ergs}\,{\rm s}^{-1}$ and a mass loss rate of
$\dot{M}_{w}=10^{-5}M_{\odot}\,{\rm yr}^{-1}$. This corresponds to a wind
velocity of ${\cal{V}}_{w}=1.78\times 10^{3}\,{\rm km}\,{\rm s}^{-1}$. The
simulation was run in a $L=5\,{\rm pc}$ box with $128^{3}$ cells. The feedback
radius was set as $r_{\rm fb}=0.25\,{\rm pc}$ ($r_{\rm fb}/\Delta x=6.4$), and
we used $N_{\rm sc}=4$ subcells. We run simulations both with and without
cooling (the latter is fully adiabatic). The results of these tests are
displayed in Figure 19.

Comparison of the profile for the Weaver et al. (1977) fully adiabatic
solution (described in Section 2 of that paper, shown in green) and our
adiabatic simulation (shown in blue) makes clear that the code accurately
captures the expected behaviour. In the realistic non-adiabatic case the shell
of swept up ambient gas is expected to cool and collapse once the cooling time
in the shell is much shorter than the lifetime of the wind. The profile of the
simulation with cooling (in red) shows that this has occurred, and that the
dense shell is at smaller radius than in the adiabatic case. From the density
slice of the non-adiabatic simulation shown in the bottom right panel of
Figure 19, it is clear that this solution matches the expected evolution of
the non-adiabatic phase of wind evolution (where the expected radius of the
shell in this case, given by Equation 1 of Paper I, is shown as a magenta
circle) and retains a high degree of spherical symmetry.

Figure 19: Wind implementation tests for a uniform background, shown for a
snapshot $10^{4}\,{\rm yrs}$ after the start of wind feedback. We show
profiles (every cell is plotted) of pressure (top left), number density (top
right), and velocity magnitude (bottom left) in radius for both the adiabatic
run (in blue) and the run with cooling (red). The Weaver et al. (1977)
solution for fully adiabatic evolution is shown in green for these same
panels. In the bottom right panel we show a slice of number density through
the $z=0$ plane of the simulation for the run with cooling. One can clearly
see the cooled, dense shell that has formed at the outer edge of the wind
bubble. The Weaver et al. (1977) solution for the radius of the shell in the
case that gas in the leading shock cools is indicated as a magenta circle.

##  Appendix B Additional Tests

In this appendix we display various tests and extensions of our theory.

###  B.1 Pure Thermal Feedback

As described in Section 3, we use a hybrid thermal/kinetic energy injection
scheme which provides for larger time-steps while at the same time avoiding
issues of discontinuity in the injected energy field near the source particle.
However, in a real star cluster, winds initially emanate from a few massive
stars in the form of bulk kinetic energy through line driving (Lucy & Solomon,
1970; Abbott, 1982; Sundqvist et al., 2014). The winds from these separate
stars then collide and shock-heat, converting the initial kinetic energy into
thermal energy within the cluster region (Krause et al., 2013; Cantó et al.,
2000). The concentrated thermal energy is then yet again converted into bulk
kinetic energy as pressure gradients produce acceleration on slightly larger
scale, as in (Chevalier & Clegg, 1985). Since the energy injection zones in
our simulations are comparable to the size of a star cluster, perhaps a more
realistic scheme would be simply to inject purely thermal energy into the
feedback region. This is also interesting to test given that thermal energy is
lost very efficiently in our model (but only at the outer edge of the bubble).

To this end we perform the same simulations described in Appendix B.2 except
adopting the value of $\dot{M}_{w}/M_{*}$ used in the main text, and now
injecting purely thermal energy rather than using the hybrid injection scheme
described in Section 3. The results of these tests are shown in Figure 20 for
the main parameters of interest ($\mathcal{R}_{b}$, $p_{r}$, and $E_{b}$), in
comparison with results for the hybrid injection scheme. The simulations with
thermal energy injection lag those with the hybrid energy injection, but only
very slightly. This is most noticeable for radial momentum, $p_{r}$, while
$\mathcal{R}_{b}$ and $E_{b}$ are both relatively insensitive to this change.
These results indicate that the bubble evolution is not strongly dependent on
the energy injection method.

Figure 20: Tests run with purely thermal energy injection, in comparison to
the hybrid scheme explained in the main text. All panels are analogous to the
top panels shown in Figure 21. Results using the hybrid injection scheme
(thermal scheme) are shown with solid (dashed) lines. The evolution of the
total radial momentum, $p_{r}$ is very slightly delayed for the thermal scheme
compared to the hybrid feedback mechanism, but reaches the same scaling by the
time the bubble radius is $2-3\mathrm{\,pc}$.

###  B.2 Reduced Mass Loss

For practical reasons, our simulations use a value for the wind mass loss rate
per unit stellar mass, $\dot{M}_{w}/M_{*}$, that is about a factor of 5 larger
than that predicted by the SB99 code for a Kroupa IMF. From Figure 1, this
results in a wind momentum injection rate a factor of $\sqrt{5}\approx 2.2$
larger than from SB99. To confirm that our theory still applies when using a
lower value for $\dot{M}_{w}/M_{*}$, we ran tests with
$\dot{M}_{w}/M_{*}=2\times 10^{-3}\,{\rm Myr}^{-1}$, similar to the SB99
level. This correspondingly reduces the momentum injection rate to
$\dot{p}_{w}/M_{*}=7.96\,{\rm km}\,{\rm s}^{-1}\,{\rm Myr}^{-1}$. With a lower
momentum injection rate, the right-hand side of Equation 14 is reduced, and
more cooling (larger $\Theta$) would be required to satisfy the condition for
efficient cooling. If this EC condition is not satisfied, energy will build
within the bubble, so that the EC predictions for radius, momentum, and energy
would underestimate the true values.

For the reduced mass-loss rate, we run three simulations ($\varepsilon_{*}=1$,
$10$, $100\%$) with $R_{\rm cloud}=20\,{\rm pc}$, $M_{\rm
cloud}=10^{5}M_{\odot}$, and $L_{\rm box}/\Delta x=128$. We note that the
choice $R_{\rm cloud}=20\,{\rm pc}$ (the low density regime) is the worst case
scenario for application of the EC theory since cooling is least efficient at
low densities, already making it harder for Equation 14 to be satisfied.

The evolution of the bubble’s effective radius ($\mathcal{R}_{b}$), the total
momentum carried by the bubble ($p_{r}$), and the total bubble internal energy
($E_{b}$) are shown in the top panels of Figure 21. Similar to the figures
shown in Section 4, we display the quantities derived from simulations using
yellow, orange, and red lines compared to the theoretical predictions of the
EC model in black. As expected, the EC theory underestimates
$\mathcal{R}_{b}$, $p_{r}$, and $E_{b}$. However, these differences are small,
especially in the case of the radial evolution.

Moreover, the temporal dependence of each quantity is still very well
described by the EC theory. This is exemplified by the bottom panels of Figure
21, which display quantities related to the ratio of the simulated values to
the values predicted by our EC theory. The fact that these quantities are
roughly constant and order-unity over the course of the simulation, especially
for higher-luminosity clusters, reflects that the EC theory still describes
evolution reasonably well.

Figure 21: Results of our test runs with smaller mass loss rates
($\dot{M}_{w}/M_{*}$). The top panels show comparisons to the EC theoretical
predictions (black) for the bubble’s effective radius (left), total radial
momentum (middle), and total energy in the bubble interior (right) for
simulations run with $\varepsilon_{*}=1\%$ (yellow), $10\%$ (orange), and
$100\%$ (red). The bottom panels show quantities related to the ratio of the
simulated and theoretically predicted values (see Section 4.4).

###  B.3 Changed Feedback Radius

For all the simulations in the main text, we keep the feedback radius, which
determines the region where wind energy is injected, constant. As a further
test of numerical robustness and convergence, we perform a set of simulations
with a feedback radius that is half that given in the main text, $r_{\rm
fb}=0.5\,{\rm pc}$. The set of models is as described in Appendix B.2 but with
the normal $\dot{M}_{w}/M_{*}$ values described in the text and now with
$L_{\rm box}/\Delta x=256$, so that $r_{\rm fb}/\Delta x\gtrsim 2$ (i.e. the
feedback region is still well resolved). We compare to the standard
simulations with $r_{\rm fb}=1\,{\rm pc}$ and $L_{\rm box}/\Delta x=256$.

The results of the comparison between $r_{\rm fb}=1\,{\rm pc}$ and $r_{\rm
fb}=0.5\,{\rm pc}$ cases is shown in Figure 22 for the standard quantities of
relevance to our theory. As expected, the results for the smaller feedback
radius simulations appear to enter the scaling regime earlier than those
simulations with larger $r_{\rm fb}$.

Figure 22: Results of tests with standard ($r_{\rm fb}=1\,{\rm pc}$, solid
curve) and and smaller ($r_{\rm fb}=0.5\,{\rm pc}$, dashed curve) feedback
radius. All panels are analogous to those shown in Figure 20. There is
excellent agreement amongst all quantities, indicating that the changed
feedback radius (higher source resolution) does not impact our results.

###  B.4 Different Turbulent Initial Conditions

When initializing the velocity fields in our simulations, we generally use the
same set of amplitudes and phases (based on the same sequence of random seeds)
to create the turbulent velocity field. With a different initial velocity
field (even for the same initial kinetic energy), the background density
structure into which the bubble expands would be different. To check whether
the results are sensitive to the specific realization of the turbulent
velocity field, we perform the same simulations as in Appendix B.2, using the
standard $\dot{M}_{w}/M_{*}$, except we initialize the velocity field with a
different random seed.

While the cloud will still have the same statistical density structure, the
gas in the immediate vicinity of the star particle will likely be different.
We expect this to result in a slightly different early-time evolution of the
wind bubble, while the same overall evolution would be followed once the
bubble has probed a significant fraction of the cloud size. In reality, the
massive stars that create these high-powered winds will preferentially form at
density maxima, but such a self-consistent treatment is left for later work.

The results of these tests are shown in Figure 23. So as not to over-emphasize
the differences at early times, we display these results on a linear (rather
than logarithmic) time scale. Comparing to the standard simulation with a
different set of turbulent phases, we see that indeed the overall evolution of
the wind-driven bubbles is insensitive to the specific density structure
within the turbulent cloud.

Figure 23: Results of our test runs comparing different initial realizations
of the turbulent velocity field, which creates different cloud density
structure. All panels are analogous to those shown in Figure 20 except here we
use a linear scale for time. This choice is so as to not overemphasize
differences in the early evolution which are due to differences in the density
structure very near to the source particle. We show the original
initialization as solid lines and the alternate initialization as dashed
lines. There is clearly very good agreement at later times between the
evolution, in spite of different density structures. Figure 24: Comparisons of
our theory (black curves) to our simulations (colored curves) for cases with
cloud mass $M_{\rm cloud}=5\times 10^{4}M_{\odot}$ and (from top to bottom
row) $R_{\rm cloud}=20\,{\rm pc}$, $10\,{\rm pc}$, $5\,{\rm pc}$, and
$2.5\,{\rm pc}$. Columns from left to right show the bubble’s effective
radius, $\mathcal{R}_{b}$, the total radial momentum, $p_{r}$, and the total
energy in the bubble interior $E_{b}$. All curves are calculated exactly
analogously to those in Figure 5, Figure 4, and Figure 6. Figure 25: Same as
Figure 24 for cases with $M_{\rm cloud}=5\times 10^{5}M_{\odot}$.

##  Appendix C Results for Different Mass Clouds

In Table 1 we lay out the full range of simulations that we ran. In the main
body of the text we only displayed and discussed results for the cases with
$M_{\rm cloud}=10^{5}\,M_{\odot}$. In this appendix we display results from
the runs with cloud masses of $M_{\rm cloud}=5\times 10^{4}M_{\odot}$ and
$5\times 10^{5}M_{\odot}$, given in Figure 24 and Figure 25, respectively.
Results for both $128^{3}$ and $256^{3}$ resolution are shown.

## References

  * Abbott (1982) Abbott, D. C. 1982, ApJ, 259, 282, doi: 10.1086/160166
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Cantó et al. (2000) Cantó, J., Raga, A. C., & Rodríguez, L. F. 2000, ApJ, 536, 896, doi: 10.1086/308983
  * Castor et al. (1975) Castor, J., McCray, R., & Weaver, R. 1975, ApJ, 200, L107, doi: 10.1086/181908
  * Chevalier & Clegg (1985) Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44, doi: 10.1038/317044a0
  * Dale & Bonnell (2008) Dale, J. E., & Bonnell, I. A. 2008, MNRAS, 391, 2, doi: 10.1111/j.1365-2966.2008.13802.x
  * Dale et al. (2012) Dale, J. E., Ercolano, B., & Bonnell, I. A. 2012, MNRAS, 424, 377, doi: 10.1111/j.1365-2966.2012.21205.x
  * Dale et al. (2013) Dale, J. E., Ngoumou, J., Ercolano, B., & Bonnell, I. A. 2013, MNRAS, 436, 3430, doi: 10.1093/mnras/stt1822
  * Dale et al. (2014) —. 2014, MNRAS, 442, 694, doi: 10.1093/mnras/stu816
  * Dwarkadas & Rosenberg (2013) Dwarkadas, V. V., & Rosenberg, D. L. 2013, High Energy Density Physics, 9, 226, doi: 10.1016/j.hedp.2012.12.003
  * El-Badry et al. (2019) El-Badry, K., Ostriker, E. C., Kim, C.-G., Quataert, E., & Weisz, D. R. 2019, MNRAS, 490, 1961, doi: 10.1093/mnras/stz2773
  * Fielding et al. (2020) Fielding, D. B., Ostriker, E. C., Bryan, G. L., & Jermyn, A. S. 2020, ApJ, 894, L24, doi: 10.3847/2041-8213/ab8d2c
  * Fierlinger et al. (2016) Fierlinger, K. M., Burkert, A., Ntormousi, E., et al. 2016, MNRAS, 456, 710, doi: 10.1093/mnras/stv2699
  * Fukushima et al. (2020) Fukushima, H., Yajima, H., Sugimura, K., et al. 2020, MNRAS, 497, 3830, doi: 10.1093/mnras/staa2062
  * Garcia-Segura et al. (1996a) Garcia-Segura, G., Langer, N., & Mac Low, M. M. 1996a, A&A, 316, 133
  * Garcia-Segura et al. (1996b) Garcia-Segura, G., Mac Low, M. M., & Langer, N. 1996b, A&A, 305, 229
  * Geen et al. (2020a) Geen, S., Bieri, R., Rosdahl, J., & de Koter, A. 2020a, arXiv e-prints, arXiv:2009.08742.  https://arxiv.org/abs/2009.08742
  * Geen et al. (2020b) Geen, S., Pellegrini, E., Bieri, R., & Klessen, R. 2020b, MNRAS, 492, 915, doi: 10.1093/mnras/stz3491
  * Geen et al. (2015) Geen, S., Rosdahl, J., Blaizot, J., Devriendt, J., & Slyz, A. 2015, MNRAS, 448, 3248, doi: 10.1093/mnras/stv251
  * Gronke & Oh (2018) Gronke, M., & Oh, S. P. 2018, MNRAS, 480, L111, doi: 10.1093/mnrasl/sly131
  * Grudić et al. (2020) Grudić, M. Y., Kruijssen, J. M. D., Faucher-Giguère, C.-A., et al. 2020, arXiv e-prints, arXiv:2008.04453.  https://arxiv.org/abs/2008.04453
  * Haid et al. (2018) Haid, S., Walch, S., Seifried, D., et al. 2018, MNRAS, 478, 4799, doi: 10.1093/mnras/sty1315
  * Harper-Clark & Murray (2009) Harper-Clark, E., & Murray, N. 2009, ApJ, 693, 1696, doi: 10.1088/0004-637X/693/2/1696
  * Harris et al. (2020) Harris, C. R., Jarrod Millman, K., van der Walt, S. J., et al. 2020, arXiv e-prints, arXiv:2006.10256.  https://arxiv.org/abs/2006.10256
  * Howard et al. (2017) Howard, C. S., Pudritz, R. E., & Harris, W. E. 2017, MNRAS, 470, 3346, doi: 10.1093/mnras/stx1363
  * Hoyer et al. (2017) Hoyer, S., Hamman, J., Fitzgerald, C., et al. 2017, Pydata/Xarray: V0.10.0, v0.10.0, Zenodo, doi: 10.5281/zenodo.1063607
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jones & Williams (2017) Jones, B. D., & Williams, J. R. 2017, Engineering Computations, 34, 1204, doi: 10.1108/EC-02-2016-0052
  * Kim & Ostriker (2017) Kim, C.-G., & Ostriker, E. C. 2017, ApJ, 846, 133, doi: 10.3847/1538-4357/aa8599
  * Kim & Ostriker (2018) —. 2018, ApJ, 853, 173, doi: 10.3847/1538-4357/aaa5ff
  * Kim et al. (2017) Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, ApJ, 834, 25, doi: 10.3847/1538-4357/834/1/25
  * Kim et al. (2018) Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2018, ApJ, 859, 68, doi: 10.3847/1538-4357/aabe27
  * Kim et al. (2019) —. 2019, ApJ, 883, 102, doi: 10.3847/1538-4357/ab3d3d
  * Kim et al. (2020) Kim, J.-G., Ostriker, E. C., & Filippova, N. 2020, arXiv e-prints, arXiv:2011.07772.  https://arxiv.org/abs/2011.07772
  * Koo & McKee (1992a) Koo, B.-C., & McKee, C. F. 1992a, ApJ, 388, 93, doi: 10.1086/171132
  * Koo & McKee (1992b) —. 1992b, ApJ, 388, 103, doi: 10.1086/171133
  * Koyama & Inutsuka (2002) Koyama, H., & Inutsuka, S.-i. 2002, ApJ, 564, L97, doi: 10.1086/338978
  * Krause et al. (2013) Krause, M., Fierlinger, K., Diehl, R., et al. 2013, A&A, 550, A49, doi: 10.1051/0004-6361/201220060
  * Krause et al. (2016) Krause, M. G. H., Charbonnel, C., Bastian, N., & Diehl, R. 2016, A&A, 587, A53, doi: 10.1051/0004-6361/201526685
  * Krause & Diehl (2014) Krause, M. G. H., & Diehl, R. 2014, ApJ, 794, L21, doi: 10.1088/2041-8205/794/2/L21
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Krumholz et al. (2019) Krumholz, M. R., McKee, C. F., & Bland -Hawthorn, J. 2019, ARA&A, 57, 227, doi: 10.1146/annurev-astro-091918-104430
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Lancaster et al. (2021) Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021, ApJ
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Leroy et al. (2018) Leroy, A. K., Bolatto, A. D., Ostriker, E. C., et al. 2018, ApJ, 869, 126, doi: 10.3847/1538-4357/aaecd1
  * Lucy & Solomon (1970) Lucy, L. B., & Solomon, P. M. 1970, ApJ, 159, 879, doi: 10.1086/150365
  * Mackey et al. (2015) Mackey, J., Gvaramadze, V. V., Mohamed, S., & Langer, N. 2015, A&A, 573, A10, doi: 10.1051/0004-6361/201424716
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Ntormousi et al. (2011) Ntormousi, E., Burkert, A., Fierlinger, K., & Heitsch, F. 2011, ApJ, 731, 13, doi: 10.1088/0004-637X/731/1/13
  * Ostriker & McKee (1988) Ostriker, J. P., & McKee, C. F. 1988, Reviews of Modern Physics, 60, 1, doi: 10.1103/RevModPhys.60.1
  * Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Rahner et al. (2017) Rahner, D., Pellegrini, E. W., Glover, S. C. O., & Klessen, R. S. 2017, MNRAS, 470, 4453, doi: 10.1093/mnras/stx1532
  * Raskutti et al. (2016) Raskutti, S., Ostriker, E. C., & Skinner, M. A. 2016, ApJ, 829, 130, doi: 10.3847/0004-637X/829/2/130
  * Raskutti et al. (2017) —. 2017, ApJ, 850, 112, doi: 10.3847/1538-4357/aa965e
  * Reback et al. (2020) Reback, J., McKinney, W., Jbrockmendel, et al. 2020, pandas-dev/pandas: Pandas 1.0.3, v1.0.3, Zenodo, doi: 10.5281/zenodo.3509134
  * Ressler et al. (2020) Ressler, S. M., Quataert, E., & Stone, J. M. 2020, MNRAS, 492, 3272, doi: 10.1093/mnras/stz3605
  * Roe (1981) Roe, P. L. 1981, Journal of Computational Physics, 43, 357, doi: 10.1016/0021-9991(81)90128-5
  * Rogers & Pittard (2013) Rogers, H., & Pittard, J. M. 2013, MNRAS, 431, 1337, doi: 10.1093/mnras/stt255
  * Schroeder (1991) Schroeder, M. 1991, Fractals, chaos, power laws. Minutes from an infinte paradise
  * Silich & Tenorio-Tagle (2013) Silich, S., & Tenorio-Tagle, G. 2013, ApJ, 765, 43, doi: 10.1088/0004-637X/765/1/43
  * Steigman et al. (1975) Steigman, G., Strittmatter, P. A., & Williams, R. E. 1975, ApJ, 198, 575, doi: 10.1086/153636
  * Stone & Gardiner (2009) Stone, J. M., & Gardiner, T. 2009, New A, 14, 139, doi: 10.1016/j.newast.2008.06.003
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137, doi: 10.1086/588755
  * Sundqvist et al. (2014) Sundqvist, J. O., Puls, J., & Owocki, S. P. 2014, A&A, 568, A59, doi: 10.1051/0004-6361/201423570
  * Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253, doi: 10.1086/191823
  * Tan et al. (2020) Tan, B., Oh, S. P., & Gronke, M. 2020, arXiv e-prints, arXiv:2008.12302.  https://arxiv.org/abs/2008.12302
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Walch et al. (2012) Walch, S. K., Whitworth, A. P., Bisbas, T., Wünsch, R., & Hubber, D. 2012, MNRAS, 427, 625, doi: 10.1111/j.1365-2966.2012.21767.x
  * Wall et al. (2020) Wall, J. E., Mac Low, M.-M., McMillan, S. L. W., et al. 2020, arXiv e-prints, arXiv:2003.09011.  https://arxiv.org/abs/2003.09011
  * Wareing et al. (2017) Wareing, C. J., Pittard, J. M., & Falle, S. A. E. G. 2017, MNRAS, 470, 2283, doi: 10.1093/mnras/stx1417
  * Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377, doi: 10.1086/155692
  * Wünsch et al. (2008) Wünsch, R., Tenorio-Tagle, G., Palouš, J., & Silich, S. 2008, ApJ, 683, 683, doi: 10.1086/589967