© 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²¹ cm⁻²) are X-rays, cosmic rays (CRs), and decaying of radionuclides. Deep in the densest parts of molecular clouds, characterised by densities n ≳ 10⁵–10⁶ cm⁻³ 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 star-forming regions, the situation is complicated by the effects of stellar X-rays (Glassgold et al. 1997), the possible exclusion of low-energy 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 X-rays to ionise the circumstellar gas depends on total fluxes and hardness of the spectra. For example, hard X-rays of energies 1, 6 and 20 keV are absorbed by column densities of 2.5 × 10²², 4.5 × 10²⁴ and 1.3 × 10²⁵ cm⁻², respectively (Igea & Glassgold 1999).

In the shielded regions near to the disc mid-plane, where X-rays 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, short-lived radionuclides (SLRs; mostly²⁶Al with half-life 7.4 × 10⁵ yr) contribute comparatively more than long-lived radionuclides (LLRs; mostly⁴⁰K with half-life 1.3 × 10⁹ yr), but decay faster. Assuming the²⁶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⁻¹⁹ s⁻¹ in the mid-plane of a disc for ≈1 Myr, while LLRslead to ionisation rates of only ≈ 10⁻²² s⁻¹, but on timescales much longer than the disc lifetimes (Umebayashi & Nakano 2009; Cleeves et al. 2013). However, since the average Galactic abundance of²⁶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²⁶Al-rich environment, as in the case of the Sun, is relatively low (Gounelle 2015).

