# The Impact of Cosmic Ray Injection on Magnetic Flux Tubes in a Galactic Disk

Roark Habegger University of Wisconsin-Madison Astronomy Department
University of Wisconsin-Madison Physics Department  Ellen G. Zweibel
University of Wisconsin-Madison Astronomy Department  University of Wisconsin-
Madison Physics Department  Sherry Wong University of Wisconsin-Madison
Astronomy Department  University of Wisconsin-Madison Computer Sciences
Department

(Received 2022-Nov-07; Revised 2023-Apr-03; Accepted 2023-Apr-21)

###### Abstract

In a seminal paper, Parker (1966) showed the vertical stratification of the
interstellar medium (ISM) is unstable if magnetic fields and cosmic rays
provide too large a fraction of pressure support. Cosmic ray acceleration is
linked to star formation, so Parker’s Instability and its nonlinear outcomes
are a type of star formation feedback. Numerical simulations have shown the
instability can significantly restructure the ISM, thinning the thermal gas
layer and thickening the magnetic field and cosmic ray layer. However, the
timescale on which this occurs is rather long $(\sim 0.4\,\mathrm{Gyr})$.
Furthermore, the conditions for instability depend on the model adopted for
cosmic ray transport. In this work, we connect the instability and feedback
problems by examining the effect of a single, spatially and temporally
localized cosmic ray injection on the ISM over $(\sim 1\,\mathrm{kpc}^{3})$
scales. We perform cosmic ray magnetohydrodynamic simulations using the
Athena++ code, varying the background properties, dominant cosmic ray
transport mechanism, and injection characteristics between our simulation
runs. We find robust effects of buoyancy for all transport models, with
disruption of the ISM on timescales as short as $100$ Myr when the background
equilibrium is dominated by cosmic ray pressure.

Galaxy Structure (622) - Cosmic Rays (329) - Magnetohydrodynamical simulations
(1966)

††journal: The Astrophysical Journal††software: Athena++ (Stone et al., 2020;
Jiang & Oh, 2018), MatPlotLib (Hunter, 2007), NumPy (van der Walt et al.,
2011; Harris et al., 2020), AstroPy (Astropy Collaboration et al., 2013, 2018)

##  1 Motivation

In 1966, Eugene Parker showed the interstellar medium’s (ISM’s) stratification
was unstable when non-thermal pressure support by magnetic field and cosmic
rays was introduced (Parker, 1966). The free energy source for the instability
is the gravitational potential energy of gas held above its natural scale
height. Parker argued that undular perturbations of the field lines, which on
average are tangent to the Galactic plane, would allow gas to slide into the
magnetic troughs, while the crests were free to buoyantly rise. The end result
would be the production of cosmic ray inflated magnetic lobes rising
perpendicular to the galactic disk and pockets of dense gas in the Galactic
plane.

In order for the instability to occur, the equation of state (EoS) of the
medium must be soft enough to allow compression into pockets. Parker assumed
the cosmic ray pressure would remain constant along field lines and that the
gas would follow an isothermal EoS. Under these conditions, any pressure
support from magnetic fields and/or cosmic rays renders the system unstable.

The Milky Way’s thermal energy density, magnetic energy density, and cosmic
ray energy density are all on the order of $\sim
10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-3}$; e.g. (Ferrière, 2001; Draine,
2011). In general, the far infrared - radio correlation (Bell, 2003) suggests
that star formation, cosmic ray acceleration, and magnetic field generation
are linked. This correlation means the Parker instability connects star
formation to the disruption of the ISM and the production of outflows. So, the
Parker instability is a form of star formation feedback and could play a key
role in galactic evolution. It could mediate the formation of dense clouds,
but also heat and disrupt them (Bustard & Zweibel, 2021). It could launch mass
from the disk as a wind or fountain (Everett et al., 2008; Uhlig et al., 2012;
Hanasz et al., 2013; Girichidis et al., 2016; Ruszkowski et al., 2017; Wiener
et al., 2017; Farber et al., 2018; Zhang, 2018; Hopkins et al., 2018; Chan et
al., 2019; Hopkins et al., 2021). Perhaps less directly related to feedback,
it could contribute to cosmic ray escape from the disk, and could influence
the galactic dynamo by growing the component of the field perpendicular to the
disk (Hanasz & Lesch, 2000; Hanasz et al., 2004, 2009).

While the non-thermal pressure in cosmic rays and aforementioned buoyancy
mechanism provide a way for cosmic rays to affect galactic evolution and the
structure of the ISM, there is evidence against the Parker instability being a
dominant effect in galaxies. The original Parker instability model has been
revised to include more physical processes and more realistic setups
(Mouschovias, 1974; Zweibel & Kulsrud, 1975; Asseo et al., 1978; Giz & Shu,
1993; Heintz & Zweibel, 2018), and analyzed in the nonlinear regime with
numerical simulations (Heintz et al., 2020; Tharakkal et al., 2022). Overall,
these results lead to more stringent conditions for instability and show the
instability requires a longer time to vertically restructure the ISM;
typically over $400\,\mathrm{Myr}$ (Heintz et al., 2020; Tharakkal et al.,
2022). This long timescale means the effects of the instability could be
overshadowed by other processes with shorter timescales.

One key difference between the different analyses lies in the treatment of
cosmic ray transport and the choice of cosmic ray EoS. A stiffer adiabatic
index ($\gamma_{g}$ (gas) or $\gamma_{c}$ (cosmic rays)) makes the system more
stable; a softer EoS is destabilizing. Parker’s original assumption,
$\gamma_{c}=0$, makes the system very unstable (Heintz & Zweibel, 2018). Of
course, not all cosmic ray effects can be expressed through an EoS, diffusion
being an obvious exception. We defer a detailed discussion of transport to
Section 2.3, but the main choices are between the self confinement model (also
known as the streaming model), confinement by extrinsic turbulence, which is
diffusive, and the limit of zero diffusion (also known as advection).

As a result of the EoS dependence on transport, the chosen transport model
determines whether the cosmic rays stabilize or destabilize the system
(Zweibel, 2017; Heintz & Zweibel, 2018) In particular, cosmic rays are
stabilizing under the advection model. This result is puzzling, given the
simple physical basis for cosmic ray buoyancy, and partly motivates the work
described in this paper.

The fundamental question connecting research on the Parker instability to
research on star formation is, “Under what conditions do cosmic rays change
the ISM of a galaxy?” To probe this question, we examine how a single
localized injection of cosmic rays (e.g., particles accelerated in the shock
waves created by supernova explosions) disrupts the ISM. We take this route
instead of a traditional Parker instability seeded by small velocity
perturbations because it directly connects our setup to star formation
feedback. Localized cosmic ray injection has been studied previously as a
driver for the galactic dynamo (Hanasz & Lesch, 2000; Hanasz et al., 2004,
2009). It produces structure similar to the Parker instability, with an
extended magnetic lobe and gas condensed in parts of the disk. Focusing on the
sources of cosmic rays instead of the background led us to ask the following:
“What is the timescale on which a single cosmic ray injection changes the
vertical structure of the surrounding ISM?” We break this question into four
sub-questions concerning a single injection’s impact on the ISM:

  * •

(Q1) How does the background medium’s stability change the injection’s impact?

  * •

(Q2) How does the choice of cosmic ray transport model change the injection’s
impact?

  * •

(Q3) How do the injection’s strength and vertical location change its impact?

  * •

(Q4) How does cosmic ray injection differ from heat injection by a thermal
explosion?

By answering these questions for a range of models within an extended
parameter space, our simulations provide a detailed look at how the injection
of cosmic rays in a single event changes the ISM on $\sim 1\,\mathrm{kpc}$
length scales and over times shorter than the $\sim 400\,\mathrm{Myr}$
required for restructuring by the Parker instability (Heintz et al., 2020;
Tharakkal et al., 2022).

Overall, we show that for a stiff cosmic ray EoS, and transport by advection
or diffusion, cosmic ray injection can loft magnetic fields and cosmic rays to
heights exceeding a kiloparsec on timescales $\approx 100\,\mathrm{Myr}$.
Cosmic ray transport dominated by streaming ( i.e. a strong magnetic field)
results in buoyant rising, producing structure similar to the Parker
instability on timescales $\gtrsim 170\,\mathrm{Myr}$. While injecting cosmic
rays in a burst replicates some features of the Parker Instability, the medium
responds on a timescale that is several times faster, and the dependence on
transport is entirely different. Taken together our simulations show that
cosmic ray injection can cause significant dynamical change on a $\sim
100\,\mathrm{Myr}$ timescale, with the change being largest if diffusion is
the dominant cosmic ray transport process, the magnetic field is weak, and the
background cosmic ray pressure is large.

The paper is split into six sections including this motivation section. In
Section 2, we provide background on the Parker instability and cosmic ray
transport. In Section 3 we cover our numerical methods, our initial conditions
for magnetohydrodynamic simulations, and our parameter choices for each
simulation run. In Section 4 we present the overall results of our
simulations. In Section 5 we compare the runs and discuss the implications of
our simulations, while answering our four questions concerning various
astrophysical parameters. In Section 6 we summarize the paper and provide key
takeaways from the work. Readers primarily interested in our results should
focus on Figures 5 & 9, Table 1, and Sections 5 & 6.

In the appendices, we provide further detail on our simulations. In Appendix
A, we compare a Parker instability simulation run with our code to one
obtained previously with a significantly different code and show they give
consistent results. In Appendix B, we discuss the boundary conditions and
dimensionality of our simulations. In Appendix C, we discuss the dependence of
our simulation results on various numerical parameters and our parameter
choices for the primary simulations.

##  2 Background

###  2.1 Magnetohydrostatic Equilibrium

The Parker instability disrupts a magnetohydrostatic equilibrium in the
galactic disk Parker (1966). With a few exceptions (Asseo et al., 1978;
Boulares & Cox, 1990), the equilibrium quantities are assumed to be functions
of $z$ alone and the magnetic field lines are horizontal and straight. The
condition for equilibrium is then

$\frac{d}{dz}\left(P_{g}+P_{c}+\frac{B^{2}}{8\pi}\right)=\rho g.$ (1)

Parker solved the stability problem for a particularly simple class of
equilibria: fixed, constant gravity $g$, constant mean squared random gas
velocity $P_{g}/\rho$, and uniform ratios of nonthermal to thermal gas
pressures parameterized by two constants $\alpha$ and $\beta$333Note Parker’s
$\beta$ is not the usual plasma beta. Instead, the parameter $\alpha$ is the
inverse of the plasma beta $\beta_{\mathrm{pl}}=P_{g}P_{B}^{-1}=\alpha^{-1}$.

