Global three-dimensional simulation of Earth's dayside reconnection using a two-way coupled magnetohydrodynamics with embedded particle-in-cell model: initial results

TL;DR

This study uses a novel 3D global simulation combining magnetohydrodynamics and particle-in-cell methods to investigate Earth's dayside magnetic reconnection and flux transfer events, capturing both large-scale and kinetic phenomena.

Contribution

It introduces a coupled MHD-EPIC model for global Earth's magnetosphere simulation, revealing detailed FTE dynamics and kinetic features consistent with observations.

Findings

FTEs form quasi-periodically near the subsolar point and evolve into flux ropes.

Kinetic features such as crescent electron distributions are captured.

Lower hybrid drift instability observed with properties matching MMS data.

Abstract

We perform a three-dimensional (3D) global simulation of Earth's magnetosphere with kinetic reconnection physics to study the flux transfer events (FTEs) and dayside magnetic reconnection with the recently developed magnetohydrodynamics with embedded particle-in-cell model (MHD-EPIC). During the one-hour long simulation, the FTEs are generated quasi-periodically near the subsolar point and move toward the poles. We find the magnetic field signature of FTEs at their early formation stage is similar to a `crater FTE', which is characterized by a magnetic field strength dip at the FTE center. After the FTE core field grows to a significant value, it becomes an FTE with typical flux rope structure. When an FTE moves across the cusp, reconnection between the FTE field lines and the cusp field lines can dissipate the FTE. The kinetic features are also captured by our model. A crescent electron phase space distribution is found near the reconnection site. A similar distribution is found for ions at the location where the Larmor electric field appears. The lower hybrid drift instability (LHDI) along the current sheet direction also arises at the interface of magnetosheath and magnetosphere plasma. The LHDI electric field is about 8 mV/m and its dominant wavelength relative to the electron gyroradius agrees reasonably with MMS observations.

Figures

Equations

$$\mathbf{E}=-\mathbf{u}\times \mathbf{B}+\frac{\mathbf{J}\times \mathbf{B}}{q_e{n}_e}-\frac{\nabla p_e}{q_e{n}_e}$$

$$\frac{\partial p_e}{\partial t}+\nabla \cdot (p_e{\mathbf{u}}_e)=(\gamma -1)(-p_e\nabla \cdot {\mathbf{u}}_e)$$

$$\Delta t=c_0\cdot \min (\Delta x/v_{rms},\Delta y/v_{rms},\Delta z/v_{rms})$$

$$\mathbf{E}=-{\mathbf{u}}_{\mathbf{i}}\times \mathbf{B}+\frac{1}{q_i{n}_i}\mathbf{j}\times \mathbf{B}-\frac{1}{q_i{n}_i}\nabla p_e=-{\mathbf{u}}_{\mathbf{e}}\times \mathbf{B}-\frac{1}{q_i{n}_i}\nabla p_e$$

$$E\sim \frac{k_B{T}_i}{q_i{r}_i}$$

Full text

\journalname

JGR-Space Physics

Center for Space Environment Modeling, University of Michigan, Ann Arbor, MI 48109, USA

Department of Physics and Astronomy, West Virginia University, Morgantown, WV 26506, USA.

KTH, Stockholm, Sweden.

Los Alamos National Laboratory, Los Alamos, NM 87545, USA

\correspondingauthor

Yuxi Chenyuxichen@umich.edu

{keypoints}

A one-hour long global simulation of Earth’s magnetosphere with kinetic modeling of the dayside reconnection

Crater FTE is found at the early stage of a flux rope formation

Kinetic phenomena are found from the global simulation

1 Introduction

Magnetic reconnection between the interplanetary magnetic field (IMF) and Earth’s intrinsic dipole field is regarded as the most important mechanism for mass and energy transfer from the solar wind to the magnetosphere. Flux transfer events (FTEs) are widely considered as a phenomenon related to dayside non-steady reconnection (Russell and Elphic, 1978). An FTE is a bundle of reconnected magnetic fluxtubes created at the magnetopause and moving anti-sunward along the magnetopause. Such events are characterized by a bipolar variation of the magnetopause normal magnetic field $B_{N}$, and are usually associated with an enhancement of core field, the magnetic field component along the axial direction of the FTE. An FTE exhibits a flux-rope structure in three-dimensional space. It has been observed that the plasma inside an FTE is usually a mixture of magnetospheric and magnetosheath plasma (Daly et al., 1981), indicating that FTEs are generated by magnetic reconnection process. The diameter of an FTE can vary from several ion inertial lengths (Eastwood et al., 2016) (a few hundred kilometers) to several Earth radii (Rijnbeek et al., 1984; Hasegawa et al., 2006). In the dawn-dusk direction along the magnetopause, FTEs can extend over a long distance (Fear et al., 2008). FTEs frequently occur as a quasi-periodic process, and Rijnbeek et al. (1984) reported that the FTEs were observed about every 8 minutes during periods of southward magnetosheath magnetic field.

FTEs have been studied with various global numerical models. Compared to local simulations, a global model can offer more realistic plasma and magnetic field context. Fedder et al. (2002) used a global ideal MHD model to study the generation of FTEs. The typical magnetic field signature is captured by their model, and their simulation suggests that the FTEs are formed by non-steady reconnection along the separator at the magnetopause. Raeder (2006) performed a high resolution ideal MHD simulation with the OpenGGCM model. FTEs formed by multiple X line reconnection (Lee and Fu, 1985) with a tilted dipole field in this study. Dorelli and Bhattacharjee (2009) revisited the FTE generation mechanism with resistive MHD using the OpenGGCM model, and the authors argue that the FTEs are generated by flow vortices and the formation of new X lines is the consequence, rather than the cause of FTE formation. Sibeck et al. (2008) studied crater FTEs with the BATS-R-US MHD model. All these global simulations are based on ideal or resistive MHD codes, and the generation of FTEs relies either on ad hoc resistivity (Dorelli and Bhattacharjee, 2009) or numerical resistivity (Fedder et al., 2002; Raeder, 2006). Recently, a 2D-3V global magnetospheric hybrid-Vlasov simulation was performed to study magnetopause reconnection and FTEs by Hoilijoki et al. (2017).

