Issue 
A&A
Volume 614, June 2018



Article Number  A111  
Number of page(s)  20  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201732202  
Published online  25 June 2018 
Cosmicray ionisation in circumstellar discs
^{1}
INAFOsservatorio Astrofisico di Arcetri,
Largo E. Fermi 5,
50125
Firenze,
Italy
email: padovani@arcetri.astro.it, galli@arcetri.astro.it
^{2}
MaxPlanckInstitut für extraterrestrische Physik,
Giessenbachstr. 1,
85748
Garching,
Germany
email: ivlev@mpe.mpg.de, caselli@mpe.mpg.de
Received:
30
October
2017
Accepted:
25
March
2018
Context. Galactic cosmic rays (CRs) are a ubiquitous source of ionisation of the interstellar gas, competing with UV and Xray photons as well as natural radioactivity in determining the fractional abundance of electrons, ions, and charged dust grains in molecular clouds and circumstellar discs.
Aims. We model the propagation of various components of Galactic CRs versus the column density of the gas. Our study is focussed on the propagation at high densities, above a few g cm^{−2}, especially relevant for the inner regions of collapsing clouds and circumstellar discs.
Methods. The propagation of primary and secondary CR particles (protons and heavier nuclei, electrons, positrons, and photons) is computed in the continuous slowing down approximation, diffusion approximation, or catastrophic approximation by adopting a matching procedure for the various transport regimes. A choice of the proper regime depends on the nature of the dominant loss process modelled as continuous or catastrophic.
Results. The CR ionisation rate is determined by CR protons and their secondary electrons below ≈130 g cm^{−2} and by electronpositron pairs created by photon decay above ≈600 g cm^{−2}. We show that a proper description of the particle transport is essential to compute the ionisation rate in the latter case, since the electron and positron differential fluxes depend sensitively on the fluxes of both protons and photons.
Conclusions. Our results show that the CR ionisation rate in highdensity environments, such as the inner parts of collapsing molecular clouds or the midplane of circumstellar discs, is higher than previously assumed. It does not decline exponentially with increasing column density, but follows a more complex behaviour because of the interplay of the different processes governing the generation and propagation of secondary particles.
Key words: cosmic rays / ISM: clouds / stars: protostars / atomic processes / molecular processes
© ESO 2018
1 Introduction
Ionisation plays a crucial role in cold, dense molecular cloud cores and in circumstellar discs, where it controls the chemical properties of the gas, the coupling to the magnetic field, and (in the latter case) the generation of turbulence via the magnetorotational instability. The primary agents of ionisation in dense gas at column densities much larger than the values at which interstellar (IS) UV photons are absorbed (i.e. at visual extinctions A_{V} ≳ 1 mag or column densities N ≳ 10^{21} cm^{−2}) are Xrays, cosmic rays (CRs), and decaying of radionuclides. Deep in the densest parts of molecular clouds, characterised by densities n ≳ 10^{5}–10^{6} cm^{−3} and A_{V} ≳ 50−100 mag, CR ionisation is dominant, leading to an ionisation fraction decreasing with density (McKee 1989; Caselli et al. 2002; Walmsley et al. 2004; Maret et al. 2006). In discs around young, active stars, or in molecular clouds close to supernova remnants or massive starforming regions, the situation is complicated by the effects of stellar Xrays (Glassgold et al. 1997), the possible exclusion of lowenergy CRs by stellar winds (Cleeves et al. 2013), and the presence of local sources of accelerated particles (Lee et al. 1998; Padovani et al. 2015, 2016). The efficiency of stellar Xrays to ionise the circumstellar gas depends on total fluxes and hardness of the spectra. For example, hard Xrays of energies 1, 6 and 20 keV are absorbed by column densities of 2.5 × 10^{22}, 4.5 × 10^{24} and 1.3 × 10^{25} cm^{−2}, respectively (Igea & Glassgold 1999).
In the shielded regions near to the disc midplane, where Xrays and CRs are strongly attenuated, radioactive elements may substantially contribute to the electron fraction. In this case the ionisation rate is proportional to the abundance of the radioactive element and its decay rate. Thus, shortlived radionuclides (SLRs; mostly^{26}Al with halflife 7.4 × 10^{5} yr) contribute comparatively more than longlived radionuclides (LLRs; mostly^{40}K with halflife 1.3 × 10^{9} yr), but decay faster. Assuming the^{26}Al abundance inferred for the early solar system from the analysis of meteorites (MacPherson et al. 1995; Umebayashi & Nakano 2009), SLRs can maintain an ionisation rate of ≈ 10^{−19} s^{−1} in the midplane of a disc for ≈1 Myr, while LLRslead to ionisation rates of only ≈ 10^{−22} s^{−1}, but on timescales much longer than the disc lifetimes (Umebayashi & Nakano 2009; Cleeves et al. 2013). However, since the average Galactic abundance of^{26}Al inferred from its γray emission is about one order of magnitude lower than the meteoritic solar system value (Diehl et al. 2006), the probability for a star of being born in a^{26}Alrich environment, as in the case of the Sun, is relatively low (Gounelle 2015).
The propagation of lowenergy CRs at low column densities, which is characteristic of the diffuse envelopes of molecular clouds, is mostly determined by resonant scattering on selfgenerated magnetohydrodynamic (MHD) waves (occurring on the scale of the particle gyroradius; see e.g. Cesarsky & Völk 1978; Everett & Zweibel 2011; Morlino & Gabici 2015; Ivlev et al. 2018). However, at gas densities higher than ≈ 10^{3} −10^{4} cm^{−3} such waves cannot affect the CR transport because they are completely damped from efficient ionneutral collisions (Ivlev et al. 2018).
At any given depth in a cloud or circumstellar disc, CRs are attenuated in a way that generally depends on characteristics of the ambient medium. In a cloud threaded by a largescale magnetic field, CRs perform helical trajectories along the local field lines, i.e. CRs gyrate many times before they collide with a particle of the medium. Therefore, if the field lines are strongly twisted, the effective column density seen by a CR particle at a given point can be much larger than the lineofsight column density at that point (Padovani & Galli 2011, 2013; Padovani et al. 2013b). The decrease of the CR ionisation rate follows a powerlaw behaviour as function of the effective column density, N, in the range 10^{20} −10^{24} cm^{−2}, corresponding to effective surface densities, Σ, below a few g cm^{−2} (Padovani et al. 2009, hereafter PGG09). At higher densities, the attenuation is generally assumed to be exponential with a characteristic scale of ≈ 96 g cm^{−2} (Umebayashi & Nakano 1981).
The main goal of this paper is to determine the attenuation of CRs at moderatetohigh surface densities (larger than a few g cm^{−2}) by includingboth the energy loss and the particle production mechanisms relevant for the inner regions of circumstellar discs or collapsing clouds and adopting appropriate models for the transport of primary and secondary CR particles. We show that above ≈ 130 g cm^{−2} the CR ionisation rate rapidly becomes dominated by electronpositron pairs that are locally produced by secondary photons. As the latter are insensitive to the magnetic field and propagate isotropically, above this threshold the ionisation is controlled by the lineofsight (rather than the effective) column density.
The paper is organised as follows: in Sect. 2, we summarise two reference models for the IS spectra of CR protons and electrons; in Sect. 3, we examine the dominant energy loss mechanisms for primary and secondary CR particles operating at different energies; in Sect. 4, we describe our modelling of the propagation and attenuation of primary CRs; the generation and propagation of secondary particles is described in Sect. 5; in Sect. 6, we compute the total CR ionisation rate, focussing on high column densities; in Sect. 7, we discuss implications of our results for various astrophysical environments; and in Sect. 8, we summarise our most important findings.
2 Interstellar CR spectra
The IS differential fluxes of CR nuclei (hereafter, IS CR spectra) at high energies, E ≳ 1 GeV nuc^{−1}, are well constrained by various sets of observations from groundbased to balloon and satellite detectors (e.g. Aguilar et al. 2014, 2015). On the other hand, lowenergy IS spectra, being strongly affected by solar modulation effects (e.g. Putze et al. 2011), cannot be reliably determined by the same means. The best available estimate of the spectra (both nuclei and electrons) at energies E ≲ 500 MeV nuc^{−1} is provided by the most recent Voyager 1 observations (Cummings et al. 2016), down to ≈ 1 MeV nuc^{−1} and ≈10 MeV for CR nuclei and electrons, respectively.
In this paper, we adopt the analytical expression for the IS spectra of CR electrons and protons^{1}, as described in Ivlev et al. (2015), (1)
where we slightly modify the values of the parameters E_{0} and a, to better reproduce the most recent Voyager 1 data release (see Table 1). A simple extrapolation of the Voyager 1 data at lower energies fails to reproduce the CR ionisation rate measured in diffuse clouds from H emission (Indriolo & McCall 2012). For this reason, we consider two different models for the CR proton spectrum: a “low” spectrum, , obtained by extrapolating the Voyager 1 data (a = 0.1), and a “high”spectrum, (a = −0.8). The resulting ionisation rates and their comparison with observational data are discussed in Ivlev et al. (2015). The and proton spectra must be considered as lower and upper bounds, respectively, to the actual average Galactic CR spectrum. Although it is generally assumed that the CR intensity measured by Voyager 1 spacecraft is not affected by the solar wind modulation, one should not forget that the observed magnetic field has not changed yet to the direction expected in the IS medium (Burlaga et al. 2013), suggesting the possibility that the spacecraft may reside in a region of compressed solar wind (Fisk & Gloeckler 2014; Gloeckler & Fisk 2015). Thus, care should be taken in interpreting the Voyager 1 measurements as representative of the Galactic spectrum.
In Sect. 4 (see Fig. 4), we show that at high column densities, typical of circumstellar discs, the propagated CR proton spectrum becomes independent of the assumptions on the slope at low energy.
3 Energy losses and attenuation of CRs
In order to calculate the ionisation induced by CRs in molecular clouds or circumstellar discs, we need to study the propagation and attenuation of CR species k (including secondaries), and derive their spectra j_{k}(E, N) as function of the column density, N, along the direction of propagation, i.e. along the local magnetic field. We consider the propagation of CRs in a semiinfinite medium and, hence, assume that half of IS CRs (with the energy spectra described in Sect. 2 and isotropic pitchangle distribution) are incident on the surface. The inclination of the magnetic field with respect to the surface can be arbitrary. To simplify the presentation, we first obtain results for zero pitch angle and compute the effect of the angle averaging on the ionisation in Sect. 6. For applications to a circumstellar disc, in Sect. 7 we consider CRs coming from both sides of the disc.
We assume that hydrogen is only in molecular form and use the IS medium composition by Wilms et al. (2000), summarised in Table A.1. The mean molecular weight of the medium, Ā, is (2)
where A_{Z} and f_{Z} are the mass number and abundance with respect to the total number of particles, respectively^{2}. The column density is related to the surface density, Σ = Ām_{p}N, where m_{p} is the proton mass. Numerically, the relation between the N and Σ is given by (3)
In order to evaluate the total CR ionisation rate, a careful treatment of showers of secondary species (photons, electrons, and positrons) produced by primary CRs through processes such as pion decay, bremsstrahlung (BS), and pair production is required. In the following subsections we describe different energy loss processes for each CR species interacting with particles of a medium of given composition.
Throughout this paper, subscripts and superscripts denote CR species and the interaction processes, respectively. The expression for the partial energy loss function, , for a particle of species k depends on the type of process l: if only a small fraction of the particle kinetic energy is lost in each collision with a particle of the medium, the process can be considered as continuous and described by the loss function (4)
where is the differential cross section of the process and E^{max} is the maximum energy lost in a collision (see e.g. Appendix B.4 for Compton scattering). In the extreme case, where the entire kinetic energy is lost in a single collision or the CR particle ceases to exist after the collision, the process is called catastrophicand the loss function becomes (5)
where is the cross section of the process. Where possible, we express the total energy loss function in terms of the loss functions for collisions with atomic hydrogen (L_{k,H}) or helium (L_{k,He}).
3.1 Protons
The proton energy loss function, L_{p}, is composed by two terms: ionisation losses at low energies () and losses due to pion production above the threshold energy E^{π} = 280 MeV (). Therefore, (6)
where ε^{ion} = 2.01 and ε^{π} = 2.17 account for the presence of heavy elements in the target medium. The two terms on the righthand side of Eq. (6) are described in detail in Appendix A.1.
3.2 Electrons and positrons
The electronand positron energy loss function, L_{e}, has contributions due to ionisation losses at low energies (), BS above ≈100 MeV (), and synchrotron above E^{syn} ≈ 1 TeV (). Then, L_{e} is given by (7)
where the factors, ε^{ion} = 2.01 and ε^{BS} = 2.24, are described in Appendix A.2. We note that the synchrotron loss term in Eq. (7) does not depend on the medium composition.
3.3 Photons
Photons are generated through BS by electrons and positrons and through decay of neutral pions (produced by CR protons). In Sects. 5 and 6, we show that the reverse process of electronpositron pair production by photons is crucial for ionisation at high N, so photon propagation should be carefully treated. The photon energy loss function, L_{γ}, is a sum of three terms: photoionisation (), Compton scattering (), and pair production losses (), (8)
where ε^{C} = 2.01 and ε^{pair} = ε^{BS}; the latter is because the pair production and BS are symmetric processes (see Appendix B.3). The dominant contribution to is provided by the Kshell photoionisation of heavy species (see e.g. Draine 2011). Detailed expressions for the three terms on the righthand side of Eq. (8) are given in Appendix A.3.
3.4 Loss functions
Figure 1 shows the total proton, electron, positron, and photon energy loss functions L_{k} (E) calculated for the IS medium composition given in Table A.1. For comparison, the loss functions in a medium of pure H atoms, L_{k,H}(E), are also plotted.Data for the ionisation by protons are taken from the Stopping and Range of Ions in Matter package (Ziegler et al. 2010); for the ionisation by electrons we adopt Dalgarno et al. (1999) and the National Institute of Standards and Technology database^{3}.
We notice a change in the asymptotic behaviour L_{p} ∝ E^{α} of the proton loss function, from α = 1.28 to α = 1.08, occurring at energy E^{as} (see Eq. (A.2)). The asymptotic behaviour of the electron and positron loss function changes from α =1 to α =2, due to the transition from the BS to synchrotron (catastrophic) losses at energy E^{syn} (see Appendix A.2). As for photons, the asymptotic behaviour of the loss function is determined by the pair production (catastrophic) losses with α = 1; at low energies, where photoionisation dominates, one can see small peaks (around 1 keV) due to Kshell ionisation of heavy species of the medium.
3.5 CSDA and catastrophic regimes
The continuous slowingdown approximation (hereafter CSDA; Takayanagi 1973) has been used so far to study the propagation of lowenergy CRs in molecular clouds. It is based on two assumptions: (i) that the energy losses are continuous, and (ii) that the pitchangle scattering (with respect to the local magnetic field) is negligible. If these assumptions are justified, then the loss function L_{k} (E) (Eq. (4)) entirely determines the modification and attenuation of the spectrum of a species k with column density N.
Figure 2 presents the socalled range functions (9)
In the CSDA framework, the kinetic energy of a CR particle decreases from E_{0} to E after traversing a column density (10)
The local spectrum at that N and energy E is then related to the IS spectrum at N = 0 and energy E_{0} via (see PGG09) (11)
where μ is the cosine of the pitch angle. The factor of 1/2 in Eq. (11) takes into account that only half of the IS CRs penetrates into the semiinfinite medium. This relation directly follows from a solution of the transport equation for the CSDA regime (Ginzburg & Syrovatskii 1964; Berezinskii et al. 1990), (12)
assuming no sources. As pointed out in Sect. 3, it is sufficient to analyse the CR propagation for μ = 1 to calculate the total ionisation rate (see Sect. 6). The CSDA is a very simple and efficient approach which, of course, has certain limitations (see also Sect. 4).
The CSDA is formally no longer applicable when catastrophic losses dominate, although in some cases it may still be used (with the corresponding loss function, Eq. (5)) to qualitatively describe propagation of CRs. Strictly speaking, when both continuous and catastrophic loss processes are present, L_{k} in Eq. (12) should include only the continuous processes, while the catastrophic processes (with the cross section σ_{k}) should be described by an additional term σ_{k}j_{k} on the lefthand side. In the following, we discuss the effect of catastrophic losses on propagation of highenergy CR electrons (Sect. 4.2) and employ a transport equation for this regime to describe propagation of secondary photons (Sect. 5.1.1).
4 Propagation of CRs at high column densities
4.1 Cosmic ray protons
At energies larger than E^{π} the interaction between CR protons and the medium leads to the production of pions. Since the pion rest mass is significant, CR protons lose a nonnegligible fraction of their energy in each collision (Schlickeiser 2002). Such losses cannot be formally considered as continuous, but treating them as catastrophic is not quite correct either since many collisions are still required to reduce the energy significantly. Furthermore, the role of elastic scattering between CR protons and target nuclei becomes increasingly important when pion production dominates losses because the relative contribution of the nuclear interactions increases. This effect can be understood by considering the momentum transfer cross section, σ_{MT}(E), which is made up of various contributions depending on the energy range (see upper panel of Fig. 3). Below about 1 MeV the elastic (Coulomb) scattering dominates, while at higher energies the CR proton interacts with the target nucleus. Between 1 MeV and 10 MeV there is a transition region where the elastic scattering (σ_{MT} ∝ E^{−2}) is modified by nuclear forces (σ_{MT} ∝ E^{−1}). Finally, above ≈1 GeV the momentum transfer cross section becomes energy independent, which is a manifestation of the hard spherelike scattering^{4}. The total momentum transfer cross section can be written as a function of the protonproton momentum transfer cross section (Appendix B.1) as (13)
where ξ = 1.49 accounts for heavy species in the ambient medium (see Eq. (B.5)).
One can easily quantify the relative importance of the elastic scattering for CR protons, as compared to the their attenuation (characterised by loss function L_{p}). We introduce the scattering parameter, (14)
which is the integral ratio of the characteristic stopping range owing to energy losses, , to the characteristic column density resulting in strong scattering, ~1∕σ_{MT}. The CSDA is a good approximation as long as , otherwise scattering is no longer negligible and a gradual transition to the diffusive transport takes place.
The lower panel of Fig. 3 shows the dependence versus E. The CSDA works for E ≲ 25 MeV, where ; according to Fig. 2, this corresponds to column densities N ≲ 7 × 10^{22} cm^{−2}. The scattering of CR protons becomes increasingly important at higher energies, and the diffusive regime operates above ≈ 1 GeV, where , corresponding to N ≳ 2 × 10^{25} cm^{−2}.
Thus, a CSDA solution obtained for the local spectrum of CR protons at low energies should be combined with the diffusion solution at high energies. In Appendix C, we describe the matching procedure for the two solutions. The exact value of the transition energy E^{tr} at which the solutions are matched, 25 MeV ≲ E^{tr} ≲ 1 GeV, turns out to have a minor effect on the final results. For example, the resulting ionisation rate ζ(N) varies by less than 10% in the corresponding range of column densities of 7 × 10^{22} cm^{−2} ≲ N ≲ 2 × 10^{25} cm^{−2}.
We nowobtain the solution for the diffusive regime, assuming continuous losses. The steadystate density of CR protons in the diffusiveregime is governed by the following equation (Ginzburg & Syrovatskii 1964), (15)
Here, ℓ is the coordinate along the local magnetic field (CR path), is the number of protons per unit volume and energy, related to j_{p} (E, ℓ) via (16)
where β is the proton velocity in unit of the speed of light, c. The diffusion coefficient is (18)
where n is the total particle number density, summed over all the medium species and the factor 3 accounts for the fact that diffusion occurs in three dimensions. Then, by substituting the definition of the energy loss function for protons, Eq. (17), introducing the timelike coordinate^{5} (19)
and taking into account that dN∕dℓ = n, we reduce Equation (15) to (20)
where . Equation (20) must be solved with the boundary condition , since only half of the IS flux penetrates into the semiinfinite medium, and zero initial condition, reflecting the fact that j_{p} (E, N) → 0 for E →∞. In analogy with the solution of the problem of heat diffusion in a halfspace (Tikhonov & Samarskii 2013), the solution is (21)
where T(E, E_{0}) ≡ T(E) − T(E_{0}).
In principle, elastic scattering of CR protons leads to new source and sink terms in the general transport equation associated with efficient energy exchange with hydrogen nuclei. In Appendix D we present a detailed discussion of this effect, and show that for realistic conditions these new terms can be safely neglected. Figure 4 shows local differential fluxes of CR protons at various surface densities Σ.
Fig. 1 Total energy loss functions of primary and secondary CR particles k (protons, electrons and positrons, and photons), computed for a medium composition given in Table A.1 (L_{k}, thick black lines) and for atomic hydrogen (L_{k,H}, thin grey lines). Protons (upper panel): ionisation losses (short dashed lines) and pion production (dotted lines); the vertical dotted line shows the energy, E^{as}, at which the proton loss function changes its asymptotic behaviour from α = 1.28 to α = 1.08. Electrons and positrons (middle panel): ionisation losses (short dashed line), BS (longdashed line), and synchrotron losses with the uncertainty on the magnetic field strength (shaded area, see Appendix A.2); the vertical dotted line shows the transition energy, E^{syn}, between BS (α = 1) and synchrotron (α = 2) losses. Photons (lower panel): photoionisation losses (dashdotted line), Compton scattering (dotted line), and pair production (shortdashed line). 
Fig. 2 Total range functions, R_{k}(E), of primary and secondary CR particles (thick black lines), Eq. (9). The inset shows the total range functions multiplied by Ām_{p}, to highlight the behaviour at large surface densities. For comparison, the range functions for atomic hydrogen are also plotted (thin grey lines). 
Fig. 3 Upper panel: total momentum transfer cross section for protonnucleus collisions, σ_{MT} (E). The elastic (Coulomb) scattering dominating at lower energies crosses over to the nuclear scattering at higher energies. Lower panel: scattering parameter , Eq. (14). The vertical grey stripe indicates a continuous transition from the CSDA regime, where and the proton scattering is unimportant, to the diffusive regime, where . 
Fig. 4 Interstellar (solid lines, labelled “IS”) and local (dashed lines) differential fluxes (or spectra) of CR protons. Only half the IS proton flux penetrates into the semiinfinite medium. Labels indicate the surface density in g cm^{−2} for the local spectra. Models of the IS proton spectra, (black) and (grey), are described by Eq. (1). 
4.2 Cosmicray electrons
Energy losses of electrons by BS overcome ionisation losses above E^{BS} ≈ 500 MeV (see middle panel of Fig. 1). As pointed out by Ginzburg & Syrovatskii (1964), the energy of a photon created by BS is generally of the order of the energy of the electron that generated it. Therefore, one can approximately treat BS as a catastrophic process. The effective cross section of this process corresponds to a column density of N^{BS} ≈ 1.5 × 10^{25} cm^{−2}, i.e. ≈ 58 g cm^{−2}.
As a consequence, CSDA slightly overestimates the electron population at E ≳ E^{BS}. However, electrons at these energies yield only a minor contribution to the ionisation rate. Our numerical results show that their effect is smaller than 2% at Σ ≈30 g cm^{−2}, and becomesvanishingly small at higher Σ. Furthermore, the results presented below in Sect. 6 demonstrate that – even in the CSDA regime – the ionisation by primary CR electrons is negligible compared to the contribution from CR protons at Σ >1 g cm^{−2} for all models of IS CRs considered in this paper.
5 Generation and propagation of secondary particles
Figure 5 presents the main ionisation routes associated with various secondary particles that can be producedby CRs in circumstellar discs. Cosmicray protons and electrons generate secondary electrons by primary ionisation (). Secondary electrons with energy larger than the H_{2} ionisation potential produce further generations of ambient electrons. Cosmicray protons colliding with protons create charged pions; through muon decay, they become electrons and positrons (π^{±}→ μ^{±}→ e^{±}), producing secondary ionisations. In addition, protonproton collisions create neutral pions decaying into photon pairs (π^{0} → 2γ). The second important source of photons is BS by electrons and positrons. One should also account for electronpositron pair production by photons (γ → e^{+} + e^{−}). In the following we give the equations needed to compute the differential flux of all the secondary particles.
5.1 Photons
The regimes of propagation for photons can be different depending on their energy. As shown in the upper panel of Fig. 6, photoionisation and pair production dominate below ≈ 5 keV and above ≈50 MeV, respectively. Both processes are catastrophic, i.e. photons disappear after the interaction with nuclei. As for Compton scattering, the relative average energy lost by a photon in each interaction with electrons is written (22)
where σ^{C}(E) is the Compton cross section (see Eq. (B.20)). For x = E∕(m_{e}c^{2}) ≫ 1 we have ⟨ΔE⟩∕E^{C} ≈ 1 − 4∕(3lnx) → 1, i.e. Compton losses become catastrophic. On the other hand, for x ≪ 1 photons transfer a small fraction of their energy, ⟨ΔE⟩∕E^{C} ≈ 3x∕2, and losses are continuous. However, below E ≈ 1 keV photoionisation is the dominant process, and losses become catastrophic again. This is shown in the lower panel of Fig. 6, where we plot the quantity ⟨ΔE⟩∕E_{γ} combining all energy loss processes for photons, (23)
In order to determine whether CSDA or diffusive regime is applicable in the intermediate energy range of 5 keV ≲ E ≲ 50 MeV, we compute the scattering parameter defined by Eq. (14). We integrate the ratio of the momentum transfer cross section (Eq. (B.21)) and the loss function for Compton scattering (see Sect. sectionlinking A.3), which yields in this energy range. This clearly implies a diffusive regime for photons with dominant Compton scattering.
In the following subsections we present the governing equations for the catastrophic and diffusive regimes, and discuss how their solutions can be matched.
Fig. 5 Ionisation diagram, explaining the effect of secondary particles that are generated (directly or indirectly) by CR protons and electrons through ionisation, pion decay (π^{0}, π^{±}), BS, and pair production (pair). The secondary particles include electrons (e^{−}, due to primary and secondary ionisation), positrons and electrons (e^{±}, due to pair production and π^{±} decay; also electrons produced in secondary ionisation are included), and photons (γ, due to BS and π^{0} decay), all contributing to the respective ionisation routes. 
Fig. 6 Upper panel: components of the cross section for photons interacting with nuclei via process l: photoionisation (σ^{PI}), Compton scattering (σ^{C}), and pair production (σ^{pair}). The momentumtransfer (MT) cross section is plotted for Compton scattering, while for catastrophic (PI and pair) processes the cross section coincides with the MT cross section. The grey lines depict the corresponding loss functions near to their crossing (in arbitrary units). The vertical grey stripes, indicating the energy intervals between the crossing points of the corresponding cross sections and loss functions, separate the diffusive and catastrophic regimes of the photon transport (see text for details). Lower panel: the relative average energy lost by a photon per collision, ⟨Δ E_{γ} ⟩∕E_{γ}, vs. the photon energy, Eq. (23). 
5.1.1 Catastrophic regime
The equation of radiative transfer for the differential flux of photons j_{γ} (E, N) (photon fluxper unit energy and solid angle) is (24)
Here, μ is the cosine of the angle with respect to the direction of CR propagation and σ_{γ} = σ^{PI} + σ^{pair} is the cross section accounting for the two catastrophic processes described in the previous section, i.e. photoionisation and pair production (see Appendices B.2 and B.3). The source function of photons, S_{γ} (E, N), namely the number of photons per unit time, energy, and solid angle produced per nucleus, (25)
is the sum of the source function for π^{0} decay from protonnucleus collisions, (26)
(each proton provides two photons), and the source function for BS (27)
where the factor 2 accounts for electrons and positrons. Here, is the differential cross section for photon production by π^{0} decay (Kamae et al. 2006), and is the differential cross section of atomic hydrogen for BS (Blumenthal & Gould 1970; see Eq. (B.8)).
Neglecting any source of photon radiation external to the cloud (i.e. at N = 0), and averaging over μ, the solution of Eq. (24) gives the photon differential flux in the catastrophic regime, (28)
The factor 1/2 in Eq. (28) takes into account the fact that only the forward (backward) propagating photons produced at N^{′} < N (N^{′} > N) contribute to the local differential flux, at a given column density N.
5.1.2 Diffusive regime
The diffusive regime of photons is conceptually similar to that of CR protons and therefore is described by a similar equation (Eq. (15)), but with additional source terms owing to the photon production by neutral pion decay and BS (Eqs. (26) and (27), respectively). The diffusion equation is then given by (29)
where the number of photons per unit volume and energy is . Substituting and (30)
we reduce Eq. (29) to (31)
where is a function of column density N and timelike coordinate (32)
Similar to the catastrophic regime, we can set zero boundary condition, , while the initial condition is . The latter condition is determined by matching the diffusive and catastrophic regimes at E = E^{tr} (specified below).
Using again the analogy with the nonhomogeneous heat diffusion problem in a halfspace (Tikhonov & Samarskii 2013), we obtain the following solution for the photon differential flux: (33)
are determined by the Green’s function (36)
In Sect. 5.1, we mentioned that the catastrophic solution obtained in Sect. 5.1.1 for the high and lowenergy catastrophic regimes must be combined with Eq. (33). Similar to the treatment of CR protons (Sect. 4.1), we introduce the transition energy E^{tr} at which the two regimes should be matched. The matching criteria are determined (i) by the applicability of the diffusion approximation, which requires that , and (ii) by a continuous transition between the solutions.
By comparing Eqs. (24) and (29) we infer that, depending on the relative importance of different terms in the equations, E^{tr} may vary: as shown in Fig. 6, the matching occurs in the energy interval limited by the intersection points of the respective cross sections or loss functions. The lower panel of Fig. 6 shows that the transition to the highenergy (pair) catastrophic regime occurs at E^{tr} ≈ 50 MeV, where both and intersect with σ^{pair}(E) and , respectively. Concerning the matching with the lowenergy (PI) catastrophic regime, the intersection of and σ^{PI} (E) is at E ≈ 6 keV, whereas and cross at E ≈ 20 keV. In this case, the crossing points do not coincide, which is easy to understand: for process k, the loss function is generally related to the corresponding cross section via L^{k} ~⟨ΔE⟩σ^{k}. At such energies, ⟨Δ E⟩∕E_{γ} ~ 10^{−1} for Compton scattering, while for photoionisation ⟨ΔE⟩∕E^{PI} = 1 by definition. Hence, the lowenergy catastrophic and diffusion solutions should be matched at 6 keV ≲ E^{tr} ≲ 20 keV. The exact value of E^{tr} turns out to have a negligible effect on the resulting ionisation rate, as for CR protons.
5.2 Electrons and positrons
Electrons and positrons have two different sources. First, pairs are produced by photons with energy above 2m_{e} c^{2}, so that the electron and positron energy is 0 ≤ E ≤ E_{γ} − 2m_{e}c^{2}. The resulting source function at a given column density for a single species (electron or positron), (37)
is determined by the differential cross section (Eq. (B.11)). Second, electrons and positrons are created through decay of charged pions, generated in protonnucleus collisions at energies above E^{π}. The corresponding source function is given by (38)
where is the differential cross section for electron and positron production by π^{±} decay (Kamae et al. 2006), which we assume to have the same scaling for target heavy nuclei as that for π^{0} production (see Appendix A.1). The total source function for electrons and positrons is then (39)
As we pointed out in Sect. 4.2, the use of CSDA to describe propagation of electrons and positrons with energies above the BS threshold E^{BS} ≈ 500 MeV leads to a slight overestimation of their differential flux. However, since the contribution of e^{±} of such energies to the ionisation rate is practically negligible, CSDA can be employed. Generally, when the stopping range is comparable to (or larger than) the local column density N, the resulting spectrum of electrons and positrons is given by the convolution of the source function, (40)
where the factor 1∕2 accounts for electrons and positrons propagating in two directions and the initial energy E_{0} > E at N_{0} is related to N by (41)
If the range is small, N_{0} − N≪ N, the spectrum is localised, (42)
A check a posteriori of the energy spectra calculated with Eq. (40) for N ≳ 10^{24} cm^{−2} shows that they are accurately reproduced by Eq. (42).
5.3 Differential fluxes of secondaries
Following the diagram sketched in Fig. 5, we use the following algorithm to compute the differential fluxes of secondary CR species: In the first step, we obtain the photon flux produced by neutral pion decay and BS, the latter being generated by electrons due to primary ionisation and by electrons and positrons from charged pion decay. Next, we calculate the contribution to the electron and positron flux due to pair production by photons. Then, we employ an iterative procedure for photons, electrons, and positrons until convergence. In Fig. 7 we present the photon, electron, and positron differential fluxes (spectra) computed for typical disc (lineofsight) surface densities:
Photonspectrum: At relatively low densities, Σ ≲1 g cm^{−2}, the lowenergy part of the spectrum, below ≈0.1 MeV, is dominated by BS of secondary electrons created in primary ionisation BS(e^{−}); at energies in the range ≈0.1–100 MeV, additional BS due to electrons and positrons created by charged pion decay and pair production, BS(e^{±}), becomes important; and above ≈100 MeV, the spectrum is mostly determined by neutral pion decay, π^{0}. When Σ exceeds ≈ 100 g cm^{−2}, the spectrum is completely due to BS(e^{±}).
Electron and positron spectrum: At surface densities ≈ 1 g cm^{−2}, the spectrum below ≈100 keV is dominated by secondary electrons due to ionisation by CR protons, then CR electrons dominate up to ≈ 10 GeV, and for higher energies the contribution of electronpositron pairs becomes the most abundant component. Above ≈ 100 g cm^{−2}, the contribution of CR electrons becomes rapidly negligible; the spectrum below ≈ 1 MeV is dominated by secondary electrons produced by CR protons and by electronpositron pairs at higher energies. For Σ ≳1000 g cm^{−2}, the spectrum is entirely made of electronpositron pairs. This latter fact, along with the dominance of BS(e^{±}) in the photon spectrum, has a decisive influence on the behaviour of the ionisation rate at large Σ, as discussed in Sect. 6.
Fig. 7 Differential fluxes (energy spectra) j_{k}(E) of photons (upper row) and electrons and positrons (lower row), plotted for four characteristic values of the surface density Σ. Each plot shows partial contributions (coloured lines) to the total differential flux (black dashed line). For photons, the contributions are due to BS from electrons (e^{−}) produced by CRs (orange line) and due to π^{0} decay (blue line) and BS from electrons and positrons (e^{±}) created by pair production and charged pion decay (green line). For electrons, the contributions are due to CR electrons (solid red line), secondary electrons (e^{−}) produced by CRs (solid and dashed orange lines), and electrons and positrons (e^{±}) generated by pair production (green solid line) and decay of charged pions (dashed and dotted green lines). 
6 Ionisation at high column densities
The total ionisation rate of molecular hydrogen due to primary and secondary CR species k (primary protons and heavier nuclei, primary electrons, electronpositron pairs, and photons) is (43)
where is the ionisation cross section of H_{2} by species k and eV is the ionisation potential of H_{2}. For protons we adopt the ionisation cross sections by Rudd (1988) and Krause et al. (2015), for electrons we use results by Kim et al. (2000), while for positrons we combine the nonrelativistic cross section by Knudsen et al. (1990) with the relativistic expression for electrons. The effect of ionisation by secondary electrons produced by species k is described by a multiplicity factor, Φ_{k}(E), which is the average number of such ionisation events, (44)
Here, eV is the average energy lost by an electron per ionisation event (Glassgold et al. 2012) and is the energy loss function for species k due to ionisationby H_{2}. Since in a broad energy range (see e.g. Appendix E for protons), the multiplicity factor Φ_{k} (E) can be practically considered as a scale factor.
The contribution from charged CR species in Eq. (43) is almost entirely dominated by energies below ~ 1 GeV, assuming unmodulated spectra and ; the effect of CR modulation is studied in Sect. 7.2. The propagation of protons (as well as electrons and positrons) at these energies is described by CSDA. So far we set the pitch angle for such particles equal to zero (μ = 1, i.e. their velocities were assumed to be parallel to the local magnetic field), but in fact their local spectra should be averaged over μ. By performing the averaging, (45)
we notice that ⟨j_{k}(E, N)⟩ can be computed from the spectra for μ = 1.
The averaging over pitch angle is unimportant for electrons and positrons generated through the pair production and π^{±} decay: their contribution to ζ turns out to be negligible for N≲ 10^{25} cm^{−2}, where the direct ionisation by CR protons dominates (see below), while at higher N their spectra are well localised (see Eq. (42), i.e. their pitch angles play no role). Ionisation by CR electrons is also negligible, as we already pointed out in Sect. 4.2. Thus, we need to perform the averaging only for CR protons, i.e. Eq. (43) for these values should be computed with ⟨j_{p} ⟩^{6}.
Summing up the contribution of CR species yields the total production rate ζ = ∑_{k}ζ_{k} of molecular hydrogen ions, H. When performing the summation, we take into account the effect of heavy CR nuclei. Since the ionisation cross section scales as Z^{2} (see PGG09) and the pion production cross section as A^{0.79} (see Geist 1991), the proton ionisation rate is increased by a factor of 1.48 and the pair and photon ionisation rate by a factor of 1.30; this is the case assuming that heavy nuclei have the same attenuation as protons and that CRs have the same composition as the IS medium (see Table A.1).
Figure 8 shows the total ionisation rate and partial contributions from various species. For Σ below the transition surface density, Σ_{tr} ≈ 130 g cm^{−2}, the ionisation is mainly due to CR protons (and their secondary electrons), while at higher surface densities the contribution of electronpositron pairs produced by photon decay becomes progressively dominant. At Σ ≳ 600 g cm^{−2}, pairs fully determine the ionisation – their contribution is about a factor of 10 larger than that of CR protons – that is, the ionisation is no longer affected by the magnetic field and hence is controlled by the line of sight, rather than the effective column density. We note that in previous studies the CR ionisation rate was computed as a function of the lineofsight surface density (e.g. Umebayashi & Nakano 1981). As long as the ionisation rate is dominated by charged particles (Σ ≲Σ_{tr}), the relevant quantity is the effective surface density seen by CRs moving along magnetic field lines (Σ_{eff}). Depending on the magnetic field configuration (see e.g. Padovani et al. 2013b), the latter is generally larger or much larger than the lineofsight surface density (Σ_{los}).
The ionisation rate can be approximated by the following expression: (46)
where angle brackets denote averaging over all the directions from the transition surface towards a given position. An application of this formula is given in Sect. 7.1.
We note that direct photoionisation is always negligible since the photoionisation cross section rapidly decreases with increasing energy (see Appendix B.2). Furthermore, because of a (partially) diffusive transport of CR protons (see Sect. 4.1), their contribution, ζ_{p} (Σ), is described by a Gaussian curve at large Σ. In Fig. 8 we compare this (solid blue) curve with ζ_{p} (Σ) calculated with the pure CSDA approach.^{7} The results essentially coincide up to Σ ≈600 g cm^{−2}, where the contribution of electrons and positrons dominates.
Our important conclusion is that at large surface densities the ionisation is determined by electronpositron pairs and that the ionisation rate is not exponential anymore. The reason why the electron and positron ionisation dominates at large Σ can be inferred from Fig. 7 by comparing the behaviour of the photon, electron, and positron spectra. We see that the photon spectrum is entirely due to BS generated by e^{±}, whereas the electron and positron spectra are entirely due to the pairs produced by photons. This indicates that the feedback loop in Fig. 5 starts playing a crucial role in the ionisation process. Physically, this is because photons are able to propagate far from the source. Therefore at large Σ, where primary CRs are completely attenuated and ionisation is due to secondary particles, photons provide the only mechanism of efficient transport and ionisation (by generating pairs).
Our results are substantially different from those obtained by Umebayashi & Nakano (1981). They found the total ionisation rate decreases exponentially with a characteristic attenuation scale of about 115 g cm^{−2} for 100 g cm^{−2} ≲Σ ≲ 500 g cm^{−2}, and about 96 g cm^{−2} at larger surface densities. Conversely, we find a characteristic scale that continuously increases with surface density, from ≈ 112 g cm^{−2} to ≈ 285 g cm^{−2} in the range 100 g cm^{−2} ≲Σ ≲ 2100 g cm^{−2}, within an error lower than 10%. This difference is because Umebayashi & Nakano (1981) treated protonproton collisions above the threshold for pion production as catastrophic losses, and described highenergy Compton scattering with the CSDA approach. In addition, we consider the presence of heavy elements both in the IS CR flux and in the target medium, and perform the pitchangleaveraging for CR protons in the CSDA regime^{8}.
For computational purposes (e.g. numerical simulations and chemical models), in Appendix F we give a polynomial fit of ζ(N) valid in the range 10^{19} cm^{−2} ≤ N ≤ 10^{27} cm^{−2}. We point out that the total rate of electron production, ζ_{e}, is slightly larger than ζ because it also includes ionisation reactions of CRs with He (see Table 1 in PGG09; contributions of heavier species of the medium are negligible). For f_{He} from Table A.1, we find (47)
with an accuracy of 1% for the same range of N.
Fig. 8 Ionisation rate per H_{2} due to primary and secondary CR species, ζ, plotted vs. the surface density Σ (bottom scale) and the column density N (top scale). The black line shows the total ionisation rate. Partial contributions to ζ include ionisation due to primary CR protons and electrons (blue and red lines, respectively), ionisation due to secondary electrons created by primary CRs (orange line), and ionisation due to electrons and positrons created by charged pion decay and pair production (green line). The blue dashed line shows the proton contribution calculated with the CSDA approach. The horizontal dashed line at 1.4 × 10^{−22} s^{−1} indicates the total ionisation rate set by longlived radioactive nuclei (LLR). For comparison, we also plot the total ionisation rate per H_{2} obtained by (Umebayashi & Nakano 1981; grey dashed line). The total rate of the electron production, ζ_{e}, is approximately 1.11ζ. 
7 Applications
The ionisation fraction is a fundamental quantity for the dynamics of the IS gas, in particular during the earliest stages of star formation, from the collapse of a molecular cloud core to the formation of an accretion disc. Before the formation of a protostar, CR ionisation regulates the degree of coupling between gas and magnetic field in the densest parts of a cloud core, setting the timescale of magnetic field diffusion (see e.g. Pinto 2008) and controlling the amount of magnetic braking of collapsing rotating envelopes (Galli et al. 2006; Li et al. 2016).
In our previous studies on CR propagation (PGG09; Padovani & Galli 2011, 2013; Padovani et al. 2013a,b, 2014), we neglected the contribution of electronpositron pairs generated by photon decay and that of relativistic protons and electrons. While this approximation is appropriate for diffuse or dense molecular clouds, it becomes invalid at the high values of column densities characteristic of circumstellar discs, where CR ionisation is dominated by both relativistic and secondary particles. As shown in Sect. 6, at surface densities larger than Σ ≈130 g cm^{−2} the ionisation rate due to CR protons quickly becomes negligible, but pair production maintains ζ larger than the LLR ionisation threshold up to Σ ≈2100 g cm^{−2}.
The ionisation rate in a circumstellar disc varies considerably with radius and vertical height above the disc midplane, and is produced byseveral ionising agents, such as Galactic CRs, accelerated particles and Xrays from the central star, and short or longlived radioactive elements mixed in the gas, whose relative importance depends on the specific conditions. A careful determination of the ionisation fraction in such environment is crucial to assess the efficiency of the magnetorotational instability and the existence of the socalled dead zones with respect to mass and angular momentum transport (Gammie 1996).
In the following subsections we concentrate on the effects of Galactic CRs penetrating in discs around protostars and young stars, and limit our analysis to the disc midplane, where terrestrial planets are likely to form. Our objective is to quantify the dependence of the CR ionisation rate on the disc physical characteristics (surface density and magnetic field profiles) rather than providing an exhaustive analysis of all possible sources of ionisation. To this goal, we consider several idealised models of magnetised discs around young stars (from Shu et al. 2007) characterised by a powerlaw behaviour of the relevant properties, including the benchmark case of the (unmagnetised) minimum mass solar nebula (MMSN; Hayashi 1981). We also examine the effect of CR modulation by a stellar wind (Cleeves et al. 2013) and calculate the ionisation rate due to stellar particles (Rab et al. 2017) in the particular case of a disc around a TTauri star at a radius of 1 AU (postponing a more detailed analysis to a future study).
We stress that our results (and results by Umebayashi & Nakano 1981) do not apply to the regions of the disc dominated by MHD turbulence. In our analysis we only include the effects of largescale magnetic fields threading the disc, ignoring the scattering and diffusion of CRs due to the turbulent magnetic field. For a recent calculation of CR propagation in this regime, see RodgersLee et al. (2017).
7.1 Ionisation in magnetised circumstellar discs by Galactic CRs
We derive the CR ionisation rate in the midplane of circumstellar discs with various surface density distributions and magnetic field profiles. We choose the models defined by Shu et al. (2007), representative of low and highmass protostars (LMP and HMP), T Tauri stars (TT), and FU Orionis stars (FU Ori). These disc models are characterised by the mass of the central star, the accretion rate, and the disc age (see Table 2 in Shu et al. 2007). They are described by powerlaw surface density profiles and midplane magnetic field strength scaling with radius ϖ as (48)
with Σ_{0} = 1.36, 8.42, 33.6, and 59.5 g cm^{−2}, and B_{0} = 9.07, 8.70, 55.2, and 164 mG for TT, LMP, FU Ori, and HMP, respectively^{9}. In contrast, the unmagnetised MMSN model has a surface density of 1.7 g cm^{−2} at ϖ = 100 AU and a surface density profile proportional to ϖ^{−3∕2} (Hayashi 1981).
The effective surface density crossed by a CR propagating along a magnetic field line, inclined with respect to the disc plane, is Σ_{eff} = Σ_{disc}∕cosψ, where cos ψ = B_{z}∕B. In the models by Shu et al. (2007), the factor 1∕ cosψ is independent of the disc type and approximately equal to 3.3. As we noted in Sect. 6, below Σ_{tr} ≈ 130 g cm^{−2} the ionisation is controlled by the effective surface density measured along magnetic field lines, while above ≈ 600 g cm^{−2} becomes independent of the magnetic field configuration and hence is determined by the lineofsight surface density. Equation (46) can be used to compute the ionisation rate in the disc midplane as a function of radius ϖ. Figure 9 shows ζ for the various models. The midplane CR ionisation rate becomes dominated by LLRs inside ϖ ≈ 0.5 AU, 0.3 AU, and 0.1 AU for MMSN, HMP, and FU Ori, respectively. Only for the most evolved TT discs and for LMP, ζ is always larger than ζ_{LLR}. We note that the typical age of a TT disc is much larger than the halflife of SLRs (e.g.^{26} Al has a halflife of 0.74 Myr; Umebayashi & Nakano 2009), whose contribution to the ionisation is therefore negligible.
The results shown in Fig. 9 illustrate the contribution of unshielded Galactic CRs to the ionisation in the disc midplane (see Sect. 7.2 for the effects of stellar modulation). Among other sources of ionisation in discs around young stars, Xrays play an important role (Igea & Glassgold 1999), but the value of the Xray ionisation rate ζ_{X} at high column densities is difficult to compute because of the limitations in the Monte Carlo scattering calculations. In practice, ζ_{X} becomes uncertain above Σ_{disc} ≈ 70 g cm^{−2} (Ercolano & Glassgold 2013). For the MMSN model, this range corresponds to radii smaller than ≈ 8 AU, where Galactic CRs – if not strongly affected by the stellar wind – would mostly dominate the ionisation in the midplane.
Fig. 9 Midplane CR ionisation rate per H_{2}, ζ, plotted against radius ϖ in the TT, LMP, FU Ori, and HMP models from Shu et al. (2007) together with the standard (unmagnetised) MMSN model from Hayashi (1981). The upper and lower borders of the shaded areas correspond to the unmodulated Galactic CR spectra and , respectively. The kinks seen at the level of ζ ≈ 2 × 10^{−18} s^{−1} occur at the transition surface density Σ_{tr} = 130 g cm^{−2} in Eq. (46). The horizontal dotted line at 1.4 × 10^{−22} s^{−1} shows the value of the ionisation rate set by LLR. 
7.2 Ionisation in discs around TTauri stars
Lowenergy Galactic CRs may be prevented from penetrating an extended heliosphere (or “T Tauriosphere”) surrounding a young star (Cleeves et al. 2013). Unfortunately, little can be said about the extent and shape of this region of CR exclusion other than scaling up the properties of the Sun’s heliosphere. Cleeves et al. (2013) suggested that the T Tauriosphere may well surround the entire disc. However, the energies of CR particles mostly responsible for the ionisation at column densities above ≈ 100 g cm^{−2} exceed a few GeV. The effect of the modulation by the stellar wind at these energies is uncertain. For TT stars, Cleeves et al. (2013) estimated values for the modulation potential ϕ at a distance of 1 AU in the range ϕ = 4.8−18 GeV, leading to a reduction of the CR flux at E = 10 GeV by a factor of ≈6 and 100, respectively (ϕ is an unknown function of distance, which could be determined from detailed magnetospheric models). Moreover, the presence of a young active star may lead to increased ionisation rate, at least in the regions closest to the star becauseof thermal ionisation (Nakano & Umebayashi 1988), particle emission in flares and/or coronal mass ejection shock waves (Reames 2013, 2015), or via locally accelerated CRs in circumstellar and jet shocks (Padovani et al. 2015, 2016, 2017).
We apply our model of CR propagation to estimate the ionisation produced at a distance of 1 AU from the protostar by two different input proton spectra: a spectrum of Galactic protons modulated by TT stellar winds (Usoskin et al. 2005; Cleeves et al. 2013), and an enhanced flux of stellar protons generated by flares in an active TT star (Feigelson et al. 2002; Rab et al. 2017)^{10}. We find significant differences with respect to previous results:
(i) Figure 10 shows the CR ionisation rate for maximum and minimum modulation by a TT stellar wind at 1 AU, corresponding to ϕ = 18 GeV and ϕ = 4.8 GeV, respectively, and labelled by “GCRs (max. modul.)” and “GCRs (min. modul.)”. For completeness, the figure also shows the ionisation rate for the Galactic CR flux modulated by an average solar wind (ϕ = 1 GeV) labelled “GCRs (Sun. av. modul.)”. To facilitate the comparison with previous studies, we do not take into account inclination of the magnetic field lines with respect to the disc plane (considered in Sect. 7.1). Compared to Cleeves et al. (2013), our result for the minimum modulation model is larger by a factor of ≈ 30 at Σ ≲ 100 g cm^{−2}, while above Σ ≈ 1200 g cm^{−2} it decreases much more abruptly. The difference is even more dramatic for the maximum modulation model. We find an ionisation rate that is larger by a factor of ≈260 at Σ ≲ 100 g cm^{−2} and is decreasing faster above Σ ≈ 1400 g cm^{−2}. The discrepancy at low surface densities is largely because of our inclusion of the process of electronpositron pair creation by photon decay and also because of the adoption of the relativistic behaviour of the ionisation cross section for protons. The faster decrease of our results at high surface densities is caused by losses due to heavy elements in the medium (see Sect. 3). It is noteworthy that, in contrast to cases of unmodulated Galactic CRs, the ionisation rate for the minimum and maximum modulation is almost entirely due to relativistic protons, propagating diffusively (see Sect. 4.1). Hence, for the modulated cases we do not perform the pitchangle averaging.
(ii) For the proton flux generated in a TT flare, labelled in Fig. 10 by “SCRs (active TT)”, our results for up to Σ ≈ 300 g cm^{−2} agree with those obtained by Rab et al. (2017) within 5%. At higher surface densities the ionisation rate computed with our model decreases slowly, since electronpositron pairs increase the ionisation by orders of magnitude. It is important to remark that is still unclear what fraction of CRs generated in a flare event can be channeled into the disc through magnetic field lines, without crossing the turbulent zone, and what part may follow open field lines perpendicular to the disc (Shu et al. 2000; Feigelson et al. 2002).
Fig. 10 Midplane CR ionisation rate per H_{2}, ζ, vs. the surface density Σ (bottom scale) and column density N (top scale) plotted for several proton spectra (at a distance of 1 AU from the central star), as shown in the inset. Galactic CRs with solar mean modulation (dotted line); minimum and maximum modulation by a TT wind (dashdotdotted and long dashed line, respectively); and stellar CRs from an active TT star (dashdotted line). The horizontal dashed line shows the ionisation rate set by LLR. For reference, black triangles on the horizontal axis indicate the values of the midplane disc density (Σ_{disc} ∕2) at 1 AU for the LMP, FU Ori, MMSN, and HMP models (134, 532, 850, and 942 g cm^{−2}, respectively). 
8 Conclusions
The main result of this paper is the characterisation of the CR ionisation rate at high column densities. In particular, we showed how the CSDA fails to describe the CR proton propagation above the energy threshold for pion production (E^{π}= 280 MeV). In fact, when a CR proton interacts with a local proton to create a pion, the energy loss of the CR proton is not small anymore and there is also a certain degree of scattering. These two effects go against the main assumptions of the CSDA, namely the infinitesimal energy loss and the conservation of pitch angles. Furthermore, we carefully described the production of secondary particles, focussing on the propagation of photons created by neutral pion decay and by secondary electrons and positrons through BS. In previous studies, photon Compton losses have been treated by assuming CSDA, but we demonstrated that for this process it is crucial to use a diffusion equation. An accurate description of photon propagation is essential, since the electron and positron fluxes depend on the photon flux.
It is important to stress the main difference between the secondary particle showers that we consider here and the CR air showers in the Earth atmosphere, where the decay length of muons or pions is comparable to the scale height of the atmosphere. A big effort has been made to explain Earth air showers observed with, for example the Auger Observatory detectors and with Imaging Atmospheric Cherenkov Telescopes such as H.E.S.S. and MAGIC. Uncertainties on the hadronic cascades are the main source of error in the determination of the composition of ultrahighenergy CRs. These errors can be reduced only through detailed air shower simulations and comparisons with Large Hadron Collider and CR data (Pierog 2017). For circumstellar discs the situation is the opposite: one can assume that muons and pions immediately decay and, hence, safely consider photons and secondary electrons and positrons as direct products of the CR interaction with the local medium.
Our main conclusions are
Interstellar CR protons and particles produced by the secondary mechanisms penetrate much farther inside a circumstellar disc with respect to what has been calculated in previous studies. As a consequence, the CR ionisation rate remains above the value set by LLRs up to the surface density of ≈2100 g cm^{−2}.
PrimaryCRs are completely attenuated at high surface densities, and (secondary) photons are then the only species that propagate and efficiently create electronpositron pairs. In turn, these pairs produce efficient ionisation and, through BS, create the next generation of photons, which leads to the ionisation feedback loop ;
The total ionisation rate ζ as a function of Σ (or N) cannot be described by an exponential attenuation law. In fact, the attenuation scale continuously increases with surface density from ≈112 g cm^{−2} to ≈285 g cm^{−2} in the range 100 g cm^{2} ≲Σ ≲2100 g cm^{2}. Our results are considerably different from the exponential attenuation by Umebayashi & Nakano (1981). The difference is because of a number of improvements in our model: (i) collisions of CR protons with energies above the pion production threshold and photon Compton scattering are treated as a diffusion process; (ii) we account for the presence of heavy elements both in the CR flux and in the target medium; and (iii) we perform the averaging over initial pitch angles of CR protons in the CSDA regime.

