STOCHASTIC ACCELERATION OF SUPRATHERMAL PARTICLES UNDER PRESSURE BALANCE CONDITIONS

In interplanetary space, there is always a significant population of particles with energies greater than a few times the solar wind thermal energy, where we normally expect very little flux if the particles have a Maxwellian distribution. These suprathermal particles (STPs) form a power-law tail in the energy spectrum. Observed spectra obtained by spacecraft in various regions of the heliosphere often possess a spectral slope corresponding to an isotropic phase space distribution function f∝p⁻⁵, which Fisk & Gloeckler (2006, 2008) called a universal or common spectrum.

The source of these particles is not completely certain (Mewaldt et al. 2005). Some may come from the solar corona, but one most noticeable suprathermal component is the pickup ions. These are originally neutral atoms, mostly interstellar neutrals, penetrating through heliospheric magnetic fields into the inner heliosphere, where they are ionized by solar UV radiation or through charge exchange with the solar wind ions. The newly born pickup ions have almost zero speed in the initial reference frame and they immediately have to convect with the outgoing solar wind plasma with an embedded magnetic field. The initial speed of these pickup ions in the plasma reference frame is then equal to the solar wind speed. Subsequent rapid pitch-angle scattering by magnetic turbulence makes them an isotropic shell distribution with a steep cutoff at exactly the local solar wind speed in the plasma reference frame. Observations in interplanetary space (Gloeckler et al. 1994) often show that pickup ions can extend to much higher speeds. Clearly, they are accelerated in the interplanetary space. Similarly, one can also argue that STPs of solar or other origins should also undergo acceleration during their transport through interplanetary space.

The mechanism for the acceleration of STPs is still under debate. A few authors suggested that interplanetary shocks can accelerate pickup ions and other SPTs (e.g., Giacalone et al. 1997), but most scientists attribute the acceleration of these particles to a stochastic mechanism in the form of diffusion in momentum space by more ubiquitous electromagnetic waves or turbulence in the solar wind plasma (Isenberg, 1987; Schwadron et al. 1996; Le Roux & Ptuskin 1998; Chalov et al. 2004, 2006). Because of the observed near isotropy of the STP distribution function, it is generally accepted that the STPs are scattered in low-frequency magnetohydrodynamic (MHD) turbulence well below the nonrelativistic proton cyclotron frequency. Faraday's induction law then implies that the ratio of the electric and magnetic field strengths of MHD turbulence is given by $|\delta \textbf {\boldsymbol E}|/|\delta \textbf {\boldsymbol B}|\simeq V_A/c\le {\cal O}(10^{-4})$ (with Alfvén velocity V_A), which in interplanetary space is well below unity. In situ electric field fluctuation measurements in the solar wind (e.g., Bale et al. 2005) confirm that $|\delta \textbf {\boldsymbol E}|\ll |\delta \textbf {\boldsymbol B}|$. These small fluctuating electric fields in the wave and turbulence can accelerate or decelerate STPs randomly, depending on the relative direction between wave variation and particle motion.

In the diffusion approximation, one can establish that the particle's momentum distribution obeys the diffusion–convection transport equation, containing both spatial and momentum diffusion and convection terms. Spatial and momentum diffusion can be directly established in the perturbative quasilinear theory (Skilling 1975; Schlickeiser 1989, 2002), but other nonlinear theories can also arrive at diffusive transport (e.g., Qin et al. 2002; Zhang 2010). It has been argued that the quasilinear approximation holds in the interplanetary medium despite the fact that interplanetary turbulence observation indicates that the turbulent magnetic field has the same strength as the ordered uniform (on scales of the particles' gyroradii) background magnetic field. The underlying reason is that most of the magnetic turbulence (about 85% according to Bieber et al. 1996) is in the form of axisymmetric transverse two-dimensional (2D) turbulence that contributes very little to particle scattering (Shalchi & Schlickeiser 2004). However, the transverse 2D turbulence is responsible for efficient wandering of magnetic field lines (see Section 2.1). For the remaining 15% magnitude turbulence, predominantly in the form of slab turbulence, the ratio δB_slab/B₀ is considerably less than unity, fulfilling the quasilinear approximation condition. As an aside we note that the proper diffusion approximation to the Fokker–Planck transport equation (Schlickeiser 1988, 1989; Schlickeiser et al. 1991, 2010), taking into account finite frequency effects of forward and backward moving MHD waves, avoids the earlier claimed 90° scattering problem. In the presence of MHD turbulence (such as Alfvén or magnetosonic waves), the momentum diffusion coefficient is always tied to the particle's pitch-angle diffusion coefficient or parallel spatial diffusion coefficient. Chalov et al. (2004, 2006) considered wave generation of particle acceleration in a self-consistent manner. They found that the amount of acceleration of the pickup ions by Alfvén or magnetosonic turbulence in the interplanetary space is not enough to overcome the adiabatic cooling of particles in the expanding solar wind and produce a significant power tail beyond the initial solar wind speed.

Stochastic acceleration of charged particles can also occur when the particles go through a magnetic medium undergoing random compression or expansion. In this mechanism it is the large-scale convective electric field that accelerates or decelerates particles in a random manner. This case has been investigated by Bykov & Toptygin (1982, 1993) and later Ptuskin (1988), who derive a diffusion equation in the momentum space using a quasilinear method with an additional ensemble average over velocity fluctuations. Ptuskin (1988) studied a case in which the particle diffusion is much too large compared to the characteristic distance between consecutive regions of plasma compression or expansion and he dismissed this mechanism as an important acceleration mechanism for cosmic rays in the interstellar medium.

Recently, observation of the common f∝p⁻⁵ spectrum as drawn renewed attention to the mechanism of particle acceleration by random compressions of large-scale plasma motion. Fisk & Gloeckler (2006, 2008) and Fisk et al. (2010) suggested that this acceleration mechanism can automatically produce the f∝p⁻⁵ spectrum, but the calculation by Jokipii & Lee (2010) disagreed. Zhang (2010) took a different approach to this problem. Instead of the quasilinear approximation, he solved the Parker transport equation for particles moving through a train of compressible waves. Although a special form of the wave trains is prescribed in the model, it is found that acceleration of particles can be very fast if the particle diffusion length is greater than the size of plasma compression region but smaller than the length of the wave train. The behavior of particle acceleration is similar to diffusion in the momentum space, but with a momentum diffusion coefficient D∝p². This mechanism, if we let it go on without significant adiabatic cooling, will soon establish a f∝p⁻³ spectrum. Such a distribution function contains too much energy or pressure. Zhang (2010) and Zhang & Lee (2011) argued that the pressure built up by the accelerated particles can modify the amplitude of the plasma compression waves, such that the rates of particle acceleration by the plasma compressions and adiabatic cooling are equal. At that point, a f∝p⁻⁵-distribution is established, independent of the details of wave or particle parameters as long as the rate of particle acceleration decreases with total particle pressure.

The derivation of the f∝p⁻⁵ spectrum in Zhang (2010) and Zhang & Lee (2011) is done in the reference frame of the particles and the intensity of particles has reached a uniform spatial distribution. This assumption may not be true on the scale size of the heliosphere as the rate of adiabatic cooling of particles and the particle diffusion coefficients could depend on the radial distance from the Sun so that the particle intensity is not uniform. It is the purpose of this paper to explore the spectrum of accelerated pickup ions or STPs in a more realistic condition of particle transport in the heliosphere. We show that the f∝p⁻⁵ spectrum can be established under some plausible conditions of particle transport properties and plasma wave environment, as long as the pressure of accelerated particles has reached a level that is able to force plasma waves to reach the balance between particle acceleration and adiabatic cooling.

Stochastic particle acceleration in cosmic sources occurs in the presence of multiple small-scale turbulent MHD plasma waves propagating in different directions (Zhang & Lee 2011). Its efficiency is characterized by a momentum diffusion coefficient D that is closely related to particle pitch-angle scattering, which causes parallel spatial diffusion along the ordered guide magnetic field, described by the parallel spatial diffusion coefficient κ_∥ (Schlickeiser 1989). The general relation

$$\begin{equation} D^{(T)}(\textbf {\boldsymbol r})=gV_A^2p^2/\kappa ^{(T)}_{\parallel }(\textbf {\boldsymbol r}) \end{equation} \tag{ 1 }$$

