Finding reconnection lines and flux rope axes via local coordinates in global ion-kinetic magnetospheric simulations

Magnetic reconnection is a crucially important process for energy conversion in plasma physics, with the substorm cycle of Earth's magnetosphere and solar flares being prime examples. While 2D models have been widely applied to study reconnection, investigating reconnection in 3D is still, in many aspects, an open problem. Finding sites of magnetic reconnection in a 3D setting is not a trivial task, with several approaches, from topological skeletons to Lorentz transformations, having been proposed to tackle the issue. This work presents a complementary method for quasi-2D structures in 3D settings by noting that the magnetic field structures near reconnection lines exhibit 2D features that can be identified in a suitably chosen local coordinate system. We present applications of this method to a hybrid-Vlasov Vlasiator simulation of Earth's magnetosphere, showing the complex magnetic topologies created by reconnection for simulations dominated by quasi-2D reconnection. We also quantify the dimensionalities of magnetic field structures in the simulation to justify the use of such coordinate systems.

1 Introduction

A long-standing issue in analysing data from magnetospheric simulations has been the identification of topological features in the magnetic field, especially in relation to reconnection and flux transfer events. These features include magnetic neutral lines; separators; and, loosely, X and O lines. In a 2D configuration, the problem is tractable via the use of a flux function (see Sonnerup, 1970, for an early example, with details given in, for example, Servidio et al., 2009). Recent examples of flux function use include Hoilijoki et al. (2017), Palmroth et al. (2017) and Hoilijoki et al. (2019). A generalized flux function can be defined in certain 3D configurations (Yeates and Hornig, 2011). Global ideal magnetohydrodynamic (MHD) simulations present a relatively simple reconnection topology that is tractable via the four-field junction (FFJ) method that identifies reconnection topologies from global magnetic field connectivity (whether or not field lines connect to Earth or to the solar wind) (Laitinen et al., 2006). However, more detailed simulations, such as Vlasiator ion-kinetic simulations, produce complex reconnection topologies impossible for FFJ to parse completely (see Pfau-Kempf et al., 2020). As such, the space physics community is challenged to find and develop new tools for analysing these topologies.

Reviewing the definitions of the topological structures of reconnecting 3D magnetic fields by Parnell et al. (2010), we can summarize three different but related concepts: magnetic null points (with several subtypes), where the magnetic field B completely vanishes; separatrices (also called separator surfaces or separatrix surfaces), which are planes that delimit different magnetic domains; and separators, which are lines connecting magnetic null points, defined as intersections of paired separatrix surfaces lying at the boundary of four different magnetic domains. The behaviour of magnetic null points has been studied in much detail – in 3D as well (e.g. Greene, 1988) – for their usefulness in reconnection analysis. Detection and classification of null points can be done readily with existing methods: Olshevsky et al. (2015, 2016) use the Poincaré index method to detect and characterize null points (where B(r)=0) in local MHD and particle-in-cell (PIC) simulations, which can be applied to spacecraft data as well. Haynes and Parnell (2007) detail methods of finding null points within a cubic cell, using a tri-linear method to find common zero-crossings for all magnetic field components. Fu et al. (2015) introduced the first-order Taylor expansion (FOTE) method for reconstructing the magnetic field topology around nulls and finding null points from multi-point spacecraft data.

As, for example, Lau and Finn (1990) note that the magnetic field lines connecting specific null-point pairs – separators – are important for reconnection. For example, Dorelli and Bhattacharjee (2009) discuss topological changes in dayside separators during flux transfer event (FTE) formation in MHD simulations. Several methods for finding separator lines have been developed. For example, Komar et al. (2013) trace related separators by stepping from magnetic null points and inspecting magnetic connectivity of candidate points, which is akin to the FFJ method used by Laitinen et al. (2006). This is used by Eggington et al. (2022) to find the main magnetospheric separator through tracing magnetic connectivity; it is noted that they found it complicated to find the initial magnetic nulls for tracing. On the other hand, Olshevsky et al. (2016) show a plethora of null points found in a magnetospheric simulation, representing a challenge for mapping the pairwise connections of null points. Glocer et al. (2016) track magnetic separators in MHD simulations and compare different techniques for it, while Haynes and Parnell (2010) introduced generalized methods for constructing the topological skeletons of magnetic fields. This includes null points, separatrices and the separator lines, not tied to specific coordinate systems, with involved field line tracing procedures. Recently, Bujack et al. (2021) also introduced similar tools. Such tools are likewise being explored in Vlasiator 3D magnetospheric simulations by Bouri et al. (2023) in order to apply machine learning methods to the topic. Machine learning methods have also seen use in the analysis of plasma simulations by, for example, Bussov and Nättilä (2021).

Closely related to the concept of separators are the magnetic X and O lines. As the intersections of separatrix surfaces, the separators delimit different magnetic domains. In the classic picture of reconnection, inspecting the magnetic field lines on a plane normal to the tangent of the reconnection line displays an X topology, with the in-plane components of B changing signs at the X line (or an X point when constrained to such a plane). This holds regardless of an out-of-plane B component, as noted by Parker (1957) and Sonnerup (1974): the hyperbolic structure remains visible in the in-plane components. Likewise, the in-plane magnetic field components change signs at an O line (or O point). This can be used to extend the definition from a true null line, defined as a continuous line with exactly zero B, which is structurally unstable (Priest and Titov, 1996) and is expected to split into null points connected by a separator. As this definition of a null line is not consequently fulfilled in practice and because we would like to identify null lines whether or not a magnetic guide field is present, we use the term in-plane null line (and, more specifically, in-plane X and O lines) in this paper to mean such sets of points that satisfy the zero-crossing condition for the in-plane magnetic field components, allowing for a common category of null lines that covers both X and O lines, with and without guide field. It should be noted, however, that the notion of an in-plane null depends on the definition of the plane, and it is not clear what plane should be chosen for this decomposition. This work proposes one possible selection of a useful coordinate system via variance analysis for this purpose. Other selections of a coordinate system are possible, as shown by, for example, Genestreti et al. (2018).

