Full-particle simulations of instabilities in a thin current sheet of the magnetospheric system prior to substorm onset

Abstract

Substorm triggering was the focus of recent two-dimensional superposed-epoch analyses based on Geotail and THEMIS data. The results indicate that magnetic-field dipolarization at X∼−8 Re and magnetic reconnection at X∼−20 Re occur simultaneously at the onset. These results imply that there are physical mechanisms that widely affect both the dipole and current-sheet regions. The analyses also have found that a local B _z enhancement appears before the substorm onset and magnetic reconnection occurs at its tailward edge. We performed four 2.5-dimensional full-particle simulations with a new initial magnetic-field structure to focus on instabilities in the magnetosphere. The structure is similar to the Earth’s dipole magnetic field combined with a stretched field and current sheet on the tailward side. The simulation with the initial magnetic-field configuration shows that nodes of the magnetic field appear in the current sheet where the growth condition of tearing instability is satisfied. The features of the instability are close to those of the electron tearing mode reported in previous simulation results. Another three simulations with a local B _z enhancement, as seen in the observational results, at various locations in the current sheet were performed to explore its impacts on the evolution of the instability. A relaxation process around the enhancement generates a new node at its tailward edge if its location satisfies the growth condition. The wavelength and dominant mode of the instability can be changed by the coupling between the process and tearing mode depending on the location of the enhancement. Our simulations reveal new features associated with tearing instabilities in the magnetospheric-field configuration.

Background

Magnetospheric substorms are manifestations of a general physical process that causes drastic changes in both the magnetosphere and ionosphere. Akasu (1964) first illustrated the time development of the auroral substorm by using data collected during the International Geophysical Year. Although a number of substorm-triggering models have been proposed on the basis of observational evidence (Lui 1991), none can as yet perfectly explain the mechanism responsible for substorm triggering.

One of the triggering models is the “near-Earth neutral line” (NENL) model (Baker et al. 1996; Hones 1979), which has been used successfully to explain energy transfer during substorms from the magnetic field to particles in the magnetotail. The cross-tail current becomes thinner as the convective electric field penetrates into the magnetosphere through the southward turning of the interplanetary magnetic field, and magnetic reconnection occurs at the NENL. The earthward flow generated by the magnetic reconnection injects energetic particles into the inner magnetosphere and causes dipolarization (Shiokawa et al. 1998). Another triggering model is the “current disruption” (CD) model (Chao et al. 1977; Lui 1996), which adopts CD as the cause of the field-aligned current flowing into the Earth’s ionosphere. The process is initiated by a cross-field current instability at X∼−8 Earth radii (Re; 1 Re = 6378 km). A rarefaction wave ejected from the CD region propagates tailward and induces current-sheet (CS) thinning. A magnetic reconnection arises at X∼−20 Re as the consequence of the process.

Machida et al. (2009, 2014) performed superposed-epoch analyses with Geotail and THEMIS data to investigate the physical conditions in the magnetotail around the substorm onset. They have suggested a new substorm-triggering model, namely, the “Catapult Current Sheet Relaxation” (CCSR) model, based on the results. This model has the following characteristics: (1) magnetic-field dipolarization and magnetic reconnection occur simultaneously (within 1–2 min) with auroral breakup; (2) from 4 min before substorm onset, the CS relaxes around X∼−13 Re; and (3) a local B _z enhancement coincides with the CS relaxation at X∼−17 Re and magnetic reconnection occurs on the tailward edge of the enhancement. These results are indicative of the need for physical mechanisms that cover a wide area and span the range from the dipole-like (DL) region to the CS region.

A number of particle simulations have been performed to study the physical processes associated with magnetic reconnection in order to investigate the causes of substorm triggering (Birn et al. 2001; Hesse et al. 1999; Kuznetsova et al. 2001). Particle simulations, where the initial magnetotail-like structure is described by a slightly modified Harris with a normal magnetic field to the CS, have led to the conclusion that tearing modes may lead to magnetic reconnection under specific physical conditions (Dreher et al. 1996; Zwingmann et al. 1990). Pellat et al. (1991) indicated that electrons magnetized by the normal magnetic field to the CS suppress the tearing mode with realistic physical values. Recent simulation studies that impose the convective electric field as a boundary condition at the magnetospheric lobe have shown that tearing instability and magnetic reconnection can be driven (Liu et al. 2014; Pritchett 2005). Sitnov and Swisdak (2011) showed that the “X-point” (the point where oppositely directed magnetic-field lines reconnect) of magnetic reconnection forms on the tailward side of the B _z -hump in the CS, as was found by Machida et al. (2009; 2014), because the electric energy of the Poynting flux is locally concentrated. The observational results reported in Machida et al. (2009; 2014), however, require simulations of both the DL field and the CS region.

