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–21]. The latter may also be split in two current layers forming bifurcated TCSs [19, 22–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.

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 N N ≫ 1 moments in the past. These moments called “the nearest neighbors” or NNs are similar to the event of interest in terms of similar values of the geomagnetic indices Sym-H and AL, their time derivatives and the solar wind input parameter v B s I M F . Here B s I M F is the southward component of the Inteplanetary Magnetic Field (IMF): B s I M F = − B z I M F if B z I M F < 0 and B s I M F = 0 otherwise (The Geocentric Solar 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 D B ≫ K N N . 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 N N 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 N N moments to form an event-oriented 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 N N , 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 R q ( i ) = ∑ k = 1 5 ( G k ( i ) − G k ( q ) ) 2 / σ G k 2 . (1) where σ G k is the standard deviation of the component G k and the coordinates G 1 - G 5 are defined by the formulae: G 1 ( t ) = 〈 S y m − H * | ∝ ∫ − Π s t / 2 0 S y m − H * ( t + τ ) cos ( π τ / Π s t ) d τ (2) G 2 ( t ) = D 〈 S y m − H * | / D t ∝ ∫ − Π s t / 2 0 S y m − H * ( t + τ ) sin ( 2 π τ / Π s t ) d τ (3) G 3 ( t ) = 〈 A L | ∝ ∫ − Π s s t / 2 0 A L ( t + τ ) cos ( π τ / Π s s t ) d τ (4) G 4 ( t ) = D 〈 A L | / D t ∝ ∫ − Π s s t / 2 0 A L ( t + τ ) sin ( 2 π τ / Π s s t ) d τ (5) G 5 ( t ) = 〈 v B s I M F | ∝ ∫ 0 τ ∞ v B s I M F ( t − τ ∞ + τ ) exp [ ( τ − τ ∞ ) / τ 0 ] d τ (6)

Here S y m − H * = A ⋅ S y m − H − B ⋅ P d y n is the pressure-corrected Sym-H index [30], P d y n 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 A L 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–5) is provided by the sine and cosine kernel functions and by the exponential function in Eq. 6. The averaging scaling parameters Π s t = 12   hr and Π s s t = 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 A L index and its time derivative ( G 3 , G 4 ) determine the strength and phase of the substorm activity, because the A L 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). 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 A L index, average v B s I M F 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 storm-time model, TS07D [36]. The latter was also justified by the empirical relationship between the v B s I M F parameter and the D s t 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 storm-substorm 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 D B = 3,668,101 records of the magnetic field vector with 5 and 15-min cadence inside and outside 5 R 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 A L , Sym-H and v B s I M F time series with 5-min cadence. At every moment t q , the subset is found as K N N points i = 1 , … , K N N satisfying the condition R q ( i ) < R N N ( K N N ) , where R q ( i ) 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 N N ∼ 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 N N values were found to cause overfitting.

