# Suppressed Cosmic Ray Energy Densities in Molecular Clouds From Streaming
Instability-Regulated Transport

Margot Fitz Axen Department of Astronomy, The University of Texas at Austin,
2515 Speedway, Stop C1400, Austin, Texas 78712-1205, USA.  Stella Offner
Department of Astronomy, The University of Texas at Austin, 2515 Speedway,
Stop C1400, Austin, Texas 78712-1205, USA.  Philip F. Hopkins TAPIR, Mailcode
350-17, California Institute of Technology, Pasadena, CA 91125, USA  Mark R.
Krumholz Research School of Astronomy and Astrophysics, Australian National
University, Canberra ACT 2600, Australia  ARC Center of Excellence for
Astronomy in Three Dimensions (ASTRO-3D), Canberra ACT 2600, Australia
Michael Y. Grudić Carnegie Observatories, 813 Santa Barbara St, Pasadena, CA
91101, USA

###### Abstract

Cosmic rays (CRs) are the primary driver of ionization in star forming
molecular clouds (MCs). Despite their potential impacts on gas dynamics and
chemistry, no simulations of star cluster formation following the creation of
individual stars have included explicit cosmic ray transport (CRT) to date. We
conduct the first numerical simulations following the collapse of a
$2000M_{\odot}$ MC and the subsequent star formation including CRT using the
STARFORGE framework implemented in the GIZMO code. We show that when CR-
transport is streaming-dominated, the CR energy in the cloud is strongly
attenuated due to energy losses from the streaming instability. Consequently,
in a Milky Way like environment the median CR ionization rate (CRIR) in the
cloud is low ($\zeta\lesssim 2\times 10^{-19}\rm s^{-1}$) during the main star
forming epoch of the calculation and the impact of CRs on the star formation
in the cloud is limited. However, in high-CR environments, the CR distribution
in the cloud is elevated ($\zeta\lesssim 6\times 10^{-18}$), and the
relatively higher CR pressure outside the cloud causes slightly earlier cloud
collapse and increases the star formation efficiency (SFE) by $50\%$ to $\sim
13\%$. The initial mass function (IMF) is similar in all cases except with
possible variations in a high-CR environment. Further studies are needed to
explain the range of ionization rates observed in MCs and explore star
formation in extreme CR environments.

##  1 Introduction

Cosmic rays (CRs) are ubiquitous and responsible for driving a variety of
astrochemical and dynamic processes. In the interstellar medium (ISM), CRs
have an energy density of $\sim$ 1 eV $\rm cm^{-3}$, which is comparable to
the energy densities of the cosmic microwave background, starlight, and
galactic magnetic field (Ferrière, 2001). In the densest regions of molecular
clouds (MCs), which are optically thick to UV radiation, CRs are the primary
driver of ionization (Dalgarno, 2006). The level of CR ionization is usually
quantified by the CR primary ionization rate (CRIR), which is the rate of CR
ionizations per hydrogen atom per unit time.

The CRIR in the Solar neighborhood inferred from the CR particle flux measured
outside the heliosphere by Voyager is $\zeta\approx 1.6\times 10^{-17}\rm
s^{-1}$ (Cummings et al., 2016; Stone et al., 2019). Since CR ionizations lead
to chemical reactions that produce molecules such as $\text{OH}^{+}$,
$\text{H}_{2}\text{O}^{+}$, and $\text{H}_{3}\text{O}^{+}$, measurements of
these molecules can also be used to constrain the CRIR indirectly, since their
abundances are sensitive to it (Indriolo et al., 2015). Observations of these
molecules along different sightlines predict a range of values for the CRIR.
In the Milky Way, the CRIR varies between $\approx 10^{-17}-10^{-15}\rm
s^{-1}$, with values $\approx$ 10 times higher in the Galactic center than in
the Galactic disk (Indriolo et al., 2015; Neufeld & Wolfire, 2017). In higher-
density regions such as high-mass star-forming regions (Sabatini et al., 2020)
and low-mass dense cores (Caselli et al., 1998; Bialy et al., 2022), the CRIR
is expected to be lower due to attenuation by the dense gas. In contrast,
recent observations near supernova remnants (Ceccarelli et al., 2011) and
towards the more active star-forming regions in the center of our Galaxy
(Rivilla et al., 2022) and NGC 253 (Behrens et al., 2022; Holdship et al.,
2022) suggest orders of magnitude higher CRIRs than would be predicted from
Voyager data (see review by Padovani 2023).

Some of the observed variation in the CRIR can be explained by the presence of
non-SN CR acceleration sources (Ceccarelli et al., 2011). While supernova
remnants have long been thought to be the primary source of galactic CRs
(Blasi, 2013), converting about 10 % of their energy into CR acceleration
(Vink, 2012), evidence indicates that stellar winds from massive stars also
accelerate CRs (Parizot et al., 2004; Bykov et al., 2020). A number of star
clusters exhibit $\gamma$-ray emission, indicative of the presence of locally
accelerated high-energy CRs (Ackermann et al., 2011; Yang & Aharonian, 2017;
Yang et al., 2018; Aharonian et al., 2019; Sun et al., 2020; Liu et al., 2022;
Pandey et al., 2024). Some estimates suggest that stellar winds convert up to
$\sim$ 10 % of their energy into CR acceleration (Aharonian et al., 2019),
while others claim it may be as high as 40 % (Pandey et al., 2024). Most
recently, ionization maps of star-forming regions hint that protostellar
sources may also accelerate CRs (Ceccarelli et al., 2014; Cabedo et al., 2023;
Pineda et al., 2024), likely accelerated by the accretion or protostellar jet
shocks (Padovani et al., 2016; Gaches & Offner, 2018).

Observations of CRIR variation have motivated a variety of theoretical studies
on the impact of CR physics on galactic processes. To date, a growing number
of hydrodynamic simulations have included explicit cosmic ray transport (CRT)
to study the effect of CRs on galaxy formation and evolution. These
calculations suggest that CR pressure may dominate over thermal pressure in
the haloes of low redshift Milky Way mass galaxies, driving galactic outflows
and suppressing star formation (Jubelgas et al., 2008; Uhlig et al., 2012;
Booth et al., 2013; Hanasz et al., 2013; Salem et al., 2014; Pakmor et al.,
2016; Simpson et al., 2016; Wiener et al., 2017; Butsky & Quinn, 2018; Jacob
et al., 2018; Chan et al., 2019; Dashyan & Dubois, 2020; Hopkins et al., 2020;
Ji et al., 2020; Farcy et al., 2022; Girichidis et al., 2022, 2023; Thomas et
al., 2023). Stratified box calculations of patches of a Milky Way-like disk
that include CRT indicate that including plasma-dependant transport parameters
changes the impact of the CRs on the gas structure (Girichidis et al., 2016;
Simpson et al., 2016; Farber et al., 2018; Holguin et al., 2019). However, all
of these effects strongly depend on the prescription adopted for CRT; for
example, Commerçon et al. (2019) find that when CRT is treated as purely
diffusive, for diffusivities $D\lesssim 10^{25}\rm cm^{2}/s$ CRs are trapped
in the gas by pressure gradients that modify the gas flow. This effect is
particularly relevant near CR acceleration sources where the CR diffusivity is
suppressed due to turbulence, causing the gas to fragment rather than group
into star-forming clumps (Semenov et al., 2021). By contrast, higher
diffusivities produce much smaller dynamical effects, because CRs are less
effectively trapped and thus large CR pressures do not build up.

Despite evidence from observations and simulations that CR properties vary
with environment and change gas properties, most simulations and analytical
models of MCs to date have neglected the impact of CRs altogether (e.g.,
Padoan & Nordlund, 2011; Hennebelle & Chabrier, 2008; Hopkins, 2012; Li et
al., 2015; Lee et al., 2019). Most simulations have assumed a spatially and
temporally uniform CRIR of $\approx 10^{-17}-10^{-16}\rm s^{-1}$ to calculate
the heating and cooling processes (Glover & Clark, 2012; Bate & Keto, 2015;
Gatto et al., 2017; Smith et al., 2020; Grudić et al., 2021; Guszejnov et al.,
2021, 2022) or have neglected CR ionization altogether. However, in addition
to neglecting the impact of CRIR variation on the gas temperature, these
simulations do not take into account the possible dynamical effect of CRs on
the gas. For example, the level of ionization determines how much the gas
couples to the magnetic field and resists collapse (Fiedler & Mouschovias,
1993; Padovani et al., 2014). Simulations of protostellar collapse show that
the CR ionization rate also determines whether a protostellar disk forms and
regulates its size (Wurster et al., 2018; Kuffmeier et al., 2020). CR energy
deposition may also modify the gas turbulence if the CR pressure is high
enough (Commerçon et al., 2019). Although some studies have investigated the
effects of including CRT in MCs via post-processing, these also neglect the
dynamical effects between the CRs and the gas (Fatuzzo & Adams, 2014; Fitz
Axen et al., 2021).

No simulations to date have both included CRT and resolved the formation of
individual stars. In this study, we carry out large-scale simulations of MC
collapse that include CRT in order to follow the propagation of CRs from the
galactic CR background. We use the STARFORGE (STAR FORmation in Gaseous
Environments) framework, which is built on the GIZMO simulation code (Hopkins,
2015) and includes all relevant feedback processes, including protostellar
jets, radiative heating, stellar winds, and supernovae. The fiducial STARFORGE
MC simulations resolve mass scales down to $\approx 0.1M_{\odot}$ and so are
able to self-consistently predict the stellar initial mass function (IMF)
(Grudić et al., 2021, 2022; Guszejnov et al., 2021, 2022). In this work we
carry out STARFORGE simulations including the GIZMO single-bin CRT treatment,
which has been used extensively in the Feedback in Realistic Environments
(FIRE) suite of cosmological simulations (Chan et al., 2019; Hopkins et al.,
2020, 2021b, 2021d, 2021c, 2021a). This implementation applied to galaxy-scale
calculations is able to reproduce observables such as the CR intensity
spectrum observed by Voyager (Hopkins et al., 2021a).

In this work we explore the impact of CRs on ionization, gas temperatures,
star formation histories, and the IMF of star-forming MCs. Section 2 describes
the physics and numerical methods used in the STARFORGE simulations, the GIZMO
CRT modules, and the simulation initial conditions. Section 3 presents the
results of simulations with different recipes for CRT. Section 4 discusses how
CRT properties shape our results and compares to prior work. Finally, Section
5 summarizes our results.

##  2 Methods

###  2.1 The STARFORGE Simulations

####  2.1.1 Physics

We simulate star-forming MCs using the STARFORGE framework, which uses the
GIZMO simulation code (Hopkins, 2015). A full description of the STARFORGE
methods is included in Grudić et al. (2021); therefore, we only summarize the
key points here.

We use the Lagrangian constrained-gradient meshless finite-mass method for
magnetohydrodynamics (MHD) (Hopkins & Raives, 2016; Hopkins, 2016) and assume
the ideal MHD limit. The simulations utilize the updated FIRE-3 radiation and
chemistry modules (Hopkins et al., 2023). They co-evolve the gas, dust, and
radiation temperature from the stellar luminosity and an external radiation
field.

Sink particles represent individual stars. Once they form, they follow the
subgrid protostellar evolution model in Offner et al. (2009). The STARFORGE
framework includes stellar feedback from protostellar jets, radiation, stellar
winds, and supernovae (SNe). Feedback from SNe, radiation, and most stellar
winds is injected by distributing mass, momentum, and energy into surrounding
gas cells using the weighting scheme described in Hopkins et al. (2018).
Feedback from protostellar jets and stellar winds, where the free-expansion
radius is smaller than the size of a wind cell, is injected using the cell
spawning technique described in Grudić et al. (2021).

Protostellar jets are modeled using the jet feedback model of Cunningham et
al. (2011) in which sink particles launch a fraction of their accreted
material $f_{w}$ along their rotational axis with a velocity that is a
fraction $f_{K}$ of the Keplerian velocity defined at the protostellar radius.
We use $f_{w}=f_{K}=0.3$, which Rosen & Krumholz (2020) found reproduces
measurements of momentum injection as a function of protostellar luminosity
(Maud et al., 2015). We do not include CR acceleration from protostellar jets.

Massive main-sequence stars inject a significant amount of energy and momentum
into their environment through radiative feedback and stellar winds. Winds are
launched for stars more massive than $2M_{\odot}$. The stars lose mass
isotropically at a mass-loss rate $\dot{M}_{w}$ moving at velocity $v_{w}$.
Wind mass-loss rates and velocities are calculated from Smith (2014, see
Equation 1 of ()) and Lamers et al. (1995, see Equation 45 of ()). The
luminosity of the winds is $\dot{E}_{w}=0.5\dot{M}_{w}v_{w}^{2}$. Typical wind
luminosities are $\sim 0.1-10L_{\odot}$ (or $\sim 10^{40}-10^{42}$ erg/yr).
Stellar winds with velocities $v_{w}>1000$ km/s inject $10\%$ of their energy
into CRs, consistent with estimates made from $\gamma$-ray measurements in
massive star clusters (Aharonian et al., 2019; Pandey et al., 2024). The CR
energy is injected by including an energy source term $S_{e}=0.1\dot{E}_{w}$
and a momentum source term
$S_{F_{e}}=S_{e}c(\hat{v_{w}}\cdot\hat{\mathbb{B}})$ parallel to the local
magnetic field $\hat{\mathbb{B}}$ in the CR transport equations (Section
2.2.1). To satisfy conservation of energy the wind velocity is reduced to
$v_{w}\rightarrow\sqrt{0.9}v_{w}$. All of our calculations form sufficiently
massive stars such that CRs are injected into the domain.

The STARFORGE framework also includes SNe, which are spawned at the end of the
lifetime of main sequence stars $>8M_{\odot}$. However, all our simulations
end before SNe occur.

####  2.1.2 Initial Cloud Properties

Our initial cloud properties are similar to those used in previous STARFORGE
simulations (Guszejnov et al., 2021, 2022; Grudić et al., 2021, 2022). We
start with a cloud of gas with mass $M_{\rm cloud}=2000M_{\odot}$ and radius
$R_{\rm cloud}=3$ pc at a uniform density of $\rho_{\rm cloud}=M_{\rm
cloud}/(4\pi R_{\rm cloud}^{3}/3)$ and a uniform temperature $T_{0}=10$ K. The
cloud is placed in a diffuse medium with a density of $\rho_{\rm cloud}/1000$,
within a periodic box of size $10R_{\rm cloud}$. The initial velocity field is
initialized so that the velocity power spectrum varies with wavenumber as
$E_{k}\propto k^{-2}$. Velocities are scaled to the value set by the turbulent
virial parameter $\alpha_{\rm turb}=2$, which measures the relative importance
of turbulence and gravity (Bertoldi & McKee, 1992; Federrath & Klessen, 2012):

$\alpha_{\rm turb}=\frac{5\sigma_{\rm 1D,cloud}^{2}R_{\rm cloud}}{GM_{\rm
cloud}},$ (1)

where $\sigma_{\rm 1D,cloud}$ is the line-of-sight velocity dispersion. The
clouds have a uniform magnetic field $B_{z}$ set by the mass-to-flux ratio
$\mu=1.3$:

$\mu=c_{1}\sqrt{\frac{E_{\rm mag}}{E_{\rm grav}}},$ (2)

where $E_{\rm mag}$ and $E_{\rm grav}$ are the magnetic and gravitational
energy, respectively, and the normalization constant $c_{1}\approx 0.4$
(Mouschovias & Spitzer, 1976). We use a mass resolution of $10^{-3}M_{\odot}$
for gas cells, so each cloud initially consists of $2\times 10^{6}$ gas cells.
The minimum resolved Jeans length is $\Delta x_{J}=36$ au.

Table 1 restates the cloud properties in our simulations. Throughout this
study, we present results in terms of the cloud freefall time, defined as

$t_{\rm ff}=\sqrt{\frac{3\pi}{32G\rho_{\rm cloud}}},$ (3)

which is $t_{\rm ff}\approx 2$ Myr. We run the simulations until the clouds
are completely dispersed by feedback, which occurs at $\sim 8$ Myr (or
$4t_{\rm ff}$).

$M_{\rm cloud}(M_{\odot})$ |  $R_{\rm cloud}$ (pc) | $\alpha_{\rm turb}$ | $\mu$ |  $T_{0}$ (K) | $M_{\rm cloud}/\Delta m$ |  $\Delta x_{J}$ (au)
---|---|---|---|---|---|---
$2\times 10^{3}$ | 3 | 2 | 1.3 | 10 | $2\times 10^{6}$ | 36
Table 1: Initial conditions of the simulated clouds. Columns are the initial
cloud mass, size, virial parameter, mass to magnetic flux ratio, temperature,
mass resolution and minimum resolved Jeans length (Guszejnov et al., 2022).

###  2.2 Cosmic Ray Methodology

####  2.2.1 Cosmic Ray Transport