$P_{B}=\frac{B^{2}}{8\pi}=\alpha P_{g}\hskip 28.90755ptP_{c}=\beta P_{g}.$ (2)

While other features of the equilibrium setup were modified in later work, the
$(\alpha,\beta)$ parameterization has generally endured, and we use it to
describe the pre-injection state in our simulations. However, we use the
equilibrium setup from Giz & Shu (1993) as our initial ISM background. In this
model the gravitational profile is smooth

$\mathbf{g}(z)=-g_{*}\hat{z}\tanh\left(\frac{z}{H_{*}}\right).$ (3)

instead of the discontinuous $-g_{*}\mathrm{Sign}(z)$ profile Parker (1966)
used. The smooth function in Equation 3 poses fewer numerical difficulties
than a discontinuous profile in numerical simulations which extend above and
below the galactic plane (Heintz et al., 2020). In Equation 3, the asymptotic
vertical gravitational acceleration $g_{*}>0$ and the gravitational scale
height $H_{*}$ depend on the structure of the galactic disk’s stellar
population.

Assuming a plane-parallel hydrostatic equilibrium, Equation 1 becomes

$\frac{dP_{\mathrm{tot}}}{dz}\equiv(1+\alpha+\beta)\frac{dP_{g}}{dz}=\rho(z)g(z)$
(4)

where $P_{\mathrm{tot}}$ is the total pressure; the sum of gas, magnetic, and
cosmic ray pressure. We solve for the equilibrium under the influence of the
gravitational profile Equation 3, using an isothermal equation of state
$P_{g}=c_{s}^{2}\rho$ with constant sound speed $c_{s}$, and midplane values
$P_{g}(0)=P_{g,0},\,\rho(0)=\rho_{0}$. The solution to Equation 4 is

$\frac{P_{g}(z)}{P_{g,0}}=\frac{\rho(z)}{\rho_{0}}=f(z)=\mathrm{sech}^{\eta}{\left[\frac{z}{\eta
H}\right]}.$ (5)

The solution depends on how the scale height of the gas $H$ relates to the
gravitational scale height. The scale height of the gas is

$H=\frac{c_{s}^{2}}{g_{*}}(1+\alpha+\beta).$ (6)

The ratio of the two scale heights $\eta=H_{*}/H$ is a constant. Taking a
limit of Equation 5 as $\eta\rightarrow 0$ (or equivalently, $H_{*}\rightarrow
0$) recovers the solutions in Parker (1966).

The geometry and coordinate system for this equilibrium are illustrated in
Figure 1. The pressures, including magnetic pressure, are homogeneous in
surfaces parallel to the $xy$ plane, which is shown as a green plane in Figure
1. For our simulations, we orient the initial magnetic field in the $\hat{x}$
direction. The magnetic field is shown as a blue arrow in Figure 1. Using
Equations 2 & 5 the magnetic field is

$\mathbf{B}(z)=\hat{x}\sqrt{8\pi\alpha P_{g,0}f(z)}.$ (7)

Similarly, the cosmic ray pressure is

$P_{c}(z)=\beta P_{g,0}f(z).$ (8)

In Figure 1, the coordinates $(\hat{x},\hat{y},\hat{z})$ approximately map to
a galactic disk’s cylindrical coordinates. The $\hat{x}$ direction is parallel
or anti-parallel to the azimuthal direction $\hat{\phi}$ (depending on whether
the magnetic field is oriented clockwise or counter-clockwise around the
galactic center), the $\hat{y}$ direction is parallel or anti-parallel to the
radial direction $\hat{r}$ (again, depending on magnetic field orientation),
and the $\hat{z}$ direction is parallel to the cylindrical $\hat{z}$
direction.

Figure 1: This schematic shows our initial setup for hydrostatic equilibrium.
We have a gravitational field (red arrows) in the $\hat{z}$ direction, which
flips sign at the midplane $z=0$ (see Equation 3). The magnetic field (blue
arrow) is oriented along the $\hat{x}$ direction and creates a magnetic
pressure supporting vertical hydrostatic equilibrium (Equations 4,7). We then
inject cosmic rays, shown by an orange circle (see Section 3.2 for injection
profile), above the midplane, which is shown as a green rectangle.

###  2.2 Parker Instability

We can assess the stability of Parker’s equilibrium to ideal (energy
conserving) small amplitude perturbations with the generalized Schwarzschild
convection criterion (Newcomb, 1961; Boulares & Cox, 1990; Zweibel, 2017).
This criterion for instability in a vertically stratified atmosphere is

$-\frac{d\ln\rho}{dz}<\frac{\rho g_{*}}{\gamma_{g}P_{g}}.$ (9)

This criterion is applicable above the midplane with $z>0$. Assuming we are
well above the midplane $(z>H_{*})$ where the gravitational acceleration is
constant $(-g_{*})$, then the logarithmic derivative is the scale height $H$
of the gas. This assumption gives the same criterion as Parker’s original
gravity profile. In this limit, the criterion for instability to occur is

$\alpha+\beta>\gamma_{g}-1.$ (10)

If we include the compressibility of cosmic rays, then we need to replace
$\gamma_{g}P_{g}$ with $\gamma_{g}P_{g}+\gamma_{c}P_{c}$ in Equation 9. Using
this change and evaluating the criterion again, we get

$\alpha+\beta(1-\gamma_{c})>\gamma_{g}-1.$ (11)

The above criterion illustrates a confusing result noted in Zweibel (2017),
and again in Heintz & Zweibel (2018). Generically, one would guess increasing
non-thermal pressures makes the system more unstable. However, Equation 11
shows there is a way for increasing cosmic ray pressure $(\beta)$ to make the
system more stable. Since the cosmic ray fluid is relativistic gas, it has
$\gamma_{c}=\frac{4}{3}$. In that case, the left hand side would be smaller
with increasing $\beta$. The physical explanation for this effect is that
compressing the cosmic rays requires work.

Self confinement by streaming leads to the polytropic relation
$P_{c}\propto\rho^{2/3}$ along magnetic flux tubes. Using this value in
Equation 11 would predict a stability threshold intermediate between “classic”
Parker ($\gamma_{c}=0$) and $\gamma_{c}=4/3$. Instead, it is shown through a
modal analysis in Heintz & Zweibel (2018) that self confinement leads to a
even larger domain of instability and faster growth rates than “classic”
Parker with $\gamma_{c}=0$. In fact, the analysis used to derive Equation 9
does not apply to the streaming model, due to the relative drift between
cosmic rays and thermal gas, and the heating that accompanies streaming.
Similarly, Equation 9 does not apply to diffusion dominated cosmic ray
transport.

Figure 2: This schematic shows the parameter space of magnetic pressure
$\alpha=P_{B}/P_{g}$ and cosmic ray pressure $\beta=P_{c}/P_{g}$. For
advective transport, we show the regions of parameter space where the initial
hydrostatic equilibrium (see Section 2.2) is unstable to the Parker
instability. The gray region is stable if cosmic rays have an adiabatic index
$\gamma_{c}=4/3$. The black region is stable if the cosmic rays are an
infinitely compressible fluid with index $\gamma_{c}=0$. The S, U, B, and C
points are where the simulations, listed in Table 1, are located in this
parameter space.

We show a schematic for advection dominated transport in Figure 2. This figure
is not quantitatively predictive of our simulation setup, which has a smooth
$g_{*}$ given by Equation 3, a non-ideal perturbation, and in some cases,
transport dominated by streaming or diffusion. The schematic shows the
$(\alpha,\beta)$ space of equilibrium solutions with an adiabatic gas exponent
$\gamma_{g}=5/3$. Using Equations 10 and 11, we show different regions of
stability determined by the cosmic ray adiabatic index $\gamma_{c}$. There is
a strictly Parker unstable region (unshaded), a strictly Parker stable region
(black), and a gray region which is stable if the cosmic rays take work to
compress (assuming $\gamma_{c}=4/3$). For smaller values of $\gamma_{g}$, the
fan of stability boundaries shifts downward so the boundaries intersect on the
vertical axis where $\alpha=\gamma_{g}-1$.

To reiterate: We choose simulation parameters throughout the $(\alpha,\beta)$
parameter space to probe different background medium conditions. This short
review of stability of the background medium illustrates how our different
background medium models would be classified in a Parker instability analysis,
and provides a short review of the Parker instability and its dependence on
cosmic ray transport.

###  2.3 Cosmic Ray Transport

The stability of the system described in Section 2.2 and illustrated in Figure
2 depends on the effective compressibility of the cosmic rays, and therefore
on cosmic ray transport. Cosmic rays are transported throughout the
interstellar medium by three mechanisms: advection, diffusion, and streaming.
The effects of these different transport mechanisms have already been examined
in global galactic models (Girichidis et al., 2016; Wiener et al., 2017).

All three mechanisms are based on the assumption that the cosmic rays exchange
momentum and energy by scattering from fluctuations in the magnetic field with
a wavelength of order the cosmic ray gyroradius (so-called gyroresonant
scattering). Advection, which applies in the limit of infinitely short mean
free path and scattering by fluctuations with no preferred propagation
direction along the background magnetic field lines, is the most
straightforward: if the thermal gas has a bulk flow in a particular direction,
the cosmic rays should flow with the gas. If advection is dominant, then the
cosmic rays behave like a relativistic fluid, with negligible inertia,
following the flow of the thermal (non-relativistic) gas.

Diffusion of cosmic rays applies when the fluctuations again have no preferred
propagation direction, but the mean free path to scattering is large enough to
allow the cosmic rays to leak through the thermal gas. If the measured
residence time and vertical scale height of cosmic rays in the Milky Way are
explained by diffusive propagation, the implied diffusion coefficient is
$\kappa_{c}\sim 3\cdot 10^{28}\,\mathrm{cm}^{2}\,\mathrm{s}^{-1}$.

