EVOLUTION OF THE PROTON VELOCITY DISTRIBUTION DUE TO STOCHASTIC HEATING IN THE NEAR-SUN SOLAR WIND

Protons in the solar wind within a few solar radii of the Sun's surface undergo an increase in temperature perpendicular to the magnetic field ${T}_{p\perp }$ (Kohl et al. 1998; Cranmer et al. 1999; Esser et al. 1999). Additionally, as the solar wind expands through the heliosphere, the total proton temperature T_p does not evolve as it would for a spherically expanding adiabatic fluid, ${T}_{p}\propto {r}^{-4/3}$, but instead decreases more slowly (Hartle & Barnes 1970). It is still unclear which mechanisms are responsible for this observed behavior.

Many studies have proposed an injection of energy into the solar wind, typically across some span of radial distance, with potential sources including low-frequency Alfvénic turbulence (Matthaeus et al. 1999; Cranmer et al. 2007), ion cyclotron waves (Cranmer 2000; Hollweg & Isenberg 2002), shocked compressive modes (Bruner 1978), and magnetic reconnection, nano-flares, and other impulsive events (Parker 1987; Cargill & Klimchuk 2004; Drake et al. 2009). Some authors have suggested that the dissipation of this energy is sufficient to accelerate the solar wind through the Sun's gravitational well and then continue heating the solar wind as it progresses through the heliosphere (Matthaeus et al. 1999; Smith et al. 2001). Other models contend that constructing the energy transport equations for non-Gaussian velocity distributions, as are expected to arise in the nearly collisionless solar wind, may be sufficient to explain radial temperature profiles without the need for heating and acceleration from energy dissipation (Scudder 1992, 2015). Identifying signatures of these proposed mechanisms that can be measured in situ by spacecraft is a necessary step for determining which mechanisms are responsible for governing the solar wind's evolution.

In this work, we focus on stochastic proton heating induced by strong, low-frequency, Alfvén-wave turbulence. A body of literature concerning stochastic heating has shown that fluctuations with amplitudes above a critical threshold are able to produce perpendicular proton heating (McChesney et al. 1987; Chen et al. 2001; Johnson & Cheng 2001; Chaston et al. 2004; Voitenko & Goossens 2004; Fiksel et al. 2009). Because perpendicular ion heating is observed in the solar corona and solar wind, and because Alfvén-wave turbulence is a dominant component of solar-wind turbulence, stochastic proton heating has been proposed as a candidate for explaining the observed coronal and solar-wind temperature profiles. Chandran et al. (2010) developed a phenomenological description of the heating rates due to Alfvén-wave turbulence with scale lengths comparable to the ion gyroscale, and along with Chandran (2010), showed that such a mechanism could plausibly explain remote temperature observations in coronal holes. Chandran et al. (2011) incorporated stochastic proton heating into a steady state two-fluid kinetic model of the fast solar wind, finding that stochastic proton heating supplied a dominant contribution to the total turbulent heating rate. van der Holst et al. (2014) incorporated stochastic heating into three-dimensional, two-fluid, numerical simulations of the solar wind and found promising agreement with observations. Bourouaine & Chandran (2013) found that heating rates from Helios observations were consistent with stochastic proton heating occurring in fast-wind streams between 0.29 and 0.64 au. Chandran et al. (2013) showed that stochastic heating can explain the observed alpha-particle-to-proton temperature ratio (Kasper et al. 2013) and how this temperature ratio depends on the thermal-to-magnetic-pressure ratio β and relative drift velocity ${\rm{\Delta }}{U}_{\alpha p}$.

We extend these previous works by evaluating the effects of stochastic heating on the evolution of the proton velocity distribution in a model fast-wind stream in the inner heliosphere under the assumption that $\beta \ll 1$. We find that stochastic heating significantly alters the distribution function from an (assumed) initial Gaussian to a flattop distribution with steep tails. These non-Gaussian features are accurately described by a modified Moyal distribution. The presence or absence of this type of distribution function in low-β regions could serve as a test for the importance of stochastic heating in the solar wind. We describe the equations of our model in Section 2 and present numerical solutions to these equations in Section 3.

