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 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 δ(β − β₀) = b(p₀)·δ(p − p₀), 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 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₀ 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 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 β, (35)where (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 , 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 (37)such that the full, non-trivial 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 well-known heat kernel of three-dimensional 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) = δ(β − β₀) /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 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 (46)Via an inverse Laplace transformation with respect to the variable q, one can recover the function P(β,s,u)(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 (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 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₀), (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 (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, , 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), (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 (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₀ = const. and values 1/3 ≤ k ≤ 4 /3. Using Eq. (12) and the dimensionless quantity p≡ p/ (m_pc), 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₀, 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₀(70)

© ESO, 2014