holds, where the constant g depends on the nature and geometry of the plasma waves, and has been determined before for the case of parallel propagating shear Alfvén waves (Dung & Schlickeiser 1990) with an additional contribution from isotropically distributed fast magnetosonic waves (Schlickeiser & Miller 1998).

Recently, Zhang & Lee (2011) demonstrated that in the case of random fluctuations of the large-scale plasma flow (Bykov & Toptygin 1993; Ptuskin 1988; Jokipii & Lee 2010), maximum stochastic acceleration occurs with different momentum and spatial diffusion coefficient:

$$\begin{eqnarray} \kappa ^{(C)}(r)\simeq V_cL&=&10^{19}V_{c,7}L_{12}\;\hbox{cm}^{2}\,\hbox{s}^{-1},\nonumber\\ D^{(C)}(\textbf {\boldsymbol r})&=&{V_cp^2\over 15L}={V_c^2p^2\over 15\kappa ^{(C)}(r)}, \end{eqnarray} \tag{ 2 }$$

where L = 10¹²L₁₂ cm is the characteristic length of the plasma velocity variations and V_c = 10⁷V_{c, 7} cm s⁻¹ is the speed of compression waves, which is not necessarily the Alfvén speed. L₁₂ and V_{c, 7} are scaling factors, accounting for deviations from the most likely values L₁₂ = 1 and V_{c, 7} = 1. κ^(C)(r) is the particle spatial diffusion coefficient along the propagation direction of the compression waves.

Ignoring drifts in non-uniform or curved guide magnetic fields, the isotropic part of the distribution function f(r, p, t) of the energetic particles in spherical position coordinates then obeys the radial diffusion–convection transport equation:

$$\begin{equation} {\partial f\over \partial t}-Q(r,p,t)={1\over r^2}{\partial \over \partial r}\left[\kappa (r)r^2{\partial f\over \partial r}\right]-V(r){\partial f\over \partial r} \end{equation}$$

$$\begin{equation} +(\hbox{div}\, \textbf {\boldsymbol V}){p\over 3}{\partial f\over \partial p}+{1\over p^2}{\partial \over \partial p}\left(p^2D{\partial f\over \partial p}\right), \end{equation} \tag{ 3 }$$

where Q(r, p, t) describes sources and sinks of particles,

$$\begin{equation} \hbox{div}\, \textbf {\boldsymbol V}={1\over r^2}{\partial \over \partial r}(r^2V(r)) \end{equation} \tag{ 4 }$$

in radial flow such as the solar wind, and

$$\begin{equation} D(\textbf {\boldsymbol r})=D^{(T)}(\textbf {\boldsymbol r})+D^{(C)}(\textbf {\boldsymbol r}). \end{equation} \tag{ 5 }$$

2.1. Radial Spatial Diffusion

The radial spatial diffusion coefficient in the transport equation (3)

$$\begin{equation} \kappa (r)=\kappa ^{(T)}_{\parallel }(r)\cos ^2\psi (r)+\kappa _\perp (r)\sin ^2\psi (r)+\kappa ^{(C)} \end{equation} \tag{ 6 }$$

is determined by the respective contributions of the parallel, perpendicular, and compressional spatial diffusion coefficients, respectively. Here, ψ(r) denotes the angle of the background Parker spiral magnetic field with respect to the radial direction, which in the outer heliosphere is near π/2. In Equation (6) we assume that the compression waves propagate predominantly along the radial solar wind direction.³ In the outer heliosphere (the likely source of pickup ions) the radial diffusion coefficient is determined predominantly by the perpendicular spatial diffusion coefficient κ_⊥(r) and the compressional contribution κ^(C).

In the absence of drifts in non-uniform or curved guide magnetic fields, due to their small gyroradii the perpendicular spatial diffusion of SPTs is dominated by magnetic field line random walk from the transverse 2D magnetic turbulence. The magnetic diffusion coefficient is given by (Bian et al. 2011; Hauff et al. 2010)

$$\begin{equation} \kappa _m= {l_\perp ^2\over l_z}{K^2\over 1+K^{4/3}} \end{equation} \tag{ 7 }$$

in terms of the Kubo number (Kubo 1963; Vlad et al. 1998; Zimbardo et al. 2000)

$$\begin{equation} K={\delta B_{\perp }\over B_0}{l_z\over l_\perp }, \end{equation} \tag{ 8 }$$

where l_z and l_⊥ are, respectively, the parallel and perpendicular auto correlation lengths of the fluctuating slab and 2D magnetic field components. Expression (7) interpolates between the quasilinear (K ≪ 1) and percolation trapping (K ≫ 1) regions of magnetic field transport. Satellite observations find that the large-scale solar wind fluctuations are nearly isotropic l_z ∼ l_⊥ ∼ 10⁶ km (Tu & Marsch 1995), implying interplanetary Kubo numbers of K = 0.85 for δB_perp/B₀ = 0.85. In the small gyroradius limit, as appropriate for STPs, the particle perpendicular spatial diffusion coefficient is given by (Jokipii 1966; Rechester & Rosenbluth 1978; Bitane et al. 2010)

$$\begin{equation} \kappa _\perp (r)=v\kappa _m=v{l_\perp ^2\over l_z}{K^2\over 1+K^{4/3}}=1.2\times 10^{10}v\;\hbox{cm}^{2}\,\hbox{s}^{-1}. \end{equation} \tag{ 9 }$$

Compared to the compressional contribution κ^(C) given in Equation (2), the field line random walk contribution (9) is negligibly small for STP velocities smaller than v = 10⁴L₁₂V_{c, 7} km s⁻¹, which is well fulfilled for pickup ions adopting the likely values L₁₂ = V_{c, 7} = 1.

2.2. Radial Transport Equation

Moreover, we assume that the momentum diffusion coefficient (5) is dominated by the contribution from compressional waves $D^{(C)}(\textbf {\boldsymbol r})$ given in Equation (2). The STP transport equation (3) then reads

$$\begin{equation} {\partial f\over \partial t}-Q(r,p,t)={1\over r^2}{\partial \over \partial r}\left[\kappa ^{(C)}(r)r^2{\partial f\over \partial r}\right]-V(r){\partial f\over \partial r} \end{equation}$$

$$\begin{equation} +(\hbox{div}\, \textbf {\boldsymbol V}){p\over 3}{\partial f\over \partial p}+{A\over p^2\kappa ^{(C)}(r)}{\partial \over \partial p}\left(p^4{\partial f\over \partial p}\right), \end{equation} \tag{ 10 }$$

with the constant

$$\begin{equation} A={V_c^2\over 15}. \end{equation} \tag{ 11 }$$

As we shall demonstrate in the following, these choices for the radial and momentum diffusion naturally account for the formation of power-law momentum distribution functions of STP as observed. It is straightforward to see that any other choice with momentum-dependent spatial diffusion coefficients will lead to non-power-law STP momentum distributions as solutions of the transport equation (10), which are not observed. The observed momentum power laws of SPT already point to the required microphysical processes causing momentum diffusion and radial spatial diffusion in the outer heliosphere.

2.3. Pressure Balance

For suprathermal charged particles the nonrelativistic treatment of Zhang & Lee (2011) is appropriate. Following their approach we introduce the pressure of energetic particles

$$\begin{equation} P(r,t)={4\pi \over 3m}\int _0^\infty {\it dp}\, p^4f(r,p,t) \end{equation} \tag{ 12 }$$

and consider a uniform (in space and time) source injecting mono-energetic particles, i.e.,

$$\begin{equation} Q(r,p,t)=Q_0\delta (p-p_0). \end{equation} \tag{ 13 }$$

Multiplying Equation (3) by 4πp⁴/3m and integrating over all momenta, we obtain for the total change of pressure (from now on we will drop the index (C) and identify κ = κ^(C))

$$\begin{eqnarray} &&{{\it dp}\over dt}\equiv {\partial P\over \partial t}-{1\over r^2}{\partial \over \partial r}\left[\kappa r^2{\partial P\over \partial r}\right]+V(r){\partial P\over \partial r} ={4\pi Q_0p_0^4\over 3m}\nonumber\\ \ms &&\quad +{4\pi \over 9m}(\hbox{div}\, \textbf {\boldsymbol V}) \int _0^\infty {\it dp}\; p^5{\partial f\over \partial p}+ {4\pi A\over 3m\kappa }\int _0^\infty {\it dp}\; p^2{\partial \over \partial p}\left(p^4{\partial f\over \partial p}\right). \nonumber\\ \end{eqnarray} \tag{ 14 }$$

