Diffusion at the Earth magnetopause: enhancement by Kelvin-Helmholtz instability

. Using hybrid simulations, we examine how particles can diffuse across the Earth’s magnetopause because of ﬁnite Larmor radius effects. We focus on tangential discontinuities and consider a reversal of the magnetic ﬁeld that closely models the magnetopause under southward interplanetary magnetic ﬁeld. When the Larmor radius is on the order of the ﬁeld reversal thickness, we show that particles can cross the discontinuity. We also show that with a realistic initial shear ﬂow, a Kelvin-Helmholtz instability develops that increases the efﬁciency of the crossing process. We investigate the distribution functions of the transmitted ions and demonstrate that they are structured according to a D-shape. It accordingly appears that magnetic reconnection at the magnetopause is not the only process that leads to such speciﬁc distribution functions. A simple analytical model that describes the built-up of these functions is proposed.


Introduction
Since the open magnetosphere concept has been proposed by Dungey (1961), possible mechanisms to ensure mass loading of the magnetosphere from the shocked solar wind are highly debated. Magnetic reconnection plays a central role in this paradigm to explain the in situ observed coupling between magnetosphere and solar wind. The exact way magnetic reconnection develops is still controversial, the main reason being that large (MHD) scales have to be coupled with small (electron) scales, which is hardly achieved in numeri-Correspondence to: R. Smets (rsm@cetp.ipsl.fr) cal studies (a same kind of problem occurs in turbulence, see e.g. Frisch, 1995).
Magnetized plasmas consist of particles gyrating around convecting magnetic field lines. When magnetic reconnection occurs, plasma travels in the newly connected regions that were prohibited prior to reconnection. An alternative way to transfer matter across a magnetic boundary is diffusion. The topology of the magnetic field lines is preserved, but because of diffusion, particles can jump from one magnetic field line to another. Because of the collisionless nature of the plasma, the origin of diffusion cannot be classical viscosity, and many studies (see e.g. Winske et al., 1995;Treumann et al., 1995) have been dedicated to diffusion processes in collisionless plasma. On the whole it is commonly admitted that diffusion is less efficient than magnetic reconnection to ensure mass loading of the inner magnetosphere from the shocked solar wind (see Phan et al., 2005).
The physical process at work at the magnetopause is related to the nature of the magnetic discontinuity: in a tangential discontinuity there is no connection between upstream and downstream plasmas. In contrast, in a rotational discontinuity, a normal component of the magnetic field exists that connects both sides. Although the theoretical distinction is clear, experimental evidences are far more complicated because the measurement of a small component of the field is technically difficult. This difficulty is not due to the instruments, but to the uncertainty in determining the normal, the planarity of the discontinuity, ... Tests have been developed to check the Rankine-Hugoniot relations and demonstrated that the magnetopause sometimes is a rotational discontinuity (Paschmann et al., 1979). In other instances, the magnetopause is assumed to be a tangential discontinuity (Papamastorakis et al., 1984). Recent observations with Cluster  lead to the first statistical studies on the nature of the discontinuity with precise multipoint observations, showing that on the same day both situations can be observed.
Published by Copernicus GmbH on behalf of the European Geosciences Union.
Using hybrid simulations, we show that finite Larmor radius (FLR) effects are responsible for anomalous diffusion, due to the thickness of the magnetopause (MP), on the order of the ion Larmor radius (see e.g. Eastmann et al., 1996). Furthermore, at the Earth MP, the shear velocity between the solar wind flowing in the anti-solar direction and the quasi stagnant plasma of the inner magnetosphere can trigger the Kelvin-Helmholtz instability (KHI). The development of such KHI can amplify the diffusion process. The kinetic properties of the plasma are investigated, and compared to the classical structures observed during magnetic reconnection.
In the case of a neutral fluid, Chandrasekhar (1961) showed that KHI develops, regardless of the shear velocity magnitude. Talwar (1964) showed that when there is a component of the magnetic field tangential to the flow, the stabilizing effect of the magnetic tension results in a threshold for the instability development. The threshold as well as the growth rate of the instability is modified by the plasma compressibility (see, e.g. Miura and Pritchett, 1982;Pu and Kivelson, 1983).
Using compressible 2-D MHD code, Miura (1984) showed that up to 2% of the energy flux can cross the magnetic boundary. Miura (1987) also showed that the anomalous tangential stress associated with the instability ensures a momentum transport across the discontinuity that can reach a few percents. Still in the MHD interpretation framework, development of KHI can lead to plasma blobs (a direct consequence of mass transport, but not the only one) because of resistivity (Nykiri and Otto, 2001), electron inertial effects (Nakamura et al., 2004), or finite Larmor radius effects Huba (1996). KHI may also trigger turbulence (Matsumoto and Hoshino, 2004), and vice-versa, Shinohara et al. (2001) showed that kinetic instability (namely the low hybrid drift instability and the current sheet kink instability) can cascade in KHI.
Using 2-D hybrid simulations, Terasawa et al. (1992) studied the time evolution of a mixing index, providing insights on the total surface close to the instability where particles, initially from both sides of the discontinuity, are mixed. Thomas and Winske (1993) also used 2-D hybrid simulations to demonstrate the existence of small-scale structures (on the order of the ion Larmor radius) possibly related to such flux transfer events (see e.g. Southwood et al., 1988, for flux transfer observations). Fujimoto et al. (1994) extended this investigation and focused on mass transport, showing that mixing is quicker and larger than that due to finite Larmor radius overlapping at the interface. Thomas (1995) obtained similar results with 3-D hybrid simulations for the same geometry, but showed that a rotation of the magnetic field stabilizes the system, and leads to a less efficient mass transport.
Recently, clear observations of the KHI have been reported by Hasegawa et al. (2004) using the four spacecrafts of the Cluster 2 mission. Furthermore, this paper provided evi-dences of plasma penetration from the magnetosheath into the inner magnetosphere. Interestingly, this paper quotes: "... we could rule out local reconnection because... we did not find signatures of plasma acceleration due to magnetic stresses and D-shaped ion distribution characteristic of reconnection". Conversely, one may wonder whether Dshaped distribution functions provide unambiguous signatures of magnetic reconnection. We show in this paper that such distribution can actually be observed in the vicinity of a KHI.
In this paper, we will not study in details the development of the KHI, that is, the possible stabilization by Landau damping as demonstrated by Ganguli (1997), the vortex pairing or inverse cascad as put forward by Belmont and Chanteur (1989), or the coupling with smaller scales process (see e.g. Shinohara et al., 2001;Matsumoto and Hoshino, 2004). We use hybrid simulations of the KHI to study the mass transport across the magnetic boundary. Electrons being considered as a massless fluid, we do not resolve electron frequencies. Instabilities like Electron-Ion-Hybrid instability (see Ganguli, 1997;Ganguli, 1993, 1994) can then not be addressed. In Sect. 2, we describe the hypotheses inherent to hybrid simulations. In Sect. 3, we present the model of the initial tangential discontinuity. In Sect. 4, we discuss the development of the KHI. In Sect. 5, we examine the structure of the ion distribution function, and discuss the possible occurence of reconnection. In Sect. 6, we focus on particle dynamics to understand the origin of D-shaped ion distribution functions. In Sect. 7, we exhibit the microscopic process and the implications for spacecraft observations.

