Eikonal Equation 0.1

Iryna YEHORCHENKO

Institute of Mathematics of NAS Ukraine, 3 Tereshchenkivs’ka Str., Kyiv-4,
Ukraine

E-mail<EMAIL_ADDRESS>
Institute of Mathematics, PAN, Poland

<EMAIL_ADDRESS>
###### Abstract

We provide a review of some symmetry-related literature on the eikonal
equations $u_{\mu}u_{\mu}=0$,$u_{\mu}u_{\mu}=1$, where lower indices at
dependent variables designate derivatives, $\mu=0,1,2,..,n$ and summation is
implied over the repeated indices.

We will consider general solutions and symmetries of these equations, and
relations of these equations with reduction of higher-order PDE. Some new
results that were needed for comprehensive presentation are also adduced. In
particular, we will also consider discrete symmetries of the eikonal equations
equivalence classes of solutions and relations of the symmetry solutions and
the general solutions. We describe in detail the procedure that allowed
obtaining of the general solution using hodograph and contact transformations
of the initial equations. Some new results by the author that were needed for
comprehensive presentation are also adduced - in particular, a general
solution of the eikonal equation for an arbitrary number of space variables.

The paper is intended both as a reference material and as a manual.

##  1 Notations and Arrangements

We will consider here a class of nonlinear first-order partial differential
equations (PDE) in the following form:

$u_{\mu}u_{\mu}=f(x_{\mu},u),$ (1)

where in general we have one real-valued dependent variable
$u=u(x_{0},x_{1},\ldots,x_{n})$, a time variable $x_{0}$ and $n$ space
variables $x_{1},\ldots,x_{n}$. I would like to point out that $u$ and other
dependent variable take values in a real space, as the complex-valued
dependent variables require different treatment and a special study.

If not indicated otherwise, the lower indices of dependent variables designate
derivatives with respect to the relevant independent variables $x_{\mu}$,
$\mu=0,1,...,n$. Here we will present results for an arbitrary number of the
space variables, but will consider in detail examples for specific cases of
$n=1,2,3$.

