ANALYTIC SOLUTION FOR SELF-REGULATED COLLECTIVE ESCAPE OF COSMIC RAYS FROM THEIR ACCELERATION SITES

Traditionally, the transport of CRs is treated differently inside and outside a DSA accelerator. Outside, CRs are thought to escape as TPs with a diffusion coefficient inferable from observations. This picture cannot be correct near the accelerator because the CR transport must be in a self-confinement regime, in which the CR streaming instability diminishes their diffusion coefficient by orders of magnitude. This CR anisotropy instability, and particularly its nonresonant extensions, macroscopically driven by the CR current and the CR pressure gradient, has a potential to generate magnetic field fluctuations well in excess of the ambient magnetic field, δB ≳ B₀ (e.g., Bell 2004; Drury & Falle 1986; Bykov et al. 2009; Malkov et al. 2010). At the very least, such strong fluctuations should justify the standard DSA assumption about the Bohm diffusion regime with the mean free path (mfp) of the order of the particle gyroradius, r_g, achieved at δB ∼ B₀. Naturally, the transport is isotropic in this regime. The ISM background turbulence, on the other hand, is much weaker $\delta B^{2}/B_{0}^{2}\lesssim 10^{-5}$ at the CR relevant length scales, thus resulting in the CR mfp ≳ 10⁵r_g (for GeV particles, r_g ≳ 10¹² cm). It is important to emphasize that under these circumstances the cross-field diffusion coefficient, κ_⊥, is suppressed by a large factor compared to the diffusion along the field line, κ_∥, i.e., κ_⊥/κ_∥ ∼ (δB/B₀)⁴ (see, e.g., Drury 1983).

It follows then that there is a problem of describing particle transport between the self-confinement (accelerator vicinity) and the TP (far from accelerator) transport regimes. To circumvent this problem, the acceleration process has been treated separately from the particle escape using one of two devices: the upper cutoff momentum and the free-escape boundary (FEB; see, e.g., Reville et al. 2009; Drury 2011, for recent discussions). As the names suggest, accelerated particles escape instantly upon reaching a prescribed boundary either in momentum or in configuration space. Their escape is assumed to have no effect on the acceleration other than through the modification of the shock structure (Moskalenko et al. 2007). In a simplified visualization of the DSA as a "box" process, suggested by Drury et al. (1999), the upper cutoff and the FEB are two sides of the box in the particle phase space. There were efforts to include self-generated waves into the description of particle acceleration and escape from strongly modified CR shocks (Malkov et al. 2002; Malkov & Diamond 2006), or in numerical treatments (Galinsky & Shevchenko 2007, 2011; Fujita et al. 2011). In most approaches, however, a sudden jump in the CR diffusion coefficient in momentum space is introduced to set an upper cutoff (e.g., Ptuskin & Zirakashvili 2005; Yan et al. 2012). A similar jump in coordinate space would result in an FEB.

As long as the CR transport outside the accelerator is treated in the TP approximation, no smooth transition between the CR acceleration and their escape is provided, regardless of the scenario for the latter. However, the diffusion coefficient rises by roughly five orders of magnitude in a transition zone that should be correspondingly large, unlike an infinitely thin FEB. It is this zone where the CR escape flux and confinement time are set by self-generated waves, thus rendering the FEB a rather implausible concept. For many SNRs it is then not yet possible to conclude whether the gamma emission is leptonic or hadronic, should such a conclusion depend on the CR escape and on the subsequent illumination of adjacent clouds. In a broader sense, there is a missing link in galactic CR generation between the CR acceleration (under a very strong self-generated wave–particle scattering) and their subsequent propagation (in a very weak interstellar turbulence). The goal of this paper is to establish such a link.

Our treatment below is applicable to the following two situations. In one situation, a shock accelerates particles continuously but some of them reach far enough or diffuse across the local field to become disconnected from the shock front and have thus chances to escape. While doing so, they drive their own waves at a gyroradius scale, whose amplitudes gradually decrease outward and so does the particle density. The second situation is a clear-cut case for this paper that deals with a CRC released into the ISM with no ongoing acceleration inside the CRC. Both situations correspond to a mesoscale transport regime, intermediate between the Bohm diffusion inside the accelerator and a global, quasi-isotropic CR diffusion in the Galaxy. The latter process is, in turn, supported by a random large-scale component of the galactic magnetic field at the scale of tens of parsecs superposed on the regular, spiral-arm-aligned toroidal field (Brunetti & Codino 2000; Strong et al. 2007; Blasi & Amato 2012). In the mesoscale transport regime, however, these components are considered as an ambient field with a scale much larger than the CR gyroradius and even larger than the CRC scale height. Before we describe in the next section the physical setting of the problem, it is worthwhile to emphasize some of the qualitative results.

First of all, for sufficiently high initial CR pressure (higher than magnetic pressure) and low ambient turbulence level, δB ≪ B₀, the spreading of the CRC strongly deviates from the oft-used TP solution. Instead of being controlled by only one (diffusive) scale $l_{{\rm D}}\sim \sqrt{D_{{\rm ISM}}t}$, the structure of the expanding CRC is more complicated. Namely, it comprises three zones, of which only the outermost has the above TP scaling, but a significantly lower CR pressure than it would maintain by diffusing through the ISM with the diffusivity D_ISM. Conversely, the CR pressure in the innermost part remains strongly enhanced. These two zones are connected by a self-similar, 1/z-part of the CR pressure profile.