For magnetic field structures such as current sheets (CSs), null lines and flux ropes, it is useful to select a suitable coordinate system that can exploit the properties of the structure to simplify the problem at hand. Methods such as minimum variance analysis (MVA, Sonnerup and Scheible, 1998) are widely used in spacecraft observational studies to investigate structures in amenable coordinate systems (Paschmann and Daly, 1998). Observational methods for the construction of local coordinate systems may require some situational adjustments, as in Gosling and Phan (2013) and Hietala et al. (2018), with a hybrid MVA coordinate system. Indeed, Genestreti et al. (2018) examine 14 different variations of generating local coordinate systems for a reconnection event observed by the Magnetospheric Multiscale (MMS) mission. The magnetic field structures that are more frequently investigated (such as CSs and flux ropes) often display more variability in some dimensions than in others, allowing for some simplifications based on the dimensionality of the structures. Methods for inspecting the dimensionality and for forming local coordinate systems are reviewed by Shi et al. (2019). Particularly, the minimum directional derivative (MDD) and minimum gradient analysis (MGA) methods may be used to define local coordinate systems based on the Jacobian of the magnetic field B, which we will build upon in this work. These are detailed in Sect. 1.1.

Figure 1 A schematic view of Earth's magnetosphere under southward-IMF driving, with zoomed-in detail panels of the dayside magnetopause and the magnetotail, sketched on the Vlasiator run described by Palmroth et al. (2023). The magnetosphere proper is delimited by the magnetopause (purple surface). Magnetic field lines on the meridional plane are shown in black, with X (red, crosses; global FFJ line annotated) and O (blue, discs; marked with FTE and plasmoid) topologies marked as lines close to the meridional plane (and one X topology circling the near-Earth region with a northward magnetic field). The white and light-blue equatorial cross section follows the magnetic Equator and the tail current sheet, with the north–south magnetic field marked, respectively, as out of and into the plane.

Global magnetospheric simulations use global coordinate systems such as the Geocentric Solar Ecliptic (GSE) or Geocentric Solar Magnetospheric (GSM) systems, similarly to spacecraft observations. However, local magnetic structures are not necessarily aligned with the axis of the simulation coordinate system. In particular CSs, flux ropes and reconnection topologies can develop unconstrained by the domain boundaries and can be oriented in virtually any direction. In this situation, methods developed to analyse spacecraft data can inspire new ways to investigate data from global simulations. Palmroth et al. (2023) used justified assumptions about the geometry of the problem to obtain useful proxies from global coordinate systems, and the aim of this paper is to generalize the approach to a local coordinate system.

To review our region of interest, Fig. 1 shows a sketch of Earth's magnetosphere based on the Vlasiator simulation described by Palmroth et al. (2023). Solar wind flows in from the right, carrying with itself the interplanetary magnetic field (IMF), which is aligned purely southward in this simulation. As the solar wind encounters Earth's magnetic field, a bow shock (grey) and a heated magnetosheath (light brown) are formed (both are shown clipped along the meridional plane). The magnetopause (purple, clipped along the equatorial plane) delimits the magnetosphere proper: the volume in which Earth's magnetic field dominates the plasma behaviour. The solar wind driving draws the magnetosphere into a long tail, with northern and southern tail lobes having opposite magnetic fields and with the tail current sheet (CS, light blue) separating the lobes. Southward IMF driving is especially interesting as it allows for magnetic reconnection on the dayside magnetopause, where the opposing magnetic fields form X topologies (dark-red lines with crosses) and allow for the transfer of energy and mass across the magnetopause and a change in global magnetic topology. The global reconnection line described by Laitinen et al. (2006) is correspondingly shown as the closed, dark-red line marked as FFJ. The reconnected magnetic field lines form flux transfer events (FTEs or plasmoids) characterized by O topologies (blue lines with discs). Similarly, the opposing magnetic fields of the tail lobes may reconnect at the tail CS, forming plasmoids in the tail current sheet. Palmroth (2023) shows complex reconnection topologies that form in the tail, represented by the ellipsis-terminated red and blue lines in the tail of Fig. 1, with similar sketches of multiple reconnection lines and flux transfer events given on the dayside magnetopause.

Figure 2 Schematic drawing of the (a) MGA and (b) MDD eigenvectors $\widehat{\mathit{L}},\widehat{\mathit{M}}$ and $\widehat{\mathit{N}}$, corresponding to the largest, middle and minimum eigenvalues in a 1D CS, along the lines of Shi et al. (2019). For both eigensystems, the vector $\widehat{\mathit{L}}$ is well-defined, while the $\widehat{\mathit{M}}$ and $\widehat{\mathit{N}}$ vectors lie on the plane perpendicular to the $\widehat{\mathit{L}}$ vectors.

In this paper we introduce and discuss alternative, physically motivated local methods for finding X and O lines in the magnetic field of global 3D magnetospheric simulations by inspecting the magnetic field on a suitable local coordinate basis. To extend the global proxies used in Palmroth et al. (2023), we introduce a local coordinate system based on the work reviewed by Shi et al. (2019), with an approach similar to the one presented by Denton et al. (2010), and describe both contouring and cell-wise FOTE methods for finding null lines, both X and O, within a global magnetospheric plasma simulation based on the magnetic field and its Jacobian matrix.

1.1 Local coordinate systems and dimensionality

Let us consider two methods reviewed by Shi et al. (2019): minimum directional derivative (MDD) and minimum gradient analysis (MGA, analogous to minimum variance analysis), both defined using the Jacobian G=∇B of the magnetic field. Local, orthogonal coordinate directions can be obtained from the eigenbases of G^⊺G (MGA) and GG^⊺ (MDD). However, as these methods obtain the initial directions as eigenvectors that are defined up to a real scalar, the directions (after normalization) are defined up to a factor of ± 1.

