Magnetohydrodynamic Turbulence in the Earth’s Magnetotail From Observations and Global MHD Simulations

Magnetohydrodynamic (MHD) turbulent flows are found in the solar wind, the magnetosheath and the magnetotail plasma sheet. In this paper, we review both observational and theoretical evidence for turbulent flow in the magnetotail. MHD simulations of the global magnetosphere for southward interplanetary magnetic field (IMF) exhibit nested vortices in the earthward outflow from magnetic reconnection that are consistent with turbulence. Similar simulations for northward IMF also exhibit enhanced vorticity consistent with turbulence. These result from Kelvin-Helmholtz (KH) instabilities. However, the turbulent flows association with reconnection fill much of the magnetotail while the turbulent flows associated with the KH instability are limited to a smaller region near the magnetopause. Analyzing turbulent flows in the magnetotail is difficult because of the limited extent of the tail and because the flows there are usually sub-magnetosonic. Observational analysis of turbulent flows in the magnetotail usually assume that the Taylor frozen-in-flow hypothesis is valid and compare power spectral density vs. frequency with spectral indices derived for fluid turbulence by Kolmogorov in 1941. Global simulations carried out for actual magnetospheric substorms in the tail enable the results of the simulations to be compared directly with observed power spectra. The agreement between the two techniques provides confidence that the plasma sheet plasma is actually turbulent. The MHD results also allow us to calculate the power vs. wave number; results that also support the idea that the tail is turbulent.


INTRODUCTION
One of the primary goals of current investigations in space physics is to understand how electromagnetic energy stored in the magnetotail is transferred to plasma energy. Turbulence is a multi-scale phenomenon that mediates the transport of energy, mass, and momentum. Unsteady, but nonrandom fluctuations in the magnetic and electric fields and flows characterize turbulence (e.g., Karimabadi et al., 2013). Turbulent spectra have been observed in many space plasmas. For example, turbulent fluctuation spectra have been found in the solar wind (Coleman, 1968;Matthaeus and Goldstein 1982;Roberts et al., 1987a;Roberts et al., 1987b;Tu et al., 1989;Marsch and Tu, 1990;Sahraoui et al., 2010), in the magnetosheath (Zimbardo et al., 2010;Li et al., 2020) and in the magnetosphere (Borovsky et al., 1997;Lui 2001;Weygand et al., 2006;Weygandf et al., 2007). In the magnetosphere, turbulence exists over a wide range of scales from large scale magnetohydrodynamic (MHD) flows to kinetic dissipation scales (see Zimbardo et al., (2010) for a review of turbulence studies in the geospace environment including both the magnetosheath and the magnetosphere). Emphasis in this paper is on our understanding of the MHD turbulence, its consequences for transport and dynamics, and, specifically, its relationship to magnetic reconnection, with a focus on the magnetotail. In particular, several studies have presented evidence that turbulence in the plasma sheet is an important mechanism for energizing plasma in the magnetotail (Borovsky et al., 1997;Angelopoulos et al., 1999;Chang, 1999;Klimas, 2000;Borovsky and Funsten 2003;Weygand et al., 2007).
In general, turbulence in plasmas can be thought of as resulting from oscillations in velocity or magnetic field driven by nonlinear processes at large scales (Kadomtsev, 1965) and as vorticity in fluid motion where the inertial forces in the vortices are larger than the forces that are damping the eddies (Leung and Gibson, 2004). The turbulence in any medium transfers energy from the largest scales to small dissipation scales, but, in some circumstances, can involve an inverse transfer from small scales to large scales (Frisch and Kolmogorov, 2001). Observations of turbulent fluctuations in the solar wind have been discussed for decades. For example, one of the primary goals of the Parker Solar Probe mission (Nature 2019) is to find the source of the Alfvénic turbulence in the solar corona. Compared to the solar wind, there are relatively few observational studies of turbulence in the magnetosphere. Such studies as there are include observations of the fluctuations in the magnetic field, plasma, and electric field measurements associated with power in both the inertial range and the dissipation range. A schematic of this process is shown in Figure 1. The large scales, or energy containing scales, will drive turbulence to an inertial range that, in fluid turbulence, was shown by Kolmogorov (1941) to have a power spectral index as a function of wave number that is -5/3. In this range of scales, turbulent eddies will break-down until dissipation causes the spectrum to steepen as viscosity damps the eddies. The beginning of the dissipative range is the Taylor scale (Kolmogorov 1941) scale and it is at this (kinetic scale) that heating occurs as the fluctuations are damped. In space, one generally measures time series, so that the power spectra are usually shown as a function of frequency but the validity of such a representation depends on the nature of the background flow and the properties of the fluctuations and often is determined by the validity of the Taylor frozen-in-flow hypothesis (Taylor, 1938), which is discussed in more detail below.
In this review, we first consider in Observations of Turbulence in the Magnetotail the observational evidence that turbulence exists in the tail. The approach used for studies of solar wind turbulence is not directly applicable to the magnetotail because contrary to the case in the solar wind, boundary effects are often important and the fluctuations may not satisfy conditions required of a random stationary process (Matthaeus and Goldstein, 1982;Perri and Balogh, 2010) which justifies construction of power spectra. In the super-Alfvénic solar wind flow, one can usually assume the validity of the Taylor frozen-in-flow hypothesis (Taylor, 1938) which justifies the transformation from frequency spectra to wave number. The Taylor hypothesis states that frequency ω relates linearly to the wave vector k (i.e., ω k v, where v is the plasma velocity) when the magnetic field evolves on a timescale longer than the time it takes it to flow past the spacecraft. In the magnetotail, the flow rarely exceeds the Alfvén speed and thus time and spatial scales are difficult to separate and determining the power spectrum of fluctuations as a function of wavenumber is challenging. Therefore, it is difficult to determine unambiguously the spectra of observed fluctuations. The flow speeds in the magnetosheath are usually higher and use of the Taylor hypothesis is frequently valid.
Consequently, we include a discussion of a variety of analysis methods in our discussion of fluctuations in the tail. In addition, analysis techniques used for studies of magnetic reconnection in the magnetotail can be applied to analyze the nature of the observed flows.
A number of generation mechanisms including flow shear instabilities such as the Kelvin-Helmholtz instability and flows from reconnection have been proposed (Matthaeus et al., 1984;Montgomery 1987;Angelopoulos et al., 1993;Borovsky et al., 1997;Lui 2001;Antonova and Ovchinnikov, 2002;Vörös et al., 2003). In Magnetohydrodynamic Simulations of Turbulence in the Magnetotail, we review work on evaluating the mechanisms by using MHD simulations. Finally, in Some Unsolved Questions About Turbulence That can be Addressed Using Modeling, we FIGURE 1 | Schematic of a typical power spectral density plot for turbulent solar wind. The lowest frequencies are at the energy containing scales and the dissipation range is at the highest frequencies. The correlation scale (left most dashed gray line) separates energy containing scales and the inertial range while the Kolmogorov scale (right most dashed gray line) is between the inertial range and the dissipation range. The third scale on the schematic is the Taylor scale (middle dashed gray line), the scale at which eddies start to damp out, (Adapted from Goldstein et al. (1995)).
consider several unsolved problems associated with turbulence in the tail and discuss future simulations to address them.