Perhaps the most systematic quasi-linear (QL) derivation of equations for the coupled evolution of CRs and Alfvén waves is given by Skilling (1975). We use them in a simplified one-dimensional (along the ambient field, z-coordinate, Figure 1) form equivalent to that used by, e.g., Bell (1978) and, with an interesting "beyond QL" interpretation, by Drury (1983),

$$\begin{equation} \frac{d}{dt}P_{{\rm {\rm CR}}}\left(p\right)=\frac{\partial }{\partial z}\frac{\kappa _{{\rm B}}}{I}\frac{\partial P_{{\rm CR}}}{\partial z} \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{d}{dt}I=-C_{{\rm A}}\frac{\partial P_{{\rm CR}}}{\partial z}-\Gamma I. \end{equation} \tag{ 2 }$$

Here, C_A is the Alfvén velocity, and the time derivative is taken along the characteristics of unstable Alfvén waves, forward propagating with respect to the flow of speed U:

$$\begin{equation} \frac{d}{dt}=\frac{\partial }{\partial t}+\left(U+C_{{\rm A}}\right)\frac{\partial }{\partial z}. \end{equation} \tag{ 3 }$$

Equation (1) above is essentially a well-known convection–diffusion equation, written for the dimensionless CR partial pressure P_CR instead of their distribution function f(p, t). We have normalized it to the magnetic energy density $\rho C_{{\rm A}}^{2}/2$:

$$\begin{equation} P_{{\rm CR}}=\frac{4\pi }{3}\frac{2}{\rho C_{{\rm A}}^{2}}vp^{4}f, \end{equation} \tag{ 4 }$$

where v and p are the CR speed and momentum, and ρ is the plasma density. The total CR pressure is normalized to dln p, similarly to the wave energy density I:

$$\frac{\langle \delta B^{2}\rangle }{8\pi }=\frac{B_{0}^{2}}{8\pi }\int I\left(k\right)d\ln k=\frac{B_{0}^{2}}{8\pi }\int I\left(p\right)d\ln p.$$

Equation (2) is a wave kinetic equation in which the energy transferred to the waves equals the difference between the total work done by the particles, (U + C_A)∇P_CR, and the work done on the fluid, U∇P_CR (Drury 1983). This interpretation of the wave generation indicates that it operates in a maximum efficiency regime. A formal QL derivation of this equation assumes that the particle momentum p is related to the wave number k by the "sharpened" resonance condition kp = eB₀/c instead of the conventional cyclotron resonance condition kp_|| = eB₀/c (note that here k = k_∥; Skilling 1975). We have included only the linear wave damping Γ and we will return to the possible role of nonlinear saturation effects later. We assume that ∂P_CR/∂z ⩽ 0 at all times, so that only the forward propagating waves are unstable. The latter inequality is ensured by the formulation of initial value problem in each of the following two settings mentioned earlier in this section. In the first setting, the initial distribution of the CRC is symmetric with respect to z = 0, so we can consider their escape into the half-space z > 0. The second setting is when the cloud is limited from the left by a reflecting wall (CD). The appropriate boundary condition is ∂P_CR/∂z = 0 at z = 0 in both cases.

Restricting our consideration to the case of coordinate-independent damping rate Γ, we obtain the following ("QL") integral of the system of Equations (1) and (2):

$$\begin{equation} P_{{\rm CR}}(z,t)=P_{{\rm CR}0}(z^{\prime })-\frac{\kappa _{{\rm B}}}{C_{{\rm A}}}\frac{\partial }{\partial z}\ln \frac{I(z,t)}{I_{0}(z^{\prime })}. \end{equation} \tag{ 5 }$$

Here, P_CR0(z) and I₀(z) are the initial distributions of the CR partial pressure and the wave energy density, respectively, and z' = z − (U + C_A)t.

Using the QL integral, the system of Equations (1) and (2) can be reduced to one nonlinear convection–diffusion equation for, e.g., wave intensity I(z, t). We will use dimensionless variables measuring the distance z in units of a_∥, which is the initial size of the CRC along the field line. The time unit is then a/C_A. Note, however, that, as the acceleration is assumed to be inactive, the particle momentum p enters the problem only as a parameter but, the initial scale height of the CRC a generally depends on p. It is also convenient to introduce a new variable for the wave energy density

$$\begin{equation} W=\frac{C_{{\rm A}}a\left(p\right)}{\kappa _{{\rm B}}\left(p\right)}I \end{equation} \tag{ 6 }$$

and similar relations for I₀ and W₀. The equation for W then takes the following form:

$$\begin{equation} \frac{\partial W}{\partial t}-\frac{\partial }{\partial z}\frac{1}{W}\frac{\partial W}{\partial z}=-\frac{\partial }{\partial z}\mathcal {P}_{0}\left(z\right). \end{equation} \tag{ 7 }$$

We have neglected the wave propagation along the characteristics given by Equation (3) and the wave damping, assuming that these processes are slower than the CRC diffusion: U, C_A, aΓ ≪ κ_B/aI ∼ cr_g/aI. When the linear damping is important (e.g., in the case of Goldreich–Shridhar cascade, Farmer & Goldreich 2004), it can be easily incorporated into the current treatment by a simple change of variables indicated in Section 4. The function $\mathcal {P}$ is defined similarly to W above,

$$\begin{equation} \mathcal {P}=\frac{C_{{\rm A}}a}{\kappa _{{\rm B}}\left(p\right)}P_{{\rm CR}}, \end{equation} \tag{ 8 }$$