The ionisation rate for Σ ≲ 130 g cm^{−2} is determined by CR protons (and their secondary electrons), while for Σ ≳ 600 g cm^{−2} it is completely controlled by secondary photons that create electronpositron pairs (producing local ionisation). Therefore, ζ(Σ) is a function of the effective surface density (measured along the magnetic field lines) at Σ ≲ 130 g cm^{−2} and of the lineofsight surface density at Σ ≳ 600 g cm^{−2}, since the photon propagation is unaffected by the magnetic field.

We show that ζ(N) can be described analytically by a powerlaw dependence in the range 10^{19} cm^{−2} ≤ N ≤ 10^{25} cm^{−2}; for practical purposes, we also give a polynomial fit in the whole range 10^{19} cm^{−2} ≤ N ≤ 10^{27} cm^{−2} (where N is related to Σ by Eq. (3)). The implementation of this fitting formula in numerical simulations and astrochemical codes is straightforward.
We applied our method to the propagation of CRs in magnetised circumstellar discs around young stars (Shu et al. 2007) where the ionisation fraction (which depends on the CR ionisation rate) is a key parameter that controls the coupling of the gas to disc magnetic field, the efficiency of the MRI instability, and the occurrence of dead zones. Our results can be easily incorporated in disc models together with the effects of other sources of ionisation (most importantly, Xrays) not considered in our analysis. However, a better understanding of the process of exclusion of Galactic CRs by stellar winds is needed for disc surface densities below ≈ 150 g cm^{−2}, where the CR ionisation is largely due to protons with energies below ≈ 1 GeV (which are strongly affected by stellar modulation).
Finally we checked how the secondary CR particles, composition of the medium, and averaging over the initial pitch angles affect the ionisation rate. Two different input proton spectra were considered: an IS CR proton flux modulated by TT stellar winds (Cleeves et al. 2013) and a local stellar proton flux generated in a flare event of an active TT (Feigelson et al. 2002; Rab et al. 2017). We found as follows:

While stellar winds are able to devoid the IS spectrum of lowenergy protons (below ≈1 GeV), the highenergy part of the spectrum is responsible for the production of electronpositron pairs through photon decay. The pair ionisation (along with the adoption of the relativistic ionisation cross section for protons) keeps the ionisation rate at Σ ≲ 100 g cm^{−2} much larger than previously calculated by a factor of 30 and 260 for the minimum and maximum modulation model, respectively. Furthermore, our value of ζ(Σ) calculated for these models decreases much faster at Σ ≳ 1200 g cm^{−2} and ≳1400 g cm^{−2}, respectively, because of the larger energy losses (determined by the medium composition).

For typical ages of TT discs, the ionisation by SLRs – if initially present – is negligible. The ionisation plateau is set by LLRs, and the CR ionisation dominates up to ≈1580 g cm^{−2} and ≈1820 g cm^{−2} for the maximum and minimum modulation model, respectively.

For ionisation due to stellar particles created in a TT flare, our results are comparable to previous calculations below ≈300 g cm^{−2} within 5%, while at higher surface densities electronpositron pairs increase the ionisation rate by orders of magnitude.
In this paper we developed a model of ionisation at high densities, above a few g cm^{−2}, particularly relevant for the inner regions of collapsing clouds and circumstellar discs. We calculated dependencies ζ(Σ), representing several characteristic energy spectra of CRs. Apart from an extreme (and also poorly constrained) case of ionisation due to enhanced flux of stellar protons, the obtained dependencies can be considered as fairly universal and applicable to any relevant environment. The principal limitation of our results is that they cannot be generally used to describe ionisation in regions dominated by MHD turbulence, which may lead to a diffusive transport of CRs (essentially dependent on properties of the turbulence). We plan to systematically investigate the effect of MHD turbulence in a separate paper.
Acknowledgements
The authors wish to thank the referee, Christopher McKee, for his careful reading of the manuscript and insightful comments that considerably helped to improve the paper and Elena Amato for valuable discussions. MP acknowledges funding from the European Unions Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement No 664931. PC acknowledges support from the European Research Council (ERC; project PALs 320620).
Appendix A Energy loss functions
In this appendix we discuss the individual contributions to the energy loss functions for (primary and secondary) CR particles interacting with a local IS medium composed by a mixture of hydrogen, helium, and heavier species according to Wilms et al. (2000), see Table A.1.
Assumed composition of IS medium.
A.1 Protons
The main contribution to L_{p} at low energies is due to ionisation losses that are proportional to the atomic number Z of the target species (BetheBloch formula, see Hayakawa 1969), so that , while for molecular hydrogen . The total ionisation loss function reads (A.1)
where ε^{ion} = 2.01.
At higher energies, above a threshold E^{π} = 280 MeV, we add energy losses due to pion production, as given by Schlickeiser (2002) and Krakau & Schlickeiser (2015), (A.2)
where β = v∕c is the ratio between the proton speed and the speed of light, and the asymptotic energy E^{as} = 200 GeV. The factor is a phenomenological correction to the pion production cross section for heavier target species (Geist 1991). Pion losses become dominant for E ≳ 1 GeV, fully determining the propagation of highenergy CRs at high column densities. The total pion production loss function reads (A.3)
where ε^{π} = 2.17.
A.2 Electrons and positrons
Ionisation losses for electrons have the same correction factor as protons, , see Eq. (A.1).
BS losses dominate at E ≳ 100 MeV. We take into account that and that the differential BS cross section is proportional to Z(Z + 1); see Appendix B.3. This yields (A.4)
where ε^{BS} = 2.24.
Synchrotron losses dominate at energies above E^{syn} ≈ 1 TeV and do not depend on the composition. Following Schlickeiser (2002), is (A.5)
where we have assumed a relation between the magnetic field strength, B, and the gas number density, n, given by Crutcher (2012) (A.6)
with B_{0} ≈ 10 μG, n_{0} = 150 cm^{−3} and κ ≈ 0.5–0.7. The value of κ recommended by Crutcher (2012) is κ = 0.65, but we assume κ = 0.5 (Nakano et al. 2002; Zhao et al. 2016) to benefit from the removal of the dependence of on n.
For positrons we adopt the same total loss function as for electrons, and therefore use the same notation L_{e} for both species.
A.3 Photons
Photoionisation and pair production are catastrophic processes. Their loss functions are proportional to the corresponding cross sections σ^{PI} (see Appendix section linking B.2) and σ^{pair} (see Appendix B.3). Compton effect is a continuous loss process; its energy loss function (Eq. (4)) is determined by the Compton differential cross section, d σ^{C}∕dE_{e} (Eq. (B.18)), and the maximum kinetic energy transferred to the recoiling electron, (Eq. (B.17)).
Since Compton losses are proportional to Z, the correction due to heavy elements is the same as ionisation losses. Hence, with ε^{C} = 2.01 (see Appendix B.4).
Appendix B Cross sections
B.1 Elastic protonnucleus collisions
The differential cross section for protonproton collision in the centreofmass reference system is given in Jackson & Blatt (1950) as the sum of three terms: one for elastic (Coulomb) scattering, d σ^{E} ∕dΩ, one for nuclear scattering, dσ^{N}∕dΩ, and an interference term that can be neglected. In a normalised form, the first two terms are written
where α = e^{2}∕ℏc is the finestructure constant, β = v∕c is determined by the relative velocity of the two protons, δ_{0}(E, ϑ) is the nuclear phase shift (Breit et al. 1939; Jackson & Blatt 1950), and ϑ is the scattering angle. The differential cross section in the centreofmass reference frame is then (B.3)
where r_{p} = e^{2}∕m_{p}c^{2} is the classical proton radius. The momentum transfer cross section is written (B.4)
To account for the collisions between CR protons and target heavy nuclei, has to be multiplied by the correction factor ξ given by (B.5)
The factor A_{Z}∕(A_{Z} + 1) on the righthand side accounts for the efficiency of momentum transfer from a CR proton to a nucleus with the mass number A_{Z} (Landau & Lifshitz 1969); in the first term we take into account that for collisions with H_{2} the cross section is .
B.2 Photoionisation
The photoionisation cross section, σ^{PI}, accounting for medium composition, is written as (B.6)
The cross sections for different species on the righthand side are given by Yeh & Lindau (1985) and Yeh (1993)^{11}, and are matched to the asymptotic behaviour (Draine 2011). Being expressed in terms of α and the Bohr radius a_{0} = ℏ^{2}∕m_{e}e^{2}, the asymptotic cross section is written (B.7)
where I_{H} = 13.6 eV is the ionisation energy of atomic hydrogen (valid for energies much larger than Z^{2} I_{H}).
B.3 Bremsstrahlung and pair production
The differential cross section for BS of electrons on hydrogen atom is given by the BetheHeitler formula (B.8)
where x = E_{γ}∕(E_{e} + m_{e}c^{2}) and (B.9)
The functions ϕ_{1}(Δ) and ϕ_{2}(Δ) are tabulated in Table 2 of Blumenthal & Gould (1970). A simple analytical fit is written as (B.10)
where c_{1} = 3∕2 and c_{2} = 4∕3. This formula has the correct behaviour both for Δ ≪ 1 and Δ ≫ 1.
For heavier species, the differential cross section for BS is a factor Z(Z + 1) larger than that of atomic hydrogen. This factor comes from the fact that BS takes place in the nuclear Coulomb field and in the field of atomic electrons. Consequently, BS losses are proportional to Z(Z + 1) rather than Z^{2} (Wheeler & Lamb 1939; Hayakawa 1969). The differential BS cross section for H_{2} is a factor of 2 larger than that of atomic hydrogen (Gould 1969). Equation (B.8) holds for relativistic energies. For lower energies we used the parameterisations given by Koch & Motz (1959) and Sacher & Schönfelder (1984). We note that the differential cross section is divergent for E_{γ} → 0.
The differential cross section for pair production by a photon in the field of a nucleus is closely related to that for BS, since the Feynman diagrams are variants of one another. For H nuclei, (B.11)
where y = (E_{e} + m_{e}c^{2})∕E_{γ} and (B.12)
The differential cross section is clearly symmetric for y ↔ 1 − y. The total pair production cross section, (B.13)
has the asymptotic limit for E_{γ} →∞, (B.14)
As for BS, it holds and .
B.4 Compton scattering
The differential cross section of Compton scattering for atomic hydrogen, expressed in terms of the incident photon energy E_{γ} and scattering angle ϑ, is given by theKleinNishina formula (Hayakawa 1969) (B.15)
where r_{e} = e^{2}∕m_{e}c^{2} is the classical electron radius. Here x = E_{γ}∕(m_{e}c^{2}) and are normalised energies before and after scattering, related by (B.16)
The kinetic energy transferred to the recoiling electron is ; its maximum value, (B.17)
corresponds to ϑ = π. The differential cross section is straightforwardly derived from Eq. (B.15): substituting ϑ(E_{γ}, E_{e}) and using , which follows from Eq. (B.16), we get (B.18)
The cross section for Compton scattering for atomic hydrogen is obtained by integrating Eq. (B.15) over the solid angle, (B.19)
where is the Thomson cross section. Asymptotically, for x ≪ 1 and for x ≫ 1. We have and , so the total Compton cross section is given by (B.20)
The momentum transfer cross section for atomic hydrogen, (B.21)
has the following analytic form for x > 10^{−3}: (B.22)
For x < 10^{−3} it tends to . Similar to Eq. (B.20), we obtain .
Appendix C Matching CSDA and the diffusive regimes for protons
The two regimes of propagation must be matched at the transition energy E^{tr} (see lower panel of Fig. 3): for E ≥ E^{tr} the diffusion solution is given by Eq. (21), which yields the matching spectrum j_{p} (E^{tr}, N) for the CSDA regime operating at lower energies. For brevity, below we omit particle indices and introduce the auxiliary function R(E, E_{0}) ≡ R(E_{0}) − R(E), determined by the range function, Eq. (9). Two solutions are possible in the CSDA regime, depending on how N compares to the transition column density N_{tr} = R(0, E^{tr}):
(i) N ≤ N^{tr}: the energy range is divided into two parts, E ≤ E_{*} and E_{*} < E < E^{tr}, with E_{*} determined from R(E_{*}, E^{tr}) = N. For E ≤ E_{*} the local flux is given by the attenuated IS spectrum, (C.1)
with E_{0} from R(E, E_{0}) = N, whereas for E_{*} < E < E^{tr} it is governed by the matching spectrum, (C.2)
with ΔN = R(E, E^{tr}). We note that for E = E_{*} we have E_{0} = E^{tr} in Eq. (C.1) and ΔN = N in Eq. (C.2). Since j(E, 0) = j_{IS}(E), the solution is continuous at E = E_{*}.
(ii) N > N^{tr}: the IS spectrum is completely attenuated, so Eq. (C.2) is valid for all E < E^{tr}.
Appendix D Losses due to elastic protonnucleus collisions
Elastic collisions of CR protons with nuclei of the medium are accompanied by energy exchange. As this process is most efficient for particles of equal mass, let us examine the effect of protonproton collisions and consider for simplicity the CSDA regime governed by Eq. (12). The energy exchange leads to a sink term − σ^{pp}(E)j_{p}(E, N) (to be added to the righthand side), describing a depopulation of CR energy state E due to elastic collisions with hydrogen nuclei. There is also a source term that consists of two contributions: , due to depopulation of higher energy CR state E + E^{′} (accompanied by exchange of energy E^{′} with a hydrogen nucleus), and , due to hydrogen nuclei that acquire energy E after collisions with CR protons; it is naturally assumed that such collisions result in dissociation of molecular hydrogen. By introducing the differential cross section of protonproton collisions, d σ^{pp}∕dΔE, which is a function of CR energy E and energy exchange ΔE, we have (D.1)
while is given by the same expression with arguments (E + E^{′}, E) for the crosssection.
The inclusion of these additional terms in Eq. (12) generally results in a complicated integrodifferential equation forCR protons. These terms (negligible for nonrelativistic energies, where ionisation losses dominate) may play a role for relativistic protons. The interaction with the medium is then mostly due to nuclear scattering, which is characterised by hardspherelike cross section (see upper panel of Fig. 3). In this case d σ^{pp}∕dΔE ≈ σ^{pp}∕E, i.e. the differential cross section does not depend on ΔE and is determined by a constant σ^{pp} (equal to mb). By substituting the resulting source and sink terms in Eq. (12), we obtain the following transport equation for CR protons (μ = 1): (D.2)
An approximate solution of this equation for high energies can be factorised, (D.3)
where is a solution of (homogeneous) Eq. (12) and is an unknown effective cross section, describing the cumulative effect of elastic nuclear collisions and depending on the form of . To obtain , we notice that the loss function in the highenergy regime is dominated by the pion production and, according to Eq. (A.2), can be roughly approximated by , where σ^{π} ≈ 32 mb is the pion production cross section (neglecting a weak logarithmic energy dependence) and δ_{π} ≈ 0.3 is the energy fraction lost in a single collision. Assuming , we get (D.4)
and . Since δ_{π} σ^{π} is of the order of σ^{pp} and ν ≈ 2.7 (relativistic part of ), we conclude that the argument of the exponential in Eq. (D.3) is much smaller than that in Eq. (D.4), i.e. the contribution of elastic collisions of CR protons can be safely neglected.
Appendix E Ionisation by CR protons at low column densities
Consider CR ionisation at relatively low N, where ionisation is still the main loss mechanism. Our numerical results show that, for model and of IS proton spectra , the contribution of CR electrons to ζ(N) can be neglected at N ≳ 10^{19} cm^{−2} and N ≳ 3 × 10^{21} cm^{−2}, respectively. Therefore, starting from these column densities we are only interested in the propagation of CR protons.
The upper panel of Fig. 1 shows that ionisation dominates losses for nonrelativistic protons, and for 10^{5} eV ≲ E ≲ 5 × 10^{8} eV the loss function is very well approximated by a single powerlaw dependence, (E.1)
where A =1.77 × 10^{−10} eV cm^{2} and s ≈ 0.82 (energy is ineV). The propagation and attenuation of such protons occurs in the CSDA regime and is governed by Eq. (12). Furthermore, from Fig. 2 we infer that these energies correspond to a very broad range of column densities, 10^{19} cm^{−2} ≲ N ≲ 10^{25} cm^{−2}. By substituting Eq. (E.1) in Eq. (12), we derive the following general solution valid for these N: (E.2)
where the function Ψ(x) is determined by matching j_{p}(E, 0) with the IS spectrum . For instance, for a powerlaw spectrum we get (E.3)
With the derived local spectrum, it is straightforward to obtain ζ(N). We substitute Eq. (E.3) in Eq. (43) and notice that the cross section of ionisation by nonrelativistic protons obeys a powerlaw scaling for E ≳ 10^{5} eV, with b ≈ s, i.e. is practically independent of E and hence Φ_{p} (E) ≈ const in Eq. (43). Then integration over E yields the following dependence: (E.4)
where q = (ν + b − 1)∕(1 + s) and c_{1,2} are constants. Equation (E.4) is obtained assuming 1 + s − b > 0, and is valid as long as q > 0, i.e. for ν > 1 − b. For b ≈ s ≈ 0.82 we obtain q ≈ 0.55ν − 0.10, which is valid for ν ≳ 0.2. We note that here ν represents the nonrelativistic part of , e.g. for model (ν = 0.8) we have q ≈ 0.34. The lower bound of applicability of Eq. (E.4) is determined by the actual form of , as mentioned above, while the upper bound is N ≈ 10^{25} cm^{−2} (or Σ ≈40 g cm^{−2}).
Appendix F Polynomial fit of the CR ionisation rate at any column density
For practicalpurposes, the total CR ionisation rate (of H_{2}) can be parameterised with the following fitting formula: (F.1)
It is applicable for column densities 10^{19} cm^{−2} ≤ N ≤ 10^{27} cm^{−2} with a maximum error of 6% and an average accuracy of 2%. Table F.1 gives two sets of coefficients, c_{k}, for both models and , since at low column densities the ionisation rate depends on the lowenergy CR proton spectrum. Figure F.1 shows a loglog plot of ζ versus N for the two models.
Fig. F.1 Total CR ionisation rate ζ per H_{2} due to primary and secondary CR species plotted vs. the column density N. The horizontal dashed line at 1.4 × 10^{−22} s^{−1} indicates thetotal ionisation rate set by LLR. Models (black) and (grey) are described by Eq. (1). 
According to Eq. (46), ζ(N) is a function of the effective column density as long as Σ ≲ 130 g cm^{−2} (N ≲ 3 × 10^{25} cm^{−2}). The excess over this transition value should be calculated as the lineofsight column density.
References
 Aguilar, M., Aisa, D., Alvino, A. et al. 2014, Phys. Rev. Lett., 113, 121102 [Google Scholar]
 Aguilar, M., Aisa, D., Alpat, B. et al. 2015, Phys. Rev. Lett., 114, 171103 [CrossRef] [PubMed] [Google Scholar]
 Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [Google Scholar]
 Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., Ginzburg, V. L., & Ptuskin, V. S. 1990, in Astrophysics of Cosmic Rays, ed. V. L. Ginzburg (NorhtHolland: Amsterdam) [Google Scholar]
 Breit, G., Thaxton, H. M., & Eisenbud, L. 1939, Phys. Rev., 55, 1018 [NASA ADS] [CrossRef] [Google Scholar]
 Burlaga, L. F., Ness, N. F., Gurnett, D. A., & Kurth, W. S. 2013, ApJ, 778, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Blumenthal, G. R., & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237 [Google Scholar]
 Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002 ApJ, 572, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Cesarsky C. J., & Völk, H. J. 1978, A&A, 70, 367 [NASA ADS] [Google Scholar]
 Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Crutcher, R. M. 2012, ARA&A, 50, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Cummings, A. C., Stone, E. C., Heikkila, B. C. et al. 2016, ApJ, 831, 18 [Google Scholar]
 Dalgarno, A., Yan, M., & Liu, W. 1999, ApJS, 125, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Diehl, R., Halloin, H., Kretschmer, K. et al. 2006, A&A, 449, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press) [Google Scholar]
 Ercolano, B., & Glassgold, A. E. 2013, MNRAS, 436, 3446 [NASA ADS] [CrossRef] [Google Scholar]
 Everett, J. E., & Zweibel, E. G., 2011, ApJ, 739, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Feigelson, E. D., Garmire, G. P., & Pravdo, S. H. 2002, ApJ, 572, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Fisk, L. A., & Gloeckler, G. 2014, ApJ, 789, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Fraschetti, F., Drake, J. J., Cohen, O., & Garraffo, C. 2018, ApJ, 853, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Galli, D., Lizano, S., Shu, F. H., & Allen, A. 2006, ApJ, 647, 374 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Geist, W. M. 1991, NuPhA, 525, 149c [NASA ADS] [Google Scholar]
 Ginzburg, V. L., & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays (New York: Macmillan) [Google Scholar]
 Glassgold, A. E., Najita, J., & Igea, J. 1997, ApJ, 480, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Glassgold, A. E., Galli, D., & Padovani, M. 2012, ApJ, 756, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Gloeckler, G., & Fisk, L. A. 2015, ApJ, 806, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Gould, R. J. 1969, Phys. Rev. 185, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Gounelle, M. 2015, A&A, 582, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayakawa, S. 1969, Cosmic Ray Physics (New York: WileyInterscience) [Google Scholar]
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl, 70, 35 [CrossRef] [Google Scholar]
 Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848 [NASA ADS] [CrossRef] [Google Scholar]
 Indriolo N., & McCall B. J. 2012, ApJ, 745, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Ivlev, A. V., Padovani, M., Galli, D., & Caselli, P. 2015, ApJ, 812, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Ivlev, A. V., Dogiel, V. A., Chernyshov, D. O. et al. 2018, ApJ, 855, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Jackson, J. D., & Blatt, J. M. 1950, Rev. Mod. Phys., 22, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Kamae, T., Karlsson, N., Mizuno, T., Abe, T., & Koi, T. 2006, ApJ, 647, 692 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, Y.K., Santos, J. P., & Parente, F. 2000, Phys. Rev. A, 62, 052710 [NASA ADS] [CrossRef] [Google Scholar]
 Koch, H. W., & Motz, J. W. 1959 Rev. Mod. Phys. 31, 920 [NASA ADS] [CrossRef] [Google Scholar]
 Krakau, S., & Schlickeiser, R. 2015, ApJ, 802, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Krause, J., Morlino, G., & Gabici, S. 2015, Proceedings of the 34th International Cosmic Ray Conference (ICRC2015), 518 [Google Scholar]
 Knudsen, H., BrunNielsen, L., Charlton, M., & Poulsen, M. R. 1990, J. Phys. B 23, 3955 [NASA ADS] [CrossRef] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1969, Mechanics, 2nd edn. (Oxford: Pergamon Press) [Google Scholar]
 Lee, T., Shu, F. H., Shang, H., Glassgold, A. E., & Rehm, K. E. 1998, ApJ, 506, 898 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Z.Y., Krasnopolsky, R., & Shang, H. 2016, in From Interstellar Clouds to StarForming Galaxies: Universal Processes?, Proceedings of the International Astronomical Union, IAU Symp., 315, 118 [NASA ADS] [Google Scholar]
 MacPherson, G. J., Davis, A. M., & Zinner, E. K. 1995, Meteoritics, 30, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Maret, S., Bergin, E. A., & Lada, C. J. 2006, Nature, 442, 425 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 McKee, C. F. 1989, ApJ, 345, 782 [NASA ADS] [CrossRef] [Google Scholar]
 Morlino, G., & Gabici, S. 2015, MNRAS, 451, L100 [Google Scholar]
 Nakano, T., & Umebayashi, T. 1988, Progr. of Theor. Phys. Suppl. 96, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., & Galli, D. 2011, A&A, 530, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., & Galli, D. 2013, in Advances in Solid State Physics, Cosmic Rays in StarForming Environments, eds. D. F. Torres, & O. Reimer, 34, 61 [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 (PGG09) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2013a, A&A, 549, C3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Hennebelle, P., & Galli, D. 2013b, A&A, 560, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Galli, D., Hennebelle, P., Commerçon, B., & Joos, M. 2014, A&A, 571, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Hennebelle, P., Marcowith, A., & Ferrière, K. 2015, A&A, 582, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2017, PPCF, 59, 014002 [NASA ADS] [Google Scholar]
 Pierog, T. 2017, EPJ Web Conf., 145, 18002 [CrossRef] [Google Scholar]
 Pinto, C., Galli, D., & Bacciotti, F. 2008, A&A, 484, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Putze, A., Maurin, D., & Donato, F. 2011, A&A, 526, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rab, Ch., Güdel, M., Padovani, M. et al. 2017, A&A, 603, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reames, D. V. 2013, Space Sci. Rev., 175, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Reames, D. V. 2015, Space Sci. Rev., 194, 303 [Google Scholar]
 RodgersLee, D., Taylor, A. M., Ray, T. P., & Downes, T. P. 2017, MNRAS, 472, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Rudd, M. E. 1988, Phys. Rev. A, 38, 12 [CrossRef] [PubMed] [Google Scholar]
 Sacher, W., & Schönfelder, V. 1984, ApJ, 279, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Berlin, Heidelberg: SpringerVerlag) [CrossRef] [Google Scholar]
 Shu, F. H., Najita, J. R., Shang, H., & Li, Z.Y. 2000, in Protostars and Planets IV, eds. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: Univ. Arizona Press), 789 [Google Scholar]
 Shu, F., Galli, D., Lizano, S., Glassgold, A. E., & Diamond, P. H. 2007, ApJ, 665, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Takayanagi, M. 1973, PASJ, 25, 327 [Google Scholar]
 Tikhonov, A. N., & Samarskii, A. A. 2013, Equations of Mathematical Physics (Dover Books on Physics) [Google Scholar]
 Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617 [NASA ADS] [Google Scholar]
 Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin, I. G., AlankoHuotari, K., Kovaltsov, G. A., & Mursula, K. 2005, J. Geophys. Rev., 110, 12108 [Google Scholar]
 Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wheeler, J. A., & Lamb, W. E. 1939, Phys. Rev., 55, 858 [NASA ADS] [CrossRef] [Google Scholar]
 Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, B., Caselli, P., Li, Z.Y., et al. 2016, MNRAS, 460, 2050 [NASA ADS] [CrossRef] [Google Scholar]
 Yeh, J. J. 1993, Atomic Calculations of Photoionization Cross Sections and Asymmetry Parameters (Langhorne, PE: Gordon & Breach Science Publishers) [Google Scholar]
 Yeh, J. J., & Lindau, I. 1985, Atomic Data and Nuclear Data Tables, 32, 1 [Google Scholar]
 Ziegler, J. F., Ziegler, M. D., & Biersack, J. P. 2010, NIMPB, 268, 1818 [Google Scholar]