Our study extends the STARFORGE simulation suite by including the GIZMO module
for “single energy bin” CRT. This evolves the CR total energy and flux
directly, integrating over the CR spectrum, so is effectively representing the
$\sim 0.5-10\,$GeV protons, which contain most of the CR energy and pressure.
This CR treatment has been used extensively in the FIRE suite of cosmological
simulations (Chan et al., 2019; Hopkins et al., 2020, 2021b, 2021d, 2021c,
2021a, 2022a, 2022b). The GIZMO implementation treats CRs as a relativistic
fluid ($\gamma_{\rm cr}=4/3$) obeying an appropriate set of fluid-like
equations derived from the collisionless Vlasov equation. We use the treatment
of CRs first formulated in Hopkins et al. (2022a), which integrates a set of
moment equations using a reduced-speed-of-light (RSOL) approximation. This
method allows for anisotropic CR transport with streaming, advection,
diffusion, and energy losses. We refer readers to Hopkins et al. (2022a) for a
detailed derivation of the CRT equations and provide only the necessary
details here.

We evolve the total CR energy $e_{\rm cr}$ and the CR energy flux $F_{e,\rm
cr}$; note that $F_{e,\rm cr}$ is a scalar quantity describing the flux
parallel to the local magnetic field; the cross-field flux is assumed to be
zero. The expressions solved for $e_{\rm cr}$ and $F_{e,\rm cr}$ are given by

$D_{t}e_{\rm cr}+\nabla\cdot(F_{e,\rm cr}\hat{\mathbb{B}})=\tilde{S}_{\rm
sc}+S_{e}-\Gamma_{\rm cr}-\Gamma_{\rm st}-\mathds{P}_{\rm
cr}:\nabla\mathbb{u},$ (4) $D_{t}F_{e,\rm cr}+c^{2}\,\hat{\bf
B}\cdot(\nabla\cdot\mathds{P}_{\rm cr})=-\bar{\nu}_{s}\left[F_{e,\rm
cr}-3\chi_{e}v_{\rm st}(e_{\rm cr}+P_{\rm cr})\right]+S_{F_{e}},$ (5)

where $\mathbb{u}$ is the gas velocity, $\hat{\bf B}$ is the local magnetic
field direction, $\mathds{P}_{\rm cr}=3P_{\rm cr}\mathds{D_{\rm cr}}$ is the
CR pressure tensor with a scalar pressure $P_{\rm cr}\approx e_{\rm cr}/3$,
$\mathds{D_{\rm cr}}$ is the Eddington tensor,
$\bar{\nu}_{s}=\bar{\nu}_{s}^{+}+\bar{\nu}_{s}^{-}$ is the total scattering
rate from forward and backward scattering rates in the Alfvèn frame (related
to the CR diffusion coefficient $D_{\parallel}=v_{\rm cr}^{2}/3\bar{\nu}_{s}$
discussed in Section 3.1.3), $v_{\rm cr}\approx c$ is the velocity of an
individual CR particle, and $\chi_{e}$ is a closure variable (discussed
further below). The streaming velocity $v_{\rm st}$ is approximated as the ion
Alfvén speed $v_{\rm Ai}=B/\sqrt{4\pi\rho_{\rm ion}}$, where $\rho_{\rm
ion}=\chi\rho_{\rm gas}=m_{\rm ion}nf_{\rm ion}$, $\chi=f_{\rm ion}m_{\rm
ion}/m_{p}$ is the mass-weighted ionization fraction, $\rho_{\rm gas}$ is the
total gas density, and $m_{\rm ion}$ and $f_{\rm ion}$ are the mean ion mass
and ionization fraction by number respectively. 111Formally, the streaming
speed is $[(\nu_{+}-\nu_{-})/(\nu_{+}+\nu_{-})]\,v_{\rm Ai}$ where the
prefactor represents the fractional asymmetry in the forward/backward
scattering rates. But for $\sim$ GeV in the limit of interest where the
streaming speed is large (the case we explore here), this is almost always
expected to be $\sim 1$.,222The mean ion mass $m_{\rm ion}$ is approximated as
$m_{\rm ion}=m_{p}$, which is equivalent to setting $\chi=f_{\rm ion}$. This
is not accurate in MCs where the dominant charge carrier is $\rm HCO^{+}$
($m_{\rm ion}=29m_{p}$) (Krumholz et al., 2020); however, since the CRT is
calculated assuming $m_{\rm ion}=m_{p}$, for consistency we make this
assumption throughout this study for analysis. We discuss the implications of
this in Section 4.3.4. The $\mathds{P}_{\rm cr}:\nabla\mathbb{u}$ term in
Equation 4 represents the adiabatic compression of CRs with the gas. The term
$\tilde{S}_{\rm sc}=(\bar{\nu}_{s}/c^{2})3\chi_{e}v_{\rm Ai}^{2}(e_{\rm
cr}+P_{\rm cr})$ represents diffusive reacceleration. The $S_{e}$ and
$S_{F_{e}}$ terms represent sources of CR energy and momentum respectively,
which in our simulations are only CR injection from stellar winds. The
$\Gamma_{\rm cr}$ and $\Gamma_{\rm st}=-(\bar{\nu}_{s}/c^{2})v_{\rm
st}F_{e,\rm cr}$ terms represent collisional and streaming instability energy
losses respectively, which are discussed further in Section 2.2.3. Finally,
the $-\bar{\nu}_{s}F_{e,\rm cr}$ and $\bar{\nu}_{s}3\chi_{e}v_{\rm st}(e_{\rm
cr}+P_{\rm cr})$ terms in Equation 5 represent scattering and streaming,
respectively, where the latter arises if the scattering is asymmetric.

The Eddington tensor $\mathds{D_{\rm
cr}}=\chi_{e}\mathds{I}+(1-3\chi_{e})\hat{\mathbb{B}}\hat{\mathbb{B}}$ is
defined in terms of the variable $\chi_{e}=(1-\langle\mu_{f}^{2}\rangle)/2$,
where a closure relation is used for $\langle\mu_{f}^{2}\rangle$ to provide a
complete system of equations. Hopkins et al. (2022a) assumes a closure
relation analogous to that adopted in radiation hydrodynamics (Levermore,
1984), which is a function that smoothly interpolates between the diffusive
and streaming limits:

$\langle\mu_{f}^{2}\rangle\approx\frac{3+4\langle\mu_{f}^{1}\rangle^{2}}{5+2(4-3\langle\mu_{f}^{1}\rangle^{2})^{1/2}},$
(6)

where $\langle\mu_{f}^{1}\rangle=F_{e,\rm cr}/v_{\rm cr}e_{\rm cr}$. In the
diffusive limit, $\langle\mu_{f}^{1}\rangle=0$ and
$\langle\mu_{f}^{2}\rangle=1/3$, so $\chi_{e}=1/3$. In the free streaming
limit, $\langle\mu_{f}^{1}\rangle=1$ and $\langle\mu_{f}^{2}\rangle=1$, so
$\chi_{e}=0$.

####  2.2.2 Coupled Gas and Cosmic-Ray Equations of Motion

Since CRs exchange energy and momentum with the gas via scattering and Lorentz
forces (Zweibel, 2013, 2017; Hopkins et al., 2022a), the default GIZMO
momentum and energy equations must be modified to account for the influence of
CRs. Spatial differences in the CR pressure tensor $\nabla\cdot\mathds{P}_{\rm
cr}$ produce a net CR current. The current parallel to $\hat{\mathbb{B}}$ is
resisted only by exchanging momentum with the gas; therefore, the force
balances the scattering + streaming term in Equation 5. The perpendicular
current is redirected by Lorentz forces that exert an equal and opposite force
on the gas, which is $\approx-\nabla_{\bot}\cdot\mathds{P}_{\rm cr}$ (Zweibel,
2017; Hopkins et al., 2022a)where for convenience we have introduced the
shorthand notation
$\nabla_{\parallel}\cdot\mathds{Q}\equiv\hat{\mathbb{B}}(\nabla\cdot\mathds{Q})$
and
$\nabla_{\perp}\cdot\mathds{Q}\equiv\nabla\cdot\mathds{Q}-\nabla_{\parallel}\cdot\mathds{Q}$
for any quantity $\mathds{Q}$. Using Equation 5, the final modified momentum
equation can be written as

$\displaystyle D_{t}(\rho\mathbb{u})+\ldots+\nabla\cdot\mathds{P}_{\rm cr}=$
$\displaystyle\hat{\mathbb{B}}(\nabla\cdot\mathds{P}_{\rm
cr})+\left(\frac{\hat{\mathbb{B}}}{c^{2}}\right)\bar{\nu}_{s}\left[F_{e,\rm
cr}-3\chi_{e}v_{\rm st}(e_{\rm cr}+P_{\rm cr})\right],$ (7)

where $\ldots$ refers to all non-CR terms. This has a pressure gradient term
and a source term that vanishes when the flux is in steady state.

The modifications to the energy equation balance the CR energy gains and
losses. These include CR energy losses due to work done on the gas outlined in
Equation 7 along with the four gain and loss terms shown in Equation 4 (energy
gained from diffusive reacceleration, CR energy injection, energy lost due to
streaming, and collisional and non-collisional losses). Collecting all these
terms, the final modified energy equation is

$D_{t}(e_{\rm tot})+\ldots=\mathbb{u}\cdot D_{t}(\rho\mathbb{u})|_{\rm
cr}+\Gamma_{\rm st}+S_{\rm gas}-\tilde{S}_{\rm sc}-S_{e},$ (8)

where $e_{\rm tot}=e_{\rm kin}+e_{\rm mag}+e_{\rm rad}+e_{\rm int}$ is the sum
of the gas kinetic, magnetic, radiation, and internal energies,
$D_{t}(\rho\mathbb{u})|_{\rm cr}$ is the CR term on the right side of Equation
7, and $S_{\rm gas}$ is the gas heating term due to CR collisional energy
losses (described in Section 2.2.3).

####  2.2.3 Energy Loss Processes

While they propagate, CRs lose energy due to interactions with non-
relativistic particles. Here we adopt an estimate for combined hadronic,
Coulomb, and ionization losses. The hadronic/Coulomb losses are from Völk et
al. (1996) and Ensslin et al. (1997) and are updated in Guo & Oh (2008). The
total CR energy loss rate due to these processes is $\Lambda_{\rm
cr}=\Lambda_{\rm coul}+\Lambda_{\rm had}+\Lambda_{\rm ion}$,

$\Lambda_{\rm cr}=7.51\times 10^{-16}(1+0.22\tilde{n}_{e}+0.125f_{\rm
neut})\left(\frac{n_{\rm H}}{\rm cm^{-3}}\right)\rm s^{-1},$ (9)

where $n_{\rm H}$ is the hydrogen number density, $\tilde{n}_{e}$ is the
number of free electrons per hydrogen nucleus, and $f_{\rm neut}=(1-f_{\rm
ion})$ is the neutral gas fraction. The energy loss term in Equation 4 is then
$\Gamma_{\rm cr}=\Lambda_{\rm cr}e_{\rm cr}$. Following Guo & Oh (2008), we
assume 1/6 of the hadronic losses and all non-hadronic losses are thermalized,
adding a volumetric heating component to the gas

$\Lambda_{\rm gas}=7.51\times 10^{-16}\left(0.17+0.22\tilde{n}_{e}+0.125f_{\rm
neut}\right)\left(\frac{n_{\rm H}}{\rm cm^{-3}}\right)\rm s^{-1}.$ (10)

The gas heating rate in Equation 8 is then $S_{\rm gas}=\Lambda_{\rm
gas}e_{\rm cr}$.

The CRs also lose energy due to the streaming instability. As they stream, CRs
excite Alfvèn waves caused by anisotropy in the CR distribution. This drives
the CRs to stream down the CR pressure gradient and effectively travel at the
ion-Alfvèn speed. The CR energy is transferred to the waves, which are quickly
thermalized and heat the gas (Kulsrud & Pearce, 1969; Skilling, 1975). CR
streaming losses are parameterized in Equation 4 by the term $\Gamma_{\rm
st}$. We describe the magnitude and impact of these losses in more detail in
Section 4.1.1.

####  2.2.4 Reduced Speed of Light Approximation

Explicitly integrating Equations 4 and 5 imposes a Courant-type time-step
limiter $\Delta t\leq C\Delta x/c$, which is numerically prohibitive because
$c$ is significantly faster than any other speed in the problem. Therefore,
analogously to radiation-hydrodynamics (RHD), we adopt an RSOL approximation,
which has been done in many previous CR studies (Jiang & Oh, 2018; Chan et
al., 2019; Hopkins et al., 2020). Following Hopkins et al. (2022a), we use an
RSOL approximation for Equations 4 and 5 equivalent to replacing
$c^{-1}D_{t}e_{\rm cr}\rightarrow\tilde{c}^{-1}D_{t}e_{\rm cr}$ and
$c^{-1}D_{t}(F_{e,\rm cr}/c)\rightarrow\tilde{c}^{-1}D_{t}(F_{e,\rm cr}/c)$ in
the original CR transport equation:

$\displaystyle\frac{1}{\tilde{c}}D_{t}e_{\rm
cr}+\nabla\cdot\left(\frac{F_{e,\rm cr}}{c}\hat{\mathbb{B}}\right)$
$\displaystyle=\frac{1}{c}\left(\tilde{S}_{\rm sc}+S_{e}-\Gamma_{\rm
cr}-\Gamma_{\rm st}-\mathds{P}_{\rm cr}:\nabla\mathbb{u}\right),$ (11)
$\displaystyle\frac{1}{\tilde{c}}D_{t}\left(\frac{F_{e,\rm
cr}}{c}\right)+(\nabla\cdot\mathds{P}_{\rm cr})$
$\displaystyle=\left(\frac{-\bar{\nu}_{s}}{c^{2}}\right)\left[F_{e,\rm
cr}-3\chi_{e}v_{\rm st}(e_{\rm cr}+P_{\rm
cr})\right]+\left(\frac{1}{c^{2}}\right)S_{F_{e}}.$ (12)

We use a value of $\tilde{c}=300\rm km/s$ in our simulations, which is
justified through a low-resolution convergence study presented in the
Appendix.

The STARFORGE simulations also use an RSOL approximation in the radiation flux
equation (Grudić et al., 2021). The default STARFORGE RSOL for radiation is
either $30\rm km/s$ (Grudić et al., 2021, 2022) or $90\rm km/s$ (Grudić et
al., 2023); however, in these runs we set $\tilde{c}=300\rm km/s$ for
radiation for consistency since it does not increase the computational cost.
We note that adopting a higher RSOL is primarily important for accurately
capturing the optically thick, strong-radiation field limits characteristic of
massive star formation (Skinner & Ostriker, 2013).

####  2.2.5 Cosmic Ray Diffusion Coefficient - Constant Coefficient Case

The key ingredient in our CR evolution equations that we have not yet
specified is the value of the diffusion coefficient $D_{\parallel}$; the value
of this parameter is, as discussed in Section 1, crucial to the outcome. For
our fiducial simulations we adopt a diffusion coefficient $D_{\parallel}$ that
is constant spatially and temporally.

We derive the diffusion coefficient following Sampson et al. (2022), where the
authors use an ensemble of MHD turbulence simulations to quantity how the
basic parameters describing ISM turbulence affect the transport of streaming
CRs. They make the assumption that the total scattering rate $\bar{\nu}_{s}$
is extremely large (driven by the CR streaming instability) and also
asymmetric. In these limits, $v_{\rm st}\rightarrow v_{\rm Ai}$ (Section
2.2.1) so that streaming dominates the CR transport and the microphysical
diffusion coefficient $D_{\parallel,\rm micro}=c^{2}/3\bar{\nu}_{s}$ arising
from CR pitch angle scattering is negligibly small in comparison to the
effective diffusion coefficient induced by random walk of magnetic field lines
(and related processes) due to unresolved turbulence. Under this assumption,
they derive a macroscopic diffusion coefficient $D_{\parallel,\rm FLW}$ that
quantifies this process; the model is effectively a Large Eddy Simulation
closure for CR transport.

Since $D_{\parallel,\rm FLW}$ does not arise from CR pitch angle scattering,
it and $D_{\parallel,\rm micro}$ cannot be formally added together in
Equations 4 and 5 (there is no scattering rate $\bar{\nu}_{s}$ that can be
derived from the two that can then be used to obtain a total diffusion
coefficient using $\bar{\nu}_{s}=c^{2}/3D_{\parallel}$). However, as long as
Equation 5 is in steady state ($D_{t}F_{e,\rm cr}\rightarrow 0$) all terms
appear as they should. In our simulations, Equation 5 reaches equilibrium on a
timescale $\Delta t_{\rm F}\sim(c/\tilde{c})\bar{\nu}_{s}^{-1}\sim 10$ yr
(Hopkins et al., 2021b), which is much shorter than the macroscopic timescales
we consider. Finally, since $D_{\parallel,\rm FLW}$ is derived under the
assumption that $D_{\parallel,\rm FLW}\gg D_{\parallel,\rm micro}$ (and we
show in Section 2.2.6 that this is a valid assumption), we arrive at
$D_{\parallel}\sim D_{\parallel,\rm FLW}$.

