key: cord-0335703-3kf2404r
authors: Gaches, Brandt A. L.; Walch, Stefanie; Lazarian, Alex
title: CRAFT (Cosmic Ray Acceleration From Turbulence) in Molecular Clouds
date: 2021-08-06
journal: nan
DOI: 10.3847/2041-8213/ac1b2f
sha: 82baf052fcc12d2171dc53dd1b8053993ddbd252
doc_id: 335703
cord_uid: 3kf2404r

Low-energy cosmic rays, in particular protons with energies below 1 GeV, are significant drivers of the thermochemistry of molecular clouds. However, these cosmic rays are also greatly impacted by energy losses and magnetic field transport effects in molecular gas. Explaining cosmic ray ionization rates of $10^{-16}$ s$^{-1}$ or greater in dense gas requires either a high external cosmic ray flux, or local sources of MeV-GeV cosmic ray protons. We present a new local source of low-energy cosmic rays in molecular clouds: first order Fermi-acceleration of protons in regions undergoing turbulent reconnection in molecular clouds. We show from energetic-based arguments there is sufficient energy within the magneto-hydrodynamic turbulent cascade to produce ionization rates compatible with inferred ionization rates in molecular clouds. As turbulent reconnection is a volume-filling process, the proposed mechanism can produce a near-homogeneous distribution of low-energy cosmic rays within molecular clouds.

Molecular clouds are immersed in a bath of cosmic rays (CRs), i.e. energetic charged particles that are accelerating and propagating in our galaxy (Schlickeiser 2002) . Low-energy particles, particularly protons with energies between 1 MeV to 1 GeV, influence the thermochemistry of molecular gas in regions which are well-shielded from ultraviolet radiation (see re-views by Dalgarno 2006; Padovani et al. 2020) . The chemistry of cold molecular gas is regulated through ion-neutral reactions. Ion-neutral chemistry is largely initiated following the ionization of H 2 and subsequent production of H + 3 H 2 + H + 2 → H + 3 + H.

Deuterium chemistry is also regulated in a similar manner through the production of H 2 D + from HD. The CR ionization rate (CRIR), ζ, is inferred primarily through molecular line features from ions, such as H + 3 , OH + and H n O + absorption and emission from species such as HCO + , DCO + and N 2 H + . These observations have shown that the CRIR in diffuse molecular gas spans 10 −16 < ζ < 10 −15 s −1 (e.g. Indriolo & McCall 2012; Indriolo et al. 2015; Neufeld & Wolfire 2017) . In the dense gas, observations infer the CRIR in a wider range from 10 −17 < ζ < 10 −15 s −1 (e.g. Caselli et al. 1998; Favre et al. 2017; Barger & Garrod 2020) . Recently, observations and astrochemical models have inferred an CRIR towards the protocluster OMC-2 FIR 4 of approximately ζ ≈ 10 −14 s −1 (Ceccarelli et al. 2014; Favre et al. 2018) .

One-dimensional models of transport through molecular gas conflict with the heightened ionization rates inferred in shielded, dense gas. These models ubiquitously predict a declining CRIR with column density, ζ(N ), whether due to energy-losses 1 or magnetic-field effects 2 (Padovani et al. 2009; Morlino & Gabici 2015; Schlickeiser et al. 2016; Ivlev et al. 2018; Phan et al. 2018; Silsbee & Ivlev 2019; Fujita et al. 2021) . A dense gas CRIR similar to that inferred in diffuse gas may necessitate a source of low-energy CRs within the gas. Some of these sources have already been posited, such as protostellar jets (Padovani et al. 2016) , protostellar accretion shocks (Padovani et al. 2016; Gaches & Offner 2018) , HII regions Padovani et al. 2019 ) and embedded stellar winds (Yang & Wang 2020) . These are localized sources of CRs acceleration that are expected to induce rather inhomogeneous distribution of CRs within molecular clouds.

