1. Introduction
Magnetic reconnection is responsible for explosive energy release in space and laboratory plasmas (Yamada, Kulsrud & Ji Reference Yamada, Kulsrud and Ji2010; Gonzalez & Parker Reference Gonzalez and Parker2016). By converting energy stored in the magnetic field to plasma energy, reconnection is able to accelerate electrons to high energies (e.g. Øieroset et al. Reference Øieroset, Lin, Phan, Larson and Bale2002; Drake et al. Reference Drake, Shay, Thongthai and Swisdak2005; Fu et al. Reference Fu, Khotyaintsev, Vaivads, Retinò and André2013; Oka et al. Reference Oka, Birn, Egedal, Guo, Ergun, Turner, Khotyaintsev, Hwang, Cohen and Drake2023). The energisation of electrons during magnetic reconnection has primarily been studied in the Lagrangian guiding centre framework (Northrop Reference Northrop1963), where the electron motion around the guiding centre is reduced to a magnetic moment
$\mu =m_ev_\perp ^2/(2B)$
, which is assumed to be conserved along collisionless orbits (e.g. Drake et al. Reference Drake, Swisdak, Che and Shay2006a
; Dahlin et al. Reference Dahlin, Drake and Swisdak2014,Reference Dahlin, Drake and Swisdak2015). Here,
$m_e$
is the electron mass and
$v_\perp$
is the velocity component perpendicular to the local magnetic field
$\boldsymbol{B}$
that has the magnitude
$B=|\boldsymbol{B}|$
. Under this assumption, one can derive a differential equation for the evolution of a single guiding centre’s energy
$\epsilon$
, as
\begin{equation} \frac {{{\rm d}}\epsilon }{{\rm d}t}=\mu \frac {\partial B}{\partial t} + q_e\left (v_\parallel \boldsymbol{b}+\boldsymbol{v}_g+\boldsymbol{v}_c\right )\cdot \boldsymbol{E}, \end{equation}
where
$\boldsymbol{b}=\boldsymbol{B}/B$
,
$q_e$
is the electron charge,
$v_\parallel$
is the velocity component parallel to the magnetic field,
$\boldsymbol{v}_g$
and
$\boldsymbol{v}_c$
are respectively the magnetic field gradient and curvature drifts, and
$\boldsymbol{E}$
is the electric field (e.g. Northrop Reference Northrop1963; Dahlin, Drake & Swisdak Reference Dahlin, Drake and Swisdak2014). The standard procedure is then to sum (1.1) over an ensemble of guiding centres in a local region to derive a fluid-like equation describing the evolution of the total electron kinetic energy density
$\mathcal{E}_{\textrm{LGCM}}$
,
\begin{equation} \frac {{\rm d}\mathcal{E}_{\textrm{LGCM}}}{{\rm d}t}= J_{e\parallel }E_\parallel + \frac {p_{e\perp }}{B}\left (\frac {\partial B}{\partial t}+\boldsymbol{u}_{E}\cdot \nabla B\right )+\left (p_{e\parallel }+m_en_eu_{e\parallel }^2\right )\boldsymbol{u}_E\cdot \boldsymbol{k}, \end{equation}
where
$p_{e\perp }$
and
$p_{e\parallel }$
are respectively the electron pressure components perpendicular and parallel to the local
$\boldsymbol{B}$
,
$J_{e\parallel }$
is the parallel component of the electron current density
$\boldsymbol{J}_e=q_en_e\boldsymbol{u}_e$
, where
$n_e=\int d^3v f_e$
is the electron density, with the electron velocity distribution function (VDF)
$f_e=f_e(\boldsymbol{r},\boldsymbol{v},t)$
;
$\boldsymbol{u}_e=n_e^{-1}\int d^3v \boldsymbol{v} f_e$
is the bulk electron velocity;
$u_{e\parallel }$
is the component of
$\boldsymbol{u}$
parallel to
$\boldsymbol{B}$
;
$\boldsymbol{u}_E=\boldsymbol{E}\times \boldsymbol{B}/B^2$
is the E-cross-B-drift; and
$\boldsymbol{k}=\boldsymbol{b}\cdot \nabla \boldsymbol{b}$
is the magnetic curvature. We will use the notation from TenBarge, Juno & Howes (Reference TenBarge, Juno and Howes2024) to refer to the first, second and third terms on the right-hand-side as
$W_\parallel$
,
$W_{{beta-gradB}}$
and
$W_{{curv0}}$
, respectively. The first term is related to acceleration by parallel electric fields, the second to betatron acceleration and the third to Fermi acceleration due to the curvature drift. We will refer to this model as the Lagrangian guiding centre model (LGCM). Equation (1.2) has been used extensively in both spacecraft data analysis (e.g. Eriksson et al. Reference Eriksson, Vaivads, Alm, Graham, Khotyaintsev and André2020; Zhong et al. Reference Zhong2020; Jiang et al. Reference Jiang2021) and numerical studies (e.g. Dahlin et al. Reference Dahlin, Drake and Swisdak2014, Reference Dahlin, Drake and Swisdak2015, Reference Dahlin, Drake and Swisdak2016; Li et al. Reference Li, Guo, Li and Li2015) to quantify electron energisation during magnetic reconnection.
However, recent work by TenBarge et al. (Reference TenBarge, Juno and Howes2024) has shown that there are some often overlooked nuances in (1.2). In particular, it is derived in a Lagrangian framework, but often directly applied to an Eulerian grid. In their study, TenBarge et al. (Reference TenBarge, Juno and Howes2024) used a continuum Vlasov–Maxwell solver to simulate single x-point reconnection and investigate the Lagrangian model of electron energisation. Their results showed that while the right-hand side of (1.2) provides a qualitative approximation of the electron fluid energisation as quantified by
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
when summed over the whole simulation domain (differences up to around
$20\,\%$
), it fails to reproduce
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
on a local, point-by-point, level. The authors then derived a complete Eulerian fluid decomposition of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
, and expressed it in several different ways, one of which is
\begin{align} \begin{split} \boldsymbol{J}_e\cdot \boldsymbol{E} &= J_{e\parallel }E_\parallel + \boldsymbol{u}_E\cdot \nabla p_{e\perp }+(p_{e\parallel }-p_{e\perp })\boldsymbol{u}_E\cdot \boldsymbol{k}+\boldsymbol{u}_E\cdot \left (\nabla \cdot \boldsymbol{\Pi }_e^a\right )+m_en_e\boldsymbol{u}_E\cdot \frac {{\rm d}\boldsymbol{u}_e}{{\rm d}t} \\ &:= W_\parallel + W_{\textrm{diam}} + W_{\textrm{curv1}} + W_{\textrm{agyro}} + W_{\textrm{pol}}, \end{split} \end{align}
where
$\boldsymbol{\Pi }^a_e = \mathbf{P}_e - p_{e\perp }\mathbf{I} - (p_{e\parallel } - p_{e\perp }) \boldsymbol{b}\boldsymbol{b}$
is the (traceless) agyrotropic part of the electron pressure tensor defined as
$\mathbf{P}_e=m_e\int d^3v (\boldsymbol{v}-\boldsymbol{u}_e)(\boldsymbol{v}-\boldsymbol{u}_e) f_e$
, and
$\mathbf{I}$
is the identity matrix. Here,
$W_{\textrm{diam}}$
is related to the diamagnetic drift,
$W_{\textrm{curv1}}$
to the Eulerian curvature drift,
$W_{\textrm{agyro}}$
to the agyrotropic drift and
$W_{\textrm{pol}}$
to the polarisation drift. We will refer to this model as the Eulerian fluid model (EFM). Equation (1.3) was derived without approximations, and it was found to accurately describe
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
on a global and local level. For the case studied by TenBarge et al. (Reference TenBarge, Juno and Howes2024), the authors found that the energisation term related to the diamagnetic drift (
$W_{\textrm{diam}}$
) was, by far, the most dominant on system scales, while all terms could have significant local contributions to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. We note in passing that other Eulerian decompositions of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
, corresponding to different groupings of the different terms, have previously been studied, such as the fluid compression/shear decomposition by Li et al. (Reference Li, Guo, Li and Birn2018). However, in that specific case, the electrons are assumed to be well magnetised such that the agyrotropic contribution is neglected.
The work of TenBarge et al. (Reference TenBarge, Juno and Howes2024) raises several interesting questions. First, as previously mentioned, it is well known that magnetic structures such as magnetic islands play an important role in accelerating electrons to high energies through Fermi-like processes (Drake et al. Reference Drake, Swisdak, Che and Shay2006a ). Such islands are formed when reconnection is triggered at multiple locations along the current sheet (e.g. Karimabadi, Daughton & Quest Reference Karimabadi, Daughton and Quest2005; Drake et al. Reference Drake, Swisdak, Schoeffler, Rogers and Kobayashi2006b ). Therefore, it is of interest to study the EFM also in the multiple x-point case, and compare it to the LGCM. On a related note, we know that the relative importance of the LGCM terms depends on the guide field strength (Dahlin, Drake & Swisdak Reference Dahlin, Drake and Swisdak2016), so it is natural to ask how the EFM terms depend on the guide field. We thus pose the question: which of the EFM terms are most important when reconnection leads to the formation of magnetic islands, and is there a guide field dependence?
Second, the inability of the LGCM to reproduce
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
begs the question what is the reason for the erroneous LGCM prediction? Is it related to a false assumption of
$\mu$
-conservation, is it due to the Eulerian/Lagrangian framework inconsistency discussed by TenBarge et al. (Reference TenBarge, Juno and Howes2024) or is it something else? Related to this, what happens to the electrons in the regions where the LGCM and EFM predict significantly different
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
?
Finally, one important feature of the LGCM is that each
$W$
term can be interpreted as having a different effect on the electron velocity distribution function. The betatron term (
$W_{\textrm{beta-gradB}}$
) is associated with perpendicular energisation, while the curvature (or Fermi) term (
$W_{\textrm{curv0}}$
) is associated with parallel energisation and ultimately the formation of power-law tails if repeated interactions are possible. There is a pedagogical value in these types of simple interpretations as they can aid in developing physical intuition. With this in mind, it is tempting to ask the question: can we make analogous interpretations of the fluid terms in (1.3)? That is, can we associate local measurements of each EFM
$W$
-term with the generation of specific VDF features?
In the present paper, we answer the aforementioned questions using particle-in-cell (PIC) simulations of multi-x-point reconnection. We find that
$W_{\textrm{diam}}$
is the dominant EFM contributor to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
on large scales for low and intermediate guide fields,
$B_g\in \{0,0.2\}$
, where
$B_g$
is normalised to the reconnecting field component. For large guide fields (
$B_g = 1$
),
$W_\parallel$
instead becomes the largest contributor (albeit still comparable to
$W_{\textrm{diam}}$
). Locally,
$W_{\textrm{agyro}}$
and
$W_{\textrm{curv1}}$
can contribute significantly to electron energisation, particularly in the lower
$B_g$
regime. We show that the assumption of
$\mu$
-conservation made by the LGCM is invalid in the low
$B_g$
case near the x-points and at the centre of the current sheets. This causes significant errors in the LGCM energisation estimate. Moreover, we find that local deviations from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
are produced by the LGCM even when
$\mu$
is conserved. This is due to the fact that the LGCM, unlike the EFM, does not describe the energy flux through the grid. Finally, we argue that we cannot easily associate local measurements of the different EFM energisation terms with the formation of features in the electron VDF.
2. Numerical set-up
In this study, we use the OSIRIS PIC-code (Fonseca et al. Reference Fonseca2002) to study electron energisation during multi-x-point magnetic reconnection in two spatial and three velocity dimensions. We use a Harris current sheet (Harris Reference Harris1962) set-up with periodic boundary conditions in
$x$
and perfectly conducting boundaries in
$y$
. We let the reconnection process grow from numerical noise. For all runs, we use an unperturbed magnetic field profile
$B_x(y)=B_\infty \tanh {(y/\lambda )}$
with an asymptotic magnetic field strength
$B_\infty$
corresponding to
$\omega _{ce\infty }/\omega _{pe0}=1/3$
. Here, the electron cyclotron frequency
$\omega _{ce}=|q_e| B/m_e$
is evaluated at
$B=B_\infty$
, and
$\omega _{pe0}=\sqrt {q_e^2 n_0/(m_e\epsilon _0)}$
is the reference electron plasma frequency, using the Harris sheet density perturbation
$n_0$
. The complete density profile is
$n_e=n_0\text{sech}^2(y/\lambda )+n_\infty$
, where we use
$n_\infty /n_0=0.2$
and a current sheet thickness of
$\lambda /d_{e0}=1.25$
, where
$d_{e0}=c/\omega _{pe0}$
is the electron inertial length. The narrow current sheet thickness was chosen to be consistent with Dahlin et al. (Reference Dahlin, Drake and Swisdak2014) and for reconnection to trigger quickly. We use two plasma components, one electron population and one ion population, with a uniform ion-to-electron temperature ratio of
$T_i/T_e=1$
.
As our baseline numerical resolution, we use a time step
$0.15\omega _{pe0}^{-1}$
and cell size
$\Delta x=\Delta y=0.25d_{e0}$
. The size of the simulation domain is
$L_x\times L_y=51.2d_{i0}\times 12.8d_{i0}$
, where
$d_{i0}=d_{e0}\sqrt {m_i/m_e}$
is the ion inertial length. Each cell is initialised with 625 macroparticles of each species. We use guide field values of
$B_g/B_{\infty }\in \{0,0.2,1\}$
, and low mass ratios
$m_i/m_e\in \{25,100\}$
similar to Dahlin et al. (Reference Dahlin, Drake and Swisdak2014) and TenBarge et al. (Reference TenBarge, Juno and Howes2024). We note that neither of the aforementioned studies investigated the
$B_g=0$
case, to avoid the problem of demagnetised electrons breaking the LGCM assumption of
$\mu$
conservation.
Henceforth, unless otherwise specified (e.g. in figure axes with explicit normalisation), we use the following normalisations in our calculations. Time is normalised by
$\omega _{pe0}^{-1}$
, lengths by
$d_{e0}$
,
$\boldsymbol{E}$
and
$c\boldsymbol{B}$
by
$m_ec^2/(ed_{e0})$
, where
$c$
is the speed of light and
$e$
is the elementary charge. Compound quantities such as
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
are normalised accordingly.
3. Results
3.1. Global energisation and parametric dependence
We will start by investigating the energisation models on large scales. The temporal evolution of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
is shown in figure 1, for the
$m_i/m_e=25$
,
$B_g=0$
run. The evolution is qualitatively similar for all runs. As reconnection starts and magnetic islands start to form (figure 1
a), patches of positive
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
are found at the x-points and the edges of the magnetic islands. The total contribution of the different energisation terms in (1.2) and (1.3) over the whole simulation domain for this time are shown in figures 2(a) and 2(b), marked by the first vertical line. We find that the curvature/Fermi term
$W_{\textrm{curv0}}$
is the dominant contributor to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
in the LGCM, while the diamagnetic term
$W_{\textrm{diam}}$
dominates in the EFM. The dominance of these two terms is consistent with the previous studies by Dahlin et al. (Reference Dahlin, Drake and Swisdak2014) and TenBarge et al. (Reference TenBarge, Juno and Howes2024).
Figure 1. Snapshots of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
for the
$m_i/m_e=25$
,
$B_g=0$
run at three different times
$t\omega _{ci\infty }\in \{8.7,11.3,16.0\}$
, where
$\omega _{ci\infty }=\omega _{ce\infty }m_e/m_i$
. The black lines are contours of the vector potential component
$A_z$
.
As the system develops further (figure 1
b), the islands contract, coalesce and move with the embedding plasma flow. Locally,
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
varies depending on the local island dynamics. The largest values of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
are found in narrow regions at the magnetic island edges. These regions are associated with local contraction of the island, and Fermi-like processes are expected to occur there (Drake et al. Reference Drake, Swisdak, Che and Shay2006a
; Dahlin et al. Reference Dahlin, Drake and Swisdak2014). Similar regions of negative
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
correspond to local expansion of the island. Consequently, islands that are convecting have a bipolar
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
signature, where the leading edge (locally expanding) is negative, and the trailing edge (locally contracting) is positive. At
$t\omega _{ci\infty }=11.3$
(second vertical lines in figures 2
a and 2
b), the LGCM finds that the Fermi-related
$W_{\textrm{curv0}}$
is indeed large, while
$W_\parallel$
gives a slight negative contribution to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
together with
$W_{\textrm{beta-gradB}}$
. In the case of the EFM,
$W_{\textrm{diam}}$
is still large, but now there are also significant net contributions from
$W_{\textrm{curv1}}\gt 0$
,
$W_\parallel \lt 0$
and
$W_{\textrm{agyro}}\lt 0$
, showing that multiple fluid drifts contribute to the net energisation. The same general picture remains valid as the islands continue to grow and merge (figure 1
c).
Figure 2. Time evolution of energisation terms summed over the whole domain (arbitrary units) for different
$m_i/m_e\in \{25,100\}$
and
$B_g/B_\infty \in \{0,0.2,1\}$
as labelled in the panels. The left column shows the guiding centre model results from (1.2) and the right column the corresponding fluid model results from (1.3). The three vertical dashed lines in panels (a) and (b) correspond to the times of the panels in figure 1. Note that we plot a shorter time interval for the
$m_i/m_e=100$
simulations in panels (g–j). The y-axis of panel (g) has been limited to better show the deviation between
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
and the LGCM sum. All data have been smoothed using a moving average with a window of
$\pm 1\omega _{ci\infty }^{-1}$
to reduce noise.
Looking at the temporal evolution of the
$W$
-terms in figures 2(a) and 2(b), we find that the LGCM generally provides a qualitative approximation of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
(compare the purple and the dashed black curves), particularly during the time interval
$t\omega _{ci\infty }\in (7,18)$
when
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
is large. However, at times, the relative difference between the LGCM sum and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
exceeds
$80\ \%$
. Integrating over the duration of the simulation reduces the relative difference to approximately
$10\ \%$
. In the EFM case, the agreement between the sum of the
$W$
-terms and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
is almost perfect. These results are consistent with those from TenBarge et al. (Reference TenBarge, Juno and Howes2024) in the single x-point case.
The remaining panels in figure 2 show the results of the same analysis applied to four other simulations with varying mass-ratio and guide field. Figures 2(c) and 2(d) show the case with moderate guide field
$B_g=0.2$
, at the same mass ratio
$m_i/m_e=25$
. There, the LGCM performs somewhat better than in the
$B_g=0$
,
$m_i/m_e=25$
case discussed previously, and the betatron
$W_{\textrm{beta-gradB}}$
term is slightly less important than before. In the EFM case, we find that the increased
$B_g$
almost removes the curvature contribution entirely and the agyrotropic contribution also becomes less important. This is expected since including an out-of-plane magnetic field component reduces the magnetic curvature across the current sheet and, consequently, the curvature drift. In addition, the guide field contributes to keeping the electrons magnetised (e.g. Swisdak et al. Reference Swisdak, Drake, Shay and McIlhargey2005; Scudder et al. Reference Scudder, Karimabadi, Daughton and Roytershteyn2015), reducing the agyrotropic pressure carried by demagnetised electrons. Increasing the guide field to
$B_g=1$
(figures 2
e and 2
f) further improves the accuracy of the LGCM to a very good agreement with
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. The relative contribution of
$W_{\textrm{curv0}}$
decreases, while
$W_\parallel$
, which for lower
$B_g$
was a negative contributor to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
, becomes large and positive. Analogously, in the EFM (figure 2
f),
$W_{\textrm{diam}}$
decreases and
$W_\parallel$
becomes the largest energisation term, albeit still comparable to
$W_{\textrm{diam}}$
. The increased importance of
$W_\parallel$
in the LGCM for increasing guide fields was previously reported by Dahlin et al. (Reference Dahlin, Drake and Swisdak2016).
Increasing the mass ratio to
$m_i/m_e=100$
does not change the results significantly. In the
$B_g=0$
case, the LGCM results (figure 2
g) are qualitatively very similar to the
$m_i/m_e=25$
case (figure 2
a), with
$W_{\textrm{curv0}}$
being the dominant positive contributor, while
$W_\parallel$
and
$W_{\textrm{beta-gradB}}$
are negative. In the EFM case (figure 2
h), the relative contribution of
$W_{\textrm{curv1}}$
is noticeably smaller than in the
$m_i/m_e=25$
case, while the other terms are less affected by the
$m_i/m_e$
change. Finally, in the
$B_g=0.2$
case (figures 2
i and 2
j), the increased mass ratio does not affect the relative importance of terms in either of the models (cf. figures 2
c and 2
d).
In summary, we conclude that the domain-integrated LGCM is generally dominated by
$W_{\textrm{curv0}}$
and the EFM is generally dominated by
$W_{\textrm{diam}}$
when
$B_g\leqslant 0.2$
. Increasing the guide field to
$B_g=1$
leads to
$W_\parallel$
becoming larger than (or comparable to)
$W_{\textrm{curv0}}$
and
$W_{\textrm{diam}}$
. The only significant
$m_i/m_e$
dependence is found for the Eulerian
$W_{\textrm{curv1}}$
, which becomes less important in the
$B_g=0$
case as
$m_i/m_e$
is increased from 25 to 100. While the EFM accurately reproduces the total
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
over the simulation domain for all runs, we find that the LGCM estimate of electron energisation agrees better with
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
for increasing values of
$B_g$
. Since the LGCM assumes
$\mu$
-conservation, and a larger
$B_g$
leads to stronger magnetisation, this result hints that the deviation between the LGCM and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
is related to the presence of non-adiabatic electrons.
Similarly to TenBarge et al. (Reference TenBarge, Juno and Howes2024), we find that the local and volume-integrated pictures of electron energisation are very different. Next, we will investigate the spatial distribution of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
and the energisation terms for the
$m_i/m_e=25$
case in more detail.
3.2. Local differences between the LGCM and EFM descriptions
The spatial structures of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
and the
$W$
-terms of (1.2) and (1.3) are shown in figure 3. The left column shows the results of the
$B_g=0$
simulation at
$t\omega _{ci\infty }=16$
. As was seen in figure 1, we find the strongest
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
at the edges of the magnetic islands (figures 3
a and 3
b) and we will focus our attention to these regions.
In the LGCM,
$W_{\textrm{curv0}}$
(figure 3
c, red) is the dominant energisation term due to the strong magnetic curvature associated with the island. In the Eulerian view, the strong diamagnetic drift associated with the pressure gradient at the island edges results in significant
$W_{\textrm{diam}}$
(figure 3
d, magenta), which is the main contributor to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. Approaching the island from the left, the sum of the energisation terms in (1.2) and (1.3) (black dashed lines in figures 3
c and 3
d) are both in good agreement with
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
(purple) until around
$x/d_{i0}\approx 9.8$
, marked by the dash-dotted vertical line. After this point,
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
starts to decrease. This is well captured by the EFM, which reveals a decreasing
$W_{\textrm{diam}}$
while
$W_{\textrm{curv1}}\approx -W_{\textrm{agyro}}$
. In contrast, the LGCM deviates from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
by predicting further energsation deeper into the island as a result of a large
$W_{\textrm{curv0}}$
. Analogous results (with opposite signs) are obtained for the right edge of the island where
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
has a negative peak.
Figure 3. Spatial dependence of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
and the
$W$
terms inside magnetic islands for
$B_g=0$
(left column) and
$B_g=0.2$
(right column). (a)
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
for the whole simulation domain. (b)
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
for the magnetic island boxed in panel (a). The black lines are contours of
$A_z$
. The dashed green lines show the y-range in which the data in panels (c) and (d) are averaged. (c) The different terms in the guiding centre model (1.2) as a function of
$x$
. (d) Same as panel (c) for the fluid model (1.3). The vertical dash-dotted line indicates where the LGCM and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
start to qualitatively deviate. The data in panels (c) and (d) have been averaged in
$y$
over the interval marked by the dashed green lines in panel (b), and smoothed in
$x$
by a moving mean of window size
$\pm 0.5d_{e0}$
. (f–h) Same format as panel (a–d) for the
$B_g=0.2$
case.
In the
$B_g=0.2$
case (right column of figure 3), the LGCM works better than in the
$B_g=0$
case, but it still predicts energisation slightly deeper into the island than the measured
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. In both cases,
$W_{\textrm{curv0}}$
is the dominant term. The EFM remains accurate, with the main difference being that
$W_{\textrm{curv1}}$
is less important than in the
$B_g=0$
case. The local differences between the models brings us to the second question posed in the introduction: why does the LGCM fail inside the islands and what actually happens to the electrons at these locations?
Since the LGCM assumes adiabatic electrons, we start by scrutinising this assumption. One important parameter which can be used to quantify whether the motion of a particle is adiabatic is the so-called adiabaticity parameter
$\kappa =\sqrt {r_{\textrm{curv}}/r_{ge}}$
, where
$r_{\textrm{curv}}=1/|\boldsymbol{k}|$
is the magnetic curvature radius and
$r_{ge}=v_{\perp } m_e/(eB)$
the electron gyroradius (Büchner & Zelenyi Reference Büchner and Zelenyi1989). An electron with
$\kappa \lt 1$
experiences strong magnetic field changes over its gyro-orbit and its magnetic moment is not conserved. We note that one may define a generalised adiabaticity parameter, accounting for all spatial variations of
$\mathbf{B}$
. We do that in Appendix A for reference, but since the alternative definition does not affect the following discussion qualitatively, we will continue using the more widely used definition, given above.
In figures 4(a), 4(b) and 4(c), we present
$|\boldsymbol{J}_e\cdot \boldsymbol{E}|$
, the sum of the right-hand side of (1.2),
$|W_{\textrm{LGCM}}|$
and
$|W_{\textrm{LGCM}}/\boldsymbol{J}_e\cdot \boldsymbol{E}|$
for the
$B_g=0$
case. The largest relative discrepancies between the LGCM and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
are found at the x-points and in the centre of the magnetic islands (figure 4
c), where
$|W_{\textrm{LGCM}}|$
is up to a few orders of magnitude larger than
$|\boldsymbol{J}_e\cdot \boldsymbol{E}|$
. These locations are regions where the thermal
$\kappa$
(using
$v_{te \perp }=\sqrt {2T_{e\perp }/m_e}$
in
$r_{ge}$
) is very small (
$\kappa \ll 1$
), as shown in figure 4(d). Since
$\kappa \ll 1$
implies non-adiabatic electron motion, we conclude that the deviations in these regions can likely be attributed to the erroneous assumption of
$\mu$
conservation made by the LGCM. The fact that electrons can exhibit non-adiabatic behaviour during magnetic reconnection is well known (Zenitani & Nagai Reference Zenitani and Nagai2016), and the failure of the LGCM near the current sheet centre and at the x-points in the low-
$B_g$
limit is not surprising. In other regions, the error does not correlate well with
$\kappa$
. One such example is the separatrices of the larger reconnection sites. There, the LGCM deviates from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
even though
$\kappa \gg 1$
.
Figure 4. Deviations from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
in the LGCM for (a–d)
$B_g=0$
and (e–h)
$B_g=0.2$
. Panels (a) and (b) show respectively
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
and the sum of the guiding centre model terms in (1.2),
$W_{\textrm{LGCM}}$
. Only values larger than
$10^{-5}$
have been included. (c)
$|W_{\textrm{LGCM}}/\boldsymbol{J}_e\cdot \boldsymbol{E}|$
using the data in panels (a) and (b). The thick black contour shows the
$\kappa =1$
threshold, where the bounded regions contain
$\kappa \lt 1$
. The colourbar is saturated. (d) The adiabaticity parameter
$\kappa$
computed using the perpendicular thermal speed. Blue regions correspond do
$\kappa \lt 1$
, i.e. where the motion of thermal electrons is not adiabatic, and red regions to
$\kappa \gt 1$
. The colourbar has been limited to the range
$\log _{10}(\kappa )\in [-0.5,0.5]$
to highlight the adiabatic/non-adiabatic transition. (e–h) Same format as panel (a–d) for the
$B_g=0.2$
case.
In the
$B_g=0.2$
case, the presence of the out-of-plane component leads to a reduced magnetic curvature (i.e. a larger
$r_{\textrm{curv}}$
), and we generally find
$\kappa \gg 1$
(figure 4h). Small values of
$\kappa$
are still found near the x-points and
$\kappa$
occasionally drops below 1. It should be noted that the plotted
$\kappa$
is for thermal electrons and electrons with
$v_\perp \gt v_{te\perp }$
have smaller
$\kappa$
. A large fraction of the electrons might therefore be non-adiabatic when the thermal
$\kappa \approx 1$
. The fact that there are still significant deviations even though
$\kappa \gt 1$
(figures 3
g and 4
h) suggests that there is another issue at play other than the lack of
$\mu$
-conservation. This suspicion is reinforced by the fact that the LGCM locally deviates from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
even in the
$B_g=1$
case (see Appendix B), where
$\kappa \gg 1$
everywhere.
We conclude that the non-adiabatic motion of electrons near x-points and the current sheet centre, particularly in the
$B_g=0$
case, is a large source of error in the LGCM. However, this is not the entire picture, and there is still the issue of the transition between the Lagrangian and Eulerian frameworks when (1.2) is evaluated on an Eulerian grid. This problem was discussed by TenBarge et al. (Reference TenBarge, Juno and Howes2024), and in the following section, we will expand on their discussion to pin-point the issue and illustrate how this affects the local electron energisation estimate.
3.3. Eulerian versus Lagrangian energisation
Specifically, the problem with evaluating (1.2) on an Eulerian grid (assuming
$\mu$
is conserved) is that it gives the impression that it provides an Eulerian view on electron energisation, i.e. that it describes the evolution of the total electron energy within a fixed volume element. However, (1.2) is derived in a Lagrangian framework from the single guiding centre equation (1.1), which describes the rate of energy gain a single guiding centre experiences at some given time. To derive (1.2) from (1.1), one takes a small volume (such that all charges experience the same fields) and adds up (1.1) for all guiding centres within the volume. Thus, what the quantity
${\rm d}\mathcal{E}_{\textrm{LGCM}}/{\rm d}t$
in (1.2) describes is the rate of energy gain per unit volume of the specific guiding centres that are currently in the volume element. Importantly, this means that the electron energy flux into or out of the volume is irrelevant for
${\rm d}\mathcal{E}_{\textrm{LGCM}}/{\rm d}t$
. This detail distinguishes the Lagrangian description from the Eulerian description, which is concerned with the evolution of the total electron energy content within the volume. The Eulerian electron energy equation therefore contains divergence terms related to the electron energy flux:
\begin{equation} \frac {\partial \mathcal{E}_{\textrm{EFM}}}{\partial t} + \nabla \cdot \left (\boldsymbol{K}+\boldsymbol{H}+\boldsymbol{q}\right ) = \boldsymbol{J}_e\cdot \boldsymbol{E}, \end{equation}
where
$\mathcal{E}_{\text{EFM}}=m_en_eu_e^2/2+\text{trace}(\mathbf{P}_e)/2:=K+U$
is the total electron energy density with
$K$
and
$U$
being respectively the bulk kinetic and thermal contributions,
$\boldsymbol{K}=K\boldsymbol{u}_e$
is the kinetic energy flux density,
$\boldsymbol{H}=U\boldsymbol{u}_e+\mathbf{P}_e\cdot \boldsymbol{u}_e$
is the enthalpy flux density, and
$\boldsymbol{q}$
is the heat flux density (e.g. Helander & Sigmar Reference Helander and Sigmar2005). The energy equation describes the evolution of the electron energy density inside a fixed volume element, a quantity that is fundamentally different from the aforementioned Lagrangian
${\rm d}\mathcal{E}_{\textrm{LGCM}}/{\rm d}t$
.
As evident from (3.1),
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
is a source term and it can be balanced by the time derivative and divergence (i.e. energy flux density) terms. The importance of the divergence terms on smaller scales is readily seen in figure 5, where the spatial distribution of the different components of (3.1) are shown in panels (a)–(h). The large
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
found at the edges of magnetic islands is mainly balanced by the divergence terms in the energy equation (see figures 5
a and 5
h, where
$\boldsymbol{Q}=\boldsymbol{K}+\boldsymbol{H}+\boldsymbol{q}$
), where
$\nabla \cdot \boldsymbol{H}$
dominates (figure 5
c) over
$\nabla \cdot \boldsymbol{K}$
(figure 5
b) and
$\nabla \cdot \boldsymbol{q}$
(figure 5
d). The time derivative terms (figure 5
e–g) tend to have the opposite sign to the divergence terms, and the internal energy term (figure 5
f) dominates over the kinetic energy term (figure 5
g). The large contribution of the energy flux terms, particularly
$\nabla \cdot \boldsymbol{H}$
, is consistent with in situ studies of magnetic reconnection (Eastwood et al. Reference Eastwood2020; Fargette et al. Reference Fargette, Eastwood, Waters, Øieroset, Phan, Newman, Stawarz, Goldman and Lapenta2024). We conclude that, on a local level, most of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
is balanced by the electron energy flux densities, with the dominant term being the enthalpy flux density. This is a key reason why the LGCM locally deviates from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
, even when the assumption of magnetic moment conservation is valid.
Figure 5. The contribution of the different terms in (3.1) to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
in the
$B_g=0$
,
$m_i/m_e=25$
simulation. (a–h) Spatial profiles for the different terms at
$t\omega _{ci\infty }=16$
. The top row of panels show the three terms in (3.1), where
$\boldsymbol{Q}=\boldsymbol{K}+\boldsymbol{H}+\boldsymbol{q}$
. The first and second columns of panels show the different divergence and time derivative terms, respectively. The divergence terms in panel (a–d) have been smoothed using a 5-point moving average to reduce noise. (i) Time dependence of the different terms summed over the whole simulation domain.
On system scales, however, the picture is very different. For our choice of boundary conditions, the system is closed in the sense that there is no net particle energy flux into or out of the simulation domain. Therefore, when we integrate the divergence terms of (3.1) over the whole simulation, the corresponding surface terms vanish, leaving us with
$\boldsymbol{J}_e\cdot \boldsymbol{E}\approx \partial \mathcal{E}_{\textrm{EFM}}/\partial t$
, as seen in figure 5(i). This is why the domain-integrated LGCM can accurately reproduce
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
when
$\mu$
-conservation applies (figure 2
e). In that case, (1.2) describes the evolution of the total energy of all guiding centres, and since they are confined to the simulation domain, this is the same thing as the evolution of the total electron energy in the simulation domain. So, when we look at a closed system in its entirety, indeed
$\int d^3x(d\mathcal{E}_{\textrm{LGCM}}/dt) = \int d^3x(\partial \mathcal{E}_{\text{EFM}}/\partial t)=\int d^3x\boldsymbol{J}_e\cdot \boldsymbol{E}$
. The difference between the domain-integrated LGCM and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
observed in the
$B_g=0$
case (figure 2
a) can thus be attributed to the lack of
$\mu$
conservation due to the small
$\kappa$
(figure 4
d). The fact that the LGCM equation is locally invalid in the
$B_g=0$
case explains why the domain-integrated LGCM both overestimates and underestimates
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
in the
$B_g=0$
case (figure 2
a), and subsequently why the time-integrated energisation more accurately reproduces
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
.
Until this point, we have followed TenBarge et al. (Reference TenBarge, Juno and Howes2024) and uncritically assumed that, since
${\rm d}\mathcal{E}_{\textrm{LGCM}}/{\rm d}t$
accurately describes the energisation of individual guiding centres, it must equal the work done on the electrons as quantified by
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. In other words, since we can write (1.2) as the scalar product of a net current density
$\boldsymbol{J}_{e,{\textrm{LGCM}}}$
and
$\boldsymbol{E}$
,
$\boldsymbol{J}_{e,\textrm{LGCM}}$
should be equal to the bulk fluid current density
$\boldsymbol{J}_e$
. This, however, is evidently not the case. As already mentioned by TenBarge et al. (Reference TenBarge, Juno and Howes2024), the guiding centre model is a particle model that is gyroperiod integrated, and it consequently does not describe certain bulk fluid properties such as the diamagnetic drift and the corresponding currents. Therefore, in general,
$\boldsymbol{J}_{e,\textrm{LGCM}}\neq \boldsymbol{J}_e$
, since
$\boldsymbol{J}_e$
describes the total electron current. In hindsight, it is therefore quite clear that we should never have expected
${\rm d}\mathcal{E}_{\textrm{LGCM}}/{\rm d}t$
to accurately model
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
to begin with, even when
$\mu$
is conserved.
Finally, we note that while (1.2) is derived in the Lagrangian framework, the right-hand terms are expressed in terms of bulk fluid quantities which, in our case, are known on the points defined by the Eulerian grid. It is perfectly valid to evaluate (1.2) on an Eulerian grid (if
$\mu$
is conserved); one must just keep in mind that a local measurement of
$d\mathcal{E}_{\textrm{LGCM}}/dt$
describes the energy evolution of the guiding centres that are presently in the local volume, not the evolution of the electron energy within the volume. In their paper, TenBarge et al. (Reference TenBarge, Juno and Howes2024) state, regarding the LGCM results, that ‘[…] neither globally nor locally does the sum of the energisation mechanisms agree well with
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
, suggesting that the Lagrangian description, unsurprisingly, does not work well in an Eulerian simulation’, which may give the impression that there is an inherent conflict between their Eulerian solver and the LGCM. We want to stress that the inherently Eulerian nature of a continuum Vlasov–Maxwell solver does not affect the applicability of (1.2). The deviations from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
observed by TenBarge et al. (Reference TenBarge, Juno and Howes2024) on a global level can likely be attributed to the demagnetisation of electrons near the x-point, as their choice of
$B_g=0.1$
is not large enough to keep the electrons well magnetised (see our
$B_g=0.2$
and
$B_g=1$
results in figures 2
c and 2
e). The local deviations are likely, as discussed above, a result of the combination of non-conserved
$\mu$
and energy fluxes.
4. Discussion
4.1. Interpretation of EFM energisation terms
With a better understanding of the differences between the LGCM and EFM, we are in a position to investigate the last question posed in the introduction, namely: can we associate local measurements of each EFM
$W$
-term with the generation of specific VDF features in a manner analogous to the betatron and Fermi terms of the LGCM? For example, one might be tempted to infer that
$W_{\textrm{diam}}\gt 0$
leads to the formation of an anisotropic thermal feature in the electron VDF. Below, we provide some arguments as to why this is not necessarily the case.
Take the dominant
$W_{\textrm{diam}}=\mathbf{u}_E\cdot \nabla p_{e\perp }$
term as an example. This term describes the scalar product of the diamagnetic drift
$\boldsymbol{u}_{\textrm{diam}}=-\nabla p_{e\perp } \times \boldsymbol{B}/(q_en_eB^2)$
with
$q_en_e\boldsymbol{E}$
. Unlike the curvature and magnetic gradient drifts that are present in the LGCM, the diamagnetic drift does not need to correspond to the drift of individual electrons. Rather, the pressure gradient introduces a local anisotropy in the momenta carried by gyrating electrons. If we apply a homogeneous electric field directed along
$\boldsymbol{J}_e$
over such a pressure gradient, the individual electrons will experience no net energy gain over a gyroperiod even though locally
$\boldsymbol{J}_e\cdot \boldsymbol{E}\gt 0$
. What does happen, however, is that all electrons (and therefore also the pressure gradient) start
$\boldsymbol{E}\times\boldsymbol{B}$
-drifting in the direction of the higher pressure. Thus, on a fluid level, the work done by the electric field on the diamagnetic bulk flow goes into pushing the fluid along the pressure gradient at the
$\boldsymbol{E}\times \boldsymbol{B}$
-velocity, as evident by the relation
$q_en_e\boldsymbol{u}_{\textrm{diam}}\cdot \boldsymbol{E}=\boldsymbol{u}_E\cdot \nabla p_{e\perp }$
.
To illustrate how the above scenario affects the local energy evolution of the plasma, we consider an Eulerian volume element that is initially in the high-pressure region. At the start, the uniform plasma is
$\boldsymbol{E}\times \boldsymbol{B}$
-drifting across the volume element, resulting in
$\nabla \cdot \boldsymbol{Q}=\partial \mathcal{E}_{\textrm{EFM}}/\partial t=\boldsymbol{J}_e\cdot \boldsymbol{E}=0$
. When the pressure gradient enters the Eulerian element, however, the presence of the diamagnetic drift yields a positive
$\boldsymbol{J}_e\cdot \boldsymbol{E}=W_{\textrm{diam}}$
. At the same time, there is a net outward flux of high-pressure plasma that is being replaced by inflowing low-pressure plasma. In terms of the energy equation, this corresponds to
$\boldsymbol{J}_e\cdot \boldsymbol{E}\gt 0$
,
$\nabla \cdot \boldsymbol{Q}\gt 0$
and
$\partial \mathcal{E}_{\textrm{EFM}}/\partial t\lt 0$
. This is exactly the situation we find at the island edge around
$x/d_{i0}=10$
in figure 5. We therefore end up with the perhaps counterintuitive conclusion that in the EFM, a positive
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
can lead to a decreasing energy density, due to the work being diverted into energy fluxes. We also note that the
$\boldsymbol{J}_e\cdot \boldsymbol{E}\gt 0$
region is not the source of the bulk flow, since the entire plasma is drifting, and it would do so even in the absence of the pressure gradient.
So far, we have considered an
$\boldsymbol{E}$
-field that is homogeneous over the electron gyro-orbits, such that no single electron gains energy during a gyroperiod in the field. If, in contrast, there are
$\boldsymbol{E}$
-field gradients, then the individual electrons can gain a net energy by the
$\boldsymbol{E}$
-field over their gyroperiod, leading to the formation of some interesting VDF feature. Importantly, however,
$W_{\textrm{diam}}$
describes the local energisation due to the local electric field and the local fluid current; it does not distinguish between situations where individual electrons gain energy and situations where they do not. It is therefore not possible to conclude from a local measurement of
$W_{\textrm{diam}}$
whether or not individual electrons are energised, or how the VDF will change.
For completeness, a few words on what differentiates the LGCM from the EFM in this context. The reason the LGCM is able to provide insights into the evolution of the VDF (assuming
$\mu$
conservation) is twofold. First, the LGCM follows the evolution of collections of guiding centres, not the evolution of the total energy content within fixed grid cells as the EFM. This means that it exclusively provides information about particle energisation, and no information about energy fluxes. Second, the LGCM is already integrated over the gyroperiod, meaning that the effects of electric field gradients (or lack thereof) on the electron gyro-scales are already accounted for. There are, of course, also limitations to the electric field scales that allow for
$\mu$
-conservation (Stephens, Brzozowski & Jenko Reference Stephens, Brzozowski and Jenko2017), but here we assume that the LGCM is applied validly.
4.2. Extrapolating to the physical mass ratio
The results presented in §§ 3.2 and 3.3 were obtained in the
$m_i/m_e=25$
case. Similar results were also obtained for
$m_i/m_e=100$
. Specifically, values of
$\kappa \lt 1$
and large errors in the guiding centre model at the x-points and near the current sheet centre were also found in the
$m_i/m_e=100$
,
$B_g=0$
case. How these results change for more realistic mass ratios remains to be investigated. It is possible that an increased mass ratio could affect the relative importance of the
$W$
-terms. However, as mentioned in the above paragraphs, the EFM energisation terms do not necessarily reflect the evolution of the VDF, and any application of a decomposition of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
into the EFM
$W$
-terms must be designed accordingly. Moreover, accurately computing the EFM terms using in situ data is a challenging task since, in addition to the
$p_\perp$
and
$p_\parallel$
required by the LGCM, the EFM also requires accurate measurements of pressure gradients in
$W_{\textrm{diam}}\propto \nabla p_{e\perp }$
and
$W_{\textrm{agyro}}\propto \nabla \cdot \boldsymbol{\Pi }^a_e$
.
One important result from our analysis is that the LGCM can produce significant errors due to the erroneous assumption of
$\mu$
-conservation when
$\kappa \lt 1$
. In our simulations, we found regions of low
$\kappa$
at the x-points, near the island edges, and inside the islands for
$B_g=0$
(figure 4
d). Increasing the guide field to
$0.2$
increased
$\kappa$
, particularly inside the islands, but low-
$\kappa$
regions still remained. Since the reconnecting current sheet is on electron scales, increasing the mass ratio should not affect
$\kappa$
near the current sheet centre or in the electron diffusion region. Indeed, from spacecraft observations, it is well known that electrons demagnetise in these regions (e.g. Burch et al. Reference Burch2016; Torbert et al. Reference Torbert2018; Cozzani et al. Reference Cozzani2019). However, since the magnetic islands grow into ion-scale structures, the associated curvature radius should increase for increasing mass ratios. We might therefore expect that the LGCM problem of demagnetised electrons at magnetic islands should be less of a concern for physical mass ratios. However, the islands still grow from electron scales, so this argument would only be valid in the later stages of reconnection. How well the electron magnetic moment is conserved may also depend on several other parameters like the plasma beta and the current sheet type (Harris versus force-free). These parameters can therefore affect the validity of the LGCM. In any case, we know from spacecraft data that electron scale structures with strong gradients can be present in physical reconnection outflows in geospace (e.g. Zhou et al. Reference Zhou2019; Leonenko et al. Reference Leonenko, Grigorenko, Zelenyi, Malova, Malykhin, Popov and Büchner2021). It is therefore important to be aware of the fact that the LGCM assumes
$\mu$
conservation, and to test that assumption before applying the model.
5. Conclusions
In the present paper, we build on the recent work by TenBarge et al. (Reference TenBarge, Juno and Howes2024) to study the Eulerian fluid model (EFM; (1.3)) and the Lagrangian guiding centre model (LGCM; (1.2)) of electron energisation during multiple x-point reconnection using particle-in-cell simulations. We set out to answer the series of questions summarised below.
-
(i) Which EFM energisation terms are most important during multi-x-point reconnection and are they affected by the guide field (
$B_g$
) or the ion-to-electron mass ratio? -
(ii) What is the reason for the different predictions by the LGCM and EFM, and what does that tell us about the electrons?
-
(iii) The energisation terms in the LGCM are often interpreted as leading to specific changes in the electron velocity distribution function (VDF). Can we make similar interpretations for the EFM energisation terms?
(i) We find that the energisation related to the diamagnetic drift (
$W_{\textrm{diam}}$
) is generally, on large scales, the dominant contributor to
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
in the EFM (see figure 2). For large guide fields (
$B_g=1$
), the contribution of parallel electric fields,
$W_\parallel$
, becomes the most important term (albeit still comparable to
$W_{\textrm{diam}}$
). The increased importance of
$W_\parallel$
for larger
$B_g$
has previously been observed for the LGCM (Dahlin et al. Reference Dahlin, Drake and Swisdak2016). The only clear mass-ratio dependence is found for the Eulerian curvature term (
$W_{\textrm{curv1}}$
) in the
$B_g=0$
case, where it becomes less important as
$m_i/m_e$
is increased from 25 to 100. Locally, other energisation terms can contribute significantly to electron energisation (see figure 3).
(ii) We show that the assumption of magnetic moment conservation made by the LGCM is invalid in the low
$B_g\leqslant 0.2$
case near the x-points and at the centre of the current sheets where the magnetic field is strongly curved. This causes significant errors in the LGCM energisation estimate. For
$B_g=1$
, the electrons remain magnetised, and we find a good agreement between the simulation domain integrated LGCM and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. Local deviations from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
are produced by the LGCM in all simulations, even when
$\mu$
is conserved. This is due to the fact that the LGCM and the EFM describe two fundamentally different quantities. The LGCM describes the evolution of the total energy of the guiding centres that are currently within a given volume element. In contrast, the EFM describes the evolution of the total electron energy content within the volume element. The latter includes the effects of energy fluxes through the volume element, whereas the former does not. This has the consequence that the two models can both be correct, while simultaneously giving seemingly contradictory local energisation rates if there is a net energy flux into or out of the volume element under consideration. By examining the different terms of the electron energy equation in figure 5, we find that the energy flux terms are important in balancing
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
on a local level, explaining the difference between the two models when
$\mu$
is conserved.
(iii) Finally, we argue that, unlike the LGCM, we cannot associate local measurements of the different EFM energisation terms with the formation of structures in the electron VDF. In essence, this is due to the fact that the EFM includes the local effect of energy fluxes and fluid drifts which need not correspond to the drift of individual particles.
In summary, the EFM describes a complete decomposition of
$\boldsymbol{J}_e$
into fluid drifts, and it therefore accurately reproduces
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. Local measurements of each constituent energisation term, however, need not describe locations where the plasma is either heated or accelerated. In contrast, the LGCM accurately describes the energisation of particles if it is used validly, i.e. when
$\mu$
is conserved, but it does not describe the local evolution of the bulk fluid energy. Local deviations between the energisation measures of the two models, even when
$\mu$
is conserved, can be attributed to the contrasting Eulerian and Lagrangian perspectives, and does not mean that one of the models is wrong. A thorough understanding of what one wishes to describe is therefore necessary before one decides to use either of the energisation models.
Acknowledgments
We are grateful to J. M. TenBarge for helpful correspondence at an early stage of this work. The computations and data handling were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.
Editor Nuno Loureiro thanks the referees for their advice in evaluating this article.
Funding
This project received funding from the Knut and Alice Wallenberg Foundation (Grant No. KAW 2022.0087) and the Swedish Research Council (Dnr. 2021-03943). K.S. acknowledges financial support from the Adlerbertska Research Foundation.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Generalised adiabaticity estimate
The standard adiabaticity parameter
$\kappa =\sqrt {r_{\textrm{curv}}/r_{ge}}$
uses the magnetic field curvature radius (Büchner & Zelenyi Reference Büchner and Zelenyi1989),
$r_{\textrm{curv}}=1/|\boldsymbol{b}\cdot \nabla \boldsymbol{b}|$
, to quantify the adiabaticity, and does not account for the scales of magnetic field gradients. A generalised length scale that incorporates all magnetic field variations can be constructed as
$L_{\nabla \boldsymbol{B}}^{-1}=\Vert \nabla \boldsymbol{B}/B\Vert$
, where
$\Vert \cdot \Vert$
denotes the spectral norm. This length scale captures both curvature and gradient effects. We can then define an adiabaticity measure analogous to
$\kappa$
as
$\kappa ^*=\sqrt {L_{\nabla \boldsymbol{B}}/r_{ge}}$
. Figure 6 shows that the regions where
$\kappa$
predicts non-adiabatic electrons (
$\kappa \lt 1$
) are confined within regions of
$\kappa ^*\lt 1$
, as expected. Because the magnetic field is strongly curved, the two quantities are similar, and the
$\kappa ^*\lt 1$
regions only extend slightly outside the
$\kappa \lt 1$
regions in the
$B_g=0$
case (figure 6
b, c). In the
$B_g=0.2$
case, the generalised
$\kappa ^*$
picks up potentially non-adiabatic regions at island boundaries (
$\kappa ^*\approx 1)$
, which are not captured by
$\kappa$
. These are regions with strong magnetic field gradients.
Figure 6. Difference between
$\kappa$
and
$\kappa ^*$
for
$B_g/B_\infty =0$
(left column) and
$B_g/B_\infty =0.2$
(right column). (a) Ratio between the LGCM sum and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. The thick black and purple contours correspond to
$\kappa =1$
and
$\kappa ^*=1$
, respectively. (b)
$\kappa$
calculated using
$r_{\textrm{curv}}$
. The thick black contours correspond to
$\kappa =1$
. (c)
$\kappa ^*$
calculated using
$L_{\nabla \boldsymbol{B}}$
. The thick purple contours correspond to
$\kappa ^*=1$
. (d–f) Same format as panel (a–c). All colourmaps are saturated.
In summary, other length and time scales can also affect the adiabaticity of electrons (Stephens et al. Reference Stephens, Brzozowski and Jenko2017). However,
$\kappa$
is sufficient to show that electrons can be non-adiabatic in both the low and intermediate guide field cases, and that
$\kappa \lt 1$
correlates strongly with large deviations between the LGCM sum and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
.
Appendix B. Local energisation and demagnetisation for
$B_g=1$
In this appendix, we present the data corresponding to figures 3 and 4 for the
$B_g=1$
case in figure 7.
Figure 7. (a–d) Spatial dependence of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
and the
$W$
terms inside magnetic islands for
$B_g=1$
. Same format as figure 3. (e–h) Deviation between the LGCM sum and
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
for
$B_g=1$
. Same format as figure 4. Note that we have increased the range of the colourmap in panel (h) compared to figure 4, as
$\log _{10}(\kappa )\gt 0.5$
everywhere in this case.
Figure 7(c) shows that the sum of the LGCM energisation terms deviates significantly from
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
in the
$B_g=1$
case. The EFM (figure 7
d), however, provides an accurate estimate of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
. We note that the reversed polarity of
$\boldsymbol{J}_e\cdot \boldsymbol{E}$
compared with the
$B_g=0$
and
$B_g=0.2$
cases (see figure 3) is due to the fact that the island in the
$B_g=1$
case is convecting in the opposite (
$-x$
) direction.
As shown in figure 7(g), the LGCM tends to provide an energisation estimate that is much less than
$|\boldsymbol{J}_e\cdot \boldsymbol{E}|$
, particularly inside the magnetic islands. These deviations cannot be explained by electron demagnetisation as quantified by
$\kappa$
, since
$\kappa \gg 1$
throughout the whole simulation domain (figure 7
h).
Open access