IS spectra of heavier nuclei (of given abundances) are also described by Eq. (1).
For a solar composition (e.g. Anders & Grevesse 1989), the mean molecular weight and resulting total ionisation rate (see Sect. 6) vary by less than 2%.
Other contributions at lower energies are described in Jackson & Blatt (1950).
According to Eq. (A.2), L_{p}(E) increases asymptotically faster than linearly, and therefore T(E) remains finite.
The effect of the averaging depends on the column density. For moderate values of N, the proton ionisation rate can be approximated by a powerlaw dependence, ζ_{p}(N) ∝ N^{−q} (see Appendix E and Fig. F.1 in Appendix F); then Eq. (45) yields ⟨ζ_{p}(N)⟩ = ζ_{p}(N)∕(1 + q). For higher column densities, the attenuation is (roughly) exponential, ζ_{p}(N) ∝ exp(−N∕N_{c}), with the characteristic scale N_{c}; then ⟨ζ_{p}(N)⟩≈ ζ_{p}(N)∕(1 + N∕N_{c}).
The CSDA curve in Fig. 8 decreases more steeply than the diffusive curve at larger Σ, since the CSDA formalism implies the existence of a certain terminal column density beyond which CR protons cannot penetrate. The latter directly follows from Eq. (9) taking into account that L_{p}(E) at high energies increases faster than linearly (Eq. (A.2)).
We note, however, that the rectilinear propagation of stellar protons assumed by Rab et al. (2017) may lead to an overestimate of the proton flux incident on the disc at 1 AU (Fraschetti et al. 2018).
See also https://vuo.elettra.eu/services/elements/WebElements.html.
All Tables
All Figures
Fig. 1 Total energy loss functions of primary and secondary CR particles k (protons, electrons and positrons, and photons), computed for a medium composition given in Table A.1 (L_{k}, thick black lines) and for atomic hydrogen (L_{k,H}, thin grey lines). Protons (upper panel): ionisation losses (short dashed lines) and pion production (dotted lines); the vertical dotted line shows the energy, E^{as}, at which the proton loss function changes its asymptotic behaviour from α = 1.28 to α = 1.08. Electrons and positrons (middle panel): ionisation losses (short dashed line), BS (longdashed line), and synchrotron losses with the uncertainty on the magnetic field strength (shaded area, see Appendix A.2); the vertical dotted line shows the transition energy, E^{syn}, between BS (α = 1) and synchrotron (α = 2) losses. Photons (lower panel): photoionisation losses (dashdotted line), Compton scattering (dotted line), and pair production (shortdashed line). 

