Hybrid Simulation of Solar-Wind-Like Turbulence

We present 2.5D hybrid simulations of the spectral and thermodynamic evolution of an initial state of magnetic field and plasma variables that in many ways represents solar wind fluctuations. In accordance with Helios near-Sun high-speed stream observations, we start with Alfvénic fluctuations along a mean magnetic field in which the fluctuations in the magnitude of the magnetic field are minimized. Since fluctuations in the radial flow speed are the dominant free energy in the observed fluctuations, we include a field-aligned v∥(k⊥)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$v_{\|}(k_{\perp })$\end{document} with an k−1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$k^{ -1}$\end{document} spectrum of velocity fluctuations to drive the turbulent evolution. The flow rapidly distorts the Alfvénic fluctuations, yielding spectra (determined by spacecraft-like cuts) transverse to the field that become comparable to the k∥\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$k_{\|}$\end{document} fluctuations, as in spacecraft observations. The initial near constancy of the magnetic field is lost during the evolution; we show this also takes place observationally. We find some evolution in the anisotropy of the thermal fluctuations, consistent with expectations based on Helios data. We present 2D spectra of the fluctuations, showing the evolution of the power spectrum and cross-helicity. Despite simplifying assumptions, many aspects of simulations and observations agree. The greatly faster evolution in the simulations is at least in part due to the small scales being simulated, but also to the non-equilibrium initial conditions and the relatively low overall Alfvénicity of the initial fluctuations.


Introduction
Not long after the discovery of the solar wind, strong evidence was found for an Alfvénic wave-like character to the observed large amplitude magnetic and velocity fluctuations (Belcher and Davis, 1971), but also spectral evidence for a turbulent cascade in these fluctuations (Coleman, 1968). Subsequent work showed clearly that the fluctuations are much more highly Alfvénic closer to the Sun, and become less so as the spectra of the fluctuations became dominated on progressively larger scales by a fluid-like "Kolmogoroff" regime with a −5/3 slope (for an extensive review, see Bruno and Carbone, 2013). Since the Alfvén wave mode is a normal mode of a magnetohydrodynamic (MHD) plasma, and the spectra seemed fluid-like, it was natural to develop MHD models that attempted to reproduce the evolution. Such models provided plausible explanations for spectral evolution and stream structure of the fluctuations (e.g., Roberts et al., 1991). However, many observations have provided insight into non-fluid, kinetic phenomena associated with both higher time (and thus spatial) resolution and more detailed investigation of the distributions of the particles. While there has been considerable work on nonthermal features of distributions (beams, haloes, temperature anisotropies, etc.; see Marsch, 2006), of most relevance here are the spectral steepenings seen at scales near and smaller than ion kinetic scales (Smith et al., 2006), and the systematic dependence of the temperature anisotropies as a function of the plasma β (the ratio of kinetic to magnetic pressure) (see, e.g., Chen et al., 2016, for recent discussion and references). As illustrated in the references just cited, there is considerable work on idealized simulations of these kinetic effects. The point of this paper is to see to what extent such effects arise naturally out of a generic model of solar wind fluctuations.
Thus, this paper presents a simulation of a simple representation of the solar wind at scales below and into the onset of kinetic effects. The case will be 2D for ease of computation, but with the mean field in the plane of the simulation to allow for Alfvénic wave propagation. Note that prior MHD simulations (e.g., Stribling, Roberts, and Goldstein, 1996) have shown that, at least in that case, 2D and 3D simulations of solar wind turbulence produce similar evolution despite the lack of a degree of freedom. More recent hybrid simulations, (e.g. Franci et al., 2018) also show substantial agreement between 2D and 3D cases. Observational evidence shows a preference for field-aligned wave-vectors for fluctuations close to the Sun (Ruiz et al., 2011), and thus, in our initial state, pure, transverse, slab Alfvén waves propagate parallel to a mean magnetic field taken to be in the x-direction, B 0x , with a −5/3 spectrum somewhat into the kinetic regime; this is reminiscent of large amplitude Alfvén waves at Helios with nearly radial field ( Figure 1, left panel, in which the mean field is nearly radial; see also Figure 2, left for a 2D representation of the transverse magnetic component, B z , at an early stage). The slope of the power spectrum was not taken to be −1 (see Figure 1) to see if the shear flow, introduced to drive the turbulence, could flatten the turbulent fluctuation spectrum. Also, the magnetic field magnitude, |B|, is made nearly constant, as observed (see the very low power in the field magnitude in Figure 1, left), using the procedure of Roberts (2012) that systematically alters the phases of the field components. Note that the assumption of frozen-in flow allows temporal spectra in the solar wind to be interpreted as equivalent to spatial spectra in the simulation.
A number of recent studies are complementary to this one, and provide considerably more detail on some aspects of the solar wind evolution. Turbulence studies in 2D usually have the mean magnetic field out of the plane of the simulation (e.g., Servidio et al., 2014), and thus largely ignore the effects of the restoring force of the mean field. (An exception is the study of Verscharen et al., 2012.) Studies using hybrid or other kinetic methods (e.g. Grošelj et al., 2018;Vasquez, Markovskii, and Chandran, 2014;Cerri, Servidio, and Califano, 2017) provide insights into kinetic effects such as ion anisotropy; they differ from Smoothed spectra from a very low cross-helicity region near 1 AU. High cross-helicity (which declines with distance from the Sun, and is high in pure coronal hole flows) is associated with flatter component spectra and greater field magnitude constancy. The magnetic and velocity component spectral evolution and cross-helicity dependence are presented by Roberts (2010); the correlation with |B| fluctuation level is new here. The data are 40.5 s resolution, from Helios 2, 1976, days 106, 23, and 29 for the respective panels.