OBSERVATIONS OF TURBULENCE IN THE MAGNETOTAIL
A simplistic method suggesting the presence of turbulent fluctuations in a time series is to examine the spectral index in the inertial range of the power spectra. Kolmogorov (1941) used dimensional analysis to argue that the spectral index for fully developed fluid turbulence should be -5/3 for spectra of power vs. wavenumber. Perhaps surprisingly, this value is frequently observed in the fluctuating magnetic fields of the solar wind and magnetosheath (also see, Podesta et al. (2007)). On the other hand, Kraichnan (1965) found a value of -3/2 for ideal isotropic incompressible MHD turbulence. The difference between the neutral fluid values and magnetized fluid value comes from the number of degrees of freedom in the fluid. In both derivations, the rate of energy transfer from the driving scale of the spectra to the dissipation range of the spectrum was held constant. Whether or not the energy transfer rate is constant is important for differentiating Kraichnan and Kolmogorov type turbulence from intermittent turbulence. In intermittent turbulence, the energy transfer rate may not be constant and the turbulence may not yet be fully developed.
Within the plasma sheet, Borovsky et al. (1997) used ISEE-2 data and found for the plasma flow velocities that the slope of the power spectral density vs. frequency was -0.8 to -2.0 while that for the magnetic fieldit was between −1.6 and −3.0. Borovsky et al. (1997) used a "random sweeping model to approximate the conversion from wavenumber to frequency. (Even though these values include the theoretical value they do not confirm the presence of turbulent flows in the plasma sheet). Several different phenomena can explain these values including waves and/or driving phenomena, or it could be that the time series is neither stationary nor fully developed. Borovsky and Funsten (2003) suggested that this range of spectral indices could result from boundary effects or a combination of driving mechanisms each with different spectral indices. Weygand et al. (2005) found spectral indices in Cluster plasma sheet magnetic field data that were closer to −2 for the inertial range but did not take into account the speed of the flows. Vӧrӧs et al. (2004), also using Cluster magnetic field data, obtained a value of −2.6, but it is not clear if this value applies to the inertial range, dissipation range, or somewhere in between. Ergun et al. (2014) used MMS data and found a clear spectral index of −5/3 within the magnetic field inertial range, but a spectral index closer −3/2 in the electric field inertial range data. Chaston et al. (2012) found a similar value within the inertial range for electric field power spectra using THEMIS plasma sheet electric field data. Overall, the studies show that for slow speed flow the Taylor hypothesis is not valid but may be valid for high speed flows like those associated with magnetic reconnection. More work on this is needed. Kozak et al. (2018) used the magnetic field measurements from four spacecraft of the Cluster mission for the analysis of turbulent processes in Earth's magnetotail. They obtained power-law scaling of the generalized diffusion coefficient indicating the presence of super-diffusion processes. Prior to the dipolarization, the spectral index was in the range between −1.68 and −2.08 while in the dipolarization on larger time scales the index was between −2.2 and −1.53. However, when they examined the data in the dipolarization on smaller time scales they found that the range was −2.89 to −2.35. They also report that the break in the spectra occurs at approximately the average proton gyrofrequency.
The wide range of spectral indices from the magnetic field and electric field data suggests the presence of intermittent turbulence within the plasma sheet. One method of determining the presence of intermittent turbulence is to look for non-self-similar scaling of the fluctuation probability distribution function. A number of studies have demonstrated that non-self-similar scaling of probability distribution functions exists within the plasma sheet Weygand et al., 2006;Stawarz et al., 2015].
To determine if there is non-self-similar scaling, one constructs probability distribution functions of the fluctuations across a range of times and calculates the kurtosis (i.e., fourth moment) for each distribution. If the kurtosis systematically decreases with increasingly temporal separation in the time series, then most likely the turbulence is intermittent (Weygand et al., 2006). All three of these studies found that magnetic field fluctuations observed with Cluster and THEMIS exhibited non-self-similar scaling of probability distributions during geomagnetic active periods within the plasma sheet. Both Stawarz et al. (2015) and Weygand et al. (2005) used single spacecraft observations to show non-self-similar scaling, which requires that the Taylor hypothesis applies, and suggests that the fluctuations are evolving slowly with respect to the time required for the plasma to flow past the spacecraft and are consequently frozen in the flow. Weygand et al. (2006) took this method one step further and avoided assuming the Taylor hypothesis by using pairs of Cluster spacecraft separated in space to show that non-self-similar scaling is present within the plasma sheet magnetic field. This non-self-similar scaling demonstrated the presence of intermittent turbulence in the plasma sheet. Thus much of the prior work suggests that turbulence is present in the plasma sheet, but that it may not be fully developed or may be intermittent.
Assuming that turbulence is present in the plasma sheet, we can determine the three fundamental associated scale lengths, viz., the correlation scale, the Taylor scale, and the Kolmogorov scale (Goldstein et al., 1995;Weygand et al., 2007). The correlation scale is the energy-containing scale and the scale at which the inertial range turbulent cascade begins to exhibit a power spectrum. The Taylor scale is the scale at which the turbulent eddies begin to damp out, and the Kolmogorov scale is the point at which dissipation begins ( Figure 1). A series of studies by Weygand et al., 2007;Weygand et al., 2009;Weygand et al., 2010) used Cluster spacecraft pairs over many years combining many different intervals to produce twodimensional cross correlation maps of the magnetic field fluctuations within the plasma sheet. An example of such a map is shown in Figure 2. From these maps, we can determine correlation scales and Taylor scales. The correlation scale is the 1/e folding distance on the correlation vs. spacecraft separation curve and the Taylor scale is the radius of curvature of the two point cross correlation function of the magnetic field fluctuations at the origin Weygand et al., 2007). Weygand et al., 2007;Weygand et al., 2009;Weygand et al., 2010) found that the correlation scale and Taylor scale varied with their angle with respect to the mean magnetic field direction for geomagnetically quiet and moderate conditions. The correlation scale and Taylor scale were longest along the mean magnetic field (about 16,400 km and 2900 km, respectively) and shortest perpendicular to the field (9200 km and 1100 km, respectively). In the plasma sheet, the resulting twodimensional correlation maps were similar to the twodimensional correlation maps of quasi two-dimensional turbulence observed in the slow solar wind Weygand et al., 2011), suggesting the presence of quasi twodimensional turbulence within the plasma sheet. Physically, the correlation scale is associated with the approximate thickness of the plasma sheet and/or driving scales and the Taylor scale is expected to be about the same size or larger than an ion gyroradius in the plasma sheet (∼700 km). Weygand et al. (2009) suggested that the difference in the Taylor scale size with magnetic field direction ( Figure 2) might be related to dispersive and dissipative effects. Similar differences were noted in the solar wind in the same paper.
The correlation scale (λ CS ) and the Taylor scale (λ TS ) derived from the two-dimensional correlation maps can be employed to determine an effective magnetic Reynolds number (R eff (λ CS / λ TS ) 2 ) for the plasma sheet. The effective magnetic Reynolds number (i.e., Lundqvist number) is an important parameter to help validate numerical MHD models and it suggests where dissipation scale or, for Hall effects, dispersion become important. Weygand et al. (2009 obtained effective magnetic Reynolds numbers between 10 and 110 for the plasma sheet. These values are similar to Lundqvist numbers for magnetospheric MHD models, but atre significantly smaller than 1600 reported in Vӧrӧs et al. (2006). El-Alaoui et al., (2010) estimated the magnetic Reynolds number from an MHD simulation to be between 100 and 1000 except near sites of reconnection where the Reynolds number is less than 10.
The statistical studies used to investigate observations of turbulence can also be applied to simulations, which enables us to compare more directly observations that indicate turbulence with the simulations. The focus of this review is on global MHD calculations of turbulence in the magnetotail but we should mention additional modeling studies using other techniques. Several calculations have modeled particle motion in turbulent fields (e.g. Taktakishvili et al. (2001), Greco et al., (2002), Zimbardo et al. (2003) We discuss the MHD simulation results in more detail in Magnetohydrodynamic Simulations of Turbulence in the Magnetotail.

MAGNETOHYDRODYNAMIC SIMULATIONS OF TURBULENCE IN THE MAGNETOTAIL
The MHD equations written in primitive form illustrate the nonlinear terms responsible for MHD turbulence.
FIGURE 2 | Correlation contour plot for the plasma sheet. Correlations calculated for one quadrant were mirrored into the other quadrants. Rperp is the separation perpendicular to the mean magnetic field direction and Rpara is the separation parallel to the magnetic field direction. The color bar gives the value of the cross correlation coefficient (Weygand et al., 2009 where ρ is the mass density, v is the plasma flow velocity, B is the magnetic field vector, P is the thermal pressure, µ o is the permeability of free space, and η is the resistivity. In the magnetotail at the largest scales, well described by MHD, flows driven on the scale of the entire system, break up into structures that cascade to smaller scales in an (energy cascade) . Borovsky and Funsten (2003) argue that plasma sheet turbulence is due to vortices, or eddies as shown in two dimensions in Matthaeus and Montgomery (1980). In the plasma sheet, turbulence produces intense mixing (Antonova and Ovchinnikov, 1999;Antonova and Ovchinnikov, 2002). A number of local (e.g. Nykyri and Otto, 2001;Nykyri et al., 2006b) and global MHD simulation studies have shown vortices forming at the magnetopause (e.g., Hwang et al., 2011;El-Alaoui et al., 2012;Hwang et al., 2012;Sorathia et al., 2019) and in the tail (e.g., Ashour-Abdalla et al., 1999;Fairfield et al., 2000;White et al., 2001;Ashour-Abdalla et al., 2002;Slinker et al., 2003;Hasegawa et al., 2004;Walker et al., 2006;Collado-Vega et al., 2007;Claudepierre et al., 2008;Hwang et al., 2011;Sorathia et al., 2019). However, the expected behavior in three dimensions has yet to be studied in any depth. Vortices observed at the magnetopause occur for both northward and southward interplanetary magnetic field (IMF) and have been interpreted as nonlinear Kelvin-Helmholtz waves (e.g. Hwang et al., 2011;El-Alaoui et al., 2012;Hwang et al., 2012;Sorathia et al., 2019). The vortices reported in the tail were, in general, not associated with boundary oscillations. For example, Ashour-Abdalla et al. (2002), Walker et al. (2006) found large-scale vortices in the central plasma sheet during simulations of prolonged intervals with southward IMF.
The results from idealized global MHD simulations of the solar wind-magnetosphere ionosphere system driven by simplified solar wind and IMF conditions (constant solar wind with either southward or northward IMF), exhibit field fluctuations with spectral properties similar to the observations (El-Alaoui et al., 2012). An interesting feature revealed by these simulations is that the fluctuation energy is transported to small regions of high dissipation as described by, for example Wan et al., (2012). However, MHD cannot reveal the processes causing the dissipation. The statistical properties of the observed fluctuations indicate that localized regions of highdissipation are formed (see, for example, a study carried out in the solar wind by Greco et al., (2009)). The simulations show that these fluctuations are associated with reconnection (El-Alaoui et al., 2009;El-Alaoui et al., 2010). However, it is unknown whether large-scale turbulence enhances, diminishes, breaks up, or otherwise affects the microprocesses involved in magnetic reconnection. Investigating this theoretically requires particle-in-cell (PIC) simulations. However, local PIC simulations of the reconnection region do not include the substantial energy input from large-scale turbulence. It is thus important to include the large scale driving as well as the microscopic dissipation in the same calculation as was done, for example, in Wu et al., (2013).
We have recently introduced a method to couple the large scale drivers to the local kinetic scales. In this approach (Walker et al., 2019), the MHD results provide the initial state and the driving boundary conditions for a particle-in-cell simulation of a substantial portion of the magnetosphere. The use of the implicit moment method as implemented in the iPic3D code (Markidis et al., 2010) allows one to consider larger domain sizes within the kinetic approach. The approach has been shown to introduce correctly the physics of reconnection at electron scales, replicating the presence of electron crescent-shaped distributions (Lapenta et al., 2017) and of unsteady reconnection processes feeding further into the turbulence cascade (Ashour-Abdalla et al., 2016) and resulting in electron (Ashour-Abdalla et al., 2015) and ion  energization. The kinetic level of description also reveals instabilities both near the reconnection site (Lapenta, 2008;Walker et al., 2018) and in the outflow (Divin et al., 2015;Lapenta et al., 2015) that further drive (Pucci et al.,2017;Lapenta et al., 2018) and impact energy exchanges Lapenta et al., 2020a). For instance, Lapenta et al. (2020b) argue that turbulent acceleration is responsible for the formation of power law tails in the distribution functions of energetic electrons.

Turbulence for idealized conditions
El-Alaoui et al., (2010) examined results from a 3D MHD simulation of the magnetosphere with nominal solar wind parameters and a southward IMF. They demonstrated that flows in the plasma sheet were consistent with turbulence. This global MHD simulation required very small grid spacing (<0.13 R E ) to resolve flow vortices and turbulence in the plasma sheet. The simulation was run with southward IMF (5 nT) for 240 min followed by a northward IMF of the same magnitude. They found that fluctuations and eddies occurred under these steady driving conditions. In the tail, the regions of high grid resolution were large, including the near earth and midmagnetotail regions. Snapshots taken during the southward IMF interval are in Figure 3. Three variables are superimposed in areas on the equatorial plane. The color contours display the B Z component of the magnetic field, the black arrows show flows in the plane and the white isocontours give the locations of the last closed field lines. Localized reconnection can be identified by flow reversals, from earthward to tailward, associated with reversals in the B Z component of the magnetic field. The B Z component was complex during this interval with filament like structure where it was large and positive at several locations. A spacecraft encountering this would see dipolarization front like signatures in the magnetic field. Vortices are apparent earthward of the reconnection (panel a). Nested within the larger vortices are smaller vortices (panel b) and within the smaller vortices we find even small vortices (panel c). This pattern of nested vortices exists throughout the magnetotail earthward of the reconnection sites throughout the simulation. Figure 4 shows the power spectral density (PSD) in the tail ( Figure 4A) with the Kolmogorov [1941] spectrum superimposed. The slope changes to -3 at ∼30 mHz, which we interpret as the dissipation regime. A qualitatively similar change in slope is also observed in the solar wind (e.g., Alexandrova,  2008; Chen et al., 2010;Sahraoui et al., 2010) and in the magnetosphere (e.g., Vörös et al., 2005;Nykyri et al., 2006a;Matthaeus et al., 2008).
We computed inertial and dissipative range PSD slopes at forty-nine locations in the equatorial magnetotail. All of the PSDs at these locations exhibited the three frequency ranges defined above: the driving (energy containing) range, the inertial range, and the dissipative range. We constructed histograms of the spectral slopes for the inertial ranges ( Figure 4B) and for dissipative ranges ( Figure 4C) for southward IMF. We found that the PSDs have a median value of −1.77 while in the dissipative range the median value is −3.9 although the distribution is broad. In a statistical distribution of observed spectra Weygand et al. (2005) found a peak in the distribution of slopes at a spectral index of −2.0 and a weak secondary peak at a slope of −1.6. There is a trend toward more negative slopes for more tailward locations.
MHD simulations require a source of dissipation for reconnection to occur. One way to provide dissipation is to add an extra term of the form ηJ where η is a resistivity and J is the current density in Ohm's law (El-Alaoui, 2001;El-Alaoui et al., 2009). Dissipation in the simulation contributes to plasma sheet turbulence in two ways. On a large scale, it leads to reconnection that drives the turbulent flows and on small scales it dissipates the energy. It is important to have this term in the MHD code. For instance, without this term reconnection does not occur in the tail .
While the turbulent flows in these calculations are related to reconnection, turbulent flows in the magnetotail associated with the Kelvin-Helmholtz (KH) instability also have been reported. Observations of changes in the magnetopause position consistent with KH have been discussed by numerous authors (e.g., Sckopke et al., 1981;Song et al., 1988;Chen et al., 1993;Fairfield et al., 2000). Intervals with northward IMF provide an opportunity to investigate boundary oscillations in the absence of strong flows from plasma sheet reconnection. For instance, several studies have reported vortices forming along the flank boundaries from MHD simulations (Li et al., 2009;El-Alaoui et al., 2010;Sorathia et al., 2019). They indicated that the KH instability was likely the source but noted that reconnection also was occurring at high latitudes in the polar cusp region (Hwang et al., 2012). Li et al. (2009) argued that the combined processes form a cold dense boundary layer. All of the papers indicated that the process was turbulent. A typical example of the resulting flows based on the El-Alaoui et al. (2010) simulation is shown in Figure 5. Power spectral densities at different locations for the northward IMF interval are shown in Figure 6. The PSD power law indices in the inertial range had more variation than in the southward case but were still consistent with turbulence. The PSD power law indices in the dissipative range also varied widely but were generally more negative than in the inertial range.
The pattern of vorticity in the tail was much more extensive for the southward IMF. This suggests that reconnection driven flows are more important than the KH instability for driving the tail into a turbulent state. However, more work needs to be done to quantify the relative contributions of reconnection and KH.
Our ability to correctly resolve the turbulent cascade resulted in large part from recent improvements in the resolution of MHD simulations [e.g., Guild et al., 2008;El-Alaoui et al., 2009). The success in simulating the overall form of turbulent spectrum (Figure 4,6) gives us confidence That the details of dissipation may not significantly affect that overall turbulent spectrum. El-Alaoui et al., (2013) investigated the properties of fluctuations during a February 7, 2009 substorm in which the WIND spacecraft provided solar wind data. The simulation of this substorm suggests that the configuration of the tail and its evolution is very complex. The simulation results and THEMIS observations in the magnetotail were remarkably similar to dipolarizations and strong flows (see Figures  3,4 of El-Alaoui et al., (2013) for a detailed comparison between the observed magnetic field and flows and the simulation). As has been reported in previous substorm simulation studies (e.g. Ashour-Abdalla et al., (2015)) a large dipolarization grows by accreting smaller earthward-moving dipolarization fronts (DFs). In this case the dipolarizations were associated with a strong channel of earthward flow that formed a large vortex near earth. Figure 7 shows two snapshots of the flows, Bz and the thermal pressure in the plasma sheet (maximum pressure surface) at two times. The plots on the left are from 0359UT and those on the right are from 0405UT. Dipolarizations from reconnection at about X −15RE have moved into the inner tail. The flows in these channels, combined with lower velocity return flows, form several large vortices in the plasma sheet. The panels in Figure 8 show successive blowups of Bz and the flows at 407UT. Within the larger vortices, we find smaller nested vortices like those found in the generic simulation described above, although the vortex structure is more complex in the event driven case (compare Figure 3 with Figures 7,8). These vortices form the low frequency end of the turbulent cascade.

MHD turbulence in simulated substorms
The MHD simulation and the observations show that the magnetotail mesoscale structures (e.g., dipolarizations, flow channels, localized reconnection sites, and flow vortices) during a substorms are in a state of rapid change. In Figure 9, we have compared the spectra of Bz from times between 0200 UT and 0500 UT, at z −3 R E and y −2 R E for three x values (−10 R E , −12 R E , and −14 R E ) with observations at THEMIS P5. The driving, inertial and dissipative scales appear in both the observations and the simulation. The slope of the spectrum in the MHD results at high frequencies becomes nearly flat as the MHD time-step size (typically between 0.01 and 0.1 s) is approached. The slopes in the upper right corner of the figure come from least squares fits to the results. The simulation (black) power levels in the driving and inertial ranges frequently overlap those observed (red). The spectral indices in the inertial ranges are similar, with an observed index at P5 of −1.74 and values from −1.85 to −1.80 from the simulation. The frequency between the inertial range and the dissipative range appears to be about 30 mHz. This feature is similar to that found in the generic cases. The numerical experiments with generic simulations have shown that the overall level of dissipation in the simulation controls the location of this breakpoint (El-Alaoui et al., 2013) but does not affect the driving and inertial ranges. Although the resistivity in the MHD model gives the observed frequency where the slope changes, it does not tell us which physical mechanisms dissipate energy at the high-frequencies.
Probability distribution functions (PDFs) enable us to characterize the turbulent flows. In particular, the fourth moment of the PDF (the kurtosis) provides a measure of non-Gaussian nature of the wings of the PDF. If the kurtosis of the PDF decreases with increasing lag (τ) the turbulence is said to be intermittent. In plasma flows, intermittent turbulence occurs if dissipation is localized to specific regions of space. Figure 10 shows the PDFs from this simulation as black curves and PDFs constructed from THEMES P5 observations in red using three lags (10, 200 and 900s). The black dotted line gives values from a Gaussian. The numbers at the top of the panels give the kurtosis. A Gaussian distribution has a kurtosis of three while the larger values reflect energy in the wings of the distributions. The kurtosis tends to decrease for increasing lags consistent with the expectations of intermittency. Similar results in the solar wind were found in Wu et al. (2013).
The vorticity in both the generic case ( Figure 3) and the event simulation (Figure 7) suggest a strong interaction between the reconnection generated flows and the tail becoming turbulent. As discussed in Observations of Turbulence in the Magnetotail, the Taylor hypothesis may not be applicable to the tail under all flow conditions. One advantage of simulations is that we can calculate the power in k space. We computed spatial power spectra in the region between −10 and −22R E in x and −4 to 8R E in y at the maximum pressure surface . The fluctuations in this box were transformed to k space by using a Fourier transform. This yields the power at the allowed wavenumbers in the square domain. These are given by k x 2πn/L and k y 2πm/L where L is the length of the box on a side and n and m are integers 0, 1, 2, 3, . . .). Most of the power is at low wavenumbers as expected ( Figure 11, bottom), since the lowest wavenumbers are dominated by large scale features of the magnetotail including the dipole field and the effects of the box's edges.
To display the variations in turbulent power as a function of time, we summed over n and m, excluding the lowest wavenumbers (0, 1, and 2). The mesoscale features with sizes less than or equal to 3 R E are included, however. We then computed the sum of the power in the remaining modes between 0340 and 0412 UT. Because dipolarizations are characterized by B z increases and the flow is primarily earthward (V x ), these two parameters are used in the analysis. The plot on bottom left of the figure is a time before a dipolarization front forms while that on the right is after formation. The total power is much higher at the second time (0353 UT). These spectra show that the additional power extends across the k-spectrum, showing that the turbulent power quickly cascades across the wavenumber range. The top panels of Figure 11 give the summed power vs. time. The maxima in the Bz and Vx power correspond with the evolution of the mesoscale structures in the magnetotail. The largest feature is a major increase in Vx power starting at 0349 UT and peaking at 0353 UT. An increase in B z power starts at about the same time and peaks a few minutes later. The start of this increase corresponds to the formation of a narrow, intense flow There is considerable work from both the data analysis and the simulations that support the conclusion that localized reconnection-driven earthward flows can generate turbulent fluctuations in the tail. Observations of localized reconnection regions with strong high-speed outflow support the idea that reconnection is an important process driving turbulence in the plasma sheet (Vörös et al., 2006;Angelopoulos et al., 1999). Huang et al., (2012) and Osman et al., (2015) have analyzed the magnetic field fluctuations observed by Cluster during a period of strong earthward plasma sheet flow (∼1200 km/s). The magnetic field fluctuations are consistent with turbulence. Localized dissipation may drive these fast flows by enabling localized reconnection. The resulting outflow jets can initiate a turbulent cascade. For example, when a high-speed flow such as a bursty bulk flow (BBF) reaches the near-Earth region, that flow is diverted leading to large vortices that are several R E across [e.g., Vörös et al., 2006]. The February 7, 2009 MHD simulation of the magnetosphere shows this explicitly. In the MHD results strong narrow flow channels form (Wiltberger et al., 2000;Ashour-Abdalla et al., 2002;Walker et al., 2006;Ge et al., 2011;Birn and Hesse, 2013) at the driving scales for turbulent vortices (e.g., El-Alaoui et al., 2013). In turn, the turbulent eddies can feedback on the reconnection process as well, leading to a very complex interplay between turbulence and dissipation (Matthaeus and Lamkin, 1986;El-Alaoui et al., 2012;Donato et al., 2013). The existence of turbulent eddies may contribute to the patchiness of the reconnection occurring in MHD simulations . Borovsky and Funsten (2003) investigated dissipation of the turbulence and argued for small-scale eddy dissipation where the vortices dissipate energy due to reconnection. Overall the results from all of these studies strongly support the idea that fast earthward flows drive large scale vortices that initiate turbulence.

SOME UNSOLVED QUESTIONS ABOUT TURBULENCE THAT CAN BE ADDRESSED USING MODELING
The simulation and data studies discussed make a strong case that flows associated with reconnection are turbulent. However, the dissipation in the MHD models used in this review depends on an additional term in Ohm's law. The extra term has free parameters so the fluid models provide little information about the physics of the dissipation region. One challenge in modeling turbulent flows in the tail is to include the physics of dissipation on a global scale. We need a model that can extend our basic understanding of kinetic turbulence and the transition between fluid and kinetic effects. Such a model would be a combined simulation that would resolve macroscopic motions and at the same time resolve ion and electron scales. It would cover the MHD turbulent cascade at large scales (inertial range) as well as Hall (dispersive) and kinetic (dissipation) ranges. Such simulations need to capture energetic exchanges at all scales. Fluid turbulence transfers the energy of large-scale flows to small-scale fluctuations and heat. As turbulence develops, the large-scale flow breaks up into vortices that cascade down to smaller and smaller scales where energy dissipates as heat. This represents a cascade of energy from the large energy-containing scales, through intermediate scales to small scales where heating occurs. In the magnetotail, the energy at the largest scales comes from the interaction between the magnetosphere and the solar wind, while at the smallest scales the dissipation is due to wave particle interactions and reconnection occurring on length scales of the ion Larmor radius or smaller. In the magnetotail, turbulence is complicated by the coupling between velocity and magnetic field fluctuations. Furthermore, turbulence in the kinetic regime involves specific normal modes of the plasma with spatial scales associated with specific waves, e.g., lower hybrid, or ion cyclotron waves.
We have recently developed (Walker et al., 2019) a model that combines simulations that resolve the macroscopic evolution and simulations that resolve ion and electron scales. The new method covers the MHD turbulent cascade at large scales (inertial range) as well as Hall (dispersive) and kinetic (dissipation) ranges. In the application to the study of turbulence, these simulations capture energetic exchanges at all scales. Fluid turbulence transfers the energy of large-scale flows to small-scale fluctuations and heat. As turbulence develops, the large-scale flow breaks up into vortices that cascade down to smaller and smaller scales where energy dissipates as heat. This represents a cascade of energy from the large energy-containing scales, through intermediate scales to small scales where heating occurs. At the magnetopause, the energy at the largest scales comes from the interaction between the magnetosphere and the solar wind, while at the smallest scales the dissipation is due to wave particle interactions and reconnection occurring on length scales of the ion Larmor radius or smaller. In the magnetotail, turbulence is complicated by the coupling between velocity and magnetic field fluctuations. Furthermore, turbulence in the kinetic regime FIGURE 11 | Time history of total small and mesoscale spectral power in the vicinity of flow channels and dipolarizations. For the region between −10 and −22R E in x and −4 to 8R E in y, the sum of the power in all x and y for mode numbers greater than two was computed (top two panels). The top panel shows the result for B z , in units of nT 2 , and the panel below it shows the result for V x , in units of (km/s) 2 . The lower two panels show the power, in units of (R E · km/s) 2 , vs. wavenumber at two times, 0344 UT on the left and 0353 UT on the right. (After El-Alaoui et al. (2016)).
Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2021 | Volume 8 | Article 620519 involves specific normal modes of the plasma with spatial scales associated with specific waves, e.g., lower hybrid, or ion cyclotron waves. It has long been recognized that the plasma distributions in the tail have non-thermal tails (a κ distribution). Recently, Lapenta et al. (2020b) have used this approach to model the acceleration of tail plasma to form this high-energy tail distribution. They argue that the main acceleration mechanism is turbulent acceleration. A similar conclusion was reached by Ergun et al. (2020) based on observations from the Magnetospheric Multi-Scale (MMS) mission.
The MMS observations come from four closely spaced satellites that provide the opportunity to investigate the turbulence in wave number space without the necessity of invoking the Taylor hypothesis. Combined with the latest kinetic simulation, these data (e.g., Ergun et al., 2018;Li et al., 2020;Ergun et al., 2020) offer great potential for understanding the importance of turbulence for acceleration in the tail and provide the community with a way to better understand the mechanisms of dissipation.
Specific issues that need to be examined include the scale of the dissipation region and its distribution in space. Fluid turbulence can develop into a state with dissipation localized in small regions of sharp gradients. How does this dissipation work in the plasma sheet when kinetic physics is included? How does dissipation occur in different regions? Is dissipation found near structures like thin current sheets (as in Wu et al., 2013), dipolarization fronts and the separatrices the result of distinct or the same processes? How do the flows become turbulent (e.g. velocity shears as in Ruffolo et al. (2020))? How important are turbulent flows for plasma heating by energy transfer from magnetic energy to particle and wave energy?
The magnetotail provides a unique laboratory for studying turbulence and reconnection. In the solar wind, turbulence has been investigated in depth for a long time (e.g. Coleman, 1968;Roberts et al., 1987a;Roberts et al., 1987b) but it is driven from the Sun over very large distances and long times but there is also evidence of reconnection driven turbulence in the solar wind (Gosling et al., 2010). In contrast in the magnetotail plasma sheet, the driving forces are well confined and comparatively shortterm, developing in a few minutes over just tens of thousands of km.