Issue 
A&A
Volume 567, July 2014



Article Number  A50  
Number of page(s)  18  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201423614  
Published online  09 July 2014 
Theoretical study of ionization profiles of molecular clouds near supernova remnants^{⋆}
Tracing the hadronic origin of GeV gamma radiation
^{1}
RuhrUniversität Bochum, Fakultät für Physik & Astronomie,
Institut für Theoretische Physik IV,
44780
Bochum,
Germany
email:
florian.schuppan@rub.de
^{2}
Universität Regensburg, Fakultät für Mathematik, 93040
Regensburg,
Germany
Received:
11
February
2014
Accepted:
18
March
2014
Context. Over the past few years, signatures of supernova remnants have been detected in gamma rays, particularly those that are known to be associated with molecular clouds. Whether these gamma rays are generated by cosmicray electrons emitting bremsstrahlung or experiencing inverse Compton scattering, or by cosmicray protons interacting with ambient hydrogen is usually not known. The detection of hadronic ionization signatures in spatial coincidence with gammaray signatures can help to unambiguously identify supernova remnants as sources of cosmicray protons.
Aims. Our central aim is to develop a method to investigate whether the gamma rays are formed by cosmicray protons. To achieve this goal, we derived the positiondependent cosmicrayinduced and photoinduced ionization rates.
Methods. To calculate hadronic signatures from cosmicrayinduced ionization to examine the origin of the observed gamma rays from molecular clouds associated with supernova remnants, we solved analytically the transport equation for cosmicray protons propagating in a molecular cloud, including the relevant momentumloss processes, and determined the proton flux at any penetration depth into the cloud.
Results. Because the solution of the transport equation is obtained for arbitrary source functions, it can be used for a variety of supernova remnants. We derived the corresponding theoretical ionization rate as a function of the penetration depth from the positiondependent proton flux, and compared it with photoinduced ionization profiles for available Xray data in a case study with four supernova remnants associated with molecular clouds. Three of the remnants show a clear dominance of the hadronically induced ionization rate, while for one remnant, Xray emission seems to dominate by a factor of 10.
Conclusions. This is the first derivation of positiondependent profiles for cosmicrayinduced ionization with an analytic solution for arbitrary cosmicray source spectra. The cosmicrayinduced ionization has to be compared with photoionization for strong Xray sources. For this purpose, measurements of Xray spectra from supernova remnant shocks in the subkeV to keV domain are necessary for a proper comparison. For sources dominated by cosmicrayinduced ionization (e.g., W49B), the ionization profiles can be used in the future to map the spatial structure of hadronic gamma rays and rotationvibrational lines induced by cosmicray protons. With instruments such as ALMA for the line signatures and CTA for the gammaray detection, this correlation study will help to identify sources of hadronic cosmic rays.
Key words: astroparticle physics / radiation mechanisms: nonthermal / cosmic rays / ISM: clouds / ISM: supernova remnants / gamma rays: ISM
Appendices are available in electronic form at http://www.aanda.org
© 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 gammaray 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 highenergy 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 lowenergy 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 crosssection 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 highenergy gamma rays and lowenergy 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 lowenergy CR protons, other possible ionization sources need to be ruled out, in particular UV photons and Xrays. 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). Xrays, 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 CRinduced ionization to enable one to attribute detectable ionization signatures to CRs (here, the term CRinduced 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 columndensity 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 gammaray 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 SNRMC systems, the positiondependent 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 positiondependent 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. Cosmicray 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 SNRMC systems is considered, and, using observed gammaray spectra, the proton fluxes at the gammaray emission regions are constructed assuming that the detected gamma radiation is caused by the decay of neutral pions formed in protonproton 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 gammaray 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 lowenergy CR protons are able to efficiently ionize the MCs, the fluxes are extrapolated toward lower energies under the assumption that they follow the same powerlaw 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 volumefilling 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, momentumdependent 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 (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 momentumdependent diffusion coefficient, Δ is the Laplace operator, b(p) is the momentumloss 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 (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_{0},p_{0},t_{0}) can be obtained by convolving the fundamental solution (2) with the source function (3)This differential CR proton number density can also be applied to many other astrophysical situations with scalar, momentumdependent diffusion and any type of momentum losses, such as stellar winds, CR diffusion in the interstellar medium, or gammaray bursts.
Two specific momentumloss processes in SNRMC 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 momentumloss timescale independent of the CR particle momentum). The loss rate for Coulomb collisions (Lerche & Schlickeiser 1982) can be given in the approximated form (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_{p}c, corresponding to a kinetic energy of 10 MeV, and a maximum momentum of 0.86 m_{p}c, 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 crosssection with increasing momentum (Rudd et al. 1985) and, on the other hand, because they can form pions by inelastic protonproton scattering (the momentum of ~0.86 m_{p}c corresponds to the threshold kinetic energy required for protons to produce pions by protonproton 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 (5)with (6)The loss rate for adiabatic deceleration at a shock velocity V can be expressed as (7)Assuming a constant, spherical expansion of the shock front, the shock velocity reads (8)where V_{0} is a constant, and R is the shock radius. Hence, the divergence of the shock velocity becomes (9)and the loss rate for adiabatic deceleration (7) yields (10)with (11)
Fig. 1 Comparison of the dynamical timescales for Coulomb losses, adiabatic deceleration, and diffusion with the typical values n_{H} = 100 cm^{3}, t_{age} = 10 kyr, and a diffusion length l = 2 pc. 

Open with DEXTER 
Note that for a constant shock velocity the quantity V_{0}/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 (longdashed, blue line for a young SNR and longshortdashed, orange line for a middleaged SNR), as well as the diffusion timescale (shortdashed, red line), as functions of the particle momentum p, are depicted. The various timescales are given by , , and T_{diff} = l^{2}/ (2D(p)), where l is the diffusion length, and the diffusion coefficient for the SNRMC systems is D(p) = 10^{28}·(p/ (m_{p}c))^{4/3} cm^{2} s^{1} with a powerlaw index of 4/3. Despite the specific choice of this powerlaw index, any powerlaw 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 nonrelativistic 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 nonrelativistic treatment, the assumption D(p) ∝ p^{4/3} leads to an overestimate of the diffusion effect for particles with momenta exceeding the nonrelativistic 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, IroshnikovKraichnan spectra can be applied as well, by adapting the powerlaw index of the diffusion coefficient. When GoldreichSridhar 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 middleaged SNRs, losses from Coulomb collisions for p ≲ 0.2 m_{p}c 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_{p}c, 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_{p}c lose all their kinetic energy due to adiabatic deceleration losses before reaching a penetration depth of 2 pc, while for middleaged SNRs, this only happens for lowmomentum particles with p ≲ 0.15 m_{p}c. The maximum penetration depth of the latter is governed by Coulomb losses.
Using both Eqs. (6) and (11), the total momentumloss rate can be expressed as (12)Moreover, the source function Q(r_{0},p_{0},t_{0}) is constructed such that it reflects the geometry of the SNRMC 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 gammaray emission and for which data samples from spectral measurements in the Xray 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 (13)where Q_{norm} denotes a normalization constant and Q_{p}(p_{0}) is the spectral shape of the lowenergy CR protons in terms of the particle momentum p_{0}, 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 zaxis is normal to this face and points into the cloud. The volume of the emission region, , 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.
Spectral parameters of the CR proton sources used to match the observed gammaray emission.
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_{p}c, 0.86 m_{p}c ] ≤ p_{0}(14)where a = a_{ad}/a_{cc}, k is the powerlaw index of the diffusion coefficient (for the data analysis, k = 4/3), erf(·) and denote the error function and the hypergeometric function, respectively, and For details, see Appendix B. In a next step, the undetermined normalization constant Q_{norm} and the momentumdependent spectral function Q_{p} are modeled such that they are adapted to the CR proton fluxes, derived from gammaray observations, in the SNRMC 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. CRinduced ionization profiles
The CRinduced 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), (15)where φ accounts for secondary ionization (i.e., ionization events induced by secondary electrons released during primary ionization events), d^{3}N_{p}/ (dE_{p} dA dt) is the differential CR proton flux, and is the direct ionization crosssection of molecular hydrogen by protons. The primary CR proton flux can be obtained with the help of gammaray observations and the differential CR proton number density (14). Treating the observed gammaray fluxes from SNRMC 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 gammaray 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 gammaray fluxes match the observational data. This is done by means of the protonproton interaction model given in Kelner et al. (2006), which describes gammaray emission for primary energies above 100 GeV, and with the delta distribution approximation for the pionproduction crosssection 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 gammaray 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 gammaray emission model, the speed of light c, yielding (16)The spectral shape function Φ_{p}(E_{p}) that is used to model the observed gammaray emission as described in more detail in Schuppan et al. (2012) is usually given by a set of typical expressions for SNRs detected at gammaray energies:

(a)
a single powerlaw in the kinetic energy Φ_{p}(E_{p}) = (E_{p}/E_{norm})^{− α},

(b)
a single powerlaw in the kinetic energy with an exponential cutoff Φ_{p}(E_{p}) = (E_{p}/E_{norm})^{− α}·exp^{(} − E_{p}/E_{c}^{)},

(c)
a double powerlaw in the kinetic energy Φ_{p}(E_{p}) = (E_{p}/E_{norm})^{− α}·^{(}1 + E_{p}/E_{br}^{)}^{α − αh},

(d)
a double powerlaw in the kinetic energy with an exponential cutoff Φ_{p}(E_{p}) = (E_{p}/E_{norm})^{− α}·^{(}1 + E_{p}/E_{br}^{)}^{α − αh}·exp^{(} − E_{p}/E_{c}^{)},
where , and E_{norm}, E_{c}, and E_{br} are normalization, cutoff, and break energies, respectively, α is the effective powerlaw index for energies below the break energy, and α_{h} is the effective powerlaw index above the break energy. Matching the different expressions for the CR proton fluxes to the observed gammaray 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 bestfit 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^{3}f_{V}/ 3 (to calculate the edge length l_{c} of the cubic emission regions), where f_{V} is a volumefilling factor (i.e., the homogeneous fraction (with hydrogen density n_{H}) of the volume of the SNRMC 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^{3} for W49B, V = 1.88 × 10^{4}pc^{3} for W44, V = 8.38 × 10^{3}pc^{3} for 3C 391, and V = 3.99 × 10^{4}pc^{3} for CTB 37A, and homogeneous hydrogen densities of n_{H} = 100 cm^{3} are used. Since the break energies and the cutoff energies of all these objects are large compared to the maximum kinetic energy considered, only the single powerlaw (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(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}/ (2 D(p)) is the corresponding diffusion timescale, and the relativistic particle speed, , 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}/ (2 D(p)), which is the mean time for a particle with momentum p to propagate the distance  r  in the case of standard threedimensional Gaussian diffusion, or via t = T_{rel}(  r  ,p) =  r  /v_{rel}(p). Moreover, the ionization rate is evaluated only along the zaxis (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 (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 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.
SNRMC systemEarth distances and radial SNR shock extensions.
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 (19)
The crosssection for direct ionization of molecular hydrogen by protons, , can be given via the parametric expression (Rudd et al. 1985) (20)with the lowenergy component and the highenergy component , where , a_{0} = 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_{p}c and p_{max} = 0.86 m_{p}c 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 GaussKronrod quadrature, but is more reliable and stable, particularly for rapidly changing integrands in onedimensional 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.
Bestfit parameters of the Xray flux models and their χ^{2} per degree of freedom.
5. Photoionization profiles
In order to determine whether there is a dominant process generating the ionization signatures of SNRMC 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 Xrays (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. Xrays, however, are capable of traversing larger columns of matter, with their soft (lowenergy) component losing energy faster than their hard (highenergy) component. In fact, hard Xrays might penetrate even farther into MCs than CRs (Tielens 2010). In this section, the Xray fluxes in the SNRMC 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 BeerLambert 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 Xray 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 Xrays. 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 Xray observations of SNRs associated with MCs provide nearEarth Xray flux data, which is inadequate for the derivation of the Xray fluxes at the MCs because, on the one hand, there are energydependent 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^{2}attenuation of the photon fluxes from the source along the path toward the instrument. Thus, the Xray fluxes at the cloud surfaces, F_{X  ray, cs}, can be derived from the Xray fluxes at the instrument, F_{X  ray, in}, by (21)where d is the distance between the SNRMC 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 Xray fluxes from interactions with traversed matter between the emission region and the detector is already taken into account in the data samples of the Xray observations F_{X  ray, in}. Since the photoinduced ionization rate is given, as in the case of CRinduced 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 Xray fluxes to the observational data samples. The parametric model function for the photon fluxes of two of the four SNRMC systems studied here, W44 and 3C 391, is chosen to be a triple powerlaw (22)where a, a_{1}, and a_{2} are dimensionless spectral indices, b_{1} and b_{2} denote break energies, respectively, and F_{0} is a constant. Owing to their significantly different Xray data samples, different parametric models for CTB 37A and W49B are used. These are, for CTB 37A, (23)and for W49B, (24)Applying these simpler parametric descriptions leads to an improved χ^{2}value per degree of freedom since a lesser number of fitting parameters is used. The bestfit parameters are given in Table 3. The extremely low value of χ^{2} per degree of freedom for the object W44 results from the large error bars of the Xray data, cf. Fig. 2b. For the objects W44, 3C 391, and CTB 37A, certain Xray data points were excluded because they are caused by line emission originating along the path between the nearEarth measuring instrument and the SNRMC 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 Xray data samples was performed by minimizing χ^{2} per degree of freedom based on the NelderMead method (Nelder & Mead 1965). This method is a nonlinear optimization technique that does not require the derivatives of χ^{2} in order to minimize χ^{2} with respect to the optimization parameters. It is particularly useful in multidimensional parameter space optimization. This routine is implemented as a C++ code via the GNU Scientific Library.
Since the total photoabsorption crosssection, 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 lowenergy photons on the photoionization rate, three descriptions of the lowenergy photon fluxes are considered:

(I)
extrapolation of the parametric model function used to fit theobserved Xray 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 cutoff below the lowest photon energy covered by the observational data, that is, no extrapolation,
where the lowest photon energies for which Xray 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 lowenergy photon fluxes are shown along with the observational Xray data samples in Fig. 2.
Fig. 2 Parametric models (I) (red, solid lines), (II) (blue, longdashed lines), and (III) fitted to Xray flux data samples (black dots and error bars) for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d). The Xray data samples used are taken from Ozawa et al. (2009) (Suzaku) for a), Harrus et al. (2006) (XMMNewton) for b), Chen & Slane (2001) (ASCA) for c), and Sezer et al. (2011) (Suzaku) for d). 

