Multiscale Nature of the Magnetotail Reconnection Onset

Mining of substorm magnetic field data reveals the formation of two X‐lines preceded by the flux accumulation at the tailward end of a thin current sheet (TCS). Three‐dimensional particle‐in‐cell simulations guided by these pre‐onset reconnection features are performed, taking also into account weak external driving, negative charging of TCS and domination of electrons as current carriers. Simulations reveal an interesting multiscale picture. On the global scale, they show the formation of two X‐lines, with stronger magnetic field variations and inhomogeneous electric fields found closer to Earth. The X‐line appearance is preceded by the formation of two diverging electron outflow regions embedded into a single diverging ion outflow pattern and transforming into faster electron‐scale reconnection jets after the onset. Distributions of the agyrotropy parameters suggest that reconnection is provided by ion and then electron demagnetization. The bulk flow and agyrotropy distributions are consistent with MMS observations.

modifying the initial tail current sheet equilibrium based on isotropic plasma models (Harris, 1962;Schindler, 1972): We change the ratio between ion and electron bulk flow velocities and embed such a TCS into a tenuous plasma background to attenuate the electrostatic field toward TCS well outside it. The charged TCS theory (Yoon & Lui, 2004b) is used to adjust both the initial plasma equilibrium and the boundary conditions. We also take into account a weak driving electric field, which appears because of the interaction SITNOV ET AL.   Sitnov, Stephens, et al., 2019) and the equatorial magnetic field distribution with gray dots showing the projections of the spacecraft coordinates for nearest neighbors selected at this moment (2:30 UT). (c) The total current density j y at the moment ω 0i t = 5.8 with overplotted magnetic field lines (isocontours of the electromagnetic vector potential A y ). (d) The electric field component E z , which reflects the thin current sheet (TCS) negative charging effect. (e)-(h) Distributions of the electron and ion current densities j ey and j iy as well as bulk flow velocities v ey and v iy (normalized by the effective Alfvén speed  n m , where B 0 is the magnetic field outside CS at x = 0) that reveal the electron-dominated (ED) TCS effect. The electric field is normalized by B 0 v A /c (where c is the speed of light), while the current densities by q 0 n 0 v A (where q 0 is the proton charge). of the magnetosphere with the solar wind (e.g., Nishimura & Lyons, 2016) and which is critically important for onsets provided by the electron demagnetization.
The new DM-guided PIC simulations are aimed to clarify roles of electron and ion demagnetization and pre-onset TCS features in the onset mechanism. To quantify demagnetization, we use the Q-parameter (Swisdak, 2016), which reflects agyrotropy of ion and electron species and hence deviations of their orbits from conventional magnetic drifts (e.g., Pritchett et al., 1991). To clarify roles of the finite plasma background, negative charging and electron current domination, the base run (hereafter Run 1) is compared with a benchmark run (Run 2) with the same background but without bulk flow velocity shifts. The simulation results are compared with post-onset DM reconstructions of substorms, and local observations of the magnetotail reconnection, including MMS data.

