© ESO, 2014

1. Introduction

Supernova remnants (SNRs) are the main candidates for accelerating Galactic cosmic rays (CRs) (see, e.g., Baade & Zwicky 1934; Ackermann et al. 2013). Particularly interesting are SNRs associated with molecular clouds (MCs) that are illuminated by CR protons accelerated at the SNR shock front, leading to the emission of gamma rays. Currently, there are about ten of these objects known to emit photons at GeV energies with typical values of 100 MeV ≲ E_γ ≲ 10 TeV detected with Fermi/LAT (see, e.g., Abdo et al. 2009, 2010a) and H.E.S.S. (see, e.g., Aharonian et al. 2008). Generally, there are three processes that might lead to the production of these gamma rays: Bremsstrahlung or inverse Compton scattering from CR electrons, or the decay of neutral pions that are formed via inelastic scattering of CR protons on ambient hydrogen. If only the detected gamma-ray spectrum is available, the processes causing the gamma radiation cannot be determined definitely from these observations since they neither allow distinguishing a leptonic from a hadronic (proton) scenario, nor do they help telling apart bremsstrahlung from inverse Compton scattering.

Although sometimes estimates for the magnetic field strength or the seed photon field for inverse Compton scattering, or deductions of average hydrogen densities from observations can be made, there is still plenty of room for interpretation regarding the origin of the gamma radiation. One way to unambiguously recognize pion production as the source of the gamma rays is the detection of neutrinos formed by the decay of charged pions, which are generated via inelastic scattering of CR protons on ambient hydrogen at a fixed ratio to the formation of neutral pions (Particle Data Group 2010). A first evidence for astrophysical neutrinos was recently announced by the IceCube collaboration (IceCube Collaboration 2013). But at this point, the observed flux of extraterrestrial neutrinos is diffuse, and the statistics does not allow statements on possible point sources. Another way to indirectly explore whether the origin of the detected gamma radiation is of hadronic nature is given by the idea that if the gamma rays are induced by high-energy CR protons (E_p ≳ 1 GeV), one also expects CR protons of lower energy, which, because of the threshold energy for pion production, do not contribute to the gamma flux (Becker et al. 2011). These low-energy protons are quite efficient in ionizing the MCs, even more efficient than CR electrons (Padovani et al. 2009), because of two effects. First, the ionization cross-section for electrons is lower than the one for protons. Second, the lower rest mass of the electrons causes them (a) to be deflected by electromagnetic fields to a higher degree; and (b) to lose their energy more rapidly than CR protons in interactions with matter. Therefore, the CR electrons cannot penetrate an MC deep enough to contribute significantly to the overall ionization rate in its interior. Hence, among all the scenarios proposed to explain the observed gamma radiation, only the hadronic one would imply effective ionization of MCs, particularly in their inner regions. A correlation study of both high-energy gamma rays and low-energy ionization signatures in spatial coincidence might therefore provide a strong support for a hadronic scenario for the generation of gamma rays.

To be able to attribute ionization signatures to low-energy CR protons, other possible ionization sources need to be ruled out, in particular UV photons and X-rays. UV photons are typically absorbed quite effectively already at low column densities of traversed matter, therefore, they play a crucial role in the total ionization only at the outer regions of the clouds (Tielens 2010). X-rays, however, can penetrate the clouds considerably deeper and are also quite effective in their ionization. Consequently, photoionization is the process inside MCs that is required to be exceeded by CR-induced ionization to enable one to attribute detectable ionization signatures to CRs (here, the term CR-induced ionization denotes ionization induced only by primary CR protons and the corresponding secondary electrons). These signatures can be distinguished in terms of their variation in the column-density dependence in the cloud. In this study, the ionization rates induced by both processes are derived as functions of the penetration depth into the cloud, thus providing a crucial clue to the dominant processes that contribute to the formation of gamma rays from SNRs associated with MCs.

The paper is organized as follows: first, the relation between gamma-ray observations and a potential underlying CR proton spectrum is described in Sect. 2. In Sect. 3, the transport equation for the propagation of the CR protons into and inside the MCs is introduced, and both the general analytic solution and a solution for specific choices of the diffusion coefficient, momentum losses, and source function are given. For four specific SNR-MC systems, the position-dependent ionization rates induced by CR protons are derived, using the solution of the transport equation, in Sect. 4. Section 5 deals with the photon fluxes both at the surfaces and inside the clouds of these systems and the corresponding position-dependent photoionization rates. The ionization rates induced by photons and CR protons are compared in Sect. 6. Finally, in Sect. 7, we give a brief overview of how ionization rates affect MCs and how they can be detected. We also provide conclusions from the results from Sect. 6 and give an outlook for future work.

2. Cosmic-ray proton spectrum

The most established processes for the acceleration of CRs are Fermi acceleration (Fermi 1949) of first and second order, in particular diffusive shock acceleration (see, e.g., Axford et al. 1977; Bell 1978a,b; Jokipii 1987; Schlickeiser 1989a,b). Within specific parameter domains, these can account for typical CR proton spectra from isolated sources. Therefore, it is not possible to determine which acceleration process is actually at work at a certain source. Usually, this poses a profound problem for the construction of theoretical models for specific objects. Nevertheless, it can be circumvented when the CR proton spectrum directly at the interaction region is derived from observational data.

In this work, a hadronic scenario for the generation of gamma rays in the SNR-MC systems is considered, and, using observed gamma-ray spectra, the proton fluxes at the gamma-ray emission regions are constructed assuming that the detected gamma radiation is caused by the decay of neutral pions formed in proton-proton interactions of CR protons with ambient hydrogen and that homogeneous volumes enclosed by the SNR shock fronts are the source regions of isotropic CR proton emission and that the MC surfaces are in contact with the shock fronts. For a parametric description of these proton fluxes, their parameters are chosen in such a way that the resulting gamma-ray emission is compatible with the observations. Models for the interaction of CR protons with ambient matter are discussed in Kelner et al. (2006) and Kamae et al. (2006), leading to parametric descriptions of the CR proton injection fluxes into the cloud with constrained kinetic proton energies E_p ≳ 280 MeV. This lower boundary comes from the minimal energy required for the formation of neutral pions, which will decay into the detected gamma rays for a hydrogen target at rest. However, since only the low-energy CR protons are able to efficiently ionize the MCs, the fluxes are extrapolated toward lower energies under the assumption that they follow the same power-law behavior as in the energy ranges described by observations. The a priori unknown shape of the CR proton spectrum at low energies is a major problem, but the extrapolation performed here is appropriate for several reasons. First, it spans less than one order of magnitude in terms of the particle momentum and, therefore, it is close to the momentum domain covered by observational data. Second, there is no apparent change in the observed spectral shape of the CR protons toward lower energies. Third, observations of ionization signatures from Indriolo et al. (2009) indicate that a CR proton flux that increases toward lower proton energies is required in order to explain these detections. Thus, the uncertainty associated with the extrapolation is reasonably small compared to other uncertainties in studies of SNRs near MCs, such as the average hydrogen density of the cloud or the volume-filling factor of the volume enclosed by the SNR shock front. The propagation of the CR protons inside the MCs is modeled in terms of a transport equation containing a general injection spectrum, momentum-dependent scalar diffusion, and ionization, Coulomb and adiabatic deceleration losses. The MCs are of homogeneous hydrogen density and of arbitrary shape.

3. Transport equation