Integrating the two last terms with respect to p by parts, with the use of f(r, p = 0, t) = f(r, p = ∞, t) = 0, then provides

$$\begin{equation} {{\it dp}\over dt}-{4\pi Q_0p_0^4\over 3m}=\gamma P, \end{equation} \tag{ 15 }$$

with

$$\begin{equation} \gamma ={10A\over \kappa }-{5\over 3}\hbox{div}\, \textbf {\boldsymbol V}={1\over 3}\left[{2V_c^2\over \kappa }-5\hbox{div}\, \textbf {\boldsymbol V}\right]. \end{equation} \tag{ 16 }$$

Equation (15) can be readily integrated over time, yielding with P(t = 0) = 0

$$\begin{equation} P(r,t)={4\pi Q_0p_0^4\over 3m\gamma }[e^{\gamma t}-1]. \end{equation} \tag{ 17 }$$

It was noted by Zhang & Lee (2011) that for γ > 0, in the presence of a large enough initial stochastic acceleration, the pressure in energetic particles (17) becomes infinitely large as t → ∞. In order to avoid such singular behavior they argue that it is necessary to require γ = 0 for a steady-state system, to which we refer as pressure balance condition

$$\begin{equation} {2V^2_c\over \kappa }-5\hbox{div}\, \textbf {\boldsymbol V}=0. \end{equation} \tag{ 18 }$$

This pressure balance condition is achieved by limiting the particle scattering and acceleration rates through a nonlinear feedback mechanism of particle pressure on plasma waves.

In spherical flows with

$$\begin{equation} V(r)=V_0r^{\alpha } \end{equation} \tag{ 19 }$$

with V₀ = constant and α ∈ [0, 1), this provides a radially increasing spatial diffusion coefficient of the form

$$\begin{equation} \kappa (r)={2V_c^2r^{1-\alpha }\over 5(2+\alpha)V_0}. \end{equation} \tag{ 20 }$$

Such a radial variation of the spatial diffusion coefficient is quite reasonable in the heliospheric magnetic field as the inferred radial dependences κ(r)∝r^ξ with ξ ∈ [1.1, 1.4] from the intensities of galactic and anomalous cosmic rays (Fujii & McDonald 2005) indicate.

2.4. Particle Transport Under Pressure Balance Condition

The radial particle transport equation (10) under the pressure balance condition (20) then reads

$$\begin{equation} {\partial f\over \partial t}-Q(r,p,t)={2V_c^2\over 5(2+\alpha)V_0r^2}{\partial \over \partial r}\left[r^{3-\alpha }{\partial f\over \partial r}\right]-V_0r^{\alpha }{\partial f\over \partial r} \end{equation}$$

$$\begin{equation} +{(2+\alpha)V_0p\over 3r^{1-\alpha }}{\partial f\over \partial p}+{(2+\alpha)V_0\over 6r^{1-\alpha }p^2}{\partial \over \partial p}\left(p^4{\partial f\over \partial p}\right). \end{equation} \tag{ 21 }$$

Zhang & Lee (2011) have solved a spatially non-diffusive transport equation for the particles for the case of a constant solar wind velocity (α = 0), demonstrating that due to the pressure balance coupling of momentum diffusion and adiabatic deceleration, the particle distribution at large momenta settles with a distribution f(p)∝p⁻⁵, which is commonly seen in STP (Fisk & Gloeckler 2006, 2008). While the spatial diffusion can be neglected for very low energy particles in the region where its density variation is small within the particle diffusion length, here we extend our study to include stationary solutions of the full spatially dependent transport equation (21) to investigate the influence of the spatially increasing spatial diffusion coefficient (20) on the momentum spectrum of accelerated particles under pressure balance conditions. Because the spatial diffusion coefficient increases with radial distance and particle energy, its effect cannot be ignored at large enough radii r and momentum p.

For separable source distributions Q(r, p) = q₁(r)q₂(p) the stationary transport equation under pressure balance condition is

$$\begin{equation} {\cal L}_rf(r,p)+{\cal L}_pf(r,p)=-{6r^{1-\alpha }q_1(r)\over (2+\alpha)V_0}q_2(p) \end{equation} \tag{ 22 }$$

with the momentum-independent spatial operator

$$\begin{equation} {\cal L}_r={\beta \over r^{1+\alpha }}{d\over dr}r^{3-\alpha }{d\over dr}- {6r\over 2+\alpha }{d\over dr} \end{equation} \tag{ 23 }$$

with the constant

$$\begin{equation} \beta =\left({6\over 2+\alpha }\right)^2{A\over V_0^2} \end{equation} \tag{ 24 }$$

and the position-independent momentum operator

$$\begin{equation} {\cal L}_p=2p{d\over {\it dp}}+{1\over p^2}{d\over {\it dp}}p^4{d\over {\it dp}}. \end{equation} \tag{ 25 }$$

The operators ${\cal L}_r$ and ${\cal L}_p$ can be understood as the rate of particle transport in r and p, respectively, normalized by half the rate of adiabatic cooling $\hbox{div }\textbf {\boldsymbol V}/6=(2+\alpha)V_0/6r^{1-\alpha }$.

As momentum boundary conditions we use finiteness at p → 0 and vanishing particle density as p → ∞, i.e.,

$$\begin{equation} f(r,p=0)=\hbox{finite},\;\;\; f(r,p\rightarrow \infty)=0. \end{equation} \tag{ 26 }$$

As spatial boundary condition we use either free-escape at the inner (r₀) and outer (R) radii, i.e.,

$$\begin{equation} f(r=r_0,p)= f(r=R,p)=0, \end{equation} \tag{ 27 }$$

or reflecting boundary conditions at r₀ and R, i.e.,

$$\begin{equation} {\partial f\over \partial r}\vert _{r=r_0}={\partial f\over \partial r}\vert _{r=R}=0, \end{equation} \tag{ 28 }$$

respectively.

According to the scattering time method (Wang & Schlickeiser 1987; Schlickeiser 2002, Ch. 14), Equation (22) is solved by the convolution

$$\begin{equation} f(p,r)=\int _0^\infty du\, M(r,u)H(p,u), \end{equation} \tag{ 29 }$$

where H(p,u) satisfies the momentum boundary conditions (26) and

$$\begin{equation} {\partial H\over \partial u}={\cal L}_pH \end{equation} \tag{ 30 }$$

with

$$\begin{equation} H(p,u=\infty)=0,\;\; H(p,u=0)=q_2(p). \end{equation} \tag{ 31 }$$

The age distribution M(r, u) satisfies the spatial boundary conditions (27) and

$$\begin{equation} {\partial M\over \partial u}={\cal L}_rM, \end{equation} \tag{ 32 }$$

with

$$\begin{equation} M(r,u=\infty)=0,\;\; M(r,u=0)={6r^{1-\alpha }q_1(r)\over (2+\alpha)V_0}. \end{equation} \tag{ 33 }$$

3.1. Leaky-box Equations

The spatial operator (23) is of Sturm–Liouville form and thus has the complete eigenfunction system M_i(r). The age distribution M(r, u) then can be expanded in this orthonormal system as

$$\begin{equation} M(r,u)=\sum _ic_iM_i(r)e^{-\lambda _iu}, \end{equation} \tag{ 34 }$$

where λ_i denotes the real eigenvalues of the spatial operator, being determined by the spatial boundary conditions (27), and the coefficients c_i have to be calculated from the initial condition (33) using orthonormality relations of M_i:

$$\begin{equation} M(r,u=0)={6r^{1-\alpha }q_1(r)\over (2+\alpha)V_0}=\sum _ic_iM_i(r). \end{equation} \tag{ 35 }$$

Inserting the expansion (34) into the convolution (29) leads to

$$\begin{equation} f(p,r)=\sum _ic_iM_i(r)R_i(p), \end{equation} \tag{ 36 }$$

where each of the functions (see Schlickeiser 2002, Ch. 14.3.2)

$$\begin{equation} R_i(p)=\int _0^\infty duH(p,u)e^{-\lambda _iu} \end{equation} \tag{ 37 }$$

obeys the leaky-box equations:

$$\begin{equation} {\cal L}_pR_i(p)-\lambda _iR_i(p)=-q_2(p). \end{equation} \tag{ 38 }$$

