THE FORMATION OF ROTATIONAL DISCONTINUITIES IN COMPRESSIVE THREE-DIMENSIONAL MHD TURBULENCE

Since the beginning of the space age, the solar wind has been regarded as an excellent natural laboratory for studying the plasma turbulence. The endeavoring research over the past decades has revealed that discontinuities are ubiquitous in the solar wind (Burlaga 1969, 1971; Smith 1973; Solodyna et al. 1977; Tsurutani & Smith 1979; Neugebauer et al. 1984; Lepping & Behannon 1986; Tsurutani et al. 1994, 1996, 2005a, 2005b; Lee et al. 1996; Horbury et al. 2001; Knetter et al. 2004; Neugebauer 2006; Vasquez et al. 2007; Li 2008; Lin et al. 2009; Sonnerup et al. 2010; Teh et al. 2011; Haaland et al. 2012; Malaspina & Gosling 2012; Paschmann et al. 2013; Wang et al. 2013). These discontinuities appear as large and rapid changes in the properties of the plasma and magnetic field and are identified as statically advected tangential discontinuities (TDs), or propagating rotational discontinuities (RDs). The TDs are characterized by small normal components of the magnetic field, large variations of magnetic field intensity, and density jumps, while the RDs have large normal components of the magnetic field, but small variations of magnetic field intensity and density (Hudson 1970).

Using measurements of the magnetic field from Pioneer 6, Burlaga (1969) observed the discontinuities in the magnetic field direction with special emphasis on their distribution in time. Based on interplanetary field measurements made by the vector helium magnetometers on board Pioneer 10 and Pioneer 11, Tsurutani & Smith (1979) investigated a possible dependence of the occurrence rate and the properties of the discontinuities on radial distance between 1 and 8.5 AU. Lepping & Behannon (1986) surveyed the data from the Mariner 10 primary mission to study the characteristics of the discontinuities in the interplanetary magnetic field at heliographic distances of 1.0, 0.72, and 0.46 AU, and found an ${r}^{-1.3\pm 0.4}$ dependence for the daily average number of discontinuities per hour. With Ulysses magnetic field and plasma data obtained at radial distances ranging between 1 and 5 AU from the Sun and at high heliographic latitudes, Tsurutani et al. (1994) discovered two regions where the occurrence rate of interplanetary discontinuities is high: in stream–stream interaction regions and in Alfvén wave trains. To determine the normals of the discontinuities, Horbury et al. (2001) explored the discontinuities measured by three spacecraft, WIND, IMP 8, and Geotail, together with the solar wind velocity measured at Geotail, and obtained quite different distributions of the discontinuity types. With magnetic field data from the ACE spacecraft, Vasquez et al. (2007) extended the survey of discontinuity properties to small spread angles of the field vectors across the discontinuity, and found that solar wind discontinuities are far more abundant at smaller than at larger spread angles. Using measurements from the WIND spacecraft, Wang et al. (2013) studied the intermittent structures in solar wind turbulence, which are identified as being mostly RDs and rarely TDs based on the technique described by Smith (1973) and Tsurutani & Smith (1979). Paschmann et al. (2013) carried out a comprehensive study of directional discontinuities and Alfvénic fluctuations in the solar wind on the basis of Cluster data.