In the text 
Fig. 2 Total range functions, R_{k}(E), of primary and secondary CR particles (thick black lines), Eq. (9). The inset shows the total range functions multiplied by Ām_{p}, to highlight the behaviour at large surface densities. For comparison, the range functions for atomic hydrogen are also plotted (thin grey lines). 

In the text 
Fig. 3 Upper panel: total momentum transfer cross section for protonnucleus collisions, σ_{MT} (E). The elastic (Coulomb) scattering dominating at lower energies crosses over to the nuclear scattering at higher energies. Lower panel: scattering parameter , Eq. (14). The vertical grey stripe indicates a continuous transition from the CSDA regime, where and the proton scattering is unimportant, to the diffusive regime, where . 

In the text 
Fig. 4 Interstellar (solid lines, labelled “IS”) and local (dashed lines) differential fluxes (or spectra) of CR protons. Only half the IS proton flux penetrates into the semiinfinite medium. Labels indicate the surface density in g cm^{−2} for the local spectra. Models of the IS proton spectra, (black) and (grey), are described by Eq. (1). 

In the text 
Fig. 5 Ionisation diagram, explaining the effect of secondary particles that are generated (directly or indirectly) by CR protons and electrons through ionisation, pion decay (π^{0}, π^{±}), BS, and pair production (pair). The secondary particles include electrons (e^{−}, due to primary and secondary ionisation), positrons and electrons (e^{±}, due to pair production and π^{±} decay; also electrons produced in secondary ionisation are included), and photons (γ, due to BS and π^{0} decay), all contributing to the respective ionisation routes. 