The diffusion coefficient of interest here is that measured on the minimum
size scale $l_{\mathrm{min}}$ that is resolved by our simulations, since
processes on larger scales are presumably resolved directly. Since GIZMO is a
Lagrangian method, the spatial resolution varies throughout the simulation
according to the local conditions. Therefore, we estimate $l_{\rm min}$ based
on the dynamic range of the cloud turbulent power spectrum. The minimum
resolved length scale for a turbulent sphere is given by $l_{\rm min}=R_{\rm
cloud}/k_{\rm max}$, where $k_{\rm max}=N_{\rm 1D}/N_{\rm cells,eff}$ is the
maximum wavenumber of the turbulence, $N_{\rm 1D}=(N_{\rm
gas}/(4\pi/3))^{1/3}$ is the 1D resolution of the simulation and $N_{\rm
cells,eff}$ is the minimum number of cells required to resolve a turbulent
mode. For our simulations $N_{\rm gas}=2\times 10^{6}$ and $N_{\rm
cells,eff}\approx 10$, which give $l_{\rm min}=0.38$ pc.

Sampson et al. (2022) present empirical fits for the parallel and
perpendicular diffusion coefficient as a function of plasma parameters to be
used as sub-grid recipes in simulations. We use the initial cloud magnetic
field value $B\approx 7\mu\rm G$ and average gas density $\rho_{\rm
cloud}\approx 18M_{\odot}/\rm pc^{3}$. For the ionization fraction, we adopt
$\chi=10^{-7}$, which is within the range of observed values for Milky-Way
molecular clouds (Caselli et al., 1998; Williams et al., 1998). For the
velocity dispersion, we adopt the observed velocity dispersion for turbulent
molecular clouds (McKee & Ostriker, 2007)

$\sigma_{\rm 1D}=0.9\left(\frac{l_{\rm min}}{1\rm pc}\right)^{0.56}\rm
km\,s^{-1},$ (13)

which gives $\sigma_{\rm 1D}\approx 0.53$ km/s. Finally, we assume Gaussian
diffusion (where the CR distribution function follows $\partial f/\partial
t=D_{\rm FLW}\nabla^{2}f$), which is predicted for the characteristics of our
problem. Following Sampson et al. (2022), we obtain values for the parallel
and perpendicular diffusion coefficients of $D_{\parallel,\rm FLW}\approx
8.33\times 10^{25}\rm cm^{2}/s$ and $D_{\bot,\rm FLW}\approx 4.3\times
10^{22}\rm cm^{2}/s$. Since $D_{\bot,\rm FLW}\ll D_{\parallel,\rm FLW}$, we
assume $D_{\bot,\rm FLW}=0$ and adopt $D_{\parallel,\rm FLW}$ as the constant
value of the diffusion coefficient.

We note that these derived diffusion coefficients are much lower than the
average galactic diffusion coefficient $D_{\rm obs}\sim 10^{28}\rm cm^{2}/s$
derived from CR data (Evoli et al., 2019) for two reasons. First, we resolve
CR streaming explicitly in the CR transport equations, while the observed
diffusion coefficient arises from both scattering and unresolved streaming
($D_{\rm obs}=D_{\parallel}+v_{\rm st}\ell_{\rm cr}$). Additionally, $D_{\rm
obs}$ is measured in the ionized ISM where $D_{\parallel,\rm FLW}\ll
D_{\parallel,\rm micro}$. The magnitude of $D_{\parallel,\rm FLW}$ is orders
of magnitude smaller when $f_{\rm ion}\sim 1$ compared to where $f_{\rm
ion}\sim 10^{-7}$ as in MCs (Sampson et al., 2022); in contrast when $f_{\rm
ion}\sim 1$, $D_{\parallel,\rm micro}$ arises from a different pitch-angle
scattering mechanism and is much larger. In principle, $D_{\parallel,\rm
micro}$ arises from CR pitch angle scattering off both the Alfvèn waves self-
generated by the streaming instability and the extrinsic ISM turbulence
($\bar{\nu}_{s}=\bar{\nu}_{s}^{\rm sc}+\bar{\nu}_{s}^{\rm et}$). Sampson et
al. (2022) assume $\bar{\nu}_{s}^{\rm sc}\gg\bar{\nu}_{s}^{\rm et}$, and that
$\bar{\nu}_{s}^{\rm sc}$ is so large that the corresponding diffusion
coefficient $D_{\parallel,\rm micro}\sim D_{\rm sc}$ is negligibly small; we
show in Section 2.2.6 that this is a reasonable assumption in MCs. In
contrast, in higher ionization environments CR scattering is dominated by
resonant interactions with extrinsic turbulence (Xu & Lazarian, 2017; Krumholz
et al., 2020), so $\bar{\nu}_{s}^{\rm et}\gg\bar{\nu}_{s}^{\rm sc}$ and
$D_{\parallel,\rm micro}\sim D_{\rm obs}\sim D_{\rm et}$ (Krumholz et al.,
2020). In short, $D_{\rm obs}$ cannot be compared to either the value
$D_{\parallel,\rm FLW}$ derived in this section or the value of
$D_{\parallel,\rm micro}$ derived in Section 2.2.6.

In addition to our derived value of $D_{\parallel,\rm FLW}$, we explore the
impact of varying $D_{\parallel,\rm FLW}$ by an order of magnitude. We also
run simulations with a non-uniform CR diffusion coefficient, which we describe
in Section 2.2.6. Our full suite of simulations, including the diffusion
coefficient used for each, is described in Table 2.

####  2.2.6 Cosmic Ray Diffusion Coefficient - Variable Coefficient Case

In addition to running simulations with a constant CR diffusion coefficient,
we also perform one simulation with a CR diffusion coefficient that varies
spatially and temporally. Unlike the diffusion coefficient described in
Section 3.1.3, which assumes CRT is dominated by streaming and the CR
scattering rate $\bar{\nu}_{s}$ is large ($D_{\parallel,\rm micro}\ll
D_{\parallel,\rm FLW}$), this simulation calculates the microphysical
diffusion coefficient $D_{\parallel,\rm micro}$ self-consistently based on the
local gas and CR properties and neglects field line wandering
($D_{\parallel,\rm FLW}\rightarrow 0$). This has the advantage that it
predicts whether CRs should be in the streaming dominated limit rather than
beginning with that postulate. It should also correctly model ‘ballistic’ CR
transport where $\bar{\nu}_{s}$ is small enough that CR scattering is entirely
negligible for the entire CR trajectory and the CRs free stream down field
lines through the cloud.

We use a model for the microphysical diffusion coefficient included in GIZMO
which encapsulates CR pitch angle scattering off both the Alfvèn waves self-
generated by the streaming instability and scattering off extrinsic ISM
turbulence ($\bar{\nu}_{s}=\bar{\nu}_{s}^{\rm sc}+\bar{\nu}_{s}^{\rm et}$)
(Hopkins et al., 2022b). For the extrinsic turbulence scattering component
$\bar{\nu}_{s}^{\rm et}$, we use the default model from Hopkins et al.
(2022b). We find in our simulations that in low-ionization environments such
as starburst galaxies and MCs, efficient ion-neutral damping prevents
extrinsic ISM turbulence from cascading down to the scales of CR gyroradii
(Krumholz et al., 2020), and so this contribution is largely negligible.

Instead, CRT is dominated by streaming along field lines at the streaming
speed $v_{\rm st}$, which is the ion-Alfv̀en speed and is modeled directly in
the streaming terms in Equations 4 and 5. In steady state, the wave damping
rate $\Gamma_{\rm damp}$ induced by ion-neutral damping and other various
damping mechanisms must be matched by the growth rate of the streaming
instability $\Gamma_{\rm SI}$ (Hopkins et al., 2022b). The bulk CR propagation
speed is $v_{D}=v_{\rm st}+v_{\rm diff}=v_{\rm Ai}+v_{\rm diff}$, where
$v_{\rm diff}$ is the diffusive velocity component that arises from scattering
off the CR-induced Alfvèn waves and $D_{\rm sc}\propto v_{\rm diff}$ is the
corresponding diffusion coefficient. Since $\Gamma_{\rm SI}\propto v_{D}$, an
increase in $\Gamma_{\rm damp}$ allows the CRs to propagate faster ($v_{\rm
diff}$ and thus $D_{\rm sc}$ increase) until once again $\Gamma_{\rm
damp}=\Gamma_{\rm SI}$. In other words, increased wave damping induces a
finite scattering rate off the Alfvèn waves $\bar{\nu}_{s}^{\rm sc}$ with a
diffusion coefficient $D_{\rm sc}=c^{2}/3\bar{\nu}_{s}^{\rm sc}$.

In general, the damping rate $\Gamma_{\rm damp}$ includes multiple damping
mechanisms including ion-neutral, dust, linear Landau and non-linear Landau
damping (with rates given in Hopkins et al. 2022b). In practice, here ion-
neutral damping dominates and other sources of damping are negligible. The
efficiency of ion-neutral damping depends on the local gas and CR properties,
and can create orders-of-magnitude variations in the CR diffusivity. The ion-
neutral damping rate is $\Gamma_{\rm IN}=\nu_{\rm in}/2$, where $\nu_{\rm in}$
is the collision frequency between ions and neutrals; in a hydrogen-helium
plasma this is approximated as

$\Gamma_{\rm IN}\approx 10^{-9}f_{\rm neutral}\sqrt{\frac{T_{\rm gas}}{1000\rm
K}}\left(\frac{\rho_{\rm gas}}{10^{-24}\rm
g/cm^{3}}\right)\left(\frac{m_{p}}{m_{\rm ion}}\right)\rm s^{-1},$ (14)

where $m_{p}$ is the proton mass, $m_{\rm ion}$ is the mean ion mass, $f_{\rm
neutral}=(1-f_{\rm ion})$ is the neutral gas fraction by number, $T_{\rm gas}$
is the gas temperature, and $\rho_{\rm gas}\approx\mu m_{p}n_{\rm gas}$ is the
total gas density with mean molecular weight $\mu$ (Kulsrud & Pearce, 1969;
Bustard & Zweibel, 2021; Hopkins et al., 2022b). The growth rate of the
streaming instability is (Kulsrud, 2005)

$\Gamma_{\rm SI}\approx\alpha\Omega_{\rm ion}\left(\frac{n_{\rm CR}}{n_{\rm
ion}}\right)\left(\frac{v_{D}}{v_{\rm Ai}}-1\right),$ (15)

where $n_{\rm ion}\approx f_{\rm ion}\rho_{\rm gas}/\mu m_{p}$ is the ion
number density, $n_{\rm CR}\approx\epsilon_{\rm cr}/E_{\rm cr}$ is the CR
number density for the CR energy density $\epsilon_{\rm cr}$ and
characteristic energy of a CR particle $E_{\rm cr}$, $\Omega_{\rm
ion}=eB/m_{\rm ion}c$ is the non-relativistic ion frequency, and
$\alpha\approx 1$ is a constant that depends on the details of the numerical
integrals over non-gyroresonant interactions. Setting $\Gamma_{\rm
SI}=\Gamma_{\rm IN}$ and using $f_{\rm ion}\approx n_{\rm ion}/n_{\rm gas}$
and $E_{\rm cr}\sim 1\rm GeV$,

$\displaystyle v_{\rm diff}\approx 300v_{\rm Ai}f_{\rm ion}f_{\rm
neut}\sqrt{\frac{T_{\rm gas}}{1000\rm K}}\left(\frac{\rho_{\rm
gas}}{10^{-24}\rm g/cm^{3}}\right)^{2}$ $\displaystyle\times\left(\frac{B}{\mu
G}\right)^{-1}\left(\frac{\epsilon_{\rm cr}}{\rm eVcm^{-3}}\right)^{-1}.$ (16)

The diffusion coefficient $D_{\rm sc}$ is related to the diffusive component
of the drift velocity as $D_{\rm sc}=v_{\rm diff}\ell_{\rm cr}$, where we
define the CR gradient scale length $\ell_{\rm cr}=e_{\rm cr}/\nabla e_{\rm
cr}$. The CR energy gradient $\nabla e_{\rm cr}$ is computed in GIZMO using a
slope-limited, second order least-squares gradient estimator (Hopkins, 2015).
Using typical cloud properties in MCs ($f_{\rm ion}\sim 10^{-7}$, $T\sim 10$
K, $\mu\approx 2.33$, $B\sim 10\mu G$, $n_{\rm gas}\sim 100\rm cm^{-3}$),
$v_{\rm diff}\sim 0.05v_{\rm Ai}(\epsilon_{\rm cr}/\rm eVcm^{-3})^{-1}$, so
$D_{\rm sc}\sim 5\times 10^{26}\,(v_{\rm A,ideal}/10\,\rm km\,s^{-1})$$(f_{\rm
ion}/10^{-7})^{-1/2}\,(\epsilon_{\rm cr}/\rm eVcm^{-3})^{-1}$ $(\ell_{\rm
cr}/\rm pc)\,cm^{2}/s$. The streaming component of the effective diffusion
coefficient at larger scales (which we resolve directly in our simulations) is
$v_{\rm Ai}\ell_{\rm cr}\sim 10^{28}\,(v_{\rm A,ideal}/10\,\rm
km\,s^{-1})$$(f_{\rm ion}/10^{-7})^{-1/2}(\ell_{\rm cr}/\rm pc)\,cm^{2}/s$.
Thus, although the streaming instability necessitates CR scattering, the CR
transport is still dominated by streaming down field lines unless
$\epsilon_{\rm cr}\ll 1\rm eV/cm^{3}$.

For our simulations, we set a lower limit of $8.33\times 10^{24}\rm cm^{2}/s$
and an upper limit of $10^{34}\rm cm^{2}/s$ for $D_{\parallel,\rm micro}$. The
upper limit is rarely used and makes little difference to our results, but the
lower limit prevents regions outside of the cloud from using very low values
of $D_{\parallel,\rm micro}$ that would cause the CRs to quickly lose all
their energy before propagating into the cloud.Equation 14 shows that the ion-
neutral damping rate outside the cloud can be low because of the high
ionization fraction and low gas density; however, it is likely that in this
regime there are other damping mechanisms which we are not accounting for
(Hopkins et al., 2022b) because of our focus on the denser molecular gas.

Finally, we note that this model neglects $D_{\parallel,\rm FLW}$; none of our
simulations consider both field-line wandering diffusion ($D_{\parallel,\rm
FLW}$) and pitch-angle scattering diffusion ($D_{\parallel,\rm micro}$), since
this is not easy to model mathematically. In this section, we have shown that
the assumption that streaming dominates the CR transport is valid ($v_{\rm
st}\gg v_{\rm diff}$); however, the assumption that $D_{\parallel,\rm
micro}\ll D_{\parallel,\rm FLW}$ used in Section 3.1.3 is not always the case.
For the remainder of this study we refer only to the diffusion coefficient
$D_{\parallel}$ rather than distinguishing between the two.

####  2.2.7 Initial CR Configuration

The simulations assume an initially uniform background $\epsilon_{\rm
cr},_{\rm med}$ for the CR energy density in the diffuse medium outside the
cloud. Integrating the CR flux spectrum in the solar neighborhood
appropriately for the energy density and ionization rate gives values of
$\epsilon_{\rm cr}\approx 1\rm eV/cm^{3}$ and $\zeta\approx 1.6\times
10^{-17}\rm s^{-1}$ respectively (Cummings et al., 2016; Stone et al., 2019);
therefore, GIZMO assumes a CRIR of $1.6\times 10^{-17}\rm s^{-1}$ per $\rm
eV/cm^{3}$. We choose our fiducial value of $\epsilon_{\rm cr},_{\rm med}=1\rm
eV/cm^{3}$ to be consistent with these measurements and previous measurements
of CR ionization rates of $\zeta\approx 10^{-17}\rm s^{-1}$ in star-forming
regions (Dalgarno, 2006). Inside the cloud, we start with 10 % of the medium
CR energy density ($\epsilon_{\rm cr},_{\rm cloud}=0.1\rm eV/cm^{3}$) to
reflect probable attenuation in the dense gas (Padovani et al., 2009, 2020).
However, low resolution tests show that the simulations progress similarly
independently of whether the initial value is $\epsilon_{\rm cr},_{\rm
cloud}=1\rm eV/cm^{3}$ or $0.1\rm eV/cm^{3}$.

Analytic calculations suggest that the CR energy density should be lower in
MCs where the gas density is higher (Padovani et al., 2009); however, analysis
of linewidth observations along sightlines towards diffuse clouds have
suggested a range of CRIRs that exceed this value (Indriolo & McCall, 2012;
Indriolo et al., 2015). The mean value in the Milky Way is $\zeta\approx
10^{-16}\rm s^{-1}$, an order of magnitude higher than the canonically assumed
value (Indriolo et al., 2015). Additionally, higher ionization rates towards
more active star-forming regions in the center of our galaxy, other galaxies
(Rivilla et al., 2022; Behrens et al., 2022; Holdship et al., 2022) and
towards supernova remnants (SNRs) (Ceccarelli et al., 2011) have been observed
(see Padovani (2023) for a review). To reflect these uncertainties, we also
perform one simulation with initial values of $\epsilon_{\rm cr},_{\rm
med}=10\rm eV/cm^{3}$ and $\epsilon_{\rm cr},_{\rm cloud}=1\rm eV/cm^{3}$.

