# Geometry, Dissipation, Cooling, and the Dynamical Evolution of Wind-Blown
Bubbles

Lachlan Lancaster Department of Astronomy, Columbia University, 550 W 120th
St, New York, NY 10025, USA  Center for Computational Astrophysics, Flatiron
Institute, 162 5th Avenue, New York, NY 10010, USA  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  Institute for Advanced
Study, 1 Einstein Drive, Princeton, NJ 08540, USA  Chang-Goo Kim Department of
Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544,
USA  Jeong-Gyu Kim Division of Science, National Astronomical Observatory of
Japan, Mitaka, Tokyo 181-0015, Japan  Korea Astronomy and Space Science
Institute, Daejeon 34055, Republic Of Korea  Department of Astrophysical
Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA  Greg L.
Bryan Department of Astronomy, Columbia University, 550 W 120th St, New York,
NY 10025, USA  Center for Computational Astrophysics, Flatiron Institute, 162
5th Avenue, New York, NY 10010, USA

(Received XXX; Revised XXX; Accepted XXX)

###### Abstract

Bubbles driven by energy and mass injection from small scales are ubiquitous
in astrophysical fluid systems and essential to feedback across multiple
scales. In particular, O stars in young clusters produce high velocity winds
that create hot bubbles in the surrounding gas. We demonstrate that the
dynamical evolution of these bubbles is critically dependent upon the geometry
of their interfaces with their surroundings and the nature of heat transport
across these interfaces. These factors together determine the amount of energy
that can be lost from the interior through cooling at the interface, which in
turn determines the ability of the bubble to do work on its surroundings. We
further demonstrate that the scales relevant to physical dissipation across
this interface are extremely difficult to resolve in global numerical
simulations of bubbles for parameter values of interest. This means the
dissipation driving evolution of these bubbles in numerical simulations is
often of a numerical nature. We describe the physical and numerical principles
that determine the level of dissipation in these simulations; we use this,
along with a fractal model for the geometry of the interfaces, to explain
differences in convergence behavior between hydrodynamical and magneto-
hydrodynamical simulations presented here. We additionally derive an
expression for momentum as a function of bubble radius expected when the
relevant dissipative scales are resolved and show that it still results in
efficiently-cooled solutions as postulated in previous work.

ISM, Stellar Winds, Star forming regions