Equation (36) indicates that the steady-state solution M(r, u) can be expressed as a (infinite) sum of products of the functions R_i(p) and c_iM_i(r). Each function R_i(p) satisfies a leaky-box equation (38), where the corresponding spatial eigenvalues λ_i enter in the form of a catastrophic escape loss rate.

3.2. Solution of the Leaky-box Equations

Inserting the momentum operator (25) then provides for the leaky-box equation (38):

$$\begin{equation} 2p{{\it dR}_i\over {\it dp}}+{1\over p^2}{d\over {\it dp}}p^4{dR_i\over {\it dp}}-\lambda _iR_i(p)=-q_2(p). \end{equation} \tag{ 39 }$$

The self-adjoint form of this equation reads

$$\begin{equation} {d\over {\it dp}}\left(p^6{dR_i\over {\it dp}}\right)-\lambda _ip^4R_i(p)=-p^4q_2(p), \end{equation} \tag{ 40 }$$

which for arbitrary injection spectra q₂(p) is solved by

$$\begin{equation} R_i(p)=\int _0^\infty {\it dp}_0\, p_0^4q_2(p_0)G(p,p_0). \end{equation} \tag{ 41 }$$

The Green's function G(p, p₀) obeys

$$\begin{equation} {d\over {\it dp}}\left(p^6{dG\over {\it dp}}\right)-\lambda _ip^4G=-\delta (p-p_0), \end{equation} \tag{ 42 }$$

and with the momentum boundary condition (26) is given by

$$\begin{equation} G(p,p_0)={(pp_0)^{-5/2}\over 2\mu _i} \left\lbrace \begin{array}{@{}ll}(p/p_0)^{\mu _i} & {\rm for}\; p\le p_0 \\ (p/p_0)^{-\mu _i} & {\rm for}\; p\ge p_0, \end{array}\right. \end{equation} \tag{ 43 }$$

with

$$\begin{equation} \mu _i=\sqrt{{25\over 4}+\lambda _i}\ge {5\over 2}. \end{equation} \tag{ 44 }$$

Inserting the Green's function (43) into Equation (41) we find

$$\begin{eqnarray} R_i(p)&=&{1\over 2\mu _ip^{5/2}}\Bigl [p^{-\mu _i}\int _0^p{\it dp}_0\, p_0^{{3\over 2}+\mu _i}q_2(p_0)\nonumber\\ \ms &&+p^{\mu _i}\int _p^\infty {\it dp}_0\, p_0^{{3\over 2}-\mu _i}q_2(p_0)\Bigr ]. \end{eqnarray} \tag{ 45 }$$

In the case of a mono-energetic injection q₂(p₀) = Q₀δ(p − p_I), we readily find the broken power-law solution:

$$\begin{equation} R_i(p,p_I)={Q_0\over 2\mu _ip_I} \left\lbrace \begin{array}{@{}ll} (p/p_I)^{\mu _i-(5/2)} & {\rm for}\; p\le p_I \\ (p/p_I)^{-\mu _i-(5/2)} & {\rm for}\; p\ge p_I, \end{array}\right. \end{equation} \tag{ 46 }$$

so that the spectral index of accelerated particles $R_i(p\ge p_I)\propto p^{-s_i}$ with momentum p ⩾ p_i is given by

$$\begin{equation} s_i=\mu _i+{5\over 2}={5\over 2}\left[1+\sqrt{1+{4\lambda _i\over 25}}\right], \end{equation} \tag{ 47 }$$

which approaches the value s_i = 5 from above, provided that the spatial eigenvalues or escape loss rates λ_i ≪ (25/4) = 6.25 are much smaller than the critical value 6.25. If such small eigenvalues are possible, a suitable explanation of the common p⁻⁵-momentum spectrum has been established.

3.3. Interlude

Equations (36), (45), and (46) represent the full formal solution of the stationary particle transport equation for arbitrary separable particle source distributions and general radial solar wind variations of the form (19). The resulting solution is the superposition of an infinite number of momentum power-law distributions, whose spectral indices are determined by the eigenvalues λ_n of the spatial propagation equation. Each power-law contribution is weighted by the respective expansion coefficient c_n, depending on the spatial source distribution of STPs. The power law resulting from the smallest spatial eigenvalue λ₁ will dominate the solution at sufficiently high-momentum values because the spectral indices from the contribution of larger spatial eigenvalues will be greater than s₁.

The determination of the eigenvalues λ_n and the expansion coefficients c_n require solving the eigenequation

$$\begin{equation} {\cal L}_rM_n(r)+\lambda _nM_n(r)=0, \end{equation} \tag{ 48 }$$

subject to the spatial boundary conditions (27) and/or (28).

Throughout this work we restrict our analysis to the case of a constant solar wind velocity (α = 0), whereas radially varying solar wind velocities will be investigated in future work. As we will demonstrate, the case of a constant solar wind velocity allows us the full analytic calculation of the eigenvalues and the expansion coefficients.

In particular, we will investigate the influence of the chosen spatial boundary conditions on the eigenvalues and thus on the power-law spectral indices of the accelerated STPs. Of particular interest is the calculation of the smallest possible eigenvalue λ₁ for different spatial boundary conditions, which is considered in Section 5. If this is much smaller than the value 6.25, the high-momentum spectrum of accelerated particles equals the p⁻⁵-spectrum.

In the following we investigate the solution of the eigenequation (48) for the constant solar wind velocity case under four different spatial boundary conditions:

(1) free-escape boundary conditions M_n(r₀) = M_n(R) = 0 at the inner (r₀) and (R) and outer radii;

(2) reflecting boundary conditions $(dM_n/dr)_{r_0}= (dM_n/dr)_R=0$ at the inner and outer radii;

(3) inner free-escape (M_n(r₀) = 0) and outer reflecting ((dM_n/dr)_R = 0) boundary conditions; and

(4) inner reflecting ($(dM_n/dr)_{r_0}=0$) and outer free-escape (M_n(R) = 0) boundary conditions.

The boundary conditions in case D closely represent the real boundary conditions of the heliosphere. At the inner boundary, the magnetic field is strong, so the majority of particles coming in from the outer heliosphere are mirrored back, suggesting that a reflective boundary would be most appropriate. At large radial distances, the magnetic field is too weak to contain these particles and the radial geometry also dilutes the particle density, making an escape boundary most likely.

For constant solar wind speed (α = 0) the steady-state transport equation (21)

$$\begin{eqnarray} -Q(r,p)&=&{V_c^2\over 5V_0r^2}{\partial \over \partial r}\left[r^3{\partial f\over \partial r}\right]-V_0{\partial f\over \partial r}\nonumber\\ \ms &&+{2V_0p\over 3r}{\partial f\over \partial p}+{V_0\over 3rp^2}{\partial \over \partial p}\left(p^4{\partial f\over \partial p}\right) \end{eqnarray} \tag{ 49 }$$

is controlled by the four closely related timescales for

(1) spatial convection

$$\begin{equation} \tau _C={R\over V_0}, \end{equation} \tag{ 50 }$$

(2) adiabatic deceleration

$$\begin{equation} \tau _A={3R\over 2V_0}={3\over 2}\tau _C, \end{equation} \tag{ 51 }$$

(3) stochastic acceleration

$$\begin{equation} \tau _S={3R\over V_0}=3\tau _C, \end{equation} \tag{ 52 }$$

and (4) spatial diffusion

$$\begin{equation} \tau _D=5\left({V_0\over V_c}\right)^2{R\over V_0}={3\over \beta _0}\tau _C, \end{equation} \tag{ 53 }$$

with the constant

$$\begin{equation} \beta _0={9A\over V_0^2}={3V_c^2\over 5V_0^2}={0.6\over M_A^2}, \end{equation} \tag{ 54 }$$

where M_A = V₀/V_c is the Mach number of the solar wind speed. The parameter β₀, determined by the Mach number of the solar wind speed, therefore has the strongest influence on the resulting particle spectrum f(p, r). τ_D is longer than τ_C in the supersonic wind, but these timescales can be comparable in the subsonic wind. This makes the spatial diffusion important in the outer heliosphere, particularly in the heliosheath.

For the case of a constant solar wind speed α = 0, the radial spatial diffusion coefficient

$$\begin{equation} \kappa (r)={\lambda (r)v\over 3}={V_c^2r\over 5V_0} \end{equation} \tag{ 55 }$$

