Intermittent turbulence in the heliosheath and the magnetosheath plasmas based on Voyager and THEMIS data

Turbulence is complex behavior that is ubiquitous in space, including the environments of the heliosphere and the magnetosphere. Our studies on solar wind turbulence including the heliosheath, and even at the heliospheric boundaries, also beyond the ecliptic plane, have shown that turbulence is intermittent in the entire heliosphere. As is known, turbulence in space plasmas often exhibits substantial deviations from normal Gaussian distributions. Therefore, we analyze the fluctuations of plasma and magnetic field parameters also in the magnetosheath behind the Earth's bow shock. Based on THEMIS observations, we have already suggested that turbulence behind the quasi-perpendicular shock is more intermittent with larger kurtosis than that behind the quasi-parallel shocks. Following this study, we would like to present a detailed analysis of intermittent anisotropic turbulence in the magnetosheath depending on various characteristics of plasma behind the bow shock and now also near the magnetopause. In particular, for very high Alfvénic Mach numbers and high plasma beta we have clear non-Gaussian statistics in the directions perpendicular to the magnetic field. On the other hand, for directions parallel to this field the kurtosis is small and the plasma is close to equilibrium. However, the level of intermittency for the outgoing fluctuations seems to be similar to that for the ingoing fluctuations, which is consistent with approximate equipartition of energy between the oppositely propagating Alfvén waves. We hope that the difference in characteristic behavior of these fluctuations in various regions of space plasmas can help to detect some complex structures in space missions in the near future.

1 Introduction

Turbulence is complex behavior that is ubiquitous in space, including the solar wind, interplanetary and interstellar media, as well as planetary and interstellar shocks (e.g., Bruno and Carbone, 2016). These shocks are usually collisionless and processes responsible for the plasma are substantially different from ordinary gases; see, e.g., Kivelson and Russell (1995) and Burgess and Scholer (2015). That is, the necessary coupling in plasma is usually provided by nonlinear structures at various scales, possibly exhibiting fractal or multifractal self-similarity properties (e.g., Burlaga, 1995; Macek, 2006). In addition, dissipation (so-called quasi-viscosity) could often result from wave damping or other processes related to electric current structures. The mechanism of complexity of space and astrophysical plasmas is still a challenge to turbulence problems (Chang, 2015).

2 Multifractal model

In our view, we should still rely on phenomenological models of intermittent turbulence, which can grasp multiplicative processes leading to complex behavior of the plasma in a simple way. As we have often argued (e.g., Macek, 2006, 2007; Macek and Wawrzaszek, 2009), the most useful concept for such a phenomenological study is a topological object, namely the generalized two-scale weighted Cantor set – an example of multifractals – as described, for example, by Falconer (1990). The turbulence model based on this set is sketched here in Fig. 1, as taken from Macek (2007). We see that at each step of construction of the generalized Cantor set, one needs to specify two scales l₁ and l₂ ($l_{\mathrm{1}}+l_{\mathrm{2}}\le \mathrm{1}$) associated with probability measures p and 1−p. In fact, fractals and multifractals could be considered a convenient mathematical language useful for understanding dynamics of turbulence, as already postulated by Mandelbrot (1982). In fact, in this review we will provide some arguments that this surprisingly simple mathematical rule provides a very efficient tool for phenomenological analysis of complex turbulent media.

Figure 1 Two-scale weighted Cantor set model for asymmetric solar wind turbulence (Macek, 2007).

Figure 2 The singularity multifractal spectrum f(α) vs. the singularity strength α with some general properties: (1) the maximum value of f(α) is D₀; (2) f(D₁)=D₁; and (3) the line joining the origin to the point on the f(α) curve where α=D₁ is tangent to the curve, as taken from Ott (1993).

Moreover, for the two-scale weighted Cantor set model, the singularity multifractal spectrum shown in Fig. 2 can easily be calculated (e.g., Ott, 1993). In particular, the width of this universal function, Δ, is obtained analytically by the following equation:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }\equiv {\mathit{\alpha }}_{\text{max}}-{\mathit{\alpha }}_{\text{min}}=D_{-\mathrm{\infty }}-D_{\mathrm{\infty }}=\left|{\displaystyle \frac{\log (\mathrm{1}-p)}{\log l_{\mathrm{2}}}}-{\displaystyle \frac{\log \left(p\right)}{\log l_{\mathrm{1}}}}\right|.\end{array}$$