In the text 
Fig. 6 Upper panel: components of the cross section for photons interacting with nuclei via process l: photoionisation (σ^{PI}), Compton scattering (σ^{C}), and pair production (σ^{pair}). The momentumtransfer (MT) cross section is plotted for Compton scattering, while for catastrophic (PI and pair) processes the cross section coincides with the MT cross section. The grey lines depict the corresponding loss functions near to their crossing (in arbitrary units). The vertical grey stripes, indicating the energy intervals between the crossing points of the corresponding cross sections and loss functions, separate the diffusive and catastrophic regimes of the photon transport (see text for details). Lower panel: the relative average energy lost by a photon per collision, ⟨Δ E_{γ} ⟩∕E_{γ}, vs. the photon energy, Eq. (23). 

In the text 
Fig. 7 Differential fluxes (energy spectra) j_{k}(E) of photons (upper row) and electrons and positrons (lower row), plotted for four characteristic values of the surface density Σ. Each plot shows partial contributions (coloured lines) to the total differential flux (black dashed line). For photons, the contributions are due to BS from electrons (e^{−}) produced by CRs (orange line) and due to π^{0} decay (blue line) and BS from electrons and positrons (e^{±}) created by pair production and charged pion decay (green line). For electrons, the contributions are due to CR electrons (solid red line), secondary electrons (e^{−}) produced by CRs (solid and dashed orange lines), and electrons and positrons (e^{±}) generated by pair production (green solid line) and decay of charged pions (dashed and dotted green lines). 