So far there are still debates regarding the origin and nature of the discontinuities in the solar wind. Nonlinear wave steepening has been suggested as the cause of its formation from both experimental observations and theoretical studies (Cohen & Kulsrud 1974; Malara & Elaoufir 1991; Tsurutani et al. 1994, 1996, 1997; Medvedev & Diamond 1996; Vasquez & Hollweg 1996, 1998a, 1998b; Medvedev et al. 1997; Buti et al. 1998, 2001; Tsurutani & Ho 1999; Tsubouchi & Matsumoto 2005; Tsubouchi 2009; Marsch & Verscharen 2011). By analyzing interplanetary magnetic field data measured by the Ulysses, Tsurutani et al. (1994, 1996, 1997) first showed that Alfvén waves in the high-speed streams are phase-steepened, with leading edges that form an RD. Furthermore, they found that the Alfvén waves have constant magnetic field magnitudes, and the magnetic field vectors swing back and forth in arc, which means that the observed Alfvén waves are essentially spherical waves, i.e., the magnetic perturbation vector rotates on the surface of a sphere rather than planar polarization. By numerically calculating the evolution of an initially parallel-propagating, elliptically polarized wave train in a cold plasma, Cohen & Kulsrud (1974) were the first to investigate this possibility, and showed that this wave evolves into a constant-B solution with RDs that rotate the field by exactly 180°. Vasquez & Hollweg (1996) continued this study, but conducted a 1.5D hybrid numerical simulation study of the evolution of obliquely propagating, linearly polarized Alfvén wave trains. They found that large-amplitude ${dB}/{B}_{0}\gt 1$ wave trains steepen and produce RDs that always rotate the field by $\lt 180^\circ$. Vasquez & Hollweg (1998a) also presented 2.5D numerical simulations of a small group of nonplanar Alfvén waves to show the generation of embedded RDs. It should be noted that in their models the formation of RDs probably occurs relatively near the Sun, where most Alfvénic fluctuations originate.

There is another model suggesting that MHD turbulence dynamically generates these discontinuities as the solar wind flows outward. Recently, numerical simulations have been done to investigate this assertion (Greco et al. 2008, 2009, 2010; Servidio et al. 2011; Zhdankin et al. 2012a, 2012b). Through analyses of MHD simulation data, Greco et al. (2008) examined the relationship between discontinuities identified by classical methods, and coherent structures identified by using intermittency statistics. They found that the two methods produce remarkably similar distributions of waiting times, and in fact identify many of the same events. Greco et al. (2009) further examined the link between intermittent turbulence and MHD discontinuities, directly comparing simulations of MHD turbulence with statistical analysis from ACE solar wind data. Their results support the notion that some solar wind discontinuities are consequences of intermittent turbulence. In direct numerical simulations of MHD turbulence with an imposed uniform magnetic field, Zhdankin et al. (2012a) investigated the statistical properties of magnetic discontinuities and concluded that the discontinuities observed in the solar wind can be reproduced by MHD turbulence. However, these works conducted statistical studies, and did not give a clear illustration of how discontinuities, especially RDs, are formed in MHD turbulence.

In the present study, we utilize a compressible three-dimensional (3D) MHD model to illustrate and analyze the formation of RDs in the turbulence. By checking the magnetic field and plasma properties, it is found that the RD is produced by the steepening of a moderate-amplitude Alfvén wave with nonuniform propagating speed. The paper is organized as follows. In Section 2, a general description of the numerical MHD model is given. Section 3 describes the results of the numerical simulation and its analysis. Section 4 is reserved for the summary and discussion.

The description of the plasma is given by a compressible 3D MHD, which involves a fluctuating flow velocity ${\boldsymbol{v}}(x,y,z,t)$, magnetic field ${\boldsymbol{b}}(x,y,z,t)$, density $\rho (x,y,z,t)$, and temperature $T(x,y,z,t)$. A uniform guide field B₀ is assumed in the z-direction, so the total magnetic field is ${\boldsymbol{B}}={{\boldsymbol{B}}}_{0}+{\boldsymbol{b}}$.

