2. Classic theory of discontinuities

At every layer, the downstream and upstream physical quantities are linked by the fundamental conservation laws: mass, momentum, energy and magnetic flux (Landau & Lifshitz Reference Landau and Lifshitz1987). The simplest case occurs whenever the number of conservation laws is equal to the number of parameters characterizing the plasma state. When this condition is met, the possible downstream states are uniquely determined as a function of the upstream state, regardless of the (non-ideal) physics at play within the layer. In particular, it is possible to describe pressure variations without any closure equation. In this case, the jumps of all quantities are determined by a single scalar parameter (namely the ‘shock parameter’ in neutral gas).

We refer hereafter to the ‘classic theory of discontinuities’ (CTD) as for the theory corresponding to this condition, which is used both for neutral media and (magnetized) plasmas. The CTD is characterized by the following simplifying assumptions: a stationary layer, 1-D variations, and isotropic pressure on both sides. For plasmas, the additional assumption of an ideal Ohm's law on both sides is considered (Belmont et al. Reference Belmont, Rezeau, Riconda and Zaslavsky2019).

In CTD the conservation laws provide a system of jump equations between the upstream and downstream physical quantities, namely the Rankine–Hugoniot conditions in neutral media and generalized Rankine–Hugoniot conditions in plasmas. The sets of equations used to compute the linear modes in hydrodynamics and MHD are similar to this system of jump equations, simply because the hydrodynamics and MHD models rely on the same conservation laws as Rankine–Hugoniot and generalized Rankine–Hugoniot, respectively. A direct consequence is that many properties are shared by the solutions of the two types of systems: linear modes and discontinuities. For a neutral medium, the linear sound wave solution corresponds to the well-known sonic shock solution, while for a magnetized plasma, the two magnetosonic waves correspond to the two main types of MHD shocks: fast and slow. However, an additional discontinuity solution, the intermediate shock, has no linear counterpart. The intermediate shock presents a reversal of the tangential magnetic field through the discontinuity, which is not observed either in the fast or in the slow mode. Furthermore, a non-compressional solution exists in both types of systems, represented by the shear Alfvén mode for linear MHD, and by the ‘rotational discontinuity’ solution for the generalized Rankine–Hugoniot system.

Focusing on magnetized plasma physics, CTD leads to distinguish compressive and rotational discontinuities. An important feature of these solutions is that the compressional and rotational solutions are mutually exclusive: the shock solutions are purely compressional, without any rotation of the tangential magnetic field (this is called the ‘coplanarity property’), while the rotational discontinuity does imply such a rotation but without any variation of the magnetic field amplitude and without any compression of the particle density (figure 1). This distinction persists whatever the fluxes along the discontinuity normal, even when the normal components $u_n$ and $B_n$ of the velocity and the magnetic field are arbitrarily small. The only exception is the ‘tangential discontinuity’ when both normal fluxes are strictly zero. This solution would correspond, for the magnetopause, to the case without any connection between the solar wind and magnetosphere. It appears as a singular case since the tangential discontinuity, with $B_n=0$, is not the limit of any of the general solutions with $B_n \ne 0$. While the limit always implies two solutions, one purely rotational and the other purely compressional, the singular solution $B_n=0$ only provides one solution where the two characters can coexist.

Figure 1. Cartoon showing the different variations of $B$ between a rotational discontinuity (a,c) and a compressive one (b,d). The top panel shows in three dimensions the variation of $B$ inside the magnetopause plane; the bottom panel shows the hodogram in this tangential plane: a circular arc for the rotational discontinuity and a radial line for shocks.

In the solar wind, discontinuities are routinely observed and several authors have performed statistics for a long time to determine the proportion of the different kinds of discontinuities, mainly focusing on the tangential and rotational ones. They conclude that in most cases tangential discontinuities (i.e. with $B_n$ small enough to be barely measurable) are the most ubiquitous (see Colburn & Sonett Reference Colburn and Sonett1966 for a pioneering work in this domain and Neugebauer Reference Neugebauer2006; Paschmann et al. Reference Paschmann, Haaland, Sonnerup and Knetter2013; Liu et al. Reference Liu, Fu, Cao, Wang, He, Guo, Xu and Yu2022, and references therein, for more recent contributions). In these studies, rotational discontinuities are identified only when $B_n$ is large enough. However, many discontinuities present features that are typical of both rotational and tangential discontinuities and are classified as ‘either’ of the two. Extending these studies in the range of small $B_n$, where all discontinuities are not necessarily ‘tangential discontinuities’ in the CTD sense, requires the study of the quasi-tangential case.

3. The Earth's magnetopause

Thanks to in-situ observations, the Earth's magnetopause has a pivotal role in testing the discontinuity theories. Indeed, the Earth's magnetopause boundary exhibits, over its entire surface, both a rotation of the magnetic field (Sonnerup & Ledley Reference Sonnerup and Ledley1974) and a density variation (Otto Reference Otto2005) since it is the junction of two media, the magnetosheath and the magnetosphere where the magnetic field and the density are different (Dorville et al. Reference Dorville, Belmont, Rezeau, Grappin and Retinò2014).

As stated above, the usual paradigm is that the magnetopause is always a tangential discontinuity and that it becomes ‘open’ only exceptionally at a few points where the boundary departs from one-dimensionality due to magnetic reconnection. Does it mean that it justifies the very radical hypothesis of a magnetopause nearly completely impermeable to mass and magnetic flux, with strictly null $B_n$ and $u_n$ and quasi-independent plasmas on both sides (apart from the normal pressure equilibrium)? From a theoretical point of view, it is clear that the singular limit from $B_n \simeq 0$ to $B_n=0$ remains to be solved. From an experimental point of view, if the components $B_n$ and $u_n$ are known to be always very small, the observations can hardly distinguish between $B_n \simeq 0$ and $B_n=0$ because of the uncertainties, due to the fluctuations and the limited accuracy in determining the normal direction (Haaland et al. Reference Haaland, Sonnerup, Dunlop, Balogh, Georgescu, Hasegawa, Klecker, Paschmann, Puhl-Quinn, Rème, Vaith and Vaivads2004; Dorville et al. Reference Dorville, Haaland, Anekallu, Belmont and Rezeau2015b; Rezeau et al. Reference Rezeau, Belmont, Manuzzo, Aunai and Dargent2018).

The results of the present paper will question the above paradigm. We will show theoretically and experimentally that CTD fails at the magnetopause and that rotation and compression can actually coexist with finite $B_n$ and $u_n$, even in the 1-D case. Such a paradigm change may be reminiscent of a similar improvement in the theoretical modelling of the magnetotail in the 70's studies (Coppi, Laval & Pellat Reference Coppi, Laval and Pellat1966; Galeev Reference Galeev1979; Coroniti Reference Coroniti1980 and references therein). In that case the authors demonstrated that even a very weak component of the magnetic field across the current layer was sufficient to completely modify the stability properties of the plasma sheet, so that the finite value of $B_n$ had to be taken into account, contrary to the pioneer versions of the tearing instability theories.

4. The role of pressure

In CTD the separation between the compressional and rotational properties of the discontinuities comes from only two equations projected on the tangential plane. These equations are the momentum equation and the Faraday/Ohm's law, that read