MGA produces a set of basis vectors where the eigenvector that corresponds to the largest eigenvalue (λ₁) is aligned with the vector that has maximal variation (${\widehat{\mathit{L}}}_{\text{MGA}}$), the second one (λ₂) corresponds to the intermediate variation (${\widehat{\mathit{M}}}_{\text{MGA}}$), and the third (λ₃) corresponds to the least variation (${\widehat{\mathit{N}}}_{\text{MGA}}$). On the other hand, the MDD method produces a set of basis vectors where the eigenvector corresponding to the largest eigenvalue shows the direction of the displacement which produces the largest variation in B (${\widehat{\mathit{L}}}_{\text{MDD}}$). Correspondingly, for the eigenvectors corresponding to the intermediate and least eigenvalues we use ${\widehat{\mathit{M}}}_{\text{MDD}}$ and ${\widehat{\mathit{N}}}_{\text{MDD}}$. Denton et al. (2010) employ the MDD eigenbasis in their method, while we use both MDD and MGA eigenvectors, as described later in Sect. 2.1.

Both of these eigensystems have the same eigenvalues, but the eigenvectors differ and are not necessarily aligned with each other. We note here that, for a 1D structure, both of these eigensystems have only one well-defined eigenvector. For example, in a 1D Harris (1962)-type CS defined by

$$\begin{array}{}\text{(1)}& \mathit{B}\left(z\right)=B_{\mathrm{0}}\tanh \left(z\right)\widehat{\mathit{x}}\end{array}$$

for some constant B₀, the eigensystems are shown in Fig. 2, with ${\widehat{\mathit{L}}}_{\text{MGA}}$ being aligned with $\widehat{\mathit{x}}$ and ${\widehat{\mathit{L}}}_{\text{MDD}}$ being aligned with $\widehat{\mathit{z}}$.

MDD also gives us a way to define the local dimensionality of a structure (Rezeau et al., 2018) from the eigenvalues of the matrix GG^⊺, with, for example, a purely 1D structure varying only as a function of a single coordinate. The quantities D₁, D₂ and D₃ describe the one-, two- and three-dimensionalities of the magnetic field and are obtained from the ratios of square roots of the eigenvalues ${\mathit{\lambda }}_{\mathrm{1}}\ge {\mathit{\lambda }}_{\mathrm{2}}\ge {\mathit{\lambda }}_{\mathrm{3}}$ of GG^⊺.

Definitions for D₁, D₂ and D₃ are given below, modified to use the square roots of eigenvalues, as suggested by Shi et al. (2019) (as the square roots of the eigenvalues are then in units of T m⁻¹):

$$\begin{array}{}\text{(2)}& {\displaystyle }{D}_{\mathrm{1}}={\displaystyle \frac{\sqrt{{\mathit{\lambda }}_{\mathrm{1}}}-\sqrt{{\mathit{\lambda }}_{\mathrm{2}}}}{\sqrt{{\mathit{\lambda }}_{\mathrm{1}}}}},\text{(3)}& {\displaystyle }{D}_{\mathrm{2}}={\displaystyle \frac{\sqrt{{\mathit{\lambda }}_{\mathrm{2}}}-\sqrt{{\mathit{\lambda }}_{\mathrm{3}}}}{\sqrt{{\mathit{\lambda }}_{\mathrm{1}}}}},\text{(4)}& {\displaystyle }{D}_{\mathrm{3}}={\displaystyle \frac{\sqrt{{\mathit{\lambda }}_{\mathrm{3}}}}{\sqrt{{\mathit{\lambda }}_{\mathrm{1}}}}}.\end{array}$$

These quantities are defined to lie in the range [0,1], and their sum is 1. For D₁≈1, the magnetic field structure is primarily 1D, such as with a CS with B≈B(z) for a direction z normal to the CS. Correspondingly, for D₂≈1, the structure is primarily a function of two coordinates and so on. These measures allow us to quantify whether or not the local 2D treatment for neutral lines is well-founded, and we return to the topic of dimensionality in Sect. 4.3.

Figure 3 Neutral point classification in a well-defined local coordinate system. On the plane where B_L=0 and on the subset of that plane where B_N=0, the O points have $\partial B_{\mathrm{N}}/\partial L<\mathrm{0}$, and the X points have $\partial B_{\mathrm{N}}/\partial L>\mathrm{0}$. Note that this classification assumes a right-handed coordinate system.

1.2 Vlasiator

In this study, the methods discussed above will be applied to a global simulation of Earth's magnetosphere performed with the hybrid-Vlasov code Vlasiator. Vlasiator is a supercomputer-scale 6D (3D+3V – three spatial dimensions and three velocity-space dimensions) ion-hybrid plasma simulation (von Alfthan et al., 2014; Palmroth et al., 2018), solving the Vlasov equation for ion components and treating electrons as a charge-neutralizing, massless fluid. Previously, Vlasiator has been used to study reconnection in 5D (2D+3V) simulations for the global magnetosphere (e.g. Hoilijoki et al., 2017; Palmroth et al., 2017) and locally in the magnetopause in full 6D (Pfau-Kempf et al., 2020). Recently, computational advances (Ganse et al., 2023) have expanded the capabilities of the model to include fully 6D global simulations (Palmroth et al., 2023; Horaites et al., 2023; Grandin et al., 2023).