We model the evolution of the proton velocity distribution f within a thin, open magnetic flux tube centered on a radially oriented background magnetic field. We start with the gyroaveraged kinetic equation (Kulsrud 1983's Equation (37)) and drop the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ terms, leaving

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+{v}_{\parallel }\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}f-\displaystyle \frac{1}{2}{v}_{\perp }{v}_{\parallel }({\rm{\nabla }}\cdot \hat{{\boldsymbol{b}}})\displaystyle \frac{\partial f}{\partial {v}_{\perp }}\\ &&\qquad +\left(\displaystyle \frac{1}{2}{v}_{\perp }^{2}{\rm{\nabla }}\cdot \hat{{\boldsymbol{b}}}+\displaystyle \frac{q}{m}{E}_{\parallel }\right)\displaystyle \frac{\partial f}{\partial {v}_{\parallel }}=0,\end{eqnarray} \tag{ 1 }$$

where q is the proton charge, m the proton mass, $\hat{{\boldsymbol{b}}}={\boldsymbol{B}}/B,{\boldsymbol{B}}$ is the magnetic field, and ${v}_{\parallel }$ and ${v}_{\perp }$ are the proton velocity components parallel and perpendicular to $\hat{{\boldsymbol{b}}}$. We assume that the proton distribution function is gyrotropic. The proton number density n and average radial flow velocity U are given by $n=\int {d}^{3}{\boldsymbol{v}}f$ and ${nU}=\int {d}^{3}{\boldsymbol{v}}{v}_{\parallel }f$.

We define the reduced distribution function

$$\begin{eqnarray}&&g({v}_{\perp })=2\pi {\displaystyle \int }_{-\infty }^{\infty }{{dv}}_{\parallel }f({v}_{\perp },{v}_{\parallel }).\end{eqnarray} \tag{ 2 }$$

We assume that the solar-wind outflow is supersonic and that, as a consequence, the average parallel proton velocity is approximately U at each ${v}_{\perp }$:

$$\begin{eqnarray}&&2\pi {\displaystyle \int }_{-\infty }^{\infty }{{dv}}_{\parallel }f({v}_{\perp },{v}_{\parallel }){v}_{\parallel }={Ug}.\end{eqnarray} \tag{ 3 }$$

Upon integrating Equation (1) over ${v}_{\parallel }$ and making use of Equation (3), we obtain an equation for the time evolution of g:

$$\begin{eqnarray}&&\displaystyle \frac{\partial g}{\partial t}+\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}({Ug})-\displaystyle \frac{1}{2}{v}_{\perp }U({\rm{\nabla }}\cdot \hat{{\boldsymbol{b}}})\displaystyle \frac{\partial g}{\partial {v}_{\perp }}=0.\end{eqnarray} \tag{ 4 }$$

We assume a steady state and that the magnetic field is nearly radial, which enables us to rewrite Equation (4) in the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}({Ug})+\displaystyle \frac{{v}_{\perp }U}{2B}\displaystyle \frac{\partial B}{\partial r}\displaystyle \frac{\partial g}{\partial {v}_{\perp }}=0.\end{eqnarray} \tag{ 5 }$$

To incorporate stochastic heating, we add a perpendicular-kinetic-energy diffusion term,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {E}_{\perp }}{D}_{E}\displaystyle \frac{\partial }{\partial {E}_{\perp }}=\displaystyle \frac{1}{{v}_{\perp }}\displaystyle \frac{\partial }{\partial {v}_{\perp }}\displaystyle \frac{{D}_{E}}{{m}^{2}{v}_{\perp }}\displaystyle \frac{\partial }{\partial {v}_{\perp }},\end{eqnarray} \tag{ 6 }$$

to the right-hand side of Equation (5):

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}({Ug})+\displaystyle \frac{{v}_{\perp }U}{2B}\displaystyle \frac{\partial B}{\partial r}\displaystyle \frac{\partial g}{\partial {v}_{\perp }}=\displaystyle \frac{1}{{v}_{\perp }}\displaystyle \frac{\partial }{\partial {v}_{\perp }}\displaystyle \frac{{D}_{E}}{{m}^{2}{v}_{\perp }}\displaystyle \frac{\partial g}{\partial {v}_{\perp }}.\end{eqnarray} \tag{ 7 }$$

For the diffusion coefficient D_E, we use the expression in Equation (17) from Chandran et al. (2010) multiplied by an exponential suppression term,

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{E}}{{m}^{2}{v}_{\perp }}\sim \delta {v}_{\rho }^{3}{{\rm{\Omega }}}_{p}\mathrm{exp}\left[-\displaystyle \frac{{c}_{2}}{\epsilon }\right],\end{eqnarray} \tag{ 8 }$$

