Numerical Modeling of Suprathermal Electron Transport in the Solar Wind: Effects of Whistler Turbulence with a Full Diffusion Tensor

The electron VDF in the solar wind consists of a Maxwellian core, a suprathermal halo, a field-aligned component strahl, and an energetic superhalo that deviates from the equilibrium. Whistler wave turbulence is thought to resonantly scatter the observed electron velocity distribution. Wave–particle interactions that contribute to Whistler wave turbulence are introduced into a Fokker–Planck kinetic transport equation that describes the interaction between the suprathermal electrons and the Whistler waves. A recent numerical approach for solving the Fokker–Planck kinetic transport equation has been extended to include a full diffusion tensor. Application of the extended numerical approach to the transport of solar wind suprathermal electrons influenced by Whistler wave turbulence is presented. Comparison and analysis of the numerical results with observations and diagonal-only model results are made. The off-diagonal terms in the diffusion tensor act to depress effects caused by the diagonal terms. The role of the diffusion coefficient on the electron heat flux is discussed.


Introduction
The electron velocity distribution function (VDF) in the solar wind differs significantly from a thermal equilibrium state and reveals an enhanced suprathermal tail (Feldman et al. 1974;Feldman et al. 1975;Marsch et al. 1982). Typically, the electron VDF can be modeled as being composed of four components: a Maxwellian "core," "halo," "strahl," and "superhalo" in ascending electron energy (Ergun et al. 1998;Wang et al. 2012). These components are significantly different from each other in their underlying distribution functions and relative number densities. The strahl is field-aligned anisotropic with respect to the magnetic field, while the remaining three components are nearly isotropic. At 1 au the Maxwellian core comprises the majority of the total number density (∼95%), and the halo, strahl, and superhalo together form the remaining < ∼ 5%. Although suprathermal electrons (the halo, strahl, and superhalo) comprise a minute fraction of the total electrons, they are chiefly responsible for most of the heat flux transported away from the Sun due to their high energy and mobility (Štverák et al. 2009).
Observations and theoretical studies of the solar wind electron velocity distribution have been implemented for decades (Montgomery et al. 1968;Feldman et al. 1975;Pilipp et al. 1987aPilipp et al. , 1987bMcComas et al. 1992). The most frequently adopted fitting model is the dual Maxwellian-kappa model which mainly possesses three ingredients: a Maxwellian core, a kappa halo, and a drifting-kappa strahl (Lazar et al. 2017). Suprathermal electrons are believed to originate in the solar corona (Pierrard et al. 1999;Štverák et al. 2008;Che & Goldstein 2014) and to escape from the solar corona along open interplanetary magnetic field lines. The strahl component is found in either the parallel or antiparallel magnetic field direction or in both in certain environments (Pilipp et al. 1987a;Anderson et al. 2012). The radial evolution of the pitch-angle width of the strahl component is contrary to an adiabatic expanding model (for instance, Lemons & Feldman 1983) which predicts that strahl electrons would be focused into a narrow beam (<1°at 1 au) due to magnetic moment conservation as the magnetic field strength weakens. Helios-1 observations from 0.3-1 au showed strahl widths from 5°to 60° (Pilipp et al. 1987a;Pilipp et al. 1987b). Based on ACE data, Anderson et al. (2012) reported that the strahl widths exihibit a range of at least 5°-90°. Ulysses observations showed that the average strahl pitch-angle width increases with heliocentric distance from ∼1 to 2.5 au (Hammond et al. 1996). Using Cassini electron measurements, Graham et al. (2017) extended the radial evolution of strahl width to a heliocentric distance of about 9 au and found that there is a constant rate of broadening of strahl pitch-angle distributions with heliocentric distance between ∼1 and 5.5 au, and that beyond this distance the strahl is likely to be completely scattered, to form part of the halo. Besides the strahl width, observations of the radial evolution of solar wind electrons from 0.3 to 4 au demonstrate that the relative number density (to the total number density) of the Maxwellian core remains unchanged with increasing heliocentric distance, while the number density of strahl electrons decreases, and halo electrons increases (Maksimovic et al. 2005;Štverák et al. 2009;Tao et al. 2016). In the latest observations from PSP, Halekas et al. (2020) found that near perihelion (∼0.17 au) the strahl is narrower and dominates the suprathermal fraction of the distribution. They also found that the halo almost disappears and its relative number density is much smaller than that at larger heliocentric distance (McComas et al. 1992), and is considerably smaller even than those reported at 0.3 au (Maksimovic et al. 2005;Štverák et al. 2009). This finding suggests that at ever closer distances, the solar wind electron VDF probably comprises only a Maxwellian core and a strahl.
The opposite radial evolution characteristic of the strahl and halo populations, as well as their both lying in the same energy range (10 2 -10 3 eV) suggests that they are essentially the same suprathermal component and that probably only the strahl component exists in the corona, and the halo electrons are pitch-angle scattered from the strahl by some mechanism as they propagate in the solar wind. The solar wind plasma is tenuous, strahl electrons are energetic, and Coulomb collisions are too infrequent to scatter suprathermal electrons efficiently into a halo that matches observations (Hammond et al. 1996;Vocks et al. 2005). At 1 au the mean free path is comparable to the typical length scale of the system (Štverák et al. 2008). It is necessary therefore to include wave-particle interactions. Vocks & Mann (2003), Vocks et al. (2005), and Vocks et al. (2008) considered the formation of the halo and strahl in the electron VDF in the solar corona and solar wind as a consequence of Whistler turbulence wave-particle interactions for electrons with a very broad energy range, even as high as hundreds of keV. They show that the quiet solar corona is able to produce suprathermal electrons by Whistler wave turbulence and that such an electron component should be present in the solar wind. By contrast, Pierrard et al. (2011) added a diffusion coefficient derived from Whistler wave turbulence to the exospheric model of Maksimovic et al. (1997), and found that the nonthermal tails of solar wind electron VDF emerge from an original Maxwellian distribution by wave-particle interactions associated with Whistler wave turbulence.
In the above models, electrons as a whole are subjected to certain scattering mechanisms collectively as they propagate in the solar wind. Alternatively, Kim et al. (2015) and Yoon et al. (2015), by assuming that the local solar wind electron VDF is a superposition of individual components (Maxwellian core, halo, and superhalo), proposed an asymptotic scattering theory for solar wind electrons in local equilibrium with plasma wave turbulence. The Maxwellian core does not experience collisionless scattering, the halo electrons only interact with Whistler wave turbulence, and the superhalo electrons interact only with Langmuir wave turbulence. Results from Yoon et al. 's (2013) asymptotic scattering theory were comparable with observations at 1 au. Boldyrev & Horaites (2019) presented an extended kinetic theory for strahl electrons scattered by multiple mechanisms, including both Coulomb collisions and wave-particle interactions of Whistler waves. They solve the drift-kinetic equation and obtain an analytic solution for the suprathermal electrons along the Parker-spiralled magnetic field lines; however, they simplified the kinetic transport equation and only considered the effect of the pitch-angle scattering term.
The latest PSP observations present a great challenge for previous models of solar wind acceleration and scattering (Halekas et al. 2020). Tang et al. (2020, Paper I) introduced a wave-particle interaction diffusion term into the kinetic transport equation that describes the interaction of only superathermal electrons with Whistler waves and developed a numerical method to solve an advection-diffusion-like Fokker-Planck kinetic equation in 3D phase space under the assumption that the diffusion tensor was purely diagonal. This numerical method was used to investigate the resonant waveparticle interaction of the suprathermal electrons in the solar wind. In the present paper, we extend the numerical approach to include the full kinetic diffusion tensor, including the offdiagonal terms whose magnitude is comparable to the diagonal ones (Tang et al. 2020) and hence are important to the kinetic transport equation. The temporal and radial evolution of the suprathermal electrons shows that resonant wave-particle interactions due to Whistler turbulence can significantly pitch-angle scatter a highly anisotropic field-aligned suprathermal electron VDF injected near the Sun into a nearly isotropic distribution at 1 au.
The layout of the present paper is as follows. Section 2 briefly reviews the suprathermal electron transport model and describes the extended kinetic transport equation with a full diffusion tensor. Section 3 presents numerical results and discussions of the electron VDFs and macroparameters. The influence of the diffusion coefficient on the radial evolution of the electron heat flux is also discussed. The summary and conclusions are given in Section 4.

