Collisionless magnetic reconnection: analytical model and PIC simulation comparison

Magnetic reconnection is believed to be responsible for various explosive processes in the space plasma including magnetospheric substorms. The Hall effect is proved to play a key role in the reconnection process. An analytical model of steady-state magnetic reconnection in a collisionless incompressible plasma is developed using the electron Hall MHD approximation. It is shown that the initial complicated system of equations may split into a system of independent equations, and the solution of the problem is based on the Grad-Shafranov equation for the magnetic potential. The results of the analytical study are further compared with a two-dimensional particle-in-cell simulation of reconnection. It is shown that both methods demonstrate a close agreement in the electron current and the magnetic and electric field structures obtained. The spatial scales of the acceleration region in the simulation and the analytical study are of the same order. Such features like particles trajectories and the in-plane electric field structure appear essentially similar in both models.

However, another mechanism is developed using MHD with the Hall effect invoking (Sonnerup, 1979;Terasawa, 1983;Hassam, 1984), which causes a Petschek-like configuration due to the generation of dispersive waves (Mandt et al., 1994;Rogers et al., 2001) and does not require anomalous resistivity.It turns out that the contribution of the Hall effect appears close to the X-line, at length scales in order of the proton inertial length (skin depth), which is defined as l p =c/ω p , where c is the speed of light, ω p = 4π ne 2 /m p is the proton plasma frequency, and m p is the proton mass.Inside this region, the Hall effect decouples the proton(ion) and electron motions, tears protons off from the magnetic field lines, while the magnetic field remains frozen into the electron fluid.
Electrons demagnetize much closer to the X-line, inside the so-called electron diffusion region (EDR), because of inertia or non-gyrotropic pressure (Vasyliunas, 1975).A scaling analysis estimates δ∼l e , where δ is the thickness of the EDR, l e =c/ω e is the electron inertial length (skin depth), and ω e = 4π ne 2 /m e is the electron plasma frequency.Studies within the frame of the Geospace Environment Modeling (GEM) Magnetic Reconnection Challenge project (Birn et al., 2001) have shown that the Hall reconnection rate is insensitive to the dissipation mechanism activated in the EDR.The rate of magnetic reconnection is nearly independent of V. Semenov et al.: Collisionless magnetic reconnection the strength of dissipation in hybrid simulations with resistivity (Mandt et al., 1994), two-fluid simulations (Ma and Bhattacharjee, 1996;Biskamp et al., 1997), and particle-incell (PIC) simulations (Shay and Drake, 1998;Hesse et al., 1999).In addition, the Hall reconnection rate is independent of the magnitude of the electron mass (Shay and Drake, 1998;Hesse et al., 1999;Pritchett, 2001;Ricci et al., 2002) and the system size (Shay et al., 1999;Huba and Rudakov, 2004).The dependence on the electron dissipation and the system size seems to be found in studies of forced reconnection and double tearing mode reconnection (Grasso et al., 1999;Wang et al., 2001;Porcelli et al., 2002;Bhattacharjee et al., 2005).Summarizing, the results obtained by all types of numerical simulations, the essential feature of fast reconnection is the presence of the Hall effect, which provides the rate of the steady-state reconnection to be approximately a constant of the order of 0.1.