We use the global Vlasiator simulation described by Palmroth et al. (2023) as a prototype case for the discussed methods. The simulation is of Earth's global magnetosphere, with an approximate conducting-sphere inner boundary condition that does not include ionospheric coupling. The simulation uses the GSE coordinate system (centred on Earth, with the x axis pointing towards the Sun, while the z axis is normal to the ecliptic in the northward direction) and no dipole tilt; thus, the dipole moment is aligned with the z axis. The simulation is set up with a predefined mesh refinement with 1000 km maximum spatial resolution in the tail CS and on the dayside magnetopause, with the simulation domain spanning $[-\mathrm{111},\mathrm{50}]R_{\mathrm{E}}$ in the X direction and $[-\mathrm{58},\mathrm{58}]R_{\mathrm{E}}$ in the Y and Z directions, with driving provided by the fast (U_sw = 750 km s⁻¹) solar wind at the temperature of T_sw = 500 kK and the density of n_sw = 1 cm⁻³ at the inflow boundary and a purely southward magnetic field B_IMF = $(\mathrm{0},\mathrm{0},-\mathrm{5})$ nT. Other boundaries employ copy conditions. The initial condition for the simulation is comprised of the velocity distribution functions initialized to Maxwellian distributions, with moments at the inner boundary set to a density of 1 cm⁻³ and a temperature of 5 × 10⁵ K with zero velocity; these distributions are smoothly tapered to solar wind values until the geocentric distance of 15.7 R_E (Horaites et al., 2023).

Despite the solar wind and IMF being constant, the kinetic physics involved in the simulation still produce a remarkably dynamic environment. Palmroth et al. (2023) described the dynamics of tail reconnection in the Vlasiator run, showing the interplay of plasma instabilities and reconnection on the tail CS. Notably, they presented a complex picture of small, independent plasmoids being generated in a bursty fashion, followed by a tail-wide reconfiguration through the merging of smaller plasmoids and their reconnection lines. Notably, the simulation setup is an idealized case, and with more realistic solar wind and IMF inputs, including a flow-aligned IMF component, dipole tilt and magnetosphere–ionosphere coupling, we expect to see even more complex behaviour, including, for example, tilted CSs (see e.g. Shen et al., 2008).

2 Methods: applying the local coordinate systems to Vlasiator

2.1 Local LMN coordinates

If the eigenvalues are not well-separated, the directions obtained from MGA and MDD may be ambiguous (see Shi et al., 2019). For example, for the CS with dimensionality 1, MGA obtains the field-aligned direction (above and below the CS), while MDD obtains the normal direction in relation to the CS, with two other directions being potentially ambiguous. The directions obtained from MGA and MDD are, however, approximately orthogonal at the CS. Since we are especially interested in CS-like structures, we choose to use the unit vectors ${\widehat{\mathit{L}}}_{\text{MGA}}$ and ${\widehat{\mathit{L}}}_{\text{MDD}}$ (described in Sect. 1.1) as the basis on which to build our coordinate system.

The electric current density $\mathit{J}\propto \mathrm{\nabla }\times \mathit{B}$ can be used to define a right-handed coordinate system by selecting ${\widehat{\mathit{L}}}_{\text{MGA}}$ as a primary vector $\widehat{\mathit{L}}$, ${\widehat{\mathit{L}}}_{\text{MDD}}$ as a secondary vector and J as a tertiary vector so that $\widehat{\mathit{L}}\times \mathit{J}\parallel {\widehat{\mathit{L}}}_{\text{MDD}}$. For example, on an idealized CS (Eq. 1), we see that $\widehat{\mathit{L}}\parallel \widehat{\mathit{x}}$, $\widehat{\mathit{N}}\parallel \widehat{\mathit{z}}$, and $\mathit{J}\parallel \widehat{\mathit{L}}\times \widehat{\mathit{N}}$. $\widehat{\mathit{L}}$ and ${\widehat{\mathit{L}}}_{\text{MDD}}$ are exactly orthogonal in this case, as in Fig. 2.

To summarize the definition of our local LMN coordinate system $\widehat{\mathit{L}},\widehat{\mathit{M}},\widehat{\mathit{N}}$, we orthogonalize the set of basis vectors with the following:

1. $\widehat{\mathit{L}}={\widehat{\mathit{L}}}_{\text{MGA}}$.

2. $\widehat{\mathit{N}}={\widehat{\mathit{L}}}_{\text{MDD}}-({\widehat{\mathit{L}}}_{\text{MDD}}\cdot \widehat{\mathit{L}})\widehat{\mathit{L}}$, with $\widehat{\mathit{N}}$ multiplied by −1, if necessary so that $(\widehat{\mathit{N}}\times \widehat{\mathit{L}})\cdot \mathit{J}>\mathrm{0}$.

3. Lastly, $\widehat{\mathit{M}}=\widehat{\mathit{N}}\times \widehat{\mathit{L}}$.

The above results in a right-handed coordinate system with $\widehat{\mathit{M}}$ oriented in the same direction as J (but not necessarily parallel). This still only requires data from G and allows for a classification of a null line in relation to X and O in this LN plane via inspection of $\partial B_N/\partial L$ (see Fig. 3).

Figure 4 Absolute value of the dot products of the primary and secondary vectors $|{\widehat{\mathit{L}}}_{\text{MGA}}\cdot {\widehat{\mathit{L}}}_{\text{MDD}}|$ used in the construction of the LMN systems, shown on the planes Y=0 and Z=0. We note that the vectors are nearly orthogonal, especially in the vicinity of the magnetopause and the tail current sheet, where we wish to employ the orthogonalized coordinate system.

The sign of $\widehat{\mathit{L}}$ and $\widehat{\mathit{N}}$ in the global coordinates remains unfixed: multiplying both $\widehat{\mathit{L}}$ and $\widehat{\mathit{N}}$ by −1 does not affect the generality of the system. For now, this is left to be defined depending on the specific application.

