# The observable properties of cool winds from galaxies, AGN, and star
clusters - II. 3D models for the multiphase wind of M82

Yuxuan Yuan1 *  |, Mark R. Krumholz1,2 *  |, Crystal L. Martin3
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2601, Australia
2ARC Centre of Excellence for Astronomy in Three Dimensions (ASTRO-3D),
Canberra, ACT 2601, Australia
3Department of Physics, University of California, Santa Barbara, Santa
Barbara, CA 93106, USA
E-mail<EMAIL_ADDRESS>(YY)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Galactic winds are a crucial player in galaxy formation and evolution, but
observations of them have proven extraordinarily difficult to interpret,
leaving large uncertainties even in basic quantities such as mass outflow
rates. Part of this uncertainty arises from the relatively simplistic models
to which complex wind observations are often fit, which inevitably discard
much of the available information. Here we present an analysis of the wind of
the nearby dwarf starburst galaxy M82 using a semi-analytic model that is able
to take advantage of the full three-dimensional information present in
position-position-velocity data cubes measured in the H i 21 cm line, the CO
$J=2\to 1$ line, and the H$\alpha$ line. Our best-fitting model produces
position-dependent spectra in good agreement with the observations, and shows
that the total wind mass flux in the atomic and molecular phases is $\approx
10$ M⊙ yr-1 (corresponding to a mass loading factor of $\approx 2-3$), with
less than a factor of two uncertainty; the mass flux in the warm ionised phase
is more poorly constrained, and may be comparable to or smaller than this. At
least over the few kpc off the plane for which we trace the outflow, it
appears to be a wind escaping the galaxy, rather than a fountain that falls
back. Our fits require that clouds of cool gas entrained into the wind expand
only modestly, suggesting they are magnetically confined. Finally, we
demonstrate that attempts to model the wind using simplifying assumptions such
as instantaneous acceleration and a constant terminal wind speed can yield
significantly erroneous results.

###### keywords:

galaxies: ISM - galaxies: general - methods: analytical

††pubyear: 2021††pagerange: The observable properties of cool winds from
galaxies, AGN, and star clusters - II. 3D models for the multiphase wind of
M82–LABEL:lastpage

##  1 Introduction

In 1963, Lynds and Sandage conducted the first H$\alpha$ observations of the
central region of the starburst galaxy M82 (Lynds & Sandage, 1963) and
discovered that the gas above and below the nucleus has been accelerated to a
velocity much larger than the local escape speed. This paper marked the
discovery of galactic winds, and shows that galaxies are not “isolated
islands” as had been thought in the eighteenth century (Wright, 1750). Modern
observations reveal that galactic winds are ubiquitous across nearly all
starburst galaxies and active galactic nuclei (AGN), at both low and high
redshift (Veilleux et al., 2005, 2020; Heckman & Thompson, 2017), and these
winds typically contain multiple phases within them. The hot phase,
characterised by high temperature ($\sim 10^{8}$ K) and speed, is widely
observed in X-ray emission, and carries most of the energy budget in galaxies
with strong outflows. Cool phase outflows ($\sim 10-10^{4}$ K), on the other
hand, are much denser than the hot component and carry most of the mass.

Galactic winds have a significant impact on the life-cycle of galaxies. They
can expel gas and quench star formation in the galactic nuclei within tens of
million years, making them a candidate to explain galaxies’ low star formation
efficiency and baryon fractions (Hopkins et al., 2011; Alatalo et al., 2011).
They carry mass, metals and dust from the ISM to the circumgalactic medium
(CGM) or even intergalactic medium (IGM; Nelson et al. 2018; Tumlinson et al.
2011; Muratov et al. 2017), making the IGM the largest reservoir of baryons
and metals across the universe (Oppenheimer & Davé, 2006). This provides a
possible solution to the long-standing “missing baryons problem”. However,
some of the mass entrained into winds may also form a galactic fountain,
recycling back to galaxies and triggering further star formation (Bregman,
1980; Melioli et al., 2015; Marasco & Fraternali, 2017; Werk et al., 2019; Li
& Tonnesen, 2020).

In the first widely accepted galactic wind model, proposed by Chevalier &
Clegg (1985), the explosion energy produced by supernovae (SNe) drives a hot,
adiabatically-expanding, fast outflow. This model describes the hot phase, but
cannot explain the cool winds that are widely observed. Theoretical models for
these began to appear several decades later, and there is still significant
uncertainty about the origin of the cool phase. In one scenario, cool gas in
the outflow launching region is rapidly shredded and shock heated into the hot
winds, but this initial hot outflow then undergoes radiative cooling beyond a
cooling radius of several kpc, resulting in the observed cool phase (Thompson
et al., 2016; Schneider et al., 2018). A second possibility is that magnetic
fields (e.g., McCourt et al., 2015; Banda-Barragán et al., 2016, 2018a),
radiative cooling (e.g., Scannapieco & Brüggen, 2015; Gronke & Oh, 2018, 2020;
Schneider et al., 2020; Kim et al., 2020a), or a combination of both (e.g.,
Cottle et al., 2020; Li et al., 2020; Kim et al., 2020b; Banda-Barragán et
al., 2021) allow cold gas to survive being shocked and entrained by a hot
outflow. Yet a third scenario is that cool gas clouds are accelerated via
radiation pressure (e.g., Krumholz & Thompson, 2013; Coker et al., 2013;
Thompson et al., 2015; Crocker et al., 2018) or cosmic ray pressure (e.g. Mao
& Ostriker, 2018; Brüggen & Scannapieco, 2020; Crocker et al., 2021b, a) from
the central source, and are therefore accelerated gently and not shocked at
all. In the latter two scenarios cool gas clouds maintain their phase
structure along their trajectories, while in the first they disappear and re-
form later.

Although numerical simulations provide the most accurate descriptions of
galactic winds, they can survey at most a very limited portion of parameter
space. This makes them unsuited to the task of extracting information from the
observations that are “sampled” from a potentially much larger portion of
parameter space. Hence analysis of the kinematic and thermodynamic structure
of observed winds requires using analytic or semi-analytic models, which are
capable of generating synthetic observations for comparison to observations
with high enough computational efficiency to allow parameter fitting. However,
these kinds of models often resort to highly simplified or heuristic
prescriptions to minimize computational cost. For example, many models (e.g.,
Steidel et al., 2010; Prochaska et al., 2011; Scarlata & Panagia, 2015;
Lochhaas et al., 2018) describe winds as fully-filled spherically symmetric
expanding shells, and adopt density, velocity and covering fraction
prescriptions that are simple functions of radius from the center. Carr et al.
(2018) adopt a more realistic biconical geometry, but still rely on a
simplified kinematic structure inherited from Scarlata & Panagia (2015).
Although these models are able to roughly reproduce stacked, spatially-
unresolved spectra, they inevitably discard much of the information hidden
within the spatially-resolved spectra available for nearby sources. Moreover,
these models solely focus on a single phase of galactic wind, and hence cannot
produce a unified picture of a multi-phase outflow. Prediction of the spatial
and velocity structure of a multi-phase wind with enough complexity to allow
useful model fitting, but with enough simplicity for that fitting to be
computationally tractable, remains a great challenge for modern observational
studies of galactic winds.

This challenge led Krumholz et al. (2017, hereafter KTOM) to propose a semi-
analytical model that gives a reasonable balance between physical complexity
and numerical efficiency. The KTOM model can generate spatially resolved
synthetic spectra based on a physical model for wind acceleration that depends
on a few physical parameters and prescriptions, suitable for comparison to the
resolved observations of wind. Unlike earlier semi-analytic models, it allows
a wide range of wind geometries and driving mechanisms, and can model emission
a variety of phases. This model is implemented as an extension of the open-
source code Derive the Energetics and SPectra of Optically Thick Interstellar
Clouds (despotic; Krumholz 2014).

In this paper we provide a first test of this model against real observations,
by using it to constrain the physical properties of the cool phases of the
outflows in the dwarf starburst galaxy M82. We aim to extract much more
detailed information about the properties of the wind than would have been
possible with previous methods, taking full advantage of the three-dimensional
position-position-velocity data cubes in multiple tracers that are available
for one of the most well-observed galactic winds. We derive basic parameters
such as mass outflow rate and wind geometry for each of the phases
independently, and highlight the relationships between them.

The structure of the paper is as follows: In Section 2, we briefly review the
KTOM model, the data to which we will apply it, and the methods we use to fit
to the observed data. We then present the results from measurements in Section
3. We next discuss our results in Section 4 and finally summarise our
conclusions in Section 5.

##  2 Methods

We first introduce the observational data to which we will fit in Section 2.1.
In Section 2.2 we summarise the KTOM model and explain how varying its
parameters affects predicted spectra. We describe our fitting method and
compare different combinations of model choices in Section 2.3.

###  2.1 Observations of the M82 outflow

Messier 82 (M82) hosts one of the most well-studied galactic winds. M82 is an
edge-on (inclination angle $\sim 80\degree$) starburst galaxy (O’Connell &
Mangano, 1978; O’Connell et al., 1995), located only 3.6 Mpc away (1 arcsec
$\sim$ 17.5 pc) (Freedman et al., 1994), making it a promising candidate for
the study of outflow around the minor axis. Hubble Space Telescope (HST)
imaging reveals that M82 hosts very young massive clusters, concentrated in
its central 500 pc (30 arcsec) nucleus. M82 is characterised by its intense
central starburst, triggered by a tidal interaction with its neighbour M81.

The M82 wind has been observed across the spectrum, from a hot phase observed
in soft X-rays ($\sim 10^{7}$ K) (Watson et al., 1984; Bregman et al., 1995;
Strickland et al., 1997; Lopez et al., 2020), to warm ionised gas traced by
H$\alpha$ ($\sim 10^{4}$ K) (McKeith et al., 1995; Martin, 1998; Westmoquette
et al., 2009), to a cool neutral phase seen in H i 21cm ($\sim 5000$ K)
(Martini et al., 2018) and a cold molecular phase traced by CO ($\lesssim 100$
K) (Walter et al., 2002; Leroy et al., 2015). Recent observations suggest that
the latter two phases exist largely on the edge of a central hot outflow, with
H$\alpha$ lying at the interface between the hot outflow and surrounding cold
gas.

Figure 1: Velocity-integrated intensity (top row) and first and second moment
maps (middle and bottom row) for the observations of the CO $2\to 1$ (Leroy et
al., 2015, left) and H i 21 cm lines (Martini et al., 2018, center and right);
the center and right columns both show the same data, and differ only in that
the central column is zoomed in so that the field of view shown matches that
for the CO in the left column, while the right column shows the full extent of
the 21 cm data. We have oriented our coordinates so that the major and minor
axis of the galaxy defined in Martini et al. (2018) lie along the $x$ and $y$
axes, respectively; arrows in the top left panel show the orientation on the
sky. We have also shifted the first moment so that a value of zero corresponds
to the systemic velocity of 211 km s-1. The white region in the centre of 21
cm map is masked because it is dominated by galaxy emission; other regions
shown in white indicate non-detections. Orange squares indicate the apertures
over which we compute average spectra for use in our fitting procedure; see
main text for details.

####  2.1.1 Neutral phases: H i and CO

In this paper we use the CO $J=2\to 1$ data obtained by Leroy et al. (2015) to
trace the cold molecular phase of the outflow, and the H i 21 cm data obtained
by Martini et al. (2018) to trace the cool neutral phase. The first data set
consists of a position-position-velocity (PPV) cube that covers a $\approx
2.5\times 2.5$ kpc region around the galaxy and the second covers an $\approx
8\times 8$ kpc region. We refer readers to the original papers for full
details on the observations. To show the overall kinematic structure of the
cool wind, we plot the velocity-integrated intensity and first and second
moment maps of these two data sets in Figure 1; we show two views of the H i
data, one zoomed in on the field that overlaps the CO data, and one zoomed out
to show the full extent of the 21 cm data. We have oriented our coordinates so
that the major and minor axes of the galaxy, as defined in Martini et al.
(2018), lie along the $x$ and $y$ axes, respectively, and we set $v=0$ to
correspond to the systemic velocity of 211 km s-1. We can clearly see evidence
for rotation both in the midplane and in the outflow region of the first
moment map. We also see that the H i outflow is much more extended than the
molecular outflow, which suggests they have different spatial distributions
and possibly different mass fluxes.

We next select representative regions over which to compute averaged spectra
to which we will compare our model. Our motivation for doing so is three-fold:
first, averaging obviously increases the signal to noise ratio, yielding a
cleaner fit. Second, while the semi-analytic model we use to generate
predicted spectra as we vary our parameters is quite fast, it is not fast
enough for it to be practical for us fit the full PPV cube at its native
spatial resolution; we must reduce the size of the data for computational
reasons. Third, while the KTOM model we use is significantly more complex than
previous galactic wind models, it is nonetheless obviously an
oversimplification compared to reality. We therefore wish to pick regions that
generally capture the rich information in the observation as a whole, but
average out the data enough that our fitting is not overly biased by “small
scale” spectroscopic features at certain positions, and gives a reasonable fit
to the entire observation.

In accord with this principle, we extract spectra at a broad range of
positions within two rectangular regions around the minor axis of the disc. We
define two extraction regions; the first, which we apply to the $2.5\times
2.5$ kpc field within which the H i and CO fields of view overlap, starts at
projected distances 0.8 kpc off the disc ($|y|>0.8$ kpc); this geometric cut
is to ensure that the emission to which we are fitting is dominated by the
outflow and not the disc. The second extends from 2.5 kpc to 8 kpc, and covers
the region within which we have only H i data. In both regions, we stack
spectra within square apertures with a size of $\pm 10^{\prime\prime}$ at each
position, and we omit apertures where emission is undetected over any part of
the aperture. We show the square regions over which we average in Figure 1,
with the inner extraction region shown in the left and central columns, and
the outer region in the right column. This procedure yields 30 (31)
independent spectra for the H i data set at central (full) region, and 18 for
the CO (central region only).

####  2.1.2 Warm ionised phase: H$\alpha$

Figure 2: Integrated intensity in the H$\alpha$ line (Kennicutt et al., 2008).
We have oriented the image to match the major and minor axes of the galaxy
defined in Martini et al. (2018), the $x$ and $y$ axes, respectively. Black
squares indicate the apertures over which we compute average spectra for use
in our fitting procedure (Martin, 1998); see main text for details.

We used the M82-2 spectrum shown in Figure 1b of Martin (1998) to measure the
kinematics of the warm ionised phase. The 3$\aas@@fstack{\prime}$5 long slit
was oriented along the outflow axis. Very high spectral resolution of 11.5
km/s was obtained using an echelle grating, and a narrowband filter blocked
all orders except the one covering the H$\alpha$ line. This resolution reveals
double peaked line profiles along the length of the slit, consistent with
emission from the surface of an expanding cone (see (Heckman et al., 1990)).
Figure 2d of Martin (1998) shows the velocity separation as a function of
position along the slit.

We extracted a spectrum for each lobe at a location near the position of
maximum velocity separation. The black squares shown in Figure 2 mark these
locations. These echellograms have an outdated format, so it was necessary to
convert them to the FITS standard.111The code used to convert old IRAF images
to FITS can be downloaded from https://github.com/mjfitzpatrick/imh2fits. The
spectra were not flux calibrated, so we scaled their intensity to the surface
brightness measured in the same aperture of H$\alpha$ image from Kennicutt et
al. (2008).

###  2.2 Galactic wind line radiative transfer model

We next summarise the most relevant features of the KTOM wind model to which
we will fit the data, leaving full details to the original paper and Thompson
& Krumholz (2016), from which it follows. For convenience we collect all of
the parameters we introduce in this section in Table 1.

