Modeling MMS Observations at the Earth's Magnetopause with Hybrid Simulations of Alfvénic Turbulence

The observational data set consists of 54 subintervals of high-resolution "burst" data collected by MMS on the duskside magnetopause on 2015 September 8 over a ∼1.5 hr long interval in which a continuous train of KH waves were observed. Such subintervals are the same as those examined by Stawarz et al. (2016) and were manually selected such that they were all located within individual KH vortices. The statistics (e.g., spectra) presented in this study are taken to be the ensemble average based on these 54 intervals. This choice allows us to investigate the properties of the turbulent plasma avoiding the periodic compressed current sheets associated with the large-scale KH wave, which would skew the statistics if included in the analysis. As a consequence of such selection, the power spectra (see Section 3.1) represent the ensemble average of turbulent fluctuations inside each of the vortices and no enhancement in the power corresponding to the scale of the KH vortices is expected to be present at the examined frequencies. In this KH event, the magnetosheath magnetic field is northward and the field strength is nearly constant across the magnetopause. The "guide field" is taken to be the average field ${{\boldsymbol{B}}}_{0}=\langle {\boldsymbol{B}}\rangle$ within each subinterval. The magnetic field is dominated by the ${{\boldsymbol{B}}}_{z\mathrm{GSE}}$ component and so the perpendicular plane associated with the turbulence is nearly the x_GSE−y_GSE plane.

The data set contains magnetic field measurements from the fluxgate (Russell et al. 2016) and the search coil (Le Contel et al. 2016) magnetometers at 1/128 and 1/8192 s cadences respectively, electric field measurements from the electric field double probes (Ergun et al. 2016; Lindqvist et al. 2016) at 1/8192 s cadence, and ion and electron particle moments, as measured by the fast plasma investigation (Pollock et al. 2016), at 0.15 and 0.03 s cadences respectively. Prior to averaging, the magnetic field was normalized to B₀, the electric field to B₀V_A, the ion and electron bulk velocities to V_A, and the electron density to n₀, where the Alfvén velocity, V_A, and the background density, n₀, are defined based on average parameters from each subinterval.

Aspects of this KH event, including magnetic reconnection (Eriksson et al. 2016a, 2016b; Li et al. 2016; Vernisse et al. 2016; Sturner et al. 2018), turbulence (Stawarz et al. 2016; Sorriso-Valvo et al. 2019), waves (Wilder et al. 2016), and mass transport (Nakamura et al. 2017a, 2017b) properties, have been examined in a number of previous studies. Eriksson et al. (2016a) provide an overview of the properties of the overall KHI within this event, which is characterized by an average ion inertial length of ${d}_{{\rm{i}}}\sim 65$ km and a KH wavelength of the order of 300 d_i.