Naturally, this quantity Δ is just the difference between the maximum and minimum dimensions related to the regions in the phase space with the least dense and most dense probability densities, and hence it has been proposed by Macek (2007) and Macek and Wawrzaszek (2009) as a degree of multifractality. Moreover, since this parameter Δ exhibits a deviation from a strict self-similarity, it can also be used as a degree of intermittency, as explained in Frisch (1995, Chapter 8). One can expect that the solar wind Δ will reveal various nonlinear phenomena, including nonlinear pressure pulses related to magnetosonic waves, as argued by Burlaga et al. (2003, 2007).

The other parameter A describing the multifractal scaling is the measure of asymmetry of the spectrum as defined by Macek and Wawrzaszek (2009):

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }A\equiv {\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}-{\mathit{\alpha }}_{\text{min}}}{{\mathit{\alpha }}_{\text{max}}-{\mathit{\alpha }}_{\mathrm{0}}}},\end{array}$$

where α=α₀ is the point at which the spectrum has its maximum, f(α₀)=1. In particular, in a simpler case when A=1 ($l_{\mathrm{1}}=l_{\mathrm{2}}=\mathrm{0.5}$), the one-scale p model is recovered (e.g., Meneveau and Sreenivasan, 1987), and for a monofractal the function in Fig. 2 is reduced to a point.

In principle, for experimental time series one can recover the multifractal spectrum and fit to either the well-known p model or the more general two-scale weighted Cantor set model. For Voyager data this can be done in the following way. That is, the generalized multifractal measures p(l) depending on scale l can be constructed using magnetic field strength fluctuations (Burlaga, 1995). Normalizing a time series of daily averages B(t_i), where $i=\mathrm{1},\mathrm{\dots },N={\mathrm{2}}^n$ for $j={\mathrm{2}}^{n-k}$, $k=\mathrm{0},\mathrm{1},\mathrm{\dots },n$,

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }p(x_j,l)\equiv {\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}+(j-\mathrm{1})\mathrm{\Delta }t}^{j\mathrm{\Delta }t}B\left(t_i\right)=p_j\left(l\right)\end{array}$$

is calculated with the successive average values 〈B(t_i,Δt)〉 of B(t_i) between t_i and t_i+Δt, for each Δt=2^k (Macek et al., 2011, 2012). When time series are obtained onboard spacecraft, it is usually possible to relate time dependence to space dependence by using the Taylor (1938) hypothesis. Because the average solar wind speed v_sw is much greater that the velocity of the space probe, we can argue that p(x_j,l) can be regarded as the probability that at a position x=v_swt, at time t, a given magnetic flux will be transferred to a spatial scale l=v_swΔt.

In this way Burlaga (1995) has shown that in the inertial range the average value of the qth moment of B at various scales l scales as

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\langle B^q\left(l\right)\rangle \sim l^{\mathit{\gamma }\left(q\right)},\end{array}$$

where the exponent γ is related to the generalized dimension, $\mathit{\gamma }\left(q\right)=(q-\mathrm{1})(D_q-\mathrm{1})$. If for a certain range of spatial scales l corresponding to a given interval of time Δt we have a straight line on a logarithmic scale using these slopes for each real q, the values of D_q can be determined with Eq. (4).

Alternatively, as explained by Macek and Wawrzaszek (2009), the multifractal function f(α) vs. scaling index α shown in Fig. 2, which exhibits universality of the multifractal scaling behavior, can be obtained using the Legendre transformation. It is worth noting, however, that we obtain this multifractal universal function directly from the slopes for a given scale range using this direct method in various situations (see Macek and Wawrzaszek, 2009; Macek et al., 2011, 2012, 2014).

3 Heliosheath turbulence

The schematic of the heliospheric boundaries is shown in Fig. 3. Voyager 1 entered the heliosheath after crossing the termination heliospheric shock at 94 AU in 2004, while Voyager 2 crossed this shock at 84 AU in 2007. It is generally accepted that after crossing the heliopause in 2012, the last boundary separating the heliosphere from the nearby interstellar medium, the Voyager 1 has ultimately left the heliosphere, while the crossing of the heliopause by Voyager 2 is expected in the very near future.

Figure 3 Schematic of the heliospheric boundaries (credit: NASA/Walt Feimer).

3.1 Heliosheath data