In this section, the relevant transport equation for the CR protons is presented, the components are explained and discussed, and the solution is given in analytic form. The transport equation describing the temporal evolution of primary CR protons inside an MC located in the vicinity of an SNR reads $\frac{\mathit{\partial }{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{-}\mathit{D}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{\Delta }{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial p}}(\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{·}{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)})\mathrm{=}\mathit{Q}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}\mathit{,}$(1)where n_p(r,p,t) is the differential number density of protons at the point r with momentum p at time t, D(p) is the scalar momentum-dependent diffusion coefficient, Δ is the Laplace operator, b(p) is the momentum-loss rate, and Q(r,p,t) is a general source function for a supernova shock front that accelerates CR protons. Note that this equation does not hold for highly relativistic momenta, since in that case the diffusion coefficient can no longer be described by a scalar quantity. Furthermore, convection and advection of particles as well as diffusion in momentum space are not considered. The Green’s function of this transport equation is calculated in detail in Appendix A using Laplace transformations and Duhamel’s principle (Duhamel 1838; Courant & Hilbert 2008), yielding $\begin{array}{ccc}\mathit{G}\mathrm{\left(}{r}\mathit{,p,t}\hspace{0.17em}\mathrm{\right|}\hspace{0.17em}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathrm{)}& \mathrm{=}& \frac{\mathrm{\Theta }\mathrm{(}{\mathit{p}}_{\mathrm{0}}\mathrm{-}\mathit{p}\mathrm{)}\mathit{\delta }\left(\mathit{t}\mathrm{+}{\mathrm{\int }}_{{\mathit{p}}_{\mathrm{0}}}^{\mathit{p}}\mathit{b}\mathrm{(}{\mathit{p}}^{\mathrm{\prime }}{\mathrm{)}}^{-1}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}{\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{·}{\left(\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\\ & & \end{array}$(2)where Θ(·) is the Heaviside step function and δ(·) is the Dirac distribution. The general solution n_p(r,p,t) for an arbitrary source function Q(r₀,p₀,t₀) can be obtained by convolving the fundamental solution (2) with the source function ${\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}\mathrm{=}\$\mathit{G}\mathrm{(}{r}\mathit{,p,t}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathrm{\left)}\mathit{Q}\mathrm{\right(}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{d}{\mathit{t}}_{\mathrm{0}}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{d}{\mathit{p}}_{\mathrm{0}}\mathit{.}$(3)This differential CR proton number density can also be applied to many other astrophysical situations with scalar, momentum-dependent diffusion and any type of momentum losses, such as stellar winds, CR diffusion in the interstellar medium, or gamma-ray bursts.

Two specific momentum-loss processes in SNR-MC systems are considered: Coulomb collision losses (i.e., ionization or excitation of interstellar matter) and adiabatic deceleration (i.e., attenuation of the CR number density due to the expansion of the shock front with a momentum-loss timescale independent of the CR particle momentum). The loss rate for Coulomb collisions (Lerche & Schlickeiser 1982) can be given in the approximated form ${\mathit{b}}_{\mathrm{cc}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\mathrm{5}\mathrm{\times }{\mathrm{10}}^{-19}{\mathit{Z}}^{\mathrm{2}}(\frac{{\mathit{n}}_{\mathrm{H}}}{{\mathrm{cm}}^{-3}})(\frac{\mathit{p}}{\mathit{Mc}}{)}^{-2}(\mathrm{11.3}\mathrm{+}\mathrm{2}\ln (\frac{\mathit{p}}{\mathit{Mc}}))\mathrm{eV}\mathrm{c}{\mathrm{m}}^{-1}\mathit{,}$(4)where Z is the charge number of the incoming CR particle, n_H is the hydrogen density of the medium, M is the particle rest mass, and c denotes the speed of light. Here, only protons (i.e., Z = 1 and M = m_p) of a minimum momentum of 0.15 m_pc, corresponding to a kinetic energy of 10 MeV, and a maximum momentum of 0.86 m_pc, corresponding to a kinetic energy of 280 MeV, are considered. Protons of lower momentum can be neglected since they do not penetrate the cloud deep enough in order to provide an observable contribution to the ionization rate. Protons of higher momentum can be neglected on the one hand, because of the rapid decrease of the ionization cross-section with increasing momentum (Rudd et al. 1985) and, on the other hand, because they can form pions by inelastic proton-proton scattering (the momentum of ~0.86 m_pc corresponds to the threshold kinetic energy required for protons to produce pions by proton-proton interactions), which renders them unavailable for the ionization of the MC. Imposing these limits on the CR proton momentum, the logarithmic contribution in Eq. (4) is small compared to the constant and is therefore, as an approximation, disregarded. Then, the loss rate for Coulomb collisions simplifies to ${\mathit{b}}_{\mathrm{cc}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{a}}_{\mathrm{cc}}{\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)}^{-2}$(5)with ${\mathit{a}}_{\mathrm{cc}}\mathrm{=}\mathrm{5.65}\mathrm{\times }{\mathrm{10}}^{-18}\left(\frac{{\mathit{n}}_{\mathrm{H}}}{{\mathrm{cm}}^{-3}}\right)\mathrm{eV}\mathrm{c}{\mathrm{m}}^{-1}\mathit{.}$(6)The loss rate for adiabatic deceleration at a shock velocity V can be expressed as ${\mathit{b}}_{\mathrm{ad}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\hspace{0.17em}\mathit{p}\mathrm{\nabla }\mathrm{·}{V}\mathit{.}$(7)Assuming a constant, spherical expansion of the shock front, the shock velocity reads ${V}\mathrm{=}{\mathit{V}}_{\mathrm{0}}\hspace{0.17em}{{e}}_{\mathit{R}}\mathit{,}$(8)where V₀ is a constant, and R is the shock radius. Hence, the divergence of the shock velocity becomes $\mathrm{\nabla }\mathrm{·}{V}\mathrm{=}\frac{\mathrm{1}}{{\mathit{R}}^{\mathrm{2}}}{\mathit{\partial }}_{\mathit{R}}({\mathit{R}}^{\mathrm{2}}{\mathit{V}}_{\mathrm{0}})\mathrm{=}\mathrm{2}\frac{{\mathit{V}}_{\mathrm{0}}}{\mathit{R}}\mathit{,}$(9)and the loss rate for adiabatic deceleration (7) yields ${\mathit{b}}_{\mathrm{ad}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{a}}_{\mathrm{ad}}\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)$(10)with ${\mathit{a}}_{\mathrm{ad}}\mathrm{=}\mathrm{2.085}\mathrm{\times }{\mathrm{10}}^{-2}\left(\frac{{\mathit{V}}_{\mathrm{0}}\mathit{/}\mathit{R}}{{\mathrm{s}}^{-1}}\right)\mathrm{eV}\mathrm{c}{\mathrm{m}}^{-1}\mathit{.}$(11)

Fig. 1 Comparison of the dynamical timescales for Coulomb losses, adiabatic deceleration, and diffusion with the typical values nH = 100 cm-3, tage = 10 kyr, and a diffusion length l = 2 pc.

Note that for a constant shock velocity the quantity V₀/R can be identified with the reciprocal age of the SNR, t_age, thus, providing a reasonable estimate of the divergence of the shock velocity, which is otherwise extremely difficult to obtain from observations. In Fig. 1, the dynamical timescales for both Coulomb losses (solid, black line) and adiabatic deceleration (long-dashed, blue line for a young SNR and long-short-dashed, orange line for a middle-aged SNR), as well as the diffusion timescale (short-dashed, red line), as functions of the particle momentum p, are depicted. The various timescales are given by ${\mathit{T}}_{\mathrm{cc}}\mathrm{=}\mathrm{(}\mathrm{d}\mathit{p}\mathit{/}\mathrm{d}\mathit{t}{\mathrm{)}}_{\mathrm{cc}}^{-1}\mathit{p}$, ${\mathit{T}}_{\mathrm{ad}}\mathrm{=}\mathrm{(}\mathrm{d}\mathit{p}\mathit{/}\mathrm{d}\mathit{t}{\mathrm{)}}_{\mathrm{ad}}^{-1}\mathit{p}$, and T_diff = l²/ (2D(p)), where l is the diffusion length, and the diffusion coefficient for the SNR-MC systems is D(p) = 10²⁸·(p/ (m_pc))^4/3 cm² s^-1 with a power-law index of 4/3. Despite the specific choice of this power-law index, any power-law dependence of the diffusion coefficient on the particle momentum allows for analytic evaluation of the solution of the transport equation with the momentum losses considered here. The general form of the diffusion coefficient for a Kolmogorov spectrum in a statistically isotropic random magnetic field (see, e.g., Ptuskin 2006) is D(p) ∝ vp^1/3, where v is the CR particle speed. In the non-relativistic limit, v ∝ p and, therefore, D(p) ∝ p^4/3, yielding a reasonable approximation of the diffusion coefficient for low momenta. Since the kinetic energies of the CR protons considered in this work also exceed values that justify the non-relativistic treatment, the assumption D(p) ∝ p^4/3 leads to an overestimate of the diffusion effect for particles with momenta exceeding the non-relativistic regime. This results in an overestimated ionization rate at large penetration depths, because the time a CR particle takes to reach a certain penetration depth decreases with an increasing diffusion coefficient, while at small penetration depths, the resulting ionization rate is underestimated. Apart from Kolmogorov spectra, Iroshnikov-Kraichnan spectra can be applied as well, by adapting the power-law index of the diffusion coefficient. When Goldreich-Sridhar spectra provide the adequate description, the diffusion coefficient splits into a parallel and perpendicular component relative to the mean magnetic field orientation. The values chosen for the average hydrogen densities of the MCs and for the diffusion coefficient in Fig. 1 are the same as those used in the data analysis in Sect. 4. For middle-aged SNRs, losses from Coulomb collisions for p ≲ 0.2 m_pc occur on shorter timescales than adiabatic deceleration losses. Thus, they dominate momentum losses of the CR protons in this regime, while at higher momenta, the losses from adiabatic deceleration prevail. For young SNRs, this transition is shifted toward lower momenta of p ≈ 0.1 m_pc, because the dynamical timescale for adiabatic deceleration is directly proportional to the age of an SNR, and, therefore, the momentum loss from adiabatic deceleration for those SNRs dominates for all momenta covered by this study. The diffusion timescale indicates how long a particle with given momentum p takes to diffuse 2 pc, ignoring the effect of momentum losses. Figure 1 indicates that for young objects, particles with a momentum p ≲ 0.5 m_pc lose all their kinetic energy due to adiabatic deceleration losses before reaching a penetration depth of 2 pc, while for middle-aged SNRs, this only happens for low-momentum particles with p ≲ 0.15 m_pc. The maximum penetration depth of the latter is governed by Coulomb losses.

Using both Eqs. (6) and (11), the total momentum-loss rate can be expressed as $\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{a}}_{\mathrm{ad}}\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)\mathrm{+}{\mathit{a}}_{\mathrm{cc}}{\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)}^{-2}\mathrm{·}$(12)Moreover, the source function Q(r₀,p₀,t₀) is constructed such that it reflects the geometry of the SNR-MC system and the energy dependence of the incident CR proton flux. This flux is derived from observations for four specific SNRs associated with MCs that show gamma-ray emission and for which data samples from spectral measurements in the X-ray energy range exist. Under the assumptions that the CR protons are constantly and isotropically emitted from the homogeneous region enclosed by the SNR shock front and that the MC is located directly adjacent to the emission volume, the source function is modeled in the specific form $\begin{array}{ccc}\mathit{Q}\mathrm{\left(}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}& \mathrm{=}& {\mathit{Q}}_{\mathrm{norm}}\mathrm{·}{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{\mathit{p}}_{\mathrm{0}}\mathrm{\right)}\mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathrm{-}\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}]\\ & & \mathrm{·}[\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)\mathrm{-}\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)]\\ & & \end{array}$(13)where Q_norm denotes a normalization constant and Q_p(p₀) is the spectral shape of the low-energy CR protons in terms of the particle momentum p₀, describing emission that is constant over a time interval of one second from a cubic emission volume with edge length l_c, seen by an observer located at the center of the face of the emission volume that is the boundary surface of the SNR shock front and the MC, with a coordinate system such that the positive Cartesian z-axis is normal to this face and points into the cloud. The volume of the emission region, ${\mathit{l}}_{\mathrm{c}}^{\mathrm{3}}$, is modeled as cubic for numerical feasibility and adapted to the spherical emission volume used in the modeling process of the gamma rays in Sect. 4.

Performing the convolution in Eq. (3) for the source function (13) yields the differential CR proton number density at all positions r inside the MC with z ≥ 0 at any time t ≥ 0 and for all particle momenta p ∈ [ 0.15 m_pc, 0.86 m_pc ] ≤ p₀$\begin{array}{ccc}{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{Q}}_{\mathrm{norm}}\mathrm{·}\mathrm{(}{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}\mathit{p,t}\mathrm{\right)}{)}^{\mathrm{3}}\mathrm{)}\mathrm{·}{\mathit{Q}}_{\mathit{p}}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}\mathit{p,t}\mathrm{\right)})}{\mathrm{8}\mathrm{·}\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{·}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}\mathit{p,t}\mathrm{\right)}{)}^{\mathrm{2}}}\mathrm{·}\underset{\mathit{h}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{x,y,z}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}}{\mathrm{􏽙}}\left[\sum _{\mathit{j}\mathrm{=}\mathrm{0}}^{\mathrm{1}}\mathrm{erf}\left(\frac{\mathrm{(}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{j}}\mathrm{·}{\mathit{h}}^{\mathrm{)}}\sqrt{{\mathit{a}}_{\mathrm{cc}}\mathrm{(}\mathrm{3}\mathrm{+}\mathit{k}\mathrm{)}}}{\sqrt{\mathrm{4}{\mathit{D}}_{\mathrm{0}}\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}\mathit{c}{\mathrm{[}{\mathrm{(}{p}^{\mathrm{\prime }}\mathrm{)}}^{\mathrm{3}\mathrm{+}\mathit{k}}{\hspace{0.17em}}_{\mathrm{2}}{\mathit{F}}_{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{2}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{-}\mathit{a}\mathrm{·}{\mathrm{(}{p}^{\mathrm{\prime }}\mathrm{)}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{]}}}_{p^{\mathrm{\prime }}\mathrm{=}\mathit{p}\mathit{/}\mathrm{\left(}{\mathit{m}}_{\mathrm{p}}\mathit{c}\mathrm{\right)}}^{p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}\mathit{p,t}\mathrm{\right)}}}}\right)\right]\mathit{,}\end{array}$(14)where a = a_ad/a_cc, k is the power-law index of the diffusion coefficient (for the data analysis, k = 4/3), erf(·) and ${}_{\mathrm{2}}\mathit{F}_{\mathrm{1}}\mathrm{(}\mathrm{·}\hspace{0.17em}\mathit{,}\mathrm{·}\hspace{0.17em}\mathrm{;}\mathrm{·}\hspace{0.17em}\mathrm{;}\mathrm{·}\mathrm{)}$ denote the error function and the hypergeometric function, respectively, and $\begin{array}{ccc}{p}_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}\mathit{p,t}\mathrm{\right)}\mathrm{=}({\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)}^{\mathrm{3}}\mathrm{·}\exp \left(\frac{\mathrm{3}{\mathit{a}}_{\mathrm{ad}}\mathit{t}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)\mathrm{+}\frac{{\mathit{a}}_{\mathrm{cc}}}{{\mathit{a}}_{\mathrm{ad}}}\left(\exp \left(\frac{\mathrm{3}{\mathit{a}}_{\mathrm{ad}}\mathit{t}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)\mathrm{-}\mathrm{1}\right){)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{.}& & \end{array}$For details, see Appendix B. In a next step, the undetermined normalization constant Q_norm and the momentum-dependent spectral function Q_p are modeled such that they are adapted to the CR proton fluxes, derived from gamma-ray observations, in the SNR-MC systems outside the clouds, yielding initial boundary values of the differential CR proton number density at the cloud surfaces. Then, Eq. (14) can be used to determine the CR proton fluxes inside the MCs that are necessary for calculating the ionization rates.

4. CR-induced ionization profiles

The CR-induced ionization rate of molecular hydrogen in MCs is computed for four specific objects: W49B, W44, 3C 391, and CTB 37A, following Padovani et al. (2009), ${\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}{r}\mathit{,t}\mathrm{\right)}\mathrm{=}\left(\mathrm{1}\mathrm{+}\mathit{\phi }\right){\mathrm{\int }}_{{\mathit{E}}_{\min }}^{{\mathit{E}}_{\max }}\frac{{\mathrm{d}}^{\mathrm{3}}{\mathit{N}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,}{\mathit{E}}_{\mathrm{p}}\mathit{,t}\mathrm{\right)}}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}\hspace{0.17em}\mathrm{d}\mathit{A}\hspace{0.17em}\mathrm{d}\mathit{t}}{\mathit{\sigma }}_{\mathrm{ion}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{E}}_{\mathrm{p}}\mathit{,}$(15)where φ accounts for secondary ionization (i.e., ionization events induced by secondary electrons released during primary ionization events), d³N_p/ (dE_p dA dt) is the differential CR proton flux, and ${\mathit{\sigma }}_{\mathrm{ion}}^{{\mathrm{H}}_{\mathrm{2}}}$ is the direct ionization cross-section of molecular hydrogen by protons. The primary CR proton flux can be obtained with the help of gamma-ray observations and the differential CR proton number density (14). Treating the observed gamma-ray fluxes from SNR-MC systems as isotropic emission from the regions enclosed by the SNR shock fronts, the primary CR proton fluxes in these regions are constructed, assuming that the gamma-ray emission is a result of the formation and subsequent decay of neutral pions from interactions of the CR protons with ambient hydrogen, such that the corresponding gamma-ray fluxes match the observational data. This is done by means of the proton-proton interaction model given in Kelner et al. (2006), which describes gamma-ray emission for primary energies above 100 GeV, and with the delta distribution approximation for the pion-production cross-section at lower energies (see, e.g., Mannheim & Schlickeiser 1994; Aharonian & Atoyan 2000). In order to study the transport of these CR protons into the interior regions of the clouds, one uses the CR proton fluxes outside the MCs as initial boundary values for the differential CR proton number densities (14) and determines the function Q_p(p) and the normalization constant Q_norm in Eq. (13) for the objects under consideration. To this end, the CR proton fluxes outside the MCs derived from gamma-ray observations are expressed in terms of a (dimensionless) spectral shape function Φ_p(E_p) and a normalization constant a_p, and linked to the momentum source factor of Eq. (13) and Q_norm by simple multiplication with the CR proton kinetic energy E_p and division by the particle speed used in the Kelner gamma-ray emission model, the speed of light c, yielding ${\mathit{Q}}_{\mathrm{norm}}\mathrm{·}{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{a}}_{\mathrm{p}}\mathrm{·}{\mathrm{\Phi }}_{\mathrm{p}}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)}\mathrm{·}\frac{{\mathit{E}}_{\mathrm{p}}\mathrm{\left(}\mathit{p}\mathrm{\right)}}{\mathit{c}}\mathrm{=}(\frac{{\mathrm{d}}^{\mathrm{3}}{\mathit{N}}_{\mathrm{p}}}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}\hspace{0.17em}\mathrm{d}\mathit{A}\hspace{0.17em}\mathrm{d}\mathit{t}}{)}_{\mathrm{obs}}\mathrm{·}\frac{{\mathit{E}}_{\mathrm{p}}\mathrm{\left(}\mathit{p}\mathrm{\right)}}{\mathit{c}}\mathrm{·}$(16)The spectral shape function Φ_p(E_p) that is used to model the observed gamma-ray emission as described in more detail in Schuppan et al. (2012) is usually given by a set of typical expressions for SNRs detected at gamma-ray energies:

• (a)

  a single power-law in the kinetic energy Φ_p(E_p) = (E_p/E_norm)^{− α},

• (b)

  a single power-law in the kinetic energy with an exponential cut-off Φ_p(E_p) = (E_p/E_norm)^{− α}·exp⁽ − E_p/E_c⁾,

• (c)

  a double power-law in the kinetic energy Φ_p(E_p) = (E_p/E_norm)^{− α}·⁽1 + E_p/E_br^{)α − α_h},

• (d)

  a double power-law in the kinetic energy with an exponential cut-off Φ_p(E_p) = (E_p/E_norm)^{− α}·⁽1 + E_p/E_br^{)α − α_h}·exp⁽ − E_p/E_c⁾,

where ${\mathit{E}}_{\mathrm{p}}\mathrm{=}({\mathit{p}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}\mathrm{+}{\mathit{m}}_{\mathrm{p}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{4}}{)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathrm{p}}{\mathit{c}}^{\mathrm{2}}$, and E_norm, E_c, and E_br are normalization, cut-off, and break energies, respectively, α is the effective power-law index for energies below the break energy, and α_h is the effective power-law index above the break energy. Matching the different expressions for the CR proton fluxes to the observed gamma-ray emission of the four specific objects considered here via the Kelner model reveals that for W49B and W44, the description (c) and for 3C 391 and CTB 37A the description (b) reproduces the observational data best. The corresponding best-fit parameters are shown in Table 1. Note that in order to obtain values for the CR proton flux normalization constant a_p, one has to make assumptions regarding the volumes V = 4πR³f_V/ 3 (to calculate the edge length l_c of the cubic emission regions), where f_V is a volume-filling factor (i.e., the homogeneous fraction (with hydrogen density n_H) of the volume of the SNR-MC complex where the pions are formed, cf. Table 1), and R is the radius of the SNR shock front (cf. Table 2), as well as the hydrogen densities of these emission regions. Here, emission volumes of V = 21.4 pc³ for W49B, V = 1.88 × 10⁴pc³ for W44, V = 8.38 × 10³pc³ for 3C 391, and V = 3.99 × 10⁴pc³ for CTB 37A, and homogeneous hydrogen densities of n_H = 100 cm^-3 are used. Since the break energies and the cut-off energies of all these objects are large compared to the maximum kinetic energy considered, only the single power-law (a), which is the E_p/E_br → 0 limit and the E_p/E_c → 0 limit of the parametric descriptions (b)−(d), is applied in calculating the ionization rates. After having determined the source term (13) for each object, the solution of the transport equation, the CR proton number density (14), is used in order to account for the propagation of the CR protons into the MCs and, thus, to calculate the ionization rates as functions of the penetration depth into the cloud as follows: the CR proton fluxes inside the MCs are related to the corresponding differential CR proton number densities by derivating the CR proton number density with respect to the CR proton kinetic energy and by multiplication with the effective particle speed v$\frac{{\mathrm{d}}^{\mathrm{3}}{\mathit{N}}_{\mathrm{p}}({r}\mathit{,p}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)}\mathit{,t})}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}\hspace{0.17em}\mathrm{d}\mathit{A}\hspace{0.17em}\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{p}}({r}\mathit{,p}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)}\mathit{,t})}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}}\mathrm{·}\mathit{v}\mathit{.}$(17)The effective particle speed is a combination of the diffusion speed v_diff = Δs/T_diff = 2D(p) / Δs, where Δs denotes the distance over which a particle was diffusing and T_diff = (Δs)²/ (2 D(p)) is the corresponding diffusion timescale, and the relativistic particle speed, ${\mathit{v}}_{\mathrm{rel}}\mathrm{=}\mathit{p}\hspace{0.17em}\mathit{c}\mathit{/}({\mathit{p}}^{\mathrm{2}}\mathrm{+}{\mathit{m}}_{\mathrm{p}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}{)}^{\mathrm{1}\mathit{/}\mathrm{2}}$, which is determined by initial values at the collision regions of the shock fronts and the MCs. If one of these speeds is dominant, it can be chosen as an approximation of the effective particle speed. It is shown in Sect. 6 that v_diff dominates over v_rel at all relevant penetration depths for the specific assumed diffusion coefficient. Here, in order to circumvent the problem that the time that has passed since the injection of the CR particles into the clouds until the measurement, while restricted to 0 <t<t_age, is generally not known, one couples the propagation time t to the position and the momentum of the particles via t = T_diff( | r | ,p) = | r | ²/ (2 D(p)), which is the mean time for a particle with momentum p to propagate the distance | r | in the case of standard three-dimensional Gaussian diffusion, or via t = T_rel( | r | ,p) = | r | /v_rel(p). Moreover, the ionization rate is evaluated only along the z-axis (cf. Sect. 3). The primary CR proton flux, as a function of this specific penetration depth into the MC at time t = T(z,p) after the injection reads $\frac{{\mathrm{d}}^{\mathrm{3}}{\mathit{N}}_{\mathrm{p}}(\mathit{z,p}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)})}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}\hspace{0.17em}\mathrm{d}\mathit{A}\hspace{0.17em}\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{p}}(\mathit{x}\mathrm{=}\mathit{y}\mathrm{=}\mathrm{0}\mathit{,z,p}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)}\mathit{,t}\mathrm{=}\mathit{T}\mathrm{\left(}\mathit{z,p}\mathrm{\right(}{\mathit{E}}_{\mathrm{p}}\mathrm{\left)}\mathrm{\right)})}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}}\mathrm{·}\mathit{v}$(18)with T = T_diff or T = T_rel and v = v_diff or v = v_rel, respectively. Note that the ionization rate at a certain penetration depth is not generated instantaneously, but builds up when the CR particles of the highest momentum considered have had the time to reach this penetration depth with their effective speed, and saturates once the CR particles of the lowest momentum arrive. This implies that the ionization rate should only be calculated up to the penetration depth ${\mathit{z}}_{\max \mathit{,}\mathrm{diff}}\mathrm{=}\sqrt{\mathrm{2}\hspace{0.17em}\mathit{D}\mathrm{\left(}{\mathit{p}}_{\min }\mathrm{\right)}\hspace{0.17em}{\mathit{t}}_{\mathrm{age}}}$  or z_max,rel = v_rel(p_min) /t_age, respectively, which correspond to the distances the particles with minimum momentum move within the ages of the SNRs.

After having computed the constrained CR proton fluxes (18) for the MCs, the ionization rate of molecular hydrogen as a function of the penetration depth into the cloud z can, therefore, be represented as the following integral $\begin{array}{ccc}{\mathit{\zeta }}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \left(\mathrm{1}\mathrm{+}\mathit{\phi }\right){\mathrm{\int }}_{{\mathit{p}}_{\min }}^{{\mathit{p}}_{\max }}\frac{\mathrm{d}{\mathit{n}}_{\mathrm{p}}(\mathit{z,p}\mathrm{\left(}{\mathit{E}}_{\mathrm{p}}\mathrm{\right)}\mathit{,t}\mathrm{=}\mathit{T}\mathrm{\left(}\mathit{z,p}\mathrm{\right(}{\mathit{E}}_{\mathrm{p}}\mathrm{\left)}\mathrm{\right)})}{\mathrm{d}{\mathit{E}}_{\mathrm{p}}}\mathrm{·}\mathit{v}(\mathit{p}\mathrm{)}\\ & & \end{array}$(19)