Simulation model
The calculations are carried out using a hybrid code (ions treated as particles and a massless electron fluid) with two spatial dimensions and three velocity dimensions. Details of the simulation algorithm are given by Winske and Quest (1986). Distances are normalized to the ion inertial length (c/ω P i ), time is normalized to the inverse of the ion gyrofrequency (ω −1 Ci ), mass is normalized to the ion mass, and density and magnetic field are normalized to their asymptotic values (the same on each side of the discontinuity), n 0 and B 0 , respectively. The size of the simulation box is X m in the x-direction and Y m in the y-direction. We use n x =256 grid points in the x-direction, and n y =128 grid points in the y-direction, X m =80 and Y m =40 (size of the box in xand y-direction), 100 macroparticles per cell and a time step t=0.005.
In the calculation of the electric and magnetic field, we neglected the displacement curent in the Maxwell-Faraday equation. Such an assumption prevents the development of high frequency modes. Ohm's law is then needed to define the electric field. Taking into acccount the Hall term, the electron pressure term and a resistivity, we have where V i is the ion fluid velocity, J the current density, P e the scalar electron pressure, η the resistivity, and n e =n i is the electron (ion) density (noted n hereinafter). We also need a closure equation for the electron fluid. As the development of the KHI is a low frequency process, we used an isothermal equation P e =n e k B T e with a constant value of the electron temperature T e , k B being the Boltzman constant.
We used a small resistivity in these calculations, mainly because of its stabilizing effect on the results, and chose η=0.0001. With this value, the diffusion time is τ D =40 000. The Alfvén time τ A being 2, the growth time of the resistive tearing mode is τ r = √ 2τ D τ A =400. We will see in the next section that the simulations are all performed over a total time T =400. Furthermore, KHI starts its linear growth phase at t∼100, and the penetration process starts at the very begining of the simulations. We can thus suspect that there is no magnetic reconnection due to anomalous resistivity that leads to the observed structures. This point will be detailed in Sect. 5 where we discuss the Walèn test.
The last classical hypothesis concern the boundary conditions. Because of the nature of the initial topology, we chose periodic conditions in the x-direction. This assumption cannot be made in the y-direction, because of the reversal of the magnetic field direction, and the bulk flow velocity. At Y =0 and Y =Y m , we assumed that particles are specularly reflected, and E X =E Z =d Y E Y =0.