The propagation of low-energy CRs at low column densities, which is characteristic of the diffuse envelopes of molecular clouds, is mostly determined by resonant scattering on self-generated 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³ −10⁴ cm⁻³ such waves cannot affect the CR transport because they are completely damped from efficient ion-neutral 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 large-scale 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 line-of-sight column density at that point (Padovani & Galli 2011, 2013; Padovani et al. 2013b). The decrease of the CR ionisation rate follows a power-law behaviour as function of the effective column density, N, in the range 10²⁰ −10²⁴ cm⁻², corresponding to effective surface densities, Σ, below a few g cm⁻² (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⁻² (Umebayashi & Nakano 1981).

The main goal of this paper is to determine the attenuation of CRs at moderate-to-high surface densities (larger than a few g cm⁻²) 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⁻² the CR ionisation rate rapidly becomes dominated by electron-positron 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 line-of-sight (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⁻¹, are well constrained by various sets of observations from ground-based to balloon and satellite detectors (e.g. Aguilar et al. 2014, 2015). On the other hand, low-energy 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⁻¹ is provided by the most recent Voyager 1 observations (Cummings et al. 2016), down to ≈ 1 MeV nuc⁻¹ 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¹, as described in Ivlev et al. (2015), $$j_k^{\text{IS}}(E)=C\frac{E^a}{{(E+E_0)}^b} {\text{eV}}^{-1} {\text{s}}^{-1} {\text{cm}}^{-2} {\text{sr}}^{-1},$$(1)

where we slightly modify the values of the parameters E₀ 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${}_3^{+}$ emission (Indriolo & McCall 2012). For this reason, we consider two different models for the CR proton spectrum: a “low” spectrum, $\mathcal{L}$, obtained by extrapolating the Voyager 1 data (a = 0.1), and a “high”spectrum, $\mathcal{H}$ (a = −0.8). The resulting ionisation rates and their comparison with observational data are discussed in Ivlev et al. (2015). The $\mathcal{L}$ and $\mathcal{H}$ 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 semi-infinite medium and, hence, assume that half of IS CRs (with the energy spectra described in Sect. 2 and isotropic pitch-angle 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 $$\overline{A}={\sum }_Z{A}_Z{f}_Z=2.35,$$(2)

where A_Z and f_Z are the mass number and abundance with respect to the total number of particles, respectively². The column density is related to the surface density, Σ = Ām_pN, where m_p is the proton mass. Numerically, the relation between the N and Σ is given by $${\frac{N}{\text{cm}}}^{-2}=2.55\times {10}^{23}{\frac{\Sigma }{\text{g} \text{cm}}}^{-2}.$$(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, $L_k^l$, 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 $$L_k^l(E)={\displaystyle {\int }_0^{E^{\text{max}}}{E}^{\prime }}\frac{\text{d}{\sigma }_k^l(E,E^{\prime })}{\text{d}{E}^{\prime }}\text{d}{E}^{\prime },$$(4)

where $\text{d}{\sigma }_k^l/\text{d}{E}^{\prime }$ 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 $$L_k^l(E)=E{\sigma }_k^l(E),$$(5)

where ${\sigma }_k^l$ is the cross section of the process. Where possible, we express the total energy loss function $L_k={\sum }_l{L}_k^l$ 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 ($L_p^{\text{ion}}$) and losses due to pion production above the threshold energy E^π = 280 MeV ($L_p^{\pi }$). Therefore, $$L_p(E)={\epsilon }^{\text{ion}}{L}_{p,\text{H}}^{\text{ion}}(E)+{\epsilon }^{\pi }{L}_{p,\text{H}}^{\pi }(E),$$(6)

where ε^ion = 2.01 and ε^π = 2.17 account for the presence of heavy elements in the target medium. The two terms on the right-hand 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 ($L_e^{\text{ion}}$), BS above ≈100 MeV ($L_e^{\text{BS}}$), and synchrotron above E^syn ≈ 1 TeV ($L_e^{\text{syn}}$). Then, L_e is given by $$L_e(E)={\epsilon }^{\text{ion}}{L}_{e,\text{H}}^{\text{ion}}(E)+{\epsilon }^{\text{BS}}{L}_{e,\text{H}}^{\text{BS}}(E)+L_e^{\text{syn}}(E),$$(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 electron-positron 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 ($L_{\gamma }^{\text{PI}}$), Compton scattering ($L_{\gamma }^{\text{C}}$), and pair production losses ($L_{\gamma }^{\text{pair}}$), $$L_{\gamma }(E)=L_{\gamma }^{\text{PI}}(E)+{\epsilon }^{\text{C}}{L}_{\gamma ,\text{H}}^{\text{C}}+{\epsilon }^{\text{pair}}{L}_{\gamma ,\text{H}}^{\text{pair}},$$(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 $L_{\gamma }^{\text{PI}}$ is provided by the K-shell photoionisation of heavy species (see e.g. Draine 2011). Detailed expressions for the three terms on the right-hand 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³.

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 K-shell ionisation of heavy species of the medium.

3.5 CSDA and catastrophic regimes

The continuous slowing-down approximation (hereafter CSDA; Takayanagi 1973) has been used so far to study the propagation of low-energy CRs in molecular clouds. It is based on two assumptions: (i) that the energy losses are continuous, and (ii) that the pitch-angle 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 so-called range functions $$R_k(E)={\displaystyle {\int }_0^E\frac{\text{d}{E}^{\prime }}{L_k(E^{\prime })}}.$$(9)

In the CSDA framework, the kinetic energy of a CR particle decreases from E₀ to E after traversing a column density $$N=R_k(E_0)-R_k(E).$$(10)

The local spectrum at that N and energy E is then related to the IS spectrum at N = 0 and energy E₀ via (see PGG09) $$L_k(E)j_k(E,N/\mu )=\frac{1}{2}{L}_k(E_0)j_k^{\text{IS}}(E_0),$$(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 semi-infinite medium. This relation directly follows from a solution of the transport equation for the CSDA regime (Ginzburg & Syrovatskii 1964; Berezinskii et al. 1990), $$\mu \frac{\partial j_k}{\partial N}-\frac{\partial }{\partial E}\left(L_k{j}_k\right)=0,$$(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 σ_kj_k on the left-hand side. In the following, we discuss the effect of catastrophic losses on propagation of high-energy 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 non-negligible 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⁻²) is modified by nuclear forces (σ_MT ∝ E⁻¹). Finally, above ≈1 GeV the momentum transfer cross section becomes energy independent, which is a manifestation of the hard sphere-like scattering⁴. The total momentum transfer cross section can be written as a function of the proton-proton momentum transfer cross section (Appendix B.1) as $${\sigma }_{\text{MT}}=\xi {\sigma }_{\text{MT}}^{\text{pp}},$$(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, $$R_{\text{sc}}(E)={\displaystyle {\int }_0^E\frac{{\sigma }_{\text{MT}}(E^{\prime })}{L_p(E^{\prime })}}\text{d}{E}^{\prime },$$(14)

which is the integral ratio of the characteristic stopping range owing to energy losses, ${\displaystyle {\int }_0^E\text{d}}{E}^{\prime }/L_p$, to the characteristic column density resulting in strong scattering, ~1∕σ_MT. The CSDA is a good approximation as long as $R_{\text{sc}}(E)\ll 1$, 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 $R_{\text{sc}}$ versus E. The CSDA works for E ≲ 25 MeV, where $R_{\text{sc}}\ll 1$; according to Fig. 2, this corresponds to column densities N ≲ 7 × 10²² cm⁻². The scattering of CR protons becomes increasingly important at higher energies, and the diffusive regime operates above ≈ 1 GeV, where $R_{\text{sc}}>1$, corresponding to N ≳ 2 × 10²⁵ cm⁻².

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²² cm⁻² ≲ N ≲ 2 × 10²⁵ cm⁻².

We nowobtain the solution for the diffusive regime, assuming continuous losses. The steady-state density of CR protons in the diffusiveregime is governed by the following equation (Ginzburg & Syrovatskii 1964), $$D_p\frac{{\partial }^2{N}_p}{\partial {\mathcal{l}}^2}=\frac{\partial }{\partial E}\left(\frac{\text{d}{E}_p}{\text{d}t}{N}_p\right).$$(15)

Here, ℓ is the coordinate along the local magnetic field (CR path), $N_p(E,\mathcal{l})$ is the number of protons per unit volume and energy, related to j_p (E, ℓ) via $$N_p(E,\mathcal{l})=\frac{4\pi j_p(E,\mathcal{l})}{\beta c},$$(16)

and $$\frac{\text{d}{E}_p}{\text{d}t}\equiv -n\beta c{L}_p(E),$$(17)

where β is the proton velocity in unit of the speed of light, c. The diffusion coefficient is $$D_p(E)\approx \frac{\beta c}{3n{\sigma }_{\text{MT}}(E)},$$(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 time-like coordinate⁵ $$T(E)=\frac{1}{3}{\displaystyle {\int }_E^{\infty }\frac{\text{d}{E}^{\prime }}{{\sigma }_{\text{MT}}(E^{\prime })L_p(E^{\prime })}},$$(19)

and taking into account that dN∕dℓ = n, we reduce Equation (15) to $$\frac{\partial F_p}{\partial T}=\frac{{\partial }^2{F}_p}{\partial N^2},$$(20)

where $F_p=\beta L_p{j}_p$. Equation (20) must be solved with the boundary condition $F_p(T,N=0)=\beta L_p{j}_p^{\text{IS}}/2$, since only half of the IS flux penetrates into the semi-infinite 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 half-space (Tikhonov & Samarskii 2013), the solution is $$\begin{array}{lll}{j}_p(E,N)\hfill & \hfill & =\frac{N}{12\sqrt{\pi }\beta (E)L_p(E)}\hfill \\ \hfill & \hfill & \times {\displaystyle {\int }_E^{\infty }\frac{\beta (E_0)j_p^{\text{IS}}(E_0)}{{\sigma }_{\text{MT}}(E_0)T^{3/2}(E,E_0)}}\exp \left[-\frac{N^2}{4T(E,E_0)}\right]\text{d}{E}_0,\hfill \end{array}$$(21)

where T(E, E₀) ≡ T(E) − T(E₀).

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 (Lk, thick black lines) and for atomic hydrogen (Lk,H, thin grey lines). Protons (upper panel): ionisation losses (short dashed lines) and pion production (dotted lines); the vertical dotted line shows the energy, Eas, 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 (long-dashed 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, Esyn, between BS (α = 1) and synchrotron (α = 2) losses. Photons (lower panel): photoionisation losses (dash-dotted line), Compton scattering (dotted line), and pair production (short-dashed line).

Fig. 2 Total range functions, Rk(E), of primary and secondary CR particles (thick black lines), Eq. (9). The inset shows the total range functions multiplied by Āmp, 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 proton-nucleus collisions, σMT (E). The elastic (Coulomb) scattering dominating at lower energies crosses over to the nuclear scattering at higher energies. Lower panel: scattering parameter $R_{\text{sc}}(E)$, Eq. (14). The vertical grey stripe indicates a continuous transition from the CSDA regime, where $R_{\text{sc}}\ll 1$ and the proton scattering is unimportant, to the diffusive regime, where $R_{\text{sc}}\gtrsim 1$.

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 semi-infinite medium. Labels indicate the surface density in g cm−2 for the local spectra. Models of the IS proton spectra, $\mathcal{L}$ (black) and $\mathcal{H}$ (grey), are described by Eq. (1).

4.2 Cosmic-ray 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²⁵ cm⁻², i.e. ≈ 58 g cm⁻².

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⁻², 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⁻² 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. Cosmic-ray protons and electrons generate secondary electrons by primary ionisation ($p_{\text{CR}}+{\text{H}}_2\to p_{{\text{CR}}^{\prime }}+{\text{H}}_2^{+}+e^{-}$). Secondary electrons with energy larger than the H₂ ionisation potential produce further generations of ambient electrons. Cosmic-ray protons colliding with protons create charged pions; through muon decay, they become electrons and positrons (π^±→ μ^±→ e^±), producing secondary ionisations. In addition, proton-proton collisions create neutral pions decaying into photon pairs (π⁰ → 2γ). The second important source of photons is BS by electrons and positrons. One should also account for electron-positron 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 $${\frac{\langle \Delta E\rangle }{E}|}^{\text{C}}=\frac{1}{E{\sigma }^{\text{C}}(E)}{\displaystyle {\int }_0^{E_e^{\text{max}}}{E}_e}\frac{\text{d}{\sigma }^{\text{C}}(E,E_e)}{\text{d}{E}_e}\text{d}{E}_e=\frac{L_{\gamma }^{\text{C}}(E)}{E{\sigma }^{\text{C}}(E)},$$(22)

where σ^C(E) is the Compton cross section (see Eq. (B.20)). For x = E∕(m_ec²) ≫ 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, $${\frac{\langle \Delta E\rangle }{E}|}_{\gamma }=\frac{{\sigma }^{\text{PI}}+L_{\gamma }^{\text{C}}/E+{\sigma }^{\text{pair}}}{{\sigma }^{\text{PI}}+{\sigma }^{\text{C}}+{\sigma }^{\text{pair}}}.$$(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 $R_{\text{sc}}(E)$ defined by Eq. (14). We integrate the ratio of the momentum transfer cross section ${\sigma }_{\text{MT}}^{\text{C}}(E)$ (Eq. (B.21)) and the loss function for Compton scattering $L_{\gamma }^{\text{C}}(E)$ (see Sect. sectionlinking A.3), which yields $R_{\text{sc}}\approx 3\times {10}^4$ 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 momentum-transfer (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 $$\mu \frac{\partial j_{\gamma }(E,N)}{\partial N}+{\sigma }_{\gamma }(E)j_{\gamma }(E,N)=S_{\gamma }(E,N).$$(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, $$S_{\gamma }=S^{{\pi }^0}+S^{\text{BS}},$$(25)

is the sum of the source function for π⁰ decay from proton-nucleus collisions, $$S^{{\pi }^0}(E,N)=2{\epsilon }^{\pi }{\displaystyle {\int }_{E+E^{\pi }}^{\infty }{j}_p}(E_p,N)\frac{\text{d}{\sigma }_{\text{H}}^{{\pi }^0}(E_p,E)}{\text{d}E}\text{d}{E}_p,$$(26)

(each proton provides two photons), and the source function for BS $$S^{\text{BS}}(E,N)=2{\epsilon }^{\text{BS}}{\displaystyle {\int }_E^{\infty }{j}_e}(E_e,N)\frac{\text{d}{\sigma }_{\text{H}}^{\text{BS}}(E_e,E)}{\text{d}E}\text{d}{E}_e,$$(27)

where the factor 2 accounts for electrons and positrons. Here, $\text{d}{\sigma }_{\text{H}}^{{\pi }^0}/\text{d}E$ is the differential cross section for photon production by π⁰ decay (Kamae et al. 2006), and $\text{d}{\sigma }_{\text{H}}^{\text{BS}}/\text{d}E$ 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 $j_{\gamma }^{\text{cat}}$ in the catastrophic regime, $$j_{\gamma }^{\text{cat}}(E,N)\text{}=\text{}\frac{1}{2}{\displaystyle {\int }_0^{\infty }{S}_{\gamma }}(E,N^{\prime })\text{\hspace{0.17em}d}{N}^{\prime }\text{}{\displaystyle {\int }_0^1\exp }\left[-\frac{{\sigma }_{\gamma }(E)|N-N^{\prime }|}{\mu }\right]\frac{\text{d}\mu }{\mu }.$$(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, $j_{\gamma }^{\text{cat}}$ 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 $$D_{\gamma }\frac{{\partial }^2{N}_{\gamma }}{\partial {\mathcal{l}}^2}=\frac{\partial }{\partial E}\left(\frac{\text{d}{E}_{\gamma }}{\text{d}t}{N}_{\gamma }\right)+4\pi n{S}_{\gamma }(E,N),$$(29)

where the number of photons per unit volume and energy is $N_{\gamma }(E,N)=4\pi j_{\gamma }(E,N)/c$. Substituting $\text{d}{E}_{\gamma }/\text{d}t=-nc{L}_{\gamma }^{\text{C}}(E)$ and $$D_{\gamma }(E)\approx \frac{c}{3n{\sigma }_{\text{MT}}^{\text{C}}(E)},$$(30)

we reduce Eq. (29) to $$\frac{\partial F_{\gamma }}{\partial T}=\frac{{\partial }^2{F}_{\gamma }}{\partial N^2}-3{\sigma }_{\text{MT}}^{\text{C}}{L}_{\gamma }^{\text{C}}{S}_{\gamma },$$(31)

where $F_{\gamma }=L_{\gamma }^{\text{C}}{j}_{\gamma }$ is a function of column density N and time-like coordinate $$T(E)=\frac{1}{3}{\displaystyle {\int }_E^{E^{\text{tr}}}\frac{\text{d}{E}^{\prime }}{{\sigma }_{\text{MT}}^{\text{C}}(E^{\prime })L_{\gamma }^{\text{C}}(E^{\prime })}}.$$(32)

Similar to the catastrophic regime, we can set zero boundary condition, $F_{\gamma }(E,0)=0$, while the initial condition is $F_{\gamma }(0,N)=F_{\gamma }^{\text{cat}}(E^{\text{tr}},N)$. The latter condition is determined by matching the diffusive and catastrophic regimes at E = E^tr (specified below).

Using again the analogy with the non-homogeneous heat diffusion problem in a half-space (Tikhonov & Samarskii 2013), we obtain the following solution for the photon differential flux: $$j_{\gamma }(E,N)=\frac{{\mathcal{F}}_1(E,N)+{\mathcal{F}}_2(E,N)}{L_{\gamma }^{\text{C}}(E)},$$(33)

where $$\begin{array}{lll}{\mathcal{F}}_1(E,N)\hfill & =\hfill & L_{\gamma }^{\text{C}}(E^{\text{tr}}){\displaystyle {\int }_0^{\infty }{j}_{\gamma }^{\text{cat}}}(E^{\text{tr}},N^{\prime })G[N,N^{\prime },T(E)]\text{d}{N}^{\prime },\hfill \\ {\mathcal{F}}_2(E,N)\hfill & =\hfill & {\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_E^{E^{\text{tr}}}{S}_{\gamma }}}(E^{\prime },N^{\prime })\hfill \\ \hfill & \hfill & \times G[N,N^{\prime },T(E)-T(E^{\prime })]\text{d}{E}^{\prime }\text{d}{N}^{\prime },\hfill \end{array}$$

are determined by the Green’s function $$\begin{array}{lll}G(N,N^{\prime },x)\hfill & =\hfill & \frac{1}{2\sqrt{\pi x}}\{\exp \left[-\frac{{(N-N^{\prime })}^2}{4x}\right]\hfill \\ \hfill & \hfill & -\exp \left[-\frac{{(N+N^{\prime })}^2}{4x}\right]\}.\hfill \end{array}$$(36)

In Sect. 5.1, we mentioned that the catastrophic solution obtained in Sect. 5.1.1 for the high- and low-energy 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 ${\sigma }_{\text{MT}}^{\text{C}}{j}_{\gamma }\gg |\partial j_{\gamma }/\partial N|$, 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 high-energy (pair) catastrophic regime occurs at E^tr ≈ 50 MeV, where both ${\sigma }_{\text{MT}}^{\text{C}}(E)$ and $L_{\gamma }^{\text{C}}(E)$ intersect with σ^pair(E) and $L_{\gamma }^{\text{pair}}(E)$, respectively. Concerning the matching with the low-energy (PI) catastrophic regime, the intersection of ${\sigma }_{\text{MT}}^{\text{C}}(E)$ and σ^PI (E) is at E ≈ 6 keV, whereas $L_{\gamma }^{\text{C}}(E)$ and $L_{\gamma }^{\text{pair}}(E)$ 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⁻¹ for Compton scattering, while for photoionisation ⟨ΔE⟩∕E|^PI = 1 by definition. Hence, the low-energy 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², so that the electron and positron energy is 0 ≤ E ≤ E_γ − 2m_ec². The resulting source function at a given column density for a single species (electron or positron), $$S^{\text{pair}}(E,N)={\epsilon }^{\text{pair}}{\displaystyle {\int }_{E+2m_e{c}^2}^{\infty }{j}_{\gamma }}(E_{\gamma },N)\frac{\text{d}{\sigma }_{\text{H}}^{\text{pair}}(E_{\gamma },E)}{\text{d}E}\text{\hspace{0.17em}d}{E}_{\gamma },$$(37)

is determined by the differential cross section $\text{d}{\sigma }_{\text{H}}^{\text{pair}}/\text{d}E$ (Eq. (B.11)). Second, electrons and positrons are created through decay of charged pions, generated in proton-nucleus collisions at energies above E^π. The corresponding source function is given by $$S^{{\pi }^{\pm }}(E,N)={\epsilon }^{\pi }{\displaystyle {\int }_{E+E^{\pi }}{j}_p}(E_p,N)\frac{\text{d}{\sigma }_{\text{H}}^{{\pi }^{\pm }}(E_p,E)}{\text{d}E}\text{d}{E}_p,$$(38)

where $\text{d}{\sigma }_{\text{H}}^{{\pi }^{\pm }}/\text{d}E$ 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 π⁰ production (see Appendix A.1). The total source function for electrons and positrons is then $$S_{e^{\pm }}=S^{\text{pair}}+S^{{\pi }^{\pm }}.$$(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, $$j_{e^{\pm }}(E,N)=\frac{1}{2L_e(E)}{\displaystyle {\int }_N^{\infty }{S}_{e^{\pm }}}(E_0,N_0)L_e(E_0)\text{\hspace{0.17em}d}{N}_0,$$(40)

where the factor 1∕2 accounts for electrons and positrons propagating in two directions and the initial energy E₀ > E at N₀ is related to N by $$|N_0-N|={\displaystyle {\int }_E^{E_0}\frac{\text{d}{E}^{\prime }}{L_e(E^{\prime })}}.$$(41)

If the range is small, |N₀ − N|≪ N, the spectrum is localised, $$j_{e^{\pm }}(E,N)=\frac{E}{L_e(E)}{S}_{e^{\pm }}(E,N).$$(42)

A check a posteriori of the energy spectra calculated with Eq. (40) for N ≳ 10²⁴ cm⁻² 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 (line-of-sight) surface densities:

Photonspectrum: At relatively low densities, Σ ≲1 g cm⁻², the low-energy 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, π⁰. When Σ exceeds ≈ 100 g cm⁻², the spectrum is completely due to BS(e^±).

Electron and positron spectrum: At surface densities ≈ 1 g cm⁻², 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 electron-positron pairs becomes the most abundant component. Above ≈ 100 g cm⁻², the contribution of CR electrons becomes rapidly negligible; the spectrum below ≈ 1 MeV is dominated by secondary electrons produced by CR protons and by electron-positron pairs at higher energies. For Σ ≳1000 g cm⁻², the spectrum is entirely made of electron-positron 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) jk(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, electron-positron pairs, and photons) is $${\zeta }_k(N)=4\pi {\displaystyle {\int }_{I_{{\text{H}}_2}}^{\infty }{j}_k}(E,N)[1+{\Phi }_k(E)]{\sigma }_k^{\text{ion}}(E)\text{\hspace{0.17em}d}E,$$(43)

where ${\sigma }_k^{\text{ion}}(E)$ is the ionisation cross section of H₂ by species k and $I_{{\text{H}}_2}=15.44$ eV is the ionisation potential of H₂. 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 non-relativistic 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, $${\Phi }_k(E)=\frac{L_k^{\text{ion}}(E)}{\langle E_e^{\text{ion}}\rangle {\sigma }_k^{\text{ion}}(E)}.$$(44)

Here, $\langle E_e^{\text{ion}}\rangle \approx 37$ eV is the average energy lost by an electron per ionisation event (Glassgold et al. 2012) and $L_k^{\text{ion}}(E)$ is the energy loss function for species k due to ionisationby H₂. Since $L_k^{\text{ion}}(E)\propto {\sigma }_k^{\text{ion}}(E)$ 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 $\mathcal{L}$ and $\mathcal{H}$; 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, $$\langle j_k(E,N)\rangle ={\displaystyle {\int }_0^1{j}_k}\left(E,\frac{N}{\mu }\right)\text{\hspace{0.17em}d}\mu \equiv N{\displaystyle {\int }_N^{\infty }\frac{j_k(E,N^{\prime })}{{N^{\prime }}^2}}\text{\hspace{0.17em}d}{N}^{\prime },$$(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²⁵ cm⁻², 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 ⟩⁶.

Summing up the contribution of CR species yields the total production rate ζ = ∑_kζ_k of molecular hydrogen ions, H${}_2^{+}$. When performing the summation, we take into account the effect of heavy CR nuclei. Since the ionisation cross section scales as Z² (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⁻², the ionisation is mainly due to CR protons (and their secondary electrons), while at higher surface densities the contribution of electron-positron pairs produced by photon decay becomes progressively dominant. At Σ ≳ 600 g cm⁻², 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 line-of-sight 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 line-of-sight surface density (Σ_los).

The ionisation rate can be approximated by the following expression: $$\zeta \simeq \{\begin{array}{ll}\zeta ({\Sigma }_{\text{eff}})\hfill & \text{for }{\Sigma }_{\text{eff}}\le {\Sigma }_{\text{tr}}\hfill \\ \langle \zeta ({\Sigma }_{\text{tr}}+{\Sigma }_{\text{los}})\rangle \hfill & \text{for }{\Sigma }_{\text{eff}}>{\Sigma }_{\text{tr}}\hfill \end{array},$$(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.⁷ The results essentially coincide up to Σ ≈600 g cm⁻², where the contribution of electrons and positrons dominates.

Our important conclusion is that at large surface densities the ionisation is determined by electron-positron 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 $e^{\pm }\stackrel{\text{BS}}{\to }\gamma \stackrel{\text{pair}}{\to }{e}^{\pm }$ 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⁻² for 100 g cm⁻² ≲Σ ≲ 500 g cm⁻², and about 96 g cm⁻² at larger surface densities. Conversely, we find a characteristic scale that continuously increases with surface density, from ≈ 112 g cm⁻² to ≈ 285 g cm⁻² in the range 100 g cm⁻² ≲Σ ≲ 2100 g cm⁻², within an error lower than 10%. This difference is because Umebayashi & Nakano (1981) treated proton-proton collisions above the threshold for pion production as catastrophic losses, and described high-energy 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 pitch-angle-averaging for CR protons in the CSDA regime⁸.

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¹⁹   cm⁻² ≤ N ≤ 10²⁷ cm⁻². 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 $${\zeta }_e(N)\approx 1.11\zeta (N),$$(47)

with an accuracy of 1% for the same range of N.

Fig. 8 Ionisation rate per H2 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 long-lived radioactive nuclei (LLR). For comparison, we also plot the total ionisation rate per H2 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 electron-positron 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⁻² 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⁻².

The ionisation rate in a circumstellar disc varies considerably with radius and vertical height above the disc mid-plane, and is produced byseveral ionising agents, such as Galactic CRs, accelerated particles and X-rays from the central star, and short- or long-lived 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 so-called 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 mid-plane, 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 power-law 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 T-Tauri 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 large-scale 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 Rodgers-Lee et al. (2017).

7.1 Ionisation in magnetised circumstellar discs by Galactic CRs

We derive the CR ionisation rate in the mid-plane 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 high-mass 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 power-law surface density profiles and mid-plane magnetic field strength scaling with radius ϖ as $${\Sigma }_{\text{disc}}={\Sigma }_0{\left(\frac{\varpi }{100 \text{AU}}\right)}^{-3/4},B_z=B_0{\left(\frac{\varpi }{100 \text{AU}}\right)}^{-11/8},$$(48)

with Σ₀ = 1.36, 8.42, 33.6, and 59.5 g cm⁻², and B₀ = 9.07, 8.70, 55.2, and 164 mG for TT, LMP, FU Ori, and HMP, respectively⁹. In contrast, the unmagnetised MMSN model has a surface density of 1.7 g cm⁻² 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⁻² the ionisation is controlled by the effective surface density measured along magnetic field lines, while above ≈ 600 g cm⁻² becomes independent of the magnetic field configuration and hence is determined by the line-of-sight surface density. Equation (46) can be used to compute the ionisation rate in the disc mid-plane as a function of radius ϖ. Figure 9 shows ζ for the various models. The mid-plane 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 half-life of SLRs (e.g.²⁶ Al has a half-life 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 mid-plane (see Sect. 7.2 for the effects of stellar modulation). Among other sources of ionisation in discs around young stars, X-rays play an important role (Igea & Glassgold 1999), but the value of the X-ray 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⁻² (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 mid-plane.

Fig. 9 Mid-plane CR ionisation rate per H2, ζ, 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 $\mathcal{H}$ and $\mathcal{L}$, 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 T-Tauri stars

Low-energy 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⁻² 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)¹⁰. 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⁻², while above Σ ≈ 1200 g cm⁻² 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⁻² and is decreasing faster above Σ ≈ 1400 g cm⁻². The discrepancy at low surface densities is largely because of our inclusion of the process of electron-positron 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 pitch-angle 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⁻² 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 electron-positron 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 Mid-plane CR ionisation rate per H2, ζ, 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 (dash-dot-dotted and long dashed line, respectively); and stellar CRs from an active TT star (dash-dotted 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 mid-plane 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 ultra-high-energy 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⁻².

• PrimaryCRs are completely attenuated at high surface densities, and (secondary) photons are then the only species that propagate and efficiently create electron-positron pairs. In turn, these pairs produce efficient ionisation and, through BS, create the next generation of photons, which leads to the ionisation feedback loop $\gamma \stackrel{\text{pair}}{\to }{e}^{\pm }\stackrel{\text{BS}}{\to }\gamma$;

• 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⁻² to ≈285 g cm⁻² 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⁻² is determined by CR protons (and their secondary electrons), while for Σ ≳ 600 g cm⁻² it is completely controlled by secondary photons that create electron-positron pairs (producing local ionisation). Therefore, ζ(Σ) is a function of the effective surface density (measured along the magnetic field lines) at Σ ≲ 130 g cm⁻² and of the line-of-sight surface density at Σ ≳ 600 g cm⁻², since the photon propagation is unaffected by the magnetic field.

• We show that ζ(N) can be described analytically by a power-law dependence in the range 10¹⁹ cm⁻² ≤ N ≤ 10²⁵ cm⁻²; for practical purposes, we also give a polynomial fit in the whole range 10¹⁹ cm⁻² ≤ N ≤ 10²⁷ cm⁻² (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, X-rays) 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⁻², 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 low-energy protons (below ≈1 GeV), the high-energy part of the spectrum is responsible for the production of electron-positron 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⁻² 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⁻² and ≳1400 g cm⁻², 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⁻² and ≈1820 g cm⁻² 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⁻² within 5%, while at higher surface densities electron-positron 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⁻², 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łodowska-Curie 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.

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 (Bethe-Bloch formula, see Hayakawa 1969), so that $L_{p,Z}^{\text{ion}}=Z{L}_{p,\text{H}}^{\text{ion}}$, while for molecular hydrogen $L_{p,{\text{H}}_2}^{\text{ion}}=2L_{p,\text{H}}^{\text{ion}}$. The total ionisation loss function reads $$L_p^{\text{ion}}(E)=\left(2f_{{\text{H}}_2}+{\sum }_{Z\ge 2}{f}_{Z}Z\right)L_{p,\text{H}}^{\text{ion}}(E)={\epsilon }^{\text{ion}}{L}_{p,\text{H}}^{\text{ion}}(E),$$(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), $$L_{p,Z}^{\pi }(E)\approx 2.57\times {10}^{-17}\frac{A_Z^{0.79}}{\beta }{\left(\frac{E}{\text{GeV}}\right)}^{1.28}{\left(\frac{E+E^{\text{as}}}{\text{GeV}}\right)}^{-0.2} \text{eV} {\text{cm}}^2,$$(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 $A_Z^{0.79}$ 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 high-energy CRs at high column densities. The total pion production loss function reads $$L_p^{\pi }(E)=\left(2f_{{\text{H}}_2}+{\sum }_{Z\ge 2}{f}_Z{A}_Z^{0.79}\right)L_{p,\text{H}}^{\pi }(E)={\epsilon }^{\pi }{L}_{p,\text{H}}^{\pi }(E),$$(A.3)

where ε^π = 2.17.

A.2 Electrons and positrons

Ionisation losses for electrons have the same correction factor as protons, $L_e^{\text{ion}}(E)={\epsilon }^{\text{ion}}{L}_{e,\text{H}}^{\text{ion}}(E)$, see Eq. (A.1).

BS losses dominate at E ≳ 100 MeV. We take into account that $L_{e,{\text{H}}_2}^{\text{BS}}=2L_{e,\text{H}}^{\text{BS}}$ and that the differential BS cross section is proportional to Z(Z + 1); see Appendix B.3. This yields $$L_e^{\text{BS}}(E)=\left(2f_{{\text{H}}_2}+{\sum }_{Z\ge 2}{f}_Z\frac{Z(Z+1)}{2}\right)L_{e,\text{H}}^{\text{BS}}(E)={\epsilon }^{\text{BS}}{L}_{e,\text{H}}^{\text{BS}}(E),$$(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), $L_e^{\text{syn}}(E)$ is $$L_e^{\text{syn}}(E)\approx 5.0\times {10}^{-14}{\left(\frac{E}{\text{TeV}}\right)}^2\text{\hspace{0.17em}eV} {\text{cm}}^2,$$(A.5)

where we have assumed a relation between the magnetic field strength, B, and the gas number density, n, given by Crutcher (2012) $$B=B_0{\left(\frac{n}{n_0}\right)}^{\kappa },$$(A.6)

with B₀ ≈ 10 μG, n₀ = 150 cm⁻³ 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 $L_e^{\text{syn}}$ 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, $E_e^{\text{max}}$ (Eq. (B.17)).

Since Compton losses are proportional to Z, the correction due to heavy elements is the same as ionisation losses. Hence, $\text{d}{\sigma }^{\text{C}}/\text{d}{E}_e={\epsilon }^{\text{C}}\text{d}{\sigma }_{\text{H}}^{\text{C}}/\text{d}{E}_e$ with ε^C = 2.01 (see Appendix B.4).

Appendix B Cross sections

B.1 Elastic proton-nucleus collisions

The differential cross section for proton-proton collision in the centre-of-mass 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 $$\begin{array}{lll}\frac{\text{d}{\tilde{\sigma }}^{\text{E}}(E,\vartheta )}{\text{dΩ}}\hfill & =\hfill & {\csc }^4\left(\frac{\vartheta }{2}\right)+{\sec }^4\left(\frac{\vartheta }{2}\right)\hfill \\ \hfill & \hfill & -{\csc }^2\left(\frac{\vartheta }{2}\right){\sec }^2\left(\frac{\vartheta }{2}\right)\cos \left[2\frac{\beta }{\alpha }\ln \tan \left(\frac{\vartheta }{2}\right)\right],\hfill \\ \frac{\text{d}{\tilde{\sigma }}^{\text{N}}(E,\vartheta )}{\text{dΩ}}\hfill & =\hfill & 4\frac{{\beta }^2}{{\alpha }^2}{\sin }^2{\delta }_0,\hfill \end{array}$$

where α = e²∕ℏc is the fine-structure constant, β = v∕c is determined by the relative velocity of the two protons, δ₀(E, ϑ) is the nuclear phase shift (Breit et al. 1939; Jackson & Blatt 1950), and ϑ is the scattering angle. The differential cross section in the centre-of-mass reference frame is then $$\frac{\text{d}{\sigma }^{pp}(E,\vartheta )}{\text{d}\text{Ω}}=\frac{r_p^2}{{\beta }^4}\left(\frac{\text{d}{\tilde{\sigma }}^{\text{E}}}{\text{d}\vartheta }+\frac{\text{d}{\tilde{\sigma }}^{\text{N}}}{\text{d}\vartheta }\right),$$(B.3)

where r_p = e²∕m_pc² is the classical proton radius. The momentum transfer cross section is written $${\sigma }_{\text{MT}}^{pp}(E)=2\pi {\displaystyle {\int }_0^{\pi }\frac{\text{d}{\sigma }^{pp}}{\text{d}\text{Ω}}}(1-\cos \vartheta )\sin \vartheta \text{d}\vartheta .$$(B.4)

To account for the collisions between CR protons and target heavy nuclei, ${\sigma }_{\text{MT}}^{pp}$ has to be multiplied by the correction factor ξ given by $$\xi =f_{{\text{H}}_2}+{\sum }_{Z\ge 2}{Z}^2\frac{A_Z}{A_Z+1}{f}_Z=1.49.$$(B.5)

The factor A_Z∕(A_Z + 1) on the right-hand 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₂ the cross section is $2{\sigma }_{\text{MT}}^{pp}$.

B.2 Photoionisation

The photoionisation cross section, σ^PI, accounting for medium composition, is written as $${\sigma }^{\text{PI}}(E)=f_{{\text{H}}_2}{\sigma }_{{\text{H}}_2}^{\text{PI}}(E)+{\sum }_{Z\ge 2}{f}_Z{\sigma }_Z^{\text{PI}}(E).$$(B.6)

The cross sections for different species on the right-hand side are given by Yeh & Lindau (1985) and Yeh (1993)¹¹, and are matched to the asymptotic behaviour (Draine 2011). Being expressed in terms of α and the Bohr radius a₀ = ℏ²∕m_ee², the asymptotic cross section is written $${\sigma }^{\text{PI}}(E)=\frac{2^8}{3Z^2}\alpha \pi a_0^2{\left(\frac{E}{Z^2{I}_{\text{H}}}\right)}^{-3.5},$$(B.7)

where I_H = 13.6 eV is the ionisation energy of atomic hydrogen (valid for energies much larger than Z² I_H).

B.3 Bremsstrahlung and pair production

The differential cross section for BS of electrons on hydrogen atom is given by the Bethe-Heitler formula $$\frac{\text{d}{\sigma }_{\text{H}}^{\text{BS}}(E_e,E_{\gamma })}{\text{d}{E}_{\gamma }}=\frac{\alpha r_e^2}{E_{\gamma }}\left\{[1+{\left(1-x\right)}^2]{\varphi }_1(\Delta )-\frac{2}{3}\left(1-x\right){\varphi }_2(\Delta )\right\},$$(B.8)

where x = E_γ∕(E_e + m_ec²) and $$\Delta =\left(\frac{m_e{c}^2}{4\alpha E_{\gamma }}\right)\frac{x^2}{(1-x)}.$$(B.9)

The functions ϕ₁(Δ) and ϕ₂(Δ) are tabulated in Table 2 of Blumenthal & Gould (1970). A simple analytical fit is written as $${\varphi }_{1,2}(\Delta )\approx 8\left[\ln \left(\frac{1}{2\alpha (1+\Delta )}\right)+\frac{c_{1,2}-\Delta }{1+2\Delta }\right],$$(B.10)

where c₁ = 3∕2 and c₂ = 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² (Wheeler & Lamb 1939; Hayakawa 1969). The differential BS cross section for H₂ 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, $$\frac{\text{d}{\sigma }_{\text{H}}^{\text{pair}}(E_{\gamma },E_e)}{\text{d}{E}_e}=\frac{\alpha r_e^2}{E_{\gamma }}\left\{[y^2+{(1-y)}^2]{\varphi }_1(\delta )+\frac{2}{3}y(1-y){\varphi }_2(\delta )\right\},$$(B.11)

where y = (E_e + m_ec²)∕E_γ and $$\delta =\left(\frac{m_e{c}^2}{\alpha E_{\gamma }}\right)\frac{1}{y(1-y)}.$$(B.12)

The differential cross section is clearly symmetric for y ↔ 1 − y. The total pair production cross section, $$\begin{array}{lll}{\sigma }_{\text{H}}^{\text{pair}}(E_{\gamma })\hfill & =\hfill & {\displaystyle {\int }_0^{E_{\gamma }}\frac{\text{d}{\sigma }_{\text{H}}^{\text{pair}}}{\text{d}{E}_e}}\text{d}{E}_e\hfill \\ \hfill & =\hfill & \alpha r_e^2{\displaystyle {\int }_0^1\left\{[y^2+{(1-y)}^2]{\varphi }_1(\delta )+\frac{2}{3}y(1-y){\varphi }_2(\delta )\right\}}\text{d}y,\hfill \end{array}$$(B.13)

has the asymptotic limit for E_γ →∞, $${\sigma }_{\text{H}}^{\text{pair}}(E_{\gamma })\to \alpha r_e^2\left[\frac{2}{3}{\varphi }_1(0)+\frac{1}{9}{\varphi }_2(0)\right]=20.6 \text{mb}.$$(B.14)

As for BS, it holds ${\sigma }_{{\text{H}}_2}^{\text{pair}}=2{\sigma }_{\text{H}}^{\text{pair}}$ and ${\sigma }_Z^{\text{pair}}=Z(Z+1){\sigma }_{\text{H}}^{\text{pair}}$.

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 theKlein-Nishina formula (Hayakawa 1969) $$\frac{\text{d}{\sigma }_{\text{H}}^{\text{C}}(E_{\gamma },\vartheta )}{\text{d}\text{Ω}}=\frac{1}{2}{r}_e^2{\left(\frac{x^{\prime }}{x}\right)}^2\left(\frac{x^{\prime }}{x}+\frac{x}{x^{\prime }}-{\sin }^2\vartheta \right),$$(B.15)

where r_e = e²∕m_ec² is the classical electron radius. Here x = E_γ∕(m_ec²) and $x^{\prime }={E^{\prime }}_{\gamma }/(m_e{c}^2)$ are normalised energies before and after scattering, related by $$\frac{1}{x^{\prime }}-\frac{1}{x}=1-\cos \vartheta .$$(B.16)

The kinetic energy transferred to the recoiling electron is $E_e=E_{\gamma }-{E^{\prime }}_{\gamma }$; its maximum value, $$E_e^{\text{max}}=\frac{2x^2}{1+2x}{m}_e{c}^2,$$(B.17)

corresponds to ϑ = π. The differential cross section is straightforwardly derived from Eq. (B.15): substituting ϑ(E_γ, E_e) and using $\text{d}\cos \vartheta /\text{d}{E}_e=-{(E_{\gamma }-E_e)}^{-2}$, which follows from Eq. (B.16), we get $$\frac{\text{d}{\sigma }_{\text{H}}^{\text{C}}(E_{\gamma },E_e)}{\text{d}{E}_e}=\frac{2\pi }{{(E_{\gamma }-E_e)}^2}\frac{\text{d}{\sigma }_{\text{H}}^{\text{C}}(E_{\gamma },E_e)}{\text{d}\text{Ω}}.$$(B.18)

The cross section for Compton scattering for atomic hydrogen is obtained by integrating Eq. (B.15) over the solid angle, $$\begin{array}{lll}{\sigma }_{\text{H}}^{\text{C}}(E_{\gamma })\hfill & =\hfill & \frac{3}{4}{\sigma }_{\text{T}}\{\frac{1+x}{x^2}\left[\frac{2(1+x)}{1+2x}-\frac{\ln (1+2x)}{x}\right]\hfill \\ \hfill & \hfill & +\frac{\ln (1+2x)}{2x}-\frac{1+3x}{{(1+2x)}^2}\},\hfill \end{array}$$(B.19)

where ${\sigma }_{\text{T}}=\frac{8}{3}\pi r_e^2$ is the Thomson cross section. Asymptotically, ${\sigma }_{\text{H}}^{\text{C}}\approx {\sigma }_{\text{T}}$ for x ≪ 1 and ${\sigma }_{\text{H}}^{\text{C}}\approx \frac{3}{8}{\sigma }_{\text{T}}\ln (2x)/x$ for x ≫ 1. We have ${\sigma }_{{\text{H}}_2}^{\text{C}}=2{\sigma }_{\text{H}}^{\text{C}}$ and ${\sigma }_Z^{\text{C}}\propto Z$, so the total Compton cross section is given by $${\sigma }^{\text{C}}=\left(2f_{{\text{H}}_2}+{\sum }_{Z\ge 2}{f}_{Z}Z\right){\sigma }_{\text{H}}^{\text{C}}=2.01{\sigma }_{\text{H}}^{\text{C}}.$$(B.20)

The momentum transfer cross section for atomic hydrogen, $${\sigma }_{\text{H},\text{MT}}^{\text{C}}(E)=2\pi {\displaystyle {\int }_0^{\pi }\frac{\text{d}{\sigma }_{\text{H}}^{\text{C}}}{\text{d}\text{Ω}}}(1-\cos \vartheta )\sin \vartheta \text{d}\vartheta$$(B.21)

has the following analytic form for x > 10⁻³: $$\begin{array}{lll}{\sigma }_{\text{H},\text{MT}}^{\text{C}}(E_{\gamma })\hfill & =\hfill & \frac{3}{4}{\sigma }_{\text{T}}[\frac{2}{{(1+2x)}^2}+\frac{2x-\ln (1+2x)}{x^2}\hfill \\ \hfill & \hfill & -\frac{2x(3+x)-(3+4x)\ln (1+2x)}{x^4}].\hfill \end{array}$$(B.22)

For x < 10⁻³ it tends to ${\sigma }_{\text{H},\text{MT}}^{\text{C}}=2{\sigma }_{\text{T}}$. Similar to Eq. (B.20), we obtain ${\sigma }_{\text{MT}}^{\text{C}}\approx 2{\sigma }_{\text{H},\text{MT}}^{\text{C}}$.

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₀) ≡ R(E₀) − 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, $$j(E,N)=j_{\text{IS}}(E_0)\frac{L(E_0)}{L(E)},$$(C.1)

with E₀ from R(E, E₀) = N, whereas for E_* < E < E^tr it is governed by the matching spectrum, $$j(E,N)=j(E^{\text{tr}},N-\Delta N)\frac{L(E^{\text{tr}})}{L(E)},$$(C.2)

with ΔN = R(E, E^tr). We note that for E = E_* we have E₀ = 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 proton-nucleus 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 proton-proton 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 right-hand 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: $S_p^{(1)}$, due to depopulation of higher energy CR state E + E^′ (accompanied by exchange of energy E^′ with a hydrogen nucleus), and $S_p^{(2)}$, 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 proton-proton collisions, d σ^pp∕dΔE, which is a function of CR energy E and energy exchange ΔE, we have $$S_p^{(1)}(E,N)={\displaystyle {\int }_0^{\infty }\frac{\text{d}{\sigma }^{pp}(E+E^{\prime },E^{\prime })}{\text{d}\Delta E}}{j}_p(E+E^{\prime },N)\text{d}{E}^{\prime },$$(D.1)

while $S_p^{(2)}(E,N)$ 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 integro-differential equation forCR protons. These terms (negligible for non-relativistic 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 hard-sphere-like 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 ${\sigma }_{\text{MT}}^{pp}\approx 3$ mb). By substituting the resulting source and sink terms in Eq. (12), we obtain the following transport equation for CR protons (μ = 1): $$\frac{\partial j_p}{\partial N}-\frac{\partial }{dE}\left(L_p{j}_p\right)=2{\sigma }^{pp}{\displaystyle {\int }_0^{\infty }\frac{j_p(E+E^{\prime })}{E+E^{\prime }}}\text{d}{E}^{\prime }-{\sigma }^{pp}{j}_p.$$(D.2)

An approximate solution of this equation for high energies can be factorised, $$j_p(E,N)\approx e^{-{\sigma }_{\text{eff}}^{pp}N}{j}_{p^{\prime }}(E,N),$$(D.3)

where $j_{p^{\prime }}(E,N)$ is a solution of (homogeneous) Eq. (12) and ${\sigma }_{\text{eff}}^{pp}$ is an unknown effective cross section, describing the cumulative effect of elastic nuclear collisions and depending on the form of $j_p^{\text{IS}}(E)$. To obtain ${\sigma }_{\text{eff}}^{pp}$, we notice that the loss function in the high-energy regime is dominated by the pion production and, according to Eq. (A.2), can be roughly approximated by $L_p^{\pi }(E)\approx {\delta }_{\pi }{\sigma }^{\pi }E$, 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 $j_p^{\text{IS}}\propto E^{-\nu }$, we get $$j^{\prime }(E,N)=e^{-(\nu -1){\delta }_{\pi }{\sigma }^{\pi }N}{j}_p^{\text{IS}}(E).$$(D.4)