Figure 2
The B z component of the field at time steps T = 10 (left) and 40 (right). The initial state has straight contours, and the final state is relatively isotropic. The range of values is about −1.5 to 1.5, where the mean magnetic field is 1.0. the present study primarily in having idealized, typically isotropic initial fluctuations, no net magnetic-velocity correlation (zero cross-helicity), and no limitation on the magnetic field magnitude fluctuations. The present study also takes a more global view, looking in less detail at more aspects of the problem in one study. Sheared flows in the simulation have the velocity, V, parallel to B 0 ; V is initially purely a function of the transverse variable, y, and has a −1 spectral slope. This shear mimics the large free energy available in sheared flows near the Sun (cf. DeForest et al., 2018, Figure 12). The simulation starts with equal and isotropic (particle) ion and (fluid) electron temperatures with a plasma β (ratio of ion thermal pressure to magnetic pressure) of about 1/3. The 2.5D hybrid code is described below. The initial condition is not an equilibrium, so it evolves quickly.
Previous, related, work studied the evolution of an initial spectrum of Alfvénic fluctuations in the solar wind using 1D hybrid codes (e.g., Liewer, Velli, and Goldstein, 2001;Ofman, Gary, and Viñas, 2002;Xie, Ofman, and Viñas, 2004;Maneva, Viñas, and Ofman, 2013), and using 2.5D hybrid codes (e.g., Ofman and Viñas, 2007;Ofman, 2010;Ozak, Ofman, and Viñas, 2015;Maneva, Ofman, and Viñas, 2015). These studies investigated the dissipation of the left hand polarized waves and the resonant anisotropic heating of the protons and alpha particles in solar wind-like plasma. The present study differs from these in its use of sheared velocity fields to mimic solar wind turbulence driving, in the use of highly Alfvénic initial waves, in the analysis methods, and in the nature of questions asked about the results.

Hybrid Simulation Model and Parameters
We employ a 2.5D hybrid model where the protons are treated as particles using the particlein-cell (PIC) method, while the electrons are modeled as a massless neutralizing fluid. Three components of the ion velocities and fields are calculated in the 2D spatial domain. The 2.5D hybrid code is an extension of the 1.5D hybrid code initially developed by Winske and Omidi (1993), generalized to 2.5D and parallelized by Ofman and Viñas (2007). Cartesian coordinates are used, with the x-direction defined along the background uniform magnetic field B 0x . The model uses the following normalizations: the distance unit is the proton inertial length, ρ i = c/ω pp , where ω pp is the proton plasma frequency, the time is in units of the inverse proton gyrofrequency, −1 p , where p = eB 0 /m p c, and the velocities are in units of the Alfvén speed, V A = B/(4πρ) 0.5 , where ρ is the mass density. Further details of the hybrid model equations and normalizations are given in previous studies (e.g., see Winske and Omidi, 1993;Ofman and Viñas, 2007;Ofman, 2010;Ofman, Viñas, and Maneva, 2014;Ofman, Viñas, and Roberts, 2017).
In the present model we use a 256 × 256 computational grid with 512 particles per cell, and we set the plasma β = 1/3. The initial velocity distribution function of the protons is Maxwellian, the initial V x spectrum of fluctuations has a 1/k spectrum, and the transverse magnetic field initially has a k −5/3 spectrum, where k is the normalized wavenumber in units of the inverse ion inertial length.