where $\delta {v}_{\rho }$ is the rms amplitude of the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity fluctuations at scale $\rho ={v}_{\perp }/{{\rm{\Omega }}}_{p}$, and ${{\rm{\Omega }}}_{p}={{qB}}_{0}/{mc}$ is the proton gyrofrequency. The effectiveness of stochastic heating depends strongly on the stochasticity parameter $\epsilon \equiv \delta {v}_{\rho }/{v}_{\perp }$. When $\epsilon \ll 1$, a proton's gyromotion is only weakly perturbed by the gyroscale fluctuations and the proton's magnetic moment is nearly conserved. When $\epsilon \sim 1$, the gyroscale electric-field fluctuations strongly distort a proton's gyromotion, which enables the time-varying electrostatic potential to cause perpendicular proton heating. Chandran et al. (2010) incorporated the exponential suppression term $\mathrm{exp}(-{c}_{2}/{\epsilon }_{p})$ into their expression for the stochastic heating rate of the entire distribution, where ${\epsilon }_{p}=\epsilon {| }_{{v}_{\perp }={w}_{p}}$ and ${w}_{p}={(2{k}_{B}{T}_{\perp }/m)}^{1/2}$. In Equation (8), we incorporate an analogous but velocity-dependent exponential suppression term directly into the velocity-diffusion coefficient. The value of c₂, taken to be 0.2 in this work, modifies the amplitude threshold above which stochastic heating becomes effective. This choice of c₂ is motivated by values found in test-particle simulations (Xia et al. 2013) and inferred from Helios 2 observations (Bourouaine & Chandran 2013). We note that D_E was derived in the low-β limit, and therefore the model constructed below is also restricted to solar-wind streams with small β.

To solve Equation (7), we first specify radial profiles for the magnetic field B₀, solar-wind velocity U, and diffusion coefficient ${D}_{E}({v}_{\perp },r)$. We adopt the magnetic-field profile used by Hollweg & Isenberg (2002),

$$\begin{eqnarray}&&{B}_{0}=\left[\displaystyle \frac{6}{{x}^{6}}+\displaystyle \frac{1.5}{{x}^{2}}\right]{\rm{Gauss}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&x=\displaystyle \frac{r}{{R}_{{\rm{s}}}},\end{eqnarray} \tag{ 10 }$$

and the density profile of Feldman et al. (1997) with the ${x}^{-2}$ term added by Chandran & Hollweg (2009),

$$\begin{eqnarray}&&n=\left(\displaystyle \frac{3.23\times {10}^{8}}{{x}^{15.6}}+\displaystyle \frac{2.51\times {10}^{6}}{{x}^{3.76}}+\displaystyle \frac{1.85\times {10}^{5}}{{x}^{2}}\right){\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ 11 }$$

Upon multiplying Equation (7) by ${v}_{\perp }$ and integrating over ${v}_{\perp }$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}\left(\displaystyle \frac{{nU}}{{B}_{0}}\right)=0,\end{eqnarray} \tag{ 12 }$$

which can also be deduced by simply noting that the magnetic flux ${B}_{0}A$ and "proton flux" ${nUA}$ through a narrow flow/flux tube of cross sectional area A are both independent of r in steady state. From Equation (12) we obtain

$$\begin{eqnarray}&&U=9.25\times {10}^{12}\ \displaystyle \frac{{B}_{0}}{n}\ \mathrm{cm}\;{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 13 }$$

where the numerical coefficient has been chosen so that U extrapolates to a value of 750 km s⁻¹ at 1 au. Profiles for B₀ and U are plotted in Figure 1, and n is plotted in panel (a) of Figure 3.

Zoom In Zoom Out Reset image size

We define $\delta {v}_{\lambda }$ to be the rms amplitude of the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity at a scale length λ measured perpendicular to ${{\boldsymbol{B}}}_{0}$. Based on theoretical predictions (Boldyrev 2006; Chandran et al. 2015), numerical simulations (Boldyrev et al. 2011; Perez et al. 2012), and solar-wind observations (Podesta 2009), we set

$$\begin{eqnarray}&&\delta {v}_{\lambda }=\delta {v}_{0}{\left(\displaystyle \frac{\lambda }{{L}_{0}}\right)}^{1/4}\end{eqnarray} \tag{ 14 }$$

