Identification of the weak-to-strong transition in Alfvénic turbulence from space plasma

The theory of anisotropic MHD turbulence has been widely accepted and adopted in plasma systems, ranging from clusters of galaxies, the interstellar medium, accretion disks, to the heliosphere Yan & Lazarian (2002); Brunetti & Lazarian (2007); Bruno & Carbone (2013). One of the most crucial predictions of the theory is an Alfvénic transition from weak to strong MHD turbulence when energy cascades from large to small scales Goldreich & Sridhar (1995); Howes et al. (2011). The self-organized process from weak to strong MHD turbulence is the cornerstone of understanding the energy cascade in the complete picture of MHD turbulence.

The critical balance model is an attractive model for describing physical behaviors of incompressible MHD (IMHD) turbulence Goldreich & Sridhar (1995). When $\tau_{A}\ll\tau_{nl}$ (referred to as weak MHD turbulence), weak interactions among the counter-propagating wave packets transfer energy to higher $k_{\perp}$, whereas no energy cascades to higher $k_{\parallel}$, where $\tau_{A}=1/(k_{\parallel}V_{A})$ is the linear Alfvén wave time, $\tau_{nl}=1/(k_{\perp}\delta V_{\perp})$ is the nonlinear time, $V_{A}$ is the Alfvén speed, $\delta V_{\perp}$ is the fluctuating velocity perpendicular to the background magnetic field ($\mathbf{B}_{0}$), and $k_{\perp}$ and $k_{\parallel}$ are wavenumbers perpendicular and parallel to $\mathbf{B}_{0}$ Galtier et al. (2000). As turbulence cascades to smaller scales, nonlinearity strengthens until reaching the critical balance $(\tau_{A}\approx\tau_{nl})$ at the transition scale ($\lambda_{CB}$). On scales smaller than $\lambda_{CB}$, Alfvén wave packets are statistically destroyed in one $\tau_{A}$. In addition to the first order interactions of counter-propagating waves, all higher orders of interactions can contribute, creating strong MHD turbulence Goldreich & Sridhar (1995); Mallet et al. (2015). In compressible MHD (CMHD), small-amplitude fluctuations can be decomposed into three eigenmodes (namely, Alfvén, fast, and slow modes) in homogeneous plasma Cho & Lazarian (2003); Makwana & Yan (2020); Zhu et al. (2020); Chaston et al. (2020); Zhao et al. (2021). Alfvén modes decoupled from CMHD turbulence, linearly independent of fast and slow modes, show similar properties to those in IMHD turbulence, e.g., the Kolmogorov spectrum and the scale-dependent anisotropy Cho & Lazarian (2003); Zhao et al. (2022). Additionally, numerical simulations have confirmed that the Alfvénic weak-to-strong transition occurs in both IMHD turbulence and Alfvén modes decomposed from CMHD turbulence Verdini & Grappin (2012); Meyrand et al. (2016); Makwana & Yan (2020).

However, such a transition has not been confirmed from observations. In this study, we present evidence for the Alfvénic weak-to-strong transition and estimate the transition scale $\lambda_{CB}$ in Earth’s magnetosheath using data from four Cluster spacecraft Escoubet et al. (2001). Earth’s magnetosheath offers a representative environment for studying plasma turbulence, given that most astrophysical and space plasmas with finite plasma $\beta$ are compressible, where $\beta$ is the ratio of the plasma to magnetic pressure.

Here we present an overview of fluctuations observed by Cluster-1 in geocentric-solar-ecliptic (GSE) coordinates during 23:00-10:00 Universal Time (UT) on 2-3 December in Fig. 1. During this period, four Cluster spacecraft flew in a tetrahedral formation with relative separation $d_{sc}\approx 200km$ (around 3 proton inertial length $d_{i}\approx 74km$) on the flank of Earth’s magnetosheath around [1.2, 18.2, -5.7] $R_{E}$ (Earth radius). We choose this time interval to study the Alfvénic weak-to-strong transition because the fluctuations satisfy the following criteria. Firstly, the background magnetic field ($\mathbf{B}$) measured by the Fluxgate Magnetometer (FGM) Balogh et al. (1997) and the proton bulk velocity ($\mathbf{V}_{p}$) measured by the Cluster Ion Spectrometry (CIS) Rème et al. (2001) are relatively stable in Figs. 1a-b. We cross-verify the reliability of plasma data, based on the consistency between the proton density ($N_{p}$) measured by CIS and electron density measured by the Waves of High frequency and Sounder for Probing of Electron density by Relaxation (WHISPER) Décréau et al. (1997) in Fig. 1c.