The cross-section for direct ionization of molecular hydrogen by protons, ${\mathit{\sigma }}_{\mathrm{ion}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}\mathit{p}\mathrm{\right)}$, can be given via the parametric expression (Rudd et al. 1985) ${\mathit{\sigma }}_{\mathrm{ion}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{l}}\mathrm{\left(}\mathit{p}\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{h}}\mathrm{\left(}\mathit{p}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathrm{l}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{h}}\mathrm{\left(}\mathit{p}\mathrm{\right)}}$(20)with the low-energy component ${\mathit{\sigma }}_{\mathrm{l}}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{a}}_{\mathrm{0}}^{\mathrm{2}}\hspace{0.17em}\mathit{C}\hspace{0.17em}\mathit{x}\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathit{D}}$ and the high-energy component ${\mathit{\sigma }}_{\mathrm{h}}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{a}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{[}\mathit{A}\mathrm{·}\ln \mathrm{(}\mathit{x}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}\mathit{B}\mathrm{]}\mathit{x}\mathrm{(}\mathit{p}{\mathrm{)}}^{-1}$, where $\mathit{x}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{}^{\mathrm{(}}\mathit{m}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathit{/}{\mathit{I}}_{\mathrm{H}}{}^{\mathrm{)}}{}^{\mathrm{(}}[\mathrm{1}\mathrm{+}(\mathit{p}\mathit{/}\mathrm{\left(}{\mathit{m}}_{\mathrm{p}}\mathit{c}\mathrm{\right)}{)}^{\mathrm{2}}{]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}$, a₀ = 5.29 × 10^-9 cm is the Bohr radius, the parameters A = 0.71, B = 1.63, C = 0.51, as well as D = 1.24, m_e and m_p are the electron and proton masses, respectively, and I_H = 13.598 eV is the ionization potential of atomic hydrogen. The integration boundaries p_min = 0.15 m_pc and p_max = 0.86 m_pc correspond to the minimum and maximum CR proton energies E_p,min = 10 MeV and E_p,max = 280 MeV. In general, the dimensionless quantity φ, which accounts for secondary ionization, depends on the energies of the primary CR protons, but for the energy range of the protons covered in this work, φ = 0.6 is an adequate constant approximation (Cravens & Dalgarno 1978). The secondary electrons have, on average, rather low kinetic energies that they lose rapidly in interactions with matter, electromagnetic radiation and the magnetic fields inside the MCs. This allows to consider the effect of secondary ionization as an instantaneous phenomenon, i.e., to neglect the propagation of the secondary electrons. The large number of primary CR particles moreover justifies the treatment of the corresponding secondary ionization via a constant. The integral (19) is solved numerically for each object, using the adaptive Monte Carlo algorithm VEGAS introduced by Lepage (1978). While in standard Monte Carlo routines the sample points are chosen according to a preassigned distribution function, VEGAS iteratively concentrates, in an adaptive manner, the distribution of the sample points on those intervals where the integrand contributes the most, starting from a uniform distribution of the sample points. This algorithm is slower than, e.g., a Gauss-Kronrod quadrature, but is more reliable and stable, particularly for rapidly changing integrands in one-dimensional integrations, which are performed here. The VEGAS routine is implemented as a C++ algorithm via the GNU Scientific Library (Galassi 2009). The results of these integrations for v = v_diff and t = T_diff as well as v = v_rel and t = T_rel are presented in Sect. 6.

5. Photoionization profiles

In order to determine whether there is a dominant process generating the ionization signatures of SNR-MC systems, one has to derive the ionization profiles for all eligible processes. To this end, the profiles emerging from photoionization are studied in addition to those induced by CR protons.

Photoionization is primarily induced by Extreme UV (EUV) radiation (13.6 eV ≤ E ≤ 200 eV) and X-rays (E> 200 eV). EUV radiation is quite effective in ionizing molecular hydrogen, but its penetration depth into MCs is short because it is rapidly absorbed. X-rays, however, are capable of traversing larger columns of matter, with their soft (low-energy) component losing energy faster than their hard (high-energy) component. In fact, hard X-rays might penetrate even farther into MCs than CRs (Tielens 2010). In this section, the X-ray fluxes in the SNR-MC systems W49B, W44, 3C 391, and CTB 37A are derived at the MC surfaces using observational data, assuming that the fluxes are composed of synchrotron radiation of CR electrons emitted at the SNR shock fronts, and for the transport into the MCs via the Beer-Lambert law. From these fluxes, the photoinduced ionization rates are calculated as functions of the penetration depth, and compared with those induced by CRs.

5.1. Photon fluxes at the MC surfaces

The X-ray fluxes in the vicinity of MCs are derived assuming isotropic radiation emission from the SNR shock fronts, which are known to emit synchrotron radiation of CR electrons in the form of X-rays. It is furthermore assumed that the SNR shock fronts coincide with the surfaces of the clouds, in accordance with the geometry used for the calculation of the ionization induced by CR protons. Note that X-ray observations of SNRs associated with MCs provide near-Earth X-ray flux data, which is inadequate for the derivation of the X-ray fluxes at the MCs because, on the one hand, there are energy-dependent losses from interactions of the photons with the interstellar medium in between the emission region and the detector, and on the other hand, there is a geometrical 1/d²-attenuation of the photon fluxes from the source along the path toward the instrument. Thus, the X-ray fluxes at the cloud surfaces, F_{X - ray, cs}, can be derived from the X-ray fluxes at the instrument, F_{X - ray, in}, by ${\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}\mathit{,}\hspace{0.17em}\mathrm{cs}}\mathrm{=}{\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}\mathit{,}\hspace{0.17em}\mathrm{in}}\mathrm{·}{\left(\frac{\mathit{d}}{\mathit{R}}\right)}^{\mathrm{2}}\mathit{,}$(21)where d is the distance between the SNR-MC system and Earth, and R is the radial extension of the SNR shock front. The values of the parameters d and R for the objects considered are given in Table 2. The attenuation of the X-ray fluxes from interactions with traversed matter between the emission region and the detector is already taken into account in the data samples of the X-ray observations F_{X - ray, in}. Since the photoinduced ionization rate is given, as in the case of CR-induced ionization (cf. Eq. (15)), in terms of an integral expression that includes the photon flux as a function of the photon energy, a continuous description of the observational data is required. This is achieved by fitting a continuous parametric model function for the X-ray fluxes to the observational data samples. The parametric model function for the photon fluxes of two of the four SNR-MC systems studied here, W44 and 3C 391, is chosen to be a triple power-law ${\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}\mathit{,}\hspace{0.17em}\mathrm{in}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{0}}\mathrm{·}{\left(\frac{\mathit{E}}{\mathrm{1}\hspace{0.17em}\mathrm{keV}}\right)}^{\mathrm{-}\mathit{a}}\mathrm{·}{\left(\mathrm{1}\mathrm{+}\frac{\mathit{E}}{{\mathit{b}}_{\mathrm{1}}}\right)}^{\mathrm{-}{\mathit{a}}_{\mathrm{1}}}\mathrm{·}{\left(\mathrm{1}\mathrm{+}\frac{\mathit{E}}{{\mathit{b}}_{\mathrm{2}}}\right)}^{\mathrm{-}{\mathit{a}}_{\mathrm{2}}}\mathit{,}$(22)where a, a₁, and a₂ are dimensionless spectral indices, b₁ and b₂ denote break energies, respectively, and F₀ is a constant. Owing to their significantly different X-ray data samples, different parametric models for CTB 37A and W49B are used. These are, for CTB 37A, ${\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}\mathit{,}\hspace{0.17em}\mathrm{in}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{0}}\mathrm{·}{\left(\frac{\mathit{E}}{\mathrm{1}\hspace{0.17em}\mathrm{keV}}\right)}^{\mathrm{-}\mathit{a}}\mathrm{·}{\left(\mathrm{1}\mathrm{+}\frac{\mathit{E}}{{\mathit{b}}_{\mathrm{1}}}\right)}^{\mathrm{-}{\mathit{a}}_{\mathrm{1}}}\mathit{,}$(23)and for W49B, ${\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}\mathit{,}\hspace{0.17em}\mathrm{in}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{E}\mathit{<}\mathrm{3.5}\hspace{0.17em}\mathrm{keV}\\ & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{E}\mathrm{\ge }\mathrm{3.5}\hspace{0.17em}\mathrm{keV}\mathit{.}\end{array}$(24)Applying these simpler parametric descriptions leads to an improved χ²-value per degree of freedom since a lesser number of fitting parameters is used. The best-fit parameters are given in Table 3. The extremely low value of χ² per degree of freedom for the object W44 results from the large error bars of the X-ray data, cf. Fig. 2b. For the objects W44, 3C 391, and CTB 37A, certain X-ray data points were excluded because they are caused by line emission originating along the path between the near-Earth measuring instrument and the SNR-MC systems (Ozawa et al. 2009; Harrus et al. 2006; Chen & Slane 2001; Sezer et al. 2011).

The fitting of the parametric models to the observational X-ray data samples was performed by minimizing χ² per degree of freedom based on the Nelder-Mead method (Nelder & Mead 1965). This method is a non-linear optimization technique that does not require the derivatives of χ² in order to minimize χ² with respect to the optimization parameters. It is particularly useful in multi-dimensional parameter space optimization. This routine is implemented as a C++ code via the GNU Scientific Library.

Since the total photoabsorption cross-section, which is relevant for the calculation of photoionization, is highest at low photon energies, the treatment of the EUV domain is crucial and consequently has to be included in this study. To evaluate the influence of low-energy photons on the photoionization rate, three descriptions of the low-energy photon fluxes are considered:

• (I)

  extrapolation of the parametric model function used to fit theobserved X-ray flux values toward lower photon energies,

• (II)

  a constant photon flux equal to the value of the parametric model function at the lowest photon energy covered by the observational data,

• (III)

  a hard cut-off below the lowest photon energy covered by the observational data, that is, no extrapolation,

where the lowest photon energies for which X-ray data are available are 0.5 keV for W49B as well as 3C 391, 0.56 keV for W44, and 0.3 keV for CTB 37A. The parametric models and the different descriptions of the low-energy photon fluxes are shown along with the observational X-ray data samples in Fig. 2.

Fig. 2 Parametric models (I) (red, solid lines), (II) (blue, long-dashed lines), and (III) fitted to X-ray flux data samples (black dots and error bars) for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d). The X-ray data samples used are taken from Ozawa et al. (2009) (Suzaku) for a), Harrus et al. (2006) (XMM-Newton) for b), Chen & Slane (2001) (ASCA) for c), and Sezer et al. (2011) (Suzaku) for d).

Fig. 3 CR-induced ionization profiles for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d), originating from diffusion-driven motion (black, solid lines) and relativistic, kinematic motion (cyan, long-dashed lines), respectively.

Fig. 4 Ionization profiles for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d), induced by CRs (black, solid lines), as well as X-ray spectra extrapolated toward lower energies by a power-law (red, short-dashed lines), by a constant (gray, long-short-dashed lines), and without extrapolation (purple, long-short-short-dashed lines), respectively.

Fig. 5 Total ionization profiles (blue, solid lines) and ionization profiles induced by X-ray spectra extrapolated toward lower energies by a power-law alone (red, dashed lines), respectively, for W49B a), W44 b), 3C 391c), and CTB 37A d).

5.2. Photon fluxes inside MCs and photoinduced ionization profiles