In the text 
Fig. 8 Ionisation rate per H_{2} due to primary and secondary CR species, ζ, plotted vs. the surface density Σ (bottom scale) and the column density N (top scale). The black line shows the total ionisation rate. Partial contributions to ζ include ionisation due to primary CR protons and electrons (blue and red lines, respectively), ionisation due to secondary electrons created by primary CRs (orange line), and ionisation due to electrons and positrons created by charged pion decay and pair production (green line). The blue dashed line shows the proton contribution calculated with the CSDA approach. The horizontal dashed line at 1.4 × 10^{−22} s^{−1} indicates the total ionisation rate set by longlived radioactive nuclei (LLR). For comparison, we also plot the total ionisation rate per H_{2} obtained by (Umebayashi & Nakano 1981; grey dashed line). The total rate of the electron production, ζ_{e}, is approximately 1.11ζ. 

In the text 
Fig. 9 Midplane CR ionisation rate per H_{2}, ζ, plotted against radius ϖ in the TT, LMP, FU Ori, and HMP models from Shu et al. (2007) together with the standard (unmagnetised) MMSN model from Hayashi (1981). The upper and lower borders of the shaded areas correspond to the unmodulated Galactic CR spectra and , respectively. The kinks seen at the level of ζ ≈ 2 × 10^{−18} s^{−1} occur at the transition surface density Σ_{tr} = 130 g cm^{−2} in Eq. (46). The horizontal dotted line at 1.4 × 10^{−22} s^{−1} shows the value of the ionisation rate set by LLR. 

In the text 
Fig. 10 Midplane CR ionisation rate per H_{2}, ζ, vs. the surface density Σ (bottom scale) and column density N (top scale) plotted for several proton spectra (at a distance of 1 AU from the central star), as shown in the inset. Galactic CRs with solar mean modulation (dotted line); minimum and maximum modulation by a TT wind (dashdotdotted and long dashed line, respectively); and stellar CRs from an active TT star (dashdotted line). The horizontal dashed line shows the ionisation rate set by LLR. For reference, black triangles on the horizontal axis indicate the values of the midplane disc density (Σ_{disc} ∕2) at 1 AU for the LMP, FU Ori, MMSN, and HMP models (134, 532, 850, and 942 g cm^{−2}, respectively). 

In the text 
Fig. F.1 Total CR ionisation rate ζ per H_{2} due to primary and secondary CR species plotted vs. the column density N. The horizontal dashed line at 1.4 × 10^{−22} s^{−1} indicates thetotal ionisation rate set by LLR. Models (black) and (grey) are described by Eq. (1). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.