Parameter | Description | Possible range / values
---|---|---
Continuous parameters
$\log(\dot{M}/\mathrm{M}_{\odot}\;\mathrm{yr}^{-1})$ | Wind mass flux | $(-1,3)$
$\log\mathcal{M}$ | Mach number in the wind launching region | $(0,4)$
$\phi$ | Inclination of outflow central axis on the plane of the sky | $(-\pi/2,\pi/2)$
$\theta_{\rm in}$ | Inner opening angle of wind cone | $(0,\pi/2)$
$\theta_{\rm out}$ | Outer opening angle of wind cone | $(\theta_{\rm in},\pi/2)$
$\tau_{0}$ | Mean optical depth of wind launching region (radiatively-driven winds only) | $(1/\Gamma,300)$
$u_{h}$ | Ratio of hot gas speed to escape speed (hot gas-driven winds only) | $(1,50)$
Discrete parameters
$y$ | Scaling of cloud size with radius: constant area, intermediate, constant solid angle | $y=1,a,a^{2}$
$m$ | Scaling of gravitational potential with radius: point or isothermal potential | $m=1,a$
$f_{p}$ | Wind driving mechanism: ideal, radiation pressure, hot gas | $f_{p}=1,1-\exp(-e^{x}\tau_{0}/y),(1-u/u_{h})^{2}$
Table 1: Summary of parameters for the KTOM model. The top part of the table
shows continuously-variable parameters, while the bottom shows discrete
parameters. Full details on the definitions of the parameters are given in
Section 2.2.

####  2.2.1 Physical model

KTOM describe a momentum-driven wind launched from a spherical region with
enclosed mass $M_{0}$, radius $r_{0}$, and escape speed
$v_{0}=\sqrt{2GM_{0}/r_{0}}$; the gas in this region has a mean surface
density $\Sigma_{0}$, and is driven outward by a mechanism that injects
momentum at a constant rate $\dot{p}$. The ratio of $\dot{p}$ to the momentum
flux provided by gravity defines the generalised Eddington ratio

$\Gamma=\frac{\dot{p}}{4\pi GM_{0}\Sigma_{0}}.$ (1)

The KTOM model applies to a medium for which $\Gamma<1$, and thus a uniform
medium would not be ejected into a wind. However, for a turbulent medium with
Mach number $\mathcal{M}$, the surface density will be lognormally
distributed, with a width that depends on $\mathcal{M}$, and material for
which the surface density relative to the mean $x=\log\Sigma/\Sigma_{0}$ is
sufficiently small will be driven outward.

Mass that enters the wind begins at radius $r_{0}$ with velocity $v=0$, but
then accelerates under the action of the driving force. The resulting velocity
evolves as

$2u\frac{du}{da}=\frac{1}{a^{2}}\left(y\Gamma e^{-x}f_{p}-m\right),$ (2)

where $u=v/v_{0}$, $a=r/r_{0}$, and $y$, $m$, and $f_{p}$ are functions that
describe, respectively, how clouds of material entrained into the wind change
their area as they move outward, how the gravitational potential varies with
radius, and how the momentum deposition mechanism operates. We refer to the
solution $u(a)$ to Equation 2 as the wind acceleration law; KTOM provide
acceleration laws for several different choices of these functions, all of
which we will consider here:

  * •

Cloud expansion: we consider $y=1$, $y=a$, and $y=a^{2}$; the first of these
possibilities corresponds to clouds that maintain constant area as they are
accelerated, the third to clouds that maintain constant solid angle, and the
second to an intermediate level of expansion.

  * •

Gravitational potential: the acceleration law depends on the shape of the
gravitational potential experienced by the wind, parameterised by the
effective enclosed mass as a function of radius $m$ in Equation 2. KTOM
consider both $m=1$ (constant enclosed mass), corresponding to a point
potential where gravity is dominated by mass in the acceleration region, and
$m=a$ (enclosed mass increasingly linearly with radius), corresponding to an
isothermal potential.

  * •

Driving mechanism: we consider three ways a wind can be driven. The first is
an idealised mechanism that simply deposits momentum at a fixed rate
proportional to the cloud cross-sectional area; for this case $f_{p}=1$. The
second is a wind driven by direct radiation pressure, in which case KTOM show
that $f_{p}=1-\exp(-e^{x}\tau_{0}/y)$, where $\tau_{0}=\kappa\Sigma_{0}$ and
$\kappa$ is the flux-mean specific opacity of the gas to radiation. The final
possibility we consider is a cool wind driven by the momentum deposited by a
hot outflow travelling at (dimensionless) speed $u_{h}$, for which
$f_{p}=(1-u/u_{h})^{2}$.

Each of these choices affects the shape of the spectra produced by the wind.
For example, as the cloud expansion law moves from constant area to constant
solid angle, winds intercept more momentum from the central engine and thus
accelerate more rapidly, pushing more emission to higher velocities; moving
from a point to an isothermal potential, thereby decelerating gas more
strongly via gravity, has the opposite effect. Different driving mechanisms
lead to different shapes for line wings. We explore these differences in more
detail below.

In addition to the acceleration law and the Mach number $\mathcal{M}$ of the
medium from which the wind is launched, wind spectra are affected by two other
primary factors: the wind mass flux and geometry. The mass flux $\dot{M}$ acts
primarily to set the normalisation of the predicted spectra, although the
relationship is somewhat more complex than this due to radiative transfer
effects. The mass flux is determined by the surface density $\Sigma_{0}$ and
by the value of $\Gamma$, which determines the minimum value of $x$ for which
the right-hand side of Equation 2 is positive. We refer to this value as
$x_{\rm crit}$, and its value is given by $x_{\rm crit}=\ln\Gamma$ for ideal
or hot gas-driven winds, and implicitly by the solution to $\Gamma e^{-x_{\rm
crit}}[1-\exp(-e^{x_{\rm crit}}\tau_{0})]=1$ for radiatively-driven winds.
Following KTOM, we choose to eliminate the dependence of $\dot{M}$ on
$\Sigma_{0}$ by instead writing the wind mass flux in terms of the star
formation rate $\dot{M}_{*}$, as

$\dot{M}=f_{A}\frac{\zeta_{M}}{\epsilon_{\rm ff}}\dot{M}_{*},$ (3)

where $f_{A}$ is the fraction of the area over which the wind can escape
(i.e., where it is not blocked by dense material such as the remainder of
M82’s disc – see below),

$\zeta_{M}=\frac{1}{2}\left[1-\mbox{erf}\left(\frac{-2x_{\rm
crit}+\sigma_{x}^{2}}{2\sqrt{2}\sigma_{x}}\right)\right]$ (4)

is the fraction of mass for which the right-hand side of Equation 2 is
positive, $\sigma_{x}$ is the width of the density distribution (determined by
$\mathcal{M}$ – see KTOM, equations 11 and 12), and $\epsilon_{\rm ff}$ is the
star formation rate per free-fall time. We adopt $\epsilon_{\rm ff}=0.01$ in
this work, based on extensive observations showing that all star-forming
systems fall near this value (see the review by Krumholz et al. 2019 and
references therein).

The final physical parameter for the wind is its geometry. Observed galactic
winds are often biconical (e.g., Shopbell & Bland-Hawthorn, 1998), with the
outflow blocked from escaping the galaxy in the disc plane, and blowing out in
directions normal to the plane. Moreover, cooler phases are often inferred to
lie on the edge of a central region dominated by hot gas. For this reason, we
consider winds with a “cone sheath” geometry, whereby the wind is confined
between two cones with outer opening angle $\theta_{\rm out}\in(0,\pi/2]$ and
inner opening angle $\theta_{\rm in}\in[0,\theta_{\rm out})$; note that we
allow $\theta_{\rm in}=0$ and $\theta_{\rm out}=\pi/2$, so we allow the
possibility for a wind that is a fully filled cone or even a sphere covering a
full $4\pi$ sr. The inner and outer opening angle determine the fraction of
the unit sphere subtended by the wind: $f_{A}=\cos\theta_{\rm
in}-\cos\theta_{\rm out}$. The central axis of the wind cone is inclined at an
angle $\phi\in(-\pi/2,\pi/2)$ relative to the plane of the sky, where $\phi=0$
corresponds to the central axis of the outflow cone lying exactly in the plane
of the sky. The values of these parameters also manifest in the spectral
shape: if $\theta_{\rm in}$ is close to $\theta_{\rm out}$, this implies that
the wind is confined to a narrow sheath, which manifests as spectra that are
more narrowly peaked, since the line of sight passes through a smaller range
of radii. Similarly, if $\phi=0$ then the near and far sides of the wind are
symmetric relative to the line of sight, implying symmetry in the red and blue
sides of the spectrum; tipping the wind cone, $\phi\neq 0$, produces an
asymmetry between the radius at which our line of sight passes through the
near and far sides of the wind cone, which in turn induces an asymmetry
between the red and blue sides of the spectrum.

Thus the physical model of the wind is parameterised by a total of three
dimensional quantities – $r_{0}$, $v_{0}$, and $\dot{M}_{*}$ – and three
dimensionless functions – $y$, $m$, and $f_{p}$ – and five or six
dimensionless numbers $\Gamma$, $\mathcal{M}$, $\theta_{\rm out}$,
$\theta_{\rm in}$, $\phi$, and, depending on the choice of driving mechanism,
either $\tau_{0}$ or $u_{h}$. Most of these we will fit, as described below,
but some can be fixed by other observations. We take $r_{0}=250$ pc, the
observed radius of the star-bursting centre of M82 (Kennicutt, 1998), and the
escape speed from this region to be $v_{0}=170$ km s-1 (Greco et al., 2012),
corresponding to a dynamical mass $M_{0}=8.2\times 10^{8}$ M⊙. M82’s star
formation rate is $\dot{M}_{*}=4.1$ M⊙ yr-1, which we derive by correcting the
value inferred by Kennicutt (1998) from their assumed Salpeter (1955) IMF to a
Chabrier (2005) IMF.

In principle the shape of the gravitational potential is also measurable;
however, in practice this has proven difficult. Sofue et al. (1992) find that
the gas in the central 2 kpc of M82 is well-fit by Keplerian rotation,
implying a point-like potential, while Greco et al. (2012) find that the
stellar rotation curve suggests a nearly-isothermal potential at radii
$\lesssim 1$ kpc, giving way a point-like potential (constant enclosed mass –
see their Figure 5) beyond $\approx 2-4$ kpc. The differences between the gas
and stellar velocity profiles are likely a result of the differential response
of these components to either the tidal perturbation from M81 or the starburst
itself. Moreover, both of these measurements are in-plane, and may have
limited applicability to the out-of-plane potential, which is what matters for
our purposes. The equipotential surfaces are almost certainly flattened, so
that, to the extent that the potential follows Greco et al.’s model of a
transition from isothermal to point-like at $\approx 2-4$ kpc in the plane,
this transition should occur closer to the galactic centre in the vertical
direction; however, since we do not know how flattened the potential is, we do
not know how much closer. Consequently, we leave the potential shape as an
unknown parameter to be fit.

####  2.2.2 Chemical and excitation model

In addition to the physical model, we require a model for the abundance and
excitation state of the atoms and molecules that produce the observable
emission. We have three different emission tracers – H i, CO, and H$\alpha$ –
which we discuss in turn.

We assume that H i emission comes from material where all the hydrogen is
atomic, and thus the mass per emitting atom is $\mu_{\rm
H~{}\textsc{i}}=2.3\times 10^{-24}$ g, as expected for the standard cosmic
ratio of H to He. Given the low Einstein $A$ of the 21 cm transition, and the
fact that H i is efficiently thermalised by the Wouthuysen-Field effect
(Wouthuysen, 1952; Field, 1959), we assume that the emitting atoms are in
local thermodynamic equilibrium (LTE). We adopt a temperature of 5000 K,
though this choice has no practical effect as long as the temperature is much
larger than the $\ll 1$ K temperature difference between the two states.
Finally, the KTOM model makes a distinction between “correlated” and
“uncorrelated” winds; the former corresponds to cases where the each of the
clouds entrained into the wind is individually optically thick (even if the
wind as a whole might be optically thin because its filling factor is small),
while the the latter to cases where the individual wind clouds are optically
thin (even if the wind as a whole might be optically thick due to the
superposition of many clouds along the line of sight). One can analogise the
correlated case to a wind made up of opaque dust grains, and the uncorrelated
case to a wind made up of transparent water droplets. Given the typically very
low optical depth of warm H i, we adopt the uncorrelated case for the 21 cm
emission.

For CO we adopt a mass per CO molecule $\mu_{\rm CO}=2.1\times 10^{-20}$ g,
corresponding to a CO abundance of $1.1\times 10^{-4}$ CO molecules per H
atom. We also assume LTE and adopt a temperature of $T=50$ K. These choices
are somewhat more uncertain and significant than the equivalent ones for H i,
and we discuss their implications and the uncertainties they induce below. We
also assume that CO is in the correlated case, given the typically high
optical depth of CO-emitting clouds. It is worth noting that the KTOM model
does not assume a fixed value of $\alpha_{\rm CO}$; instead, the value of
$\alpha_{\rm CO}$ is computed self-consistently from the wind acceleration law
using the large velocity gradient approximation.

Finally we come to H$\alpha$, which is the most complex case because it can be
produced in two distinct phases. One is in fully ionised gas, and the other is
in predominantly neutral gas that is subject to an ionising flux, and
recombines at the ionised-neutral interface. The KTOM model is intended to
treat primarily the former, for this component we have mass per H+ ion
$\mu_{\rm H^{+}}=2.3\times 10^{-24}$ g, which gives an energy emission rate
per unit volume $\Lambda n_{\rm H}^{2}$, where $n_{\rm H}$ is the number
density of H nuclei, and the coefficient $\Lambda=3.9\times 10^{-25}$ erg cm3
s-1. We further assume that the H$\alpha$ line is optically thin. Because the
emission rate varies as the square of density, calculation of the H$\alpha$
emission requires that we adopt a clumping factor $c_{\rho}=\langle n_{\rm
H}^{2}\rangle/\langle n_{\rm H}\rangle^{2}$, where the angle brackets indicate
a volume average.222In the optical literature it is more common to define the
filling factor $f=\langle n\rangle^{2}/\langle n^{2}\rangle=1/c_{\rho}$. The
two are obviously equivalent. The emissivity is proportional to the value of
$c_{\rho}$. We discuss our treatment of $c_{\rho}$ below.

For the second component, the emission will be proportional to the ionising
flux to which the neutral gas is subjected. We have no direct knowledge of the
ionising radiation field present in the wind of M82, so we make the simplest
possible assumption that it is constant, in which case the emissivity at a
given position and velocity is simply proportional to atomic hydrogen mass,
and thus to the H i 21 cm signal, in that voxel. The scaling between this
signal and the emitted H$\alpha$ is unknown, and thus becomes an additional
parameter to be fix.

Given these choices, we can compute a model-predicted spectrum of the wind in
H i 21 cm and CO $J=2\to 1$ emission at an arbitrary position using equations
56 and 78 of KTOM, which are implemented in despotic (Krumholz, 2014). We use
the same procedure for H$\alpha$, except that we add the emission predicted by
despotic to an emission term that is proportional to the 21 cm signal at the
corresponding position and velocity.

###  2.3 Fitting method

We next outline our strategy for fitting the model to the data. We consider a
series of wind models described by all possible combinations of driving
mechanism, potential and expansion law. For each combination, we optimise the
set of physical parameters under a Bayesian framework. We then identify the
combination that gives the best fit between model and data. We carry out this
exercise independently for the Northern and Southern sides of the wind, since
the wind is observed to be asymmetric between these two hemispheres (Leroy et
al., 2015; Martini et al., 2018).