Typical FTEs are associated with an enhancement of the field strength at the center of a flux rope. On the other hand, the so-called crater FTEs show more complicated structure: the center field is surrounded by two ‘trenches’ and the field strength usually show a dip just at the center (LaBelle et al., 1987; Owen et al., 2008). Typical FTEs are more frequently observed than crater FTEs (Zhang et al., 2010). The generation mechanism of crater FTEs has been explored with both numerical simulations (Sibeck et al., 2008) and analytic models (Zhang et al., 2010). Zhang et al. (2010) proposed that crater FTEs are the initial stage of typical FTEs based on hundreds of events selected from THEMIS observations. The structure of the core field can be even more complicated, for example, Eriksson et al. (2016) found a tripolar core field flux rope at the magnetopause.

It is widely accepted that the formation of FTEs is related to the dayside magnetopause reconnection, which is a kinetic process for collisionless plasma. Therefore it is important to include proper kinetic effects into the numerical model in order to produce FTEs in a physical way. The MHD with embedded PIC (MHD-EPIC) model developed by Daldorff et al. (2014) makes it feasible to use a kinetic model to study reconnection and FTEs with realistic magnetospheric configuration for the first time. Because of the small kinetic scales inside the magnetosheath, for example, the ion inertial length $d_{i}$ is about 60km $\sim 1/100\,R_{E}$, we have to artificially increase the kinetic scales in the present 3D global simulation. As shown by our companion paper (Tóth et al., 2017, submitted companion paper), this scaling has on significant effect on the large scale structures, while the kinetic phenomena occur at linearly increased scale. Since the kinetic scale physics is included in our global model, the reconnection related kinetic phenomena, like the crescent shape electron phase space distribution, the Larmor electric field and the lower hybrid drift instability (LHDI), are all captured by the model. The crescent distribution was first found by Hesse et al. (2014) from 2D local simulation, then observed by the Magnetospheric Multiscale (MMS) mission recently (Burch et al., 2016). It is formed by the magnetosheath electrons reaching the stagnation point and accelerated by the Hall electric field (Bessho et al., 2016; Shay et al., 2016). This special distribution has been proposed as an indicator of the magnetic reconnection location (Hesse et al., 2014). The Larmor electric field is potentially another signature that can help to identify the location of reconnection site (Malakit et al., 2013). It is on the magnetosphere side, normal to the current sheet and pointing away from the X line. The lower hybrid drift instability (LHDI) develops along the current direction (Daughton, 2003; Roytershteyn et al., 2012), and it has been observed recently by MMS satellites (Graham et al., 2016). LHDI was considered as a potential source to create anomalous resistivity for reconnection (Huba et al., 1977), but previous research (Mozer et al., 2011) has suggested the related resistivity may be not large enough. However a recent 3D simulation showed LHDI may still play an important role near the diffusion region because of the presence of turbulence (Price et al., 2016) .

In the following sections we will describe the MHD-EPIC model, the simulation setup, and then discuss the simulation results.

2 Model description

The MHD-EPIC model has been successfully applied to investigate the interaction between the Jovian wind and Ganymede’s magnetosphere, where the ion inertial length is large compared to the size of its magnetosphere (Tóth et al., 2016). In this paper, the same model is applied to study Earth’s magnetosphere, which is more challenging because of the small kinetic scale. The MHD-EPIC model two-way couples the BATS-R-US (Powell et al., 1999; Tóth et al., 2008) MHD code and the implicit particle-in-cell code iPIC3D (Markidis et al., 2010) through the Space Weather Modeling Framework (SWMF) (Tóth et al., 2005, 2012). A general description of the these models and the simulation setup is provided in this session.

2.1 Global MHD model: BATS-R-US

In order to make the MHD model as complete as possible, the Hall term and the electron pressure gradient term are included in the generalized Ohm’s law, and a separate electron pressure equation is solved. The generalized Ohm’s law we use is:

[TABLE]

where $q_{e}$, $n_{e}$ and $p_{e}$ are the charge per electron, electron number density and electron pressure, respectively. The electron pressure is obtained from:

[TABLE]

where $\gamma=5/3$ is the adiabatic index, and $\mathbf{u_{e}}=\mathbf{u}-\mathbf{J}/(q_{e}n_{e})$ is the electron velocity.

From the numerical prospective, it is not trivial to incorporate the Hall term into the MHD equations. The Hall MHD equations support the whistler mode wave, which is dispersive and the characteristic speed is inversely proportional to the wave length. Since the shortest wave length that can be resolved in a numerical system is twice the cell size, the fastest whistler wave speed is proportional to $1/\Delta x$. For an explicit time integration scheme, the time step is limited by the CFL condition, which leads to a time step approximately proportional to $1/(\Delta x)^{2}$ for Hall MHD. In order to use a reasonably large time step, a semi-implicit time discretization is employed (Tóth et al., 2012). The semi-implicit scheme treats the stiff terms, which is the Hall term here, and other terms separately. Excluding the Hall term, the rest of the equations are updated with an explicit scheme, and the time step is only limited by the fast magnetosonic wave speed. The Hall term is handled by an implicit solver after the explicit update has been done.