Therefore, we performed four 2.5-dimensional full-particle simulations with a new initial magnetic-field structure in order to focus on instabilities in the magnetosphere. The structure is similar to the Earth’s dipole magnetic field combined with a stretched magnetic field and CS on the tailward side. The first simulation does not have a local B _z enhancement, as seen in the observational results (Machida et al. 2009, 2014), in its initial condition. We discuss the process leading to the appearance of instabilities, as well as the growth rate and wavelength of the modes resulting from such inhomogeneous initial magnetospheric-field conditions. We also conducted three other simulations with a local B _z enhancement at various locations in the CS. Based on the simulation results, we discuss the effects of the presence of the B _z enhancement on the resulting instabilities.

Simulation settings

The 2.5-dimensional full-particle simulations were performed in the presence of a characteristic magnetic-field structure corresponding to a DL field region that is connected to the magnetotail region. We have set periodic boundary conditions in both the x (Sun–Earth) and z (north–south) directions of the geocentric solar magnetospheric coordinate system. Our “2.5-dimensional” setup means that two dimensions cover real space, while the particle-velocity and electric/magnetic-field vectors have three dimensions. To quantify the particle motions, the Lorentz equation is solved by using the Buneman–Boris method. Maxwell’s equations are solved based on the fast Fourier transform (FFT) method (Birdsall and Langdon 1985). In this paper, we refer to the positive (negative) x direction as “earthward” (“tailward”), in accordance with the terrestrial magnetotail picture.

Figure 1 a shows the magnetic-field structure, which is defined by the following vector potentials:

$$ A_{y}(x,z)=A_{y0}(x,z)+A_{y1}(x,z), $$

((1))

Fig. 1

Initial configurations of the simulations. a Initial condition of the magnetic-field line selected in the lower left of the simulation system and b in the overall system of Run(0). c Hot-particle distribution (n _h) determined by the numerical “quiet start” method

$$\begin{array}{*{20}l} A_{y0}(x,z)=& B_{0}d\ln\left[\frac{\cosh\left(\frac{z}{d}\right)+p\cos\left\{2\pi\left(-\frac{x}{2\pi d} +\frac{1}{2}\right)\right\}}{\sqrt{1-p^{2}}}\right]\\ & \qquad (-\pi < x/d < 0), \end{array} $$

((2))

$$\begin{array}{*{20}l} A_{y0}(x,z)=&B_{0}d\ln\left[\frac{\cosh\left(\frac{z}{d}\right)+p}{\sqrt{1-p^{2}}}\right]\\ &\qquad (-4\pi < x/d < -\pi), \end{array} $$

((3))

$$ A_{y1}(x,z)=-\epsilon_{1}B_{0}D\exp\left(-\frac{z^{2}}{d^{2}}\right) \qquad (-\alpha/d < x/d < 0), \\ $$

((4))

$$\begin{array}{*{20}l} A_{y1}(x,z)=&-\epsilon_{1} B_{0}\left\{\!-C\left(\frac{x+\alpha}{d}\right)+D\!\right\}\exp\left\{E\left(\!\frac{x+\alpha}{d}\!\right)\right.\\ &\quad-\left.\frac{z^{2}}{d^{2}}\right\} (-4\pi < x/d < -\alpha/d), \end{array} $$

((5))

where d is the CS half-thickness. We use p=0.4 (0≤p<1, this value controls the roundness of the DL magnetic field; when p=0, the magnetic-field structure is the Harris solution, for example), so that the maximum B _z field in the DL region is about a half of the magnetic field in the lobe region (B ₀=1). Other parameters are ε ₁=0.03 (the ratio of A _y1 to A _y0), C=10d/3, E=3/4, D=C/E=40d/9, and α/d=π−1/E∼1.8. The main field given by A _y0 corresponds to Earth’s DL magnetic field (−π<x/d<0) and the magnetic field of the CS in the magnetotail (−4π<x/d<−π). The field defined by A _y1 corresponds to the finite northward-directed component of the magnetic field in the CS; its effects are described in this section.

The DL field of A _y0 is represented by a modification of the “Kelvin’s cats eyes (KCE)” magnetic field (Schindler 2007), which satisfies the magnetohydrostatic (MHS) conditions in two-dimensional (2D) space. This KCE field is obtained by setting $f(\zeta) = \{\exp (\mathrm {i}\zeta)+p\}/\sqrt {1-p^{2}} \,\,\,\, (\zeta = x + {\mathrm i}z)$ in the analytical solution of the Grad–Shafranov equation, i.e.,

$$\begin{array}{@{}rcl@{}} A_{y}(x,z) = \ln\frac{1+\left|f(\zeta)\right|^{2}}{2\left|\frac{{\mathrm d}f}{{\mathrm d} \zeta}\right|}, \end{array} $$

((6))

under the condition that particles follow a Maxwell velocity distribution. The KCE field has the following functional form:

$$\begin{array}{@{}rcl@{}} A_{y}(x,z) = \ln\left(\frac{\cosh(z)+p\cos(x)}{\sqrt{1-p^{2}}}\right). \end{array} $$