Results
This section will show a number of ways in which the simulations correspond to reality, and a couple of instructive "misses" as well. Figure 3 summarizes the spectral information in the simulation found by treating cuts through the simulation as time series that would be measured by spacecraft. To obtain good statistics in the spectra, each 256 point line across the box is treated as a member of an ensemble, and all 256 lines are treated as multiple passes through the plasma that can be averaged. All the spectra in Figure 3 are calculated in this way, and given the mean field in the x-direction, it is possible to make very clear parallel and perpendicular spectra using cuts along and transverse to the field (x) direction. The terms "perpendicular" and "parallel" components in the spectra always refer to k-space directions.

Spectral Evolution with Slopes and Breaks Similar to the Solar Wind
Panel (a) of Figure 3 shows the perpendicular dependence of the velocity spectrum in black; it has the k −1 dependence mentioned above. The parallel spectra of magnetic (blue) and velocity (green) fields are shown along with a black −5/3 ideal spectral line. Finally, the parallel spectrum of the magnetic field magnitude has the very low values produced by the construction mentioned above.
Panel (b) shows an early time in the simulation. Already the magnetic magnitude spectrum (purple) is rising toward the component spectra. The cyan trace shows that the perpendicular spectrum of |B|, initially zero, has quickly attained similar levels to the parallel. The shearing of the fields by the flow (see Figure 2) has also generated k-parallel magnetic component fluctuations (red). These fluctuations reflect the flatness of the driving parallel flow spectrum, so there is at least some hint that a possible origin of flat component spectra could be the transverse structure of the flow. On the other hand, the k-parallel magnetic fluctuations (blue) show no sign of flattening.
The parallel magnetic component fluctuations (but not the velocity) show an enhancement at the beginning of the kinetic scales in Panel (b). This enhancement changes some- Figure 4 The evolution of the run in σ c -σ r space. The evolution starts from the highly Alfvénic state represented by the red dots (T = 0) and spreads to progressively lower σ c , reaching the black dot distribution at T = 99.
what, but persists for about 20 time steps. This is similar to the observed peaks in such spectra that are often attributed to the excitation of ion cyclotron waves (Jian et al., 2009). The lack of a peak in the "compressional" spectrum in the observations is reflected in the lack of a bump in the |B| spectrum of the simulation. Future investigation of this peak will check to see if these two phenomena are similar in other ways.
Moving to later times, in Panel (c) the peak is gone and all the component spectra are developing a steeper slope above the kρ i = 1 point. The spectrum follows the dashed line (more visible in Panel (d)) which has a −2.8 slope; this is fairly typical of this range in observations (Smith et al., 2006). The power in |B| slowly continues to approach the component power (it can never reach it), and it takes on the −5/3 spectrum at low k that is observed (Figure 1, above). The final panel shows a general decline of the spectra with a −5/3 slope at low wavenumbers and −2.8 at higher k. It may not be significant here, since such small scales are not wholly reliable in the simulation, but the green and black velocity spectra show enhancements at high k, reminiscent of observations (Roberts, 2010). Although the inertial scale in the simulation does not move to lower wavenumber as it does in the solar wind, the spectral break from inertial range to kinetic scales appears to be moving to lower k in the simulation, simply due to the continuing cascade of power out of the inertial scales.