The typical solar wind condition at 1AU with purely southward IMF is used as the boundary condition to drive the magnetosphere: $\mathbf{B}=(0,0,-5)$ nT, mass density $\rho=5\,\mathrm{amu/cm^{3}}$, ion pressure $p_{i}=3.45\times 10^{-3}$ nPa, and solar wind velocity $\mathbf{u}=(-400,0,0)$ km/s. Electron pressure $p_{e}=8p_{i}=2.76\times 10^{-2}$ nPa is used, so that after crossing the shock, where the ions are heated by converting bulk into thermal energy while the electron thermal energy changes adiabatically, the ion-electron pressure ratio is about $p_{i}/p_{e}\sim 2.5$. Wang et al. (2012) shows that the temperature ratio $T_{i}/T_{e}$ in the solar wind varies from $0.1\sim 2$, and the ratio is about $4\sim 12$ inside the magnetosheath. The $T_{i}/T_{e}$ ratio, which is the same as $p_{i}/p_{e}$, used in the simulation is close to but slightly smaller than the typical observed ratio. We use $T_{i}/T_{e}=1/8$, because our numerical experiments show that the electrons can be numerically heated by PIC if colder solar wind electrons are used as boundary condition. A magnetic dipole with 30116 nT field strength at the magnetic equator is used. Its magnetic axis is aligned with the z axis. The total magnetic field $\mathbf{B}$ is split into the intrinsic dipole field $\mathbf{B_{0}}$ and the deviation $\mathbf{B_{1}}$. A three-dimensional block-adaptive Cartesian grid is used to cover the whole magnetosphere: $-224\,R_{E}<x<32\,R_{E}$, $-128\,R_{E}<y<128\,R_{E}$ and $-128\,R_{E}<z<128\,R_{E}$. Since we focus on the dayside dynamics in this paper, the mesh along the dayside magnetopause is refined to high resolution with $\Delta x=1/16\,R_{E}$ (see Figure 1). 59 million cells are used in total. At the inner boundary $r=2.5\,R_{E}$, the density is fixed as $28\,\mathrm{amu/cm^{3}}$, the pressure and the magnetic field $\mathbf{B_{1}}$ have zero gradient, the radial velocity is zero, while the tangential velocity is calculated from the ionosphere electrodynamics model developed by Ridley et al. (2004).

2.2 Implicit particle-in-cell model: iPIC3D

The semi-implicit particle-in-cell code iPIC3D was developed by Markidis et al. (2010). The advantage of iPIC3D over explicit particle-in-cell codes is that iPIC3D is linearly unconditionally stable, so that iPIC3D can handle larger time step and larger cell size than explicit PIC. Compared to the explicit PIC method, the cell size of iPIC3D is chosen based on the scale of interest instead of the Debye length, and the time step of iPIC3D is not limited by the plasma frequency or the speed of light, but the accuracy condition, which requires ${v}_{rms}\Delta t/\Delta x<1$ on all grid nodes for all species, where ${v}_{rms}$ is the root mean square of macro-particle velocities. In order to make the simulation as efficient as possible while keeping the accuracy condition satisfied, we use an adaptive time step:

[TABLE]

calculated for each grid nodes and the minimum is taken over the whole PIC mesh. The root mean square velocity $v_{rms}$ is similar to the thermal velocity but contains the effect of bulk velocity. $c_{0}$ is a coefficient that should be smaller than 1. $c_{0}=0.4$ is used for the simulation in this paper.

Since the focus of this paper is the dayside magnetopause reconnection, the embedded PIC box is placed near the sub-solar magnetopause, where reconnection happens under purely southward IMF. In the GSM coordinates, the region inside $8\,R_{E}<x<12\,R_{E}$ and $-6\,R_{E}<y,z<6\,R_{E}$ (see Figure 1) is solved by iPIC3D. The PIC region covers the magnetopause and it is just inside the bow shock. Since the size of the ion diffusion region is the same order as the ion inertial length, such kinetic scale should be resolved in order to capture reconnection kinetic physics. However, the ion inertial length $d_{i}=c/\omega_{pi}$ is about 60km $\sim 1/100\,R_{E}$ for a typical magnetosheath density of $20\,\mathrm{amu/cm^{3}}$. The length is so small that it is extremely difficult to resolve even for a 3D global MHD model, not to mention the PIC code. Scaling up the kinetic length helps to reduce computational resources. In the present simulation, all the fluid values, including density, pressure, velocity, IMF and dipole field strength, hence the derived values like the sound speed, Alfven velocity and plasma beta, are realistic so that the global structure of the magnetosphere is comparable to the real situation. On the other hand, the ion inertial length is scaled up 16 times to about $1/6\,R_{E}$ in the magnetosheath by artificially increasing ion mass per charge by a scaling factor of 16. Since all the quantities are normalized in the numerical model, there are several ways to understand or interpret the scaling. One way is treating the scaling as changing the charge of ions and electrons. Compared with the original system, we reduce the charge by a factor of 16 while all the other basic physical quantities, like mass per ion, number density, and temperature remain realistic. From the perspective of ideal magnetohydrodynamics, the scaled system is exactly equivalent to the original one. For a particle-in-cell code, the reduction of charge per ion reduces the electromagnetic force on an ion and therefore increases the gyroradius and gyroperiod by a factor of 16. But the gyroradius and the gyroperiod are still several orders smaller than the global spatial and temporal scale, for example the distance from Earth to the magnetopause and the time for the plasma moving from the subsolar point to the cusp, respectively. How the scaling changes the structure of reconnection is discussed in details in our companion paper by Tóth et al. (2017, submitted companion paper). We also apply a reduced ion-electron mass ratio $m_{i}/m_{e}=100$, which is sufficiently large to separate the electron and ion scales. We choose $\Delta x=1/32\,R_{E}$ as the PIC grid resolution so that $\Delta x/d_{i}\sim 5$ and $\Delta x/d_{e}\sim 0.5$. This resolution keeps a balance between the computational cost and the requirement of resolving kinetic scales. 216 particles per cell per species are used and there are about 9 billion particles in total inside the domain initially. Our numerical experiments suggest smoothing the electric field $\mathbf{E}$ and the current density $\mathbf{j}$ can help to suppress the numerical noise (Tóth et al., 2017, submitted companion paper).