In total, we run six different simulations varying the CR diffusion
coefficient and initial CR energy density. These are described in Table 2. We
note that run M2e3_noCRT does not explicitly model CR transport but instead
assumes a fixed CR energy density of $\epsilon_{\rm cr}=1\rm eV/cm^{3}$ (or a
CRIR of $\zeta\approx 1.7\times 10^{-17}\,\rm s^{{}_{1}}$), which is used only
to calculate the heating and cooling of the gas.

Simulation | $M_{\rm cloud}[M_{\odot}]$ | $\epsilon_{\rm cr},_{\rm med}[\rm eV/cm^{3}]$ | $D_{\rm||}[\rm cm^{2}/s]$
---|---|---|---
M2e3_E1_D8e25 | $2\times 10^{3}$ | 1 | $8.33\times 10^{25}$
M2e3_noCRT | $2\times 10^{3}$ | 0 | N/A
M2e3_E10_D8e25 | $2\times 10^{3}$ | 10 | $8.33\times 10^{25}$
M2e3_E1_D8e24 | $2\times 10^{3}$ | 1 | $8.33\times 10^{24}$
M2e3_E1_D8e26 | $2\times 10^{3}$ | 1 | $8.33\times 10^{26}$
M2e3_E1_DVAR | $2\times 10^{3}$ | 1 | Variable
Table 2: Details of the variations we use in the initial conditions and
physics of our simulations. Columns are the simulation name, cloud mass,
initial CR energy density (outside the cloud), and CR diffusion coefficient.

##  3 Results

###  3.1 Cloud Evolution

####  3.1.1 Overview

We begin with a discussion of the cloud evolution for the fiducial cloud
including CRT, M2e3_E1_D8e25, and a simulation with the same conditions, which
assumes a uniform CRIR but no explicit CRT, M2e3_noCRT. Figure 1 shows the
time evolution of the projected gas density, density-weighted gas temperature,
density-weighted 1D velocity dispersion, and density-weighted CR energy
density of the M2e3_E1_D8e25 run, with the positions of stars superimposed in
blue. The 1D velocity dispersion for a cell is calculated from the $z$
component of the gas velocity $v_{z}$ as $\sigma_{1\rm
D}=\sqrt{(v_{z}-\bar{v_{z}})^{2}}$, where $\bar{v_{z}}$ is the mean value of
$v_{z}$ for all cells along the line of sight. The evolution of the gas in the
calculation including CRT follows a similar pattern to that of previous
STARFORGE simulations (Grudić et al., 2022; Guszejnov et al., 2021, 2022). At
early times before star formation starts, the initial turbulence in the cloud
leads to the formation of dense filaments, which host gravitationally unstable
cores. By 1.9 Myr ($\approx t_{\rm ff}$), some of the cores have collapsed and
formed the first stars (column 1).

Figure 1: Projected density (top row), density-weighted temperature (second
row), density-weighted 1D velocity dispersion (third row) and density-weighted
CR energy density (bottom row) of the M2e3_E1_D8e25 simulation at four
different times between $t_{\rm ff}-2.5t_{\rm ff}$ ($\sim 2-5$ Myr). The SFE
at each time is given in the bottom right corner. Sink particles are plotted
in blue, sized by their mass.

At $\sim 2.9$ Myr ($\approx 1.5t_{\rm ff}$), the rate of star formation peaks
(column 2). The forming protostars launch high-velocity outflows, which slow
star formation and stir turbulence in the cloud. By 4.3 Myr ($\approx
2.2t_{\rm ff}$), massive stars ($>8M_{\odot}$) launch stellar winds and
produce bubbles of ionizing radiation. This creates a central cavity, which
causes the surrounding gas to expand (column 3). By $2.5t_{\rm ff}$, feedback
has completely disrupted the cloud and the cloud is dispersed.

The bottom row of Figure 1 follows the evolution of the CR energy density. The
cloud boundary is clearly visible because the CR energy losses happen at a
faster rate in the cloud material than in the external medium. However, the CR
energy density is smoother and exhibits less filamentary substructure than the
gas density. The evolution of the CR energy density inside the cloud is shown
quantitatively in Figure 2, which plots the evolution of the median (solid)
and mean (dashed) CR energy density as a function of time for various gas
densities. Inside the cloud (blue, purple, and black lines), the CR energy
density is lower and declines more steeply than outside the cloud (cyan);
however, it is largely independent of gas density. The only significant
variations before $\sim 1t_{\rm ff}$ are caused by artificial trapping in the
dense gas (Appendix A), which creates slightly elevated median values in the
high density gas (black) (Hopkins et al., 2020). The uniformity of the CR
distribution indicates that CRs are losing most of their energy before
penetrating into the filaments where stars are forming.

In all of our simulations massive stars reach the minimum mass for stellar
wind launching between $\sim 1.5-2t_{\rm ff}$. Except for in the M2e3_noCRT
simulation, the massive stars inject CR energy and momentum as they lose mass
(Section 2.1.1), causing the mean CR energy density in the cloud to rise. Many
stars undergo accretion which rapidly increases their mass; consequently, the
CR energy injected by each source is nearly constant or increasing. Figure 2
shows that by $\sim 1.75-2t_{\rm ff}$ there is enough CR injection to sharply
increase the median CR energy density in the cloud, especially in the denser
gas where the stars are forming (black). Once the stellar winds rarefy the gas
the massive stars are no longer embedded, which causes the CR energy density
in the dense gas to decline.

Finally, the dip in the median CR energy density at low densities (cyan) at
$\sim 2.5t_{\rm ff}$ occurs because the cloud gas disperses faster than the
CRs diffuse into this gas. The energy density dips to match the values at
$n\gtrsim 500-10^{3}\rm cm^{-3}$ (blue and purple), reflecting the CR levels
in the regions previously at these gas densities. Once the background CRs
stream into the low density dispersing gas, the median CR energy density rises
again.

Figure 2: Median (solid) and mean (dashed) CR energy density (left axis) and
the corresponding CRIR (right axis) for gas in the indicated density range or
lower limit for the M2e3_E1_D8e25 simulation. The initial number density
inside the cloud is $n_{\rm cloud}\approx 300\rm cm^{-3}$; therefore, we
choose an upper cutoff of 30 $\rm cm^{-3}$ to represent gas outside the cloud.
The cutoffs of $500\rm cm^{-3}$ and $10^{3}\rm cm^{-3}$ are chosen from
various definitions to define the cloud boundary, while $10^{4}\rm cm^{-3}$ is
a typical density of filamentary substructure and star formation.

For comparison, Figure 3 shows the evolution of the gas properties of the
M2e3_noCRT run. In this simulation a uniform CR heating rate is used which is
equivalent to a CR energy density of $\epsilon_{\rm cr}=1\rm eV/cm^{3}$
($S_{\rm gas}=\Lambda_{\rm gas}e_{\rm cr}=\Lambda_{\rm gas}\epsilon_{\rm
cr}/n_{\rm cr}$ with $\Lambda_{\rm gas}$ given in Equation 8), while in the
M2e3_E1_D8e25 run $\epsilon_{\rm cr}<1\rm eV/cm^{3}$ throughout most of the
cloud. Figure 3 shows that despite the differences in the CRIRs between the
two runs, the effect of including CRT on the gas structure and overall star
formation is minimal. Differences in the gas structure past $\sim$ 4 Myr can
largely be attributed to variations in stellar feedback from the different
sink particle distributions.

Figure 3: Same as Figure 1 except for the M2e3_noCRT simulation. In this
simulation, a uniform CR energy density of $\epsilon_{\rm cr}=1\rm eV/cm^{3}$
($\zeta\approx 2\times 10^{-17}\rm s^{-1}$) is used in the heating and cooling
functions.

The effect of including CRT on the gas temperature and velocity dispersion is
also minimal. Figure 4 shows the evolution of the median temperature (top) and
velocity dispersion (bottom) as a function of time at various gas densities
for all simulations. The velocity dispersion for a cell is calculated from the
three components of the gas velocity $v_{x}$, $v_{y}$, and $v_{z}$ as
$\sigma_{3\rm
D}=\sqrt{(v_{x}-\bar{v_{x}})^{2}+(v_{y}-\bar{v_{y}})^{2}+(v_{z}-\bar{v_{z}})^{2}}$,
where $\bar{v_{x}}$, $\bar{v_{y}}$, and $\bar{v_{z}}$ are the mean values of
$v_{x}$, $v_{y}$, and $v_{z}$ respectively for all cells. Figure 4 shows the
temperature in the M2e3_E1_D8e25 simulation (black) is actually a couple
degrees colder than the M2e3_noCRT simulation (cyan) due to the lower CR
heating. Similarly, the median velocity dispersion for the two simulations is
virtually identical until $t\approx t_{\rm ff}$, indicating CRs are not
significantly affecting the gas dynamics in the M2e3_E1_D8e25 run. Beyond
$t\approx 1-1.5t_{\rm ff}$, the cloud begins to disperse due to feedback and
differences in these quantities are generally due to minor variations in the
stellar mass distribution and cannot be attributed to CRs.

Figure 4: Median temperature(top) and median velocity dispersion (bottom)
versus time for all gas with number density $n>1000\rm cm^{-3}$ for all
simulations.

####  3.1.2 Star Formation in a High CR Environment

We now compare the fiducial M2e3_E1_D8e25 CR simulation with the
M2e3_E10_D8e25 simulation, which follows a cloud embedded in an environment
with a 10 times higher initial CR energy density. The top two rows of Figure 5
show the projected CR energy density for all CRT simulations at $t=1.5t_{\rm
ff}$ (approximately the time of peak star formation). The CR energy density in
the M2e3_E10_D8e25 cloud declines at a similar rate to that of the
M2e3_E1_D8e25 cloud. Despite the higher CR environment, the CR energy density
inside the M2e3_E10_D8e25 cloud is relatively low, smaller than the typical 1
eV cm-3 generally assumed for clouds. In both runs, by $t=1.5t_{\rm ff}$ the
CR energy density in the cloud has declined by over an order of magnitude from
its initial value, with values of $\epsilon_{\rm cr}\lesssim 10^{-2}\rm
eV/cm^{3}$ and $\epsilon_{\rm cr}\lesssim 10^{-1}\rm eV/cm^{3}$ for the
M2e3_E1_D8e25 and M2e3_E10_D8e25 clouds, respectively.

Figure 5: Top: Density-weighted CR energy density for all simulations at
$t=1.5t_{\rm ff}$ ($\approx 3$ Myr). Sink particles are plotted in blue, sized
by their mass. The SFE at this time for each simulation is shown in the bottom
right corner. For the M2e3_noCRT panel, a CR energy density of $\epsilon_{\rm
cr}=1\rm eV/cm^{3}$ is shown, corresponding to the assumed ionization rate of
$\zeta\sim 1.7\times 10^{-17}\rm s^{-1}$. Bottom: Projected gas density for
all simulations at $t=1.5t_{\rm ff}$ ($\approx 3$ Myr). Sink particle
locations are indicated by blue stars, sized by their mass. The current SFE
for each simulation is shown in the bottom right corner. Despite significant
differences in the CR treatments, the star formation activity and distribution
is qualitatively similar.

The higher CR energy does, however, affect the rate at which the cloud
collapses and the time star formation commences. The bottom two rows of Figure
5 show the projected gas density for all simulations at $t=1.5t_{\rm ff}$
($\approx 3$ Myr). At this time, the M2e3_E10_D8e25 cloud (top right) is more
compact and there are more stars compared to the M2e3_E1_D8e25 and M2e3_noCRT
clouds (top left and middle, respectively). This occurs because the relatively
high CR energy density outside the cloud acts as an additional pressure on the
cloud. While this relative pressure imbalance does not appreciably affect the
gas evolution in the M2e3_E1_D8e25 simulation, it is sufficient to accelerate
the dynamical evolution of the cloud in the M2e3_E10_D8e25 simulation.
Consequently, by $1.5t_{\rm ff}$ more than twice as much of the cloud mass in
the M2e3_E10_D8e25 run has formed (or been accreted onto) stars.

Although the higher CR background causes a faster initial collapse than in the
fiducial M2e3_E1_D8e25 simulation, the other gas properties are qualitatively
similar. The top panel of Figure 4 shows that the gas in the M2e3_E10_D8e25
simulation cools to a slightly lower median temperature than the M2e3_noCRT
simulation (cyan) once the internal CR energy density declines. However, the
cloud temperature is higher after $\sim 1t_{\rm ff}$ because the elevated star
formation causes stronger and earlier heating. The bottom panel shows that the
median gas velocity dispersion is also higher after $\sim 0.5t_{\rm ff}$, a
consequence of the earlier star formation and subsequent impact of feedback
within the cloud. Overall, due to the strong CR attenuation in the dense gas,
star formation in the high CR environment proceeds relatively similarly to
star formation occurring in a typical CR environment.

####  3.1.3 Impact of the CR Diffusion Coefficient on Star Formation

We now turn to a comparison with our simulations with different CR diffusion
coefficients. The second row of Figure 5 shows the projected CR energy density
at $t=1.5t_{\rm ff}$ for the low diffusion coefficient M2e3_E1_D8e24, high
diffusion coefficient M2e3_E1_D8e26, and variable diffusion coefficient
M2e3_E1_DVAR simulations. The CR energy density in the cloud is even lower for
the M2e3_E1_D8e24 (bottom left) run than for the M2e3_E1_D8e25 run, because
the CRs scatter more, and consequently suffer greater attenuation by the time
they are able to random-walk into the cloud interior. In the M2e3_E1_D8e26 run
(bottom middle), the CRs scatter less frequently and lose comparatively little
energy, so the cloud and medium reach an equilibrium value of $\epsilon_{\rm
cr}\approx 1\rm eV/cm^{3}$ (second row middle panel of Figure 5).

The right panel in the second row of Figure 5 shows that the CR energy density
in the M2e3_E1_DVAR simulation evolves qualitatively similarly to the
M2e3_E1_D8e24 simulation. Computing the microphysical diffusion coefficient
$D_{\parallel,\rm micro}$ typically gives lower average values than the value
of $D_{\parallel,\rm FLW}$ used in the M2e3_E1_D8e25 simulation. Figure 6
(left) shows the median diffusion coefficient at different gas densities in
the M2e3_E1_DVAR simulation. Outside the cloud (cyan line), the low gas
densities and high ionization result in inefficient ion-neutral damping
(Equation 14) so the diffusion coefficient is typically equal to the minimum
allowed value. Inside the cloud (purple and blue lines), the values of
$D_{\parallel}$ increase as the cloud collapses because the higher densities
create more efficient ion-neutral damping, which damps the Alfvèn waves and
decreases the CR scattering rate. The right panel of Figure 6 shows a 2D
histogram of the amplitude of the CR diffusion coefficient as a function of
gas density at $t=0.5t_{\rm ff}$. Before $\sim 1t_{\rm ff}$ the diffusion
coefficient throughout most of the cloud is between the M2e3_E1_D8e25 value
and the M2e3_E1_D8e24 value except at the highest gas densities. At densities
$n<1000\rm cm^{-3}$, a large fraction of the simulation space computes our
lower cutoff value $D_{\parallel}=8.33\times 10^{24}\rm cm^{2}/s$, which we
impose due to the additional contribution of multiple wave damping mechanisms
at lower gas densities (Section 2.2.6). For most of the cloud collapse, the
high scattering rate causes the CRs to lose energy as in the M2e3_E1_D8e24
simulation.

Figure 6: Left: Median CR diffusion coefficient for gas cells in the indicated
density range or lower limit for the M2e3_E1_DVAR simulation. Dashed
horizontal lines show the constant values of the diffusion coefficient used in
the M2e3_E1_D8e25 (black), M2e3_E1_D8e24 (red), and M2e3_E1_D8e26 (brown)
simulations. The drop in the median value at $\approx 1.75t_{\rm ff}$ is due
to a numerical error in setting the lower cutoff; however, since this occurs
after the cloud has collapsed and the CRs have lost most of their energy it
does not affect our results. Right: Total gas mass as a function of gas
density and CR diffusion coefficient $D_{\rm sc}$ (not including resolved
streaming, Section 2.2.6) for the M2e3_E1_DVAR simulation at $t=0.5t_{\rm
ff}$.

For all three of these variations, the effect of CRT on the cloud evolution
and gas properties is minimal. The bottom row of Figure 5 shows that the gas
density at $t=1.5t_{\rm ff}$ looks virtually identical to the gas density of
the M2e3_E1_D8e25cloud. As expected, the temperature of the M2e3_E1_D8e24 and
M2e3_E1_DVAR simulations shown in Figure 4 (red and magenta, respectively) is
lower than for the M2e3_E1_D8e25 simulation. The temperature in the
M2e3_E1_D8e26 simulation is similar to that of the M2e3_E1_D8e25 simulation,
because both have a value of $\epsilon_{\rm cr}\approx 1\rm eV/cm^{3}$
throughout the cloud. The velocity dispersion for all three is comparable,
since star formation proceeds similarly.