when λ is in the inertial range, where L₀ is the outer-scale correlation length. We take L₀ to be proportional to the radius of the flux tube in which the turbulence is embedded. Specifically, we use Equation (33) from Chandran & Hollweg (2009), setting ${L}_{0}=5000{[7.5{\rm{G}}/{B}_{0}(r)]}^{1/2}\;{\rm{km}}$.

We assume that Equation (14) can be extrapolated, at least approximately, all the way to the proton-gyroradius scale ρ, which yields

$$\begin{eqnarray}&&\delta {v}_{\rho }=\delta {v}_{0}{\left(\displaystyle \frac{\rho }{{L}_{0}}\right)}^{1/4}.\end{eqnarray} \tag{ 15 }$$

This assumption neglects the possible back reaction of the heating process on the turbulent power spectrum. For example, if stochastic heating drains energy from the turbulent cascade at scales $\sim \rho$, it could reduce $\delta {v}_{\rho }$ below the value in Equation (15). Chandran (2010) modeled this back reaction by including a reduction factor in the expression for $\delta {v}_{\rho }$. In the present paper, we find that even without such a reduction factor the ion heating rate remains less than the cascade power except within a narrow radial interval near the innermost radius r₀ in our numerical solutions (see Equation (23) and Figure 4), which is dominated by artifacts associated with our imposition of a Gaussian velocity distribution at r₀ (see discussion below). We have thus refrained from multiplying the right-hand side of Equation (15) by a reduction factor in order to keep our model as simple as possible.

For the radial profile of $\delta {v}_{0}$, we use the reflection-driven turbulence model of Chandran & Hollweg (2009), converting their expression for the outward-propagating Heinemann–Olbert variable ${g}_{{\rm{HO}}}$ into the velocity amplitude

$$\begin{eqnarray}&&\delta {v}_{0}=\displaystyle \frac{{g}_{a}}{2}\left(\displaystyle \frac{{\eta }^{1/4}}{1+{\eta }^{1/2}}\right){\left(\displaystyle \frac{{v}_{A}}{{v}_{{Aa}}}\right)}^{\chi /2},\end{eqnarray} \tag{ 16 }$$

where ${v}_{A}={B}_{0}/{(4\pi {mn})}^{1/2}$ is the Alfvén speed, $\eta \equiv n/{n}_{a}$, and ${v}_{{Aa}},{g}_{a}$, and n_a are values of ${v}_{A},{g}_{{\rm{HO}}}$, and n at the Alfvén critical point, which is at $r=11.1\ {R}_{{\rm{S}}}$ in our model. The factor χ models the reduction of the efficiency of wave reflection for waves with periods below ∼1 hr. Both values of χ used in Chandran & Hollweg (2009), $\chi =1$ and 0.65, lead to qualitatively similar heating rates and velocity distributions in our model. In this work, we only present results from the $\chi =0.65$ model, which gives turbulence amplitudes closer to observed values. The value of g_a is set to $7.2\times {10}^{7}\;\mathrm{cm}\;{{\rm{s}}}^{-1}$ to match constraints from Helios 2 (Marsch et al. 1982) and Ultraviolet Coronagraph Spectrometer (UVCS) observations (Esser et al. 1999). Profiles for $\delta {v}_{0}$ and $\delta {v}_{\rho }$ are shown in Figure 1. With the above assumptions, we can rewrite Equation (8) as

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{E}}{{m}^{2}{v}_{\perp }}={{\rm{\Omega }}}_{p}^{1/4}\delta {v}_{0}^{3}{\left(\displaystyle \frac{{v}_{\perp }}{{L}_{0}}\right)}^{3/4}\mathrm{exp}\left(\displaystyle \frac{-{c}_{2}{v}_{\perp }^{3/4}{({{\rm{\Omega }}}_{p}{L}_{0})}^{1/4}}{\delta {v}_{0}}\right).\end{eqnarray} \tag{ 17 }$$

We solve Equation (7) numerically using the Crank–Nicholson finite-difference method. The reduced distribution function $g({v}_{\perp },r)$ is evaluated at 10⁵ points logarithmically spaced between ${r}_{0}=4\ {R}_{{\rm{S}}}$ and ${r}_{{\rm{f}}}=30\ {R}_{{\rm{S}}}$. The value of r₀ is chosen so that the solar wind is supersonic and weakly collisional. Based upon UVCS observations (Kohl et al. 1998; Esser et al. 1999), we take ${T}_{\perp 0}$ to be $2\times {10}^{6}\;{\rm{K}}$. We also take $g({v}_{\perp },{r}_{0})$ to be a Gaussian. The initial density, ${n}_{0}=2.51\times {10}^{4}\;{\mathrm{cm}}^{-3}$, is found through evaluation of Equation (11) at r₀.