The typical magnetic field strength in the magnetosheath is about 30 nT, and the corresponding ion gyro-frequency is $\Omega_{ci}=0.0286$Hz and $\Omega_{ce}=2.86$Hz with scaled charge-mass ratio. As mentioned above, the time step of iPIC3D is determined by the accuracy condition (Eq. 3). From the simulation, we find the maximum thermal speed of electrons inside the PIC domain is about 2500km/s, which leads to a time step of $\Delta t\sim 0.03s\sim 10^{-3}\Omega_{ci}^{-1}\sim 0.1\Omega_{ce}^{-1}$ with cell size $\Delta x=1/32\,R_{E}$. Therefore, the time step is small enough to resolve the gyro-motion of both electrons and ions.

2.3 Coupling between BATS-R-US and iPIC3D

BATS-R-US and iPIC3D are coupled through the Space Weather Modeling Framework (SWMF). These two models are compiled together to generate a single executable file. Both models can run simultaneously on specified processors and the information exchange is paralleled and handled by the Message Passing Interface (MPI). The details of the two-way coupling has been described by Daldorff et al. (2014).

In the simulation presented in this paper, we run the Hall MHD code first with the local time stepping scheme to reach a steady state. Then BATS-R-US sends the information, including density, velocity, pressure and magnetic field, to iPIC3D. iPIC3D initializes the electric field based on the Ohm’s law. The Maxwellian distributed particles are generated according to the fluid information so that iPIC3D and BATS-R-US have consistent density, velocity and pressure at the same position. After the PIC initialization, the MHD and PIC models update independently with their own time steps. The coupling frequency between these two models can be set to a value that is independent of the MHD or PIC time step. During the coupling, iPIC3D calculates moments of the particle distribution function, such as the density, velocity and pressure, and overwrites the MHD cells overlapped with PIC region. In return the MHD model provides electromagnetic field as well as particle boundary conditions for iPIC3D. For the particle boundary, iPIC3D removes the particles in the boundary cells, and re-generates new particles based on the fluid variables obtained from MHD. Between the two coupling time points, iPIC3D uses the latest information obtained from BATS-R-US as boundary condition during each iteration. In the simulation presented here, the time step for BATS-R-US and iPIC3D are around $\Delta t_{MHD}=0.015$ s and $\Delta t_{PIC}=0.032$ s, respectively. The coupling time interval is set to a small value $\Delta t_{couple}=0.005$ s so that MHD and PIC are coupled every time step. We note that the time step of PIC is even larger than the MHD, because the MHD time step is limited near the magnetic poles due to the high Alfven speed, while these regions are outside the PIC domain.

We used to generate particles in only one ghost cell layer (Daldorff et al., 2014) as particle boundary condition. Our numerical experiments suggest using more layers (5 layers specifically in this paper) as the particle boundary, while the electromagnetic field boundary is still only enforced at the outermost layer, is helpful to smoothly transit from PIC to MHD. The MHD cells overlapped with the PIC particle boundary are not overwritten by PIC. Similar technique has been used to implement open boundary condition for stand-alone PIC simulations (Peng et al., 2015).

We run the simulation on 6400 processors for 170 hours to model one hour simulation time on Blue Water supercomputer (Bode et al., 2012). iPIC3D and BATS-R-US use about $80\%$ and $15\%$ of the simulation time, respectively. The coupling and other overheads use the remaining $5\%$.

2.4 Energy conservation

Even though the PIC region is not a closed system, therefore mass and energy flow into and out of the region, it is still important to check the energy variation during the simulation to make sure the PIC model does not suffer from numerical heating or cooling. The normalized energy changes are shown in Figure 2. Throughout the simulation, the total energy $E_{t}$ variation is less than $3\%$. The small variation suggests that the numerical heating or cooling are insignificant. The initial condition for iPIC3D is under MHD equilibrium, but not necessarily under Vlasov equilibrium. The electromagnetic field energy $E_{EM}$ and kinetic energy of each species normalized by the initial total energy are also shown in Figure 2. During the first several minutes, energy is transferred from the particles to the electromagnetic field. After 200s, the ion and electron energy decreases about $5\%$, while the electromagnetic field energy increases from 0.3 to about 0.36. This is the transition from the MHD steady state to a PIC preferred solution. Further change of these energies are gradual and small. $E_{EM}$ is mainly magnetic field energy, which is about 3 orders larger than electric field energy.

3 Results

3.1 Overview

The iPIC3D code is initialized from a steady Hall MHD state, which is shown in Figure 1. The steady state is obtained from the Hall MHD run by using a local time stepping scheme, and a reconnection X line already exists near the equatorial plane along the dayside magnetopause. Since the local time stepping scheme is diffusive in this case, the reconnection signature near the X line is weak, for example, the Hall magnetic field strength is only about 1 nT. The PIC code inherits the magnetic field topology and starts evolving based on Maxwell’s equations and the motion of the macro-particles. An overview of the evolution of the dayside magnetopause is shown in Figure 3, which contains the Hall magnetic field $B_{y}$ and the field lines at the meridional plane inside the PIC box. At $t=70$ s, $B_{y}$ has already increased to about 8 nT. The Hall field extends far away from the X line with roughly the same field strength for each branch. 15s later, south of the existing reconnection point, another X line starts to form at around $x=10.2\,R_{E}$ and $z=-1\,R_{E}$. At $t=145$ s, both X lines can be seen clearly, and a flux rope like structure forms between the two X lines. The top X line has moved to about $z=0.5$. The bottom X line is almost steady so far, but it will move southward later. At $t=325$ s, the top and bottom X lines reach about $z=1.8$ and $z=-3.5$, respectively, and the center of the flux rope is moving southward with the bottom X line. Since the flux rope is moving away from the top X line, the current sheet between them becomes unstable and a secondary flux rope is generated (rightmost panel of Figure 3). During the one hour simulation, flux ropes form near the subsolar point and move toward poles quasi-periodically. More details about reconnection and flux ropes, both macroscopic and microscopic scales, are discussed in the following sub-sections.