###  3.2 Star Formation Properties

####  3.2.1 Star Formation History of the Cloud

We now turn to various star formation metrics that can be obtained from the
simulations and compared to observations. The star formation efficiency (SFE)
is defined as the percentage of the initial cloud mass that has been converted
to stars:

${\rm SFE}(t)=M_{*}(t)/M_{0},$ (17)

where $M_{*}(t)$ is the total stellar mass at time $t$. Figure 7, left, shows
the evolution of the SFE for all runs. Qualitatively, the SFE evolves
similarly in all cases. Once star formation starts, the SFE undergoes a period
of rapid growth as the cloud collapses and the densities increase. The rate of
star formation peaks at $\approx 1.5t_{\rm ff}$. At this time, feedback begins
to disperse the cloud, causing the slope of the SFE to decrease. By $\approx
2t_{\rm ff}$, the SFE has plateaued and the cloud is dispersing.

Figure 7 shows that the inclusion of CRT does not have a significant impact on
the star formation history of the cloud for any of our simulations with
$\epsilon_{\rm cr}=1\rm eV/cm^{3}$. Up until $\approx 1.5t_{\rm ff}$, the SFE
looks very similar for all of these runs, and minor variations after this can
be attributed to differences in star formation feedback. The SFE in the
M2e3_E10_D8e25 cloud (blue lines) begins the period of steep rise slightly
earlier than the other clouds due to the external CR pressure on the cloud
discussed in Section 3.1.2; however, once star formation starts it proceeds at
roughly the same rate for all of the clouds once the background CR level has
declined. For all of our clouds, variations in the early evolution of the SFE
are smaller than variations caused by different initial turbulent realizations
with the same physics (Guszejnov et al., 2022).

To evaluate the significance of the variation between runs, we compare our
results to a suite of 103 simulations of $M=2,000$ M⊙ clouds evolved with
different turbulent seeds and no CRT. 100 of these clouds were run with
$E_{\rm mag}/E_{\rm grav}=0.01$ and were presented in Grudić et al. (2023),
while the other three were run with $E_{\rm mag}/E_{\rm grav}=0.1$ but are
otherwise identical. We take the mean value of the final SFE from the three
higher-magnetization clouds plus our M2e3_noCRT run as the ‘standard’ final
SFE with no CRT, and assume the standard deviation in the final SFE calculated
from the suite of 100 lower magnetization clouds is representative of the
variation. The $\pm 1\sigma$ region ($\rm SFE\approx 7.42\pm 1.77$), is
indicated as a grey shaded area in the left panel of Figure 7.

The second column of Table 3 shows the SFE at the final time $4t_{\rm ff}$.
The final SFE of our M2e3_noCRT cloud is the lowest. All of our $\epsilon_{\rm
cr}=1\rm eV/cm^{3}$ clouds have slightly higher values, which could be due to
the extra confining CR pressure; the M2e3_E1_D8e26 and M2e3_E1_DVAR are within
1$\sigma$ while the M2e3_E1_D8e25 and M2e3_E1_D8e24 clouds are within
2$\sigma$. The M2e3_E10_D8e25 cloud has a more than 2$\sigma$ greater final
SFE.

Figure 7: Left: Evolution of the star formation efficiency (SFE) versus time
for all simulations. Shaded regions are the $1\sigma$ (dark gray) and
$2\sigma$ (light gray) final SFE values obtained from a set of 100 2,000 M⊙
cloud simulations without CRT that vary the initial turbulent seed. Right:
Median mass for $M_{*}\geq 0.1M_{\odot}$ versus time for all runs. Shaded
regions are the $1\sigma$ (dark gray) and $2\sigma$ (light gray) final $M_{\rm
med}$ values obtained from a set of previous simulations run without CRT,
varying the initial turbulent seed.

####  3.2.2 Final Stellar Mass Distribution

Figure 7, right, shows the evolution of the median mass of all sink particles
with $M>0.1M_{\odot}$. As in the left panel, the shaded region shows the mean
of the final values from the four $E_{\rm mag}/E_{\rm grav}=0.1$ simulations
(including the M2e3_noCRT run) with the $1\sigma$ region from the 100 $E_{\rm
mag}/E_{\rm grav}=0.01$ simulations. There is no significant difference in the
evolution of $M_{\rm med}$ between the CRT and M2e3_noCRT simulations, except
for the M2e3_E1_D8e25 simulation. The median mass for that run experiences a
sharp jump between $\sim 1.25-1.5t_{\rm ff}$ and then rises again sharply at
$\sim 1.8t_{\rm ff}$. While the earlier jump may be due to statistical
variation, the latter may be correlated with the CR injection by stellar winds
discussed in Section 3.1, which was strongest in our M2e3_E1_D8e25 simulation.
It occurs at approximately the same time as the sharp increase in the median
CR energy density in the cloud from stellar wind injection shown in Figure 2.
The additional CR energy may affect the gas temperature and pressure, which
would effect the stellar mass distribution. A larger statistical study is
required to draw firm conclusions.

The final value of $M_{\rm med}$ at $t=4t_{\rm ff}$ is shown in the third
column of Table 3. Except for the M2e3_E1_D8e25 simulation, all of our
$\epsilon_{\rm cr}=1\rm eV/cm^{3}$ runs are within $1\sigma$. The
M2e3_E10_D8e25 simulation is slightly higher, within $2\sigma$.

At the final time, the $\epsilon_{\rm cr}=1\rm eV/cm^{3}$ simulations have
each formed $\sim 120-170$ stars. Figure 8 shows the stellar mass
distributions at $4t_{\rm ff}$ (left) and their normalized cumulative
distribution functions (right). The shaded region represents the $1\sigma$
region in each mass bin from the 100 $E_{\rm mag}/E_{\rm grav}=0.01$
simulations for a mass distribution with 136 stars and a high mass slope
$\alpha=-1$ in $dN/d\ln M$, shown by the gray dashed line. Qualitatively, the
IMF for all of the $\epsilon_{\rm cr}=1\rm eV/cm^{3}$ simulations appears
reasonably well described by a Chabrier IMF (Chabrier, 2005). Prior work
indicates that protostellar jets are responsible for setting the mass scale of
the IMF (Guszejnov et al., 2021). The high mass slope is flatter than the
canonical Chabrier value of $\alpha=-1.35$ but is consistent with previous
STARFORGE simulations, which find $\alpha=-1$, possibly due to some missing
physics, such as disk mediated accretion and fragmentation, that enhance the
number of massive stars (Grudić et al., 2022; Guszejnov et al., 2022). The
peak of the M2e3_E1_D8e25 simulation is higher than that of the other
simulations because of the higher median mass (Figure 7). Except for this, we
find no qualitative difference between the IMFs of the $\epsilon_{\rm cr}=1\rm
eV/cm^{3}$ CRT simulations and the IMF of the M2e3_noCRT simulation.

The M2e3_E10_D8e25 simulation produced twice as many stars as the other
simulations ($\sim$ 220). Figure 8 shows that the number of stars formed with
$M\lesssim 0.1M_{\odot}$ is larger than the other simulations, resulting in a
visual peak in the IMF at $\sim 0.1M_{\odot}$ rather than $\sim 0.3M_{\odot}$.
However, the median mass of the M2e3_E10_D8e25 simulation above $0.1M_{\odot}$
is higher than all of the other simulations except the M2e3_E1_D8e25
simulation (Figure 7) because this simulation also formed more high-mass
stars. Indeed, the high-mass slope of the M2e3_E10_D8e25 simulation in Figure
8 appears flatter than $\alpha=-1$ and the distribution is well outside the
$1\sigma$ region at $M>1M_{\odot}$. Higher resolution simulations with a
larger statistical sample of stars are required to confirm this trend.

Figure 8: Left: Stellar initial mass function (IMF) for all runs at $t=4t_{\rm ff}$. The shaded region at $M<0.1M_{\odot}$ indicates the low-mass incompleteness region (Grudić et al., 2021). Dashed lines show the Chabrier IMF (Chabrier, 2005) with the standard slope of $\alpha=-1.35$ (black) and the slope obtained from previous STARFORGE simulations, $\alpha=-1$ (gray) (Grudić et al., 2022; Guszejnov et al., 2022). Shaded region shows the $1\sigma$ values in each mass bin obtained from a set of 100 2000 $M_{\odot}$ cloud simulations without CRT. Right: Normalized CDF of the IMF for all simulations. Simulation Name | SFE (%) | $M_{\rm med}(M>0.1M_{\odot})$ | $M_{\rm max}$
---|---|---|---
M2e3_E1_D8e25 | 9.53 | 0.50 | 29.3
M2e3_noCRT | 7.35 | 0.38 | 30.6
M2e3_E10_D8e25 | 12.90 | 0.44 | 22.5
M2e3_E1_D8e24 | 9.30 | 0.37 | 23.3
M2e3_E1_D8e26 | 9.17 | 0.38 | 22.7
M2e3_E1_DVAR | 8.19 | 0.39 | 19.1
Table 3: Final values for the SFE, median mass above the incompleteness
region, and maximum mass for all simulations at $4t_{\rm ff}$.

###  3.3 Ionization

As shown in Figure 5, we find that in a Milky Way like environment, for low to
moderate values of the CR diffusion coefficient the CRIR is roughly uniform
throughout the cloud. The median value ranges between $\zeta\approx 4\times
10^{-20}-2\times 10^{-19}\rm s^{-1}$ during the main star forming epoch of the
simulations ($t\approx 1-1.5t_{\rm ff}$). Figure 9 shows a phase plot of the
total simulation mass as a function of the CR energy density and gas density
at $1.5t_{\rm ff}$. The CRIR is shown on the right, which we calculate from
the CR energy density assuming an ionization rate of $1.6\times 10^{-17}\rm
s^{-1}$ per $\rm eV/cm^{3}$ (Cummings et al., 2016). The blue line shows the
mass-weighted mean CR energy density and corresponding CRIR at each gas
density. The M2e3_E1_D8e25, M2e3_E1_D8e24 and M2e3_E1_DVAR clouds have a CRIR
that strongly peaks at a value of $\sim 0.5-1\times 10^{-19}$s-1, where some
regions with $n<10^{4}$ reach $\zeta\approx 10^{-18}$s-1. The M2e3_E1_D8e26
run has a nearly uniform CRIR of $\zeta\sim 10^{-17}$s-1 throughout the cloud.

The M2e3_E10_D8e25 simulation has elevated ionization rates throughout the
cloud compared to the M2e3_E1_D8e25 simulation. The median value of the CRIR
ranges from $\approx 1-6\times 10^{-18}\rm s^{-1}$ between 1-1.5 $t_{\rm ff}$.
Although the average CR energy density declines over time in both the
M2e3_E10_D8e25 and M2e3_E1_D8e25 clouds, the CR distribution is more variable
throughout the M2e3_E10_D8e25 cloud. Figure 9 shows that the typical cloud gas
in the M2e3_E1_D8e25 and M2e3_E10_D8e25 simulations has a characteristic
ionization rate of $\zeta\approx 10^{-19}\rm s^{-1}$ and $2\times 10^{-18}\rm
s^{-1}$, respectively, relatively independent of gas density. The blue line
indicates that the M2e3_E1_D8e25 cloud is strongly peaked at this value,
especially at $n\gtrsim 10^{4}\rm cm^{-3}$. In contrast, the M2e3_E10_D8e25
cloud has an average CRIR between $\approx 10^{-18}-10^{-17}\rm s^{-1}$
depending on gas density, indicating that the distribution is not strongly
peaked and has a range of values above and below the median CRIR. There are
regions with $\zeta\gtrsim 10^{-17}\rm s^{-1}$ at all gas densities.

Figure 9: The total mass of cells as a function of gas density and CR energy
density for all simulations at $t=1.5t_{\rm ff}$ ($\approx 3$ Myr). The blue
line shows the mass-weighted mean CR energy density at each gas density. On
the right hand is the corresponding CR ionization rate. The horizontal blue
line in the M2e3_noCRT panel is the MW CR energy density ($\epsilon_{\rm
cr}=1\rm eV/cm^{3}$) assumed when CRT is not used.

Figure 10 shows a phase plot of the ionization fraction for all runs at
$t=1.5t_{\rm ff}$. We compute the ionization fraction following the convention
used in GIZMO from the electron fraction $f_{e}$, neutral hydrogen fraction
$f_{\rm H}$, and HII fraction $f_{\rm HII}$ as MAX($f_{e},(1-f_{\rm H})$,
$f_{\rm HII}$). At this time, ionization produced by stellar feedback is still
minimal, so the degree of ionization reflects the local CR energy density. The
pale yellow line in the M2e3_noCRT panel shows the contribution assuming a
fixed CRIR of $1.6\times 10^{-17}\rm s^{-1}$, which is used in the M2e3_noCRT
simulation. The dark patches in the upper left are produced by ionization
caused by radiative feedback, while the area in the lower left shows the
contribution from CR ionization. Since the CRIR is low for most of the runs,
the contribution from CRs to the overall gas ionization is minor compared to
the ionization caused by radiative feedback once star formation is underway.

Figure 10: Total mass of cells as a function of gas density and ionization
fraction for all simulations at $t=1.5t_{\rm ff}$ ($\approx 3$ Myr). The
yellow regions in each panel correspond to CR ionization, while the dark
patches are from radiative feedback. The pale yellow line in the M2e3_noCRT
panel shows the contribution assuming a fixed CRIR of $1.6\times 10^{-17}\rm
eV/cm^{3}$, which declines with gas density as $1/\sqrt{n_{H}}$. The
descending line in all CRT simulations with ionization fraction $<1e-8$ and
density $<1e3$ is due to a known bug in GIZMO where one particle per MPI rank
is assigned $\epsilon_{\rm cr}=0$. Given the fact that we use $2\times 10^{6}$
cells and only $\sim 200$ MPI ranks, this bug impacts very few cells and thus
does not have an effect on our results.

##  4 Discussion

###  4.1 Propagation of CRs in Molecular Clouds

In this section we discuss the implications of our results for the propagation
of CRs within MCs.

####  4.1.1 Role of the CR-Driven Streaming Instability

We find that for all of the simulations except the M2e3_E1_D8e26 simulation,
the CR energy density in the cloud is significantly lower than it is outside
the cloud. The CR energy density in the cloud declines quickly as the cloud
collapses and is not replenished by CRs from outside the cloud. To understand
this result, we consider the relevant timescales for CR propagation and
compare them to the loss timescales. The following expressions describe the
timescale for each transport mechanism to propagate CRs a distance $L$ (e.g.,
Habegger et al., 2023):

$\tau_{\rm adv}=L/v_{\rm adv}\approx 1L_{\rm pc}v_{\rm km}^{-1}\rm Myr,$ (18)
$\tau_{\rm str}=L/v_{\rm st}\approx 1L_{\rm
pc}B_{\mu}^{-1}\rho_{-21}^{1/2}f_{-7}^{1/2}\rm kyr,$ (19) $\tau_{\rm
diff}=L^{2}/D_{||}\approx 3L_{\rm pc}^{2}D_{26}^{-1}\rm kyr,$ (20)

for advection, streaming, and diffusion, respectively, where we assume $v_{\rm
adv}=|\mathbb{u}|$ (the gas velocity), $v_{\rm st}=v_{\rm Ai}$, and $L_{\rm
pc},v_{\rm km},B_{\mu},\rho_{-21},f_{-7}$ and $D_{26}$ are the length scale,
advection velocity, root-mean-square magnetic field, gas density, ionization
fraction, and the diffusion coefficient in units of 1 pc, 1 km/s, $1\mu G$,
$10^{-21}\rm g/cm^{3}$, $10^{-7}$, and $10^{26}\rm cm^{2}/s$, respectively.

While they propagate, the CRs lose energy through streaming and collisional
losses. If the amount of time for CRs to lose their energy through either one
of these mechanisms is considerably shorter than the time for the CRs to
propagate through the cloud, then the CRs will quickly lose energy, and thus
the CR energy density will be low inside the cloud. The energy loss timescale
from collisions is obtained from Equation 9:

$\tau_{\rm cat,loss}=\Lambda_{\rm cr}^{-1}=\left(7.51\times
10^{-16}\right)^{-1}\left(\frac{n_{H}}{\rm cm^{-3}}\right)^{-1}\rm s.$ (21)

The streaming loss term is calculated consistently in the simulations as
$\Gamma_{\rm st}=-(\bar{\nu}_{s}/c^{2})v_{\rm st}F_{e,\rm cr}$, where
$\Gamma_{\rm st}=\Lambda_{\rm st}e_{\rm cr}$ and $\Lambda_{\rm st}$ is the
streaming energy loss rate. The streaming loss timescale can be approximated
as (Chan et al., 2019):