The main aim of our Voyager studies is to look at the measure of multifractal scaling in the heliosheath. Because in the distant heliosphere the magnetic fields have mainly azimuthal components, one can use the magnitude of the magnetic fields $\left|\mathit{B}\right|$ to estimate the probability measures and using straight lines according to Eq. (4) in a certain scale range, as with those seen in Figs. 1 and 2 of the paper by Macek et al. (2014), and in this way we can calculate the multifractal singularity spectrum. The results using the data gathered onboard both the Voyager 1 and 2 spacecraft immersed in the heliosheath are presented in Fig. 4, case (a) at 94–97 AU for the year 2005 and (c) at 105–107 AU for the year 2008 for Voyager 1, and Voyager 2 in case (b) at 85–88 AU for the year 2008 and (d) at 88–90 AU for the year 2009, respectively (Macek et al., 2012, Fig. 5).

Figure 4 The singularity spectrum f(α) as a function of a singularity strength α. The values are calculated for the weighted two-scale (continuous lines) model and the usual one-scale (dashed lines) p model with the parameters fitted using the magnetic fields (diamonds) measured by Voyager 1 in the heliosheath at various heliocentric distances of (a) 94–97 AU and (c) 105–107 AU, and by Voyager 2 at (b) 85–88 AU and (d) 88–90 AU, respectively, taken from Macek et al. (2012).

Table 1 The measure of chi-square fitting for the weighted two-scale and one-scale p models.

Figure 5 The parameter Δ quantifying multifractality in the heliosphere as a function of the distances from the Sun together with a periodic function shown by dotted (dashed) lines during different phases of the solar cycle (SSN is the sunspot number) in the heliosphere (Macek et al., 2011) and in the heliosheath (Macek et al., 2014). The heliospheric termination shock (TS) and the heliopause (HP) crossings by Voyager 1 are indicated.

Figure 6 The degree of multifractality Δ and asymmetry A in the heliosphere as a function of the distances from the Sun. The termination shock (TS) and the heliopause (HP) crossings by Voyager 1 are also indicated (cf. Macek, 2012).

It seems that the two-scale weighted Cantor set model fits the data better than the classical p model. To support this result in a more quantitative way, we have used the weighted χ², consisting of a sum of squares of differences between the spectrum obtained from data and the model, each normalized to unit variance (e.g., Press et al., 1992). This measure of fit quality in the two cases is shown in Table 1. We see that the values obtained for the two-scale model are at least one order in magnitude smaller than that for the standard p model. This means that the generalized Cantor set model is in fact substantially better.

We have also calculated the degree of multifractality Δ, as given in Eq. (1), for Voyager 1 in the heliosheath depending on the heliospheric distances during different phases of the solar cycle, the minimum (MIN), maximum (MAX), declining (DEC), and rising (RIS) phases, which is now demonstrated in Fig. 5 (left panel). We clearly see that this multifractality measure obtained for the scale range of Δt from 2 to 16 days decreases steadily with the heliospheric distance and is modulated by the solar activity following the sunspot numbers (SSN, indicated at the bottom panel), as taken from Macek et al. (2011, Fig. 2). In the right panel we show Δ calculated in the heliosheath for two various scaling ranges from 2 to 16 and from 4 to 32 days (cf. Macek et al., 2014). The crossings of the termination shock (TS) and the heliopause (HP) by Voyager 1 are indicated by vertical dashed lines. We see that in the heliosheath the degree of multifractality basically still follows the periodic dependence fitted inside the heliosphere (Macek et al., 2011, 2012). It is worth noting that after crossing the heliopause at ∼122 AU, the value of Δ suddenly drops to zero, and nonmultifractal (nonintermittent) smoothly varying magnetic fields are observed by Voyager 1, as indicated in the right panel of Fig. 5; see Macek et al. (2014). This would mean that the entire heliosphere with turbulent plasma inside is immersed in a relatively quiet ambient very local interstellar medium.