Initial topology
We focused on tangential discontinuities with a 180 • rotation angle for the magnetic field. There is thus no normal component of velocity and magnetic field, and the only constraints for the tangential components are to satisfy the pressure balance. Such a geometry corresponds to the flanks of the MP when the interplanetary magnetic field (IMF) is southward. By definition, the rotation of the magnetic field is in a plane tangential to the discontinuity which we refer to as the XZ plane. As illustrated in Fig. 1, the magnetic field is in the +Z direction in side 1 of the discontinuity, and in the −Z direction in side 2. At the discontinuity, the rotation is such that the magnetic field is in the +X direction. The velocity profile is also represented, varying from −V 0 in side 1 to +V 0 in side 2.
We use a fluid equilibrium and consider the following profile P 0 (y) = P 1 + (P 2 − P 1 )χ (2) (with =(y − Y 0 )/λ, Y 0 =Y m /2, and χ=(1+ tanh )/2) for the magnetic field magnitude B 0 (y), the density n 0 (y), and the β i (y) and β e (y) parameter (β being the ratio between kinetic to magnetic pressure for ions and electrons, respectively). The tangential magnetic field is analytically described by with α=π(1−χ ). We first perform a reference run where the magnetic field magnitude, density, β i and β e parameters are unifom. We and β e =1.0. One remarks there is a gradient on the magnetic field direction, but not on its magnitude. As discussed later, ions are hence also magnetized in the field reversal region. Except when indicated, these values will also be used in the subsequent runs. We did not consider any kinetic equilibrium because none of them are satisfactory in our case. Some provided an analytical solution for tangential discontinuities with no velocity shear and no magnetic shear (Channell, 1976;Mottez, 2003). Lee and Kan (1979) provided a semi-analytical solution of tangential discontinuities including a velocity shear and a small magnetic shear. Roth and DeKeyser (1996) described a large (though incomplete) class of solutions given as convolution between exponential and error functions. The common problem of these equilibriums is that the proposed solution may not be unique. The numerical equilibrium itself requires some time (even short) to adjust, either because of numerical uncertainties (inherent to semi-analytical solutions), or because the given equilibrium is not the most likely. This adjustment time is inherent to the problem and cannot be avoided. As a result, we used a fluid equilibrium and monitored the ion distribution functions to ensure that the equilibrium is reached before physical processes of interest occur.
To monitor the time evolution of the distribution function at the beginning of the simulation, we rebuilt the 10 independent components of the heat flux Q. We focus on the heat flux Q because fine structures play an increasing role in the www.ann-geophys.net/25/271/2007/ Ann. Geophys., 25, 271-282, 2007 high order moments of the distribution functions. The integrale to build Q is made on the cells located at Y =Y 0 , corresponding to an integral conducted along the x-direction, at the center of the discontinuity. Being at the center of the discontinuity, the distribution function strongly depends upon the kinetic equilibrium. Figure 2 depicts the time evolution of the XXY component of the heat flux Q. This component of the heat flux is the only one that oscillates at the begining of the simulation, the others being approximately constant. The Q XXY value oscillates with a period about 6, that is close to the gyroperiod (2π). After about 3 gyroperiods, the value of Q XXY is close to zero. Though not represented here, this level is constant until T =400. Calculating the Q XXY value at different Y positions, the level of these oscillations proves to be always lower, and reaches the noise level with the same caracteristic time. This suggests that a kinetic equilibrium is reached at t∼20. As we will see hereinafter, the linear growth phase of the KHI starts at t∼100. Hence, there is no overlapping between the phase during which the equilibrium is reached, and the begining of the KHI development.
Another requirement of this study is to identify the position of the boundary between side 1 and side 2, at whatever stage of the simulation. To do so, Nykiri and Otto (2001) calculated the z-component of the magnetic vector potential. In the present case, because of the rotation of the magnetic field, a simple definition can be used: we define the MP location as the location where the z-component of the magnetic field changes sign. Figure 3 depicts the z-component of the magnetic field at (a) t=120 (b) t=200 and (c) t=260. Fig-Fig. 3. Color-coded z-component of the magnetic field at (a) t=120, (b) t=200, and (c) t=260 in the x-y-plane. ure 3a is the begining of the growing phase of the KHI, where the interface starts to ripple. Figure 3b shows the end of the linear growth phase of the KHI. In this figure, two wavelengths are noticeable at the interface. During the growing Ann. Geophys., 25,[271][272][273][274][275][276][277][278][279][280][281][282]2007 www.ann-geophys.net/25/271/2007/ phase, the wave fronts steepen until saturation. Thereafter, the two waves merge yielding new waves. The number of waves is not very clear in Fig. 3c. Fourier analysis presented later in the paper will help to picture the time evolution. After t∼260, a new growing phase start, with the build up of 2 new waves fronts that grow until a new saturation. As discussed at the end of Sect. 1, inverse cascad occurence colsely depends on the size of the simulation box.