We propose a new source of low-energy CRs in dense magnetized molecular clouds: particles accelerated within zones of turbulent reconnection. Turbulence is known to be part and parcel of the molecular cloud dynamics (McKee & 1 The dominant energy losses for 1 MeV -1 GeV protons are due to ionizing atomic and molecular material and the production of pions. 2 Magnetic-field effects such as screening, mirroring and streaming instabilities.

Ostriker 2007) and it is known to be accompanied by fast turbulent reconnection (Lazarian & Vishniac 1999; Eyink et al. 2011 Eyink et al. , 2013 . The latter process is known to induce the acceleration of energetic particles (de Gouveia dal Pino & Lazarian 2005a; Kowal et al. 2011; Lazarian et al. 2020) . We explain in Section 2 our simplified model of CR production in reconnection zones in magneto-hydrodynamic turbulence. Section 3 presents the results of these calculations and discusses the broad implications of the mechanism.

We assume a spherical, magnetized, hierarchical cloud with magneto-hydrodynamic (MHD) turbulence. MHD turbulence is injected at some large scale, L, with a velocity dispersion, σ 0 , and cascades through the cloud. Figure 1 gives a schematic of the proposed mechanism. At a given length scale, , the turbulent linewidth is given by the linewidth-size relation (Larson 1981; McKee & Ostriker 2007; Heyer & Dame 2015) σ

where we use β = 0.5. The density of the gas is calculated by assuming the gas can be prescribed by the virial parameter, α V , defined by

where G is the gravitational constant. The density is thus

(4)

Finally, the magnetic field is calculated using the empirical fit from Crutcher & Kemball (2019)

where n = ρ/(2.33×m H ) is the number density, B 0 = 10 µG, n 0 = 300 cm −3 and κ = 0.65. The energy within the MHD turbulence cascade depends on the Alfvén Mach number,

In the ideal-MHD case, the dissipiation rate of the specific energy per unit mass is given by (Lazarian et al. 2020 

The dissipation rate of energy per unit volume is then ε = ρ . A fraction of this energy, f CR goes into CR acceleration, such that

where we take f CR = 0.01 as an estimated lower limit of the acceleration efficiency. We assume particles are accelerated within the turbulent reconnection regions via a first-order Fermi process (de Gouveia dal Pino & Lazarian 2005b). Following Khiali et al. (2015) , the CRs are isotropically injected with an exponentially suppressed power law

where we take γ = 2 and E 0 = 10 GeV 3 . Increasing E 0 negligibly impacts our main results, due to the weak dependence of the CRIR on super-GeV CRs. Further, changing γ between 2 and 3/2 produces no qualitative changes in the results, nor quantitative variations over an order of magnitude. The normalization factor, Q 0 , is calculated by assuming

where E min = 13.6 eV and E max = 100 GeV. These bounds have a minor impact on the over-all results of the work. Determining the injection and maximum energies requires particlein-cell calculations of the CR acceleration and injection within molecular cloud reconnection zones. However, even if energy losses are ignored, the necessary acceleration timescale from E max = eBv 2 A δt to accelerate protons up to 100 GeV exceeds molecular cloud lifetimes for much of the parameter space.

The CR proton spectrum from the reconnection zones is a balance of injection and energylosses. The steady-state energy-loss solution (Longair 2011) for the number density of protons within the reconnection region, N p (E), is

where dE dt is calculated using a prescribed loss function, L(E)

and M s is the sonic Mach number, M s = σ(l)/c s , c s = k b T /µm H , µ is the mean molecular weight and T = 10 K, and v CR (E) is the relativistic velocity of the CR. We utilize the loss function given in Padovani et al. (2009) . Turbulent reconnection is an essential part of the turbulent cascade (Lazarian & Vishniac 1999 ) and a volume-filling process. This induces in CR acceleration and we model the resulting CR number density accelerated by turbulent reconnection at length scale by assuming the CRs diffuse from the reconnection zones and undergo energy losses. We assume an energy- 

using δ = 0.5 and different values of D 0 . The diffusion length scale is defined by

4 The process of CRs diffusion in MHD turbulence is pretty complicated with different components of MHD modes acting very differently on CRs (see Yan & Lazarian 2004) and its effects for the diffusion parallel and perpendicular to the mean magnetic field is also not trivial (see Lazarian & Yan 2014) . However, for the same of simplicity, in this paper we adopt the simplest possible assumptions about the diffusion. This assumption is further justified by Lazarian & Xu (2021) which presented a new non-resonant scattering process in turbulent magnetized media.

The energy-loss scale, or the range, R(E) is given by the stopping column, (Padovani et al. 2009 )

We then define a transport length scale, T , −1

Finally, we use the volume-filling fraction of the reconnection zones

where V rec is the volume of a sheet-like reconnection zone,

l A = M −3 A for super-and trans-Alfvénic turbulence and l A = M 2 A for sub-Alfvénic turbulence (Lazarian et al. 2020) . Following Lazarian & Vishniac (1999) the reconnection zone width is

The volume of a region of radius, is V l = 4 3 π 3 . The number density of transported CRs is

Either number density, N p (E) or N T (E), can be converted to a flux by

The resulting CRIR due to turbulent reconnection at length scale, , is

where σ(E) is the total ionization cross section. We use the empirical fit from Rudd et al. (1985) .

For the following results, our canonical cloud is virialized, α V = 1, and turbulence is injected at a scale of 1 pc with a turbulent linewidth of 1 km s −1 . This results in the cloud primarily being sub-Alfvénic (see e.g. Crutcher & Kemball 2019). Figure 2 shows the steady-state flux, j p (E), and the transported flux, j T (E), as a function of scale and CR energy, E. Both Q(E) and j p (E) are weakly dependent on the length scale, . The CRIR associated with j p (E) are of order ζ ≈ 7 × 10 −10 − 10 −9 s −1 (not shown in the figure) . These CRIRs are too high to be physical, and highlight the necessity of treating the diffusion of the CRs throughout the rest of the cloud structures.

The final CR spectrum, j T (E), shows significant variation with both D 0 and . The flux at low energies is dramatically decreased due to energy losses from ionizations, Coulomb interactions and pion production while at high energies the flux is suppressed by diffusion. Our model predicts reconnection will seed the cloud with protons of energies E ≈ 10 6 − 10 10 eV. The resulting spectrum is greatly sensitive to the diffusion coefficient, D 0 . For diffusion coefficients D 0 = 10 29 and 3 × 10 28 cm 2 s −1 , the resulting spectrum is relatively flat. However, for lower values of D 0 , the spectrum is only flat for reconnection driven by the smallest scales of turbulence. We find a significant change in behavior between high D 0 and low D 0 values. For D 0 = 3 × 10 28 and 10 29 cm 2 s −1 , we find the flux typically increases with the turbulence driving length scale, . For D 0 = 10 28 and 3 × 10 27 cm 2 s −1 the flux decreases with increasing length scales. Conversely, if the particles travel ballistically, the transport length T ≈ R and the CRIR increases dramatically to ζ ≈ 10 −14 s −1 for the fiducial model. This CRIR is far outside the observed range within the Milky Way, except in sightlines towards the galactic center (Indriolo et al. 2015) . Therefore, in the framework of CR acceleration by magnetic reconnection, low energy CRs must travel diffusely or the acceleration efficiency must be f CR 0.01 to be consistent with observations. Figure 3 shows the resulting CRIR, ζ( ) as a function of length scale, , for different values of D 0 . We highlight in the blue box the range of inferred values of ζ for molecular clouds in the Milky Way. We find that for values of D 0 > 10 28 cm 2 s −1 , our model is able to produce (or overproduce) the CRIR inferred in molecular clouds. For low values of D 0 , the CRIR is below even the effective minimum ionization rate ≈ 10 −19 in molecular clouds due to radioactive nuclei decay (Adams et al. 2014) . Table 1 shows the power-spectrum averaged CRIR, defined as Note-The power-spectrum averaged cloud cosmic-ray ionization rates (s −1 ), ζ W , (Eq. 22) for different values of D 0 . The value in the parenthesis indicates the power. The fiducial model uses the parameters α V = 1, β = 0.5, σ 0 = 1 km s −1 and L = 1 pc. The rest of the rows delineate models with a specific parameter variation: "Less Bound" corresponds to α V = 2, "Strong Turbulence" to σ V = 2.5 km s −1 and "Kolmogorov Turbulence" to β = 0.33. The ↓ represents ionization rates below the radionuclide ionization rate (Adams et al. 2014) , ζ RN ≈ 10 −19 .

where k = 2π/ and W (k) is the isotropic kinetic energy turbulence spectrum

and σ(k) = σ 0 k 0 k β . We find that for all cases, increasing D 0 (and thus allowing CRs to propagate more easily throughout the cloud) systematically increases the CRIR. For the "less bound" clouds, the ionization rate increases due to the decreased average density, and hence CRs lose less energy through the cloud. Similarly, for both the "strong turbulence" model and "Kolmogorov turbulence" model, for which the turbulence strength is not increased, the produced CRIR is increased due to the enhanced turbulent power throughout the driving scales.

Our "strong turbulence" model represents regions of significant driving, such as in regions of enhanced star formation feedback (e.g. Offner & Liu 2018) (e.g., nearby protostar jets, high-mass stars and supernovae) or in the Galactic Center (Kauffmann et al. 2017) . Due to the strength of the turbulence, there is a significantly enhanced produced CRIR, far exceeding that observed in Solar neighborhood clouds. However, CRIRs on the order of 10 −14 s −1 are observed through H + 3 absorption towards the Galactic Center (Indriolo et al. 2015) .

Most of the clouds in the Milky Way are not entirely virialized, and exhibit virial parameters greater than 1 (Heyer & Dame 2015) . Therefore, our model predicts that within these clouds, reconnection within the MHD turbulence produces enough MeV -GeV protons to sustain CRIRs, ζ > 10 −16 s −1 .

This mechanism directly correlates the CRIR and the properties of the magnetohydrodynamic turbulence within molecular clouds, along with the transport physics of low-energy CRs. Therefore, it may be possible to verify this mechanism with co-spatial observations of the ionization rate in dense gas, the magnetic field strength and the turbulence properties through observations of molecular ions and dust polarization maps. However, inferring the CRIR from such observations will rely on understanding the diffusion coefficient. Conversely, if the CRIR is dominated in dense gas by our proposed mechanism, it may be possible to infer the properties of the magneto-hydrodynamic turbulence from the CRIR through backwards modelling.

It is worth discussing the great uncertainty in the diffusion coefficient. Within the Milky Way, cosmic-ray transport studies and observations have indicated an average diffusion coefficient between D 0 = 10 28 -3 × 10 28 cm 2 s −1 (Evoli et al. 2020) . However, regarding the dense gas, studies have shown a spread over several orders of magnitude, from D 0 = 10 27 -10 30 cm 2 s −1 (Dogiel et al. 2015; Owen et al. 2021 ). As such, it is even more paramount to understand what is constraining the CR transport within molecular gas, and how the local environment changes the diffusion coefficient. Therefore, a widespread is observed, although dwarf galaxies appear to necessitate a higher ionization rate. If this is the case, our model predicts that dwarf galaxies would have CRIRs significantly greater than Milky Way-type galaxies.

We have proposed a novel mechanism for the production of low-energy CRs in molecular clouds through Fermi acceleration in regions undergoing turbulence reconnection: CRAFT. Since the MHD turbulence cascades across a wide range of scales, and since both the turbulence and the reconnection are volume-filling processes, we expect this will produce an approximately isotropic and homogenous floor to the CRIR. Historically, there has been a contradiction between the constant CRIRs used in astrochemical models (Röllig et al. 2007 ) and theoretical calculations, which have ubiquitously shown that energy losses would produce steep gradients with a low CRIR towards the cloud's center. Furthermore, observations indicate the ionization rate in dense gas is not significantly lower than more diffuse regions. Our results would instead show that a properly chosen constant CRIR may be actually appropriate when modeling the dense gas in molecular clouds.

Proceedings of the National Academy of Science

Cosmic Ray Astrophysics Schlickeiser

We thank the anonymous referee for their useful comments improving this work. This work was funded by theww ERC starting grant No. 679852 'RADFEEDBACK'. AL acknowledges the support by NASA TCAN AAG1967 and NSF AST 1816234. B.A.L.G would also like to thank his canine office mate, Mojo Gaches, for providing constant support during this work during the Coronavirus pandemic, although due to constantly sleeping on the job is not a coauthor.