and, again, the index 0 refers in Equation (7) to the initial CR distribution P_CR0(z). Thus, according to Equation (5), for the dimensionless CR partial pressure we have

$$\begin{equation} \mathcal {P}=\mathcal {P}_{0}-\frac{\partial }{\partial z}\ln \frac{W}{W_{0}}. \end{equation} \tag{ 9 }$$

This quantity is governed by the equation

$$\begin{equation} \frac{\partial \mathcal {P}}{\partial t}=\frac{\partial }{\partial z}\frac{1}{W}\frac{\partial \mathcal {P}}{\partial z}, \end{equation} \tag{ 10 }$$

but its solution can be written down using the integral given by Equation (9), after the solution to Equation (7) is obtained.

While letting C_A → 0 in the wave–particle collective propagation, we utilize the finiteness of C_A in determining the boundary condition for Equation (7) at z = 0 as follows. Returning from the nonlinear diffusion Equation (7) to Equation (2), one may see that the wave energy does not actually "diffuse" but it is generated locally by particles that diffuse. Had the neglected wave diffusion spread the waves to the point of their stability, viz., z = 0, they would be convected away from it by virtue of C_A > 0. Recalling the symmetry (reflection) boundary condition ∂P_CR/∂z = 0 at z = 0 (no wave generation), we thus set a fixed boundary value W(0, t) = W₀, where W₀ is an initial (small) wave noise. It is instructive to investigate the case in which the initial noise is the same throughout the entire half-space z > 0 (this limitation can be straightforwardly relaxed), so that the waves will be generated entirely by escaping particles, thus emphasizing the self-confinement. The second boundary condition is set by W → W₀ for z → ∞, and all t < ∞. Note that in general, W₀ = W₀(p). These conditions determine the boundary value problem given by Equation (7) completely. However, we are interested primarily in diffusion of CRs outside the region of their initial localization, viz., at z > 1, where the source term in this equation vanishes. Therefore, the problem given by Equation (7) can be split into two separate problems, one in 0 ⩽ z ⩽ 1 and another in 1 < z < ∞ domain with the following junction conditions:

$$\begin{equation} \frac{\partial }{\partial z}\ln W\left|_{1-}^{1+}=W\right|_{1-}^{1+}=0. \end{equation} \tag{ 11 }$$

These are the continuity conditions for both the wave energy density and pressure across the edge of the initial CR localization.

The nonlinear boundary value problem in the region z > 1 given by Equations (7) and (11) describes the CR propagation outside the region of their initial localization. This problem can be solved exactly using the following self-similar substitution:

$$\begin{equation} W=\frac{1}{t^{\alpha }}w\left(\zeta \right),\;\;\;{\rm where}\;\zeta =\frac{z}{t^{\beta }}. \end{equation} \tag{ 12 }$$

Submitting this to Equation (7) with $\mathcal {P}_{0}=0$ yields α = 2β − 1. The boundary condition at infinity (W → W₀ ≠ 0, z → ∞), on the other hand, requires α = 0, so that the equation for w reads

$$\begin{equation} \frac{d}{d\zeta }\frac{1}{w}\frac{dw}{d\zeta }+\frac{\zeta }{2}\frac{dw}{d\zeta }=0, \end{equation} \tag{ 13 }$$

with $\zeta =z/\sqrt{t}$. It is interesting to observe that this equation has a simple special solution w = 2/ζ² = 2t/z². This solution describes a stationary particle distribution $\mathcal {P}\propto 1/z$ that is obviously related to the Bell (1978) asymptotic particle distribution in self-excited waves far away from a shock. In both cases, this is a singular limit of the problem as it implies a zero background turbulence level, W₀ = 0, and thus the number of particles in z > 0 half-space is infinite (Lagage & Cesarsky 1983; Drury 1983). Physically, this solution is different from what we find below in that it requires a permanent source of CRs at the origin so that their flux steadily drives waves that linearly grow in time.⁸ As we shall see, this simple special solution is an important singular limit of a more general and completely regular solution.

The general solution to Equation (13) can be found in quadratures by swapping ζ and w as dependent and independent variables and introducing an auxiliary function V(w) by the following substitution:

$$\begin{equation} \frac{dw}{d\zeta }=-\sqrt{2}w^{3/2}V\left(w\right). \end{equation} \tag{ 14 }$$

The equation for V(w) takes the following form:

$$\begin{equation} w\frac{d}{dw}w\frac{dV}{dw}+\frac{1}{4}\left(\frac{1}{V}-V\right)=0. \end{equation} \tag{ 15 }$$

This equation can be easily integrated, so that the function w(ζ) can be subsequently found from Equation (14). The first integral of Equation (15) is as follows:

$$\begin{equation} \left(\frac{\partial V}{\partial \ln w}\right)^{2}=R(V)\equiv \frac{1}{4}(V^{2}-2\ln V-q), \end{equation} \tag{ 16 }$$

where q is an integration constant. From Equation (14) and from the boundary condition w → W₀ at ζ → ∞ we infer that V(W₀) = 0, so that from Equation (16) we find

$$\begin{equation} w=W_{0}e^{\int _{0}^{V}\frac{dV^{\prime }}{\sqrt{R\left(V^{\prime }\right)}}}. \end{equation} \tag{ 17 }$$

The function w(ζ) is then determined by the following relation for ζ(V), Equation (14):