3.2 Evolution of FETs

The meridional cut of the first two FTEs formed in the simulation are already shown in Figure 3. When we go beyond the 2D view, more complicated but completed structures arise. The flux ropes colored with the ion velocity z component $u_{iz}$ at different times are shown in Figure 4. At $t=100\,\mathrm{s}$, a short flux rope appears near the subsolar point. It is labeled as FTE-A. This flux rope extends from $y\sim-1\,R_{E}$ to $y\sim 1\,R_{E}$ in the dawn-dusk direction. It suggests that next to the primary X line near $z=0$, another X line starts to form south of the subsolar point. We have checked a series of 2D $x-z$ plane cuts, and found that the signature of reconnection, like the ion jets, at the second X line is clear at $y=0$, but appears very weak far away from the Sun-Earth line, for example at $y=0.78\,R_{E}$ or $y=-0.78\,R_{E}$. At $t=150$ s, the flux rope has extended significantly in both dawn and dusk directions. Along the flux rope, the ion velocity varies. Close to the dusk side (positive y), the reconnection at the second X line produces fast northward ion jet flow to slow down the southward flow from the primary X line, so that the flux rope moves relatively slowly. Close to the dawn side (negative y), the reconnection at the second X line is not strong enough to offset the southward flow ejected from the primary X line. The varying ion velocity leads to an inclined flux rope. At $t=240$ s, the flux rope is even more tilted because of the varying ambient ion jet velocity. A new small flux rope, FTE-B in Figure 4, is generated at $t=320$ s above FTE-A. FTE-A bifurcates at $y\sim-2.5$ and the new branch extends along the dawn-northward direction. FTE-A keeps moving southward while FTE-B is growing. At $t=540$ s, a large portion of FTE-A, except for the dawn part, already moves to the southern edge of the PIC domain ($z=-6$). FTE-B elongates significantly along the dawn-dusk direction. It is twisted at the dawn side so that the axial direction is almost parallel to the z-axis. At the dusk side, FTE-B connects to a newly formed flux rope FTE-C. At $t=660$ s, FTE-B and FTE-C have merged and become indistinguishable. These 3D plots suggest:1) flux ropes arise from multiple X line reconnection and can grow in time along the dawn-dusk direction, 2) the pole-ward moving velocity varies along a flux rope and makes them tilted, and 3) two flux ropes can merge and form a new long rope.

Since the PIC code is two-way coupled with the MHD model, the flux ropes can smoothly move out of the PIC region. Figure 5 shows a series of $j_{y}$ and field lines of FTE-A in the meridional plane after it leaves the PIC domain. FTE-A moves southward along the magnetopause after it is generated near the subsolar point. At $t=600\,\mathrm{s}$, the flux rope is already close to the southern cusp. There is strong axial current $j_{y}\sim 0.02\mu$A/m2 near the center of the flux rope. As FTE-A moves toward the cusp, $j_{y}$ inside the flux rope decreases in intensity, which indicates the dissipation of the magnetic helicity, as we can see at $t=660\,\mathrm{s}$. When the FTE reaches the center of the cusp ($t=720\,\mathrm{s}$), the field lines at the leading edge of the FTE and the cusp field lines are anti-parallel and creates a narrow and short current sheet with negative $j_{y}$ around $x\sim 4\,R_{E}$ and $z\sim-9.5\,R_{E}$. The ion velocity $u_{iz}$ at $x=4\,R_{E}$ in Figure 6 shows a jump around $z=-9.5\,R_{E}$. The narrow current sheet and the velocity jump imply that reconnection occurs between the flux rope field lines and the cusp field lines. At $t=840\,\mathrm{s}$, after FTE-A leaves the cusp, the signature of the flux rope becomes very weak: even though the magnetic field is still perturbed, the $j_{y}$ component is close to zero near the center and no ‘O’ line can be found. Finally, the remnant of the flux rope completely disappears as it moves toward the tail. Here we show the dissipation of FTEs in the meridional plane. But FTEs were observed along the distant tail magnetopause ($x=-67\,R_{E}$) on the dusk flank (Eastwood et al., 2012). One possibility to explain the conflict is that these survived FTEs may bypass the cusps and move along the flank from the dayside to tail magnetopause.

3.3 Magnetic field signature

Since the most widely used indicator of FTEs in satellite data is the magnetic field signature, we discuss how the flux rope magnetic field looks like along a virtual satellite trajectory. A series of meridional cuts are shown in Figure 7 to illustrate the magnetic field evolution. At $t=290\,\mathrm{s}$, north of the FTE-A event, there is an X line at about $z=1\,R_{E}$ surrounded by the quadrupolar Hall magnetic field $B_{y}$. As expected, the two branches on the magnetosheath side with amplitude of $\sim 30$ nT are stronger than the other two on the magnetosphere side with amplitude of $\sim 10$ nT. Near the X line, the magnetosheath and magnetosphere are separated by a current sheet accompanied with very weak magnetic field. 30 s later, another X line near $z=0$ arises, and an O line forms between the two X lines. Around the edge of the O lines, the azimuthal component of the magnetic field grows, while the $B_{y}$ component is still very weak just near the center. We note that the strong field on the magnetosheath side of the flux rope is mainly contributed by the $B_{z}$ component because of the accumulation of the inflow magnetic field lines. The reconnection at the northern X line is stronger than the southern one, so the ion jet around the O line is moving southward with a slow speed less than $100\,\mathrm{km/s}$. Inside the O line, the pressure starts increasing. $100\,\mathrm{s}$ later, the pressure at the center of the flux rope already reaches about $1.3\,\mathrm{nPa}$ while the core field is still small. At $t=540\,\mathrm{s}$, the O line structure continues to grow as the two X lines move northward and southward, respectively. We can see the core field $B_{y}$ at the center of the O line has grown to a significant value of $\sim 30\,\mathrm{nT}$ now, while the center pressure drops to $\sim 1.0\,\mathrm{nPa}$. The converging jets from the two X lines are comparable and the flux rope is almost steady. $180\,\mathrm{s}$ later, the core field grows to $\sim 40\,\mathrm{nT}$ and the corresponding pressure drops to about $0.8\,\mathrm{nPa}$. The whole structure at this stage is moving northward driven by the ion jet generated by the southern X line. To demonstrate the scaling factor has weak influence on the global structures, we perform another simulation with ion inertial length increased by a factor of 32. The simulation results are shown in Figure 8. The FTE in Figure 8 shows similar dynamic process as the event in Figure 7: the core field grows gradually and the ion pressure is anti-correlated with the core field strength. The FTEs in Figure 8 and Figure 7 also have comparable sizes.