(4.1)\begin{gather} \rho\frac{{\rm d} \boldsymbol{u}}{{\rm d} t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{P}_i +\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{P}_e=\boldsymbol{J}\times \boldsymbol{B}, \end{gather}

(4.2)\begin{gather}\boldsymbol{\nabla}\times \boldsymbol{E}={-}\frac{\partial \boldsymbol{B}}{\partial t}, \end{gather}

where

(4.3)\begin{equation} \boldsymbol{E}={-}\boldsymbol{u}\times\boldsymbol{B}+\frac{1}{ne} \boldsymbol{J}\times\boldsymbol{B}-\frac{1}{ne}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{P}_e, \end{equation}

where $\boldsymbol B$ is the magnetic field and $\boldsymbol u$ is the flow velocity in a reference frame where the layer is steady.

Considering 1-D gradients along the normal direction $\boldsymbol n$, neglecting the non-ideal terms in Ohm's law and integrating across the layer, these two equations, projected on the tangential plane, give

(4.4)\begin{gather} \rho_{2} u_{n2} \boldsymbol u_{t2}-B_{n2} \boldsymbol B_{t2}/\mu_0 = \rho_{1} u_{n1} \boldsymbol u_{t1}-B_{n1} \boldsymbol B_{t1}/\mu_0 , \end{gather}

(4.5)\begin{gather}B_{n2} \boldsymbol u_{t2}-u_{n2} \boldsymbol B_{t2}=B_{n1} \boldsymbol u_{t1}-u_{n1} \boldsymbol B_{t1} . \end{gather}

Due to the divergence free equation, the values $B_{n1}$ and $B_{n2}$ are equal and will be written as $B_{n}$ without index in the following. Similarly, $\rho _1 u_{n1}$ and $\rho _2 u_{n2}$ are equal because of the continuity equation and will be simply noted $\rho u_n$ in the following. Here, the indices $n$ and $t$ indicate the projection along the normal and in the tangential plane, respectively, while indices 1 and 2 indicate the two sides of the discontinuity. It is important to note that, in CTD, the pressure divergence terms do not appear in (4.4) because of the assumption in this theory that the pressure is isotropic on both sides so that their integration gives terms of the form $(p_2-p_1)\boldsymbol n$, with no component in the tangential plane.

We see that all terms in these two equations are proportional to $B_n$ or $u_n$, so that any non-ideal term, even small, can become dominant when these two quantities tend to zero (if these non-ideal terms do not tend to zero at the same time). As the distinction between compressional and rotational character fully relies on this system of equations, this evidences the necessity of investigating the quasi-tangential case for resolving the usual singularity of the tangential discontinuity. We note that the left- and right-hand sides of (4.5) can be put equal to zero by choosing the ‘De Hoffmann–Teller’ tangential reference frame where the electric field is zero (Belmont et al. Reference Belmont, Rezeau, Riconda and Zaslavsky2019). However, this choice, even if it can simplify some calculations, is not necessary here. Finally, the variables $\boldsymbol u_{t}$ can be eliminated from the system by a simple linear combination of the two equations, leading to

(4.6)\begin{equation} (u_{n2} -u_{n0})\boldsymbol{B}_{t2} = (u_{n1} -u_{n0})\boldsymbol{B}_{t1}, \end{equation}

where

(4.7)\begin{equation} u_{n0}=\frac{B_{n}^2} {\mu_0 \rho u_n} = \mathrm{cst}. \end{equation}

Equation (4.6) leads to the distinction between shocks, where the tangential magnetic fields on both sides are collinear (but with different modules), and rotational discontinuities, where the terms inside the brackets must be equal to zero. Rotational discontinuities correspond to a propagation velocity equal to the normal Alfvén velocity and imply $u_{n1}=u_{n2}=u_{n0}$, and therefore, an absence of compression of the plasma.

As previously stated, the separation between the compressional and rotational characters mainly derives, in CTD, from the assumption of isotropic pressures on both sides, which prevents the pressure divergences having tangential components. When the isotropic hypothesis is relaxed (Hudson Reference Hudson1971), the set of conservation equations is no longer sufficient to determine a unique downstream state for a given upstream one. As a consequence, the global result depends on the non-ideal processes occurring within the layer. In addition to anisotropy effects, finite Larmor radius (FLR) effects can be expected to break the gyrotropy of the pressure tensor around $\boldsymbol B$ in the case of thin boundaries between different plasmas. This means that the main effect that explains departures from CTD comes from the tangential component of the divergence of the pressure tensor, which must be taken into account in the momentum equation. On the other hand, the non-ideal effects related to the generalized Ohm's law are negligible, at least in the examples shown in this paper. The possible types of discontinuities in an anisotropic plasma have been discussed in several papers a long time ago (Abraham-Shrauner Reference Abraham-Shrauner1967; Lynn Reference Lynn1967; Chao Reference Chao1970; Neubauer Reference Neubauer1970), and the present paper improves the analysis in light of the new experimental possibilities given by the MMS measurements.

When the dynamics drives the conditions for the pressure tensor to become anisotropic (and a fortiori in the non-gyrotropic case) the $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {P}$ term comes into play linking upstream and downstream quantities. Considering the ‘simple’ anisotropic case, i.e. keeping the gyrotropy around $\boldsymbol B$, it has been shown (Hudson Reference Hudson1971) that the $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}$ term then just introduces a new coefficient:

(4.8)\begin{equation} \alpha = 1-\frac{p_{{\parallel}}-p_{{\perp}}}{ B^2 / \mu_0}. \end{equation}

This coefficient has been interpreted as a change in the Alfvén velocity $V'^2_{An}=\alpha V_{An}^2$, but it appears more basically as a change in (4.6):

(4.9)\begin{equation} (u_{n2} -\alpha_2 u_{n0})\boldsymbol{B}_{t2} = (u_{n1} -\alpha_1 u_{n0})\boldsymbol{B}_{t1}. \end{equation}

This equation shows that, in this simple anisotropic case, coplanar solutions still exist ($\boldsymbol {B}_{t2}$ and $\boldsymbol {B}_{t1}$ are collinear), but that whenever $\alpha _2$ is not equal to $\alpha _1$, the equivalent of the rotational discontinuity now implies compression, i.e.

(4.10)\begin{equation} u_{n2} \ne u_{n1} \quad \text{if} \ \alpha_2 \ne \alpha_1, \end{equation}

since $u_{n2}=\alpha _2 u_{n0}$ and $u_{n1}=\alpha _1 u_{n0}$. The variation of $u_n$ explains why the modified rotational discontinuity can be ‘evolutionary’ (Jeffrey & Taniuti Reference Jeffrey and Taniuti1964), the nonlinear steepening being counter-balanced at equilibrium by non-ideal effects for a thickness comparable with the characteristic scale of these effects.

There is actually no additional conservation equation available that would allow the jump of the anisotropy coefficient $\alpha$ to be determined. Consequently, there is no universal result that gives the downstream state as a function of the upstream one, regardless of the microscopic processes going on within the layer. This remains valid for the full anisotropic case, with non-gyrotropy. As soon as the ion Larmor radius $\rho _i$ and the ion inertial length $d_i$ are not fully negligible with respect to the characteristic scale $L$ of the layer, kinetic effects, and in particular FLR effects, which make the pressure tensor non-gyrotropic, must be taken into account to describe self-consistently the internal processes. Then, the effect of the divergence of the pressure tensor is no longer reduced to adding a coefficient $\alpha$ since its tangential component is no longer collinear with $\boldsymbol B_t$. Such an effect has already been reported and analysed in the context of magnetic reconnection (Aunai et al. Reference Aunai, Retinò, Belmont, Smets, Lavraud and Vaivads2011; Aunai, Hesse & Kuznetsova Reference Aunai, Hesse and Kuznetsova2013) and in kinetic modelling of purely tangential layers (Belmont, Aunai & Smets Reference Belmont, Aunai and Smets2012; Dorville et al. Reference Dorville, Belmont, Aunai, Dargent and Laurence2015a). It has also been investigated in the case of linear modes where they are responsible for the transition from shear Alfvén into kinetic Alfvén waves (Hasegawa & Uberoi Reference Hasegawa and Uberoi1982; Belmont & Rezeau Reference Belmont and Rezeau1987; Cramer Reference Cramer2001). On the other hand, it has never been introduced in the context of quasi-tangential discontinuities.

If a simple anisotropy preserving gyrotropy around $\boldsymbol {B}$ can be straight fully taken into account for modelling the pressure tensor and using it in fluid equations, introducing non-gyrotropy does not lead to a general and simple modelling for the pressure tensor. It would demand a priori a full kinetic description or, at least, some expansions assuming that these effects are small enough (see Braginskii (Reference Braginskii1965) for the pioneer work in this field, Passot & Sulem (Reference Passot and Sulem2006) and references therein). Several papers have investigated the changes in rotational discontinuities when such non-ideal effects are introduced (Lyu & Kan Reference Lyu and Kan1989; Hau & Sonnerup Reference Hau and Sonnerup1991; Hau & Wang Reference Hau and Wang2016). These theoretical papers used different analytical models based on different simplifying assumptions. Contrary to these papers, we will not use such kind of assumptions. Instead, we will just analyse the observed magnetic hodograms, and show that their shape is incompatible with a gyrotropic pressure.

5. The magnetopause normal

When studying the magnetopause with in situ measurements, the most basic geometric characteristic to be determined is the normal to its surface (which may vary during the crossing). An accurate determination of the magnetopause normal is actually a fundamental condition for determining reliable estimates of the normal components of both the magnetic and the mass fluxes. Moreover, having a good estimation of the normal direction is also necessary to determine the speed of the structure and its thickness. Quantitatively speaking, to determine the normal component of the magnetic field sufficiently well (assuming that $B_n/|\boldsymbol {B}|\sim 2\,\%$), an accuracy of the normal should be of the order of $\delta \theta < 1^\circ$. In the literature, a good accuracy of determination of the normal is considered to be of the order of 5 $\%$ (Denton et al. Reference Denton, Sonnerup, Russell, Hasegawa, Phan, Strangeway, Giles, Ergun, Lindqvist, Torbert, Burch and Vines2018).

Beyond determining the normal direction, some ‘reconstruction methods’ can be used to provide a more global view of the large-scale structure around the spacecraft. Although these methods have proven to provide remarkable results (Hasegawa et al. Reference Hasegawa, Sonnerup, Klecker, Paschmann, Dunlop and Rème2005; De Keyser Reference De Keyser2008; Denton et al. Reference Denton, Torbert, Hasegawa, Dors, Genestreti, Argall, Gershman, Le Contel, Burch, Russell, Strangeway, Giles and Fischer2020) they will not be used here (the first two studies assume the Grad–Shafranov equations to be valid, implying stationary MHD, and are therefore not appropriate to investigate the non-MHD effects such as the FLR effects).