$\tau_{\rm str,loss}=\Lambda_{\rm st}^{-1}=\left(\frac{1}{\gamma_{\rm
cr}-1}\right)|\hat{B}\cdot\hat{\nabla}\epsilon_{\rm
cr}|^{-2}\left(\frac{\ell_{\rm cr}}{v_{\rm st}}\right)\rm s,$ (22)

where $\hat{\nabla}\epsilon_{\rm cr}$ is the unit vector of the CR energy
density gradient $\nabla\epsilon_{\rm cr}$ and $\ell_{\rm cr}=\epsilon_{\rm
cr}/\nabla\epsilon_{\rm cr}$ is the CR gradient scale length.

The advection timescale (Equation 18) reflects the speed of the non-
relativistic gas and thus represents the gas crossing time in the simulation.
Since our simulations use a comoving Lagrangian frame which moves with the
gas, it is independent of the choice of $\tilde{c}$. The four other timescales
($\tau_{\rm str}$, $\tau_{\rm diff}$, $\tau_{\rm cat,loss}$ and $\tau_{\rm
str,loss}$), however, represent CR propagation and losses away from the gas
and thus these timescales are a factor $\Gamma_{\rm RSOL}=c/{\tilde{c}}$
slower due to the RSOL approximation (Hopkins et al., 2022a). For our value of
$\tilde{c}=300$ km/s, $\Gamma_{\rm RSOL}=1000$. This caveat does not influence
the conclusions of this section because we are concerned with the relative
importance of these timescales and all of the timescales scale proportionally.
We discuss the implications of the choice of $\tilde{c}$ further in Section
4.3 and in Appendix A.

Figure 11 (bottom) shows the five timescales as a function of gas density for
the M2e3_E1_D8e25 (solid), and M2e3_E1_D8e26 (dashed) simulations at
$0.5t_{\rm ff}$, where the relevant gas quantities used to calculate the
transport and loss timescales are displayed in the top panel. To reflect the
behavior of the simulations we scale all timescales except for $\tau_{\rm
adv}$ by a factor of $\Gamma_{\rm RSOL}$. For the M2e3_E1_D8e25 simulation
transport due to streaming dominates over advective and diffusive transport
throughout the cloud. This is especially true at the highest densities
($n\gtrsim 100\rm cm^{-3}$), where the ionization fraction is $f_{\rm
ion}<10^{-7}$, resulting in $v_{\rm Ai}>1000$ km/s. The high streaming
velocity also produces efficient energy losses from the streaming instability.
Streaming losses dominate for gas densities $n<10^{4}\rm cm^{-3}$ and
collisional energy losses dominate at higher densities. Although the
magnitudes of the transport and loss timescales vary with the CR and gas
properties, all simulations except the M2e3_E1_D8e26 simulation have
significant streaming instability energy losses, which cause the CRs to lose a
significant amount of energy before propagating far in the cloud. To provide
additional evidence for this result, we run low-resolution tests where we
disable collisional losses and streaming transport and losses separately (see
Appendix B).

Figure 11: Top: Median values of the CR energy density $\epsilon_{\rm cr}$,
RMS magnetic field $B_{\rm rms}$, ionization fraction $f_{\rm ion}$, CR
gradient scale length $\ell_{\rm cr}=\epsilon_{\rm cr}/\nabla\epsilon_{\rm
cr}$, and streaming velocity $v_{\rm st}$ as a function of density as used in
the calculation of the streaming timescale $\tau_{\rm str}$ (Equation 19) and
the streaming loss timescale $\tau_{\rm str,loss}$ (Equation 22) at $0.5t_{\rm
ff}$. The vertical gray dotted lines correspond to gas densities
$n=1,10,100\rm cm^{-3}$. Bottom: Advection timescale $\tau_{\rm adv}$,
streaming timescale $\tau_{\rm str}$, diffusion timescale $\tau_{\rm diff}$,
streaming loss timescale $\tau_{\rm str,loss}$, and collisional loss timescale
$\tau_{\rm cat,loss}$ from Equations 18, 19, 20, 21 and 22, respectively, as a
function of gas density for the M2e3_E1_D8e25 (solid) and M2e3_E1_D8e26
(dashed) simulations at $0.5t_{\rm ff}$. The transport timescales $\tau_{\rm
adv}$, $\tau_{\rm str}$ and $\tau_{\rm diff}$ assume a characteristic length
scale of $L=3$ pc. The values of the streaming timescales $\tau_{\rm str}$ and
$\tau_{\rm str,loss}$ and the advection timescale $\tau_{\rm adv}$ are
calculated using the median values of gas and CR properties at each gas
density. All timescales except $\tau_{\rm adv}$ are scaled by $\Gamma_{\rm
RSOL}=c/\tilde{c}$. The vertical gray dotted lines correspond to gas densities
$n=1,10,100\rm cm^{-3}$.

The diffusive transport timescale is inversely proportional to the CR
diffusion coefficient. For the M2e3_E1_D8e26 simulation, diffusive transport
is comparable to streaming transport. Physically, the CRs scatter less and so
effectively propagate more efficiently. The streaming timescale is also
slightly higher due to the higher ionization fraction from the higher CR
energy density in the cloud. Finally and most significantly, the streaming
loss timescale is larger for the M2e3_E1_D8e26 simulation because the CR
distribution is more uniform. The cumulative effect is that in the
M2e3_E1_D8e26 simulation, CRs are able to propagate through the cloud without
losing a significant amount of energy.

Our simulations show that streaming instability energy losses dominate over
collisional losses for $\sim$ GeV CRs, which has important physical and
observational implications. The CR energy lost to the streaming instability
goes into exciting Alfvèn waves, which are quickly damped by ion-neutral
collisions that thermalize the energy. This means that the CR streaming
instability provides a stronger source of gas heating then CR ionization.
Furthermore, the CR energy is not lost to $\gamma$-ray production, especially
not at the cloud boundary where streaming instability energy losses are
strongest. This supports $\gamma$-ray observations from MCs that are in
agreement with the solar neighborhood CRIR of $\zeta\approx 10^{-17}\rm
s^{-1}$ (Krumholz et al., 2023) even though linewidth observations suggest
they have elevated CRIRs by approximately an order of magnitude (Indriolo et
al., 2015). In other words, our simulations support the conclusion that the
high CRIR values observed in Milky Way MCs do not come from an elevated total
CR background. Instead, they may result from additional sources of low-energy
ionizing CRs ($\sim$ MeV), which we discuss more in Section 4.2.1.

On the other hand, $\gamma$-ray measurements from starburst galaxies indicate
that they are proton-calorimeters, with all of their CR energy lost to
catastrophic collisions. Our results suggest that in these galaxies the
streaming speeds in MCs must be slower than in a Milky-Way environment because
otherwise the CRs lose all of their energy to streaming losses. The ionization
fraction is higher on average in starbursts than in Milky Way type galaxies,
with typical values from $\chi\approx 10^{-7}-10^{-4}$ (Krumholz et al.,
2020). In regions of higher ionization the ion-Alfvèn speed is lower and so
the CRs can get ‘trapped’ and lose their energy to hadronic losses rather than
streaming losses.

####  4.1.2 Exclusion of CRs from Molecular Clouds

Section 4.1.1 discusses how streaming instability energy losses cause the CRs
initialized in the cloud to quickly lose energy once the cloud starts
collapsing. However, this does not explain why CRs from outside the cloud do
not propagate into the cloud and offset the energy losses. In this section, we
discuss the effects that may be producing this ‘exclusion’ of CRs from MCs,
which has been noted in many previous studies (Kulsrud & Pearce, 1969;
Skilling & Strong, 1976; Cesarsky & Volk, 1978; Chandran, 2000; Padoan &
Scalo, 2005; Silsbee et al., 2018; Padovani & Galli, 2011; Padovani et al.,
2013; Everett & Zweibel, 2011; Bustard & Zweibel, 2021).

The primary effect keeping CRs from propagating into the cloud is that
streaming energy losses are even stronger at the cloud boundary than inside
the cloud. Figure 11 shows that at $0.5t_{\rm ff}$ the streaming loss
timescale of the M2e3_E1_D8e25 simulation drops sharply at the outer cloud
boundary region ($n\sim 1-10\rm cm^{-3}$) and reaches a minimum at the inner
cloud boundary region ($n\sim 10-100\rm cm^{-3}$) because of the sharp
gradient in the CR energy density. This indicates that CRs lose a significant
amount of energy as they propagate through the boundary and into the cloud. We
emphasize that this outcome is not due to the fact that we started the
simulations with $\epsilon_{\rm cr},_{\rm cloud}=0.1\epsilon_{\rm cr},_{\rm
med}$, and it occurs even in tests beginning with $\epsilon_{\rm cr},_{\rm
cloud}=\epsilon_{\rm cr},_{\rm med}$. Rather, the gradient in the CR energy
density always develops at the cloud boundary because of the drop in the
ionization fraction and resulting increase in the ion-Alfvèn velocity there
(Figure 11).

The degree of CR exclusion from the cloud interior also depends on the
diffusion coefficient. Across the entire cloud boundary region ($n\sim
1-100\rm cm^{-3}$), the diffusive transport timescale is comparable to that of
streaming transport. Physically, the CRs have a longer path length in the
M2e3_E1_D8e26 simulation, and gradients in the CR energy density are wiped out
faster. For the M2e3_E1_D8e26 simulation, the streaming loss timescale is
roughly uniform across the cloud boundary, because no sharp gradients in the
CR energy density develop.

The CR exclusion at the cloud boundary may also be enhanced by magnetic field
effects. Figure 12 shows slices of the gas density (top), CR energy density
(middle), and root-mean-squared (RMS) magnetic field $B_{\rm rms}$ (bottom) at
$0.5t_{\rm ff}$ through the middle of the domain for the M2e3_E1_D8e25 (left),
M2e3_noCRT (middle) and M2e3_E1_D8e26 (right) simulations. The magnetic field
in the cloud is up to an order of magnitude higher than the magnetic field in
the ambient medium, setting up a strong magnetic gradient. However, the
magnetic field behavior at the cloud boundary is complex: there is a region
with slightly elevated magnetic field that is a few $\mu G$ higher than the
medium value at $n\sim 1\rm cm^{-3}$, followed by a pocket of low-magnetic
field that coincides with the boundary region ($n\sim 1-100$cm-3). This
magnetic behavior is similar in all three simulations and coincides with the
phase transition from ionized to neutral gas. In the simulations with CRT, the
highly magnetized ring aligns with a peak in the CR energy density (shown by
the yellow rim). In the M2e3_E1_D8e26 simulation, a ring of enhanced CR energy
density is faintly visible, but the CR distribution is mostly uniform
throughout the cloud and medium. The CRs have a longer path length in the
M2e3_E1_D8e26 simulation, which allows them to more easily scatter over
magnetic field asymmetries into the cloud.

Figure 12: Gas number density (top), CR energy density (middle) and RMS
magnetic field (bottom) for the M2e3_E1_D8e25 (left), M2e3_noCRT (middle), and
M2e3_E1_D8e26 (right) simulations at $0.5t_{\rm ff}$. For the M2e3_noCRT
panel, a CR energy density of $\epsilon_{\rm cr}=1\rm eV/cm^{3}$ is shown,
corresponding to the assumed ionization rate of $\zeta\sim 1.7\times
10^{-17}\rm s^{-1}$. The gray dotted contours correspond to gas densities
$n=1,10,100\rm cm^{-3}$, also marked in Figure 11.

We note that in our simulations this ‘ring’ effect may be due to the
initialization of the CR energy and the varying cell size across the cloud
boundary. If the CRs were able to diffuse freely in this region then the
excess may be rapidly smoothed out; instead the enhancement persists even in
the M2e3_E1_D8e26 simulation. Still, the spatial correlation between the CRs
and the magnetic field deserves consideration due to the number of studies
that have looked at the effects of large-scale non-uniformity of the magnetic
field on CR propagation. Magnetic field lines converge in regions undergoing
gravitational collapse which results in a proportional increase in the number
density of CRs, a phenomenon known as “focusing”. On the other hand, the CRs’
pitch angles increase in response to the growing field and more particles are
reflected back, known as “mirroring.” A variety of studies have investigated
the competing effects of magnetic mirroring and focusing on CR propagation in
MCs (Kulsrud & Pearce, 1969; Skilling & Strong, 1976; Cesarsky & Volk, 1978;
Chandran, 2000; Padoan & Scalo, 2005; Silsbee et al., 2018) and in dense cores
(Padovani & Galli, 2011; Padovani et al., 2013; Fatuzzo & Adams, 2014). Which
effect dominates depends on the assumed magnetic geometry and details of the
CR microphysics.

Recent studies have suggested that although magnetic fields may partially
screen CRs from MCs, the effect is at most modest. Padovani & Galli (2011) and
Padovani et al. (2013) find that mirroring dominates over focusing for dense
cores but only causes a reduction in the CR energy density by a factor of 2-3
compared to the case where these effects are neglected. Silsbee et al. (2018)
found that the two effects cancel except in magnetic “pockets,” local minima
in the field strength, where mirroring again dominates by trapping CRs in the
pocket; however, even in these regions the reduction in the CRIR is less than
an order of magnitude. We conclude that although the CR enhancement along the
magnetic field gradient at the boundary may be physical rather than numerical,
the primary effect excluding CRs from the cloud interior is streaming
instability energy losses.

Past numerical studies have identified CR exclusion from MCs; however, most of
them have attributed the cause to effects other than streaming losses.
Skilling & Strong (1976) and Cesarsky & Volk (1978) found a drop in CR energy
density for CRs of $E<300\rm MeV$ due to ionization losses, thereby creating a
CR gradient that amplified the streaming instability at the boundaries of
clouds. Although we have a gradient in the CR energy density at the cloud
boundary due to the streaming instability, our single-bin approximation of CRs
approximates $\sim$ GeV CRs which do not experience the strong ionization
losses present in $\sim$ MeV CRs, so this effect is not present here.

Our findings are consistent with a few prior studies. Everett & Zweibel (2011)
considered a similar simulation configuration in which ionization losses are
also subdominant to streaming losses. They showed that even a small CR
pressure gradient amplifies the growth of the CR streaming instability.
However, they found only a small drop in CRs at the boundaries of clouds
rather than the almost complete exclusion we find. Bustard & Zweibel (2021)
considered more complex cloud geometries and a more advanced CRT model. Their
results show gradients in the CR energy density at the cloud boundary similar
to those in our simulations.

Most prior studies (Kulsrud & Pearce, 1969; Skilling & Strong, 1976; Cesarsky
& Volk, 1978; Chandran, 2000; Padoan & Scalo, 2005; Silsbee et al., 2018;
Padovani & Galli, 2011; Padovani et al., 2013) solve a set of simplified CR
propagation equations that model a steady-state, nearly isotropic CR
distribution, and they do not model the evolution of the gas and magnetic
field. Even more importantly, they make the approximation that CR pitch-angle
scattering due to the streaming instability is negligible and only model CR
streaming along field lines. This is a reasonable assumption on microphysical
timescales such as the CR gyro time or ion-neutral damping timescale, where
streaming instability energy losses are insignificant (Kulsrud & Pearce,
1969); however, on the longer timescales of MC evolution comparable to CR
transport timescales, CRs are expected to lose a large fraction of their
energy (Wentzel, 1971). Everett & Zweibel (2011) and Bustard & Zweibel (2021)
both model CRs interacting with diffuse clouds; thus, they use more advanced
moment-based CRT models and consider $\sim$ GeV CRs to more accurately capture
streaming transport and scattering. However, Everett & Zweibel (2011) do not
model the evolution of the cloud properties. Bustard & Zweibel (2021) provides
the most accurate point of comparison to our simulations; however, their two-
moment CRT implementation is only accurate for steady state, close-to-
isotropic cases. None of these studies considered the full effect of streaming
energy losses on the further cloud evolution, which we have shown is an
important component of CR physics that should be considered in future studies.

####  4.1.3 Are Streaming Instability Energy Losses Unavoidable in Molecular
Clouds?

The results of our simulations suggest that streaming instability energy
losses strongly limit the CR energy density in MCs. In this section, we
discuss whether this is an unavoidable outcome of CR physics in MCs.

We first consider the situation where the CRT is streaming dominated, which is
the case in all except the M2e3_E1_D8e26 simulation. In this limit, the
Eulerian equation for the CR energy density solved by many other studies
(e.g., Chan et al., 2019) converges to $\partial\epsilon_{\rm cr}/\partial
t=-\epsilon_{\rm cr}\nabla\cdot\mathbb{v_{\rm st}}-\gamma_{\rm cr}|v_{\rm
st}|\epsilon_{\rm cr}/\ell_{\rm cr}$. This equation limits to $\epsilon_{\rm
cr}\rightarrow 0$ if the streaming velocity increases going into the cloud
($\nabla\cdot\mathbb{v_{\rm st}}$ is positive) and/or $|\epsilon_{\rm
cr}\nabla\cdot\mathbb{v_{\rm st}}|<|\gamma_{\rm cr}v_{\rm st}\epsilon_{\rm
cr}/\ell_{\rm cr}|$. The only way for $\epsilon_{\rm cr}$ to stabilize at a
higher value is if $|\epsilon_{\rm cr}\nabla\cdot\mathbb{v_{\rm
st}}|>|\gamma_{\rm cr}v_{\rm st}\epsilon_{\rm cr}/\ell_{\rm cr}|$ and $v_{\rm
st}$ decreases going into the cloud. However, this is opposite of the expected
behavior, because the ionization fraction decreases in the cloud and the
streaming velocity increases (Figure 11). Consequently, streaming losses are
unavoidable if streaming is the dominant transport mechanism.