Streaming, or self confinement, applies when the magnetic fluctuations are
generated by the cosmic rays themselves as a result of directional anisotropy
in the frame of the fluctuations. It can be shown that in the short mean free
path regime, the streaming direction is down the cosmic ray pressure gradient.
This instability appears in the kinetic theory for cosmic ray transport, and
it drives the cosmic rays into a bulk flow at the Alfven speed
$\mathbf{v}_{A}\equiv B/\sqrt{4\pi\rho}$ (Kulsrud & Pearce, 1969; Wentzel,
1969; Skilling, 1975). The streaming instability also heats the gas: the
cosmic rays transfer energy to hydromagnetic waves through gyro-resonance,
which then dissipate energy into the thermal gas (Kulsrud & Pearce, 1969;
Wentzel, 1969; Skilling, 1975). This heating appears as a source term in the
equation for the thermal energy density and the cosmic ray energy density. The
heating rate from the streaming instability is $\mathbf{v}_{A}\cdot\gradient
P_{c}$ (see Zweibel (2017) for discussion and references).

Each of these transport methods drives a different characteristic response in
the ISM when cosmic rays are injected. Advection and streaming both drive
steep fronts of gas and cosmic ray energy. Whereas advection will result in
the cosmic ray fluid moving at the flow speed, bringing the gas with it,
streaming allows the cosmic rays to move ahead of the gas if the Alfven speed
is faster than the sound speed. In the streaming case, a front of cosmic rays
will move through the gas while heating the gas according to the Alfvenic
heating term $\mathbf{v}_{A}\cdot\gradient P_{c}$. Diffusion, in contrast to
both streaming and advection, drives a smoother flow of cosmic rays while also
exerting a force on the gas along the cosmic ray pressure gradient.

We provide the following useful timescales for each transport mechanism to
move cosmic rays through a distance $L$ of the ISM:

$\tau_{\mathrm{adv}}=L\cdot\left(\frac{\gamma_{g}P_{0}}{\rho_{0}}\right)^{-1/2}\approx
65L_{\mathrm{kpc}}T_{4}^{-1/2}\rm{Myr},$ (12)
$\tau_{\mathrm{str}}=L\cdot\left(\frac{B_{0}}{\sqrt{4\pi\rho_{0}}}\right)^{-1}\approx
35L_{\mathrm{kpc}}\rho_{-24}^{1/2}B_{\mu}\rm{Myr},$ (13)
$\tau_{\mathrm{diff}}=\frac{L^{2}}{\kappa_{c}}\approx
30L_{\mathrm{kpc}}^{2}\kappa_{28}^{-1}\rm{Myr},$ (14)

where $L_{\mathrm{kpc}}$, $T_{4}$, $B_{\mu}$, $\rho_{-24}$, and $\kappa_{28}$
are the length scale, gas temperature, gas density (assuming full ionization),
and cosmic ray diffusivity in units of $1\,\mathrm{kpc}$,
$10^{4}\,\mathrm{K}$, $1\mu\mathrm{G}$,
$10^{-24}\,\mathrm{g}\,\mathrm{cm}^{-3}$, and
$10^{28}\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$. In our simulations, the gas flows
are generally subsonic, so the advection timescale is an imperfect measure:
the correct timescale will be larger. For advection dominated simulations, we
set the simulation parameters such that the streaming and diffusion timescales
are much larger than an order unity multiple of the advection timescale in
Equation 12. However, when the advection and streaming timescales are similar
according to Equations 12 & 13, the streaming transport dominates the dynamics
in our simulations because of this overestimated flow speed. We only use these
timescales to determine the dominant transport mechanism, defaulting to
streaming when the streaming and advection timescales are similar (i.e.
Simulation U \- see Table 1).

##  3 Setup

###  3.1 Basic Equations and Numerical Methods

Our simulations are based on fluid equations for the thermal gas - cosmic ray
- magnetic field system Breitschwerdt et al. (1991); Zweibel (2017). These
equations accommodate cosmic ray transport by both diffusion and streaming;
for the latter, it is assumed that the scattering waves propagate down the
cosmic ray pressure gradient. Enforcing this feature creates numerical
difficulties near extrema in cosmic ray pressure. Although this can be handled
through a smoothing length (Sharma et al., 2010), we adopt the two moment
formulation of streaming transport developed by Jiang & Oh (2018) (see also
Thomas et al. (2021)) as implemented in the Athena++ code (Stone et al.,
2020), which modifies the equations presented in Breitschwerdt et al. (1991)
through a time dependent evolution equation for the cosmic ray flux.

The implementation of Athena++ from Jiang & Oh (2018) solves the following
equations:

$\partialderivative{\rho}{t}+\divergence{\left(\rho\mathbf{v}\right)}=0,$ (15)
$\partialderivative{\mathbf{B}}{t}-\curl{\left(\mathbf{v}\crossproduct\mathbf{B}\right)}=0,$
(16)
$\partialderivative{\big{(}\rho\mathbf{v}\big{)}}{t}+\divergence{\left(\rho\mathbf{v}\mathbf{v}+\left(P_{g}+\frac{1}{2}B^{2}\right){\bf\it\mathbb{1}}-\mathbf{B}\mathbf{B}\right)}\\\
=\rho\mathbf{g}+\hat{{\bf\it\sigma}}_{c}\dotproduct\left(\mathbf{F}_{c}-\gamma_{c}\mathbf{v}E_{c}\right),$
(17)
$\partialderivative{E}{t}+\divergence{\left(\mathbf{v}\left(E+P_{g}+\frac{1}{2}B^{2}\right)-\mathbf{B}\left(\mathbf{B}\dotproduct\mathbf{v}\right)\right)}\\\
=\rho\mathbf{g}\dotproduct\mathbf{v}-(\gamma_{c}-1)\left(\mathbf{v}+\mathbf{v}_{s}\right)\dotproduct\gradient{E_{c}},$
(18)
$\partialderivative{E_{c}}{t}+\divergence{\mathbf{F}_{c}}=(\gamma_{c}-1)\left(\mathbf{v}+\mathbf{v}_{s}\right)\dotproduct\gradient{E_{c}},$
(19)
$\frac{1}{V_{m}^{2}}\partialderivative{\mathbf{F}_{c}}{t}+(\gamma_{c}-1)\gradient{E_{c}}=-\hat{{\bf\it\sigma}}_{c}\dotproduct\left(\mathbf{F}_{c}-\gamma_{c}\mathbf{v}E_{c}\right).$
(20)

The dynamical variables are gas density $\rho$, bulk flow velocity
$\mathbf{v}$, gas pressure $P_{g}$, magnetic field $\mathbf{B}$, the combined
internal and kinetic energy density $E=P/(\gamma_{g}-1)+\rho v^{2}/2$, the
cosmic ray energy density $E_{c}=P_{c}/(\gamma_{c}-1)$, and the cosmic ray
energy flux $\mathbf{F}_{c}$. We include the effects of a gravitational
acceleration $\mathbf{g}$, given by Equation 3. We use $\gamma_{c}=4/3$ and
set the modified speed of light parameter to $V_{m}=0.3c$. For more discussion
on convergence according to this parameter, see Appendx C.

###  3.2 Initial Conditions

The initial profiles of density, gas pressure, magnetic field, and cosmic ray
pressure are described in Section 2. We use midplane density and pressure
values of $\rho_{0}=10^{-24}\,\mathrm{g}\,\mathrm{cm}^{-3}$ and
$P_{g,0}=10^{-12}\,\mathrm{g}\,\mathrm{cm}^{-3}$. We use an asymptotic
gravitational acceleration $g_{*}=4\cdot
10^{-9}\,\mathrm{cm}\,\mathrm{s}^{-2}$ and a gravitational scale height
$H_{*}=100\,\mathrm{pc}$ (Heintz et al., 2020). This height is smaller to what
is used in Parker instability simulations (see App. A and Rodrigues et al.
(2016); Heintz et al. (2020)). Additionally, this height decreases our initial
profile’s height ratio $\eta$ when compared to Milky Way and Parker
instability simulations (Rodrigues et al., 2016; Heintz et al., 2020). This
adjustment means our injection can occur in a region where the gravity has
nearly reached the asymptotic value $g_{*}$. As the flux tube of interest
rises, it will move in a constant gravitational acceleration.

The simulation grid is $200\times 50\times 300\,\mathrm{cells}$ (ordered with
respect to $(\hat{x},\hat{y},\hat{z})$) . The cell size is $10\times 20\times
10\,\mathrm{pc}^{3}$, giving a total simulation volume $2\times 1\times
3\,\mathrm{kpc}^{3}$. The third dimension $(\hat{z})$ is off-center of the
midplane ($z=0$), extending from $z=-1\,\mathrm{kpc}$ to $z=+2\,\mathrm{kpc}$.
The cosmic rays are injected above the midplane ($z_{0}>0$), so this off-
centering focuses the computational resources on the injection and the
resulting flows. We do not use any adaptive or static mesh refinement. We
adopt outflow boundary conditions in the $\hat{x}$ and $\hat{y}$ directions
and vacuum (diode) boundary conditions in the $\hat{z}$ direction to minimize
the effect of domain size. See Appendix B for additional discussion.

The equilibrium is setup to within single-precision floating point numerical
error, which is less than the density and pressure floors of our simulations.
There are some waves created by the interaction of the equilibrium with the
vertical boundaries where there is an abrupt transition to vacuum. These waves
are small and have no long term effect on the dynamics of our simulation.

The cosmic ray energy density injection profile is a 3D Gaussian function:

$\varepsilon_{c}(x,y,z)=\frac{E_{\mathrm{SN}}}{\left(2\pi
r_{\mathrm{SN}}^{2}\right)^{3/2}}\cdot\exp\left[\frac{-1}{2r_{\mathrm{SN}}^{2}}\left|\mathbf{r}-\mathbf{r}_{0}\right|^{2}\right].$
(21)

The parameters of the injection are its position $\mathbf{r}_{0}$, radius
$r_{\mathrm{SN}}$, and total energy injected $E_{\mathrm{SN}}$. We add this
perturbation onto the background cosmic ray pressure profile in Equation 5.
Integrating the profile, Equation 21, over volume gives $E_{\mathrm{SN}}$ as
the total energy injected. The injection occurs at different heights above the
midplane depending on the other parameters under consideration (see Table 1).

