Exclusion of cosmic rays from molecular clouds by self-generated electric fields
DDraft version October 1, 2020
Preprint typeset using L A TEX style emulateapj v. 01/23/15
EXCLUSION OF COSMIC RAYS FROM MOLECULAR CLOUDS BY SELF-GENERATED ELECTRIC FIELDS
Kedron Silsbee , Alexei V. Ivlev Max-Planck-Institut f¨ur Extraterrestrische Physik, 85748 Garching, Germany
Draft version October 1, 2020
ABSTRACTIt was recently discovered that in some regions of the Galaxy, the cosmic ray (CR) abundanceis several orders of magnitude higher than previously thought. Additionally, there is evidence thatin molecular cloud envelopes, the CR ionization may be dominated by electrons. We show that forregions with high, electron-dominated ionization, the penetration of CR electrons into molecular cloudsis modulated by the electric field that develops as a result of the charge they deposit. We evaluatethe significance of this novel mechanism of self-modulation and show that the CR penetration can bereduced by a factor of a few to a few hundred in high-ionization environments, such as those foundnear the Galactic center. INTRODUCTION
Understanding the transport of cosmic rays (CRs) indense gas is one of the big open questions of astrophysics.Low-energy CRs govern the evolution of molecular cloudsand the formation of stars (Caselli & Ceccarelli 2012;Padovani et al. 2020), being the dominant source of ion-ization (McKee 1989; Keto & Caselli 2008; Neufeld &Wolfire 2017) and UV emission (Prasad & Tarafdar 1983)above a column density of ∼ cm − . These pro-cesses affect both the chemistry (Keto & Caselli 2008;Keto et al. 2014) and thermodynamics (Galli et al. 2002;Glassgold et al. 2012a; Ivlev et al. 2019) of the clouds.Furthermore, the level of ionization governs the degree towhich the gas is coupled to the magnetic field (Shu et al.1987). This has profound implications for the existenceand size of disks around young stars (Zhao et al. 2016,2018).For a long time the CR ionization rate ζ was thoughtto likely be on the order of 10 − s − , based on mea-surements of the CR abundance near Earth (Spitzer &Tomasko 1968). More recently there have been measure-ments of ζ in nearby molecular clouds (Indriolo & Mc-Call 2012), suggesting ζ as high as 10 − s − towardssome clouds. In some environments the CR ionizationrates can be orders of magnitude higher still. In theCentral Molecular Zone (CMZ) of the Galaxy, Le Pe-tit et al. (2016) and Oka et al. (2019) estimate ζ of1 − × − s − and 2 × − s − respectively. Yusef-Zadeh et al. (2007) suggest a rate of 5 × − s − in theSagittarius C region. There is also evidence of extremelyhigh CR abundance near young stars (Ceccarelli et al.2014; Ainsworth et al. 2014).It is not known whether CR protons or electrons arethe primary source of ionization. As shown in Padovaniet al. (2018), if the spectra of electrons and protons mea-sured by the Voyager probes (Stone et al. 2019) are ex-trapolated down to lower energies, then ζ is dominatedby electrons at column densities lower than 2 × cm − .If the CR electron and proton spectra have the slope ap-propriate for acceleration in strong shocks, then ioniza-tion is dominated by electrons unless protons dominate e-mail: [email protected]: [email protected] the total CR energy by factors of tens.The commonly used free-streaming model for the CRtransport in clouds and disks (Padovani et al. 2009,2018) holds that they propagate along local magneticfield lines without substantial pitch-angle scattering, andlose energy due to interactions with the gas in the cloud(losses are dominated by ionization for non-relativisticparticles). The fundamental effect completely neglectedin such models (also those including CR scattering onmagnetic disturbances) is an inevitable net deposition ofcharge within the cloud. Low-energy CRs are absorbedin the cloud, becoming thermalized charged particles. Soas to maintain charge balance within the cloud, the ther-mal plasma must transport a net current. Because theplasma has a finite conductivity, this implies the pres-ence of a long-range electric field which acts to modulatethe penetration of CRs.This mechanism of CR modulation – which has notbeen considered thus far, to our knowledge – is the topicof the present Letter. We show that the self-generatedelectric field is strong enough to have a large effect forreasonable parameters, provided CR electrons dominatethe ionization. LINEAR REGIME
Let us calculate the steady-state electric potential in acloud in the limit that the incoming CR flux is not mod-ulated by the electric field. We approximate the cloud asa slab of a weakly ionized cool gas with uniform density n , embedded in a warm infinitely conducting mediumfilled with CRs. The magnetic field lines are assumedto be straight, but may enter the cloud at an arbitraryangle with respect to the surface. The column density N relevant to the CR propagation is defined by integratingthe density along the magnetic field . The distance z ismeasured in this direction, too, and set to 0 at the cloudcenter, where N is half of the total cloud value N cl . Fortypical diffuse clouds (even with the extreme ionizationimplied by our model), the plasma conductivity paral-lel to the magnetic field is higher by at least 6 orders ofmagnitude than the perpendicular conductivity. Hence,the current due to charge deposition by CRs in the re-gion between 0 and z within the cloud must be simplybalanced by the parallel plasma current at position z . a r X i v : . [ a s t r o - ph . H E ] S e p Keeping in mind that the cloud is bombarded by CRsfrom both sides, the CR current J CR at column N isgiven by the integral J CR ( N ) = 2 πq CR (cid:90) dµµ (cid:90) E ext ([ N cl − N ] /µ ) E ext ( N/µ ) d E j i ( E ) , (1)where q CR = ± e is the charge of a CR particle and µ is the cosine of the pitch angle. The integration limitsare determined by the extinction energy E ext ( N ) – thelowest energy of a particle that can penetrate to columndepth N . The external (initial) spectrum of CRs, j i ( E ),is assumed to be isotropic and given by j i ( E ) = j (cid:18) E E (cid:19) a s − cm − eV − sr − . (2)To determine E ext ( N ), we must introduce the ionizationloss function. This is given in Padovani et al. (2018) as L ≡ d E dN = L (cid:18) E E (cid:19) s , (3)and then E ext ( N ) = E (cid:18) NN (cid:19) s , (4)where N = (1 + s ) − E /L . Performing the integral inEquation (1) yields J CR ( N ) = 2 π (1 + s ) q CR j E (1 − a )(1 + 2 s + a ) (cid:104) ( ˜ N cl − ˜ N ) − a s − ˜ N − a s (cid:105) , (5)where ˜ N ≡ N/N . For a = 1: J CR ( N ) = πq CR j E s ln (cid:18) N cl − NN (cid:19) . (6)If we are interested in column densities between 10 and 10 cm − , the relevant energies are from 5 to300 KeV for electrons, and from 100 KeV to 8 MeVfor protons. The loss functions on these intervals arewell approximated by Equation (3) with s = 0 . L =2 . × − eV cm , E = 10 KeV for electrons, and s = 0 . L = 3 . × − eV cm , E = 10 MeV forprotons. Even though s may change a little, dependingon the energy range, these variations have negligible im-pact on our results. Therefore, in this Letter we employthe above values for numerical calculations, while keep-ing s explicitly in the analytical results. As in Padovaniet al. (2018), we use the number density of all gas par-ticles, rather than of hydrogen atoms, and assume thehydrogen to be molecular.Let us denote with E the electric field component par-allel to the magnetic field. The value of E ( N ) in thecloud is obtained from the the steady-state condition J CR ( N ) + J pl ( E ) = 0 , (7)where J pl = σE is the plasma current along the magneticfield, determined by the corresponding electric conduc-tivity (Braginskii 1965), σ = 5 × s − (cid:18) T
50 K (cid:19) / . (8) Equation (8) assumes a fully ionized plasma. The ex-pected ionization fraction in the outer layers of a cloudis in excess of 10 − (Neufeld & Wolfire 2017). Since theelectron-neutral collision cross section is lower by 5 or 6orders of magnitude than the electron-ion cross sectionfor such conditions, the neutrals have a negligible effecton the parallel conductivity.Consider the CR spectra measured from the Voyagerprobes (Stone et al. 2019) and extrapolated to lower en-ergies as in Padovani et al. (2018), and a cloud with N cl = 6 × cm − , n = 60 cm − , and T = 50 K (Draine2011). From Equations (5), (7), and (8) we derive themagnitude of the electric potential energy e | φ ( N ) | forCR protons ( p ) and electrons ( e ), and compare thesewith the respective extinction energies E ext ( N ). We ob-tain that eφ p is completely negligible at any N , while e | φ e ( N ) | is just a factor of 30 less than E ext ,e ( N ). Sincethe electron spectrum dominates the ionization at lowcolumn density, this implies that if the CR abundancewere increased by a factor of 30, then ζ would be sig-nificantly affected by the electron charge buildup. Suchan increased spectrum would result in a total ionizationrate of somewhat less than 10 − s − at a column densityof 10 cm − . This value is within the range estimatesmade by Indriolo & McCall (2012), suggesting that thiseffect may play a role, even in local molecular clouds.The importance of the electric field (at a given ζ )is substantially higher if the ionization is dominatedby electrons. Consider locations adjacent to a strongshock which is acting as a source of CRs, so one canexpect F ( p ) ∝ p − for the particle density in momen-tum space . Assuming an electron to proton ratio of χ ,and column densities such that the ionization is dom-inated by non-relativistic particles (up to a few times10 cm − for electrons, and a few times 10 cm − forprotons), we obtain the energy spectra given by Equation(2) with j ,e /j ,p = χm p /m e . Then Equation (6) showsthat the deposited charge is dominated by electrons if χ > m e /m p . From Equation 31 of Silsbee & Ivlev (2019),we find that the ratio of the ionization rates at a givencolumn density is ζ e /ζ p ∼ χ ( m p /m e ) s/ (1+ s ) (for simplic-ity, we assume the same s for electrons and ions and set L ,e /L ,p ∼ m e /m p ). From this we conclude that, if χ isgreater than a few percent, then both the ionization andthe charge deposition are dominated by electrons.Studies of particle acceleration are still uncertain asto the value of χ – there is evidence that it is less than1% in quasi-parallel shocks (Park et al. 2015). On theother hand, there is recent evidence (Spitkovsky et al.2019) that quasi-perpendicular shocks in fact preferen-tially accelerate electrons. Hence, it is reasonable to as-sume that there are regions where more than a few % ofthe electrons have been produced in quasi-perpendicularshocks, and therefore for the remainder of this Letter,we consider regions in which CR electrons dominate theionization rate. HIGH-FLUX LIMIT FOR CR ELECTRONS
We now approach the problem from a different per-spective. Instead of assuming the electric field to be asmall perturbation on the propagation of CRs, we con-sider it to be the dominant effect and treat ionizationlosses as a perturbation. To be more precise, we assumethat at every position of interest within the cloud, theelectric potential φ satisfies e | φ | (cid:29) E ext .Let us consider, as before, a slab of uniform gas witha constant angle between the magnetic field and the sur-face. Now the distance coordinate z , measured along themagnetic field, is set to be 0 at one edge of the cloud.We posit that the absolute value of the potential as afunction of z over some range of distances is given by | φ ( z ) | = E e (cid:18) zz (cid:19) f , (9)with the length scale z and exponent 0 < f < φ ( z ) isnot valid near the center of the cloud, so we restrict ourattention to z much smaller than the cloud size. Thisallows us to solve for E , the electric field componentparallel to the magnetic field, as a function of position: E ( z ) = E (cid:18) zz (cid:19) f − , (10)where E = f E / ( ez ).Instead of using variables E and µ , for our prob-lem below we find it more convenient to work withthe “parallel” and “perpendicular” energies, E (cid:107) = E µ and E ⊥ = E (1 − µ ). The local spectrum per unit E (cid:107) and E ⊥ can be conveniently calculated from the localdensity in the momentum space. According to the Li-ouville theorem, the CR density in momentum spaceis conserved along the phase trajectories, that is tosay F ( p , r ) = F i ( p + 2 me | φ ( r ) | ). Combining thiswith a general relation j ( E , µ ) = p F ( p, µ ), we find F ( p, z ) = j i ( E + e | φ | ) / [2 m ( E + e | φ | )]. Then, notingthat 2 πp ⊥ dp ⊥ dp (cid:107) = πm (cid:112) m/ E (cid:107) d E (cid:107) d E ⊥ , we multiply F ( p, z ) with the physical velocity (cid:112) E /m and the pre-factor πm (cid:112) m/ E (cid:107) , which yields the spectrum expressedin new variables, j ( E (cid:107) , E ⊥ , z ) = π (cid:115) EE (cid:107) j i ( E + e | φ | )( E + e | φ | ) , (11)where we use E = E (cid:107) + E ⊥ for brevity.As the first step, we equate the current of CRs whichis absorbed beyond position z due to the losses and theplasma current along the magnetic field. Integrating overthe initial CR distribution at the cloud edge, we obtain σE ( z ) e = (cid:90) ∞ e | φ ( z ) | d E (cid:107) (cid:90) E cr ⊥ ( E (cid:107) )0 d E ⊥ j ( E (cid:107) , E ⊥ , , (12)with σ from Equation (8). Particles with E ⊥ = 0will have zero kinetic energy at the turning point, andwill therefore be stopped and contribute to the chargebuildup. Particles with relatively large E ⊥ have enoughtransverse energy that they are accelerated back to thecloud edge before their energy is damped. Hence, forparticles with a given E (cid:107) there is a critical value of E ⊥ ,denoted E cr ⊥ , which determines their trapping inside thecloud. As discussed in Appendix A, the dynamics of aparticle with initial transverse energy E ⊥ in the presenceof losses are determined by a dimensionless number M = eE turn nL ( E ⊥ ) , (13) where E turn is the parallel electric field at the turn-ing point. The critical initial transverse energy E cr ⊥ for E (cid:107) = e | φ ( z ) | corresponds to M cr ≈ .
6: For
M < M cr ,particles are stopped by the losses near the turning point,otherwise they return back to the cloud edge. UsingEquations (9) and (10), we find that the electric fieldat the turning point is E turn = E (cid:18) E E (cid:107) (cid:19) − ff . (14)Combining Equations (13) and (14), we find E cr ⊥ = E (cid:18) M cr nL eE (cid:19) s (cid:18) E (cid:107) E (cid:19) − fsf . (15)Plugging Equation (15) into Equation (12) and approxi-mating that E cr ⊥ (cid:28) E (cid:107) (which is verified in Section 3.2),we can evaluate the integral in Equation (12) under con-dition f (1 + as ) >
1. This yields σE e (cid:16) z z (cid:17) − f = πj E (cid:18) M cr nL eE (cid:19) s × sff (1 + as ) − (cid:16) z z (cid:17) f (1+ as ) − s . (16)Matching powers of z , we find f = 1 + s s + as . (17)This allows us to solve for z : nz N = (cid:32) (1 + s ) s − f s +1 M s cr asσnπe j N (cid:33) s s . (18)We note that the condition f (1 + as ) > a >
0. We finally derive e | φ | E ext = (cid:18) N nz (cid:19) s s + as (cid:18) NN (cid:19) s (1+ s − a )(1+ s )(1+ s + as ) , (19)naturally, invariant with respect to the choice of E . Werequire e | φ | / E ext (cid:29) N , depending on the sign of the slope. Infact, however, the slope is very small: ≈ .
13 ( − .
03) for a = 1 ( a = 2). For this reason, as a practical matter, overthe range of column density of interest the solution eitherapplies everywhere, or applies nowhere – depending onthe magnitude of nz /N , which is the chief parametercharacterizing the effect of self-generated field.Now we can calculate the ionization rate in the high-flux limit. Again, we assume that the regular (ionization)losses play no role in determining the local CR spectrum,which is determined purely by the external spectrum andthe electric potential. The primary CR ionization rate ofH at position z is given by ζ φ ( z ) = 2 (cid:15) (cid:90) ∞ (cid:90) ∞ d E (cid:107) d E ⊥ j ( E (cid:107) , E ⊥ , z ) L ( E ) , (20)where j ( E (cid:107) , E ⊥ , z ) is given by Equation (11) and (cid:15) is themean energy lost per primary ionization event (Silsbee& Ivlev 2019), which we take to be 58 eV. We obtain ζ φ = 4 πBj L E (cid:15) (cid:18) e | φ | E (cid:19) − ( a + s − , (21)where B ≡ B (2 − s, a + s −
1) is the beta function (seeAppendix B). By comparing Equation (21) with Equa-tion (31) of Silsbee & Ivlev (2019), which describes the“regular” CR ionization rate ζ ( N ), we obtain ζ φ ζ ≈ . (cid:18) e | φ | E ext (cid:19) − ( a + s − , (22)where e | φ | / E ext is given by Equation (19) and the pre-factor is accurate within 3% for 1 ≤ a ≤
2, see Equation(B3).We point out that the sign of a + s − ζ ( N ), Equation 33of Silsbee & Ivlev (2019). In case a + s − < a + s − > Magnitude of the effect
Using Equations (3), (8), and (21), we rewrite Equa-tion (19) in terms of the physical parameters: e | φ | E ext ≈ . . a ) T − . N . ( ζ − /n ) . , (23)with T in units of 250 K, N in units of 10 cm − , ζ − (evaluated at same column density) in units of10 − s − , and n in units of 30 cm − . Equation (23)is accurate to within 2.5% for 1 ≤ a ≤ a = 1 (appropriate for acceleration by strong shocks)and a = 2 (for weaker shocks with compression ratioof 2) (Blandford & Ostriker 1978). In both cases, wechoose j so that (with the electric field included) theprimary ionization rate at N = 10 cm − is equal to4 × − s − , based on the values of 1 − × − s − reported in Le Petit et al. (2016) for the CMZ.The top panel of Figure 1 shows a comparison of e | φ ( N ) | with E ext ( N ). We use T = 250 K and n =30 cm − , based on the observations in Le Petit et al.(2016). The bottom panel shows a comparison of theionization rate calculated from Equation (21) with thatcalculated ignoring electric fields (Equation 31 of Silsbee& Ivlev (2019)). At a representative column density of10 cm − , e | φ | / E ext ≈ a = 1,and ≈
11 for a = 2, leading to reductions in the ion-ization rate by factors of about 2.7 and 40, respectively.Note that for T ≈
50 K, suggested by Figure 9 of Bisbaset al. (2015) for our values of ζ , the reduction would beabout 5.3 and 200, respectively.Limits of the non-relativistic formulation are reachedat higher column densities, where e | φ | (cid:38) mc . In Ap-pendix C, the calculations presented in Equations (11)–(19) are redone in the ultra-relativistic regime, assum-ing that outside the cloud the CR density in momentum N , cm , s a = 1: electric fielda = 1: regular lossesa = 2: electric fielda = 2: regular losses E n e r g y , e V ext e | |, a = 1 e | |, a = 2 Fig. 1.—
The top panel shows a comparison between the electricenergy e | φ | [Equation (19)] and the extinction energy E ext [Equa-tion (4)] for CR electrons. The black and red curves are for anelectron spectrum with a = 1 and 2, respectively (normalized suchthat ζ = 4 × − at N = 10 cm − ). The bottom panel showsthe corresponding ionization rate, plotted without and with takinginto account the self-generated electric field [Equation 31 of Silsbee& Ivlev (2019) and Equation (21) of this Letter, respectively]. space has the same power-law slope for both relativis-tic and non-relativistic particles. It is also assumed thatthe critical kinetic energy is still non-relativistic near theturning point – this premise is shown to be valid for N substantially higher than 10 cm − , i.e., well applicablefor molecular clouds. We find that the electric potential(19) is modified in the ultra-relativistic regime as φ rel ( N ) = φ ( N ) (cid:18) NN rel (cid:19) − s (1+ a )(1+ s )(1+ s + as )(1+2 s +2 as ) , (24)Here N rel is the column density at which e | φ ( N ) | fromEquation (19) is equal to ≈ . mc , see Equation (C6).For Figure 1, N rel = 2 − × cm − and the exponentvaries between − .
26 and − .
22 for 1 ≤ a ≤
2. Hence,relativistic effects only lead to a minor modification ofthe self-generated electric field, and the CR modulationremains essentially unchanged.
Notes on the Derived Solution
Here we verify important assumptions made to derivethe above results, and briefly discuss some immediateimplications.
Anisotropy
We can check the assumption made after Equation(15), that E cr ⊥ (cid:28) E (cid:107) . Setting e | φ | ≈ E (cid:107) , we rewrite Equa-tion (15) as E cr ⊥ E (cid:107) ≈ (cid:18) M cr f (1 + s ) (cid:19) s (cid:18) E ext e | φ | (cid:19) s +1 . (25)For the parameters in Figure 1, the anisotropy due toCR deposition is expected to be less than a few %. Thismeans that the excitation of the streaming instability bypenetrating CRs will be suppressed compared to a caseof no electric field (Morlino & Gabici 2015; Ivlev et al.2018) (where the pitch angle anisotropy could be of orderunity). Electric Fields from Alfv´en Waves
We have implicitly assumed that the deposition of CRsis the only source of a large-scale electric field parallel tothe magnetic field. In fact, an electric field could alsobe produced by Alfv´en waves present in the cloud asa result of turbulence. In ideal MHD, such fields areperpendicular to the local magnetic field, and thus dono work on CRs. It is worth noting though that thestrength of this electric field, associated with turbulentmotions with the velocity u , is on the order of uB/c ∼ − statV cm − for typical parameters. This is ∼ times stronger than the field from the modulation effect,and therefore it could be significant if there were even avery small deviation from orthogonality.In particular, Bian et al. (2010) and Klimushkin &Mager (2014) suggest that under realistic conditions,Alfv´en waves are able to generate a small parallel electricfield, but its strength is highly uncertain. The authorsdeveloped a model under which the ratio of parallel toperpendicular electric field is roughly the squared ratioof the ion gyroradius to the wavelength. Estimates ofthe lower wavelength bound for Alfv´en waves, given inAppendix C of Kulsrud & Pearce (1969), suggest thatwaves are cut off around the ambipolar damping scale of ∼ cm for our conditions. For typical magnetic fields,the ion gyroradius is a few times 10 cm, so the parallelelectric field arising from such mechanism is smaller thanthe perpendicular field by some 19 orders of magnitude.We are not aware of a model which shows a significantparallel electric field in the long-wavelength limit. How-ever, as both the cutoff scale of Alfv´en waves and themagnitude of parallel electric field are subject to signifi-cant uncertainties, this could be an interesting avenue offuture work. Joule Heating
The electric current J pl = σE induced in the gas dueto CR deposition represents an additional source of gas heating. The rate of the resulting Joule heating is H J = σE . (26)This should be compared to the rate of regular gas heat-ing by CRs, given by H CR = η(cid:15)ζ φ n, (27)where η is an efficiency factor of order 40% (Glassgoldet al. 2012b). As shown in Appendix D, their ratio is H J H CR = Q (cid:18) e | φ | E ext (cid:19) − s + s , (28)where Q varies between ≈ ≈
10 for 1 ≤ a ≤ ≤ a ≤ a , Joule heating becomes subdominant inthe limit of very strong self-modulation. CONCLUSION
We propose a novel mechanism of CR self-modulation,which can substantially reduce the penetration of CRelectrons into molecular clouds. The penetration is lim-ited by the electric fields generated due to the depositionof those same electrons. If the electron spectrum is pro-duced by acceleration in strong shocks, the ionizationrate can be reduced by a factor of a few in the high-ionization environments found in our Galaxy, such asCMZ. The reduction becomes much stronger for steeperspectra, appropriate for weaker shocks. Hence, the highionization rates near the Galactic center could imply evenhigher CR energy densities than previously thought. Theeffect is more pronounced at lower gas number densi-ties, where direct measurements of the ionization can bemade. The ionization rate in denser regions will thereforebe much higher than would be predicted from measure-ments coupled with conventional models of CR trans-port. Furthermore, the electric current induced in thegas due to the CR deposition represents an additionalheating source. We show that the resulting Joule heatingcould be of similar magnitude to the regular gas heatingby CRs.
REFERENCESAinsworth, R. E., Scaife, A. M. M., Ray, T. P., et al. 2014, ApJLetters, 792, L18Bian, N. H., Kontar, E. P., & Brown, J. C. 2010, A&A, 519, A114Bisbas, T. G., Papadopoulos, P. P., & Viti, S. 2015, ApJ, 803, 37Blandford, R. D., & Ostriker, J. P. 1978, ApJ Letters, 221, L29Braginskii, S. I. 1965, Reviews of Plasma Physics, 1, 205Caselli, P., & Ceccarelli, C. 2012, A&A Review, 20, 56Ceccarelli, C., Dominik, C., L´opez-Sepulcre, A., et al. 2014, ApJLetters, 790, L1Draine, B. T. 2011, Physics of the Interstellar and IntergalacticMedium (Princeton: Princeton University Press)Galli, D., Walmsley, M., & Gon¸calves, J. 2002, A&A, 394, 275Glassgold, A. E., Galli, D., & Padovani, M. 2012a, ApJ, 756, 157—. 2012b, ApJ, 756, 157Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91Ivlev, A. V., Dogiel, V. A., Chernyshov, D. O., et al. 2018, ApJ,855, 23Ivlev, A. V., Silsbee, K., Sipil¨a, O., & Caselli, P. 2019, ApJ, 884,176Keto, E., & Caselli, P. 2008, ApJ, 683, 238Keto, E., Rawlings, J., & Caselli, P. 2014, MNRAS, 440, 2616Klimushkin, D. Y., & Mager, P. N. 2014, Ap&SS, 350, 579 Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445Le Petit, F., Ruaud, M., Bron, E., et al. 2016, A&A, 585, A105McKee, C. F. 1989, ApJ, 345, 782Morlino, G., & Gabici, S. 2015, MNRAS, 451, L100Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163Oka, T., Geballe, T. R., Goto, M., et al. 2019, ApJ, 883, 54Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A,614, A111Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, Space ScienceReviews, 216, 29Park, J., Caprioli, D., & Spitkovsky, A. 2015, Phys. Rev. Lett.,114, 085003Prasad, S. S., & Tarafdar, S. P. 1983, ApJ, 267, 603Shu, F. H., Adams, F. C., & Lizano, S. 1987, Ann. Rev. Astron.Astrophys., 25, 23Silsbee, K., & Ivlev, A. V. 2019, ApJ, 879, 14Silsbee, K., Ivlev, A. V., Padovani, M., & Caselli, P. 2018, ApJ,863, 188Spitkovsky, A., Xu, R., & Tsiolis, V. 2019, in AAS/High EnergyAstrophysics Division, AAS/High Energy AstrophysicsDivision, 107.10
Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971Stone, E. C., Cummings, A. C., Heikkila, B. C., & Lal, N. 2019,Nature Astronomy, 3, 1013 Yusef-Zadeh, F., Muno, M., Wardle, M., & Lis, D. C. 2007, ApJ,656, 847Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018,MNRAS, 473, 4868Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050APPENDIX
APPENDIX A: CALCULATION OF E CR ⊥ ( E (cid:107) ) As stated in the main text, we consider the limit e | φ | (cid:29) E ext . In this case only particles with sufficientlysmall E ⊥ can be trapped in a cloud, as they are slowednearly to a stop at the turning point, thus sufferingstrong ionization losses. To estimate E cr ⊥ ( E ), we use theloss function given by Equation (3). We are dealing withnon-relativistic particles, such that their velocities are v = (cid:112) E /m , where m is the particle mass.We note that for a particle moving parallel to an elec-tric field with strength E , there is a critical energy E cr such that the drag force nL ( E cr ) due to energy losses iscompensated by the acceleration from the electric field, E cr = E (cid:18) nL eE (cid:19) s . (A1)We assume that if E cr is reached after turnaround, thisoccurs in a short distance from the turning point, andwe can therefore calculate E cr ⊥ assuming a constant elec-tric field. This assumption is verified at the end of ourcalculation.In a constant electric field E , the equations of motionfor the parallel and transverse velocities are: m ˙ v (cid:107) = − eE − nL ( E ) v (cid:107) v , m ˙ v ⊥ = − nL ( E ) v ⊥ v . (A2)Normalizing the velocity to the initial transverse velocity,˜ v = v/v ⊥ i , and time to τ = mv ⊥ i nL ( E ⊥ i ) , (A3)we then arrive at the equations˙˜ v (cid:107) = − M − ˜ v (cid:107) ˜ v − s − , ˙˜ v ⊥ = − ˜ v ⊥ ˜ v − s − , (A4)containing a single dimensionless number M = eEnL ( E ⊥ i ) . (A5)To distinguish between the local and initial values, herewe identify the latter with the subscript i . The numericalsolution of Equations (A4) shows that for M < M cr ≈ . M/M cr ≤ . z turn satisfies∆ z < . v ⊥ i τ = 0 . M E ⊥ i / ( eE ). Noting that the turn-ing point occurs approximately where e | φ ( z ) | = E (cid:107) i andusing Equations (9) and (10), we find z turn = f E (cid:107) i / ( eE ),and ∆ zz turn < . Mf E ⊥ i E (cid:107) i . (A6)As shown in Equation (25), in the limit e | φ | (cid:29) E ext wehave E cr ⊥ / E (cid:107) (cid:28) z/z turn (cid:28)
1. Thus, the assumption that trapped particles are stopped near theturning point is justified.
APPENDIX B: DERIVATION OF EQUATIONS (21) AND(22)
We substitute L ( E ) and j ( E (cid:107) , E ⊥ , z ), given by Equa-tions (3) and (11), to Equation (20) and, normalizingenergies by e | φ | , readily obtain Equation (21) with thepre-factor proportional to the following double integral: I ( a, s ) = (cid:90) ∞ (cid:90) ∞ dpdq ( p + q ) / − s √ p ( p + q + 1) a +1 . (B1)We replace the integration variables p, q by x , y andrewrite the integral in polar coordinates with r = (cid:112) x + y and tan θ = y/x . Integrating over θ between0 and π/ r = (cid:112) t/ (1 − t ) yields I ( a, s ) = 2 B (2 − s, a + s − , (B2)expressed via the beta function.Equation (22) is derived by comparing Equation (21)with Equation 31 of Silsbee & Ivlev (2019), which de-scribes the “regular” CR ionization rate ζ ( N ). We notethat Equation 31 should be multiplied by 2 π , due to adifferent normalization of the CR spectrum in Silsbee &Ivlev (2019), and d denotes s in our Letter. Also, unlikeSilsbee & Ivlev (2019), we assume that the magnetic fieldhas constant strength, so there is no magnetic mirroring(Silsbee et al. 2018). The integral I f , entering Equation31 and given by Equation 32 of Silsbee & Ivlev (2019),can also be expressed via the beta function, by substitut-ing x s = t/ (1 − t ) for the integration variable. UsingEquation (4), we finally obtain ζ φ ζ = 2( a + 2 s ) B (2 − s, a + s − B (cid:16) s , a + s − s (cid:17) (cid:18) e | φ | E ext (cid:19) − ( a + s − . (B3)The pre-factor of ( e | φ | / E ext ) is a slowly varying functionof a , equal to ≈ . ≤ a ≤ APPENDIX C: THE ULTRA-RELATIVISTIC REGIME
We calculate the electric potential in the limit e | φ | (cid:29) mc , assuming that the kinetic energy at the turningpoint is still less than mc , so the loss function of Equa-tion (3) can be used.Consider interstellar CRs with the density in momen-tum space having the same power-law slope for both rela-tivistic and non-relativistic energies. The correspondingkinetic energy spectrum reads j i ( E ) = j (cid:18) mc E E + 2 mc E (cid:19) a . (C1)For non-relativistic particles, Equation (C1) is reducedto the spectrum of Equation (2), adopted in the paper.In the ultra-relativistic regime, the spectrum becomes j i ( E ) = j , rel ( E / E ) a with j , rel = (2 mc / E ) a j . Thisallows us to easily extend our calculations to the ultra-relativistic case.Following the same logic as in the main text, but using E = pc instead of E = p / (2 m ), we find j ( E (cid:107) , E ⊥ , z ) = 2 π E ⊥ j i ( E + e | φ | )( E + e | φ | ) , (C2)with E = (cid:113) E (cid:107) + E ⊥ , in lieu of Equation (11). WhileEquation (12) remains unchanged, the calculation of theinitial E cr ⊥ proceeds differently. We assume that thecritical energy near the turning point is non-relativistic.Then its value is still given by Equation (15). Next,we note that the perpendicular momentum p ⊥ is a con-served quantity (neglecting losses). Setting p ⊥ / (2 m ) atthe turning point equal to the RHS of Equation (15),we find the critical value of the initial transverse energy E ⊥ = cp ⊥ , E cr ⊥ = (cid:112) mc E (cid:18) M cr nL eE (cid:19) s (cid:18) E (cid:107) E (cid:19) − f rel2 sf rel , (C3)to be substituted in Equation (12). As before, we ap-proximate E cr ⊥ (cid:28) E (cid:107) and obtain an equation analogousto Equation (16), which yields f rel = 1 + s s + 2 as , (C4)and nz , rel N = (cid:34) (1 + s ) s − f s +1rel M s cr asσnπe j , rel N (cid:18) E mc (cid:19)(cid:35) s s . (C5)Equation (C5) is similar to Equation (18) where param-eters a , f , and j are replaced with the respective ultra-relativistic values, and the extra factor E / (2 mc ) origi-nates from the square-root factor in Equation (C3).Using Equation (C5), we obtain an ultra-relativistic re-lation e | φ | rel / E ext versus N . By comparing this with thenon-relativistic relation, Equation (19), we derive Equa-tion (24) where N rel is the column such that e | φ ( N rel ) | mc = 2 (cid:34) (cid:18) s + 2 as s + as (cid:19) s +1 (cid:35) a ≡ ψ ( a, s ) . (C6)For 1 ≤ a ≤
2, we have ψ ( a, s ) ≈ . N max where the assumption breaks down,we use Equation (13) with E cr ⊥ = mc , which givesthe electric field E max at that turning point. Substi-tuting this to E max = f rel ( n/N max ) | φ max | , which fol-lows from Equations (9) and (10), we obtain e | φ max | =( M cr /f rel ) N max L ( mc ). Next, we introduce the columndensity N ∗ ≈ × cm − at which the electron extinc-tion energy in Equation (4) is equal to mc . Combiningthe two equations, we derive e | φ max | mc = M cr f rel (1 + s ) N max N ∗ . (C7) Finally, by virtue of Equations (9) and (C6) we write e | φ max | /mc = ψ ( N max /N rel ) f rel . Equating with Equa-tion (C7) gives N max N rel = (cid:18) ψf rel (1 + s ) M cr N ∗ N rel (cid:19) − f rel . (C8)Note that N max is comparable to, or larger than N ∗ .For the conditions illustrated in Figure 1, N rel ≈ × cm − (2 × cm − ) for a = 1 (2), resulting in N max ≈ × cm − (2 × cm − ). Thus, N max (cid:29) N rel and our assumption is well justified for molecularclouds. APPENDIX D: JOULE HEATING
We calculate the ratio of Joule heating H J , given byEquation (26), to regular gas heating H CR by CRs, givenby Equation (27). Substituting Equation (21) into Equa-tion (27), we obtain H CR = 4 πηBj L E n (cid:18) e | φ | E (cid:19) − ( a + s − . (D1)Inserting E = f | φ | /z in Equation (26) and keeping inmind that E /L = (1 + s ) N , we can then write theratio as H J H CR = (1 + s ) f πηB σne j N (cid:18) N nz z z (cid:19) (cid:18) e | φ | E (cid:19) a + s +1 . (D2)Next, by virtue of Equations (4) and (18) this can bewritten as H J H CR = Q (cid:18) N nz (cid:19) s − s (cid:18) NN (cid:19) a + s s (cid:16) z z (cid:17) (cid:18) e | φ | E ext (cid:19) a + s +1 , (D3)where Q = (1 + s ) − s f − s M s cr ηasB , (D4)is a function of a and s (for given η ). Then, inserting z /z = ( eφ/ E ) − /f with f from Equation (17), we find H J H CR = Q (cid:18) N nz (cid:19) (cid:18) NN (cid:19) s (1+ s − a )(1+ s )2 (cid:18) e | φ | E ext (cid:19) s (1+ s − a )1+ s s − s . (D5)We notice that, using Equation (19), the first two fac-tors in the brackets can be expressed via e | φ | / E ext . Thisfinally yields H J H CR = Q (cid:18) e | φ | E ext (cid:19) − s + s . (D6)For parameters of Figure 1, Q varies between about 7and 10 for 1 ≤≤
We calculate the ratio of Joule heating H J , given byEquation (26), to regular gas heating H CR by CRs, givenby Equation (27). Substituting Equation (21) into Equa-tion (27), we obtain H CR = 4 πηBj L E n (cid:18) e | φ | E (cid:19) − ( a + s − . (D1)Inserting E = f | φ | /z in Equation (26) and keeping inmind that E /L = (1 + s ) N , we can then write theratio as H J H CR = (1 + s ) f πηB σne j N (cid:18) N nz z z (cid:19) (cid:18) e | φ | E (cid:19) a + s +1 . (D2)Next, by virtue of Equations (4) and (18) this can bewritten as H J H CR = Q (cid:18) N nz (cid:19) s − s (cid:18) NN (cid:19) a + s s (cid:16) z z (cid:17) (cid:18) e | φ | E ext (cid:19) a + s +1 , (D3)where Q = (1 + s ) − s f − s M s cr ηasB , (D4)is a function of a and s (for given η ). Then, inserting z /z = ( eφ/ E ) − /f with f from Equation (17), we find H J H CR = Q (cid:18) N nz (cid:19) (cid:18) NN (cid:19) s (1+ s − a )(1+ s )2 (cid:18) e | φ | E ext (cid:19) s (1+ s − a )1+ s s − s . (D5)We notice that, using Equation (19), the first two fac-tors in the brackets can be expressed via e | φ | / E ext . Thisfinally yields H J H CR = Q (cid:18) e | φ | E ext (cid:19) − s + s . (D6)For parameters of Figure 1, Q varies between about 7and 10 for 1 ≤≤ a ≤≤