Over the years, several methods have been developed with the purpose to precisely determine the normal direction (see, e.g. Haaland et al. Reference Haaland, Sonnerup, Dunlop, Balogh, Georgescu, Hasegawa, Klecker, Paschmann, Puhl-Quinn, Rème, Vaith and Vaivads2004; Shi et al. Reference Shi, Tian, Bai, Hasegawa, Degeling, Pu, Dunlop, Guo, Yao, Zong, Wei, Zhou, Fu and Liu2019). The most common is the minimum variance (MVA) introduced with the first measurements of the magnetic field in space (Sonnerup & Cahill Reference Sonnerup and Cahill1967; Sonnerup & Scheible Reference Sonnerup and Scheible1998). This method, which requires single spacecraft measurements, provides a global normal, i.e. a single normal vector for each entire time series across the boundary. The tool is based on the assumption that the boundary is a perfectly 1-D and stationary layer crossing the spacecraft. Other notable examples are the generic residue analysis technique (Sonnerup et al. Reference Sonnerup, Haaland, Paschmann, Dunlop, Rème and Balogh2006), which consists of a generalization of the MVA to other parameters than $\boldsymbol {B}$, and the BV method (Dorville et al. Reference Dorville, Belmont, Rezeau, Grappin and Retinò2014), which combines magnetic field and velocity data. Even though these methods can give an accurate normal determination (Dorville et al. Reference Dorville, Haaland, Anekallu, Belmont and Rezeau2015b), they provide, like MVA, a global normal and, thus, they cannot provide the necessary basis for investigating the variations of the magnetopause normal within the structure and test the possible departures from mono-dimensionality. Let us finally recall that waves and turbulence, which are always superimposed to the laminar magnetopause profiles, bring strong limitations in the normal direction accuracy for all methods, in particular these global ones.

In this context, multi-spacecraft missions have represented a fundamental step in increasing the accuracy of the normal determination, allowing one to determinate the gradients of the measured fields. A notable example is the minimum directional derivative (MDD, Shi et al. Reference Shi, Shen, Pu, Dunlop, Zong, Zhang, Xiao, Liu and Balogh2005) method. This tool generally uses the magnetic field data, but it must be kept in mind that it is not based on specific properties of this field. The MDD technique is a so-called ‘gradient based method’ since the calculation of the normal is based on the experimental estimation of the dyadic tensor $\boldsymbol {G} = \boldsymbol {\nabla } \boldsymbol {B}$. This tensor gradient can be obtained from multi-spacecraft measurements using the reciprocal vector method (Chanteur Reference Chanteur1998). The MDD method consists in diagonalizing the matrix $\boldsymbol {L}=\boldsymbol {G}\boldsymbol {\cdot } \boldsymbol {G}^{\rm T}$, finding the normal direction as the eigenvector corresponding to the maximum eigenvalue. Moreover, the gradient matrix can also be used for estimating the dimensionality of the boundary from the ratio between the eigenvalues. A way of finding a quantitative determination of this dimensionality was proposed in Rezeau et al. (Reference Rezeau, Belmont, Manuzzo, Aunai and Dargent2018).

For the vector $\boldsymbol {B}$, the MDD method makes use only of the spatial derivatives $\partial _i \boldsymbol {B}$, which are accessible at each time step thanks to the four-point measurements today available with multi-spacecraft space missions. In this sense, it is the opposite of the MVA method, which makes use only of the temporal variances of the $\boldsymbol {B}$ components. It therefore allows for an instantaneous determination of the normal at any point inside the layer, while MVA can only provide a single normal for a full crossing. In addition, contrary to MVA, MDD does not make any assumption about the geometry of the layer (1-D variations or not), and about the physical properties of the vector used. Indeed, it can be applied to the magnetic field data but also to any other vector since the property $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}=0$ is not used.

However, due to waves and turbulence, the magnetopause can present locally 2-D properties that are insignificant for the profiles we are looking for. For this reason, we will focus here on intervals where the magnetopause is mainly one dimensional, discarding the crossings in which local 2-D features are observed. The intervals considered as one dimensional are those for which $\lambda _{{\rm max}}\gg \lambda _{{\rm int}}$. Here $\lambda _{{\rm max}}$ and $\lambda _{{\rm int}}$ are defined as the highest and the intermediate eigenvalues of the matrix ${\boldsymbol{G}}$. In this limit, the ordering between $\lambda _{{\rm int}}$ and $\lambda _{{\rm min}}$ (i.e. the smaller eigenvalue) is not relevant in defining the intervals. Specifically, we use the parameter, $D_1 = (\lambda _{{\rm max}}-\lambda _{{\rm int}})/\lambda _{{\rm max}}$, which enables us to quantify this mono-dimensionality of the magnetopause as a function of time.

A more recent tool proposed to study the magnetopause is the hybrid method presented in Denton et al. (Reference Denton, Sonnerup, Hasegawa, Phan, Russell, Strangeway, Giles, Gershman and Torbert2016, Reference Denton, Sonnerup, Russell, Hasegawa, Phan, Strangeway, Giles, Ergun, Lindqvist, Torbert, Burch and Vines2018), in which the orientation of the magnetopause is obtained through a combination of the MDD and MVA methods, resulting in an improved accuracy of the normal direction.

The only limitation to the MDD accuracy comes from the uncertainty of the spatial derivatives that it uses. In particular, the local gradient matrix is calculated through the reciprocal vector technique (Chanteur Reference Chanteur1998), which assumes linear variations between the spacecraft. Because small-scale waves and turbulence are always superimposed on the magnetopause profiles being searched for, this assumption cannot be well respected without some filtering. This filtering actually leads to introducing part of the temporal information on the variations, but it still allows keeping local information inside the layer whenever one filters only the scales sufficiently smaller than those associated to the full crossing duration. The quality of the filtering is therefore the biggest challenge to complete for getting accurate results. For instance, simple Gaussian filters done independently on the four spacecraft would provide insufficient accuracy: this can be observed by the fact that, when doing it, the relation $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {B}=0$ is violated in the result. In the following section, it is shown how the MDD method can be included in a fitting procedure of the four spacecraft simultaneously and where this relation can be imposed as a constraint. We also show that, when no constraint is added, this procedure justifies the use of MDD with data that are filtered independently.

6. A new tool

The tool we present here, namely GF2 (gradient matrix fitting), has been derived from the MDD method. The digit 2 indicates that in the version of the tool that we use here the data are fitted with a 2-D model (it can be shown that fitting with a 1-D model is mathematically equivalent to the standard MDD technique used with smoothed data). Differently from the original method, we assume that the structure under investigation can be fitted locally (i.e. in each of the small sliding window used along the global crossing), by a 2-D model. This does not imply that the magnetopause is assumed globally two dimensional. As for MDD, the instantaneous gradient matrix ${\boldsymbol{G}}$ is obtained from the data using the reciprocal vector's technique (Chanteur Reference Chanteur1998). When performing the 2-D fit in each sliding window, we then impose some physical constraints, which could be checked only a posteriori with the classic MDD method.

The model $\boldsymbol {G}_{{\rm fit}}$ is obtained as follows:

(6.1)\begin{equation} \boldsymbol{G}_{{\rm fit}}=\boldsymbol{e}_0 \boldsymbol{B}_{e0}'+\boldsymbol{e}_1 \boldsymbol{B}_{e1}'. \end{equation}

Here we define $\boldsymbol {e}_0$ and $\boldsymbol {e}_1$ as two unit vectors in the plane perpendicular to the direction of invariance and $\boldsymbol {B}_{e0}'$ and $\boldsymbol {B}_{e1}'$ as the variation of the magnetic field along these two directions.

By performing the fit, we impose $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {B}=0$ (as used in MVA but ignored in standard MDD). In the model, this can be written as

(6.2)\begin{equation} \boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{B}_{e0}'+\boldsymbol{e}_1 \boldsymbol{\cdot}\boldsymbol{B}_{e1}'=0. \end{equation}

In order to fit the experimental $\boldsymbol {G}$ by the model $\boldsymbol {G}_{{\rm fit}}$, the following quantity has to be minimised:

(6.3)\begin{align} D_2& =\operatorname{{\rm Tr}} [(\boldsymbol{G}_{{\rm fit}}-\boldsymbol{G})\boldsymbol{\cdot}(\boldsymbol{G}_{{\rm fit}}-\boldsymbol{G})^{\rm T}]\nonumber\\ & =\boldsymbol{B}_{e0}'^2-2\boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{G}\boldsymbol{\cdot} \boldsymbol{B}_{e0}'+\boldsymbol{B}_{e1}'^2-2\boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{G}\boldsymbol{\cdot} \boldsymbol{B}_{e1}'+\operatorname{{\rm Tr}}(\boldsymbol{G} \boldsymbol{G}^{\rm T}). \end{align}