Naturally, the multifractal spectrum can be related to nonlinear Alfvén waves, associated with discontinuities, or mirror mode structures due to some plasma instabilities, or possibly current sheets (Borovsky, 2010; Tsurutani et al., 2011a, b) generated upstream of the termination shock, as discussed in our previous paper (Macek and Wawrzaszek, 2013). In this way we have applied the multifractal model (Macek, 2007; Macek and Szczepaniak, 2008) to solar wind turbulence in the entire heliosphere (Szczepaniak and Macek, 2008; Macek and Wawrzaszek, 2009; Macek et al., 2011, 2012), also beyond the ecliptic plane (Wawrzaszek and Macek, 2010; Wawrzaszek et al., 2015), and even at the heliospheric boundaries (Burlaga et al., 2013; Macek et al., 2014), and have shown that turbulence could often be intermittent. By the way, it would be difficult to argue that there is an asymmetry in these spectra for the Voyager 1 data, but there are some deviations from the symmetric spectrum for Voyager 2. In summary, the values of the degree of intermittency calculated from our two-scale weighted Cantor set model are presented in Fig. 6 (cf. Macek, 2012).

As is known, turbulence in space and astrophysical plasmas exhibits deviations from normal distributions, and these higher moments are often considered signatures of intermittency. In particular, kurtosis – the fourth moment of the probability density function – is often used as a measure of intermittency (Bruno et al., 2003; Bruno and Carbone, 2013).

Figure 7 Schematic of the THEMIS mission (credit: NASA/ESA).

4 Magnetosheath turbulence

Naturally, nonlinear structures responsible for turbulence have already been identified in planetary environments, in the solar wind, and also in the magnetosheath (e.g., Alexandrova, 2008). In particular, the magnetic fluctuations using Wind (Lion et al., 2016) and Cluster multi-spacecraft have been analyzed at ion scales (Yordanova et al., 2008; Roberts et al., 2016; Perrone et al., 2016, 2017). In addition, some results on very high-resolution data on electron scales have recently been provided by the Magnetospheric Multiscale (MMS) mission (Yordanova et al., 2016; Chasapis et al., 2017). Moreover, on the basis of kinetic simulations by Karimabadi et al. (2014), one can suggest some interesting relationships of turbulent processes near shocks with reconnection processes. But in spite of progress in MHD simulations, including Hall effects, the physical mechanisms of turbulent behavior are still not sufficiently clear.

Table 2 List of selected interval samples (mm.dd.hh.MM).

Various space missions provide unique observational data, which help to understand phenomena in our environment in space. In particular, the THEMIS mission was launched by NASA in 2007 in order to resolve macroscale phenomena occurring during substorms (Sibeck and Angelopoulos, 2008), as schematically presented in Fig. 7. In addition, for the first time THEMIS data were used for analysis of turbulence at the terrestrial bow shock. That is, we have suggested that turbulence behind the quasi-perpendicular shock is more intermittent with larger kurtosis than that behind the quasi-parallel shocks (Macek et al., 2015).

In this review paper, besides turbulence in the heliosheath, as has already been discussed in Sect. 3, now in Sect. 4 we continue our study in the entire magnetosheath also near the magnetopause. However, since it would be difficult to obtain the full multifractal spectrum using the THEMIS data, at present we only examine how the degree of multifractality resulting in deviation from the normal distribution, which is also a level of intermittency, depends on the characteristics of the solar wind and magnetospheric plasmas (Macek et al., 2017). The data under study are briefly described in Sect. 4.1. In Sect. 4.2 we present the results of our analysis, showing in particular that at high Alfvénic Mach numbers turbulence becomes clearly intermittent. The importance of this intermittent behavior for space plasmas is underlined in Sect. 5.

4.1 Magnetosheath data

We analyze various time samples acquired during the long period between 2008 and 2015 from the THEMIS mission consisting of a quintet (A, B, C, D, and E) of space probes (Sibeck and Angelopoulos, 2008), as listed in Table 2. We have selected the following 24 intervals in the magnetosheath (without any evident large-scale static plasma structures): 11 samples measured after crossing the bow shock, denoted by BS, and 13 samples obtained before leaving the magnetosheath, i.e., near the magnetopause, denoted by MP. The time resolution here is 3 s and these samples taken in the Geocentric Solar Ecliptic (GSE) reference system are all (except for no. 11) longer than 4 h. Naturally, the length of each sample depends on the orbit of a particular probe immersed in the magnetosheath during some periods of time. Please note that the timescales in the magnetosheath are much shorter than that in the heliosheath.

Various characteristic plasma parameters, namely the Alfvén Mach number, M_A, the plasma parameter beta, β, and the magnetosonic Mach number, M_ms, are calculated in the solar wind upstream: first before crossing the bow shock (before entering the magnetosheath) and next in the magnetosphere (before crossing the magnetopause). The plasma β is the ratio of the thermal pressure p to the magnetic pressure B²∕(2μ₀ρ), where ρ=mN is the mass density for ions of mass m and the number density N (μ₀ denotes the permeability of free space).