The velocity grid has 10³ linearly spaced points between ${v}_{\perp ,0}=0$ and ${v}_{\perp ,{\rm{f}}}=10{w}_{0}$, where ${w}_{0}=1.82\times {10}^{7}\;\mathrm{cm}\;{{\rm{s}}}^{-1}$ is the thermal speed of the Gaussian distribution at r₀. We set $\partial g/\partial {v}_{\perp }=0$ at ${v}_{\perp }=0$ and $g({v}_{\perp ,{\rm{f}}})=0$ at ${v}_{\perp }={v}_{\perp ,{\rm{f}}}$. The choice of velocity and radial resolutions is sufficient to remove the spurious decaying oscillations that arise for under-resolved Crank–Nicholson solutions. To check that our outer boundary condition does not significantly alter our results, we recalculate our numerical solution with ${v}_{\perp ,{\rm{f}}}=15{w}_{0}$ and $20{w}_{0}$, retaining the same resolution in ${v}_{\perp }$. These extended models (not shown) produce virtually the same reduced distribution functions as the model with ${v}_{\perp ,{\rm{f}}}=10{w}_{0}$. We have verified that our numerical method conserves particles.

We solve Equation (7) with D_E given by Equation (17) and with D_E set equal to zero for comparison. We present contour plots of $g({v}_{\perp },r)$ in Figure 2 for the ${D}_{E}=0$ (left panel) and ${D}_{E}\ne 0$ (right) cases. When ${D}_{E}=0$, particles shift from high to low ${v}_{\perp }$ as a result of magnetic moment conservation. When D_E is given by Equation (17), the narrowing of g is arrested and reversed between 7 and $10\;{R}_{{\rm{S}}}$. At larger distances, particles do lose energy as they flow away from the Sun, but at a much slower rate than in the ${D}_{E}=0$ case.

Zoom In Zoom Out Reset image size

In Figure 3, we plot radial profiles of the density

$$\begin{eqnarray}&&n(r)={\displaystyle \int }_{0}^{{{\rm{v}}}_{\perp },{\rm{f}}}{{dv}}_{\perp }{v}_{\perp }g({v}_{\perp },r),\end{eqnarray} \tag{ 18 }$$

perpendicular temperature

$$\begin{eqnarray}&&{T}_{\perp }(r)=\displaystyle \frac{m}{2{k}_{B}n(r)}{\displaystyle \int }_{0}^{{{\rm{v}}}_{\perp },{\rm{f}}}{{dv}}_{\perp }{v}_{\perp }^{3}g({v}_{\perp },r),\end{eqnarray} \tag{ 19 }$$

and excess kurtosis

$$\begin{eqnarray}&&\kappa (r)=\displaystyle \frac{n(r){\displaystyle \int }_{0}^{{{\rm{v}}}_{\perp },{\rm{f}}}{{dv}}_{\perp }{v}_{\perp }^{5}g({v}_{\perp },r)}{{\left({\displaystyle \int }_{0}^{{{\rm{v}}}_{\perp },{\rm{f}}}{{dv}}_{\perp }{v}_{\perp }^{3}g({v}_{\perp },r)\right)}^{2}}-3.\end{eqnarray} \tag{ 20 }$$

The density profiles for the ${D}_{E}=0$ and ${D}_{E}\ne 0$ cases are the same and in agreement with the input density profile, Equation (11), as they must be since Equation (11) follows from Equations (7), (9), and (13) via Equation (12). The perpendicular temperature ${T}_{\perp }$ initially decreases, before increasing to a peak of $2.5\times {10}^{6}\;{\rm{K}}$ at $8\ {R}_{{\rm{S}}}$. The temperature then falls monotonically to $1.2\times {10}^{6}\;{\rm{K}}$ at $30\ {R}_{{\rm{S}}}$. As foreshadowed by the left panel of Figure 2, ${T}_{\perp }$ for the ${D}_{E}=0$ case rapidly falls from r₀ outward, reaching $6.6\times {10}^{3}\;{\rm{K}}$ at $30\ {R}_{{\rm{S}}}$.