####  2.3.1 Parameter fitting

For the cases of H i and CO, we fit five parameters, $\phi$,
$\theta_{\mathrm{in}}$, $\theta_{\mathrm{out}}$, $\mathcal{M}$, and $\dot{M}$,
for ideal momentum driven winds, with an additional parameter $\tau_{0}$ for
radiation driven and $u_{h}$ for hot gas driven winds, respectively; note that
we do not need to treat $\Gamma$ or $x_{\rm crit}$ as independent parameters,
because they can be deduced from the other parameters via Equation 3 and
Equation 4. The case of H$\alpha$ is handled slightly differently, and we
defer a discussion of it to Section 2.3.3. We adopt uniform priors for
$\sin(\phi)$ ($\phi\in(-\pi/2,\pi/2)$), $\theta_{\mathrm{in}}\in[0,\pi/2)$,
$\theta_{\mathrm{out}}\in(\theta_{\rm in},\pi/2)$,
$\log\mathcal{M}\in(0,4)$333In principle we could adopt an informative prior
for $\mathcal{M}$, since the velocity dispersions and approximate temperatures
of the various gas phases are measureable quantities. However, we shall see
below that all other quantities are essentially uncorrelated with
$\mathcal{M}$, so there is little point to doing so., and $\log\dot{M}/{\rm
M}_{\odot}\,{\rm yr}^{-1}\in(-1,3)$. We note that the range of mass outflow
rate is broad enough to encompass essentially all previously reported results
for mass outflow rates (e.g., Muratov et al., 2015; Fluetsch et al., 2019;
Roberts-Borsani & Saintonge, 2019; Roberts-Borsani et al., 2020; Roberts-
Borsani, 2020). We also adopt flat priors on $\tau_{0}$ subject to the
constraint $\tau_{0}>1/\Gamma$ (since a wind is only launched at all if this
condition is met), and flat priors on $u_{h}\in(1,50)$, corresponding to hot
wind speeds from $170-8,500$ km s-1.

In addition to these model parameters, we allow the zero of velocity to vary
with position. The reason we do so is that there is considerable evidence that
the outflow from M82 is rotating (Leroy et al., 2015), whereas our KTOM model
does not include rotation. We approximate that the effect of rotation is
simply to shift the zero point of the spectrum as a function of position, and
include this effect by leaving the zero point at each position as a free
parameter to be optimised.

We compute the likelihood of model given the observations on a spectrum-by-
spectrum basis, so the total likelihood function is given by the product of
the likelihoods over all sky positions and velocity bins,

$\mathcal{L}(f_{\mathrm{obs}}|\mathbf{Q})=-\frac{1}{2}\sum_{i}\sum_{j}\frac{\left[f_{i}(v_{j}-v_{0,i};\mathbf{Q})-f_{i,\mathrm{obs}}(v_{j})\right]^{2}}{\sigma_{ij}^{2}}$
(5)

where $\mathbf{Q}$ is the vector of parameters being fit ($\phi$, $\theta_{\rm
in}$, $\theta_{\rm out}$, $\mathcal{M}$, $\dot{M}$, $\tau_{0}$ or $u_{\rm h}$
for radiatively-driven or hot gas-driven models), $f_{i}$ and $f_{i,\rm obs}$
represent the predicted and observed spectra at the $i$th position (i.e., at
the $i$th box shown in Figure 1), the spectra are measured at velocities
$v_{j}$, and $\sigma_{ij}$ is the observational uncertainty at position $i$
and velocity bin $j$. At each position $i$, we evaluate the sum using the
velocity shift $v_{0,i}$ that maximises the likelihood at that position. The
posterior PDF can then be calculated via Bayes’ theorem

$P(\mathbf{Q}|f_{\mathrm{obs}})\propto\mathcal{L}(f_{\mathrm{obs}}|\mathbf{Q})P_{\mathrm{prior}}(\mathbf{Q})$
(6)

We calculate the posterior probability distribution function, from which we
determine the best-fitting parameters with their uncertainties, using the
affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler emcee
(Foreman-Mackey et al., 2013). We repeat this calculation for each possible
combination of driving mechanism, potential and expansion law, using 40 MCMC
walkers, each iterating over 500 steps to sample the posterior distribution of
parameters. In general we find that the distribution of walkers needs $\sim$
50 iterations to stabilise, so we burn in the first 100 iterations and derive
the posterior PDFs from the rest.

####  2.3.2 Which combination of driving mechanism, potential and expansion
law?

Having optimised physical parameters for all combinations of the driving
mechanism, potential, and expansion law, denoted as (dm, p, ex), we are now in
a position to distinguish which combination gives the best fit to the observed
data. In order to decide this, for each combination (dm, p, ex), we find the
largest value of the likelihood function returned by any of the MCMC sample
points, denoted as $\hat{\mathcal{L}}_{\mathrm{(dm,p,ex)}}$, and compute the
corresponding Akaike information criteria (AIC) (see Sharma 2017)

$\mathrm{AIC}_{\mathrm{(dm,p,ex)}}=2k-2\ln\hat{\mathcal{L}}_{\mathrm{(dm,p,ex)}}$
(7)

where $k=5,$ 6 and 6 are the number of parameters for ideal, radiation and hot
gas driven winds, respectively: 3 parameters to describe the cone sheath
geometry, 1 parameter for the Mach number, 1 parameter to describe the mass-
outflow rate, and 1 extra parameter to describe the optical depth or hot gas
velocity for the radiation-driven or hot gas cases, respectively. The
corresponding Akaike weight for each combination then is

$\displaystyle w(\mathrm{dm,p,ex})$
$\displaystyle=\frac{\mathrm{e}^{-\Delta_{\mathrm{(dm,p,ex)}}/2}}{\sum\limits_{\rm(dm,p,ex)}\mathrm{e}^{-\Delta_{\rm(dm,p,ex)}/2}}$
(8) $\displaystyle\Delta_{\mathrm{(dm,p,ex)}}$
$\displaystyle=\mathrm{AIC}_{\mathrm{(dm,p,ex)}}-\min\left(\mathrm{AIC}_{\mathrm{(dm,p,ex)}}\right)$

which gives the relative probability for each combination of driving
mechanism, potential, and expansion law.

####  2.3.3 Modifications for H$\alpha$

As discussed in Section 2.2.2, H$\alpha$ is substantially more complex than H
i or CO, both because it potentially arises from two distinct phases – fully
ionised gas and the surfaces of neutral clouds – and because the former of
these has an emissivity that scales as the square of the local volume density,
and thus involves a clumping factor $c_{\rho}$. Our method for handling the
emission from H i cloud surfaces is straightforward; as discussed in Section
2.2.2, we assume that this emission is proportional to the H i mass at every
PPV voxel, and the constant of proportionality between the 21 cm signal and
the H$\alpha$ signal simply becomes another parameter to be fit by the MCMC.

The clumping factor $c_{\rho}$ requires more care; since it simply applies an
overall scaling factor to the signal, it is essentially degenerate with the
mass outflow rate. For this reason, if we were to treat $c_{\rho}$ as a fit
parameter, our derived mass outflow rate would depend sensitively on our prior
for $c_{\rho}$. For this reason, we choose to carry separate fits using a
range of fixed $c_{\rho}$ values that span the plausible range, and allow us
to determine a corresponding plausible range of mass outflow rates. We discuss
the results of this experiment below.

##  3 Results

In this section, we fit our model to each of the observational data sets
described in Section 2.1, using the fitting method described in Section 2.3.
We examine the H i data first in Section 3.1, since they are the most
straightforward conceptually, then move on to the CO data in Section 3.2 and
the H$\alpha$ data in Section 3.3.

H i 21 cm (Martini et al., 2018)
---
Models |  | Best fit parameters
Driver | Potential | Expansion | $w$ | $\phi$ | $\theta_{\mathrm{in}}$ | $\theta_{\mathrm{out}}$ | $\log\dot{M}$ | $\log\mathcal{M}$ | $\tau_{0}$ | $u_{h}$
Northern hemisphere
Ideal | Point | Area | 0.52 | $\mathbf{4.73_{-5.10}^{+4.21}}$ | $\mathbf{30.08_{-2.64}^{+3.54}}$ | $\mathbf{85.96_{-6.15}^{+2.76}}$ | $\mathbf{0.61_{-0.11}^{+0.14}}$ | $\mathbf{\gtrsim 1.57}$ | - | -
Intermediate | 0 | $-37.55_{-2.96}^{+63.91}$ | $38.99_{-8.07}^{+1.97}$ | $69.41_{-13.78}^{+2.55}$ | $1.39_{-0.24}^{+0.04}$ | $\gtrsim 1.92$ | - | -
Solid angle | 0 | $40.37_{-91.88}^{+16.76}$ | $48.23_{-11.60}^{+4.98}$ | $68.80_{-6.53}^{+4.05}$ | $1.61_{-0.27}^{+0.06}$ | $1.41_{-0.71}^{+0.27}$ | - | -
Isothermal | Area | 0 | $3.91_{-6.85}^{+8.71}$ | $11.91_{-3.36}^{+5.96}$ | $75.84_{-8.02}^{+7.70}$ | $1.21_{-0.07}^{+0.05}$ | $\gtrsim 2.84$ | - | -
Intermediate | 0.02 | $4.16_{-4.22}^{+10.39}$ | $14.58_{-6.12}^{+41.44}$ | $77.84_{-8.63}^{+8.34}$ | $0.71_{-0.25}^{+0.17}$ | - | - | -
Solid angle | 0 | $-0.52_{-40.44}^{+38.77}$ | $30.53_{-22.24}^{+21.60}$ | $72.22_{-15.41}^{+13.37}$ | $0.36_{-1.71}^{+1.65}$ | - | - | -
Radiation | Point | Area | 0.22 | $\mathbf{4.46_{-6.25}^{+6.19}}$ | $\mathbf{28.10_{-18.63}^{+7.04}}$ | $\mathbf{82.68_{-11.67}^{+5.51}}$ | $\mathbf{0.63_{-0.15}^{+0.16}}$ | $\mathbf{\gtrsim 1.50}$ | $\mathbf{76.17_{-23.73}^{+93.99}}$ | -
Intermediate | 0 | $2.84_{-39.36}^{+31.42}$ | $25.85_{-17.48}^{+23.03}$ | $72.75_{-16.06}^{+12.81}$ | $0.44_{-1.34}^{+1.62}$ | - | $64.30_{-38.31}^{+119.05}$ | -
Solid angle | 0 | $-2.52_{-36.65}^{+37.70}$ | $23.80_{-15.23}^{+19.36}$ | $67.20_{-12.52}^{+15.30}$ | $0.63_{-1.99}^{+1.57}$ | - | $81.87_{-53.75}^{+121.21}$ | -
Isothermal | Area | 0 | $7.27_{-6.56}^{+12.93}$ | $14.53_{-5.24}^{+15.59}$ | $76.91_{-9.60}^{+8.36}$ | $1.22_{-0.08}^{+0.07}$ | $3.08_{-1.23}^{+0.92}$ | $72.16_{-29.63}^{+107.31}$ | -
Intermediate | 0.01 | $3.96_{-3.94}^{+8.63}$ | $27.75_{-17.74}^{+6.47}$ | $80.07_{-8.25}^{+6.65}$ | $0.67_{-0.17}^{+0.14}$ | $1.41_{-0.99}^{+1.11}$ | $102.83_{-42.36}^{+83.02}$ | -
Solid angle | 0 | $16.17_{-58.68}^{+31.19}$ | $28.27_{-12.49}^{+14.78}$ | $63.98_{-4.99}^{+5.55}$ | $0.88_{-0.50}^{+0.62}$ | $2.03_{-0.53}^{+0.34}$ | $41.17_{-13.49}^{+14.20}$ | -
Hot gas | Point | Area | 0.23 | $\mathbf{3.10_{-5.37}^{+5.11}}$ | $\mathbf{54.51_{-24.74}^{+7.24}}$ | $\mathbf{78.17_{-6.32}^{+6.53}}$ | $\mathbf{0.57_{-0.13}^{+0.14}}$ | $\mathbf{\gtrsim 1.38}$ | - | $\mathbf{19.84_{-7.85}^{+7.47}}$
Intermediate | 0 | $-36.08_{-3.52}^{+78.08}$ | $38.63_{-4.44}^{+2.15}$ | $69.41_{-4.15}^{+2.93}$ | $1.33_{-0.19}^{+0.06}$ | $\gtrsim 1.62$ | - | $15.64_{-4.15}^{+7.47}$
Solid angle | 0 | $-35.48_{-11.73}^{+87.18}$ | $42.33_{-12.97}^{+7.20}$ | $69.33_{-6.51}^{+4.87}$ | $1.43_{-0.13}^{+0.09}$ | $\lesssim 0.96$ | - | $15.92_{-4.49}^{+8.27}$
Isothermal | Area | 0 | $2.96_{-6.73}^{+5.54}$ | $14.64_{-4.76}^{+15.43}$ | $70.28_{-8.67}^{+12.73}$ | $1.16_{-0.09}^{+0.10}$ | $\gtrsim 2.93$ | - | $13.18_{-3.80}^{+8.40}$
Intermediate | 0 | $2.80_{-3.96}^{+4.08}$ | $29.81_{-10.14}^{+4.41}$ | $79.38_{-9.57}^{+6.65}$ | $0.70_{-0.14}^{+0.13}$ | $\gtrsim 1.57$ | - | $24.63_{-11.85}^{+4.05}$
Solid angle | 0 | $-40.03_{-10.19}^{+88.70}$ | $42.71_{-9.32}^{+8.82}$ | $69.55_{-6.66}^{+6.15}$ | $1.50_{-0.18}^{+0.10}$ | $1.28_{-0.94}^{+0.55}$ | - | $14.53_{-4.01}^{+7.57}$
Southern hemisphere
Ideal | Point | Area | 0.85 | $6.09_{-36.60}^{+5.77}$ | $42.86_{-32.34}^{+7.57}$ | $65.26_{-6.51}^{+12.53}$ | $0.62_{-0.10}^{+0.08}$ | $1.65_{-0.46}^{+0.50}$ | - | -
Radiation | Point | Area | 0.15 | $7.74_{-42.15}^{+4.60}$ | $32.23_{-25.81}^{+11.20}$ | $68.33_{-7.67}^{+9.41}$ | $0.63_{-0.08}^{+0.08}$ | $1.71_{-0.39}^{+0.46}$ | $55.88_{-20.82}^{+18.29}$ | -
Large scale Northern hemisphere
Ideal | Point | Area | $0.65$ | $-0.25_{-14.60}^{+22.61}$ | $22.35_{-14.93}^{+36.49}$ | $84.40_{-9.37}^{+4.27}$ | $0.79_{-0.15}^{+0.11}$ | $1.97_{-1.08}^{+0.70}$ | - | -
Ideal | Isothermal | Area | 0.01 | $2.95_{-11.66}^{+11.98}$ | $17.08_{-9.38}^{+8.06}$ | $67.72_{-8.27}^{+8.03}$ | $1.63_{-0.13}^{+0.10}$ | $1.73_{-0.37}^{+0.47}$ | - | -
Ideal | Isothermal | Intermediate | 0.04 | $3.88_{-18.59}^{+22.68}$ | $18.81_{-12.01}^{+9.02}$ | $77.50_{-13.54}^{+8.64}$ | $0.99_{-0.19}^{+0.18}$ | $1.19_{-0.70}^{+1.02}$ | - | -
Radiation | Point | Area | 0.20 | $3.98_{-17.91}^{+27.68}$ | $21.33_{-13.66}^{+27.89}$ | $82.69_{-8.62}^{+5.53}$ | $0.79_{-0.14}^{+0.11}$ | $2.06_{-0.99}^{+0.66}$ | $78.76_{-39.33}^{+108.75}$ | -
Hot gas | Point | Area | 0.08 | $-5.67_{-6.41}^{+8.36}$ | $55.08_{-7.58}^{+6.73}$ | $77.79_{-6.83}^{+6.27}$ | $0.77_{-0.13}^{+0.12}$ | $1.70_{-0.55}^{+0.72}$ | - | $21.25_{-7.67}^{+6.48}$
Large scale Southern hemisphere
Ideal | Point | Area | $0.57$ | $-11.21_{-5.54}^{+34.96}$ | $40.95_{-24.80}^{+20.50}$ | $81.43_{-22.75}^{+6.79}$ | $0.78_{-0.17}^{+0.16}$ | $1.47_{-0.52}^{+0.72}$ | - | -
Ideal | Point | Intermediate | 0.29 | $-12.71_{-22.05}^{+13.89}$ | $39.85_{-10.25}^{+4.45}$ | $75.66_{-19.13}^{+11.86}$ | $0.90_{-0.14}^{+0.19}$ | $0.41_{-0.32}^{+0.86}$ | - | -
Radiation | Point | Area | 0.13 | $-10.66_{-7.39}^{+17.36}$ | $40.01_{-8.12}^{+2.81}$ | $80.41_{-20.75}^{+7.72}$ | $0.79_{-0.16}^{+0.17}$ | $1.68_{-0.58}^{+0.74}$ | $67.83_{-27.90}^{+93.86}$ | -