Figure 8 Alfvén Mach number (a), total plasma beta (b), and magnetosonic Mach number (c) near the bow shock (BS, red circles) and the magnetopause (MP, white triangles).

All three of these plasma parameters vs. sample number are depicted in Fig. 8. We see that the Alfvén Mach numbers can vary substantially with the limiting value of about 25 ($\mathrm{5}\le M_{\mathrm{A}}\le \mathrm{25}$), and that in most cases β is below 5 (only three cases are above 10). However, the magnetosonic Mach numbers are rather moderate: $\mathrm{3.6}\le M_{\text{ms}}\le \mathrm{7.5}$.

Figure 9 Kurtosis of the increments of the Elsässer vectors, δz^±, in the magnetosheath vs. the Alfvén Mach number for the components $\mathit{\delta }{z}_{\parallel }^{\pm }$ (parallel to the mean field vector B_o, panels a and b) and two components perpendicular to the field: (1) $\mathit{\delta }{z}_{⟂\mathrm{1}}^{\pm }$ (in the plane containing B_o and the X axis in the GSE system, panels c and d), and (2) $\mathit{\delta }{z}_{⟂\mathrm{2}}^{\pm }$ (perpendicular to the plane, panels e and f) near the bow shock (BS, red circles), with averages marked by a continuous line, and near the magnetopause (MP, white triangles), marked by dashed lines, as observed by THEMIS for samples listed in Table 2.

Figure 10 Kurtosis of the increments of the Elsässer vectors, δz^±, in the magnetosheath vs. the total plasma beta β for the components $\mathit{\delta }{z}_{\parallel }^{\pm }$ (parallel to the mean field vector B_o, panels a and b) and two components perpendicular to the field: (1) $\mathit{\delta }{z}_{⟂\mathrm{1}}^{\pm }$ (in the plane containing B_o and the X axis in the GSE system, panels c and d), and (2) $\mathit{\delta }{z}_{⟂\mathrm{2}}^{\pm }$ (perpendicular to the plane, panels e and f) near the bow shock (BS, red circles) and near the magnetopause (MP, white triangles), as observed by THEMIS for samples listed in Table 2.

Figure 11 Kurtosis of the increments of the Elsässer vectors, δz^±, in the magnetosheath vs. the magnetosonic Mach number for the components $\mathit{\delta }{z}_{\parallel }^{\pm }$ (parallel to the mean field vector B_o, panels a and b) and two components perpendicular to the field: (1) $\mathit{\delta }{z}_{⟂\mathrm{1}}^{\pm }$ (in the plane containing B_o and X axis in the GSE system, panels c and d), and (2) $\mathit{\delta }{z}_{⟂\mathrm{2}}^{\pm }$ (perpendicular to the plane, panels e and f) near the bow shock (BS, red circles) and magnetopause (MP, white triangles), as observed by THEMIS for cases listed in Table 2.

4.2 Results for the magnetosheath

Using the values of plasma and magnetic fields shown in Figs. 1 and 2 of the paper by Macek et al. (2017), we can calculate the Elsässer variables, ${\mathit{z}}^{\pm }=\mathit{V}\pm {\mathit{V}}_{\mathrm{A}}$, where the characteristic Alfvénic velocity is given by ${\mathit{V}}_{\mathrm{A}}=\mathit{B}/({\mathit{\mu }}_{\mathrm{0}}\mathit{\rho }{)}^{\mathrm{1}/\mathrm{2}}$ (Elsasser, 1950). It is worth noting that the sign is taken here relative to the local average magnetic field B_o, which certainly depends on the timescale τ responsible for turbulence (Kiyani et al., 2013), as recently noted by Gerick et al. (2017). Because the time period during which this average background magnetic field is calculated, say dτ, should be substantially larger than the timescale of turbulence τ, we have taken d=10. By the way, in turbulence the dependence of statistical moments on spatial scales is often considered. For example, based on spacecraft measurements in the solar wind, one can estimate spatial scales by using the Taylor (1938) hypothesis (e.g., Macek and Wawrzaszek, 2009). However, in the magnetosheath the solar wind velocity is substantially reduced, and this approach can be somewhat less certain (Mangeney et al., 2006), especially for some plasma parameters (Perri et al., 2017). Therefore, it is better to analyze directly time samples obtained onboard several space probes, as is the case with the THEMIS mission.