The numerical data set was obtained from a high-resolution 2D DNS of plasma turbulence performed with the hybrid (particle-in-cell ions and mass-less fluid electrons) code CAMELIA (Franci et al. 2018a). The temporal and spatial units are the inverse ion gyrofrequency, ${{\rm{\Omega }}}_{{\rm{i}}}^{-1}$, and the ion inertial length, d_i. The magnetic field is expressed in units of the magnitude of the ambient magnetic field, B₀, the ion bulk velocity in units of the Alfvén velocity, V_A, and the ion density in units of the background density, n₀, as for the respective observational quantities. The initial conditions are the same as those in Franci et al. (2015a, 2015b). These consist of a homogeneous plasma in the (x, y) simulation plane, embedded in a uniform out-of-plane ambient magnetic field, ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{z}$. The simulation plane corresponds to the plane perpendicular to ${{\boldsymbol{B}}}_{0}$ in the observations, which is nearly the x_GSE−y_GSE plane, as noted above. The plasma is perturbed by in-plane balanced Alfvénic fluctuations, i.e., magnetic and ion bulk velocity fluctuations with zero mean cross helicity. These are the superposition of modes with the same amplitude and random phases in the range $-0.2\leqslant {k}_{x,y}{d}_{{\rm{i}}}\leqslant 0.2$, so that ${k}_{\perp }^{\max }{d}_{{\rm{i}}}\sim 0.3$, which we will refer to as the "injection scale." The parallel wavevector with respect to ${{\boldsymbol{B}}}_{0}$ is by construction ${{\boldsymbol{k}}}_{\parallel }=0$, so that ${\boldsymbol{k}}\equiv {{\boldsymbol{k}}}_{\perp }$. This, however, still allows for parallel propagation with respect to the local mean field, as a consequence of the bending of the magnetic field lines due to the in-plane magnetic fluctuations. Our initialization does not intend to be representative of the real environment explored by MMS since the KHI, which is the driving mechanism of the observed turbulence, is not included in our DNS. We do, however, reproduce similar plasma conditions by setting the following physical parameters accordingly to their observational values: ion beta β_i = 0.42, electron beta ${\beta }_{{\rm{e}}}=0.065$, amplitude of the initial fluctuations ${{\boldsymbol{B}}}^{\mathrm{rms}}/{{\boldsymbol{B}}}_{0}\sim 0.14$, and injection scale ${k}_{\perp }^{\max }{d}_{{\rm{i}}}\sim 0.3$. The numerical parameters are 16384 particle-per-cell (ppc), 4096 × 4096 grid points, box size ${L}_{x}={L}_{y}=256\,{d}_{{\rm{i}}}$ (of the same order of the KH wavelength for the observed event, see above), spatial resolution Δx = Δy = d_i/16, time resolution ${\rm{\Delta }}t=0.01\,{{\rm{\Omega }}}_{{\rm{i}}}^{-1}$ for the particle advance and ${\rm{\Delta }}{t}_{B}={\rm{\Delta }}t/10$ for the magnetic field advance, and resistivity $\eta =2\times {10}^{-4}\,4\pi {V}_{{\rm{A}}}{c}^{-1}{{\rm{\Omega }}}_{{\rm{i}}}^{-1}$. The optimal value for η has been chosen empirically, based on our previous simulations with similar parameters. This prevents an accumulation of the energy in magnetic fluctuations at small scales due to numerical noise and assures a satisfactory conservation of the total energy (Franci et al. 2015a). The development of the turbulent cascade is qualitatively the same as in Franci et al. (2017). The simulation has been run until $350\,{{\rm{\Omega }}}_{{\rm{i}}}^{-1}$, when the rms value of the current density has reached a plateau. At this time, a turbulent quasi-stationary state has been achieved, characterized by magnetic field structures at all scales (Figure 1(a)) and small-scale current structures between vortices (Figure 1(b)). We perform our analysis at $300\,{{\rm{\Omega }}}_{{\rm{i}}}^{-1}$, without performing any time average. We have verified that spectra and other properties do not change significantly between that time and the end of the simulation.

Zoom In Zoom Out Reset image size