††journal: ApJ††software:  scipy (Virtanen et al., 2020), numpy (Harris et
al., 2020), matplotlib (Hunter, 2007), adstex
(https://github.com/yymao/adstex) ††thanks: Simons Fellow

##  1 Introduction

In the modern theory of galaxy evolution and star formation, energetic
feedback from stars (McKee & Ostriker, 2007; Ostriker & Kim, 2022; Walch et
al., 2015; Krumholz et al., 2018; Krumholz, 2014; Sun et al., 2023; Hopkins et
al., 2014; Kruijssen et al., 2019; Chevance et al., 2023) and central black
holes or active galactic nuclei (AGN, Faucher-Giguère & Quataert, 2012;
Zakamska & Greene, 2014; King & Pounds, 2015) is essential in regulating the
flow of mass and energy, and hence the growth of a galaxy and its stellar
population (Somerville & Davé, 2015; Naab & Ostriker, 2017).

Feedback plays multiple roles on multiple scales. AGN and ‘super-bubble’
feedback from repeated supernovae (SNe) drive outflows from galaxies overall
(Tomisaka & Ikeuchi, 1986; Mac Low & McCray, 1988; McCray & Kafatos, 1987; Koo
& McKee, 1992; Faucher-Giguère & Quataert, 2012; Sharma et al., 2014; Krause &
Diehl, 2014; Cicone et al., 2014; King & Pounds, 2015; Yadav et al., 2017; Kim
et al., 2017; El-Badry et al., 2019; Orr et al., 2022; Gentry et al., 2017,
2019), set the energy and mass- and metal- loading of galactic winds
(Schneider et al., 2020; Fielding & Bryan, 2022; Smith et al., 2024; Kim &
Ostriker, 2018; Kim et al., 2020a, b; Pandya et al., 2020), and regulate the
energy balance of the galaxy’s circum-galactic medium (CGM, Fielding et al.,
2017; Faucher-Giguère & Oh, 2023; Pandya et al., 2023). Expanding remnants
from individual and clustered SNe are also believed to bear primary
responsibility for driving turbulence on large scales in the ISM (McKee &
Ostriker, 1977; Elmegreen & Scalo, 2004; Agertz et al., 2009; Ostriker &
Shetty, 2011; Kim et al., 2011, 2013; Hennebelle & Iffrig, 2014; Girichidis et
al., 2016; Ostriker & Kim, 2022). Radiation feedback from stars locally acts
to disperse the molecular gas around star-forming regions (Dale et al., 2005,
2012; Walch et al., 2012; Haid et al., 2018; Grudić et al., 2021, 2022; Kim et
al., 2016, 2018, 2021; Geen et al., 2015, 2016, 2017; He et al., 2019) but
also provides Far Ultraviolet (FUV) radiation which is essential to the
equilibrium thermal balance of the surrounding ISM (Ostriker et al., 2010; Kim
et al., 2023b; Ostriker & Kim, 2022; Kim et al., 2023a). Kinetic feedback from
stellar winds (Lancaster et al., 2021a, b; Fierlinger et al., 2016; Krause et
al., 2013; Harper-Clark & Murray, 2009; Rosen et al., 2014; Lopez et al.,
2011; Rogers & Pittard, 2013; Wareing et al., 2017; Dale et al., 2014) and
proto-stellar jets (Cunningham et al., 2011; Offner et al., 2016; Offner &
Chaban, 2017; Murray et al., 2018) can act to regulate the accretion of mass
onto individual stars (Rosen et al., 2020; Guszejnov et al., 2021), as well as
help to disrupt the parent giant molecular clouds (GMCs) in which these stars
form (Lancaster et al., 2021c; Dale, 2015).

In many of the above situations, high velocity flows shock to create a
‘bubble’ of very hot, very diffuse gas that is ‘blown’ into the surrounding
medium, so-called Wind-Blown Bubbles (WBBs). While in this work we will
generally explore the parameter space of mechanical feedback relevant to WBBs
from stellar winds, much of the material presented here is relevant to
feedback from AGNs, super-bubbles at late stages, or proto-stellar jets with
appropriate changes to the relevant cooling physics and pressures to be
considered. It should be noted that the case of an individual SN remnant (as
well as most of the evolution of a super-bubble) is quite distinct from WBB
dynamics, in that the majority of the work imparted by feedback is done by
individual shocks propagating into surrounding gas, rather than by application
of a continuous pressure force by the hot bubble interior. Because SN feedback
is impulsive, the strict resolution requirements we discuss in this work are
not necessary to accurately model the momentum injection by SNe in the ISM.

It is clear from recent work that the effectiveness of feedback from WBBs is
strongly dependent upon the exchange of energy between the hot gas in a
bubble’s interior and the cold gas in its radiative shell111We have implicitly
specified to a given region of the parameter space relevant to the expansion
of WBBs; we will explain further in Section 2 that this is generally the
regime of greatest interest. (Gentry et al., 2017, 2019; El-Badry et al.,
2019; Lancaster et al., 2021a, b). In particular, such energy exchange,
mediated by either micro-(conduction) or macro-(turbulent) dissipation, can
potentially lead to much of the feedback energy being radiated away.
Additionally, energy exchange across this interface can have drastic
consequences for the mass and density of the hot phase, which in turn has
strong implications for observational signatures of these interfaces (Chu &
Mac Low, 1990; Parkin & Pittard, 2010; Rogers & Pittard, 2014; Krause et al.,
2014; Toalá & Arthur, 2018). Understanding the exchange of energy across such
interfaces is therefore a fundamental problem to the theory of star formation
and galaxy evolution.

As we review in Section 2, correctly resolving this energy exchange should
entail resolving the scales appropriate to the dissipative physics that
mediate this exchange. In some situations, including shear-driven turbulent
boundary layers, numerical simulations indicate that it may not be necessary
to fully resolve the relevant dissipative scales in order for the cooling rate
to be converged (Fielding et al., 2020; Tan et al., 2021).

However, in the present work we will show through detailed analysis of three-
dimensional hydrodynamical (HD) and magneto-hydrodynamical (MHD) simulations
that this property does not hold for the case of simulated WBBs. In
particular, in Section 2 we will expand upon the fluid structure analysis of
Lancaster et al. (2021a) relevant for such bubbles to show exactly how the
geometry of and dissipation at the WBB’s interface determines its dynamical
evolution. We combine this with an extensive discussion of the physics
relevant to setting the geometry and dissipative properties of this interface,
and how they should set the transfer of energy when they are properly resolved
in simulations. We quantify the numerical requirements for the relevant
physical scales to be resolved, and demonstrate that this is difficult to
achieve even for very high numerical resolution. Finally, we present a
consistency condition on cooling in our simulations which can act to determine
the degree of dissipation in the ‘unresolved’ case.

In Section 3 we describe the numerical methods used in our simulations and
present the details of the main suite of simulations that we perform. In
Section 4 we use these simulations to extensively validate the physical
picture described in Section 2 and demonstrate how the resolution dependence
of the bulk properties of the WBBs can be explained through the resolution
dependence of the fractal geometry of the bubbles and dissipation at their
interfaces governed through the consistency condition laid out in Section 2.5.
Finally, in Section 5 we discuss this work in the context of past results both
in regards to WBBs and the mixing layer literature as well as avenues for
future work. We summarize our main conclusions in Section 6.

##  2 Background and Theory

In this section we review the theory of WBBs relevant to the present work and
present several novel interpretations and results. In Section 2.1 we provide a
formal definition of the problem under consideration. In Section 2.2 we review
the results presented in Appendix A, in particular how the fluid structure in
the shocked wind determines the dynamical evolution of the bubble. In Section
2.3 and Section 2.4 we discuss the processes and physical scales relevant for
setting the boundary conditions on the fluid flow interior to the bubble. In
Section 2.5 we consider how a consistency condition that the boundary layer
and bubble interior must mutually satisfy, and how this is affected by the
limited numerical resolution of the simulations we present in Section 4.
Finally, in Section 2.6 we derive the momentum evolution expected for a WBB
where the boundary conditions on its surface are determined by physically
resolved turbulent dissipation.

###  2.1 Problem Definition

We consider the action of a constant mass-loss rate, $\dot{M}_{w}$, and
mechanical wind luminosity, $L_{w}$, source term into a general background
medium. Specifically, the equations of motion describing this problem are:

$\frac{\partial\rho}{\partial
t}+\mathbf{\nabla}\cdot\left(\rho\mathbf{v}\right)=\dot{M}_{w}K(|\mathbf{x}|)\,,$
(1) $\frac{\partial\left(\rho\mathbf{v}\right)}{\partial
t}+\mathbf{\nabla}\cdot\left[\left(P+\frac{B^{2}}{8\pi}\right)\mathbb{I}+\rho\mathbf{v}\otimes\mathbf{v}-\frac{\mathbf{B}\otimes\mathbf{B}}{4\pi}\right]\\\
=\dot{p}_{w}(|\mathbf{x}|)\mathbf{\hat{r}}K(|\mathbf{x}|)\,,$ (2)
$\frac{\partial E}{\partial
t}+\mathbf{\nabla}\cdot\left[\mathbf{v}\left(E+P+\frac{B^{2}}{8\pi}\right)-\frac{\mathbf{B}\left(\mathbf{B}\cdot\mathbf{v}\right)}{4\pi}\right]=\\\
(\mathcal{G}-\mathcal{L})+L_{w}K(|\mathbf{x}|)\,.$ (3)

Equations 1, 2, and 3 are the equations for the conservation of mass,
momentum, and energy respectively. Here, $\rho$ is the mass-density of the
fluid, $\mathbf{v}$ is the fluid velocity, $P$ is its thermal pressure,
$K(|\mathbf{x}|)$ is some spherically symmetric kernel function with compact
support222$K(|\mathbf{x}|)\neq 0$ only in some small region and $\int
K(|\mathbf{x}|)d\mathbf{x}=1$. centered on the origin, $\mathbf{B}$ is the
magnetic field in units with a magnetic permeability of unity, $\mathbb{I}$ is
the identity matrix, $\otimes$ denotes the tensor product, and $E$ is the
total energy density of the fluid given by

$E=e+\frac{1}{2}\rho v^{2}+\frac{B^{2}}{8\pi}\,,$ (4)

where $e$ is the internal energy density, which we will assume to be given by
an ideal gas law $e=P/(\gamma-1)$; we assume $\gamma=5/3$. The source terms
$\dot{M}_{w}K(\mathbf{x})$, $\dot{\mathbf{p}}_{w}(\mathbf{x})K(\mathbf{x})$,
and $L_{w}K(\mathbf{x})$ in Equation 1, Equation 2, and Equation 3,
respectively, denote the input of mass, momentum, and energy from the central
source, which we assume to be constant in time. $\dot{M}_{w}$ and $L_{w}$
together define the wind velocity

$V_{w}=\sqrt{\frac{2L_{w}}{\dot{M}_{w}}}\,.$ (5)

and wind momentum input rate

$\dot{p}_{w}=\dot{M}_{w}V_{w}=2L_{w}/V_{w}=\sqrt{2L_{w}\dot{M}_{w}}$ (6)

We assume that these source terms are “turned on” at some time $t=0$, before
which they are zero and after which they are constant.

Finally, $\mathcal{G}$ and $\mathcal{L}$ in Equation 3 denote thermal heating
and cooling of the gas, respectively, which depend on the local (global)
distribution of the gas temperature and density for optically thin (thick)
radiative cooling and heating. Generally we will assume optically-thin cooling
and heating (Kim et al., 2023b). The above equations are supplemented and made
complete by the ideal induction equation and divergence-free constraint on the
magnetic field.

Under the action of the above source terms a hot ‘bubble’ naturally develops
around the source resulting from the injected mass and energy; this is
commonly referred to as a Wind-Blown Bubble (WBB). At the edge of the bubble
will be an interface with cooler gas. For stellar wind bubbles (and
potentially other applications as well), the bubble will be surrounded by gas
which is photo-ionized. Depending on the regime in which feedback is taking
place this photo-ionized region may or may not be “trapped” by the shell of
the WBB (Peters et al., 2010; Geen & de Koter, 2022; Geen et al., 2023, 2024),
in which case the WBB shell would consist of both photo-ionized ‘warm’ gas and
cool, neutral gas. In the present work, we shall not explicitly include the
effects of photo-ionized gas, but the analysis here may be generalized to
allow for it.

Depending on the exact values of $L_{w}$, $\dot{M}_{w}$, and properties of the
stresses (pressure, density, velocity, magnetic field) in the surrounding
medium, different regimes explaining the expansion of this bubble will be
relevant. Here we focus on the intermediate stage of evolution, which applies
to most high-energy stellar wind sources. This stage has previously been
referred to as the ‘pressure-driven snowplow’ (PDS) stage (Weaver et al.,
1977). However, as argued in Lancaster et al. (2021a, b) (see also El-Badry et
al., 2019), in fact the dynamics differ from the original PDS conception due
to cooling at the interface between the hot gas and the shell. In contrast,
Weaver et al. (1977) assumed that the only energy loss from the bubble’s hot
interior is due to $PdV$ work done in expanding the shell.

What we specifically mean by the intermediate stage is that we assume the
following:

  * 1.

We have left the ‘free-expansion’ stage where the bubble expands ballistically
at $V_{w}$. This occurs very early for all parameter values of interest.

  * 2.

We are dealing with a ‘fast-wind’ in the language of Koo & McKee (1992), where
cooling in the wind bubble interior is negligible.

  * 3.

The exterior stresses outside the bubble and shell are negligible compared to
the pressure in the bubble interior. This implies that the bubble expands
super-sonically (and super-Alfvénically) into the background medium.

  * 4.

The outer shell of the WBB, mostly consisting of gas that was initially shock-
heated when swept up by the bubble, has had enough time to cool back to
ambient temperature.

With these assumptions the bubble naturally has 4 distinct regions, as
described in Weaver et al. (1977): (i) the Free Wind, (ii) the Shocked Wind,
(iii) the shell of the WBB, and (iv) the background medium. The presence of
the shocked wind arises from our first assumption above: the bubble is
expanding slower than $V_{w}$ so the wind must shock. As we will discuss
below, the exact relative volume of (ii) as opposed to (i) will be determined
by the properties of the flow in (ii) and boundary conditions. The structure
of such a bubble, discussed further below, is laid out schematically in Figure
1.

We will define the volume of the free and shocked winds as $V_{\rm bub}$, and
define its volume-equivalent radius as

$R_{\rm bub}\equiv\left(\frac{3V_{\rm bub}}{4\pi}\right)^{1/3}\,.$ (7)

By equating the $dV_{\rm bub}/dt$ derived from a time derivative of the above
with the $dV_{\rm bub}/dt$ derived from considering the motion of the bubble’s
surface it is straightforward to show that

$\frac{dR_{\rm bub}}{dt}=\left\langle W_{\rm bub}\right\rangle\frac{A_{\rm
bub}}{4\pi R_{\rm bub}^{2}}\,.$ (8)

where $\left\langle W_{\rm bub}\right\rangle=(dV_{\rm bub}/dt)/A_{\rm bub}$ is
the expansion velocity of the bubble boundary averaged over the whole surface
(this average is denoted by the brackets $\left\langle\cdot\right\rangle$) and
$A_{\rm bub}$ is the area of the surface of the WBB. This states that the rate
of change of the bubble’s linear dimension is its average outward velocity at
its surface ‘modulated’ by the extent to which the bubble is aspherical,
denoted here by the ratio of the bubble’s surface area to its surface area if
it were a perfect sphere. It was shown in Lancaster et al. (2021b) that, for
fractal bubbles, $4\pi R_{\rm bub}^{2}/A_{\rm
bub}\approx\left\langle\mathbf{\hat{r}}\cdot\mathbf{\hat{n}}\right\rangle$,
the average value of the dot product of the radial and normal vectors over the
surface. This would imply that $\frac{dR_{\rm bub}}{dt}\approx\left\langle
W_{\rm
bub}\right\rangle/\left\langle\mathbf{\hat{r}}\cdot\mathbf{\hat{n}}\right\rangle$,
the factor in the denominator corrects for the reduced surface velocity when
$A_{\rm bub}\gg 4\pi R_{\rm bub}^{2}$.

Figure 1: A schematic representation of the interface of a wind-blown bubble
with the shell consisting of the ambient material it has swept up, in the thin
shell limit. Three aspects of this interface are emphasized. Left: A contact
discontinuity between the ambient material and the wind material, with no
movement of mass or energy across this surface. The hot gas velocity then must
match the velocity of the surface $v_{\rm hot}=\dot{R}_{\rm bub}$. Right: The
inclusion of conductive heat transport across this interface causes mass-
loading of the bubble interior. The flow is usually strongly subsonic
resulting in the hot gas velocity still matching the surface velocity. Middle:
Turbulence in the hot gas, seeded by hydrodynamical instabilities and
inhomogeneities in the background, modifies the geometry of the bubble’s
surface and causes mixing between the hot wind and the surrounding cold gas,
changing the nature of energy transport. This turbulent mixing and subsequent
cooling can be drastic enough to draw in gas much faster than the bubble
itself is expanding, acting as a new boundary condition on the velocity of the
hot gas. The relative volume between the free and shocked wind (the placement
of the reverse shock) can change depending on the boundary conditions.

With the assumptions stated above we have now restricted ourselves to the case
of a bubble with an adiabatic interior (a “fast wind”) and a radiative
exterior (a shell that has cooled). We will see in the next section that the
evolution of such a bubble depends crucially on the physics at the interface
between these two regions.

###  2.2 Bubble Structure Highlights

In Appendix A we expand upon the Appendix of Lancaster et al. (2021a) to
derive the fluid structure in the wind bubble’s interior without assuming a
purely radial flow. In particular we assume only that the structure of the
shocked wind is determined by the properties of a steady, adiabatic, isobaric
flow and discuss the justification for these assumptions and the time-scales
over which they apply.

In the context of a general, aspherical wind bubble and taking into account
the ‘dilation’ of streamlines through a deformed, messy, shocked wind region,
we show that the only relevant stress communicated from the wind bubble
interior to the shell is the thermal pressure force, locally $P_{\rm
hot}dA_{\rm bub}$, where $P_{\rm hot}$ is the (isobaric) pressure of the
shocked wind and $dA_{\rm bub}$ is an area element of the bubble surface.
Locally, this force is in the direction normal to the bubble (or shell)
surface, $\mathbf{\hat{n}}$. We show that the radial pressure force exerted by
the bubble is

$\displaystyle F_{r}$ $\displaystyle\approx$ $\displaystyle 4\pi R_{\rm
bub}^{2}P_{\rm hot}$ (9) $\displaystyle=$
$\displaystyle\frac{3}{4}\dot{p}_{w}\left(\frac{R_{\rm bub}}{R_{\rm
f}}\right)^{2}\,$ (10)

(see Equation A20 and Equation A19, respectively). The second expression shows
that a larger force corresponds to a larger radius of the bubble, $R_{\rm
bub}$, relative to that of the free wind shock, $R_{\rm f}$. As we discuss in
Appendix A, this equation is missing a term due to the Reynolds stress which
is only relevant in the extreme cooling limit where $F_{r}=\dot{p}_{w}$ and
$R_{\rm bub}/R_{\rm f}\approx 1$. We note that the derivation fully accounts
for the ‘foldedness’ of the wind bubble (see Figure 16 of Lancaster et al.,
2021b, and also Equation A17 here). Even though the bubble surface may be
highly folded
($\left\langle\mathbf{\hat{r}}\cdot\mathbf{\hat{n}}\right\rangle\ll 1$ and
$4\pi R_{\rm bub}^{2}\ll A_{\rm bub}$), the net radial force in Equation 9 may
still be written only in terms of the volume equivalent bubble radius $R_{\rm
bub}$.

The interior pressure, as we see from Equation 9, determines the bulk dynamics
of the bubble. As we discuss further in Appendix A, equating the enthalpy flux
through the interior shock, $15L_{w}/16$, with the total enthalpy flux through
the bubble’s surface, $5P_{\rm hot}A_{\rm bub}\left\langle v_{\rm
out}\right\rangle/2$, allows us to write (in the context of a steady bubble)

$P_{\rm hot}=\frac{3}{8}\frac{L_{w}}{\left\langle v_{\rm out}\right\rangle
A_{\rm bub}}\,,$ (11)

where

$\left\langle v_{\rm out}\right\rangle\equiv A_{\rm bub}^{-1}\int_{\rm A_{\rm
bub}}\left(\mathbf{v}-\mathbf{W}\right)\cdot\mathbf{\hat{n}}dA\,.$ (12)

Here, $\mathbf{v}$ is the fluid velocity333When $\dot{R}_{\rm
bub}\ll|\mathbf{v}|$ near the surface, the bubble’s exterior surface is moving
very slowly compared to the velocity at which gas is being advected into the
surface and we may consider $\mathbf{v}-\mathbf{W}\approx\mathbf{v}$ in
Equation 12, and it is dotted into the normal to the surface
$\mathbf{\hat{n}}$. The above assumes a clear definition of the bubble’s
surface, which requires additional care in the case that the surface has a
fractal structure. In all physical cases, this fractal structure will be
regulated or smoothed-out at some small scale determined by dissipative
physics. However, as we discuss in Appendix B, one can give a meaningful
definition to both $A_{\rm bub}$ and $\left\langle v_{\rm out}\right\rangle$
measured on scales larger than this dissipative scale, provided that both are
measured using the same ‘ruler.’444We shall analyze measurement scale
dependencies of $A_{\rm bub}$ from our simulations in Section 4.4. Our
measurements of $\left\langle v_{\rm out}\right\rangle$ will be at the
smallest scale available in the simulation, $\Delta x$. In Appendix A we
discuss how, if Equation 11 is not satisfied at a given time, the energy
dynamics of the bubble interior will act to restore it by losing or gaining
energy, causing a change in pressure.

Using Equation 9 and Equation 11 we can write a momentum equation for the
evolution of the shell surrounding the bubble as

$\dot{p}_{r}=\frac{3L_{w}}{8}\frac{4\pi R_{\rm bub}^{2}}{\left\langle v_{\rm
out}\right\rangle A_{\rm bub}}\,.$ (13)

The momentum equation of Lancaster et al. (2021a) for the shell surrounding a
bubble (Equation 18 of that paper with $\alpha_{R}=1$) is

$\frac{dp_{r}}{dt}=\frac{d}{dt}\left(\frac{4\pi}{3}\bar{\rho}R_{\rm
bub}^{3}\frac{dR_{\rm bub}}{dt}\right)=\alpha_{p}\dot{p}_{w},$ (14)

where $\bar{\rho}$ is the mean density in the environment555Equation 14 has
solution $R_{\rm
bub}=\left[3\alpha_{p}\dot{p}_{w}t^{2}/(2\pi\bar{\rho})\right]^{1/2}$, and the
limit of maximally efficient cooling (‘EC’) takes $\alpha_{p}\rightarrow 1$..
Comparing these and using $\dot{p}_{w}=2L_{w}/V_{w}$, we can identify the
momentum enhancement factor in terms of the bubble parameters as

$\alpha_{p}=\frac{3}{4}\frac{V_{w}/4}{\left\langle v_{\rm
out}\right\rangle}\frac{4\pi R_{\rm bub}^{2}}{A_{\rm bub}}\,.$ (15)

We can see from the above that whether or not a bubble appears to be
“momentum-driven” (i.e. $\alpha_{p}\approx{\rm const.}$) depends solely on how
$4\pi R_{\rm bub}^{2}/A_{\rm bub}$ and $\left\langle v_{\rm
out}\right\rangle/V_{w}$ scale in time with respect to one another; these
factors are driven by geometry666The $4\pi R_{\rm bub}^{2}/A_{\rm bub}$
dependence in Equation 15 can be thought of straightforwardly as the result of
enthalpy flux being dependent on the total area (increase in $A_{\rm bub}$
results in further losses) whereas the radial momentum is dependent only on
the radially outward component. and dissipative processes respectively. We
will next discuss the physics that govern each of these quantities.

###  2.3 The Bubble-Shell Interface

In Figure 1 we have laid out three different ways of thinking about the
boundary between the bubble’s interior and its shell: a contact discontinuity
(ideal fluid dynamics), a turbulent, radiative mixing layer (TRML) which
simply emphasizes the importance of hydrodynamical instabilities in modifying
the geometry of the interface and enhancing micro-physical dissipation, or a
conductive evaporation or condensation front (non-ideal fluid dynamics, with
cooling). In reality, the boundary between the wind-blown bubble and its shell
may combine all of these features. Ignoring for a moment the diffusion of
mass, there is a ‘contact surface’ that separates the injected mass from the
surrounding mass. This surface is distorted and amplified in area by dynamical
instabilities as well as the larger scale inhomogeneities in the background
medium. Across this distorted surface, heat dissipation acts to carry energy
from the hot wind to the cold shell. As we can see from Equation 15, the
details of the physics of this interface, through its geometry (given by $4\pi
R_{\rm bub}^{2}/A_{\rm bub}$) and the rate at which energy is diffused across
the surface (which determines $\left\langle v_{\rm out}\right\rangle$) have a
profound impact on the dynamics of the bubble.

We will show in Section 4.4 that the geometry of these interfaces is well
described by a fractal model. In the next section we discuss in depth the
various diffusive physics that can set $\left\langle v_{\rm
out}\right\rangle$.

Figure 2: Dissipative scales as a function of temperature (left) at
$P/k_{B}=10^{5}\,{\rm K}\,{\rm cm}^{-3}$ and pressure (right) using the values
at $T_{\rm pk}$. Left panel: $\ell_{\rm cool}$ from turbulence-mediated mixing
as defined in Equation 24 with $L=4\,{\rm pc}$ and $v_{t}(L)=300\,{\rm km/s}$;
$\lambda_{F}$ as defined in Equation 18 with temperature dependent $\kappa$
(see Equation 19). The vertical gray dashed lin indicates the temperature
which yields the shortest cooling time $\approx 1.6\times 10^{4}\,{\rm K}$.
Right panel: The diffusive scales at the temperatures responsible for mixing
as a function of pressure. Each length scale is displayed relative to the
bubble radius for the ‘efficiently cooled’ solution (inverse of Equation A4
with $R_{\rm f}\to R_{\rm bub,EC}$ and $\dot{p}_{w}$ corresponding to our
fiducial simulations). The turbulent diffusive scale $\ell_{\rm cool}$ is
shown with $v_{t}(L=4\,{\rm pc})=75,\,150,\,300,\,\&\ 600\,{\rm km/s}$ which
correspond to the lighter to darker blue lines, respectively.
$\lambda_{F}(T_{\rm pk})$ is shown in brown (see Equation 18) and the
resolution of our highest resolution simulations is shown in black. The rough
range in pressures spanned by our simulations are shown as a grey band.

###  2.4 Mechanisms and Physical Scales for Dissipation

The setting of $\left\langle v_{\rm out}\right\rangle$ on the smallest scales
has to be regulated by balancing the transport of heat through dissipative
processes and the rate at which energy can be cooled away. The latter is
straightforwardly given by a characteristic cooling time

$t_{\rm
cool}(T;P)\equiv\frac{P}{n^{2}\Lambda(T)}=\frac{(k_{B}T)^{2}}{P\Lambda(T)}\,,$
(16)

where $\Lambda(T)$ is the cooling function777In the simulations presented
here, and in reality, $\Lambda$ is strictly speaking not only a function of
$T$ but also of, primarily, the electron fraction $x_{e}$. However, we ignore
that dependence here since the phases we care most about are fully ionized.,
$n$ is the number density of the gas, $k_{B}$ is Boltzmann’s constant, and $T$
is the temperature. Generally we will consider isobaric layers (though this
does not hold strictly in our simulations) so $t_{\rm cool}(T;P)$ may be
regarded solely as a function of temperature (at given pressure $P$). The
minimum cooling time as a function of pressure determines the temperature of
gas the cools most rapdily, which we will refer to as $T_{\rm pk}$

$T_{\rm pk}\equiv{{\rm arg}\,{\rm min}}_{T}t_{\rm cool}(T)\,.$ (17)

Since the dependence of $t_{\rm cool}$ on pressure does not change the
temperature dependence, $T_{\rm pk}$ is independent of pressure and is equal
to $1.59\times 10^{4}\,{\rm K}$ for our adopted cooling function (which
assumes solar metallicity).

We next discuss two forms of physical dissipation and their relevant spatial
and velocity scales. Where appropriate, we discuss how magnetic fields alter
these means of dissipation. The dissipative scales defined below are functions
of the pressure of the medium, $P$, and the temperature, $T$, which may vary
throughout the medium, although the bubble interior and boundary layer are
roughly isobaric. In Figure 2 we show these scales as functions of both
temperature (left panel) and pressure (right panel) for given choices of the
parameters of the dissipative physics, defined and explained below.

####  2.4.1 Thermal Conduction

The first dissipation mechanism is thermal conduction mediated by individual
particle collisions over a mean free path (Spitzer & Härm, 1953; Spitzer,
1962), parameterized by the (temperature-dependent) thermal conductivity
$\kappa(T)$, which produces a heat flux $-\kappa\nabla T$. In balancing heat
transport through conduction against cooling in the energy equation, one can
derive a characteristic length scale $\lambda_{F}$, which we will refer to as
the ‘Field Length’ 888Note that this differs by order unity from the length at
which thermal instability is suppressed by thermal conduction (Field,
1965).(Field, 1965; Zel’Dovich & Pikel’Ner, 1969; Kim & Kim, 2013; Tan et al.,
2021):

$\lambda_{F}(T)=\sqrt{\frac{\kappa(T)T}{n^{2}\Lambda(T)}}=\sqrt{\frac{\kappa(T)Tt_{\rm
cool}(T)}{P}}\,.$ (18)

Below we take $\kappa$ as given by either

$\frac{\kappa_{SP}(T)}{{\rm erg}/{\rm s}/{\rm cm}/{\rm
K}}=\begin{cases}2.5\times 10^{3}\left(\frac{T}{\rm K}\right)^{1/2}&T<T_{\rm
eq}\\\ 6\times 10^{-7}\left(\frac{T}{\rm K}\right)^{5/2}&T\geq T_{\rm
eq}\end{cases}$ (19)

which is a combination of Spitzer (high temperatures) and Parker (low
temperatures, Parker, 1953) conductivity, where the point of switching between
the two $T_{\rm eq}=6.45\times 10^{4}\,{\rm K}$ is the point at which they are
equal.

In the left panel of Figure 2, we show $\lambda_{F}$ as a function of
temperature $T$ at $P/k_{B}=10^{5}{\rm K\ cm^{-3}}$. The cooling function used
here, $\Lambda(T)$, is that from Kim & Ostriker (2017). Since $t_{\rm
cool}\propto P^{-1}$, it is clear that $\lambda_{F}\propto P^{-1}$ if one
considers a fixed temperature. In the right panel of Figure 2, we show
$\lambda_{F}$ relative to the bubble radius $R_{\rm bub}$, selecting
$\lambda_{F}(T_{\rm pk})$. The bubble radius is calculated as a function of
$P_{\rm hot}$ assuming the efficiently cooled solution (Equation A4 with
$R_{\rm f}\to R_{\rm bub,EC}$) with $\dot{p}_{w}=4.8\times
10^{4}M_{\odot}/{\rm km}/{\rm s}/{\rm Myr}$. Really, we should choose the
value at which $\lambda_{F}$ is minimized in the range of temperatures probed
by the conductive interface, roughly $5000\,{\rm K}<T<10^{8}\,{\rm K}$, which
is slightly different from $T_{\rm pk}$ due to the temperature dependence of
$\kappa(T)T$ in the numerator of Equation 18. However, the dependence of
$\lambda_{F}$ on $T$ is so strongly influenced by $\Lambda(T)$ that this
temperature is essentially $T_{\rm pk}$, so we will use the latter for
simplicity.

Note that, in the presence of a magnetic field, $\mathbf{B}$, thermal
conduction across field lines is expected to be suppressed, due to suppression
of electron transport across field lines. In the context of WBBs expanding
into a magnetized medium, the process of ‘magnetic draping’ of the field lines
around the bubble will virtually guarantee that the magnetic field is mostly
parallel to the bubble’s surface, suppressing conduction. Furthermore, recent
plasma physics simulations and experiments have demonstrated that thermal
conduction can be strongly suppressed due to plasma instabilities in
high-$\beta$, weakly collisional and turbulent plasmas, as may be the case in
some regimes of these mixing layers (Roberg-Clark et al., 2016; Komarov et
al., 2018; Drake et al., 2021; Meinecke et al., 2022). Additionally, as
derived in Cowie & McKee (1977), there is a maximum heat flux that can be
carried by electrons, $q_{\rm max}\sim(3/2)\rho c_{s,{\rm iso}}^{3}$ (where
$c_{s,{\rm iso}}$ is the isothermal sound speed), which limits the flux when
the temperature gradient scale length, $T/|\nabla T|$, is comparable to the
mean free path of electrons. All of these effects can decrease $\kappa$ by an
order of magnitude, making $\lambda_{F}$ a factor of 3 or more smaller.

We can derive a rough estimate of the value of the velocity that would be
induced by resolved thermal conduction by assuming that all cooling occurs in
a layer of thickness $\lambda_{F}(T_{\rm pk})$. We assume the cooling proceeds
at a rate per unit area of $\lambda_{F}n^{2}\Lambda(T_{\rm pk})$ and balance
this against advection of enthalpy from the hot gas into the layer at a rate
$5Pv_{\rm cond}/2$. This gives us

$v_{\rm cond}\approx\frac{2}{5}\sqrt{\frac{\kappa(T_{\rm pk})\Lambda(T_{\rm
pk})}{k_{B}^{2}T_{\rm pk}}}\approx\frac{\lambda_{F}}{t_{\rm cool}}\,,$ (20)

which is independent of pressure.

We should note that there is another scale associated with the action of
thermal conduction, which comes from the balancing of the advection term with
the conductive heat transport term in a one-dimensional energy equation of a
conductive cooling layer (e.g. Tan et al., 2021; Tan & Oh, 2021). This
‘evaporation length’ scale is given by

$\lambda_{\rm evap}=\frac{2\kappa T}{5Pv}$ (21)

where $v$ is the gas velocity. Though the value of this length scale depends
on the details of the evaporative mass loss (which is subject to the energy
lost in the cooling layer – see section 5.3 of El-Badry et al., 2019), this
scale is generally much larger than $\lambda_{F}$. It is therefore generally
possible to resolve an evaporative mass flow from the surface even if the
details of cooling in the surface are unresolved in a given simulation.

####  2.4.2 Turbulent Dissipation

The second dissipation mechanism is turbulent diffusivity, which has been
explored extensively in recent literature (Fielding et al., 2020; Tan et al.,
2021; Lancaster et al., 2021a). In particular, we imagine a turbulent cascade
in the hot gas with energy containing scale $L$ and velocity normalization
$v_{t}(L)$, following a spectrum

$v_{t}(\ell)=v_{t}(L)\left(\frac{\ell}{L}\right)^{p}\,,$ (22)

with $p=1/3$ for subsonic, Kolmogorov turbulence. We define the eddy turn-over
time at scale $\ell$ as

$t_{\rm eddy}(\ell)\equiv\frac{\ell}{v_{t}(\ell)}\,.$ (23)

The mixing that leads to cooling is expected to take place on a scale
$\ell_{\rm cool}$ such that $t_{\rm eddy}(\ell_{\rm cool})=t_{\rm cool}$,
where Equation 16 gives $t_{\rm cool}$. Combining Equation 22 and Equation 23
with this definition we obtain the cooling length

$\ell_{\rm cool}\equiv L\left(\frac{t_{\rm cool}}{t_{\rm
eddy}(L)}\right)^{\frac{1}{1-p}}\,.$ (24)

In this picture, eddies on scales larger than $\ell_{\rm cool}$ mix too slowly
to contribute significantly to energy losses. At scales $\lesssim\ell_{\rm
cool}$, the surface area of the interface is rapidly enhanced on a time-scale
shorter than $t_{\rm cool}$ (since $t_{\rm eddy}(\ell)<t_{\rm cool}$ for all
$\ell<\ell_{\rm cool}$). In particular, this enhancement in area is so rapid
that it does not matter what the micro-physical means of energy transport
across the surface is; the area will always become large enough to carry the
heat being mixed by the turbulence. In the language of the combustion
literature (Tan et al., 2021), we can think of mixing on the scale $\ell_{\rm
cool}$ as occurring in the ‘well-stirred reactor’ limit, at which point an
effective turbulent diffusivity may be applied. Either way, this results in
the $\ell_{\rm cool}$-scale eddy being ‘mixed’ to the intermediate
temperatures where cooling is rapid.

The velocity scale associated with turbulent dissipation is built in:

$v_{\rm turb}\approx v_{t}(\ell_{\rm cool})=v_{t}(L)\left(\frac{t_{\rm
cool}}{t_{\rm eddy}(L)}\right)^{\frac{p}{1-p}}\,.$ (25)

We can see that this is very similar to Equation 20 if we remember that
$v_{t}(\ell_{\rm cool})=\ell_{\rm cool}/t_{\rm cool}$.

In the presence of a magnetic field, turbulent heat diffusion is also altered
due to the modification of instabilities which drive the turbulence (Das &
Gronke, 2023; Zhao & Bai, 2023). Roughly speaking, instabilities will be
suppressed if the velocity scale associated with the instability (by which we
mean the length scale associated with the instability, divided by its linear
theory growth time) is less than the Alfvèn speed in the medium
(Chandrasekhar, 1961). This likely has drastic consequences for the WBB as,
due to compression of the magnetic field, the shell of the WBB will be
strongly magnetically dominated. This is likely the cause of the large
differences in fractal structure (Section 4.4) and turbulent properties
(Section 4.7) between our HD and MHD simulations.

In this section and the previous one we have defined two mechanisms of
physical heat dissipation and the associated length scales relevant to
resolving them. These length scales are functions of temperature, as shown in
the left panel of Figure 2. Considering one dissipation mechanism alone, it is
relevant to resolve the minimum such length scale in the relevant temperature
range for the problem. This is essentially determined by where the cooling
rate is maximized, due to the strong temperature dependence of $\Lambda(T)$.
However, considering a realistic case where both dissipation mechanisms are
acting, it is only important to resolve the larger of the two length scales,
as the associated dissipation mechanism will dominate the heat transport. For
example, in the right panel of Figure 2 we can see that, for
$v_{t}(L)=300\,{\rm km/s}$, the turbulent dissipation scale, $\ell_{\rm
cool}$, is always greater than length scale associated with thermal
conduction, $\lambda_{F}$. In this case, it is only necessary to resolve
$\ell_{\rm cool}$, which is can be orders of magnitude larger than
$\lambda_{F}$ at lower pressures (later in the bubble evolution). If instead
we were simulating a bubble where $\ell_{\rm cool}<\lambda_{F}$, as is the
case for the $v_{t}(L)=75\,{\rm km/s}$ case in Figure 2, it would then be
necessary to resolve $\lambda_{F}$, though we will see that this is often not
the case.

####  2.4.3 Numerical Dissipation

When the appropriate dissipative scales are resolved (the maximum of
$\ell_{\rm cool}$ or $\lambda_{F}$ for a given condition of the gas) physical
dissipation controls the pace at which energy can be transported across the
interface and therefore sets $\left\langle v_{\rm out}\right\rangle$. As we
can see from comparing $\Delta x\approx 0.1\,{\rm pc}$ at our highest
resolution to the relevant dissipative scales in the right panel of Figure 2,
we are only just near resolving the turbulent dissipation at the latest times
in our simulations for relevant choices of the turbulent dissipative
parameters (see Section 4.7). The question then arises: how is $\left\langle
v_{\rm out}\right\rangle$ set when the relevant physical dissipative scales
are not properly resolved?

In keeping with the previous two subsections, let us first suppose that
$\left\langle v_{\rm out}\right\rangle$ is set by dissipation, this time of a
numerical nature. In numerical dissipation the relevant dissipative velocity
scale associated with dissipation across contact discontinuities (as is
relevant here) is roughly the local gas velocity relative to the grid. For
example, consider the advection equation in one-dimension,
$\partial_{t}q+a\partial_{x}q=0$ with $q$ some scalar field and $a$ the
advection speed. If we solve this equation on a fixed grid with spatial
discretization $\Delta x$, using a first-order upwind scheme one can
analytically show that the actual equation being solved (to second order
accuracy) is $\partial_{t}q+a\partial_{x}q=C_{\rm num}\partial_{x}^{2}q$ with
numerical diffusivity $C_{\rm num}=\Delta xa(1-|c|)/2$ with $c$ the Courant
number (see e.g. Chapter 5 of Toro (2009) or Chapter 20 of Press et al.
(2002)). From $C_{\rm num}$ we can read off the relevant spatial and velocity
scale for the dissipation: $\Delta x$ and $a$, the resolution and the
advection velocity. While it is much more difficult to derive an exact
numerical diffusivity for a higher-order method, let alone a system as
complicated as the MHD equations, the same basic principle still holds. That
is, for diffusion across contact discontinuities, the numerical diffusive
velocity scale, $v_{\rm num}$, is comparable to the local velocity of the gas
relative to the grid:

$v_{\rm num}\approx v_{\rm gas}\,,$ (26)

in analogy to Equation 20 and Equation 25.

In the previous sections we posited that the local gas velocity is set by the
relevant diffusive velocity at the interface $v_{\rm gas}\approx v_{\rm diff}$
with $v_{\rm diff}=v_{\rm cond},\,v_{\rm turb}$. However, here, this provides
no constraint, as the relevant diffusive velocity is the local gas velocity.
We will see in the next section that a consistency condition on cooling at the
wind bubble’s interface will lead to a meaningful estimate of the numerical
diffusive velocity as

$v_{\rm num}\approx\frac{\Delta x}{t_{\rm cool}}\,.$ (27)

###  2.5 Relation Between $\left\langle v_{\rm out}\right\rangle$ and Cooling

As we discuss in Appendix A, the assumption of an entirely ‘steady’ flow in
the bubble interior necessitates that most of the mass and energy in the
bubble is advected out of it. As discussed in Section 2.4, this transfer can
be accomplished by the diffusive flux associated with a turbulent boundary
layer in the physical case; numerical diffusion can alternatively play this
role if the turbulence is not sufficiently resolved. Whether physically or
numerically driven, mixing of the energy from the interior with mass from the
shell can produce temperatures where cooling is very strong, so the energy
mixed into the shell can be largely radiated away. This loss of energy is
essentially a consistency condition on the assumption of a ‘steady’ interior:
if energy were not lost, there would be no movement of energy out of the hot
gas bubble, as is required by the steady flow assumption.

If we assume for simplicity that this cooling is dominated by gas at $T_{\rm
pk}$, given in Equation 17, and that the volume of gas at this temperature is
$V_{\rm pk}$ then the cooling rate of this gas is

$\dot{E}_{\rm cool,pk}=\left(\frac{P(T_{\rm pk})}{k_{B}T_{\rm
pk}}\right)^{2}\Lambda(T_{\rm pk})V_{\rm pk}=\frac{P(T_{\rm pk})V_{\rm
pk}}{t_{\rm cool}(T_{\rm pk})}$ (28)

where $\Lambda(T)$ is the cooling function (we have ignored heating here),
$P(T_{\rm pk})$ is the pressure of the gas in $V_{\rm pk}$, and $t_{\rm cool}$
is the cooling time (Equation 16). If we assume that the cooling layer is in
pressure equilibrium with the hot gas, as seems reasonable from resolved
simulations (El-Badry et al., 2019; Fielding et al., 2020), then we have
$P(T_{\rm pk})=P_{\rm hot}$.

We can define the thickness of the cooling region as

$w_{\rm cool}\equiv\frac{V_{\rm pk}}{A_{\rm bub}}\,.$ (29)

Then, with $\dot{E}_{\rm cool,pk}\approx L_{w}$ from the assumption of strong
cooling, and using Equation 11 and Equation 28, we have

$\left\langle v_{\rm
out}\right\rangle\approx\frac{3}{8}\frac{\sqrt{\Lambda(T_{\rm pk})V_{\rm
pk}L_{w}}}{A_{\rm bub}k_{B}T_{\rm pk}}\propto\left(\frac{w_{\rm cool}}{A_{\rm
bub}}\right)^{1/2}\,.$ (30)

The last proportionality has ignored factors having to do with the exact
properties of the cooling (as well as the source luminosity), assuming them to
be roughly constant.

It is worth discussing what the above relation is telling us about both the
bubble and boundary layer dynamics and thermodynamics. Suppose that, at a
given pressure, $\left\langle v_{\rm out}\right\rangle$ were too high for
Equation 30 to hold. This incongruence can be seen as a failure from either
the assumption of cooling matching the energy flux being provided
($\dot{E}_{\rm cool,pk}\approx L_{w}$), or in the context of a steady state
bubble (Equation 11). From the perspective of the bubble, $\left\langle v_{\rm
out}\right\rangle$ being too large means that the bubble is losing energy at
its outer surface faster than it gains energy through the wind shock, which
will lead to pressure loss. If we take $\left\langle v_{\rm out}\right\rangle$
as given, steady-state can be re-established when de-pressurization reaches
the point at which Equation 11 is again satisfied (which would also involve an
increase in $R_{\rm f}$ from Equation A4). In our picture, this decrease in
pressure also affects the rate of cooling in the interface layer, with the new
lower pressure driving down cooling from Equation 28. From the perspective of
the cooling layer, $\left\langle v_{\rm out}\right\rangle$ being too large
means that it is being provided energy faster than it can cool it away. If we
enforce, in this thought experiment, fixed $\left\langle v_{\rm
out}\right\rangle$ and isobaric conditions, the build up of energy resulting
from insufficient cooling would cause an expansion of the cooling layer (an
increase in $w_{\rm cool}$). This would continue until the cooling layer was
once again in equilibrium according to Equation 30.

We emphasize that in the resolved case, Equation 30 is not a prediction of the
value of $\left\langle v_{\rm out}\right\rangle$ but rather a consistency
condition on the relationship between $\left\langle v_{\rm out}\right\rangle$,
the properties of the cooling, and the geometry of the cooling layer ($w_{\rm
cool}$ & $A_{\rm bub}$). In fact, if we replace $P(T_{\rm pk})$ in Equation 28
with Equation 11 and assume $\dot{E}_{\rm cool,pk}=L_{w}$ as above we have

$\left\langle v_{\rm out}\right\rangle=\frac{3}{8}\frac{w_{\rm cool}}{t_{\rm
cool}}\sim\frac{w_{\rm cool}}{t_{\rm cool}}\,.$ (31)

Comparing this with Equation 20 and Equation 25 we see that, with the
appropriate choices of $w_{\rm cool}=\lambda_{F}$ and $\ell_{\rm cool}$
respectively, this condition approximately holds.

In the unresolved case, Equation 30 must still hold, and therefore acts as a
condition on $\left\langle v_{\rm out}\right\rangle$ that effectively
contributes to setting its value if $w_{\rm cool}$ is numerically rather than
physically controlled. In fact, if we assume $w_{\rm cool}\approx\Delta x$ in
Equation 31, we get a similar relation to the resolved cases: $\left\langle
v_{\rm out}\right\rangle\sim\Delta x/t_{\rm cool}$, as foreshadowed by
Equation 27.

In the present simulations, the enhanced cooling in the shell can also lead to
a pressure dip if $w_{\rm cool}$ is not fully resolved. We will return to the
measured scaling behavior of $\left\langle v_{\rm out}\right\rangle$ in the
context of our simulations in Section 4.5.

Figure 3: Scaling relations relevant to the theoretical value for $\alpha_{p}$
given in Equation 34. Left: Evolution of $\alpha_{p}$ with $R_{\rm bub}$ for
$v_{t}(L)=10\,\,{\rm km}\,{\rm s}^{-1}$ and three values of the excess fractal
dimension, $d$ (other relevant parameters given in the text). Areas that are
dis-allowed by energy and momentum conservation are marked as shaded regions;
the upper limit is specific to the parameter choices made for this plot.
Right: Exponents that appear in Equation 34 as a function of $d$ for $p=1/3$
and $p=1/2$ in solid and dashed lines respectively.

###  2.6 Momentum Evolution with Turbulent Dissipation

We now wish to give a model for how the bubble should evolve in the situation
that dissipative physics is controlled by turbulent mixing at small scales, as
in Section 2.4.2, rather than numerics. We specify here to the case of
turbulent mixing since, as demonstrated by the right panel of Figure 2, we
expect turbulent heat transport to dominate over conduction in most real WBBs.
To this end, we can assume $\left\langle v_{\rm out}\right\rangle$ is equal to
a resolved, turbulent dissipation velocity $v_{t}(\ell_{\rm cool})$ (see
Equation 25). We can also replace $A_{\rm bub}$ with a fractal model for its
evolution, given by

$A_{\rm bub}(\ell)=4\pi R_{\rm bub}^{2}\left(\frac{R_{\rm
bub}}{\ell}\right)^{d}\,,$ (32)

where $d$ is the so-called ‘excess fractal dimension.’ This is discussed in
depth in Fielding et al. (2020); Lancaster et al. (2021a) and we return to it
in Section 4.4.

Substituting for $\left\langle v_{\rm out}\right\rangle$ and using $A_{\rm
bub}(\ell_{\rm cool})$ in Equation 11, we arrive at

$P_{\rm hot}=\frac{3L_{w}}{32\pi R_{\rm bub}^{2}v_{t}(L)}\left(\frac{R_{\rm
bub}}{L}\right)^{-d}\left(\frac{t_{\rm cool}}{t_{\rm
eddy}(L)}\right)^{\frac{d-p}{1-p}}\,.$ (33)

Since $t_{\rm cool}$ depends on $P_{\rm hot}$ (assuming the cooling layer is
isobaric; we will return to this below), we can rearrange the above to solve
for the dependence of $P_{\rm hot}$ on the other parameters. Doing this, and
setting $\alpha_{p}=4\pi R_{\rm bub}^{2}P_{\rm hot}/\dot{p}_{w}$ we can derive

$\alpha_{p}=\frac{3V_{w}}{16v_{t}(L)}\left[\frac{3L_{w}\Lambda(T_{\rm
pk})}{32\pi Lv_{t}^{2}(L)(k_{B}T_{\rm pk})^{2}}\right]^{\frac{p-d}{1-2p+d}}\\\
\times\left(\frac{R_{\rm bub}}{L}\right)^{\frac{d-2p+dp}{1-2p+d}}\,.$ (34)

It should be noted that the above relation puts no restrictions on physically
allowable values of $\alpha_{p}$ or the pressure. In reality, $\alpha_{p}\geq
1$ and has an upper limit imposed by energy conservation from the Weaver et
al. (1977) solution.

If we choose $p=1/3$ for Kolmogorov turbulence and $d=1/2$ as seems reasonable
in our hydrodynamical simulations (see Section 4.4), we have

$\alpha_{p}=\frac{3}{16}\frac{V_{w}}{v_{t}(L)}\left[\frac{3L_{w}\Lambda(T_{\rm
pk})}{32\pi Lv_{t}^{2}(L)(k_{B}T_{\rm pk})^{2}}\right]^{-1/5}\,;$ (35)

i.e. $\alpha_{p}$ is independent of the bubble size $R_{\rm bub}$, and also
(based on the $-1/5$ power) insensitive to the details of the cooling
function. As we can see in the right panel of Figure 3, the physically
motivated choice of $p$ and $d$ that results in a $\alpha_{p}$ which is
constant in time is somewhat serendipitous: large changes from $d=1/2$ result
in large changes in the exponent of $\alpha_{p}$ with $R_{\rm bub}$.
Fortunately, for reasonable values for the HD case, near $d=1/2$, this
exponent does not change drastically.

We can re-write Equation 35 in a more physically motivated way as

$\alpha_{p}=\frac{3}{16}\frac{V_{w}}{v_{t}(L)}\left[\frac{32\pi}{3}\frac{P_{\rm
hot}L^{3}(t_{\rm cool}/t_{\rm eddy}(L))}{L_{w}t_{\rm
eddy}(L)}\right]^{1/5}\,.$ (36)

Consider the last fraction in the square bracket term. The numerator is a
measure of the amount of thermal energy that is able to be stored in the
shocked gas: $P_{\rm hot}L^{3}$ is the thermal energy in a turbulent eddy and
$t_{\rm cool}/t_{\rm eddy}(L)$ is a measure of the fraction of energy that can
be retained before being lost to dissipation and cooling. The denominator is a
measure of how much energy is being provided to the bubble. If the energy that
the bubble can retain (numerator) is large compared to the rate at which
energy is provided to the shocked wind (denominator), then the momentum
enhancement is higher.

If we take quantities relevant to the high-resolution HD simulations that we
present below ($L=4\,{\rm pc}$, $v_{t}(L)=600\,{\rm km/s}$, $T_{\rm
pk}=1.59\times 10^{4}\,{\rm K}$, $\Lambda(T_{\rm pk})=2\times 10^{-22}\,{\rm
erg/s}\,{\rm cm}^{3}$, $L_{w}=4.875\times 10^{37}\,{\rm erg/s}$,
$V_{w}=3230\,{\rm km/s}$, $\overline{n}_{\rm H}=100\,{\rm cm}^{-3}$) the
quantity in square brackets above is roughly $1/5$ and the prefactor in front
is roughly unity, implying $\alpha_{p}\approx 0.2$ and constant in time. The
fact that $\alpha_{p}<1$ is simply because we have not explicitly enforced
momentum conservation in the derivation given above. Physically, however,
$\alpha_{p}\geq 1$ and is bounded from above by the corresponding energy-
driven solution from Weaver et al. (1977). These bounds are shown as grey
regions in the left-panel of Figure 3 for the $L_{w}$ and $V_{w}$ as quoted
above and an assumed mean background density that is identical to our
simulations.

The fact that the derived $\alpha_{p}<1$ for the parameters quoted above
simply means that turbulent mixing over a fractal surface is likely more than
capable of dispersing enough energy for the bubble’s evolution to approach the
momentum-conserving solution ($\alpha_{p}=1$) as posited in Lancaster et al.
(2021a). Motivated by the fact that turbulence is likely reduced in the, more
realistic, MHD case we instead plot the $\alpha_{p}$ evolution with $R_{\rm
bub}$ according to Equation 34 for an assumed $v_{t}(L)=10\,\,{\rm km}\,{\rm
s}^{-1}$, slightly above the observed Alfvén speeds in our MHD simulations.
Though the turbulence in the hot gas in our MHD simulations lies far above
this value (see Figure 15), this choice is likely more representative of
turbulence in resolved MHD mixing layers. This choice results in moderate
$\alpha_{p}$ enhancement above unity. Interestingly, in this formalism, fixing
other parameters and reducing the fractal dimension of the bubble results in a
higher $\alpha_{p}$, but one that decreases with time ($R_{\rm bub}$).

Figure 4: Visualization from our fiducial simulation at a resolution of
$N_{\rm res}=512$ approximately $t=0.80\,{\rm Myr}$ after the wind has been
turned on. All smaller panels depict slices through the $z=0$ midplane of the
simulation (clockwise from top-right): the gas temperature, the black line is
a $T=10^{6}\,{\rm K}$ iso-contour; the net cooling rate in the gas; loci in
the slice-plot of the thermal phases we use here; the thermal pressure of the
gas; the magnitude of the velocity field. The large panel in the top left
shows a three-dimensional rendering of the $T=10^{6}\,{\rm K}$ iso-temperature
surface within the simulation. The color of the surface corresponds to the
radius of that part of the surface and is included for clarity of viewing the
geometry of the surface. Light blue corresponds to surface elements far from
the wind source and bright red corresponds to surface elements close to the
wind source.  Figure 5: Visualization exactly analogous to Figure 4 but now
for the N512_mhd simulation. We have additionally replaced the phase splitting
slice in the bottom right panel with a schematic diagram that is similar to
the middle panel of Figure 1 but for the MHD case. We can see that structure
in the interface is suppressed by the presence of the magnetic field.

##  3 Methods

To provide a test of the theory of WBB evolution discussed above, we run
three-dimensional (magneto)hydrodynamic simulations, detailed below.

###  3.1 Numerical Methods

We use the Athena code (Stone et al., 2008). The numerical methods employed in
these simulations are nearly identical to those used in Lancaster et al.
(2021b); in particular we use the same ‘hybrid/subcell method’ for injecting
wind mass (tracked as a passive scalar field $\rho_{\rm wind}$) and energy
onto the grid and the same basic setup, which we will detail further below.
However, unlike Lancaster et al. (2021b), we employ the chemistry and
radiative cooling module of Kim et al. (2023b), which tracks the non-
equilibrium photochemistry of molecular, atomic, and ionized hydrogen as well
as the equilibrium chemistry of carbon- and oxygen-bearing species (${\rm C}$,
${\rm C}^{+}$, ${\rm CO}$, ${\rm O}$, and ${\rm O}^{+}$), and their effects on
the thermal heating and cooling of the medium. For the heating and ionization
of neutral gas, we adopt uniform FUV background radiation of Draine (1978) and
cosmic ionization rate per hydrogen $\xi_{\rm cr}=2\times 10^{-16}\,{\rm
s}^{-1}$ (Indriolo et al., 2015). For cooling by helium and metals at high
temperature ($T>2\times 10^{4}\,{\rm K}$), we use the tabulated collisional
ionization equilibrium (CIE) cooling function from Gnat & Ferland (2012). Our
hydrogen cooling is fully dynamical, dependent on the local radiation field
and ionization state of the gas and therefore not assumed to be CIE.

Our simulations are all run using the linearized Roe Riemann solver for flux
calculations, the piecewise linear method (PLM) for second-order spatial
reconstruction, and the unsplit second-order van Leer integrator for temporal
evolution (Stone et al., 2008, and references therein). We also include the
‘H-correction’ of Sanders et al. (1998) to prevent carbuncle instability when
the difference in signal speeds in different grid-aligned directions is
greater than $20~{}\,{\rm km}\,{\rm s}^{-1}$. We do not include gravity in
these simulations, and we use open/outflow boundary conditions.

We run simulations of both hydrodynamics (HD) and magneto-hydrodynamics (MHD).
In our MHD simulations, we employ Athena’s default constrained transport (CT)
method for evolving the magnetic field, which we will refer to as the CT-
Contact scheme (Gardiner & Stone, 2005). However, we found that this method is
subject to numerical instabilities in regions of our simulation with strong
source terms such as the bubble’s edge (where cooling is strong) and the wind
injection region. In order to suppress these numerical instabilities we
implement the upwind constrained transport (UCT) HLL scheme of Londrillo & Del
Zanna (2004) as described in Mignone & Del Zanna (2021). As the UCT-HLL method
is more diffusive999This diffusion only applies to the magnetic field, not to
the other fluid variables. than the CT-Contact scheme, we only employ it in
the regions where these instabilities were observed: in the wind-feedback
region and at the wind bubble interface, which we identify by regions with a
strong gradient in the magnetic field $\delta B/B>10$. We leave investigation
as to the origin of these numerical instabilities to future work.

Table 1: Simulation List.

Name | Hydro Solver | $N_{\rm res}$ | $\Delta x\,[{\rm pc}]$ | $\dot{M}_{w}/M_{*}\,[{\rm Myr}^{-1}]$
---|---|---|---|---
N128_hydro | HD | 128 | $0.31$ | $2.965\times 10^{-3}$
N256_hydro | HD | 256 | $0.16$ | $2.965\times 10^{-3}$
N512_hydro | HD | 512 | $0.08$ | $2.965\times 10^{-3}$
N128_mhd | MHD | 128 | $0.31$ | $2.965\times 10^{-3}$
N256_mhd | MHD | 256 | $0.16$ | $2.965\times 10^{-3}$
N512_mhd | MHD | 512 | $0.08$ | $2.965\times 10^{-3}$
N128_hydro_Mwmod | HD | 128 | $0.31$ | $10^{-2}$
N256_hydro_Mwmod | HD | 256 | $0.16$ | $10^{-2}$

Note. — Description of the parameters of our simulations.

###  3.2 Simulation Description

In all of the simulations detailed here, we first set up a turbulent, multi-
phase background before initializing a constant mechanical luminosity,
$L_{w}$, and mass loss rate, $\dot{M}_{w}$, wind at the center of the domain.
Here we choose parameters of our background medium to be representative of a
GMC with a total mass of $M_{\rm cl}=10^{5}\,M_{\odot}$ and a radius of
$R_{\rm cl}=20\,{\rm pc}$. In particular, we initialize our simulations within
a box of side length $L_{\rm box}=2R_{\rm cl}=40\,{\rm pc}$ with a given
number of resolving elements per grid-aligned direction $N_{\rm res}\equiv
L_{\rm box}/\Delta x$. We set the density everywhere to $n_{{\rm
H},0}=86.25\,\,{\rm cm}^{-3}$ and initialize a Gaussian, random-phase
turbulent velocity field with power spectrum $|v_{k}|^{2}\propto k^{-4}$ for
wavenumbers in the range $2\leq kL_{\rm box}/(2\pi)\leq 64$. The normalization
of the velocity field is chosen so that the total kinetic energy in the domain
is twice the binding energy, $2W_{\rm cl}=(6/5)GM_{\rm cl}^{2}/R_{\rm cl}$, of
a uniform-density spherical cloud of mass $M_{\rm cl}$ and radius $R_{\rm
cl}$. The initial one-dimensional velocity dispersion is therefore $\sim
3\,\,{\rm km}\,{\rm s}^{-1}$. In our simulations with magnetic fields we
additionally set a uniform magnetic field everywhere of
$\mathbf{B}_{0}=B_{0}\mathbf{\hat{z}}$ with $B_{0}=13.5\,\mu{\rm G}$. This
value was chosen based on setting the dimensionless mass-to-flux ratio equal
to 2 as in Kim et al. (2021).

We then allow the simulation to evolve, without any further turbulent driving,
until the total kinetic energy has become equal to $W_{\rm cl}$. The thermal
pressure is initially set everywhere to $P/k_{B}=6000\,{\rm K}\,\,{\rm
cm}^{-3}$. While this is above the equilibrium pressure at these densities,
the simulation domain quickly relaxes to its equilibrium pressure during the
initial turbulent evolution.

After the turbulent evolution, which lasts approximately $0.7\,{\rm Myr}$, the
wind is turned on. For all simulations, we assume a stellar mass of
$5000\,M_{\odot}$ with stellar wind luminosity per unit stellar mass of
$L_{w}/M_{*}=9.75\times 10^{33}\,{\rm erg}\,{\rm s}^{-1}\,M_{\odot}^{-1}$. We
consider two different values of the mass-loss rate per unit stellar mass:
$\dot{M}_{w}/M_{*}=10^{-2}\,{\rm Myr}^{-1}$ and $2.965\times 10^{-3}\,{\rm
Myr}^{-1}$. The first corresponds to values used in the simulations of
Lancaster et al. (2021b) while the latter is more representative of the value
derived from Starburst99 for a solar metallicity population with a Kroupa IMF
(Leitherer et al., 1999, 1992). These choices of mass loss rate correspond to
wind-velocities of $V_{w}=1759$ and $3230\,\,{\rm km}\,{\rm s}^{-1}$ and the
momentum input rate of $\dot{p}_{w}=8.795\times 10^{4}$ and $4.788\times
10^{4}M_{\odot}\,{\rm km}\,{\rm s}^{-1}\,{\rm Myr}^{-1}$, respectively, for
the given value of $L_{w}/M_{*}$. The higher wind velocity is our fiducial
value. The simulations with $\dot{M}_{w}/M_{*}=10^{-2}\,{\rm Myr}^{-1}$ and
$V_{w}=1759\,\,{\rm km}\,{\rm s}^{-1}$ will be referred to as our ‘mass loss
rate modified’ simulations, using the addendum Mwmod. The full parameter suite
of the simulations is detailed in Table 1.

For all simulations we inject mass, momentum, and energy using the hybrid
feedback method detailed in Lancaster et al. (2021b). Our effective injection
radius, $r_{\rm fb}/\Delta x=3.2,\,3.2,$ and $6.4$ for the simulations with
$N_{\rm res}=128,\,256,$ and $512$ respectively. Pittard et al. (2021)
demonstrated in one-dimensional simulations that it was necessary for $r_{\rm
fb}$ to be less than a maximum value, $r_{\rm inj,max}=\sqrt{\dot{p}_{w}/4\pi
P_{\rm amb}}$, in order for a WBB to be successfully inflated, with $P_{\rm
amb}$ the ambient pressure in the medium. For our default parameter values
$r_{\rm inj,max}=55\,{\rm pc}$, which is larger than the domain size of our
simulations. Thus, all of our simulations lie in the regime $r_{\rm fb}\ll
r_{\rm inj,max}/10$ recommended by Pittard et al. (2021) in order to properly
capture energy-driven bubbles.

Table 2: Definitions of gas phases.

Phase | Temperature Condition
---|---
WHIM | $10^{4}\,{\rm K}<T<10^{6}\,{\rm K}$
Warm Neutral Gas | $6085\,{\rm K}<T<10^{4}\,{\rm K}$
Thermally Unstable Gas | $181\,{\rm K}<T<6085\,{\rm K}$
Cold Neutral Gas | $T<181\,{\rm K}$

Note. — Temperature Definitions for colder phases. The wind phases are
specified in the text.

###  3.3 Measurement Details

Throughout the following section we will refer to various definitions of
thermal phases of the gas in our simulations. These are generally similar to
those outlined in Lancaster et al. (2021b), with the exception that we re-
label ‘Ionized Gas’ as ‘Warm-Hot Ionized Medium’ (WHIM) and make a few small
changes to the temperature ranges given our new cooling module; these are
given in Table 2. The wind phases are as in Lancaster et al. (2021b) where the
shocked wind is defined as all gas hotter than $10^{6}\,{\rm K}$ and with a
radial velocity less than half the wind velocity. The free wind is defined as
all gas with a radial velocity that is greater than half $V_{w}$ as well as
all gas contained within the feedback region.

In Figure 4 and Figure 5 we provide a rendering of our highest resolution
fiducial simulation: N512_hydro and N512_mhd respectively. It is clear that
the bubble’s geometry is very complicated and highly ‘folded’. It is also
clear from the top-right panel that the $T=10^{6}\,{\rm K}$ iso-temperature
surface encloses nearly all of the WBB’s volume.

We are particularly concerned in this work with properties of the gas near the
bubble’s surface. In order to precisely locate this surface within our
simulations we follow Fielding et al. (2020) and Lancaster et al. (2021b) by
using the marching_cubes algorithm from scikit-image (Van der Walt et al.,
2014). Specifically, we measure iso-temperature surfaces at
$\log_{10}\left(T/{\rm K}\right)=3.75-6$ in logarithmic steps of $0.25$ (10
temperatures) and using step sizes of $\ell/\Delta
x=1,\,2,\,4,\,8,\,16,\,32,\,\&\,64$. We can only perform such measurements as
post-processing on simulation snapshots, which we output from the simulations
every $\sim 0.01\,{\rm Myr}$ (100 per simulation). A three-dimensional
rendering of such an iso-temperature surface is given in the left-panel of
Figure 4 for $T=10^{6}\,{\rm K}$ and a spacing of $\ell/\Delta x=1$ for the
$N_{\rm res}=512$ simulation.

After using the marching_cubes algorithm to define the iso-temperature
surfaces we calculate the total area of the surface $A_{\rm bub}(T,\ell)$ and
interpolate the fluid variables to each facet of the surface that is returned
by the algorithm. We then calculate the component of the fluid velocity normal
to each facet of the surface, $v_{\rm out}$, and the dot product of the facet
normal and the radial vector, $\mathbf{\hat{n}}\cdot\mathbf{\hat{r}}$. The
quantities $\left\langle v_{\rm out}\right\rangle$ and
$\left\langle\mathbf{\hat{r}}\cdot\mathbf{\hat{n}}\right\rangle$, presented
below, correspond to the area-weighted average over all the faces of the iso-
temperature surfaces. Unless it is stated otherwise, we use the iso-
temperature surfaces defined by $T=10^{6}\,{\rm K}$ and measured on the
resolution scale $\ell=\Delta x$ for all analyses detailed below.

We will also be concerned with the distribution of gas in the simulation as a
function of temperature. To investigate this in detail we take each snapshot
of the simulations and compute the volume-weighted, pressure-weighted, and
cooling-weighted histograms in temperature in 200 logarithmically spaced
temperature bins between $10^{1.5}\,{\rm K}<T<10^{8.5}\,{\rm K}$. For these
measurement we restrict to gas that has a wind mass fraction $f_{\rm
wind}=\rho_{\rm wind}/\rho>10^{-4}$ and a radial velocity $v_{r}<V_{w}/2$ to
avoid contamination from the background and free wind respectively. However,
the former constraint does somewhat bias measurement in our MHD simulations,
where the momentum from the forward-shock of the wind bubble is carried away
from wind-polluted material by MHD waves.

Further details on measurements of specific physical quantities are presented
as needed in the remainder of the paper.

##  4 Results

Figure 6: Comparison of the total radial momentum carried by WBBs in two of
our simulations, N256_hydro (solid) and N256_hydro_Mwmod (dashed). These
simulations have the same energy input rate, $L_{w}$, but differ in their
input wind mass-loss rates, $\dot{M}_{w}$, and therefore have different wind
velocities, $V_{w}$, and input momentum rates, $\dot{p}_{w}$ (as indicated in
the key). For each simulation blue lines indicate the total measured radial
momentum as a function of time whereas black lines indicate the momentum that
would result from $p_{r}=\dot{p}_{w}t$, (i.e. $\alpha_{p}=1$). Despite
different input $\dot{p}_{w}$, the resulting momentum carried by each bubble
is quite similar. Figure 7: As in Figure 6 but now comparing HD (blue) and MHD
simulations (green). In each case, darker curves indicate higher numerical
resolution. Notably, the HD and MHD simulations have the opposite behaviour as
a function of resolution.

###  4.1 Radial Momentum

One of the most important aspects of WBBs is the dynamical effect they have on
their surroundings, exemplified by the bulk momentum that they inject (Keller
et al., 2022), so it is natural to first investigate this aspect of our
simulations. Here we measure the total radial momentum carried by the wind
bubble by computing evolution of the total radial momentum within the
simulation domain, and then subtracting out the evolution of the radial
momentum from a simulation that is identical except it does not include wind
feedback. This comparison to the no-feedback simulation is necessary since,
due to our outflow boundary conditions and the turbulent conditions, the total
radial momentum in the simulation is non-zero even when there is no feedback.
We use this method instead of the method of Lancaster et al. (2021b), which
measured total radial momentum in gas that included some fraction of wind
material, as we found that in our MHD simulations momentum is able to be
carried out of wind-polluted material by MHD waves, which has been observed
previously (e.g. Offner & Liu, 2018).

In Lancaster et al. (2021a) we posited that turbulent mixing in WBBs followed
by radiative cooling was strong enough to reduce the bubble evolution to
nearly the ‘momentum conserving’ regime. That is, the momentum communicated to
the surroundings by the bubble is very close to that put in by the wind,
$p_{r}\sim\dot{p}_{w}t$ or $\alpha_{p}\sim 1$. The simulations of Lancaster et
al. (2021b) showed measured values of $\alpha_{p}$ between 1 and 4 (see their
Fig. 8). Additionally, these simulations showed that the measured energy
retention fraction was consistent with the prediction for the measured values
of $\alpha_{p}$. In Lancaster et al. (2021a, b) we further argued that in the
case of a radial flow in the bubble, the speed at the edge of the bubble,
$\sim(3/16)V_{w}/\alpha_{p}$, should match the effective diffusion velocity
that carries material into the shell. Here, from Equation 15, we have taken a
further step to see that the value of $\alpha_{p}$ is actually determined by
the interplay of the bubble’s geometry and dissipation at its surface in time.
The quasi-steady state assumption imposes the global constraint
$\alpha_{p}\approx(R_{\rm bub}/R_{\rm f})^{2}$, which has a lower limit of
unity, but there is no intrinsic reason why the ratio should be exactly 1 or
why it should be constant in time.

As we can see upon inspection of Figure 6, this is indeed not the case. Here
we show the total radial momentum evolution in our default HD simulations with
two different input wind momenta, $\dot{p}_{w}$. As described in Section 3,
this difference is mediated by a difference in $\dot{M}_{w}$, while both
simulations have the same mechanical wind luminosity, $L_{w}$. Though the
input momenta, $\dot{p}_{w}$, are different by about a factor of 2 between
these simulations, the momentum they impart to background are the same to
within 10%. That is, $\alpha_{p}$ is a factor $\sim 2$ larger in the
simulation with lower momentum injection rate (having lower $\dot{M}_{w}$,
higher $V_{w}$). As we will discuss further in Section 4.3, the agreement in
$\dot{p}_{r}$ is largely due to the fact that momentum input to the shell is
governed by dissipation processes at the bubble’s surface (which set $P_{\rm
hot}$ through through Equation 11), and these dissipation processes remain
largely unchanged between these two simulations. For fixed $\left\langle
v_{\rm out}\right\rangle A_{\rm bub}/R_{\rm bub}^{2}$, Equation 13 says that
$\dot{p}_{r}\propto L_{w}$, i.e. the same value for the two models shown,
consistent with the simulation result of Figure 6. At the same time, Equation
15 says that $\alpha_{p}\propto V_{w}$; this implies a larger value for
$\alpha_{p}$ in the fiducial simulation, which has higher $V_{w}$ and lower
$\dot{p}_{w}$.

The second interesting issue we would like to introduce in Figure 7 is the
opposite resolution dependence of the total radial momentum in our HD and MHD
simulations. While the HD simulations have decreasing momentum output with
increased resolution, the MHD simulations have the opposite behaviour. As we
will discuss below, and should be clear from Equation 15, this must be solely
due to differences in the relative scaling of dissipation and geometry between
the two types of simulation.

Figure 8: A validation of the driving dynamics of the bubble laid out in
Section 2.2. Line colors in all panels are as in Figure 7 and all comparisons
are shown as a function of $R_{\rm bub}$. Top Panel: The rate of change of the
total radial momentum carried by the WBBs divided by the radial force given in
Equation 9. Bottom Panel: Same as top panel, now dividing by Equation 10.
Figure 9: Left Panel: The hot gas pressure as a function of time as measured
in the simulations (solid lines), showing excellent agreement with the
measurement derived from Equation 11 (dashed lines). The MHD simulations are
left out for clarity, but they also show excellent agreement with Equation 11.
Right Panel: We show the bubble’s radius as a function of time in the highest
resolution HD (blue) and MHD (green) runs. For comparison we also show, as
dashed lines, the radial evolution derived for the constant-momentum input
solution of Lancaster et al. (2021a) with $\alpha_{p}=3.75$ in black, as well
as the constant cooling fraction solution of El-Badry et al. (2019) with
$\theta=0.9,\,0.95$ in red and yellow respectively.

###  4.2 Validity of the Momentum Equation

In Section 2.2 we argued that it is valid to describe the instantaneous force
exerted on the bubble’s surface using Equation 9. In Lancaster et al. (2021b),
this was indeed shown to be a very good approximation in the HD simulations
presented there. Before presenting the details of the bubble’s structure we
would like to validate this result in the current simulations.

In the top panel of Figure 8 we show the rate of change of the total radial
momentum carried by the bubbles (the time derivative of the curves in Figure
7) divided by the force as computed from Equation 9. Specifically, the
pressure is the volume-averaged value within the shocked wind region and
$R_{\rm bub}$ is as computed from Equation 7. We can see that the momentum
rate of change is always within an order unity factor from the prediction with
smooth variations over time (represented here by $R_{\rm bub}$) and is an
especially good prediction for the MHD simulations. In the bottom panel of
Figure 8 we show a comparison of the same force to the second prediction given
in Equation 10, which shows even better agreement.

The fact that these predictions are slightly under-estimating the true force
could be due, in part, to the fact that we are ignoring the contribution of
the Reynolds stress term, which should be quite sub-dominant. Further under-
estimates could be due to the slight inaccuracy of the ‘effective area’
assumption behind Equation 9. This point is consistent with the fact that the
equation is a better approximation in the MHD simulations (top panel of Figure
8), where the bubble’s surface is less folded and one might expect $A_{\rm
bub}$ to be closer to $4\pi R_{\rm bub}^{2}$. Regardless, the agreement seen
here is quite remarkable considering the extremely complex bubble geometry
depicted in Figure 4.

In the left panel of Figure 9 we show $P_{\rm hot}/k_{B}$, as a function of
time along with the prediction of Equation 11 with $\left\langle v_{\rm
out}\right\rangle$ and $A_{\rm bub}$ measured on iso-temperature surfaces with
$T=10^{6}\,{\rm K}$. It is clear from the correspondence between the measured
pressure and the predicted pressure that Equation 11 is an excellent
description of the bubble dynamics. This validates the picture discussed in
Appendix A and Section 2.2: that the bubble’s interior is nearly steady and
that its pressure is determined by the dissipation (leading to $\left\langle
v_{\rm out}\right\rangle$) and the geometry ($A_{\rm bub}$) of its interface
with surrounding gas, as exemplified in Equation 11 and Equation 15.
Additionally, the fact that the prediction shown in Figure 8 is given based on
measurements from an iso-temperature surface at $T=10^{6}\,{\rm K}$, which
lies at the bubbles’ furthest outer edges (see the iso-contour in the top
right panel of Figure 4), indicates that the steady-state bubble interior
structure is a good approximation.

It is worth emphasizing that Figure 8 and the left panel of Figure 9 depict
completely independent comparisons reflecting different aspects of bubble
dynamics. Figure 8 shows that the measured radial force is equal to the bubble
pressure multiplied by an effective area computed from the bubble volume
(which is much smaller than the actual area of the bubble’s surface). The left
panel of Figure 9 shows that the measured flow of thermal energy out of the
bubble’s (highly complex) bounding surface is consistent with the steady rate
at which energy is being pumped in by the wind at its center.

It is also clear from the left panel of Figure 9 that the pressure in the hot
gas is not a converged quantity with resolution, even if we can consistently
predict its value based on measurements at the bubble’s interface. In
particular, at a given time, higher resolution HD simulations have a lower
bubble pressure. This resolution dependence then must be manifest in the
resolution dependence of the quantities measured at the bubble’s interface,
which we explore below.

Finally, in the right panel of Figure 9 we show the radial evolution of the
bubbles in our highest resolution simulations in both HD and MHD as functions
of time. For comparison we show the radial evolution expected for momentum-
driven bubbles (Steigman et al., 1975; Lancaster et al., 2021a)

$R_{\rm
bub}=\left(\frac{3\alpha_{p}\dot{p}_{w}t^{2}}{2\pi\overline{\rho}}\right)^{1/4}\,,$
(37)

where, in Figure 9, we use $\alpha_{p}=3.75$. We also show the radial
evolution expected for a bubble which cools a set fraction, $\theta$, of its
energy in its interface (El-Badry et al., 2019)

$R_{\rm
bub}=\left(\frac{125L_{w}(1-\theta)t^{3}}{154\pi\overline{\rho}}\right)^{1/5}\,,$
(38)

where in Figure 9 we use $\theta=0.9$ and $0.95$.

It is clear that the radial evolution of the bubbles in the right panel of
Figure 9 are better represented by the ‘momentum-driven’ solution (which
requires cooling to increase as time goes on (Lancaster et al., 2021a)) than
by the constant $\theta$ solutions. We investigate the details of energy loss
in the bubbles in the next section.

Figure 10: Summary characteristics of interior and interface properties of the
bubble as a function of resolution and equivalent area $4\pi R_{\rm bub}^{2}$,
both with (bottom panels) and without (top panels) magnetic fields. The
resolution of each simulation is indicated by the legends in panels (e) and
(f). Panels (a) & (b): Hot gas pressure, which shows at later times the same
opposing trends with resolution as in Figure 7. For reference we show the
pressure for an efficiently-cooled bubble (Equation A4 with $R_{\rm f}\to
R_{\rm bub}$) Panels (c) & (d): The outward velocity averaged over the
bubble’s surface, $\left\langle v_{\rm out}\right\rangle$ (see Equation 12),
decreasing as resolution increases. These values are all much smaller than the
expected value for unimpeded radial flow $V_{w}\left(R_{\rm f}/R_{\rm
bub}\right)^{2}/4$, which is nearly constant with time and in the range
$200-300\,\,{\rm km}\,{\rm s}^{-1}$ in the simulations over the time range
shown. Panels (e) & (f): The area of the bubble surface, $A_{\rm bub}$,
increasing as resolution increases. We additionally include the $A_{\rm
bub}=4\pi R_{\rm bub}^{2}$ 1-to-1 line in black to emphasize that the true
bubble surface areas are much larger than their sphere-equivalent surface
areas.

###  4.3 The Interface

In Figure 10 we show key properties of the bubble’s interior and interface as
a function of its sphere-equivalent area, $4\pi R_{\rm bub}^{2}$. In the left
hand panels of Figure 10 we can see the difference in HD (top) vs. MHD
(bottom) bubble dynamics that we saw in Figure 7: the opposite trend with
resolution of the bubble’s pressure (and hence radial force through Equation
9; see Figure 8).

Since Figure 9 shows that Equation 11 is an accurate predictor of the bubble’s
interior pressure, we can attribute the difference in resolution dependence of
$P_{\rm hot}$ to a difference in resolution dependence of the quantities
$\left\langle v_{\rm out}\right\rangle$ and $A_{\rm bub}$ defined at the
interface, which we show in the middle and right panels of Figure 10. In
particular, since $P_{\rm hot}\propto(\left\langle v_{\rm out}\right\rangle
A_{\rm bub})^{-1}$ and consistently decreases at higher resolution in the HD
simulations, the increase of $A_{\rm bub}$ must be steeper than the decrease
of $\left\langle v_{\rm out}\right\rangle$ at smaller $\Delta x$. In the HD
simulations, we have found dependence on resolution varying roughly as
$\left\langle v_{\rm out}\right\rangle\propto(\Delta x)^{1/4}$ and $A_{\rm
bub}\propto(\Delta x)^{-1/2}$: thus, increased surface area wins over
decreased outflow velocity normal to the surface, which in these simulations
ultimately is due to numerical dissipation. In contrast, the relative scalings
with resolution in the MHD simulations nearly compensate for one another:
$\left\langle v_{\rm out}\right\rangle\propto(\Delta x)^{1/4}$ and $A_{\rm
bub}\propto(\Delta x)^{-1/4}$. In the end, however, the decreased dissipation
marginally wins over the increased bubble surface area, resulting in the
increased pressure and therefore increased momentum input to the surroundings.

The question remains as to the origin of this dependence on numerical
resolution, $\Delta x$. In addition, we would like to understand why
$\left\langle v_{\rm out}\right\rangle$ decreases and $A_{\rm bub}/(4\pi
R_{\rm bub}^{2})$ increases for larger bubble radius (a proxy for later time
in Figure 10). We investigate the scalings of $A_{\rm bub}$ in Section 4.4 and
$\left\langle v_{\rm out}\right\rangle$ in Section 4.5, below.

###  4.4 Scaling of $A_{\rm bub}$

As has been discussed extensively in the literature, dynamically unstable
interfaces can lead to fractal, self-similar structures (Fielding et al.,
2020; Tan et al., 2021; Lancaster et al., 2021a, b). We will show in this
section that the scaling of $A_{\rm bub}$ with time and resolution can be
entirely explained through a fractal model for the geometry of these WBBs. We
will leave discussion as to the origin of this fractal structure for Section
5.

There are several ways to think about fractal structure. Traditionally in
mathematics (Mandelbrot, 1983; West et al., 1997) and in the astrophysical
literature (Meneveau & Sreenivasan, 1987; Shibata & Tanuma, 2001; Federrath et
al., 2009), describing something as a fractal means parameterizing how some
aspect of its structure, say $f$, varies as a function of scale, say $\ell$,
as a power law in that scale, $f(\ell)\propto\ell^{a}$ for some $a$. For
example, how the measurement of the size of the boundary of some region varies
as a function of the physical scale on which you measure it: like measuring
the coastline of England with a mile long ruler (not able to describe the
structure of small bays and inlets) versus a foot long ruler (maybe not able
to describe every rock). For our current purposes, and generally in the
mixing-layer literature, we are interested in the area of the WBB interface
and how this varies with scale. We can describe this scaling as (see e.g.
Lancaster et al., 2021b; Fielding et al., 2020)

$d_{\ell}\equiv-\frac{d\log A_{\rm bub}(\ell)}{d\log\ell}\,.$ (39)

We see that $d_{\ell}=0$ corresponds to an interface with no finer-scale
structure (area does not get larger when measured at smaller scales) and
$d_{\ell}>0$ corresponds to the formation of smaller scale structure (area
increases as scale on which it is measured decreases).

In mathematics, fractal structures can be infinitely self-similar,
continuously forming structure on smaller and smaller scales. Physical
structures, in contrast, have to be regulated and eventually become smooth on
some scale. This implies that $d_{\ell}$ is not constant as a function of
$\ell$ but must vary, as it must return to $0$ at the scale on which structure
becomes smooth. In the context of real WBBs this is at least the mean free
path of individual particles, at which point the idea of a boundary between
fluid elements becomes meaningless, and is more likely to have to do with
dissipative physics. In our ideal (M)HD simulations, however, this scale is
the resolution scale101010It is interesting to note that, in the mathematical
world of ideal fluid dynamics, there is no regulatory scale, and one can form
infinitely self-similar structures..

Figure 11: Fractal dimension of the WBB interfaces measured using Equation 39
in the $N_{\rm res}=512$ versions of both the HD (blue) and MHD (green)
simulations as a function of $R_{\rm bub}$. Horizontal black lines indicate
the fractal dimensions used for the scaling shown in Figure 12: $d=1/2,\,1/4$.
The vertical dashed line indicates when the magnetic pressure of the bubble’s
shell comes into rough equilibrium with the bubble’s interior thermal pressure
in the MHD simulation.

In Figure 11 we show the fractal dimension as measured using Equation 39 at
$\ell=8\Delta x$ in the $N_{\rm res}=512$ versions of the HD and MHD
simulations. Specifically, this is calculated as the logarithmic difference of
$A_{\rm bub}$ measured on $\ell=4\Delta x$ and $16\Delta x$. We choose this
scale as it is the scale that seems to be least affected by the ends of the
“fractal dynamic range” so that it is to the best extent probing the true
fractal dimension. We see from Figure 11 that the fractal dimension is
relatively constant as a function of radius/time in the HD simulations. In
contrast, the MHD simulation has a clear decrease in its fractal dimension
over time. If we attribute some part of the fractal structure of these bubbles
to the growth of MHD instabilities at the interface, this decrease in the
fractal structure can be understood as a suppression of these instabilities as
the bubble expands, sweeping up the background magnetic field and causing the
field strength to grow at the interface. This process is schematically shown
in the bottom right panel of Figure 5. Indeed, this period of decreasing $d$
coincides with the bubble’s thermal pressure coming into equilibrium with the
magnetic pressure of its shell, indicated by the vertical dashed line in
Figure 11. We indicate the time when the shell magnetic pressure first comes
within 20% of the bubble’s thermal pressure with a vertical line in Figure 11.

Figure 12: Demonstration of the fractal nature of the bubble’s interfaces. The
surface area is shown as a function of linear dimension, $R_{\rm bub}$,
additionally scaling with resolution, in comparison to power law relations.
See text.

Instead of considering the largest scales as fixed and investigating the
fractal nature of the interface by examining the structure on smaller and
smaller scales, we can also do the reverse. That is, we can consider the
smallest scales as fixed (the resolution or dissipative scale) and investigate
how the structure of the interface changes as it becomes larger. Investigating
the fractal nature of the bubble in this way is a distinct feature of problems
like the WBB (as opposed to cloud-crushing or TRML studies) as the outer scale
actually varies in time. We can then alternatively define the excess fractal
dimension as

$d_{R}\equiv\frac{d\log A_{\rm bub}}{d\log R_{\rm bub}}-2\,.$ (40)

In the above we assume that $A_{\rm bub}$ is measured on some fixed small
scale. In practice we use the resolution scale of the simulations. This
definition includes the $-2$ since we would expect that the surface area of a
3-dimensional structure scales with the square of its linear dimension;
$d_{R}$ then quantifies the amount by which it grows more than this.

In the context of Equation 40 we can show the areas of the bubble at fixed
linear scale, $R_{\rm bub}$, as a pure geometrical measurement of the fractal
qualities of the bubbles. However, at higher resolution there is a larger
dynamic range on which to form structure. We can account for this using
Equation 39 assuming a given $d_{R}$ and backing out the factor needed to
multiply the low-resolution simulations by in order to match the high
resolution simulations: $\left(512/N_{\rm res}\right)^{d_{R}}$.

This comparison is shown in Figure 12 where we take $d_{R}=1/2$ for the HD
simulations and $d_{R}=1/4$ for the MHD simulations for this resolution
correction. These choices are generally consistent with the mean values shown
in Figure 11. The fact that the lines at differing resolution match very well
with this correction points to the fact that it is simply the bubble’s fractal
structure that is causing this scaling.

We additionally show as black lines the scalings of $A_{\rm bub}$ with $R_{\rm
bub}$ for the same values of $d_{R}$ assumed in the resolution scaling.
Specifically, we show

$A_{\rm bub}\left(\frac{N_{\rm res,hr}}{N_{\rm res}}\right)^{d_{R}}=C4\pi
R_{\rm bub}^{2}\left(\frac{R_{\rm bub}}{\Delta x_{\rm hr}}\right)^{d_{R}}$
(41)

where $N_{\rm res,hr}=512$ and $\Delta x_{\rm hr}\approx 0.08\,{\rm pc}$ is
the resolution in the highest resolution simulations. We have set
$C=0.56,\,0.80$ and $d_{R}=1/2,\,1/4$ for the HD and MHD simulations,
respectively. The fact that these lines match very well with the simulations
indicate the self-consistency of the fractal picture and how it explains the
behaviour of $A_{\rm bub}$ with resolution and scale. Furthermore, the fact
that the black line is an over-estimate of the true $A_{\rm bub}$ in the MHD
simulations at late times is consistent with Figure 11 where we see a decrease
in the fractal dimension, below $1/4$, at late times.

In conclusion, we find that the fractal structure of the bubble’s surface,
with excess dimension $d\approx 1/2$ in the HD case and $d\sim 0.25$ in the
MHD case, is able to account for the scaling of $A_{\rm bub}$ with numerical
resolution seen in Figure 10e,f.

Figure 13: An empirical argument for the scaling of $\left\langle v_{\rm
out}\right\rangle$. Line-styles are as specified in Figure 7. Left Panel: The
measured value of $\left\langle v_{\rm out}\right\rangle$ compared to that
predicted by Equation 43. The one-to-one relation is shown for reference in
black. Middle Panel: The width of the cooling layer, $w_{\rm cool}$, as a
function of the simulation’s numerical resolution, $\Delta x$. As black lines
we show scaling with resolution of $w_{\rm cool}=\Delta x$ and $w_{\rm
cool}=\Delta x/2$. Right Panel: The factor by which the pressure dips from the
bubble interior to the pressure at peak cooling $f_{\rm dip}=P(T_{\rm
pk})/P_{\rm hot}$ as a function of $\left\langle v_{\rm out}\right\rangle$.

###  4.5 Value and Scaling of $\left\langle v_{\rm out}\right\rangle$

In Section 2.4 we discussed the dissipative scales that need to be resolved in
order for $\left\langle v_{\rm out}\right\rangle$ to be set by physical
processes. When these processes are not well resolved and cooling is strong
enough so that the local value of $v_{\rm out}$ (the velocity of the fluid
into the surface, relative to the surface) is large compared to the velocity
of the surface, $W_{\rm bub}$, then the local numerical dissipation is set by
the local grid velocity, $v_{\rm out}$. But, as we motivated in Section 2.4,
it is the dissipation which should be setting $\left\langle v_{\rm
out}\right\rangle$ in the first place.

In Section 2.5 we described how cooling provides a constraint on the value of
$\left\langle v_{\rm out}\right\rangle$. In particular, through its role in
setting $P_{\rm hot}$ (through Equation 11), $\left\langle v_{\rm
out}\right\rangle$ sets the pressure in the cooling layer (assuming the
cooling layer is iso-baric with the bubble interior) and hence sets the
cooling rate. Requiring a given cooling rate then requires a specific value of
$\left\langle v_{\rm out}\right\rangle$.

For our numerical simulations, in fact the cooling layer is not perfectly iso-
baric with the hot bubble interior, but instead has a pressure dip:

$f_{\rm dip}\equiv\frac{P(T_{\rm pk})}{P_{\rm hot}}\,.$ (42)

Accounting for this in Equation 30 gives us

$\left\langle v_{\rm out}\right\rangle_{\rm pred}=\frac{3\sqrt{\Lambda(T_{\rm
pk})V_{\rm pk}L_{w}}}{8A_{\rm bub}k_{B}T_{\rm pk}}f_{\rm dip}\,,$ (43)

where we have included the subscript for clarity of comparison in Figure 13.
If we remove the factors that are explicitly dependent on the details of
cooling we have

$\left\langle v_{\rm out}\right\rangle_{\rm pred}\propto f_{\rm
dip}\sqrt{\frac{w_{\rm cool}}{A_{\rm bub}}}$ (44)

where we have used $V_{\rm pk}=w_{\rm cool}A_{\rm bub}$ with $w_{\rm cool}$
the thickness of the cooling layer.

In the left hand panel of Figure 13 we show $\left\langle v_{\rm
out}\right\rangle$ measured from the simulations against the prediction of
Equation 43, where we take $V_{\rm pk}$ as the volume of gas with temperatures
$10^{4}<T/{\rm K}<4.8\times 10^{4}$. This is chosen to match the rough
temperature range over which the highest rate of cooling occurs (see Appendix
D). We measure $T_{\rm pk}$ and $P(T_{\rm pk})$ using the cooling and pressure
weighted temperature histograms as detailed in Section 3.3, and
$\Lambda(T_{\rm pk})$ measured by interpolating the collisional ionization
equilibrium (CIE) cooling curve of Kim & Ostriker (2017).

Using this CIE $\Lambda(T)$ is not strictly faithful to the non-equilibrium
cooling that we are employing in these simulations. While we could employ the
equilibrium cooling curve from Kim et al. (2023b) for the non-photoionized
case (e.g. top panel of Figure 3 there), this is negligibly different than Kim
& Ostriker (2017) for the temperature ranges of interest here. Additionally,
the true cooling in our simulations will differ from the equilibrium curve of
Kim et al. (2023b) due to the combination of turbulent and numerical mixing of
the electron fraction, $x_{e}$, across the cooling layer. The CIE cooling
curve is therefore a useful simplification. Using a single value of $T_{\rm
pk}$ is of course a great simplification itself, since there is a wide range
of temperatures where cooling occurs in these simulations. It is less obvious,
but using a single value of $P(T_{\rm pk})$ is also strictly inaccurate since,
in these simulations, there are large pressure fluctuations in the cooling
layer.

With these caveats in mind, the correspondence seen in the left panel of
Figure 13 acts as a test of Equation 44, that is, the scaling of $\left\langle
v_{\rm out}\right\rangle$ and its dependence on the width of the cooling layer
$w_{\rm cool}$, the area of the bubble’s surface, $A_{\rm bub}$, and the
pressure dip between the hot gas bubble and the cooling layer, $f_{\rm dip}$.
In particular, this shows that the resolution and time dependence of
$\left\langle v_{\rm out}\right\rangle$ can be explained through the
resolution and time scaling of these quantities.

We investigate the resolution scaling of $w_{\rm cool}$ (which is independent
of time in all of our simulations) in the middle panel of Figure 13. This
width appears to scale roughly linearly with $\Delta x$. This is not
unexpected; in the strong-cooling limit (McCourt et al., 2018; Fielding et
al., 2020; Tan et al., 2021) the cooling region will always occupy roughly a
one-cell thickness volume over the interface unless the dissipative scales are
resolved. The fact that $w_{\rm cool}<\Delta x$ illustrates how quickly this
transition occurs. Indeed, for a wider definition of the “cooling region” that
includes hotter, not as rapidly cooling gas we would find $w_{\rm cool}>\Delta
x$ and scaling linearly with it. Together this means that the edge of the
bubble always has approximately one cell at intermediate temperatures, but
this is not always at the most rapidly cooling temperatures.

Finally, we’d like to investigate a relationship between the pressure at the
cooling interface and $\left\langle v_{\rm out}\right\rangle$ or equivalently
$P_{\rm hot}$. In the right panel of Figure 13 we show the fractional dip in
pressure from the wind-bubble interior to the peak-cooling gas, $P(T_{\rm
pk})/P_{\rm hot}$, as a function of $\left\langle v_{\rm out}\right\rangle$.
While $f_{\rm dip}$ is strongly dependent on resolution itself, it is clear
from this panel that $\left\langle v_{\rm out}\right\rangle$ and $f_{\rm dip}$
follow a relation to one another that is roughly independent of resolution but
somewhat dependent on whether the medium is magnetized; thus, there may in
principle be dependence on the numerical solver (MHD vs. HD). As we discuss in
Appendix C, this difference likely has to do with the details of the Riemann
problem being solved at the interface and how these differ between the HD and
MHD cases. With a fully resolved turbulent layer, this pressure dip should
tend to unity (i.e. disappear, see Fielding et al., 2020). In any case, this
relationship follows the regulatory pattern we set out at the beginning of
this section: if $\left\langle v_{\rm out}\right\rangle$ increases not only
does $P_{\rm hot}$ decrease (through the bubble structure relationship laid
out in Appendix A) but the pressure dip also gets more extreme, leading to an
even smaller $P(T_{\rm pk})$ and thus much decreased cooling. The opposite is
true for decreased values of $\left\langle v_{\rm out}\right\rangle$.

If we accept this final relation as empirical, then it closes the loop on our
causality argument for the dynamics of the bubble: $\left\langle v_{\rm
out}\right\rangle$ adjusts so that the pressure in the cooling gas produces a
fixed $\dot{E}_{\rm cool,pk}$ given the volume of cooling gas $V_{\rm
pk}=w_{\rm cool}A_{\rm bub}$. This adjustment of $\left\langle v_{\rm
out}\right\rangle$ is not only dependent on how $\left\langle v_{\rm
out}\right\rangle$ relates to $P_{\rm hot}$ but also how it relates to $f_{\rm
dip}$. This means the detailed value of $\left\langle v_{\rm
out}\right\rangle$ has a non-trivial dependence on the details of the Riemann
solver and cooling, which determine $f_{\rm dip}$. With $\left\langle v_{\rm
out}\right\rangle$ determined, $P_{\rm hot}$ then sets the dynamics of the
bubble. In particular, the fact that the whole dynamics problem comes down to
a question of energy regulation through cooling explains why the resulting
momentum is independent of $\dot{p}_{w}$ but has the same value as long as
$L_{w}$ is the same, as we saw in Figure 6. That is, after $P_{\rm hot}$
adjusts such that the cooling rate matches the input energy rate in Equation
28, the same $P_{\rm hot}$ sets the expansion rate through Equation 9 (see
Figure 8).

Given the complexity of the picture laid out above, and the fact that in the
end it is dependent upon an empirical relationship between $\left\langle
v_{\rm out}\right\rangle$ and $f_{\rm dip}$, it is still quite possible that
physical turbulent dissipation plays some role in the setting of $\left\langle
v_{\rm out}\right\rangle$ as well. We will return to this in Section 4.7.

Our detailed analysis has produced significant insight into bubble dynamics
and the reasons for resolution dependence in numerical simulations, but a
number of mysteries remain. In particular, the relationship between the exact
fraction of energy lost in the cooling layer, the fractional dip in the
pressure, $f_{\rm dip}$, and $\left\langle v_{\rm out}\right\rangle$ remains
uncertain. These issues are likely explained by the detailed nature of
numerical dissipation at the WBB interfaces and how this numerical dissipation
is driven or modified by the presence of strong cooling and a strong mean
velocity. As these details are not dependent upon the three-dimensional,
dynamical nature of the simulations at hand, they are more likely to be made
clear through the investigation of much simpler, one-dimensional simulations,
where the parameters of the problem can be more easily tuned by hand. We leave
this detailed investigation for future work.

Finally, we have found an empirical relation for $\left\langle v_{\rm
out}\right\rangle$ based solely on $R_{\rm bub}$ and $\Delta x$ as

$\left\langle v_{\rm out}\right\rangle\propto R_{\rm bub}^{-1/2}\Delta
x^{1/4}\,.$ (45)

The success of this empirical scaling is shown in Figure 14. We will
additionally note that there is no hidden dependence in the above scaling
relations on $\dot{p}_{w}$, $V_{w}$, or $\dot{M}_{w}$: the “Mwmod” simulations
outlined in Table 1 have nearly identical values of $\left\langle v_{\rm
out}\right\rangle$ to their unmodified counter-parts.

We emphasize that this empirical scaling of $\left\langle v_{\rm
out}\right\rangle$ should be understood as describing simulations where the
cooling length in the boundary layer is not resolved, so instead diffusion is
numerical. Additionally, in the resolved case we would expect $\left\langle
v_{\rm out}\right\rangle$ to depend on turbulence and/or the true thermal
conductivity, as discussed in Section 2.6.

Figure 14: The mean outflow velocity $\left\langle v_{\rm out}\right\rangle$
as a function of $R_{\rm bub}$ for each of the simulations at different
resolution. Both HD (top panel) and MHD (bottom panel) are shown. The
empirical scaling relation of Equation 45 with a chosen norm (equal to
$302,\,350\,{\rm km/s}\,{\rm pc}^{1/4}$ in the HD and MHD simulations
respectively) is shown as the black lines. Colored line styles for different
resolution are as indicated in Figure 7.

###  4.6 A Geometry-Driven Solution

Combining the arguments from Section 4.4 and Section 4.5 with Equation 15 we
can derive a direct relationship between the momentum enhancement factor,
$\alpha_{p}$, the geometric properties of the interface, and the numerical
properties of diffusion across the interface. In particular, accepting a
certain fractal dimension, $d$, that describes the bubble’s surface

$A_{\rm bub}\propto R_{\rm bub}^{2+d}\Delta x^{-d}\,,$ (46)

along with Equation 45 and Equation 15, we arrive at

$\alpha_{p}\propto V_{w}R_{\rm bub}^{1/2-d}\Delta x^{d-1/4}\,.$ (47)

If we then apply $d=1/2$ for the HD simulations and $d=1/4$ for the MHD
simulations, as seems reasonable from Section 4.4, we arrive at

$\alpha_{p,{\rm HD}}\propto V_{w}\Delta x^{1/4}$ (48) $\alpha_{p,{\rm
MHD}}\propto V_{w}R_{\rm bub}^{1/4}\,.$ (49)

If we instead apply $d=1/5$, as seems appropriate at late times in the MHD
simulations (see Figure 11), then we have $\alpha_{p,{\rm MHD}}\propto
V_{w}R_{\rm bub}^{3/10}\Delta x^{-1/20}$. We emphasize again that the scalings
in Equation 48 and Equation 49 are intended as an interpretation of our
current simulations, in which numerical diffusion affects the value of
$\left\langle v_{\rm out}\right\rangle$ as empirically characterized in
Equation 45. For the fully resolved case, there could be no dependence of
$\alpha_{p}$ on $\Delta x$.

These relations elucidate several mysteries:

  * 1.

The “momentum-conserving”-like behavior of the HD simulations, given the lack
of $R_{\rm bub}$ dependence in Equation 48, as well as the very weak
dependence on resolution.

  * 2.

As all of the above relations are linearly dependent on $V_{w}$, the resulting
output momentum is $\dot{p}_{r}\propto\alpha_{p}\dot{p}_{w}\propto L_{w}R_{\rm
bub}^{1/2-d}\Delta x^{d-1/4}$, this explains the lack of $\dot{p}_{w}$ (or
$V_{w}$) dependence in Figure 6. In fact, this is explicitly not a “momentum-
driven” solution, but an energy driven solution with losses to cooling that
are large and increase over time. We should note that there is likely some
implicit additional dependence of $\left\langle v_{\rm out}\right\rangle$ on
$L_{w}$ (and possibly $V_{w}$) given the causal structure of how $\left\langle
v_{\rm out}\right\rangle$ is set, which we have not accounted for here. This
dependence is apparent in the $\alpha_{p}$ dependence on $L_{w}$ in Lancaster
et al. (2021b) (see Figure 8).

  * 3.

The different momentum scaling of the MHD simulations, which have explicit
$R_{\rm bub}$ dependence, and their opposite resolution dependence behaviour
in Figure 6 from the HD models. This is apparent especially at late times
where the fractal dimension is much reduced and the scaling of $\alpha_{p}$
with $\Delta x$ becomes a negative power.

Figure 15: Turbulence in the hot, shocked wind. Line styles are as in Figure
7. Left Panel: We show the second order structure function of the turbulent
velocity field, $v_{t}$ (details in text), on a scale $\ell=4\,{\rm pc}$,
which is roughly the energy containing scale of the turbulence across time in
each of the simulations, as a function of the bubble’s radius. We only show
data for where the turbulence at this scale is well resolved: $R_{\rm
bub}>4\,{\rm pc}$ and $R_{\rm bub}>20\,\Delta x$. Right Panel: $\left\langle
v_{\rm out}\right\rangle$ vs. the turbulent velocity, adjusted to match the
predicted turbulent velocity on the resolution scale. We assume a power-law
scaling of the turbulent velocity field of $p=1/2$, which is consistent with
the strucutre functions in our simulations, likely due to the fact that the
turbulent dynamic range is not very well resolved, further details in text. We
show the 1-to-1 relationship between $\left\langle v_{\rm out}\right\rangle$
and the measure of turbulence as a black line.

###  4.7 Relation Between $\left\langle v_{\rm out}\right\rangle$ and
Turbulence

In the TRML literature it has often been shown that the amount of cooling in
the interface layers is straightforwardly dependent upon the amount of
turbulence in the mixing layer (Fielding et al., 2020; Tan et al., 2021; Das &
Gronke, 2023). Indeed, in Lancaster et al. (2021a) and Lancaster et al.
(2021b) the amount of turbulence in the hot, shocked wind gas is causally tied
to the amount of mixing, and therefore dissipation, at the wind bubble’s
interface. While we will discuss nuances in this connection in Section 5.1, it
is worth reporting the relationship between turbulence and dissipation in the
simulations at hand.

We quantify turbulence in our simulations by measuring the second-order
structure function of the non-radial velocity field in the hot, shocked wind.
We do this in exactly the same manner as described in Section 4.5 of Lancaster
et al. (2021b) except that we restrict solely to the $T>10^{6}\,{\rm K}$
gas111111In Lancaster et al. (2021b) all gas with $v_{r}<V_{w}/2$ and
$T>2\times 10^{4}\,{\rm K}$ was used.. To give a brief description of the
process here: we randomly choose 400 points in the shocked wind volume and
calculate the average norm of the difference in the non-radial velocity field
as a function of distance, $\ell$, away from the chosen point. We then use all
400 samplings of the square deviation as a function of distance to compute the
mean over the samples and then take the square-root. This gives us
$v_{t}(\ell)$, the second-order structure function as a function of scale.

In the left panel of Figure 15 we show $v_{t}(4\,{\rm pc})$, which we found to
be approximately the energy containing scale for the turbulence across all of
the simulations, as a function of $R_{\rm bub}$. We see that, while the
turbulent velocity is relatively well converged in the MHD simulations, it
steadily increases with resolution in the HD simulations121212We note that the
turbulent velocity measured here is larger than in Lancaster et al. (2021b)
likely due to both (i) our measurement being restricted to higher temperature
gas and (ii) $V_{w}$ being a factor of $\sim 1.8$ larger in these
simulations.. This is likely due, at least in part, to decreased numerical
viscosity at higher resolution, allowing for a better resolved dynamic range
in the HD simulations. In contrast, small-scale turbulence is likely
suppressed in the MHD simulations due to the stabilization of instabilities
from the magnetic field. If this small-scale suppression is well-resolved in
even the lowest resolution simulations, then we would not expect a better
resolved dynamic range, and therefore increased turbulence, at higher
resolution.

There is an additional effect that can have a significant impact on the
resolution dependence of turbulence in the HD simulations. As we can see in
Figure 7, in the HD simulations we have lower $\alpha_{p}$ at higher
resolution. Through Equation 10, this implies a smaller ratio of $R_{\rm
bub}/R_{\rm f}$, which in turn implies more direct contact between the free-
wind and dense clumps in the background cold gas. This process is directly
observable in Figure 4. This direct contact between the dense gas and the free
wind creates oblique bow shocks which seed vorticity, and hence turbulence, in
the shocked wind material. Indeed, such bow-shocks were posited as a main
source of turbulence in the hot gas in Lancaster et al. (2021b). However, we
still expect turbulent mixing at the interface, seeded by instabilities, even
in the absence of this effect.

It is also interesting to compare the turbulence and dissipation at the
bubble’s surface, quantified by $\left\langle v_{\rm out}\right\rangle$. To
this end we take the energy-containing scale turbulence, which we believe to
be well resolved, and extrapolate it to the expected turbulent velocity at the
resolution scale using $v_{t}(\Delta x)=v_{t}(L)(\Delta x/L)^{1/3}$, the
expected scaling for Kolmogorov turbulence. We can see from the right panel of
Figure 15, where we compare $\left\langle v_{\rm out}\right\rangle$ with this
scaled turbulent velocity, that they have a close to linear relation with one
another, as one would expect from recent work in TRMLs (Fielding et al., 2020;
Lancaster et al., 2021a; Tan et al., 2021; Das & Gronke, 2023). It is
encouraging that this relationship appears more linear in the higher
resolution simulations. Thus, while our current simulations are not fully
resolved, this result suggests that dissipation mediated by turbulence could
be achieved with an enhancement in resolution. The right panel of Figure 2
also indicates that at least for the pressures present at late time, we are
close to being able to resolve $\ell_{\rm cool}$.

##  5 Discussion

###  5.1 Relation to Past Work

A key thesis of this work, laid out in Section 2 and Section 4, is that the
dynamics of WBBs are intimately dependent upon the movement of energy across
their bubble-shell interfaces. In particular, we have emphasized how this
depends both on the multi-scale geometry of these interfaces and the small-
scale dissipation processes that act across them. These two effects acting
together make for a particularly difficult physical problem.

If we ignore geometry, considering one-dimensional problems, the role of
dissipation across such interfaces has long been solved under the further
assumptions of time-steady, isobaric flows (Zel’Dovich & Pikel’Ner, 1969;
Cowie & McKee, 1977; Inoue et al., 2006; Kim & Kim, 2013). This is
unfortunately not the picture that applies to WBBs. The full complications of
the problem were circumvented in Weaver et al. (1977) by ignoring the
contribution of cooling at the interface, deriving a solution that required
all conductive heat flux to be balanced by evaporation of mass from the shell
into the hot, wind bubble.

El-Badry et al. (2019) ameliorated this picture by running spherically
symmetric, radiative, one-dimensional simulations of super-bubble evolution
with both explicit Spitzer conductivity ($\propto T^{5/2}$) and a constant
diffusivity at low temperatures, $\kappa_{\rm mix}$ (meant to account
approximately for turbulent mixing) and showed that a predictable, constant
fraction of the wind energy is lost to cooling within the interface layer,
$\dot{E}_{\mathrm{cool}}\propto\kappa_{\rm mix}^{1/2}$ (see Eq. 36 of El-Badry
et al., 2019). This essentially solves the problem from a ‘purely diffusive’,
one-dimensional picture, albeit by appealing, in an effective way, to a
mechanism of diffusion that is inherently three-dimensional: turbulence. We
note that the problem, of theoretical interest no doubt, that appeals to only
Spitzer-like heat dissipation in the one-dimensional context but includes
cooling in the dynamical WBB picture has not been fully resolved (in the
numerical sense, see Appendix A of El-Badry et al., 2019).

Turbulence spans both sides of the divide we have drawn here between
dissipation and geometry, and it is worth emphasizing which aspects of
turbulence lie on either side in the context of these interfaces. One can
treat turbulence purely as a micro-physical dissipation, as it is treated in
El-Badry et al. (2019) and countless other astrophysical applications (e.g.
Smagorinsky, 1963; Röpke et al., 2007). However, the root of this ‘effective’
micro-physical dissipation is itself geometry: turbulence is stretching-
enhanced diffusion, it only acts to increase the surface area over which true
micro-physical dissipation may act. As we discussed in Section 2.4.2, the
scale on which turbulence effectively becomes micro-physical for these
interfaces, is the cooling length, $\ell_{\rm cool}$, at which point cooling
effectively ‘stands still’ as the turbulent cascade efficiently mixes
material. Above this scale, larger scale turbulence acts to significantly
enhance the surface area over which energy transport takes place, while being
too slow at mixing to be of actual importance in mediating that energy
transport. This is the geometric role that it serves. The distinction between
these two regimes is essentially about the ratio of cooling time to the eddy-
turnover time on the largest scales (quantified by $\xi$ in Fielding et al.
(2020) and the Damköhler number, ${\rm Da}$ in Tan et al. (2021)).

While El-Badry et al. (2019) simply treated $\kappa_{\rm mix}$ as a parameter,
Lancaster et al. (2021a) linked its expected value to predictions from the
theory of Fielding et al. (2020) discussed in Section 2.4.2. Folded into this
theory, and the model for diffusivity given in Lancaster et al. (2021a), is
the geometric role that turbulence plays, accounted for using an additional
factor for the increase in fractal area as a bubble expands. While, in this
case, the energy lost to cooling would no longer necessarily be a constant
fraction of the energy input rate, modifying the WBB expansion solution, the
evaporative flow and temperature of the hot medium could still be computed
following El-Badry et al. (2019).

The true complexity of this problem is realized when we allow resolved
dissipation to act in concert with complicated geometry. Here the geometry
changes the temperature and velocity structure of the flow in three
dimensions, changing the nature of the conductive heat transport. In
particular, while either condensation or evaporation is allowed exclusively in
one-dimensional problems, in this complicated three-dimensional picture, both
could occur simultaneously.

The largest further complication to this picture is the presence of magnetic
fields, which alter the turbulent cascade (Schekochihin, 2022; Grete et al.,
2023), reduce the strength of physical conduction processes (Spitzer & Härm,
1953; Spitzer, 1962; Meinecke et al., 2022), and alter the geometric
enhancement of surfaces, as we have demonstrated here. While the MHD picture
of TRMLs has begun to be investigated in recent numerical works (Das & Gronke,
2023; Zhao & Bai, 2023), the picture of the roles of physical versus numerical
dissipation in these simulations still appears unclear.

In the present work, we have emphasized that the resolution of the relevant
micro-physical diffusive scales (outlined in Section 2.4) is important to
accurately capture the dynamics of WBBs. We have presented a detailed study of
the origin of the numerical diffusivity in these simulations, as well as why
it has a weak scaling with resolution and leads to momentum-driven-like
behavior (summarized in Section 4.6). In particular, the decrease in
dissipation at higher resolution is accomplished by a near equal increase in
surface area, making the net result the same in our MHD simulations. At the
same time, from the results of Section 4.7, it is likely that turbulent
dissipation may play some role in the setting of $\left\langle v_{\rm
out}\right\rangle$ in our simulations.

These results may also provide an explanation for the lack of resolution
dependence in simulations of Kelvin-Helmholz mixing with cooling, even when
$\ell_{\rm cool}$ is not resolved. In particular, for the Kelvin-Helmholz
problem in which there is not a strong mean velocity field flowing into the
interface (unlike the present situation), we would expect the numerical
dissipation velocity to scale as $\Delta x^{1/2}$, as was demonstrated for
one-dimensional problems in Tan et al. (2021). Similarly, as is demonstrated
in Fielding et al. (2020) and here, we would expect the area of the fractal
interface to scale as $\Delta x^{-1/2}$. These effects perfectly balance one
another in this case. We emphasize that the Kelvin-Helmholz problem is quite
different from the present case, where there _is_ a strong flow into the
mixing layer, with a comparable effective velocity of dissipation.

Finally, what does the numerical dependence demonstrated here tell us about
our ability to use numerical simulations to predict the evolution of WBBs? In
Lancaster et al. (2021a, b) we described the root of the “momentum-
driven”-like behaviour present in both Lancaster et al. (2021b) and the
simulations presented here as due to extreme cooling, with $\alpha_{p}$
approaching its lower bound of unity. We argued that the strong cooling was
enabled by turbulence-driven mixing of hot gas into denser material at the
interface with the shell.

What we have demonstrated in this work is that the turbulence, while likely
playing a role in the enhancement of the fractal structures in our
simulations, is not the only, or necessarily dominant, source of diffusivity.
This is clearly evidenced by the resolution dependence of the total momentum
carried by the bubble, as laid out in Figure 6, as well as the resolution
dependence of $\left\langle v_{\rm out}\right\rangle$ (laid out in Section
4.5), which plays a major role in setting the bubble’s interior pressure.

While numerical dissipation is playing a role in setting $\dot{p}_{r}$ in the
current simulations, given the available resolution, this does not change the
overall physical picture for WBB evolution outlined in Lancaster et al.
(2021a). As discussed in Section 2.4 and above, in the real world turbulence
should play a major role in both enhancing the surface area of the interface
and in the transport that leads to intermingling of hot and cool gas at small
scales. The main difference between Lancaster et al. (2021a) and the picture
suggested by this work is that this evolution need not be “momentum-driven” in
the $\alpha_{p}=1$ or even $\alpha_{p}=const.$ sense, but rather could have
$\alpha_{p}$ that varies with time, as suggested by Equation 15.

However, as we showed in Section 2.6, we can use Equation 15 to derive the
momentum scaling that we would expect in WBBs where dissipation at the
interface is mediated through turbulent mixing (described in Equation 25), and
the overall bubble geometry is given by a fractal geometry. If we apply a
Kolmogorov scaling for the turbulence and an excess fractal dimension of
$d=1/2$ (as seems relevant for our hydro simulations) this results in an
$\alpha_{p}$ which is constant in time (“momentum-driven”) with an order unity
value for the parameters relevant to the simulations described here. This
validates and lends further meaning to the physical picture laid out in
Lancaster et al. (2021a).

###  5.2 Prospects for Future Work

While the argument laid out in Section 2.6 resolves the issue of the evolution
of WBBs whose interiors can be thought of as quasi-steady (details in Appendix
A) and whose turbulent properties and fractal geometry are known, several open
questions remain. First among these is when these conditions may apply. Since
the relevant speeds involved in setting the interior and boundary conditions –
the sound speed and flow speeds within the bubble and the turbulent velocities
in the mixing layer – are large compared to the bubble expansion velocity, a
quasi-steady state should be easy to establish.

The question of the nature of energy transport and geometry at turbulent
interfaces has some solid theoretical basis (Fielding et al., 2020; Tan et
al., 2021). However, some open questions remain as to the details of the
resolution (in)dependence of the results obtained from these simulations and
the interpretation of the scaling of the total cooling with the relevant
physical parameters ($t_{\rm cool}$, $v_{t}(L)$, etc.). Understanding the
nature of the resolution dependence of these simulations is, in our view,
essential to a confident theoretical model for dissipation in these systems.
To settle this issue it is essential that future numerical work spatially
resolves the dynamics of these interfaces while including explicit, physically
motivated, diffusive properties and demonstrating that physical diffusion is
dominant over numerical diffusion.

This problem is made even harder by the ubiquitous presence of magnetic fields
in astrophysical fluids, which these interfaces are likely not exempt from.
While this problem has begun to be explored numerically (Das & Gronke, 2023;
Zhao & Bai, 2023), questions of numerical dependence still remain. The
magnetic field not only changes the geometry of these interfaces, as we have
demonstrated here, but also reduces the strength of thermal conduction
(Roberg-Clark et al., 2016; Komarov et al., 2018; Drake et al., 2021; Meinecke
et al., 2022). If magnetic fields act to significantly dampen the turbulent
diffusivity, through altering the cascade at these interfaces by stabilizing
small-scale instabilities, as seems likely from our study as well as Zhao &
Bai (2023) and Das & Gronke (2023), this makes the resolution requirements
even stricter. Numerical studies that resolve the dynamics and dissipative
processes, especially in the presence of magnetic fields, will likely be quite
costly but are absolutely necessary.

Even if we have an understanding of the small-scale diffusive processes
mentioned above, we still need a way of accurately accounting for them in
larger-scale simulations where we have little hope of resolving the
appropriate scales. While one way forward is to potentially include a modified
diffusivity that accounts for this unresolved physics (El-Badry et al., 2019;
Lancaster et al., 2021a), one still needs to ensure that the numerical
diffusivity is less than this new, modified diffusivity. One promising avenue
for ensuring this is to explicitly separate the high and low entropy
components of the fluid in a subgrid manner (Smith et al., 2024; Weinberger &
Hernquist, 2023; Butsky et al., 2024). The validity of this approach is
subject to the details of the coupling between these subgrid components, which
must be tuned to match the understanding gained from resolved interface
simulations.

##  6 Conclusions

We end with a summary of our main conclusions:

  * 1.

We demonstrate that the evolution of a WBB is governed by boundary conditions
at its surface (see Section 2.2 and Appendix A). These boundary conditions are
a combination of geometry, embodied by the area of the bubble’s surface,
$A_{\rm bub}$, and dissipation, embodied by the interface’s ability to absorb
high-entropy gas, $\left\langle v_{\rm out}\right\rangle$, as seen in Equation
15. Increased $A_{\rm bub}$ gives more surface area over which energy can be
dissipated, reducing the degree to which momentum is enhanced above that input
by the wind, $\alpha_{p}$. At the same time decreased $\left\langle v_{\rm
out}\right\rangle$ limits losses at each point on the surface, increasing
$\alpha_{p}$.

  * 2.

We demonstrate that solutions with constant $\alpha_{p}$ are not necessarily
“momentum-dirven” (see Figure 6). Our past work suggested, with evidence from
numerical simulations, that turbulent mixing at the bubble’s interface leads
to efficient enough cooling that these WBBs always enter a $\alpha_{p}\approx
1$ “momentum-driven” limit, Figure 6 demonstrates that this is not actually a
“momentum-driven” evolution, as is also apparent from a fair interpretation of
Equation 15. In Section 4 we showed how this apparent “momentum-driven”-like
behaviour is simply determined by the relative scaling of $A_{\rm bub}$ and
$\left\langle v_{\rm out}\right\rangle$ both in time and with resolution,
$\Delta x$. In particular, differences in these scaling relations between the
HD and MHD simulations, due mostly to differences in $A_{\rm bub}$ from
differing fractal structures, can explain the opposite convergence behaviour
apparent in these simulations, as seen in Figure 7. This is summarized in
Section 4.6.

  * 3.

In analyzing the scaling of $A_{\rm bub}$, Section 4.4 demonstrated, using
three different measurement techniques, how the surfaces of the simulated WBBs
are well described as fractals with an excess fractal dimension of
$d=1/2,\,1/4$ for the HD and MHD simulations respectively. The fractal
dimension of the MHD WBBs is also clearly time-dependent and this is likely
related to the growth of magnetic pressure in the shell of the WBB as the
bubble evolves. This strong magnetic pressure both halts dynamical
instabilities such as Rayleigh-Taylor and Kelvin-Helmholtz but also creates a