Figure 4 shows the global magnetosphere as reproduced by the Vlasiator simulation, with Earth at the origin, the Sun on the lower right (in the direction of the positive x axis), and Earth's dipole field aligned to have the North Pole on the positive z axis. The figure displays the absolute value of the dot product ${\widehat{\mathit{L}}}_{\text{MGA}}\cdot {\widehat{\mathit{L}}}_{\text{MDD}}$, a measure of non-orthogonality for the primary and secondary vectors of our LMN coordinate system on the Y=0 and Z=0 planes. In regions such as the magnetopause (X≈10 R_E at the subsolar point) and the bow shock (X≈15 R_E), hosting boundaries and CSs that are expected to be 1D structures, we see that the non-orthogonality is low. The magnetotail lobes and parts of the magnetosheath display consistent alignment of the primary and secondary vectors, indicating that the most-varying component of B varies most in the same direction, e.g. when approaching the polar regions along the magnetic field. This shows that care must be taken when considering these LMN coordinates, especially in the lobes and inner magnetosphere. The transition between the CS and the lobes in particular displays a complex structure.

While the above local coordinate system is defined where there is some variation in B, such that $|{\widehat{\mathit{L}}}_{\text{MGA}}\cdot {\widehat{\mathit{L}}}_{\text{MDD}}|<\mathrm{1}$, and a non-zero electric current density J, our interest lies in reconnection-supporting CS structures, which naturally fulfil these conditions.

2.2 Contouring method

2.2.1 Finding reconnection-supporting CSs

In Palmroth et al. (2023), the CS centre was defined by the zero-crossing of the Earth-centred radial magnetic field component, B_r=0. This selects the crossover plane between the northern and southern lobes in the tail, with B_r=0 accounting for some radial component closer to the flanks. B_x=0 gives similar results in this case.

Figure 5 A Vlasiator tail CS (t = 1267 s). The CS centre plane is located by a contouring algorithm for B_L=0, forming a wavy sheet, which is coloured by the current density $\left|\mathit{J}\right|$. For reference, the current density is also shown on the $X=-\mathrm{15}{R}_{\mathrm{E}}$ plane, indicating a match with the actual location of the CS. The local LMN basis vectors are shown on the CS (red: $\widehat{\mathit{L}}$, green: $\widehat{\mathit{M}}$, blue: $\widehat{\mathit{N}}$) at arbitrary locations, and the B vectors are shown above and below the CS on a $Y=-\mathrm{9}{R}_{\mathrm{E}}$ slice (magenta).

The local coordinate system, as given above, allows us to define an LMN system at each spatial cell of the simulation, analogously to the definition at each time in a time series for spacecraft observations. Finding B_L=0 surfaces shows the centres of symmetric CSs where the locally most-varying component of B changes its sign. While such a definition for the centre of a general current sheet is not strictly correct, it is sufficient for our purposes. These B_L=0 surfaces we deem to be the most interesting on the basis of supporting reconnection. However, with contouring algorithms, regions of consistent orientation of the LMN system need to be selected. As an example, the LMN coordinates of a Vlasiator tail CS, as shown in Fig. 5, have been homogenized by requiring $\widehat{\mathit{L}}\cdot \widehat{\mathit{x}}>\mathrm{0}$ and have been focused onto the CS by requiring $\left|\mathit{J}\right|$ > 2 nA m⁻². Figure 5 demonstrates, firstly, correct extraction of the CS midplane and, secondly, details of the local LMN basis vectors: the $\widehat{\mathit{L}}$ vector tracks the orientation of the lobe fields, the $\widehat{\mathit{N}}$ vector tracks the CS normal, and the $\widehat{\mathit{M}}$ vector is in-plane and pointing in the general direction of the current.

Figure 6 (a) Null lines obtained from B_r=0 and B_z=0 on the Vlasiator tail CS – local zoom-in of the simulation within the domain presented in Fig. 5. (b) LMN-contoured null lines. Both sheets are coloured by positive or negative B_z or B_N, respectively, with the null lines coloured by $-{\partial }_r{B}_z$ and ∂_LB_N, respectively. The strong CS flapping leads the former to miss a plasmoid on the flap, as shown by the field lines (magenta).

2.2.2 Null lines

In Palmroth et al. (2023), the CS centre proxy B_r=0 serves as a basis for finding null lines on the sheet via the additional constraint of B_z=0. An example of these proxies and the flapping of the CS is shown in Fig. 6a. The sign of B_z serves as a proxy for field topology (given as sheet colour). We note that strong flapping of the CS causes the CS normal component to deviate strongly from the Z direction at points. Using the LMN basis allows us to take into consideration the local behaviour of the CS, improving upon the proxy model.

Finding the surfaces with B_L=0 and further intersecting those with the surfaces of B_N=0 yields a set of connected lines, as shown in Fig. 6b. For example, the contouring operator used by VisIt (Childs et al., 2012) finds these lines with topological connectivity. The left-hand side panel of Fig. 6 shows an example of the $B_r=B_z=\mathrm{0}$ contouring missing a plasmoid at a prominent CS flap, whereas the LMN method correctly finds the axis of the plasmoid. We note that the overall structure of the tail X and O lines is still well-described by the previous method outside of the region with large-amplitude flapping.

Above, and in general, we assume the existence of the normal component B_N for the detection of in-plane nulls. The potential degenerate case with B_N=0 over the entire CS (a tangential discontinuity) is hardly encountered in practice.

To use the LMN coordinate system and the components of magnetic field $B_{\mathrm{L}},B_{\mathrm{M}},B_{\mathrm{N}}$ for contouring requires the LMN coordinate systems to be consistently oriented: if the ambiguity in the sign of L and N coordinates happens to flip between neighbouring cells, the B_L and B_N components will naturally change signs as well, as illustrated in Fig. 7. This is a source of false detections, which can be fixed by re-orienting neighbouring coordinate systems. In the magnetotail, for example, the choice we make is for L to point sunward so that, if the constructed LMN system has an initially tailward L, we multiply the L and N vectors by −1 in these cells. Automatic determination of consistently oriented neighbourhood charts could be considered in a future update. However, it may not be possible to provide a consistent orientation for the LMN coordinates in an entire domain, such as a global magnetospheric simulation. As long as there is a consistent orientation, the simple contouring operation works for current-sheet-like structures and produces topologically connected neutral lines.