Analytical study
Analytical studies of magnetic reconnection in the frame of Hall MHD (HMHD) meet difficulties that arise due to the EDR physics.However, as far as the reconnection rate does not depend on the mechanism of dissipation, one can restrict EDR contribution studies to the size of the EDR rather than its internal structure.Numerical simulations (Shay et al., 1998) confirmed the thickness of the EDR to be of the order of l e for anti-parallel Hall reconnection.The length of the EDR was thought to be ∼10 l e (Shay et al., 1999;Huba and Rudakov, 2004), but later results suggested that it expands to the edge of the proton dissipation region, i.e., ∼10 l p (Daughton et al., 2006).This contradiction seems to be resolved in recent PIC simulations which indicate that the EDR develops into a two-scale structure along the outflow direction.The length of the electron current layer is found to be sensitive to the proton/electron mass ratio, approaching 0.6 l p for a realistic electron mass value.In addition, an elongated outflow electron jet is formed in the outflow region, and its length extends to 10 s of the proton skin depth (Shay et al., 2007).
Analytical studies of the problem may be simplified by using the electron Hall MHD (EHMHD) approximation.In the nearest vicinity of the stagnation point, at a length scale of the order of l p , the proton velocities are small compared to the electron velocities.Hence, one may consider the electric current in this region as the electron current only (Biskamp, 2000), where V e is the electron bulk velocity.For a detailed EHMHD-analysis of the problem, see Uzdensky and Kulsrud (2006).These authors have shown that for quasistationarity and translational symmetry assumed, the magnetic field and electron velocity can be expressed in terms of just a single one-dimensional function.This function is a magnetic potential of the in-plane magnetic field (poloidal flux function).
In addition, the authors have found out that, neglecting the ion current, one gets the Grad-Shafranov equation for this potential.
In the work of Korovinskiy et al. (2008) an analytical model of self-consistent steady-state collisionless magnetic reconnection in an incompressible plasma was developed based on the Grad-Shafranov equation for the magnetic potential.We outline this model in next paragraphs and review the main results derived; in-depth study is provided in cited paper.
Firstly, we chose a coordinate system as follows: The Xaxis coincides with the magnetic field direction at infinity (in the upper semiplane), the Y-axis is directed along the X-line, and the Z-axis is perpendicular to both of them.We assume homogeneity in the Y direction, so all quantities are assumed to be independent of Y .We avoid the description of the EDR internal processes and consider only its size, which we suppose to be of the order of l e in its cross section (Z direction) and ηl p along it (X direction), where η is a coefficient of the order of 1. Outside the EDR the plasma is supposed to be nonresistive; furthermore, it is assumed to be quasi-neutral and incompressible.
The two-fluid description of our problem is determined by the following equations, ∇ × E = 0, (5) Here, Eq. ( 2) is the equation of the proton motion, where P p is the scalar proton gas pressure, and V p is the proton bulk velocity; Eq. ( 3) is the Ohm law, where P e is the scalar electron gas pressure; Eq. ( 4) is the Ampère law; Eq. ( 5) is the Faraday law; Eq. ( 6) is the Gauss law; and Eq. ( 7) is the mass conservation law for each particle species.
In our steady-state 2.5 D case, the electric field E y must be a constant, according to Faraday's law (Eq.5).So, we define where is the reconnection rate which we assume to be small, 1, and Here, B 0 is the magnetic field value above the X-line at the upper boundary of the examined region and V A is the corresponding proton Alfvén velocity.
To resolve the system (2-7), we introduce dimensionless quantities: The magnetic field strength B=B/B 0 , the proton and electron bulk velocities Ṽp,e =V p,e /V A , the electric field strength Ẽ=E/E A , the gas pressure Pp,e =P p,e /P 0 , and the length scales r=r/ l p .Here P 0 =B 2 0 /4π, and r=(x, y, z).Secondly, we introduce an electric potential ˜ via Ẽ=− ∇ ˜ .Omitting the tildes, we rewrite Eqs.(2-7) for normalized quantities, bearing in mind the EHMHD approximation (1), Note that has a linear dependence on Y coordinate, so that ∂ /∂y=− .Under this point of view, we can present as a sum of two terms, namely (x, y, z)=φ(x, z)− y.Using an effective potential φ eff ≡φ−P e , we eliminate quantity P e from the Ohm law (10).
We also introduce a magnetic potential A(x, z), where e y is the unit vector and ⊥ denotes the XZ plane.
At last, we note that accordingly to the Ampère law (Eq.11), the magnetic field B y is the stream function for the electron in-plane velocity (Biskamp, 2000), Bearing in mind the EHMHD approximation (1), we note that Eq. ( 10) looks very similar to the condition of magnetohydrostatic equilibrium, In 2.5 D models, expression (16) yields the famous Grad-Shafranov equation (Schindler, 2007), where ⊥ is the 2 D Laplace operator.Analogously, we obtain where G(A) is an unknown modelling function.The other equations of the system (9-13) take the following form V py (r) = r r 0 Here, Eq. ( 19) is the equation for the out-of-plane magnetic field B y , where k is the quadrant number and ds f l is an elementary displacement along the projection of the magnetic field line onto the XZ plane; Eq. ( 20) is the equation for the effective electric potential φ eff ; Eq. ( 21) is the Bernoulli equation for the in-plane motion of protons, where 2 is a total pressure and C tr is a constant along the trajectory; Eq. ( 22) is the continuity equation, where V p⊥ is the proton in-plane velocity; and Eq. ( 23) is the equation for the out-of-plane proton velocity V py , where ds tr is an elementary displacement along the projection of the proton trajectory onto the XZ plane.
Thus, the initial complicated system splits, and the solution of the problem bases on the solution of the Grad-Shafranov equation for the magnetic potential (18).A scaling of the problem allows us to make use of the boundary layer approximation ∂/∂x ∂/∂z.Under this approximation, the Laplace Eq. ( 18) has a following solution with the boundary condition B z (x, 0).The unknown function G(A) has a simple physical meaning, namely it is the main part of the electric potential, while its derivative is the out-of-plane electron velocity.This allows us to preset function G(A) manually in order to obtain the solution of the problem (Korovinskiy et al., 2008), or get it from some other source, e.g., from PIC simulations.
Note that Eqs.(18-23) do not contain any dissipation so solution obtained is, strictly speaking, nonapplicable inside the EDR.Indeed, accordingly to Eq. ( 19) the magnetic field B y tends to infinity at the origin, where the in-plane magnetic field goes to zero.Therefore, we must interrupt the calculation of B y at the EDR boundary and have to consider this region separately.Advantageously, the EDR is very thin and comparatively short, the function B y (x, z) is smooth, and B y (0, 0)=0 due to the symmetry condition.Therefore, it is justified to neglect EDR contribution in B y .
Equation ( 19) claims also that the extreme values of B y are located at the separatrices of the in-plane magnetic field, as well as extreme values of |∇G(A)| and |V ey | are.Using the condition for solvability of Eq. ( 24), we obtain estimation of the extremum electric field, max |E|∼10E A .Note that outside the EDR, the contribution of ∇P e is negligible small, so φ≈φ eff .As for the electron velocity V ey , the Ampère law yields where δ is the EDR width measured in electron skin depths l e and V Ae = m p /m e V A is the electron Alfvén velocity.At last, the proton in-plane motion obeys the Bernoulli Eq. ( 21).Under the simple assumption =const, the drop of the electric potential accelerates protons up to V A .More realistically, the total pressure is fixed by its distribution at the upper boundary of the EHMHD domain, so that = (x).Thus, the system of Eqs.(18-25) expresses a self-consistent solution of our problem based on the modelling function G(A) with the boundary conditions B z (x, 0) and (x, z max ).To check the effectiveness of the model we take these functions from the PIC simulation of reconnection and then compare our analytical solution and the numerical model.