If not indicated otherwise, the lowercase Greek letters used as indices are
assumed to run from $0$ to $n$, where $n$ is an arbitrary natural number,
usually $n>1$ (the case $n=1$ needs special consideration and will be
discussed separately. Lowercase Greek letters used as indices (e.g. $\mu$ will
designate the $\mu$-th independent variable $x_{\mu}$ or used as an index at a
symbol for a dependent variable the partial derivative with respect to
$x_{\mu}$).

Also if not indicated otherwise, the lowercase Latin letters used as indices
are assumed to run from $1$ to $n$, where $n$ is an arbitrary natural number.
Usually $n>1$ (the case $n=1$ needs special consideration and will be
discussed separately). The lowercase $n$ will designate the number of the
independent space variables, or when used as an index will designate the
$n$-th independent variable $x_{n}$ or the partial derivative with respect to
$x_{n}$. Other lowercase Latin letters used as indices (e.g. $k$ will
designate the $k$-th independent variable $x_{k}$ or used as an index at a
symbol for a dependent variable the partial derivative with respect to
$x_{k}$).

The following summation arrangements will be assumed (if not specifically
indicated otherwise): repeated lowercase Latin indices will mean summation
from $1$ to $n$ in the Euclidean space:

$x_{k}x_{k}=x_{1}x_{1}+\ldots+x_{n}x_{n}.$

A summation arrangement using an expression on indices will be used as follows
($m$ takes values from 1 to $n-k$):

$x_{k+m}w_{m}=x_{k+1}w_{1}+\ldots+x_{n}w_{n-k}.$

The repeated lowercase Greek indices will mean summation from $0$ to $n$ in
the Minkowski space:

$x_{\mu}x_{\mu}=x_{0}x_{0}-x_{1}x_{1}-\ldots-x_{n}x_{n}.$

The indices at dependent variables that are used to indicate partial
derivatives will mean the following:

$\displaystyle u_{k}=\frac{\partial u}{\partial x_{k}},\qquad
u_{kl}=\frac{\partial u}{\partial x_{k}\partial x_{l}},$ (2)

with a similar formulae applicable to the lowercase Greek indices at dependent
variables, and to higher-order partial derivatives.

A letter with an index below

$\mathop{u}\limits_{m}$

means the set of all partial derivatives of the function $u$ of the order $m$,
with respect to all its independent variables.

$g_{\mu\nu}$ is a metric tensor

$g_{\mu\nu}={\rm diag}\\{1,-1,\ldots,-1\\}$

in an $(n+1)$-dimensional Minkowski space.

$\delta_{ab}$ is a metric tensor

$\delta_{ab}={\rm diag}\\{1,1,\ldots,1\\}$

in an $n$-dimensional Euclidean space.

We will use the following notations in presentation of first-order
differential operators:

$\displaystyle\partial_{\mu}=\frac{\partial}{\partial
x_{\mu}},\qquad\partial_{u}=\frac{\partial}{\partial
u},\qquad\partial_{u_{\mu}}=\frac{\partial}{\partial{u_{\mu}}}.$

We assume that all functions considered are sufficiently smooth, and all
dependent and independent variables take values in the real space.

There is a serious problem with equation names in the mathematics in general.
The name ”eikonal equation” might be used for something different than (1),
and the equations from the class (1) or with specific values of $f(x_{\mu},u)$
may be called also ”Hamilton equations” for $f(x_{\mu},u)=1$. E.g. in [1] the
same equation (1) for $f=0$, $n=2$ was called a Hamilton-Jacobi equation.

We will call the Hamilton-Jacobi equation the following equation:

$u_{t}-u_{a}u_{a}=f(t,x_{a},u).$ (3)

##  2 Introduction

The idea of this paper emerged from the need of the author to have a
comprehensive reference material to be cited in the papers presenting new
results, and from the lack of anything suitable in the existing literature. It
looks like many useful formulae from the author’s manuscript notebooks and not
actually accessible literature, e.g. [2], are worth to be published for the
mathematical community in general. It is worth noting that even a Wikipedia
article ”Eikonal Equation” [3] as of 30.11.2022 does not contain any
information on the symmetry and on the known parametric general solution of
the eikonal equation for low dimensions.

The intended audience of this paper is wider than ”the symmetry analysis
people”, so many detailed arrangements and calculations that are obvious for
the above group are written explicitly. The results that can be attributed to
other authors are cited as such. However, I failed to find relevant citations
for many results and calculations, even though they are not difficult and are
probably of a student-exercise level, but needed in a reference paper. Such
results will be presented without citations, and I will greatly appreciate any
information for the relevant references if available. The list of the results
in this paper I believe at the moment of writing to be new will be adduced at
the end. The literature cited cannot be regarded as a comprehensive list, I
mentioned only was immediately relevant for the presentation.

##  3 Symmetry of the eikonal equations with 0 and 1

Full group classification for the equations of the form

$u_{a}u_{a}=F(t,u,u_{t})$ (4)

was performed in [4].

There are two eikonal equations, each having remarkable symmetry properties:

$u_{\mu}u_{\mu}=0,$ (5)

and

$u_{\mu}u_{\mu}=1.$ (6)

The maximal Lie symmetry algebra of continuous transformations for the
equation (5) may be described by the following basis operators:

$\displaystyle
X=2g_{\mu\nu}c^{\mu}(u)x_{\nu}x_{\beta}\partial_{\beta}-g_{\mu\nu}c^{\beta}(u)x_{\mu}x_{\nu}\partial_{\beta}+g_{\mu\nu}b_{\mu\beta}(u)x_{\nu}\partial_{\beta}+$
(7)
$\displaystyle+d(u)x_{\beta}(u)\partial_{\beta}+a^{\beta}(u)\partial_{\beta}+\eta(u)\partial_{u};$

and here $c^{\mu}$, $b_{\mu\beta}(u)$, $d(u)$, and $a^{\beta}$ are arbitrary
sufficiently smooth functions.

The maximal symmetry algebra of continuous transformations for the equation
(6) may be described by the following basis operators:

$\displaystyle\partial_{\mu},\qquad\partial_{u},\qquad
J_{ab}=x_{a}\partial_{b}-x_{b}\partial_{a},$ (8) $\displaystyle
J_{0a}=x_{0}\partial_{b}+x_{b}\partial_{0},\qquad
J_{u0}=u\partial_{0}+x_{0}\partial_{u},\qquad
J_{ua}=u\partial_{a}-x_{a}\partial_{u},$ $\displaystyle
D=x_{0}\partial_{0}+x_{a}\partial_{a}+u\partial_{u},\qquad
K_{a}=2x_{a}D+s^{2}\partial_{a},\qquad K_{u}=2uD+s^{2}\partial_{u},$
$\displaystyle
K_{0}=2x_{0}D-s^{2}\partial_{0}\qquad(s^{2}=x_{\mu}x_{\mu}-u^{2}).$

Here we used the forms of the symmetry operators listed in [4]. In (7), (8) we
assume that $n$ is arbitrary, though to our knowledge these results were
strictky proved only for the case $n=3$.

The maximal Lie symmetry of the equation (8) for $n=3$ was found in [5], and
for (7) where $n=3$ it was found in [6]. Lie symmetry for (7), $n=2$ was
described in [1]. There exist some papers with the titles like ”Symmetries of
the Eikonal equation” (see e.g. [7]) where authors are obviously not
acquainted with the previous literature on the subject and basics of the Lie
symmetry analysis, do not cite many available papers on low-dimensional
eikonal equations and either repeat previously known results or present not
substantiated and not comprehensive results.

So, I see the need to give here a proof of the results on the maximal Lie
symmetry of the equations (5) and (6) stated by the formulae in (7), (8).

A proof would require going to the basics of the Lie symmetry analysis
procedures (see e.g. [8]).

An equation $F(x,u,\underset{1}{u},\ldots,\mathop{u}\limits_{k})=0$ (we
consider here one scalar functions) is invariant under a first-order
differential operator

$\displaystyle
X=\xi^{\mu}(x_{\nu},u)\partial_{x_{\mu}}+\eta(x_{\nu},u)\partial_{u}$ (9)

such that

$\mathop{X}\limits_{k}F\big{|}_{F=0}=0,$ (10)

where $\mathop{X}\limits_{k}$ is a $k$-th prolongation of the operator $X$.

A proof that some set of first-order differential operators is indeed a
maximal Lie invariance algebra of the equation being considered is to prove
that a linear combination of the operators from this set is a maximal operator
satisfying the condition (10).

For our particular proof we need only first prolongation of the operator $X$
(9) (this is the operator acting not only on independent and dependent
variables, but also at the derivatives, and its coefficients are fully
determined by the coefficients of the operator $X$:

$\displaystyle\mathop{X}\limits_{1}=X+\zeta^{\mu}(x_{\mu},u)\partial_{u_{\mu}},$
(11)

where the coefficients

$\displaystyle\zeta^{\mu}(x_{\mu},u)=\eta_{\mu}+\eta_{u}u_{\mu}-\xi^{\nu}_{\mu}u_{\nu}-\xi^{\nu}_{u}u_{\mu}u_{\nu}.$
(12)

Application of the operator (11),(12) to the equations (5) and (6) gives the
following expression:

$\displaystyle\mathop{X}\limits_{1}u_{\mu}u_{\mu}=2\zeta^{\mu}u_{\mu}=2(\eta_{\mu}u_{\mu}+\eta_{u}u_{\mu}u_{\mu}-\xi^{\nu}_{\mu}u_{\nu}u_{\mu}-\xi^{\nu}_{u}u_{\mu}u_{\mu}u_{\nu})=$
(13) $\displaystyle
2(\eta_{\mu}u_{\mu}-\xi^{\nu}_{\mu}u_{\nu}u_{\mu})+2u_{\mu}u_{\mu}(\eta_{u}-\xi^{\nu}_{u}u_{\nu})$

From the condition that (13) has to vanish (even before transition to the
relevant manifolds determined by equations (5) and (6)) we get for both
equations the following conditions for the coefficients of the operator $X$
from coefficients at $u_{\mu}u_{\nu}$:

$\displaystyle\xi^{\nu}_{\mu}+\xi^{\nu}_{\mu}=0$ (14)
$\displaystyle\xi^{\nu}_{\nu}=\xi^{\mu}_{\mu}\qquad{\rm if}\qquad\nu\neq\mu.$

No summation over the repeated indices is implied in (14).

Transition to the manifold of equation (5) gives more conditions in addition
to (14):

$\displaystyle\eta_{\mu}=0.$ (15)

From (14) and (15) we get explicit forms of coefficients of the operator $X$
(7).

Transition to the manifold of equation (6) gives more conditions in addition
to (14):

$\displaystyle\eta_{\mu}=\xi^{\mu}_{u},\qquad\eta_{u}=\xi^{\mu}_{\mu},$ (16)

no summation over the repeated indices is implied in (16).

From the equations (14) and (16) we get coefficients of the operator $X$ and
the basis of the maximal symmetry algebra of continuous transformations (8).
Conditions (14) imply presence of the Poincaré and conformal algebras. For (6)
these algebras appear in the space of one more variable than e.g. for the wave
equation $\Box u=0$ in $n$ space dimensions - the dependent variable is
included in the same way as the independent variables. The maximal Lie
invariance algebra for (5) appears to be infinitely dimensional, and
coefficients of its operators include arbitrary functions of the dependent
variable $u$. It is easy to check that (5) has a remarkable property - a
function of any of its solutions is also a solution.

The discrete transformations for the equations (5) and (6) are quite obvious.
But for a comprehensive presentation we need to note that these equations are
invariant under time reflection $t^{\prime}=-t$, space reflections
$x^{\prime}_{\mu}=-x_{\mu}$, dependent variable reflection $u^{\prime}=-u$,
under a permutation group for the independent space variables. A hodograph
transformation $x^{\prime}_{1}=u$, $u^{\prime}=x_{1}$ that is a local
transformation for (6) will also be a discrete transformation for this
equation.

##  4 General Solutions

In consideration of the general solutions for the eikonal equations, we will
follow the algorithm, given in [2], where general solutions for the equations
(6) for $n=3$ are adduced. We generalize the relevant formulae for an
arbitrary number of space variables.

The algorithm in [2] uses contact transformations, nonlocal transformations
involving dependent and independent variables and first derivatives. The
general form of such contact transformations is as follows:

$\displaystyle x^{\prime}_{\mu}=\varphi_{\mu}(x_{\mu},u,u_{\mu})$
$\displaystyle u^{\prime}=\Phi(x_{\mu},u,u_{\mu})$ (17) $\displaystyle
u^{\prime}_{\mu}=\Phi_{\mu}(x_{\mu},u,u_{\mu}).$

Here $x^{\prime}_{\mu},u^{\prime}$ are new transformed variables,
$u^{\prime}_{\mu}$ are transformed first derivatives, $\varphi_{\mu}$, $\Phi$
and $\Phi_{\mu}$ are sufficiently smooth functions on their variables. Note
that indices at $\Phi_{\mu}$ are just indices, not derivatives.

We will first present in detail the procedure for construction of the general
solution for the simplest PDE of the type (6) for only two space variables:

$u_{1}^{2}+u_{2}^{2}=1.$ (18)

Note that if we apply the partial derivative operators to (18), we get the
system

$\displaystyle u_{1}u_{11}+u_{2}u_{12}=0,$ (19) $\displaystyle
u_{1}u_{21}+u_{2}u_{22}=0,$ (20)

from which we can see linear dependence of the above equations considered as
linear equations for the first derivatives of $u$, $u_{1}$ and $u_{2}$.

So, the determinant of the matrix of the second derivatives of $u$

$\displaystyle U={\left(\begin{matrix}u_{11}&u_{12}\\\ u_{21}&u_{22}\\\
\end{matrix}\right)}$

is equal to zero. However, here there will be two cases to be treated
separately:

  * •

${\rm rank}(U)=0$,

  * •

${\rm rank}(U)=1$.

In the first case, we have $u_{1}=u_{2}=const$, and $u$ linearly dependent on
$x_{1}$ and $x_{2}$:

$u=c_{1}x_{1}+c_{2}x_{2}+c,$ (21)

where $c_{1}$, $c_{2}$ and $c$ are arbitrary constants such that
$c_{1}^{2}+c_{2}^{2}=1$.

In the second case we obtain an implicit general solution depending on a
parameter by means of contact transformations.

We apply to (18) the following contact transformations in the case if
$u_{x_{1}}\neq 0$:

$\displaystyle H=x_{1}u_{x_{1}}-u,\qquad u_{x_{1}}=y_{1},\qquad x_{2}=y_{2};$
$\displaystyle H_{y_{1}}=x_{1},$ (22) $\displaystyle H_{y_{2}}=-u_{x_{2}};$

and respectively

$\displaystyle u=y_{1}H_{y_{1}}-H,\qquad H_{y_{1}}=x_{1},\qquad y_{2}=x_{2};$
$\displaystyle u_{x_{1}}=y_{1},$ (23) $\displaystyle u_{x_{2}}=-H_{y_{2}}.$

The resulting equation is a linearisable ordinary differential equation (ODE)

$H^{2}_{y_{2}}+y^{2}_{1}=1.$ (24)

We can linearize it as follows:

$H_{y_{2}}=\sqrt{1-y^{2}_{1}}$ (25)

A general solution of (25) may be taken as follows:

$H=-y_{2}\sqrt{1-y^{2}_{1}}-\Psi(y_{1}),$ (26)

where $\Psi(y_{1})$ is a sufficiently smooth arbitrary function on $y_{1}$.
Further it is convenient to introduce a parameter $\tau$ to be used instead of
$y_{1}$.

Transformations (23) allow going back to expressions for the function $u$, and
we obtain the following general solution for (18):

$\displaystyle u=\tau H_{\tau}-H=x_{1}\tau+x_{2}\sqrt{1-\tau^{2}}+\Psi(\tau),$
(27) $\displaystyle
x_{1}=H_{\tau}=x_{2}(\sqrt{1-\tau^{2}})^{-1}-\Psi_{\tau}(\tau),$

where $\Psi(\tau)$ is a sufficiently smooth arbitrary function on $\tau$.

We may rewrite the expressions (27) in a more convenient form:

$\displaystyle u=x_{1}\tau+x_{2}\sqrt{1-\tau^{2}}+\Psi(\tau),$ (28)
$\displaystyle x_{1}-x_{2}(\sqrt{1-\tau^{2}})^{-1}+\Psi_{\tau}(\tau)=0$

In our calculations that used the contact transformations, we assumed
$u_{x_{1}}\neq 0$, or, more generally, the rank of the matrix of the second
derivatives of the function $u$

$\displaystyle{\left(\begin{matrix}u_{11}&u_{12}\\\ u_{21}&u_{22}\\\
\end{matrix}\right)}$

is equal to one.

Let us consider the case of equation (6) with an arbitrary number of space
variables.

Similarly to (19), we can apply the partial derivative operators to (6) and
conclude linear dependence of the resulting equations considered as linear
equations for the first derivatives of $u$, and that the rank of the matrix of
the second derivatives of the function $u$ may take values from 0 to $n$.

In the case when the rank of the matrix of the second derivatives of the
function $u$ is equal to 0, all first derivatives $u_{\mu}$ are constants, and
$u$ is linearly dependent on $x_{\mu}$:

$u=c_{\mu}x_{\mu}+c,$ (29)

where $c_{\mu}$ and $c$ are arbitrary constants such that $c_{\mu}c_{\mu}=1$.

In the case when the rank of the matrix of the second derivatives (we will use
the notation $U$ for this matrix) of the function $u$ is equal to $n$, we can
apply a procedure very similar to what was used above for the equation (18),
and by means of the contact transformations (see [2] where it is considered
for $n=3$)

$\displaystyle H=x_{a}u_{x_{a}}-u,\qquad y_{a}=u_{x_{a}},\qquad x_{0}=y_{0};$
$\displaystyle H_{y_{a}}=x_{a},$ (30) $\displaystyle H_{y_{0}}=-u_{x_{0}};$

obtain an $n$-parameter general solution (${\rm rank}(U)=n$):

$\displaystyle
u=-x_{a}\tau_{a}+x_{0}\sqrt{1+\tau_{d}\tau_{d}}+\Psi(\tau_{d}),$ (31)
$\displaystyle
x_{a}-x_{0}\tau_{a}(\sqrt{1+\tau_{d}\tau_{d}})^{-1}-\Psi_{\tau_{d}}(\tau_{a})=0$

where $\Psi(\tau_{a})$ is a sufficiently smooth arbitrary function on $n$
parameters $\tau_{a}$.

Further general solutions for every rank $k=1,n-1$ will be written depending
on $k$ parameters.

In each such case we may introduce $n-k$ new functions $w_{m}$, were $m$ are
indices taking values from 1 to $n-k$. Without any limitation of generality,
we can assume

$\displaystyle u_{0}=\sqrt{1+u_{b}u_{b}+w_{m}w_{m}},$ (32) $\displaystyle
u_{k+m}=w_{m}(u_{b}),$

where the indices $m$ run from 1 to $n-k$, and the indices and symbols of
derivatives $b$ run from 1 to $k$, symbols of derivatives $k+m$ run from $k+1$
to $n$. Thus, we can take for every rank $k=1,n-1$ of the matrix $U$ $n-k$ new
functions $u_{k+m}=w_{m}(u_{b})$ depending on the first $k$ derivatives
$u_{b}$. General solutions depending on $k$ parameters may be taken as
follows:

$\displaystyle
u=-x_{b}\tau_{b}+x_{0}\sqrt{1+\tau_{d}\tau_{d}+w_{m}w_{m}}+w_{m}x_{k+m}+\Psi(\tau_{d}),$
(33) $\displaystyle
x_{b}-x_{0}(\sqrt{1+\tau_{d}\tau_{d}+w_{m}w_{m}})^{-1}-w_{m\tau_{b}}x_{k+m}-\Psi_{\tau_{b}}(\tau_{d})=0$

The expressions in (33) imply a maybe unusual summation over $m$ from 1 to
$n-k$.

Proof of the fact that (33) is a rank $k$ general solution of (6) for $n$
space variables can be done by means of a contact transformation

$\displaystyle H=x_{b}u_{x_{b}}-u,\qquad y_{k+m}=u_{x_{k+m}},\qquad
x_{0}=y_{0};\qquad x_{b}=y_{b};$ $\displaystyle
H_{y_{b}}=x_{b}H_{y_{0}}=-u_{x_{0}},$ (34) $\displaystyle
H_{y_{k+m}}=-u_{x_{k+m}}$

applied similarly to the procedures in [2].

A general solution for (5) and for arbitrary number of independent variables
was found in [9].

##  5 A Note of the General Solution and the Known Symmetry Exact Solutions
for Eikonal Equations

A well-known symmetry solution - the radial solution of the equation (6)

$u=\sqrt{x_{\mu}x_{\mu}}$ (35)

has the rank $n$ and corresponds to the following solution represented by the
general solution formula (31)

$\displaystyle u=x_{a}\tau_{a}-x_{0}\sqrt{1+\tau_{b}\tau_{b}},$ (36)
$\displaystyle x_{a}-x_{0}\tau_{a}(\sqrt{1+\tau_{b}\tau_{b}})^{-1}=0,$

The parameters $\tau_{a}$ will take values

$\tau_{a}=\frac{x_{a}}{\sqrt{x_{\mu}x_{\mu}}}.$ (37)

The formula (31) includes arbitrary functions, so all solutions (31) cannot be
equivalent to symmetry solutions. In general (not always), formulas for the
general solutions give wider classes of solutions than the Lie symmetry
method.

As to the zero-rank general solution (29), all such solutions are symmetry
solutions. For more symmetry solutions of low-dimensional eikonal equations
see e.g. [6] and [10].

##  6 Some Relations Among the Eikonal Equations and the Hamilton-Jacobi
Equation

As was established by Sophus Lie [11], all scalar first-order partial
differential equations can be transformed into one another by some suitable
contact transformations.

A hodograph transformation

$u=y_{0},\qquad x_{0}=w,\qquad x_{a}=y_{a}$ (38)

will get from a system of two eikonal equations

$\displaystyle u_{\mu}u_{\mu}=0,$ (39) $\displaystyle v_{\mu}v_{\mu}=0,$

a system of an eikonal equation without a time variable (with only $n$ space
variables) and a Hamilton-Jacobi equation

$\displaystyle w_{y_{a}}w_{y_{a}}=1,$ $\displaystyle
v_{y_{a}}v_{y_{a}}=2v_{y_{0}},$ (40)

##  7 New Results, Conclusions and Further Work

The new results of this paper are a proof of maximal Lie symmetries for the
eikonal equations (5) and (6) and a general solution for (6) for arbitrary
number of space variables.

The eikonal equation is important in many fields, and in symmetry analysis it
often appears in reduction conditions for many higher-order equations, in
particular, for all equations invariant under the orthogonal and Poincaré
groups ([12], [13]). Its general solution is a first step to find general
solutions of reduction conditions.

There are also multiple papers searching for numerical solutions of the
boundary problems for eikonal equations. On my opinion, it would be
interesting to look for these using the general solutions.

Further work will include study of the coupled eikonal equations in higher
dimensions (for the case of 2 space variables see [14]), relevant conditional
and hidden symmetries and looking for higher-order equations with the same or
similar symmetry properties as the eikonal equation.

##  8 Acknowledgements

The first and foremost my acknowledgements go to the Armed Forces of Ukraine
and to the Territorial Defence Forces of Ukrainian Regions due to whom I am
alive and is able to work.

Please remember that Russia is an aggressor country and still plans to kill
all Ukrainians not going to be their subordinates, and Russian scientists
contribute to killings despite words about peace from a tiny portion of them;
if even by trying to persuade naïve people that Russia is a civilised country.

I would like also to thank the Institute of Mathematics of the Polish Academy
of Sciences for their hospitality and grant support, to the National Academy
of Sciences of the USA and the National Centre of Science of Poland for their
grant support.

Research was supported by Narodowe Centrum Nauki, grant number
2017/26/A/ST1/00189.

## References

  * [1] Boyer C.P., Kalnins E.G., and Miller Jr. W., Symmetry and separation of variables for the Hamilton–Jacobi equation $W^{2}_{t}-W^{2}_{x}-W^{2}_{y}=0$, J. Math. Phys. 19, 200–211 (1978).
  * [2] Fushchych W.I, Zhdanov R.Z., Revenko I.V., Compatibility and solution of non-linear d’Alembert-Hamilton equations, Preprint 90.39, Institute of Mathematics of Acad.Sci. Ukr. SSR, Kyiv, 1990, 65 p.
  * [3] Wikipedia: Eikonal equation, Accessed on 30.11.2022.
  * [4] Popovych R.O., Yehorchenko I.A., Group Classification of Generalized Eikonal Equations. Ukrainian Mathematical Journal 53, 1841–-1850 (2001).
  * [5] Fushchich W.I., Shtelen W.M., The symmetry and some exact solutions of the relativistic eikonal equation. Lett. Nuovo Cimento 34, 498–-502 (1982).
  * [6] Fushchych W.I., Serov N.I., The symmetry and some exact solutions of the nonlinear many-dimensional Liouville, d’Alembert and eikonal equation, J. Phys. A 16, 3645–3656 (1983).
  * [7] Ghanam R., Thompson G., Symmetries of the Eikonal equation, Communications in Nonlinear Science and Numerical Simulation 60, 137–144 (2018).
  * [8] Olver P.J., Application of Lie groups to differential equations, New York, Springer Verlag, 1987.
  * [9] Collins C.B., Null field solutions of the wave equation and certain generalizations, J. Math. Phys. 24, 22-28 (1983).
  * [10] Fedorchuk V.M., Fedorchuk V.I., On classification of symmetry reductions for the eikonal equation. Symmetry 8(6), 51 (2016).
  * [11] Lie Sophus, Vorlesungen über continuierliche Gruppen mit geometrischen und anderen Anwendungen. BG Teubner, 1893.
  * [12] Fushchych W.I., Zhdanov R.Z., Revenko I.V., On the general solution of the d’Alembert equation with a nonlinear eikonal constraint and its applications, J. Math. Phys. 36, 7109–7127 (1995).
  * [13] Yehorchenko I.A., Ansatzes and exact solutions for nonlinear Schrödinger equations, arXiv:1412.1889 [math-ph] (2014).
  * [14] Yehorchenko I.A., General solution for a coupled system of eikonal equations in two space variables, arXiv:1712.01948 (2017).


# Electron and proton energization in 3D reconnecting current sheets in
semirelativistic plasma with guide magnetic field

Gregory R. Werner Center for Integrated Plasma Studies, Physics Department,
390 UCB, University of Colorado, Boulder, CO 80309, USA  Dmitri A. Uzdensky
Center for Integrated Plasma Studies, Physics Department,
390 UCB, University of Colorado, Boulder, CO 80309, USA

###### Abstract

Using 3D particle-in-cell simulation, we characterize energy conversion, as a
function of guide magnetic field, in a thin current sheet in semirelativistic
plasma, with relativistic electrons and subrelativistic protons. There,
magnetic reconnection, the drift-kink instability (DKI), and the flux-rope
kink instability all compete and interact in their nonlinear stages to convert
magnetic energy to plasma energy. We compare fully 3D simulations with 2D in
two different planes to isolate reconnection and DKI effects. In zero guide
field, these processes yield distinct energy conversion signatures: ions gain
more energy than electrons in 2D$xy$ (reconnection), while the opposite is
true in 2D$yz$ (DKI), and the 3D result falls in between. The flux-rope
instability, which occurs only in 3D, allows more magnetic energy to be
released than in 2D, but the rate of energy conversion in 3D tends to be
lower. Increasing the guide magnetic field strongly suppresses DKI, and in all
cases slows and reduces the overall amount of energy conversion; it also
favors electron energization through a process by which energy is first stored
in the motional electric field of flux ropes before energizing particles.
Understanding the evolution of the energy partition thus provides insight into
the role of various plasma processes, and is important for modeling radiation
from astrophysical sources such as accreting black holes and their jets.

magnetic reconnection — plasmas — acceleration of particles

## References

  *

##  1 Introduction

In many plasma environments, magnetic energy is converted to plasma energy at
thin current sheets (CSs), likely via magnetic reconnection (e.g., Zweibel &
Yamada, 2009). Reconnection occurs in diverse regimes and, even in 2D,
generates such complexity that we owe much of our understanding to numerical
simulation. Unfortunately, 3D simulation studies are far costlier than 2D,
motivating efforts to understand 3D CS evolution in terms of 2D processes
whenever possible.

A thin CS, considered in 2D [specifically in “2D$xy$”, the $x$-$y$ plane
perpendicular to the current $j_{z}$ supporting a magnetic field $\hat{\bf
x}B_{0}\textrm{tanh}(y/\delta)$ reversing over thickness $2\delta$], will
typically be most unstable to tearing, leading, in its nonlinear stage, to
reconnection and rapid magnetic energy release. In 2D$yz$, the same CS may
endure the drift-kink instability (DKI), which can also release magnetic
energy rapidly (e.g., Pritchett et al., 1996)—if the separation between
electron and ion scales is small (Daughton, 1999; Hesse & Birn, 2000; Scholer
et al., 2003), as in relativistic pair plasma (Zenitani & Hoshino, 2005, 2008;
Yin et al., 2008; Liu et al., 2011; Cerutti et al., 2014; Guo et al., 2014,
2015). In 3D, these instabilities can compete and interact. Although 2D-like
reconnection may dominate DKI in 3D (Kagan et al., 2013; Sironi & Spitkovsky,
2014; Guo et al., 2014, 2015; Cerutti et al., 2014; Werner & Uzdensky, 2017;
Guo et al., 2021), Werner & Uzdensky (2021) highlighted how the DKI can
dramatically slow reconnection. Furthermore, in 3D only, flux ropes formed by
reconnection can decay via the magnetohydrodynamic kink instability (Kagan et
al., 2013; Markidis et al., 2014; Zhang et al., 2021; Schoeffler et al.,
2023), releasing additional magnetic energy (Werner & Uzdensky, 2021).

Using particle-in-cell (PIC) simulation, this work characterizes energy
conversion as a function of guide magnetic field, in 3D CSs in magnetically-
dominated, collisionless electron-proton plasma in the semirelativistic
regime. In semirelativistic plasma, where electrons are ultrarelativistic and
protons subrelativistic, relativistic effects diminish the scale separation
between (e.g., gyroradii of) electrons and ions (protons), making DKI
consequential. CSs in this regime likely play important roles in releasing
magnetic energy in astrophysical systems such as jets, coronae, and accretion
flows around black holes; characterizing electron and ion energization is
crucial for modeling flaring reconnection-powered emission from accreting
black holes (Chael et al., 2018; Dexter et al., 2020; Ressler et al., 2020,
2023; Scepi et al., 2022; Hankla et al., 2022).

Until now (to our knowledge), no 3D PIC studies of electron-ion reconnection
(or CS evolution) have characterized energy conversion in the nonradiative
semirelativistic regime, although several 2D reconnection studies have
(Melzani et al., 2014b, a; Guo et al., 2016; Rowan et al., 2017, 2019; Werner
et al., 2018; Ball et al., 2018, 2019).

After describing the simulation setup (§2), we will focus on the different
energy components, discussing the release of magnetic energy $U_{Bxy}$ (§3),
the increase in guide field energy $U_{Bz}$ (§4), the electric field energy
$U_{E}$ (§5), and the electron and ion energies $U_{e}+U_{i}=U_{\rm ptcl}$
(§6). We summarize in §7.

##  2 The simulations

Simulations were run with the Zeltron electromagnetic PIC code (Cerutti et
al., 2013), initialized with semirelativistic electron-proton plasma in a
standard Harris sheet configuration with uniform background plasma (similar to
Werner et al., 2018)—initial physical and numerical parameters are listed in
Table 1. The simulation frame (with zero initial electric field) is the zero-
momentum frame, where ions drift slowly and electrons carry the current
supporting the field reversal.

Table 1:  Initial plasma and numerical simulation parameters.  $\mu\equiv m_{i}/m_{e}=1836$ | ion/electron mass ratio
---|---
$B_{0}$ | asymptotic upstream (reversing) magnetic field
$n_{b}=n_{be}=n_{bi}$ | background electron/ion density
$T_{b}=18.36m_{e}c^{2}=0.01m_{i}c^{2}$ | background electron/ion temperature
$\sigma_{i}\equiv B_{0}^{2}/(4\pi n_{b}m_{i}c^{2})=0.5$ | background ion cold magnetization
$\sigma_{e}\equiv B_{0}^{2}/(4\pi n_{b}m_{e}c^{2})=\mu\sigma_{i}=918$ | background electron cold magnetization
$\delta=0.33\sigma_{i}\rho_{i0}$ | Harris CS half thickness ($\rho_{i0}\equiv m_{i}c^{2}/eB_{0}$)
$\mathbf{B}(y)=\hat{\bf z}B_{gz}-\hat{\bf x}B_{0}\textrm{tanh}(y/\delta)$ | magnetic field profile
$-0.6c\hat{\bf z}$, $+0.0004c\hat{\bf z}$ | electron, ion bulk drift velocities in CS
$\eta\equiv n_{d,\rm max}/n_{b}=5$ | CS overdensity
$n_{d}(y)=\eta n_{b}\textrm{cosh}^{-2}(y/\delta)$ | CS electron/ion density profile
$T_{d}=51m_{e}c^{2}=0.28m_{i}c^{2}$ | CS electron/ion (co-moving) temperature
$\beta_{be}=\beta_{bi}\equiv 8\pi n_{b}T_{b}/B_{0}^{2}=0.04$ | upstream plasma beta for electrons/ions
$v_{A}=0.57c$ | upstream Alfvén velocity for $B_{gz}=0$
$L_{y}\equiv L=55.3\sigma_{i}\rho_{i0}$ | system size (in $y$)
$L_{x}/L=1$ | system $x/y$ aspect ratio (for 2D$xy$ and 3D)
$L_{z}/L=\textrm{max}(1,B_{gz}/B_{0})$ | system $z/y$ aspect ratio (for 2D$yz$ and 3D)
$T=20L/c$ | simulation duration
$\Delta y=\Delta z=0.036\sigma_{i}\rho_{i0}$ | grid cell size in $y$ and $z$
$\Delta x=0.058\sigma_{i}\rho_{i0}$ | grid cell size in $x$
$\Delta t=0.64\Delta y/c$ | timestep (the 3D Courant-Friedrichs-Lewy maximum)
$M_{\rm ppcs}=10$ (3D), 160 (2D) | macroparticles per cell per species (4 species)
$M_{\rm ppc\phantom{s}}=40$ (3D), 640 (2D) | total macroparticles per cell
Table 2:  Initial length scales normalized to $\sigma_{i}\rho_{i0}=\sigma_{e}\rho_{e0}$, where $\rho_{s0}=m_{s}c^{2}/eB_{0}$, $\sigma_{s}=B_{0}^{2}/(4\pi n_{bs}m_{s}c^{2})$ $\Delta y=\Delta z$ | 0.036 | grid cell size in $y$, $z$
---|---|---
$\Delta x$ | 0.058 | grid cell size in $x$
$\rho_{be}$ | 0.060 | upstream electron gyroradius ($B_{gz}=0$)
$\lambda_{De}=\lambda_{Di}$ | 0.14 | upstream electron/ion Debye length
$d_{e}$ | 0.25 | upstream relativistic electron skin depth
$\delta$ | 0.33 | CS half-thickness
$\rho_{bi}$ | 0.35 | upstream ion gyroradius ($B_{gz}=0$)
$\sigma_{e}\rho_{e0}$ | 1.006 | twice the energized electron gyroradius, if all magnetic energy went to electrons
$d_{i}$ | 1.406 | upstream nonrelativistic ion skin depth
$L$ | 55.306 | system size in $x$, $y$, and $z$ (except: if $B_{gz}/B_{0}>1$, then $L_{z}/L=B_{gz}/B_{0}$)

The initial CS is perpendicular to $y$, with current in the $+z$ direction and
reversing magnetic field in the $\pm x$ directions (approaching value $B_{0}$
upstream). We add uniform guide field $B_{gz}\hat{\bf z}$, exploring several
strengths $B_{gz}/B_{0}\in\\{0,0.25,0.5,1,2,4\\}$ (except $B_{gz}/B_{0}=4$ was
not run in 3D, due to expense). We also explore three dimensionalities: 2D$xy$
(2D simulation ignoring $z$—classic 2D reconnection), 2D$yz$ (ignoring $x$,
for studying DKI), and full 3D. Also: in one case, 2D$yz$ and $B_{gz}=0$, we
ran an ensemble of 3 simulations, identical except for differing randomization
of initial particles.

Table 2 lists length scales of interest in terms of
$\sigma_{i}\rho_{i0}=\sigma_{e}\rho_{e0}$, where $\rho_{s0}\equiv
m_{s}c^{2}/(eB_{0})$; $\sigma_{e}\rho_{e0}$ is the gyroradius of an
“energized” electron with energy $2B_{0}^{2}/(8\pi n_{be})$ in perpendicular
field $B_{0}$.

The system size, containing a single CS, is $L\equiv
L_{y}=55.3\sigma_{i}\rho_{i0}$ and $L_{x}=L_{z}=L$ [except when
$B_{gz}>B_{0}$, we increase $L_{z}=(B_{gz}/B_{0})L_{y}$; and 2D simulations
lack either the $x$ or $z$ dimension]. Boundary conditions are periodic in $x$
and $z$, and conducting/particle-reflecting in $y$. The size $L$ is large with
respect to initial kinetic scales and CS thickness [e.g., $L/(2\delta)\approx
83$], and $L>20\sigma_{i}\rho_{i0}$ puts it in the “large system regime” as
defined by Werner et al. (2016) (for relativistic pair plasma in 2D$xy$).

We empirically found that using anisotropic grid cells with $\Delta
x/1.6=\Delta y=\Delta z=\sigma_{i}\rho_{i0}/28$ avoids numerical instabilities
while yielding the best computational performance. We similarly found that 160
macroparticles/species/cell were needed in 2D to avoid unphysical electron-ion
energy exchange at late times, whereas, in 3D, 10 macroparticles/species/cell
sufficed (with 4 species: background and CS electrons and ions).

No initial perturbation is used to kickstart reconnection, because such a
perturbation (if uniform in $z$) can significantly suppress 3D effects (Werner
& Uzdensky, 2021).

##  3 Qualitative evolution and transverse magnetic energy released

(a)(b)

Figure 1:  (a) $U_{Bxy}(t)$ for selected simulations ($B_{gz}$ is labeled for
each). Three (differently-seeded) 2D$yz$ simulations with $B_{gz}=0$
demonstrate the enhanced stochastic variability of DKI. (b) $|\Delta U_{Bxy}|$
(large symbols) and $\Delta U_{\rm ptcl}$ (small symbols) versus $B_{gz}$,
normalized to $U_{Bxy0}$, the initial transverse magnetic energy.

Over time, transverse magnetic energy $U_{Bxy}$ is released to other forms
($U_{Bxy}$ excludes energy in guide field $B_{z}$, which cannot be
released—see §4). Figure 1(a) shows the evolution of $U_{Bxy}(t)$ for selected
simulations. Figure 1(b) shows the net loss in $U_{Bxy}$ over $T=20L/c$ as a
function of $B_{gz}$, for all simulations (large symbols). Ultimately this
energy $\Delta U_{Bxy}$ goes almost entirely to particles ($\Delta U_{\rm
ptcl}$, small symbols; cf. §6), with the small difference $|\Delta
U_{Bxy}|-\Delta U_{\rm ptcl}$ comprising mainly $\Delta U_{Bz}$, increased
guide field energy (see §4).

In 2D$xy$ (Fig. 1a, blue triangles), magnetic energy decreases rapidly due to
magnetic reconnection until saturation after 5–10$\>L/c$. During that time we
observe the familiar behavior of plasmoid-dominated reconnection. Reconnection
pumps magnetic energy and flux from the upstream into plasmoids (magnetic
islands), continually growing the total plasmoid volume and the magnetic
energy therein. Apart from merging, plasmoids are permanent (stable) in 2D.
Eventually (after 5–10$\>L/c$), in a closed system, the largest plasmoid
becomes so large ($\sim\\!\\!L$) that it halts reconnection—even though
considerable transverse magnetic energy may remain upstream as well as in the
plasmoid.

In 2D$yz$ (Fig. 1, red diamonds), tearing and reconnection are forbidden, but
the nonlinear DKI can rapidly release magnetic energy. We observe
qualitatively similar behavior—the rippling of the CS, sometimes leading to a
catastrophic folding-over on itself—as reported in nonrelativistic electron-
ion plasma with artificially low $m_{i}/m_{e}$ and in relativistic pair plasma
(Ozaki et al., 1996; Pritchett & Coroniti, 1996; Pritchett et al., 1996; Zhu &
Winglee, 1996; Zenitani & Hoshino, 2007; Cerutti et al., 2014; Werner &
Uzdensky, 2021). The nonlinear DKI effectively thickens the CS, significantly
slowing but not fully halting energy conversion. A DKI mode may enter the
nonlinear stage when its rippling amplitude exceeds its rippling wavelength,
but—like waves breaking in the ocean—whether and where this occurs has a
substantial element of randomness (see Figs. 28, 29, 35, in Werner & Uzdensky,
2021). While a large statistical study is beyond the scope of this work, this
stochastic variation is shown via the three 2D$yz$ simulations with $B_{gz}=0$
(red diamonds and thin dashed lines) in Fig. 1(a).

In 3D (Fig. 1, magenta circles), $U_{Bxy}$ can be released by reconnection and
DKI. In addition, flux ropes (3D plasmoids) formed by reconnection may decay
via the flux-rope kink instability, releasing still more magnetic energy than
in 2D (Fig. 1b). Because kinking flux ropes decay, reconnection will not halt
as in 2D$xy$, although CS thickening (due to DKI and flux-rope instability)
slows both DKI and reconnection; all the 3D simulations in Fig. 1(a) show
continued but slow energy conversion beyond $20\>L/c$ (Werner & Uzdensky,
2021).

In Fig. 1(b) we see that, for $B_{gz}=0$, all dimensionalities yield
comparable $|\Delta U_{Bxy}|\sim 0.2$–$0.3U_{Bxy0}$. While a small guide field
$B_{gz}/B_{0}=1/4$ has only a small effect in 2D$xy$ and 3D, it strongly
inhibits energy conversion in 2D$yz$.

Increasing $B_{gz}$ has two main effects, in all cases: it slows energy
conversion, and decreases the net magnetic energy release $|\Delta U_{Bxy}|$.
A guide field of $B_{gz}/B_{0}\gtrsim 1/2$ stabilizes the DKI in 2D$yz$ (cf.
Zenitani & Hoshino, 2008; Schoeffler et al., 2023), preventing any release of
$U_{Bxy}$. For 2D$xy$ and 3D, increasing $B_{gz}/B_{0}$ from 0 to 1 slows
energy conversion and reduces $|\Delta U_{Bxy}|$ by roughly half; by
$B_{gz}/B_{0}=2$, only about $0.1U_{Bxy0}$ is released (this depends on the
aspect ratio: for $L_{y}\gtrsim L_{x}$, we expect $|\Delta
U_{Bxy}|/U_{Bxy0}\sim 0.1L_{x}/L_{y}$). Previous 2D$xy$ studies showed that
$B_{gz}$ slows reconnection by reducing the Alfvén speed $v_{A,x}$ in the
outflow direction (Liu et al., 2014; Werner & Uzdensky, 2017); also, $B_{gz}$
resists plasmoid compression, halting reconnection at smaller $|\Delta
U_{Bxy}|$. These trends hold in 3D as well, except that increasing
$B_{gz}/B_{0}$ from 0 to $1/4$ releases more energy, possibly because $B_{z}$
suppresses the DKI and prevents interference with reconnection, or possibly
this is stochastic variation; we leave exploration of this subtrend to future
work.

Overall, after $20\>L/c$, $\Delta U_{Bxy}$ is about 15% higher in 3D than in
2D$xy$ because 3D evolution accesses all the 2D dissipation channels plus the
flux-rope instability; the difference between 2D$xy$ and 3D narrows slightly
with higher $B_{gz}$ as guide field increases uniformity in $z$. We note that
with more time, $|\Delta U_{Bxy}|$ would increase in 3D, but not in 2D$xy$.

##  4 Energy gain in guide field

(a) |  | (b)
---|---|---
(c) |  | (d)
Figure 2:  (a) Net gain in $U_{Bz}$ versus $B_{gz}$, in 2D$xy$ (blue
triangles), 3D (magenta circles), and 2D$yz$ (red diamonds showing zero gain).
(b) Net gain in $U_{E}(B_{gz})$. (c) For 2D$xy$ only, $U_{E}(t)$ vs. time, for
different $B_{gz}$, normalized to the ultimately-released magnetic energy [not
shown: in 2D$yz$ and 3D, $U_{E}(t)$ remains relatively small]. $U_{E}$
dissipates only after most of $U_{Bxy}$ has dissipated, as indicated by the
circles placed where $\Delta U_{Bxy}(t)=0.75\Delta U_{Bxy}(20L/c)$. (d) For
2D$xy$, the plasmoid $E_{y}$ as a function of $B_{gz}$, estimated from Eqs.
(3)–(4) (solid line), compared with the maximum $E_{y}$ measured in
simulations (triangles). (To reduce sensitivity to noise, the maximum was
measured over time of the 99th percentile of $|E_{y}|$.)

In our simulations, the guide-field energy $U_{Bz}$ can only increase from its
initial value; the $B_{z}$-flux through any $x$-$y$ plane is conserved and
$B_{z}$ is initially uniform, yielding the lowest $U_{Bz}$ for the fixed flux.
Figure 2(a) shows the net gain $\Delta U_{Bz}$ versus $B_{gz}$ for each
dimensionality. In the 2D$yz$ case, $\Delta U_{Bz}\approx 0$, and for
$B_{gz}=0$ neither 2D$xy$ nor 3D show noticeable gain. However, for
$B_{gz}>0$, reconnection outflows compress plasma and guide field in flux
ropes, increasing $U_{Bz}$. The gain $\Delta U_{Bz}$ increases with
$B_{gz}/B_{0}\lesssim 0.5$ because there is more $B_{z}$ to be compressed;
however, stronger $B_{gz}$ rivals the plasma in resisting compression (with
adiabatic index 2, greater than the adiabatic index 4/3 or 5/3 of relativistic
or nonrelativistic plasma), and $\Delta U_{Bz}$ peaks around $0.04\,U_{Bxy0}$
at $B_{gz}/B_{0}\sim 0.5$ before decreasing. The similarity between 2D$xy$ and
3D suggests reconnection occurs to similar extents—though slower in
3D—resulting in similar compression of flux ropes, with energy remaining in
$U_{Bz}$ even as flux ropes decay in 3D.

We note that the contribution of the Hall quadrupole $B_{z}$ field is
negligible on global energy scales, because it is localized to small regions
around X-lines; in contrast, the $U_{Bz}$ in plasmoids becomes globally
significant because it grows with plasmoids, which approach the system size.

##  5 Energy gain and subsequent loss in electric field

The net gain in electric field energy $U_{E}$ is negligible in all cases—less
than $0.004\,U_{Bxy0}$ (Fig. 2b). For 2D$yz$ and 3D, $U_{E}(t)$ is small at
all times. However (see Fig. 2c), in 2D$xy$ with strong guide field,
$U_{E}(t)$ builds up substantially over time before declining (after most
reconnection has occurred); this effect increases with $B_{gz}$. For
$B_{gz}/B_{0}=2$, 17% of the ultimately-released magnetic energy resided at
one time (around $t=5L/c$) in $U_{E}$, and for $B_{gz}/B_{0}=4$, over 30%.

This effect—by which substantial $U_{Bxy}$ is first converted to $U_{E}$
before energizing particles—has not been reported before, to our knowledge. We
find that (when $B_{gz}/B_{0}\gtrsim 0.5$) $U_{E}$ is dominated by $E_{y}$ in
plasmoids; this is simply the motional electric field $E_{m,y}\approx
v_{x}B_{z}/c$ due to $B_{z}$ in a plasmoid moving with relativistic speed
$v_{x}$ along the layer. This effect should appear in any 2D reconnection
simulation with guide field and relativistic outflows, even in pair plasma.
After the last major plasmoid merger, as reconnection is ending, only one
large but stationary plasmoid remains, and thus the motional electric field
ultimately decays.

We can estimate how $E_{m,y}\approx v_{x}B_{z}/c$ scales with guide field,
assuming that $B_{z}\sim B_{gz}$ and that plasmoid velocity $v_{x}$ scales
with $v_{A,x}$, the Alfvén velocity projected along $x$. Since the typical
reconnection electric field is $E_{\rm rec}\sim 0.1v_{A,x}B_{0}/c$ (for fast
collisionless reconnection), $E_{m,y}/E_{\rm rec}\sim 10B_{gz}/B_{0}$, which
exceeds unity for even small $B_{gz}$. In absolute terms,
$v_{A,x}=c[B_{0}^{2}/(4\pi h+B_{0}^{2}+B_{gz}^{2})]^{1/2}$, where $h$ is the
relativistic plasma enthalpy density including rest-mass (Liu et al., 2014;
Sironi et al., 2016; Werner & Uzdensky, 2017), and thus

$\displaystyle E_{m,y}(B_{gz})$ $\displaystyle\sim$
$\displaystyle\sqrt{\frac{B_{gz}^{2}B_{0}^{2}}{4\pi
h+B_{0}^{2}+B_{gz}^{2}}}\rightarrow\left\\{\begin{array}[]{c@{\textrm{ for
}}l}E_{m,y}\propto B_{gz}&B_{gz}^{2}\ll B_{0}^{2}+4\pi h\\\ E_{m,y}\sim
B_{0}&B_{gz}^{2}\gtrsim B_{0}^{2}+4\pi h\end{array}\right..$ (3)

(In this paper, $4\pi h\approx 4\pi
n_{b}m_{i}c^{2}=B_{0}^{2}/\sigma_{i}=2B_{0}^{2}$.)

This estimate of $E_{m,y}(B_{gz})$ roughly explains the $B_{gz}$-dependence of
$E_{y}$ in our 2D$xy$ simulations; Fig. 2(d) shows the maximum observed
$E_{y}$ compared with

$\displaystyle E_{y,\rm max}(B_{gz})$ $\displaystyle=$
$\displaystyle\sqrt{E_{y0}^{2}+[\alpha E_{m,y}(B_{gz})]^{2}},$ (4)

which adds in quadrature the motional field $E_{m,y}(B_{gz})$ given by Eq. (3)
and $E_{y0}\equiv E_{y,\rm max}(0)\approx 0.14B_{0}$, an estimate of $E_{y}$
from other sources (which we measure in the simulation with $B_{gz}=0$ and
assume it remains similar for $B_{gz}>0$), with empirical factor
$\alpha\approx 0.6$ accounting for, e.g., $v_{x}<v_{A,x}$.

Notably, when $U_{E}$ dissipates at the end of reconnection, it preferentially
energizes electrons, as described next.

##  6 Energy gain in electrons and ions

(a) |  | (b)
---|---|---
(c) |  | (d)
Figure 3:  The net gain in (a) electron and (b) ion energy, versus $B_{gz}$,
normalized by initial transverse magnetic energy $U_{Bxy0}$, over the full
$T=20\>L/c$. (c) The electron gain fraction $q_{e}$ versus $B_{gz}$. (d)
$q_{e}$ for each third (in terms of $\Delta U_{Bxy}$) of each simulation (the
smallest markers, shifted left, show $q_{e}$ in the first third; the largest,
shifted right, show the last third). 2D$yz$ cases with $B_{gz}/B_{0}>0.5$,
hence negligible energization, are omitted.

In all dimensionalities, increasing guide field generally suppresses both
electron and ion energization relative to the available magnetic energy
$U_{Bxy0}$ (Fig. 3a,b). In 2D$yz$, both electron and ion net gains ($\Delta
U_{e}$ and $\Delta U_{i}$) drop precipitously with $B_{gz}$. In 2D$xy$ and 3D,
energy gains also fall, but not to zero, and $\Delta U_{i}$ falls much faster
than $\Delta U_{e}$. As $B_{gz}/B_{0}$ increases from 0.25 to 2, $\Delta
U_{e}$ falls from $0.09U_{Bxy0}$ to $0.05U_{Bxy0}$ in 2D$xy$, and from
$0.16U_{Bxy0}$ to $0.09U_{Bxy0}$ in 3D; in contrast, $\Delta U_{i}$ falls from
about $0.17U_{Bxy0}$ to just $0.02U_{Bxy0}$ (almost an order of magnitude!),
and is nearly identical in 2D$xy$ and 3D. The additional magnetic energy
released in 3D (compared with 2D$xy$) goes mostly to electrons. At low
$B_{gz}$, this may be due to the DKI and/or the flux-rope kink instability; at
higher $B_{gz}$, DKI is quenched, suggesting that the flux-rope kink
preferentially energizes electrons. The electron energy gain fraction
$q_{e}\equiv\Delta U_{e}/(\Delta U_{e}+\Delta U_{i})$ grows with $B_{gz}$
(Fig. 3c), even though $\Delta U_{e}$ decreases. In 2D$xy$, $q_{e}\approx
0.35$ at $B_{gz}=0$, increasing to $q_{e}\approx 0.5$ (equipartition) around
$B_{gz}/B_{0}=0.75$, and reaching $q_{e}>0.8$ by $B_{gz}/B_{0}=4$. In 3D,
$q_{e}(B_{gz})$ is higher than in 2D$xy$ by about 0.1 (though this gap narrows
as $B_{gz}$ increases), reaching equipartition near $B_{gz}/B_{0}=0.25$, and
$q_{e}\approx 0.7$ at $B_{gz}/B_{0}=1$. In 2D$yz$, $q_{e}$ is $\sim 0.2$
higher than in 3D, but we emphasize that, for $B_{gz}/B_{0}\geq 0.5$, $\Delta
U_{e}$ in 2D$yz$ is tiny compared with 2D$xy$ and 3D (hence hard to measure).

The 3D case with $B_{gz}=0$ is exceptional, converting less energy than higher
$B_{gz}/B_{0}=0.25$, but the difference is within the stochastic variation
observed in 2D$yz$ for $B_{gz}=0$ (cf. Figs. 1b and 3a,b).

Similar values of $q_{e}$ were previously measured in 2D$xy$ PIC reconnection
simulations in the (nearly) semirelativistic, plasmoid-dominated regime. For
$B_{gz}=0$, $q_{e}\simeq 0.25$–$0.35$ was found in similar simulations (Rowan
et al., 2017, 2019; Werner et al., 2018; Zhang et al., 2021) and also in
nonrelativistic laboratory experiment (e.g., Yamada et al., 2014). Similar
increases in $q_{e}$ with $B_{gz}$ were also reported in 2D: e.g., Rowan et
al. (2019) found $q_{e}$ increasing to equipartition around
$B_{gz}/B_{0}\approx 1$ and up to 0.92 for $B_{gz}/B_{0}=6$; with reduced mass
ratio $\mu=25$, Melzani et al. (2014a) found $q_{e}$ increasing to 0.67 at
$B_{gz}/B_{0}=1$. There have, however, been no 3D PIC simulations in a
quantitatively comparable regime (to our knowledge) that have measured
$q_{e}$.

Although a detailed analysis of energization mechanisms is beyond the scope of
this Letter, the time-dependence of $q_{e}$ offers some insight. Figure 3(d)
plots $q_{e}$ separately for each third (in terms of $\Delta U_{Bxy}$, not
time) of each simulation. In 2D$yz$, we observe no significant time-dependence
of $q_{e}$. For 2D$xy$ and strong $B_{gz}$, however, the last third (when
$U_{E}$ dissipates; see Fig. 2c) has much higher $q_{e}$ than the first two
thirds. Therefore, we conclude that energization by the motional electric
field $E_{m,y}$ differs from reconnection energization mechanisms that operate
earlier (or without $B_{gz}$), and yields higher $q_{e}$. In 3D, for strong
$B_{gz}$, we see the last two-thirds have higher $q_{e}$ than the first third,
suggesting that $U_{E}$ builds up a little in 3D, but is promptly dissipated.
That is, energization-via-$E_{m,y}$ operates continuously throughout 3D
reconnection, but only at the end of 2D$xy$ reconnection due to the high
integrity of plasmoids. Besides being useful in their own right, measurements
of $q_{e}$ can thus help distinguish between different energization
mechanisms.

##  7 Summary

Using PIC simulation, we evolve a thin current sheet in 3D semirelativistic
plasma (see Table 1), investigating energy conversion for different guide
magnetic fields $B_{gz}$. The 3D simulations are compared with 2D$xy$
(reconnection) and 2D$yz$ (drift-kink) simulations with the same initial
conditions.

We quantify the net changes in transverse and guide magnetic, electric,
electron kinetic, and ion kinetic energies. We find:

  * •

Increasing $B_{gz}$ slows energy conversion and reduces the total released
magnetic energy.

  * •

Increasing $B_{gz}$ suppresses ion energization more than electron
energization, increasing $q_{e}\equiv\Delta U_{e}/(\Delta U_{e}+\Delta
U_{i})$. In 3D, $\Delta U_{e}/U_{Bxy0}$ falls from 0.16 to 0.09 as
$B_{gz}/B_{0}$ increases from 0.25 to 2, while $\Delta U_{i}/U_{Bxy0}$ falls
from 0.17 to 0.02; $q_{e}$ rises from roughly $0.5$ to $0.8$.

  * •

In 3D only, the flux-rope kink instability further energizes particles (mostly
electrons), releasing the reconnected transverse magnetic energy stored in
flux ropes, but not their guide-field energy (which is increased by
compression during reconnection).

  * •

Relativistic reconnection in strong guide field transfers energy to transient
motional electric fields $E_{m,y}$ in flux ropes before preferentially
energizing electrons. In 2D$xy$, $U_{E}$ grows substantially, and dissipates
only at the end of reconnection; in 3D, $U_{E}$ dissipates promptly, possibly
due to flux-rope decay.

  * •

Different plasma processes operating in thin CSs yield different electron/ion
energy partitions. We roughly distinguished electron energy fractions $q_{e}$
due to four mechanisms: DKI, classic 2D reconnection, dissipation of
$E_{m,y}$, and flux-rope kink instability. Their relative roles (especially
for DKI and reconnection) in 3D may be sensitive to the CS configuration.

We have thus characterized the conversion of magnetic energy to electron and
ion kinetic energy in 3D CSs, addressing a fundamental plasma-physics question
and obtaining new insight into the interplay of multiple dissipation channels
yielding different energization of electrons versus ions. Importantly, we
measured, from first-principles simulation, electron and ion heating fractions
critical for subgrid modeling of electron energization. These results will
help connect general-relativistic magnetohydrodynamic simulations (e.g.
Ripperda et al., 2019; Dexter et al., 2020; Scepi et al., 2022) with
observations of flaring emission from black hole accretion flows.

Acknowledgments This work was supported by NSF grants AST-1806084 and
AST-19033350 and by NASA grants 80NSSC20K0545 and 80NSSC22K0828. The 3D
simulations used computer time provided by the U.S. Department of Energy’s
(DOE) Innovative and Novel Computational Impact on Theory and Experiment
(INCITE) Program, and in particular used resources from the Argonne Leadership
Computing Facility, a U.S. DOE Office of Science user facility at Argonne
National Laboratory, which is supported by the Office of Science of the U.S.
DOE under Contract No. DE-AC02-06CH11357. The 2D simulations were run on the
Frontera supercomputer at the Texas Advanced Computing Center (TACC) at The
University of Texas at Austin.

## References

  * Ball et al. (2018) Ball D, Sironi L, & Özel F, 2018, “Electron and Proton Acceleration in Trans-relativistic Magnetic Reconnection: Dependence on Plasma Beta and Magnetization,” ApJ, 862(1), 80. DOI:10.3847/1538-4357/aac820 arXiv:1803.05556
  * Ball et al. (2019) Ball D, Sironi L, & Özel F, 2019, “The Mechanism of Electron Injection and Acceleration in Transrelativistic Reconnection,” ApJ, 884(1), 57. DOI:10.3847/1538-4357/ab3f2e arXiv:1908.05866
  * Cerutti et al. (2013) Cerutti B, Werner GR, Uzdensky DA, & Begelman MC, 2013, “Simulations of Particle Acceleration beyond the Classical Synchrotron Burnoff Limit in Magnetic Reconnection: An Explanation of the Crab Flares,” ApJ, 770, 147. DOI:10.1088/0004-637X/770/2/147 arXiv:1302.6247
  * Cerutti et al. (2014) Cerutti B, Werner GR, Uzdensky DA, & Begelman MC, 2014, “Three-dimensional Relativistic Pair Plasma Reconnection with Radiative Feedback in the Crab Nebula,” ApJ, 782, 104. DOI:10.1088/0004-637X/782/2/104 arXiv:1311.2605
  * Chael et al. (2018) Chael A, Rowan M, Narayan R, Johnson M, & Sironi L, 2018, “The role of electron heating physics in images and variability of the Galactic Centre black hole Sagittarius A*,” MNRAS, 478(4), 5209–5229. DOI:10.1093/mnras/sty1261 arXiv:1804.06416
  * Daughton (1999) Daughton W, 1999, “Two-fluid theory of the drift kink instability,” J. Geophys. Res., 104(A12), 28701–28708. DOI:10.1029/1999JA900388
  * Dexter et al. (2020) Dexter J, Jiménez-Rosales A, Ressler SM, Tchekhovskoy A, Bauböck M, de Zeeuw PT, Eisenhauer F, von Fellenberg S, Gao F, Genzel R, Gillessen S, Habibi M, Ott T, Stadler J, Straub O, & Widmann F, 2020, “A parameter survey of Sgr A* radiative models from GRMHD simulations with self-consistent electron heating,” MNRAS, 494(3), 4168–4186. DOI:10.1093/mnras/staa922 arXiv:2004.00019
  * Guo et al. (2014) Guo F, Li H, Daughton W, & Liu YH, 2014, “Formation of hard power laws in the energetic particle spectra resulting from relativistic magnetic reconnection,” Phys. Rev. Lett., 113, 155005. DOI:10.1103/PhysRevLett.113.155005
  * Guo et al. (2021) Guo F, Li X, Daughton W, Li H, Kilian P, Liu YH, Zhang Q, & Zhang H, 2021, “Magnetic Energy Release, Plasma Dynamics, and Particle Acceleration in Relativistic Turbulent Magnetic Reconnection,” ApJ, 919(2), 111. DOI:10.3847/1538-4357/ac0918 arXiv:2008.02743
  * Guo et al. (2016) Guo F, Li X, Li H, Daughton W, Zhang B, Lloyd-Ronning N, Liu YH, Zhang H, & Deng W, 2016, “Efficient production of high-energy nonthermal particles during reconnection in a magnetically dominated ion-electron plasma,” ApJ Lett., 818(1), L9
  * Guo et al. (2015) Guo F, Liu YH, Daughton W, & Li H, 2015, “Particle acceleration and plasma dynamics during magnetic reconnection in the magnetically dominated regime,” ApJ, 806(2), 167
  * Hankla et al. (2022) Hankla AM, Scepi N, & Dexter J, 2022, “Non-thermal emission from the plunging region: a model for the high-energy tail of black hole X-ray binary soft states,” MNRAS, 515(1), 775–784. DOI:10.1093/mnras/stac1785 arXiv:2206.12018
  * Hesse & Birn (2000) Hesse M, & Birn J, 2000, “Near- and Mid-tail Current Flow During Substorms: Small- and Large-Scale Aspects of Current Disruption,” Geophysical Monograph Series, 118, 295. DOI:10.1029/GM118p0295
  * Kagan et al. (2013) Kagan D, Milosavljević M, & Spitkovsky A, 2013, “A Flux Rope Network and Particle Acceleration in Three-dimensional Relativistic Magnetic Reconnection,” ApJ, 774, 41. DOI:10.1088/0004-637X/774/1/41 arXiv:1208.0849
  * Liu et al. (2011) Liu W, Li H, Yin L, Albright BJ, Bowers KJ, & Liang EP, 2011, “Particle energization in 3D magnetic reconnection of relativistic pair plasmas,” Physics of Plasmas, 18(5), 052105. DOI:10.1063/1.3589304 arXiv:1005.2435
  * Liu et al. (2014) Liu YH, Daughton W, Karimabadi H, Li H, & Peter Gary S, 2014, “Do dispersive waves play a role in collisionless magnetic reconnection?,” Phys. Plasmas, 21(2), 022113. DOI:10.1063/1.4865579
  * Markidis et al. (2014) Markidis S, Lapenta G, Delzanno GL, Henri P, Goldman MV, Newman DL, Intrator T, & Laure E, 2014, “Signatures of secondary collisionless magnetic reconnection driven by kink instability of a flux rope,” Plasma Phys. Contr. F., 56(6), 064010\. DOI:10.1088/0741-3335/56/6/064010 arXiv:1408.1144
  * Melzani et al. (2014a) Melzani M, Walder R, Folini D, Winisdoerffer C, & Favre JM, 2014a, “The energetics of relativistic magnetic reconnection: ion-electron repartition and particle distribution hardness,” Astron. Astrophys., 570, A112
  * Melzani et al. (2014b) Melzani M, Walder R, Folini D, Winisdoerffer C, & Favre JM, 2014b, “Relativistic magnetic reconnection in collisionless ion-electron plasmas explored with particle-in-cell simulations,” Astron. Astrophys., 570, A111
  * Ozaki et al. (1996) Ozaki M, Sato T, Horiuchi R, & Complexity Simulation Group, 1996, “Electromagnetic instability and anomalous resistivity in a magnetic neutral sheet,” Phys. Plasmas, 3(6), 2265–2274. DOI:10.1063/1.871908
  * Pritchett & Coroniti (1996) Pritchett PL, & Coroniti FV, 1996, “The Role of the Drift Kink Mode in Destabilizing Thin Current Sheets,” J. Geomagn. Geoelectr., 48(5), 833–844. DOI:10.5636/jgg.48.833
  * Pritchett et al. (1996) Pritchett PL, Coroniti FV, & Decyk VK, 1996, “Three-dimensional stability of thin quasi-neutral current sheets,” J. Geophys. Res., 101(A12), 27413–27430. DOI:10.1029/96JA02665
  * Ressler et al. (2023) Ressler SM, White CJ, & Quataert E, 2023, “Wind-fed GRMHD simulations of Sagittarius A*: tilt and alignment of jets and accretion discs, electron thermodynamics, and multiscale modelling of the rotation measure,” MNRAS, 521(3), 4277–4298. DOI:10.1093/mnras/stad837
  * Ressler et al. (2020) Ressler SM, White CJ, Quataert E, & Stone JM, 2020, “Ab Initio Horizon-scale Simulations of Magnetically Arrested Accretion in Sagittarius A* Fed by Stellar Winds,” ApJ Lett., 896(1), L6. DOI:10.3847/2041-8213/ab9532 arXiv:2006.00005
  * Ripperda et al. (2019) Ripperda B, Porth O, Sironi L, & Keppens R, 2019, “Relativistic resistive magnetohydrodynamic reconnection and plasmoid formation in merging flux tubes,” MNRAS, 485(1), 299–314. DOI:10.1093/mnras/stz387 arXiv:1810.10116
  * Rowan et al. (2017) Rowan ME, Sironi L, & Narayan R, 2017, “Electron and Proton Heating in Transrelativistic Magnetic Reconnection,” ApJ, 850(1), 29. DOI:10.3847/1538-4357/aa9380 arXiv:1708.04627
  * Rowan et al. (2019) Rowan ME, Sironi L, & Narayan R, 2019, “Electron and Proton Heating in Transrelativistic Guide Field Reconnection,” ApJ, 873(1), 2. DOI:10.3847/1538-4357/ab03d7 arXiv:1901.05438
  * Scepi et al. (2022) Scepi N, Dexter J, & Begelman MC, 2022, “Sgr A* X-ray flares from non-thermal particle acceleration in a magnetically arrested disc,” MNRAS, 511(3), 3536–3547. DOI:10.1093/mnras/stac337 arXiv:2107.08056
  * Schoeffler et al. (2023) Schoeffler KM, Grismayer T, Uzdensky D, & Silva LO, 2023, “High-energy synchrotron flares powered by strongly radiative relativistic magnetic reconnection: 2D and 3D PIC simulations,” MNRAS, 523(3), 3812–3839. DOI:10.1093/mnras/stad1588 arXiv:2303.16643
  * Scholer et al. (2003) Scholer M, Sidorenko I, Jaroschek CH, Treumann RA, & Zeiler A, 2003, “Onset of collisionless magnetic reconnection in thin current sheets: Three-dimensional particle simulations,” Phys. Plasmas, 10(9), 3521–3527. DOI:10.1063/1.1597494
  * Sironi et al. (2016) Sironi L, Giannios D, & Petropoulou M, 2016, “Plasmoids in relativistic reconnection, from birth to adulthood: first they grow, then they go,” MNRAS, 462, 48–74. DOI:10.1093/mnras/stw1620 arXiv:1605.02071
  * Sironi & Spitkovsky (2014) Sironi L, & Spitkovsky A, 2014, “Relativistic Reconnection: An Efficient Source of Non-thermal Particles,” ApJ Lett., 783, L21. DOI:10.1088/2041-8205/783/1/L21 arXiv:1401.5471
  * Werner & Uzdensky (2017) Werner GR, & Uzdensky DA, 2017, “Nonthermal Particle Acceleration in 3D Relativistic Magnetic Reconnection in Pair Plasma,” ApJ Lett., 843, L27. DOI:10.3847/2041-8213/aa7892 arXiv:1705.05507
  * Werner & Uzdensky (2021) Werner GR, & Uzdensky DA, 2021, “Reconnection and particle acceleration in three-dimensional current sheet evolution in moderately magnetized astrophysical pair plasma,” J. Plasma Phys., 87(6), 905870613. DOI:10.1017/S0022377821001185 arXiv:2106.02790
  * Werner et al. (2018) Werner GR, Uzdensky DA, Begelman MC, Cerutti B, & Nalewajko K, 2018, “Non-thermal particle acceleration in collisionless relativistic electron-proton reconnection,” MNRAS, 473, 4840–4861. DOI:10.1093/mnras/stx2530 arXiv:1612.04493
  * Werner et al. (2016) Werner GR, Uzdensky DA, Cerutti B, Nalewajko K, & Begelman MC, 2016, “The Extent of Power-law Energy Spectra in Collisionless Relativistic Magnetic Reconnection in Pair Plasmas,” ApJ Lett., 816, L8. DOI:10.3847/2041-8205/816/1/L8 arXiv:1409.8262
  * Yamada et al. (2014) Yamada M, Yoo J, Jara-Almonte J, Ji H, Kulsrud RM, & Myers CE, 2014, “Conversion of magnetic energy in the magnetic reconnection layer of a laboratory plasma,” Nat. Commun., 5, 4774. DOI:10.1038/ncomms5774
  * Yin et al. (2008) Yin L, Daughton W, Karimabadi H, Albright BJ, Bowers KJ, & Margulies J, 2008, “Three-Dimensional Dynamics of Collisionless Magnetic Reconnection in Large-Scale Pair Plasmas,” Phys. Rev. Lett., 101(12), 125001. DOI:10.1103/PhysRevLett.101.125001
  * Zenitani & Hoshino (2005) Zenitani S, & Hoshino M, 2005, “Three-Dimensional Evolution of a Relativistic Current Sheet: Triggering of Magnetic Reconnection by the Guide Field,” Phys. Rev. Lett., 95(9), 095001. DOI:10.1103/PhysRevLett.95.095001 arXiv:astro-ph/0505493
  * Zenitani & Hoshino (2007) Zenitani S, & Hoshino M, 2007, “Particle Acceleration and Magnetic Dissipation in Relativistic Current Sheet of Pair Plasmas,” ApJ, 670, 702–726. DOI:10.1086/522226 arXiv:0708.1000
  * Zenitani & Hoshino (2008) Zenitani S, & Hoshino M, 2008, “The Role of the Guide Field in Relativistic Pair Plasma Reconnection,” ApJ, 677, 530–544. DOI:10.1086/528708 arXiv:0712.2016
  * Zhang et al. (2021) Zhang Q, Guo F, Daughton W, Li H, & Li X, 2021, “Efficient Nonthermal Ion and Electron Acceleration Enabled by the Flux-Rope Kink Instability in 3D Nonrelativistic Magnetic Reconnection,” Phys. Rev. Lett., 127(18), 185101. DOI:10.1103/PhysRevLett.127.185101 arXiv:2105.04521
  * Zhu & Winglee (1996) Zhu Z, & Winglee RM, 1996, “Tearing instability, flux ropes, and the kinetic current sheet kink instability in the Earth’s magnetotail: A three-dimensional perspective from particle simulations,” J. Geophys. Res., 101(A3), 4885–4898. DOI:10.1029/95JA03144
  * Zweibel & Yamada (2009) Zweibel EG, & Yamada M, 2009, “Magnetic Reconnection in Astrophysical and Laboratory Plasmas,” ARA&A, 47, 291–332. DOI:10.1146/annurev-astro-082708-101726