We set a moving time window with a five-hour length and a five-minute moving step. The selection of a five-hour length ensures that we obtain measurements at low frequencies (large scales) while the mean magnetic field ($\mathbf{B}_{0}$) within the moving time window is approaching the local mean field at the selected largest spatial scale. The uniform $\mathbf{B}_{0}$ is independent of the transformation between real and wavevector space; however, it differs from the theoretically expected local mean field at each scale. To assess such differences, Supplementary Fig. 1 shows that the local mean field at different scales is closely aligned with $\mathbf{B}_{0}$ most of the time, suggesting that $\mathbf{B}_{0}$ approximating the local mean field is acceptable. To further address this limitation of the mode decomposition method, which relies on a perturbative treatment of fluctuations in the presence of a uniform background magnetic field Cho & Lazarian (2003), we provide results obtained using various time window lengths, all of which show similar conclusions (Supplementary Fig. 2).

Fig. 1d shows spectral slopes of the trace magnetic field and proton velocity power calculated by fast Fourier transform (FFT) with three-point centered smoothing in each time window. These spectral slopes at spacecraft-frame frequency $f_{sc}\approx[0.001Hz,0.1f_{ci}]$ are close to $-5/3$ or $-3/2$ (the proton gyro-frequency $f_{ci}\approx 0.24Hz$), suggesting turbulent fluctuations are in a fully-developed state. The remaining magnetosheath fluctuations with spectra close to $f_{sc}^{-1}$ are typically populated by uncorrelated fluctuations Hadid et al. (2015, 2018) and are beyond the scope of the present paper. Fig. 1e shows the average proton plasma $\beta_{p}$ is around $1.4$. Finally, Fig. 1f shows that the turbulent Alfvén number $M_{A,turb}\equiv\delta V_{p}/V_{A}\approx\delta B/(2B_{0})\approx 0.33$, suggesting that fluctuations include substantial Alfvénic components and satisfy the small-amplitude fluctuation assumption (the nonlinear terms ($\delta V_{p}^{2}$, $\delta B^{2}$) are weaker than the linear terms ($V_{A}\delta V_{p}$, $B_{0}\delta B$)). Nevertheless, the average magnetic compressibility $C_{\parallel}(f_{sc})=\frac{|\delta B_{\parallel}(f_{sc})|^{2}}{|\delta B_{\parallel}(f_{sc})|^{2}+|\delta B_{\perp}(f_{sc})|^{2}}\approx 0.34$, indicating fluctuations are a mixture of Alfvén and compressible magnetosonic (fast and slow) modes Sahraoui et al. (2020), where $\delta B_{\parallel}$ and $\delta B_{\perp}$ are the fluctuating magnetic field parallel and perpendicular to $\mathbf{B}_{0}$.

Due to the homogeneous and stationary state of the turbulence (Supplementary Fig. 3), we can utilize frequency-wavenumber distributions of Alfvénic power, i.e. magnetic power $P_{B_{A}}(k_{\perp},k_{\parallel},f_{sc})$ and proton velocity power $P_{V_{A}}(k_{\perp},k_{\parallel},f_{sc})$, to investigate the structure of turbulence. Alfvénic fluctuations are extracted based on their incompressibility and fluctuating directions perpendicular to $\mathbf{B}_{0}$ (see Methods). The extraction is taken at each time window. To distinguish spatial and temporal evolutions without any spatiotemporal hypothesis, we determine wavevectors by combining the singular value decomposition method Santolík et al. (2003) (to obtain $\mathbf{\hat{k}}_{SVD}$) and multispacecraft timing analysis Pincon & Glassmeier (2008) (to obtain $\mathbf{\hat{k}}_{A}$). It is worth noting that $\mathbf{\hat{k}}_{A}$ is not completely aligned with $\mathbf{\hat{k}}_{SVD}$. Namely, $\mathbf{\hat{k}}_{A}$ may deviate from $\mathbf{\hat{k}}_{SVD}$ by angle $\eta$ (Supplementary Fig. 4). Thus, we present the results under $\eta<10^{\circ}$, $\eta<15^{\circ}$, $\eta<20^{\circ}$, $\eta<25^{\circ}$, and $\eta<30^{\circ}$. Given the marginal impact of different choices of $\eta$ (Supplementary Figs. 5 and 6), spectral results are displayed by taking the data set under $\eta<30^{\circ}$ as an example without loss of generality. The fluctuations, in this event, count for 42$\%$ of total Alfvénic fluctuations.