PIC simulation
The explicit particle-in-cell code P3D (Zeiler et al., 2002) is used for a simulation of 2.5 D reconnection.In brief, the P3D is an electromagnetic full particle code; the Boris algorithm (Birdsall and Langdon, 1991) is used for the numerical solution of the equation of motion.The electromagnetic field solver uses leapfrog scheme to advance fields in time.For the initial condition, we take conventional a Harris neutral current sheet (Harris, 1962), where with a background plasma density n b =0.2 and a half-width of the initial current sheet λ=0.4 l p .
The magnetic field is normalized to its maximum value in the lobes and the density is normalized to its current sheet maximum.A moderate initial GEM-type perturbation (Birn et al., 2001) is added to ignite reconnection where L x =L z =38.4 l p are the sizes of the computational box and the intensity of perturbation is 0 =0.3.A quasistationary state is achieved at t=15 −1 p , where −1 p is the inverse proton gyrofrequency (see Fig. 7), and the simulation parameters at t=20 are further taken as a reference to be compared with the analytical study.The mass ratio is m p /m e =64 and the temperature ratio is T p /T e =3/2.
Open boundary conditions for the fields are implemented at the exhaust boundaries to allow a free outflow of the plasma (Divin et al., 2007;Pritchett, 2001).
A perfect electric conductor (PEC) boundary closes the simulation box at z=±19.2.Under the boundary conditions adopted, not more than 15% of magnetic flux and particles escape through the outflow boundary by t=20.In the following section, the results of our simulations are presented.