Figure 7 Schema of LN sign ambiguity and its effect. The signs of L and N vectors are not fixed by the construction and may vary arbitrarily. On the left is a correct evaluation of a B_N=0 (vertical dashed line) value between two points A and B, while on the right there is an incorrect evaluation of a B_N=0 value (vertical dashed line) between points B and C. The LN directions are locally correct but not necessarily consistent between neighbouring cells.

When there is no consistent orientation available, naïve contouring does not work in an off-the-shelf manner. For example, in the case of the upstream solar wind, where Vlasiator shows very little variation in IMF, the LMN coordinate system is dominated by numerical precision artefacts, leading to inconsistent orientation showing up as noise in these contouring methods.

2.3 Cell-wise FOTE method

An alternative that does not require choices on orientation is to use the local data of B and G to construct first the local coordinate system and secondly a linear expansion of the components B_L and B_N at the cell centre. This sidesteps the ambiguity in the sign of the LN directions. Solving for the planes of B_L=0 and B_N=0 and their intersection is then straightforward in this linear approximation: $B_L+\mathrm{\nabla }{B}_L\cdot \mathit{r}=\mathrm{0}$ is solved via an exact, single-step gradient descent along the negative gradient:

$$\begin{array}{}\text{(5)}& \mathit{r}=-{\displaystyle \frac{B_{\mathrm{L}}}{\left|\mathrm{\nabla }{B}_{\mathrm{L}}\right|}}\cdot {\widehat{\mathit{e}}}_{\mathrm{\nabla }{B}_{\mathrm{L}}},\end{array}$$

where ${\widehat{\mathit{e}}}_{\mathrm{\nabla }{B}_{\mathrm{L}}}$ is the unit vector in the direction of the gradient. This defines the plane B_L=0 with the normal direction ${\widehat{\mathit{e}}}_{\mathrm{\nabla }{B}_{\mathrm{L}}}$ and a point r relative to the cell centre. Similarly, we find B_N=0 and the intersection line of these two planes.

Figure 8 (a) Komar et al. (2013) 150° case reproduction: FOTE null lines (translucent blue cells) compared to analytical separator (magenta) and B field line stubs (light grey, from Y=0 plane) after slight propagation; note that the analytical separator is contained almost everywhere by the cells marked by the FOTE method. Bottom panels (b, c): LN slices at M=0 with magnetic field lines of in-plane components in blue, B_L=0 in black and B_N=0 in magenta. Background colour shows in-plane magnetic field magnitude. (b) This figure shows a detection of an X line in the vacuum superposition case (origin at cell marked by arrow b). (c) This figure shows an unexpected but correct detection of an X line in the vacuum superposition case (cell marked by arrow c).

The line where $B_{\mathrm{L}}=\mathrm{0}=B_{\mathrm{N}}$ would then again signify a null line, but the linear approximation can only be considered to be good within a small neighbourhood, that is, in the cell. It is then also straightforward to check whether or not this line intersects with the cell being analysed in the simulation coordinate system – if yes, we can consider the cell to contain a neutral line. We minimize a signed distance function (SDF; see e.g. Keinert et al., 2013) to the cell along the line, given by the definitions adapted from Korndörfer et al. (2015):

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{d}(\mathit{r},\mathit{b})& {\displaystyle }={\text{abs}}^{\mathrm{e}}\left(\mathit{r}\right)-\mathit{b},\text{(7)}& {\displaystyle }{f}_{\text{box}}(\mathit{r},\mathit{b})& {\displaystyle }=\left|\right|{\max }^{\mathrm{e}}(\mathit{d},\mathbf{0})\left|\right|+\max \left({\min }^{\mathrm{e}}\right(\mathit{d},\mathbf{0}\left)\right),\end{array}$$

where op^e denotes a component-wise operator (op) returning a vector of the same dimensionality, r is the query point, b is the upper corner of a rectangular cuboid centred at the origin and aligned with the coordinate system basis, and f_box is the SDF of the query point with respect to the cuboid. The SDF describes the distance of a point to the surface of an object, with positive values indicating that the point lies outside the object and negative values indicating that the point lies inside. Therefore, an SDF gives a continuous measure of the distance of the line with the cell, with negative values indicating the neutral line intersecting the cell and positive values indicating no intersection. This allows some evaluation and acceptance of near misses. A discriminating value of SDF≈0.25…0.36 is found to produce uninterrupted null lines in the tail, as expected based on field topology. Figure 8a shows examples of the FOTE method in use for Vlasiator data.

Further, we may evaluate the derivative ∂_LB_N at the cell and again classify the possible neutral line with it. Eigenvalues of G might as well be used for this purpose, as for degenerate nulls in Parnell et al. (1996), but we opt to use a straight derivative for a clear physical meaning. The main disadvantage for now is that this method does not provide a ready-made topology, only sets of cells that are considered to contain a segment of an X or O line.

Figure 9 Komar et al. (2013) 150° case reproduction: contoured null lines (blue tubes) compared to analytical separator (magenta) and B field line stubs (light grey, from Y=0 plane) along the lines of Hu et al. (2009) and Komar et al. (2013) after slight propagation.

3 Validation with an ideal case

We take the southward-IMF initial configurations of Hu et al. (2009) and Komar et al. (2013) for a test case that has null lines; the mentioned studies also provide analytical solutions for the magnetosphere using the methods by Yeh (1976). For the magnetic field configuration, we use as the constant IMF background magnetic field ${\mathit{B}}_{\text{IMF}}=(\mathrm{0},\mathrm{32.5},-\mathrm{56.3})$ nT for a field strength of 65 nT at a 150° clock angle and an unscaled terrestrial dipole aligned with the z axis. The domain size is $(\pm \mathrm{2}\times {\mathrm{10}}^{\mathrm{5}}\mathrm{km}{)}^{\mathrm{3}}$, and the maximum spatial resolution is 3125 km. The current-free initial condition is not suitable to apply our method to since we use the current density vector to define the handedness of our coordinate basis. However, we choose to propagate the initial condition for a small time for the global plasma dynamics (7.048 ms). The background plasma is initialized at a density of 1 cm⁻³ as a radial taper of temperature and velocity from the inner boundary values of 5 MK and zero velocity to the solar wind values of 0.5 MK and 750 km s⁻¹ to break the exact current-free condition. This enables the use of our methods to extract the X line in the case with IMF clock angle θ = 150° (mostly southward IMF). To guard against current-free edge cases, we can use, as a fall-back case, the MGA maximum and medium variance coordinates for L and N directions instead.