$$\begin{equation} \zeta =-\frac{1}{\sqrt{2W_{0}}}\int \frac{dV}{V\sqrt{R\left(V\right)}}e^{-\frac{1}{2}\int _{0}^{V}\frac{dV^{\prime }}{\sqrt{R\left(V^{\prime }\right)}}}. \end{equation} \tag{ 18 }$$

Using the identity 1/V = V − 2dR/dV and integrating Equation (18) by parts, after some manipulations the last equation simplifies considerably:

$$\begin{equation} \zeta =\sqrt{\frac{2}{w}}(2\sqrt{R}+V). \end{equation} \tag{ 19 }$$

It is useful to introduce an auxiliary constant V₁ related to the constant q and being the smaller of the two roots of R(V), which is R(V₁) = 0. The requirement of a real zero of the function R(V) ensures mapping two "copies" of the domain of V onto the full range 0 < ζ < ∞, whereby ζ → ∞ for V → 0. Indeed, ζ(V) in Equation (19) diverges at V = 0 as ${\sim }\sqrt{-\ln V}$. However, increasing V from V = 0 to V = V₁ does not cover the full range 0 < ζ < ∞ yet, as may be seen from Equation (19); ζ decreases from ∞ only to $\zeta _{1}\equiv V_{1}\sqrt{2/w\left(V_{1}\right)}>0$. To continue the integral curve to ζ = 0, it is necessary to switch branches of $\sqrt{R}$ at V = V₁ in Equations (17) and (19) and so continue the integral in Equation (17) back to V < V₁ (along the second copy of the V-domain). Decreasing then V from V = V₁ to V = V₀ ≡ exp (− q/2) brings the integral curve to the point ζ = 0, w = w_max = w(ζ = 0) (Figure 2). The explicit formal representation of the general solution of Equation (13) is given in Appendix A along with a numerical example (Figure 3). Equations (17) and (19) determine the solution w(ζ), where the integration constant V₀ (or q) should be obtained from the matching condition, Equation (11). This can be done by considering Equation (7) inside the region of initial CR localization, i.e., 0 < z < 1, where its right-hand side is nonzero.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given an initial distribution $\mathcal {P}_{0}\left(z\right)$, Equation (7) can always be integrated numerically in the finite domain 0 < z < 1, so that the full solution will be obtained from the boundary conditions and from the results of the previous section. However, according to Equation (10), already for t ≳ W₀, where W₀ ≪ 1, the initial profile $\mathcal {P}_{0}\left(z\right)$ will be redistributed over the unity interval in such a way as to approach a quasi-steady state in which the flux $W_{0}^{-1}\partial W/\partial z$ in Equation (7) through z = 0 will be balanced by the source integral, $\mathcal {P}_{0}\left(0\right)$ (recall that $\mathcal {P}_{0}\left(1\right)=0$). The particle flux through the z = 1 boundary decays as t^−1/2 with time (since |∂w/∂ζ| < ∞ at ζ → 0 +; see the preceding section or Equation (23) below). Therefore, for the self-similar stage of the cloud relaxation we can write

$$\begin{equation} \ln \frac{W}{W_{0}}=\int _{0}^{z}\mathcal {P}_{0}dz-B\left(t\right)z, \end{equation} \tag{ 20 }$$

where the "integration constant" B depends slowly (in the above sense) on time. Using the first matching condition in Equation (11), i.e., the CR pressure continuity, we can specify B(t) as follows: B∝t^−1/2 → 0 as t → ∞ (see Equation (23) below). This determines the self-similar (outer, z > 1) solution by the second matching condition in Equation (11):

$$\begin{equation} w\left(0\right)\equiv \lim _{t\rightarrow \infty }w\left(\frac{1}{\sqrt{t}}\right)=w_{{\rm max}}=W_{0}e^{\Pi }. \end{equation} \tag{ 21 }$$

Here we have introduced the integrated partial pressure as follows:

$$\begin{equation} \Pi =\int _{0}^{1}\mathcal {P}_{0}dz. \end{equation} \tag{ 22 }$$

The function B(t) and thus the internal solution W(z, t) in Equation (20) may be obtained using this equation, the first matching condition in Equation (11), and Equation (14),

$$\begin{equation} B(t)={-}\frac{1}{\sqrt{t}}\frac{\partial }{\partial \zeta }\ln w\left|_{\zeta =1/\sqrt{t}}=\sqrt{2w/t}V\right|_{\zeta =1/\sqrt{t}}. \end{equation} \tag{ 23 }$$

Note that the particle pressure at z > 1 is completely determined by the turbulence level w, Equation (9).

Now we can determine the integration constant q introduced in the preceding section. From Equations (21) and (A3) we obtain the following equation for q:

$$\begin{equation} \int _{0}^{V_{1}}dV/\sqrt{R\left(V\right)}+\int _{V_{0}}^{V_{1}}dV/\sqrt{R\left(V\right)}=\Pi , \end{equation} \tag{ 24 }$$

where R(V₁) = 0, R(V₀) = −V₀/2, while R(V) is given by Equation (16). In the most interesting case q ≈ 1, it is convenient to use the constant ε² = (q − 1)/2 in place of q. In this case, Π ≫ 1, ε ≪ 1, V₁ ≈ 1 − ε, and