Comparison of the results
A plot of the electric current j ey at x=0 is presented in Fig. 1 at the left, where the EDR is a well-recognizable region where the electron current dominates over the proton one.The EDR half-width δ comes up to 3/4 l p , i.e., δ≈6 l e under the used mass ratio.According to our estimation (26), the analytical model gives max |V ey |≈7 V A , and the value obtained in the PIC simulation is 6 V A (see Fig. 1,right).The electron/proton current ratio is j e /j p ≈11 in the origin and it decreases moving away from the X-line.As far as the EHMHD assumption is j e j p , we restricted the region of the analytical model by the value x=4, where j e /j p ≈5. Analogously, the upper boundary of the modelling region is restricted by the value z max =30 l e .This value corresponds to z max =4 in the PIC simulation (m p /m e =64) and z max =0.7 in the analytical modelling (m p /m e ≈1840).At last, the EDR half-length reaches 2 l p .
The PIC simulation provides us dG/dA≡V ey (A) and B z (x, 0) (see Fig. 2).
As for the total pressure (x, z max ), it turns out to be a linearly increasing but weakly varying quantity, (0, z max )=0.61 and (4, z max )=0.66.The last parameter of the analytical model is the reconnection rate ≡E y .Its value obtained in the simulation is 0.2.
The electron trajectories obtained from the analytical study and the PIC simulation are compared in Fig. 3.The magnetic field separatrix mapped by the electric current is clearly visible in both cases.In fact, this picture shows a classical Hall current structure (Sonnerup, 1979), observed in the magnetosphere (e.g.Alexeev et al., 2005) and in laboratory experiments (e.g.Cothran et al., 2005).The electron jet in X direction is visible as well, in agreement with results of other authors (Daughton et al., 2006;Shay et al., 2007).The dependence of the velocity of this jet of X is plotted in Fig. 4. The analytical model underestimates the electron velocity (approximately 50%) as compared to that of the PIC simulation.
The proton velocities demonstrate a better agreement, with acceleration up to 1.5±0.1 V A in both models.The proton trajectories and the magnetic field structure are presented in Fig. 5.At last, the electric fields E z shown in Fig. 6, are in good qualitative agreement as well.The localization of the extremum of E z corresponds to the jump of the electric potential across the separatrices as predicted.

Conclusions
One can see that the plasma characteristics obtained in both models are qualitatively equal.As for numeric values, the analytical model is not precise everywhere.While some values predicted are quite accurate (e.  others differ noticeably (e.g., V ex , E z ).We attribute this difference to the simplifications adopted.Namely, while outside the EDR P e is a weakly varying quantity, indeed, the situation is completely different inside the EDR.The electron pressure turns there to be an anisotropic tensor and ∇ Pe becomes the dominant term in the Ohm law.This term is responsible for the freezing-out of electrons in a thin and very stretched (∼10 l p ) layer mapping the X axis, called external EDR, where electron jets develop (see Fig. 6).Though this effect is completely out of the scope of our analytical study.Nevertheless, the analytical solution obtained demonstrates all essential Hall reconnection features and a close qualitative agreement with results of the PIC simulation.
Besides, this solution claims that a powerful mechanism of electron acceleration in the X-line direction is required.Accordingly to the estimation (26), it must accelerate electrons up to the electron Alfvén velocity inside the EDR and on the separatrices.At the downstream edge of the EDR, these accelerated electrons are deflected by the Lorentz force in X-direction and then get decelerated in the outflow region, pulling protons there.

Fig. 2 .
Fig. 2. PIC simulation data: Function dG/dA≡V ey (A) (on the left); magnetic field B z (x, 0) (on the right).The point A=0 on the left picture corresponds to the magnetic field separatrices.Simulation data are shown by the thin line, approximation is shown by the thick curve.

Fig. 3 .
Fig. 3. Electron trajectories in XZ plane: analytical model (on the left) and PIC simulation (on the right).

XFig. 6 .
Fig. 6.Distribution of electric field E z in XZ plane: analytical model (top) and PIC simulation (bottom).

Fig. 7 .
Fig. 7. PIC simulation: reconnected magnetic flux and reconnection electric field E R =E y histories.