Using the photon fluxes obtained at the surfaces of the MCs, the corresponding fluxes inside the clouds, which undergo an attenuation due to interactions with cloud matter, as functions of the photon energy and the penetration depth, are calculated. This is done by means of the Beer-Lambert formula, ${\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}}\mathrm{\left(}\mathit{z,E}\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}\mathit{,}\hspace{0.17em}\mathrm{cs}}\mathrm{·}\exp (\mathrm{-}\mathit{\tau }\mathrm{\left(}\mathit{z,E}\mathrm{\right)})\hspace{0.17em}\mathit{,}$(25)which relates the photon flux F_{X - ray} at a certain penetration depth and energy to the initial photon flux (21) at the cloud surface. The quantity τ(z,E) = σ_pa(E) n_Hz is the optical depth for a homogeneous cloud, σ_pa(E) is the total photoabsorption cross-section, and n_H is the hydrogen density. It should be remarked that the actual photon paths are longer than the direct ones, for which deviations caused by scattering are neglected. Therefore, the photoinduced ionization signatures calculated here are upper limits, leading to conservative estimates with respect to the detectability of CR-induced ionization signatures, since under the influence of scattering, the photons lose more energy at lower penetration depths, which leaves less energy for ionization at large penetration depths. The photoinduced ionization rate, as a function of the penetration depth, following Maloney et al. (1996), yields ${\mathit{\zeta }}_{\mathrm{X}\mathrm{-}\mathrm{ray}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{i}}}{{\mathit{I}}_{{\mathrm{H}}_{\mathrm{2}}}}{\mathrm{\int }}_{{\mathit{E}}_{\min }}^{{\mathit{E}}_{\max }}{\mathit{F}}_{\mathrm{X}\mathrm{-}\mathrm{ray}}\mathrm{\left(}\mathit{z,E}\mathrm{\right)}\mathrm{·}\mathit{E}\mathrm{·}{\mathit{\sigma }}_{\mathrm{pa}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{E,}$(26)where f_i is the fraction of the absorbed photon energy leading to ionization, I_H2 = 15.4 eV is the ionization potential of molecular hydrogen, and the integral is the total energy deposition by the absorbed photons per hydrogen nucleus. The parametrization used for a description of the total photoabsorption cross-section per hydrogen nucleus is also provided in Maloney et al. (1996) as an empirical, broken power-law fit to experimental data $\begin{array}{ccc}& & {\mathit{\sigma }}_{\mathrm{pa}}\mathrm{\left(}\mathit{E}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{·}{\left(\frac{\mathit{E}}{\mathrm{1}\hspace{0.17em}\mathrm{keV}}\right)}^{\mathrm{-}\mathit{\gamma }}\\ & & \\ & & \mathrm{with}\\ & & {\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\mathrm{6.0}\mathrm{\times }{\mathrm{10}}^{-23}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}\mathit{\gamma }\mathrm{=}\mathrm{3}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\mathrm{0.014}\hspace{0.17em}\mathrm{keV}\mathit{<}\mathit{E}\mathrm{\le }\mathrm{0.5}\hspace{0.17em}\mathrm{keV}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\mathrm{2.6}\mathrm{\times }{\mathrm{10}}^{-22}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}\mathit{\gamma }\mathrm{=}\mathrm{8}\mathit{/}\mathrm{3}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{0.5}\hspace{0.17em}\mathrm{keV}\mathit{<}\mathit{E}\mathrm{\le }\mathrm{7.0}\hspace{0.17em}\mathrm{keV}\mathit{,}\\ & & {\mathit{\sigma }}_{\mathrm{0}}\mathrm{=}\mathrm{4.4}\mathrm{\times }{\mathrm{10}}^{-22}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{\gamma }\mathrm{=}\mathrm{8}\mathit{/}\mathrm{3}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{7.0}\hspace{0.17em}\mathrm{keV}\mathit{<}\mathit{E}\hspace{0.17em}\mathit{.}\end{array}$(27)The breaks in the above expression come from contributions of heavy elements in the MCs to the total cross-section for which there is a certain threshold energy for shell absorption, namely oxygen K-shell absorption at 0.5 keV and iron K-shell absorption at 7 keV. In contrast to the ionization rate induced by CRs, the photoinduced ionization rate is not dominated by direct ionization events. The fraction of the absorbed photon energy yielding ionization is approximately f_i ≈ 40% of the deposed photon energy (Voit 1991; Dalgarno et al. 1999), so a photon with E = 1 keV absorbed by the cloud medium induces a total of 26 ionization events. This allows a description of the total photoionization rate in terms of the photon energy deposed in the cloud matter. The lower integration boundary is set at the photon energy required to produce an ion pair, E_min = 37 eV, in order to obtain an upper limit of photoinduced ionization, but it should be noted that no X-ray measurements for E ≲ 0.1 keV exist. The upper integration boundary is fixed at 10 keV, but a higher value would not change the results significantly because of the rapid decrease of both the total photoelectric energy deposition cross-section and the photon flux with increasing photon energy. The ionization profiles induced by photons (26) are calculated numerically for each object as well as each parametric description of the low-energy photon flux (I)−(III) using the VEGAS integration technique and shown in the following section.

6. Results

The CR- and photoinduced ionization profiles for W49B, W44, 3C 391, and CTB 37A computed in Sects. 4 and 5.2, respectively, are shown and compared. The cases of purely diffusion-driven motion as well as purely relativistic, kinematic motion of the CR protons and their implications for the CR-induced ionization rate are shown in Fig. 3.

For the specific diffusion coefficient assumed here, the diffusion-driven motion dominates the relativistic, kinematic motion for all considered particle momenta. This can be seen from the lower overall magnitude of the ionization rate for kinematic motion, indicating that the kinematically driven particles take longer to reach a certain penetration depth and consequently lose more energy at low penetration depths. This in turn means that fewer particles are capable of ionizing the molecular hydrogen at large penetration depths. This relation between the diffusion and the relativistic speeds, v_diff>v_rel, arises because the magnetic field inhomogeneities that determine the diffusion coefficient effectively boost the CR particles to diffusion speeds exceeding their relativistic speeds. Since the diffusion-driven and the relativistic, kinematic motions occur simultaneously, the effective particle speed is a combination of both. However, because the former is dominant, it is used as an approximation for the effective CR particle speed.

In Fig. 4, the photoionization profiles for all three low-energy extrapolations of the X-ray data are presented, and compared with the ionization profiles induced by CRs. The power-law extrapolation (I) toward the UV photon energy domain leads to overestimates of the photon fluxes and, therefore, to reasonable upper limits of the photoinduced ionization rates, while the cut-off X-ray fluxes (III) yield lower limits of the photoionization rates. Particularly at very short penetration depths, the ionization rates differ distinctly, because on the one hand, at low photon energies the total photoabsorption cross-section is large, which means that low-energy photons are absorbed very efficiently and, on the other hand, the parametric descriptions of the photon fluxes at low energies are different. With increasing penetration depth, the different photoionization profiles become asymptotically equivalent, which has two reasons. On the one hand, at high photon energies the total photoabsorption cross-section is small, and on the other hand, the parametric descriptions of the photon fluxes at high energies are identical.

Among the four objects studied in this work, W49B shows the highest CR-induced ionization rate, more than two orders of magnitude higher than the corresponding ionization rate found for 3C 391, which has the second-highest CR-induced ionization rate of the objects studied here. The rather flat photoinduced ionization profile for W49B results from the flat spectral shape of the photon flux (cf. Fig. 2). The ionization in this cloud is clearly dominated by CRs. In CTB 37A, the influence of photons on the ionization significantly exceeds the influence of CR protons. Figure 4d indicates that this behavior does not depend on the extrapolation of the X-ray flux toward lower photon energies, because even without any extrapolation (III), the photoinduced ionization rate dominates, except for the very surface region of the MC. For both W44 and 3C 391, the CR-induced ionization dominates up to a certain penetration depth, 6.5 pc for W44 and 3.2 pc for 3C 391, while at larger penetration depths photons dominate the ionization. Since the parametric descriptions (I) of the low-energy photon fluxes used for these calculations are overestimates of the real photon fluxes, the actual photoionization rates are lower. Moreover, it should be kept in mind that the effects of magnetic fields, that is, magnetic focusing and magnetic mirroring, on the CR-induced ionization rates have not been taken into account in this approach. Padovani & Galli (2011) considered these effects in their studies, and came to the conclusion that the influence of magnetic fields, on average, decreases the ionization rate induced by CR protons by a factor of 3 to 4. For the rather diffuse clouds examined here, the decrease of the CR-induced ionization rate would drop to a factor of 3 or less. This would not change the overall trend of which ionization process is dominant in the studied SNR-MC systems, because the ionization profiles induced by CR protons and photons differ by at least a factor of 10 at all relevant penetration depths for both W49B and CTB 37A. For W44 and 3C 391, they differ by approximately a factor of 10 for a wide range of penetration depths even if extrapolation (I) is used. For W44, 3C 391, and CTB 37A, volume-filling factors of 1 are applied. While this assumption is common, a volume-filling factor ≪1 is more realistic. Taking this into consideration, the CR proton flux normalizations would increase, and so would the CR-induced ionization rates. Hence, for W44 and 3C 391, an observational distinction between CRs and photons as the dominant source of ionization seems feasible. For the calculation of the CR proton normalization for W49B, a volume-filling factor of 0.06 was applied (Abdo et al. 2010b). Changing this factor to 1 would yield a CR-induced ionization rate that would still be the highest among all the objects considered and would also be orders of magnitude higher than the corresponding ionization rate induced by photons. It should be pointed out that the total ionization rates very close to the surfaces of the MCs, i.e., in the outer ~0.1 pc, strongly depend on the estimates used to model the EUV regime. In these regions, the contribution of EUV radiation to the total ionization rate is highest (cf. the rapid decline of the photoinduced ionization rate at very low penetration depths in Fig. 4, particularly the difference between the different extrapolations of the X-ray fluxes seen in Figs. 4b−4d). Since the ionization profiles induced by different processes cannot be detected separately, it is reasonable to compare the total ionization rates with the ionization rates expected from photoinduced ionization alone. Detecting an ionization rate that clearly exceeds the theoretical photoinduced ionization rate (enhanced ionization rate) would imply a significant contribution from low-energy CR protons, indicating that high-energy CR protons are the driving force for gamma-ray emission. In the following, for a conservative comparison, only the extrapolation of the X-ray fluxes toward lower photon energies by the power-law (I) is discussed. The total ionization profiles analyzed here are simply the sum of the photoinduced ionization profiles, with the power-law extrapolation (I), and the CR-induced ionization profiles, with v = v_diff as the effective CR particle speed and the corresponding diffusion time t = T_diff. A comparison of the total ionization profiles and the photoinduced ionization profiles for each object is shown in Fig. 5. A huge difference between the CR- and photoinduced ionization profiles implies that the total ionization rate almost coincides with the ionization profile from the dominant process in the MC. For W49B, the total ionization rate practically coincides with the CR-induced ionization rate (cf. Fig. 4a), while for CTB 37A, it coincides with the ionization rate induced by photons (cf. Fig. 4d). The impact of the calculated ionization rates on the chemistry inside the MCs, with a focus on the detectability of ionization rates, is discussed in the following section.

7. Conclusions and outlook

We calculated the position-dependent H₂ ionization rates of molecular clouds near supernova remnants, induced by photons and cosmic-ray protons, respectively, and compared them to examine a potential hadronic formation scenario of the observed gamma radiation from these systems. The cosmic-ray proton fluxes were constructed by means of the analytic solution of the transport equation for cosmic-ray protons inside the clouds, while the photon fluxes were obtained using the Beer-Lambert law. In addition to the systems consisting of a supernova remnant and a molecular cloud that we studied here, the solution of the transport equation can also be used for many other astrophysical objects with scalar, momentum-dependent diffusion and any type of momentum losses, such as stellar winds, cosmic-ray diffusion in the interstellar medium, or gamma-ray bursts.

Enhanced ionization rates compared with the ionization rates induced by photons alone are expected if the gamma-ray emission seen from systems containing supernova remnants and molecular clouds originates from the decay of neutral pions that are formed in proton-proton interactions of the cosmic-ray protons with ambient hydrogen. The method presented here is particularly promising for objects of low X-ray fluxes and steep, incident cosmic-ray proton spectra, when the spectra can be extrapolated toward lower energies by single power-laws. For the SNR-MC system W49B, the ionization rate induced by cosmic-ray protons significantly exceeds the ionization rate induced by photons. In contrast, for CTB 37A, the total ionization rate is dominated by photoinduced ionization, while for the objects W44 and 3C 391, the results presented here favor cosmic-ray protons as the dominant source of ionization. The results shown in Figs. 4 and 5 are a definite motivation for further studies of enhanced ionization rates attributable to cosmic-ray sources. But since the ionization rate is not directly detectable, one has to derive it for instance from observational spectroscopic data. One method to estimate the ionization rate is to examine the absolute intensity of molecular rotation-vibrational lines, as stated in Becker et al. (2011). It is also possible to analyze the ratio of the abundance of certain molecules or molecular ions (see, e.g., O’Donnell & Watson 1974; Black & Dalgarno 1977; Indriolo et al. 2009). Computing rotation-vibrational lines requires solid estimates on several physical parameters of the molecular clouds, such as composition, temperature, ionization degree, electron density, and a few others. Many of these parameters are quite difficult to estimate, which means that predictions for rotation-vibrational line intensities are inaccurate. Hence, the molecular ion H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ is especially well-suited for observational examinations of ionization rates because only very few cloud parameters are involved in calculating its formation rate (Goto et al. 2013). Particularly in W49B and 3C 391, very high ionization rates are found (10^-12s^-1 at depths of 1 pc into clouds with a hydrogen density of n_H = 100 cm^-3). The resulting ionization degrees are on the order of 10^-4 at depths of 1 pc and 10^-3 or higher at the cloud surfaces. However, a penetration depth of 1 pc at a density of 100 cm^-3 corresponds to a column density of 3 × 10²⁰ cm^-2, so these regions will also be significantly affected by stellar Far UV radiation (E< 13.6 eV), which has an impact on the chemistry and temperature in the clouds. Accordingly, it is also possible that these outer layers of the clouds are photodissociation regions. This needs to be taken into account in the astrochemistry models for the prediction of molecular rotation-vibrational lines to test the theoretical ionization rates by means of observations.

Some of the molecular species linked to high ionization rates are most effectively and hence solely detectable in absorption. This means that the molecular clouds in which these molecules are to be found have to be somewhat close to Earth, and the sightlines toward these clouds need to feature a bright source in the background for illumination such that the light is absorbed by the molecules in order to trace them. This might restrain the target-selection process of systems composed of a supernova remnant associated with molecular clouds to nearby complexes, that is, a few kpc from Earth, such as W44. These are particularly interesting because the resolution available with current and next-generation telescopes will allow studying their substructure. The currently most promising experiments for the suggested correlation study are the Cherenkov Telescope Array (CTA) and Fermi/LAT for the detection of the gamma rays, and the Atacama Large Millimeter/submillimeter Array (ALMA) for the detection of the molecular lines. CTA will reach a spatial resolution of ~0.03° (Wagner & CTA Consortium 2010; CTA Consortium & Hermann 2011). While this is still lower than for submm, IR, or far-IR (FIR) telescopes, CTA measurements can be used to determine regions of enhanced gamma-ray emission that then can be observed with IR telescopes to search for a correlation predicted for hadronic emission scenarios. The spatial resolution for a selection of relevant telescopes for the suggested correlation study is given in Table 4. Since different observations indicate high ionization rates through the measurement of H$\begin{array}{c}\mathrm{+}\\ \mathrm{3}\end{array}$ (see, e.g., Indriolo et al. 2010), it is possible to find the ionization-induced signatures that are spatially correlated with the GeV gamma-ray emission, which will help pinpointing supernova remnants as sources of cosmic-ray protons.

Apart from the specific supernova remnants discussed above, there are currently six additional ones that are associated with molecular clouds. However, for most of these no suitable X-ray data are currently available, because the X-ray fluxes in the corresponding energy ranges are either below the detection limit or are not aligned with the gamma-ray emission regions. For these systems, precise measurements of X-ray spectra by experiments such as Swift, XMM-Newton, Suzaku, or Chandra, particularly for energies E ~ 0.1 keV, can help to estimate the photoionization rate.

To obtain even more accurate theoretical predictions for ionization profiles induced by cosmic rays, additional components of both the cosmic-ray spectrum and the cloud matter can be taken into account. Considering also heavier CR particles in addition to protons for calculating the ionization rate of molecular hydrogen, the model can be applied by adequately adapting the Coulomb-loss coefficient and the ionization cross-section, for example by using the Bethe-Bloch approximation (Bethe 1933; Padovani et al. 2009). The relative abundances of cosmic-ray helium or heavier nuclei at the acceleration regions are poorly known, however. In addition to the ionization of molecular hydrogen in the clouds, the model can also be used to calculate the ionization rates of other species such as helium by simply replacing the ionization cross-section. This further improves estimates of the ionization degree in the molecular clouds required for the chemistry models.

To conclude, observable ionization signatures are expected and, combined with gamma-ray observations, correlations that can only be caused by a common primary component (see, e.g., Becker 2008) can be revealed.

Online material

Appendix A: General solution of the transport equation

In order to find a general solution of the transport Eq. (1) using Green’s method, first the fundamental solution G(r,p,t | r₀,p₀,t₀), that is, the solution for the Dirac source distribution $\begin{array}{ccc}\mathit{Q}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}\mathrm{=}{\mathit{Q}}_{\mathrm{0}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{p}\mathrm{-}{\mathit{p}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathit{,}& & \end{array}$is determined.The fundamental transport equation is then given by $\frac{\mathit{\partial G}}{\mathit{\partial t}}\mathrm{-}\mathit{D}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{\Delta }\mathit{G}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial p}}(\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{·}\mathit{G})\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{p}\mathrm{-}{\mathit{p}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(28)In terms of the function R(r,p,t) = b(p)·G(r,p,t), Eq. (28) becomes $\frac{\mathit{\partial R}}{\mathit{\partial t}}\mathrm{-}\mathit{D}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{\Delta }\mathit{R}\mathrm{-}\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\frac{\mathit{\partial R}}{\mathit{\partial p}}\mathrm{=}\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{p}\mathrm{-}{\mathit{p}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(29)With the new coordinate $\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{=}\mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}\mathit{c}}{\mathrm{3}\hspace{0.17em}{\mathit{a}}_{\mathrm{ad}}}\ln \left({\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)}^{\mathrm{3}}\right)|\begin{array}{c}\mathit{p}\\ {\mathit{p}}_{\mathrm{0}}\end{array}\mathit{,}$(30)induced by ∂_β = − b(p)∂_p, and δ(β − β₀) = b(p₀)·δ(p − p₀), Eq. (29) yields $\frac{\mathit{\partial R}}{\mathit{\partial t}}\mathrm{-}\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{\Delta }\mathit{R}\mathrm{+}\frac{\mathit{\partial R}}{\mathit{\partial \beta }}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(31)This partial differential equation can be solved by first applying a Laplace transformation L[·] with respect to the time variable t$L\mathrm{\left[}\mathit{R}\mathrm{\right]}\mathrm{\equiv }G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\exp \mathrm{(}\mathrm{-}\mathit{st}\mathrm{)}\mathit{R}\mathrm{\left(}{r}\mathit{,\beta ,t}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{t}\mathit{.}$(32)Note that the lower integration limit in the definition of the Laplace transformation (32) fixes t₀ = 0 without loss of generality. The Laplace transform of Eq. (31) reads $\begin{array}{ccc}\mathrm{-}\mathit{R}\mathrm{(}{r}\mathit{,\beta ,t}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{+}\mathit{s}\hspace{0.17em}G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{-}\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{\Delta }G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{+}{\mathit{\partial }}_{\mathit{\beta }}G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathit{.}& & \end{array}$(33)Since R is related linearly to the differential CR proton number density, R ∝ dN/ dV, where N is the (finite) number of protons, and V is the volume within which the protons are distributed, this function vanishes at t = 0, because the distribution of CR protons at this time is restricted to the boundary r = r₀ of the MC. Then, Eq. (33) reduces to the inhomogeneous, linear partial differential equation of first order in β and second order in r, $\mathrm{\Delta }G\mathrm{-}\frac{\mathrm{1}}{\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}\mathrm{[}\mathit{s}\mathrm{+}{\mathit{\partial }}_{\mathit{\beta }}\mathrm{]}G\mathrm{=}\mathrm{-}\frac{\mathit{\delta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}}{\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(34)This equation can be solved explicitly by using Duhamel’s principle (Duhamel 1838; Courant & Hilbert 2008), which is a general method to find solutions of inhomogeneous, linear partial differential equations in terms of the solutions of the Cauchy problems for the corresponding homogeneous partial differential equations, that is, by interpreting the inhomogeneity as a boundary value condition in a higher-dimensional space labeled with an additional auxiliary variable. This principle is applied here on a linear partial differential equation with a product-separable inhomogeneity of single-variable factors with respect to the variables r and β, $(\mathit{Ô}{r}\mathrm{+}\mathit{Ô}\mathit{\beta ,s})G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathit{F}}_{\mathrm{1}}\mathrm{\left(}{r}\mathrm{\right)}{\mathit{F}}_{\mathrm{2}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathit{,}$(35)where $\mathit{Ô}{r}\mathrm{=}\mathrm{\Delta }\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{Ô}\mathit{\beta ,s}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}\mathrm{[}\mathit{s}\mathrm{+}{\mathit{\partial }}_{\mathit{\beta }}\mathrm{]}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{F}}_{\mathrm{1}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}^{\mathrm{3}}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}\mathrm{)}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{F}}_{\mathrm{2}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{=}\frac{\mathit{\delta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}}{\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}\mathrm{·}$(36)From the structure of the left-hand side of Eq. (35), it directly follows that the solution of the corresponding homogeneous equation is product-separable G_hom = S(r)P(β,s). An embedding of the original (r,β,s)-space into the extended, higher-dimensional (r,β,s,u)-space, where $\mathit{u}\mathrm{\in }{R}_{\mathrm{0}}^{\mathrm{+}}$, implies that there is a family of functions S_u(r): = S(r,u) and P_u(β,s): = P(β,s,u), fulfilling the relations Ô_rS(r,u) = ∂_uS(r,u) and Ô_β,sP(β,s,u) = ∂_uP(β,s,u) with the boundary conditions $\mathit{S}\mathrm{(}{r}\mathit{,u}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}{\mathit{F}}_{\mathrm{1}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{,}\hspace{0.17em}\hspace{0.17em}\mathit{P}\mathrm{(}\mathit{\beta ,s,u}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}{\mathit{F}}_{\mathrm{2}}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathit{,}\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}\mathit{S}\mathrm{(}{r}\mathit{,u}\mathrm{=}\mathrm{\infty }\mathrm{)}\mathrm{=}\mathit{P}\mathrm{(}\mathit{\beta ,s,u}\mathrm{=}\mathrm{\infty }\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}$(37)such that the full, non-trivial solution of Eq. (35) is given by the integral $G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathit{S}\mathrm{\left(}{r}\mathit{,u}\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{\beta ,s,u}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{u}\mathit{.}$(38)Technically, this integral represents the “summation” over the entire family of homogeneous solutions in the extended coordinate space with boundary conditions compatible with the inhomogeneity. Within this setting, Eq. (34) decouples into the simpler set of partial differential equations $\frac{\mathit{\partial S}\mathrm{\left(}{r}\mathit{,u}\mathrm{\right)}}{\mathit{\partial u}}\mathrm{=}\mathrm{\Delta }\mathit{S}\mathrm{\left(}{r}\mathit{,u}\mathrm{\right)}$(39)and $\frac{\mathit{\partial P}\mathrm{\left(}\mathit{\beta ,s,u}\mathrm{\right)}}{\mathit{\partial u}}\mathrm{+}\frac{\mathrm{1}}{\mathit{D}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}\mathrm{[}\mathit{s}\mathrm{+}{\mathit{\partial }}_{\mathit{\beta }}\mathrm{]}\mathit{P}\mathrm{\left(}\mathit{\beta ,s,u}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(40)which has to be solved in order to determine the fundamental solution G via the integral in Eq. (38). The solution of Eq. (39) is the well-known heat kernel of three-dimensional Euclidean space $\mathit{S}\mathrm{\left(}{r}\mathit{,u}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{4}\mathit{\pi u}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\exp (\mathrm{-}\frac{\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}\mathit{u}})\mathrm{·}$(41)A solution of Eq. (40) can directly be found after performing a Laplace transformation with respect to the variable u$L\mathrm{\left[}\mathit{P}\mathrm{\right]}\mathrm{\equiv }P\mathrm{\left(}\mathit{\beta ,s,q}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\exp \mathrm{(}\mathrm{-}\mathit{qu}\mathrm{)}\mathit{P}\mathrm{\left(}\mathit{\beta ,s,u}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{u}\mathit{.}$(42)Then, Eq. (40) becomes $(\mathit{qD}\mathrm{\left(}\mathit{\beta }\mathrm{\right)}\mathrm{+}\mathit{s})P\mathrm{+}\frac{\mathit{\partial }P}{\mathit{\partial \beta }}\mathrm{=}\mathit{\delta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathit{.}$(43)Note that the boundary condition P(β,s,u = 0) = δ(β − β₀) /D(β) was already fixed in (37). Using an integrating factor of the form $\exp \left(\mathrm{-}\mathit{s\beta }\mathrm{-}\mathit{q}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}\right)$, Eq. (43) can be rewritten as $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial \beta }}\left[\exp (\mathit{s\beta }\mathrm{+}\mathit{q}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }})\mathrm{·}P\right]\mathrm{=}\mathit{\delta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\exp (\mathit{s\beta }\mathrm{+}\mathit{q}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }})& & \end{array}$(44)and solved by simple integration, leading to $P\mathrm{\left(}\mathit{\beta ,s,q}\mathrm{\right)}\mathrm{=}(\mathrm{\Theta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathrm{+}\mathit{C}\mathrm{\left(}\mathit{s,q}\mathrm{\right)})\exp (\mathrm{-}\mathit{s}\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathit{q}{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }})\mathit{,}$(45)where Θ(·) is the Heaviside step function and C(s,q) is an integration constant with respect to β. Since only momentum-loss processes are considered, there are no particles with momenta larger than their initial momentum p>p₀, corresponding to β<β₀, at any time. Therefore, because the function S(r,u) is independent of β, P(β,s,q) must vanish for β<β₀ implying C(s,q) = 0. Then, the solution of Eq. (43) reads $P\mathrm{\left(}\mathit{\beta ,s,q}\mathrm{\right)}\mathrm{=}\mathrm{\Theta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\exp (\mathrm{-}\mathit{s}\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathit{q}{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }})\mathit{.}$(46)Via an inverse Laplace transformation with respect to the variable q, one can recover the function P(β,s,u)$\mathit{P}\mathrm{\left(}\mathit{\beta ,s,u}\mathrm{\right)}\mathrm{=}{L}^{-1}\mathrm{\left[}P\mathrm{\right]}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{i}}{\mathrm{\int }}_{\mathit{c}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{i\infty }}^{\mathit{c}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{i\infty }}\exp \mathrm{\left(}\mathit{qu}\mathrm{\right)}P\mathrm{\left(}\mathit{\beta ,s,q}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{q}\mathit{.}$(47)Since P is well-defined and finite everywhere, one is free to choose the real-valued constant c = 0. Therefore, from using the relation $\mathit{\delta }\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{i}}{\mathrm{\int }}_{\mathrm{-}\mathrm{i\infty }}^{\mathrm{+}\mathrm{i\infty }}\exp (\mathit{w}\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{0}}\mathrm{)})\hspace{0.17em}\mathrm{d}\mathit{w}\mathit{,}$(48)it directly follows that $\mathit{P}\mathrm{\left(}\mathit{\beta ,s,u}\mathrm{\right)}\mathrm{=}\mathrm{\Theta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\exp (\mathrm{-}\mathit{s}\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)})\mathit{\delta }\left(\mathit{u}\mathrm{-}{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}\right)\mathit{.}$(49)Substituting the functions (41) and (49) into the integral in Eq. (38) leads to $\begin{array}{ccc}G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Theta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\exp (\mathrm{-}\mathit{s}\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)})}{\mathrm{(}\mathrm{4}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{u}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\exp \left(\mathrm{-}\frac{\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}\mathit{u}}\right)\mathit{\delta }\left(\mathit{u}\mathrm{-}{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}\right)\mathrm{d}\mathit{u}\mathit{.}& & \end{array}$(50)Because $\mathrm{0}\mathrm{\le }{{}^{\mathrm{\int }}\mathit{\beta }_{\mathrm{0}}}_{\mathit{\beta }}^{}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}\mathit{<}\mathrm{\infty }$, integration with respect to the variable u yields $G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Theta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\exp (\mathrm{-}\mathit{s}\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)})}{{\left(\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\exp \left(\frac{\mathrm{-}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}}\right)\mathrm{·}$(51)In order to obtain R(r,β,t) from G(r,β,s), another inverse Laplace transformation, now with respect to the variable s, is performed,$L^{-1}\mathrm{\left[}G\mathrm{\right]}\mathrm{\equiv }\mathit{R}\mathrm{\left(}{r}\mathit{,\beta ,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathrm{i}}{\mathrm{\int }}_{\mathit{c}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{i\infty }}^{\mathit{c}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{i\infty }}\exp \mathrm{\left(}\mathit{st}\mathrm{\right)}G\mathrm{\left(}{r}\mathit{,\beta ,s}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{s,}$(52)where again c = 0 is chosen, resulting in $\mathit{R}\mathrm{\left(}{r}\mathit{,\beta ,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Theta }\mathrm{(}\mathit{\beta }\mathrm{-}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}\mathit{\delta }\mathrm{(}\mathit{t}\mathrm{-}\mathit{\beta }\mathrm{+}{\mathit{\beta }}_{\mathrm{0}}\mathrm{)}}{{\left(\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\exp \left(\frac{\mathrm{-}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\int }}_{{\mathit{\beta }}_{\mathrm{0}}}^{\mathit{\beta }}\mathit{D}\mathrm{\left(}{\mathit{\beta }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\beta }}^{\mathrm{\prime }}}\right)\mathrm{·}$(53)The Green’s function becomes $\mathit{G}\mathrm{(}{r}\mathit{,p,t}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\frac{\mathrm{\Theta }\mathrm{(}{\mathit{p}}_{\mathrm{0}}\mathrm{-}\mathit{p}\mathrm{)}\mathit{\delta }\left(\mathit{t}\mathrm{+}{\mathrm{\int }}_{{\mathit{p}}_{\mathrm{0}}}^{\mathit{p}}\mathit{b}\mathrm{(}{\mathit{p}}^{\mathrm{\prime }}{\mathrm{)}}^{-1}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}{\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{·}{\left(\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\exp \left(\frac{\mathrm{-}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}}\right)\mathrm{·}$(54)The general solution n_p(r,p,t) of the transport Eq. (1), with a momentum-loss rate given by Eq. (12) and for an arbitrary source function Q, can be obtained by convolving the fundamental solution G(r,p,t | r₀,p₀) with a source term Q(r₀,p₀,t₀), ${\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}\mathrm{=}\$\mathit{G}\mathrm{(}{r}\mathit{,p,t}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathrm{\left)}\mathit{Q}\mathrm{\right(}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathrm{d}{\mathit{t}}_{\mathrm{0}}\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{\mathit{r}}_{\mathrm{0}}\hspace{0.17em}\mathrm{d}{\mathit{p}}_{\mathrm{0}}\mathit{.}$(55)This differential CR proton number density can also be applied to many other astrophysical situations with scalar, momentum-dependent diffusion and any type of momentum losses, such as stellar winds, CR diffusion in the interstellar medium, or gamma-ray bursts.

Appendix B: Specific source function for SNR-MC systems

The source function Q(r₀,p₀,t₀) is modeled for four specific SNRs associated with MCs showing gamma-ray emission for which data samples from spectral measurements in the X-ray energy range exist. Here, the source spectrum is assumed to be of the specific form $\begin{array}{ccc}\mathit{Q}\mathrm{\left(}{{r}}_{\mathrm{0}}\mathit{,}{\mathit{p}}_{\mathrm{0}}\mathit{,}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}& \mathrm{=}& {\mathit{Q}}_{\mathrm{norm}}\mathrm{·}{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{\mathit{p}}_{\mathrm{0}}\mathrm{\right)}\mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathrm{-}\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}]\mathrm{·}[\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)\mathrm{-}\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)]\mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{z}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathrm{)}\mathrm{-}\mathrm{\Theta }\mathrm{\left(}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}]\\ & & \end{array}$(56)where Q_norm denotes a normalization constant, Q_p(p₀) is the spectral shape of the low-energy CR protons in terms of the particle momentum p₀, and l_c characterizes the extent of the emission region. A source function of this type describes emission that is constant over a period of time (normalized to the unit time interval of one second) from a cubic emission volume, seen by an observer located at the center of a face of the emission volume that coincides with the cloud surface, with a coordinate system such that the positive Cartesian z-axis is normal to this face and points into the cloud. The cubic geometry is chosen over the more physical, spherical geometry for numerical feasibility. The volume of the cube-shaped emission region, ${\mathit{l}}_{\mathrm{c}}^{\mathrm{3}}$, is adapted to the spherical emission volume used in the modeling process of the gamma rays in Sect. 4. For the specific source function (56), the integrations with respect to r₀, p₀, and t₀, that have to be performed in order to determine the differential CR proton number density (55), $\begin{array}{ccc}{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}& \mathrm{=}& {\mathit{Q}}_{\mathrm{norm}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{R^{\mathrm{3}}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}\frac{\mathrm{\Theta }\mathrm{(}{\mathit{p}}_{\mathrm{0}}\mathrm{-}\mathit{p}\mathrm{)}\mathrm{·}{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{\mathit{p}}_{\mathrm{0}}\mathrm{\right)}\mathrm{·}\mathit{\delta }\left(\mathit{t}\mathrm{+}{\mathrm{\int }}_{{\mathit{p}}_{\mathrm{0}}}^{\mathit{p}}\mathit{b}\mathrm{(}{\mathit{p}}^{\mathrm{\prime }}{\mathrm{)}}^{-1}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}{\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{·}{\left(\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}\exp \left(\frac{\mathrm{-}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}}\right)\\ & & \mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathrm{-}\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}]\mathrm{·}[\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)\mathrm{-}\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)]\mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{z}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathrm{)}\mathrm{-}\mathrm{\Theta }\mathrm{\left(}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}]\\ & & \end{array}$(57)can be done separately. The specific expressions for the actual spectral shapes Q_p(p₀) and the normalization constants Q_norm used for the astrophysical objects of interest are of no relevance for these integrations. Since the Green’s function G(r,p,t | r₀,p₀) is independent of t₀, only the time-dependent factor of the source function, Q_time(t₀) = Θ(t₀) − Θ(t₀ − 1) has to be integrated with respect to t₀, yielding ${\mathit{I}}_{\mathrm{time}}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{\infty }}[\mathrm{\Theta }\mathrm{\left(}{\mathit{t}}_{\mathrm{0}}\mathrm{\right)}\mathrm{-}\mathrm{\Theta }\mathrm{(}{\mathit{t}}_{\mathrm{0}}\mathrm{-}\mathrm{1}\mathrm{)}]\hspace{0.17em}\mathrm{d}{\mathit{t}}_{\mathrm{0}}\mathrm{=}\mathrm{1.}$(58)The spatial integration $\begin{array}{ccc}{\mathit{I}}_{\mathrm{space}}& \mathrm{=}& \int R^{\mathrm{3}}\exp \left(\frac{\mathrm{-}\mathrm{(}{r}\mathrm{-}{{r}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{4}{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}}\right)\mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathrm{-}\mathrm{\Theta }\mathrm{(}{\mathit{x}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}]\mathrm{·}[\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)\mathrm{-}\mathrm{\Theta }\left({\mathit{y}}_{\mathrm{0}}\mathrm{-}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\right)]\\ & & \mathrm{·}[\mathrm{\Theta }\mathrm{(}{\mathit{z}}_{\mathrm{0}}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathrm{)}\mathrm{-}\mathrm{\Theta }\mathrm{\left(}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}]\hspace{0.17em}{\mathrm{d}}^{\mathrm{3}}{\mathit{r}}_{\mathrm{0}}\end{array}$(59)results in a product of error functions $\begin{array}{ccc}{\mathit{I}}_{\mathrm{space}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{8}}{\left(\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{·}\underset{\mathit{h}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{x,y,z}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}}{\mathrm{􏽙}}\left[\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{\mathrm{1}}\mathrm{erf}\left(\frac{{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{j}}\mathrm{·}\mathit{h}}{\sqrt{\mathrm{4}{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}}}\right)\right]\mathrm{·}& & \end{array}$(60)Then, one finds the following momentum integral $\begin{array}{ccc}{\mathit{I}}_{\mathrm{mom}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{\mathit{p}}_{\mathrm{0}}\mathrm{\right)}}{\mathit{b}\mathrm{\left(}\mathit{p}\mathrm{\right)}}\mathrm{·}\mathrm{\Theta }\mathrm{(}{\mathit{p}}_{\mathrm{0}}\mathrm{-}\mathit{p}\mathrm{)}\mathit{\delta }\left(\mathit{t}\mathrm{+}{\mathrm{\int }}_{{\mathit{p}}_{\mathrm{0}}}^{\mathit{p}}\mathit{b}\mathrm{(}{\mathit{p}}^{\mathrm{\prime }}{\mathrm{)}}^{-1}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\right)\mathrm{·}\underset{\mathit{h}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{x,y,z}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}}{\mathrm{􏽙}}\left[\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{\mathrm{1}}\mathrm{erf}\left(\frac{{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{j}}\mathrm{·}\mathit{h}}{\sqrt{\mathrm{4}{\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}}}\right)\right]\hspace{0.17em}\mathrm{d}{\mathit{p}}_{\mathrm{0}}\mathit{.}& & \end{array}$(61)In order to solve this integral, first, one has to explicitly evaluate the integrals in the arguments of the Dirac distribution and the error functions, respectively. Starting with the integral ${{}^{\mathrm{\int }}}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{/}\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}$, keeping the solution as general as possible, the momentum dependence of the diffusion coefficient is assumed to be $\mathit{D}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{0}}{\left(\frac{\mathit{p}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)}^{\mathit{k}}$(62)with D₀ = const. and values 1/3 ≤ k ≤ 4 /3. Using Eq. (12) and the dimensionless quantity p≡ p/ (m_pc), one obtains ${\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\frac{\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{D}}_{\mathrm{0}}\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}\mathit{c}}{{\mathit{a}}_{\mathrm{cc}}}\hspace{0.17em}{\mathrm{\int }}_p^{p_{\mathrm{0}}}\frac{{\left(p^{\mathrm{\prime }}\right)}^{\mathrm{2}\mathrm{+}\mathit{k}}}{\mathit{a}\mathrm{·}{\mathrm{\left(}{p}^{\mathrm{\prime }}\mathrm{\right)}}^{\mathrm{3}}\mathrm{+}\mathrm{1}}\hspace{0.17em}\mathrm{d}{p}^{\mathrm{\prime }}\mathit{,}$(63)where a = a_ad/a_cc. An analytic solution of this integral can be found in Gradshteyn & Ryzhik (1965, formula 3.914 number 5), yielding ${\mathrm{\int }}_{\mathit{p}}^{{\mathit{p}}_{\mathrm{0}}}\frac{\mathit{D}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{b}\mathrm{\left(}{\mathit{p}}^{\mathrm{\prime }}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{D}}_{\mathrm{0}}\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}\mathit{c}}{{\mathit{a}}_{\mathrm{cc}}}\hspace{0.17em}\left[\frac{{\left(p^{\mathrm{\prime }}\right)}^{\mathrm{3}\mathrm{+}\mathit{k}}{\hspace{0.17em}}_{\mathrm{2}}{\mathit{F}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{2}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{-}\mathit{a}\mathrm{·}{\left(p^{\mathrm{\prime }}\right)}^{\mathrm{3}}\mathrm{)}}{\mathrm{3}\mathrm{+}\mathit{k}}\right]\begin{array}{c}{p}_{\mathrm{0}}\\ p^{\mathrm{\prime }}\mathrm{=}p\end{array}\mathrm{·}$(64)Here, $\begin{array}{ccc}\mathrm{2}{\mathit{F}}_{\mathrm{1}}\mathrm{\left(}a\mathit{,}b\mathrm{;}c\mathrm{;}z\mathrm{\right)}\mathrm{\equiv }\frac{\mathrm{\Gamma }\mathrm{\left(}c\mathrm{\right)}}{\mathrm{\Gamma }\mathrm{\left(}b\mathrm{\right)}\mathrm{\Gamma }\mathrm{(}c\mathrm{-}b\mathrm{)}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\frac{{\mathit{t}}^{b\mathrm{-}\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{t}{\mathrm{)}}^{c\mathrm{-}b\mathrm{-}\mathrm{1}}}{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{t}\hspace{0.17em}z{\mathrm{)}}^a}\hspace{0.17em}\mathrm{d}\mathit{t}& & \end{array}$denotes the hypergeometric function and $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{\left(}c\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{t}}^{c\mathrm{-}\mathrm{1}}\exp \mathrm{(}\mathrm{-}\mathit{t}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{t}& & \end{array}$the complete Gamma function. Substituting this and ${\mathrm{\int }}_{{\mathit{p}}_{\mathrm{0}}}^{\mathit{p}}\mathit{b}\mathrm{(}{\mathit{p}}^{\mathrm{\prime }}{\mathrm{)}}^{-1}\hspace{0.17em}\mathrm{d}{\mathit{p}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}\mathit{c}}{\mathrm{3}\hspace{0.17em}{\mathit{a}}_{\mathrm{ad}}}\ln \left(\frac{{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{p}_{\mathrm{0}}^{\mathrm{3}}}{{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{p}^{\mathrm{3}}}\right)$(65)into Eq. (61), one finds $\begin{array}{ccc}{\mathit{I}}_{\mathrm{mom}}& \mathrm{=}& \frac{{\mathit{m}}_{\mathrm{p}}\mathit{c}}{\mathit{b}\mathrm{\left(}p\mathrm{\right)}}{\mathrm{\int }}_p^{\mathrm{\infty }}\mathit{\delta }(\mathit{t}\mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}\mathit{c}}{\mathrm{3}\hspace{0.17em}{\mathit{a}}_{\mathrm{ad}}}\ln \left(\frac{{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{p}_{\mathrm{0}}^{\mathrm{3}}}{{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{p}^{\mathrm{3}}}\right))\mathrm{·}{\mathit{Q}}_{\mathit{p}}\mathrm{\left(}{p}_{\mathrm{0}}\mathrm{\right)}\\ & & \end{array}$(66)The zeros of the argument of the Dirac distribution, as a function of the momentum p₀, are given by $p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}\mathrm{=}(p^{\mathrm{3}}\mathrm{·}\exp \left(\frac{\mathrm{3}{\mathit{a}}_{\mathrm{ad}}\mathit{t}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)\mathrm{+}\frac{{\mathit{a}}_{\mathrm{cc}}}{{\mathit{a}}_{\mathrm{ad}}}\left(\exp \left(\frac{\mathrm{3}{\mathit{a}}_{\mathrm{ad}}\mathit{t}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}}\right)\mathrm{-}\mathrm{1}\right){)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(67)Hence, the Dirac distribution can be written as $\begin{array}{ccc}\mathit{\delta }(\mathit{t}\mathrm{-}\frac{{\mathit{m}}_{\mathrm{p}}\mathit{c}}{\mathrm{3}\hspace{0.17em}{\mathit{a}}_{\mathrm{ad}}}\ln \left(\frac{{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{p}_{\mathrm{0}}^{\mathrm{3}}}{{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}{p}^{\mathrm{3}}}\right))\mathrm{=}\frac{\mathrm{(}{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}{)}^{\mathrm{3}}\mathrm{)}}{{\mathit{m}}_{\mathrm{p}}\mathit{c}\mathrm{·}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}{)}^{\mathrm{2}}}\mathrm{·}\mathit{\delta }(p_{\mathrm{0}}\mathrm{-}{p}_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)})\mathit{.}& & \end{array}$(68)Subsequently, the momentum integral becomes $\begin{array}{ccc}{\mathit{I}}_{\mathrm{mom}}& \mathrm{=}& \\ & & \mathrm{·}\underset{\mathit{h}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{x,y,z}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}}{\mathrm{􏽙}}\left[\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{\mathrm{1}}\mathrm{erf}\left(\frac{\mathrm{(}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{j}}\mathrm{·}{\mathit{h}}^{\mathrm{)}}\sqrt{{\mathit{a}}_{\mathrm{cc}}\mathrm{(}\mathrm{3}\mathrm{+}\mathit{k}\mathrm{)}}}{\sqrt{\mathrm{4}{\mathit{D}}_{\mathrm{0}}\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}\mathit{c}{\mathrm{[}{\mathrm{(}{p}^{\mathrm{\prime }}\mathrm{)}}^{\mathrm{3}\mathrm{+}\mathit{k}}{\hspace{0.17em}}_{\mathrm{2}}{\mathit{F}}_{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{2}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{-}\mathit{a}\mathrm{·}{\mathrm{(}{p}^{\mathrm{\prime }}\mathrm{)}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{]}}}_{p^{\mathrm{\prime }}\mathrm{=}p}^{p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}}}}\right)\right]\mathrm{·}\end{array}$(69)Note that $\mathrm{\Theta }(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}\mathrm{-}p)\mathrm{=}\mathrm{1}$ for all values of $p_{\mathrm{0}}^{\mathrm{zero}}$ and p. Then, combining of (58), (60) and (69) leads to the differential CR proton number density of a cubic emission source region with edge length l_c for all positions r inside the MC, with z ≥ 0, at any time t ≥ 0 and for all particle momenta p∈ [ 0.15, 0.86 ] ≤ p₀$\begin{array}{ccc}{\mathit{n}}_{\mathrm{p}}\mathrm{\left(}{r}\mathit{,p,t}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{Q}}_{\mathrm{norm}}\mathrm{·}\mathrm{(}{\mathit{a}}_{\mathrm{cc}}\mathrm{+}{\mathit{a}}_{\mathrm{ad}}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}{)}^{\mathrm{3}}\mathrm{)}\mathrm{·}{\mathit{Q}}_{\mathit{p}}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)})}{\mathrm{8}\mathrm{·}\mathit{b}\mathrm{\left(}p\mathrm{\right)}\mathrm{·}(p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}{)}^{\mathrm{2}}}\\ & & \mathrm{·}\underset{\mathit{h}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{x,y,z}\mathrm{+}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}}{\mathrm{􏽙}}\left[\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{\mathrm{1}}\mathrm{erf}\left(\frac{\mathrm{(}{\mathit{l}}_{\mathrm{c}}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{j}}\mathrm{·}{\mathit{h}}^{\mathrm{)}}\sqrt{{\mathit{a}}_{\mathrm{cc}}\mathrm{(}\mathrm{3}\mathrm{+}\mathit{k}\mathrm{)}}}{\sqrt{\mathrm{4}{\mathit{D}}_{\mathrm{0}}\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}\mathit{c}{\mathrm{[}{\mathrm{(}{p}^{\mathrm{\prime }}\mathrm{)}}^{\mathrm{3}\mathrm{+}\mathit{k}}{\hspace{0.17em}}_{\mathrm{2}}{\mathit{F}}_{\mathrm{1}}{\mathrm{(}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{2}\mathrm{+}\frac{\mathit{k}}{\mathrm{3}}\mathrm{;}\mathrm{-}\mathit{a}\mathrm{·}{\mathrm{(}{p}^{\mathrm{\prime }}\mathrm{)}}^{\mathrm{3}}\mathrm{)}}^{\mathrm{]}}}_{p^{\mathrm{\prime }}\mathrm{=}p}^{p_{\mathrm{0}}^{\mathrm{zero}}\mathrm{\left(}p\mathit{,t}\mathrm{\right)}}}}\right)\right]\mathrm{·}\end{array}$(70)