Figure 8a shows a comparison between the FOTE method and the analytical solution. The cells detected by the FOTE method enclose the analytical separator nearly everywhere, but they show, in addition, a few extra null line segments further away. Figure 8b and c show details of the X-line detection by plotting L–N cross sections of the in-plane nulls. The origins of these cross sections are placed in the cells with FOTE-detected in-plane nulls, with the cross sections taken at M=0. The LMN basis is the local basis from the detection cell, extrapolated to the immediate neighbourhood. The slices show projections of the magnetic field lines in the LN plane and zero-crossing lines for B_N and B_L. The cross sections confirm these cells do contain in-plane null lines, as per our definition, but these are not described by the analytical solution for the separator. Figure 8 shows that these are indeed valid detections with an X topology. Notably, these detections are not present when inspecting the initial state with an MGA-only coordinate system and are likely a result of the initial velocity and temperature gradients causing minor perturbations.

Figure 10 (a) The tail CS from Vlasiator simulation (EGI, t = 1267 s). Similar format to Fig. 6 LMN plot, with contoured X and O lines. (b) FOTE detections of X and O topologies (outlined cubic cells) compared with the contoured detections of panel (a) (black curves). Panels (c) and (d) show details from points marked by c and d, specifically a slice of the local L–N plane with in-plane magnetic field lines from streamline plotting of the in-plane magnetic field components and with zero contours of these components being annotated. The slices show that point c is correctly identified as an O topology, and point d is correctly identified as an X topology. The background colour is the strength of the in-plane magnetic field.

Figure 9 shows a comparison between the contouring method and the analytical solution. Here, the LMN coordinate system has been regularized so that the L direction is chosen to fulfil $\widehat{\mathit{L}}\cdot \widehat{\mathit{z}}>\mathrm{0}$. For this particular case, the LMN contouring of the superposition field produces prominent conical artefacts (not shown) where the coordinate system is not oriented consistently with this particular choice, and the data to be contoured are filtered to lie within one cell of a FOTE null line as a way to choose compact neighbourhoods of regularized orientation. Less stringent filtering can be used for other, more local cases (see, e.g. Sect. 4.1). Good agreement with the analytical solution is also seen in this case.

4 Results in Vlasiator

4.1 CS null lines in detail

Figure 10a and b show a section of the Vlasiator tail CS (a subset of Fig. 5). Firstly, panel (a) shows the in-plane null lines as found by the contouring method, coloured with ∂_LB_N (with red colours signifying X topologies and blue lines signifying O topologies). Evidently, the method struggles at the site of strong current sheet flapping. The CS is coloured by B_N, positive values are coloured blue, and negative values are coloured beige, extending the north–south classification in Fig. 6a to a flapping CS. The LMN coordinate bases are shown at some locations on the sheet as red–blue–green arrows. Further, two in-plane null line locations marked as (c) and (d) are chosen as examples (noting that these coincide with the flap seen in Fig. 6a), and for these locations the local LN plane magnetic field is plotted below in the corresponding panels (c) and (d), showing that these null lines are indeed correctly found and classified.

To validate the observed null lines with both methods in the case of the tail CS, Fig. 10b shows both contouring and FOTE detections at work, with the contoured in-plane null lines shown as black curves and the FOTE-detected null-line containing cells as wireframes coloured with ∂_LB_N. The methods agree well, even if the FOTE method may have some stray detections. These strays are found to be a minor issue, showing up on the current sheet as isolated cells flagged as containing null lines, while the true detections form coherent structures.

Figure 11 Dawn flank of the magnetosphere, with FOTE detections as outlined cells, using the β^∗ parameter on the vertical slice and on a magnetopause proxy at ${\mathit{\beta }}^{\ast }=\mathrm{0.5}$.

We may note that the local coordinate systems spanned by our LMN basis describe well the in-plane X and O lines in the simulation. However, this leads to the observation that the X and O line axes might not be aligned with the $\widehat{\mathit{M}}$ direction, as also observed by Pathak et al. (2022) with MMS data. Figure 10 demonstrates this in the Vlasiator tail CS, with X and O lines necessarily breaking their alignment from the $\widehat{\mathit{M}}$ direction when forming into loops, for example at $X=-\mathrm{17}{R}_{\mathrm{E}}$, Y=3 R_E. The effect of this misalignment on reconnection in Vlasiator simulations is an ongoing subject of study.

4.2 Null lines in global Vlasiator 6D magnetosphere

To extend our analysis to the global Vlasiator magnetosphere, we use the generic FOTE method that is not limited by construction of consistently oriented neighbourhoods. The FOTE method reveals a sinuous and alternating neutral line structure on the magnetopause flanks, shown in Fig. 11. The O lines are marked with blue cells signifying dayside FTEs extending onto the flanks as a part of the sinuous pattern. The figure uses a magnetopause proxy from the β^∗ parameter of Xu et al. (2016), defined as

$$\begin{array}{}\text{(8)}& {\mathit{\beta }}^{\ast }={\displaystyle \frac{P_{\text{thermal}}+P_{\text{dynamic}}}{P_B}}.\end{array}$$