Evolution in the Cross-Helicity/Residual Energy Plane
'The normalized cross-helicity (or Alfvénicity) of the fluctuations is given by σ c = 2 δb · δv /( (δb) 2 + (δv) 2 ) and the "residual energy" by: σ r = ( δv 2 − δb 2 )/( δb 2 + δv 2 ) where the brackets represent averages, here taken over 50 grid points, although the results are essentially the same with a wide range of averages. The δ's indicate that a running mean has been subtracted out over the averaging interval. Here, v is the wind velocity and b is the magnetic field in Alfvén speed units that make b and v able to be compared directly as in the simulation. Wicks et al. (2013, and the references therein) showed that the evolution in a typical solar wind parcel is from the Alfvénic state with well correlated and equipartitioned velocity and magnetic fluctuations (in red in Figure 4 for the simulation) to a state with a spread similar to that in black dots in Figure 4 but shifted down toward a state of lower kinetic than magnetic energy. This same evolution can be seen in Figure 1 (right) where in the "evolved" state shown, the velocity spectrum is lower than the magnetic. (The Figure 5 A plot showing the initial values of T ⊥ /T vs. β in contours, and the same quantities for the end of the run in gray-scale point densities. These show very similar evolution to that seen in the inner heliosphere by Matteini et al. (2007). The line shows the empirical fit to Helios data for high-speed wind by Marsch, Ao, and Tu (2004). The relationship of this line to linear instability boundaries is shown in the paper of Hellinger et al. (2006). continuation of this evolution is shown by Roberts (2010); the velocity remains below the field.) Thus the simulation captures only part of the dynamics. One missing ingredient could be the expansion of the plasma, since different conservation laws govern the kinetic and magnetic evolution in that case. However, previous studies have shown a residual energy evolution away from zero, and often to the observed negative values, when either the mean field is out of the 2D simulation plane (Papini et al., 2019;Franci et al., 2015) or when the mean field in the plane goes to zero (Roberts et al., 1992). Thus at least part of the problem here is that the one-dimensional k-space perpendicular to the mean field is not sufficient to capture the residual energy evolution. Figure 5 shows the proton temperature anisotropy versus the parallel plasma beta measured at the beginning (contours) and end (grayscale) of the simulation. The behavior of the solar wind in this space has been studied by many authors in the context of the possible limiting of the evolution of the states by instability (e.g., firehose, mirror) thresholds (see Chen et al., 2016, and the references therein). Note that there is no expansion in the simulation, and the time scale is much faster than in the solar wind, but what we find is very like what is seen in the inner heliosphere (Matteini et al., 2007). The plot shows the same trend as the empirical fit to Helios data shown, but it is offset; our initial condition did not match the fit line, and this should be changed in a subsequent study. The results in this study are consistent with linear instability boundaries; see the plots and discussion by Hellinger et al. (2006).

Kinetic Energy, Magnetic Energy, and Cross-Helicity 2D Evolution
in k-Space Figure 6 shows 2D power spectra of the fluctuations. Initially the magnetic field only contains parallel propagating Alfvén waves, whereas the velocity also contains the shear modes. Spectrally, shear couples to the wave modes, producing modes that fill in the k-space, eventually becoming nearly isotropic. There is little change after T = 60. While there is a very slight perp dominance (reminiscent of the "evolved" solar wind cf. Ruiz et al., 2011), the effect here is not as dramatic as in some other cases (Roberts et al., 1992;Ghosh et al., 1998). The initial wave-like spectrum is just that of the parallel propagating Alfvén waves, and thus all the cross-helicity is concentrated where the wave modes are defined (see Figure 7). The correlated region spreads, and the shear modes introduce an opposite correlation (unclear why; shear should contribute zero, on average). The memory of the initial geometry persists for some time, but eventually the distribution is nearly isotropic in σ c as well as in power. The ideal MHD conservation of cross-helicity may be responsible for the slow decline of the correlation (still dominated by the initial values) through the evolution.

Conclusions
In many respects the simulation presented here follows the solar wind evolution, although on a much faster time scale. The initial conditions are chosen to match various observational constraints: flat velocity spectrum, low |B| fluctuations; and slab-like and highly Alfvénic magnetic and associated velocity fluctuations. Spectra, including that of the magnitude of the magnetic field, evolve to very close to the observed slopes, and the thermal evolution is similar as well. The fluctuations in the magnitude of the magnetic field increase as they do in observations. We find a (transient) enhancement in the magnetic field spectrum near the proton inertial length, again as observed. There is (weak) evidence in the simulations that flat magnetic spectra could be generated by velocity shear with similar spectra. The slow decline in cross-helicity is also consistent with the evolution of highly-Alfvénic regions in the solar wind (Roberts, 2010). The disparate time scales and the equipartition of magnetic and velocity energies have testable explanations in terms of more Alfvénic and larger-scale initial states, and expansion during the evolution. Thus the simple state assumed here provides an interesting first approximation to the conditions early in the turbulent evolution of solar wind fluctuations. Observations from Parker Solar Probe and Solar Orbiter should  , 5, 40, and 99 (top left, right; bottom left, right); σ c = −1 (red); σ c = 1 (purple). make the considerations here more precise, and provide additional, more stringent, tests of the ideas. Guided by this exploratory simulation, future simulations should have higher resolution, three dimensions in space, higher cross-helicity for the shear component, and a thermal anisotropy consistent with observations. All these conditions are possible with current computing resources, and should get us closer to a full representation of solar wind turbulence.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.