In each simulation, we use the same injection radius
$r_{\mathrm{SN}}=50\,\mathrm{pc}$. This radius is large for a single supernova
shock, but reasonable for our default energy injection from $\sim 10$
supernovae (see Section 3.3). This larger radius allows our simulation (with
$10-20\,\mathrm{pc}$ resolution) to reliably sample the injection profile. To
avoid sampling errors (aliasing) which will change the total injected energy,
we center the injection in the $xy$ plane on a cell center. The grid has cell
faces along $x=0$ and $y=0$ planes, so we choose to center the injection at
$x_{0}=5\,\mathrm{pc}$ and $y_{0}=10\,\mathrm{pc}$. As a result of this
placement, the peak of the profile occurs in the center of a computational
cell. The height $z_{0}$ and total energy $E_{\mathrm{SN}}$ are varied across
our simulations (see 1)

For reference, if the pressure of $10^{51}\,\mathrm{erg}$ injected cosmic rays
were uniformly distributed throughout the volume of the injection tube (radius
$50\,\mathrm{pc}$, length $2\,\mathrm{kpc}$) is $\sim 0.73\times 10^{-12}$
dynes cm-2, corresponding roughly to a doubling of background cosmic ray
pressure. In this strictly 1D situation, with no horizontal pressure gradient,
the tube would expand radially. This expansion would reduce its density and
cause it to float upward buoyantly, with a characteristic rise time $\sim
100\,\mathrm{Myr}$ (for our adopted gravitational field and scale height).

Parameters |  |  |  | Timescales
---|---|---|---|---
$\rho_{0}=10^{-24}\,\mathrm{erg}\,\mathrm{cm}^{-3}$ | $H_{*}=100\,\mathrm{pc}$ |  | $V_{m}=0.3c$ | $\tau_{\mathrm{adv}}=76\,\mathrm{Myr}\left(\frac{\gamma_{g}}{5/3}\right)^{-1/2}$
$P_{0}=10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-3}$ | $g_{*}=4\cdot 10^{-9}\,\mathrm{cm}\,\mathrm{s}^{-2}$ |  | $\kappa_{c,\perp}=3\cdot 10^{24}\frac{\mathrm{cm}^{2}}{\mathrm{s}}$ | $\tau_{\mathrm{str}}=69\,\mathrm{Myr}\left(\alpha\right)^{-1/2}$
$c_{s}=10^{6}\,\mathrm{cm}\,\mathrm{s}^{-1}$ | $r_{\mathrm{SN}}=50\,\mathrm{pc}$ |  | $\gamma_{c}=4/3$ | $\tau_{\mathrm{diff}}=10\,\mathrm{Myr}\left(\frac{\kappa_{c}}{10^{28}\frac{\mathrm{cm}^{2}}{\mathrm{s}}}\right)^{-1}$
$T=1.21\cdot 10^{4}\,\mathrm{K}\,(\bar{m}/m_{p})\,$ | $B_{0}=5.06\,\mu\mathrm{G}\sqrt{\alpha}$ |  |  |

Run Name | $\alpha$ | $\beta$ | $\gamma_{g}$ | $H\left(\mathrm{pc}\right)$ | $\eta$ | $z_{0}\left(\mathrm{pc}\right)$ | Weight $\left(M_{\odot}\right)$ | $E_{\mathrm{SN}}\left(\mathrm{erg}\right)$ | CR | $\kappa_{c,\parallel}\left(\mathrm{cm}^{2}\,\mathrm{s}^{-1}\right)$ | Question
---|---|---|---|---|---|---|---|---|---|---|---
S | $0.1$ | $0.1$ | $5/3$ | $97$ | $1.03$ | $45$ | $1.28\cdot 10^{4}$ | $10^{51}$ | Adv. | $3\cdot 10^{24}$ | Q1
U | $1.0$ | $1.0$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51}$ | Str. | $3\cdot 10^{24}$ | Q1,Q3,Q4
Uiso | $1.0$ | $1.0$ | $1.1$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51}$ | Str. | $3\cdot 10^{24}$ | Q1
Cadv | $0.1$ | $1.9$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51}$ | Adv. | $3\cdot 10^{24}$ | Q1,Q2
Bstr | $1.9$ | $0.1$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51}$ | Str. | $3\cdot 10^{24}$ | Q1,Q2
Cdiff | $0.1$ | $1.9$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51}$ | Diff. | $3\cdot 10^{28}$ | Q1,Q2
Bdiff | $1.9$ | $0.1$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51}$ | Diff. | $3\cdot 10^{28}$ | Q2
Uheavy | $1.0$ | $1.0$ | $5/3$ | $243$ | $0.41$ | $25$ | $3.27\cdot 10^{4}$ | $10^{51}$ | Str. | $3\cdot 10^{24}$ | Q3
Ublast | $1.0$ | $1.0$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{52}$ | Str. | $3\cdot 10^{24}$ | Q3
Utherm | $1.0$ | $1.0$ | $5/3$ | $243$ | $0.41$ | $105$ | $2.24\cdot 10^{4}$ | $10^{51\,*}$ | Str. | $3\cdot 10^{24}$ | Q4

Table 1: The first table shows the parameters held fixed across our
simulations, as well as derived timescales. The second table shows the
parameters we vary for each simulation run we perform. The first column shows
each simulation run name. The next columns show background medium parameters:
magnetic pressure to gas pressure ratio $\alpha$, cosmic ray pressure to gas
pressure ratio $\beta$, thermal adiabatic index $\gamma_{g}$, gas scale height
$H$, and profile exponent $\eta$ (see Equation 5). The second set of columns
shows injection parameters: initial height of the injection $z_{0}$, the
weight in the column above the injection, and the volume integrated injection
energy $E_{\mathrm{SN}}$ (see Equation 21). The third set of columns shows
dominant cosmic ray transport mechanism and the parallel cosmic ray diffusion
coefficient $\kappa_{c,\parallel}$. The final column shows which questions the
simulation applies to (see Sections 1 & 3.3).

###  3.3 Explanation of Parameters

We present the results of ten simulations of cosmic ray injection. The initial
parameters for each simulation are shown in Table 1. We choose the parameters
for the simulations to help us answer the four questions, presented at the end
of Section 1 and labeled Q1, Q2, Q3, and Q4.

We first considered what values of $\alpha$ and $\beta$ are reasonable in the
solar neighborhood. Measurements suggest the pressures (thermal, cosmic ray,
and magnetic) are nearly equal (Ferrière, 2001). Therefore, we choose a base
run with $\alpha=\beta=1$. Simulations with these values are named with a
prefix U. In Figure 2, this is on a stability boundary when cosmic ray
streaming and diffusion are not included. However, each simulation in Table 1
includes transport by diffusion (even when it is not dominant), and most
include streaming. Therefore, the stability implied by Figure 2 does not
strictly apply, suggesting this equal pressure background medium is Parker
unstable. To compare this background with a strictly Parker stable background,
we run a simulation with $\alpha=\beta=0.1$. We use a prefix ‘S’ to refer to
this combination of values. Our first two simulations, labelled U and S, allow
us to address Q1: how the stability of the background medium changes the
effect of cosmic ray injection. For both these simulations, we use a gas
adiabatic constant consistent with a monatomic ideal gas $\gamma_{g}=5/3$, a
negligible cosmic ray diffusion constant $\kappa_{c}=3\cdot
10^{24}\,\mathrm{cm}^{2}\,\mathrm{s}^{-1}$, and an injection energy equivalent
to the estimated cosmic ray injection energy of $\sim 10$ supernovae (Caprioli
& Spitkovsky, 2014), such as might be expected from an association of coeval
massive stars.

Our third simulation, Uiso, has the same parameters as U except for the
adiabatic gas constant $\gamma_{g}=1.1$, which brings the system closer to the
isothermal case originally considered by Parker (1966) and is expected to be
more unstable. This simulation also addresses Q1.

To answer Q2, how cosmic ray transport affects the behavior of the system post
injection, we need a way to differentiate the dominant cosmic ray transport
mechanism. We do this by considering each transport timescale, given by
Equations 12, 13, 14. It turns out that the flows in our simulations are
subsonic and our advection timescale assumes propagation at the sound speed.
So, when streaming and advection timescales are close, we assume streaming
dominates. Simulations U and S, are instances of streaming and advection
dominance, respectively.

To better isolate the impact of each transport mechanism, we probe two other
points in the $(\alpha,\beta)$ plane of Figure 2. We run a simulation with a
large magnetic field $\alpha=1.9$ and low amount of background cosmic rays
$\beta=0.1$, so it is in the top left corner of Figure 2. This simulation is
labelled Bstr, and streaming dominates the transport of the injected cosmic
rays because the Alfven speed is much higher than the flow speed and the
diffusion rate. We also run a simulation Cadv with $\alpha=0.1$ and
$\beta=1.9$, placing it in the bottom right hand corner of Figure 2. With this
simulation, we probe a medium which is Parker stable due to the stiff equation
of state with large non-thermal pressure where streaming and diffusion are
subdominant transport mechanisms relative to advection.. The simulation Cdiff
has the same values of $\alpha$ and $\beta$, but uses a diffusion constant
$\kappa=3\cdot 10^{28}\,\mathrm{cm}^{2}\,\mathrm{s}^{-1}$ which is close to
estimated Milky Way values.

While these three simulations allow us to examine how streaming, advection,
and diffusion affect the response of the medium to cosmic ray injection, they
do not completely eliminate the influence of the background state; while the
nonthermal pressures are equal in all three, the role of magnetic tension in
Bstr is larger than in Cadv or Cdiff due to the stronger background magneic
field. Therefore, to determine whether the streaming cosmic rays or the strong
magnetic field was more critical, we ran Bdiff. That simulation had the same
initial parameters as Bstr, but with a higher diffusion coefficient and
streaming terms removed. Although diffusion and advection do not have an ‘off-
switch’ in the implementation from Jiang & Oh (2018), streaming transport
does, and we make use of it when running Bdiff.

Bstr and Cadv also address Q1 because they are at different points in the
$(\alpha,\beta)$ plane. They have different background cosmic ray and magnetic
pressures from simulations U or S, allowing us to learn how that aspect of the
background medium changes the impact of cosmic ray injection.

Question 3 probes the effects of injection properties, so we adjust parameters
related to the injection. Using simulation U as a control case, we compare
both of the following simulations to it, and only change single parameters.
First, we use simulation Uheavy with the cosmic ray injection at a lower
height in the disk, meaning it has more weight above it. Second, we use a
simulation, Ublast, with the strength of $10\%$ of $100$ supernovae to
understand how different amounts of energy injection change the results. Other
than the weight above the injection and the injection energy, we keep the
variables the same as simulation U, which we compare each of these cases
against.