The observational and numerical spectral properties of the electromagnetic and plasma fluctuations are shown in Figures 2 and 3, respectively, using red dashed–dotted lines for MMS and black solid ones for the DNS. Four vertical dashed lines mark the particle characteristic scales, i.e., from small to large wavenumbers, the ion inertial length, d_i, the ion gyroradius, ${\rho }_{{\rm{i}}}=\sqrt{{\beta }_{{\rm{i}}}}{d}_{{\rm{i}}}$, the electron inertial length, d_e, and the electron gyroradius, ${\rho }_{{\rm{e}}}=\sqrt{{\beta }_{{\rm{e}}}}{d}_{{\rm{e}}}$. In the top panels, the 1D power spectra of each field are compared. The MMS spectra have been converted from frequency, f, to perpendicular wavenumber with respect to the ambient magnetic field, k_⊥, using the "frozen flow" Taylor hypothesis (Taylor 1938). Under the assumptions that the average perpendicular ion bulk flow velocity computed over the 54 subintervals, ${U}_{0,\perp }=170\,\mathrm{km}\,{{\rm{s}}}^{-1}$, is much larger than its parallel counterpart, ${U}_{0,\parallel }=33\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Stawarz et al. 2016) and that k_⊥ ≫ k_∥ (consistent with Alfvénic turbulence), such a hypothesis yields ${k}_{\perp }=2\pi f/{U}_{0,\perp }$ (Bourouaine et al. 2012). The difference between this estimate of k_⊥ and the longitudinal wavenumber estimated as 2π f/U₀ is here only 2%, indicating that the spacecraft is essentially sampling k_⊥ as the plasma flows past it. Therefore, in the remaining analysis we will utilize the estimate of k_⊥ obtained under the assumption that k_⊥ ≫ k_∥, as this is the component of the wavevector most compatible with the DNS. The DNS spectra have not been rescaled in amplitude and their portions affected by numerical noise are drawn as dotted lines. In the bottom panels we compare the values of the local spectral index, α, obtained by performing many power-law fits on small intervals in ${k}_{\perp }$. For each field, horizontal power laws with characteristic slopes are drawn as a reference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 2(a), we analyze the fluctuations of the magnetic field, ${\boldsymbol{B}}$. The MMS spectrum exhibits a triple power-law-like behavior over four decades in wavenumber, with a first steepening around the ion scales and a second one around the electron scales. It is closely followed by the DNS spectrum over two full decades where the numerical noise is negligible, i.e., in the wavenumber interval $0.3\lesssim {k}_{\perp }{d}_{{\rm{i}}}\lesssim 30$. The bottom panel shows that the first power-law slope is compatible with the −5/3 prediction for the inertial-range turbulent cascade (Goldreich & Sridhar 1995). A transition is observed in correspondence of k_⊥ d_i ∼ 3, consistently with previous DNS with a similar plasma beta (Franci et al. 2015a, 2015b, 2016). A second, steeper (α ∼ −3.2), power law is observed at sub-ion scales, extended for over a decade and a half in MMS, up to k_⊥ d_i ∼ 100, and over slightly less than a decade in the DNS, until the resolution limit is reached. A second transition is observed in the data just before the electron scales are met, followed by a third, much narrower interval where the instrumental noise becomes important, so that no clear power law can be inferred, despite a hint of a further steepening.

Three distinct, clearer, power-law intervals are observed in the MMS spectrum of the electric field, ${\boldsymbol{E}}$ (Figure 2(b)). The agreement with the DNS spectrum is particularly remarkable, since in the hybrid model ${\boldsymbol{E}}$ is computed from the generalized Ohm's law involving the other fields and their derivatives, so it is more sensitive to both physical parameters (e.g., ion and electron beta) and small-scale numerical noise (e.g., from the density and the current density). The first interval is compatible with a power law with slope −3/2, followed by a flattening to ∼−0.8 at sub-ion scales. Both these spectral indices are consistent with theoretical predictions from the generalized Ohm's law (Franci et al. 2015a). In correspondence with the electron scales, the MMS spectrum shows a second transition and a steepening, hinting at a third power-law interval with slope ∼ − 2.8, although less extended than a decade, for 300 ≲ k_⊥d_i ≲ 1000.

The spectral properties of the density fluctuations are compared in Figure 3(a). MMS measures both the ion, ${n}_{{\rm{i}}}$, and the electron density fluctuations, n_e, with two different instruments and resolutions. In the hybrid model, since the electrons are treated as a mass-less fluid, these two quantities are assumed to be the same. Therefore, here we only show the observed electron density, which has a higher time resolution, corresponding to the DNS spatial resolution. In the inertial range, the MMS and DNS density fluctuations exhibit a very different spectral behavior, due to the fact that the latter is initialized with no density fluctuations and this keeps the large-scale compressibility lower. The MMS spectrum exhibits a power-law behavior with a slope compatible with -1, consistent with previous spacecraft observation in the solar wind (Šafránková et al. 2015). At ion and sub-ion scales, a small difference in the level of fluctuations between MMS and the DNS is observed, of the same order of the one observed for the electron velocity spectra. Nevertheless, a good agreement is recovered for the slope below the ion scales, which is consistent with −2.8. This confirms that the level of compressibility in the inertial range does not impact significantly the nature of the fluctuations in the kinetic range (Cerri et al. 2017). The shape of the MMS spectrum of n_i (not shown) is the same as for n_e in the whole range of resolved scales: the slope is the same, i.e., ∼−1 at large scales and ∼−2.8 at small scales, while the level of fluctuations is slightly lower, likely due to the MMS FPI density estimates, and in better agreement with, but still slightly overestimating, the DNS spectrum at sub-ion scales.