$$\begin{eqnarray*} \Pi =2\int _{0}^{V_{1}}\frac{dV}{\sqrt{R}} &\approx & 2^{3/2}\int _{0}^{1-\varepsilon }\frac{dV}{\sqrt{\left(V-1\right)^{2}-\varepsilon ^{2}}} \nonumber\\ &\approx & 2^{3/2}\ln \frac{2}{\varepsilon }+\mathcal {O}\left(1\right), \end{eqnarray*}$$

so that the turning point of the solution at V = V₁ approaches the critical point V = 1, where R has a minimum (Figure 2)

$$V_{1}=1-2e^{-\Pi /2\sqrt{2}}.$$

For Π ≪ 1, we can write instead, $V_{1}\approx V_{0}(1+V_{0}^{2}/2)$ where V₀ = exp (− q/2).

To conclude this section, the self-similar expansion of the CRC described by Equations (7)–(10) is controlled by two parameters. One parameter is the background turbulence level W₀(p) in the media into which the cloud expands. The second parameter is the integrated pressure of the cloud, Π. Although the initial wave energy density inside the cloud (z < 1) is likely to be higher than W₀, we have adopted the background value W₀ also inside the cloud for simplicity, as this should not influence the self-similar CR propagation outside the cloud. In addition, efficient wave generation by a dense expanding CRC renders the initial value for the wave energy density inside the CRC unimportant.

The most concise characterization of the escape may be given by looking at the following two parts of the (conserved) integral particle pressure:

$$\Pi =\Pi _{0}+\Pi _{1}={\rm const},$$

where

$$\Pi _{0}=\int _{0}^{1}\mathcal {P}dz\;{\rm and}\;\Pi _{1}=\int _{1}^{\infty }\mathcal {P}dz$$

refer to the regions inside and outside of the initial CRC, respectively. Recall that at t = 0, Π₀ = Π and Π₁ = 0, while at t = ∞, an opposite CR distribution is reached: Π₀ = 0 and Π₁ = Π. Note that Equation (21) specifies the total work ultimately done by particles on waves, while they diffuse from 0 < z < 1 to 1 ⩽ z < ∞. As Π is conserved, it is natural to define the CR confinement time as the time at which Π₀ (Π₁) drops (raises) to a half of Π. Substituting W from Equation (12) into Equation (9), we obtain

$$\Pi _{1}\left(t\right)=\left.\ln \frac{w}{W_{0}}\right|_{\zeta =1/\sqrt{t}}.$$

For the half-life time t = t_1/2 we may use the following equations:

$$\begin{equation} \Pi _{1}(t_{1/2})=\Pi _{0}(t_{1/2})=\frac{1}{2}\Pi. \end{equation} \tag{ 25 }$$

In the simple TP case (Π ≪ 1), the half-life time amounts to (see Appendix B)

$$\begin{equation} t_{1/2}\approx \frac{1}{4\sigma }W_{0} \end{equation} \tag{ 26 }$$

with σ ≈ 0.23. In the opposite case Π ≫ 1, for t_1/2 we obtain

$$\begin{equation} t_{1/2}\approx t_{1}=1/\zeta _{1}^{2}\equiv \frac{w_{1}}{2V_{1}^{2}}\approx \frac{w_{1}}{2} \end{equation} \tag{ 27 }$$

(see Equation (19)). Note that for Π ≫ 1, V₁ ≈ V_1/2 (see Appendix B), and for the half-life time we obtain

$$\begin{equation} t_{1/2}=\frac{1}{2}W_{0}e^{\Pi /2}. \end{equation} \tag{ 28 }$$

It is clear that the nonlinear delay factor exp (Π/2) slows down the escape considerably, compared to the TP solution.

Another important characteristic of the particle escape is the spatial width of their distribution. Formally, the self-similar solution is scale-invariant ∼2/ζ, for Π ≫ 1 over the most part of the spatial distribution of the partial pressure (see below). Therefore, to characterize the width we use the point ζ_1/2, related to the half-time t_1/2, as follows:

$$\zeta _{1/2}=\frac{1}{\sqrt{t_{1/2}}}=\sqrt{\frac{2}{W_{0}}}e^{-\Pi /4}.$$

In physical coordinate z, this point moves outward as $z_{1/2}=\zeta _{1/2}\sqrt{t}$ starting from z = 1 at t = t_1/2. The expansion rate of a dense CRC (Π ≫ 1) is thus exponentially low.

Considering the spatial distribution of the spreading CRC in detail, it is convenient to start with the region far away from the source, where W → W₀. This asymptotic behavior corresponds to the case V ≪ 1. Introducing an auxiliary function U(ζ),

$$\begin{equation} U=-\frac{\partial }{\partial \zeta }\ln w=\sqrt{2w}V, \end{equation} \tag{ 29 }$$

that is related to the particle pressure through $\mathcal {P}=U\left(\zeta \right)/\sqrt{t}$, and using Equations (16) and (19), we obtain for U the following expression:

$$\begin{equation} U\left(\zeta \right)=\sqrt{2W_{0}}V_{0}e^{-W_{0}\zeta ^{2}/4}. \end{equation} \tag{ 30 }$$

Therefore, the asymptotic CR pressure depends on the following two parameters: the ISM background diffusivity $W_{0}^{-1}\left(p\right)\equiv W^{-1}\left(p,z=\infty \right)$ (see Equation (10)) and on the CR source pressure Π(p), but only through the parameter V₀ = exp (− q/2). The latter grows linearly with Π ≪ 1 but it saturates when Π ≫ 1:

$$\begin{equation} V_{0}\approx \left\lbrace \begin{array}{@{}ll}\frac{1}{\sqrt{2\pi }}\Pi , & \Pi \ll 1\\ \exp \left(-\frac{1}{2}-4e^{-\Pi /\sqrt{2}}\right), & \Pi \gg 1. \end{array}\right. \end{equation} \tag{ 31 }$$

For practical calculations, it is more convenient to combine both asymptotic regimes into the following simple interpolation formula:

$$\begin{equation} V_{0}=\left[(\sqrt{2\pi }/\Pi )^{5/2}+\left(\sqrt{e}+\frac{3}{2}e^{-\Pi /\sqrt{2}}\right)^{5/2}\right]^{-2/5} \end{equation} \tag{ 32 }$$

which not only recovers both Π ≪ 1 and Π ≫ 1 regimes but also reproduces the transition zone accurately. The latter property is achieved by the choice of numerical parameters of the interpolation, i.e., the powers 5/2, 3/2, and −2/5. The quality of the interpolation may be seen from Figure 4, where the formula in Equation (32) is shown against exact numerical points calculated from Equation (24). The partial pressure $\mathcal {P}\left(z,t\right)$ is then given by

$$\begin{equation} \mathcal {P}=U/\sqrt{t}=\sqrt{2W_{0}/t}V_{0}\left(\Pi \right)e^{-W_{0}z^{2}/4t}. \end{equation} \tag{ 33 }$$

Note that for strong sources (Π ≫ 1) the CR density at large distances becomes independent of the source strength. This is the regime of an efficient ablative escape suppression in which a small (and Π-independent!) number of leaking particles leave behind enough Alfvén fluctuations to limit the leakage of particles remaining in the source to exactly the rate given by Equation (33).

Zoom In Zoom Out Reset image size

The spatial distribution of the CRs becomes even more universal closer to the origin where it falls off as $\mathcal {P}\approx 2/z$ before it turns into an innermost flat-top part of the entire distribution for ζ² < D_NL(p). We have introduced the "self-confinement" diffusion coefficient D_NL as

$$\begin{equation} D_{{\rm NL}}=\frac{2}{V_{0}^{2}}D_{{\rm ISM}}e^{-\Pi }=\frac{2}{V_{0}^{2}w_{{\rm max}}} \end{equation} \tag{ 34 }$$

with $D_{{\rm ISM}}=W_{0}^{-1}$. To obtain this intermediate (D_NL < ζ² < D_ISM) part of the CR spatial distribution, we expand the analytic solution given by Equations (17) and (19) in small V₁ − V ≪ 1. Using the expression for the self-similar CR pressure U from Equation (29) we obtain the following expansion near ζ = ζ₁ ≡ ζ(V₁):

$$\begin{equation} U=\frac{2}{\zeta }\left[1-\frac{3\epsilon }{2}\left(\frac{2}{w_{1}}\right)^{2^{-3/2}}\zeta ^{-\sqrt{2}}-\frac{\epsilon }{2}\left(\frac{w_{1}}{2}\right)^{2^{-3/2}}\zeta ^{\sqrt{2}}\dots \right], \end{equation} \tag{ 35 }$$

where $w_{1}=2V_{1}^{2}/\zeta _{1}^{2}$, and ≡ (q − 1)/2. Note that while the solution has a branching point in auxiliary variable V at V = V₁, it is a regular and single-valued function of the physical variable ζ, throughout the entire half-space ζ > 0, including the branching point of ζ(V) at V = V₁, as it should be. We also notice that while ζ decreases starting from ζ₁, the solution grows slower than 1/ζ, leveling off toward the origin. The above expansion, however, becomes inaccurate for smaller ζ and we alter it below.

Now we turn to this innermost part of the distribution where the particle diffusion coefficient is most strongly decreased. It is convenient to expand the solution in a series in ξ = V/V₀ − 1 ≪ 1. Using the representation of w given by Equation (A2), we obtain

$$w=w_{{\rm max}}\left[\frac{\xi -a+\sqrt{\xi ^{2}-2a\xi +b^{2}}}{b-a}\right]^{-2b},$$

where $a=(1-V_{0}^{2})/(1+V_{0}^{2})$ and $b=V_{0}/\sqrt{\vphantom{A^A}\smash{\hbox{{1+V_{0}^{2}}}}}$. Then, using the expression (A4) for small ζ and Equation (29), we obtain for $U(\zeta)=\mathcal {P}\sqrt{t}$ the following simple result:

$$\begin{equation} U\approx \frac{\sqrt{2w_{{\rm max}}}V_{0}}{1+w_{{\rm max}}\zeta ^{2}/4}\approx \sqrt{2w_{{\rm max}}}V_{0}e^{-w_{{\rm max}}\zeta ^{2}/4}, \end{equation} \tag{ 36 }$$

valid for w_maxζ² ≲ 1. Note that we have rewritten, with the same accuracy, the denominator in the standard "diffusion" form, exp (− w_maxζ²/4), which shows that the CR diffusivity is diminished by a factor W₀/w_max = exp (− Π) ≪ 1 compared to its background level $W_{0}^{-1}$. If, however, w_max ≳ aC_A/κ_B (see Equation (6)), a better approximation would be to simply replace the diffusion coefficient by its Bohm value, as the neglected nonlinear wave interactions (see Section 4) should render δB ∼ B₀. The result in Equation (36) should then be replaced by