implies the radially increasing scattering mean free path

$$\begin{equation} \lambda (r)={1\over 5M_A^2}{V_A\over v}r, \end{equation} \tag{ 56 }$$

which for an adopted Mach number M_A = 2 yields the small ratio

$$\begin{equation} {\lambda (r)\over R}={V_c\over 20v}{r\over R}, \end{equation} \tag{ 57 }$$

so that for all STPs with velocities v ≫ 0.05 V_c the spatial diffusion approximation is justified.

Moreover, for a constant solar wind speed α = 0, the spatial operator (23) reads

$$\begin{equation} {\cal L}_r={\beta _0\over r}{d\over dr}r^3{d\over dr}-3r{d\over dr}. \end{equation} \tag{ 58 }$$

The initial condition (33) becomes

$$\begin{equation} M(r,u=0)={3rq_1(r)\over V_0}. \end{equation} \tag{ 59 }$$

With the ansatz (34) for the age distribution the eigenfunctions obey the eigenequation:

$$\begin{eqnarray} &&{\cal L}_rM_i(r)+\lambda _iM_i(r)\nonumber\\ &&\quad=\beta _0r^2{d^2M_i(r)\over dr^2}+3(\beta _0-1)r{dM_i(r)\over dr}+\lambda _iM_i(r)=0, \nonumber\\ \end{eqnarray} \tag{ 60 }$$

or

$$\begin{equation} r^2{d^2M_i\over dr}-2\eta r{dM_i\over dr}+\Lambda _iM_i=0 \end{equation} \tag{ 61 }$$

with

$$\begin{equation} \eta ={3\over 2}\left({1\over \beta _0}-1\right)={5\over 2}\left[M_A^2-{3\over 5}\right]>0 \end{equation} \tag{ 62 }$$

and

$$\begin{equation} \Lambda _i=\lambda _i/\beta _0. \end{equation} \tag{ 63 }$$

4.1. Self-adjoint Form

The self-adjoint form of the eigenequation (61) reads

$$\begin{equation} {d\over dr}\left(r^{-2\eta }{dM_i\over dr}\right)+\Lambda _ir^{-2(\eta +1)}M_i=0, \end{equation} \tag{ 64 }$$

providing with the boundary conditions (27)–(28) the orthonormality relations

$$\begin{equation} \int _{r_0}^Rdr\, r^{-2(\eta +1)}M_n(r)M_m(r)=j_n\delta _{m,n} \end{equation} \tag{ 65 }$$

with

$$\begin{equation} j_n=\int _{r_0}^Rdr\, r^{-2(\eta +1)}M_n^2(r). \end{equation} \tag{ 66 }$$

4.2. Normal Form

With M_i(r) = r^ηE(r) the normal form of Equation (61) reads

$$\begin{eqnarray} &&r^2{d^2E\over dr^2}+\left[\Lambda _i-\eta (\eta +1)\right]E\nonumber\\ &&\quad =r^2{d^2E\over dr^2}+{1\over 4}E+\left[\Lambda _i-\left(\eta +{1\over 2}\right)^2\right]E=0. \end{eqnarray} \tag{ 67 }$$

Equation (67) is readily solved by E(r)∝r^k with

$$\begin{equation} \left(k-{1\over 2}\right)^2=\left(\eta +{1\over 2}\right)^2-\Lambda _i. \end{equation} \tag{ 68 }$$

Obviously, for each specified boundary conditions (A)–(D), we have to consider the two cases of (1) small eigenvalues Λ_i < (η + (1/2))², corresponding to λ_i < β₀[η + (1/2)]², and (2) large eigenvalues Λ_i < (η + (1/2))², corresponding to λ_i ⩾ β₀[η + (1/2)]². The first case is particularly important because it implies the p⁻⁵-spectrum at high momenta if the smallest eigenvalue λ₁ is much smaller than 6.25. We consider each case in turn in the next section and in the Appendix, respectively.

In this case we set

$$\begin{equation} \psi ^2= (\eta +(1/2))^2-\Lambda _i>0. \end{equation} \tag{ 69 }$$

Equation (68) then provides

$$\begin{equation} E(r)=r^{1/2}[a_1r^{\psi }+a_2r^{-\psi }]. \end{equation} \tag{ 70 }$$

5.1. Cases A and D

The free-escape boundary condition at r = R gives

$$\begin{equation} a_2=-a_1R^{2\psi }, \end{equation} \tag{ 71 }$$

so that

$$\begin{equation} E(r)=a_1r^{1/2}\sinh \left[\psi \ln \left({R\over r}\right)\right]. \end{equation} \tag{ 72 }$$

5.1.1. Case A

The free-escape condition at the inner boundary E(r₀) = 0 in this case then demands

$$\begin{equation} \sinh [\psi \ln (R/r_0)]=0, \end{equation} \tag{ 73 }$$

providing as the only possibility ψ = 0, which implies according to Equation (72) the unusable trivial solution E_i(r) = 0. Therefore, in this case of the boundary conditions no small eigenvalue Λ_i exists.

5.1.2. Case D

The solution (72) implies

$$\begin{equation} M(r)=a_1M_1(r)=a_1r^{\eta +(1/2)}\sinh [\psi \ln (R/r)], \end{equation} \tag{ 74 }$$

yielding

$$\begin{eqnarray} {dM\over dr}&=&a_1r^{\eta -(1/2)}[(\eta +(1/2))\sinh [\psi \ln (R/r)]\nonumber\\ &&-\psi \cosh [\psi \ln (R/r)]]. \end{eqnarray} \tag{ 75 }$$

The inner reflecting boundary conditions $(dM/dr)_{r_0}=0$ in this case then provide the condition

$$\begin{equation} \tanh [\psi \ln (R/r_0)]={2\psi \over 1+2\eta } \end{equation} \tag{ 76 }$$

for the smallest first eigenvalue Λ₁ in this case. For ψ > η + (1/2), no solutions of Equation (76) exist, because the tanh function is limited to values less than unity. Consequently, the value of ψ has to be smaller than η + (1/2), which according to Equation (69) implies that the eigenvalue Λ₁ is nonnegative.

Moreover, for a nontrivial solution of the transcendental equation (76) the slope ln (R/r₀) of the tanh  function at ψ = 0 has to be larger than the slope (1/(1/2) + η) of the line (2ψ/1 + 2η), i.e.,

$$\begin{equation} \ln (R/r_0)>{2\over 1+2\eta }={2\over 5M_A^2-2}. \end{equation} \tag{ 77 }$$

We therefore find for this boundary condition case one small eigenvalue Λ₁, provided that Equation (77) is fulfilled.

With y₀ = ln (R/r₀) and K = ψy₀ we can write Equation (76) as

$$\begin{equation} {\tanh K\over K}={2\over (1+2\eta)y_0}. \end{equation} \tag{ 78 }$$

For values of

$$\begin{equation} y_0=\ln (R/r_0)\gg {2\over 5M_A^2-2} \end{equation} \tag{ 79 }$$

the solution of Equation (78) can be approximated as

$$\begin{equation} K=(\eta +(1/2))y_0-\epsilon ,\quad \epsilon \ll 1 \end{equation} \tag{ 80 }$$

with

$$\begin{eqnarray} && 1-{\epsilon \over (\eta +(1/2))y_0}=\tanh [(\eta +(1/2))y_0-\epsilon ]\nonumber\\ &&\quad \simeq 1-2e^{-(1+2\eta)y_0-2\epsilon }\simeq 1-2e^{-(1+2\eta)y_0}, \end{eqnarray} \tag{ 81 }$$

yielding

$$\begin{equation} \epsilon \simeq (1+2\eta)y_0e^{-(1+2\eta)y_0}. \end{equation} \tag{ 82 }$$

Equation (80) then becomes

$$\begin{equation} K\simeq (\eta +(1/2))y_0[1-2e^{-(1+2\eta)y_0}]. \end{equation} \tag{ 83 }$$

With the definition of K = ψy₀ we obtain as the only small solution of the transcendental equation (76)

$$\begin{equation} \psi _1\simeq (\eta +(1/2))[1-2e^{-(1+2\eta)y_0}]. \end{equation} \tag{ 84 }$$

Combination with Equation (69) then provides for the small eigenvalue