Table 2: Best-fitting parameters and Akaike weights obtained for various
combinations of driving mechanism (Column 1), potential (Column 2) and
expansion law (Column 3) for the warm neutral phase outflow. The top part of
the table shows results for the Northern hemisphere; here we show all possible
combinations of driving mechanism, potential, and expansion law, and highlight
the cases with the highest Akaike weights $w$ (Column 4). The bottom part of
the table is for the Southern hemisphere, and here we omit models with
$w<0.01$, which are ruled out with $>99\%$ confidence, for brevity. The fit
quantities shown in Columns 5 - 11 are the inclination angle $\phi$ (Column
5), inner and outer opening angle of the wind $\theta_{\rm in}$ and
$\theta_{\rm out}$ (Columns 6 and 7), logarithm of mass outflow rate
$\log\dot{M}$ measured in $M_{\odot}$ yr-1 (Column 8), logarithm of Mach
number $\mathcal{M}$ (Column 9), wind optical depth $\tau_{0}$ (Column 10;
radiatively-driven models only), and hot gas dimensionless velocity $u_{h}$
(Column 11; hot gas-driven models only). For most fit quantities, values are
specified with the 50th percentile value given first, and the 84th to 50th and
16th to 50th percentile ranges listed as the superscript and subscript values,
respectively; thus for example the entry $4.73^{+4.21}_{-5.10}$ for $\phi$ in
the (Ideal, Point, Area) case indicates a 50th percentile value
$\phi=4.73^{\circ}$, and a 16th to 84th percentile range $\phi=-0.37^{\circ}$
to $8.94^{\circ}$. The exception is $\mathcal{M}$, where in most cases our
posterior probability distribution is bounded above only by the prior we
impose; for these cases we report only a 16th percentile lower limit. Figure
3: Corner plot showing the 1D and 2D histograms of the posterior PDFs of all
parameters (see Table 1) describing the Northern side of the warm neutral
phase outflow traced by H i emission. The model fit shown is for an ideal wind
in a point gravitational potential, for clouds that maintain constant area as
they flow outward; this combination of driving, potential, and cloud expansion
yields the highest Akaike weight as compared to other combinations.

###  3.1 Warm neutral phase

As discussed in Section 2.1.1, for the H i data we define four distinct
regions: the Northern and Southern Hemispheres of the wind, and for both an
inner region that overlaps with the CO and H $\alpha$ data, and an outer
region where only H i is detected. We fit each region independently. For
simplicity we first discuss our results for the Northern inner region, using
this region as an example to explain and explore our method in Section 3.1.1
\- Section 3.1.4. We then present results for the the inner Southern region in
Section 3.1.5, and for the outer extraction region where only H i is detected
in Section 3.1.6. Since we are fitting the two hemispheres separately, we
divide the values of $\dot{M}$ returned by our MCMC fits by a factor of two,
since we are fitting only half the available solid angle at a time; all values
of $\dot{M}$ listed below include this factor of two division.

####  3.1.1 Fitting results

We summarise the best-fit parameters, uncertainties, and Akaike weights
returned by our pipeline as applied to the central region of the Northern side
of the outflow in Table 2. As discussed in Section 2.3.1, we fit each possible
combination of driving mechanism ($f_{p}$), gravitational potential ($m$), and
expansion law ($y$) separately, so we report a best fit for each combination.
However, from the Akaike weights shown in the table, we find that by far the
best fit to the data comes from a point gravitational potential and clouds
that maintain constant area; the Akaike weights corresponding to the three
possible driving mechanisms (ideal, radiative, and hot gas) are similar enough
that our fits do not allow us to distinguish between them. We show the joint
and marginal posterior PDFs for our fit parameters for the ideal, point
potential, constant area case in Figure 3; the results for the radiatively-
driven and hot gas-driven cases are qualitatively identical. We see that the
fit parameters are well constrained and largely non-degenerate, with most of
them converged into a tiny island of parameter space, with the exception of
$\mathcal{M}$, for which we only obtain a lower limit that
$\log\mathcal{M}\gtrsim 1.4$.

Our best-fitting geometric parameter values are largely consistent with and
extend the results of previous observational studies. The near-zero
orientation angle $\phi$ indicates the edge-on nature of the M82 and the large
opening angles $\theta_{\rm in}$ and $\theta_{\rm out}$ implies the geometry
of the outflow is approximately a conical sheath, as has been proposed before.
The mass outflow rate at the nuclear region is about $4$ M⊙ yr-1,
corresponding to a mass-loading factor of $\eta_{\rm
wn}=\dot{M}/\dot{M_{*}}\approx 1.0$, which is also within the range of values
derived in previous literature (Strickland & Heckman, 2009). However, our
uncertainty range is far small than in previous work – even taking the extreme
ranges of all of our acceptable models, our fit requires $\dot{M}$ to lie
between $3.0$ and $6.2$ M⊙ yr-1 with 68% confidence; we have therefore
constrained the mass outflow rate to $\approx 50\%$. We remind readers that
this is the mass outflow rate for the atomic phase on the Northern side of the
wind only; we return to other phases and to the Southern hemisphere below.

####  3.1.2 Model validation

Before accepting the best-fitting parameters we have obtained as definitive,
it is important to establish that the our parametric model gives a reasonable
description of that data. To this end, we compare our model-predicted spectra
to the observed ones in Figure 4. In this figure, the grey band shows the
observed spectrum (with $1\sigma$ errors) at each of our sample positions,
while the blue and orange lines show predicted spectra computed using
parameters from the single highest-likelihood model found by the MCMC (for the
ideal, point potential, constant area case), and from 10 randomly-selected
walkers at the final iteration of the MCMC, respectively. We see that the
theoretical spectra give a fairly good fit to the observations at all
positions, both in strength and shape.

Having shown that our model can provide a reasonable fit to the spectra at our
chosen positions, we next investigate whether our fits give a good match to
the entire 2D map, or, in other words, whether our model can reproduce the
data that are not included in the fits. In Figure 5 we compare the observed
integrated intensity and intensity-weighted second moment maps to model maps
generated using our highest-likelihood combination of parameters (blue line in
Figure 4). In this figure, we mask the region $z=0-0.6$ kpc to avoid
contamination from the disc. Outside the masked region, we see that the
predicted and observed maps show reasonably similar morphology and absolute
value, illustrating that our model successfully reproduces the bulk structure
of the outflow. We do see that our model predicts significantly more limb-
brightening than is observed in the integrated-intensity map, and a sharper
decrease in second moment as one moves away from the central axis. This may be
due to our oversimplification of the true geometry, as we assume there is a
hard cutoff of wind material at the inner and outer boundary. At a given
height, this feature of our model implies increasing line of sight distance
and decreasing velocity range as the impact parameter ($x$) approaches the
inner cone, which overestimates the velocity-integrated intensity and
underestimates the second moment. This artificial boundary effect is not
present in real galaxies, as the geometry is more likely to resemble a
biconical frustum (for which the “base” in the plane of the galaxy is a disc
of finite area, rather than a point), and the boundary of the volume occupied
by H i is much less sharp than in our idealised model.

Figure 4: Predicted versus observed H i 21 cm spectra. Each panel shows the
spectrum at one of our sample positions, expressed as brightness temperature
$T_{B}$ as a function of velocity $v$, with $v=0$ corresponding to the
systemic velocity of M82. The panels are arranged in the same pattern as in
Figure 1, i.e., the upper left panel here is the spectrum at the position
corresponding to the upper leftmost of the orange boxes indicated in Figure 1;
legends in each panel give the $(x,y)$ coordinates of the box centre, in units
of kpc. Blue lines show the predicted spectra for the set of parameters that
gives the largest posterior probability found by our MCMC fit, while orange
lines show spectra predicted using the parameters of 10 random walkers at the
last iteration of MCMC sampling. For comparison, we also show the observed
spectrum with its $1\sigma$ errors (grey band). Figure 5: Velocity-integrated
brightness temperature (top row) and brightness temperature-weighted second
moment map (bottom row) for the observations (left), theoretical predictions
using our best-fitting model (middle) and the residual between the two
normalised by the observations (right) for the Northern side H i 21 cm line
data. The black lines show our best-fitting angles $\theta_{\rm in}$ and
$\theta_{\rm out}$ for the inner and outer angles of the outflow cone. The
grey central region is a mask to block out $0-0.6$ kpc around the disc
midplane, where disc emission is significant. Figure 6: Same as Figure 4, but
now showing how the predicted spectra vary in response to changes in the fit
parameters. The blue line and grey region are the same as in Figure 4, and
show the highest-likelihood model and the observations, respectively. Other
colours show predicted spectra produced by starting from the highest-
likelihood parameters and increasing $\phi$ by $10^{\circ}$ (orange),
decreasing $\theta_{\rm in}$ and $\theta_{\rm out}$ by $10^{\circ}$ (green),
and doubling $\dot{M}$ (red).

####  3.1.3 Understanding parameter constraints

We next examine how model-predicted spectra vary as we change parameters, with
the goal of both justifying the tight constraints returned by our MCMC fits,
and gaining insight into which features of the observed spectra provide
constraints on particular parameters. In Figure 6, we show predicted spectra
produced by taking the single highest likelihood model (blue line in Figure 4)
and changing the parameters on which our model gives tight constraints one by
one. We first increase the inclination angle $\phi$ by 10 degrees, thereby
tilting the outflow axis relative to the plane of the sky (orange lines in
Figure 6). We see that the increased value of $\phi$ gives rise to an
asymmetry in the spectra. The reason for this is simple: positive velocities
correspond to the far side of outflow while negative velocities correspond to
the near side. A small value of $\phi$ means we are observing both near and
far side of wind at similar distances from the nuclear region, so only optical
depth effects, which are negligible in the case of H i, could induce asymmetry
between the negative and positive velocity sides of the spectra. However, for
non-zero $\phi$, our line of sight passes through the the near and far sides
of the outflow at two different distances from the nucleus. The near and far
sides therefore have different emissivities, leading to a sharp change in the
spectrum near zero velocity. Such a sharp feature is not observed, which is
why our model prefers an orientation angle $\phi$ close to zero.

The green line in Figure 6 shows the effect of decreasing both the inner and
outer opening angles $\theta_{\rm in}$ and $\theta_{\rm out}$ by 10∘. This
change makes the profile narrower for lines of sight close to the central axis
of the wind (central columns in Figure 6), and both narrower and more intense
near $v=0$ for lines of sight farther from the central axis (leftmost and
rightmost columns). The physical origin of this behaviour is that, as we
narrow the outflow cone at fixed $\dot{M}$, we both reduce the amount of
material for which the velocity vector is well-aligned to the line of sight
(thus depressing the high-velocity wings of the line) and create denser
outflows, increasing the intensity overall, particularly along heavily limb-
brightened lines of sight. Thus both the shape of the line wings and the
amount of variation in intensity from one line of sight to another serve to
constrain the opening angle.

Finally, we see that, as expected, the mass outflow rate acts as primarily as
a normalisation of the spectra, affecting the overall intensity while leaving
the shape of the profile mostly unchanged.444The shapes of the spectra are not
completely insensitive to $\dot{M}$ because higher $\dot{M}$ requires higher
$\Gamma$, i.e., that the momentum flux be closer to the Eddington value; this
in turn manifests as somewhat different wind acceleration laws. However, this
is a second-order effect. With increased $\dot{M}$, the generated spectra
overestimate the observed ones at nearly all velocities and positions. This
strong effect explains why our model is able to produce such strong
constraints on $\dot{M}$.

Figure 7: Same as Figure 4, but now showing how the predicted spectra vary in
response to the choices of potential and expansion law, for ideal winds. We
show the best-fitting spectra for each individual combination of potential and
expansion law. The spectrum produced with the combination that yields the
highest Akaike weight is shown by the thickest line and the observed spectrum
with 1$\sigma$ errors is shown by the grey region, as in Figure 4. Spectra
adopting point and isothermal gravitational potentials are shown as solid and
dashed lines, respectively, while those adopting constant area, intermediate,
and constant solid angle expansion laws are shown in blue, yellow and red
line, respectively.

####  3.1.4 Understanding constraints on the potential and expansion law

We have seen that our fitting not only provides strong constraints on wind
parameters for fixed combinations of driving mechanism ($f_{p}$),
gravitational potential ($m$), and expansion law ($y$), comparison between the
models clearly favours particular potentials and expansion laws. As a
reminder, the former describes the shape of the gravitational potential, and
thus the rate of gravitational deceleration, as the wind moves out. The latter
dictates how clouds’ cross-sectional areas, and thus the amount of momentum
per unit time they intercept from the wind driving mechanism, change as they
flow outward.

To understand why we are able to constrain these parameters, it is helpful to
begin by examining Table 2. In this table, we see that some models yield
unrealistically large orientation angles $\phi\sim\pm 40^{\circ}$ with a huge
uncertainty, while others yield $\phi\approx 0$, consistent with the edge-on
nature of M82. The reason for this is simple: we observe emission near zero
velocity in all spectra, which implies either that there is some gas that
truly has near-zero velocity at all positions, or that the zero-velocity
emission is coming from gas whose velocity vector lies entirely in the plane
of the sky. The first option is only possible for certain combinations of
cloud expansion and potential: examining the acceleration law, Equation 2, we
see that acceleration of the wind can remain close to zero at all radii only
for the combinations $y=1$, $m=1$ (constant area, point potential), $y=1$,
$m=a$ (constant area, isothermal potential), and $y=a$, $m=a$ (intermediate
area, isothermal potential) – and it is precisely these combinations that
produce $\phi\approx 0$. For all other combinations, the presence of emission
at zero velocity forces the MCMC fit to select models that have substantial
amounts of material launched along the plane of the sky, which, for a conical
sheath geometry, in turn requires $\phi\gg 0$ or $\phi\ll 0$. However, this
compensation worsens the fit to the remainder of the spectrum, yielding a
lower Akaike weight overall.