For Q4, which compares the effects of cosmic rays to direct energization of
thermal gas, we run a simulation, Utherm, with the same initial parameters as
the Parker unstable simulation U, except the injection is in thermal pressure
instead of in cosmic ray pressure.

Figure 3: After injecting cosmic rays on a magnetic flux tube, we track that
tube’s movement. For each simulation shown, the solid black line is maximum
height of the flux tube at a given time. The green dotted line is the average
mass in the region $x<|250\,\mathrm{pc}|$, as a fraction of the initial
average mass. The purple solid, dashed, and dash-dotted lines are the gas,
magnetic, and cosmic ray pressure, respectively, averaged in the same region.
The second $y$ axis on the right side applies to the lines of mass and
pressures, with each quantity $Q\in\\{\langle M\rangle,\langle
P_{g}\rangle,\langle P_{B}\rangle,\langle P_{c}\rangle\\}$ normalized by
either $P_{0}=10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-3}$ or the initial value
(for mass $\langle M\rangle$).

##  4 Results

We focus on a flux tube enclosing the initial injection. Since we use the
ideal MHD equations and because there is effectively no perpendicular cosmic
ray diffusion, the gas and cosmic rays stay on the flux tube. All the
simulations share some common features. The cosmic rays in the injection
sphere are overpressured by more than a factor of 20. This creates a low
density cavity within a few Myr and launches a pressure driven flow away from
the injection sphere. The cavity itself is buoyant because of reduced density
and the outflow reduces the density along a progressively longer segment of
the tube. The result is a rising magnetized arch which becomes ever more
buoyant as gas is drained from the tube by a gravity driven downflow. The only
forces that counter the rise are magnetic tension due to field line curvature,
the inertia of the overlying gas, and adiabatic cooling of both the thermal
and cosmic ray gases as they expand. The configuration has strong up-down
asymmetry, because the effects which slow the flux tube’s rise increase near
the midplane.

For each simulation, we examine physical quantities near the center of the
tube $(x=\pm 250\,\mathrm{pc})$. This restriction allows us to focus on the
dynamics in the center of the flux tube, where the initial cosmic ray
injection occurred. Additionally, it minimizes the impact of the boundary
conditions on our results for early times in the simulation (See Appendix B
for a discussion of boundary conditions). In Figure 3 we show the maximum
height along the flux tube as a solid black line, the average gas pressure as
a solid purple line, the average magnetic pressure as a dashed purple line,
the average cosmic ray pressure as a dash-dotted purple line, and the average
mass as a dotted green line. The average flux tube height in the $\pm
250\,\mathrm{pc}$ region is approximately the same as the maximum height.
These quantities help us understand how the injection evolves through time and
how it causes the flux tube to change.

The most striking part of the results in Figure 3 are the plots for
simulations Cadv, Cdiff, Ublast, and Utherm. In these four simulations, the
flux tubes rise over $400\,\mathrm{pc}$ in less than $100\,\mathrm{Myr}$. The
rise is accompanied by decreases in mass and pressure in the central region of
the flux tube.

The other simulations do not exhibit such a quick, drastic rise in flux tube
height, nor decreases in pressure and mass in the central region of the flux
tube. Simulations S and Uheavy end up being stable to the perturbation
(meaning the flux tube rise is limited). S even finds a new equilibrium after
an initial rise. Simulations U, Uiso, and Bstr all eventually begin to
buoyantly rise, but only rising to $200\,\mathrm{pc}$, at a time of
$150\,\mathrm{Myr}$ after the cosmic ray injection. These simulations each
have streaming as the dominant cosmic ray transport - instead of creating a
violent disruption in the center of the magnetic flux tube, streaming puts
energy into heating the entire length of the magnetic flux tube. However, the
stronger injection in Ublast still launches material rapidly when compared to
U, even though streaming is dominant.

To more clearly illustrate the dynamics in our simulations, we show 2D cuts
from simulation Cdiff in Figure 4. The first column shows gas density, the
second column shows vertical momentum, and the third column shows horizontal
momentum. The green lines are magnetic field lines along the flux tube
analyzed in Figure 3. The first row is a snapshot at $4\,\mathrm{Myr}$, the
second row at $60\,\mathrm{Myr}$, and the third row at $120\,\mathrm{Myr}$.
Once the cosmic rays diffuse, they create an over-pressured flux tube which
begins to move upward (middle column, top row). After the flux tube bends, gas
begins to fall down the curved lines (right column, middle row). After some
time, this process builds until a large up-flow in the center and down-flow
along the magnetic field lines. The final density plot (left column, bottom
row) also shows how the flux tube is able to lift some mass out of the disk.

We show the cosmic ray pressure distributions for simulations Cadv, Bstr, and
Cdiff in Figure 5. The first column shows simulation Cadv, the second column
shows Bstr, and the third column shows Cdiff. The first and second rows show
the same time dumps for each simulation (top row is $3\,\mathrm{Myr}$, middle
row is $50\,\mathrm{Myr}$). The final row shows the later development of the
injection for each simulation. The Bstr simulation does not grow as much
vertically because magnetic tension holds the flux tube down for most of the
simulation. Only once the field lines are able to bend does the flux tube
begin to rise, because gas is able to leave the flux tube at a faster rate.
This effect is also shown in the second row of plots in Figure 3. Regardless
of cosmic ray transport by streaming or diffusion, the strong magnetic field
simulations (Bstr and Bdiff) exhibit slowed or negligible rise of the flux
tube. Streaming appears to be the most effective in the large magnetic field
case because the gas gets heated in addition to being overpressured.

Figure 4: Plots of gas density (first column), vertical momentum (second
column), and horizontal ($\hat{x}$) momentum (third column) at a $y=0$ cut of
simulation Cdiff. In this simulation, cosmic ray diffusion is the dominant
transport mechanism. Early in the simulation, the cosmic rays quickly diffuse
along the magnetic field lines. This diffusion leaves behind a flux tube with
increased average cosmic ray pressure. This region begins to rise buoyantly,
shown in the middle row. The final row shows gas falling back toward the disk
around the injection, while the location continues to launch gas upward at
$x=0$. The normalization constants are
$P_{0}=10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-3}$ and
$\rho_{0}c_{s}=10^{-18}\,\mathrm{g}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}$. The
green lines are magnetic field lines; the middle line is integrated from the
$x=-1\,\mathrm{kpc}$ bound at the height of the initial injection. The other
lines are initialized at $\pm r_{\mathrm{SN}}$ and $\pm 2r_{\mathrm{SN}}$ of
the injection height. Figure 5: Plots of cosmic ray pressure at a $y=0$ cut of
simulations Cadv (first column), Bstr (second column), and Cdiff (third
column). The normalization constant is
$P_{0}=10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-3}$. The green lines are magnetic
field lines. The dashed white lines are contours of equal cosmic ray pressure.
In the first two rows, cuts are all from the same time snapshots, but the
final row shows simulation Bstr at a later time snapshot than the other two
simulations.

###  4.1 A Caveat on Simulation _Utherm_

In simulation Utherm, instead of injecting cosmic rays, we injected energy as
thermal heating. This thermal energy injection had the same magnitude,
$10^{51}\,\mathrm{erg}$, as the cosmic ray injection in simulation U.

While this simulation highlights the effectiveness of thermal energy
injection, our simulations lack radiative cooling. This injected gas should
cool by emitting radiation while it expands. We estimate a radiative cooling
time scale to determine the effect of radiative cooling on the injection,
using Equation 34.4 from Draine (2011):

$t_{\mathrm{cool}}=1.1\cdot
10^{5}\,\mathrm{yr}\left(\frac{T}{10^{6}\,\mathrm{K}}\right)^{1.7}\left(\frac{n_{\mathrm{H}}}{\mathrm{cm}^{3}}\right)^{-1}$
(22)

For simulation Utherm, the peak temperature in the center of the injection is
$10^{5.04}\,\mathrm{K}$, compared to the background $T_{0}\sim
10^{4}\,\mathrm{K}$ Therefore, the radiative cooling time for this injection
is $t_{\mathrm{cool}}\approx 2.2\,\mathrm{kyr}$. This time scale is so short
that this thermal injection should dissipate before causing a disruption in
the ISM like we see in simulation Utherm. While temperatures are even higher
in actual supernova remnants ($\sim 10^{7}\,\mathrm{K}$ for the Sedov-Taylor
phase (Draine, 2011)), our injection represents a long term average impact of
those remnants. Radiative cooling will limit the impact of the thermal
injection, as compared to a cosmic ray injection. A similar effect is seen in
models of supernova heated gas launched from galactic disks, which tend to be
fountains rather than winds (Shapiro & Field, 1976; Rosen & Bregman, 1995;
Bustard et al., 2018, 2020).

##  5 Discussion

To the questions posed in Section 1, we now have the following answers:

  * •

(Q1) A purely Parker stable medium limits the effect of cosmic ray injection.
An isothermal medium is more prone to disruption by a cosmic ray injection.
See Figure 6.

  * •

(Q2) Cosmic ray advection and diffusion drive changes in the ISM on $\sim 100$
Myr time scales, through explosive launching and buoyancy, respectively. Large
magnetic field strength (which also implies streaming dominated transport)
delays the flux tube’s rise, but eventually causes buoyant rising. In the
large magnetic field case, streaming is more disruptive than diffusion. See
Figure 7.

  * •

(Q3) Stronger injections drive more explosive flows, and injections closer to
the midplane take longer to launch the flux tube. See Figure 8.

  * •

(Q4) Thermal injection drives buoyant rising of the flux tube on a shorter
time scale than cosmic ray injection. Cosmic rays decrease the average mass
along the flux tube at a slower rate. Eventually, cosmic ray injections
overtake thermal injections in height. However, simulation Utherm
overestimates the efficacy of thermal injection because it lacks cooling (see
Section 4.1). See Figures 8, 9.

In the following subsections, we provide more complete explanations for these
answers.

###  5.1 Dependence on Background Medium (Q1)

We explored how changes in the background medium affected the evolution of a
cosmic ray injection. The results are shown in Figure 6, which plots the
average mass and height along the flux tube of each simulation.

Figure 6: Plot of maximum tube height (black) and average mass (green) within
$x=\pm 250\,\mathrm{pc}$ of the injection location for simulations related to
Q1. Note the vertical axis has been extended compared to Figure 3. Solid lines
show simulation S, dashed lines show simulation U, dash-dotted lines show
simulation Uiso, dotted lines show simulation Bstr and dots show simulation
Cadv. The streaming dominated simulations, U and Bstr, rise late. The
advection dominated simulations, Cadv and S, have different results: the
medium with more cosmic ray pressure is significantly disrupted. Simulation
Uiso with $\gamma_{g}=1.1$ allows more mass loss later in the simulation, when
compared with simulation U, which had $\gamma_{g}=5/3$.