Numerical Model in an Expanding Solar Wind Background
Details of the numerical model have been given in Tang et al. (2020), and we briefly review the basic features. With increasing heliocentric distance, the relative number density of strahl electrons decreases and that of halo electrons increases, while the relative number density of the core remains unchanged (Maksimovic et al. 2005;Štverák et al. 2009;Tao et al. 2016). Beyond about 5.5 au, the strahl is almost completely scattered into a halo (Graham et al. 2017). PSP observations show that near perihelion (∼0.17 au), the strahl is narrower and dominates the suprathermal fraction of the distribution, and the halo is almost absent (Halekas et al. 2020). We assume that the halo and strahl components are the same electron population, the halo being the result of pitchangle scattering of the strahl as it propagates in the solar wind. At an inner boundary of 0.1 au, it is assumed that there are only a Maxwellian core and a highly field-aligned strahl component. These different components interact with different scattering mechanisms individually since the distinct components may have different origins (Pierrard et al. 1999;Maksimovic et al. 2005). We simply assume, for the present, that a core Maxwellian distribution exists, does not experience scattering, and is maintained throughout the solar wind . We then follow the radial evolution of suprathermal electrons (at an inner boundary it is a strahl only) subject to certain scattering mechanisms in an expanding plasma and magnetic field background. Therefore the electron VDF comprises two parts, where f c is the Maxwellian core and f s the superimposed suprathermal electron distribution. At the inner boundary of 0.1 au, the suprathermal electron distribution is only the strahl, while at heliocentric distances >0.1 au the suprathermal electron distribution may include both halo and strahl electrons. The Maxwellian core VDF is isotropic, where n c and T c denote the number density and temperature of Maxwellian core electrons, respectively.
A narrow cold beam represents the strahl component injected at the inner boundary of 0.1 au and can be expressed as: where n s0 is the number density of suprathermal electron and v d is the drifting velocity in the solar wind frame. T s0 = T c0 ∼ 5 × 10 5 K is the temperature of the suprathermal electrons which corresponds to the temperature of the Maxwellian core at the inner boundary. The value of n s0 is chosen such that the calculated number density of suprathermal electrons n s at 1 au equals one twentieth of the Maxwellian core electrons since observations show that the suprathermal electrons comprise approximately only 5% of the total electron number density. The drifting velocity v d was chosen such that in the solar wind frame the suprathermal electron beam has a kinetic energy of 100 eV. At the outer boundary of r ∞ = 3 au, no electrons enter from infinity: The radially expanding background magnetic field corresponds to that of Adhikari et al. (2017), as adapted from Weber & Davis (1967), where r is in au, r a ∼ 0.05 au, ω a = 2.9 × 10 −6 rad s −1 , U = 400 km s −1 and B a = 2.08 × 10 3 nT, so that B(r = 1) = 0.45 nT is the magnetic field at 1 au. The rest of the parameters of the expanding background plasma are adopted from Horaites et al. (2015Horaites et al. ( , 2017 and Tang et al. (2020): where n c0 and T c0 are the Maxwellian core electron number density and temperature at 1 au, respectively, r is in au, and k B is the Boltzmann's constant.

The Full Form of the Kinetic Transport Equation with
Diffusion Tensor The kinetic equation for the suprathermal electron VDF f (r, v, μ, t) in a radial magnetic field and a constant radial flow is: where r is the heliocentric distance, v and μ the velocity magnitude and the cosine of pitch angle in the solar wind frame, and m is the electron mass. E ∥ represents the parallel electric field along the magnetic field and (δf/δt) sc the scattering term including wave-particle interactions, Coulomb collisions, and any other possible scattering mechanisms. For heliocentric distances 0.1 au, the solar wind is too tenuous for Coulomb collisions to produce the scattering necessary to match the suprathermal electron observations (Vocks et al. 2005). The scattering mechanisms must involve wave-particle interaction rather than collisions (Pagel et al. 2007;Saito & Gary 2007). Of the possible wave-particle interactions, the most likely candidate is resonant scattering associated with Whistler wave turbulence (Graham et al. 2017). We follow the assumption in Paper I that suprathermal electrons resonate only with Whistler wave turbulence (Pierrard et al. 2001;Kim et al. 2015;Yoon 2015). The scattering term hence becomes: which is chosen to be of the form (Schlickeiser 1989), The diffusion tensor for nonrelativistic electrons is expressed as (Steinacker &Miller 1992 andPierrard et al. 2011) 1 is the normalization constant and s = 3/2 ∼ 2 the spectral index of Whistler waves.
In order to avoid the singularity as D μv at μ tends to 0, we have assumed an empirical form of resonance broadening near μ=0 both in Paper I and the present paper and treat this region as follows: ⎨ ⎩ This ensures that the diffusion coefficients do not blow up as μ goes to 0. However, following Tang et al. (2020), we can alternatively use a spectral index for the Whistler wave pdf of s = 2 in the calculation. This choice avoids the singularity problem with D μv since μ goes to 0; the denominator of D μμ disappears when s = 2.
In Paper I, we showed that the diffusion coefficients (11a) are much larger than the realistic values obtained from observation.
In order to use them in our calculation, we multiply the diffusion coefficients by a small constant fraction f r , which is an artificial method to reduce the diffusion coefficients to possibly more realistic values. Furthermore, we also showed that D μμ > D vμ > D vv , showing that pitch-angle scattering dominates. In our previous paper, we kept only D μμ and D vv , but ignored the off-diagonal term D vμ despite its being larger than D vv . In the present paper, we retain the full diffusion tensor (11a) in our calculation. The inclusion of the off-diagonal terms in the kinetic transport equation gives rise to new effects on the transport of the suprathermal electron VDF as shown in Section 3. The offdiagonal diffusion terms D vμ , unlike the diagonal terms D μμ and D vv which are monotonic with respect to the variables r, v, and μ, are odd functions about μ = 0, as shown in Figure 1.
By introducing a new distribution function Y = Y(r, v, μ) such that f = Y/v 2 r 2 , we obtain an advection-diffusion-like kinetic transport equation with the full form diffusion tensor in the new variable Y in the three-dimensional phase space (r, v, μ): where the constant fraction f r was omitted for the sake of notational simplicity. The last three terms on the left-hand side are equations of characteristics, along which electrons move in the phase space (r, v, μ). There is one additional term in the third and fourth terms on the left-hand side that comes directly from D vv and D vμ , respectively. The third term in general is positive, and causes electrons to move in a direction of increasing velocity which therefore leads to an acceleration effect. The value of the third term depends on the pitch angle, and determines the direction in μ along which an electron moves. Figure 2 shows the effects of including the off-diagonal terms into the characteristics in phase space (r, v, μ). Graph A shows the characteristics for the full diffusion tensor, while graph B shows the characteristics when only the diagonal diffusion terms are included (Tang et al. 2020). All characteristics start near the inner boundary. Black lines represent characteristics starting with μ > 0, while red lines with μ < 0. Apparently, by including the off-diagonal terms in the diffusion tensor, the characteristics are first bent toward the central (r,v) plane where μ = 0, which is the new effect caused by off-diagonal diffusion terms. Then the characteristics bend toward the top of the (r,v) plane where μ ≈ 1 again and move to large distances. The first and second terms on the right-hand side give rise to diffusion in velocity space (v, μ), which corresponds to heating and pitch-angle scattering of the electron distribution function. The third and fourth terms are off-diagonal terms that were ignored in Tang et al. (2020). We show in the next section that the off-diagonal terms introduce additional advection in the (v, μ) plane, and ultimately depress the diffusion caused by the diagonal terms.

The Electron Velocity Distribution Functions
The kinetic transport Equation (12) yields the radial evolution of the distribution function f s when approaching steady state. The total electron VDF f t comprises the Maxwellian core f c and the suprathermal electron f s which is the solution of the kinetic transport Equation (12). In our simulations, the fractional constant f r = 2.5 × 10 −5 is used as in case 1 of Paper I. The total VDF f t at 1 au is shown in Figure 3 in which the left graph displays the result of the full diffusion tensor. The suprathermal electron distribution is a combination of a halo and strahl since the strahl has been pitch-angle scattered by Whistler waves. The observed solar wind electron VDF is given by Wang et al. (2012; Figure 5) and by Yoon et al. (2013; Figure 3). Our constructed electron VDF (Figure 3, left) resembles those observed results in the energy range of the Maxwellian core and halo (and strahl). This similarity with the typical form of the observed electron VDF indicates that waveparticle interactions associated with Whistler wave turbulence are able to scatter a cold beam of electrons (strahl) to form the observed distribution.
The right graph of Figure 3 is the electron VDF obtained from the diagonal diffusion tensor presented in Paper I. The two electron VDFs are basically similar but the new distribution does not extend to a velocity as high as ∼10 8 m s −1 . As discussed in Paper I, particle energization is due to velocity diffusion caused by Whistler turbulence scattering of electrons. Evidently, electron energization is less effective when the full diffusion tensor is included. Figure 4 shows the evolution of the electron VDF with heliocentric distance in the two-dimensional phase space (v ∥ , v ⊥ ). The left panel shows results using the full diffusion tensor. The injected suprathermal electron VDF is thoroughly pitch-angle scattered from a cold beam and evolves gradually into an isotropic distribution at about 0.1 au where the scattering of Whistler wave-particle interactions is strong. The suprathermal electron distribution at 0.1 au is almost isotropic and different from the cold beam distribution injected on the inner boundary. This is because as a steady-state solution is approached, the inner boundary distribution of the suprathermal electrons is different from the initial cold beam (t = 0) since scattered and focused electrons can return to their point of origin. The energy and angular distribution of those electrons returning to the Sun (μ < 0 at r = r 0 ) is calculated by the code. At large distance as far as ∼1-2 au where wave-particle scattering is weak, sunward electrons, although the proportion is minute,  persist. In addition, although D vv is the smallest element of the diffusion tensor and has the weakest effect, it nonetheless acts to heat the electron distribution function which can be seen from the suprathermal electrons experiencing diffusion in velocity space, becoming broader in v.
The right panel of Figure 4 shows results from using a diagonal diffusion tensor, i.e., Equation (13) in Paper I. A distinct feature of the full diffusion tensor solutions is that the electron VDFs are cooler (narrower) than those corresponding to the diagonal-only diffusion tensor solution at all heliocentric distances. The off-diagonal terms D vμ and D μv , surprisingly, depress the diffusion of the electron VDF in the (v ∥ , v ⊥ ) plan. Second,the inclusion of the off-diagonal terms leads to fewer antisunward propagating electrons in the right half plane. Thus the effect of pitch-angle scattering is weakened as well. Consequently, D vμ and D μv reduce the pitch-angle scattering and heating compared to the diagonal tensor-only solutions. Figure 5 shows the cuts of a total solar wind electron VDF parallel and perpendicular to the magnetic field. Štverák et al. (2009) presented the radial evolution of the parallel cut of observed electron VDFs in the slow and fast solar wind (see their Figures 6 and 10). The correspondence of Figure 5 with observational results is good in that a similar radial evolution tendency is exhibited. Sunward-propagating electrons are apparent as an enhancement in the v ∥ < 0 region which can also be seen in Figure 4. They apparently do not belong to the Maxwellian core, but were scattered from the antisunward propagating strahl by wave-particle interactions with Whistler turbulence. The emergence of sunward-propagating electrons shows the effectiveness of Whistler scattering of a highly anisotropic field-aligned cold beam into an almost isotropic distribution. The parallel cut of the solar wind electron VDF is asymmetric which means that although the highly anisotropic field-aligned cold beam is scattered into an almost isotropic distribution, there is an imbalance between inward-and outward-propagating electrons. Since the suprathermal electrons streaming in the antisunward direction are still fieldaligned, they can be regarded as part of the strahl.
In Paper I, we showed the radial evolution tendency of the parallel and perpendicular cuts through the solar wind electron VDF for cases with different fractional constants f r in Figures 8  and 11. In the present paper, we use the same fractional constant f r = 2.5 × 10 −5 , as in case 1. By comparing the new graphs with those in Paper I, we find that although the full diffusion tensor results are calculated with the same fractional constant, they resemble more closely case 2 in which the fractional constant is smaller (1 × 10 −5 ). The wave-particle interaction in case 2 is weaker than in case 1, illustrating that the off-diagonal terms D vμ and D μv reduce the effect of waveparticle scattering.

The Plasma Parameters
Once the steady state is approached, we have constructed the (total) solar wind electron VDF f t at different heliocentric distances. By taking moments of the electron VDF, various plasma parameters (for instance, electron number density, temperature, and especially, heat flux), pitch-angle distribution, and anisotropy can be calculated. The appropriate definitions are given in Paper I.
In Figures 6 and 7 black lines show the differential anisotropy ξ and pitch-angle distributions (PADs) of the suprathermal electrons from 0.1 to 3 au for a full diffusion tensor solution. The differential anisotropy increases with increasing electron radial distance and kinetic energy, approaching ∼0.8. It is evident that the suprathermal electron VDF at smaller distances is more isotropic than that at larger distances which reveals that the initially highly field-aligned cold beam at the inner boundary is efficiently scattered by Whistler waves at smaller distances. As heliocentric distance  increases, the scattering effect of Whistler waves is weaker, hence the suprathermal electrons are not scattered sufficiently by Whistler waves to maintain quasi-isotropy, and their VDFs return to anisotropic again.
Red lines in Figures 6 and 7 represent the corresponding data of the diagonal-only tensor model of Paper I. It can be seen that the anisotropy ξ of the full tensor result increases faster than that of the diagonal-only tensor. Therefore, we find that the full tensor results are more anisotropic, which is a consequence of the off-diagonal diffusion terms depressing pitch-angle scattering. Figure 8 shows the radial evolution of the solar wind electron number density and bulk velocity. Here n c represents the Maxwellian core number density in Equation (6). n s,diag is the suprathermal electron number density calculated using the diagonal diffusion terms only, while n s,full uses the full diffusion tensor. At 1 au, both n s,diag and n s,full are about 20 times smaller than n c , consistent with observations. u s is exactly the bulk velocity of the total electrons since the mean velocity of a Maxwellian distribution disappears in the solar wind frame. The heliocentric evolution of the electron number density does not change greatly after including the off-diagonal  tensor terms. However, including the off-diagonal terms in the model reduces the electron bulk velocity.
In Figure 9(A), the black solid line shows the suprathermal electron temperature T s,t as a function of distance, while dotted and dashed-dotted black lines are the parallel and perpendicular temperature that satisfy T s,t = (1/3)(T s,∥ + 2T s,⊥ ). Figures 9B and  9C compare temperatures for the full diffusion tensor case and the diagonal-only case, showing that temperatures for the full diffusion tensor are slightly smaller than that for the diagonalonly case. Although T s,t becomes smaller, it is still higher than that shown in Figure 1(b) of Pierrard et al. (2016). However, the temperature anisotropy A = T ⊥ /T ∥ for the full diffusion tensor case becomes larger than that for the diagonal diffusion-termsonly case. Furthermore, the temperature anisotropy for the full diffusion tensor case exceeds 1 at very small heliocentric distances and decreases to less than 1 with increasing distance.   Scime et al. (2001). Right: comparison of total heat fluxes calculated from the full diffusion tensor model (black line) and the diagonal-only model (red line). The two curves are not offset from one another but illustrate that the amplitudes are different by a factor of about 2.

Solar Wind Electron Heat Flux
Another parameter of great interest both theoretically and observationally is the solar wind electron heat flux q e . In contrast to fluid models, the electron heat flux is not assumed a priori, but calculated everywhere from the third moment of the electron velocity distribution. The calculated suprathermal electron heat flux is the solar wind electron heat flux since the total solar wind VDF is assumed to be composed of a Maxwellian core and the scattered suprathermal distribution, and the heat flux of a Maxwellian distribution is zero. Figure 10 displays the solar wind electron heat flux radial evolution. Scime et al. (2001) show that the solar wind electron heat flux varies as q e = 7.4 × r −2.9 μ Watts/m 2 from Ulysses observations which is shown by the blue dashed line. The calculated heat fluxes of our numerical method decrease as r −2.2 , which is flatter than observed. The right graph in Figure 10 compares the radial evolution of electron heat flux between the diagonal-only model (red line) and full diffusion tensor model (black line). The two heat fluxes follow the same r −2.2 radial scaling but the heat flux from the full diffusion tensor model is smaller by a factor of about 2 than that from the diagonal-only model. Thus the inclusion of the off-diagonal terms just slightly reduces the heat flux.
We check the possible radial dependence of the electron heat flux on the assumed value of the diffusion coefficients since the full diffusion tensor is used and therefore makes this test reliable. We first increase the fractional constant f r from 2.5 × 10 −5 to 7.5 × 10 −5 to model larger diffusion. The results are shown in the left panel of Figure 11. The index of radial evolution for the electron heat flux decreases from r −2 and can be divided into two stages. The radial evolution index decreases to −2.6 within 1 au and achieves a value of −2.1 beyond 1 au. This result suggests some dependence of the electron heat flux radial evolution index on the diffusion coefficients. The larger the value of the diffusion coefficients, the steeper the radial evolution and hence, the smaller the index. Figure 11 (right) shows the comparison of the electron heat flux calculated using different fractional constants. We find that a larger f r produces a larger heat flux intensity. We ascribe this to the stronger D vv that increases the heat flux intensity. On the left-hand side of the transport Equation (12), the second term inside the derivative of velocity v is related to D vv . This term results from the conversion of the transport equation into a conservative form by introducing f = Y/v 2 r 2 , i.e., the equation solved in our numerical code. The positive value of the derivative of the velocity generates a larger electron velocity, as we can see from the characteristics in Figure 2A. This is the motion of the electrons as a whole, yielding electrons with higher velocity. A larger f r (a larger D vv ) results in a more intense motion of this kind for all electrons and hence a larger heat flux intensity. Second, as the fractional constant f r increases from 2.5 × 10 −5 to 10 × 10 −5 , the two heat fluxes have the same radial evolution indices. This suggests that the electron heat flux radial evolution index will not keep decreasing as diffusion is further strengthened, but approaches a final value. The comparison also shows that as the fractional constant is continuously increased, the intensity of the heat flux increases as well and, although it is still smaller than observed, is approaching the observational intensity gradually, which suggests that a stronger diffusion produces a larger heat flux. In order to produce the observed heat flux, the diffusion terms should be suitably strong. Matching the observed intensity of the heat flux is probably not that important since it can be adjusted by increasing the number density at the inner boundary.
Besides using the third moment of the electron VDF to calculate the heat flux, we use a second indirect method to calculate the electron heat flux based on the lower-order moments. This serves as a useful check. The energy flux Q, heat flux q, bulk velocity u, scalar pressure p, and viscosity tensor π are related via (Zank 2014): á ñ . If we ignore the viscous transport of energy π · u and choose a coordinate system such thatê z is along the magnetic field, we obtain the formula for the parallel heat flux (also see Meyer-Vernet 2007), Figure 12 compares the solar wind electron heat fluxes calculated using the two methods, assuming a larger fractional constant f r = 7.5 × 10 −5 . The black line (labeled old) represents Figure 11. Left: the radial evolution of solar wind electron heat flux using the full diffusion tensor with larger values of the fractional constant f r = 7.5 × 10 −5 . Right: the comparison of electron heat fluxes calculated with different fractional constants. The black lines represent the fractional constant f r = 2.5 × 10 −5 , 7.5 × 10 −5 and 10 × 10 −5 , respectively, and the blue dashed line shows the observational result obtained by Scime et al. (2001).
the heat flux calculated from the third-moment method used before, while the red line (labeled new) illustrates the heat flux calculated from Equation (14). The two lines overlap each other showing that the two approaches for calculating the heat flux produce the same results. The convergence of the two approaches indicates the accuracy of the third-order moment approach for the electron heat flux and the numerical method for electron transport.

Summary and Conclusion
We have completed development of the previous numerical method proposed in Paper I for solving the Fokker-Planck kinetic transport equation for solar wind suprathermal electrons interacting with Whistler wave turbulence. Specifically, the diffusion tensor has been extended from a diagonal only to a full general form. The solar wind electron VDF is assumed to be composed of two components: a Maxwellian core and suprathermal electrons which interact with Whistler wave turbulence. By adopting an expanding background of constant velocity radial flow and magnetic field, we obtained the transport equation for suprathermal electrons that describes their spatial and temporal evolution in the solar wind. By solving the equation with the full diffusion tensor, we follow the transport of the suprathermal electrons, and hence the spatial evolution of the solar wind electron VDF. The numerical result demonstrates that wave-particle interactions due to Whistler wave turbulence play an important role in the formation of the solar wind electron VDF and hence the radial evolution of macroparameters, such as the electron heat flux.
A highly field-aligned cold beam of electrons, representing the suprathermal electrons, is injected at the inner boundary. Numerical results are summarized as follows.
1. The full form of the wave-particle interaction tensor associated with Whistler wave turbulence affects the electron pitch-angle scattering and heating, which can be inferred from the full form of the kinetic transport Equation (12). Specifically, D μμ describes pitch-angle scattering and D vv electron heating.
2. The constructed solar wind electron VDFs at 1 au and their radial evolution extended to 4 au resemble observations, which validates the assumption and treatment that the electron VDF is composed of two individual components: the Maxwellian core electron distribution function and the suprathermal electron VDF that is subject to wave-particle interactions throughout the solar wind. 3. The radial evolution of the differential anisotropy and the pitch-angle distribution, together with the electron VDF, shows that the initially highly field-aligned anisotropic suprathermal electrons are significantly pitch-angle scattered and modestly heated by Whistler waves. The effect of Whistler wave turbulence decreases with increasing distance, being unable to scatter the suprathermal electrons sufficiently strongly. The suprathermal electron cannot maintain quasi-isotropy with increasing heliocentric distance and the electron VDF returns to anisotropy far from the Sun. 4. The suprathermal electron temperature slightly increases with increasing distance, which is consistent with observations shown by Pierrard et al. (2016). The temperature of the full diffusion tensor is slightly smaller than that of the diagonal case, but it is still higher than observed. By contrast, the temperature anisotropy of the full diffusion tensor model becomes larger than that for the diagonal diffusion-termsonly case, exceeds 1 at very small heliocentric distances, and decreases to less than 1 with increasing distance. By contrast, the temperature anisotropy of the diagonal-only model is smaller than the unit (see Paper I). 5. A second method to calculate the electron heat flux, besides taking the third moment of the electron VDF directly, using other lower-order moments, is adopted. A comparison between the two methods is shown in Figure 12. The convergence of the two approaches indicates the accuracy of the third-order moment approach for the electron heat flux and the numerical method for electron transport. 6. The off-diagonal terms of the diffusion tensor D μv and D vμ act to depress diffusion in both pitch-angle and velocity space caused by the diagonal diffusion terms D μμ and D vv . Compared with prior simulation results in Paper I, the suprathermal electron temperature and heat flux intensity are slightly smaller, while the number density and bulk velocity scarcely change.
We note that our numerical results yield a strong anisotropy at ∼1 au due to relatively weak scattering of electrons by Whistler wave turbulence, which is contrary to most observations. The diffusion tensor corresponding to Whistler wave turbulence used in the present paper rapidly decreases with distance and f r is a constant. Unfortunately, there is little information available to guide our choice in how the energy density of Whistler wave turbulence changes with increasing heliocentric distance. Thus, it is entirely possible that the diffusion tensor should not decrease as fast, and of course f r is an empirical parameter that attempts to ensure that the scattering of electrons is not overwhelmingly strong based on the theoretical model of Whistler turbulence that seems to overestimate the levels. Second, although electrons in the halo/ strahl energy range appear to best resonate with Whistler waves, it is possible that we should not be neglecting other waves. It is also possible that electrons propagate in a marginal stability regime, thereby introducing the possibility that the Figure 12. Comparison of solar wind electron heat flux using the two calculation methods. The black line (old) represents the heat flux calculated from the third-moment method used before, while the red line (new) illustrates the heat flux calculated from Equation (14). The blue dashed line shows the fitting result obtained by Scime et al. (2001). electron beam alternates between a scattered state and a beam state with the concomitant alternating between self-generation of Whistler turbulence and decay of turbulence (this is the essence of a stochastic growth theory of electron beam propagation advocated by, e.g., Robinson & Cairns 1993). This level of sophistication is beyond the scope of the current paper, however.
The numerical results can be considered to describe the formation and evolution of the solar wind halo and strahl, which is consistent with the consensus that halos are pitchangle scattered from the strahl as electrons propagate in the solar wind. We find that the value of the diffusion coefficients associated with Whistler wave turbulence affects the intensity and radial evolution property of electron heat flux. The observations of heat flux consequently can be used as criteria for determining the required wave-particle interaction of whistler turbulence in theoretical models. The intensity of the calculated heat flux, which is an order of magnitude smaller than observed, can be increased by increasing the halo electron number density at the inner boundary. However, the index describing the radial variation of the heat flux is sensitive to the value of the diffusion tensor. In this sense, observations can place constraints on the index, which provides some insight into physical values of the diffusion tensor and possibly f r .
This work was partially supported by the NSF EPSCoR project OIA-1655280 "Connecting the Plasma Universe to Plasma Technology in AL: The Science and Technology of Low-Temperature Plasma." B.T. acknowledges the support of a NASA Earth and Space Science Fellowship Program grant HELIO19R-0001. G.P.Z. acknowledges the partial support of a NASA Parker Solar Probe contract SV4-84017.