To ensure the reliability of wavenumber determination, we establish a minimum threshold of $k>1/(100d_{sc})$ and $k_{\parallel}>10^{-5}km^{-1}$. Consequently, our observations exclude the ideal two-dimensional (2D) case ($k_{\parallel}=0$), where $\tau_{A}$ is infinity, indicating persistent strong nonlinearity Galtier et al. (2003); Nazarenko (2007). Nevertheless, quasi-2D (small $k_{\parallel}$) modes are present, as $k_{\parallel}$ is much smaller than $k_{\perp}$ at small wavenumbers (Fig. 2). These quasi-2D fluctuations satisfies $\tau_{A}<\tau_{nl}$, as shown in Fig. 3c, and exhibit weak nonlinearity. This weak turbulent state occurs since $\delta B_{A}^{2}(k_{\perp},k_{\parallel})/B_{0}^{2}$ is very low, where Alfvénic magnetic energy density at $k_{\perp}$ and $k_{\parallel}$ is calculated by $\delta B_{A}^{2}(k_{\perp},k_{\parallel})=\sum_{k_{\perp}=k_{\perp}}^{k_{\perp}\rightarrow\infty}\sum_{k_{\parallel}=k_{\parallel}}^{k_{\parallel}\rightarrow\infty}\int_{0}^{\infty}P_{B_{A}}(k_{\perp},k_{\parallel},f_{sc})df_{sc}$.

Two-dimensional wavenumber distributions of magnetic energy are calculated by

$\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})=D_{B_{A}}(k_{\perp},k_{\parallel})/D_{B_{A},max}$ is normalized by the maximum magnetic energy in all $(k_{\perp},k_{\parallel})$ bins, displayed by the spectral image and contours in Fig. 2. Compared to the isotropic dotted curves, $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ is prominently distributed along the $k_{\perp}$ direction, suggesting a faster perpendicular cascade. This anisotropic behavior is more pronounced at higher wavenumbers, consistent with previous simulations and observations Cho & Vishniac (2000); He et al. (2011); Makwana & Yan (2020); Zhao et al. (2022).

Moreover, $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ is compared with the modeled 2D theoretical energy spectra based on strong turbulence Yan & Lazarian (2008); Goldreich & Sridhar (1995)

where the injection scale $L_{0}\approx[4.6\times 10^{4},8.1\times 10^{4}]km$ is approximately estimated by the correlation time $T_{c}\approx[1300,2300]s$ and rms perpendicular fluctuating velocity $\delta V_{p\perp}\approx M_{A,turb}V_{A}\approx 35kms^{-1}$. In Fig. 2, $\hat{I}_{A}(k_{\perp},k_{\parallel})$ is normalized $I_{A}(k_{\perp},k_{\parallel})$ by a constant value (one-third of the maximum magnetic energy in all $(k_{\perp},k_{\parallel})$ bins), displayed by color contours with black dashed curves. The 2D distribution $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ shows two different properties: (1) At $k_{\perp}<2\times 10^{-4}km^{-1}$, $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ is mainly concentrated at $k_{\parallel}<7\times 10^{-5}km^{-1}$ cascading along $k_{\perp}$ direction, suggesting that little energy cascade parallel to the background magnetic field, consistent with energy distributions in weak MHD turbulence Galtier et al. (2000). (2) At $k_{\perp}>2\times 10^{-4}km^{-1}$, $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ starts to distribute to higher $k_{\parallel}$, and both wavenumber distributions and intensity changes of $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ are almost consistent with $\hat{I}_{A}(k_{\perp},k_{\parallel})$. It indicates that $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ captures some theoretical characteristics of strong MHD turbulence Goldreich & Sridhar (1995). Besides, $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ is in good agreement with the Goldreich-Sridhar scaling $k_{\parallel}\propto k_{\perp}^{2/3}$ Goldreich & Sridhar (1995). This result further confirms that the properties of $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ at $k_{\perp}>2\times 10^{-4}km^{-1}$ are closer to those in strong MHD turbulence. The change in $\hat{D}_{B_{A}}(k_{\perp},k_{\parallel})$ from purely stretching along the $k_{\perp}$ direction to following the Goldreich-Sridhar scaling $k_{\parallel}\propto k_{\perp}^{2/3}$ reveals a possible transition in the energy cascade.