We illustrate this effect in Figure 7, which shows how the quality of fits
varies with the different choices of potential and expansion law, for ideal
winds. There are six lines, corresponding to each possible combination of
potential (indicated by solid versus dashed) and expansion law (indicated by
colour); each line corresponds to the set of parameters with the highest
likelihood recovered by the MCMC for that combination of expansion law and
potential. We see that although all models produce near zero-velocity
emission, the models where the wind is more strongly accelerated ($y=a$ or
$y=a^{2}$) clearly show a deficit of emission at the centre, and a line shape
function completely different to that of the observed spectra. In contrast,
the combination $y=1$, $m=1$ that has the highest Akaike weight provides a
much better comparison to the observations, both at line centre and line
wings. Our method is therefore able to distinguish between different
combinations of cloud expansion and potential with confidence.

####  3.1.5 Results for the Southern hemisphere

Thus far our discussion has been limited to the Northern side of the wind. We
use the same method to analyse the Southern hemisphere, again limiting
ourselves for now to the inner region that overlaps the CO data. We summarise
the best-fit parameters, uncertainties, and Akaike weights for the Southern
side of the outflow in Table 2. For brevity here we show results only for
models with Akaike weight $w>0.01$, since all other models are ruled out with
$>99\%$ confidence. As for Northern hemisphere, we find that the best fit to
the data comes from a point gravitational potential and clouds that maintain
constant area, and the Akaike weights corresponding to the first two possible
driving mechanisms are still similar enough that our fits do not allow us to
distinguish between them. However, unlike in the Northern hemisphere, the fit
parameters for the South converge to more than one island, though within each
island they are well constrained and largely non-degenerate. Overall we find
that our model provides a somewhat poorer, less constrained fit to the
Southern Hemisphere than to the North, likely because the wind in the South
suffers greater tidal forces from M81 that produce structures not captured in
our model; these tidal effects are evident in the data, which show significant
asymmetry about the minor axis of M82, and a significant blue-shift compared
to the Northern Hemisphere. However, for the largest island, we still recover
near-zero orientation angles, and most importantly, mass outflow rates with
values similar to those of Northern hemisphere, illustrating the robustness of
out fit. We further note that there are changes to some parameters compared to
the Northern Hemisphere, especially for the opening angles, which change by
$\Delta\theta_{\rm in}~{}10^{\circ}$ and $\Delta\theta_{\rm
out}~{}-20^{\circ}$, which shows the asymmetric geometry of the wind system.

####  3.1.6 Results for the large-scale extraction region

We have thus far presented our fits for the central region of M82 where data
for all three phases overlap. However, warm neutral phase traced by H i
emission extends well beyond the region over which the other phases are
detected. Therefore it is worthwhile to investigate this single phase outflow
to the greatest extent we can achieve; in addition to providing further
information by itself, this allows us to perform a useful consistency check,
since our inner and outer extraction regions represent non-overlapping parts
of the wind. This exercise therefore tests whether our results are
reproducible when we look at a different part of the same structure. We
therefore repeat our analysis for the Northern and Southern outer extraction
regions.

Figure 8: Same as figure Figure 3, but show the results for H i at larger
scale

As in Section 3.1.1, we summarise the fitting results in Table 2; the table
reports fits for the large-scale region on both the North and South sides, and
we again omit models with Akaike weights $<0.01$ for brevity. From the models
remaining, we again find that clouds of constant area being driven out of a
point potential provide the best fit to the data (though models with
intermediate expansion are not ruled out quite as strongly as for the inner
region), and we are unable to differentiate strongly between different
possible wind acceleration mechanisms. We show the joint and marginal
posterior PDFs for our fit parameters for the ideal, point, area case for the
Northern hemisphere in Figure 8; results for the Southern hemisphere are
qualitatively similar. We see that the fit parameters converge into two
islands of parameter space characterised by a near-zero orientation angle
$\phi$ and reasonable values of opening angles $\theta_{\rm in}$ and
$\theta_{\rm out}$, verifying that the warm neutral phase of the outflow also
displays a cone sheath geometry at large scale. The mass outflow rate is
approximately $6.2$ M⊙ yr-1, corresponding to a mass-loading factor of
$\approx 1.50$. This mass outflow rate is consistent at the $1\sigma$ level
with our previous estimate based on the inner part of the wind, which
demonstrates the robustness of our model to different scales we are fitting.
The constraints on the geometry are also consistent at the $1\sigma$ level,
but we find that our fits to the larger-scale data have much larger
uncertainties, $\sim 20^{\circ}$ for $\theta_{\rm in}$ and $\sim 10^{\circ}$
for $\phi$ and $\theta_{\rm out}$. This may be caused by the stronger tidal
torque induced by M81 on the larger scale.

Figure 9: Same as Figure 4, but showing the results for H i at larger scale.
As in Figure 4, each panel corresponds to one of the boxes shown in Figure 1.

In Figure 9 we compare predicted and observed spectra for the best-fitting
case for the Northern hemisphere, just as we did in Figure 4. We again see
that the theoretical spectra give a fairly good fit to the observations at all
positions.

The relationship between the Northern and Southern hemispheres of the wind on
the large scale are very similar to those on small scale. The Southern
hemisphere gives comparable mass outflow rate and slightly different opening
angle compared to the North.

###  3.2 Cold molecular phase

We next repeat the main steps in the analysis presented in Section 3.1 for the
cold molecular phase as traced by the CO $J=2\to 1$ observations.

CO $J=2\to 1$ (Leroy et al., 2015)
---
Models |  | Best fit parameters
Driver | Potential | Expansion | $w$ | $\phi$ | $\theta_{\mathrm{in}}$ | $\theta_{\mathrm{out}}$ | $\log\dot{M}$ | $\log\mathcal{M}$ | $\tau_{0}$ | $u_{h}$
Northern hemisphere |
Ideal | Point | Area | 0.54 | $-3.16_{-12.22}^{+9.57}$ | $11.34_{-7.85}^{+11.36}$ | $80.87_{-10.92}^{+6.30}$ | $0.30_{-0.17}^{+0.19}$ | $0.83_{-0.30}^{+0.22}$ | - | -
Radiation | Point | Area | 0.24 | $-2.50_{-12.06}^{+11.88}$ | $11.88_{-8.33}^{+11.80}$ | $77.52_{-9.92}^{+8.40}$ | $0.29_{-0.17}^{+0.22}$ | $0.84_{-0.34}^{+0.29}$ | $78.08_{-29.00}^{+99.18}$ | -
Hot Gas | Point | Area | 0.22 | $-4.59_{-9.61}^{+8.73}$ | $14.63_{-9.98}^{+9.93}$ | $80.28_{-9.17}^{+6.84}$ | $0.36_{-0.20}^{+0.15}$ | $0.74_{-0.25}^{+0.28}$ | - | $14.66_{-5.34}^{+10.09}$
Southern hemisphere
Ideal | Point | Area | 0.01 | $11.57_{-22.24}^{+18.33}$ | $26.63_{-6.76}^{+7.99}$ | $54.60_{-3.26}^{+12.44}$ | $-0.00_{-0.11}^{+0.83}$ | $1.49_{-1.40}^{+0.21}$ | - | -
Radiation | Point | Area | 0.04 | $17.97_{-28.48}^{+16.09}$ | $15.38_{-7.12}^{+11.09}$ | $55.05_{-4.06}^{+11.19}$ | $-0.13_{-0.25}^{+0.74}$ | $1.71_{-1.15}^{+0.33}$ | $61.89_{-21.75}^{+35.19}$ | -
Hot gas | Point | Area | 0.95 | $18.96_{-39.51}^{+3.39}$ | $10.72_{-6.92}^{+28.25}$ | $56.22_{-3.29}^{+19.68}$ | $-0.34_{-0.10}^{+0.98}$ | $2.07_{-1.43}^{+0.26}$ | - | $8.81_{-4.61}^{+4.11}$

Table 3: Same as table Table 2, but show the results of fits for CO. Figure
10: Same as figure Figure 3, but show the results for CO.

Similar to Section 3.1, we summarise the fitting results to the Northern side
of the outflow in Table 3. We omit models with Akaike weights $<0.01$ for
brevity, and from the models remaining, we again find that clouds of constant
area being driven out of a point potential provide the best fit to the data,
and that we are unable to differentiate strongly between different possible
wind acceleration mechanisms. We show the joint and marginal posterior PDFs
for our fit parameters for the ideal, point, area case in Figure 10. Again, we
see that the fit parameters converge into a tiny island of parameter space
characterised by a near-zero orientation angle $\phi$ and reasonable values of
opening angles $\theta_{\rm in}$ and $\theta_{\rm out}$, verifying that the
molecular phase of the outflow also displays a cone sheath geometry (Walter et
al., 2017). The mass outflow rate is approximately $2.2$ M⊙ yr-1,
corresponding to a mass-loading factor of $\eta_{\rm
cm}=\dot{M}/\dot{M_{*}}\approx 0.55$. As with H i, the mass outflow rate is
very well-constrained; considering the extremes of all three models with
acceptable Akaike weights, it lies between $1.3$ and $3.2$ M⊙ yr-1 with 68%
confidence.

Figure 11: Same as Figure 4, but showing the results for CO. As in Figure 4,
each panel corresponds to one of the boxes shown in Figure 1.

In Figure 11 we compare predicted and observed spectra for the best fitting
case, just as we did in Figure 4 for the H i data. We again see that the
theoretical spectra give a fairly good fit to the observations at all
positions, though the quality of the fits for H i is somewhat better than that
for CO. This is due to the greater complexity of modeling CO lines, which
require us to make several approximations that may fail in certain parts of
the wind. First, we rely on the large velocity gradient (LVG) approximation to
evaluate optical depths, which becomes invalid near zero velocity, where
thermal broadening becomes important. Second, we assume that the molecules are
in LTE at a fixed temperature at all positions; both assumptions likely become
invalid in low-density regions, where the collision rate is too low to
thermalise the molecules, or to efficiently radiate away the energy deposited
by shocks. The failure of this assumption may explain why our models fit some
positions better than others.

Finally, and most importantly, we caution that our need to choose a fixed
temperature for the CO molecules likely induces a systematic uncertainty that
exceeds the uncertainties quoted above. To first order, for an optically thick
line such as CO $J=2\rightarrow 1$, we expect the emissivity per unit mass to
scale close to linearly with the gas temperature, so a factor of two
uncertainty on the gas temperature corresponds roughly to a factor of two
uncertainty in the mass flux. Since the confidence intervals we obtain from
our fits are considerably narrower than this, the gas temperature likely
represents the single largest source of uncertainty remaining. Unfortunately,
it is not possible to constrain the temperature using the single CO line to
which we have access.

The relationship between the Northern and Southern hemispheres of the wind for
CO are very similar to those for H i. Specifically, while the model spectra
qualitatively reproduce the observed ones reasonably well, there are clear
tidal features produced by M81 that are not captured, and as a result the fit
is somewhat more uncertain. The mass outflow rates we infer are consistent
within the uncertainties with those found in the North, but the uncertainties
in the South as substantially larger, so that the 68% confidence interval
spans $\approx 1$ dex.

###  3.3 Warm ionised phase

We finally repeat the main steps in the analysis presented in Section 3.1 for
the warm ionised phase as traced by the H$\alpha$ observations, with the
differences in analysis described in Section 2.3.3. We consider three possible
values for the clumping factor, $c_{\rho}=10,100,1000$, spanning a range of
values measured in individual H ii regions (Kennicutt, 1984). We summarize the
fitting results to the Northern side of the outflow in Table 4, but only for
three highest-weight combinations of driving mechanism, potential and
expansion, for each of the three clumping factors. Contrary to what we find
for H i and CO 2$\to$ 1 line, the fitting for H$\alpha$ favors point
gravitational potential and clouds that maintain constant solid angle (or
intermediate cross sectional area), and we still cannot differentiate strongly
between the driving mechanisms. However, the fit parameters still converge
into a tiny island of parameter space, with a near-zero orientation angle
$\phi$ and reasonable values of opening angle $\theta_{\rm in}$ and
$\theta_{\rm out}$, verifying the cone sheath geometry of warm ionised outflow
(Martin, 1998). Most of the fit parameters are similar under the variation of
clumping factors, with the exception of the mass outflow. We see the
uncertainty of mass flux is about 10 across three clumping factors, from $5$
M${}_{\odot}{\rm yr}^{-1}$ at $c_{\rho}=10$, to $1.62$ M${}_{\odot}{\rm
yr}^{-1}$ at $c_{\rho}=100$, to $0.34$ M${}_{\odot}{\rm yr}^{-1}$ at
$c_{\rho}=1000$.

Figure 12: Predicted versus observed H$\alpha$ spectrum for the Northern
hemisphere. Each panel shows one choice of clumping factor. Blue lines show
the predicted spectra for the set of parameters that gives the largest
posterior probability found by our MCMC fit, while orange lines show spectra
predicted using the parameters of 10 random walkers at the last iteration of
MCMC sampling. For comparison, we also show the observed spectrum with its 1
$\sigma$ errors (grey region).

In Figure 12, we compare predicted and observed spectra, for the ideal, point
potential, constant solid angle case, and for three choices of clumping
factor. We see that the theoretical spectra give a fairly good fit to the
observations, particularly at line wings. All fits qualitatively reproduce the
dip near line centre that is the hallmark of cone sheath geometry and an
expansion law whereby clouds continually accelerate – this is feature in the
H$\alpha$ data, which is not present in the H i or CO, is the reason why our
model prefers constant solid angle or intermediate expansion for H$\alpha$,
and constant area for H i and CO. However, we also see that there is no clear
systematic difference between the spectra predicted for different clumping
factors, which limits our model’s capacity to distinguish differences in
clumping factor. This in turn introduces an order of magnitude uncertainty
into the measurement of mass outflow rate. We conclude that the accurate
measurement of mass outflow rate for H$\alpha$ data, and by extension for any
tracer for which the emissivity scales quadratically rather than linearly with
the local volume density, is possible only if one has access to independent
constraints on this factor.

We also apply the pipeline to the Southern hemisphere. However, we find that
the MCMC cannot distinguish a favoured driving mechanism, potential, or
expansion law with confidence, and the model-predicted spectra are a poor
match to the observations. This is likely because the observed spectrum is
severely blue-shifted due to tidal force of M81. For these reasons we do not
discuss our fitting results further for the Southern hemisphere.