At the early time when the O line just formed, for example, at $t=420\,\mathrm{s}$, the weak core field is surrounded by relatively large toroidal fields. We argue that this is an example of the so-called ’crater FTEs’ that have been observed by spacecrafts (LaBelle et al., 1987; Zhang et al., 2010). Since the O line moves slowly during its initial stage of formation, the magnetic field observed at a fixed point can not reflect its global structure. Instead, the magnetic field along the magnetopause (the red curve in the left panel of Figure 9) is shown in the right panel of Figure 9 to illustrate its magnetic field structure. Along the magnetopause, from south to north, the $B_{x}$ field, which is roughly normal to the magnetopause, reaches a local minimum of $\sim-15$ nT at $z=0$ and then quickly increases to $\sim 15$ nT at $z=1\,R_{E}$. The flux rope is bounded by the depressed magnetic field ‘trenches’ at $z=-0.2\,R_{E}$ and $z=2\,R_{E}$ as indicated by $B_{t}$. The depression results from the low magnetic field strength inside the current sheet as can be seen from the right panel of Figure 9. $B_{t}$ reaches local maximum at the same position of the $B_{x}$ peaks ($z=0\,R_{E}$ and $z=1\,R_{E}$), while the field strength decreases to about $10$ nT between the peaks. We refer to the event on 30 July 2007 observed and analyzed by Zhang et al. (2010) as a comparison. Figure 6 of Zhang et al. (2010) shows the magnetic field signature of this event. Even though the 30 July 2007 event has a large guide field (corresponding to $B_{y}$ component in our simulation), and its magnetic field around the flux rope is more steady than our simulation, the whole structure of this event is similar to what is shown in Figure 9.

As the flux rope evolves, the core field strength grows to a significant value. The magnetic field measured at a fixed position $x=10.2\,R_{E}$, $z=2.75\,R_{E}$ is shown in the right panel of Figure 10. The vertical dashed line at $t=760\,\mathrm{s}$ represents the location of the maximum $B_{t}$. Around this time, the $B_{x}$ field, which is roughly perpendicular to the magnetopause, jumps from $\sim 5$ nT to $\sim-20$ nT within about 25 s. At $t=760\,\mathrm{s}$, both the axial field $B_{y}$ and the total field $B_{t}$ reach a maximum. These features match the signatures of an FTE with typical flux rope structure (Zhang et al., 2010). During the one-hour long simulation, there are ten FTEs with significant core field moving across the southern PIC edge. The occurrence frequency is consistent with observations (Rijnbeek et al., 1984) and previous MHD simulations (Raeder, 2006).

The IMF is purely southward in our simulation and there is no uniform background guide field at the magnetopause. But a significant core field can still arise during the FTE generation and evolution as seen in Figure 7. When a flux rope is still close to the X lines, the core field may be encompassed by the Hall magnetic field generated by the reconnection, resulting in complicated guide field structure. The $B_{M}$ field at $t=540\,\mathrm{s}$ is shown in Figure 11. In order to compare with observations, the magnetic field has been transformed into a boundary normal coordinate system ($\mathbf{LMN}$), in which the $\mathbf{N}$ component points outward, normal to the magnetopause, the $\mathbf{M}$ component is determined by $\mathbf{N}\times\mathbf{Z}_{GSM}$ and the $\mathbf{L}$ component completes the right-hand coordinate system. Since the plot is shown in the meridional plane, the $\mathbf{Y}_{GSM}$ direction is anti-parallel to the $\mathbf{M}$ direction. Around the flux rope center, the guide field $B_{M}$ is negative, while the southern part of this flux rope is surrounded by positive $B_{M}$. The polarity of the positive ’Y’ shape $B_{M}$ is consistent with the Hall magnetic field generated by the X line at $z=-1\,R_{E}$. If a satellite is moving across the flux rope along the red solid line in the left panel of Figure 11, the satellite will observe a tripolar guide field structure (right panel of Figure 11). Similar structure was first observed in the solar wind (Eriksson et al., 2015), and then it was observed by the Polar satellite at the magnetopause (see Figure 1 of Eriksson et al. (2016)). The Polar event shows a large negative $B_{M}$ core field bounded by two narrow $B_{M}$ depressions in the presence of a large background guide field. There is no background guide field in our simulation and thus the right panel of Figure 11 shows a pure tripolar structure: the large negative $B_{M}$ field is surrounded by two relative small positive peaks. Despite the difference in the background guide field, the topology of $B_{M}$ obtained from our simulation is very similar to the Polar observation.

3.4 Kinetic features

We have examined the global structure of the FTEs in the previous discussion. In this subsection, we will demonstrate that the underlying kinetic physics is properly captured by our model. The Larmor electric field, identified by Malakit et al. (2013), is a localized electric field that appears on the magnetospheric side of the dayside reconnection site. The x-component of the electric field $E_{x}$ at the end of the simulation (t=3600s) is shown in Figure 12. The positive $E_{x}$ pointing towards the Sun along the magnetopause is the Hall electric field, while behind the Hall electric field, the localized negative field pointing towards the Earth is the Larmor electric field. A 1D cut through the reconnection site along the x direction is also shown in Figure 12. The Larmor field strength is -3 mV/m, the magnetospheric side ambient field is about 2 mV/m, and the nearby Hall field is about 12 mV/m. These values are reasonably close to the MMS observation by Graham et al. (2016): the Hall electric field strength is $\sim 20\,\mathrm{mV/m}$ and the Larmor field strength is about 10 mV/m (see Figure 2 of Graham et al. (2016)).