At every moment of interest t q , the K N N 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 N N ≫ 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: B ( e q ) ( ρ , ϕ , z ) = B ( e q ) ( ρ , ϕ , z ; D ) + B ( e q ) ( ρ , ϕ , z ; D T C S ) (7) 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 T C S 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 2 M angular Fourier harmonics (even and odd parity in ϕ) with the total number of N + 2 M ⋅ 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 ζ = z 2 + D 2 . 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 0 n ( s ) : ( A ϕ ) 0 n ( s ) = 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 T C S , are determined from fitting the model to data. In particular, to distinguish between thick current sheets and TCS, we impose the condition D T C S < D T C S ( max ) = 1 R 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 T C S 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 d y n 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 d y n . 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 F A C = 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 ion-scale TCS ( D T C S ≪ 1 R E ) inside a much thicker CS ( D > 2 R 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² 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 T C S ≈ 0.2 R 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 | < 1 R E and | z | < 5 R 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 N N = 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 20 R 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 20 R 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.

In this section we describe the global structure and dynamics of reconnection on the example of a relatively small and short substorm (13 February 2008 02:05–02:55 UT) considered earlier in [27] with the reconstruction parameters ( M , N ) = ( 6,8 ) , K N N = 32,000 and N F A C = 16 . The latter analysis is extended here by increasing the maximum radial distance of the spacecraft data used in the reconstruction from 31 R E to 35 R 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 ≈ − 20 R E ) and midtail ( x < − 27 R 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 single-point 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 30 R 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 ∂ E y / ∂ x = − ∂ B z / ∂ t , ∂ E y / ∂ z = ∂ B x / ∂ t , (8)

It suggests that the temporal variations of B x and B z magnetic field components determine the spatial gradients of the dawn-dusk (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–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–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 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 d B 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 d B x , z values in regions x ≈ − 20 R E and x ≈ − 31 R 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.01 v 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 E y / ( v A B 0 / c ) ∼ 0.1 [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 magnetotail equilibria with some of the observed growth phase features.

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 ( o b s ) M e r r ( N N ) = ∑ j ∈ S N N ∑ i = x , y , z w j w ( 0 ) ( r ) [ B i ( mod ) ( r ( j ) ) − B i j , o b s ] 2 , (9) where S N N is a set of K N N ( B ) magnetometer measurements of the magnetic field components B i j , o b s with ephemeris r ( j ) , corresponding to the selected set of K N N 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 w j = exp [ − ( R q ( j ) / σ R N N ) 2 / 2 ] , (10)

Here R q ( j ) is the distance (1) of the corresponding NN from the query point q and R N N 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 q ( j ) < R N N . 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 z I M F < 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 A L ( t ) . The start of the recovery phase (23:40 UT, vertical blue dashed line) corresponds to the minimum of the A L 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 - S y m − 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 ∼ − 20 R E along with the pre-exisitng X-line near x ∼ − 30 R 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 ≲ 15 R 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 ) = ( 3 R 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 I n t ( B z ) = ∫ B z d s (where d s 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 I n t ( B z ) along the tail shown in Figure 9 as functions of the arc’s 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 − 17 R E is roughly equivalent to the decrease in flux at 17 − 26 R E . If the dipolarization were provided by an increase of the lobe field reconnection at ≳ 30 R E that would be seen as a net increase of flux within ∼ 30 R 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 ′ G S M and x ″ G S M are the coordinates along the dashed lines in Figures 6, 7). These figures show processes similar to the 13 February 2008 dipolarization and shown in Figs. 12 and 13 in [27] as well as in Figures 3, 4. In particular, both dipolarizations reveal stronger variations of d B x , x ′ and B z magnetic field components in the near-Earth reconnection region (around X n ) compared to the midtail one (around X m ). The differences between the magnetic field variations in Figures 10, 11 and those in Figures 3, 4, such as for instance, different relative phases of d B x and d B z variations can be explained by the larger time difference used in case of the 13 February substorm to calculate these variations. In fact, Figures 3, 4 describe the magnetic field variations during the 25-min long 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 > − 27 R E and − 2 R E ≲ z ≲ 4 R E . The near-Earth X-line during this dipolarization forms in the center of the TCS, which extends from − 28 R E to − 9 R E (Figure 10A). It only moderately redistributes its current density (Figure 10B).

In contrast, during the second dipolarization, the (already shorter, less than ∼ 10 R 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–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].

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 ≈ − 22 R E . According to Figure 6A, similar pre-onset features in the form of a wide valley with small B z values at R ∼ 22 R 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 ψ = L B 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 a 3-D box with dimensions L x ( s ) × L y ( s ) × L z ( s ) = 80 d i × 5 d i × 20 d i , where d i = c / ω p i is the ion inertial scale and ω p i = ( 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.1 R E [82], the distance between X-lines in our run is 30 d i ∼ 3 R 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 y ( s ) = 20 d i (see, for instance, Fig. 13 in [82]) suggest that the selected value of L y = 5 d 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 = 1 d 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 y ( d r ) 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. It remains constant throughout the main part of the box length L x , being only attenuated near open boundaries. For example near the left boundary E y ( d r ) = E 0 [ tanh ( x − δ ) / δ ] 2 with δ = 0.1 L x . The actual structure of the driving electric field remains 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 / ω 0 i 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 ≲ − 24 d 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–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 ∼ − 50 d 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–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 ∼ − 15 d 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 ∼ − 23 d i ) describe north-south 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.5 B 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 ∼ − 5 d 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–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–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 quasi-neutrality n e ≈ n i ).

To solve this problem, it has recently been proposed [109] to employ the new kinetic parameter P i − D ( α ) = − Π i j ( α ) D i j ( α ) ( α = e , i ), the double contraction of deviatoric pressure tensor Π i j ( α ) = p i j ( α ) − δ i j p ( α ) (where p ( α ) = p i i ( α ) / 3 ) and traceless strain-rate tensor D i j ( α ) = ( ∂ i v j ( α ) + ∂ j v i ( α ) ) / 2 − δ i j θ ( α ) / 3 (with θ ( α ) = ∇ ⋅ v ( α ) ), which was introduced earlier in [110]. It was demonstrated [109] that the P i − 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 P i − 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–114].

In Figures 14A–D we present kinetic dissipation measures obtained in PIC simulations and averaged over the y direction 0 < y < 5 d 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 〈 P i − D ( e ) 〉 y remains irregular and not obviously positive, its 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 ∫ 0 x 〈 P i − D ( i ) 〉 y d x 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 P i − 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 j ( i ) may give trivial results for structures moving much faster than V m a x = δ 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 have interesting similarities. In particular, both simulations and data show the accumulation of positive P i − D values for electrons after integration (over x in simulations and in time in observations) seen in Figures 14B–G. Both in simulations and in observations (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.

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 = − 20 R E at the substorm onset, which is defined in our work as a transition to the A L ( 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 ≈ − 30 R 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 30 R 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 < 32 R 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 = 80 d 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–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 P i − 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.2 R 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 31 R E < R < 55 R 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.