((7))

There are two advantages in regard to the use of the KCE field. First, because the B _z component derived from Eq. (2) always becomes 0 at x/d=−π along z, we can easily define the magnetic configuration of the CS, which is connected to the DL field, by keeping the values of the B _x and B _z components constant on the tailward side of x/d=−π. Equation (3) expresses the constant B _x and B _z in the x direction. Second, periodic boundary conditions can be applied to the KCE field because B _z =0 and ∂B _z /∂x=0 at x=0. The magnetic configuration of the initial conditions in the overall simulation system is shown in Fig. 1 b. This configuration is obtained by reversing the magnetic field shown in Fig. 1 a (which displays the lower left quadrant of the simulation system) and applying this to the entire region. Figures of simulation results for the simulation system display only the same region as Fig. 1 a because the residual region had almost symmetrical structures owing to the periodic boundary conditions.

The physical parameters are defined as follows. The ion-electron mass and temperature ratio are m _i/m _e=100 and T _i/T _e=4, respectively. The ratio of the electron gyro-frequency to the electron plasma-frequency is Ω _e/ω _pe=eB ₀/(m _e ω _pe)=1. We adopt d=0.5δ _i=5c/ω _pe (δ _i: average ion inertial length in the simulation system) as the CS half-thickness. The spatial scales of the simulation system are normalized by d. This simulation box has a size of L _x /d×L _z /d=8π×8π and covers 1024 × 1024 grid points, but the region shown in our figures is restricted to half of those ranges in both the x and z directions (L _x /d×L _z /d=4π×4π). The grid lengths Δx and Δz are constrained to Δx=Δz=0.1227c/ω _pe. The average number of particles in each cell is 25, where we adopt a ratio of the total number of hot to cold particles of 3:2. The cold particles are distributed homogeneously across the overall simulation system. We have used the relation n _h (x,z)∝ exp{−2A _y (x,z)} (Toichi 1972) to determine the positions of the hot particles carrying the diamagnetic current, which leads to the characteristic magnetic-field configuration. To assign positions of particles, we applied the “quiet start method” (Birdsall and Langdon 1985) to the relation and solved it by inversion. The cumulative distribution function D(x,z) of the density function n(x,z) is divided equally with respect to the total particle number N _p, and each particle position (x _i ,z _i ) (i: particle number, 1≤i≤N _p) can be obtained by solving the inverse of the relation D(x _i ,z _i )=i/N _p. This procedure was performed numerically because it could not be done analytically because of the 2D setup. Figure 1 c shows the hot-particle density across the entire simulation system resulting from the procedure. In regard to the a velocity distribution for hot particles, we used the shifted Maxwellian, so the distribution function is

$$\begin{array}{@{}rcl@{}} f_{s} = \frac{n(x,z)}{\pi^{3/2}v_{\text{th}_{s}}^{3}}\exp\left[-\frac{{v_{x}^{2}}+\left(v_{y}-U_{\mathrm{s}}\right)^{2}+{v_{z}^{2}}}{v_{\text{th}_{s}}}\right], \end{array} $$

((8))

where U _s is the drift velocity and $\phantom {\dot {i}\!}v_{{\text {th}}_{s}}$ is the thermal velocity. U _s is calculated by using Ampere’s law, en(U _i−U _e)=μ ₀∇×B, and the charge neutrality, U _i/T _i=−U _e/T _e. Furthermore, in order to quickly induce instabilities in the CS without unwanted electric-field effects in the specific simulation system, we have set smaller particle pressures than those of the magnetic field by a factor of 3, $v_{\text {th}_{s}}^{2}\equiv 2T_{\mathrm {s}}/3m_{\mathrm {s}}$. This setting is also effective to focus on couplings between the instabilities and relaxation process when a B _z local enhancement is included; this will be discussed further in “Discussion and conclusions” section.

Figure 2 a shows the B _z component in the magnetic equatorial plane (z/d=0). The B _z component of the DL magnetic field, B _z0, which is defined by A _y0, is shown as the purple line. B _z1, defined by A _y1, expresses the small positive normal B _z component that is present in the real magnetospheric CS. It is indicated by the green line. There is a discontinuity in the B _z differential at the boundary between the DL and CS regions (x/d=−π), which we henceforth refer to as the “DCB” (DL–CS boundary). This discontinuity is formed by connecting the DL field and CS regions so that B _x is continuous at the location. The magnetic field in the CS region does not perfectly satisfy the MHS condition originally (i.e., it deviates from the Harris solution), so the transient adjustment of the CS initially appears in the E _y field. This discontinuity may be seen as unnatural, but it has been shown in a magnetohydrodynamical calculation that such a discontinuity can appear at the inner plasma-sheet boundary through penetration of the convective electric field during a substorm’s growth phase (Voigt 1986). B _z1 reaches its maximum value ∼ 0.04 near the DCB and decreases toward the tail side. Electrons magnetized by normal magnetic field B _z1 suppress the growth tearing mode (Pellat et al. 1991), so the region near the DCB is stable, while the stabilizing effect is less at the far tail. If we do not add B _z1 to the simulation system, magnetic reconnection occurs at the location of the B _z differential discontinuity (x/d=π). The B _z1 is useful to stabilize this reconnection, which is unimportant in the context of this study.