We can disregard the last term, since it is independent of the fit parameters. To impose the physical constraints, we use Lagrange multipliers, minimizing

(6.4)\begin{align} D_2& =\boldsymbol{B}_{e0}'^2-2\boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{G}\boldsymbol{\cdot} \boldsymbol{B}_{e0}'+\boldsymbol{B}_{e1}'^2-2\boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{G}\boldsymbol{\cdot} \boldsymbol{B}_{e1}'+2\lambda(\boldsymbol{e}_0{\cdot}\boldsymbol{B}_{e0}' +\boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{B}_{e1}')\nonumber\\ & =\boldsymbol{B}_{e0}'^2-2\boldsymbol{e}_0\boldsymbol{\cdot}(\boldsymbol{G}-\lambda\boldsymbol{I}) \boldsymbol{\cdot}\boldsymbol{B}_{e0}'+\boldsymbol{B}_{e1}'^2-2\boldsymbol{e}_1\boldsymbol{\cdot} (\boldsymbol{G}-\lambda\boldsymbol{I})\boldsymbol{\cdot}\boldsymbol{B}_{e1}'. \end{align}

By assuming in the first approximation that the direction of invariance $\boldsymbol {e}_2$ is known, we can choose the two vectors $\boldsymbol {e}_0$ and $\boldsymbol {e}_1$ as an arbitrary orthonormal basis for the plane of variance. For performing the minimisation, we have just to impose equal to zero the derivatives with respect to $\boldsymbol {B}_{e0}'$, $\boldsymbol {B}_{e1}'$ and $\lambda$, obtaining (6.2) and

(6.5)\begin{gather} \boldsymbol{B}_{e0}'=\boldsymbol{e}_0\boldsymbol{\cdot}(\boldsymbol{G}-\lambda\boldsymbol{I}), \end{gather}

(6.6)\begin{gather}\boldsymbol{B}_{e1}'=\boldsymbol{e}_1\boldsymbol{\cdot}(\boldsymbol{G}-\lambda\boldsymbol{I}). \end{gather}

By introducing these two equations in (6.2) we obtain

(6.7)\begin{equation} \lambda=\frac{G_{00}+G_{11}}{2}, \end{equation}

from which we get the values of $\boldsymbol {B}_{e0}'$ and $\boldsymbol {B}_{e1}'$. At this point, the matrix $\boldsymbol {G}_{{\rm fit}}$ is fully determined. We can then look for its eigenvalues and eigenvectors, as in the standard MDD method, and get the normal $\boldsymbol {n}$ and the tangential directions $\boldsymbol {t}_1$ (i.e. the one orthogonal to the direction of invariance) from this smooth fit.

The choice of the direction of invariance has actually no major influence on the determination of the normal direction, neither on the estimation of the 2-D effects. For large 2-D effects, one could choose the direction of MVA obtained by applying directly the standard MDD method to the data. Nevertheless, for almost 1-D cases (the most common situation), the spatial derivatives in the tangential directions are generally much smaller than the noise, so this result is not reliable. We simply choose here the constant M direction given by MVA, which is often considered as the direction of the X line if interpreted in the context of 2-D models of magnetic reconnection (cf. for instance, Phan et al. Reference Phan, Shay, Gosling, Fujimoto, Drake, Paschmann, Oieroset, Eastwood and Angelopoulos2013 for a typical use of this choice and Aunai et al. Reference Aunai, Hesse, Lavraud, Dargent and Smets2016, Denton et al. Reference Denton, Sonnerup, Russell, Hasegawa, Phan, Strangeway, Giles, Ergun, Lindqvist, Torbert, Burch and Vines2018, Liu et al. Reference Liu, Hesse, Cassak, Shay, Wang and Chen2018 for discussions about it).

Finally, another useful by-product of the method can be obtained: comparing the spatial derivatives and the temporal ones and using a new fitting procedure, we can compute the two components of the velocity of the structure $V_{n0}$ and $V_{t1}$ with respect to the spacecraft. Only the motion along the invariant direction then remains unknown.

6.1. Normal from ions mass flux

This tool can be easily adapted to any other vector dataset by just changing the physical constraint. In particular, we chose to study the structure using the ion mass flux data. In this case we impose mass conservation $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {\varGamma }_i=-\partial _t n_i$ (with $\boldsymbol {\varGamma }_i=n_i\boldsymbol {u}_i$). Equation (6.2) now writes

(6.8)\begin{equation} \boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{\varGamma}_{e0}'+\boldsymbol{e}_1\boldsymbol{\cdot} \boldsymbol{\varGamma}_{e1}'+\partial_t. n_i=0. \end{equation}

Therefore, when using the Lagrange multipliers, (6.4) changes to

(6.9)\begin{align} D_2& =\boldsymbol{\varGamma}_{e0}'^2-2\boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{G}\boldsymbol{\cdot} \boldsymbol{\varGamma}_{e0}'+\boldsymbol{\varGamma}_{e1}'^2-2\boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{G} \boldsymbol{\cdot}\boldsymbol{\varGamma}_{e1}'+2\lambda(\boldsymbol{e}_0\boldsymbol{\cdot} \boldsymbol{\varGamma}_{e0}'+\boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{\varGamma}_{e1}'+\partial_t n_i)\nonumber\\ & =\boldsymbol{\varGamma}_{e0}'^2-2\boldsymbol{e}_0\boldsymbol{\cdot}(\boldsymbol{G}-\lambda\boldsymbol{I}) \boldsymbol{\cdot}\boldsymbol{\varGamma}_{e0}'+\boldsymbol{\varGamma}_{e1}'^2-2\boldsymbol{e}_1\boldsymbol{\cdot} (\boldsymbol{G}-\lambda\boldsymbol{I})\boldsymbol{\cdot}\boldsymbol{\varGamma}_{e1}'+2\lambda\partial_t n_i. \end{align}

By using the same algorithm as above, the constraint can now be written as

(6.10)\begin{equation} \lambda=\frac{G_{00}+G_{11}+\partial_t n_i}{2}. \end{equation}

6.2. Dimensionality index

From this procedure, we can also derive another significant result: we can obtain an indicator of the importance of the 2-D effects in the profiles, free of the parasitic noise effects. Specifically, we can estimate the variation of the magnetic field along the normal by projecting the $\boldsymbol {G}_{{\rm fit}}$ matrix along it $\text {var}_n=|\boldsymbol {G}_{{\rm fit}}\boldsymbol {\cdot } \boldsymbol {n}|$. Consequently, if we designate the variation along $\boldsymbol t_1$ as $\text {var}_t$, we can introduce a new dimensionality index:

(6.11)\begin{equation} \mathscr{D}_{{\rm GF2}}=\frac{\text{var}_n-\text{var}_t}{\text{var}_n}. \end{equation}

This index can usefully be compared with the instantaneous index $D_1 = (\lambda _{{\rm max}}-\lambda _{{\rm int}})/\lambda _{{\rm max}}$ of Rezeau et al. (Reference Rezeau, Belmont, Manuzzo, Aunai and Dargent2018).

7. Expected accuracy and tests of the tool

