Journal cover Journal topic
Annales Geophysicae An interactive open-access journal of the European Geosciences Union
Journal topic
ANGEO | Articles | Volume 38, issue 2
Ann. Geophys., 38, 275–286, 2020
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.
Ann. Geophys., 38, 275–286, 2020
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.

Regular paper 02 Mar 2020

Regular paper | 02 Mar 2020

A multi-fluid model of the magnetopause

A multi-fluid model of the magnetopause
Roberto Manuzzo1,2, Francesco Califano2, Gerard Belmont1, and Laurence Rezeau1 Roberto Manuzzo et al.
  • 1LPP, CNRS, Ecole Polytechnique, Sorbonne Université, Univ. Paris-Sud, Observatoire de Paris, Université Paris-Saclay, PSL Research University, Paris, France
  • 2Department of Physics E. Fermi, Universitá di Pisa, Pisa, Italy

Correspondence: Francesco Califano (

Back to toptop

Observation of the solar wind–magnetosphere boundary provides a unique opportunity to investigate the physics underlying the interaction between two collisionless magnetized plasmas with different temperature, density and magnetic field topology. Their mixing across the interface as well as the boundary dynamics are affected by the development of fluid (and kinetic) instabilities driven by large-scale inhomogeneities in particle and electromagnetic fields. Building up a realistic initial equilibrium state of the magnetopause according to observations is still a challenge nowadays. In this paper, we address the modeling of the particles and electromagnetic field configuration across the Earth's magnetopause by means of a three-fluid analytic model. The model relies on one hot and one cold ion population as well as a neutralizing electron population. The goal is to create an analytic model that is able to reproduce the observations as closely as possible. Some parameters of the model are set using a fitting procedure that aims to minimize their difference with respect to experimental data provided by the Magnetospheric MultiScale (MMS) mission. All of the other profiles, concerning the electron pressure and the relative densities of the cold and hot ion populations, are calculated in order to satisfy the fluid equilibrium equations. Finally, using a new tri-fluid code, we check the stability of the large-scale equilibrium model for a given experimental case and provide proof that the system is unstable to reconnection. This model could be of interest for the interpretation of satellite results and for the study of the dynamics at the magnetosphere–solar wind boundary.

1 Introduction
Back to toptop

The solar wind–magnetosphere boundary, known as the magnetopause, is characterized by the presence of magnetic and velocity shears as well as jumps in magnetic and velocity magnitudes in addition to jumps in plasma density and temperature. These inhomogeneities are the sources of many plasma instabilities at different spatiotemporal scales (Labelle and Treumann1988) that, in turn, often trigger secondary instabilities at smaller scales. As an example, secondary instabilities such as magnetic reconnection, Kelvin–Helmholtz and/or Rayleigh–Taylor instability can efficiently develop on the shoulder of the primary instability; the Kelvin–Helmholtz instability at the low-latitude magnetopause is a prime example of such a case (see Faganello and Califano2017 and references therein).

All of these phenomena can cause significant entry of magnetosheath plasma mass (Paschmann1997), momentum (Dungey1961) and energy (Lee and Roederer1982) into the magnetosphere. The study of the magnetopause is of particular interest, as this system offers the unique opportunity to study a two-plasma large-scale interaction in conditions not achievable in the laboratory. The magnetopause physics is also of basic importance in studies addressing the Sun–Earth interaction, in particular concerning the impact of solar wind disturbances on the terrestrial environment and space-weather forecasting efforts (see for instance Baker and Lanzerotti2016). Modeling space plasmas using data provided by multi-spacecraft missions has been much developed during (and in preparation of) the “cluster era” (Büchner et al.1998). Concerning the magnetopause data, one of the key points relates to the mixing between magnetospheric and magnetosheath plasmas and the resulting non-Maxwellian shape of the distribution functions (hereafter referred to as “dfs”) observed in these regions (Bosqued et al.2001; Frey et al.2003; Phan et al.2005; Retinò et al.2005). These dfs are often reminiscent of those observed in reconnection kinetic simulations (Nakamura and Scholer2000; Tanaka et al.2008; Aunai et al.2011). Some of these dfs can be compared using simple analytic models, as in the pioneering work by Cowley and Owen (1989). As the populations originating from the two different sides of the magnetopause differ with respect to density and temperature, modeling the mixing requires the use of at least a multi-population model. In the perspective of investigating the dynamics of the magnetopause mixing layer using a three-fluid numerical simulation, the main target of this paper is to build up a three-fluid equilibrium that is as realistic as possible to initialize the simulation.

The three-fluid model that we propose in the present paper aims to distinguish between ions of magnetospheric and magnetosheath origin. It must be considered as a first step towards a four-population model, which would also allow one to distinguish between the electrons from both origins. This simplification allows nonessential details related to the presence of two electron species, which are not likely to have a major role in structuring the large-scale equilibrium, to be avoided. In addition, note that for a correct modeling of the electrons in the magnetopause vicinity, one should also split the magnetospheric electrons into at least two subpopulations: one “cold” poorly measured population, carrying the density, and one “hot” population, carrying the pressure. Finally, starting from a previously existing two-fluid code, it is relatively straightforward to create a multi-ion one-electron code where the electrons only provide an “Ohm's law”. However, a fully multi-population code (i.e., including several electron populations) requires a radically different approach that is presently under investigation and is a matter for future work.

In the past, several multi-population models trying to simulate the plasma interaction between the magnetosheath and magnetosphere have been developed. In particular, for the modeling of tangential layers (i.e., without normal magnetic field and without normal velocity, Bn=un=0), the kinetic models are by definition the most complete ones as they rely on a Vlasov equilibrium. Nevertheless, as these models are considerably more complicated than the fluid models, using the latter can be considered to be a necessary first step for modeling the magnetopause as closely to the observational data as possible. Note also that all of the equilibria built via dfs that are only single-valued functions of the particle invariants of motion (Channell1976) cannot really be considered to be “multi-population” models: they ignore the issue of accessibility (Whipple et al.1984) and, therefore, cannot distinguish between particles of magnetospheric or magnetosheath origin. Some models (Roth et al.1996; De Keyser and Roth1998) have partly accounted for the accessibility problem by considering functions of the invariants that are bi-valued. In recent years, new models have been developed by Belmont et al. (2012) and Dorville et al. (2015) that account for accessibility by considering more general multivalued dfs. Nevertheless, it must be kept in mind that all of these models involve many free parameters. In Dorville et al. (2015), the authors even emphasize that there is no general constraint for fixing the normal electric field profile in their kinetic model; thus, the normal electric field profile can be chosen in a fully arbitrary manner. Furthermore, even if a few profiles can be fixed in a more or less realistic way, all of the others depend on simple mathematical assumptions, which are largely arbitrary; therefore, the latter profiles generally still remain far from realistic.

A different approach to creating a large-scale configuration at the boundary between the magnetosheath and the magnetosphere is to run a global simulation that eventually reaches a steady state. This has been done recently utilizing hybrid and kinetic codes (Chen et al.2017; Palmroth et al.2018), although the question of the spatial resolution is then difficult to solve for investigating the internal structure of the magnetopause boundary. It will be interesting to compare the results obtained using the abovementioned methods with those obtained using the methods discussed in the present paper where the model instead relies on spacecraft data as the initial condition for the simulation.

In summary, the lack of realistic equilibria in the literature generally makes the initialization of the magnetopause studies difficult with respect to kinetic simulations. Recently the multi-population character of the medium has been taken into account by Dargent et al. (2017), although for a different purpose from that in this study. Using a particle-in-cell (PIC) code, they investigated the influence of cold magnetospheric ions on the development of magnetic reconnection. In this paper, the magnetospheric plasma includes two populations with different temperatures in order to account for the presence of the cold ions that are observed in the magnetosphere close to the magnetopause.

Multi-fluid models have been developed in various domains, although generally not for magnetopause studies. These studies address multispecies evolution involving chemical processes and collisions. They have been used to investigate planetary atmospheres (Modolo et al.2006; Ma et al.2007), the solar chromosphere (Alvarez Laguna et al.2016) and basic plasma physics problems (e.g., drift turbulence in Shumlak et al.2011).

In this paper, we present a new technique to create a three-fluid equilibrium that derives directly from satellite observations. The model assumes 1-D gradients in the normal direction and a tangential boundary (Bn=0) at the magnetopause. The magnetic and velocity shear are both taken into account in a realistic way. Profiles are chosen that best fit the data from the Magnetospheric MultiScale (MMS) mission (Burch et al.2016), for which the time-to-space conversion has been performed using recent techniques presented in Manuzzo et al. (2019). As will be shown in Sect. 4, the method provides a cold and a hot contribution that are in qualitative agreement with observations, even if the model only uses the global ion macroscopic moments as inputs.

2 Observations
Back to toptop

We use MMS data from 16 October 2015 at 13:05:34 + 40 s UT that occurred at [8.3, 8.5, −0.7] RE in GSE frame, which embeds a magnetopause crossing. In Fig. 1, we plot the experimental data that the equilibrium model attempts to reproduce. This interval shows the standard signatures of the region where the magnetospheric and magnetosheath plasmas meet (magnetopause crossing): a reversal of the magnetic field and a change in the energy distributions.

Figure 1MMS data for the event on 16 October 2015 at 13:05:34 UT + 40 s. (a) Normal and tangential (to the magnetopause plane) components of the magnetic field. (b, c) Ion and electron spectrograms. The first abscissa is the spatial coordinate normal to the magnetopause Xn (see text), and the second abscissa is time. (d, e, f) Ion distribution functions (idfs) recorded by the Fast Plasma Investigation (FPI) instruments in the magnetosphere, in the overlapping region and in the magnetosheath, respectively. These idfs are projected on the tangential plane by integration over the normal component of the velocity.


In Fig. 1a–c, data are plotted as functions of a spatial coordinate Xn, which is the projection of the spacecraft path along the direction normal to the magnetopause (units of di,MSh with di,MSh70 km). Once the function Xn(t) is known, it allows for the recovery of the spatial profiles of all quantities of interest from the corresponding temporal series of the spacecraft magnetosheath–magnetosphere crossing, even with a variable velocity. For the sake of completeness, we also give the time corresponding to each given value of Xn in the abscissa of Fig. 1c. The spatial coordinate X(t) is determined by the integration of the velocity of the spacecraft with respect to the magnetic structure. In order to determine this spatial velocity, we use a new method from Manuzzo et al. (2019), which is an extension of the spatiotemporal difference (STD) method first proposed by Shi et al. (2006). The main difference with respect to the Shi model consists of taking the possible weak nonstationarities of the magnetic structure into account , i.e., assuming t,scB=t,0XB+t,0B and not just t,scB=t,0XB (the indexes “sc” and “0” indicate the spacecraft frame and the frame where the intrinsic variation of the magnetopause structure is the smallest, respectively). This allows for some singularities in the results to be avoided.

Black points have been superimposed on the two spectrograms (Fig. 1b, c) to indicate their maxima. This allows one to more easily distinguish where the magnetosheath and the magnetospheric plasma interact, as indicated by discontinuities in the curve joining the maxima. The region where the two plasmas partially overlap in space is emphasized by a blue rectangle in Fig. 1b, which is centered at Xn∼3 di.

In Fig. 1d–f, we plot the 2-D ion distribution functions (idfs) in the plane tangential to the magnetopause. They are obtained by integration over the out-of-plane (normal) component of the velocity. Each plot is the average of five single idfs recorded within a ∼0.75 s long interval (equivalent to 0.5di). The radius of the distribution functions is 103 km s−1 and the purple full circle drawn at their centers determines the bottom limits of the energy of the FPI instrument (10 eV  53 km s−1 for ions). The direction of the local magnetic field is indicated by a white arrow. In Fig. 1e (mixing region), one can observe that the idf contains two peaks that are emphasized by the blue and red dashed circles. These two circles have a diameter equal to the magnetosheath and magnetospheric thermal velocities, respectively. The same circles are shown for the magnetosheath and magnetospheric idfs in Fig. 1f and d. In these two asymptotic media, we see that there is only one single peak. Note that the idf shown in Fig. 1d was recorded a little earlier (10:20:00 UT + 2 s) during a “clear” observation of the magnetosphere and, therefore, allows the presence of magnetosheath particles when the spacecraft is too close to the magnetopause to be avoided. Conversely, Fig. 1e shows a mixture of the magnetosheath and magnetospheric populations at the same time. However, as the two peaks are close to each other and the distributions of the two populations are partly superposed, it is not possible to clearly separate the hot/cold contributions, which is a necessary input for the multi-population model to be built by a direct fit of this region. In the next section, we will explain a new method capable of separating the two particle components – even in such complex situations. Note that the electron spectrogram of Fig. 1c shows energy maxima that lay just at the bottom limit of the instrument (≃10 eV) in the magnetospheric region (Xn≤2di). This indicates the presence of cold electrons in the magnetosphere. The role of this poorly measured cold electron population is not relevant for the magnetopause pressure equilibrium, but it could be significant in the electron bulk velocity. However the electron population parameters will be not determined by a direct fit of the data nor will those of the two ion populations (cold/hot). Instead, they will be determined by another method based on the equilibrium equations, which we will describe in the next section.

3 The three-fluid model
Back to toptop

3.1 Equilibrium equations

We present a three-fluid collisionless model here that includes two proton populations (one cold and one hot) and one electron population. The cold ion population models the ions of magnetosheath origin and disappears more and more on the magnetospheric side. Conversely, the hot population models the ions of magnetosphere origin and disappears on the magnetosheath side.

The continuity and ion momentum equations are derived from the first two moments of the Vlasov equation. We impose charge neutrality, and the displacement current is neglected. We assume isotropic pressures and adopt a polytropic closure for all populations. These equations are coupled to the electromagnetic field via the Faraday equation, and we use an Ohm's law that takes the electron pressure gradient into account but neglects electron inertia effects. The three-fluid system of equations reads as follows:


where the index α stands for “all populations”, and the index β represents “all ion populations”.

3.2 Determination of the fluid profiles

We aim to establish a tangential 1-D equilibrium to mimic the magnetopause observations previously presented as closely as possible. Assuming /t = 0, this is done in three steps.

In step 1, we impose the magnetic field B as well as the density ni, temperature Ti and velocity Ui profiles by fitting the data using a combination of hyperbolic tangents, as explained in Sect. 4. At this stage, the ion quantities are chosen without distinguishing between the cold and hot populations.

In step 2, we deduce the electron density ne and velocity Ue using Eqs. (1a) and (1b), respectively. The temperature Te is deduced from


where the total pressure, Ptot, has been imposed to be constant in order to fulfill the equilibrium conditions.

As far as Pe is concerned, we note that Pe is much smaller than Pi+PB (see Fig. 2) and that it is difficult to estimate it precisely due to experimental uncertainties. As a consequence, the constant value of Ptot used in the model is taken as the maximum of the measured total pressure Pi+PB+Pe. Note that the total pressure is not really constant due to the abovementioned uncertainties on Pe and because the experimental structure may be slightly out of equilibrium. The modeled Pe is deduced from Eq. (2), with the choice for Ptot ensuring only positive values for Pe. Finally, the electric field E is deduced from the Ohm's Law, as in Eq. (1g).

Figure 2Comparison between Pe, Pi and PB. The total pressure is indicated as Ptot. Numerical values are all normalized to the magnetosheath pressure PMSh∼3200 eV cm−3.


In step 3, we now split the global proton population into two different populations, cold and hot (hereafter referred to as “ic” and “ih”, respectively) to distinguish the magnetospheric and magnetosheath populations. The densities (nic and nih), pressures (Pic and Pih) and currents (Jic and Jih) of the two ion populations add to form the total ion density, pressure and current as follows:


The temperatures of the cold and hot ion populations, Tic and Tih, are assumed to be constant. As the global ion temperature profile Ti is known, their values are obtained from the satellite data using the following two limits:


The temperature ratio between the two populations is set by the value of the dimensionless parameter as follows:


Using Eq. (3b), the contributions of each population to density and pressure are fully determined by the Ti profile and the temperature ratio Υ:


The perpendicular currents and, thus, the corresponding velocities are fully determined by Eq. (1d). On the contrary, the parallel currents cannot be determined by the abovementioned system of equilibrium equations. Therefore, we will set them using a reasonable choice for the parameter ϕ which is equal to the ratio of the cold parallel ion current to the total parallel ion current as seen in the electron frame:


The parallel components of the hot and cold ion velocities can have opposite directions; thus ϕ is defined in the [-1,1] range, whereas Γ and Π are defined in the [0,1] range. The reasonable choice for ϕ is suggested by the data and will be discussed in more details in the next section. In general, the asymptotic values of the cold and hot ion currents are chosen in agreement with the asymptotic values of nic and nih so that all of the corresponding values of the velocities Uic and Uih have reasonable values, although one of the two densities nic or nih tends to nearly zero on each side.

In order to implement this model into a numerical simulation, a compromise is necessary because the multi-fluid code cannot deal with a population having a zero density somewhere in the domain. To avoid this problem, we introduce the parameters ϵ(c)≪1 and ϵ(h)≪1, and we modify the initialization so that the cold and hot densities tend to ϵ(h)ni and (1−ϵ(h))ni on the magnetospheric side and to (1−ϵ(c))ni and ϵ(c)ni on the magnetosheath side. The temperatures are changed according to the following equations:


where Tic and Tih indicate the observed values corresponding to the model, and TiMSph and TiMSh indicate the temperatures corresponding to the magnetospheric and magnetosheath values of Ti. A similar correction is made for the ion velocities (see the following sections).

4 Data vs. analytic profiles
Back to toptop

We apply the procedure to the case study introduced in Sect. 2. In Fig. 3, we compare the model field profiles with those obtained with the MMS data. The model profiles for the magnetic field, the ion temperature and density are obtained by a fit procedure (Fig. 3a, b, and c, respectively), whereas the others are calculated from the equilibrium equations. The fits are obtained by means of analytic functions. For a given quantity Q, the fitting functions have the following form:


where Xn is the coordinate along the direction normal to the magnetopause (as discussed in Sect. 2). The parameters aQ,j, bQ,j, cQ,j and dQ,j are the free parameters shaping the analytic profiles, and j is the component index. The maximum value of j depends on the fitted quantity: the analytic profiles are considered to be good fits of the data if they correctly shape the large-scale configuration as well as the position and the length scale of the gradients within the magnetopause layer. An example of a “good fit” is given in Fig. 3a–c. It is worth noticing that the particle boundary, observed on density and temperature, has a length scale smaller than the magnetic boundary (by a ratio ≃0.25) and that its position is considerably shifted toward the magnetosphere with respect to the center of the magnetic jump. This may indicate the presence of a boundary layer, possibly made of magnetosheath plasma, that is observed on the magnetospheric side of the magnetopause (Hasegawa2012). Such features cannot be reproduced in the framework of a magnetohydrodynamic (MHD) equilibrium model. To the best of our knowledge, these features have not even been introduced in the context of a kinetic model.

In Fig. 3b, we show the temperature profiles as obtained with our model equilibrium. The total ion population temperature Ti has been obtained by a fit, and it is superposed onto the cold ion population temperature Tic (blue curve) and the hot ion population temperature Tih (red curve). The figure has been drawn using ϵ(h)=0.35 and ϵ(c)=0.05, which determine the values of Tic and Tih via Eqs. (7a) and (7b), respectively.

One observes that the global temperature is well fitted by the model outside the mixing region, but the fit is less accurate in the 1.0diXn2.5di interval. In this interval, the real total ion temperature actually becomes larger than its magnetospheric asymptotic limit. Unfortunately this feature can not be reproduced by the present three-fluid model with hot and cold constant temperatures, as the Γ≥0 constraint forces the Ti profile to be lower than Tih everywhere (see Eq. 5a). This little deviation is acceptable as the model mainly aims to reproduce the asymptotic trends, with the observed inner region most likely being out of equilibrium.

In Fig. 3c, we show the density profiles. As explained in the previous section, the hot and cold ion contributions to the total density ni are computed by means of the Γ function, which is fixed once the global Ti profile and the temperature ratio Υ are fixed (Eq. 5a). In Fig. 3a–g, the two vertical (dashed black) lines indicate the limits of the region where 1/4Γ3/4. Note that the cold ion density rapidly falls to very low values in more or less ∼2 di (where di≃70 km is the ions inertial length measured in the magnetosheath), whereas the hot population density remains at nearly the same value over a longer interval (between 0 and 8 di).

The electron density and velocity profiles are obtained from the equilibrium equations. However, these quantities are not plotted here because their experimental counterparts are likely to be biased in the magnetosphere by the cold electron population, which is below the bottom energy threshold of the FPI instrument (as mentioned in Sect. 2). Conversely, in Fig. 3d, we plot the normal component of the electric field, which is obtained from the three-fluid model using Eq. (1g). We see that the electric field calculated by the model agrees quite well with the experimental field. As the electric field reflects the electron dynamics, this shows (independently of the electron measurements) that the electron dynamics is correctly accounted for in the model.

Figure 3Comparison between the magnetopause profiles as observed by MMS on 16 October 2015 at 13:05:34 + 40 s UT and those used for the three-fluid model equilibrium. Satellite data are represented by dashed lines, and the extrapolated profiles used in the model are represented by continuous lines. The Xn coordinate represents the spatial coordinate normal to the magnetopause Xn. In the panels we show the magnetic field (a), the temperatures (b), the density (c), the normal component of the electric field (d), the parameter Γ, Π and Φ (e), and the parallel and perpendicular components of the ion current (f and g, respectively) as computed from particle data. The two vertical (black dashed) lines highlight the transition region (1/4Γ3.4). The blue and orange lines adopted for the electric and magnetic fields represent the normal and the tangential components of the fields. The square roots of the temperatures in panel (b) are plotted in velocity units in order to make the comparison with the idfs shown in Fig. (1) easier. For the sake of clarity, the curves shown in panels (f) and (g) have been multiplied by a factor of 10 in the 0.0Xn3.0 interval.


The parallel components of the cold and hot ion currents are set by ϕ (Eq. 6). As long as there are no cold ions on the magnetospheric side and no hot ions on the magnetosheath side, the asymptotic constraints on ϕ would be


Nevertheless, owing to the compromises necessary to implement the model in the multi-population numerical simulation, the cold and hot densities actually take small but not strictly null values on both sides. To determine the corresponding parallel currents, corrections similar to Eqs. (7a) and (7b) are applied with the assumption that, on each side, the parallel velocities of the cold and hot populations, in the electron frame, are equal to each other and, therefore, equal to the global one. Under this assumption, it can be easily shown that the asymptotic values of Φ are equal to those of Γ:


Note that for the particular MMS event considered, the global ion parallel velocities are observed to be quasi-null on each side, so that the same asymptotic property holds for the velocities of the two populations.

Between the two limits above, a reasonable choice for the ϕ profile is that the length of its gradients is of the same order as the scale length of the density and temperature gradients, i.e.,  1–2 di. The position of the main gradient of ϕ is set in order to separate the magnetopause thickness into two parts, each representing the length proportional to the gyro-radii of the two populations (their ratio is ≃2).

In Fig. 3e, we show the model profiles for Γ, Π and Φ. Because of the differences in temperature between the two components, the profile in Π (concerning the pressures) noticeably differs from the profile in Γ (concerning the densities).

Finally, in Fig. 3f and g, we show the results concerning the parallel and perpendicular components of the ion currents. Once more, one observes that the global ion current is well fitted with the exception of the perpendicular current in the mixing region, which is less accurate. This is due to the previously mentioned small inaccuracy in the modeled ion temperature in this region.

5 Numerical simulations
Back to toptop

5.1 Setup

Here we give an example of a three-fluid numerical simulation with the aim of demonstrating the possibility of numerically studying the above system by starting from an equilibrium not far from a real one (not only qualitatively but also quantitatively). For the sake of simplicity, we limit the geometry (to 2-D) and the numerical code (to 3-D). A detailed numerical study relying on such an approach will be the focus of future work.

The three-fluid model introduced in this paper has been used to initialize a 2-D three-fluid numerical simulation of the interaction between the solar wind and the Earth's magnetopause. The numerical simulation is intended to mimic the MMS crossing on 16 October 2015 at 13:05:34 + 40 s UT. This simulation has been performed using a three-fluid numerical code that solves Eqs. (1a)–(1g). The code originates from a two-fluid 3-D parallel code largely used for the study of the interaction of the solar wind with the magnetosphere (see Fadanelli et al.2018 and references therein). The three-fluid code adapts the new equations to the algorithm of the two-fluid code presented in Faganello et al. (2009). It advances in time using a standard third-order Runge–Kutta algorithm (Canuto1988). It uses sixth-order explicit finite differences along the periodic y direction and a sixth-order compact finite difference scheme with spectral like resolution for spatial derivative along the inhomogeneous x direction (see Lele1992 for the significance and technical details of compact finite differences with spectral like resolution). The numerical stability is guaranteed by means of a spectral filter along the periodic y direction and a spectral-like filtering scheme along the inhomogeneous x direction. The code is parallelized along the periodic y direction (Lele1992). The code has been validated by standard numerical tests. In particular, by separately selecting the two cold and hot ion populations, we have reproduced the propagation of ion acoustic and Alfvén waves.

To initialize the simulation presented in this paper, we take the model profiles represented in Fig. 3 as the initial equilibrium, including the few modifications for the cold and hot ion density components (with respect to the basic model) for computational reasons that are discussed at the end of Sect. 3. The three-fluid code handles a normalized version of the three-fluid set of equations (Eqs. 1a1g) using the following characteristic quantities. Magnetic field and ion densities are normalized to their mean values in the magnetosheath (out of the magnetopause layer) B0=BMSh and n0=nMSh. Velocity is normalized to the corresponding Alfvén velocity V0=VA,MSth=B0/(μ0n0mi), and time is normalized to the corresponding inverse ion gyro-frequency t0=1/Ωc,MSh=mi/qiB0. The other normalization quantities are obtained from those just introduced: the lengths are normalized to l0=V0t0 (i.e., the ion inertial length), the temperatures are normalized to T0=miV02 and the electric field is normalized to E0=V0B0.

The simulation box dimensions are given by Lx=160di and Ly=20πdi, and the box is discretized using nx=800 and ny=320 grid points corresponding to dx=dy=0.2di. We have checked that the equilibrium configuration remains stable for several thousand ion cyclotron times in the absence of an initial perturbation due to the very low values of the numerical noise and the high accuracy of the numerical method.

5.2 Results

The large-scale equilibrium configuration used to initialize the simulation is unstable with respect to the reconnection mode. At t=0 we add an initial perturbation δB=×A to the equilibrium. The potential vector is given by a sum of random phase modes as follows:


where k=kx2+ky2 and ϕ[0,2π) are random phases, and a(x) is a Gaussian-like convolution profile in the inhomogeneous direction going to zero at both boundaries; a(x) is given by the following equation:


where xmp=Lx/2=60 and Lmp=1.66 are the position and the thickness of our magnetopause model, respectively, and a0 is a small amplitude such that the maximum absolute value of the initial perturbation is 10-3.

Figure 4Δ vs. the wave number my for an equilibrium magnetic field ∼tanh(x) (dimensionless). The unstable modes are those associated with positive values of Δ.


Figure 5Development of the reconnection instability. (a) The modulus of the Fourier transform of δbx along y vs. x at five fixed time instants. The plots correspond to the fastest growing mode, m=2, on a log scale. The two red dashed vertical lines indicate the space interval from Fig. 3. (b) The first five eigenmodes' growth vs. time. The orange curve corresponds to the most unstable mode, m=2, which is plotted in panel (a). (c) The growth rate values vs. ky calculated by a best fit of the slopes in panel (b); the colors also correspond to those used in (b).


The simulation is run for about 1500 ion cyclotron times. Very rapidly, the initial perturbation reorganizes and sets up the reconnection eigenmodes that are identified by their wave number in the y periodic direction (each m wave number is easily recovered by taking the Fourier transform of the perturbation along the y direction at a given time). Following the classical reconnection theory (Furth et al.1963) (but ignoring the density inhomogeneity), we have checked that our equilibrium is Δ unstable for the first five eigenmodes. We recall here that the Δ parameter depends on the equilibrium magnetic shear and the wavelength of the perturbation. It defines the instability threshold condition (Δ0). The unstable modes can be deduced by looking at Fig. 4 where we plot Δ as a function of the wave number my. We see that only the first five modes my=1,,5 correspond to a positive value of Δ, which is in agreement with the simulation where the my≥6 are stable in the linear phase (see discussion below). In Fig. 5a, we plot the profile of the fastest growing eigenmode (corresponding to m=2) of the x component of the magnetic field fluctuation δbx. The plot is along the inhomogeneous x direction at five different times (see the legend) on a log scale. The two red vertical dashed lines indicate the spatial window of the equilibrium represented in Fig. 3. This shows that after an initial transient that is needed to set up the normal mode shape, the reconnection instability develops around the region where the magnetic field reverts, 0x16 (see also Fig. 3). As the equilibrium is asymmetric, in particular regarding the cold and hot ion density that vary in a different location with respect to the point where the magnetic field inverts, the eigenmode is not symmetric considering the point where the magnetic field inverts, Xn≃6.4. To the best of our knowledge, this is the first time that the reconnection instability has been investigated in the framework of a three-fluid approach in a nonsymmetric equilibrium directly representing the large-scale configuration taken from a satellite data event. Our goal here is to present the possibility of setting up a “realistic” large-scale initial equilibrium configuration to be simulated by a three-fluid approach. The nonlinear development of the system and, in particular, the mixing efficiency will be the subject of future work. Still, in Fig. 5b, we plot the eigenmodes growth vs. time in normalized units (log scale). We see the exponential growth of the first five modes, my=1,,5. Modes with my=6,7 are instead stable. The orange curve corresponds to the most unstable mode, my=2, which is plotted in panel (a). Despite the strong inhomogeneity of the system where, as discussed before, the magnetic inversion and the density variations occur at different locations, we see very clear exponential growth with a constant growth rate. The linear phase lasts until about t≃1000, after which the nonlinear phase begins. The growth rate values of the five unstable modes are reported in Fig. 5c, confirming the theoretical trend shown in Fig. 4. In particular, my=2 is the most unstable mode. In Fig. 6, we show the shaded iso-contours of the cold ion population, Ni,c, at the beginning of the saturated phase, t=1455. We see the formation of a hole structure corresponding to the region where the cold ion density grows and eventually reaches the asymptotic magnetosheath value. To show the cold density hole, we cut along the inhomogeneous x direction at y=38 (see dashed line in the bottom frame of Fig. 6). In Fig. 7, we show the same quantities for the hot ion fluctuations. We see “complementary” behavior in the sense that a bump is now generated more or less in correspondence with the cold ion hole. However, as already discussed, it is not the goal of this paper to study the nonlinear dynamics and mixing properties of the cold and hot ion populations. Our aim here is limited to presenting a method able to obtain a “realistic” three-fluid equilibrium starting from a set of satellite data that can be used as an initial condition for the investigation of the dynamics in the framework of a three-fluid approach.

Figure 6Shaded iso-contours of the cold ion fluctuations, Ni,c-Ni,c(t=0) at t=1455Ωi-1. The bottom panel shows a plot of the same quantities vs. x at y≃39, corresponding to the horizontal dashed line in the shaded iso-contours. Numerical values are normalized to the magnetosheath density NMSh∼10 cm−3.


Figure 7Same as in Fig. 6 but for the hot ion fluctuations, Ni,h-Ni,h(t=0).


6 Conclusions
Back to toptop

The huge number of spacecraft data available today offers a lot of information about the magnetopause, especially those from the MMS mission due to their high-resolution particle data. Thus, magnetopause modeling can now be improved in view of these observations, which show that this boundary is never the simple textbook boundary generally considered. Beyond the natural asymmetry in temperature and density between the magnetosphere and magnetosheath plasmas, the first important ingredient to consider is the strong velocity shear that arises at the boundary as well as the magnetic shear which is a defining property of the magnetopause. Furthermore, the gradients concerning the particles and those concerning the magnetic field generally have different locations and show different scale lengths. Therefore, the model also has to be able to take these characteristics into account.

In this paper, for the first time, we present a three-fluid equilibrium directly derived from data using a magnetopause crossing by MMS. The derivation of the model is based on a fit of the experimental data to the most reliable data, which is completed by a “realistic” solution of the equilibrium fluid equations for the others. The relative densities of the hot and cold ion populations calculated using the equilibrium equations provide an a posteriori check of our three-fluid model. In particular, this information helps to understand the different bulk quantities observed in the ion distribution functions (see Fig. 1d).

Furthermore, a preliminary study shows that the model can be implemented in a three-fluid numerical simulation, validating the correctness of the equilibrium solution. The detailed study of the long time evolution of the magnetopause instability will be the subject of future work.

It may seem contradictory to consider the data as characteristic of some magnetopause equilibrium and observe afterward that this equilibrium is not stable and should not last for long (even if the reconnection phenomenon is never “immediate”). To justify this point, one must understand that the main characteristics that are taken into account are the asymptotic values on each side and, in particular, the velocity shear between the magnetosheath and magnetosphere. These conditions are not changed by the instability. On the contrary, the positions and the scale of the different gradients can indeed be partly modified by the instability. We think that this is one of the interesting issues that can be investigated by the time evolution observed in the simulation. However the following questions are issues that will be resolved in future work: how is the system stability impacted? (A parametric study is needed.) How does the system change in time due to nonlinear effects? Will the simulation converge toward a new more stable equilibrium state representative of the real system?

Investigating the magnetopause stability and trying to understand topics such as when and where reconnection phenomena can be triggered and how the plasmas from both sides can be mixed are still currently challenging issues to attack using numerical simulations. However, knowing that the stability of a physical system is given by the specific initial equilibrium state, it must be kept in mind that the resulting nonlinear dynamics, in particular the mixing properties, also strongly depend on the choice of the initial equilibrium. As a consequence, it is very important to initialize a simulation with a configuration that is as realistic as possible. In most of the published literature, simulations have been initialized with relatively simple configurations, such as Harris sheets or modified Harris sheets, with little relationship to the real magnetopause. Therefore, the realistic three-fluid equilibrium presented in this paper should allow for this work to be taken a step further, and the same method could be applied to other experimental cases in the future.

Data availability
Back to toptop
Data availability. 

All the data used are available from the MMS Science Data Center: (MMS Science Data Center2019).

Author contributions
Back to toptop
Author contributions. 

RM, LR and GB performed the experimental work. RM and FC developed the numerical code and carried out the simulation. All authors contributed to the final paper.

Competing interests
Back to toptop
Competing interests. 

The authors declare that they have no conflict of interest.

Back to toptop

We acknowledge ISCRA for awarding us access to the supercomputer facilities at CINECA, Italy, where the calculations were performed. The French involvement in MMS is supported by CNES and CNRS.

Financial support
Back to toptop
Financial support. 

This project (Francesco Califano) has received funding from the European Unions Horizon 2020 research and innovation programme under grant agreement no. 776262 (AIDA,, last access: 19 February 2020).

Review statement
Back to toptop
Review statement. 

This paper was edited by Minna Palmroth and reviewed by Johan De Keyser and one anonymous referee.

Back to toptop

Alvarez Laguna, A., Mansour, N., Lani, A., and Poedts, S.: MULTI-FLUID MODELING OF THE EARTH'S MAGNETOSPHERE, in: Comparative Heliophysics Program, Summer 2016 at NASA Ames Research Center, available at: (last access: February 2020), 2016. a

Aunai, N., Belmont, G., and Smets, R.: Proton acceleration in antiparallel collisionless magnetic reconnection: Kinetic mechanisms behind the fluid dynamics, J. Geophys. Res.-Space, 116, A09232,, 2011. a

Baker, D. N. and Lanzerotti, L. J.: Resource Letter SW1: Space Weather, Am. J. Phys., 84, 166–180,, 2016. a

Belmont, G., Aunai, N., and Smets, R.: Kinetic equilibrium for an asymmetric tangential layer, Phys. Plasmas, 19, 022108,, 2012. a

Bosqued, J. M., Phan, T. D., Dandouras, I., Escoubet, C. P., Rème, H., Balogh, A., Dunlop, M. W., Alcaydé, D., Amata, E., Bavassano-Cattaneo, M.-B., Bruno, R., Carlson, C., DiLellis, A. M., Eliasson, L., Formisano, V., Kistler, L. M., Klecker, B., Korth, A., Kucharek, H., Lundin, R., McCarthy, M., McFadden, J. P., Möbius, E., Parks, G. K., and Sauvaud, J.-A.: Cluster observations of the high-latitude magnetopause and cusp: initial results from the CIS ion instruments, Ann. Geophys., 19, 1545–1566,, 2001. a

Büchner, J., Kuska, J.-P., and Wiechen, H.: Numerical Modelling and Simulation for Multi-Spacecraft Data Analysis: Approaches and Examples, ISSI Scientific Reports Series, 1, 449–478, 1998. a

Burch, J. L., Torbert, R. B., Phan, T. D., Chen, L.-J., Moore, T. E., Ergun, R. E., Eastwood, J. P., Gershman, D. J., Cassak, P. A., Argall, M. R., Wang, S., Hesse, M., Pollock, C. J., Giles, B. L., Nakamura, R., Mauk, B. H., Fuselier, S. A., Russell, C. T., Strangeway, R. J., Drake, J. F., Shay, M. A., Khotyaintsev, Y. V., Lindqvist, P.-A., Marklund, G., Wilder, F. D., Young, D. T., Torkar, K., Goldstein, J., Dorelli, J. C., Avanov, L. A., Oka, M., Baker, D. N., Jaynes, A. N., Goodrich, K. A., Cohen, I. J., Turner, D. L., Fennell, J. F., Blake, J. B., Clemmons, J., Goldman, M., Newman, D., Petrinec, S. M., Trattner, K. J., Lavraud, B., Reiff, P. H., Baumjohann, W., Magnes, W., Steller, M., Lewis, W., Saito, Y., Coffey, V., and Chandler, M.: Electron-scale measurements of magnetic reconnection in space, Science, 352, aaf2939,, 2016. a

Canuto, C.: Spectral methods in fluid dynamics, Springer series in computational physics, Springer-Verlag, ISBN 978-3-5405-2205-8, 1988. a

Channell, P. J.: Exact Vlasov–Maxwell equilibria with sheared magnetic fields, Phys. Fluids, 19, 1541–1545,, 1976. a

Chen, Y., Tóth, G., Cassak, P., Jia, X., Gombosi, T. I., Slavin, J. A., Markidis, S., Peng, I. B., Jordanova, V. K., and Henderson, M. G.: Global Three-Dimensional Simulation of Earth's Dayside Reconnection Using a Two-Way Coupled Magnetohydrodynamics With Embedded Particle-in-Cell Model: Initial Results, J. Geophys. Res.-Space, 122, 10318–10335,, 2017. a

Cowley, S. W. H. and Owen, C. J.: A simple illustrative model of open flux tube motion over the dayside magnetopause, Planet. Space Sc., 37, 1461–1475,, 1989. a

Dargent, J., Aunai, N., Lavraud, B., Toledo‐Redondo, S., Shay, M. A., Cassak, P. A., and Malakit, K.: Kinetic simulation of asymmetric magnetic reconnection with cold ions, J. Geophys. Res.-Space, 122, 5290–5306,, 2017. a

De Keyser, J. and Roth, M.: Equilibrium conditions and magnetic field rotation at the tangential discontinuity magnetopause, J. Geophys. Res.-Space, 103, 6653–6662,, 1998. a

Dorville, N., B. G., Aunai, N., D. J., and Rezeau, L.: Asymmetric kinetic equilibria: Generalization of the BAS model for rotating magnetic profile and non-zero electric field, Phys. Plasmas, 22, 092904,, 2015. a, b

Dungey, J. W.: Interplanetary Magnetic Field and the Auroral Zones, Phys. Rev. Lett., 6, 47–48,, 1961. a

Fadanelli, S., Faganello, M., Califano, F., Cerri, S. S., Pegoraro, F., and Lavraud, B.: North-South Asymmetric Kelvin-Helmholtz Instability and Induced Reconnection at the Earth's Magnetospheric Flanks, J. Geophys. Res.-Space, 123, 9340–9356,, 2018. a

Faganello, M. and Califano, F.: Magnetized Kelvin–Helmholtz instability: theory and simulations in the Earth's magnetosphere context, J. Plasma Phys., 83, 535830601,, 2017. a

Faganello, M., Califano, F., and Pegoraro, F.: Being on time in magnetic reconnection, New J. Phys., 11, 063008,, 2009. a

Frey, H. U., Phan, T. D., Fuselier, S. A., and Mende, S. B.: Continuous magnetic reconnection at Earth's magnetopause, Nature, 426, 533–537,, 2003. a

Furth, H. P., Killeen, J., and Rosenbluth, M. N.: Finite-Resistivity Instabilities of a Sheet Pinch, Phys. Fluids, 6, 459–484,, 1963. a

Hasegawa, H.: Structure and Dynamics of the Magnetopause and Its Boundary Layers, Monographs on Environment, Earth and Planets, 1, 71–119,, 2012. a

Labelle, J. and Treumann, R. A.: Plasma waves at the dayside magnetopause, Space Sci. Rev., 47, 175–202,, 1988. a

Lee, L. C. and Roederer, J. G.: Solar wind energy transfer through the magnetopause of an open magnetosphere, J. Geophys. Res.-Space, 87, 1439–1444,, 1982. a

Lele, S. K.: Compact finite difference schemes with spectral-like resolution, J. Comput. Phys., 103, 16–42,, 1992. a, b

Ma, Y.-J., Nagy, A. F., Toth, G., Cravens, T. E., Russell, C. T., Gombosi, T. I., Wahlund, J.-E., Crary, F. J., Coates, A. J., Bertucci, C. L., and Neubauer, F. M.: 3D global multi-species Hall-MHD simulation of the Cassini T9 flyby, Geophys. Res. Lett., 34, L24S10,, 2007. a

Manuzzo, R., Belmont, G., Rezeau, L., Califano, F., and Denton, R. E.: Crossing of Plasma Structures by Spacecraft: A Path Calculator, J. Geophys. Rese.-Space, 124, 10119–10140,, 2019. a, b

MMS Science Data Center: Index of /mms/sdc/public/data, available at:, last access: June 2019. a

Modolo, R., Chanteur, G. M., Dubinin, E., and Matthews, A. P.: Simulated solar wind plasma interaction with the Martian exosphere: influence of the solar EUV flux on the bow shock and the magnetic pile-up boundary, Ann. Geophys., 24, 3403–3410,, 2006. a

Nakamura, M. and Scholer, M.: Structure of the magnetopause reconnection layer and of flux transfer events: Ion kinetic effects, J. Geophys. Res.-Space, 105, 23179–23191,, 2000. a

Palmroth, M., Ganse, U., Pfau-Kempf, Y., Battarbee, M., Turc, L., Brito, T., Grandin, M., Hoilijoki, S., Sandroos, A., and von Alfthan, S.: Vlasov methods in space physics and astrophysics, Living Reviews in Computational Astrophysics, 4, 1,, 2018. a

Paschmann, G.: Observational Evidence for Transfer of Plasma Across the Magnetopause, Space Sci. Rev., 80, 217–234,, 1997. a

Phan, T. D., Escoubet, C. P., Rezeau, L., Treumann, R. A., Vaivads, A., Paschmann, G., Fuselier, S. A., Attié, D., Rogers, B., and Sonnerup, B. U. Ö.: Magnetopause Processes, Space Sci. Rev., 118, 367–424,, 2005. a

Retinò, A., Bavassano Cattaneo, M. B., Marcucci, M. F., Vaivads, A., André, M., Khotyaintsev, Y., Phan, T., Pallocchia, G., Rème, H., Möbius, E., Klecker, B., Carlson, C. W., McCarthy, M., Korth, A., Lundin, R., and Balogh, A.: Cluster multispacecraft observations at the high-latitude duskside magnetopause: implications for continuous and component magnetic reconnection, Ann. Geophys., 23, 461–473,, 2005. a

Roth, M., De Keyser, J., and Kuznetsova, M. M.: Vlasov theory of the equilibrium structure of tangential discontinuities in space plasmas, Space Sci. Rev., 76, 251–317,, 1996. a

Shi, Q. Q., Shen, C., Dunlop, M. W., Pu, Z. Y., Zong, Q.-G., Liu, Z. X., Lucek, E., and Balogh, A.: Motion of observed structures calculated from multi-point magnetic field measurements: Application to Cluster, Geophys. Res. Lett., 33, l08109,, 2006. a

Shumlak, U., Lilly, R., Reddell, N., Sousa, E., and Srinivasan, B.: Advanced physics calculations using a multi-fluid plasma model, Comput. Phys. Commun., 182, 1767–1770,, 2011. a

Tanaka, K. G., Retinò, A., Asano, Y., Fujimoto, M., Shinohara, I., Vaivads, A., Khotyaintsev, Y., André, M., Bavassano-Cattaneo, M. B., Buchert, S. C., and Owen, C. J.: Effects on magnetic reconnection of a density asymmetry across the current sheet, Ann. Geophys., 26, 2471–2483,, 2008. a

Whipple, E. C., Hill, J. R., and Nichols, J. D.: Magnetopause structure and the question of particle accessibility, J. Geophys. Res.-Space, 89, 1508–1516,, 1984. a

Publications Copernicus
Short summary
We investigate the magnetopause stability and mixing using a new three-fluid model aimed at reproducing the system configuration obtained directly from satellite data. This realistic model is a basic starting point for numerical simulations; however, the realistic three-fluid equilibrium presented in this paper should allow this work to be taken a step further and could be applied to other experimental cases in the future.
We investigate the magnetopause stability and mixing using a new three-fluid model aimed at...