Fig. 2

Initial conditions of the B _z component in the equatorial plane (z/d=0). a–d show the B _z component in Run(0)–(3), respectively. Purple lines indicate the magnetic field of the DL region B _z0, green lines indicate the small northward magnetic field in the CS B _z1, red lines express the enhancement of B _z2, and black lines show the total of B _z

A time step of the simulation is ω _pe Δt=0.2, and we ran the simulation for 12,000 steps (Ω _i t=24). We call this simulation “Run(0),” since it serves as the reference simulation for the other runs in this study.

Observations of a magnetospheric CS during a substorm show that a region of local B _z enhancement appears at X∼−17 Re a few minutes prior to substorm onset, and a magnetic reconnection occurs at the tailward edge of the enhanced region at the substorm onset (Machida et al. 2014). To investigate the effect of a local B _z enhancement in the CS, we have included such an enhancement in the CS of Run(0). We use the following vector potential, A _y2, to introduce the local B _z enhancement:

$$\begin{array}{*{20}l} A_{y2}(x,z)&=-\epsilon_{2}B_{0}\gamma d\exp\left(-\frac{z^{2}}{d^{2}}\right)\\ &\qquad (-l_{x} < x/d < 0), \end{array} $$

((9))

$$\begin{array}{*{20}l} A_{y2}(x,z)&=\epsilon_{2}B_{0}\gamma d\tanh\left(\eta\ln\left[-\frac{x+l_{x}}{\gamma d}\right]\right)\\ &\qquad\exp\left(-\frac{z^{2}}{d^{2}}\right) \qquad (-4\pi < x/d < -l_{x}/d), \end{array} $$

((10))

where ε ₂=0.008 (the ratio of A _y2 to A _y0), η=4 (this parameter controls the steepness of the local B _z enhancement), and γ=0.6 (the width of the local B _z enhancement). Figure 2 b–d shows the initial B _z configurations in the equatorial plane with the local B _z enhancement (B _z2; red line), in addition to the initial conditions of Run(0). The maximum value of the local enhancement was set to approximately 10 % of the maximum of B _z0. We also set the width of the enhancement so as to constrain it to approximately one third of the width of the DL region, again based on observational results (Machida et al. 2014). The local maximum of B _z2 appears at $\phantom {\dot {i}\!}x/d=-\gamma \{(2\eta -1)/(2\eta +1)\}^{{1/2\eta -l}_{x}}$. This value asymptotically approaches to x/d≃−γ−l _x as η becomes greater than unity. We performed more than 10 simulation runs, but here we only show 3 relevant runs where the maximum location of the enhancement x _h /d has been shifted from x _h /d∼−3.6 [Run(1)] to x _h /d∼−6.8 [Run(3)] by changing the term l _x in intervals of x/d∼−1.6 (see Table 1). The overall structures of the magnetic fields in the initial conditions in Run(1)–(3) are hardly distinguishable from that in Run(0). We also focus on the effect of a local B _z enhancement on instabilities, which depends on the locations where the enhancement is seeded in the CS.

Table 1 Locations of the local B _z enhancement x _h during initial conditions

Results

Instability in the current sheet

First, we mention that the stronger magnetic-field pressure causes CS thinning. Figure 3 shows the J _y profiles along z at x/d=−8 in Run(0) at times Ω _i t=0.4 and Ω _i t=4.0 from the start of the simulation. J _y in the simulation, which is shown as the black line, initially follows the Harris solution with a half-thickness d=5c/ω _pe (blue line). Because of compression due to the pressure difference, it proceeds to follow the thinner Harris profile for d=1.6c/ω _pe (red line) as seen in Fig. 3 b. The plasma pressure in the CS is 3 times smaller than that of the magnetic field, so the effective CS half-thickness has become d ^′=1.6c/ω _pe, which is about one third of the initial CS half-thickness. The CS in other runs becomes similarly thin.

Fig. 3

Comparisons of J _y in the simulation Run(0) to the Harris solution along z. a J _y (x/d=−8) profile in the simulation (black line) follows that of the Harris solution for d=5 (blue line) at the initial stage of the simulation Ω _i t=0.4. b J _y (x/d=−8) profile in the simulation (black line) follows the Harris solution for d=1.6 (red line) at the time Ω _i t=4.0