Development of the Kelvin-Helmholtz instability
We start our analysis with the reference run described above.
As mentionned in Sect. 3, there is an initial velocity shear and reversal of the magnetic field, with no gradient in density, temperature and magnetic field magnitude. Prior to studying the ion distribution functions, we quantify the effect of the KHI on the mass transport across the discontinuity. We use the δZ value proposed by Thomas and Winske (1993), and defined by where X m is the system length in the x-direction, n 0 the initial ion density on each side of the boundary, and N T the total number of ions crossing the interface (from side 2 to side 1) defined in Sect. 3. The crossing process being symetric, the N T value computed considering particles crossing from side 1 to side 2 is statistically the same. Even if δZ is homogeneous to a distance, it has to be seen as the normalized number of particles crossing the MP from side 2 to side 1. The time evolution of this quantity is presented in Fig. 4. It should be noted here that at the begining of the simulation, particles close to the boundary may experience repeated excursions across it. To remove such behaviors from δZ calculated, the value of N T (and thus δZ) is computed only with particles which remain at least one gyroperiod on the initial side of the discontinuity prior crossing. As a result, one has δZ=0 until t=2π. In Fig. 4, δZ is strictly growing with time. Between t∼0 and t∼160, one observes a linear phase during which δZ evolves linearly with time. This demonstrates that the crossing process starts at the very begining of the simulation, when the KHI is not developed yet. Such a behavior was already obtained by Thomas and Winske (1993), and suggests that the crossing process is not directly triggered by the KHI. Between t∼160 and t∼220, a second linear phase is noticeable but with a higher slope value. This interval corresponds to the non-linear stage of the growing phase of the instability. Even though the crossing process is not directly triggered by the KHI, KHI clearly has an influence on this process and enhances it. Then, we observe a series of phases during which the slope of δZ increases and decreases.
To compare the time evolution of δZ with the stage of the KHI, we calculated the Fourier transform of the spatial series (in the x-direction) of the y-component of the magnetic field, averaged in the y-direction. In Fig. 5, we show the logarithm of this value (color-coded) as a function of time and K X mode. We also calculated the time derivative of δZ. The obtained value being quite noisy (as usual with numerical derivative calculation), we filtered it with a gaussian window of 8 points width. The obtained result is multiplied by a normalization coefficient (to fit the K X scale) and displayed with solid thick line in Fig. 5. At t∼100, the KHI develops and grows until t∼200. This can be seen in, Fig. 5 by the broadening of the B Y spectra. At t∼200, the instability reaches its saturation level after what, vortex coalescence is observed until t∼250. This phase clearly corresponds to a lowering of the spectral width of B Y . We then observe until t=400 two others saturation stages, preceded by a growing phase, and followed by a coalescence phase. One can note a good agreement between the time derivative of δZ and the spectral width of B Y .
For the sake of clarity, we introduce the following denomination for the particle orbits: we call "short crossings" particles that cross several times the discontinuity, with a time of flight between 2 crossings smaller than 2π, and "long crossings" particles that cross the discontinuity several times, with a time of flight between two crossings larger than 2π. We must also make a distinction between even number of crossings and odd number of crossing: odd values are associated to particles that change sides, whereas even values are associated to particles that end up on the same side of the discontinuity as they were initially.