The spectra of the ion bulk velocity, ${{\boldsymbol{U}}}_{{\rm{i}}}$ (Figure 3(b)), show no clear extended power-law interval in either the inertial range or the kinetic one. The level of the ion bulk velocity fluctuations strongly drops at k_⊥d_i ∼ 1, so there is not even a full decade separating the break from the injection scale. Approaching the kinetic scales, both the MMS and the DNS spectra flatten due to instrumental and numerical noise, respectively. In the range of scales not affected by noise, the two spectra are still in quite good qualitative agreement.

Figure 3(c) compares the MMS and DNS spectra of the electron bulk velocity, which in the hybrid model is computed as ${{\boldsymbol{U}}}_{{\rm{e}}}={{\boldsymbol{U}}}_{{\rm{i}}}-({\rm{\nabla }}\times {\boldsymbol{B}})/n$. By chance, the maximum wavenumber is the same, since the spatial resolution of the DNS corresponds to the observational time resolution for this field. A small discrepancy in the level of fluctuations is observed over the whole range of scales, whose origin is not clear. The DNS spectrum of ${{\boldsymbol{U}}}_{{\rm{e}}}$ in the kinetic range is just that of ${\boldsymbol{B}}$ multiplied by ${k}_{\perp }^{2}$, consistently with the fact that the ion bulk velocity fluctuations are negligible at those scales and so the current is almost entirely supported by the electron bulk motion, i.e., ${\rm{\nabla }}\times {\boldsymbol{B}}={\boldsymbol{J}}={n}_{{\rm{i}}}{{\boldsymbol{U}}}_{{\rm{i}}}-{n}_{{\rm{e}}}{{\boldsymbol{U}}}_{{\rm{e}}}\sim -n\,{{\boldsymbol{U}}}_{{\rm{e}}}$. Despite the slightly different level, the slope of the two spectra of ${{\boldsymbol{U}}}_{{\rm{e}}}$ is almost exactly the same at all scales. The only difference is a small flattening in the DNS spectrum around the ion scales, which brings the level of fluctuations at kinetic scales slightly closer to the observational one. This can be seen as a small bump in the value of the local slope, visible in the bottom panel, between the scale corresponding to d_i and the scale of the break in the magnetic field spectrum.

The slopes of the different fields in the sub-ion range, i.e., 3 ≲ k_⊥d_i  ≲  40, are not all independent. The one of the electric field, ${\alpha }_{{\boldsymbol{E}}}\sim -0.8$, is related to the one of the density and of the parallel magnetic fluctuations (not shown), ${\alpha }_{n}\sim {\alpha }_{{{\boldsymbol{B}}}_{\parallel }}\sim -2.8$, by the simple relation ${\alpha }_{{\boldsymbol{E}}}\sim {\alpha }_{n({{\boldsymbol{B}}}_{\parallel })}-2$. This comes from the electron pressure gradient term and the Hall term in the generalized Ohm's law, the latter being dominated by its first-order contribution $({\rm{\nabla }}\times {{\boldsymbol{B}}}_{\parallel })\times {{\boldsymbol{B}}}_{0}$ (Franci et al. 2015a). The kinetic-scale slope of the electron bulk velocity spectrum is instead related to the one of the perpendicular magnetic fluctuations (not shown, but similar to the total magnetic field, representing its dominant contribution), ${\alpha }_{{{\boldsymbol{B}}}_{\perp }}\sim -3.2$, according to ${\alpha }_{{{\boldsymbol{U}}}_{{\rm{e}}}}\sim {\alpha }_{{{\boldsymbol{B}}}_{\perp }}+2$. This is due to the fact that the out-of-plane current density is much larger than its in-plane component and the ion bulk velocity drops at ion scales, so ${\boldsymbol{J}}\sim {{\boldsymbol{J}}}_{\parallel }\sim {\rm{\nabla }}\times {{\boldsymbol{B}}}_{\perp }\sim {n}_{0}{{\boldsymbol{U}}}_{{\rm{e}}}$.