In this section we test the accuracy of the GF2 tool. To accomplish this, we exploit a case crossing, which will be investigated in detail in the following section. The crossing considered comes from MMS data (Burch & Phan Reference Burch and Phan2016), taking place at around 22:11 on 28th December 2015. For this study, we use data from the FluxGate magnetometers (FGM, Russell et al. Reference Russell2016), providing the magnetic field data, the electric double probe (Ergun et al. Reference Ergun, Tucker, Westfall, Goodrich, Malaspina, Summers, Wallace, Karlsson, Mack, Brennan, Pyke, Withnell, Torbert, Macri, Rau, Dors, Needell, Lindqvist, Olsson and Cully2016; Lindqvist et al. Reference Lindqvist2016), for the electric field, and dual ions and electrons spectrometer instrument (Pollock et al. Reference Pollock2016), for plasma measurements. An overview of the event is shown in figure 2, where both the magnetic field and ion bulk velocity are given in geocentric solar ecliptic (GSE) coordinates. For this event, the spacecraft are located in $[7.6, -6.7, -0.8]$ $R_E$ in GSE coordinates (where $R_E$ is the Earth's radius).

Figure 2. Main features of the crossing of the 28th December 2015. From top to bottom: (a) the magnetic field (in ${\rm nT}$), (b) the ion particle density (in m$^{-3}$), (c) ion velocity (in km s$^{-1}$), (d) total current (computed from the curlometer (Dunlop et al. Reference Dunlop, Southwood, Glassmeier and Neubauer1988), in nA m$^{-2}$), (e) the ion and ( f) electron spectrograms (energies are shown in eV). Vertical lines indicate the time interval chosen for the case study.

The temporal interval in which we observe the shear in the magnetic field and the crossing in the particle structure is about 8 s, enough to allow for high resolution for both sets of measurements. The crossing is chosen by also analysing the dimensionality of the magnetic field measurements averaged along the crossing. In particular, the dimensionality parameter defined in (6.11), denoted as $\mathscr {D}_{{\rm GF2}}$, and the one introduced in Rezeau et al. (Reference Rezeau, Belmont, Manuzzo, Aunai and Dargent2018), denoted as $D_1$, were considered. In this interval, indeed, we have $D_{1, {\rm mean}}=0.97 \pm 0.03$ while $\mathscr {D}_{{\rm GF2}}=0.89\pm 0.06$, both highlighting that the crossing exhibits 1-D features throughout the time interval. We remind here that in burst mode, the frequency of magnetic field measurements is 132 Hz while it is 6.67 Hz for ions. To conduct the following study, it is necessary to interpolate all measurements at the same times. We did it by testing two sampling frequencies: the magnetic field and the ion ones. The results obtained are consistent with the two methods. All figures shown in the paper are obtained with the sampling times of the MMS1 magnetic field. Furthermore, the crossing is observed quasi-simultaneously for the two quantities, with a large interval where the two kinds of results can be compared.

As a first test, we compare in figure 3 the normals obtained by GF2 and those by the standard MDD technique (using data smoothed in a 0.31 s time window), for both the magnetic field and the ion data. For reference, we also compare the result of ${\rm GF2}_B$ with the MVA one.

Figure 3. Comparison for the normals obtained with GF2 with respect to the MDD tools. Panel (a) shows the magnetic field and (b) the ion mass flux, measured by the four MMS spacecraft. Panels (c,d) show the magnetic and the ion normal, respectively. The continuous (respectively dashed) line correspond to the components of GF2 (respectively MDD) normal. Horizontal dotted lines indicate the MVA normal obtained along the whole interval. Vertical dashed lines correspond to the time interval boundaries for the crossing, which are different for the magnetic field and the ion mass flux.

Vertical dashed lines indicate the time interval during which all the satellites are inside the boundary. We observe that the time required for the ions flux to complete the crossing (of about 5 s) is shorter than for the magnetic field (about 8 s). To perform a quantitative analysis of the differences, we studied the angles between the different normals obtained through GF2, MDD and MVA, as shown in figure 4(a).

Figure 4. (a) Angle between the normals obtained using the state-of-the-art tools (MDD, MVA) and GF2. The subscripts $B$ and $ions$ indicate whenever the magnetic field or the ion flux measurements are used. (b) Dimensionality of the structure as a function of time; here both the $\mathscr {D}_{{\rm GF2}}$ (continuous line) and the $D_1$ (dashed line) indices are shown, for both the magnetic field (blue) and ions (red) data.

The first striking result is that all these results are quite consistent. Almost all the directions are less than 10$^\circ$ apart from each other, with an average difference of about 5$^\circ$. The major exception concerns the comparison between MVA and ${\rm GF2}_B$ during the last second of the interval where the two directions appear to be up to 35$^\circ$ apart. This can be explained by the fact that the local normals are observed (by ${\rm GF2}_B$ as well as by ${\rm MDD}_B$) to differ noticeably in this part from their averaged value and that MVA is not able to detect such a change. Looking in more detail, we can see a slight difference between the first part of the crossing (between 2 and $\sim$6.5 s), where the two normals ${\rm GF2}_B$ and ${\rm MDD}_B$ differ by less than 5$^\circ$, and the second part, where the angle between the normals can be up to 10$^\circ$ (probably due to a smaller ratio signal/noise for the gradient matrix $\boldsymbol {G}$). The normals derived from ion measurements are not much different from those derived from the magnetic field, showing that the particle and magnetic structures are approximately identical. In figure 4(b) the dimensionality of the structures is analysed as a function of time, by using both the $\mathscr {D}_{{\rm GF2}}$ and the $D_1$ (Rezeau et al. Reference Rezeau, Belmont, Manuzzo, Aunai and Dargent2018) parameters, as explained above. Even if the numerical values of the two indices are slightly different, they both indicate structures close to one-dimensionality in the first part, with a – small but significant – decrease in the second part. This increased departure from mono-dimensionality can explain the slight difference between the two parts when comparing the normals from standard MDD and GF2 techniques.

The present test does not allow us to state that GF2 is more accurate than standard MDD (this will be checked in future work by comparing the two tools in a global simulation involving realistic turbulence) but it shows at least a good agreement between the two approaches. We will see in the following that this accuracy is anyway sufficient to prove the role of FLRs.

In order to smooth the small fluctuations over the time interval and to reduce the statistical error associated with the determination of the normal, we can compare the directions averaged along the crossing time. Mean values obtained through the tools presented above are shown in table 1. Here we observe that all the averaged normals differ by less than 10$^\circ$. Specifically, we observe that the normals obtained with GF2, MDD and (Denton et al. Reference Denton, Sonnerup, Russell, Hasegawa, Phan, Strangeway, Giles, Ergun, Lindqvist, Torbert, Burch and Vines2018) are similar, with a difference of less than 1$^\circ$ (with ours being closer to that from Denton et al. Reference Denton, Sonnerup, Russell, Hasegawa, Phan, Strangeway, Giles, Ergun, Lindqvist, Torbert, Burch and Vines2018). The MVA normal, instead, differs around 6$^\circ$ from all these other normals. Finally, we also observe that the one computed with ions flux data is the most distant. This is interpreted to be due to the higher uncertainty of particles measurements.

Table 1. Magnetopause normal vectors obtained with the main tools presented above averaged in the time interval and their angle with respect to the normal obtained with GF2 using the magnetic field data (in degrees).

8. Case study

In this section we undertake a detailed analysis of the previously mentioned crossing case by employing the normal obtained using the GF2 tool. Here, we focus on the time interval between 2 and 9 s in figure 2. To mitigate the potential influence of non-unidimensionality effects, we chose to exclude the last second of the time interval studied in the preceding section for the magnetic field (where both $\mathscr {D}_{{\rm GF2}}$ and $D_1$ show that the structure is less one dimensional and where we observe that the normal is more different from the averaged one). To carry out this analysis, we study the hodogram of the magnetic field in the tangential plane. Here the tangential results are presented in a basis ($\boldsymbol {T}_1, \boldsymbol {T}_2$) chosen as

(8.1)\begin{gather} \boldsymbol{T}_1=\boldsymbol{n}_{{\rm mean}} \times \hat{\boldsymbol{b}}, \end{gather}

(8.2)\begin{gather}\boldsymbol{T}_2=\boldsymbol{n}_{{\rm mean}} \times \boldsymbol{T}_1. \end{gather}

where $\hat {\boldsymbol {b}}=\boldsymbol {B}/|\boldsymbol {B}|$ and $\boldsymbol {n}_{{\rm mean}}$ is the direction of the averaged normal in the chosen time interval. Note that the choice of the reference frame ($\boldsymbol {T}_1, \boldsymbol {T}_2$) is just a convention. The shape of the hodogram remains unaffected by this choice except for the corresponding rotation in this tangential plane. The direction $\boldsymbol {t}_1$, which characterizes the direction of the second dimension of the model in GF2 and that is also in the tangential plane, is generally different from $\boldsymbol {T}_1$.

If CTD was valid everywhere, the hodogram of the magnetic field in the tangential plane for a rotational discontinuity would correspond to a circular arc with constant radius while a shock would correspond to a radial line (as shown in figure 1). For this reason, the hodogram is a good tool to recognize the cases for which the CTD fails at describing the magnetopause. The hodogram for this case is shown in figure 5. We observe a clear ‘linear’ (although not radial) hodogram. This non-radial variation of the magnetic field although not predicted by CTD, is a striking feature of the hodogram. It cannot be explained by a departure from the 1-D assumption since we have measured that the crossing can be considered as one dimensional to a good degree of accuracy. It is therefore due to an intrinsic property of the layer itself. Also, in figure 5 we present the hodogram derived from the local normal (un-averaged, violet line). It is clear that averaging does not affect the shape of the hodogram.

Figure 5. Hodogram in the tangential plane of the magnetic field for a magnetopause crossing by MMS in 28.12.2015 from 22:12:02 to 22:12:09. See text for the significance of the arrows. Here $B_{T1}$ and $B_{T2}$ are the projections of $\boldsymbol {B}$ along the tangential directions computed as described in the text. The black line (respectively violet) is the hodogram when the $\boldsymbol {n}_{{\rm mean}}$ (respectively $\boldsymbol {n}$) value is used to define the reference frame.

To further analyse the causes of the disagreement between the hodogram of this case crossing and what is expected from CTD, we compare the different terms of the tangential momentum equation and Faraday/Ohm's law. As discussed above, indeed, these two equations are responsible for the distinction between the rotational and tangential discontinuities in CTD. This is the object of figure 6, where we plot the different terms of the two equations projected along the $\boldsymbol {n}_{{\rm mean}}$ (panels (a) and (b)) and $\boldsymbol {t}_{1,{\rm mean}}$  (panels (c) and (d)) directions obtained using the GF2 tool (averaged over the whole time interval). The influence of the averaging of the $\boldsymbol {t}_{1,{\rm mean}}$ direction on the results is discussed in Appendix A. We do not show the quantities along the direction of invariance, which are dominated by noise. The current and the gradient matrix for the pressure term are obtained via the reciprocal vector method described in Chanteur (Reference Chanteur1998).

Figure 6. Terms of Ohm's law (panels (a) and (c), units of mV m$^{-1}$) and the momentum equation (panels (b) and (d), units of $10^{-15}$ kg m s$^{-2}$), projected in the normal direction n (a,b) and in the tangential direction $t_1$ (c,d ). To reduce the noise, a running average with a time window of 0.35 s is applied to the electric field measurements. Shaded regions in panel $(d)$ represent the estimated uncertainties of the divergence of the pressure (red), the $\boldsymbol {J}\times \boldsymbol {B}$ (blue) and the classic inertial term (green). Concerning Ohm's law, we included the sum $\boldsymbol {U}\times \boldsymbol {B}-\boldsymbol {J}\times \boldsymbol {B}/nq$ to facilitate the readability (blue dashed line). Note that the terms of the tangential Faraday/Ohm's law used in the text are just the derivatives of those in (a) (apart from a $\pi /2$ rotation).

Concerning Ohm's law (figure 6, panels (a) and (c)), we see that the electric field is well counter-balanced by the $\boldsymbol {u}\times \boldsymbol {B}$ and $\boldsymbol {J}\times \boldsymbol {B}/nq$ terms (ideal and Hall terms). Outside the layer, on both sides, the ideal Ohm's law is satisfied, as assumed in CTD (this is not visible on the figure, which is a zoom on the inner part of the layer, and where the Hall term is important). It has been shown in the literature that $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_e/nq$ is not always negligible in Ohm's law and that it can even be dominant close to an electron diffusion regions. This has been predicted theoretically (Hesse et al. Reference Hesse, Neukirch, Schindler, Kuznetsova and Zenitani2011, Reference Hesse, Aunai, Sibeck and Birn2014) and observed experimentally (Torbert et al. Reference Torbert2016; Genestreti et al. Reference Genestreti, Nakamura, Nakamura, Denton, Torbert, Burch, Plaschke, Fuselier, Ergun, Giles and Russell2018), but it is not the case for events like this one. We observe that at approximately $3.5$ s, the $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_e/nq$ is not entirely negligible along the tangential direction (a similar peak can also be observed in panel $(c)$ for the term associated with the electron pressure in the momentum equation). However, during this time interval, this value is not dominant, this term being smaller than both the electric field and the $\boldsymbol {J}\times \boldsymbol {B}/nq$ components. Furthermore, this effect exhibits a local characteristic, as $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_e/nq$ is only non-negligible within a small subinterval (with respect to the magnetopause temporal width). It is therefore not likely to be indicative of proximity to a reconnection point.

Concerning the momentum equation, shown in panels (b) and (d) of figure 6, we observe that, in the normal direction (panel b), the $\boldsymbol {J}\times \boldsymbol {B}$ term is counter-balanced by the divergence of the ion pressure tensor, as expected. But, if the isotropic condition assumed in CTD was valid, we would expect the divergence of the ion pressure tensor to be zero in the tangential direction, or at least negligible with respect to the inertial term $\rho \,{\rm d}\boldsymbol {u}/{\rm d}t$. On the contrary, we observe that the $\boldsymbol {J}\times \boldsymbol {B}$ term along $\boldsymbol {t}_1$ is of the same order of magnitude as the divergence of the ion pressure tensor, and one order of magnitude larger than all the other terms. Panel $(d)$ also shows an estimation of the error on the relevant terms: $\boldsymbol {J}\times \boldsymbol {B}$, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_i$ and the classical inertial term. It is known that measurement errors are difficult to estimate, especially at small scales. In order to validate our results, however, we sought to obtain an upper bound of the error associated with the quantities of interest. For that purpose, an overestimation of the uncertainty of the measurements (acquired as the maximum during the crossing of the errors available in Fast Plasma Investigation (FPI) datasets for the pressure tensor and from the FGM nominal error for the magnetic field) was exploited. These values are propagated as a statistical (i.e. quadratic) error (by assuming that the errors on the reciprocal vectors can be neglected with respect to that of other physical quantities).

From panel $(d)$ of figure 6, we see that the $\boldsymbol {J}\times \boldsymbol {B}$ and the $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_i$ terms are pointing in opposite directions and balancing each other. If valid in the first part of the interval, this conclusion cannot be safely trusted due to measurement uncertainty, but we observe that in the middle part (particularly between 3.5 and 6 s) it is evident that the two quantities counterbalance each other while the classical inertia term $\rho \,{\rm d}\boldsymbol {u}/{\rm d}t$ is much lower with respect to the others. This proves that the tangent $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_i$ term plays a fundamental role in the magnetopause equilibrium.

This point can be emphasized also by analysing the hodogram. In figure 5 the arrows are directed along the directions of the tangential plane that are physically relevant for the problem: (i) the tangent to the hodogram (green), which indicates the total variation of $\boldsymbol B_t$; (ii) the radial direction (red), which corresponds to the plasma compression; (iii) the $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_{it}$ direction (blue), which is the direction of the divergence of the ion pressure tensor in the tangential plane, and therefore, corresponds to a term that is absent in CTD. The relative lengths of the arrows are chosen proportional to the corresponding term magnitudes. These directions are averaged in two subintervals (bold hodogram). The striking result is that the total variation is mainly determined by the non-classic term $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {P}_{it}$ and not by the radial classic one. This explains the very recurrent (even if not reported in the literature hitherto) feature that the hodograms are almost linear but not radial.

8.1. Comparison of the width of the magnetopause to relevant physical lengths

Finally, we compare the width of the magnetopause ($L$) to the two main ion-related lengths: the ion Larmor radius ($\rho _L$) and the ion inertial length ($d_i$). The magnetopause width is estimated using the normal velocity obtained from the GF2 tool. By averaging the velocity of the magnetopause in the normal direction, we can estimate $L=\boldsymbol {V}_{n, {\rm mean}}\Delta t$ (where $\Delta t$ is the time length of the full crossing). These three scales are shown in figure 7. We observe that this width is larger than the ion Larmor radius and the ion inertial length all across the crossing, but only two to five times larger, which appears sufficient to drive the observed kinetic effects.

Figure 7. Comparison of the magnetopause width ($L$) with the ion inertial length ($d_i$) and the ions Larmor radius ($\rho _L$). Vertical lines highlight the considered temporal interval.

9. Ion pressure tensor analysis

To further investigate the question of the ion non-gyrotropy with respect to the magnetic field and quantify this effect, let us now examine the properties of the ion pressure tensor and introduce a new non-gyrotropy index. For that purpose, we define the matrices $\boldsymbol {P}_\parallel =p_\parallel \boldsymbol {bb}$, where $\boldsymbol {b}=\boldsymbol {B}/|\boldsymbol {B}|$, $p_\parallel =\boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {P}_i \boldsymbol {\cdot }\boldsymbol {b}$ and $\boldsymbol {P}_\perp =\boldsymbol {P}_i-\boldsymbol {P}_\parallel$. By defining $p_1$ and $p_2$ the maximum and intermediate eigenvalues of the $\boldsymbol {P}_\perp$ matrix, we define

(9.1)\begin{equation} D_{ng, \perp}=\frac{p_1-p_2}{p_1 + p_2}. \end{equation}

In figure 8(a) this parameter is compared with the non-gyrotropy index presented in Aunai et al. (Reference Aunai, Hesse and Kuznetsova2013). The two indices define non-gyrotropy differently, Aunai et al. (Reference Aunai, Hesse and Kuznetsova2013)'s index defining non-gyrotropy as the ratio of the non-gyrotropic to the gyrotropic part of the tensor (instantaneous), while ours makes use of the 2-D modelling of the data used in GF2 (averages on sliding windows). We note how both indices are significantly different from zero, approximately of the order of 0.1 within the boundary, corresponding to clearly present, although not predominant, non-gyrotropic effects. We note a decrease in both indices outside the magnetopause, as expected, but it is worth noting also that, despite a continuous decrease, these indices remain relatively high in the time interval just preceding the crossing, in a region where the magnetic field, density and pressure tensor are almost constant. This can be understood by noting that an ion velocity gradient is observed in this interval, suggesting that the non-diagonal terms of the pressure tensor could be due there to a kind of gyroviscous effect, the non-diagonal terms of the pressure tensor (Braginskii Reference Braginskii1965) being due to FLRs (Stasiewicz Reference Stasiewicz1993). One must take care that, in this interval, the pressure tensor has low values characterized by larger relative errors, which could partially influence this result. To further analyse this question, we have estimated the uncertainties on both non-gyrotropy indices. This estimation is derived from the nominal uncertainties of the FPI dataset. The diagonal terms have higher values and lower relative errors. Concerning the time interval before the crossing that we discuss here, the diagonal terms have errors of approximately 5 %, whereas off-diagonal terms have an average relative error about 50 %. We observe on figure 6 that this way of estimating the uncertainty well encompasses the variance of the results. It confirms that, within the crossing interval, all relative errors are smaller than 10 %, as considered in the Ohm's law study (figure 6).

Figure 8. Panels ($a$) and ($b$) show the evolution of the $D_{ng, \perp }$ and $D_{ng, {\rm Aunai}}$ (Aunai et al. Reference Aunai, Hesse and Kuznetsova2013) indices, respectively, along with their estimated uncertainties. Thin lines correspond to the real-time values while thick lines to an averaged window of 1 s. (c) Evolution of the eigenvalues of the $\boldsymbol {P}_i$ matrix (averaged on the four spacecraft). The dotted line indicates the magnetopause crossing. The red dotted lines in panel ($c$) highlight the time interval studied in figure 9.

In addition, a preliminary study appears to confirm the validity of the gyroviscous interpretation. Using the theoretical expressions given in Stasiewicz (Reference Stasiewicz1989), we can compare the variations of the non-diagonal terms of the pressure tensor with the spatial derivatives of the flow velocity, and evidence a fairly good correlation (see Appendix B).

Figure 8(c) also shows the evolution of the eigenvalues of the $\boldsymbol {P}_i$ tensor, averaged on the four spacecraft. This figure shows how outside of the magnetopause the three eigenvalues tend to converge towards each other meaning that these media are close to isotropy. However, inside the magnetopause, we note a transition in the behaviour of the intermediate eigenvalue, shifting from a value close to the minor one to being closer to the major eigenvalue. The minor eigenvalue exhibits a significant deviation from the other two towards the last part of the crossing.

Focusing on the temporal interval marked by the red square in figure 8, this transition is further investigated in figure 9 where we show the ions’ distribution functions in the tangential plane (with respect to the magnetopause) for four different intervals during the crossing, highlighting the non-gyrotropy of the ions’ distribution function over time. Velocity distribution functions (VDFs) (printed using a linear 2-D interpolation on a Cartesian grid in the chosen plane using the Pyspedas library) are here averaged in the corresponding time intervals framed with the same colour as in the bottom plot where the eigenvalues of the ion pressure tensor are plotted again (the time length decreases as the density increases).

Figure 9. (a) Ions’ VDFs in the tangential plane (the $\boldsymbol {T}_1$–$\boldsymbol {T}_2$ plane) averaged in four different time periods. Velocity axes are between $-$220 and 220 km s$^{-1}$. (b) Eigenvalues of the pressure tensor (same interval as in the red dashed square of figure 8c). The four coloured boxes are used to distinguish the four time intervals.

Finally, we analysed the non-gyrotropy with respect to a generic direction, i.e. without imposing that this direction is the magnetic field direction. Specifically, we have looked at a direction, denoted as $\boldsymbol {g}$, around which the rotated matrix could be rewritten as follows:

(9.2)\begin{equation} \begin{pmatrix} P_2 & 0 & 0\\ 0 & P_1 & 0\\ 0 & 0 & P_1 \end{pmatrix}. \end{equation}

To achieve this, we employ a minimization algorithm to derive the rotation matrix $\boldsymbol {M}$ that allows us to put the pressure tensor data under a form as close as possible to this one. Results from this study are shown in figure 10 (here shown for MMS2). Figure 10(a) displays the variation of $P_1$ and $P_2$ along the crossing, where $P_2$ consistently exceeds $P_1$. In addition, we imposed an upper limit on the temporal variation of the gyrotropic direction $\boldsymbol {g}$, excluding points with significant temporal variations (indicated by the thin line). Consequently, the remaining points reflect instances where the direction of $\boldsymbol {g}$ can be considered as stable and reliable. The vector $\boldsymbol {g}$ itself is represented in panels ($b$) and ($c$), where it is clear that the direction of gyrotropy is not close to the magnetic field direction: it is close to $\boldsymbol {n}_{{\rm mean}}\times \boldsymbol {B}$, the component along $\boldsymbol {B}$ being smaller and varying. This result reminds us that at boundaries such as the magnetopause, the strong gradients can break the isotropy as much, and even more here, than the magnetic field, so that gyrotropy can be around another vector than $\boldsymbol {B}$. A similar remark had already been made in Belmont et al. (Reference Belmont, Aunai and Smets2012) concerning the modelling of a tangential discontinuity.

Figure 10. (a) Evolution of parameters $P_1$ and $P_2$. (b,c) Projections of the gyrotropy direction in two planes. The ordinate is the direction of $\boldsymbol {B}$, the abscissa is the direction of $\boldsymbol {n}_{{\rm mean}}\times \boldsymbol {B}$ for panel (b) and $\boldsymbol {n}_{{\rm mean}}$ for panel (c).

10. Dataset selection

In order to expand the results on a statistical basis, we selected a dataset of 146 crossings, chosen from the largest one reported in Nguyen et al. (Reference Nguyen, Aunai, Michotte De Welle, Jeandet, Lavraud and Fontaine2022) and Michotte De Welle et al. (Reference Michotte De Welle, Aunai, Nguyen, Lavraud, Génot, Jeandet and Smets2022). From this database, the following conditions were required in order to carry out an accurate study.

1. (i) The MMS data are in burst mode.

2. (ii) The crossing duration is between 3 and 15 s. Too short crossings do not have a sufficient number of points within the structure (ion measurements are every 0.15 s). Too long crossings may imply non-stationary structures.

3. (iii) Partial crossings are discarded. For that, we impose a density threshold less than 4 cm$^{-3}$ in the magnetosphere and larger than 15 cm$^{-3}$ in the magnetosheath.

4. (iv) Only cases presenting simultaneous crossing features in the particles and magnetic field are considered, in order to compare normals computed at the same time.

In addition to these basic requirements, we also excluded some of the selected crossings for criteria that demand a more detailed analysis of the internal structure of the boundary. First, we excluded 2-D features. The quantitative determination of the dimensionality was done with the parameters presented in Rezeau et al. (Reference Rezeau, Belmont, Manuzzo, Aunai and Dargent2018) and the dimensionality index presented in § 6.2, which are functions of the ratio between the eigenvalues of the gradient matrix. Namely, we considered only crossings with $D_1>0.9$ and $\mathscr {D}_{{\rm GF2}}>0.8$, these two parameters being averaged on the crossing interval. These parameters are calculated at each time step but, due to waves and turbulence, attention must be paid to the fact that some of these 2-D features can be only local and insignificant for the profiles we are looking for. It is the reason why we use only the averaged values. The 146 selected crossings span from September 2015 to December 2017 (included). We can observe in figure 11 that the crossings are evenly distributed in the $x$, $y$ plane. Regarding the $z$ component, there is a prevalence of cases at negative $z$.

Figure 11. Spatial distribution of the selected database of crossings on the (a) $x,y$ and (b) $x,z$ planes. The dashed grey lines represent the magnetopause location (Shue et al. Reference Shue, Chao, Fu, Russell, Song, Khurana and Singer1997).

The list of crossings can be found in supplementary materials. For each crossing, the classification and the physical quantities relevant for the study (normals, dimensionality index, non-gyrotropy index and the main characteristic lengths discussed above for the case crossing) are included.

11. Statistical study of the magnetic hodograms

The previous results about the role of the FLR effects at the magnetopause are now carried out statistically. This study aims to generalize the results obtained from the case crossing studied in the previous section and to estimate the role played by FLRs at the magnetopause.

The database described above has first been used to perform a statistical study on the hodogram shapes, to determine how often linear hodograms are observed in magnetopause crossing. Having an estimation of the percentage of crossings that do not conform to CTD allows us to gauge how frequently the assumptions made by this theory do not accurately represent the magnetopause. For this purpose, we separate the crossings into the following different classes, this classification being based on CTD distinctions and on the preceding findings.

1. (i) Linear crossings, i.e. straight lines not passing through the origin as in the above case study.

2. (ii) Radial crossings, including all linear crossings whose best fit line passes through the origin (considering uncertainty). These crossings correspond to CTD compressional discontinuities.

3. (iii) Circular crossings, when the distance from the origin is constant. These cases correspond to CTD rotational discontinuities.

4. (iv) Other crossings, whose features are not included in the previous classes. This class includes crossings with various features, e.g. circular hodograms not centred on the origin, crossing characterized by two different hodograms in two subintervals, etc., and crossings that do not have an obvious distinction between the previous classes, due to noise.

To classify each crossing, we only focus on its central time interval, where the gradients are maximum. By considering larger time intervals, the hodograms’ shape becomes more complex because the variations out of this interval are generally unrelated to the main boundary jumps. Selecting only the middle part of the crossing provides simpler and more conformal hodograms. Even if the boundary jumps are not fully completed in this part, this will not prevent comparing the experimental results with CTD predictions since this theory, when valid, is based on conservation laws for any subinterval of the discontinuity. When this theory fails to reproduce the observed properties, we can interpret those new features as coming from kinetic effects, therefore confirming the limitation of CTD to describe the magnetopause boundary. To that purpose, for each dataset, we selected the crossing temporal interval following the algorithm used in Haaland et al. (Reference Haaland, Sonnerup, Dunlop, Balogh, Georgescu, Hasegawa, Klecker, Paschmann, Puhl-Quinn, Rème, Vaith and Vaivads2004, Reference Haaland, Reistad, Tenfjord, Gjerloev, Maes, DeKeyser, Maggiolo, Anekallu and Dorville2014) and Paschmann et al. (Reference Paschmann, Haaland, Phan, Sonnerup, Burch, Torbert, Gershman, Dorelli, Giles, Pollock, Saito, Lavraud, Russell, Strangeway, Baumjohann and Fuselier2018b) to estimate the spatial scale of the magnetopause (intervals are identified as 75 $\%$ of the magnetic field $\boldsymbol {B}_L$ component variation).

The classification performed here differs from previous attempts to classify magnetopause hodograms, as seen in studies such as Sonnerup & Ledley (Reference Sonnerup and Ledley1974), Berchem & Russell (Reference Berchem and Russell1982) and Panov et al. (Reference Panov, Artemyev, Nakamura and Baumjohann2011). In these previous works, hodograms were categorized as C-shaped or S-shaped based on their form in the tangential plane. However, unlike those studies, we considered the central part of the crossing, rather than considering the entire temporal interval. Our classification of hodograms involves the following two-step process.

1. (1) Visual inspection: initially, all hodograms are visually inspected to identify the cases that are clearly not linear or circular, which are classified separately as ‘Others’. Additionally, a preliminary distinction is made between crossings with circular and linear features.

2. (2a) Analysis of hodograms with possible circular features: for these crossings, we analyse the variation of the modulus of the magnetic field in the plane, allowing for a maximum possible variation of 20 %. This accounts for factors such as turbulence and waves propagating alongside the magnetopause. Any crossings exceeding this 20 % threshold are categorized as ‘Others.’

3. (2b) Analysis of hodograms with possible linear or radial features: these crossings undergo an initial assessment to confirm their linearity. This involves examining the width-to-length ratio of the crossing, with any ratio exceeding 20 % classified as ‘Other’. Finally, the remaining crossings are classified as either radial or linear based on whether their projection passes through the origin.

From this database, we found the following distribution:

1. (i) 36.3 $\%$ (53/146) of the crossings present linear features;

2. (ii) 2.7 $\%$ (4/146) of the crossings present circular features (rotational discontinuity);

3. (iii) 15.8 $\%$ (23/146) of the crossings present radial features (compressional discontinuity);

4. (iv) 45.2 $\%$ (66/146) of the crossings could not be interpreted definitely as either of the three before (presenting more than one feature at the same time).

It follows that more than a third of the selected crossings show linear features, emphasizing that the fundamental role FLR effects have on magnetopause structure is found in a significant number of crossings.

It could be interesting to compare the above results with the several classifications that were previously published (see Liu et al. Reference Liu, Fu, Cao, Wang, He, Guo, Xu and Yu2022 and references therein). These previous classifications were not based on the analysis of the rotational and compressional properties as done here, but on the normal component of the magnetic field and its magnitude (background and variation) (Smith Reference Smith1973; Burlaga, Lemaire & Turner Reference Burlaga, Lemaire and Turner1977; Tsurutani & Smith Reference Tsurutani and Smith1979; Neugebauer & Giacalone Reference Neugebauer and Giacalone2010). For such a comparison, however, one should take care that there are important differences in the definitions: in these previous classifications in particular, any discontinuity is named ‘tangential’, whatever its other properties, as soon as the measured $B_n$ is sufficiently smaller than $B$, the threshold for this ratio being, for instance, of the order of 0.3 (Smith Reference Smith1973; Burlaga et al. Reference Burlaga, Lemaire and Turner1977; Tsurutani & Smith Reference Tsurutani and Smith1979; Neugebauer & Giacalone Reference Neugebauer and Giacalone2010; Liu et al. Reference Liu, Fu, Cao, Wang, He, Guo, Xu and Yu2022). This is of course a very different approach from the one we use here since, even when $B_n$ is small (and even if barely measurable), we consider that different kinds of discontinuities exist, with different properties.

As done for the case study above, it was possible to study on a statistical basis (i) the ratio between the width of the magnetopause and the ion Larmor radius, and (ii) the non-gyrotropy index. For both parameters, the case study appears rather typical. On average, the magnetopause was found to be approximately 6.5 times the ion Larmor radius, only slightly smaller (6.1) for linear hodograms. Similarly, the non-gyrotropy index $D_{ng, \perp }$ has an average value of 0.07, only slightly higher (0.08) for linear hodograms. The $D_{ng, {\rm Aunai}}$ index has even comparable averages for the four different classes. It therefore seems that, although non-gyrotropy has been demonstrated above to play an important role, the non-gyrotropy index alone is not decisive for predicting unequivocally the shape of the hodograms. This question should be the subject of future works.

12. A comparison between the magnetic and the particles normals

For each crossing, both the magnetic and the particles normals were computed with the GF2 tool. Thanks to the high resolution of the MMS measurements, we can measure the local fluctuations of the normals inside the magnetopause around their mean values. However, in order to compare the magnetic and ion geometries, a single average normal was used for each case. The mean normal is obtained inside the same time interval as in the previous study.

To study the differences between the two normals, we compared them via their departure from the Shue model's normal (where the magnetopause is assumed to be a paraboloid, Shue et al. Reference Shue, Chao, Fu, Russell, Song, Khurana and Singer1997). This normal was obtained using the solar wind and Interplanetary magnetic field (IMF) properties from the OMNI dataset (King & Papitashvili Reference King and Papitashvili2005). The time delay between the crossing time and the measurement time of the solar wind relevant parameters is estimated by using the propagation method used in Michotte De Welle et al. (Reference Michotte De Welle, Aunai, Nguyen, Lavraud, Génot, Jeandet and Smets2022) (which was adapted from Šafránková et al. Reference Šafránková, Nmeček, Dušík, Přech, Sibeck and Borodkova2002). The procedure for acquiring these parameters is as follows: (i) the distance from the bow shock's nose (where OMNI data are defined) to the crossing location, projected along the Earth–Sun axis, is estimated; (ii) we estimate the solar wind's propagation time ($t_{{\rm est}}$) between these two points, assuming an average solar wind velocity of 400 km s$^{-1}$; (iii) the solar wind velocity $V_{{\rm sw}}$ is then determined from the OMNI dataset, averaging over a 2-min interval centred on the crossing time adjusted by the time delay $t_{{\rm est}}$; and (iv) ultimately, a final time delay is computed based on $V_{{\rm sw}}$, which is subsequently utilized to obtain final values of solar wind and IMF parameters. The crossings for which OMNI data computed with this procedure are missing (10 out of 146) were left out of this analysis.

In figure 12(a) we plot the angle between the nominal normal and the magnetic and particles normals, respectively. In this figure we observe that most of the crossings are along the diagonal, corresponding to cases where the two normals, ionic and magnetic, are similar (82 points out of 146 are between the two thin lines, which indicate differences of ${\pm }10^\circ$).

Figure 12. (a) Comparison between the angle between the theoretical normal (Shue et al. Reference Shue, Chao, Fu, Russell, Song, Khurana and Singer1997) and the magnetic and ion normals. (b) Distribution of the angle between the magnetic and ion normals. Here the markers for each point are chosen depending on whether each crossing respects the criteria on dimensionality, stationarity and normal variance on the ions flux measurements (see Appendix C for further details). Colours in the histograms are used accordingly. Blue, green and yellow points indicate the crossing with small variance on the ions normal direction within the crossing, good one-dimensionality and good stationarity. Black points indicate the crossings respecting all criteria, red points not any criteria.

The cases are distributed throughout the plane, with many cases above 40$^\circ$, although we observe a cluster at lower angles, between zero and 30$^\circ$. The largest angles correspond to a magnetopause very far from the paraboloid shape assumed in Shue's model, which relies on the assumption of a magnetopause at (or near) equilibrium. The departures are likely to be related to surface waves on the boundary itself.

Finally, the distribution of the angles between the two normals is shown in figure 12(b). Here we evidence again that most of the cases studied (82 out of 146) are below 20$^\circ$, with the maximum of the distribution at 10$^\circ$. However, we also observe again that several cases have much larger angles, up to 90$^\circ$. The strongest departures are problematic and deserve further investigation. This appears to be due to the more complex ion structure with respect to the magnetic one. As the criteria used for the dataset selection were built from magnetic data, they are not as relevant when considering ion normals. This is evidenced in figure 12, where the colours indicate how several ion criteria are satisfied. These criteria concern respectively the dimensionality, the stationarity and the variance of the normal direction. All details are given in Appendix C. Focusing on points respecting all the criteria for the ions flux (black markers and hodogram), we observe that only a few crossings are outside the diagonal. Only two of these crossings have angles above 40$^\circ$.