$$\begin{equation} \Lambda _1\simeq (1+2\eta)^2e^{-(1+2\eta)y_0}=(1+2\eta)^2\left({R\over r_0}\right)^{-(1+2\eta)}. \end{equation} \tag{ 85 }$$

Equations (63) and (54) then yield

$$\begin{eqnarray} \lambda _1&=&\beta _0\Lambda _1=3{\big[5M_A^2-2\big]^2\over 5M_A^2}\left({r_0\over R}\right)^{5M_A^2-2}\nonumber\\ \ms &=& 15M_A^2\left(1-{2\over 5M_A^2}\right)^2e^{-(5M_A^2-2)\ln (R/r_0)}, \end{eqnarray} \tag{ 86 }$$

which for large R ≫ r₀ and Mach numbers M_A > 1 is small compared to unity.

The expansion coefficient a₁ can be determined by using the orthonormality property of M₁(r), i.e.,

$$\begin{eqnarray} j_1&=&\int _{r_0}^Rdr\, r^{-2(\eta +1)}M_1^2(r)= \int _{r_0}^Rdr\, r^{-1}\sinh ^2[\psi _1\ln (R/r)\nonumber\\ \ms &=&{1\over \psi _1}\int _0^{\psi _1\ln (R/r_0)}dx\, \sinh ^2x\nonumber\\ \ms &=&{\sinh [2\psi _1\ln (R/r_0)]\over 4\psi _1}-{1\over 2}\ln (R/r_0). \end{eqnarray} \tag{ 87 }$$

5.2. Cases B and C

The solution (70) implies

$$\begin{equation} M(r)=r^{\eta +(1/2)}[a_1r^{\psi }+a_2r^{-\psi }], \end{equation} \tag{ 88 }$$

so that the outer reflecting boundary condition (dM/dr)_R = 0 in cases B and C yields

$$\begin{equation} {dM\over dr}=a_3r^{\eta -(1/2)}\sinh \left[\psi \ln (R/r)\right]. \end{equation} \tag{ 89 }$$

5.2.1. Case B

The inner reflecting boundary condition $(dM/dr)_{r_0}=0$ in this case then can only be fulfilled for ψ = 0, so that no small eigenvalues with nontrivial solutions exist in this case either.

5.2.2. Case C

The outer reflecting boundary condition (dM/dR)_R = 0 implies for the solution (88)

$$\begin{eqnarray} M(r)&=&a_4r^{\eta +(1/2)}[(\eta +(1/2)+\psi)(R/r)^{\psi }\nonumber\\ \ms &&-(\eta +(1/2)-\psi)(R/r)^{-\psi }] =2a_4r^{\eta +(1/2)} \nonumber\\ \ms &&\times \left[(\eta +(1/2))\sinh [\psi \ln (R/r)]+\psi \cosh [\psi \ln (R/r)]\right].\nonumber\\ \end{eqnarray} \tag{ 90 }$$

The inner free-escape condition M(r₀) = 0 then provides the condition

$$\begin{equation} \tanh [\psi \ln (R/r_0)]=-{\psi \over \eta +(1/2)} \end{equation} \tag{ 91 }$$

with the only solution ψ = 0. Also in this case of boundary conditions no small eigenvalue with nontrivial eigensolution exists.

5.3. Interlude

As the most important result of this section we found that the spatial boundary condition D (reflecting inner and free-escape outer boundary) is the only case where a small eigenvalue is possible, provided that

$$\begin{equation} \ln (R/r_0)>{2\over 1+2\eta }={2\over 5M_A^2-2}. \end{equation} \tag{ 92 }$$

The eigenvalue is then given by

$$\begin{equation} \lambda _1=15M_A^2\left(1-{2\over 5M_A^2}\right)^2e^{-(5M_A^2-2)\ln (R/r_0)}, \end{equation} \tag{ 93 }$$

which for large R ≫ r₀ and Mach numbers M_A > 1 is small compared to unity and to the value 6.25. This mode gives a negligible particle escape loss.

For completeness we consider in Appendix the case of large eigenvalues Λ_i ⩾ (η + (1/2))².

6.1. Collection of Earlier Results

According to Equations (36) and (46) the full stationary solution for α = 0 and mono-energetic injection with rate Q₀ at momentum p_I reads

$$\begin{eqnarray} &&f(p,r)={Q_0\over 2p_I}\sum _{i=1}^\infty {c_i\over \mu _i}M_i(r)\nonumber\\ &&\quad \times \left[\left({p\over p_I}\right)^{\mu _i-{5\over 2}}H[p_I-p]+\left({p\over p_I}\right)^{-(\mu _i+{5\over 2})}H[p-p_I]\right] \nonumber\\ \end{eqnarray} \tag{ 94 }$$

with the step function H(x) = (1 + (x/|x|))/2. For arbitrary spatial source distributions q₁(r), the expansion coefficients are given by

$$\begin{equation} c_i={3\over V_0j_i}\int _{r_0}^Rdr\, r^{-2\eta -1}q_1(r)M_i(r) \end{equation} \tag{ 95 }$$

with

$$\begin{equation} j_i=\int _{r_0}^Rdr\, r^{-2(\eta +1)}M_i^2(r). \end{equation} \tag{ 96 }$$

For an inner reflecting (at r₀) and outer (at R) free-escape boundary condition, we found

$$\begin{equation} j_1={\sinh [2\psi _1\ln (R/r_0)]\over 4\psi _1}-{1\over 2}\ln (R/r_0), \end{equation} \tag{ 97 }$$

$$\begin{equation} M_1(r)=r^{\eta +{1\over 2}}\sinh \left[\psi _1\ln (R/r)\right] \end{equation} \tag{ 98 }$$

with

$$\begin{equation} \eta ={5\over 2}\big(M_A^2-0.6\big), \end{equation} \tag{ 99 }$$

and

$$\begin{equation} \lambda _1={3\over 5M_A^2}\Lambda _1={3\over 5M_A^2}(1+2\eta)^2\left({r_0\over R}\right)^{1+2\eta }. \end{equation} \tag{ 100 }$$

Likewise, for i ⩾ 2

$$\begin{equation} M_i(r)=r^{\eta +{1\over 2}}\sin \left[\nu _i\ln (R/r)\right] \end{equation} \tag{ 101 }$$

with

$$\begin{equation} \nu _i=(i-1)\pi \left[1+{1\over (\eta +(1/2))\ln (R/r_0)}\right], \end{equation} \tag{ 102 }$$

so that

$$\begin{equation} j_i={1\over 2}\ln (R/r_0)-{\sin \big(2\nu _i\ln {R\over r_0}\big)\over 4\nu _i}, \end{equation} \tag{ 103 }$$

$$\begin{eqnarray} &&\fontsize{9}{11}\hbox{$\dsty \lambda _i={3\over 5M_A^2}\left[\!\left(\eta +{1\over 2}\right)^2+(i-1)^2\pi ^2\left[\!1+{1\over (\eta +(1/2))\ln (R/r_0)}\!\right]^2\!\right]$}\nonumber\\ \end{eqnarray} \tag{ 104 }$$

and

$$\begin{eqnarray} &&\fontsize{8.5}{11}\hbox{$\mu _i=\left[\dsty {25\over 4}+{3\over 4(3+2\eta)}\right((1+2\eta)^2+4(i-1)^2\pi ^2+{16(i-1)^2\pi ^2\over (1+2\eta)\ln {R\over r_0}}\right)\right]^{1/2}.$}\nonumber\\ \end{eqnarray} \tag{ 105 }$$

6.2. Illustrative Choice of Parameters

In the slow solar wind we adopt the Mach number M_A = 1.35, implying the parameter value η = 3.056. This small Mach number should be regarded as a best-fit value and it corresponds to rather strong velocity fluctuations. For higher values of M_A we obtain much too small values of the first expansion coefficient a₁ compared to the others a_i, i ⩾ 2, so that the momentum spectrum of accelerated STPs approaches the p⁻⁵-power-law spectrum (shown in Figures 1 and 3) at particle momenta that are too high. Condition (79) for the existence of a small eigenvalue then demands

$$\begin{equation} \ln (R/r_0)\gg 0.28, \end{equation} \tag{ 106 }$$

which is well fulfilled for our adopted values R = 10r₀.

Zoom In Zoom Out Reset image size

In Table 1 we calculate the first five spatial eigenvalues λ_i and the associated power-law spectral indices μ_i(λ_i). The first extremely small eigenvalue Λ₁ implies μ₁ = 2.5, whereas the higher eigenvalues yield steeper spectral index values μ_i > 3.8.

6.3. Expansion Coefficients

Guided by more exact determinations (Zhang & Schlickeiser 2012) we approximate the spatial radial distribution of injected particles as a uniform outer ring distribution:

$$\begin{equation} q_1(r)=H[r-r_1]H[r_2-r] \end{equation} \tag{ 107 }$$

with r₁ = 0.5R and r₂ = 0.9R. It is convenient to scale these expansion coefficients as

$$\begin{equation} c_i={3\over V_0R^{\eta -{1\over 2}}}a_i \end{equation} \tag{ 108 }$$

with the dimensionless values

$$\begin{equation} a_i={R^{\eta -{1\over 2}}\over j_i}\int _{r_0}^Rdr\, r^{-2\eta -1}q_1(r)M_i(r). \end{equation} \tag{ 109 }$$

With Equation (94) we then obtain for the full stationary solution for α = 0 and mono-energetic injection with rate Q₀ at momentum p_I

$$\begin{equation} f(p,r)={3Q_0R\over 2p_IV_0}F(p,r) \end{equation} \tag{ 110 }$$

with

$$\begin{eqnarray} F(p,r)&=&{1\over R^{\eta +{1\over 2}}}\sum _{i=1}^\infty {a_i\over \mu _i}M_i(r)\Bigl [\left({p\over p_I}\right)^{\mu _i-{5\over 2}}H[p_I-p]\nonumber\\ &&+\left({p\over p_I}\right)^{-(\mu _i+{5\over 2})}H[p-p_I])\Bigr ]. \end{eqnarray} \tag{ 111 }$$