The intermittency properties are compared in Figure 4, showing the PDFs of the in-plane increments of a magnetic field component perpendicular to ${{\boldsymbol{B}}}_{0}$, $\delta {B}_{{\rm{y}}}({k}_{\perp })/\delta {B}_{{\rm{y}},\mathrm{rms}}$. Six different values of k_⊥d_i are explored, representing the inertial range, k_⊥d_i = 0.6 (panel (a)) and 1.2 (b), the ion-scale transition, k_⊥d_i = 2.4 (c) and 2.33 (f), and the ion kinetic range, ${k}_{\perp }{d}_{{\rm{i}}}=4.8$ (d) and 9.6 (e). For the DNS (black curves), these are simply obtained from k_⊥ = 2π/ℓ, where ℓ is the spatial lag in the simulation plane. For MMS (red curves), this only holds for panel (f), where the increments are actually computed from the spatial lag, while for (a)–(e) the temporal lag, τ, has been converted to the corresponding wavenumber using the Taylor hypothesis.

Zoom In Zoom Out Reset image size

The PDFs of magnetic fluctuations present very different shapes at different scales. At scales close to the injection scale, it is quite similar to a Gaussian distribution (dotted black curve), deviating from it only in correspondence of the largest fluctuations. The PDFs become less and less Gaussian toward smaller lags, developing extended tails at kinetic scales. For fluctuations with an amplitude smaller than 4 times the rms value, $\delta {B}_{{\rm{y}}}({k}_{\perp })\lt 4\,\delta {B}_{{\rm{y}},\mathrm{rms}}$, the numerical and observational PDFs are almost overlapped at all scales, apart from small differences in the kinetic range. In the tails, the PDFs from the DNS are usually higher than their observational counterpart (such differences are almost negligible if we consider that there the number of counts is small and concentrated in very few points in the simulation domain). Such intermittency properties are consistent with the excess kurtosis, shown in Figure 5. Both the DNS (black) and the MMS (red) kurtosis are very small at large scales, down to k_⊥ d_i ∼ 1. They increase in the kinetic range, with a similar but still different trend, such that their values at k_⊥ d_i ∼ 10 differ by a factor of ∼2. The DNS kurtosis is seen to be very sensitive to the tails of the PDFs, where the results suffer from resolution effects. Moreover, at earlier times, it assumes larger values by about 2 orders of magnitude, peaking at more than 100 in correspondence with the onset of magnetic reconnection, when strong current sheets are just about to start being disrupted.

Zoom In Zoom Out Reset image size

The cascade rate properties can be investigated by using the statistical von Kármán–Howarth/Politano–Pouquet law (de Kármán & Howarth 1938; Politano & Pouquet 1998) for the mixed third-order structure function

$$\begin{eqnarray}&&{\boldsymbol{Y}}=\left\langle \delta {\boldsymbol{U}}{\left|\delta {\boldsymbol{U}}\right|}^{2}+\delta {\boldsymbol{U}}{\left|\delta {\boldsymbol{b}}\right|}^{2}-2\delta {\boldsymbol{b}}\left(\delta {\boldsymbol{U}}\cdot \delta {\boldsymbol{b}}\right)\right\rangle ,\end{eqnarray} \tag{ 1 }$$

where the increments of the bulk velocity, $\delta {\boldsymbol{U}}$, and of the magnetic field in Alfvénic units, $\delta {\boldsymbol{b}}$, are computed between two points x and x + l and $\langle \,\rangle$ denotes the average over x. This law predicts for the MHD inertial range a linear scaling of Y_L with the separation scale l, where Y_L is longitudinal component. Figure 6 compares −Y_L/l as a function of l from MMS and the DNS. Both curves exhibit a similar profile with a plateau, but they do not show a clear linear scaling of Y_L (or a constant value of $-{Y}_{L}/l$), a signature of the inertial range. This is likely due to the fact that the energy injection takes place at scales not far enough from the ion ones and in the sub-ion range the Hall physics must be included (Papini et al. 2019). Indeed, recent simulations and observations indicate that the turbulent energy cascade partly continues via the Hall term (Hellinger et al. 2018; Bandyopadhyay et al. 2020). Figure 6, however, shows that MMS and the DNS exhibit very similar values (for the separation length in the range between about 50–200 km). Assuming isotropy, we get (at the plateau) a cascade rate of $\sim 5\times {10}^{7}\,{\rm{J}}\,{\mathrm{kg}}^{-1}\,{{\rm{s}}}^{-1}$, a value about 10 times larger than the one observed in the Earth's magnetosheath (Bandyopadhyay et al. 2018; Hadid et al. 2018).

Zoom In Zoom Out Reset image size