Fig. 3a shows the compensated spectra ($k_{\perp}^{5/3}E_{B_{A}}(k_{\perp})$), where the magnetic energy spectral density is defined as $E_{B_{A}}(k_{\perp})=\frac{\delta B_{A}^{2}(k_{\perp})}{2k_{\perp}}$, and $\delta B_{A}^{2}(k_{\perp})$ is magnetic energy density at $k_{\perp}$ (see Methods). In Zone (2), $k_{\perp}^{5/3}E_{B_{A}}(k_{\perp})$ is roughly consistent with $k_{\perp}^{5/3-2}$ (the dashed line), indicating that spectral slopes of $E_{B_{A}}(k_{\perp})$ are around $-2$. In Zone (3), on the other hand, $k_{\perp}^{5/3}E_{B_{A}}(k_{\perp})$ is almost flat, suggesting that $E_{B_{A}}(k_{\perp})$ satisfies the Kolmogorov scaling ($E_{B_{A}}(k_{\perp})\propto k_{\perp}^{-5/3}$). The sharp change in spectral slopes of $E_{B_{A}}(k_{\perp})$ from $-2$ to $-5/3$ is apparent evidence for the transition of turbulence regimes Verdini & Grappin (2012); Meyrand et al. (2016). In addition, $E_{B_{A}}(k_{\perp})\propto k_{\perp}^{-1}$ appears in a substantial portion of Zone (1), indicating the weak turbulence forcing in action Schekochihin et al. (2012); Makwana & Yan (2020).

Fig. 3b shows the variation of $k_{\parallel}$ versus $k_{\perp}$ given the same Alfvénic magnetic energy. As $k_{\perp}$ increases, $k_{\parallel}$ is relatively stable at $k_{\parallel}\approx 7\times 10^{-5}km^{-1}$ in Zone (1). In Zone (3), the variation of $k_{\perp}$ versus $k_{\parallel}$ agrees with the Goldreich-Sridhar scaling $k_{\parallel}\propto k_{\perp}^{2/3}$ (the dashed line). Fig. 3c shows $k_{\perp}-k_{\parallel}$ distributions of nonlinearity parameter ($\chi_{B_{A}}(k_{\perp},k_{\parallel})$), which is one of the most critical parameters in distinguishing between weak and strong MHD turbulence Howes et al. (2011), where $\chi_{B_{A}}(k_{\perp},k_{\parallel})$ is calculated by $\frac{k_{\perp}\delta B_{A}(k_{\perp},k_{\parallel})}{k_{\parallel}B_{0}}$ (see Methods). At the corresponding parallel and perpendicular wavenumbers in Fig. 3b, $\chi_{B_{A}}(k_{\perp},k_{\parallel})$ is much less than unity at most wavenumbers in Zone (1), whereas $\chi_{B_{A}}(k_{\perp},k_{\parallel})$ increases towards unity and follows the scaling $k_{\parallel}\propto k_{\perp}^{2/3}$ in Zone (3). These results suggest a transition from weak to strong nonlinear interactions, agreeing with theoretical expectations and simulations Howes et al. (2011); Verdini & Grappin (2012); Meyrand et al. (2016).

With the measurements of proton velocity fluctuations, we observe a similar Alfvénic weak-to-strong transition (Supplementary Fig. 7). The transition scale ($\lambda_{CB}$) is estimated by the smallest perpendicular wavenumber of strong turbulence ($k_{\perp,CB}$), where $\lambda_{CB}\approx 1/k_{\perp,CB}$. For both magnetic field and proton velocity fluctuations, $k_{\perp,CB}$ is around $3\times 10^{-4}km^{-1}$, marked by the second vertical lines in Fig. 3 and Supplementary Fig. 7. The consistency in the transition scales estimated by magnetic field and proton velocity measurements further confirms the reliability of our findings.

A notable perturbation is present in Zone (2), as a result of local enhancements of magnetic energy at $k_{\perp}\approx 1.8\times 10^{-4}km^{-1}$ (Fig. 2), leading to the simultaneous existence of strong nonlinearity ($\chi_{B_{A}}\approx 1$) and weak nonlinearity ($\chi_{B_{A}}\ll 1$) in the wave number range corresponding to those in Fig.3b. Thus, the Alfvénic weak-to-strong transition more likely occurs within a ‘region’ rather than at a critical wavenumber. Besides, we do not discuss the fluctuations in Zone (4). The deviations of data sets under $\eta<10^{\circ}$ and $\eta<15^{\circ}$ in Zone (3) of Fig. 3b are likely due to the limited data samples (Supplementary Fig. 6). The uncertainties mentioned above do not affect our main conclusions.