The MHD equations are written in the following non-dimensional form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot \rho {\boldsymbol{u}}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{u}}}{\partial t}+{\rm{\nabla }}\cdot \left[\rho {\boldsymbol{uu}}+(p+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2}){\boldsymbol{I}}-{\boldsymbol{BB}}\right]=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\rm{\nabla }}\cdot \left[{\boldsymbol{u}}(e+p+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2})-({\boldsymbol{u}}\cdot {\boldsymbol{B}}){\boldsymbol{B}}\right]={\rm{\nabla }}\cdot ({\boldsymbol{B}}\times \eta {\boldsymbol{j}}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{uB}}-{\boldsymbol{Bu}})=\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&e=\displaystyle \frac{1}{2}\rho {{\boldsymbol{u}}}^{2}+\displaystyle \frac{p}{\gamma -1}+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2},\ \ {\boldsymbol{j}}={\rm{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 5 }$$

which corresponds to the total energy density and current density, respectively. Here ρ is the mass density; ${\boldsymbol{u}}=({v}_{x},{v}_{y},{v}_{z})$ are the x, y, and z-components of velocity; p is the thermal pressure; ${\boldsymbol{B}}$ denotes the magnetic field; t is time; $\gamma =5/3$ is the adiabatic index; and η is the magnetic resistivity.

Three independent parameters, an initial mean density ${\rho }_{0}$, a characteristic length L, and a characteristic plasma speed ${v}_{0}=\delta b/\sqrt{4\pi {\rho }_{0}}$, with $\delta b=\langle {b}^{2}{\rangle }^{1/2}$, are used to normalize the MHD equations. Other variables are normalized by their combinations. The dimensionless numbers appearing in the equations are the Mach number ${M}_{s}={v}_{0}/{c}_{s}$, where ${c}_{s}=\sqrt{\gamma {p}_{0}/{\rho }_{0}}$ is the sound speed, and the magnetic Reynolds number ${R}_{m}={v}_{0}L/\eta$. Here we take M_s to be 0.5, consistent with the solar wind observations at 1 AU, and R_m to be 1000, which is limited by the available spatial resolution. The uniform guide field B₀ is two times that of the fluctuating field $| \delta {\boldsymbol{b}}|$.

We consider periodic boundary conditions in a cube with a side length of $2\pi L$ and a resolution defined by the number of grid points, which is 512³, and run a simulation from an initial state with kinetic and magnetic energy per unit mass $\langle {v}^{2}\rangle =\langle {b}^{2}\rangle =1$. The fluctuations initially populate an annulus in the Fourier k-space such that $1\leqslant k\leqslant 8$, with constant amplitude and random phases (Matthaeus et al. 1996; Dmitruk et al. 2004). The initial normalized cross helicity is set to be 0.9 such that the primordial fluctuations are highly Alfvénic. The initial density and thermal pressure are set to be uniform.

To solve the equations, we employ a splitting-based finite-volume numerical scheme. The fluid part is solved by the Godunov-type central scheme and the magnetic part by the constrained transport approach, in conjunction with the method called second-order Monotone Upstream Schemes for Conservation Laws for reconstruction, and with the approximate Riemann solvers of Harten–Lax–van Leer for calculation of the numerical fluxes (Feng et al. 2011). The explicit second-order Runge–Kutta stepping with total variation diminishing is applied in the time integration.

Figure 1 shows the distribution of the current density in the z-direction J_z in the x–z (left) and x–y  (right) planes at t = 1.25. The arrows superposed on the images are the projections of the magnetic field vectors in the x–z (left) and the x–y planes (right). The black ellipses mark the region where the identified RD is formed. From this figure, we can see that a large-scale background magnetic field is clearly present in the z-direction. As a result of the well-known anisotropic behavior of magnetic field fluctuations in an MHD with an imposed uniform guide field (Matthaeus et al. 1996), current density structures preferentially align along the guide field direction, as shown in the left panel, and become much more varying in the perpendicular cross section, as shown in the right panel. Also, J_z appears to be large in magnitude at the location where the identified RD is formed.

Zoom In Zoom Out Reset image size

In order to detect the RD, we first seek the regions with large normalized partial variance of increments (${\rm{PVI}}$) of the magnetic field vector. ${\rm{PVI}}$ in 3D space is defined as

$$\begin{eqnarray*}&&{\rm{PVI}}(x,y,z)\\ &&=\;\left(\begin{array}{l}| {\boldsymbol{B}}(x+{\rm{\Delta }}x,y+{\rm{\Delta }}y,z+{\rm{\Delta }}z)\\ -{\boldsymbol{B}}(x-{\rm{\Delta }}x,y-{\rm{\Delta }}y,z-{\rm{\Delta }}z)| \end{array}\right)/\\ &&\left(\sqrt{\begin{array}{l}\langle | {\boldsymbol{B}}(x+{\rm{\Delta }}x,y+{\rm{\Delta }}y,z+{\rm{\Delta }}z)\\ -{\boldsymbol{B}}(x-{\rm{\Delta }}x,y-{\rm{\Delta }}y,z-{\rm{\Delta }}z){| }^{2}\rangle \end{array}}\right),\end{eqnarray*}$$

where ${\rm{\Delta }}x,{\rm{\Delta }}y$, and ${\rm{\Delta }}z$ are the grid increment in the x-, y-, and z-directions, respectively. We first sample the magnetic field ${\boldsymbol{B}}$, plasma velocity ${\boldsymbol{v}}$, density ρ, and temperature T along a linear path through the region with ${\rm{PVI}}\gt 2$, and then perform along that path the minimum variance analysis (MVA; Sonnerup & Cahill 1967) using the magnetic field data to find the maximum variance direction (L), intermediate variance direction (M), minimum variance direction (N), and their corresponding eigenvalues ${\lambda }_{1}$, ${\lambda }_{2}$, and ${\lambda }_{3}$. Finally, we check the 3D plasma and field structure of the possible events.

Figure 2 shows the magnetic field and plasma parameters along the sampling interval (s) as well as B_L–B_M hodogram for an identified RD. In this figure, the magnetic field has been converted into the Alfvén velocity ${{\boldsymbol{V}}}_{{\rm{A}}}$using ${{\boldsymbol{V}}}_{{\rm{A}}}={\boldsymbol{B}}/\sqrt{\rho }$ . We can see that the sampling series of  ${\rm{PVI}}$ is rather close to 0, except near the center. The large jumps of the Alfvén velocity and plasma velocity mainly occur along the L direction, with V_L jumping from −1.0 to 1.0. Their N-components, V_N,  are nearly constant and are equal to 1.6, along with a nearly negligible jump of $| {{\boldsymbol{V}}}_{{\rm{A}}}|$, the magnitude of Alfvén speed. Also, the Alfvén speed and plasma speed are nearly in agreement over the entire interval, including the jump across the RD itself. The positive correlation between them implies that the RD propagation is anti-parallel to the magnetic field. The density ρ, temperature T, and total pressure P_t (which consists of the summed magnetic pressure and thermal pressure) all exhibit relatively constant traces throughout the interval. To be noted, the jump of V_L and the large value of V_N, together with the relatively slight change of P_t as well as ρ, corroborate that this event is an RD.

Zoom In Zoom Out Reset image size

The B_L–B_M hodogram shows that the phase-steepened rotation of the Alfvén wave has an "arc-like" shape, indicating that the simulated Alfvén wave is nearly an arc-polarized phase-steepened Alfvén wave, consistent with observations by Tsurutani et al. (1994, 1996, 1997), and Tsurutani & Ho (1999). The perturbation vector starts at the upper-left corner herein, and rotates to the upper right corner. The resulting trajectory is approximately an arc with its center at the zero point.

Figure 3 exhibits the 3D structure of the identified RD. The green lines in the left panel denote magnetic field lines and the yellow arrows in the right panel are plasma velocity vectors, which are converted into the de Hoffman–Teller (HT) frame. The light gray surface is the isosurface where ${\rm{PVI}}=4$, showing the RD, and red, green, and blue arrows displaying the x-, y-, and z-directions, respectively. This figure shows that in the HT frame, magnetic field and plasma velocity have evident normal components across the RD, and they both rotate by a certain angle. Also, the plasma velocity is field-aligned on both sides of the discontinuity, in accordance with the Alfvénic nature of RD. The identified RD appears as a thin surface and its normal inclines to the z-axis.

Zoom In Zoom Out Reset image size

To see how the identified RD is formed, Figure 4 shows the isosurfaces of B_x (Left) and V_x (Right) at different moments in time. These isosurfaces are associated with the Alfvén wave as shown below. ${\rm{PVI}}=3$ is used here to exhibit the identified RD at these four moments, and is shown as the light gray surface. From this figure, it is notable that the mutual approaching of the two isosurfaces of B_x (V_x), which describes the steepening of the Alfvén wave, leads to the formation of an RD. At t = 0.95, the isosurface with ${B}_{x}=-0.55$ (${V}_{x}=-0.45$) is positioned relatively away from the isosurface with B_x = 0.75 (V_x = 0.45), and the transition of B_x (V_x) from ${B}_{x}=-0.55$ (${V}_{x}=-0.45$) to B_x = 0.75 (V_x = 0.45) is gentle. The isosurface where ${\rm{PVI}}=3$ is small and thin. At t = 1.10, the isosurfaces with ${B}_{x}=-0.55$ (${V}_{x}=-0.45$) and B_x = 0.75 (V_x = 0.45) approach each other, and the transition between the two isosurfaces of B_x (V_x) becomes steep. Correspondingly, the isosurface where ${\rm{PVI}}=3$ grows. Afterward, that is, at t = 1.25, the two isosurfaces of B_x (V_x) are close enough to each other, and the transition between the two isosurfaces of B_x (V_x) becomes steeper. The isosurface where PVI = 3 grows larger and thicker, and the identified RD is fully grown. However, at t = 1.40, the isosurfaces with ${B}_{x}=-0.55$ (${V}_{x}=-0.45$) are far away from the isosurface with B_x = 0.75 (V_x = 0.45), and the transition between the two isosurfaces of B_x (V_x) becomes gentle again. As a result, the isosurface where ${\rm{PVI}}=3$ becomes small and thin. The identified RD starts to collapse.

Zoom In Zoom Out Reset image size

To understand the type of waves before the RD is formed and to investigate the process of the evolution of the isosurfaces mentioned above, Figure 5 exhibits the distribution of B_x (the first row), V_x (the second row), and V_A (the third row) in the neighborhood of the identified RD in the plane x–z at different moments in time. Superposed by arrows are the projections of the magnetic field (the first row), the velocity field in the HT frame (the second row), and negative Alfvén speed field (the third row) vectors in the x–z plane. The black ellipses mark the position where the identified RD is formed. Near the ellipses the magnitudes of B_x and V_x are almost the same (they are set to have the same color scales so that the same color stands for the same value), and the directions of the in-plane projection of the ${\boldsymbol{B}}$ and ${\boldsymbol{V}}$ vectors are almost identical. Hence in this vicinity we have ${\boldsymbol{V}}={\boldsymbol{B}}/\sqrt{{\rho }_{0}}$ (we recall that ${\rho }_{0}\approx 1$), which agrees with the polarity relations of an Alfvén wave. In other words, the RD is detected in a neighboring Alfvénic environment that apparently favors RD formation.

Zoom In Zoom Out Reset image size

Therefore, V_A can be regarded as the propagation velocity of the structures relative to the location where the RD forms. Hence it is significant to trace the changes of the Alfvén speed. In the third row in Figure 5 where the evolution of the Alfvén speed is shown, we can see that at t = 0.95, there is a difference of Alfvén speed across the black ellipse, which makes the layers with negative B_x (blue in Figure 4), propagate faster than that with positive B_x (green in Figure 4). At t = 1.10, there is the evidence for an approaching and squeezing of these two layers. The difference of Alfvén speed there remains. This will drive these two layers further to approach each other. At t = 1.25, their transition becomes sharp, and the difference of Alfvén speed nearly disappears. This status continues until t = 1.40, when the transition becomes gentle as a result of the faster propagation of the layers with positive B_x than that with negative B_x. It is obvious that the difference of Alfvén speed makes the Alfvén wave steepen, a process which forms the identified RD.

In this study, we use a simulation of the decaying compressive 3D MHD turbulence with an imposed uniform guide field as a test case to explore the formation of RDs in MHD turbulence. Motivated by solar wind observation at 1 AU, we consider a moderate fluctuation amplitude corresponding to $\delta b/{B}_{0}=0.5$ and high Alfvénic correlations, with the normalized cross helicity initially of 0.9. A case study is thus conducted to illustrate the origin of RDs in MHD turbulence.

The numerical simulation shows the well-known anisotropic behavior of the turbulent MHD field, with the current density structures preferentially aligning along the guide field direction and scattering in the perpendicular plane. To detect the RDs in this simulated magnetofluid, we first seek the regions with large ${\rm{PVI}}$, then conduct an MVA by sampling the parameters along a linear path, and finally check the 3D structure of the possible RD events. One clear RD is identified with sharp field rotations and little variations of the magnetic field intensity or density. At the same time, in the HT frame, the plasma velocity is nearly in agreement with the Alfvén speed, and is field-aligned on both sides of the discontinuity, satisfying the Walen relation that expresses the Alfvénic nature of an RD. The normal direction of the identified RD inclines to the z-axis, and propagates anti-parallel to the guide field.

The comprehensive information obtained by the simulation of the magnetic field and plasma parameters associated with the RD implies that the RD is produced by the steepening of the moderate-amplitude Alfvén wave with nonuniform Alfvén speed in the ambient turbulence. Before the RD is formed, the layers with negative B_x smoothly transit to the layers with positive B_x. However, there is a difference of Alfvén speed across them, which makes the layers with negative B_x chase after their counterpart with positive B_x. As they are driven by the neighboring turbulence to approach and squeeze each other, the transition between them undergoes further steepening until the difference of Alfvén speed nearly disappears. At the same time, the identified RD is formed.

In this work we investigated the formation of an RD due to the nonlinear steepening of an Alfvén wave. It should be mentioned that magnetic decreases (MD)s are associated with the phase-steepened edges of the Alfvén waves (Tsurutani et al. 2002a, 2002b, 2003, 2005a, 2005b, 2009, 2011). These MDs are referred to as the dissipation of the phase-steepened Alfvén waves due to the fact that the proton-perpendicular temperature is higher inside the MDs than outside, and the wave observations at both ion and electron scales indicate that instabilities are driven by this heating process. As argued by Tsurutani et al. (2007), the existence of MDs at the location of the RDs will confuse the situation for determining what type of discontinuity might be occurring. Although these features of the Alfvén waves are beyond the present work, future work may be devoted to this direction to see the formation of MDs, as well as their influences.

Also, the phase-steepened Alfvén waves are argued to be intermediate shocks, which are rapidly evolving spatially (Tsurutani et al. 2005a, 2005b). This fits together with the idea that these Alfvénic waves have been freshly created and are spherical in nature. Such further exploration may involve simulations considering the interaction between the phase-steepened Alfvén wave with the high density MD. Will the MD damp out the intermediate shock? How will slow shocks form? There are many interesting topics here.

Another consequence associated with the wave phase-steepening is a spreading of the wave spectral power (Tsurutani et al. 2003). The phase-steepened edges form high-frequency power. The remaining trailing edges are enhancing low-frequency power. As indicated from the results presented by Tsurutani et al. (2005b), perhaps a substantial part of interplanetary turbulence is due to this physical process. It will be interesting to see whether the wave spectral power in this simulation of decaying MHD turbulence has such a spreading.

Certainly, many RDs are generated in the simulation. Statistical studies could be done to see whether or not all RDs are produced by the nonlinear steepening of an Alfvén wave. Meanwhile, TDs are also formed. It can be inspiring to investigate their generation mechanism in MHD turbulence. Furthermore, the parameters we adopted in the simulation (e.g. Mach number, cross helicity, plasma β) may influence the forming of RDs or TDs. In the future, we plan to conduct a parameter study to investigate the possible effects induced by parameter variation on the generation of discontinuities in MHD turbulence.

This work is supported by NSFC under contract No. 41304133, 41474147, 41231069, 41222032, 41174148, 41421003, and 41204105. The numerical calculations have been completed on a computing system at Peking University. Figure 3 was partly produced by VAPOR (Clyne et al. 2007). L.P.Y., J.S.H., C.Y.T., E.M., and X.W. are also members of the ISSI/ISSI-BJ International Team-304.