Zoom In Zoom Out Reset image size

A Gaussian distribution has no excess kurtosis: ${\kappa }_{{\rm{G}}}=0$. Our ${D}_{E}\ne 0$ distribution is platykurtic, with negative excess kurtosis. For the reduced distribution, κ departs from 0 rapidly before leveling off at $\kappa \approx -0.8$ near $5\ {R}_{{\rm{S}}}$. The value of κ then gradually increases for larger radial distances, but the distribution remains highly non-Gaussian over the entire radial range under consideration in this work.

We calculate the perpendicular heating rate per unit mass, ${Q}_{\perp }$, by multiplying Equation (7) by ${v}_{\perp }^{3}/2n$ and integrating over ${v}_{\perp }$, which gives

$$\begin{eqnarray}&&\displaystyle \frac{{{Bk}}_{B}U}{m}\displaystyle \frac{d}{{dr}}\left(\displaystyle \frac{{T}_{\perp }}{B}\right)={Q}_{\perp },\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{Q}_{\perp }=\displaystyle \frac{1}{2n}{\displaystyle \int }_{0}^{\infty }{{dv}}_{\perp }{v}_{\perp }^{2}\displaystyle \frac{\partial }{\partial {v}_{\perp }}\displaystyle \frac{{D}_{E}}{{m}^{2}{v}_{\perp }}\displaystyle \frac{\partial g}{\partial {v}_{\perp }}.\end{eqnarray} \tag{ 22 }$$

We plot ${Q}_{\perp }$ in Figure 4. For this figure, we approximate Equation (22) by replacing the upper limit of integration with ${v}_{\perp ,{\rm{f}}}$.

Zoom In Zoom Out Reset image size

Beyond the immediate neighborhood of r₀, as r increases ${Q}_{\perp }$ levels off near $5\ {R}_{{\rm{S}}}$ and then decreases monotonically. Our assumption of a Gaussian distribution at r₀ is for concreteness. We do not expect this assumption to hold in the corona, and therefore our model is inaccurate near r₀. To understand the effects of the boundary condition at r₀, we recalculate the numerical solution using ${r}_{0}\in [2,3,5]\ {R}_{{\rm{S}}}$, keeping ${T}_{\perp 0}=2\times {10}^{6}\;{\rm{K}}$, and using ${n}_{0}({r}_{0})$ from Equation (11). While the resulting heating rates differ near r₀, all of the solutions quickly converge to the same radial profile. A similar convergence is seen for ${T}_{\perp }$ and κ in Figure 4, as well as for the reduced distribution function g (not shown).

We compare the perpendicular proton heating rate ${Q}_{\perp }$ with the approximate total turbulent heating rate (per unit mass)

$$\begin{eqnarray}&&{Q}_{{\rm{turb}}}=\displaystyle \frac{{z}^{-}{({z}^{+})}^{2}}{4{L}_{0}}\end{eqnarray} \tag{ 23 }$$

in the model of Chandran & Hollweg (2009), where ${z}^{+}$ (${z}^{-}$) is the rms amplitude of the Alfvén-wave-like fluctuations that are outward-propagating (inward-propagating) when viewed in the local plasma frame. Except in the immediate vicinity of r₀, the perpendicular heating rate is less than the turbulent heating rate, indicating that not all of the cascade power is dissipated by protons via stochastic heating.

To characterize the shape of the velocity distribution, we consider two types of fits to g. First, we fit g using a least-squares method to a Gaussian of the form

$$\begin{eqnarray}&&{g}_{{\rm{G}}}({v}_{\perp },r)=\displaystyle \frac{2{n}_{{\rm{G}}}(r)}{{w}_{{\rm{G}}}^{2}(r)}\mathrm{exp}\left(-\displaystyle \frac{{v}_{\perp }^{2}}{{w}_{{\rm{G}}}^{2}(r)}\right)\end{eqnarray} \tag{ 24 }$$

at each radial grid point. The density and perpendicular temperature (${T}_{{\rm{G}},\perp }={{mw}}_{{\rm{G}}}^{2}/2{k}_{B}$) profiles calculated from the fitted ${g}_{{\rm{G}}}$ are shown in panels (a) and (b) of Figure 3. The Gaussian density ${n}_{{\rm{G}}}$ agrees with the input density, while ${T}_{{\rm{G}},\perp }$ overestimates ${T}_{\perp }$ by a factor that is $\simeq 1.7$ at $r\gtrsim 5{R}_{{\rm{s}}}$.