Ion distribution functions
In this section, we build the ion distribution functions associated with particles that cross the discontinuity. Building ion distribution functions requires to define a protocol. Indeed, with particle-like simulations, distribution function results from an average on a statistically representative population. In the same way as in Fig. 4, distribution functions correspond to particles that cross the discontinuity, with an initial time of flight prior to the first crossing at least equal to 2π. Figure 6 displays the number of crossing particles colorcoded in each cell (of size X by Y ), depending on their X and Y position at t=240. The solid line indicates the position of the MP identified by the sign change of B Z . Figure 6 does not take into account "short crossing", but only particles that cross once the discontinuity. Hence, particles on sides 2 were initially on side 1 and vice-versa. One can clearly see that the crossing particles form islands located between the edges of the surface instability. Crossing particles being located in small regions, we build the distribution function associated with these 4 islands (2 on each side). Figure 7 depicts the ion distribution functions at t=240 associated with particles that experience 2 crossings. Since we have very similar distribution functions in these 4 islands (whatever the side we consider), the one presented in Fig. 7 shows an average over all islands. The left and right panel of Fig. 7 show the color coded phase space density in the plane perpendicular to the local magnetic field and in the parallel -perpendicular velocity direction, respectively. The left panel of Fig. 7 reveals that the distribution is gyrotropic. This result is not unexpected: because of the constraints on the time of flight, the distribution functions correspond to particles that experience an adiabatic motion. The hypothesis of random gyrophases is thus satisfied. Furthermore, the computed value of the temperature is approximately the same as the one used for the initialization. No evidence of heating or cooling is noticeable.
In the right panel of Fig. 7, a D-shaped distribution can be seen, most particles having a positive value of their parallel velocity. One should note that the sign of the bulk flow associated to the D-shaped distribution (positive in Fig. 7) does not depend on the side of the MP where it is observed. As mentioned earlier, the particles considered here fly more than 2π after crossing the discontinuity. If this time of flight is increased the same overall pattern is observed, with just a slightly smaller number of particles. As discussed in the introduction, the build-up of such D-shaped distributions is often associated with magnetic reconnection when such Dshaped distribution is observed (see e.g. Cowley and Shull, 1983). We already saw that the numerical resistivity is supposedly too small to let resistive tearing occur. To explore this distribution further and rule out the occurence of magnetic reconnections, we perform Walèn test (see e.g. Sonnerup et al., 1981).
This test consists in comparing the fluid velocity to the Alfvén velocity (for both x-and y-components). To do so, we build at t=240, the components of the Alfvén and fluid velocity from the magnetic and fluid structures we computed along a virtual crossing, at X=25, and between Y =20.3 and Y =26.5. This arises from the position of the MP, located (at X=25) at Y =23.4 and corresponds to 20 cells (10 in each   Figure 8 demonstrates that the Walèn test is not satisfied because there is no way to interpolate these points with a straight line having a slope equal to one. Similar calculations changing the localization and/or the stage of the KHI lead to the same conclusion. Accordingly, there is no magnetic reconnection in the present set of simulations. Starting from a tangential discontinuity, regardless of the stage of the KHI, Table 1. Parameters used for different runs. λ is the initial half thickness of the discontinuity, V 0 , B 0 and n 0 the initial asymptotic value of the fluid velocity, magnetic field magnitude and density, respectively, and β i the ionic β parameter. the discontinuity never becomes rotational. To understand the building of the D-shaped distributions, one must consider the charged particles dynamics.

Particle dynamics
With the D-shaped distribution obtained in Fig. 7, we calculate the associated mean parallel velocity noted V . A variety of runs were performed to characterize the conditions of occurence of such distributions, as well as the computed value V . They are summarized in Table 1. As indicated in Sects. 2 and 3, we have for run 1, λ=1.0, V 0 =0.5, B 0 =1.0, n 0 =1, β i =1.0, β e =1.0. We indicate for the other runs values that differs from these ones. This choice of parameters allows to explore differents values of the Alfvén velocity, of the shear velocity, of the thermal Larmor radius, and of the polarization of the discontinuity. For these 7 runs, we construct the ion distribution function as in Fig. 7. They are not depicted here, but in all cases, D-shaped structures are obtained. The efficiency of the penetration process can differ from one run to another by a factor of 2, which has (for our purpose) only consequences on the statistics to calculate V . We performed some other runs reversing the velocity shear direction and/or the direction of the asymptotic magnetic field, and obtained that it does not change the positive value of the mean parallel velocity in the D-shaped ion distribution function. The only way to achieve this is to reverse the sense of rotation of the magnetic field of the initial tangential discontinuity (being thus initially in the −X direction at the interface between side 1 and side 2). Other test runs showed that the magnitude of this mean parallel velocity only depends on the initial half thickness of the discontinuity λ and on the magnitude of the initial asymptotic magnetic field.
We computed the mean parallel velocity V and calculated the ζ value, defined as where C is the gyroperiod with the asymptotic magnetic field magnitude. This is the average distance traveled along the magnetic field by the particle (forming the core of the D-shaped ion distribution function) during a cyclotron turn. Figure 9 depicts the computed ζ value, depending on the initial λ value. It can be seen in this figure that a straight line with a slope close to one connects the various runs. Accordingly, the building of a D-shaped structure does not depend upon the magnitude of the Alfvén velocity or the shear velocity. We thus suspect the magnetic field rotation to play a role in this process. This idea is strenghtened by the fact that reversing the sense of rotation of the magnetic field changes the sign of the mean parallel velocity associated to distribution function. To examine this issue further, we examined test particle orbits. Figure 10 shows the orbit of a particle in the time varying electric and magnetic field during T =200. This particle starts initially at X∼4 and Y ∼22 on side 2. Following the shear velocity, this particle drifts in the +X direction. White triangles are drawn every T =20. Around t=160, the particle crosses the discontinuity (indicated by the black triangle), and ends up on side 1 of the discontinuity until the end of the simulation. This type of orbit is representative of particles that compose the bulk of the D-shaped structure. In other words, these particles never exhibits short crossing, and seldom long crossing. Crossing always occurs on time scale on the order of the gyroperiod, with nearly-adiabatic sequences prior to and after this crossing. One can understand that, in view of Fig. 10, the structure of the magnetic field, even if corrugated through time, has a structure in the z-direction similar to the field reversal initially imposed. Such behaviors were already reported in Smets (2000). They are not due to the curvature of the magnetic field lines (the initial magnetic field lines being straight lines as indicated in Fig. 1), but to the scale of the field reversal. In a like manner to the adiabatic parameter κ BZ defined by Büchner and Zelenyi (1989), Smets (2000) defined the κ RS parameter (RS standing for Reversal Sheet) as the square root of the ratio between the thickness of the reversal sheet and the maximum Larmor radius. This study is hence done at κ RS ∼1, and as predicted by Smets (2000), we never observe Speiser-like orbits and particle trapping at the interface, because of FLR effects. This type of non-linear dynamics has been reported by Smets (2002) using test-particle calculations in realistic magnetic topologies, obtained by MHD simulation of the KHI. But this work was done at κ RS ∼0.6 which, as demonstrated by Smets (2000), allows particle trapping.
As we already saw in Fig. 4, penetration process starts at the begining of the simulation, before the KHI develops. We can suspect that this process is not directly linked to the occurence of the KHI, but rather to the magnetic topology of the discontinuity. The KHI just appears as a speeding up factor for the penetration process.

Microscopic process
A heuristic interpretation of the penetration process can be obtained with a rough picture of the magnetic field reversal. As a matter of fact, details of the magnetic profile (hyperbolic tangent in the present case) is not necessary and knowledge of the asymptotic value of the magnetic field in each side and Fig. 11. Schematic view of the typical orbit followed by a particle crossing the magnetic field reversal.
at the center of the discontinuity is enough to drive the essential features. The associated process is depicted in Fig. 11. We considered side 1 and 2 with a magnetic field that is considered constant, in the +Z direction and −Z direction, respectively and we neglect the drift effect of the electric field: the associated fluid velocity being tangential to the discontinuity, it does not play a significant role in this process. At the interface, we consider also a constant magnetic field, in the +X direction that we call reversal sheet. Let us consider a particle, initially on side 2, with a positive parallel velocity. Starting from an initial position with positive Z value, this particle, because of its parallel velocity, goes toward the Z=0 plane. Furthermore, because of the magnetic field and its perpendicular velocity, the particle gyrates around a magnetic field line, with a constant Larmor radius. With the appropriate initial position (as chosen in Fig. 11), this particle skims the interface between side 2 and the reversal sheet. When at the interface, because of magnetic noise, the particle can step aside, in the reversal sheet. This point is marked A in Fig. 11, and for simplicity, chosen to be in the Z=0 plane. To reach side 1 at point marked B, the Larmor radius in the reversal sheet has to be equal to the half thickness of the discontinuity λ. At point A, the perpendicular velocity needed in the calculation of the Larmor radius is in the −Z direction, that is the parallel velocity in side 2. Because there is no electric field, the kinetic energy of the particle is constant, and the parallel velocity in side 1 is equal to the parallel velocity in side 2 (that turns to be a perpendicular velocity in the reversal sheet).
The above explanation is based on geometric arguments, and does not depend upon the details of the discontinuity, provided it is tangential. This suggests that the analytical profile that we chose for the magnetic field is not critical. On the other hand, it is clear from Fig. 5 that magnetic field associated to high values of K X plays an important role in this process. We computed the same figure as in Fig. 5    components of the electric field, and did not get any correlation between electric noise and enhancement of penetration process. It seems that magnetic noise is the essential source of penetration, even if ruling out the role of the electric noise would require us to perform some other runs with an appropriate low band filtering.
With the mechanism we just described, we explain easily why the crossing paticles are the one for which the computed ζ value is equal to λ. Further insights into this mechanism may be obtained from Fig. 12 that shows the initial particle velocity depending on its final value. Left panel of Fig. 12 is for parallel velocity, and right panel for perpendicular velocity. The first remark is that we have a clear symetry about the first bisector of the axis. In the left panel of Fig. 12, the largest number of particles crossing the interface is for both initial and final parallel velocity close to one, that is consistent with the mechanism discussed above. Furthermore, following the isocontour associated to the highest number of crossings, particles with a final parallel velocity larger than one have an initial parallel velocity smaller than one, and vice-versa. This is also consistent with the proposed mechanism. Going back to Fig. 11, let us consider a particle reaching point A with a parallel velocity giving a ζ value larger than λ. The Larmor radius in the reversal sheet being larger than λ, the particle will reach side 1 before achieving a complete half cyclotron turn. Thus, only a fraction of the perpendicular velocity in the reversal sheet will end-up in parallel velocity in side 1. This is the reason why crossing particles with large initial parallel velocity have a smaller one after crossing the discontinuity. This process being totally symetric, we explain in the very same way why particles with a small initial parallel velocity will have a larger one after crossing the discontinuity. In these two cases, initial parallel (perpendicular) velocity can transform in perpendicular (parallel) velocity after crossing the discontinuity, and viceversa. This clearly explains the right pattern of Fig. 12.

Discussion
The process described here can be viewed as a filter for the crossing particles. Figure 13 shows a schematic view of the distribution function, depending on the parallel velocity. The solid line represents the initial distribution function, the dashed line represents the distribution function associated with the crossing particles. For simplicity, we neglect the mean value of the bulk velocity in the initial distribution function. The vertical line represents the parallel velocity for which we have ζ =λ. Particles with a parallel velocity larger than this value are able to cross the discontinuity. The others can hardly achieve it.
The resulting D-shaped distribution is only associated with crossing particles. At a given point, there is no way to experimentally distinguish between crossing particles and noncrossing particles. This means that the observed distribution results from the superposition of a drifting maxwellian associated to the non-crossing particles, and a D-shaped distribution associated to crossing particles. Construction of realistic distribution is beyond the scope of this paper, but calculating the mixing index as suggested by Terasawa et al. (1992), we found that there exists regions where the proportion of crossing/non-crossing particles is half-half (at t=400). We claim that, even if fuzzed, D-shaped distribution should be observable in the vicinity of the MP, when the IMF is southward.
The spatial extension of the region in which such D-shaped distribution are observable around a tangential discontinuity is questionable. From Fig. 6, it appears that it is around the field reversal on a layer thickness of about 20 (even if it is not filled only by crossing particles). This layer may actually be much thicker since we have seen that the penetration process depends on the thickness of the discontinuity (only its efficiency depends on the KHI phase), ζ being a linear function of λ. The small value of δZ is clearly due to the fact that δZ=0 at t=0. As long as the thickness of the discontinuity remains constant (and we have no evidence of its enhancement) we will have penetration process. A proper answer to this question would need an asymptotic study of this behavior, on a time much larger than a few hundreds of ionic cyclotronic periods to reach a stationary state.
There also exists a density gradient at the MP. This gradient can be located at a different location from the magnetic field gradient that defines the boundary layer. We do not consider such a density gradient, and will explain why it is not a problem for our results. We computed two runs with a density gradient at the interface, one with a density ratio equal to 4, and another with a ratio equal to 10 (which is a quite realistic value at the MP, see e.g. Paschmann et al., 1978). We observe the same process of discontinuity traversal, and a D-shaped distribution. The relation between λ and ζ is still linear, but with a different value of the slope (depending on the density gradient). The diamagnetic drift associated with such a density gradient should then be considered in the ex-planation. But the important point however is that, looking at in situ observations, the density gradient can occur at a different location from the gradient of the magnetic field direction: density gradient and field reversal can then be decoupled (see e.g. Paschmann et al., 1978). Observations also show a thinner MP with a southward IMF than with a northward one, and a density gradient (considered as the MP thickness) thicker than the field reversal layer. The consequences of these remarks are that calculations with a realistic density gradient should not affect the build-up of D-shaped distribution if the density gradient is not located at the same position as the field reversal. If both density gradients and field reversal occur at the same location, this would only change the relation between λ and ζ , but not the overall structure of a D-shaped distribution.

Conclusions
The hybrid simulations performed demonstrate that in a tangential discontinuity as is the case in the flanks of the magnetosphere during southward IMF conditions, the KHI does not trigger the transport of matter across the MP but rather enhances and/or speed up it. The transmitted particles exhibit D-shaped distribution, the average parallel speed being controlled by the half thickness of the discontinuity and the asymptotic magnetic field magnitude. We hence conclude that there is no one-to-one relationship between magnetic reconnection and D-shaped distribution functions as is commonly postulated.