Regions with small values of $\left|{\partial }_{\mathrm{L}}{B}_{\mathrm{N}}\right|$ (grey) should be inspected for their validity, in general, as these signify low variability in the magnetic field. These regions may be susceptible to local, insignificant perturbations fulfilling the FOTE criterion. The alternating pattern on the magnetopause at X<20 still suggests that these are indeed real magnetic structures, with alternating positive and negative B_N regions. We note, however, that the behaviour of the magnetopause on the flank in the present simulation is not yet studied in detail, and the shown time state (1267 s) may still contain initialization artefacts deep in the tail.

Figure 12 Overview of local MDD dimensionality, shown on the planes Y=0, Z=0. The dimensionality values $D=(D_{\mathrm{1}},D_{\mathrm{2}},D_{\mathrm{3}})\in [\mathrm{0},\mathrm{1}{]}^{\mathrm{3}}$ are mapped to RGB colour channels via $\text{rgb}=D/\mathrm{2}+(\mathrm{0.5},\mathrm{0.5},\mathrm{0.5})$ to obtain a three-channel representation of the triplet. See the ternary diagram at the top left for reference: the more red there is, the more 1D the local magnetic field structure is; green dominance signifies two-dimensionality, and blue dominance signifies three-dimensionality.

4.3 Dimensionality of magnetic structures in Vlasiator

The applicability of the LMN coordinate system should be inspected in terms of dimensionality of the magnetic structures to ensure that the structures inspected are locally quasi-2D. Figure 12 shows an overview of the dimensionality measures D₁, D₂ and D₃ (Eqs. 2–4, described in Sect. 1.1) applied to a Vlasiator 6D simulation. Red colour corresponds to regions dominated by 1D variation (sheet-like structures), green corresponds to regions dominated by 2D structures (e.g. plasmoids and flux ropes), and blue corresponds to regions dominated by 3D structures. The upstream appears to be noisy as the magnetic field is nearly constant. The regions of interest for reconnection are naturally found in CSs, that is, mostly 1D structures with some embedded 2D features. These can be used to constrain the local basis analysis to compatible regions.

Figure 13 Local dimensionality of a Vlasiator simulation obtained via the MDD analysis; ternary histogram over most of the magnetosphere. Each cell outside of the inner boundary with r>5 R_E and constrained to −40$R_{\mathrm{E}}<X<\mathrm{20}{R}_{\mathrm{E}}$ and ± 40 R_E in Y and Z is shown here. The lack of completely 3D structures indicates that the in-plane approximation is usable.

An interesting finding in the dimensionality parameters in the 6D Vlasiator run can be observed in the lack of dominantly 3D structures, as shown in Fig. 13. Magnetic structures exhibiting D₃ dominance appear to be limited: the occurrence of cells is increasingly rare as D₃→1. The prominent collection of cells having D₃=D₁ is found to be dominated by cells in the lobe regions of the magnetosphere. Whether or not these features are related to the divergence-free condition ∇⋅B, intrinsic properties of the magnetic field being convected by plasma, the boundary condition or other constraints remains to be studied.

5 Discussion and conclusions

Here, we have presented two local methods to acquire X and O topologies of the magnetic field from plasma simulations of Earth's magnetosphere, enabling fast and efficient identification of possible reconnection sites and flux transfer events. The algorithms work on O(N) complexity, i.e. operating once per grid cell (and possibly using a contouring algorithm of similar complexity) without requiring field line or separator tracing such as that used by Haynes and Parnell (2010). These methods use solely the magnetic field and its Jacobian, and so the found X lines may or may not be actively reconnecting – other methods need to be used to characterize the reconnection activity at the found sites, but the present methods have the advantage of being constrained to a local neighbourhood and only requiring knowledge about the magnetic field and its derivatives.

Using the dimensionality measures introduced by Rezeau et al. (2018) shows well that the magnetopause, the bow shock and the tail CS are essentially 1D structures, with intermittent, embedded 2D features. As we have a purely southward-IMF driving in the prototype case, this is consistent with Zeiler et al. (2002) observing that non-guide-field reconnection is mostly 2D. Future simulations including an IMF B_y component will provide stronger guide field reconnection and possibly different dimensionality of detected structures.

The construction of the local coordinate system involves some choices. The first in-plane direction $\widehat{\mathit{L}}$ is given by the local MVA analogue (MGA). The normal direction $\widehat{\mathit{N}}$ is chosen to be the primary eigenvector of the MDD system, orthogonalized with respect to $\widehat{\mathit{L}}$, to ensure good behaviour in nearly 1D CSs; the vector ∇B_L could be a conceivable option to supplant the secondary MDD vector, but this vector is not readily available from the magnetic field and its Jacobian.

For the contouring method, the limitation of requiring consistently oriented local coordinates is a slight issue. It could be mitigated by automated construction of local neighbourhood charts of consistent orientation, as well as automatically finding the extents where the local coordinate charts are sensible. As presented here, the contouring method relies on manual restriction of neighbourhoods for analysis. The FOTE method is given as a generalized alternative, operating cell-wise in the full domain, with the caveat that, in contrast to the contouring method, the FOTE method does not produce topological connectivity for the in-plane null lines; that is, it can be used to mark cells containing null lines but not which neighbouring cells those lines connect to. In contrast, the contouring method outputs line topologies from the isocontours of $B_{\mathrm{L}}=\mathrm{0}=B_{\mathrm{N}}$.

In comparison to topological methods, such as the magnetic skeleton method of Haynes and Parnell (2010), the present method does not require field line tracing or knowledge on boundary flows like the method proposed by Titov et al. (2009). On the other hand, Lapenta (2021) proposed a local Lorentz transformation method using both magnetic and electric fields, the latter of which the present method does not require. The present method can acquire useful null line objects via contouring, presenting something of a middle ground between local and topological analyses, adding another tool for studies of reconnection.

The presented FOTE and contouring methods are fundamental to investigate magnetic structures in a physically meaningful coordinate system, supporting reconnection studies in quasi-2D current sheets. The quantification of dimensionalities of magnetic structures may also provide new insights into plasma processes: for example, the magnetospheric tail lobes were found to be characterized by D₁≈D₃.