Figure 4 a, b displays the B _z component in the equatorial plane (z/d=0), as well as the magnetic-field line and poloidal current, J _y , at times Ω _i t=8.4 and Ω _i t=15.2 in Run(0). The B _z value shown in Fig. 4 a is averaged along the z direction, i.e., for −0.5≤z/d≤0.5, to eliminate small perturbations. Nodes of B _z have appeared at x/d∼−5.8 and −9.4 in the CS. When we refer to a “node,” we mean a pinched point of the magnetic-field structure in the CS. Such points can be defined exactly as points where B _z goes from negative to positive along the tail axis. We have estimated the wavelength of the mode as λ/d∼4 based on the distance between the two nodes. This wavelength is also rewritten as kd ^′∼0.5, where k≡2π/λ by adopting the thinned CS thickness d ^′. This wavelength is close to that reported by Bessho and Bhattacharjee (2014).

Fig. 4

Simulation results and growth rate analysis of Run(0). a B _z (z/d=0) and b magnetic-field line and poloidal current J _y (color) at the times Ω _i t=8.4 and Ω _i t=15.2. c B _z (z/d=0) averaged during the period 0≤Ω _i t≤4.0 in Run(0). d Time developments of the power of the modes resolved for B _z in Run(0)

We have also calculated the growth rate of the modes that appeared in Run(0) by the following procedure.

• 1. Calculate the geometric average of the B _z (z/d=0) during the period 0<Ω _i t<4.0. The geometric average is defined as $\hat {B_{z}}$.

• 2. Subtract $\hat {B_{z}}$ from B _z for each time step (defined as b _z ) to remove the B _z component that originally existed in the DL region.

• 3. Take the FFT of b _z in the x direction [we define the transformed b _z as $\tilde {b_{z}}(m)$], and plot the power of $\tilde {b_{z}}(m)$, $\left (\left |\tilde {b_{z}}(m)\right |^{2}\right)$, for each m.