$$\begin{equation} U=\sqrt{2aC_{{\rm A}}/\kappa _{{\rm B}}}V_{0}\left(\Pi \right)e^{-aC_{{\rm A}}z^{2}/4\kappa _{{\rm B}}t}. \end{equation} \tag{ 37 }$$

Equations (30), (35), and (36) provide the explicit asymptotic representation of the exact implicit solution given by Equations (17) and (19). This representation is rigorous but somewhat impractical, so we provide, again, an interpolation formula that accurately describes the solution $U\left(\zeta \right)=\mathcal {P}\sqrt{t}$ in the entire range of 0 < ζ < ∞ (Figure 5),

$$\begin{equation} U=2[\zeta ^{5/3}+\left(D_{{\rm NL}}\right)^{5/6}]^{-3/5}e^{-W_{0}\zeta ^{2}/4}, \end{equation} \tag{ 38 }$$

where D_NL(Π) is defined in Equation (34) with V₀ given by the interpolation in Equation (32). The last formula recovers the limiting cases given by Equations (30), (35), and (36) except a slight modification of Equation (30) at large ζ, so that an extra 1/ζ-factor is accepted in the interests of parameterization simplicity. According to Figure 5, however, this deviation from the correct asymptotic expansion of the exponential tail of the distribution is not really noticeable.

Zoom In Zoom Out Reset image size

The overall profile of the partial pressure presented in log–log format (Figure 5) unveils the CR escape as a highly structured process. It comprises a nearly flat-top core (Equation (36)), a nearly 1/ζ pedestal (Equation (35)), and an exponential foot (Equation (30)). We already noted that the escape rate of CR in the foot saturates with the source strength for Π ≫ 1. This is easy to understand as the foot is separated from the core (source) by the flux controlling pedestal, where the CR transport is self-regulated in such a way that a fixed CR flux streams through it to the foot regardless of the strength of the CR source. The solution is shown in Figure 5 as a function of ζ using the self-similar variables $U=\mathcal {P}\sqrt{t}$ and $\zeta =z/\sqrt{t}$. Remarkably, the pedestal portion of the profile does not change with time also in the physical variables, $\mathcal {P},z$: $\mathcal {P}\approx 2/z$. The core, however, sinks in as $\mathcal {P}\propto 1/\sqrt{t}$, but in addition, it expands in z as $\sqrt{D_{{\rm NL}}t}$, to conserve the integrated pressure contained in the core at the expense of the pedestal that shrinks accordingly. The latter disappears completely at t ∼ D_ISM/D_NL, after which the further escape proceeds in a TP regime but with a smaller diffusion coefficient in the core–pedestal region, since waves persist even after most of particles have escaped. At this point, the wave damping should become important (see Section 4).

As the TP approximation is widely used in calculating the CR escape from their sources, we also show, for comparison, one such example in Figure 5. We plot the following simple TP solution

$$\begin{equation} U=\sqrt{\frac{W_{0}}{\pi }}\Pi e^{-W_{0}\zeta ^{2}/4} \end{equation} \tag{ 39 }$$

that can be obtained from Equation (38) or Equation (33) with V₀(Π) taken from Equation (31) or Equation (32) for Π ≪ 1. To compare the self-regulated escape with the TP one, we formally extend the above solution to large Π, as this is normally done within TP approaches. Such an extension clearly underestimates the nonlinear level of the CR pressure by a factor ∼Π⁻¹exp (Π/2) ≫ 1 in the core and in the most of the pedestal. By virtue of lacking self-regulation, it also overestimates the pressure in the exponential foot by a factor Π ≫ 1.

The integrated partial pressure Π is the most important parameter that regulates the CR escape. Therefore, we consider it in some detail, returning to the dimensionful variables and rewriting Π as follows (see Equations (8) and (22)):

$$\begin{equation} \Pi \simeq 3\frac{C_{{\rm A}}}{c}\frac{a\left(p\right)}{r_{{\rm g}}\left(p\right)}\frac{\bar{P}_{{\rm CR}}\left(p\right)}{B_{0}^{2}/8\pi }, \end{equation} \tag{ 40 }$$

where a(p) is the initial size of the CRC, r_g is the CR gyroradius, and $\bar{P}_{{\rm CR}}\left(p\right)$ is their average partial pressure inside the cloud. Although the ratio of Alfvén speed to speed of light, C_A/c, is typically very small (∼10⁻⁵ to 10⁻⁴), the remaining two ratios in the last equation may be fairly large. Indeed, a/r_g should amount to several c/U_ST ≫ 1, with U_ST being the shock speed at the end of the ST stage, as a should exceed the precursor size cr_g/U_ST. Indeed, a significant part of the shock downstream region filled with CR, which have been convected there over the entire history of CR production, may significantly contribute to, if not dominate, the total size of the cloud a (Kang et al. 2009). This is particularly relevant to the CR release during the decreasing maximum energy (reverse acceleration phase; Gabici 2011). In the forward acceleration regime (growing maximum energy) Bohm diffusion, the downstream CR scale height is similar to that of the upstream. A reasonable estimate for $U_{{\rm {\rm ST}}}$ is c/U_ST ≳ 10³. Finally, the CR pressure in the source should considerably exceed the magnetic pressure as both quantities are roughly in equipartition in the background ISM. Alternatively, as is usually assumed, the accelerated electrons are in equipartition with the magnetic energy inside the source, since they should have lost their excessive energy via synchrotron radiation. As electrons are thought to be involved in the acceleration at ∼10⁻² of the proton level, the last ratio in Equation (40) is then ∼10², thus giving, perhaps, an upper bound to the pressure parameter Π ∼ 10²/M_A, where M_A = U_ST/C_A.