With the spatial source distribution (107) we then obtain for the first expansion coefficient:

$${\fontsize{8.2}{11}{\begin{eqnarray} a_1&=&{R^{\eta -{1\over 2}}\over j_1}\int _{0.5R}^{0.9R}dr\, r^{-2\eta -1}M_1(r)\nonumber\\ &=&{R^{\eta -{1\over 2}}\over j_1}\int _{0.5R}^{0.9R}dr\, r^{-\eta -{1\over 2}}\sinh \left(\psi _1\ln {R\over r}\right)\nonumber\\ &=&{1\over j_1\psi _1}\int _{\psi _1\ln {10\over 9}}^{\psi _1\ln 2}dx\, e^{(2\eta -1)x \over 2\psi _1}\sinh x={1\over j_1\big[\psi ^2_1-(\eta -{1\over 2})^2\big]}\nonumber\\ &&\times \left [2^{\eta -{1\over 2}}\left[\psi _1\cosh (\psi _1\ln 2)-\left(\eta -{1\over 2}\right)\sinh (\psi _1\ln 2)\right]\nonumber\\ &&-\left({10\over 9}\right)^{\eta -{1\over 2}}\left[\psi _1\cosh \left(\psi _1\ln {10\over 9}\right)-\left(\eta -{1\over 2}\right)\sinh \left(\psi _1\ln {10\over 9}\right)\right]\right ].\nonumber\\ \end{eqnarray}}} \tag{ 112 }$$

For i ⩾ 2 we derive

$${\fontsize{8.6}{11}{\begin{eqnarray} a_i&=&{1\over j_i\nu _i}\int _{\nu _i\ln {10\over 9}}^{\nu _i\ln 2}dx\, e^{(2\eta -1)x\over 2\nu _i}\sin x\nonumber\\ &=&{1\over j_i\big[\nu _i ^2+(\eta -{1\over 2})^2\big]}\left [2^{\eta -{1\over 2}}\left[\left(\eta -{1\over 2}\right)\sin (\nu _i\ln 2)- \nu _i\cos (\nu _i\ln 2)\right]\nonumber\\ &&-\left({10\over 9}\right)^{\eta -{1\over 2}}\left[\left(\eta -{1\over 2}\right)\sin \left(\nu _i\ln {10\over 9}\right)- \nu _i\cos \left(\nu _i\ln {10\over 9}\right)\right]\right ].\nonumber\\ \end{eqnarray}}} \tag{ 113 }$$

In the last column of Table 1 we calculate the first five expansion coefficients a_i.

The stationary solution (110) for α = 0 and mono-energetic injection with rate Q₀ at momentum p_I finally reads

$$\begin{eqnarray} F(p,r)&=&{f(p,r)\over {3Q_0R\over 2p_IV_0}}\nonumber\\ &=&\left({r\over R}\right)^{\eta +{1\over 2}}\right [{2a_1\over 5}\sinh \left(\psi _1\ln {R\over r}\right) \right(H[p_I-p]\nonumber\\ &&+\left({p\over p_I}\right)^{-5}H[p-p_I]\left)\nonumber\\ &&+\sum _{i=2}^\infty {a_i\over \mu _i}\sin \left(\nu _i\ln {R\over r}\right)\left (\left({p\over p_I}\right)^{\mu _i-{5\over 2}}H[p_I-p]\nonumber\\ &&+\left({p\over p_I}\right)^{-(\mu _i+{5\over 2})}H[p-p_I]\right)\right ]. \end{eqnarray} \tag{ 114 }$$

In Figures 1 and 2, we show the momentum spectra of the accelerated STPs inside (r = 0.7R) and outside (r = 0.1R) of the adopted spatial source distribution. In each plot we show the total momentum distribution and the contribution from the higher eigenvalues i ⩾ 2. The momentum spectra are flat below p_i. Immediately above p_i, due to the influence of higher-order eigenfunctions, they exhibit a slight falloff, until at large enough momenta they flatten toward the p⁻⁵-spectrum. This is best visible in Figure 2.

Zoom In Zoom Out Reset image size

The radial profile of the accelerated STPs at p = 4p_I and p = 10p_I shown in Figures 3 and 4, respectively, first decreases because of the reflecting boundary condition at r = 0.1R, before they increase significantly at r ⩾ 0.4R toward the outer ring source distribution. At small radial distances the radial profile suffices a bit from the well-known mathematical Gibbs phenomenon of Fourier series expansion (Carslaw 1950, Ch. 9), providing negative particle densities as we have truncated the sum over i in Equation (114) at i = 10. This mathematical problem can be easily alleviated by modifying the radial boundary conditions requiring constants values M(r₀) = M(R) = c₀ at the edges and additionally $(dM/dr)_{r_0}=0$. The profile in Figure 2 should therefore be regarded as the relative radial variation of accelerated STPs.

Table 1. The Spatial Eigenvalues λ_i, Spectral Indices μ_i, and Dimensionless Expansion Coefficients a_i for η = 3.056 (M_A = 1.35) and R = 10r₀

i	λi	μi	ai
1	1.29 × 10 −6	2.5	1.17 × 10 −5
2	8.26	3.80	1.27
3	20.53	5.16	−0.28
4	40.99	6.84	−0.11
5	69.63	8.67	0.30
6	106.45	10.56	−0.31
7	151.45	12.49	0.11
8	204.64	14.44	−0.03
9	266.01	16.41	−0.17
10	335.56	18.38	0.13

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We explored the acceleration of suprathermal charged particles in the heliosphere under pressure balance conditions including for the first time the radial spatial particle diffusion and convection in the solar wind. In particular, we were interested in deriving the physical conditions for which the stationary phase space distribution of STPs approaches the power-law distribution f∝p⁻⁵, often seen in spacecraft observations. We first showed that power-law momentum distributions only result if the radial spatial diffusion and the momentum diffusion of particles are caused by the contributions from compressional waves. Inclusions of other momentum-dependent spatial diffusion coefficients will lead to non-power-law STP momentum distributions, which are not observed.

For separable source distributions in momentum and position we analytically solved the stationary particle transport equation for a radially constant solar wind speed V₀ and a momentum-independent radial spatial diffusion coefficient, by expanding the solution in terms of the eigenfunctions of the spatial operator. For either reflecting or free-escape spatial boundary conditions, these spatial eigenfunctions form a complete orthonormal basis. We demonstrated that the resulting stationary solution at any position within the finite heliosphere is the superposition of an infinite sum of power laws in momentum below and above the (assumed mono-momentum) injection momentum p_I. The smallest spatial eigenvalue determines the flattest power law, to which the full stationary solution approaches at large (p ≫ p_I) and small (p ≪ p_I) enough momenta.