Next we consider the case where streaming and diffusive transport are
comparable as in our M2e3_E1_D8e26 simulation. This simulation models an
environment with a strong macroscopic diffusivity $D_{\parallel,\rm FLW}$
where turbulence advects the magnetic field lines at a rate faster than
streaming energy losses remove energy. However, this scenario was only
realized by ad-hoc setting the diffusion coefficient to a high value and is
not consistent with the degree of field-line wandering expected in mostly-
neutral gas unless $f_{\rm ion}\lesssim 10^{-9}$ (Sampson et al., 2022).
Figure 10 shows that an ionization fraction this low is rarely achieved in our
simulations, especially near the cloud boundary, and, furthermore, is not in
line with measured values in MCs, (Caselli et al., 1998). Moreover, this limit
could not be sustained because the CRs would stream with minimal energy
losses, which would increase $f_{\rm ion}$ and thereby decrease the ion Alfvén
velocity.

It is also impossible to achieve this limit with microphysical diffusion. As
discussed in Section 2.2.6, if the streaming instability is the only source of
microphysical scattering then achieving $v_{\rm diff}\gg v_{\rm Ai}$ would
require $\epsilon_{\rm cr}\ll 1\rm eV/cm^{3}$. Adding additional sources of CR
scattering from extrinsic turbulence would decrease $D_{\parallel,\rm micro}$;
in other words, assuming scattering due to self-generated turbulence alone
provides an upper limit to $D_{\parallel,\rm micro}$. Therefore, it is not
possible to reach the diffusion-dominated transport limit unless
$\epsilon_{\rm cr}\ll 1\rm eV/cm^{3}$ in the cloud from the start.

An alternative way to reduce streaming losses is to suppress the growth of the
streaming instability by limiting the drift velocity to less than the ion-
Alfvèn velocity ($v_{D}=v_{\rm st}+v_{\rm diff}<v_{\rm Ai}$). In steady state,
the growth rate of the streaming instability (Equation 15) is matched by the
ion neutral damping rate (Equation 14), which sets the diffusive component of
the velocity $v_{\rm diff}$. However, Equation 15 shows that when
$v_{D}<v_{\rm Ai}$, $\Gamma_{\rm SI}<0$ and the streaming instability will
never grow. In our simulations, we directly model CR streaming using $v_{\rm
st}=v_{\rm Ai}$ and include $v_{\rm diff}$ using the diffusion coefficient
$D_{\parallel}$, so by construction we cannot capture a situation where
$v_{D}<v_{\rm Ai}$. However, in principle the streaming velocity could be less
than the ion Alfvèn velocity. One way to achieve this is if $v_{\rm Ai}\gtrsim
c/3$ (where the factor of 3 comes from appropriately including the moments of
the CR distribution function). Using appropriate gas quantities, ($B\sim 10\mu
G,n\sim 100\rm cm^{-3},m_{\rm ion}\sim 29m_{p}$), $v_{\rm Ai}\approx 0.5f_{\rm
ion}^{-1/2}\rm km/s$ and so achieving $v_{\rm Ai}\gtrsim c/3$ would require
$f_{\rm ion}\lesssim 10^{-9}$. For the reasons discussed above, this limit
could not be sustained.

Another way to keep $v_{D}<v_{\rm Ai}$ is through an extrinsic source of
strong, symmetric scattering modes that scatter the CRs isotropically. The
full expression for the CR streaming speed is $v_{\rm st}=v_{\rm
Ai}(\bar{\nu}_{s}^{+}-\bar{\nu}_{s}^{-})/(\bar{\nu}_{s}^{+}+\bar{\nu}_{s}^{-})$,
where $\bar{\nu}_{s}^{+}$ and $\bar{\nu}_{s}^{-}$ are the scattering modes for
forward and backward propagating waves, respectively (Hopkins et al., 2021b).
Extrinsic sources of turbulence tend to produce symmetric scattering with
$\bar{\nu}_{s}^{+}\approx\bar{\nu}_{s}^{-}$. In this limit, $v_{\rm
st}\rightarrow 0$ so the $v_{\rm st}$ terms in Equations 4 and 5 vanish and
the diffusive acceleration term $\tilde{S}_{\rm sc}$ dominates, meaning the
CRs may gain energy as they propagate.

For this model to be viable, the turbulence must be strong enough such that
$v_{\rm diff}\sim D_{\parallel,\rm micro}/\ell_{\rm cr}\ll v_{\rm Ai}$ so the
streaming instability is not excited in the first place. Using appropriate
quantities for the diffuse ISM ($B\sim\mu G,n\sim 1\rm cm^{-3},m_{\rm ion}\sim
m_{p}$), $v_{\rm Ai}\approx 2\rm km/s$ $f_{\rm ion}^{-1/2}$, and so
$D_{\parallel}\lesssim 5\times 10^{23}(\ell_{\rm cr}/\rm pc$) $f_{\rm
ion}^{-1/2}\rm cm^{2}/s$. This value is much lower than the microphysical
diffusivity in the diffuse ionized ISM ($f_{\rm ion}\approx 1$), where
simulations estimate the diffusivity at a few $\times 10^{28}\rm cm^{2}/s$
(Trotta et al., 2011). Moreover, this scattering would have to be strong
enough to compete with ion-neutral damping in MCs. We tested this scenario by
including a standard model for extrinsic turbulence scattering in our
M2e3_E1_DVAR simulation and found that this limit was never reached
(otherwise, there would have been no streaming transport or losses in that
simulation).

We conclude that any scenario in which streaming instability energy losses are
avoided in MCs requires new physics. Either our understanding of the CR
streaming instability and self confinement is wrong, or there are sources of
extrinsic turbulence that are stronger than any currently known that can
produce the necessary scattering.

###  4.2 Comparison to Observations

####  4.2.1 Ionization Within Molecular Clouds

Measured values of the CRIR in the Milky Way vary by over an order of
magnitude. The value measured directly by the Voyager spacecraft, which has
the lowest amount of uncertainty, is $\zeta\sim 2\times 10^{-17}\rm s^{-1}$
(Stone et al., 2019). However, this only probes the CRIR in the Solar
neighborhood. Indirect measurements of the CRIR obtained from $\rm H_{3}^{+}$
observations of diffuse molecular gas suggest the CRIR varies over a few
orders of magnitude within the Milky Way, with an average value of a few
$\times 10^{-16}\rm s^{-1}$ (Indriolo & McCall, 2012; Indriolo et al., 2015;
Neufeld & Wolfire, 2017; Kálosi et al., 2023).

Although most of our simulations assume a CR background consistent with that
measured by Voyager, only the M2e3_E1_D8e26 simulation reached a CR level
throughout the cloud of $\zeta\sim 10^{-17}$. Our M2e3_E10_D8e25 simulation,
which adopted an external CR energy density 10 times that measured by Voyager,
had a larger range of CRIRs throughout the cloud than the other simulations
(top right panel of Figure 9). However, it still had a significantly lower CR
energy density inside the cloud due to CR exclusion. These results suggest
that the observed elevated ionization rates are unlikely to come from spatial
variations in the galactic CR distribution pervading the clouds. For low to
moderate values of the CR diffusion coefficient, the CRs lose a significant
amount of energy due to the streaming instability. For high values of the CR
diffusion coefficient, the CRs do not scatter enough to create significant
spatial variations.

Measurements of higher density gas tracers indicate that there is tentative
evidence that the CRIR declines in higher column density gas (e.g., Padovani,
2023). Ionization maps of dense MCs using $\rm H_{2}D^{+}$ show CRIRs lower
than those in diffuse clouds ranging from a few $\times 10^{-17}-10^{-16}\rm
s^{-1}$, in agreement with this conclusion (Sabatini et al., 2023). However,
they also show that within MCs the CRIR often spans orders of magnitude as
well. Observations of $\text{HCO}^{+}$ and $\text{N}_{2}\text{H}^{+}$ along
sightlines towards the intermediate-mass protostar OMC-2 FIR 4 indicate values
of $\approx 10^{-14}\rm s^{-1}$, a couple orders of magnitude higher than the
mean CRIR (Ceccarelli et al., 2014; Fontani et al., 2017). Meanwhile,
ionization maps of the Class 0 protostar B335 show large fractions of ionized
gas ($f_{\rm ion}\approx 1-8\times 10^{-6}$) and suggest the ionization
increases towards small envelope radii, reaching values of $\sim~{}10^{-14}\rm
s^{-1}$ a few hundred au from the central protostar (Cabedo et al., 2023).

Models used to derive the CRIR using dense gas tracers have significant
uncertainties and their accuracy varies with the source properties
(Shingledecker et al., 2016). Still, none of our simulations produced regions
approaching the elevated values inferred for B335. They show that streaming
instability energy losses result in a uniform CRIR throughout the cloud,
without enhancements in denser regions.

Our study lends support to the idea that there may be internal sources of CR
acceleration within MCs associated with the protostellar shocks (Padovani et
al., 2016; Gaches & Offner, 2018; Fitz Axen et al., 2021). Krumholz et al.
(2023) estimate that protostellar sources contribute $\sim 1/4-1/3$ as much
energy to the galactic CR energy budget as stellar winds. Low mass stars form
earlier and are more numerous than high mass stars, and so would accelerate
CRs earlier in the cloud collapse than stellar winds. These CRs may explain
the observed high CRIRs, especially near protostars, as well as the overall
elevated ionization rates inferred in diffuse clouds. However, the degree of
ionization inside clouds is sensitive to CR physics and the CR diffusion
coefficient, which is highly uncertain. Thus, ionization rate measurements and
comparisons based on synthetic observations can provide helpful constraints
for models of CR acceleration and propagation within MCs.

####  4.2.2 Star Formation in a High CR environment

Many regions in the Milky Way and other galaxies have a CRIR that is higher
than the Voyager value of $\sim 10^{-17}\rm s^{-1}$. In particular, recent
observations show evidence of elevated CRIRs towards the center of galaxies.
Although the Galactic average is $\approx 2\times 10^{-16}\rm s^{-1}$
(Indriolo et al., 2015), the CRIR measured in the central parsec of the Milky
Way is $\approx 2$ orders of magnitude higher than that in the clouds outside
the galactic center (Goto, 2014; Indriolo et al., 2015; Rivilla et al., 2022).
Extragalactic measurements towards the central molecular zone of the starburst
galaxy NGC 253 suggest even higher values. Holdship et al. (2022) studied
several star-forming regions there and found CRIR values between $1-80\times
10^{-14}\rm s^{-1}$. Behrens et al. (2022) found that the central GMCs of NGC
253 had $\zeta\approx 10^{-13}-10^{-11}\rm s^{-1}$, greater than four orders
of magnitude higher than the Milky Way value.

We note that the accuracy of these measurements are in tension with
theoretical predictions. Krumholz et al. (2023) argue that these CRIR values
are inconsistent with the measured star formation rates, which determine the
amount of energy available for accelerating CRs and that the true CRIRs are
significantly lower, though still likely a factor of 10-100 greater than the
Voyager value. They are also inconsistent with other observational and
theoretical expectations, which predict an anticorrelation between CRIR and
gas density (Indriolo & McCall, 2012; Neufeld & Wolfire, 2017; Padovani,
2023). However, regions of the ISM near SNRs may have elevated CR levels, as
SNRs are thought to be the primary accelerators of Galactic CRs (Aharonian,
2013). Indeed, measurements of CRIRs targeting clouds near SNRs show elevated
rates of $\approx 100$ times the Galactic average (Ceccarelli et al., 2011;
Vaupré et al., 2014).

We found that star formation occurs more rapidly and the SFE is a few percent
higher in our M2e3_E10_D8e25 simulation, which raises the CR background level
by only an order of magnitude. Thus, we might expect larger differences in
star formation in even more extreme CR environments, such as near SNRs.
However, the details of the effect are highly sensitive to the assumed value
of the diffusion coefficient. Raising the diffusion coefficient by an order of
magnitude caused a large difference in the CR propagation within the cloud,
diminishing energy losses from the streaming instability and exclusion by
magnetic fields. Meanwhile, decreasing the diffusion coefficient by an order
of magnitude caused even more rapid CR energy losses. These results suggest a
higher diffusion coefficient and a higher CR environment acting in concert
could produce more extreme variations in star formation outcomes. Further work
is needed to investigate this parameter space.

###  4.3 Simulation Caveats

####  4.3.1 Energy-Independent Treatment of CRs

Our simulations make a number of approximations in the treatment of CR
physics. We use a single-bin model for the CR fluid that represents the total
CR energy rather than evolving a full multi-bin energy spectrum. The treatment
of the CR velocities and the energy loss function assume $\sim$ GeV CRs.
Simulations modeling a full multi-bin CR energy spectrum that extends down to
$\sim$ MeV and up to $\sim$ TeV energies would likely predict different CR
behavior within the dense gas. For example, $\sim$ MeV CRs more efficiently
ionize gas than $\sim$ GeV CRs. Collisional energy losses are anticorrelated
with gas density (see Equation 9), so multi-bin CRT simulations might better
reproduce the apparent density-dependent CRIR seen in observations (Padovani,
2023). Modeling CRs with $\gtrsim$ GeV energies would likely have a minimal
impact on our results, because these CRs are less numerous than CRs of lower
energies (Stone et al., 2019). However, these CRs produce $\gamma$-rays, so
comparisons to observed $\gamma$-ray emission could lead to additional
constraints on CRT in MCs.

We note that modeling a full CR spectrum would not change the conclusions
presented in 4.1.3 because they apply at all CR energies as long as the
streaming instability dominates the CR scattering. CRs of all energies will
stream at the ion Alfvèn speed and experience streaming energy losses. Thus,
the cloud will still develop a low CR energy density eventually, even if there
is measurable ionization and $\gamma$ ray emission before the CRs lose their
energy.

####  4.3.2 Reduced Speed of Light Approximation

We use a reduced-speed-of-light (RSOL) of 300 km/s for our simulations to make
them computationally tractable. This is higher than the value used for other
STARFORGE calculations of $\tilde{c}=30\rm km/s$ (Grudić et al., 2021, 2022)
or $\tilde{c}=90\rm km/s$ (Grudić et al., 2023), which has been thoroughly
tested for non-CRT simulations in which the RSOL affects the radiative
transfer calculation. The choice of RSOL should not have a significant impact
on our results as long as it is higher than the other relevant velocities in
the simulations such as the CR streaming speed $v_{\rm st}$ and the diffusion
speed set by the diffusivity $D_{\parallel}$. Otherwise, the CR transport and
loss processes may be artificially slowed down (Hopkins et al., 2020). As
discussed in Section 4.1.1 and shown in Figure 12, the streaming velocity is
an order of magnitude higher than 300 km/s in many regions of the cloud. In
order to justify our choice of RSOL, we run a suite of low-resolution
simulations that vary the RSOL. We find that using a higher RSOL increases the
speed of CR transport and rate of CR energy losses in the cloud but does not
significantly effect the star formation properties. We show these results and
discuss the implications further in Appendix A.

####  4.3.3 Sub-Grid Models and Resolution

Our simulations predict an IMF with a high-mass slope $\alpha=-1$ consistent
with results from previous STARFORGE calculations (Figure 8). As discussed in
Grudić et al. (2022), limitations due to the feedback models and resolution
may produce a shallower IMF than the commonly assumed value of $\alpha=-1.35$.

The STARFORGE protostellar jet feedback model, discussed in Section 2.1.1,
uses fixed parameters $f_{w}=f_{K}=0.3$, which parameterize the jet velocities
but are observationally uncertain (Cunningham et al., 2011). Variations in
these parameters which increase the jet velocity may regulate massive SF
(Guszejnov et al., 2021; Rosen et al., 2021).

We are limited in resolution in both length scale and mass (Grudić et al.,
2022; Guszejnov et al., 2022). Size scales smaller than the Jeans length are
not well resolved in the simulations, so protostellar disks are unresolved.
Disk fragmentation may affect the stellar multiplicity and mass ratio
distribution of binaries (Offner et al., 2023), so fully resolving disks may
steepen the IMF.

####  4.3.4 Calculation of the Streaming Speed

To calculate the ion Alfvèn velocity, GIZMO assumes $m_{\rm ion}=m_{p}$, which
is equivalent to setting $\chi=f_{\rm ion}$. This is not a correct assumption
for MCs, where the dominant charge carrier is $\rm HCO^{+}$ ($m_{\rm
ion}=29m_{p}$), and so overestimates $v_{\rm Ai}$ by a factor of
$\sim\sqrt{29}\approx 5$. Additionally, we do not account for the effect of
dust damping on the CR streaming speed. Dust grains can damp the amplitude of
the Alfvèn waves from the streaming instability, reducing the scattering rate
and increasing the streaming speed (Squire et al., 2021).