One of a few other ways the parameter Π can be looked at is to express it through the average acceleration efficiency $\mathcal {E}=2\bar{P}_{{\rm CR}}/\rho \bar{U}^{2}$, where $\bar{U}$ is the shock velocity, appropriately averaged over the CR acceleration history (one may expect $\bar{U}\gg U_{{\rm {\rm ST}}})$. Then $\Pi \simeq 3(\bar{U}^{2}/cC_{{\rm A}})(a/r_{{\rm g}})\mathcal {E}\simeq 3A(\bar{U}^{2}/U_{{\rm {\rm ST}}}C_{{\rm A}})\mathcal {E}$, where we have introduced a factor A = (a/r_g)U_ST/c ≳ 1. In this form, the estimate indeed boils down to the acceleration efficiency $\mathcal {E}$ with the remaining quantities ($A,\bar{U}$ and U_ST) being more accessible to specific models. Since the acceleration efficiency $\mathcal {E}$ is believed to be at least ≳ 0.1 for productive SNRs, we conclude that the control parameter Π is rather large than small. The caveat here is that it might, in some cases, be too large to limit the applicability of the above treatment. We return to this issue in the Discussion section.

Once the major control parameter is known, we can calculate the distribution function f_CR of the CR, released into the ISM, depending on the distance from the source z and time, using the parameterization in Equation (38). It is convenient to rewrite it in the form of the CR partial pressure $\hat{P}_{{\rm CR}}$ in physical units as follows:

$$\begin{equation} \frac{\hat{P}_{{\rm CR}}}{B_{0}^{2}/8\pi }\simeq 0.8\!\left(\!\frac{M_{{\rm A}}}{A}\frac{10^{4}\,{\rm yr}}{t}\frac{n}{{\rm 1\,cm^{-3}}}\frac{E}{{\rm 1\,GeV}} \! \right)^{1/2}\!\left(\!\frac{1\,\mu {\rm G}}{B} \! \right)^{3/2}U\left(\zeta \right). \end{equation} \tag{ 41 }$$

Here, M_A = U_ST/C_A, n is the plasma number density, E is the particle energy, and B is the magnetic field. The self-similar coordinate ζ can be represented, in turn, in the following way:

$$\begin{equation} \zeta \simeq 17\left(\frac{M_{{\rm A}}}{A}\frac{{\rm 1\,GeV}}{E}\frac{B}{1\,\mu {\rm G}}\frac{10^{4}\,{\rm yr}}{t}\right)^{1/2}\frac{z}{{\rm pc}}. \end{equation} \tag{ 42 }$$

Once this quantity is calculated for given z and t, the CR partial pressure in Equation (41) can be obtained from Equation (38) or from Figure 5. Note that the scale of the CR pressure is determined by the maximum of $U=U\left(0\right)=\sqrt{2W_{0}}V_{0}\exp \left(\Pi /2\right)$, where V₀ is given by Equation (32). The CR partial pressure also depends on the energy through W₀ and Π, apart from the factor $\sqrt{E}$ in Equation (41), which will be discussed below. Finally, the background diffusion coefficient of CR in physical units is related to the dimensionless diffusivity $W_{0}^{-1}$ as follows:

$$D_{{\rm ISM}}=6.6\times 10^{27}\frac{10^{-4}}{W_{0}}\frac{a}{{\rm 1\,pc}}\frac{B}{1\,\mu {\rm G}}\left({\frac{1\,{\rm cm}^{-3}}{n}}\right)^{1/2}\,{\rm cm^{2}}\,{\rm s}^{-1}.$$

As the transition from the flat-top part of the normalized (Equation (8)) partial pressure $\mathcal {P}\left(p\right)$ to the pedestal is described by a function $\mathcal {P}\approx 2\lbrace z^{5/3}+[D_{{\rm NL}}(p)t]^{5/6}\rbrace ^{-3/5}$, the pressure $\mathcal {P}(p)$ is almost p-independent (the corresponding particle momentum distribution $f_{{\rm CR}}\propto \kappa _{{\rm B}}p^{-4}/\hat{z}$, in the physical units, $\hat{z}=az$) for such momenta where D_NL(p) < z²/t, for the fixed dimensionless observation point z and time t. At p = p_br, where D_NL(p_br) = z²/t, the spectrum incurs a break. If we approximate D_NL∝p^δ near p ∼ p_br, the break will have an index δ/2, as may be seen from the above formula for $\mathcal {P}$. The index δ thus derives from the momentum dependence of the D_ISM(p) and from that of the coordinate-integrated CR pressure Π(p) (through the factor exp (− Π) in Equation (34)). Furthermore, if we represent exp (− Π)∝p^−σ and D_ISM∝p^λ at p ∼ p_br, so that δ = λ − σ, then $\mathcal {P}$ is flat at p < p_br for δ > 0 and steepens to p^−δ/2 at p = p_br. Conversely, if δ < 0, $\mathcal {P}$ raises with p as p^−δ/2 at p < p_br and it levels off at p > p_br. Note, however, that Π is momentum independent in the important case of the p⁻⁴ distribution of the initial CRC with the scale height a(p)∝κ_B. This case corresponds to the TP DSA spectrum, ∝p⁻⁴ with a scale height a estimated from the shock precursor size. The break has then an index λ/2 and it is entirely due to the D_ISM momentum dependence.