The Parker stable simulation (Simulation S) weathers the cosmic ray injection,
at least under transport by advection and streaming, reaching a new stable
equilibrium. When the background is Parker unstable according to the original
criterion (Equation 10), the injection drives significant change. For the
simulations U, Bstr, and Cadv, we can use Equation 10 to show the combined
nonthermal pressures dominate the gas compressibility
$(\alpha+\beta=2)>(\gamma_{g}-1=2/3)$. The simulations with a background
biased towards magnetic pressure (Bstr) or cosmic ray pressure (Cadv) launch
the flux tube faster than the equipartition case (U). Therefore, the initial
nonthermal pressures, given by $(\alpha,\beta)$, determine how prone the
system is to disruption by cosmic ray injection. This dependence also suggests
an ISM with non-uniform $\alpha,\beta$ would provide a complex environment for
cosmic ray injection, since different directions and positions could be more
(or less) prone to being disrupted.

The case of an isothermal-like background with $\gamma_{g}=1.1$ responds on a
slightly shorter timescale than a background medium with $\gamma_{g}=5/3$.
This conclusion comes from considering the average height achieved by the flux
tubes in simulations Uiso and U. In simulation Uiso, buoyant rising begins
after $80\,\mathrm{Myr}$ and stays significant through the end of the
simulation. For simulation U, buoyant rising becomes dominant only after
$140\,\mathrm{Myr}$. The average mass also decreases at a faster rate at late
times in simulation Uiso. Because $\gamma_{g}$ is lower, it takes less work to
compress the gas and push it off the flux tube once buoyancy kicks in. This
change allows mass to flow at a faster rate, driving the buoyancy force to
become larger.

###  5.2 Dependence on Cosmic Ray Transport (Q2)

We ran four simulations to determine the impact of cosmic ray transport:
simulation Bstr is streaming dominated, simulation Cadv is advection
dominated, Cdiff is diffusion dominated, and Bdiff is diffusion dominated with
a strong magnetic field. With these simulations, it is clear that streaming
does a poor job of launching material when compared to diffusion and
advection, which both disrupt the ISM on a short time scale $\sim
100\,\mathrm{Myr}$. However, Bdiff shows this delay is not a result of the
cosmic ray transport by streaming. Instead, the large magnetic field
(necessary for streaming to be dominant) resists any bending created by the
injection. Streaming is more effective than diffusion at disrupting the
vertical structure when magnetic tension is a dominant force. Figure 7 shows
the flux tube averaged mass and height for these three simulations.

Figure 7: Plot of maximum tube height (black) and average mass (green) within
$x=\pm 250\,\mathrm{pc}$ of the injection location for simulations related to
Q2. Note the vertical axis has been extended compared to Figure 3. Solid lines
show simulation Cdiff, dashed lines show simulation Cadv, dash-dotted lines
show simulation Bstr, and dotted lines show simulation Bdiff. The streaming
case has the flux tube rise buoyantly, but only after a long time $\gtrsim
100\,\mathrm{Myr}$. In contrast, the advection and diffusion dominated cases
rise early $\lesssim 100\,\mathrm{Myr}$ and cause a significant amount of mass
loss from the flux tube. The streaming case is primarily inhibited by the
magnetic field strength, not the streaming transport, because in simulation
Bdiff the flux tube takes an even longer time to rise.

Simulation Bstr does not start to rise buoyantly until a significant amount of
mass has been lost from the flux tube. In comparison, the buoyant rise in
simulation Cdiff happens quickly (see panels of Figure 3). The advective
simulation Cadv is mainly driven by an explosive launching, instead of
buoyancy. However, it and the diffusion case produce similar results in terms
of flux tube movement. The difference between those two simulations is more
apparent in the beginning, when the launching is different.

Of these simulations, the diffusion case is the most surprising. Initially,
one might think cosmic ray diffusion may have a smaller effect on feedback
processes than streaming and advection, because there is less time for cosmic
rays to impact the ISM (See Equation 14) when using the Milky Way value of the
diffusion coefficient. However, our simulation Cdiff shows a large cosmic ray
injection can generate enough force through the cosmic ray pressure gradient
to move mass out of the disk and bend the magnetic field.

###  5.3 Dependence on Injection Characteristics (Q3 & Q4)

We ran three simulations focused on the injection characteristics. The first,
Uheavy, placed the injection lower in the galactic disk. While the flux tube
rose later than in other simulations, this simulation gave similar results
when compared to U. Simulation Ublast examined how increasing the injection
energy to $10^{52}\,\mathrm{erg}$ would change the dynamics. This large
injection caused significant change in the ISM, causing the flux tube to rise
rapidly. This simulation suggests that future simulations, with multiple
injections, could lead to rapid disruption. In Simulation Utherm, we replaced
the cosmic ray injection with a thermal injection of similar magnitude. This
injection contained no mass, but was an over-pressure $(T\sim
10^{5}\,\mathrm{K}>T_{0}\sim 10^{4}\,\mathrm{K})$ region. This thermal energy
injection caused rapid mass flow off the flux tube, driving a buoyant rise.
However, the simulation is physically unrealistic because we do not include
radiative cooling (see Section 4.1). The tube averaged mass and height of
these simulations are shown as part of Figure 8.

Figure 8: Plot of maximum tube height (black) and average mass (green) within
$x=\pm 250\,\mathrm{pc}$ of the injection for simulations related to Q3 and
Q4. Note the vertical axis has been extended compared to Figure 3. Solid lines
show simulation U, dashed lines show simulation Uheavy, dash-dotted lines show
simulation Ublast, and dotted lines show simulation Utherm. Simulation Uheavy
does eventually rise from the disk after $200\,\mathrm{Myr}$. Simulation
Ublast drives an explosive launching faster than the thermal injection in
Utherm.

Figure 9:  Top Row, left plot: The evolution over time of each flux tube’s
height and mass. Solid lines and filled markers correspond to the tube’s peak
height. Each simulation is shown until $130\,\mathrm{Myr}$, with markers
placed at $10\,\mathrm{Myr}$ increments. Black circles show simulation Cadv,
red stars show simulation Cdiff, and blue squares show simulation Utherm.
While the thermal injection in simulation Utherm initially moves material up
faster, at late times the cosmic ray injections in simulations Cadv and Cdiff
overtake simulation Utherm. After $130\,\mathrm{Myr}$, the diffusion
transported cosmic ray injection has created the most vertically extended flux
tube. Top Row, Right Plot: For the same simulations, we show the average
momentum in the vertical direction in the region of $z>1.5\,\mathrm{kpc}$.
Cdiff creates the largest mass flux. The decrease around $140\,\mathrm{Myr}$
is when the flux tube leaves the simulation through the upper boundary. Bottom
Row: Cuts of simulations Cadv, Cdiff, and Utherm. These simulations were the
fastest at disrupting the vertical structure, with each image shown at a time
of $130\,\mathrm{Myr}$. The black lines are magnetic field lines and the white
dashed lines are density contours. These plots are cuts of our 3D simulations
at $y=0$. Simulation Cdiff reaches the highest point, with significantly
adjusted magnetic field structure.

###  5.4 Fast Disruption: Cosmic rays vs. Thermal Injection

The simulations which produced the largest changes in the shortest time for an
injection energy $10^{51}\,\mathrm{erg}$ were Cadv, Cdiff, and Utherm. These
injections had flux tubes rise faster than a traditional Parker Instability,
and caused the flux tube to lose approximately half its mass in $\lesssim
100\,\mathrm{Myr}$. In the top left panel of Figure 9 we show the evolution of
peak tube height against mass. The markers on each line denote
$10\,\mathrm{Myr}$ steps, whereas each line is made with points at a
resolution of $1\,\mathrm{Myr}$. The final marker is at $130\,\mathrm{Myr}$
for each simulation. Even though the cosmic ray injection simulations take a
longer time to start rising, they rise at a faster rate than the thermal
injection. They overtake the thermal injection after approximately
$90\,\mathrm{Myr}$, and they have more mass at that time. The density at $y=0$
of the three simulations at this final time, $130\,\mathrm{Myr}$, is shown in
the bottom panel of Figure 9. Each simulation produces a different magnetic
field structure. Considering simulation Utherm is an overestimate of thermal
injection’s impact (see cooling time argument in Section 4.1), we can focus on
the difference between Cadv and Cdiff. The left and middle plots in the lower
row of Figure 9 show that cosmic ray transport changes the both the flow of
gas around the rising flux tube and the shape of the magnetic field lines. The
top right plot of Figure 9 shows the average vertical momentum of the cells
with $z>1.5\,\mathrm{kpc}$ for each of the three simulations. Following the
red dashed line for simulation Cdiff, it is clear the diffusion dominated
simulation creates the largest mass flux from the disk. While this gas may
fall down given enough time and vertical expansion distance, we are unable to
continue to follow that gas because it leaves out of the top of our simulation
$(z=2\,\mathrm{kpc})$.

##  6 Conclusion

We ran ten simulations of cosmic ray injection into a vertically stratified
medium, using the Athena++ code. These simulations illustrated the effect of
cosmic ray injection in a galactic disk, on intermediate scales $(\sim
1\,\mathrm{kpc})$ larger than the ISM’s fine structure and smaller than the
entire galaxy. By exploring an extended parameter space, we produced the
highest resolution picture yet of localized cosmic ray injection on these
mesoscales. We also showed that cosmic ray transport dictates the impact of
cosmic rays on the ISM. Below, we provide the key points and results of this
work:

  * •

Cosmic ray diffusion of a cosmic ray injection can change the ISM’s vertical
structure and a galaxy’s magnetic field on timescales $<100\,\mathrm{Myr}$.

  * •

The large magnetic field strength necessary for cosmic ray streaming dominance
over diffusion limits the rapid disruption of the flux tube because of
increased magnetic tension. Streaming is effective at disrupting the vertical
structure, but takes a longer time than in weak magnetic field cases.

  * •

A flux tube disrupted by cosmic ray injection will rise faster at later times
than one disrupted by thermal injection, producing a larger mass flow of
material out of the galactic disk.