Fig. 4 presents $k_{\perp}$ versus $f_{rest}$ distributions of magnetic energy, where $f_{rest}$ is the frequency in the plasma flow frame. At $k_{\perp}<5\times 10^{-5}km^{-1}$, magnetic energy is concentrated at $f_{rest}\approx f_{A}$, where $f_{A}$ is Alfvén frequency (horizontal dotted lines with error bars). At $k_{\perp}>1\times 10^{-4}km^{-1}$, the range of $f_{rest}$ broadens, mostly deviating from $f_{A}$. Nevertheless, the boundary of fluctuating frequencies is roughly consistent with the scaling $f_{rest}\propto k_{\perp}^{2/3}$ (the dashed line), indicating that magnetic energy at these wavenumbers satisfies the scaling $k_{\parallel}\propto k_{\perp}^{2/3}$ due to $f_{rest}\propto k_{\parallel}$ for Alfvén modes. These results suggest that Alfvénic fluctuations with strong nonlinear interactions do not agree with linear dispersion relations but satisfy the wavenumber scaling of Alfvén modes. The change from single-frequency to broadening-frequency fluctuations with increasing $k_{\perp}$ suggests a possible transition of turbulence regimes.

The Alfvénic transition of weak to strong turbulence during cascades to smaller scales is one of the cornerstones of the modern MHD theory. Despite being proposed decades ago, evidence for confirming the existence of the Alfvénic transition is lacking. In this paper, we present direct evidence of the Alfvénic transition via different angles: e.g., the transition of energy spectra (Fig. 3a), Goldreich-Sridhar type envelope for the nonlinear parameter (Fig. 3c), and the spread of $f_{rest}$ on small scales (Fig. 4; See Tab. 1 for a summary). Our observation demonstrates that the Alfvénic transition to strong turbulence is bound to occur with the increase of nonlinearity even fluctuations on large scales are considered as ”small amplitude” ($M_{A,turb}\approx 0.33$). We want to point out that plasma parameters in the analyzed event are generic, and Alfvénic weak-to-strong transition can occur in other astrophysical and space plasma systems. The impact of our findings goes beyond the study of turbulence itself to particle transport and acceleration Schlickeiser (2002); Yan (2021), magnetic reconnection Matthaeus & Lamkin (1986); Lazarian & Vishniac (1999), star formation Crutcher (2012); Padoan et al. (2014), and all the other relevant fields (see, e.g. Zhang & Yan, 2011; Hirashita & Yan, 2009).

It is worth noting that timing analysis determines the actual wavevectors of the Alfvénic magnetic field. In contrast, the SVD method determines the best estimate of the wavevector sum in three magnetic field components Santolík et al. (2003). Thus, $\mathbf{k}_{A}$ is not completely aligned with $\hat{\mathbf{k}}_{SVD}$. Besides, we restrict our analysis to fluctuations with small angle $\eta$ between $\hat{\mathbf{k}}_{SVD}$ and $\mathbf{k}_{A}$, to ensure the reliability of the extraction process (the third step). With relaxed $\eta$ constraints, more sampling points are involved; thus, the uncertainty from limited measurements decreases. On the other hand, with relaxed $\eta$ constraints, $k_{A}$ deviates more from $k_{SVD}$, which may increase the uncertainty. This letter presents results from five data sets under $\eta<10^{\circ}$, $\eta<15^{\circ}$, $\eta<20^{\circ}$, $\eta<25^{\circ}$, and $\eta<30^{\circ}$ to investigate the effects of uncertainties introduced by the combination of the SVD method and timing analysis.

$P_{\epsilon_{A}}(k_{\perp},k_{\parallel},f_{sc})$ are obtained by averaging $P_{\epsilon_{A}}(k_{\perp},k_{\parallel},f_{sc},t)$ over effective time points in all time windows at each $f_{sc}$ and each $\mathbf{k}$, where $\epsilon=V,B$ represents the proton velocity ($V$) and magnetic field ($B$).

The nonlinearity parameter is estimated by $\chi_{B_{A}}(k_{\perp},k_{\parallel})\approx k_{\perp}\delta B_{A}(k_{\perp},k_{\parallel})/(k_{\parallel}B_{0})$, where the Alfvénic magnetic energy density is calculated by $\delta B_{A}^{2}(k_{\perp},k_{\parallel})=\sum_{k_{\perp}=k_{\perp}}^{k_{\perp}\rightarrow\infty}\sum_{k_{\parallel}=k_{\parallel}}^{k_{\parallel}\rightarrow\infty}\int_{0}^{\infty}P_{B_{A}}(k_{\perp},k_{\parallel},f_{sc})df_{sc}$, and $B_{0}$ in Alfvén speed units is around $106km/s$.