We next fit g to a modified Moyal distribution (Moyal 1955) of the form

$$\begin{eqnarray}&&{g}_{{\rm{Mo}}}=A(r)\mathrm{exp}\left(\displaystyle \frac{1}{2}\left[\displaystyle \frac{{v}_{\perp }^{2}}{{w}_{{\rm{Mo}}}^{2}(r)}-\mathrm{exp}\left(\displaystyle \frac{{v}_{\perp }^{2}}{{w}_{{\rm{Mo}}}^{2}(r)}\right)\right]\right).\end{eqnarray} \tag{ 25 }$$

As this distribution is not frequently employed in the solar-physics literature, we plot selected properties in Figure 5. In the left panel, ${g}_{{\rm{Mo}}}$ is plotted with fixed $A=1.0$ for ${w}_{{\rm{Mo}}}$ varying from 0.5 to 10. In the right panel, the root mean square velocity and excess kurtosis for both ${g}_{{\rm{G}}}$ and ${g}_{{\rm{Mo}}}$ are plotted as a function of w. In the Gaussian case, ${v}_{{\rm{rms}}}={w}_{{\rm{G}}}$, while for the modified Moyal distribution, ${v}_{{\rm{rms}}}\approx 0.877{w}_{{\rm{Mo}}}$. Both distributions have a constant $\kappa ,{\kappa }_{{\rm{G}}}=0$ and ${\kappa }_{{\rm{Mo}}}\approx -0.781$.

Zoom In Zoom Out Reset image size

At each radial grid point, we calculate the best-fit modified Moyal distribution ${g}_{{\rm{Mo}}}({v}_{\perp },r)$ for the reduced distribution $g({v}_{\perp },r)$. We then calculate the density, perpendicular temperature (${T}_{\mathrm{Mo},\perp }={{mv}}_{{\rm{rms}}}^{2}/2{k}_{B}$), and excess kurtosis ${\kappa }_{{\rm{Mo}}}$ by numerically integrating the fitted function and plot these quantities in panels (a)–(c) of Figure 3. The fitted Moyal density in panel (a) is in agreement with the input density, and overlaps with plots of n. The temperature corresponding to the best-fit Moyal distribution, ${T}_{\mathrm{Mo},\perp }$ panel (b), slightly underestimates ${T}_{\perp }$ between r₀ and $5\ {R}_{{\rm{S}}}$, beyond which ${T}_{\mathrm{Mo},\perp }$ and ${T}_{\perp }$ are in good agreement. The fitted excess kurtosis ${\kappa }_{{\rm{Mo}}}$ has a constant value of −0.781 for all radial distances, which compares well to κ calculated from g.

In Figure 6 we plot g and the associated fit functions ${g}_{{\rm{G}}}$ and ${g}_{{\rm{Mo}}}$ at eight radial distances. Near ${r}_{0}=4\ {R}_{{\rm{S}}}$, we see that the core of g has flattened, but that the tail is still Gaussian. By $5\ {R}_{{\rm{S}}}$, the Gaussian features of g have disappeared, and ${g}_{{\rm{Mo}}}$ is visually indistinguishable from g. When compared to the Gaussian ${g}_{{\rm{G}}},g$ has a much flatter core and a significantly steeper tail. This shape is a product of the rapid energy diffusion at small ${v}_{\perp }$ and negligible energy diffusion at superthermal perpendicular velocities.

Zoom In Zoom Out Reset image size

We have solved for the evolution of the reduced distribution function $g({v}_{\perp })$ in a model fast-solar-wind stream in the presence of stochastic heating under the assumption that β is small. Stochastic heating produces perpendicular proton heating, as is observed in the solar wind. It also causes g to develop significant non-Gaussian features, specifically a platykurtic (negative excess kurtosis) flattop distribution, which is well modeled by a modified Moyal distribution (Equation (25)). Detailed measurements of the proton distribution function from the upcoming Solar Probe Plus mission will provide a wealth of data that can be compared with the model results presented here, providing a test for the importance of stochastic heating in the solar-wind acceleration region.

We thank Phil Isenberg for helpful discussions. This work was supported by NSF grants AGS-1258998, AGS-1331355, and PHY-1500041 and NASA grant NNX15AI80G.