H$\alpha$ (Martin, 1998)
---
Models |  | Best fit parameters
Driver | Potential | Expansion | $w$ | $\phi$ | $\theta_{\mathrm{in}}$ | $\theta_{\mathrm{out}}$ | $\log\dot{M}$ | $A$ | $\log\mathcal{M}$ | $\tau_{0}$ | $u_{h}$
North
$c_{\rho}=10$
Ideal | Point | Intermediate | 0.30 | $-0.68_{-2.15}^{+1.52}$ | $40.88_{-5.59}^{+2.98}$ | $77.50_{-8.45}^{+8.41}$ | $0.16_{-0.21}^{+0.24}$ | $0.37_{-0.06}^{+0.08}$ | $1.80_{-0.52}^{+0.70}$ | - | -
Ideal | Point | Solid angle | 0.49 | $-10.95_{-14.09}^{+10.95}$ | $29.66_{-10.57}^{+15.14}$ | $68.65_{-12.77}^{+13.83}$ | $0.73_{-0.24}^{+0.17}$ | $0.46_{-0.07}^{+0.06}$ | $1.98_{-0.59}^{+0.59}$ | - | -
Radiation | Point | Solid angle | 0.21 | $-1.59_{-10.53}^{+1.84}$ | $28.31_{-3.21}^{+3.69}$ | $79.09_{-11.57}^{+7.90}$ | $0.67_{-0.43}^{+0.23}$ | $0.38_{-0.06}^{+0.09}$ | $1.41_{-0.87}^{+1.00}$ | $74.07_{-13.41}^{+35.31}$ | -
$c_{\rho}=100$
Ideal | Point | Intermediate | 0.32 | $-0.61_{-3.42}^{+2.77}$ | $40.30_{-12.68}^{+3.69}$ | $77.08_{-10.28}^{+9.44}$ | $-0.37_{-0.30}^{+0.22}$ | $0.42_{-0.08}^{+0.08}$ | $2.51_{-0.73}^{+0.35}$ | - | -
Ideal | Point | Solid angle | 0.53 | $-5.94_{-16.22}^{+5.91}$ | $23.64_{-4.23}^{+18.25}$ | $75.54_{-10.84}^{+9.69}$ | $0.21_{-0.34}^{+0.21}$ | $0.44_{-0.09}^{+0.06}$ | $1.87_{-0.98}^{+0.80}$ | - | -
Radiation | Isothermal | Solid angle | 0.15 | $-0.93_{-8.70}^{+2.11}$ | $28.90_{-5.00}^{+4.00}$ | $80.46_{-10.28}^{+7.13}$ | $-0.47_{-0.33}^{+0.39}$ | $0.37_{-0.08}^{+0.06}$ | $0.50_{-0.28}^{+0.56}$ | $76.57_{-13.72}^{+62.42}$ | -
$c_{\rho}=1000$
Ideal | Point | Intermediate | 0.37 | $-0.93_{-5.02}^{+2.32}$ | $41.32_{-7.62}^{+3.32}$ | $76.16_{-8.74}^{+9.15}$ | $-1.04_{-0.34}^{+0.30}$ | $0.40_{-0.07}^{+0.10}$ | $2.34_{-0.65}^{+0.46}$ | - | -
Ideal | Point | Solid | 0.53 | $-2.55_{-18.45}^{+3.29}$ | $22.26_{-3.58}^{+22.69}$ | $76.73_{-9.47}^{+7.77}$ | $-0.47_{-0.36}^{+0.28}$ | $0.42_{-0.07}^{+0.08}$ | $1.74_{-0.56}^{+0.64}$ | - | -
Radiation | Isothermal | Solid | 0.10 | $-1.06_{-13.82}^{+2.96}$ | $31.48_{-6.28}^{+6.21}$ | $78.77_{-12.24}^{+7.80}$ | $-1.46_{-0.42}^{+0.70}$ | $0.38_{-0.08}^{+0.10}$ | $0.38_{-0.23}^{+0.66}$ | $74.30_{-19.24}^{+76.51}$ | -

Table 4: Same as Table 2, but showing the results of fits to the H$\alpha$
data. The three blocks in the table show results for three different clumping
factors, $c_{\rho}=10$, 100, and 1000, as indicated.

##  4 Discussion

Having presented our fits for the multi-phase outflow of M82, we are now in a
position to come up with an overall physical picture of the cool wind in M82.
We first discuss the mass budget of different phases of the wind in Section
4.1, and present a broader overview of our model for the wind structure in
Section 4.2. We discuss the implications of our results for the origin and the
evolution of launched cool gas clouds in Section 4.3. We then compare our
self-consistent methods of identifying outflow rates and distinguishing winds
from fountains to earlier methods in Section 4.4, and discuss the implications
for the outflow energy and momentum budget in Section 4.5 We finally discuss
several caveats to our model and their possible effects on our results in
Section 4.6.

###  4.1 The M82 wind mass budget

A primary result of our analysis is that we have, for the first time, obtained
tight constraints on the mass outflow rates for the cool phases of the M82
wind. Collating the results from Table 2 for the inner extraction region,
where we have data in all tracers, and taking a simple average of the 50th,
16th, and 84th percentile values for cases with acceptable Akaike weights
(reasonable, since the results are qualitatively identical), we measure mass
outflow rates of $\dot{M}=4.0^{+1.6}_{-1.0}$ M⊙ yr-1 and $4.2^{+0.85}_{-0.79}$
M⊙ yr-1 in the Northern and Southern sides of the H i wind, respectively; if
we approximate the distributions as Gaussian, and sum the outflow rates in
both hemispheres, the 16th to 84th percentile range for the total outflow rate
is $\dot{M}=8.2^{+1.8}_{-1.3}$ M⊙ yr-1. The corresponding mass-loading factors
in the warm neutral phase, for our fiducial star formation rate
$\dot{M}_{*}=4.1$ M⊙ yr-1, are $\eta_{\rm wn}=1.0^{+0.40}_{-0.25}$ and
$1.0^{+0.21}_{-0.20}$ for the Northern and Southern winds separately, and
$\eta_{\rm wn}=2.0^{+0.44}_{-0.31}$ for the sum.

Repeating this exercise for the CO data (Table 3), the molecular outflow rates
in the Northern and Southern sides separately are $\dot{M}=2.0^{+1.1}_{-0.65}$
M⊙ yr-1 and $0.70^{+4.2}_{-0.13}$ M⊙ yr-1, respectively, and the total
molecular outflow rate is $\dot{M}=2.7^{+4.3}_{-0.66}$ M⊙ yr-1. The
corresponding mass loading factors in the cold molecular phase are $\eta_{\rm
cm}=0.48^{+0.27}_{-0.16}$ and $0.17^{+1.0}_{-0.031}$ for the Northern and
Southern sides separately, and $\eta_{\rm cm}=0.66^{+1.1}_{-0.16}$ for the
total. We remind readers that, in addition to our quoted statistical
uncertainties, we also have a factor of $\sim 2$ systematic uncertainty
arising from the unknown gas temperature in the wind.

We can therefore draw a few high-level conclusions. The total mass outflow
rate for the neutral (atomic plus molecular) phases of M82 is $\approx 10$ M⊙
yr-1, and the total mass loading factor for these phases $\approx 2-3$. We can
be confident in these numbers to better than a factor of two uncertainty. The
atomic phase carries moderately more mass than the molecular phase, but both
contribute at order unity. We are much less certain about the contribution
from the warm ionised phase traced by H$\alpha$, due to the uncertain clumping
factor. It could be as much as comparable to the neutral phases, if the
clumping factor is small, or it could be an order of magnitude less, if the
clumping factor is large. For comparison, Martin (1998) estimate a mass
outflow rate of $24$ M⊙ yr-1 for an assumed H$\alpha$ volume filling factor
(equivalent to $1/c_{\rho}$, where $c_{\rho}$ is our clumping factor) of
$0.1$. This is a factor of several larger than our estimate of $\approx 5$ M⊙
yr-1 at $c_{\rho}=10$, but given the large differences in assumed geometry, a
difference of this size is not particularly surprising.

###  4.2 Structure of the wind

We can also characterise the structure of the wind deduced from our fits more
broadly. Since our goal here is qualitative understanding rather than detailed
statistical analysis, for clarity we will show only a single result for each
CO and H i component, omitting H$\alpha$ since our fits for it are generally
poor, and using only the results from the inner region, where the CO and H i
data overlap. The models we show are those derived from the single highest-
likelihood model found by our MCMC for that component, using ideal driving,
constant area clouds, and a point-like potential, which is in all cases is one
of the models with high Akaike weight. We pause to note that our finding that
a point-like potential generally fits the data better than an isothermal
potential implies, per our discussion in Section 2.2.1, that either the gas is
more reflective of the potential shape than the stars, or that the potential
is flattened enough that most of the region that we survey, which begins $0.8$
kpc above the plane, is above the height where the isothermal-Keplerian
transition occurs.

Figure 13: Radial profiles of, from top to bottom, mean density, mean radial
velocity, momentum flux, and energy flux carried by the different wind
components. In all panels, blue shows H i and orange shows CO, solid lines
show the Northern hemisphere and dashed lines show the Southern hemisphere,
and the grey vertical line shows the wind launch radius $r_{0}=250$ pc. For
density, the black dotted line shows as $\rho\propto r^{-2}$ scaling for
comparison; the left axis is mass density, and the right axis is number
density of H nuclei, computed assuming a mean mass of $2.34\times 10^{-24}$ g
per H nucleon, as expected for Solar abundances. In the velocity panel, dark
lines show the mass flux-weighted mean and light lines show the mass-weighted
mean; see Equation 12. For the momentum and energy fluxes, the solid black
line is the sum of all components; the left axis shows total flux, and the
right axis shows flux normalised by star formation rate, so that the quantity
shown is the rate per unit mass of stars formed at which stellar feedback must
supply momentum and energy to drive the wind to that radius.

We show the radial profiles of mean density, mean velocity, radial momentum
flux, and radial energy flux in Figure 13. We derive these mean values
following equations 13 and 17 of KTOM. Specifically, for each component we use
the parameters of the best-fitting model derive the mass flux $\dot{M}$, the
fraction of area covered by the wind $f_{A}$, the dispersion $\sigma_{x}$ of
the lognormal surface density distribution $p_{M}(x)$ (equation 10 of KTOM),
the critical logarithmic surface density $x_{\rm crit}$ below which material
is ejected, and the mass fraction $\zeta_{M}$ below this surface density. From
these, plus the wind acceleration law $u_{a}(x)$ for the highest-likelihood
model (Equation 2), we define the mean density and the momentum and energy
fluxes at any specified radius $a=r/r_{0}$ by

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\frac{\dot{M}}{4\pi
f_{A}r_{0}^{2}v_{0}}\left[\frac{1}{\zeta_{M}}\int_{-\infty}^{x_{\rm
crit}}\frac{p_{M}}{a^{2}u_{a}}\,dx\right]$ (9) $\displaystyle\dot{p}$
$\displaystyle=$
$\displaystyle\dot{M}v_{0}\left[\frac{1}{\zeta_{M}}\int_{-\infty}^{x_{\rm
crit}}u_{a}p_{M}\,dx\right]$ (10) $\displaystyle\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dot{M}v_{0}^{2}\left[\frac{1}{\zeta_{M}}\int_{-\infty}^{x_{\rm
crit}}u^{2}_{a}p_{M}\,dx\right].$ (11)

Since the differential mass of material launched with surface density $x$
scales as $p_{M}(x)$, and the differential mass flux as $p_{M}(x)u_{a}(x)$,
these definitions also immediately permit us to define the mass-weighted and
mass flux-weighted mean velocities as

$\langle v\rangle_{M}=\frac{\dot{p}}{\dot{M}}\qquad\langle
v\rangle_{\dot{M}}=\sqrt{\frac{2\dot{E}}{\dot{M}}}.$ (12)

There are a few features of Figure 13 worthy of comment. First, the density
and velocity structure asymptote to the values expected for a constant-
velocity wind ($v\sim\mbox{constant}$ and $\rho\propto r^{-2}$) only at radii
of a few kpc. At smaller radii, the wind is still accelerating, and as a
result the density profile is steeper than $r^{-2}$. Second, both the mass-
weighted and mass flux-weighted mean velocities are rather small, $\lesssim
150$ km s-1, and at no radius does the mean velocity exceed the escape speed
from the launch region, $v_{0}=170$ km s-1. While there is a tail of material
extending to higher velocity (as is required to explain the observed line
wings), most of the mass moves much more slowly. The reason it is possible for
material moving below the escape speed to nevertheless escape the galaxy is
that the gas is not ballistic; instead, the driving mechanism supplies
momentum continuously rather than instantaneously, so the net force on gas
entrained into the wind is outward rather than inward. Third, this slow
acceleration of the wind is reflected in where the wind acquires it energy.
The wind acquires $50\%$ of its final, terminal kinetic energy $\gtrsim 0.5$
kpc away from the galaxy, and $25\%$ $\gtrsim 1$ kpc away.

Figure 14:  $X_{\mathrm{CO}(2\to 1)}$ or $\alpha_{\mathrm{CO}(2\to 1)}$ as a
function of position with the M82 wind; the coordinate system here is the same
as that used in Figure 1. Dashed lines indicate the inner and outer wind cone
edges projected on the sky, and the grey around the origin indicates the wind
launch region, where our model is not valid. The hatching marks the region
$|z|<0.8$ kpc that we mask in our analysis due to contamination from disc
emission; our model is fit only to data outside this region.

We can also use this model to explore the extent to which the CO “$X$” factor
(or equivalently the $\alpha$ factor) varies across the wind, which is
important for interpreting CO observations. Recall that in our model we do not
assume a light to mass conversion for CO; instead, we compute it self-
consistently from the large-velocity gradient approximation (equation 65 of
KTOM). We show the results in Figure 14. We see that there is a factor of
$\sim 2$ difference between the North and South sides of the wind; this is
likely a reflection not of a real difference, so much as a reflection of the
range of uncertainty in our model; there are other parameter options found by
the MCMC with fits that are nearly as good, but where the North-South
differences are smaller. For the North side, we generally have
$\alpha_{\mathrm{CO}(2\to 1)}\sim 1$ M⊙ pc-2 / [K km s-1], comparable to or
slightly smaller than the estimate provided by Leroy et al. (2015); or values
on the South side are a factor of $\sim 2$ smaller, but, as noted above, this
level of disagreement is not surprising given both Leroy et al.’s stated
factor of $\approx 2$ systematic uncertainty and the comparably large
uncertainty in our analysis that arises from our need to assume a gas
temperature.

A more interesting disagreement is the trend of $\alpha_{\mathrm{CO}(2\to 1)}$
with position: Leroy et al. find higher values of $\alpha_{\mathrm{CO}(2\to
1)}$ in the bright parts of the wind than in the faint parts, while we find
the opposite. It is unclear which estimate is correct. Leroy et al. derive
their conclusions from the dust infrared emission, and their analysis relies
on the assumption that the dust in the wind has the same basic emission
properties (implying similar temperatures and sizes) as the dust in the disc;
this hypothesis may or may not be correct. In our model, where the gas
temperature and level populations are assumed to be constant, the trends in
$\alpha_{\mathrm{CO}(2\to 1)}$ arise simply from the variation of the line-of-
sight velocity gradient across the wind: gas closer to the launch region and
to the central outflow axis has larger LOS velocity gradients, and hence
smaller $\alpha_{\mathrm{CO}(2\to 1)}$, because those are regions where the
wind radial velocity is changing most rapidly, and where the radial direction
is most closely aligned with our line of sight. Conversely, regions far from
the central axis and the launch point have slowly-varying wind radial
velocity, and radial vectors that are mostly in the plane of the sky. These
are basic physical effects that it seems must be present; however, it is
possible that in reality they are mitigated or outweighed by variations in the
gas temperature or excitation state with position, which we do not model.

###  4.3 Origin and evolution of cool gas clouds

As discussed in Section 1, it is unclear if the cool gas observed in galactic
winds is primarily entrained as cold gas, or if it represents hot gas that has
cooled and condensed at some distance from the galactic plane. In M82 we see
cool gas spread continuously at locations from immediately above the galactic
plane to several kpc off it, which suggest that we are observing in this
system is entrainment rather than recondensation. However, this raises the
further question of how these cold clouds are able to survive entrainment into
the hot wind without disruption.