Open with DEXTER 
Fig. 3 CRinduced ionization profiles for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d), originating from diffusiondriven motion (black, solid lines) and relativistic, kinematic motion (cyan, longdashed lines), respectively. 

Open with DEXTER 
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 Xray spectra extrapolated toward lower energies by a powerlaw (red, shortdashed lines), by a constant (gray, longshortdashed lines), and without extrapolation (purple, longshortshortdashed lines), respectively. 

Open with DEXTER 
Fig. 5 Total ionization profiles (blue, solid lines) and ionization profiles induced by Xray spectra extrapolated toward lower energies by a powerlaw alone (red, dashed lines), respectively, for W49B a), W44 b), 3C 391c), and CTB 37A d). 

Open with DEXTER 
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 BeerLambert formula, (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_{H}z is the optical depth for a homogeneous cloud, σ_{pa}(E) is the total photoabsorption crosssection, 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 CRinduced 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 (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 crosssection per hydrogen nucleus is also provided in Maloney et al. (1996) as an empirical, broken powerlaw fit to experimental data (27)The breaks in the above expression come from contributions of heavy elements in the MCs to the total crosssection for which there is a certain threshold energy for shell absorption, namely oxygen Kshell absorption at 0.5 keV and iron Kshell 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 Xray 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 crosssection 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 lowenergy 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 diffusiondriven motion as well as purely relativistic, kinematic motion of the CR protons and their implications for the CRinduced ionization rate are shown in Fig. 3.
For the specific diffusion coefficient assumed here, the diffusiondriven 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 diffusiondriven 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 lowenergy extrapolations of the Xray data are presented, and compared with the ionization profiles induced by CRs. The powerlaw 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 cutoff Xray 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 crosssection is large, which means that lowenergy 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 crosssection 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 CRinduced ionization rate, more than two orders of magnitude higher than the corresponding ionization rate found for 3C 391, which has the secondhighest CRinduced 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 Xray 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 CRinduced 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 lowenergy 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 CRinduced 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 CRinduced 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 SNRMC 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, volumefilling factors of 1 are applied. While this assumption is common, a volumefilling factor ≪1 is more realistic. Taking this into consideration, the CR proton flux normalizations would increase, and so would the CRinduced 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 volumefilling factor of 0.06 was applied (Abdo et al. 2010b). Changing this factor to 1 would yield a CRinduced 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 Xray 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 lowenergy CR protons, indicating that highenergy CR protons are the driving force for gammaray emission. In the following, for a conservative comparison, only the extrapolation of the Xray fluxes toward lower photon energies by the powerlaw (I) is discussed. The total ionization profiles analyzed here are simply the sum of the photoinduced ionization profiles, with the powerlaw extrapolation (I), and the CRinduced 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 CRinduced 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.
Spatial resolutions of the studied SNRMC systems for a selection of relevant telescopes.
7. Conclusions and outlook
We calculated the positiondependent H_{2} ionization rates of molecular clouds near supernova remnants, induced by photons and cosmicray protons, respectively, and compared them to examine a potential hadronic formation scenario of the observed gamma radiation from these systems. The cosmicray proton fluxes were constructed by means of the analytic solution of the transport equation for cosmicray protons inside the clouds, while the photon fluxes were obtained using the BeerLambert 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, momentumdependent diffusion and any type of momentum losses, such as stellar winds, cosmicray diffusion in the interstellar medium, or gammaray bursts.
Enhanced ionization rates compared with the ionization rates induced by photons alone are expected if the gammaray emission seen from systems containing supernova remnants and molecular clouds originates from the decay of neutral pions that are formed in protonproton interactions of the cosmicray protons with ambient hydrogen. The method presented here is particularly promising for objects of low Xray fluxes and steep, incident cosmicray proton spectra, when the spectra can be extrapolated toward lower energies by single powerlaws. For the SNRMC system W49B, the ionization rate induced by cosmicray 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 cosmicray 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 cosmicray 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 rotationvibrational 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 rotationvibrational 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 rotationvibrational line intensities are inaccurate. Hence, the molecular ion H is especially wellsuited 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^{12}s^{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^{20} 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 rotationvibrational 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 targetselection 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 nextgeneration 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 farIR (FIR) telescopes, CTA measurements can be used to determine regions of enhanced gammaray 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 (see, e.g., Indriolo et al. 2010), it is possible to find the ionizationinduced signatures that are spatially correlated with the GeV gammaray emission, which will help pinpointing supernova remnants as sources of cosmicray 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 Xray data are currently available, because the Xray fluxes in the corresponding energy ranges are either below the detection limit or are not aligned with the gammaray emission regions. For these systems, precise measurements of Xray spectra by experiments such as Swift, XMMNewton, 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 cosmicray 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 Coulombloss coefficient and the ionization crosssection, for example by using the BetheBloch approximation (Bethe 1933; Padovani et al. 2009). The relative abundances of cosmicray 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 crosssection. 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 gammaray observations, correlations that can only be caused by a common primary component (see, e.g., Becker 2008) can be revealed.
Acknowledgments
We would like to thank R. Schlickeiser, J. H. Black, and S. Schöneberg for helpful and inspiring discussions. We acknowledge helpful comments and suggestions from the anonymous referee, which considerably improved the quality of the manuscript. We are also grateful to M. Ozawa for providing the Xray data from Chandra for W49B. Furthermore, we thank Hongquan Su for informing us about revised observational values for the distance of W44, CTB 37A, and 3C 391, and H. Fichtner for a careful reading of the manuscript. FS also thanks M. Mandelartz for helpful advice concerning programming.
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 706, L1 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, Science, 327, 1103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010b, ApJ, 722, 1303 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2006, A&A, 457, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, A&A, 483, 509 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F. A., & Atoyan, A. M. 2000, A&A, 362, 937 [NASA ADS] [Google Scholar]
 Axford, W. I., Leer, E., & Skadron, G. 1977, ICRC, 11, 132 [NASA ADS] [Google Scholar]
 Baade, W., & Zwicky, F. 1934, PNAS, 20, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, J. K. 2008, Phys. Rep., 458, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, J. K., Black, J. H., Safarzadeh, M., & Schuppan, F. 2011, ApJ, 739, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 1978a, MNRAS, 182, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 1978b, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Bethe, H. 1933, in Handbuch der Physik (Berlin: Springer), 24, 491 [Google Scholar]
 Black, J. H., & Dalgarno, A. 1977, ApJS, 34, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, Y., & Slane, P. O. 2001, ApJ, 563, 202 [NASA ADS] [CrossRef] [Google Scholar]
 Courant, R., & Hilbert, D. 2008, Methods of Mathematical Physics, Differential Equations, Methods of Mathematical Physics (Wiley) [Google Scholar]
 Cravens, T. E., & Dalgarno, A. 1978, ApJ, 219, 750 [NASA ADS] [CrossRef] [Google Scholar]
 CTA Consortium, & Hermann, G. 2011, Nucl. Phys. B Proc. Suppl., 212, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Dalgarno, A., Yan, M., & Liu, W. 1999, ApJS, 125, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Duhamel, H. 1838, in Mémoires présentés par divers savants à l’Académie des Sciences de l’Institut de France, ed. Académie des Sciences (Paris), 5, 440 [Google Scholar]
 Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] [Google Scholar]
 Galassi, M. 2009, GNU Scientific Library: reference manual for GSL version 1.12 (Network Theory) [Google Scholar]
 Goto, M., Indriolo, N., Geballe, T. R., & Usuda, T. 2013, J. Phys. Chem. A, 117, 9919 [CrossRef] [Google Scholar]
 Harrus, I., Smith, R., Slane, P., & Hughes, J. 2006, in The Xray Universe 2005, ed. A. Wilson, ESA SP, 604, 369 [Google Scholar]
 IceCube Collaboration 2013, Science, 342 [arXiv:1311.5238] [Google Scholar]
 Indriolo, N., Fields, B. D., & McCall, B. J. 2009, ApJ, 694, L257 [NASA ADS] [CrossRef] [Google Scholar]
 Indriolo, N., Blake, G. A., Goto, M., et al. 2010, ApJ, 724, 1357 [NASA ADS] [CrossRef] [Google Scholar]
 Jokipii, J. R. 1987, ApJ, 313, 842 [NASA ADS] [CrossRef] [Google Scholar]
 Kamae, T., Karlsson, N., Mizuno, T., Abe, T., & Koi, T. 2006, ApJ, 647, 692 [NASA ADS] [CrossRef] [Google Scholar]
 Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74, 034018 [NASA ADS] [CrossRef] [Google Scholar]
 Lepage, G. P. 1978, J. Comput. Phys., 27, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Lerche, I., & Schlickeiser, R. 1982, MNRAS, 201, 1041 [NASA ADS] [CrossRef] [Google Scholar]
 Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Mannheim, K., & Schlickeiser, R. 1994, A&A, 286, 983 [NASA ADS] [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, The Computer Journal, 7, 308 [CrossRef] [Google Scholar]
 O’Donnell, E. J., & Watson, W. D. 1974, ApJ, 191, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Ozawa, M., Koyama, K., Yamaguchi, H., Masai, K., & Tamagawa, T. 2009, ApJ, 706, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, M., & Galli, D. 2011, A&A, 530, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, M., Galli, D., & Glassgold, A. E. 2009, A&A, 501, 619 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Particle Data Group 2010, J. Phys. G: Nucl. Part. Phys., 37, 5021 [Google Scholar]
 Peck, A. B., & Beasley, A. J. 2008, J. Phys. Conf. Ser., 131, 2049 [NASA ADS] [CrossRef] [Google Scholar]
 Ptuskin, V. 2006, J. Phys. Conf. Ser., 47, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Rudd, M. E., Kim, Y.K., Madison, D. H., & Gallagher, J. W. 1985, Rev. Mod. Phys., 57, 965 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 1989a, ApJ, 336, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 1989b, ApJ, 336, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Schuppan, F., Becker, J. K., Black, J. H., & Casanova, S. 2012, A&A, 541, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sezer, A., Gök, F., Hudaverdi, M., & Ercan, E. N. 2011, MNRAS, 417, 1387 [NASA ADS] [CrossRef] [Google Scholar]
 Tielens, A. 2010, The Physics and Chemistry of the Interstellar Medium (Cambridge University Press) [Google Scholar]
 Voit, G. M. 1991, ApJ, 377, 158 [NASA ADS] [CrossRef] [Google Scholar]
 Wagner, R., & CTA Consortium 2010, J. Phys. Conf. Ser., 203, 2121 [NASA ADS] [CrossRef] [Google Scholar]
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_{0},p_{0},t_{0}), that is, the solution for the Dirac source distribution is determined.The fundamental transport equation is then given by (28)In terms of the function R(r,p,t) = b(p)·G(r,p,t), Eq. (28) becomes (29)With the new coordinate (30)induced by ∂_{β} = − b(p)∂_{p}, and δ(β − β_{0}) = b(p_{0})·δ(p − p_{0}), Eq. (29) yields (31)This partial differential equation can be solved by first applying a Laplace transformation L[·] with respect to the time variable t(32)Note that the lower integration limit in the definition of the Laplace transformation (32) fixes t_{0} = 0 without loss of generality. The Laplace transform of Eq. (31) reads (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_{0} of the MC. Then, Eq. (33) reduces to the inhomogeneous, linear partial differential equation of first order in β and second order in r, (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 higherdimensional space labeled with an additional auxiliary variable. This principle is applied here on a linear partial differential equation with a productseparable inhomogeneity of singlevariable factors with respect to the variables r and β, (35)where (36)From the structure of the lefthand side of Eq. (35), it directly follows that the solution of the corresponding homogeneous equation is productseparable G_{hom} = S(r)P(β,s). An embedding of the original (r,β,s)space into the extended, higherdimensional (r,β,s,u)space, where , implies that there is a family of functions S_{u}(r): = S(r,u) and P_{u}(β,s): = P(β,s,u), fulfilling the relations Ô_{r}S(r,u) = ∂_{u}S(r,u) and Ô_{β,s}P(β,s,u) = ∂_{u}P(β,s,u) with the boundary conditions (37)such that the full, nontrivial solution of Eq. (35) is given by the integral (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 (39)and (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 wellknown heat kernel of threedimensional Euclidean space (41)A solution of Eq. (40) can directly be found after performing a Laplace transformation with respect to the variable u(42)Then, Eq. (40) becomes (43)Note that the boundary condition P(β,s,u = 0) = δ(β − β_{0}) /D(β) was already fixed in (37). Using an integrating factor of the form , Eq. (43) can be rewritten as (44)and solved by simple integration, leading to (45)where Θ(·) is the Heaviside step function and C(s,q) is an integration constant with respect to β. Since only momentumloss processes are considered, there are no particles with momenta larger than their initial momentum p>p_{0}, corresponding to β<β_{0}, at any time. Therefore, because the function S(r,u) is independent of β, P(β,s,q) must vanish for β<β_{0} implying C(s,q) = 0. Then, the solution of Eq. (43) reads (46)Via an inverse Laplace transformation with respect to the variable q, one can recover the function P(β,s,u)(47)Since P is welldefined and finite everywhere, one is free to choose the realvalued constant c = 0. Therefore, from using the relation (48)it directly follows that (49)Substituting the functions (41) and (49) into the integral in Eq. (38) leads to (50)Because , integration with respect to the variable u yields (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,(52)where again c = 0 is chosen, resulting in (53)The Green’s function becomes (54)The general solution n_{p}(r,p,t) of the transport Eq. (1), with a momentumloss 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_{0},p_{0}) with a source term Q(r_{0},p_{0},t_{0}), (55)This differential CR proton number density can also be applied to many other astrophysical situations with scalar, momentumdependent diffusion and any type of momentum losses, such as stellar winds, CR diffusion in the interstellar medium, or gammaray bursts.
Appendix B: Specific source function for SNRMC systems
The source function Q(r_{0},p_{0},t_{0}) is modeled for four specific SNRs associated with MCs showing gammaray emission for which data samples from spectral measurements in the Xray energy range exist. Here, the source spectrum is assumed to be of the specific form (56)where Q_{norm} denotes a normalization constant, Q_{p}(p_{0}) is the spectral shape of the lowenergy CR protons in terms of the particle momentum p_{0}, 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 zaxis 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 cubeshaped emission region, , 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_{0}, p_{0}, and t_{0}, that have to be performed in order to determine the differential CR proton number density (55), (57)can be done separately. The specific expressions for the actual spectral shapes Q_{p}(p_{0}) 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_{0},p_{0}) is independent of t_{0}, only the timedependent factor of the source function, Q_{time}(t_{0}) = Θ(t_{0}) − Θ(t_{0} − 1) has to be integrated with respect to t_{0}, yielding (58)The spatial integration (59)results in a product of error functions (60)Then, one finds the following momentum integral (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 , keeping the solution as general as possible, the momentum dependence of the diffusion coefficient is assumed to be (62)with D_{0} = const. and values 1/3 ≤ k ≤ 4 /3. Using Eq. (12) and the dimensionless quantity p≡ p/ (m_{p}c), one obtains (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 (64)Here, denotes the hypergeometric function and the complete Gamma function. Substituting this and (65)into Eq. (61), one finds (66)The zeros of the argument of the Dirac distribution, as a function of the momentum p_{0}, are given by (67)Hence, the Dirac distribution can be written as (68)Subsequently, the momentum integral becomes (69)Note that for all values of 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_{0}(70)
All Tables
Spectral parameters of the CR proton sources used to match the observed gammaray emission.
Bestfit parameters of the Xray flux models and their χ^{2} per degree of freedom.
Spatial resolutions of the studied SNRMC systems for a selection of relevant telescopes.
All Figures
Fig. 1 Comparison of the dynamical timescales for Coulomb losses, adiabatic deceleration, and diffusion with the typical values n_{H} = 100 cm^{3}, t_{age} = 10 kyr, and a diffusion length l = 2 pc. 

Open with DEXTER  
In the text 
Fig. 2 Parametric models (I) (red, solid lines), (II) (blue, longdashed lines), and (III) fitted to Xray flux data samples (black dots and error bars) for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d). The Xray data samples used are taken from Ozawa et al. (2009) (Suzaku) for a), Harrus et al. (2006) (XMMNewton) for b), Chen & Slane (2001) (ASCA) for c), and Sezer et al. (2011) (Suzaku) for d). 

Open with DEXTER  
In the text 
Fig. 3 CRinduced ionization profiles for the objects W49B a), W44 b), 3C 391 c), and CTB 37A d), originating from diffusiondriven motion (black, solid lines) and relativistic, kinematic motion (cyan, longdashed lines), respectively. 

Open with DEXTER  
In the text 
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 Xray spectra extrapolated toward lower energies by a powerlaw (red, shortdashed lines), by a constant (gray, longshortdashed lines), and without extrapolation (purple, longshortshortdashed lines), respectively. 

Open with DEXTER  
In the text 
Fig. 5 Total ionization profiles (blue, solid lines) and ionization profiles induced by Xray spectra extrapolated toward lower energies by a powerlaw alone (red, dashed lines), respectively, for W49B a), W44 b), 3C 391c), and CTB 37A d). 

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