Simulation Setup Constrained by TCS Features and DM Results
The major problem in constraining first-principle simulations of the magnetotail reconnection on the global scale is the extreme sparsity of in-situ observations. Recent DM studies (Sitnov et al., 2008;Stephens et al., 2019) mitigated this problem by enriching the small number of probes available at the moment of interest by observations made at other moments when the magnetosphere was in similar storm and substorm states ("nearest neighbors" or NNs is the southward component of the Interplanetary Magnetic Field). Swarms of synthetic NN probes, large enough for global reconstructions, and at the same time, small enough to reflect the state of the magnetosphere, then can be used. An example of such a swarm of probes for K NN = 32,000 NNs selected in the database of ∼4.10 6 historical records of the geomagnetic field is shown in Figure 1b for the growth phase of the February 13, 2008 substorm .
Large K NN numbers allow one to employ in the magnetic field reconstruction a sophisticated description of the magnetic field in the form of a system of basic functions Tsyganenko & Sitnov, 2007). Moreover, it becomes possible to distinguish between the evolution of the ion-scale TCS and a halo of the thicker CS with the thickness ≳1R E , where R E is the Earth's radius ( Figure 1a). With its 5 min time cadence, the DM reconstruction is found to reproduce the tail current sheet thinning and dipolarization during substorms , as well as the formation of new X-lines at the substorm onset . It also reveals that some onsets are indeed preceded by the formation of B z humps ( Figure 1b).
Being guided by this empirical pre-onset picture, we start simulations from a two-dimensional (2D) equilibrium with a B z hump (Sitnov & Schindler, 2010). Its magnetic structure is specified in the supporting information Material (hereafter SIM). As one can see from Figure 1c, it is consistent with the DM reconstructions (Figures 1a and 1b). The selected values of the CS thickness L = 1d i , where d i = c/ω pi is the ion inertial length (   2 0 4 / pi i e n m is the plasma frequency; n 0 is the plasma number density at the earthward side of the simulation box near the neutral plane z = 0), also matches the DM picture with L TCS ≈ 0.2R E .
Tail CS plasmas are usually described by isotropic particle distributions with shifted Maxwellians (e.g., Lembege & Pellat, 1982)  where the canonical momentum P yα = m α v y + (q α /c)A y and the total energy W α = m α v 2 /2 + q α ϕ (ϕ is the electrostatic potential) are the integrals of motion making any function f 0α (P yα , W α ) to satisfy Vlasov equation.
In many CS models, it is additionally assumed that electron and ion drift velocities v Dα obey the condi- where T i and T e are ion and electron temperatures, to maintain charge neutrality ϕ = 0. However, the CS thinning process (Lu et al., 2016;Pritchett & Coroniti, 2013) ). An additional background plasma population with the density n b /n 0 = δ b ≪ 1, where n 0 is the CS density maximum, eliminates the invariance of the transformation and it also shields the electrostatic field E z outside CS. As follows from (Yoon & Lui, 2004b), in the limit   0 Di v and β(x) = 1 this field can be presented in the form The three-dimensional (3D) PIC simulations were performed using the explicit massively parallel code P3D (Zeiler et al., 2002) in a 3D box with dimensions L x × L y × L z = 80d i × 5d i × 20d i that were optimized to consider a sufficiently long CS, comparable to the corresponding tail region in the DM reconstructions (80d i correspond to ∼8R E assuming that d i ∼ 0.1R E ). The box extension in the dawn-dusk direction was chosen to accommodate at least one wavelength of the largest (flapping) perturbations detected in earlier simulations (Sitnov et al., 2014). Other simulation parameters are provided in the SIM.
To allow free propagation of reconnection motions along the tail we employed a modification of open boundary conditions (Divin et al., 2007;Sitnov & Swisdak, 2011), including free reconnection electric field component at the boundary ∂E y /∂x = 0 and taking explicitly into account the electrostatic field at the left boundary x = 0: es z E being determined by 2. The newly injected equilibrium plasma at the boundary to maintain the finite pressure gradient ∂p/∂x = j y B z has the same velocity shift as the ED distributions inside the box. At the right boundary x = -L x where the CS density is small the electrostatic field is neglected. The boundary conditions in y and z directions are periodic and conducting (Sitnov & Swisdak, 2011;Sitnov et al., 2014). To mimic the solar wind driving of the tail, a weak external electric field ( ) dr y E is applied at top and bottom boundaries as is detailed in the SIM.
Figures 1c-1h provide key details of the resulting selfconsistent plasma configuration at ω i0 t = 5.8, well before the major instabilities develop. Figure 1d shows the established electrostatic field, consistent with the equilibrium theory (Yoon & Lui, 2004b), and in particular, the analytical profile (2). Figures 1e and 1f show that the current is dominated by electrons. They can be contrasted with the corresponding current distributions for the benchmark Run 2 provided in Figure S1 of the SIM and describing an ion-dominated CS. Distributions of the y component of the bulk flow velocities presented in Figures 1g and 1h are quite different from the idealized ED CS with v Di = 0 and v De = −1 that would be expected for zero background. However, the reduction of the electron bulk flow velocity outside CS is an expected effect of the background population, which is particularly well seen in Figure S1e. The increase of the (initially zero) ion velocity v Di seen in Figure 1h occurs outside TCS where the current density is small (Figures 1c and 1f).