Our results provide insight into possible solutions, since they suggest that,
when launched, cold molecular and warm neutral clouds maintain approximately
constant cross-sectional area as they propagate outwards, so their
morphologies are well-preserved through the entrainment processes. Previous
theoretical studies have not reached consensus on whether, and under what
circumstances, this is possible. Purely hydrodynamic simulations from Hopkins
& Elvis (2010) show that Kelvin-Helmholtz (KH) instabilities stretch cold
clouds and increase their cross-sectional area dramatically, while MHD studies
from (McCourt et al., 2015; Banda-Barragán et al., 2018b; Zhang et al., 2017;
Grønnow et al., 2017) show that the magnetic pressure with $\beta\sim 1$ can
suppress mixing (Gentry et al., 2019), confining cloud expansion transverse to
the direction of magnetic field, and allowing the cloud to survive being
entrained. Alternately, expansion can be limited if clouds are accelerated by
infrared radiation pressure (Huang et al., 2020), though the modest column
density of M82 suggests that infrared radiation pressure is likely unimportant
for it (Crocker et al., 2018). Our results therefore favour a magnetic
confinement scenario, and hence suggest that M82 is characterised by a strong
magnetic field with field lines tracing the trajectory of launched clouds.
This scenario is consistent with the recent measurements of magnetic field
strength and morphology in M82 (Yoast-Hull et al., 2013; Jones et al., 2019;
Buckman et al., 2020; Lopez-Rodriguez et al., 2021), which show that the
magnetic field threads from starburst core all the way to intergalactic medium
with a strength of $B\sim 300$ $\mu$G and $\beta\sim 0.7$.

###  4.4 Outflows versus winds

While our fits allow us to quantify with relatively small uncertainty how much
gas the winds expel, we are left with another important question: what is the
ultimate fate of the gas in the galactic winds in M82? Will it be driven out
to the CGM or even further into IGM, or will it be captured by the
gravitational potential and fall back as a fountain? Our model fits suggest
the former. In the KTOM model, whether a given outflow is a wind or a fountain
depends on the indices of the cloud expansion law, $y\propto a^{p}$, and the
gravitational potential, $m\propto a^{q}$. If $p\geq q$, then the outward
force on clouds grows faster than or keeps pace with gravity, so that clouds
always move outward and never fall back, while for $p<q$ gravity eventually
overwhelms acceleration, resulting in a fountain. Almost all our best fitting
models are $p=0$ (constant area) and $q=0$ (point-like potential); the only
exceptions are the large-scale H i in the South, where a model with $p=1$ and
$q=0$ (and thus still a wind) is acceptable (Akaike weight $w=0.29$), and
H$\alpha$ in the North, where all our acceptable fits again have $p\geq q$ and
thus are winds. Therefore our fits generally favour a wind over a fountain,
though only by a small amount since we have $p=q$ in most cases.

This conclusion is at odds with those of Leroy et al. (2015) and Martini et
al. (2018), though we are fitting the same data. It is therefore important to
understand why. There are two main features in the data that Leroy et al.
(2015) and Martini et al. (2018) identify as favouring a fountain. The first
is that the surface brightness of the CO and H i emission fall steeply within
$\approx 1-2$ kpc of the galactic centre, which Leroy et al. interpret as a
steep fall in the mass flux within this region (c.f. their Figure 18). Our
fits also recover this steep fall in surface brightness, but are able do so
with a constant mass flux, i.e., a wind rather than a fountain. The underlying
reason for this difference in interpretation is a difference in assumptions:
KTOM model cool clouds accelerating continuously under the action of some
source of momentum deposition, as illustrated in Figure 13, while Leroy et al.
assume that cool material is immediately accelerated to its terminal speed
near the galactic plane, and after that can only decelerate in response to
gravity or drag. The latter assumption imposes strong constraints on the
surface brightness profile, because if the velocity can only decrease with
radius, then a fixed mass flux implies that density must decrease with radius
no more sharply than $\rho\propto r^{-2}$, or, equivalently, that the surface
density must decrease with projected radius no more steeply than $r_{\rm
proj}^{-1}$. Thus Leroy et al.’s assumption of instantaneous acceleration
means that the observed fall in surface brightness with radius, which is
steeper than $r_{\rm proj}^{-1}$ near the galactic centre, can only be
interpreted as a diminution of the mass flux. If one relaxes this assumption,
and allows the acceleration to take place over a finite range of radius, then
a fit such as that shown in Figure 13 becomes possible, whereby the mass flux
remains constant and yet the surface density fall more steeply than $r_{\rm
proj}^{-1}$ because the wind is still accelerating at small radii.

It is worth emphasising two further points in this regard. First, the KTOM
model does not require a slowly-accelerating wind; the range in radius over
which the wind accelerates is controlled by the parameters of the model.
Therefore the fact that the MCMC fit prefers gradual acceleration is
significant; gradual acceleration gives a better match to the observed
spectra. Second, there is no physical reason to assume that the cool
components of winds are instantaneously accelerated. While such an assumption
has been common in the literature for reasons of simplicity, there is no
physical model for how the cool components of winds are accelerated that
produces instantaneous acceleration. Quite the opposite: if the acceleration
were instantaneous, the resulting shock would almost certainly dissociate and
ionise the cool neutral gas. More gradual acceleration, as in the model we
present here, not only fits the data just as well, it is a much more natural
physical model.

The other line of evidence that Martini et al. (2018) argue favours a fountain
rather than an outflow is the position-velocity diagram of the H i in the
Southern hemisphere along the minor axis (c.f. their Figure 3), which shows a
relatively sharp decrease in maximum H i velocity as one moves away from the
galactic centre. They interpret this decrease as evidence that the wind is
decelerating. To understand why our fits do not take this decrease as decisive
evidence for deceleration, it is helpful to examine the underlying spectra,
which we show for the Southern hemisphere along the line $x=-1.5$ kpc555We
show $x=-1.5$ kpc rather than $x=0$, i.e., exactly along the minor axis,
because this avoids a region of non-detections along the minor axis. By
contrast, at $x=-1.5$ we achieve strong detections of the signal at all
vertical distances from $-3$ to $-8$ kpc; c.f. Figure 1. in Figure 15.
Examining the figure, it is clear that there is a decrease in emission in the
negative velocity wing from $y=-3$ to $-4$ kpc, but decreases at larger radii
are much less clear. In particular, it seems possible that the lack of
emission beyond $v=-150$ km s-1 are larger radii is simply an issue of the
signal in the line wing diminishing below the detection limit, rather than an
actual decrease in the gas velocity that would be indicative of deceleration.
Much the same point can be made with regard to Martini et al.’s Figure 3, and
the general effect – that weak lines systematically underestimate the velocity
range of emission because the line wings drop below the detection threshold
faster than the line core – is a well-known bias in atomic and molecular line
observations (e.g., Yuan et al., 2020). Since our model likelihood function
properly accounts for the uncertainties on the data (c.f. Equation 5), it will
naturally account for this effect, which may be why our fits indicate that
winds are acceptable models even for the large-scale H i emission in the
Southern hemisphere.

Figure 15: H i spectra in the Southern hemisphere, at positions from
$(x,y)=(-1.5,-3)$ kpc to $(-1.5,-8)$ kpc, as indicated in the legend. Filled
bands indicate the uncertainty range.

###  4.5 Outflow energy and momentum budget

Our best-fitting model also eases some tension in the energy and momentum
budgets for models in which the wind is assumed to flow at constant speed. We
use the model proposed by Leroy et al. (2015) as an example, but emphasise
that the issues we point out are generic to constant-speed wind models. In
order for feedback to launch a wind with mass flux $\dot{M}$ at fixed speed
$v_{w}$, it must provide a momentum and kinetic energy per unit mass of stars
formed

$\displaystyle\frac{p}{M_{*}}$ $\displaystyle=$
$\displaystyle\frac{\dot{M}}{\dot{M}_{*}}v_{w}$ (13)
$\displaystyle\frac{E}{M_{*}}$ $\displaystyle=$
$\displaystyle\frac{\dot{M}}{2\dot{M}_{*}}v_{w}^{2},$ (14)

where $\dot{M}_{*}$ is the star formation rate of the stellar population
powering the wind. Leroy et al.’s model has $v_{w}=450$ km s-1, and under
their assumption of constant $v_{w}$, the high surface brightness at small
radii requires that the outflow mass flux at launch be $\approx 50$ M⊙ yr-1;
this only diminishes to $\approx 10$ M⊙ yr-1 at $2-3$ kpc above the plane.
Inserting these factors into the equations above, we find $p/M_{*}=5400$ km
s-1 and $E/M_{*}=2.5\times 10^{49}$ erg M${}_{\odot}^{-1}$.666These figures
are for our fiducial star formation rate $\dot{M}_{*}=4.1$ M⊙ yr-1; adopting
Leroy et al. (2015)’s preferred star formation rate of $\approx 10$ M⊙ yr-1
would halve them, but this would not change the qualitative conclusions we
draw. These figures are difficult to supply: for blast waves driven by single
SNe, Gentry et al. (2017) find that the terminal momentum and kinetic energy
are $\approx 3\times 10^{5}$ M⊙ km s-1 and $\approx 3\times 10^{49}$ erg,
while for the most efficient superbubbles produced by clustered SNe (those
that suffer the smallest radiative losses), these budgets rise to $\approx
3\times 10^{6}$ M⊙ km s-1 and $\approx 10^{50}$ erg per SN. For a standard IMF
that produces approximately one supernova per 100 M⊙ of stars formed, we
therefore have firm upper limits of $p/M_{*}\approx 3\times 10^{4}$ km s-1 and
$E/M_{*}\approx 10^{48}$ erg M${}_{\odot}^{-1}$ for superbubbles, and $3000$
km s-1 and $3\times 10^{47}$ erg M${}_{\odot}^{-1}$ for non-clustered SNe. The
constant velocity model for the wind fits within the momentum budget for
superbubbles, but it exceeds the energy budget by a factor of $\approx 30$;
for single SNe, it exceeds the momentum budget by a factor of 2, and the
energy budget by a factor of 100.777Our conclusions about the energy budget
are different from those of Leroy et al. (2015), who argue that their model
fits within $\approx 30\%$ of the available energy. Partly this is because
Leroy et al. adopt a mass flux of $\approx 10$ M⊙ yr-1, which their model
reaches only $\approx 3$ kpc above the plane, rather than the $\approx 50$ M⊙
yr-1 with which it starts. However, a larger issue is that Leroy et al. assume
that all of the $\approx 10^{51}$ erg provided by a SN can be converted to
kinetic energy of the outflow, and none is lost to radiation. Neither analytic
models nor simulations support such an assumption. Gentry et al. (2017) find
that even in the most energy-efficient superbubbles put $\approx 90\%$ of the
initial SN energy into radiation, while for single SNe this rises to $\approx
97\%$. Thus the largest plausible mechanical energy budget is a factor of 10
smaller than the value that Leroy et al. adopt.

The energy and momentum requirements of a KTOM wind are smaller. This in part
because relaxing the assumption of constant speed allows the mass flux at the
base to be smaller, and in part because relaxing the assumption that all wind
material travels at a single speed allows the high-velocity line wings to be
produced by a relatively small mass of low-density material that is
accelerated to high speed, while the bulk of the mass moves more slowly.
Quantitatively, we can compute the wind energy and momentum requirements
simply by evaluating Equation 10 and Equation 11 in the limit $a\to\infty$, in
which case $u_{x}\to\Gamma e^{x}-1$. Doing so using the highest-likelihood
values for ideal, point, constant area models, as in Section 4.2, we find that
our wind model requires a momentum and kinetic energy per unit mass of stars –
including both hemispheres and both the H i and CO components – of
$p/M_{*}=210$ km s-1 and $E/M_{*}=2.7\times 10^{47}$ erg M${}_{\odot}^{-1}$,
respectively. (Roughly 75% of this is in the H i component, and 25% in the
CO.) Thus our model, unlike the constant velocity model, fits comfortably
within the available energy and momentum budgets.

A broader point to take from this analysis is that models in which a galactic
wind is assumed to be ballistic or to move at constant speed – both constant
with radius, and constant with density – can be seriously misleading. These
assumptions can make the difference between interpreting observations as a
wind versus a fountain, and they can lead to order of magnitude errors in the
inferred energy budget.

###  4.6 Caveats

We end this section by noting several caveats for our results, as our
analytical model inevitably ignores some physical processes.

The first is that our model assumes that each phase in the wind is composed of
a continuous population of clouds with a temperature and chemical composition
that is independent of time and distance from the central galaxy. However, it
is possible that galactic winds evolve in thermal or chemical properties as
the they propagate outward. Depending on the nature of the change, this could
either move gas into or out of the cool phases as the outflows propagate
outward. The observed weak cooling in the wind of M82 found by Hoopes et al.
(2003) suggests that cooling of hot phases into the cool phase is unlikely to
be a strong effect in this galaxy, but it is unclear if this result applies to
outflows in general, or only to M82. Conversely, clouds may dissociate from
molecular gas into atomic gas, causing the density of molecular gas to fall
and that of atomic gas to rise, or they may ionise, causing mass to disappear
in H i and appear in tracers such as H$\alpha$. Determining how these thermal
and chemical processes affect the cool wind requires future observations of
the chemical and thermal structure of the wind of M82 and a more accurate
description of the thermal and chemical properties of the wind in our model.

The second limitation of our modelling originates from our ignorance of the
larger environment around M82, which may induce complex structures in the real
outflows that are not accounted for in our simple model. For example, we have
already seen that there is evidence that the torque from M81 induces a tidal
stream from M82, and this might alter the structure and distribution of matter
in the outflow of M82. This effect is most obvious for the warm neutral
outflow traced by H i, which shows a ‘northwest spur’ and ‘southeast spur’
induced by the drag exerted by ambient environment (Martini et al., 2018),
This ‘NW-SE’ asymmetry is in contrast with our model prediction, which is
symmetric about the minor axis.

##  5 Conclusion

###  5.1 The wind of M82

Analyses of galactic wind observations have long been limited lack of
theoretical tools to extract information from the rich, 3D data sets that we
can obtain from line observations of nearby galaxies. In this paper, we make a
first step to improve this situation by using a novel semi-analytic wind model
to constrain the properties of the multi-phase outflow of the nearby starburst
M82. We model the outflow as a continuous population of gas clouds being
momentum-driven out of a turbulent galactic disc. The kinematic structure of
the wind is then determined by model parameters describing the outflow
geometry and mass flux, and physical prescriptions for the outflow driving
mechanism, the gravitational potential, and the rates at which cool clouds
expand as they flow outward. We use this model to fit position-position-
velocity (PPV) data measured in the H i 21 cm, CO $J=2\to 1$, and H$\alpha$
lines for the wind of M82. Our main conclusions are as follows:

  * •

Our best-fitting model shows good overall agreement with the observations,
both for the sample spectra used for the fits and for the full 2D moment maps,
although several areas of tension still remain due to the simplified geometry
and physics recipes adopted by our model.

  * •

The wind has an edge-on orientation and a conical sheath geometry for all
phases, consistent with earlier visual analyses. The atomic outflow has a very
broad outer opening angle $(\gtrsim 80^{\circ})$, reaching nearly to the
equator, while the CO is still broad, but somewhat more narrowly confined
(outer opening angle $\approx 60^{\circ}$).

  * •

The total mass flux carried by the warm neutral phase is $\approx 8$ M⊙ yr-1,
while that carried by the molecular phase is $\approx 2$ M⊙ yr-1. The latter
is uncertain at the factor of $\approx 2$ level, while the former is tightly
constrained, with only $\approx 25\%$ uncertainty. These fluxes are similar to
each other and to the SFR at the order-of-magnitude level, implying a total
mass-loading factor of a few. The mass flux is the ionised component is far
less certain, due to the quadratic dependence of H$\alpha$ emissivity on the
local volume density. We tentatively conclude that the ionised gas mass flux
is no more than that in the neutral phases, but are unable to draw any
stronger conclusions.

  * •