Even though the ion inertial length is scaled up by a factor of 16 in the present simulation, the electric field strength is not sensitive to the scaling factor. Ignoring the electron inertia term, the generalized Ohm’s is:

[TABLE]

Tóth et al. (2017, submitted companion paper) shows the electron velocity $\mathbf{u_{e}}$ of the current sheet does not change with the scaling factor while the current sheet width scales. The gradient of electron pressure is inversely proportional to the scaling factor, because the pressure jump is fixed across the current sheet and the current sheet width is proportional to the scaling factor. Since the charge per ion or electron is also reduced by the same factor, the scaling does not change the electric field strength. Besides the scaling of the ion inertial length, a reduced ion-electron mass ratio $m_{i}/m_{e}=100$ is used in this study to increase electron kinetic scales (see section 2.2). The influence of the mass ratio $m_{i}/m_{e}$ has been studied in numerous papers (Shay and Drake, 1998; Hesse et al., 1999; Ricci et al., 2004; Shay et al., 2007; Lapenta et al., 2010). For the Larmor electric field , Malakit et al. (2013) specifically estimates its amplitude to be:

[TABLE]

where $k_{B}$ is the Boltzmann’s constant, $T_{i}$, $q_{i}$ and $r_{i}$ are the temperature, charge per ion and ion Larmor radius of the ions on the magnetospheric side. In the simulation, $q_{i}$ reduces by a factor of 16 and $r_{i}$ becomes 16 times larger compared to the realistic situation, while the temperature $T_{i}$ does not change. So, the scaling of inertial length should not influence the strength of the Larmor electric field. On the magnetosheath side, our simulation shows the ion temperature is about $2\times 10^{6}\,K$, magnetic field strength is about 60 nT. Substituting these values into Eq. 5 gives $E\sim 5.5\,\mathrm{nT}$. As mentioned above, the value obtained from simulation is about -3 mV/m.

The crescent shape electron phase space distribution has been observed near the electron diffusion region at the dayside magnetopause by MMS (Burch et al., 2016). The same distribution is also found in our 3D global simulation. The phase space distribution of electrons inside a cube region: $10.27\,R_{E}<x<10.33\,R_{E}$, $-0.3\,R_{E}<y<0.3\,R_{E}$ and $-2.1\,R_{E}<z<-1.9\,R_{E}$ is shown in Figure 12. The crescent distribution is found in the $V_{y}-V_{x}$ plane, corresponding to the two velocity components perpendicular to the magnetic field. The crescent hot electrons are drifting along negative y direction with a speed close to 3000 km/s. The direction of the flow is consistent with the $\mathbf{E}\times\mathbf{B}$ direction, and the velocity of the crescent particles is very close to the MMS observation (Burch et al., 2016). Slightly further away from the reconnection site, where the Larmor field appears, inside a cube $10.08\,R_{E}<x<10.14\,R_{E}$, $-0.3\,R_{E}<y<0.3\,R_{E}$ and $-2.1\,R_{E}<z<-1.9\,R_{E}$, the ion phase space distribution also presents crescent like shape as it is shown in Figure 12(c). The crescent ions drift in positive y direction because $E_{x}$ is negative. We also checked the distributions for particles inside the current sheet but far from the reconnection site, and no crescent distributions are found for either electrons or ions.

Kinetic effects along the magnetopause current direction are also captured by our 3D MHD-EPIC model. Figure 13 shows the fully developed lower hybrid drift instability (LHDI) at the end of the simulation (t=3600 s) at the $z=-3\,R_{E}$ plane. The electric field $E_{M}$ shown in Figure 13 is the $\mathbf{M}$ component in the boundary normal coordinates, and $\mathbf{M}$ is anti-parallel to the current direction. The black curve in Figure 13 separates the negative and positive $B_{z}$. We can see the LHDI appears along the magnetopause on the magnetospheric side. A closer view of $E_{M}$, as well as $B_{z}$, ion mass density $\rho_{i}$ and electron velocity $u_{ey}$ are also shown Figure 13. It is clear to see the LHDI arising near the interface of magnetosheath and magnetosphere, where there is a sharp density gradient. $B_{z}$, $\rho_{i}$ and $u_{ey}$ show sawtooth pattern at the same location. The amplitude of the LHDI electric field is about 8 mV/m, which is consistent with MMS observation (Graham et al., 2016). The dominant wave length shown in Figure 13(b) is about $0.38\,R_{E}$, and the ambient magnetosheath side electron gyroradius is about $r_{e}=0.025\,R_{E}$ with the artificially changed charge per electron mass ratio, which results in $kr_{e}\sim 0.4$, where $r_{e}=m_{e}v_{e}/(q_{e}B)$ and $v_{e}$ is defined as $v_{e}=\sqrt{2T_{e}/m_{e}}$. The value of $kr_{e}$ is also consistent with observation (Graham et al., 2016) and theory (Daughton, 2003). LHDI at different time and different location is analyzed, the value of $kr_{e}$ varies from $\sim 0.3$ to $\sim 0.5$, and $kr_{e}\sim 0.4$ is a typical value. Similar as the argument above with the Ohm’s law, the electric field strength is not sensitive to the scaling, that is why the LHDI electric field strength agrees with MMS observations. But the length scale does change with the scaling. The charge per mass of electron $q_{e}/m_{e}$ is artificially increased by a factor of 294 in the simulation, and the electron thermal velocity reduces by a factor of $\sqrt{18.36}=4.3$ for $m_{i}/m_{e}=100$. The magnetic field is realistic, hence the electron gyroradius is about 68 times larger than in reality. If we scale back the LHDI wavelength of the simulation by the same factor, it will be $\sim 35$ km. As a comparison, MMS observed 10km $\sim 13$km wavelength (Graham et al., 2016). Figure 13(f) shows the isosurfaces of $E_{M}=4$ mV/m colored by the ion velocity $u_{iz}$ viewed from the Sun. Along the magnetic field direction, the isosurfaces are cut off two or three times. The ion velocity jumps or even change directions across a cut-off region. It suggests these cut-off regions corresponding to the reconnection sites and the LHDI electric field is weak near the diffusion regions (Pritchett, 2013).