Multiscale Reconnection Onset
The state of the tail CS later in this run (ω i0 t = 55.8) is described in Figure 2, close to the formation of two X-lines at x ≈ −25d i and −65d i . It shows the development of a tearing mode with the wavelength λ PIC ∼ 40d i , comparable to the global inhomogeneity scale in X ∼ d i /0.03. Its growth becomes energetically favorable in spite the electron compressibility effect (Lembege & Pellat, 1982) due the B z hump (Sitnov & Schindler, 2010). However, earlier simulations of the B z -hump effect usually revealed the formation of a single X-line in the trail of the dipolarization front (DF) that formed from the hump (Bessho & Bhattacharjee, 2014;Pritchett, 2015;Sitnov et al., 2013Sitnov et al., , 2017. A similar process is seen in Figure S2, an analog of Figure 2 for the ion-dominated CS (see, in particular, Figure S2a). The second X-line, which forms farther in the tail in Runs 1 and 2, is provided by the external driving. It goes through electron demagnetization, and further, steady reconnection, consistent with the earlier models (Hesse & Schindler, 2001;Liu et al., 2014). In particular, Sitnov et al. (2021) in a run similar to the present ones, but with the zero background plasma and without the ED effect (referred to below as Run 0), showed that the reconnection electric field in that "midtail" region is homogeneous and its rate is close to the theoretical estimates E y ∼ 0.1 (Cassak et al., 2017).
The formation of two X-lines (Figure 2a) matches the global DM-based reconnection picture (Sitnov et al., 2021) suggesting that simulations properly describe this type of reconnection onsets. The distance SITNOV ET AL.

10.1029/2021GL093065
between X-lines in the DM picture λ DM ≳ 10R E (Sitnov et al., 2021, Figure 2 and 8) matches λ PIC . Consistent with the DM results, stronger B z perturbations are found closer to Earth (x ≳ −40d i ) suggesting that the reconnection process there is more unstable and hence the electric field is more inhomogeneous (∂E y /∂x = −∂B z /∂t). The latter is seen from Figure 2i where the E y peak exceeding the classical 0.1 value is SITNOV ET AL.  The dented structure of the B z hump seen in Figure 2a for x > −20d i is likely caused by two-stream (e.g., Yoon & Lui, 2004a) and lower hybrid drift (LHD) (Divin et al., 2015;Yoon & Lui, 2004a) instabilities. The former is caused by the velocity shift of the CS electron population relative to the background, which is small compared to the electron thermal speed    v v T T T m m but significant compared to the ion thermal speed. The latter is present already in the ion-dominated CS in case of the finite background, as is seen from Figures S2a-S2d. Note that these 3D effects in 2D ED TCS are different from the 2D evolution of 1D ED TCS (e.g., Lu et al., 2020, Figure 5), although in both pictures the electron and ion currents become comparable in density as the corresponding plasmas evolve in time ( Figure S3b). Figure 2b shows that the TCS remains negatively charged, but its charging is rather provided by the external driving, because similar E z patterns are found in Runs 0 and 2 with ion-dominated CSs ( Figure S2b and Figure 13b in Sitnov et al., 2021). At the same time, E z patterns in Runs 1 and 2 show more ion-scale structuring, compared to Run 0, including striations, different from LHD effects (e.g., Figure 2c and Divin et al., 2015).  (Phan et al., 2016;Torbert et al., 2018). While ions gain some bulk velocity (Figure 2e), the CS remains ED over the whole run, as is seen from Figure S3. The temperature distributions show that the TCS heats up (the original temperatures T e = 0.125 and T i = 0.375 in the code units) with strong gradients across the sheet matching local observations (Artemyev et al., 2017;Lu et al., 2020). Figures 2h and 2j show development of B/I and flapping motions (FM) in the reconnecting TCS. Their wavelengths are consistent with (Pritchett & Coroniti, 2013;Sitnov et al., 2014). The main distinctions, compared with the zero background case (Run 0, Figures 13a, 13e, and 13g in Sitnov et al., 2021) are smaller amplitudes and longer growth times of FM, B/I and reconnection motions. For the latter two this is consistent with the stability theory (Merkin & Sitnov, 2016;Sitnov, Birn, et al., 2019). In addition, the comparison of Figures 2j and S2j reveals more elongated shape of FM striations for ED TCS, which is consistent with the original idea of obtaining ED TCS by the coordinate transformation to the system moving with ions (Yoon & Lui, 2004b).

Electron and Ion Plasma Watersheds and Agyrotropy Picture
The most impressive new features found in this study are the distributions of the bulk flow velocities and agyrotropy parameters for electrons and ions shown in Figure 3. They reveal two patterns of diverging electron flows with significant velocities (≳v A ) marked by dashed green lines with two arrowheads in Figure 3a. We call them "watersheds" (WSs) to distinguish from flows diverging from an X-line. They form within a single ion WS (Figures 3d and 3e) in the regions where the new X-lines form later (Figures 3c and 3f). As is seen from the comparison of Figures 3a and 3b, WSs form at least 10 ion gyrotimes before the onset. They break the paradigm of the topology change as a prerequisite of the reconnection outflow formation. Such a counter-intuitive dynamics is nevertheless quite consistent with the tearing theory: It suggests that the reconnection free energy comes from the mutual attraction of the parallel current filaments in the CS (Forslund, 1969;Galeev, 1984) rather than from any topology changes (see also Bellan (2006) for a water-beading analogy). Hence, it may occur when electrons are yet magnetized by the field B z (Schindler, 1974). The demagnetization of electrons and ions shown in Figures 3g-3l in terms of the agyrotropy parameter Q (Swisdak, 2016) confirms this interpretation: Electrons reveal significant agyrotropy near the EDR only after the topology change (Figure 3i) whereas ions are demagnetized well before that moment and quite away from any X-lines or their expected locations (Figures 3j-3l). As is seen from the comparison of Figures 3a and 3d, electron WSs more closely follow the perturbed magnetic field than the ion WS, which is consistent with more inertial response of less magnetized ions to the field perturbations. After the formation of new X-lines, every electron WS (Figure 3c) is accompanied by the ion WS ( Figure 3f) and it forms a conventional electron-scale reconnection jet pattern (Figure 3c) around the EDR (Figure 3i).  Figure 4b). The comparison of the Q-parameters (Figures 4d and 4h) shows that electrons in observations are even more magnetized during their embedded reconnection event. This can be explained by both the artificially small mass ratio in simulations (m i /m e = 128) and by the significant guide field effects as discussed by Chen et al. (2019).
Multiple electron WSs are consistent with the electron tearing onset models (Liu et al., 2014, Figure 5), although the predicted wavelength λ ∼ d i may be too small. Besides, both our simulations and MMS observations predict that some WSs (all but one in observations shown in Figure 4g) are embedded into a single earthward ion flow. This is more consistent with the ion tearing driven by DF acceleration (Sitnov et al., 2017, Figure 1). Embedded WSs also resemble the so-called "electron-only" reconnection observed in the turbulent magnetosheath . However, in contrast to the latter, in our simulations the embedded electron WSs appear well before the formation of the corresponding X-lines (Figures 3a-3b), whereas in MMS observations of the tail, the corresponding B z variations inside the earthward ion flow resemble more DFs than X-lines because of very small negative B z regions (Figure 4e).