Since both streaming transport and losses scale with $v_{\rm Ai}$, these
caveats do not change the major conclusions presented throughout Section 4
emphasizing the importance of streaming losses when the CR transport is
dominated by streaming. However, they may modify when CRs are in the streaming
versus diffusion dominated regimes (Section 4.1.1). These uncertainties are
comparable to others in the chemistry module, and for given values of $v_{\rm
st}$ and $v_{\rm diff}$ our conclusions are robust.

##  5 Conclusions

In this work, we present numerical simulations of the collapse of MCs
including CRT as part of the STARFORGE project. These simulations follow the
formation of individual stars and include all major stellar feedback
mechanisms including protostellar jets, radiation, and stellar winds. We
expand on previous STARFORGE simulations by explicitly modeling CRT using a
single-energy bin model of CRs, which are evolved as a fluid using an
M1-moment based method. We explore variations in the initial CR energy density
and CR diffusion model. Our main findings are as follows:

1\. In a fiducial CR environment where $\epsilon_{\rm cr}=1\rm eV/cm^{3}$,
including CRT does not produce a strong difference in the dynamics of the gas
collapse compared to simulations without CRT (Figures 1 and 3). This is true
regardless of the model used for the diffusion coefficient (Figure 5).

2\. In a high CR environment where $\epsilon_{\rm cr}=10\rm eV/cm^{3}$, the
enhanced CR pressure outside the cloud causes faster collapse, earlier star
formation and a higher final SFE (Figure 5). The SFE in the high CR
environment increases earlier due to the earlier cloud collapse but progresses
at the same rate as clouds embedded in a typical CR environment.

3\. For low to moderate values of the CR diffusion coefficient, streaming
instability energy losses strongly attenuate the CRs (Figure 2) . However,
when the diffusion coefficient is high enough ($\gtrsim 10^{26}\rm cm^{2}/s$)
the CRs lose less energy and distribute uniformly throughout the cloud at a
level similar to that outside the cloud (Figure 5).

4\. In a CR environment where $\epsilon_{\rm cr}=1\rm eV/cm^{3}$ the evolution
of the SFE is insensitive to the CR treatment (Figure 7). The final SFE is
slightly higher for all of the CRT clouds than for the cloud with no explicit
CRT. The no-CRT simulation has a $\rm SFE\approx 7$ %, while all of the CRT
simulations have $\rm SFE=8-10$ %. The high-CR cloud has a $\rm SFE>12$ %.

5\. Variations in the CR diffusion coefficient do not have a statistically
significant impact on the IMF of the clouds (Figure 8). All simulations in an
environment with $\epsilon_{\rm cr}=1\rm eV/cm^{3}$ have a final stellar mass
distribution consistent with previous STARFORGE calculations, with a peak at
$\sim 0.2-0.3M_{\odot}$ and a high mass slope of $\alpha\approx-1$. However,
one of the CRT simulations (M2e3_E1_D8e25) has a higher median mass, possibly
due to stellar wind CR feedback that changes the gas conditions inside the
cloud. In a high-CR environment with $\epsilon_{\rm cr}=10\rm eV/cm^{3}$ the
cloud forms both more low mass ($M<0.1M_{\odot}$) and high mass
($M>1M_{\odot}$) stars.

6\. In a fiducial CR environment the CRIR is roughly uniform throughout the
cloud (Figure 9). For low to moderate values of the CR diffusion coefficient
it is up to a few orders of magnitude less than outside the cloud, with a
median value $\zeta\lesssim 2\times 10^{-19}\rm s^{-1}$ at $1-1.5t_{\rm ff}$.
In a high-CR environment simulation the CRIR in the cloud is higher than in a
fiducial environment ($\zeta\lesssim 6\times 10^{-18}\rm s^{-1}$) and there is
a larger spread of values, but is still less than the commonly assumed value
of $\sim 10^{-17}\rm s^{-1}$ throughout most of the cloud. Overall, the CR
background does not provide a significant contribution to the ionization of
the cloud compared to radiative feedback once star formation is underway
(Figure 10).

Further numerical studies exploring the impact of a high CR diffusion
coefficient in high-CR environments may provide further insights and
constraints on CRT. Furthermore, numerical simulations including locally
injected CRs from protostellar sources may explain the elevated CRIRs and
variations within clouds observed from various sources.

This research is part of the Frontera computing project at the Texas Advanced
Computing Center, and used computing award AST21002. Frontera is made possible
by National Science Foundation award OAC-1818253. This research was supported
in part by NASA ATP grant 80NSSC20K0507. This material is based upon work
supported by the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research, Department of Energy Computational
Science Graduate Fellowship under Award No. DE-SC0021110. Support for PFH was
provided by NSF Research Grants 1911233, 20009234, 2108318, NSF CAREER grant
1455342, NASA grants 80NSSC18K0562, HST-AR-15800. MRK acknowledges support
from the Australian Research Council with the Discovery Projects and Laureate
Fellowship schemes, awards DP 230101055 and FL220100020. Support for MYG was
provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51479
awarded by the Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., for NASA, under
contract NAS5-26555.

##  Appendix A Choice of the Reduced Speed of Light

We use a RSOL of $\tilde{c}=300$ km/s for both the CRs and radiation in our
simulations. This makes the simulations computationally feasible by limiting
the timestep criterion to $\Delta t\leq C_{\rm cour}\Delta x/\tilde{c}$ in all
regions of the grid. The choice of $\tilde{c}$ should be irrelevant if it is
much faster than the other velocities in the simulation, as the CR propagation
equations converge to the same solutions in steady-state independent of
$\tilde{c}$. However, the top panel of Figure 11 shows that in many regions
during the cloud evolution the CR streaming velocity is larger than 300 km/s.
Consequently, all transport and loss timescales except the advection timescale
are slowed down by a factor of $\Gamma_{\rm RSOL}=c/\tilde{c}=1000$ (Hopkins
et al., 2020), which is used to scale the timescales plotted in the bottom
panel of Figure 11.

It is not feasible to run the simulations at full resolution with
$\tilde{c}\sim 10^{3}-10^{4}\rm km/s$ in a reasonable amount of wall clock
time. Instead, we test the impact of our choice of $\tilde{c}$ by performing
six low resolution simulations ($N_{\rm gas}=2\times 10^{5}$) that vary the
value of $\tilde{c}$. Three of these start with the initial CR background
equivalent to that used in the M2e3_E1_D8e25 simulation ($\epsilon_{\rm
cr},_{\rm med}$ = 1 $\rm eV/cm^{3}$) and three start with the initial CR
energy density used in the M2e3_E10_D8e25 simulation ($\epsilon_{\rm cr},_{\rm
med}$ = 10 $\rm eV/cm^{3}$). For each initial CR configuration, we adopt
values $\tilde{c}=75,300$, and 1200 km/s. All six tests use the fiducial
diffusion coefficient $D_{\parallel}=8.33\times 10^{25}\rm cm^{2}/s$. The six
different runs are listed in Table 4.

Figure 13 shows slice plots of the gas number density (top) and CR energy
density (bottom) for the M2e3_E1_D8e25 setup at $t=t_{\rm ff}$ using
$\tilde{c}=75\rm km/s$ (left), $\tilde{c}=300\rm km/s$ (middle), and
$\tilde{c}=1200\rm km/s$ (right). Increasing $\tilde{c}$ causes faster CR
energy losses because the streaming loss timescale $\tau_{\rm str,loss}$
scales inversely with $\tilde{c}$. The biggest discrepancy occurs near the
cloud boundaries where streaming energy losses are significant, as discussed
in Section 4.1.2. The top left panel of Figure 14 shows the evolution of the
median CR energy density as a function of time for various gas densities for
the three M2e3_E1_D8e25 setup tests. This plot shows that by $1t_{\rm ff}$ the
$\tilde{c}=75\rm km/s$ run has $\epsilon_{\rm cr}\gtrsim 10^{-2}\rm eV/cm^{3}$
in the cloud, while the $\tilde{c}=1200\rm km/s$ run has $\epsilon_{\rm
cr}<10^{-3}\rm eV/cm^{3}$. The variation between different choices of
$\tilde{c}$ is slightly more than the factor of four difference in
$\Gamma_{\rm RSOL}$, which indicates that the choice of $\tilde{c}$ has an
impact on the CR energy evolution even greater than the difference expected
from timescale analysis.

Another potential consequence of using a low value of $\tilde{c}$ is that CRs
can become artificially trapped in dense gas (Hopkins et al., 2020). Since the
calculations adopt a Lagrangian frame of reference that moves with the gas,
the advection timescale $\tau_{\rm adv}$ is independent of $\tilde{c}$, but
the streaming and diffusion timescales $\tau_{\rm str}$ and $\tau_{\rm diff}$
are not. Using a low value of $\tilde{c}$ can cause CR advection to dominate
and not allow the CRs to escape through diffusion and streaming. Figure 14
shows that all three tests have slightly elevated values in gas densities
$n>10^{4}\rm cm^{-3}$ (black lines) at early times but reach a uniform CR
energy density in the cloud at slightly earlier times for a higher value of
$\tilde{c}$. Still, this difference is minor compared to the discrepancy
between the different tests for the same density range, indicating that the
enhancement caused by artificial CR trapping is subdominant compared to the
differences caused by streaming energy losses. However, this effect may
produce more significant differences in the evolution for smaller clouds.

Figure 13: Gas number density (top) and CR energy density (bottom) at $1t_{\rm
ff}$ for the runs with a reduced-speed-of-light $\tilde{c}=75$ km/s (left),
$\tilde{c}=300$ km/s (middle), and $\tilde{c}=1200$ km/s (right), starting
with the M2e3_E1_D8e25 CR distribution.

In addition to these differences in the CR transport, there is a correlation
between the value of $\tilde{c}$ used and the stellar mass distribution
formed. Figure 14 shows that the M2e3_E1_D8e25_RSOL1200 simulation does not
have a noticeable increase in the median CR energy density between $\sim
1.5-2.5t_{\rm ff}$ like the other M2e3_E1_D8e25 setup tests. The stars formed
in that test simulation did not have stellar winds powerful enough to inject a
significant amount of CR energy. To illustrate this point, the bottom row of
Figure 14 shows a plot of the evolution of the CR energy injected by stellar
winds (left) and the most massive star in the simulation (right) for our six
test runs. For a given initial CR distribution, a higher value of $\tilde{c}$
produces more rapid initial accretion onto a few massive stars, causing the
injected CR energy to rise as the stellar winds become more powerful due to
the more massive stars. However, this initial accretion phase is cut off
earlier, capping the stellar wind energy. In the M2e3_E1_D8e25_RSOL1200
simulation, no stars form that are massive enough to have stellar winds that
noticeably increase the median CR energy density in the cloud.

Despite these caveats, Figure 13 shows that the evolution of the gas density
is similar for all choices of $\tilde{c}$. This is partially because the
M2e3_E1_D8e25 setup is initialized with $\epsilon_{\rm cr}=1\rm eV/cm^{3}$,
which is not high enough to strongly influence the gas dynamics even with
minimal energy losses. Additionally, the effects described on the CRT are not
as drastic as they initially seem because at the relevant length scales our
choice of $\tilde{c}$ is large enough. Figure 11 shows the timescales for CRs
to propagate all the way through the cloud ($L=3$ pc) through streaming and
diffusion alone are $\approx 1-10$ Myr, comparable to the advection timescale
(which is a good approximation of the dynamical timescale). However, since the
advection timescale is independent of $\tilde{c}$, the gas collapse transports
CRs through the cloud and helps offset the slow transport timescales. The CR
energy density will reach a relatively smooth distribution on smaller
timescales than shown in Figure 11 because the simulation resolution is much
smaller than the cloud radius ($\Delta x_{J}\sim 36\rm au$), and CR gradients
at these scales will be wiped out quickly. Figures 2 and 11 show that by
$0.5-0.75t_{\rm ff}$ the CR energy density at $n\gtrsim 500\rm cm^{-3}$ is
roughly independent of gas density, indicating that CRs are indeed able to
propagate out of the densest gas.

Furthermore, despite the differences in the stellar accretion, differences in
the evolution of the CR distribution in the cloud do not change the SFE until
near the end of the cloud evolution. Figure 14, top right, shows the evolution
of the SFE for the six different runs. As in Section 3.2, we assume SFE
$\approx 7.42\pm 1.77$ as the $1\sigma$ region, although the exact numbers
should not be compared to our main simulations with $N_{\rm gas}=2\times
10^{6}$ as resolution causes variations in the results as well. For both the
M2e3_E1_D8e25 variations and M2e3_E10_D8e25 variations, the evolution of the
SFE begins similarly for the different values of $\tilde{c}$. However, the
time at which star formation slows and the cloud disperses occurs later for a
lower choice of $\tilde{c}$, raising the final value of the SFE.

The final values of the SFE for the six runs are shown in Table 4. If we
assume SFE $\approx 7.42\pm 1.77$, then both the M2e3_E1_D8e25_RSOL75 and
M2e3_E10_D8e25_RSOL75 tests have values of the SFE above the $1\sigma$ region.
The difference between the M2e3_E1_D8e25_RSOL75 and the M2e3_E1_D8e25_RSOL300
tests is $\approx 2\%$, while between the M2e3_E10_D8e25_RSOL75 and
M2e3_E10_D8e25_RSOL300 tests it is $\approx 3\%$. The RSOL1200 tests did not
have as large of a variation in the final SFE from the RSOL300 tests, and only
had a difference in the SFE of $\sim 0.35\%$ and $\sim 1.2\%$ for the
M2e3_E1_D8e25_RSOL1200 and M2e3_E10_D8e25_RSOL1200 tests respectively. Since
these differences are less than $1\sigma$ and we focus this study on the SFE
rather than the IMF, we conclude that using $\tilde{c}=300$ km/s is
sufficiently resolved for our simulations. However, the trends in the SFE and
stellar mass distribution suggest that future numerical simulations of CR
propagation in MCs, particularly studies exploring a high CR environment,
should carefully consider the choice of $\tilde{c}$.

Simulation Name | RSOL [km/s] | SFE (%)
---|---|---
M2e3_E1_D8e25_RSOL75 | 75 | 9.3
M2e3_E1_D8e25_RSOL300 | 300 | 7.46
M2e3_E1_D8e25_RSOL1200 | 1200 | 7.09
M2e3_E10_D8e25_RSOL75 | 75 | 10.9
M2e3_E10_D8e25_RSOL300 | 300 | 7.8
M2e3_E10_D8e25_RSOL1200 | 1200 | 6.62
Table 4: Simulation name, RSOL, and final values for the SFE for six different
tests varying $\tilde{c}$. All simulations use $N_{\rm gas}=2\times 10^{5}$
and $D_{\parallel}=8.33\times 10^{25}\rm cm^{2}/s$.

Figure 14: Top Left: Median CR energy density for gas in the indicated density
range or lower limit for the M2e3_E1_D8e25_RSOL75 (dashdot),
M2e3_E1_D8e25_RSOL300 (solid), and M2e3_E1_D8e25_RSOL1200 (dashed)
simulations. Top Right: Evolution of the SFE for the RSOL tests. Shaded
regions are the $1\sigma$ (dark gray) and $2\sigma$ (light gray) final SFE
values obtained from a set of 100 2,000 M⊙ cloud simulations without CRT that
vary the initial turbulent seed. We note that both the tests with the fiducial
RSOL of $\tilde{c}=300\rm km/s$ do not evolve exactly identically to the full
resolution analogs due to numerical differences. Bottom left: Evolution of the
total CR energy injected by stellar winds for the RSOL tests. Bottom right:
Evolution of the maximum stellar mass for the RSOL tests.

##  Appendix B Tests With Varying Physics

Figure 11 shows that for our M2e3_E1_D8e25 simulation both CR transport and
energy losses are dominated by streaming, especially at the cloud boundary. As
we discuss in Section 4.1.2, these effects cause CRs that originated in the
cloud to lose most of their energy and prevent CRs outside the cloud from
propagating in without significant energy losses (Figure 12). In order to
confirm these conclusions, we run low resolution simulations ($N_{\rm
gas}=2\times 10^{5}$) out to $\sim 1.5t_{\rm ff}$ where we vary the CR
transport and loss physics. The first one has identical transport and loss
physics to our high resolution simulations. In the second we disable all
heating/cooling interactions with the gas, and consequently all collisional
losses (Equations 9 and 10). In the final test simulation we disable both
streaming transport and energy losses (all terms including $v_{\rm st}$ in
Equations 4 and 5). All three of these simulations start with the initial CR
background and diffusion coefficient used in the M2e3_E1_D8e25 simulation
($\epsilon_{\rm cr},_{\rm med}$ = 1 $\rm eV/cm^{3}$ and