Now, following our previous work on THEMIS data, the kurtosis of the increments of the various components of both Elsässer vectors z^±, $\mathit{\delta }{\mathit{z}}^{\pm }(t,\mathit{\tau })={\mathit{z}}^{\pm }(t+\mathit{\tau })-{\mathit{z}}^{\pm }\left(t\right)$, can be calculated for any given scale τ, taken in units of time resolution (Macek et al., 2015, Eq. 1). As is known, the Alfvénic increments perpendicular to the direction of B_o and the parallel compressive (slow-mode-like) increments should provide rather different contributions to the turbulent behavior of the solar wind plasma (e.g., Bruno et al., 2003; Oughton and Matthaeus, 2005). Therefore, we have performed our calculation in the Mean Field (MF) coordinate system, as described by Bruno and Carbone (2013). That is, the direction parallel to the local mean field B_o in the GSE system is taken along the versor $\widehat{\mathit{z}}={\mathit{B}}_{\mathrm{o}}/B_o$ of the new MF reference system (the symbol $\widehat{}$ is used for a unitary vector). This allows us to calculate the parallel components of both Elsässer vectors $\mathit{\delta }{z}_{\parallel }^{+}$ and $\mathit{\delta }{z}_{\parallel }^{-}$. Next, in order to obtain two other components perpendicular to the field B_o, $\mathit{\delta }{z}_{⟂\mathrm{1}}^{+}$, $\mathit{\delta }{z}_{⟂\mathrm{1}}^{-}$ and $\mathit{\delta }{z}_{⟂\mathrm{2}}^{+}$, $\mathit{\delta }{z}_{⟂\mathrm{2}}^{-}$, we take for the latter case the axis in the direction perpendicular to the plane containing the mean field B_o and the X axis in the GSE system (taken here as positive from the Sun), which is approximately consistent with the radial component of the mean solar wind velocity (V), that is, along $\widehat{\mathit{y}}=\widehat{\mathit{z}}\times \widehat{\mathit{x}}$. The remaining transverse components $\mathit{\delta }{z}_{⟂\mathrm{1}}^{\pm }$ are along $\widehat{\mathit{x}}=\widehat{\mathit{y}}\times \widehat{\mathit{z}}$ in this plane, which completes the right-handed orthogonal MF system.

The obtained values of kurtosis of the increments of the fluctuations of the Elsässer variables for the outgoing and ingoing Alfvénic fluctuations, respectively, z⁺ and z⁻, as observed by THEMIS in the magnetosheath near the bow shock (BS, red circles) and magnetopause (MP, white triangles) vs. the Alfvén Mach number, the total plasma beta β, and the magnetosonic Mach number, corresponding to Fig. 8, are presented in Figs. 9–11 for all 24 cases listed in Table 2. The departure of the probability density functions from normal distributions for the selected four cases corresponding to Figs. 1 and 2 of the paper by Macek et al. (2017), namely near the bow shock, cases (a) and (b), and near the magnetopause, cases (c) and (d), for a given timescale τ = 9 s, are illustrated in Figs. 12 and 13, respectively. The dependence of the kurtosis on the timescale τ is depicted in the corresponding Figs. 14 and 15.

Figure 12 The probability density functions (PDFs) of the increments of the parallel (white circles) and two perpendicular (black diamonds and squares) components of the Elsässer variables, δz⁺ (a', b') and δz⁻ (a”, b”), for a given timescale τ=9 s, near the bow shock (cases (a) and (b) in Table 2).

Figure 13 The probability density functions (PDFs) of the increments of the parallel (white circles) and perpendicular (black diamonds and squares) components of the Elsässer variables, δz⁺ (c', d') and δz⁻ (c”, d”), for a given timescale τ = 9 s, near the magnetopause (cases (c) and (d) in Table 2).

Figure 14 Kurtosis for the increments of the parallel (white circles) and two perpendicular (black diamonds and squares) components of the Elsässer variables, δz⁺ (a', b') and δz⁻ (a”, b”), as a function of timescale τ near the bow shock (cases (a) and (b) in Table 2).