Discussion and Conclusion
In this study, we combined for the first time in 3D PIC simulations multiple features of the pre-onset magnetotail: Magnetic flux accumulation, ion-scale TCS formation earthward of it, solar wind driving and negative charging of the current sheet with electrons being the dominant current carriers. The first feature inferred from data-mining reconstructions agrees with earlier statistical studies (Machida et al., 2009;Wang et al., 2004) and recent remote sensing analysis (Sergeev et al., 2018). The last feature (domination of electrons as current carriers) is consistent with local observations, but the corresponding equilibrium (Yoon & Lui, 2004b) is substantially modified by two-stream instabilities. The obtained multiscale onset picture involves both electron and ion mechanisms with electron WSs being embedded into a single ion watershed, and it agrees with MMS observations. Simulations also match the global post-onset reconnection picture SITNOV ET AL.  offered by data-mining reconstructions. Thus, they reproduce the inherently multiscale nature of magnetotail reconnection.
The idea of WSs, that is, diverging plasma flows internally driving the reconnection onset in the magnetotail has been formulated earlier using global hybrid (Lin & Swift, 2002) and MHD (Siscoe et al., 2009;Tanaka et al., 2019) simulations. Triggering reconnection by non-reconnection instabilities is also popular in the solar physics (e.g., Török & Kliem, 2005), but such instabilities (e.g., kink) are different from tearing. In the magnetotail, the ion tearing instability does not change the magnetic field topology either, as electrons must be magnetized by the field B z . It may even start at MHD scales (Birn et al., 2018;Merkin et al., 2015). But it naturally triggers reconnection due to ion and electron WSs. Thus, its generalizations for more complex plasma configurations, such as in the solar corona, can provide direct link from global MHD-scale motions to microscale reconnection physics. The main advance of the present study is the fully kinetic description of the plasma watershed hierarchy from MHD to kinetic scales and its evolution from the observation-guided global tail configurations to the topology change and the quantitative description of electron and ion demagnetization.

Data Availability Statement
The data used in the study are archived on Zenodo (http://doi.org/10.5281/zenodo.4553406). 80NSSC20K1271, 80NSSC20K1787, and 80NSSC19K0396, as well as NSF Grant AGS-1744269. Simulations were made possible by the NCAR's computational and information systems laboratory (http://doi.org/10.5065/D6RX99HX), supported by the NSF, as well as the NASA high-end computing program through the NASA advanced supercomputing division at Ames Research Center.