Figure 4 c shows the geometric average of B _z (0≤Ω _i t≤4) in Run(0) resulting from procedure (1). The reason why we did not use the initial condition, B _z (x,z=0,t=0), but have used the geometric average is that the configuration of B _z is changed by the compression of the CS during the early stages of the simulation because of the pressure difference. Figure 4 d indicates the development of $\left |\tilde {b_{z}}(m)\right |^{2}$ (m=2,3,4,5,6, and 8) in Run(0). $\left |\tilde {b_{z}}(m=6)\right |^{2}$ first develops in the early stages, and this mode corresponds to the wavelength seen in Fig. 4 a, kd ^′∼0.5. The growth rate of this mode is γ/Ω _i=0.37. The magnitude of this growth rate is about a half of the asymptotic electron tearing instability (Laval et al. 1966), $\gamma /\Omega _{\mathrm {i}} = \sqrt {\pi }\left (v_{\text {th}_{\mathrm {e}}}/d'\right)\left (\rho _{\mathrm {e}}/2d'\right)^{3/2}\left [(T_{\mathrm {i}}+ T_{\mathrm {e}})/T_{\mathrm {e}}\right ]\left (1-k^{2}d'^{2}\right)/\Omega _{\mathrm {i}} = 0.75$, when substituting the following parameters $\phantom {\dot {i}\!}v_{\text {th}_{\mathrm {e}}}=0.154c$, $\phantom {\dot {i}\!}v_{{\text {th}}_{\mathrm {i}}}=0.0326c$, and d ^′=1.6c/ω _pe, which are at the time of the onset Ω _i t=6.8. This growth rate compared to the theory is also consistent with the result of “case I” in previous research by Bessho and Bhattacharjee (2014). In their simulations, the electron tearing mode showed a smaller growth rate than that of the theory by a factor of 3. We investigated the growth rates of several other tearing theories including the “ion tearing instability” (γ/Ω _i=4.03) by Schindler (1974), “single particle tearing instability” (γ/Ω _i=2.93) by Pritchett et al. (1991), and “anisotropic tearing instability” (γ/Ω _i=2.04) by Quest et al. (2010), but the above theoretical growth rate is closest to our case. Therefore, this instability is likely to be electron tearing instability.

The small node of the tearing mode at x/d∼−9.4 is absorbed by interactions with other nodes soon after Ω _i t=8.4. On the other hand, the node at x/d∼−5.8 develops into a dominant magnetic reconnection in the simulation. The earthward flow generated by the reconnection conveys positive B _z to the DL region. This corresponds to an earthward propagating dipolarization front (Birn et al. 2011; Runov et al. 2009). The tailward flow conveying negative B _z couples with another small node of the tearing mode at x/d∼−13, which appears at that location because of the mode’s periodicity, but that node cannot be seen because Fig. 4 a, b does not cover that region. The dominant mode changes in correspondence with the change of the magnetic-field structure. Because a node of the tearing mode develops to a dominant magnetic reconnection and other small nodes disappear as the simulation progresses, the dominant mode number decreases from m=6 to m=4 in Fig. 4 d. The dominant mode further decreases to m=2 during the final stage of this simulation through nonlinear processes such as the coalescence of magnetic islands.

It should be noted that there are impacts of the characteristic initial conditions on the onset of the tearing instability. The normal magnetic field B _z1 given by non-self-consistent vector potential A _y1 decreases from 0.04 (x/d∼3) to 0 toward the tailward end. The J × B force resulting from incomplete force balance due to the the vector potential conveys the B _z earthward. The nodes of the tearing mode, however, appear at the tailward region where the tearing growth condition kρ _e>1 is satisfied, so the above J × B force has less of an effect on the tearing mode in this case. By performing a three-dimensional (3D) particle-in-cell (PIC) simulation, Pritchett (2005) demonstrated that the convective electric field locally induced as a boundary condition decrements the normal field as well as the thickness of the CS locally, and once the location satisfies the growth condition kρ _e>1, tearing instability appears and develops into magnetic reconnection. On the other hand, CS thinning in our simulation does not change B _z in the CS because the magnetic flux density is not changed in the tail axis by the compression. Liu et al. (2014) conducted similar, but 2D, particle simulations to those of Pritchett (2005) and indicated that an electron tearing instability easily occurs when the anisotropy of electron T _e⊥/T _e∥ is over unity because of CS thinning. In that context, Run(0) shows a similar feature to that of Liu et al. (2014) because the electron anisotropy is T _e⊥/T _e∥=1.07 at Ω _i t=6.8, but the wavelength in Run(0) is different from their result. Therefore, the CS thinning makes it easier for the instability to arise when the electron anisotropy increases, and it induces periodic nodes similar to those seen in Bessho and Bhattacharjee (2014).

Effects of the local B _z enhancement

Run(1) to Run(3) show different time developments depending on the initial locations of the local B _z enhancement seeded in the CS. Figure 5 shows B _z (z/d=0) at Ω _i t=8.4 and Ω _i t=15.2, as well as the simulated magnetic-field structure and J _y at Ω _i t=15.2 of each run. The panels showing B _z (z/d=0) also include the B _z component in Run(0) as a dotted line for comparison purposes.

Fig. 5

Simulation results of Run(1)–(3). B _z (z/d=0) at a Ω _i t=8.4 and b Ω _i t=15.2 in Run(1)–(3). B _z profile in Run(0) is shown for comparison (dotted line). c Same as Fig. 4 b at Ω _i t=15.2 in Run(1)–(3)

The local enhancement in Run(1) cannot induce a node around its position because there is a strong stabilizing effect caused by the normal magnetic field around the location. The nodes of the tearing mode appear at almost the same locations (x/d∼−5.6,−10.1) as those of Run(0). The node of the tearing mode at x/d∼−5.6 evolves into a dominant magnetic reconnection and a dipolarization-front structure is also generated from that location [see Fig. 5 b of Run(1)]. The B _z component and the magnetic-field structure are similar to those in Run(0), except that the start of the tearing onset is slightly later than that in Run(0).

In Run(2), the enhancement generates a node of B _z on its tailward side, and we found x/d∼−4.8 at Ω _i t=8.4. The notion that the node (i.e., the reconnection X-point) is located on the tailward side of the local B _z enhancement is consistent with observations in which the NENL is formed on the tailward side of the local B _z enhancement (Machida et al. 2009, 2014). The reason for the appearance of the node can be considered to be the convection of the enhancement earthward by J × B force due to the non-self-consistency of A _y2. Similar relaxation processes were seen in Sitnov and Swisdak (2011) and Bessho and Bhattacharjee (2014), but they adopted a relatively large-scale B _z -hump in the CS. Even if the growth condition kρ _e>1 is not satisfied, an earthward moving B _z -hump generates the magnetic reconnection in their case. However, in our local B _z enhancement case, kρ _e<1 suppresses the instability as seen in Run(1). The local B _z enhancement is originally located in the region where the growth condition is satisfied in Run(2) and the instability has appeared. This enhancement can change the wavelength of the nodes of the tearing mode in addition to their locations. Three nodes can be found in Fig. 5 a of Run(2) at x/d∼−4.8,−7.7, and −10.9; the corresponding wavelength is λ/d∼3 (kd ^′∼0.7). We speculate that the a J _y -free region near the DCB at x/d∼−2 acts as an effective boundary for the tearing mode, and the wavelength (kd ^′∼0.7) is defined by the distance between the effective boundary (x/d∼−2) and the directly induced node (x/d∼−4.8). Interestingly, the B _z profile of Run(2) at Ω _i t=15.2 shows that the node at x/d∼−7.7 (but not that at x/d∼−4.8) becomes a dominant magnetic reconnection. This result reflects the fact that the growth of the instability at x/d∼−4.8 is slightly suppressed by the finite normal magnetic field, while the growth of the node at x/d∼−7.7 is not suppressed because B _z ∼0. The growth of B _z generated by the reconnection is smaller than that seen in Run(0) at the same time, since the interaction with the node at x/d∼−4.8 interferes with the growth of the instability.

Just as for Run(2), the local B _z enhancement also generates a node on its tailward side in Run(3). B _z1 reaches almost to zero at the locus where the enhancement is seeded, and the suppression on the tearing instability is less around this region. Therefore, the enhancement quickly generates a node by the relaxation process, and the developing B _z attains a greater value than that of Run(0) at Ω _i t=8.4 and Ω _i t=15.2. In Fig. 5 a, b of Run(3), the earthward flow generated by the magnetic reconnection reaches a locus close to the DL region. Enhanced plasma pressure in the DL region produced by the earthward flow decelerates the flow, and the magnetic field in the flow piles up near the DL region at x/d∼−3 [see Fig. 5 c of Run(3)]. This pileup is caused by the two-dimensionality of this simulation because particles cannot escape into the y direction. This implies that the earthward flow with B _z (dipolarization front) generated by the near-Earth reconnection cannot produce the dipolarization near the DL region in 2D particle simulations.

Figures 6 a, b shows the growths of the modes in Run(2) and Run(3) in the same manner as shown for Run(0). The m=8 mode, which is greater than that of Run(0), becomes dominant at an early stage in Run(2). This mode corresponds to the shorter wavelength (kd ^′∼0.7) discussed previously. The dominant mode number decreases in a similar way to Run(0), but the dominant modes (m=5,3, and 2) are different from those seen in Run(0) because of the effects of the local B _z enhancement. In Run(3), it can be seen that the same mode numbers (m=6,4, and 2) as those in Run(0) become dominant, but these modes grow earlier. This is because the magnetic reconnection is rapidly induced by the B _z enhancement in this run. Although not shown in this paper, the time development of the modes in Run(1) is very similar to the result obtained for Run(0) since B _z in Run(0) and Run(1) exhibits similar time variations. These results confirm that the local B _z enhancement can modify the mode of the tearing instability in this periodic-magnetosphere system.

Fig. 6

Growth rate analyses for Run(2) and Run(3). a and b show time developments of the power of the modes resolved for B _z in Run(2) and Run(3), the same as Fig. 4 d

Discussion and conclusions

The tearing growth condition, kρ _e>1, is rewritten as follows:

$$ B_{z}/B_{0}\le\frac{\rho_{\mathrm{i}}}{2d}\sqrt{\frac{T_{\mathrm{e}}}{T_{\mathrm{i}}}\cdot\frac{m_{\mathrm{e}}}{m_{\mathrm{i}}}}, $$

((11))

given that kd ^′=0.5 and ρ _i/d ^′=1, where ρ _i is the ion-gyro radius for B ₀ (Pritchett 1994). If realistic physical quantities in the magnetotail (ρ _i∼d, T _i/T _e=5, and m _i/m _e=1836) are substituted into the right-hand side of the following inequality, the condition for instability becomes B _z /B ₀≤0.005. In this classical view, satisfying this condition seems to be difficult, even if the CS becomes very thin before an expansion phase. Recent studies, however, have suggested that the thinning process of the CS caused by the solar wind can lead the CS to an unstable condition. Initial 2D equilibrium of an isotropic CS develops into the kinetic one-dimensional form of the anisotropic CS via the growing solar wind pressure, and the CS may become unstable for the tearing mode (Zelenyi et al. 2009). Anisotropy of electrons generated by CS thinning also facilitates the tearing mode (Liu et al. 2014; Quest et al. 2010). Therefore, we think that it is important to investigate the basic physical features of the tearing mode in the magnetosphere when discussing substorm-triggering mechanisms.

This paper focuses on the instabilities in a magnetosphere -like system spanning the regions from the DL region to the CS by employing 2.5-dimensional full-particle simulations. Based on a unique method, we have adopted a new initial magnetic-field condition for the simulations. Instabilities in the CS are rapidly induced because of the low plasma pressure. The reconnection electric field “ E _y ” in the equatorial plane is expressed by using the electron-fluid equation of motion as follows:

$$ E_{y} = -v_{x\mathrm{e}} B_{z} - \frac{1}{en_{\mathrm{e}}}\left(\frac{\partial P_{xy\mathrm{e}}}{\partial x} + \frac{\partial P_{yz\mathrm{e}}}{\partial z}\right) - \frac{m_{\mathrm{e}}}{e}\frac{dv_{y\mathrm{e}}}{dt}. $$

((12))

In Fig. 7, the first term (Lorentz term), second term (electron pressure term), and third term (electron inertial term) of the right-hand side at Ω _i t=7.2 and 7.6 in Run(0) are indicated by a green line, blue line, and red line, respectively. Although there are some fluctuations, the electron pressure term gives two positive peaks at x/d∼−6,−9, where the nodes of the tearing mode appear. This means that the electron dynamics is crucial to the energy dissipation in the early time of the instability because the electron pressure term dominates in diffusion regions (Hesse et al. 1999; Kuznetsova et al. 2001). Therefore, we conclude that this instability is an “electron tearing instability,” even though the initial plasma sheet is not in an equilibrium state because of the low plasma pressure.

Fig. 7

Terms of E _y given by the electron-fluid equation of motion in Run(0). a and b indicate terms of the right-hand side of Eq. (12) at Ω _i t=7.2 and 7.6 in the equatorial plane of Run(0). Green lines are the first term (Lorentz term), blue lines are the second term (electron pressure term), and red lines are the third term (electron inertial term)

The nodes of the electron tearing mode appear in Run(0) where the growth condition is satisfied. The wavelength and growth rate of the mode show similar features to those of Bessho and Bhattacharjee (2014). In the real magnetotail, this wavelength corresponds to λ∼3Re, given that the thinned CS half-thickness is d ^′∼1500 km during substorms (Asano et al. 2004). This wavelength is shorter than the average distance between the DCB and NENL by a factor of 3–4, but it is consistent with the observation of the magnetic flux rope produced by a few Re spatial scales of multiple X-lines (Eastwood et al. 2005). If it is assumed that the J _y -free region near the DCB acts as an effective boundary for the tearing mode, such multiple X-lines may appear at the specific locations with each wavelength corresponding to the condition of the magnetosphere.

Three additional simulations with the local B _z enhancement included in the CS region, as observed by Machida et al. (2009, 2014), were performed. We found that in Run(1), where the enhancement is located near the DCB, nodes of the tearing instability appeared at similar loci as in Run(0) because the normal magnetic-field component of B _z1 suppresses the instability at that site. However, the local enhancement located tailward from the DCB induces a node of B _z on its tailward edge in Run(2) and Run(3) because of the relaxation process caused by the local enhancement. This result is similar to the findings of previous research (Bessho and Bhattacharjee 2014; Sitnov and Swisdak 2011). This also can explain the observations showing magnetic reconnection on the tailward side of the B _z enhancement at the substorm onset (Machida et al. 2009, 2014). An interesting result is obtained in regard to the location of the main node of the tearing instability, which develops into the X-point of a dominant magnetic reconnection. The node, which is preferentially induced tailward of the local B _z enhancement, develops into a dominant reconnection X-point in Run(3). In Run(2), however, the induced node does not develop into a dominant reconnection, but the node located one wavelength down the tailward direction becomes a dominant X-point. This is because a greater normal magnetic field B _z1 at the earthward node than that at the tailward region impairs the growth of the instability more effectively. This result implies that, even if such a local B _z enhancement induces a node of the tearing mode during the substorm period, the nodes at one (or several) wavelength(s) tailward of the node may develop into a magnetic reconnection depending on the locus of the enhancement with respect to the background normal B _z component.

We found that the enhancement changes the wavelength and mode number particularly in Run(2) by coupling between the relaxation process and the tearing n. Rapid appearance of the instability by the CS thinning from the initial pressure difference caused such coupling. Even if such a local enhancement exists, the change of the wavelength of the tearing mode may not occur in the real magnetotail because it is possibly caused by the boundary conditions in the x direction; hence, this result requires a reassessment of the boundary conditions. That is an aspect of our future work. We have also shown that a dominant mode of the tearing instability switches because of the influence of nonlinear effects, such as the coalescence of magnetic islands, as the simulations progress.

In this study, we had to deal with several free parameters that define the form of the magnetospheric system. We note that p and ε ₂ are particularly important. As p, which defines the roundness of the DL region, decreases, the shape of the region becomes closer to that of the Harris sheet and the J _y -free region is not formed. A smaller value of ε ₂ than what was adopted may not preferentially induce nodes of the tearing instability. Nevertheless, the ratios determined by these two parameters, B _z0:B ₀=1:2 (B _z0; the maximum of B _z in the DL region, B ₀; B _x in the lobe region) and B _z0:B _z2=10:1 (B _z2; the maximum of B _z2), do not deviate significantly from the observations (Machida et al. 2009, 2014).

The local B _z enhancement has been included artificially in this study, but why and how such an enhancement exists at X∼−17 Re prior to substorm onset is still an open question. One possibility is that the enhancement is produced by an interaction between the thinning process of the CS and the earthward flows with a small B _z component that propagate from a distant neutral line. More extensive observations or sophisticated simulations are necessary to further investigate this problem. The physical mechanism related to dipolarization has not occurred in our simulations, but CD is also statistically expected to occur simultaneously with magnetic reconnection around the DCB (Machida et al. 2009, 2014; Miyashita et al. 2009). The DL region’s shape and size remain almost constant in the simulations, despite the development of a tearing instability in the CS. This implies that because the DL region acts as an effective boundary for the tearing mode, the tearing mode cannot directly affect the physical processes in the DL region. Therefore, CD is considered to be the consequence of other types of instabilities. One of the limitations is that in a 2D system, any possible couplings between a tearing mode and other instabilities cannot be modeled properly. In the future, we plan to use full 3D simulations to investigate whether such couplings occur or not.