Our simulations provided useful heuristic results by considering cosmic ray
injection in a stratified medium. These simulations are clearly limited in
their realism because galactic disks are not uniformly stratified. The ISM in
galactic disks is multiphase and clumpy. Additionally, by neglecting the
dynamics of the stellar gravitational potential, we lack forcing terms which
could change the effects of these energy injections. Future work may need to
focus on diffusion as the primary cosmic ray transport mechanism, along with
implementation of more realistic ISM conditions (multiple gas phases, galactic
rotation, etc.). The non-constant $\alpha$, $\beta$ of the injection have a
significant impact. This non-uniformity should be extended to the background
medium, instead of using the constant $\alpha,\beta$ assumption originally
introduced by Parker (1966). Variable magnetic field strength is particularly
important, as it could amplify the importance of cosmic ray streaming in low
$\beta_{\mathrm{pl}}$ regions.

In these simulations, we neglected radiative heating and cooling, but those
processes are an important consideration in the ISM. The radiative cooling
would be particularly important for the gas heated by cosmic ray streaming.
Cooling will minimize, and possibly remove, the impact of the thermal
injection in simulation Utherm. Additionally, having only a single cosmic ray
injection by supernovae in our simulation volume over a time of
$100\,\mathrm{Myr}$ is unrealistic. In future work, we aim to include multiple
injections at different times and locations. Separating the injections in
space opens up a variety of other situations which make it difficult to
isolate the buoyancy the injection creates, hence why we only considered a
single injection location in this work.

For large scale (cosmological or galactic) simulations, there needs to be some
consideration of cosmic ray injection. Even simulations which include cosmic
rays generally do not resolve their injection into the ISM. Our work shows a
cosmic ray injection can drive an upward flow in $<100\,\mathrm{Myr}$ after
their injection, and most of that upward movement actually happens in
$<20\,\mathrm{Myr}$. A resolution of $\sim 100\,\mathrm{pc}$ would be enough
to illustrate the impact of an upward flow along the magnetic field (consider
bottom row of Figure 9). A sub-grid physics module related to cosmic ray
injection could initiate this flow before letting it evolve independently. In
the $y$ direction, across the magnetic field and in the plane of the disk,
higher resolution would be necessary. The width of the flux tube barely
reaches $100\,\mathrm{pc}$ in that direction, and the overall dynamics in that
direction are minimal. The main impact the $y$ direction had in our
simulations is to allow mass above the flux tube to move out of the path of
the rising flux tube. The actual resolution of that dimension is less
significant (See Appendix B).

Our work also shows cosmic ray injection is an important part of dynamics in a
galaxy. Moving, heating, and compressing gas all have an impact on where stars
form and on galactic structure. Since the eventual rise of the flux tube
happens in a short time $(<20\,\mathrm{Myr})$, the dynamics created by cosmic
ray injection are less likely to be washed out by galactic rotation and other
dynamical processes. Our results also shed light on the respective roles of
cosmic ray compressibility and buoyancy alluded to in Section 2.2. Finally,
our results suggest that including effects of localized cosmic ray injection
in global simulations could be a necessary step in accurately modeling
galactic outflows and evolution.

##  7 Acknowledgements

We would like to thank Yan-Fei Jiang for sharing the code for cosmic rays in
Athena++, originally presented in Jiang & Oh (2018). We are very thankful for
the comments provided by an anonymous reviewer. Their suggestions
significantly improved this work. We would like to thank Evan Heintz and Chad
Bustard for helpful discussions about our simulations and the Parker
instability. We also appreciate commentary from Ryan Farber and Cassi Lochhaas
which improved the paper. This work was funded by NSF grant AST-2007323 and
NASA FINESST Grant 80NSSC22K1749.

## References

  * Asseo et al. (1978) Asseo, E., Cesarsky, C. J., Lachieze-Rey, M., & Pellat, R. 1978, ApJ, 225, L21, doi: 10.1086/182784
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bell (2003) Bell, E. F. 2003, ApJ, 586, 794, doi: 10.1086/367829
  * Boulares & Cox (1990) Boulares, A., & Cox, D. P. 1990, ApJ, 365, 544, doi: 10.1086/169509
  * Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Voelk, H. J. 1991, A&A, 245, 79
  * Bustard et al. (2018) Bustard, C., Pardy, S. A., D’Onghia, E., Zweibel, E. G., & Gallagher, J. S., I. 2018, ApJ, 863, 49, doi: 10.3847/1538-4357/aad08f
  * Bustard & Zweibel (2021) Bustard, C., & Zweibel, E. G. 2021, ApJ, 913, 106, doi: 10.3847/1538-4357/abf64c
  * Bustard et al. (2020) Bustard, C., Zweibel, E. G., D’Onghia, E., Gallagher, J. S., I., & Farber, R. 2020, ApJ, 893, 29, doi: 10.3847/1538-4357/ab7fa3
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, ApJ, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Chan et al. (2019) Chan, T. K., Kereš, D., Hopkins, P. F., et al. 2019, MNRAS, 488, 3716, doi: 10.1093/mnras/stz1895
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press)
  * Everett et al. (2008) Everett, J. E., Zweibel, E. G., Benjamin, R. A., et al. 2008, ApJ, 674, 258, doi: 10.1086/524766
  * Farber et al. (2018) Farber, R., Ruszkowski, M., Yang, H. Y. K., & Zweibel, E. G. 2018, ApJ, 856, 112, doi: 10.3847/1538-4357/aab26d
  * Ferrière (2001) Ferrière, K. M. 2001, Reviews of Modern Physics, 73, 1031, doi: 10.1103/RevModPhys.73.1031
  * Fryxell et al. (2000) Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273, doi: 10.1086/317361
  * Girichidis et al. (2016) Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19, doi: 10.3847/2041-8205/816/2/L19
  * Giz & Shu (1993) Giz, A. T., & Shu, F. H. 1993, ApJ, 404, 185, doi: 10.1086/172267
  * Hanasz et al. (2004) Hanasz, M., Kowal, G., Otmianowska-Mazur, K., & Lesch, H. 2004, ApJ, 605, L33, doi: 10.1086/420697
  * Hanasz & Lesch (2000) Hanasz, M., & Lesch, H. 2000, ApJ, 543, 235, doi: 10.1086/317077
  * Hanasz et al. (2013) Hanasz, M., Lesch, H., Naab, T., et al. 2013, ApJ, 777, L38, doi: 10.1088/2041-8205/777/2/L38
  * Hanasz et al. (2009) Hanasz, M., Wóltański, D., & Kowalik, K. 2009, ApJ, 706, L155, doi: 10.1088/0004-637X/706/1/L155
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Heintz et al. (2020) Heintz, E., Bustard, C., & Zweibel, E. G. 2020, ApJ, 891, 157, doi: 10.3847/1538-4357/ab7453
  * Heintz & Zweibel (2018) Heintz, E., & Zweibel, E. G. 2018, ApJ, 860, 97, doi: 10.3847/1538-4357/aac208
  * Hopkins et al. (2021) Hopkins, P. F., Chan, T. K., Squire, J., et al. 2021, MNRAS, 501, 3663, doi: 10.1093/mnras/staa3692
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jiang & Oh (2018) Jiang, Y.-F., & Oh, S. P. 2018, ApJ, 854, 5, doi: 10.3847/1538-4357/aaa6ce
  * Kulsrud & Pearce (1969) Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445, doi: 10.1086/149981
  * Mouschovias (1974) Mouschovias, T. C. 1974, ApJ, 192, 37, doi: 10.1086/153032
  * Newcomb (1961) Newcomb, W. A. 1961, Physics of Fluids, 4, 391, doi: 10.1063/1.1706342
  * Parker (1966) Parker, E. N. 1966, ApJ, 145, 811, doi: 10.1086/148828
  * Rodrigues et al. (2016) Rodrigues, L. F. S., Sarson, G. R., Shukurov, A., Bushby, P. J., & Fletcher, A. 2016, ApJ, 816, 2, doi: 10.3847/0004-637X/816/1/2
  * Rosen & Bregman (1995) Rosen, A., & Bregman, J. N. 1995, ApJ, 440, 634, doi: 10.1086/175303
  * Ruszkowski et al. (2017) Ruszkowski, M., Yang, H. Y. K., & Zweibel, E. 2017, ApJ, 834, 208, doi: 10.3847/1538-4357/834/2/208
  * Shapiro & Field (1976) Shapiro, P. R., & Field, G. B. 1976, ApJ, 205, 762, doi: 10.1086/154332
  * Sharma et al. (2010) Sharma, P., Colella, P., & Martin, D. F. 2010, SIAM Journal on Scientific Computing, 32, 3564, doi: 10.1137/100792135
  * Skilling (1975) Skilling, J. 1975, MNRAS, 172, 557, doi: 10.1093/mnras/172.3.557
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, The Astrophysical Journal Supplement Series, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Tharakkal et al. (2022) Tharakkal, D., Shukurov, A., Gent, F. A., et al. 2022, arXiv e-prints, arXiv:2212.03215, doi: 10.48550/arXiv.2212.03215
  * Thomas et al. (2021) Thomas, T., Pfrommer, C., & Pakmor, R. 2021, MNRAS, 503, 2242, doi: 10.1093/mnras/stab397
  * Uhlig et al. (2012) Uhlig, M., Pfrommer, C., Sharma, M., et al. 2012, MNRAS, 423, 2374, doi: 10.1111/j.1365-2966.2012.21045.x
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Wentzel (1969) Wentzel, D. G. 1969, ApJ, 156, 303, doi: 10.1086/149965
  * Wiener et al. (2017) Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 906, doi: 10.1093/mnras/stx127
  * Zhang (2018) Zhang, D. 2018, Galaxies, 6, 114, doi: 10.3390/galaxies6040114
  * Zweibel (2017) Zweibel, E. G. 2017, Physics of Plasmas, 24, 055402, doi: 10.1063/1.4984017
  * Zweibel & Bruhwiler (1992) Zweibel, E. G., & Bruhwiler, D. L. 1992, ApJ, 399, 318, doi: 10.1086/171927
  * Zweibel & Kulsrud (1975) Zweibel, E. G., & Kulsrud, R. M. 1975, ApJ, 201, 63, doi: 10.1086/153858