Figure 15 Kurtosis for the increments of the parallel (white circles) and perpendicular (black diamonds and squares) components of the Elsässer variables, δz⁺ (c', d') and δz⁻ (c”, d”), as a function of timescale τ near the magnetopause (cases (c) and (d) in Table 2).

Figures 9–11 show kurtosis for the increments of all the components of the Elsässer vectors, for all the cases listed in Table 2, but for only one scale. Even though there is no very clear dependence on these plasma parameters, one can notice that the value of kurtosis often decreases with Alfvénic Mach number along the local magnetic field and sometimes increases in the perpendicular directions. We see from Fig. 9 that near the bow shock for the outgoing fluctuations kurtosis along the magnetic field $\mathit{\delta }{z}_{\parallel }^{+}$ somewhat decreases from 16.56 ± 0.06 at lower M_A (bin: $\mathrm{5}\le M_{\mathrm{A}}\le \mathrm{15}$) to 10.28 ± 0.06 at higher M_A (bin: $\mathrm{15}<M_{\mathrm{A}}\le \mathrm{25}$, even though we have only two points in this bin). But near the magnetopause, where we have more points in the letter bin, we can observe a clear significant decrease from 12.58 ± 0.05 to 7.25 ± 0.05. We have basically a similar behavior for the ingoing fluctuations, $\mathit{\delta }{z}_{\parallel }^{-}$: a decrease from 19.69 ± 0.06 to 16.86 ± 0.06 at the BS and a clear decrease from 17.66 ± 0.05 to 6.20 ± 0.05 at the MP. Here the standard deviations of kurtosis are rather small, about 0.06, as calculated according to Press et al. (1992). However, for the transverse components $\mathit{\delta }{z}_{⟂\mathrm{1},\mathrm{2}}^{\pm }$, kurtosis seems to be more scattered and often rather increasing with the Alfvénic Mach number, seeming not only anisotropic, but also seeming to be non-gyrotropic, with differences in two perpendicular components. We can see from Fig. 10 that $\mathit{\delta }{z}_{\parallel }^{\pm }$ decrease with plasma β, approaching normal distribution for high β, when the thermal pressure dominates the plasma behavior. But we do not see any clear regularity for the dependence of $\mathit{\delta }{z}_{⟂\mathrm{1},\mathrm{2}}^{\pm }$ on β. It also seems from Fig. 11 that the value of kurtosis is not very sensitive to the magnetosonic Mach number, but admittedly the range of this parameter considered in Table 2 is rather limited: $\mathrm{3.6}\le M_{\text{ms}}\le \mathrm{7.5}$.

Additionally, for the four clearly quasi-perpendicular cases (illustrated in Figs. 1 and 2 of the work by Macek et al., 2017, cases 3, 5, 15, and 16 listed in Table 2), the dependence of the parallel and perpendicular components of kurtosis on timescale τ for both the outgoing (z⁺) and ingoing fluctuations (z⁻) is now presented in Figs. 14 and 15, taken from (Macek et al., 2017). We can see that kurtosis behind the bow shock, Fig. 14, cases (a) and (b) in Table 2, could sometimes (for β∼1) be smaller than that near the magnetopause, Fig. 15, cases (c) and (d) in Table 2. We can generally notice only small differences between z⁺ and z⁻, and therefore the outgoing and ingoing fluctuations seem to be similar, which is roughly consistent with equipartition suggested by Tu et al. (1989). On the other hand, behind the bow shock for small plasma β∼1 (when the thermal pressure and the magnetic pressure are similar in the magnetized plasma), but with a moderate Alfvénic Mach number M_A≈9, Fig. 14a' and a”, case (a) in Table 2, we see only small kurtosis with approximately Gaussian normal distribution (i.e., close to equilibrium). For similar M_A≈12 and somewhat higher plasma β∼4, Fig. 14b' and b”, case (b) in Table 2, both parallel and perpendicular components of the Elsässer vectors are active. A similar behavior is also observed near the magnetopause, Fig. 15c' and c”, case (c) in Table 2. Finally, it is worth noting that for the highest value of M_A = 25 and β = 16.5, as illustrated in Fig. 15d' and d” (case (d) in Table 2), the perpendicular $\mathit{\delta }{z}_{⟂\mathrm{1},\mathrm{2}}^{\pm }$ components of fluctuations of Elsässer vectors are much larger than the parallel $\mathit{\delta }{z}_{\parallel }^{\pm }$ components. This exhibits a clear intermittent anisotropic turbulence with non-Gaussian probability distributions in transverse directions. On the other hand, the plasma along the local magnetic field is rather close to equilibrium.

