of 12
Effective Resistivity in Relativistic Reconnection: A Prescription Based on Fully Kinetic
Simulations
Abigail Moran
1
, Lorenzo Sironi
1
,
2
, Aviad Levis
3
,
4
, Bart Ripperda
4
,
5
,
6
,
7
, Elias R. Most
8
,
9
, and Sebastiaan Selvi
1
1
Department of Astronomy, Columbia University, New York, NY 10027, USA;
abigail.moran@columbia.edu
2
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave., New York, NY 10010, USA
3
Department of Computer Science, University of Toronto, Toronto, ON M5S 2E4, Canada
4
David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada
5
Canadian Institute for Theoretical Astrophysics, 60 St. George St., Toronto, ON M5S 3H8, Canada
6
Department of Physics, University of Toronto, 60 St. George St., Toronto, ON M5S 1A7, Canada
7
Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada
8
TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
9
Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA
Received 2024 November 29; revised 2024 December 18; accepted 2024 December 18; published 2025 January 10
Abstract
A variety of high-energy astrophysical phenomena are powered by the release
via magnetic reconnection
of the
energy stored in oppositely directed
fi
elds. Single-
fl
uid resistive magnetohydrodynamic
(
MHD
)
simulations with
uniform resistivity yield dissipation rates that are much lower
(
by nearly 1 order of magnitude
)
than equivalent
kinetic calculations. Reconnection-driven phenomena could be accordingly modeled in resistive MHD employing
a nonuniform,
effective
resistivity informed by kinetic calculations. In this work, we analyze a suite of fully
kinetic particle-in-cell
(
PIC
)
simulations of relativistic pair-plasma reconnection
where the magnetic energy is
greater than the rest mass energy
for different strengths of the guide
fi
eld orthogonal to the alternating
component. We extract an empirical prescription for the effective resistivity,
/
∣∣ (∣∣
(
)
)
h
a
=+
++
JJ
Benc
pp
t
p
eff
0
11
,
where
B
0
is the reconnecting magnetic
fi
eld strength,
J
is the current density,
n
t
is the lab-frame total number
density,
e
is the elementary charge, and
c
is the speed of light. The guide
fi
eld dependence is encoded in
α
and
p
,
which we
fi
t to PIC data. This resistivity formulation
which relies only on single-
fl
uid MHD quantities
successfully reproduces the spatial structure and strength of nonideal electric
fi
elds and thus provides a promising
strategy for enhancing the reconnection rate in resistive MHD simulations.
Uni
fi
ed Astronomy Thesaurus concepts:
High energy astrophysics
(
739
)
;
Plasma astrophysics
(
1261
)
;
Magnetic
fi
elds
(
994
)
;
Magnetohydrodynamics
(
1964
)
1. Introduction
Strong magnetic
fi
elds in astrophysical compact sources
provide a reservoir of magnetic energy. This energy can be
released to the plasma
resulting in particle acceleration and
nonthermal emission
when anti-aligned
fi
eld lines annihilate
in a process called magnetic reconnection. In a number of
astrophysical sources, reconnection occurs in the relativistic
regime, where the magnetic energy exceeds the plasma rest
mass energy
(
for reviews, see M. Hoshino & Y. Lyubarsky
2012
; D. Kagan et al.
2015
; F. Guo et al.
2020
,
2024
)
.
Relativistic reconnection can power a variety of high-energy
phenomena, such as emission from black hole coronae,
magnetar
fl
ares, blazar jet
fl
ares, radio and gamma-ray emission
from pulsar magnetospheres, fast radio bursts, and
fl
ares from
supermassive black hole magnetospheres.
Magnetic reconnection refers to the breaking and reconnect-
ing of oppositely directed
fi
eld lines. This requires the
ideal
condition
+
áñ
́=
E
v
B
c
0
s
to be violated for each relevant plasma species
s
. Here,
E
and
B
are the electromagnetic
fi
elds, while
v
s
is the mean three-
velocity of species
s
.In
collisional
plasmas, the ideal condition
can be broken by resistive effects due to binary particle
collisions
encoded by the resistivity appearing in Ohm's law.
When resistive effects are not important, the magnetic
fi
eld is
frozen
into the
fl
uid, as prescribed by Alfvén's theorem
(
H. Alfvén
1943
; also known as the
fl
ux freezing theorem
)
.In
dilute astrophysical plasmas, binary collisions are rare, so the
collisional resistivity is often insuf
fi
cient to break
fl
ux freezing
on interesting time and length scales.
Reconnection occurring in the
collisionless
regime requires a
kinetic description. Since the typical separation between plasma
scales and global scales is very large, kinetic descriptions, e.g.,
employing the particle-in-cell
(
PIC
)
method, are unaffordable at
realistic scale separations. Fluid-type approaches such as
magnetohydrodynamics
(
MHD
)
, while suitable to model the
global dynamics, are by construction collisional and therefore
unable to capture collisionless effects. In fact, single-
fl
uid
resistive MHD simulations with uniform resistivity yield
reconnection rates in the plasmoid-dominated regime that are
much lower
(
by nearly 1 order of magnitude
)
than equivalent
kinetic calculations
(
J. Birn et al.
2001
; D. A. Uzdensky et al.
2010
; L. Comisso & A. Bhattacharjee
2016
;P.A.Cassaketal.
2017
)
. This discrepancy impacts the timescale of reconnection-
powered
fl
ares, e.g., in black hole magnetospheres
(
A. Bransgr-
ove et al.
2021
; A. Galishnikova et al.
2023
)
.
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
https:
//
doi.org
/
10.3847
/
2041-8213
/
ada158
© 2025. The Author
(
s
)
. Published by the American Astronomical Society.
Original content from this work may be used under the terms
of the
Creative Commons Attribution 4.0 licence
. Any further
distribution of this work must maintain attribution to the author
(
s
)
and the title
of the work, journal citation and DOI.
1
A large body of work has focused on identifying the
processes that can break the ideal condition in collisionless or
weakly collisional plasmas
here, wave
particle interactions
provide a form of effective collisionality. In pair plasmas, fast
reconnection is mediated by the off-diagonal terms of the
pressure tensor
(
N. Bessho & A. Bhattacharjee
2005
,
2007
;
M. Hesse & S. Zenitani
2007
; M. Melzani et al.
2014
;
M. Goodbred & Y.-H. Liu
2022
)
, which are also important for
electron-ion plasmas in the small, electron-scale diffu-
sion region
(
L. R. Lyons & D. C. Pridmore-Brown
1990
;
R. Horiuchi & T. Sato
1994
; H. J. Cai & L. C. Lee
1997
;
M. M. Kuznetsova et al.
1998
; J. Egedal et al.
2019
)
.
By identifying the dominant contributors to the breaking of
fl
ux freezing in collisionless plasmas, it may be possible to
write the corresponding nonideal electric
fi
eld as
η
eff
J
here,
η
eff
is some effective resistivity and
J
the electric current
density
which could be incorporated in resistive, single-
fl
uid
MHD approaches as a kinetically motivated subgrid prescrip-
tion
(
R. M. Kulsrud
1998
,
2001
; D. Biskamp & E. Schwarz
2001
; D. Uzdensky
2003
; N. Bessho & A. Bhattacharjee
2010
;
S. Zenitani et al.
2010
; B. Ripperda et al.
2019
; N. F. Loureiro
2023
)
. In general, parameterizing kinetic effects as an effective
resistivity is a nontrivial task
(
e.g., E. Hirvijoki et al.
2016
;
M. Lingam et al.
2017
)
. By means of a statistical analysis based
on PIC simulations, S. Selvi et al.
(
2023
)
identi
fi
ed the
mechanisms driving the nonideal electric
fi
eld in the general-
ized Ohm's law, for the case of relativistic pair-plasma
reconnection. The effective resistivity proposed by S. Selvi
et al.
(
2023
)
for the zero guide
fi
eld case
(
and earlier suggested
by N. Bessho & A. Bhattacharjee
2007
,
2012
)
has been shown
to successfully enhance the reconnection rate in resistive MHD
simulations
(
M. Bugli et al.
2024
)
.
As we discuss below, the form of effective resistivity
proposed by S. Selvi et al.
(
2023
)
suffers from a few
limitations, which may hamper its applicability. In this work,
rather than analyzing nonideal terms in the generalized Ohm's
law, we adopt an empirical approach. We perform a suite of
PIC simulations of relativistic pair-plasma reconnection with
varying guide
fi
eld strength, and we formulate an empirical
prescription for the effective resistivity
η
eff
, which is derived
directly from our PIC runs through a data-driven parameteriza-
tion. Our proposed model depends only on the electric current
density and the plasma number density, both of which are
readily available in resistive MHD codes. As compared to
S. Selvi et al.
(
2023
)
, the form of
η
eff
that we obtain has four
main advantages: it is written explicitly in single-
fl
uid MHD
quantities, does not depend on spatial derivatives, is coordinate
agnostic, and is valid for any guide
fi
eld. We demonstrate that
the formulation of
η
eff
we propose successfully reproduces the
spatial structure and strength of nonideal electric
fi
elds in our
PIC simulations, thus providing a promising strategy for
enhancing the reconnection rate in resistive MHD approaches.
2. PIC Simulation Setup
Our simulations are performed with the 3D PIC code
TRISTAN-MP
(
A. Spitkovsky
2005
)
. We use a 2D
x
y
domain, but we track all components of the particle velocity
and of the electromagnetic
fi
elds. Although the physics of
particle acceleration in relativistic reconnection is dramatically
different between 2D and 3D
(
H. Zhang et al.
2021
,
2023
)
, the
reconnection rate
which is de
fi
ned as the plasma in
fl
ow
velocity and is the focus of our work
is roughly the same
(
e.g., L. Sironi & A. Spitkovsky
2014
; G. R. Werner &
D. A. Uzdensky
2017
)
.
The in-plane magnetic
fi
eld is initialized in a
Harris
equilibrium
(
E. G. Harris
1962
)
,
/
ˆ
()
p
=D
Bx
By
tanh 2
in
0
,
where the direction of the in-plane
fi
eld reverses at
y
=
0 over a
thickness
Δ
=
70
c
/
ω
p
. Here,
c
/
ω
p
is the depth of the plasma
skin, and
/
w
p
=
ne m
4
p
0
2
is the plasma frequency, where
n
0
is the total number density of electron
positron pairs far from
the layer,
m
is the electron
/
positron mass, and
e
is the
elementary charge. We parameterize the
fi
eld strength in the
plane
B
0
by the magnetization
s
p
=
B
nmc
4
,
0
2
0
2
which we take to be
σ
=
50. We consider guide
fi
elds of
magnitude
B
g
/
B
0
=
0.0, 0.3, 0.6, and 1.0.
The upstream region is initialized with
n
0
=
64 particles per
cell
(
including both species
)
. We resolve the plasma skin depth
c
/
ω
p
with
fi
ve cells and evolve the simulation up to
w
-
4
500
p
1
.
In Appendix
A
, we choose
c
/
ω
p
=
20 cells and demonstrate
that our results are robust to spatial resolution. In Appendix
A
,
we also display cases that include strong synchrotron cooling
losses. For our
fi
ducial runs, the length of the domain in the
x
-
direction of plasma out
fl
ows is
L
x
=
1920
c
/
ω
p
. We use open
boundaries for
fi
elds and particles along the
x
-direction. The
box grows in the
y
-direction as the simulation progresses,
allowing for more plasma and magnetic
fl
ux to enter the
domain. At the end of the simulations, the length of our box
along the
y
-axis is comparable to
L
x
.
3. Resistivity Formulation
We begin by considering Ohm's law for resistive relativistic
single-
fl
uid MHD
(
S. S. Komissarov
2007
)
:
(·)
(
)
()
hr
G+ ́- = -
E
v
BEvvJv
cc
1
,1
e
2
where
Γ
is the bulk
fl
uid Lorentz factor,
v
is the
fl
uid three-
velocity,
ρ
e
is the electric charge density, and
η
is the
collisional resistivity. In a collisionless plasma, the replacement
of
η
by
η
eff
in Equation
(
1
)
can be regarded as the de
fi
nition of
an effective resistivity that incorporates kinetic effects in
single-
fl
uid resistive MHD. As we justify in Appendix
B
,we
can further assume that
|
ρ
e
v
|
=
|
J
|
and
Γ
;
1 and that the third
term in the square bracket is negligible, which yields
()
h
º+ ́=
*
EE
v
BJ
c
,2
eff
where
E
*
is the nonideal electric
fi
eld. The spatial structure of
the
z
-component of the nonideal electric
fi
eld,
*
E
z
, is shown in
the bottom row of Figure
1
, at a representative time after the
simulation has achieved a quasi-steady state
(
i.e., the
reconnection rate attains a quasi-steady value
)
. The
fi
gure
emphasizes that nonideal regions are generally larger for
increasing guide
fi
eld. We also present the spatial structure of
the total particle density
n
t
(
top row; in units of
n
0
)
and of the
magnetic energy density
(
middle row; in units of
/
p
B
8
0
2
)
, for
both
B
g
/
B
0
=
0
(
left column
)
and
B
g
/
B
0
=
1
(
right column
)
.
In order to determine the effective resistivity
η
eff
, we focus
on the
z
-component of Equation
(
2
)
, which dominates the
nonideal
fi
eld for the whole range of
B
g
/
B
0
we explore. The
2
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.
z
-component
*
E
z
is the only signi
fi
cant component for zero
guide
fi
eld, being 1 to 2 orders of magnitude larger than other
components; for nonzero guide
fi
elds, we still determine
η
eff
from
h
=
*
E
J
z
z
eff
, but we show that the same effective
resistivity properly describes other components, speci
fi
cally
h
=
*
E
J
x
x
eff
(
see Section
4
)
. Our new prescription for the
effective resistivity is derived using a data-driven phenomen-
ological model with two free parameters, which are bench-
marked with PIC simulations. We compare the performance of
our prescription to the resistivity model from S. Selvi et al.
(
2023
)
based on a kinetic approach
and to its extension
employing MHD quantities.
3.1. Kinetically Motivated Resistivity
S. Selvi et al.
(
2023
)
analyzed PIC simulations of relativistic
reconnection in pair plasmas and identi
fi
ed the terms that
dominate the nonideal electric
fi
eld in the generalized Ohm's
law
(
M. Hesse & S. Zenitani
2007
)
. Their analysis was
restricted to regions of electric dominance, de
fi
ned as having
>+
E
BB
zxy
222
(
which is nearly identical to the condition
|
E
|
>
|
B
|
in the case of zero guide
fi
eld
)
. They found that the
z
-component of the nonideal electric
fi
eld could be written as
()
h
==
áñ
áñ
¶á ñ
*
EJ
m
ne
u
v
vJ
,3
z
S
z
t
ez
ez
yey z
23,kin
2
where
n
t
is the total number density
(
including both electrons
and positrons
)
;
v
ez
and
u
ez
are, respectively, the mean
electron three- and four-velocity in the
z
-direction;
10
and
v
ey
;
v
y
is the mean three-velocity along
y
, which is roughly
the same for both species
(
hereafter, we call
v
y
the single-
fl
uid
y
velocity
)
.
The effective resistivity proposed by S. Selvi et al.
(
2023
)
in
Equation
(
3
)
has a few limitations:
(
i
)
it provides a satisfactory
description of the nonideal electric
fi
eld only for
B
g
=
0; and
(
ii
)
it was derived considering regions of electric dominance,
which are only a subset of the regions hosting nonideal
fi
elds
(
L. Sironi
2022
; S. R. Totorica et al.
2023
)
, where resistive
effects are important. In order to derive Equation
(
3
)
, S. Selvi
et al.
(
2023
)
used the approximation
()
()
¶áñáñ»áñ¶áñ
nv u
nu
v
,4
y
e
ey
ez
e
ez
y
ey
which is valid only in the vicinity of the center of the current
sheet. In fact, as shown in Figure
2
, the effective resistivity in
Equation
(
3
)(
hereafter,
η
S
23,kin
)
provides a reasonable descrip-
tion of the nonideal
fi
eld near the center of the layer
(
|
y
|
ω
p
/
c
1
)
, where
|
E
|
>
|
B
|
(
blue shaded area
)
, but it
signi
fi
cantly overestimates the ground truth
(
i.e., the direct
measurement of
*
E
z
from PIC runs
)
farther away from the layer
(
|
y
|
ω
p
/
c
1
)
.
For use in single-
fl
uid MHD codes, Equation
(
3
)
needs to be
rewritten using
fl
uid quantities. As we have already discussed
above, the mean three-velocity along
y
is roughly the same for
the two species,
v
ey
;
v
y
. The most reasonable approx-
imation for the ratio between the mean four- and three-
velocities of a given species is
u
ez
/
v
ez
;
γ
, where the
mean particle Lorentz factor
(
including both bulk and internal
motions
)
can be derived from the
T
00
component of the stress
energy tensor as
γ
=
T
00
/
n
t
mc
2
. This leads to a form of
Equation
(
3
)
that can be implemented in MHD:
()
hg
==áñ¶
*
EJ
m
ne
vJ
.5
z
S
z
t
yy z
23,MHD
2
As shown in Figure
3
, Equation
(
5
)
is an excellent
approximation of the kinetic form in Equation
(
3
)(
compare
top right and bottom right panels
)
. However, as anticipated in
Figure
2
, the two forms overestimate the true resistivity
(
top
left of Figure
3
)
, especially at the boundaries of the current
layer. In an earlier version of S. Selvi et al.
(
2023
)
, Equation
(
3
)
was cast in an alternative form, approximating
/
()
()
áñ
áñ
-
u
v
Jenc
1
1
,6
ez
ez
zt
2
which only holds if each species has negligible internal
motions and moves in the
z
-direction with dimensionless drift
Figure 1.
Spatial distribution of the particle number density
n
t
(
top row; in units of
n
0
)
of the magnetic energy density
(
middle row; in units of
/
p
B
8
0
2
)
and of the
z
-component of the nonideal electric
fi
eld as de
fi
ned in Equation
(
2
)(
bottom row; in units of
B
0
)
, for simulations with
B
g
/
B
0
=
0
(
left
)
and
B
g
/
B
0
=
1
(
right
)
. The
snapshots are taken at a representative time to show the nonideal electric
fi
eld and plasmoid structure after the simulations have achieved a quasi-steady state. We
de
fi
ne
x
=
0 as the edge of the simulation domain
(
so, the plot only shows a portion of the domain
)
, while
y
=
0 in the midplane.
10
At
X
-points, positrons and electrons have opposite
v
ez
and
u
ez
, but the
ratio
u
ez
/
v
ez
is the same for both species.
3
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.
speed of
|
J
z
|/
en
t
c
. This was recently rewritten by M. Bugli
et al.
(
2024
)
in the form
()
()
h
=¶+G
*
en c
mc
e
vE
1
.7
B
t
yy
z
24
2
2
While the approximation in Equation
(
6
)
leading to
Equation
(
7
)
does not generally hold, as shown by the poor
agreement between the top right and bottom left panels in
Figure
3
, Equation
(
7
)
appears to provide a remarkably good
proxy for the ground truth
(
compare top left and bottom left
)
.
While Equation
(
7
)
appears to improve upon the kinetically
motivated model by S. Selvi et al.
(
2023
)
, it loses some of the
physical motivation of Equations
(
3
)
and
(
5
)
.
While useful, the forms of effective resistivity presented in
this subsection have some undesirable properties:
(
i
)
they only
apply to the case of zero guide
fi
eld;
(
ii
)
they contain a spatial
derivative, which makes them dif
fi
cult to include in relativistic
MHD codes while maintaining causality
(
L. Del Zanna et al.
2007
)
;
(
iii
)
they only apply to the main layer and not to the
antireconnection layers in between merging plasmoids
(
which
extend along
y
and for which the relevant velocity derivative is
x
v
x
)
; and
(
iv
)
they retain a dependence on the system
geometry
(
e.g., via the
z
-component
*
E
z
)
, which makes it hard
to incorporate in global MHD simulations where current sheets
will be curved, oscillating, and generally not aligned with the
coordinate axes. In the Section
3.2
, we turn to a more agnostic
approach that avoids some of these issues.
3.2. Prescriptive Resistivity
To overcome the limitations of the model by S. Selvi et al.
(
2023
)
, we propose an empirical approach. We expect that in
regions of strong currents
as de
fi
ned below
the nonideal
electric
fi
eld should approach
|
E
*
|
(
v
in
/
c
)
B
0
, where
v
in
is the
reconnection rate
(
i.e., the in
fl
ow velocity of plasma into the
layer; see Figure
4
)
, which implies that the effective resistivity
should be
∣∣
()
h
J
v
c
B
.8
eff
in 0
A choice of
η
eff
|
J
|
1
is inapplicable in regions of small
electric currents, where the resistivity should vanish. We
therefore design a form such that
η
eff
|
J
|
p
for small current
densities, where
p
>
0 is a free parameter. More precisely, this
should occur where
|
J
|
=
en
t
c
. Adding a normalization factor
α
, this motivates choosing a form
/
∣∣
∣∣
(
)
∣∣[
(
∣∣) ]
()
h
aa
=
+
=
+
++
+
J
JJJ
B
en c
B
en c
1
.9
p
p
t
p
t
p
eff
0
11
0
1
We will determine free parameters
α
and
p
from PIC
simulations. This scales as
/
∣∣ (
)
h
μ
+
J
en c
p
t
p
eff
1
at small
currents and approaches
η
eff
=
α
B
0
/
(
2
|
J
|
)
for
|
J
|
;
en
t
c
.We
therefore expect
α
/
2
;
v
in
/
c
, as we indeed
fi
nd below
(
see also
Figure
4
)
. The condition
|
J
|
;
en
t
c
corresponds to the charge
starvation regime, i.e., all charge carriers move at near the
speed of light. This limit is indeed realized in the inner region
of the current sheet: as the bottom panel of Figure
2
shows, the
1D pro
fi
les of
J
z
and
n
t
have the same shape, suggesting a
Figure 2.
1D slice of the domain along
y
through an
X
-point, for the simulation
with zero guide
fi
eld. The top panel shows the
z
-component of the nonideal
electric
fi
eld in units of
B
0
; the middle panel shows the resistivity; and the
bottom panel shows the electric current
J
z
(
in blue
)
, the number density
n
t
(
in
orange
)
, and the drift speed
v
dr
,
z
/
c
;
J
z
/
en
t
c
(
in green
)
. In the top and middle
panels, we present in blue the ground truth obtained directly from our
simulation, while other colors show various choices for
η
eff
, as described in the
legend. Our prescription for resistivity
(
Equation
(
10
))
is shown as
η
(
α
,
p
best
)
and
η
(
α
,
p
high
)
for the values of
p
de
fi
ned in Section
4
and corresponding
α
values. Regions where
|
E
|
>
|
B
|
are shaded in blue.
4
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.
nearly constant drift velocity
J
z
/
en
t
;
0.9
c
. In fact, if we
de
fi
ne the drift velocity
v
dr
J
/
en
t
, our prescription can be
written as
/
∣∣[
( ∣ ∣) ]
()
h
a
=
+
+
Jv
B
c
1
.10
dr
p
eff
0
1
In the inner region of the current sheet, where
|
v
dr
|
;
c
(
green
line in the bottom panel of Figure
2
)
, we obtain
η
eff
|
J
|
1
,
which matches the double-peaked shape of the ground truth
(
i.e.,
E
z
/
J
z
)
in the middle panel of Figure
2
. We emphasize that
the density dependence in
v
dr
J
/
n
t
is a key ingredient of our
resistivity model
in fact, the density in the middle of the sheet
can be signi
fi
cantly larger than in the immediate upstream; see
bottom panel of Figure
2
.
In Section
5
, we provide an equivalent, more general version
of Equation
(
9
)
suitable for implementation within resistive
MHD codes.
4. Results
To determine the optimal values of
α
and
p
in Equation
(
10
)
,
we consider the
z
-component of the nonideal
fi
eld and de
fi
ne a
loss, or data-
fi
t metric
() ∣
()∣∣∣
()
å
aha
=-
**
Lp
E
pJE
,,,11
xy
z
z
z
,
eff
2
i.e., we minimize the L2 loss
(
the mean squared error
)
between
η
eff
J
z
and the measured
*
E
z
. The loss is weighted by
∣∣
*
E
z
to
ensure that the large regions with negligible nonideal
fi
elds do
not skew our
fi
ndings. We calculate the optimal parameters
α
and
p
by minimizing this loss through a simple grid search. We
create a composite domain including several time snapshots of
the PIC simulations. For each case with varying guide
fi
eld, the
snapshots
(
roughly 15 in each case
)
are equally spaced from the
time when reconnection
fi
rst attains a quasi-steady state up to
the end of our simulations,
ω
p
t
=
4500. For each snapshot, we
consider a region extending along the whole domain in
x
and
with thickness 64
c
/
ω
p
along
y
(
suf
fi
cient to enclose the largest
plasmoids
)
, centered around the current sheet. As a representa-
tive case, the loss for
B
g
/
B
0
=
1 is shown in Figure
5
. The
manifold shows that a valley of small loss, with values of
L
near the global minimum
(
shown by the red point
)
, stretches
across a wide range of
p
.
In order to determine the optimal
p
and de
fi
ne a range of
acceptable values, we adopt the following procedure. We begin
by minimizing the loss on many small, randomly selected
regions
(
hereafter,
patches
)
of the composite domain. In each
Figure 3.
A comparison between the measured nonideal electric
fi
eld
*
E
z
(
top left
)
and its reconstruction
η
eff
J
z
based on different choices of
η
eff
:
η
S23,kin
(
Equation
(
3
))
in top right,
η
B24
(
Equation
(
7
))
in bottom left, and
η
S23,MHD
(
Equation
(
5
))
in bottom right. All panels are normalized to
B
0
. Within each panel, horizontal black lines
separate different time snapshots: the
fi
rst one is taken when the reconnection rate shown in Figure
4
fi
rst settles into a steady state, and the others follow after 450,
810, and 1080
w
-
p
1
, respectively. The horizontal axis is measured with respect to the center of the portion of domain that is displayed. The derivatives in Equations
(
3
)
,
(
5
)
, and
(
7
)
are computed as numerical derivatives on cells downsampled by a factor of 2.
Figure 4.
Reconnection rate
(
i.e., the plasma in
fl
ow velocity normalized by
c
)
over time for each guide
fi
eld case, as indicated in the legend. The reconnection
rate is measured as the mean in
fl
ow velocity in the region
y
=
[
672, 672
]
c
/
ω
p
. In the same color, we show with the shaded area the acceptable ranges of
α
/
2 of our resistivity prescription, which we de
fi
ne in Section
3.2
.
5
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.
small patch, we
fi
nd that there is a clearly preferred value of
p
(
i.e., a sharp minimum of Equation
(
11
)
, as opposed to the wide
minimum we
fi
nd on global scales
)
, which we will use to de
fi
ne
a range of acceptable values of
p
. We continue adding regions
until the results converge, meaning that repeatedly selecting the
same number of random patches produces the same outcome,
regardless of which regions are chosen. We vary the patch size
depending on the guide
fi
eld strength, such that the patch is
twice larger than the typical extent of a region with signi
fi
cant
nonideal
fi
elds
(
see bottom panels in Figure
1
)
. The number of
patches and the patch size used in this step are indicated in
Table
1
. To ensure that the loss in a given patch is informative
(
which is not the case for patches with small
*
E
z
)
, we require
()
()
()
>+
**
EE
median
median
threshold,
12
zz
patch
global
where the median is computed in a given patch
(
left-hand side
)
or over the whole composite domain
(
right-hand side
)
. The
threshold indicated on the right-hand side is the difference
between the 55th and 45th percentiles of the distribution of
*
E
z
in the whole domain. We use these percentiles instead of the
standard interquartile range to ensure a more robust analysis
that includes a greater portion of the domain. We de
fi
ne
p
best
as
the value that minimizes the loss when considering the
combined area of all patches.
We then create the distribution of values of
p
that minimize
the loss in each patch. The difference
Δ
p
low
between the 50th
Figure 5.
Loss manifold computed with Equation
(
11
)
for
B
g
/
B
0
=
1, as a
function of
α
and
p
. The red point marks the
(
α
,
p
)
combination of minimum
loss. Many combinations in the dark blue valley produce losses very close to
the global minimum. A correlation between the two parameters can also
be seen.
Table 1
Best-
fi
t
p
with Upper and Lower Limits of the Acceptable Range
(
Last
Column
)
for Each Guide Field
B
g
/
B
0
Patch Dim.
(
c
/
ω
p
)
Num. Patches
p
0.0
[
80, 16
]
2500
-
+
0.00
0.0
0
1.73
0.3
[
160, 40
]
2000
-
+
9
.59
5.37
8.59
0.6
[
240, 40
]
600
-
+
1
5.4
2.1
3.8
1.0
[
120, 24
]
1500
-
+
1
8.2
6.2
5.1
Note.
We indicate the dimensions
(
along
x
and
y
, respectively
)
of the patches
used to compute the distribution of values of
p
(
second column
)
as well as the
number of patches
(
third column
)
.
Table 2
For Each Guide Field We Show the Best-
fi
t Value of
p
(
Second Column
)
and
the Lowest and Highest Acceptable Values
(
Third Column
)
B
g
/
B
0
p
best
[
p
low
,
p
high
]
α
(
p
)
0.0
0.00
[
0.00, 1.73
]
0.0369
p
+
0.3268
0.3
9.59
[
4.22, 18.2
]
0.0046
p
+
0.1295
0.6
15.4
[
13.3, 19.2
]
0.0017
p
+
0.1002
1.0
18.2
[
12.0, 23.3
]
0.0010
p
+
0.0702
Note.
These bounds are calculated as described in the text. We also show a
function that returns the optimal
α
for a given
p
within this range
(
fourth
column
)
.
Figure 6.
Best
fi
t
(
p
best
)
and upper bound
(
p
high
)
as a function of guide
fi
eld
strength.
Table 3
The Results Obtained When Varying the Size of Patches Used to Compute
p
(
Default Values Are in Table
1
)
and the Percentiles Used in the Threshold for
Patch Selection
(
Default Values Are 55
45
)
B
g
/
B
0
Patch Dim.
(
c
/
ω
p
)
Threshold per.
p
0.0
[
60, 32
]
55
45
-
+
0.00
0.00
1.26
0.0
[
80, 16
]
60
40
-
+
0.00
0.0
0
1.33
1.0
[
80, 40
]
55
45
-
+
1
7.8
5.9
5.
4
1.0
[
120, 24
]
60
40
-
+
1
8.2
4.8
4.8
Note.
As before, for
B
g
/
B
0
=
0.0, we use 2500 patches, and for
B
g
/
B
0
=
1.0,
we use 1500 patches.
Table 4
Performance L2 Loss of Various Resistivity Models as Compared to the
Measured
*
E
z
, for Different Guide Fields
B
g
/
B
0
Model
*
E
z
()
*
E
z
2
1
0.0
η
(
α
,
p
best
)
J
z
5.052
3.517
21.35
η
(
α
,
p
high
)
J
z
5.157
3.599
17.93
0.3
η
(
α
,
p
best
)
J
z
3.121
1.085
12.97
η
(
α
,
p
high
)
J
z
3.129
1.086
13.00
0.6
η
(
α
,
p
best
)
J
z
0.410
0.067
3.913
η
(
α
,
p
high
)
J
z
0.412
0.067
3.920
1.0
η
(
α
,
p
best
)
J
z
0.108
0.008
2.047
η
(
α
,
p
high
)
J
z
0.110
0.009
2.048
Note.
We vary the weight of the loss function as indicated in the last three
columns. We exclude cells where
n
t
<
1.
6
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.
and 16th percentiles of this distribution gives the lower limit on
allowed values of
p
,
[]
=-D
ppp
max
, 0
low
best
low
, while the
difference
Δ
p
high
between the 84th and 50th percentiles gives
the upper limit
p
high
=
p
best
+
Δ
p
high
. These values are
indicated for each guide
fi
eld in Table
2
. We plot
p
best
and
p
high
as a function of guide
fi
eld in Figure
6
. Table
3
demonstrates
that our
fi
ndings do not depend much on the patch size or the
threshold percentiles employed in Equation
(
12
)
.We
fi
nd that
the optimal value of
p
is robust to patch size and varies by
2%
when the threshold percentiles are altered. Similarly, the upper
and lower limits on on
p
decrease by a modest amount when
varying these parameters.
Within the range
[
p
low
,
p
high
]
, we then consider 15 evenly
spaced values of
p
. For each value, we
fi
nd the optimal
α
using
the loss function on the combined area of all patches. This
reveals that the two parameters are correlated. We interpolate to
fi
nd
α
(
p
)
and show the resulting linear
fi
ts in Table
2
. The
range of
α
/
2 corresponding to the interval
[
p
low
,
p
high
]
is
shown by the colored shaded bands in Figure
4
. As expected,
α
/
2 matches well with the measured reconnection rate, for all
guide
fi
eld cases
(
solid lines of the same color
)
.
We
fi
nally assess how well the reconstructed
η
eff
J
z
captures
the nonideal
fi
eld
*
E
z
obtained directly from our PIC
simulations. Table
4
shows the L2 loss obtained for
p
=
p
best
or
p
=
p
high
. For each of the two choices, the corresponding value
of
α
is obtained from the linear
fi
t in Table
2
.We
fi
nd that,
regardless of the weight adopted in the loss function
(
no
weight,
∣∣
*
E
z
or
(
)
*
E
z
2
)
, the L2 loss increases by less than
~
10%
when using
p
high
, as compared to choosing
p
best
(
and for the
unweighted loss of
B
g
/
B
0
=
0,
p
high
performs better than
p
best
)
.
We therefore regard all solutions within the range of
[
p
best
,
p
high
]
as acceptable.
This is also con
fi
rmed by the 2D spatial pro
fi
les shown in
Figure
7
. For all the guide
fi
elds we explore, we present the
ground truth in the left column
(
i.e., the nonideal
fi
eld
*
E
z
measured directly from our simulations
)
, the reconstruction
η
eff
(
α
,
p
best
)
J
z
in the middle column
(
here,
α
is the value
corresponding to
p
best
based on the linear
fi
t in Table
2
)
, and
the reconstruction
η
eff
(
α
,
p
high
)
J
z
in the right column
(
here,
α
is
Figure 7.
A comparison between the measured nonideal electric
fi
eld
*
E
z
(
left column
)
and its reconstruction
η
eff
J
z
based on our prescription in Equation
(
10
)
, for the
whole range of guide
fi
elds we explored. The middle column shows
η
eff
(
α
,
p
best
)
J
z
(
here,
α
is the value corresponding to
p
best
based on the linear
fi
t in Table
2
)
, while
the right column shows
η
eff
(
α
,
p
high
)
J
z
(
α
is the value corresponding to
p
high
)
. All panels are normalized to
B
0
. Within each panel, horizontal black lines separate
different time snapshots: the
fi
rst one is taken when the reconnection rate shown in Figure
4
fi
rst settles into a steady state, and the others follow after 450, 810, and
1080
w
-
p
1
respectively. The horizontal axis is measured with respect to the center of the portion of domain that is displayed
(
which is a small fraction of the composite
domain used to determine the best-
fi
t values of
α
and
p
)
.
7
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.
the value corresponding to
p
high
)
. The plot shows that the two
reconstructions are equally good for nonzero guide
fi
elds, while
for
B
g
/
B
0
=
0, the case
η
eff
(
α
,
p
high
)
J
z
seems to capture better
the longitudinal extent of nonideal regions. Most importantly,
our prescriptive resistivity performs clearly better than the
kinetically motivated models presented in Figure
3
.
Although our prescriptive resistivity was benchmarked with
the
z
-component of the nonideal
fi
eld, it can successfully model
other nontrivial components that appear for nonzero guide
fi
elds. This is shown in Figure
8
. While
*
E
y
is negligible for all
guide
fi
elds, there are distinct areas in which
*
E
x
is signi
fi
cant
for nonzero guide
fi
eld cases. We calculate
η
eff
J
x
via
Equation
(
10
)
using the same
α
and
p
from the analysis of
the
z
-component described above. From the results in Figure
8
,
we can conclude that the scalar resistivity in Equation
(
10
)
provides a satisfactory description of all components of the
nonideal electric
fi
eld across the whole range of guide
fi
elds
that we explore.
5. Discussion
We have performed a suite of PIC simulations of relativistic
pair-plasma reconnection with varying guide
fi
eld strength, and
we have formulated an empirical prescription for the effective
resistivity
η
eff
in Equation
(
9
)
or equivalently Equation
(
10
)
.
Our prescription depends on two free parameters,
α
and
p
,
which are derived directly from our PIC runs
with
α
/
2
expected to be comparable to the dimensionless reconnection
rate. As compared to the kinetically motivated model proposed
by S. Selvi et al.
(
2023
)
, the form of
η
eff
that we propose has
four main advantages: it is explicitly written in single-
fl
uid
MHD quantities, does not depend on spatial derivatives, is
coordinate agnostic, and is valid for any guide
fi
eld. It depends
only on the electric current density and the particle number
density
(
and the two free parameters discussed above
)
.We
have demonstrated that the scalar resistivity we propose
successfully describes the spatial structure and strength of all
components of the nonideal
fi
eld. It thus provides a promising
strategy for enhancing the reconnection rate in relativistic
resistive MHD approaches.
To con
fi
rm the robustness of our
fi
ndings, we demonstrate in
Appendix
A
that the form in Equation
(
10
)(
with
α
and
p
determined from our reference runs
)
provides an excellent
description of nonideal
fi
elds in independent simulations,
which either include synchrotron cooling or resolve the plasma
skin depth with 20 cells
(
as compared to 5 cells for our
reference runs
)
.
We conclude with an important remark. Our prescription in
Equation
(
9
)
can be equivalently written as
∣∣
()
∣∣
∣∣
h
=
a
-
+
*
*
*
E
en c
.13
E
E
t
B
eff
p
0
1
1
In the limit of very high
p
, the square bracket is elevated to a
very small power, yielding a contribution of order unity.
Furthermore, Figure
5
suggests that, as long as
p
is large, our
results do not signi
fi
cantly depend on its precise value. In the
Figure 8.
Same as Figure
7
, but for the
x
-component of the nonideal
fi
eld, which appears for nonzero guide
fi
eld cases. The resistivity is based on Equation
(
10
)
and
employs the same values of
α
and
p
as in Figure
7
. Although our model was not developed using the
x
-component of
E
*
, it is able to recover the basic structure of
regions with signi
fi
cant
*
E
x
. The magnitude of
*
E
x
is in units of
B
0
, and the times shown are the same as in Figure
7
.
8
The Astrophysical Journal Letters,
978:L45
(
12pp
)
, 2025 January 10
Moran et al.