Data Mining Reconstruction of Magnetotail Reconnection and Implications for Its First-Principle Modeling

Magnetic reconnection is a fundamental process providing topological changes of the magnetic field, reconfiguration of space plasmas and release of energy in key space weather phenomena, solar flares, coronal mass ejections and magnetospheric substorms. Its multiscale nature is difficult to study in observations because of their sparsity. Here we show how the lazy learning method, known as K nearest neighbors, helps mine data in historical space magnetometer records to provide empirical reconstructions of reconnection in the Earth’s magnetotail where the energy of solar wind-magnetosphere interaction is stored and released during substorms. Data mining reveals two reconnection regions (X-lines) with different properties. In the mid tail ( ∼ 30 R E from Earth, where R E is the Earth’s radius) reconnection is steady, whereas closer to Earth ( ∼ 20 R E ) it is transient. It is found that a similar combination of the steady and transient reconnection processes can be reproduced in kinetic particle-in-cell simulations of the magnetotail current sheet.


INTRODUCTION
Charged particles, electrons and ions forming space plasmas usually drift in the ambient magnetic field making plasmas frozen in that field [1]. The frozen-in condition may be broken when oppositely directed field lines approach each other so closely that particles become unmagnetized and their orbits become different from conventional drift motions. As a result, magnetic field lines may change their connectivity near so-called X-lines in the process of magnetic reconnection. This process was introduced to explain major sources of space weather disturbances on the Sun, solar flares [2,3] and coronal mass ejections (CMEs) [4]. It was also invoked by Dungey [5] to describe the structure of the Earth's magnetosphere, the plasma bubble surrounding our planet and protecting its life from the hazardous stream of high-energy particles emitted by our star. According to Dungey, reconnection takes place on the day side of the magnetospheric boundary, the magnetopause, to provide the solar wind plasma entry into the magnetosphere through the reconnected magnetic flux tubes. Then the flux tubes reconnect again on the night side, in the region where the Earth's dipole magnetic field lines are stretched in the antisunward direction forming the magnetotail. Finally, to explain a delayed explosive response of the polar regions of the magnetosphere to solar wind disturbances during substorms [6], Hones [7] proposed that the substorm explosions are powered by the unsteady reconnection in the tail due to the formation of another "near-Earth" X-line.
The magnetotail reconnection is important not only as a key element of the space weather chain. It occurs in space plasma practically in the absence of particle collisions. Similar collisionless reconnection processes are expected to occur in the solar corona during flares and CMEs, where in-situ observations are impossible [1]. They are also expected in sufficiently hot laboratory plasmas that are investigated on the way to controlled nuclear fusion [8]. Thus, the magnetotail represents a natural space laboratory for collisionless reconnection due to many dedicated missions, such as Geotail [9], Cluster [10], THEMIS [11] and MMS [12].
The magnetotail is also very interesting because it reveals different regimes of reconnection. On the one hand, it must experience steady reconnection, which was conjectured by Dungey [5] in his description of the magnetospheric convection cycle and later confirmed in observations of steady magnetospheric convection (SMC) regimes [13]. On the other hand, the magnetotail experiences unsteady reconnection during substorms [7].
Both the first-principle modeling and the empirical reconstruction of magnetotail reconnection are very difficult to perform because of its multiscale nature. It links global reconfigurations of the nightside magnetosphere to kinetic processes on the scales of ion or even electron gyroradii that provide irreversibility for global reconfigurations. As a result, its kinetic particle-in-cell (PIC) simulations describing the full dynamics of electrons and ions (largely protons) and their self-consistent electromagnetic fields [14] are usually limited to the immediate X-line vicinity [15] and the moments after the X-line formation in global magnetohydrodynamic models [16], where the reconnection onset is provided due to numerical or ad hoc plasma resistivity. Moreover, it is very difficult to take into account that the magnetotail itself becomes multiscale prior to the reconnection onset. In-situ observations suggest that it may contain thin (ion-scale) current sheets (TCS) embedded into a thicker current sheet (CS) [17][18][19][20][21]. The latter may also be split in two current layers forming bifurcated TCSs [19,[22][23][24].
The major problem in the empirical reconstruction of magnetotail reconnection, common to all in-situ space observations, is the extreme sparsity of these observations with fewer than a dozen probes available at any moment. To solve this problem, it has recently been proposed to mine data in the multi-mission database covering many years of historical spaceborne magnetometer observations [25,26]. It was found that such a data-mining (DM) method resolves the formation of embedded TCSs in the growth phase of substorms and their decay after the substorm onset. It also resolves the formation of the near-Earth X-lines during substorms [27]. Here we show that the DM approach allows one to resolve the formation of two different X-lines in the magnetotail during substorms. Moreover, it becomes possible to quantitatively assess their steadiness. We also show that PIC simulations guided by the DM reconstruction of the magnetotail reproduce the formation of X-lines and reconnection regimes similar to those found in the DM analysis.

DATA MINING METHOD
In the DM approach, the geomagnetic field is reconstructed using not only a few points of spaceborne magnetometer measurements available at the moment of interest, but also a much larger number of other measurements made at the K NN ≫ 1 moments in the past. These moments called "the Magnetospheric coordinate system (GSM) coordinate system is used throughout this paper. Its origin is at the center of the Earth; the X-axis is directed toward the Sun; the y-axis is defined as the cross product of the GSM x-axis and the magnetic dipole axis, directed positive toward dusk; the z-axis is defined as the cross product of the x-and y-axes.) The large number of NNs is at the same time much smaller than the size of the database K DB ≫ K NN . This allows one to fit with the NN subset a complex empirical magnetic field model [26], and at the same time, to make the model reconstructions sufficiently flexible to reflect the characteristic variations of the magnetosphere during storms and substorms.
This approach resembles very much the "lazy-learning" pattern recognition technique known as the K-nearest neighbor (KNN) learning [28,29]. At the same time, our DM approach differs from conventional KNN regression methods, where both finding the NNs ("mining") and regressions (model fitting) are made in the same space. Here, as is illustrated in Figure 1, we first detect NNs as a (sub)set of K NN present and historical moments in similar phases of similar substorms. Their similarity ("neighborhood") is quantified by the closeness of the corresponding global magnetospheric activity parameters and their time derivatives to their values at the moment of interest ( Figure 1A). Then we use these K NN moments to form an eventoriented subset of the original database of magnetic field observations ( Figure 1B) and to fit our magnetic field model with this subset ( Figure 1C).
The NN subset is formed by points G (i) [G 1 (t i ), . . . , G 5 (t i )], i 1, . . . K NN , in the 5-D space that are closest to the query point (moment of interest t q ) G (q) [G 1 (t q ), . . . , G 5 (t q )] by the Euclidean metric where σ G k is the standard deviation of the component G k and the coordinates G 1 -G 5 are defined by the formulae: Here Sym − H * A · Sym − H − B · P dyn is the pressurecorrected Sym-H index [30], P dyn is the solar wind dynamic pressure (in nPa) and the values of A and B are taken to be 0.8 and 13.0, respectively. The functions G 1 and G 3 in Eqs. 2, 4 describe weighted moving averages of the indices Sym-H and AL limited to their past values (see [25] for further details), while G 2 and G 4 , defined by Eqs. 3, 5, describe the corresponding smoothed time derivatives. Weighting in moving averages (2)(3)(4)(5) is provided by the sine and cosine kernel functions and by the exponential function in Eq. 6. The averaging scaling parameters Π st 12 hr and Π sst 2 hr reflect the characteristic storm and substorm scales. The parameter G 5 defined by Eq. 6 describes the integral effect of the magnetic flux accumulation in the tail during the growth phase due to the dayside reconnection. Its scale τ 0 0.5 hr is selected based on observed values of a typical growth phase duration [31]. The selected upper integration limit in Eq. 6 τ ∞ 6τ 0 corresponds to six e-folding times.
In the 5-D space of the binning parameters (2)-(6), the AL index and its time derivative (G 3 , G 4 ) determine the strength and phase of the substorm activity, because the AL index reflects the strength of the substorm electrojet [32]. These parameters may still be insufficient to capture the substorm growth phase, which is characterized by the accumulation of the magnetic flux in the tail lobes without any significant electrojet enhancement. To take this effect into account, we involve in the analysis the solar wind electric field parameter through the binning variable G 5 . Furthermore, many substorms occur at the moments of the storm activity, which may substantially modify the substorm evolution of the magnetosphere [33]. To take these effects into account, we further extend the binning space at the expense of the parameters G 1 and G 2 reflecting the storm-time index Sym-H and its time derivative (to distinguish between main and recovery storm phases).
The conjecture that the substorm dynamics of the magnetosphere is coherent and hence the distribution of its magnetic field can be determined by a few control parameters had been formulated many years ago (e.g., [34], and refs. therein). The event-oriented subset of the database formed using the selected set of K NN nearest neighbors. It is used to fit the magnetic field model. Gray dots show the projections of the spacecraft coordinates on the equatorial plane, where the normal magnetic field component B z is color coded. (C) The resulting magnetic field lines (black) and the equatorial field B z (color-coded) assuming zero tilt angle (adapted from [27]).
Later, the singular spectrum analysis of substorms [35] revealed that the mean-field dynamics of the magnetosphere can be described as a motion on a folded 2-D surface in a 3-D state space formed by the average AL index, average vB IMF s parameter and its average time derivative. An increase of the dimensionality through the Sym-H index and its time derivative, to take magnetic storms into account and to distinguish between their main and recovery phases, is consistent with the original DM-based stormtime model, TS07D [36]. The latter was also justified by the empirical relationship between the vB IMF s parameter and the Dst index, a 1-h time resolution analog of Sym-H [37]. Independent description of storms and substorms assuming their common solar wind driver is consistent with recent analysis of the stormsubstorm relationship using a multivariate information-theoretic approach [38]. The further increase of state space dimensionality (e.g., using higher time derivatives of storm and substorm indices as well as the solar wind input parameter) is also possible (e.g., [25]). However, it was found [39] that the effect of higher dimensions often resembles the second-order phase transition fluctuations that require a probabilistic description of the magnetospheric states [40].
The database consists of K DB 3, 668, 101 records of the magnetic field vector with 5 and 15-min cadence inside and outside 5R E , respectively, in archived data from IMP-8, Geotail, Polar, GOES-08, GOES-09, GOES-10, GOES-12, Cluster, THEMIS, Van Allen Probes and MMS missions covering more than two decades (1995-2017) of observations. The KNN subsets are selected using AL, Sym-H and vB IMF s time series with 5-min cadence. At every moment t q , the subset is found as K NN points i 1, . . . , K NN satisfying the condition R (i) q is defined by Eq. 1, as is illustrated in Figure 1A. Since the resulting magnetic field geometry is determined by the instantaneous KNN swarm of virtual probes, its time resolution is largely determined by the global parameter cadence. This is seen, for instance, from rapid substorm dipolarizations reproduced by the KNN method in [27] ( Fig. 8i), when the B z field increases from 3 to 10 nT in 5 min over a significant part of the magnetotail. It was found [26,27] that the use of NN subsets with K NN ∼ 32, 000 and the magnetic field model parameters specified below provides both sufficient selectivity of the model, which allows one to distinguish different substorm phases, and the high spatial resolution to resolve the distinctive features of the magnetospheric morphology in these phases, such as TCSs (and their buildup and decay), flux accumulation regions and X-lines. Smaller K NN values were found to cause overfitting.
At every moment of interest t q , the K NN subset of the database, whose elements neighbor t q in the state and input space of the magnetosphere by the metric (1), is used to fit the geomagnetic field model SST19 [26]. Since K NN ≫ 1, its architecture can be made quite complex and flexible (compared, for instance, with the event-oriented models using only a few points of data available at the moment of interest [41,42]) to capture key features of the substorm current system. In fact, we only assume that the magnetic field is formed by two major current systems inside the magnetosphere, equatorial and field-aligned currents, whose contributions are presented as sums of basis functions with the corresponding amplitude coefficients (more general 3-D expansions of the magnetic field using radial basis functions were considered in [43]). Moreover, to describe the multi-scale structure of the equatorial currents, including the formation of embedded and bifurcated TCSs, these currents are described by two independent expansions: where (ρ, ϕ, z) are cylindrical coordinates in a system with the origin at the center of the Earth and the z axis normal to the equatorial plane. They represent the magnetic field of thick and thin current sheets with the same structure determined by the approximate solution for the magnetic field of an arbitrary distribution of equatorial currents [44] with different thickness parameters D and D TCS to be derived from the fitting with the NN subset. Each expansion is a finite-sum approximation of an integral solution of the Ampère's equation for the magnetic field of an infinitely thin CS (D 0) above and below the equatorial plane z 0 by separation of variables. Tsyganenko and Sitnov [44] showed that the sum consists of N azimuthally symmetric radial expansions and 2M angular Fourier harmonics (even and odd parity in ϕ) with the total number of N + 2M · N elements. The basis functions of the solution for the vector potential with an infinitely thin CS contain factors like exp(−k|z|). Their regularization comes from assuming the finite CS half-thickness D and it can be provided by replacing the function |z| by the smooth function ζ . The radial expansions include Bessel functions and they can be exemplified by the azimuthal component A ϕ of the vector potential corresponding to the azimuthally symmetric group of basis functions B (s) 0n : (A ϕ ) (s) 0n J 1 (k n ρ)exp(−k n z 2 + D 2 √ ), where J 1 is the Bessel function of the first order, k n n/ρ 0 , and ρ 0 is the radial scale, corresponding to the largest mode in the radial expansion.
The parameters ρ 0 , N and M are fixed because they determine the adopted resolution of the expansions in Eq. 7. Other parameters, such as the weights of individual radial and azimutal harmonics, as well as the CS thickness parameters D and D TCS , are determined from fitting the model to data. In particular, to distinguish between thick current sheets and TCS, we impose the condition D TCS < D (max) TCS 1R E . The latter value is intermediate between the observed thick and thin current sheet values [19] and it does not significantly constrain the specific values of D and D TCS inferred from data. Thus, the spatial resolution of such an expansion is determined by the number of terms in expansions (7) and can be increased to any desired level, commensurate with the data density. To take the global scaling of currents due to variations of the solar wind dynamic pressure P dyn into account, each amplitude coefficient in expansions (7) is further expanded in two parts, one of which is constant and another is a linear function of P dyn . The equatorial expansion has several other nonlinear parameters to take into account global deformations of the tail CS along the dawn-dusk direction and arising from the finite dipole tilt angle, which are described in [44].
Another major group of currents are the field-aligned currents (FACs), connecting the ionosphere with the magnetopause and the tail CS. It is described in SST19 using a similar system of finite current elements [45], sufficiently flexible to reproduce the spiral FAC structure at low latitudes [46] whose night-side part is likely associated with the Harang discontinuity [47]. Each element of the FAC system is described as the magnetic field of two deformed conical surfaces corresponding to Region 1 (R1) and Region 2 (R2) FACs [48]. The size of each system is an adjustable parameter, while their azimuthal distribution is controlled by the relative contributions of two groups of basis functions with odd and even symmetry due to factors sin(lϕ) and cos(lϕ), (l 1, 2, . . .). The first group represents the main part of the FAC system, in which the dusk-side currents have the same magnitude but opposite direction to those at dawn, while the second group has an even distribution of currents with respect to the noon-midnight meridian plane, which allows one to model the azimuthal rotation of the FACs.
Originally two groups of such FAC elements were proposed in [44] to describe R1 and R2 systems in their DM-based storm-time model, TS07D [36]. Later, it was proposed [45] to use more elements similar to the original TS07D FAC basis functions, shifted in latitude to describe more complex FAC distributions. Eventually, Stephens and coauthors [26] showed that the FAC system can be described with many details important for substorm reconstructions, including the Harang discontinuity and the substorm current wedge [49], with the following set of elements. It consists of N FAC 16 basis functions with the first two Fourier harmonics (l 1, 2) for R1 and R2, as well as their latitude-shifted clones. Each element in equatorial and FAC expansions is independently shielded (has its own subsystem of Chapman-Ferraro-type currents at the magnetopause).
Thus, the resulting DM algorithm, which links the SST19 model with KNN binning, represents a typical "gray box" model combining empirical algorithms with physics-based constraints [50]. As is shown in [26,27], it reproduces the multiscale CS thinning process with the formation of an ionscale TCS (D TCS ≪ 1R E ) inside a much thicker CS (D > 2R E ), which takes place in the substorm growth phase and causes stretching of the tail magnetic field lines in the antisunward direction. In particular, Figs. 11a-11c in [27] show the current distribution in the equatorial plane during the growth phase of the 13 February 2008 substorm discussed below. The corresponding current distributions in the meridional plane presented in their Figs. 12a-12c reveal the multiscale CS structure with an ion-scale TCS embedded into a thick CS halo. The peak TCS current density ∼8 nA/m 2 is consistent with in-situ Cluster observations (see, for example, Figs. 2-4 and 9 in [19]). Further quantitative analysis made in [27] showed that while the TCS thickness in DM reconstructions remains approximately constant D TCS ≈ 0.2R E (Fig. 8e), consistent with Cluster and THEMIS observations [19,20,51], their strength (measured as the TCS contribution to the total tail current) changes drastically with the substorm phase ( Fig. 12d). At the same time, the contribution of the TCS to the total tail current is relatively small (≈ 1/6). It is worth noting that, as is seen from the comparison of Figs. 10a and 11a in [27], the extended TCS forms earthward of the flux accumulation region (B z hump).
The DM SST19 algorithm also describes the magnetic field dipolarization in the expansion phase (see, for example, Fig. 4 in [26]) with the formation of a substorm current wedge seen as a curl of the difference between the expansion and growth phase magnetic field distributions ( [26], Fig. 10). The disappearance of TCS after the dipolarization can be seen, for instance, from the comparison of Fig. 12d in [27] with other panels in Fig. 12. It can also been seen from their Fig. 8f, where the relative strengths of thin and thick CSs are quantified by integrating the current density over the regions |z| < 1R E and |z| < 5R E .
The new DM reconstruction has a characteristic property of machine learning algorithms [29,52]: given more data in the database it may provide more details about the magnetospheric structure and evolution. In particular, as is shown in [27], with adding to the database first two years of the MMS mission data (2016-2017), it becomes possible for the same SST19 model with the parameters (M, N) (6, 8) and K NN 32, 000 to resolve more details of the magnetotail structure and evolution. In particular, the 2017 MMS data help resolve the X-lines forming largely beyond 20R E , where the pre-MMS database had a substantial drop in the occurrence rate distribution ( [27], Fig. 1). This is seen in particular, from the comparison of the SST19 validation using THEMIS data beyond 20R E in Fig, 2e of [26] with THEMIS validation in a similar region in Fig. S6 of [27]. The former reveals clear signatures of overfitting while the latter does not. In spite of a relatively small total number of the new MMS data, they fill the main gap in the existing database distribution ( [27], Fig. 1) and thus become particularly important in solving the overfitting problem. The latter analysis is extended here by increasing the maximum radial distance of the spacecraft data used in the reconstruction from 31R E to 35R E (largely, due to IMP8 data). Figure 2 shows the equatorial magnetic field distribution at the moment 02:40 UT in the expansion phase of this substorm. It reveals the formation of two X-lines X n and X m in the near-Earth (x ≈ − 20R E ) and midtail (x < − 27R E ) regions, respectively. They are seen as earthward parts of the B z 0 contours in the distribution of the equatorial north magnetic field component B z in Figure 2 and they are additionally marked by blue arrows. The tailward parts of the B z 0 contours represent O-lines. This global X-line reconstruction is quite unique. In fact, because of the extreme sparsity of in-situ space observations, such reconstructions were not available before the machine learning era. Earlier, Nagai et al. [53,54] described the location and the dawn-dusk extension of X-lines using singlepoint observations. More recently, the reconstructions of the X-line vicinity were made by processing multi-probe MMS data with Grad-Shafranov [55] and polynomial [56] techniques. However, these were still very local reconstruction, largely limited to the size of the MMS tetrahedron (<30 km). Here we demonstrate for the first time how the DM approach based on the KNN algorithm resolves simultaneously two X-lines in the near-Earth and midtail regions.
The formation of transient near-Earth X-lines, which was proposed by Hones [7] as a mechanism of substorms, has been discussed since that time in many studies, including correlated multi-probe and remote sensing analyses (see, for instance [57,58], and references therein). At the same time, persistent reconnection in the midtail around 30R E follows from THEMIS and ARTEMIS statistics of traveling compression regions [59,60]. However, neither the co-existence of the second, midtail X-line X m with X n nor its relatively steady reconnection, as suggested by Dungey's convection cycle [5], have ever been demonstrated. Here we not only resolve two X-lines in the tail but also propose a method to quantify their steadiness.
This can be done using the Faraday's law, which in the 2-D picture of reconnection takes the form It suggests that the temporal variations of B x and B z magnetic field components determine the spatial gradients of the dawndusk (reconnection) electric field. If the magnetic field varies slowly, the corresponding reconnection electric field is broadly distributed over the whole reconnection region. This justifies the concept of the reconnection rate, one of the key global parameters characterizing steady reconnection regimes [61][62][63]. Reconnection can also be unsteady with the electric field being localized in space and the magnetic field changing in time consistent with (8). For example, Sitnov and Swisdak [64] showed reconnection regimes with the electric field localized near dipolarization fronts (DFs) [65][66][67] with their values strongly exceeding the steady reconnection values. Localization of the dawn-dusk component of the electric field near DFs was later confirmed by Cluster [68] and THEMIS [69] observations. The noon-midnight meridional maps of magnetic field lines presented in Figures 3, 4 reveal interesting distinctions of magnetic reconnection in the mid tail region and closer to Earth (X m and X n vicinities) in this substorm. As is seen from the comparison of solid and dashed field lines in these figures, reconnection near X n is accompanied by strong changes of the in the noon-midnight meridional plane with overplotted magnetic field lines (solid lines for the moment t 2 and dashed lines for the moment t 0 ) for the 13 February 2008 substorm. The approximate location of the X-lines X n and X m is marked by gray arrows. Magnetic field lines start from the ionosphere at 60 + with 2 + step in latitude. White disks 1R E < r < 9R E in panels (B) and (C) mask magnetic field variations in the inner magnetosphere. Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 644884 6 magnetic field geometry, especially, earthward of that X-line, while near X m the geometry barely changes, which is seen particularly well in the lobe region. The color-coded variations of the z-and x-components of the magnetic field between moments t 0 02 : 30 UT and t 2 02 : 55 UT dB x,z B x,z (t 2 ) − B x,z (t 0 ) in the same noon-midnight meridional plane in Figures 3, 4 quantify these steady and unsteady reconnection regimes.
The difference in dB x,z values in regions x ≈ − 20R E and x ≈ − 31R E seen from Figures 3, 4 suggests that the reconnection process near X m is more steady-state than near X n . The unsteady nature of the near-Earth reconnection is particularly well seen from B z variations earthward of X n in Figure 4. Moreover, the analysis of the equatorial B z (x) profiles with the 5-min cadence provided in Fig. 8i of [27] shows that the main part of the B z changes earthward of X n shown in Figure 4 occurs in the 5-min interval between 02 : 35 UT and 02 : 40 UT. Furthermore, the analysis of the magnetic flux redistribution in the lobes made in [27] gives an estimate of the average electric field in the steady-state reconnection region E y ∼ 0.01v A B 0 /c for B 0 40 nT and v A 1, 000 km/s, consistent with the theoretical estimates that impose the upper limit for the reconnection rate [61, 62, and refs. therein] or ≈ 0.2 [63]. Therefore, one can expect the reconnection near X m to be steady and its electric field homogeneous in space, whereas near X n to be more transient and structured. Below we show that a similar combination of steady and unsteady reconnection regions can be reproduced in PIC simulations of weakly driven during the growth and expansion phases, respectively. The moment of time, when each B z -profile was sampled, is specified by the corresponding colored text in the format "DOY-hour-minute".
Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 644884 magnetotail equilibria with some of the observed growth phase features.

6 AUGUST 2017 SUBSTORM: COMPLEX RECONNECTION PICTURE RESOLVED USING ADVANCED DM METHOD
In this section we consider another substorm event with a more complex structure of dipolarizations and X-lines. It occurred on 6 August 2017 and is also interesting because of a possible X-line crossing detected by the MMS mission. Its signatures were the B z reversal ( Figure 5C), the ions bulk flow reversal and the large dawnward electron bulk flow velocity (not shown). At the same time, at the moment of the reversal the |B x | and B y magnetic field components were relatively large ( Figures 5A,B) so that the total magnetic field exceeded 10 nT. We reconstruct this event using an advanced version of the KNN algorithm where the statistical weights of NNs depend on their proximity to the event of interest (e.g., [29]). In this algorithm, the model magnetic field B (mod) is determined by minimizing the RMS of its deviation from observations B (obs) where S NN is a set of K (B) NN magnetometer measurements of the magnetic field components B j,obs i with ephemeris r (j) , corresponding to the selected set of K NN nearest neighbors; w (0) is the original weighting factor, which is a function of the real-space distance r from Earth, introduced in [44] to mitigate the spatial inhomogeneity of observations, especially at geosynchronous orbit. A distinctive feature of the weighted KNN algorithm is that each term in the sum in Eq. 9 has now an additional weighting factor Here R (j) q is the distance (1) of the corresponding NN from the query point q and R NN is the radius of the sphere containing NNs in the binning space (G 1 , . . . , G 5 ). When σ ≫ 1, all distance-modulated weights w j ≈ 1 and NNs are not weighted. In contrast, for σ < 1, the new weighs w j are well modulated within the sphere R (j) q < R NN . This increases the statistical weight of measurements that were made at the more similar state and input conditions of the magnetosphere, according to the metric (1).
The weighted KNN approach is shown to result in better sensitivity of the model to variations of the magnetospheric state (e.g., storm or substorm phase) by using effectively much smaller numbers of NNs without overfitting [70]. Below we provide the DM reconstruction results with the parameters similar to those used in the previous section and with the weighting factor σ 0.5. Validation results for this event using the MMS1 probe data are presented in the left panels of Figure 5 and they show a reasonable agreement, especially for the B z component, where it does not exceed ∼2 nT.
The reconstruction summary for this substorm in the format used earlier in [26,27] is presented in the right panels of Figure 5. Following [26], we consider the growth phase starting from the first point with B IMF z < 0 in the 5-min cadence series (vertical red dashed line corresponding to t 04:00 UT). The onset time 04:20 UT (vertical orange dashed line) is selected because of the strong change of the negative slope of the AL(t). The start of the recovery phase (23:40 UT, vertical blue dashed line) corresponds to the minimum of the AL index. The recovery phase is postulated to end when the AL > −25 nT, in accordance with [26,27].
Figures 5H-J show weak storm activity: small and constant values of -Sym − H * and low-latitude field-aligned currents FAC R2. In the expansion phase (yellow zone) the amplitude of the TCS ( Figure 5L) decreases, consistent with the earlier DM analyses [26,27] (with small variations of the thickness parameters D and D_{TCS}, according to Figure 5K), while the equatorial magnetic field in the near-Earth tail ( Figure 5M) increases making the magnetic field more dipole-like. At the substorm onset, the evolution of the equatorial magnetic field along the midnight meridian (red lines in Figures 5N, O) reveals wavy perturbations similar to the tearing mode (e.g., Fig. 6.2.9 in [71]). However, their wavelength is rather macroscopic, in contrast to the electron-or ion-scale tearing modes discussed earlier in theory and kinetic simulations of the magnetotail reconnection onset ( [72,73] and refs. therein). Further it results in the formation of new X-lines ( Figure 5O) whose structure and evolution are better seen in Figures 6, 7.
An interesting feature of this event is that the magnetic field dipolarization in the expansion phase has two sub-phases: 04:20-04: 35 UT and 04:40-04:55 UT. Indeed, Figures 6, 7, which describe the evolution of the equatorial magnetic field and current, reveal two successive reconfigurations developing in the premidnight and postmidnight sectors. Figure 6 shows that during the first dipolarization a new X-line forms at x ∼ − 20R E along with the pre-exisitng X-line near x ∼ − 30R E . As it is shown in Figure 8, the used NN subsets are sufficiently extended over the tail to resolve both X-lines. As is seen from Figures 6D-F, the equatorial current during this dipolarizaation becomes bifurcated.
The second dipolarization described in Figure 7 causes stronger and more global changes of the near-Earth magnetic field (regions R(15R E in Figures 7A-C). It also causes not only the formation of new flux ropes in the postmidnight sector but the azimuthal extension of the region of the depressed or even reversed equatorial magnetic field. According to Figures 7D-F, this is accompanied by the reduction of the equatorial current density. To quantify these processes, we integrated the equatorial field B z over arcs similar to dashed blue lines in Figure 7A from the dawn to dusk magnetopause boundaries. Each arc represents a part of the circle with the center (x, y) (3R E , 0) (the shift is used to avoid integration over whole circles within the magnetopause). As was argued in [27], the distribution along the tail of the corresponding integral parameter Int(B z ) B z ds (where ds is the arc length element) may be a good proxy of the magnetic flux evolution in the closed field line region of the magnetotail. The distributions of Int(B z ) along the tail shown in Figure 9 as functions of the arc's Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 644884 8 most tailward value of x indicate that the main flux redistribution in this substorm is provided by the second dipolarization. They also suggest that the main part of the near-Earth dipolarization is provided by the redistribution of magnetic flux in the closed field line region. Indeed, the area under these curves is now magnetic flux. According to Figure 9B, the increase in flux during the second dipolarization in the region 10 − 17R E is roughly equivalent to the decrease in flux at 17 − 26R E . If the dipolarization were provided by an increase of the lobe field reconnection at a30R E that would be seen as a net increase of flux within ∼ 30R E .
To further investigate two dipolarizations occuring during this substorm, we present in Figures 10, 11 the meridional cuts of the cross current and in-plane magnetic field components B x′,x″ and B z in the planes marked by dashed lines in Figures 6, 7 (x ′ GSM and x ″ GSM are the coordinates along the dashed lines in Figures 6, 7).  interval covering the whole expansion phase of that substorm, whereas Figures 10, 11 describe 15-min long partial dipolarizations that constitute the more complex tail reconfiguration during the 6 August 2017 substorm.
As one can seen from the comparisons of Figures 10A,B, the first dipolarization is relatively weak, and it does not cause any significant flux redistribution, according to Figure 9B. During the first dipolarization, the magnetic field variations near X n ( Figures  10C,D) are confined to the region x > − 27R E and −2R E (z(4R E . The near-Earth X-line during this dipolarization forms in the center of the TCS, which extends from −28R E to −9R E ( Figure 10A). It only moderately redistributes its current density ( Figure 10B).
In contrast, during the second dipolarization, the (already shorter, less than ∼ 10R E in the radial extension) TCS disappears ( Figures 11A,B), the near-Earth X-line forms at its tailward end and these processes are associated with a significant flux redistribution shown in Figure 9B (compare yellow and green curves).
It is important to note that these processes of the tail thinning and dipolarization often occur under weak variations of the lobe magnetic field. Its weak variations in the growth phase were reported in [74][75][76] and they are seen in Figure 10C as well as in [27] (Figs. 12, 16, S4 and S13). Even rapid dipolarization processes shown in Figures 4, 11C are accompanied by more gradual lobe field variations, consistent with other data analyses [58,74].

KINETIC SIMULATIONS OF MAGNETOTAIL RECONNECTION GUIDED BY EMPIRICAL RECONSTRUCTIONS
In order to understand the physical mechanisms of the formation of several X-lines in the magnetotail and their different reconnection regimes revealed in the DM analysis, we performed PIC simulations of the tail current sheet equilibrium sharing some of the observed pre-onset tail features. In particular, the reconstruction of the February 13 event discussed above in ( [27], Fig. 8h) suggests that the near-Earth reconnection is preceded by the formation of a flux accumulation region near x ≈ − 22R E . According to Figure 6A, similar pre-onset features in the form of a wide valley with small B z values at R ∼ 22R E and the enhanced B z ridge earthward of that valley took place prior to the 6 August 2017 substorm. To take these features into account, the PIC simulations start from a 2-D equilibrium with a B z hump described by the vector potential A (0) [0, −ψ(x, z), 0], where ψ LB 0 ln(β(x)cosh{z/[Lβ(x)]}), L is the characteristic current sheet thickness parameter, and the x-axis points from Earth to Sun. Its variation along the tail is determined by the function β(x) exp[ε 1 g(ξ)], with ξ x/L, ε 1 ≪ 1 and −g(ξ) ξ + (α/ε 2 ){1 + tanh[ε 2 (ξ − ξ 0 )]}, which provides a region of accumulated magnetic flux near ξ ξ 0 . This is seen from the magnetic field profile B z (x, z 0) ε 1 B 0 {1 + αcosh −2 [ε 2 (ξ − ξ 0 )]} having a characteristic hump. The  corresponding class of isotropic plasma equilibria was first proposed in [77] based on the 2-D generalization [78] of the 1-D Harris model [79] to describe spontaneous onset of the ion tearing instability. The PIC simulations were performed using an open boundary modification [64,80] of the explicit massively parallel code P3D [81] in where d i c/ω pi is the ion inertial scale and ω pi (4πe 2 n 0 /m i ) 1/2 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). The choice of such a relatively long in x and narrow in the y-direction box was motivated by the available computer resources and the necessity to cover a large portion of the tail containing both X-lines resolved by the DM method and described in the previous section. In particular, with d i ∼500 km ∼ 0.1R E [82], the distance between X-lines in our run is 30d i ∼ 3R E , that is only 3-4 times smaller than in the DM reconstruction. At the same time, our previous simulations of similar equilibria with shorter in x and wider in y boxes, up to L (s) y 20d i (see, for instance, Fig. 13 in [82]) suggest that the selected value of L y 5d i with periodic boundaries in the y-direction is sufficient to reproduce major structuring in that direction, including ballooning/ interchange and flapping motions.
The plasma parameters include the mass ratio m i /m e 128, ion-to-electron temperature ratio T i /T e 3 and the effective Alfvén speed v A B 0 / 4πn 0 m i √ c/15 where c is the speed of light. The equilibrium magnetic field parameters are ε 1 0.03, ε 2 0.2, α 3, and ξ 0 −30 with the CS thickness parameter L 1d i . The magnetic and electric fields are normalized, respectively, by B 0 and v A B 0 /c. The coordinates are normalized by d i and velocity components by v A . The simulation grid has 2560 × 160 × 640 cells with ≈ 230 particles per cell corresponding to n n 0 . The magnetic field configuration at the early stage of the run is shown in Figure 12A.
In contrast to earlier simulations ( [83], and refs. therein) that described spontaneous onset regimes, here we drive the system by imposing a weak external electric field E (dr) y at top and bottom boundaries. This setup resembles earlier simulations of the externally driven electron tearing [73], and the whole setup is therefore a combination of the earlier ion and electron tearing modeling efforts. Still, in contrast to earlier setups with localized in x driving fields [73,84,85] and similar to [64], we do not assume any localization of the driving electric field along the tail.  Color-coded distributions of the x'-and z-components of the magnetic field variation between moments t 0 04 : 20 UT and t 2 04 : 35 UT dB x′ ,z B x′ ,z (t 2 ) − B x′ ,z (t 0 ) in the same meridional plane with overplotted magnetic field lines (solid lines for the moment t 2 and dashed lines for the moment t 0 ). x′ GSM is the coordinate along the dashed lines in Figure 6.
Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 644884 insufficiently resolved in observations and it can only be conjectured from global MHD simulations (e.g. [86]). In this situation, the assumption of the homogeneous electric field appears to be the most plausible ad hoc assumption. The driving field amplitude E 0 smoothly increases in a half of the ion gyrotime 1/ω 0i at the beginning of the run and then remains constant with E 0 0.05. The external driving first results in the CS thinning and stretching, which are seen particularly well in the tailward part of the box ( Figure 12B). It also causes the buildup of the plasma pressure in the region x( − 24d i (not shown), consistent with previous studies of the driven reconnection regimes (e.g., Fig. 9 in [73]). This makes the CS configuration more similar to empirical reconstructions with extended TCS, such as for instance in Figures 10A, 11A (see also [27], Figs. 12b). At some point, the first X-line X ′ m forms in the "tailward" part of the simulation box ( Figure 12C). However, the second X-line X ′ n forming later in the left ("earthward") part of the box ( Figure 12D) is not the secondary X-line caused by the tearing instability of the reconnection exhausts (e.g. [87]), because it also forms in the absence of any primary X-lines [88][89][90] or when the primary X-line shows no reconnection signatures [72,82]. X ′ n rather forms because of the flux starvation effect created by the earthward-moving DF in its trailing part. As it was shown in [72,88,89], the DF appears from the original B z hump due to its spontaneous acceleration and further localization in x.
It is very interesting that the magnetic field perturbations shown in Figure 12E strongly resemble the DM reconstructions of substorm dipolarizations shown in Figures 4, 10D, 11D with much stronger bipolar B z perturbations around the near-Earth X-line compared to the midtail region. This suggests that reconnection near X ′ n is unsteady, in contrast to the steady midtail reconnection process at X ′ m . This conclusion is further confirmed in our simulations by the analysis of the electric field and plasma parameters. Figure 13A shows that the electric field distributions around the X-lines are indeed drastically different. Around X ′ m (x ∼ − 50d i ) the distribution of E y (x, z) is homogeneous and its value E y (x, z) ≈ 0.1 is consistent with the theoretical estimates [61][62][63]. These are strong indications of the steady reconnection process. In particular, the broad distribution of the electric field E y over a large area in the plane (x, z) justifies the concept of the reconnection rate, measured by E y , as a global parameter, which characterizes the reconnection process as a whole. In contrast, near X ′ n (x ∼ − 15d i ) the reconnection electric field strongly varies in space. However, not all these variations are related to unsteady reconnection. In particular, the sign-alternating variations of E y near the O-line (x ∼ − 23d i ) describe northsouth flapping motions of the CS as a whole, which are also seen in Figure 13F as strong variations of the magnetic field B x ( ∼ 0.5B 0 ) without noticeable B z variations in the same region ( Figure 13G). The properties of non-reconnection flapping and ballooning/interchange motions (seen in the region x ∼ − 5d i in Figure 13G) in this run with a relatively small extension in the y-direction are similar to non-reconnection motions investigated with larger in y boxes in PIC simulations of spontaneous reconnection onset regimes [82], where they are compared with the corresponding magnetotail observations and other kinetic simulations. At the same time, earthward of X n ′ and near the DF, the electric field is structured in the y-direction due to ballooning/ interchange perturbations that are best seen in variations of the magnetic field B z ( Figure 13G). All in all, the electric field associated with the earthward motion of the DF is highly localized near its B z peak and its value strongly exceeds the steady-state reconnection limit 0.2 [63]. Note, that such strong values of the reconnection electric field were reported before in simulations of the ion tearing instability ( [82], Fig. 5) and interchange-driven reconnection ( [92], Fig. 11). Thus, the kinetic reconnection picture in our PIC simulation, which combines steady and unsteady reconnection regions, is quite consistent with the empirical DM-based reconstructions described in the previous section. Moreover, kinetic simulations reveal its features that cannot be captured from the empirical geomagnetic field analysis, because they represent spontaneous or small-scale plasma modes or they are not reflected in the magnetic field data at all. The examples of the first group of such phenomena are flapping and ballooning/interchange motions seen in Figures 13A,B. They are indeed observed in the tail [93][94][95][96][97], although their relation to substorms and their reconnection modes remains a topic of ongoing discussions [98]. Another example is DFs, with their ion-scale leading edges and fast (v x ∼ v A ) earthward propagation (e.g. [93]), forming out of relatively stationary and macroscopic B z -humps (compare, for instance, Figures 12A,E).
In Figures 13B-D we present another group of signatures, which cannot be resolved using the DM analysis. Figure 13B shows the electric field directed toward the neutral plane z 0 and arising in ion and sub-ion-scale TCS due to different motions of electrons and ions on those scales [85,[99][100][101][102]. Similar effects of the electric field directed toward a negatively charged TCS were shown in PIC simulations ( [101], Fig. 8) and in observations ( [23], Fig. 9). Figures 13C,D show plasma signatures that are usually associated with the electron diffusion region (EDR) in steady reconnection regimes: The first shows super-Alfvénic dawnward electron flows [103] that have been found one of the key distinctive EDR features in recent MMS observations of the magnetotail reconnection [104]. The second reveals non-gyrotropic electron motions that are quantified using the agyrotropy parameter Q e √ proposed by [91] and shown later in MMS observations as a distinctive EDR signature [105]. Finally, in Figure 14 we present the kinetic dissipation parameters for the unsteady part of this run and compare them with similar parameters inferred from MMS observations. In contrast to the steady-state reconnection area near X ′ m , the unsteady reconnection region near X ′ n does not reveal impressive EDR signatures, such as the super-Alfvénic dawnward electron flows or the agyrotropy enhancement (left parts of Figures 4B,C). This is because the main process in this region is the formation and fast earthward movement of a DF and the resulting dipolarization of the magnetic field configuration [72,88,89]. Moreover, many key aspects of these processes can be described by ideal MHD models [106,107], whereas the DF formation and acceleration processes are shown to resemble the ion tearing instability [64,72] supported by the ion Landau dissipation [108]. However, quantifying the latter in simulations and observations is a challenging problem because the conventional single-fluid measure, the Joule heating rate cannot distinguish between ion and electron Landau dissipation in collisionless magnetospheric plasmas. Indeed, the energy conversion rates in the frame moving with ions or electrons j · E ′ e,i (where j j i + j e , E ′ e,i E + v e,i × B/c, j e,i are the electron/ion currents in the laboratory frame of reference and v e and v i are the electron and ion bulk flow velocities) are same for ion and electron species j · E ′ e ≈ j · E ′ i assuming plasma quasineutrality n e ≈ n i ).
To solve this problem, it has recently been proposed [109] to employ the new kinetic parameter Pi − D (α) −Π (α) ij D (α) ij (α e, i), the double contraction of deviatoric pressure tensor ii /3) and traceless strainrate tensor , which was introduced earlier in [110]. It was demonstrated [109] that the Pi − D parameters represent direct analogs of the MHD Joule heating as an entropy variation measure and that they have different distributions for electrons and ions. It was shown that in the regions with ion Pi − D (i) peaks, at the leading part of the DF, ion distributions show signatures of multi-flow motion, including ions reflected from the DF. Such multi-flow ion motions have indeed been detected at DFs in Cluster, THEMIS, and MMS observations [111][112][113][114].
In Figures 14A-D we present kinetic dissipation measures obtained in PIC simulations and averaged over the y direction 0 < y < 5d i , along with the corresponding profile of the magnetic field B z shown here to provide the global context for this local investigation. As one can see from Figures 14A,B, while the linear distribution of the electron dissipation parameter 〈Pi − D (e) 〉 y remains irregular and not obviously positive, its FIGURE 13 | Steady and unsteady reconnection regions in weakly driven magnetoatail at the moment ω 0i t 45.8 corresponding to Figure 12E. (A-D) The distributions in the plane y 2.5d i of the electric field components E y and E z , the electron bulk flow velocity −V ey and the electron agyrotropy parameter Q e [91] marking the localization of the electron diffusion region (in the latter case, to reduce noise in simulation outputs, the original numerical distributions are averaged over 20 × 20 grid cells). (E-G) The distributions of the electric field E y and the magnetic field components B x and B z in the neutral plane z 0.
integration along the tail reveals its persistent accumulation upstream of the DF structure (red line in Figure 14B). The increase starts from the X ′ n vicinity with another buildup near the corresponding O-line. The ion dissipation parameter is even more impressive: Already its average over the y-coordinate reveals a peak near the DF ( Figure 14C), and when integrated along the tail x 0 〈Pi − D (i) 〉 y dx builds up near the DF and remains elevated farther in the tail (red line in Figure 14D). Figures 14E-H show the dissipation parameters similar to those in Figures 14A-D but now derived from MMS observations of a DF on 6 July 2017, a relatively rare case of a slow moving DF with the ion bulk flow speed smaller than 200 km/s. The four-probe sub-ion-scale MMS observations of the electromagnetic field and plasma parameters provide the unique opportunity to measure the kinetic dissipation parameters Pi − D for both electrons and ions. At the same time, even with the MMS capability of calculating higher moments of the plasma distribution, the assessment of the kinetic dissipation parameters remains a challenging problem. In particular, even in the MMS burst mode with the sampling time δt 0.15 s [115] and probe spacing δr(20 km, any velocity gradient estimates necessary for calculation of the tensor D (i) ij may give trivial results for structures moving much faster than V max δr/δt ≈ 133 km/s. Thus, MMS data is only appropriate so far to study the kinetic dissipation in relatively slow moving DFs.
In spite of these caveats, simulation and observation results presented in Figure 14 (Figures 14B,F) the electron dissipation builds up behind the DF, upstream of the ion dissipation buildup, in the regions with relatively small values of the magnetic field, while for ions the dissipation starts accumulating at or even before the DF.

Error Analysis of Empirical Reconstructions
In this study we provided a DM reconstruction of magnetic reconnection in the Earth's magnetotail associated with its dipolarizations during substorms. A direct validation of this reconstruction can only be provided using a limited number of in-situ observations available at the moment of interest. This is an unavoidable feature of the DM method as a data discovery x 0 〈Pi − D (i) 〉 y dx (red line) at the moment ω 0i t 45.8, which is also described in Figures 12E, 13. In all panels the magnetic field B z profile is shown by blue lines. Panels tool, which extracts from data the information (e.g., on the global structure of the magnetotail), which cannot be obtained by other methods. We simply have no real constellations of ∼ 3 · 10 4 probes to comprehensively validate our results. Still, the 13 February 2008 substorms were validated by all five THEMIS probes (Figs. S6-S7 in [27]), while for the 6 August 2017 event, the MMS1 validation results are presented in Figure 5. Moreover, the uncertainty of the DM method caused by averaging over the NN bins can be quantified by comparing the original values of the parameters G 1 -G 5 with their NN means. For the 13 February 2008 reconstruction such information was provided in Fig. 19 of [27]. For the 6 August 2017 substorm we provide it in Figure 15. This figure shows in particular that during the dipolarization intervals considered in Section 4 and shown by vertical dashed lines, the maximum deviation of the binning parameters averaged over their NN bins from their original values defined by Eqs. 4-6 does not exceed ∼10% (the largest deviation is seen for 〈AL| at the end of the second dipolarization interval). This means that statistical errors of the presented reconstruction of the magnetic field during this substorm are much smaller compared to major variations of the binning parameters. Therefore the presented DM-based picture of magnetotail reconfigurations should indeed reflect the characteristic features of magnetic reconnection during substorms. Consistent with the analysis of the 13 February 2008 substorms [27], we have found that the relatively strong deviations of the binning parameters from their means over NNs take place for the solar wind parameter 〈vB IMF s and the AL index in the recovery phase. This suggests that the solar wind and the magnetosphere after substorms are less coherent (perhaps turbulent) and hence less reproducible, compared to the evolution of the magnetosphere during growth and expansion phases.
An important source of uncertainty in the present NN approach may be the instrument errors and combining probes from different epochs. Fortunately, the accuracy of magnetic field measurements critical for our investigation (with a few nT accuracy necessary to resolve the B z magnetic field in the tail) was sufficiently high. In particular, the IMP8 magnetometer was good to 0.3 nT [116] and later missions had largely better instruments (e.g. [117][118][119]) with a few caveats. Significant errors (up to 7 nT) were found for some geosynchronous missions [120] and they were mitigated using inter-spacecraft calibration. The errors in the external magnetic field (difference between the measured and dipole magnetic field values) may also be large in the inner magnetosphere because of the spacecraft attitude uncertainty and large values of the dipole field there [121]. However, this is not an issue in the magnetotail.

Implications for Local Reconnection Models and Tearing Stability
The concept of magnetic reconnection was introduced to explain explosive energy release and rapid changes of magnetic field topology associated with solar flares [2,3], magnetospheric substorms [7,108,122] and laboratory current disruptions ( [123], and refs. therein). But its theory turned out to be built mainly on models of steady-state reconnection regimes ( [63,[124][125][126][127], and refs. therein). The few exceptions include the tearing instability theory [87,108,122,128], and catastrophe models of coronal mass ejections and solar flares [129,130].
At the same time, the description of transition from the slow evolution of the tail to its rapid reconfiguration has long been complicated by the almost universal tearing stability of the tail current sheet provided by magnetization of electrons due to nonzero northward magnetic field B z [131,132]. As a result, the tail can be unstable when electrons become unmagnetized, under the condition B z /B 0 (kρ 0e , where B 0 is the field outside CS, k is the wave vector and ρ 0e is the thermal electron gyroradius in the field B 0 [73,122,133]. The resulting electron tearing instability is enabled by the free energy of the mutual attraction of the parallel electric current filaments and the electron Landau dissipation of unmagnetized electrons. In PIC simulations, the corresponding electron-demagnetization mediated reconnection (EDMR) onset used to be reproduced due to stretching and thinning of a CS by the external electric field [73,101,134]. It is important that after the electron tearing instability phase (or in its absence in simulations with spatially localized driving [84,85]) the reconnection process becomes quasi-steady ( [83], and refs. therein), consistent with regimes found earlier in kinetic simulations with non-self-consistent setups using 1-D CS equilibria with an imposed X-line perturbation ( [127], and refs. therein).
In 1974 Schindler [108] hypothesized that the tail could become unstable even with magnetized electrons if the CS is sufficiently thin to demagnetize ions and provide their Landau dissipation. The corresponding tearing instability must be much faster compared to the electron tearing. However, later it was found [135] that magnetized electrons change the free energy of the tearing mode, and eventually Lembege and Pellat [131] showed that the corresponding sufficient stability condition coincides with the Wentzel-Kramers-Brillouin (WKB) approximation π(B z /B 0 )(kL z , which allows one to consider stability neglecting the CS variations along the tail with the scale L x ∼ L z (B 0 /B z ) (L z is the CS half-thickness) making the ion tearing impossible. A missing key for ion tearing destabilization was found relatively recently when it was discovered [77] that the stability condition derived by Lembege and Pellat [131] is only valid for constant B z values. If B z changes along the tail, the stability condition takes the form π(B z /B 0 )C 2 d (kL z , where the parameter C d VB z /(πL z ) is determined by the flux tube volume per unit magnetic flux V dl/B. In particular, in the presence of a flux accumulation area with the tailward gradient of B z , the parameter C d > 1 and a room for instability arises. The corresponding instability had indeed been found in PIC simulations with ad hoc configurations having B z (x) profiles with a hump [72,88,89]. Since electrons remained initially magnetized by the field B z , the instability was essentially the ion tearing. It first led to the formation of an earthward-moving dipolarization front (DF), in whose wake new X-lines formed due to the flux starvation process [89]. The resulting iondemagnetization dominated reconnection (IDMR) onset did not require any external driving and could be considered as spontaneous or "internally driven" by the DF formation and evolution processes.
Despite this clarity in the tearing stability theory and consistent simulation results, until now, the role of EDMR and IDMR regimes in the actual magnetotail dynamics remained unclear. In particular, it is unknown if/when the driving (ultimately due to the solar wind) is sufficiently strong to squeeze the CS down to electron scales and to provide EDMR with the subsequent steady reconnection, and when (if any) B z humps form to provide IDMR.
The present study provides interesting implications for the magnetotail stability and reconnection onset mechanisms. Our DM reconstructions suggest that both steady and unsteady reconnection regimes are possible in the magnetotail during substorms. At the same time, our PIC simulations guided by empirical reconstructions suggest that both IDMR and EDMR regimes are possible in the tail. Moreover, the former resembles the unsteady reconnection, while the latter becomes eventually steady, consistent with the classical fast and steady reconnection models ( [62] and refs. therein).

Role of Thin Current Sheets
The use in Section 5 of isotropic plasma equilibria with shifted Maxwellian distributions for ions and electrons, inherited from the 1962 Harris solution [79], to explain the reconnection features found in our DM reconstructions may be questioned in view of another discovery in the DM analysis of substorms, namely the buildup of extended TCSs in the substorm growth phase and their decay in the expansion phase [26,27] (see also Figures 10A,B,  11A,B of the present study).
The analysis of 2-D isotropic equilibrium models [136] suggests that they require strongly stretched magnetic field configurations (with sufficiently large values of the ratio B 0 /B z ) to explain the formation of the ion-scale TCS sufficiently far from the Earth. Large values of B 0 /B z are required to maintain the force balance between the magnetic field line tension and the pressure gradient 1/L x ∼ (B z /B 0 )/L z , where L x is the inhomogeneity scale of the TCS, L z is its half-thickness and B 0 is the lobe field [137]. Modeling TCSs with L x ≫ L z (B 0 /B z ) might require more sophisticated equilibria with anisotropic and agyrotropic particle distributions (e.g. [136], and refs. therein).
Indeed, three of four substorm events on 13 February 2008 considered in [27] had relatively small values of B 0 /B z ∼ 10 (according to their Figs. 15b-15d), whereas their aspect ratios L x /L z often exceeded 50 ( Fig. 16 in [27]). That finding was consistent with signatures of the multiscale structure of the magnetotail inferred from local observations of the pre-onset CSs [19,20,23,51].
However, this is not the case for the event considered in Section 3, whose specific features (the B z hump and the ion-scale TCS earthward of it) guided our PIC simulations. In that first substorm of the 13 February 2008 series, the ratio B 0 /B z reaches 70 in the late growth phase (yellow line in Fig. 15a, corresponding to 02:25 UT). Thus, the specific substorm event, considered in Section 3 of our DM analysis is close to the isotropic force balance state and it can be consistently described by 2-D isotropic CS equilibrium models of the class [78]. Moreover, the specific parameters used in our simulations correspond to L x /L z ∼ B 0 /B z ≈ 33 and D TCS ∼ 0.1R E and they are quite close to similar TCS parameters of the first substorm in the 13 February 2008 series in its late growth phase (02:25 UT): L x /L z ∼ 25, B 0 /B z 20 − 70 and L z 0.2R E .
One can also provide more general arguments why the isotropic 2-D models can still be used in the local stability analysis of the realistic magnetotail. First, statistical studies show that the tail plasmas away from the dipole region are weakly anisotropic [138,139]. At the same time, the DM reconstructions demonstrate that the current of the embedded TCS in the late growth phase may be small, compared to the total current, as is seen, for instance, from Figure 5L (this is the case for all four 13 February 2008 events as is seen from Fig. 8f in [27]). This suggests that the embedded TCS features and underlying non-isotropic plasma properties may only serve to provide the formation of the ion-scale TCSs sufficiently far from Earth, where their local stability properties can still be realistically reproduced by PIC simulations with isotropic equilibria and open x-boundaries. This is consistent with the results of statistical studies based on Geotail data [140], which suggest that the near-Earth X-line mainly forms near the tailward edge of the TCS. This appears to be the case during the second dipolarization in the 6 August 2017 event (Figure 11), although this is likely not the case during the first dipolarization when the near-Earth X-line forms in the middle of a very long TCS ( Figure 10). Besides, even if the initial TCS is relatively short because of the corresponding force balance, the simulations performed in Section 5 suggest that it becomes more stretched and closer to empirical TCS reconstructions due to the external driving. To conclude, while some substorm dipolarizations certainly require a generalization of the isotropic plasma approximation, as it was outlined in [136], others can still be described using the conventional class of isotropic CS models [78,79].

CONCLUSION
In this study, we investigated for the first time the magnetotail reconnection picture using modern data-mining methods, which allow us to employ for the reconstruction not only the magnetic field measurements available at the moment of interest but also other events in the historical database when the magnetosphere was in similar global states (substorm phases). The DM reconstruction revealed two distinctly different regions of magnetic reconnection with weak and strong changes of the magnetic field geometry. For both the 13 February 2008 and the 6 August 2017 substorms considered in our study the near-Earth X-line appears near x −20R E at the substorm onset, which is defined in our work as a transition to the AL(t) index evolution with a strong negative slope (dashed vertical orange line in Figure 5G-M). This result is consistent with the original conjecture of Hones [7], later single-and multi-probe studies of the near-Earth X-lines [57,58], as well as with the plasmoid statistics [141]. In both events, the near-Earth X-line first appears in the pre-midnight sector (Figures 2, 6B), which is consistent with the earlier statistical investigations using Geotail [142] and Cluster [143] data.
In addition to earlier investigations, our DM reconstruction reveals that the near-Earth X-line (X n ) often co-exists with another more distant midtail X-line (X m ) located at x ≈ − 30R E . In spite of the fact that its location is near the edge of the main cloud of historical magnetometer measurements [44], the analysis of data in the corresponding NN bins (Figures 2, 8) shows that the selected NN subsets provide sufficiently broad radial coverage of data to resolve both X-lines. The finding of the midtail X-line is consistent with another group of earlier observations suggesting persistent reconnection in the midtail around 30R E , which was inferred from THEMIS and ARTEMIS statistics of traveling compression regions [59,60]. However, the coexistence of near-Earth and midtail X-lines has never been demonstrated before.
Moreover, the DM analysis shows that reconnection regimes at near-Earth and midtail X-lines are different. The near-Earth X-line appears at the substorm onset and then disappears from that region or reappears in another near-Earth region, e.g., in the postmidnight sector (compare Figures 7A,B or Figures 2, 10 in [27]). In contrast, the midtail X-line, after its appearance within the reconstruction validity region (here R < 32R E ) in the late growth phase remains relatively stable and only gradually approaching the Earth (Figures 6A-C, 7A-C). Furthermore, the analysis of the magnetic field changes in the meridional plane (Figures 3, 4), which according to the Faraday's law 8) quantifies the steadiness of the reconnection process, suggests that the latter is relatively steady near X m and transient at X n .
To understand the physical mechanisms of the formation of several X-lines in the magnetotail and their different reconnection regimes, we performed 3-D PIC simulations of a relatively long (L x 80d i ) tail CS region with open boundaries in the Sun-Earth direction. A new aspect of simulations was the combination of the initial TCS configuration having a region of the flux accumulation (B z hump) with a relatively weak and homogeneous external driving. The formation of the flux accumulation regions prior to unsteady reconnection in the near-Earth tail is found in the DM reconstruction of both substorm events (Fig. 8h in [27], as well as Figures 6A, 7A), consistent with earlier statistical results [144,145]. Recently, it has been inferred from remote-sensing observations of 30-100 keV energy electrons precipitating from the tail CS during the substorm growth phase [146]. This feature is also interesting because the corresponding region with the tailward B z gradient (earthward of the B z hump) has been found in the tail stability theory [77] to be the only mechanism of destabilization of the ion tearing mode [108]. The second feature, the external driving was used before to reproduce the tail reconnection onset through the electron tearing instability [73]. It was also used in forced reconnection models [84,85].
The reconnection picture in PIC simulations, guided by the DM reconstructions, is found to be surprisingly consistent with the empirical picture of the magnetotail reconnection. It also reveals two reconnection areas with distinctly different reconnection regimes, whose steadiness can now be checked using the explicit distributions of the electric field in the meridional plane ( Figure 13A). It is found that farther in the tail, the reconnection process is steady and it reveals many signatures of the sustained collisionless reconnection process with the region of agyrotropic electron motion in its center. The corresponding dusk component of the electric field is broadly distributed in the meridional plane and hence it becomes effectively a global parameter of this reconnection regime. Its value E y ≈ 0.1 matches earlier theoretical estimates for this regime supported by local PIC simulations [61][62][63]. At the same time, the evolution of the B z hump is found to result in an unsteady reconnection process with the peak electric field near the dipolarization front exceeding the steady reconnection rate limit by more than an order of magnitude, the result, which is consistent with earlier PIC simulations of local unsteady reconnection regions [82,92]. The analysis of kinetic dissipation parameters in the unsteady reconnection region shows that the ion dissipation parameter Pi − D (i) peaks near the DF and it is further accumulated upstream of the propagating front. The electron dissipation is largely accumulated behind the DF near new X-and O-lines. Similar ion and electron dissipation parameters are inferred from MMS observations. Both empirical and first-principle pictures of magnetotail reconnection still need further refinement. The present DM approach provides an empirical picture on the magnetotail on the time scales greater than 5 min and on the spatial scales larger than ∼ 0.2R E for the TCS thickness and a few R E in the equatorial plane. On these scales, the magnetic field dipolarization is likely a cumulative effect of the smaller-scale processes, such as multiple DFs (e.g. [58,66,147]). These cumulative effects are not yet reproduced in PIC simulations. On the other hand, the midtail X-lines are found close to the gap region 31R E < R < 55R E in historical data [44,148]. Thus, a better resolution of the midtail reconnection picture requires more measurements in that gap region. PIC simulations were made in a relatively thin CS, whose non-Harris properties, such as its negative charging and multiscale structure, are only partially captured now due to the external driving. In simulations with thicker CSs and broader B z humps, as well as more realistic values of the parameters m i /m e and c/v A , one can expect stronger negative charging effects, slower growth of DFs and subsequent reconnection, as well as weaker electron dissipation. A further improvement of the tail reconnection and stability picture is also required to better reproduce less stretched embedded TCS.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: http://doi.org/10.5281/zenodo.4383387.

ACKNOWLEDGMENTS
We acknowledge interesting and useful discussions of the obtained results with Nikolai Tsyganenko, Rumi Nakamura, Ferdinand Plaschke, Slava Merkin and Shin Ohtani. We thank the many spacecraft and instrument teams and their PIs who produced the data sets we used in this study, including the Cluster, Geotail, Polar, IMP-8, GOES, THEMIS, Van Allen Probes and MMS, particularly their magnetometer teams. We also thank the SPDF for the OMNI database for solar wind values, which is composed of data sets from the IMP-8, ACE, WIND, and Geotail missions, and also the WDC in Kyoto for the Geomagnetic indices. Simulations were made possible by the NCAR's Computational and Information Systems Laboratory (doi: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.