Figure 10: Top row: Vertical profiles of horizontally averaged gas density
(solid red line), cosmic ray pressure (dash-dotted blue line), and magnetic
pressure (dotted green line). The initial vertical profile for the quantities
is shown as well (dashed black line). The initial profile is the same for each
quantity because we normalize by the initial midplane value of each quantity,
and each physical variable is proportional to the other (see Section 2).
Bottom row: Gas density at corresponding time dump to the vertical profile
above each image. The total mass in the simulation is conserved, with a loss
of $\sim 0.1\%$ through the diode boundary conditions in the vertical
$(\hat{z})$ directions.

##  Appendix A Comparison with Parker Instability Simulations

To test our numerical methods and initial conditions, we ran a 2D simulation
matching one in Heintz et al. (2020). Those simulations used the FLASH code
(Fryxell et al., 2000), along with the streaming transport method in Sharma et
al. (2010). Matching those parameters (originally based on Rodrigues et al.
(2016)), we use $\rho_{0}=6.76\cdot 10^{-25}\,\mathrm{g}\,\mathrm{cm}^{-3}$,
$P_{g,0}=8.19\cdot 10^{-13}\,\mathrm{erg}\,\mathrm{cm}^{-3}$, $g_{*}=2\cdot
10^{-9}\,\mathrm{cm}\,\mathrm{s}^{-2}$, $\eta=2$, $H=250\,\mathrm{pc}$,
$\alpha=10/3$, $\beta=1.25$, and $V_{m}=0.01c$.

We find that that our simulation method is consistent with that of Heintz et
al. (2020), and our results are illustrated in Figure 10. This figure shows
horizontally $(\hat{x})$ averaged profiles of the gas density, cosmic ray
pressure, and magnetic pressure at several time dumps, compared with the
original profile for the quantities. Below the averaged profiles, we show the
gas density at each selected time dump. The magnetic and cosmic ray pressures
decrease more slowly away from the plane as the magnetic field becomes bent,
while gas is compressed towards the midplane. Gravity pulls the gas down along
the magnetic flux tubes, which have turned perpendicular to the disk in
several locations. Our simulation matches expected behavior (see Figure 11 of
Heintz et al. (2020), a similar plot) and evolves on a similar time scale to
Tharakkal et al. (2022). These profiles illustrate the transition from a
linear growth regime to a nonlinear regime examined in depth by Tharakkal et
al. (2022). Additionally, these 2D simulations illustrate the cosmic ray
method implemented in Athena++ by Jiang & Oh (2018) is useful in studying the
Parker instability, since it is in reasonable agreement with other numerical
simulations.

##  Appendix B Boundary Conditions and Dimensionality

In the direction of stratification $(\hat{z})$, we use vacuum (also known as
diode) boundary conditions. These conditions make it impossible for inflow to
occur, since the boundary cells are set to the density and pressure floor of
the numerical simulation. We also avoid setting up a steep discontinuity at
this boundary by extending the simulation several scale heights in the
vertical direction. Because of this extension, the cells near the boundary are
already almost at the density and pressure floors when the simulation starts.
Then, any dynamical activity beyond that initial equilibrium profile will
propagate out of the simulation.

Figure 11: Plots matching Figure 3, but comparing simulations Bstr and Cdiff
to identical simulations with periodic boundary conditions, BstrPer and
CdiffPer. For each simulation, the solid black line is maximum height of the
flux tube at a given time. The green dotted line is the average mass in the
region $x<|250\,\mathrm{pc}|$, as a fraction of the initial average mass. The
purple solid, dashed, and dash-dotted lines are the gas, magnetic, and cosmic
ray pressure, respectively, averaged in the same region. The second $y$ axis
on the right side applies to the lines of mass and pressures, with each
quantity $Q\in\\{\langle M\rangle,\langle P_{g}\rangle,\langle
P_{B}\rangle,\langle P_{c}\rangle\\}$ normalized by either
$P_{0}=10^{-12}\,\mathrm{erg}\,\mathrm{cm}^{-3}$ or the initial value (for
mass $\langle M\rangle$).

In directions in the plane of the disk ($\hat{x}$ and $\hat{y}$) we use
outflow boundary conditions. While these allow inflow from boundary cells
after the crest of a wave passes the boundary, there is very little error if
the gas is quickly moving out of the simulation at that boundary. Our current
problem satisfies this assumption because the flow is either static, or large
(when the perturbation pushes gas down a magnetic flux tube and toward the
boundary). Ideally, we would again use vacuum boundary conditions to stop
inflow. However, we cannot use a vacuum in this direction since this would
create a significant discontinuity near the midplane of the simulation, where
there are gas densities and pressures above the floor values. Therefore,
outflow boundary conditions are the best ones for treating a single, spatially
and temporally isolated, injection.

With the goal of minimizing numerical errors, simulations often exploit
periodic boundary conditions. In the current work, we wanted to avoid tying
down the flux tubes at the boundaries of the simulation in this study, because
that can have a stabilizing effect (Zweibel & Bruhwiler, 1992). However, to
make sure our simulation results were not significantly different when using
periodic boundaries, we ran two simulations, BstrPer and CdiffPer, with
periodic boundary conditions but the same initial conditions as Bstr and
Cdiff, respectively. The flux tube dynamics of those simulations are shown in
Figure 11. Simulations Cdiff and CdiffPer show very similar results. For
simulations Bstr and BstrPer, the difference is larger: the periodic boundary
conditions tie down the magnetic flux tube at the boundaries, making it more
difficult for the tube to rise after $100\,\mathrm{Myr}$. Overall, the trends
are the same, but the timing is delayed when using periodic boundaries.

We also examined the effect of dimensionality on our simulations. Initial 2D
simulations differed from our results in 3D simulations significantly, with
tubes taking much longer to rise in 2D simulations. This delay is caused by
flux tubes above our injection being unable to move out of the way of the
rising flux tube where the injection took place. This effect is also apparent
in Figure 12, where we examine different resolutions and numerical sizes in
the third dimension, $\hat{y}$. For $8$ cells in the $\hat{y}$ direction and
resolution $\Delta y=62.5\,\mathrm{pc}$, the tube slowly rises before bursting
upward. This trend is similar to our 2D simulations, whose initial rise takes
longer. In these simulations, the rise only happens quickly because we used a
large injection energy $E_{\mathrm{inj}}\approx 10^{51}\,\mathrm{erg}$ similar
to simulation Ublast. For $16$ cells and a resolution $\Delta
y=31.25\,\mathrm{pc}$, we see a smooth well-behaved flux tube rise. By
increasing the width of the box to $32$ and $64$ numerical cells with the same
resolution as the $16$ cell run, we get clear convergence of results. From
these simulations, we chose a width of $\pm 0.5\,\mathrm{kpc}$ and resolution
of $\Delta y=20\,\mathrm{pc}$, resulting in $50$ numerical cells.

Figure 12: Height of flux tube in simulations with large injection energy
$E_{\mathrm{inj}}\approx 10^{51}\,\mathrm{erg}$, similar to simulation Ublast.
These simulations were run to examine the impact of the third dimension,
$\hat{y}$, on our results. The red solid line shows a simulation with only $8$
numerical cells in the $\hat{y}$ direction, the blue dashed line used $16$
numerical cells, the green X markers are from a simulation with $32$ numerical
cells, and the purple plus markers are from a simulation with $64$ numerical
cells. The simulations with more room in the $\hat{y}$ direction are well
converged at a length of $\pm 0.5\,\mathrm{kpc}$. We chose to use $50$
numerical cells in our production runs, with a resolution of $\Delta
y=20\,\mathrm{pc}$.

Figure 13: Convergence of key numerical parameters. Top row, first plot:
Average cell-by-cell numerical error between $\mathrm{CFL}=0.2$ and
$\mathrm{CFL}=0.1$. Top row, second plot: Average cell-by-cell numerical error
between $\mathrm{CFL}=0.15$ and $\mathrm{CFL}=0.1$. Top row, third plot:
Average cell-by-cell numerical error between a run with a third order Runge-
Kutta time integrator and a run with a fourth order Runge-Kutta time
integrator. Middle row: Average cell-by-cell numerical error between runs with
varying $V_{m}$. Each has error measured with respect to a $V_{m}=c$
simulation. Bottom row: Plot of vertical momentum $p_{z}$ in simulations with
different resolutions $dx=dz$.

##  Appendix C Simulation Convergence

To mitigate the significance of purely numerical parameters we performed
several convergence tests. We used 2D simulations to narrow our choice of time
integrator, Courant-Friedrich-Lewy (CFL) number, resolution, and the reduced
speed of light constant $V_{m}$ (see Equation 20). Using 2D simulations
allowed us to save computational resources while exploring the convergence of
these simulations. We illustrated in Section B that 2D simulations would be
less accurate because they limit the movement of magnetic flux tubes, but most
of the fast, dynamic flow is still in the $xz$-plane. For numerical
parameters, the motion in that plane is where we need be concerned. All the
convergence tests here used a weak injection of $10^{50}\,\mathrm{erg}$ with
other parameters equal to those of simulation U. With the weaker injection and
2D restricted motion, the flux tube begins to rise buoyantly by
$200\,\mathrm{Myr}$.

The time integrator choice of 3rd order Runge-Kutta and CFL number ($0.2$)
were well converged, as the average cell-by-cell error is under $10\%$ in the
top row of Figure 13. The error associated with the modified speed of light is
larger: in the middle row, the cell-by-cell error for various values of
$V_{m}$ are shown, each compared to a simulation with $V_{m}=c$. From these,
we pick $V_{m}=0.3c$ because there is not a huge increase in accuracy by going
to higher values of $V_{m}$. Increasing $V_{m}$ beyond $0.3c$ would decrease
the timestep below $0.1\,\mathrm{kyr}$, increasing computational resource
requirements without an equivalent increase in accuracy. Finally, in the
bottom row of Figure 13, we show the vertical momentum $p_{z}$ for simulations
with varying resolution. Instead of comparing these simulations via
interpolation of the high resolution simulations, we focus on how similar the
dynamics are between each simulation. Each simulation had an aspect ration of
$1$, with resolution $\Delta x=\Delta z$. Any higher resolution than
$15\,\mathrm{pc}$ appears to reproduce a similar vertical momentum on the flux
tube. Also, the structures are similar between resolutions
$15.62\,\mathrm{pc}$, $7.81\,\mathrm{pc}$, and $3.91\,\mathrm{pc}$. From
these, we decided to run our perfomance simulations with a resolution $\Delta
x=\Delta z=10\,\mathrm{pc}$.