We find that the atomic and molecular clouds retain near constant area as the
flow outward, rather than expanding to intercept more momentum. This suggests
a picture where strong magnetic fields in M82 thread the launched gas clouds
and prevent them from expanding and subsequently being destroyed by Kelvin-
Helmholtz instability as they move outward along their trajectories.

  * •

At least over the region covered by the H i and CO data cubes, the outflow is
most consistent with being a wind that will escape the galaxy, rather than a
fountain that will fall back. Attempts to deduce whether the outflow is a wind
or a fountain that make the simplifying assumption that the outflow is
instantaneously accelerated to its final velocity, rather than being allowed
to accelerate continuously, likely produce misleading conclusions.

###  5.2 Future prospects

A natural question to follow this work is to what extent similar detailed
model-fitting can be used in other systems, and for which systems it is likely
to be the most fruitful. One obvious point to draw from our analysis is that
the cleanest, most unambiguous results – and, in M82, the best fits – come
from H i observations. These have the advantage over molecular data that there
is no need to adopt a gas temperature or to assume that the molecules are in
local thermodynamic equilibrium, and has the advantage over collisionally-
excited lines such as H$\alpha$ (or C ii) that there is no need to worry about
a clumping factor. Unfortunately, due to the weakness of the 21 cm line,
measurements of H i emission from winds are scarce. The largest extant sample
of extrplanar H i is the HALOGAS sample of 15 galaxies (Marasco et al., 2019).
However, MeerKAT is expected to increase this number dramatically over the
next few years (Maddox et al., 2021). These observations are very promising
targets for the methods we have deployed here, and application of the KTOM
model to them should similarly allow constraints on the neutral gas wind mass
flux to better than a factor of two.

Another possibly powerful approach is to use molecules, but reduce the
uncertainty about the temperature and excitation state by observing multiple
transitions and/or multiple isotopologues. Within galactic discs, measuring
spatially-resolved measurements in a sufficiently large set of such lines
allows strong constraints on the position-dependent gas temperature and degree
of thermal equilibration (e.g., Sharda et al., 2021); in principle the same
methods could be used to derive these quantities for galactic winds, thereby
removing our single largest source of systematic uncertainty.

A third approach that we have not explored here, but that can be used directly
within the context of the KTOM model, is absorption line measurements. Such
measurements also provide a method to constrain outflows that is not dependent
on assumed gas temperatures or clumping factors, and surveys have begun to
gather samples of reasonable size (e.g., Schroetter et al., 2016, 2019). Here
the primary systematic uncertainty is likely to be the abundances of the
absorbing species, coupled to the fact that such measurements will establish
the outflow rate only in the warm ionised phase, not in other phases.
Unfortunately the set of galaxies for which large-scale absorption
measurements exist and the set for which atomic and molecular measurements
exists is essentially non-overlapping. Thus it is likely to be some time
before we can obtain strong constraints on the outflow rates in both the
neutral and warm ionised phases for the same galaxy.

## Acknowledgements

We thank A. Leroy and P. Martini for providing us with access to their H i and
CO data. MRK acknowledges support from the Australian Research Council through
Future Fellowship award FT180100375. This research was undertaken with the
assistance of resources and services from the National Computational
Infrastructure (NCI), grant jh2, which is supported by the Australian
Government. We also gratefully acknowledge the contribution of Mike
Fitzpatrick who wrote a python script to byte swap the echelle IRAF images,
making them readable on modern computers.

## Data availability

The source code and data used in this paper are available from
https://github.com/yyx319/multiphase-outflow-M82.

## References

  * Alatalo et al. (2011) Alatalo K., et al., 2011, ApJ, 735, 88
  * Banda-Barragán et al. (2016) Banda-Barragán W. E., Parkin E. R., Federrath C., Crocker R. M., Bicknell G. V., 2016, MNRAS, 455, 1309
  * Banda-Barragán et al. (2018a) Banda-Barragán W. E., Federrath C., Crocker R. M., Bicknell G. V., 2018a, MNRAS, 473, 3454
  * Banda-Barragán et al. (2018b) Banda-Barragán W. E., Federrath C., Crocker R. M., Bicknell G. V., 2018b, MNRAS, 473, 3454
  * Banda-Barragán et al. (2021) Banda-Barragán W. E., Brüggen M., Heesen V., Scannapieco E., Cottle J., Federrath C., Wagner A. Y., 2021, MNRAS,
  * Bregman (1980) Bregman J. N., 1980, ApJ, 236, 577
  * Bregman et al. (1995) Bregman J. N., Schulman E., Tomisaka K., 1995, ApJ, 439, 155
  * Brüggen & Scannapieco (2020) Brüggen M., Scannapieco E., 2020, ApJ, 905, 19
  * Buckman et al. (2020) Buckman B. J., Linden T., Thompson T. A., 2020, MNRAS, 494, 2679
  * Carr et al. (2018) Carr C., Scarlata C., Panagia N., Henry A., 2018, ApJ, 860, 143
  * Chabrier (2005) Chabrier G., 2005, The Initial Mass Function: From Salpeter 1955 to 2005. p. 41, doi:10.1007/978-1-4020-3407-7_5
  * Chevalier & Clegg (1985) Chevalier R. A., Clegg A. W., 1985, Nature, 317, 44
  * Coker et al. (2013) Coker C. T., Thompson T. A., Martini P., 2013, ApJ, 778, 79
  * Cottle et al. (2020) Cottle J., Scannapieco E., Brüggen M., Banda-Barragán W., Federrath C., 2020, ApJ, 892, 59
  * Crocker et al. (2018) Crocker R. M., Krumholz M. R., Thompson T. A., Clutterbuck J., 2018, MNRAS, 478, 81
  * Crocker et al. (2021a) Crocker R. M., Krumholz M. R., Thompson T. A., 2021a, MNRAS, p. arXiv:2006.15821
  * Crocker et al. (2021b) Crocker R. M., Krumholz M. R., Thompson T. A., 2021b, MNRAS, 502, 1312
  * Field (1959) Field G. B., 1959, ApJ, 129, 536
  * Fluetsch et al. (2019) Fluetsch A., et al., 2019, MNRAS, 483, 4586
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Freedman et al. (1994) Freedman W. L., et al., 1994, ApJ, 427, 628
  * Gentry et al. (2017) Gentry E. S., Krumholz M. R., Dekel A., Madau P., 2017, MNRAS, 465, 2471
  * Gentry et al. (2019) Gentry E. S., Krumholz M. R., Madau P., Lupi A., 2019, MNRAS, 483, 3647
  * Greco et al. (2012) Greco J. P., Martini P., Thompson T. A., 2012, ApJ, 757, 24
  * Gronke & Oh (2018) Gronke M., Oh S. P., 2018, MNRAS, 480, L111
  * Gronke & Oh (2020) Gronke M., Oh S. P., 2020, MNRAS, 494, L27
  * Grønnow et al. (2017) Grønnow A., Tepper-García T., Bland-Hawthorn J., McClure-Griffiths N. M., 2017, ApJ, 845, 69
  * Heckman & Thompson (2017) Heckman T. M., Thompson T. A., 2017, arXiv e-prints, p. arXiv:1701.09062
  * Heckman et al. (1990) Heckman T. M., Armus L., Miley G. K., 1990, ApJS, 74, 833
  * Hoopes et al. (2003) Hoopes C. G., Heckman T. M., Strickland D. K., Howk J. C., 2003, ApJ, 596, L175
  * Hopkins & Elvis (2010) Hopkins P. F., Elvis M., 2010, MNRAS, 401, 7
  * Hopkins et al. (2011) Hopkins P. F., Quataert E., Murray N., 2011, MNRAS, 417, 950
  * Huang et al. (2020) Huang X., Davis S. W., Zhang D., 2020, ApJ, 893, 50
  * Jones et al. (2019) Jones T. J., et al., 2019, ApJ, 870, L9
  * Kennicutt (1984) Kennicutt R. C., 1984, ApJ, 287, 116
  * Kennicutt (1998) Kennicutt R. C., 1998, ApJ, 498, 541
  * Kennicutt et al. (2008) Kennicutt R. C., Lee J. C., Funes J. G., J. S., Sakai S., Akiyama S., 2008, ApJS, 178, 247
  * Kim et al. (2020a) Kim C.-G., et al., 2020a, arXiv e-prints, p. arXiv:2006.16315
  * Kim et al. (2020b) Kim C.-G., et al., 2020b, ApJ, 903, L34
  * Krumholz (2014) Krumholz M. R., 2014, MNRAS, 437, 1662
  * Krumholz & Thompson (2013) Krumholz M. R., Thompson T. A., 2013, MNRAS, 434, 2329
  * Krumholz et al. (2017) Krumholz M. R., Thompson T. A., Ostriker E. C., Martin C. L., 2017, MNRAS, 471, 4061
  * Krumholz et al. (2019) Krumholz M. R., McKee C. F., Bland -Hawthorn J., 2019, ARA&A, 57, 227
  * Leroy et al. (2015) Leroy A. K., et al., 2015, ApJ, 814, 83
  * Li & Tonnesen (2020) Li M., Tonnesen S., 2020, ApJ, 898, 148
  * Li et al. (2020) Li Z., Hopkins P. F., Squire J., Hummels C., 2020, MNRAS, 492, 1841
  * Lochhaas et al. (2018) Lochhaas C., Thompson T. A., Quataert E., Weinberg D. H., 2018, MNRAS, 481, 1873
  * Lopez-Rodriguez et al. (2021) Lopez-Rodriguez E., Guerra J., Asgari-Targhi M., Schmelz J. T., 2021, arXiv e-prints, p. arXiv:2102.03362
  * Lopez et al. (2020) Lopez L. A., Mathur S., Nguyen D. D., Thompson T. A., Olivier G. M., 2020, ApJ, 904, 152
  * Lynds & Sandage (1963) Lynds C. R., Sandage A. R., 1963, ApJ, 137, 1005
  * Maddox et al. (2021) Maddox N., et al., 2021, A&A, 646, A35
  * Mao & Ostriker (2018) Mao S. A., Ostriker E. C., 2018, ApJ, 854, 89
  * Marasco & Fraternali (2017) Marasco A., Fraternali F., 2017, MNRAS, 464, L100
  * Marasco et al. (2019) Marasco A., et al., 2019, A&A, 631, A50
  * Martin (1998) Martin C. L., 1998, ApJ, 506, 222
  * Martini et al. (2018) Martini P., Leroy A. K., Mangum J. G., Bolatto A., Keating K. M., Sandstrom K., Walter F., 2018, ApJ, 856, 61
  * McCourt et al. (2015) McCourt M., O’Leary R. M., Madigan A.-M., Quataert E., 2015, MNRAS, 449, 2
  * McKeith et al. (1995) McKeith C. D., Greve A., Downes D., Prada F., 1995, A&A, 293, 703
  * Melioli et al. (2015) Melioli C., Brighenti F., D’Ercole A., 2015, MNRAS, 446, 299
  * Muratov et al. (2015) Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691
  * Muratov et al. (2017) Muratov A. L., et al., 2017, MNRAS, 468, 4170
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 477, 450
  * O’Connell & Mangano (1978) O’Connell R. W., Mangano J. J., 1978, ApJ, 221, 62
  * O’Connell et al. (1995) O’Connell R. W., Gallagher John S. I., Hunter D. A., Colley W. N., 1995, ApJ, 446, L1
  * Oppenheimer & Davé (2006) Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265
  * Prochaska et al. (2011) Prochaska J. X., Kasen D., Rubin K., 2011, ApJ, 734, 24
  * Roberts-Borsani (2020) Roberts-Borsani G. W., 2020, MNRAS, 494, 4266
  * Roberts-Borsani & Saintonge (2019) Roberts-Borsani G. W., Saintonge A., 2019, MNRAS, 482, 4111
  * Roberts-Borsani et al. (2020) Roberts-Borsani G. W., Saintonge A., Masters K. L., Stark D. V., 2020, MNRAS, 493, 3081
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Scannapieco & Brüggen (2015) Scannapieco E., Brüggen M., 2015, ApJ, 805, 158
  * Scarlata & Panagia (2015) Scarlata C., Panagia N., 2015, ApJ, 801, 43
  * Schneider et al. (2018) Schneider E. E., Robertson B. E., Thompson T. A., 2018, ApJ, 862, 56
  * Schneider et al. (2020) Schneider E. E., Ostriker E. C., Robertson B. E., Thompson T. A., 2020, ApJ, 895, 43
  * Schroetter et al. (2016) Schroetter I., et al., 2016, ApJ, 833, 39
  * Schroetter et al. (2019) Schroetter I., et al., 2019, MNRAS, 490, 4368
  * Sharda et al. (2021) Sharda P., et al., 2021, MNRAS,
  * Sharma (2017) Sharma S., 2017, ARA&A, 55, 213
  * Shopbell & Bland-Hawthorn (1998) Shopbell P. L., Bland-Hawthorn J., 1998, ApJ, 493, 129
  * Sofue et al. (1992) Sofue Y., Reuter H. P., Krause M., Wielebinski R., Nakai N., 1992, ApJ, 395, 126
  * Steidel et al. (2010) Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289
  * Strickland & Heckman (2009) Strickland D. K., Heckman T. M., 2009, ApJ, 697, 2030
  * Strickland et al. (1997) Strickland D. K., Ponman T. J., Stevens I. R., 1997, A&A, 320, 378
  * Thompson & Krumholz (2016) Thompson T. A., Krumholz M. R., 2016, MNRAS, 455, 334
  * Thompson et al. (2015) Thompson T. A., Fabian A. C., Quataert E., Murray N., 2015, MNRAS, 449, 147
  * Thompson et al. (2016) Thompson T. A., Quataert E., Zhang D., Weinberg D. H., 2016, MNRAS, 455, 1830
  * Tumlinson et al. (2011) Tumlinson J., et al., 2011, Science, 334, 948
  * Veilleux et al. (2005) Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 43, 769
  * Veilleux et al. (2020) Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, arXiv e-prints, p. arXiv:2002.07765
  * Walter et al. (2002) Walter F., Weiss A., Martin C., Scoville N., 2002, AJ, 123, 225
  * Walter et al. (2017) Walter F., et al., 2017, ApJ, 835, 265
  * Watson et al. (1984) Watson M. G., Stanger V., Griffiths R. E., 1984, ApJ, 286, 144
  * Werk et al. (2019) Werk J. K., et al., 2019, ApJ, 887, 89
  * Westmoquette et al. (2009) Westmoquette M. S., Gallagher J. S., Smith L. J., Trancho G., Bastian N., Konstantopoulos I. S., 2009, ApJ, 706, 1571
  * Wouthuysen (1952) Wouthuysen S. A., 1952, AJ, 57, 31
  * Wright (1750) Wright T., 1750, An original theory or new hypothesis of the universe : founded upon general phaenomena of the visible creation; and particularly the Via the laws of nature, and solving by mathematical principles : the Lactea …compris’d in nine familiar letters from the author to his friendand : illustrated with upward of thirty graven and mezzotinto plates …, doi:10.3931/e-rara-28672.
  * Yoast-Hull et al. (2013) Yoast-Hull T. M., Everett J. E., Gallagher J. S. I., Zweibel E. G., 2013, ApJ, 768, 53
  * Yuan et al. (2020) Yuan Y., Krumholz M. R., Burkhart B., 2020, MNRAS, 498, 2440
  * Zhang et al. (2017) Zhang D., Thompson T. A., Quataert E., Murray N., 2017, MNRAS, 468, 4801