Even though there is no clear regularity in Figs. 14 and 15 showing dependence of the kurtosis on scale τ, it seems that kurtosis near the bow shock (Fig. 14) is rather similar to that near the magnetopause (Fig. 15). That is, based on Figs. 9 and 10, it seems that near the bow shock (circles) the intermittency seems to decrease with the Alfvénic Mach number M_A and decrease with the plasma beta β near the magnetopause (triangles). We also see some difference between z⁺ and z⁻ in Fig. 14b' and b”, behind the bow shock (BS), and some scatters near the magnetopause (MP) for small scales (Fig. 15c”, for τ<30 s), and therefore we could consider the other cases in Table 2. In fact, generally speaking, we have verified that the level of intermittency for the outgoing fluctuations z⁺ is usually similar to that for the ingoing fluctuations z⁻, which exhibits approximate equipartition of energy between these oppositely propagating Alfvén waves.

5 Conclusions

Using our weighted two-scale Cantor set model, which is a convenient tool to investigate the asymmetry of the multifractal spectrum, we confirm the characteristic shape of the universal multifractal singularity spectrum. In fact, as seen in Fig. 4, f(α) is a downward concave function of scaling indices α. We show that the degree of multifractality for magnetic field fluctuations of the solar wind falls steadily with the distance from the Sun and seems to be modulated by the solar activity also in the heliosheath. Moreover, we have considered the multifractal spectra of fluctuations of the interplanetary magnetic field strength before and after crossing of the heliospheric termination shock by Voyager 1 and 2 near 94 and 84 AU from the Sun, respectively.

Further, we have provided important evidence that the large-scale magnetic field fluctuations reveal the multifractal structure not only in the outer heliosphere, but also in the entire heliosheath, even near the heliopause. Naturally, the evolution of the multifractal distributions should be related to some physical (MHD) models, as suggested by Burlaga et al. (2003, 2007). The driver of the multifractality in the heliosheath could be the solar variability on scales from hours to days, fast and slow streams or shock interactions, and other nonlinear structures discussed by Macek and Wawrzaszek (2013). In our view, any accurate physical model must reproduce the multifractal spectra. In particular, the observed nonmultifractal scaling after the heliopause crossing suggests nonintermittent behavior in the nearby interstellar medium, consistent with the smoothly varying interstellar magnetic field reported by Burlaga and Ness (2014). We have identified the scaling region of fluctuations of the interplanetary magnetic field.

In fact, using our two-scale model based on the weighted Cantor set, we have examined the universal multifractal spectra before and after crossing by Voyager 1: the termination shock at 94 AU and before crossing the heliopause at distances of about 122 AU from the Sun. Moreover, inside the heliosphere we observe the asymmetric spectrum, which becomes more symmetric in the heliosheath. We confirm that multifractality of magnetic field fluctuations embedded in the solar wind plasma for large scales decreases slowly with the heliospheric distance, demonstrating that this quantity is still modulated by the solar cycles further in the heliosheath, and even in the vicinity of the heliopause, possibly approaching a uniform nonintermittent behavior in the nearby interstellar medium. We propose this change in behavior as a signature of the expected crossing of the heliopause by Voyager 2 in the near future.

Regarding the magnetosheath, we have shown that turbulence for small scales is intermittent in the entire magnetosheath, in regions near the bow shock, and even near the magnetopause. In particular, we have found that near the magnetopause at very high Alfvénic Mach numbers M_A and high plasma β the probability density functions of compressive fluctuations parallel to the local average magnetic field should be nearly normal and close to equilibrium with small kurtosis, while in the transverse Alfvénic turbulence, resulting from nonlinear interactions, is non-gyrotropic with large kurtosis for the Elsässer variables. These fluctuations are more intermittent than that at the lower Alfvénic Mach numbers and plasma beta behind the bow shock. On the other hand, the level of intermittency for the outgoing fluctuations (z⁺) seems to be approximately similar to that for the ingoing fluctuations (z⁻). In view of the space investigation in the near future, including the THOR mission (e.g., Vaivads et al., 2016), we expect that the difference in characteristic behavior of these fluctuations in various regions of the magnetosheath will be able to help in identifying some new complex structures in space plasmas.