Because of the pressure balance condition, intrinsically relating the momentum and spatial diffusion coefficients with the rate of adiabatic cooling, the most sensitive parameters determining the solution are (1) the Mach number of the solar wind M_A = V₀/V_c, where V_c denotes the speed of the solar wind compressional waves, (2) the location of the inner and outer spatial boundaries r₀ and R, and (3) the location and extent of the adopted outer ring spatial source distribution. By studying in detail four different spatial boundary conditions (reflecting and/or free-escape at the inner and outer boundaries r₀ and R), we demonstrated that only for the case of a reflecting inner and a free-escape outer boundary does one small eigenvalue λ₁ < (15M²_A/4)[1 − 0.4M⁻²_A]² exist, yielding a flattest power-law index μ₁ = 5, so that the power-law distribution f∝p⁻⁵ results at sufficiently large momentum values. The other three spatial boundary conditions imply steeper momentum spectra.

We illustrated our results for the favorite reflecting inner and free-escape outer boundary by adopting a uniform outer ring spatial source distribution (between 0.5R and 0.9R) and a Mach number of M_A = 1.35. The calculated momentum spectra inside and outside this source distribution indeed show the flattening toward the p⁻⁵-spectrum at high momenta. The radial profiles of STPs first slightly decrease because of the reflecting boundary condition at r = 0.1R before they increase significantly at r ⩾ 0.4R. By comparing with observed radial gradients of STPs it should be possible to determine the free parameters of our solution.

Our explanation for the f∝p⁻⁵ power-law momentum spectrum of STPs is based on the pressure balance concept, implying that the timescale of effective turbulent momentum diffusion of particles is kept in balance with the timescale of adiabatic cooling in the expanding solar wind. This is reasonable in the supersonic part of the solar wind, where adiabatic deceleration dominates all other loss mechanisms. However, f∝p⁻⁵ momentum spectra are also frequently observed in the heliosheath where adiabatic cooling is negligible. In the heliosheath, particle losses due to charge exchange and particle escape from active acceleration regions are major loss mechanisms that can also limit the development of the power-law slope (see Zhang & Schlickeiser 2012 for details). It is the balance between stochastic acceleration and (any) loss mechanisms that may result in a near p⁻⁵-momentum spectrum.

The work of T.A. and R.S. was partially supported by the Deutsche Forschungsgemeinschaft through grants Schl 201/23-1 and Schl 201/25-1. M.Z. acknowledges partial support by NASA grants NNX08AP91G, NNX09AG29G, and NNX09AB24G.

In this case we set

$$\begin{equation} \nu ^2=\Lambda _i-(\eta +(1/2))^2\ge 0. \end{equation} \tag{ A1 }$$

A.1. Cases A and D

For the free-escape boundary conditions at R in these two cases, we obtain

$$\begin{equation} E(r)=c_1r^{1/2}\sin [\nu \ln (R/r)]. \end{equation} \tag{ A2 }$$

A.1.1. Case A

Free-escape at the inner boundary r₀ in this case then requires

$$\begin{equation} \nu \ln (R/r_0)=n\pi ,\; \; n=1,2,3,\ldots , \end{equation} \tag{ A3 }$$

implying

$$\begin{eqnarray} M(r)&=&r^{\eta }E(r)=\sum _{n=1}^\infty c_nM_n(r)\nonumber\\ &=&r^{\eta +(1/2)}\sum _{n=1}^\infty c_n\sin \left(n\pi {\ln (R/r)\over \ln (R/r_0)}\right). \end{eqnarray} \tag{ A4 }$$

Likewise, Equation (A3) provides for the eigenvalues according to Equations (A1) and (63)

$$\begin{eqnarray} \lambda _n&=&\beta _0\left[(\eta +(1/2))^2+{n^2\pi ^2\over \ln ^2(R/r_0)}\right]\nonumber\\ &=&{3\over 5M_A^2}\left[{\big(5M_A^2-2\big)^2\over 4}+{n^2\pi ^2\over \ln ^2(R/r_0)}\right]. \end{eqnarray} \tag{ A5 }$$

It remains to determine the expansion coefficients c_n in Equation (A4) from the initial condition (59), i.e.,

$$\begin{equation} \sum _{n=1}^\infty c_nM_n(r)={3rq_1(r)\over V_0}. \end{equation} \tag{ A6 }$$

With solution (A4) used in Equation (66) we obtain

$$\begin{equation} j_n={1\over 2}\ln (R/r_0). \end{equation} \tag{ A7 }$$

Applying the orthonormality condition (65) to the initial condition then yields for arbitrary but specified special source distributions the expansion coefficients

$$\begin{eqnarray} c_n&=&{6\over V_0\ln (R/r_0)}\int _{r_0}^Rdr\, r^{-(1+2\eta)}q_1(r)M_n(r)\nonumber\\ &=&{6\over V_0\ln (R/r_0)}\int _{r_0}^Rdr\, r^{-\eta -(1/2)}q_1(r) \sin \left[n\pi {\ln (R/r)\over \ln (R/r_0)}\right].\nonumber\\ \end{eqnarray} \tag{ A8 }$$

A.1.2. Case D

Instead of Equation (74) we find here the solution

$$\begin{equation} M(r)=a_1r^{\eta +(1/2)}\sin [\nu \ln (R/r)], \end{equation} \tag{ A9 }$$

implying

$$\begin{eqnarray} {dM\over dr}&=&a_1r^{\eta -(1/2)}[(\eta +(1/2))\sin [\nu \ln (R/r)]\nonumber\\ &&-\nu \cos [\nu \ln (R/r)] ]. \end{eqnarray} \tag{ A10 }$$

The inner reflecting boundary $(dM/dr)_{r_0}=0$ then yields the transcendental equation

$$\begin{equation} \tan [\nu \ln (R/r_0)]={\nu \over \eta +(1/2)}, \end{equation} \tag{ A11 }$$

which has an infinite number of solutions. For (η + (1/2))ln (R/r₀) ≫ 1 the solutions are given by

$$\begin{equation} \nu _n=(n-1)\pi \left[1+{1\over (\eta +(1/2))\ln (R/r_0)}\right],\,\, n=2,3,\ldots \end{equation} \tag{ A12 }$$

A.2. Cases B and C

Instead of Equation (89) we obtain here

$$\begin{equation} {dM\over dr}=a_3r^{\eta -(1/2)}\sin \left[\nu \ln (R/r)\right]. \end{equation} \tag{ A13 }$$

A.2.1. Case B

The inner reflecting boundary condition $(dM/dr)_{r_0}=0$ is then fulfilled by the eigenvalues

$$\begin{equation} \nu \ln (R/r_0)=n\pi , \;\; n=1,2,3,\ldots , \end{equation} \tag{ A14 }$$

providing the same eigenvalue equation (A5) as in case A. Here the eigensolutions read

$$\begin{eqnarray} M(r)&=&\sum _{n=1}^\infty c_nM_n(r)\nonumber\\ &=&r^{\eta +(1/2)}\sum _{n=1}^\infty c_n \left [(\eta +(1/2))\sin \left[n\pi {\ln (R/r)\over \ln (R/r_0)}\right]\nonumber\\ &&+{n\pi \over \ln (R/r_0)}\cos \left[n\pi {\ln (R/r)\over \ln (R/r_0)}\right]\right ]. \end{eqnarray} \tag{ A15 }$$

The expansion coefficients c_n are given by using the eigensolutions (A15) in Equation (A8).

A.2.2. Case C

Instead of Equation (91) we obtain here

$$\begin{eqnarray} M(r)&=&a_5r^{\eta +(1/2)}[(\eta +(1/2))\sin [\nu \ln (R/r)]\nonumber\\ &&+\nu \cos [\nu \ln (R/r)]], \end{eqnarray} \tag{ A16 }$$

so that the inner free-escape boundary condition M(r₀) = 0 yields the transcendental equation

$$\begin{equation} \tan [\nu \ln (R/r_0)]=-{\nu \over \eta +(1/2)}, \end{equation} \tag{ A17 }$$

which has an infinite number of solutions for ν. For ln (R/r₀) ≫ 2/(1 + 2η) one finds approximately

$$\begin{equation} \nu _n\simeq (n-1)\pi \left[1-{2\over (1+2\eta)\ln (R/r_0)}\right],\; n=2,3,\ldots \end{equation} \tag{ A18 }$$