3.5 Comparison with Hall MHD

For comparison, we also run a pure Hall MHD simulation with the same setup except the PIC region is removed and the MHD grid resolution around the dayside magnetopause is refined to $1/32\,R_{E}$, which is the resolution used by PIC in the MHD-EPIC run. Even for Hall MHD, resolving the ion inertial length is necessary in order to capture the Hall effect correctly. Due to the small kinetic scale inside the magnetosheath, scaling the ion inertial length is also required for a global Hall MHD simulation since Hall MHD is also computationally expensive as we will see. We note that the ion inertial length in the pure Hall MHD simulation is also scaled up by a factor of 16 so it can be resolved by $1/32\,R_{E}$ cell. Hall MHD is reasonably optimized by using semi-implicit scheme to overcome the time step imposed by the whistler mode wave and speed up the simulation. It still takes 6400 cores running about 67 hours to model one hour because of the high resolution and the stiffness of the Hall term. As a comparison, the MHD-EPIC simulation (170 hours on 6400 cores) is about 2.5 times more expensive. Hall MHD produces the Hall magnetic field near the X line and generates flux ropes in a way similar to MHD-EPIC. But Hall MHD can not reproduce the kinetic features, neither the crescent particle distributions nor the LHDI.

4 Summary and conclusion

We have performed a one-hour long high-resolution global simulation with the MHD-EPIC model to study dayside reconnection and FTEs. Our simulation is the first attempt to investigate the FTEs and reconnection with kinetic physics resolved in a realistic magnetopause environment. Although the kinetic scale is artificially increased to reduce the computational cost, the model still captures the kinetic features very well. MMS observations, like the crescent particle phase space distribution and LHDI, are reproduced in our model. The FTEs from the simulation also agree well with spacecraft observations. The key results from the present simulation are:

• •

When an FTE arises, its cross section is small and it is short in the dawn-dusk direction. During its growth, the cross section increases and the FTE extends along the dawn-dusk direction.

• •

An FTE forms near the subsolar point and moves toward the poles under steady southward IMF conditions. When the FTE reaches the cusp, reconnection happens between the FTE magnetic field and the cusp magnetic field lines, thus dissipating the FTE. The signature of FTE is weak behind the cusps.

• •

FTE is flanked by two reconnection sites during its formation, and the converging ion jets around the FTE are found.

• •

The present simulation confirms that the ’crater FTEs’ magnetic field signature can be found at the early stage of an FTE formation when the axial magnetic field is still weak. A strong core field may develop as the FTE evolves, and the Hall magnetic field may provides the initial seed core field. Therefore a fully developed FTE has the typical strong core field structure.

• •

A tripolar guide field structure is found from our simulation.

• •

The Larmor electric field is found near the reconnection site on the magnetospheric side, and its amplitude is about -3 mV/m.

• •

A crescent electron phase space distribution is found near the reconnection site where the Hall electric field reaches its maximum. A similar distribution is also found for ions at the place where the Larmor electric field appears.

• •

The lower hybrid drift instability (LHDI) appears at the interface of the magnetosheath plasma and magnetosphere plasma. The LHDI electric field peak strength is about 8 mV/m, and a typical ratio between its wavelength and the electron gyroradius is about $kr_{e}\sim 0.4$. The simulation agrees with the MMS observations and theory.

Compared to the models relying on ad hoc resistivity or numerical resistivity to generate FTEs or investigate reconnection process, our 3D MHD-EPIC model makes one significant step forward by incorporating a self-consistent kinetic description of reconnection into a global MHD model. While the kinetic scales are increased by artificially reducing the charge per mass for both ions and electrons, all the other parameters are realistic. The scaling changes the size of kinetic features, for example the wavelength of LHDI, but other values, like the strength of Larmor electric field or LHDI electric field, are not modified by the scaling. Another artificial change is the solar wind electron pressure. It is set to a value 8 times larger than the ion pressure so that $p/p_{e}\sim 2.5$ inside the magnetosheath while the ratio is usually about $4\sim 12$ from observation (Wang et al., 2012). The artificially increased electron pressure can help to stabilize the simulation, and it does not deviate significantly from the observed values. We plan to improve this in the future studies.

The MHD-EPIC model offers a powerful tool to study magnetospheric physics. The PIC code only covers the dayside magnetopause in the present simulation. As a natural extension, it can be elongated to cover the bow shock so that the kinetic processes associated with the bow shock can be modeled. Another future application is covering the tail reconnection site with another PIC region, so that both dayside and tail reconnections are handled by a kinetic code and we can study substorm in a more realistic way.

Acknowledgements.

This work was supported by the INSPIRE NSF grant PHY-1513379, NSF strategic capability grant AGS-1322543, NASA grant NNX16AF75G, NASA grant NNX16AG76G, and the Space Hazards Induced near Earth by Large, Dynamic Storms (SHIELDS) project DE-AC52-06NA25396, funded by the U.S. Department of Energy through the Los Alamos National Laboratory Directed Research and Development program. Computational resources supporting this work were provided on the Blue Waters super computer by the NSF PRAC grant ACI-1640510, on the Pleiades computer by NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, and from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. The SWMF code (including BATS-R-US and iPIC3D) is publicly available through the csem.engin.umich.edu/tools/swmf web site after registration. The output of the simulations presented in this paper can be obtained by contacting the first author Yuxi Chen.