General relativistic force-free electrodynamics with a discontinuous
Galerkin-finite difference hybrid method
Yoonsoo Kim ,
1
,*
Elias R. Most ,
1,2
William Throwe ,
3
Saul A. Teukolsky ,
1,2,3
and Nils Deppe
4,3
1
Theoretical Astrophysics 350-17,
California Institute of Technology
, Pasadena, California 91125, USA
2
Walter Burke Institute for Theoretical Physics,
California Institute of Technology
,
Pasadena, California 91125, USA
3
Cornell Center for Astrophysics and Planetary Science,
Cornell University
, Ithaca, New York 14853, USA
4
Department of Physics,
Cornell University
, Ithaca, New York 14853, USA
(Received 1 April 2024; accepted 14 May 2024; published 10 June 2024)
Relativistic plasmas around compact objects can sometimes be approximated as being force-free. In this
limit, the plasma inertia is negligible and the overall dynamics is governed by global electric currents. We
present a novel numerical approach for simulating such force-free plasmas, which allows for high accuracy
in smooth regions as well as capturing dissipation in current sheets. Using a high-order accurate
discontinuous Galerkin method augmented with a conservative finite-difference method, we demonstrate
efficient global simulations of black hole and neutron star magnetospheres. In addition to a series of
challenging test problems, we show that our approach can
—
depending on the physical properties of the
system and the numerical implementation
—
be up to
10
× more efficient than conventional simulations,
with a speedup of
2
–
3
× for most problems we consider in practice.
DOI:
10.1103/PhysRevD.109.123019
I. INTRODUCTION
Compact objects such as neutron stars and black holes
can feature some of the strongest magnetic fields in the
Universe. Under these conditions, the environments sur-
rounding them can be filled with a highly conducting
plasma. The plasma dynamics of these magnetospheres
are thought to be responsible for several observable
transients in the radio
[1
–
3]
and x-ray
[4
–
8]
bands.
While the description of emission processes fundamen-
tally necessitates modeling the relevant kinetic scales
[9,10]
, the available energy budget as well as the presence
of any dissipative or emitting region inside the magneto-
sphere is a result of the bulk dynamics. It is this latter
aspect that our present work aims to advance. Since these
scenarios are highly nonlinear, their effective description
requires numerical approaches.
The global dynamics of the plasma is usually modeled
under several simplifying assumptions. In a very strongly
magnetized magnetosphere, the inertia of the plasma can
approximately be neglected
[11,12]
. In this force-free
electrodynamics (FFE) state, the evolution of the system
is governed largely by bulk currents, obtained via an
effective closure of the Maxwell equations. It is important
to point out that the main assumption
—
neglecting plasma
inertia
—
can break down, e.g., during shock formation, as
well as the absence of physically meaningful dissipation in
reconnection regions. In this regime, the closest extension
of force-free electrodynamics is magnetohydrodynamical
(MHD) models, retaining a single-component plasma
rest-mass density. MHD studies of relativistic magneto-
spheres are not commonly employed (see, e.g., Ref.
[13]
for a notable exception in neutron star magnetospheres).
Instead, most studies adopt an FFE approach.
Recent examples include applications to magnetar quakes
[14]
, nonlinear steepening of Alfv ́
en waves
[15,16]
, mag-
netar giant flares
[17
–
20]
, outbursts from gravitational
collapse of a neutron star
[21]
(see also Refs.
[22,23]
for
related studies in electrovacuum), and black hole and
neutron star magnetospheres
[24
–
28]
. Apart from isolated
compact objects, force-free electrodynamics has also been
employed in the context of jets from massive black hole
mergers
[29,30]
and potential electromagnetic precursor to
gravitational wave events involving merger of compact
objects
[31
–
37]
.
Several approaches have been adopted in the literature
for numerically solving the FFE equations. Most com-
monly, either unlimited finite-difference
[25,27,29,38,39]
or conservative finite-volume schemes
[34,40
–
46]
have
been employed. These methods are robust, can easily
capture strong gradients inside the magnetosphere, and
work well with commonly employed mesh refinement
techniques
[47]
. However, they come with a major draw-
back. Properly capturing wave solutions over long inte-
gration times, e.g., Alfv ́
en waves
[15]
, requires a large
number of grid points, especially when less accurate
*
Corresponding author: ykim7@caltech.edu
PHYSICAL REVIEW D
109,
123019 (2024)
2470-0010
=
2024
=
109(12)
=
123019(23)
123019-1
© 2024 American Physical Society
versions of finite-difference/volume schemes are being
used. This prohibitively increases computational costs,
especially for applications such as compact binary mag-
netospheres in which scale separations can span two
orders of magnitude.
On the other hand, spectral-type methods such as the
pseudospectral method offer exponential convergence for
smooth solutions, providing a maximum in accuracy over
computational cost. Several studies have made use of
spectral schemes to solve the FFE equations
[48
–
50]
.
One serious limitation of spectral methods is the appear-
ance of unphysical oscillations (Gibbs phenomenon) near a
discontinuity or a large gradient e.g. current sheets, which
are naturally present in compact object magnetospheres.
Remedying these numerical instabilities requires special
treatments such as filtering or a limiting procedure
[51]
.In
addition, globally spectral methods are not easily paral-
lelizable, making it difficult to simulate physical scenarios
with large scale separations.
To counteract this shortcoming, popular approaches in
the literature focus on spectral element methods in which
the computational domain is divided into nonoverlapping
spectral elements, communicating only with the directly
neighboring elements through the element boundaries. This
approach allows for highly parallelizable implementations
while retaining the exponential convergence property for
smooth solutions. A concrete example of this approach, a
discontinuous Galerkin (DG) method, is gaining its pop-
ularity in computational fluid dynamics and astrophysics
[e.g.
[52
–
61]
], as well as in FFE
[62,63]
.
While a DG scheme naturally permits a discontinuity at
the element boundary, without special care to suppress
unphysical oscillations, it suffers from the same fate as
globally spectral methods described above. Several strat-
egies have been proposed in the DG literature, which are
frequently referred to as
limiters
. Common types of DG
limiters are implemented as direct manipulations on spec-
tral coefficients, addition of artificial viscosity, or a flat-
tening correction of the solution with respect to its average
value within an element. We refer the reader to
[53,64]
and
references therein for the available types of DG limiters and
related discussions.
DG limiters currently are not particularly accurate or
reliable compared with corresponding finite-volume or
finite-difference techniques, especially for curved meshes
or relativistic applications (see, e.g., Ref.
[64]
). A recently
developed alternative strategy is to supplement the DG
evolution with a more robust subelement discretization,
which has been mostly chosen to be finite volume [e.g.
[53
–
55,65
–
69]
]. Motivated by the idea of the
a posteriori
finite-volume limiting approach of Ref.
[65]
, the discon-
tinuous Galerkin-finite difference (DG-FD) hybrid method
was introduced by
[64]
(see also Refs.
[58,70]
for appli-
cations to relativistic fluid dynamics simulations).
In this paper, we present a new numerical scheme and
code for general-relativistic FFE simulations based on a
discontinuous Galerkin discretization. Our motivation is
twofold. First, we explore the suitability of the DG-FD
hybrid approach to enable large-scale, parallel yet accurate
numerical simulations, especially of compact binary mag-
netospheres. Second, since FFE on physical grounds has
very localized regions of nonsmoothness such as current
sheets, these simulations serve as an ideal testbed to
calibrate and assess the usefulness of the DG-FD hybrid
approach. Our hybrid scheme also incorporates previously
developed implicit-explicit time integration schemes
[71]
(see Refs.
[34,35]
for applications to the FFE system),
which allows us to enforce a set of algebraic constraints
present in the FFE system. This joint approach achieves
high-order convergence in smooth regions while capturing
discontinuous features such as magnetic reconnection
points and current sheets.
This article is organized as follows. In Sec.
II
, we briefly
review Maxwell
’
s equations in general relativity and
introduce the formulation we adopt in this work. We also
discuss our strategy for maintaining the force-free con-
ditions in simulations. In Sec.
III
, we describe the numeri-
cal implementation of spatial discretization, time stepping,
and the discontinuous Galerkin-finite difference hybrid
solver. We present results from a set of test problems in
Sec.
IV
, and conclude with a discussion of result in Sec.
V
.
In this paper, we adopt geometrized (
c
¼
G
¼
1
)
Heaviside-Lorentz units, where electric and magnetic fields
have been rescaled by
1
=
ffiffiffiffiffiffi
4
π
p
compared to Gaussian units.
We use the abstract index notation using latin indices
(
a; b;
) for spacetime tensors, but reserve
f
i; j; k;
...
g
for
spatial tensors. We follow the sign convention of the Levi-
Civita tensor from
[72]
,
ε
abcd
¼
ffiffiffiffiffiffi
−
g
p
½
abcd
, where
g
is the
determinant of spacetime metric and
½
abcd
¼
1
with
½
0123
¼þ
1
is the flat-space antisymmetric symbol.
1
II. GENERAL RELATIVISTIC FORCE-FREE
ELECTRODYNAMICS
We begin by outlining the mathematical description used
to numerically study magnetospheric dynamics. This
includes a general relativistic formulation of electrodynam-
ics in a curved spacetime, which is then specialized to the
force-free case: general relativistic force-free electrody-
namics (GRFFE).
A. Electrodynamics
The dynamics of electric and magnetic fields is governed
by the Maxwell equations. In covariant form, they are
given by
∇
b
F
ab
¼
J
a
ð
1
Þ
1
Note that an opposite sign convention is sometimes adopted
in the literature [e.g.,
[32,73]
].
KIM, MOST, THROWE, TEUKOLSKY, and DEPPE
PHYS. REV. D
109,
123019 (2024)
123019-2
∇
b
F
ab
¼
0
ð
2
Þ
where
F
ab
and
F
ab
¼
ε
abcd
F
cd
=
2
are the electromagnetic
field tensor and its dual, and
J
a
is the electric 4-current
density.
For the standard
3
þ
1
decomposition of the spacetime
metric
ds
2
¼
−
α
2
dt
2
þ
γ
ij
ð
dx
i
þ
β
i
dt
Þð
dx
j
þ
β
j
dt
Þ
;
ð
3
Þ
where
α
is lapse,
β
i
is the shift vector, and
γ
ij
is the spatial
metric, the normal to spatial hypersurfaces is given by
n
a
¼ð
1
=
α
;
−
β
i
=
α
Þ
;n
a
¼ð
−
α
;
0
Þ
:
ð
4
Þ
In terms of the normal vector
n
a
, electromagnetic field
tensor
F
ab
and its dual
F
ab
can be decomposed as
F
ab
¼
n
a
E
b
−
n
b
E
a
−
ε
abcd
B
c
n
d
;
ð
5
Þ
F
ab
¼
−
n
a
B
b
þ
n
b
B
a
−
ε
abcd
E
c
n
d
;
ð
6
Þ
where
n
a
E
a
¼
n
a
B
a
¼
0
:
ð
7
Þ
E
a
¼ð
0
;E
i
Þ
and
B
a
¼ð
0
;B
i
Þ
are electric and magnetic
fields in the frame of an Eulerian observer. One can read off
E
a
and
B
a
from
F
ab
using the following relations
E
a
¼
F
ab
n
b
;
ð
8
Þ
B
a
¼
−
1
2
ε
abcd
n
b
F
cd
¼
−
F
ab
n
b
:
ð
9
Þ
While analytically complete, Maxwell equations cannot
be directly evolved numerically, as any violation of the
divergence constraints [Eqs.
(1)
and
(2)
with
a
¼
0
]
will break strong hyperbolicity of the system
[74,75]
.
This can be avoided by either using constrained transport
approaches
[76]
or extending the system using effective
Lagrange multipliers
[77]
. We here adopt the latter
approach. The extended (or augmented) Maxwell equa-
tions
[73,78]
are
∇
a
ð
F
ab
þ
g
ab
ψ
Þ¼
−
J
b
þ
κ
ψ
n
b
ψ
ð
10
Þ
∇
a
ð
F
ab
þ
g
ab
φ
Þ¼
κ
φ
n
b
φ
ð
11
Þ
∇
a
J
a
¼
0
ð
12
Þ
where auxiliary scalar fields
ψ
and
φ
propagate diver-
gence constraint violations of electric field and magnetic
field.
κ
ψ
and
κ
φ
are damping constants, leading to an
exponential damping of the constraints in the character-
istic timescales
κ
−
1
ψ
;
φ
.
Performing a standard
3
þ
1
decomposition of the
extended Maxwell
’
s equations
(10)
–
(12)
using the normal
vector
n
a
and the spatial projection operator
h
a
b
≡
δ
a
b
þ
n
a
n
b
, we get
ð
∂
t
−
L
β
Þ
E
i
−
ε
ijk
ð
3
Þ
D
j
ð
α
B
k
Þþ
αγ
ij
D
j
ψ
¼
−
α
J
i
þ
α
KE
i
;
ð
13a
Þ
ð
∂
t
−
L
β
Þ
B
i
þ
ε
ijk
ð
3
Þ
D
j
ð
α
E
k
Þþ
αγ
ij
D
j
φ
¼
α
KB
i
;
ð
13b
Þ
ð
∂
t
−
L
β
Þ
ψ
þ
α
D
i
E
i
¼
−
ακ
ψ
ψ
þ
α
q;
ð
13c
Þ
ð
∂
t
−
L
β
Þ
φ
þ
α
D
i
B
i
¼
−
ακ
φ
φ
;
ð
13d
Þ
ð
∂
t
−
L
β
Þ
q
þ
D
i
ð
α
J
i
Þ¼
α
qK;
ð
13e
Þ
where
D
i
¼
h
a
i
∇
a
is the spatial covariant derivative,
K
is
the trace of extrinsic curvature,
q
¼
−
n
μ
J
μ
is the electric
charge density measured by an Eulerian observer, and
J
i
¼
h
i
a
J
a
is the spatial electric current density. Here we also
defined the spatial Levi-Civita tensor associated with the
spatial metric as
ε
abc
ð
3
Þ
≡
n
d
ε
dabc
:
ð
14
Þ
The Lie derivative along the shift vector applied to a spatial
vector
E
i
is
L
β
E
i
¼
β
j
∂
j
E
i
−
E
j
∂
j
β
i
;
and same for
B
i
on the left-hand side of Eq.
(13)
, while it is
simply a directional derivative (e.g.
L
β
ð
q
Þ¼
β
i
∂
i
q
) when
applied to a scalar variable.
Evolution equations
(13)
can be cast into conservative
form
∂
t
U
þ
∂
j
F
j
¼
S
;
ð
15
Þ
with evolved variables
U
¼
ffiffiffi
γ
p
2
6
6
6
6
6
6
4
E
i
B
i
ψ
φ
q
3
7
7
7
7
7
7
5
≡
2
6
6
6
6
6
6
4
̃
E
i
̃
B
i
̃
ψ
̃
φ
̃
q
3
7
7
7
7
7
7
5
;
ð
16
Þ
fluxes
GENERAL RELATIVISTIC FORCE-FREE ELECTRODYNAMICS
...
PHYS. REV. D
109,
123019 (2024)
123019-3
F
j
¼
2
6
6
6
6
6
6
6
6
4
−
β
j
̃
E
i
þ
α
ð
γ
ij
̃
ψ
−
ε
ijk
ð
3
Þ
̃
B
k
Þ
−
β
j
̃
B
i
þ
α
ð
γ
ij
̃
φ
þ
ε
ijk
ð
3
Þ
̃
E
k
Þ
−
β
j
̃
ψ
þ
α
̃
E
j
−
β
j
̃
φ
þ
α
̃
B
j
̃
J
j
−
β
j
̃
q
3
7
7
7
7
7
7
7
7
5
;
ð
17
Þ
and source terms
S
¼
2
6
6
6
6
6
6
6
6
4
−
α
ffiffiffi
γ
p
J
i
−
̃
E
j
∂
j
β
i
þ
̃
ψ
ð
γ
ij
∂
j
α
−
αγ
jk
Γ
i
jk
Þ
−
̃
B
j
∂
j
β
i
þ
̃
φ
ð
γ
ij
∂
j
α
−
αγ
jk
Γ
i
jk
Þ
̃
E
k
∂
k
α
þ
α
̃
q
−
α
̃
ψ
ð
K
þ
κ
ψ
Þ
̃
B
k
∂
k
α
−
α
̃
φ
ð
K
þ
κ
φ
Þ
0
3
7
7
7
7
7
7
7
7
5
;
ð
18
Þ
where
Γ
i
jk
are the Christoffel symbols associated with the
spatial metric. A prescription for the electric current density
J
i
(Ohm
’
s law) needs to be supplied to close the system.
B. Force-free limit
In the magnetospheres of neutron stars and black holes,
we expect copious production of electron-positron pairs
[79]
. The resulting plasma will be highly conductive,
effectively screening electric field components parallel to
the magnetic field. In addition, the magnetization of the
plasma will be very high, allowing us to consider the limit
in which the Lorentz force density vanishes and the plasma
becomes force-free.
The force-free conditions are given as
F
ab
J
b
¼
0
;
ð
19
Þ
F
ab
F
ab
¼
0
;
ð
20
Þ
F
ab
F
ab
>
0
:
ð
21
Þ
In terms of
E
i
,
B
i
,
q
and
J
i
, these conditions are
qE
i
þ
ε
ijk
ð
3
Þ
J
j
B
k
¼
0
;
ð
22
Þ
E
a
B
a
¼
E
i
B
i
¼
0
;
ð
23
Þ
B
2
−
E
2
>
0
;
ð
24
Þ
where
E
2
¼
E
a
E
a
¼
E
i
E
i
and
B
2
¼
B
a
B
a
¼
B
i
B
i
. The
first condition
(22)
corresponds to the vanishing Lorentz
force density, and the second one
(23)
shows the screening
of electric field along magnetic field lines. The third
condition
(24)
is called magnetic dominance, and violation
of this constraint flags the breakdown of force-free electro-
dynamics; characteristic speeds associated with Alfv ́
en
modes become complex and Maxwell equations are no
longer hyperbolic
[80]
. Physically,
E
2
≈
B
2
means that the
plasma drift speed approaches the speed of light, beyond
which the FFE approximation breaks down.
The force-free conditions also give constraints on the
electric current density. Equation
(22)
gives
J
i
in the form
J
i
¼
q
ε
ijk
ð
3
Þ
E
j
B
k
B
2
þ
ð
J
l
B
l
Þ
B
2
B
i
;
ð
25
Þ
which leaves the parallel component
J
l
B
l
undetermined.
The first term on the right hand side of
(25)
, the drift
current, is perpendicular to both electric and magnetic
fields and shows that electric charge moves collectively
with the drift velocity
v
d
¼
ε
ijk
ð
3
Þ
E
j
B
k
=B
2
.
Requiring Eq.
(23)
to always be satisfied, we obtain a
closed form expression of the parallel current
J
l
B
l
as
[43,81]
J
l
B
l
¼
ε
ijk
ð
3
Þ
ð
B
i
D
j
B
k
−
E
i
D
j
E
k
Þ
−
2
E
i
B
j
K
ij
;
ð
26
Þ
which reduces to
J
l
B
l
¼
B
i
ð
∇
×
B
Þ
i
−
E
j
ð
∇
×
E
Þ
j
ð
27
Þ
in the special relativistic limit
[12]
.
The parallel current Eq.
(26)
contains the spatial deriv-
atives of
E
and
B
, the dynamical variables that we evolve.
Including these derivatives in the source terms changes the
principal part of the Maxwell PDE system, and the resulting
system of equations is not strongly hyperbolic
[80]
.
A straightforward way to keep the force-free conditions
satisfied in numerical simulations is to algebraically impose
Eqs.
(23)
and
(24)
in the time evolution
[25,27,29,46,48
–
50]
. This commonly employed approach exactly ensures
the force-free conditions, but reduces the numerical accu-
racy to first-order convergence in time.
As we aim to implement a higher-order numerical
scheme for GRFFE, we consider an alternative strategy.
We adopt the driver term approach first implemented in
[31,82]
and applied in later studies [e.g.
[34,35]
]. In this
method, a stiff relaxation term is added to the electric
current density
J
i
to continuously damp the violation of the
force-free conditions. We adopt the following electric
current density prescription
[35]
J
i
¼
q
ε
ijk
ð
3
Þ
E
j
B
k
B
2
þ
η
E
j
B
j
B
2
B
i
þ
R
ð
E
2
−
B
2
Þ
B
2
E
i
;
ð
28
Þ
where
R
ð
x
Þ
≡
max
ð
x;
0
Þ
is the rectifier function and
η
is a
relaxation parameter. The parallel current consists of the
terms in the square bracket in Eq.
(28)
, each being
proportional to the violation of the force-free conditions
(23)
and
(24)
. They are coupled to the evolution of electric
field and drive the solution to the force-free limit with the
KIM, MOST, THROWE, TEUKOLSKY, and DEPPE
PHYS. REV. D
109,
123019 (2024)
123019-4
characteristic damping timescale
η
−
1
. The limiting case
η
→
∞
corresponds to the ideal force-free limit.
A caveat to the FFE simulations with a parallel electric
current is that the energy loss from an Ohmic dissipation
J
i
E
i
is removed out from and no further tracked in
simulations; therefore, total electromagnetic energy is
not conserved.
2
While numerical dissipation will also
contribute to the energy loss, the amount of energy
dissipation in current sheets [corresponding to the rectifier
term in Eq.
(28)
] dominates, albeit likely at a different rate
compared to a full kinetic reconnection model [e.g.
[9,83]
].
III. NUMERICAL IMPLEMENTATION
In this section, we describe the details of our numerical
scheme and its implementation. We present our method of
spatial discretization in Sec.
III A
, time integration in
Sec.
III B
, and the adaptive discontinuous Galerkin-finite
difference hybrid solver in Sec.
III C
. Our numerical
scheme described here is implemented in the open source
numerical relativity code
S
p
ECTRE
[84]
.
A. Domain decomposition and spatial discretization
The computational domain typically used in astrophysics
or numerical relativity simulations is simple enough to be
decomposed into a set of nonoverlapping deformed cubes.
We divide the domain into these deformed cubes, which are
called subdomain elements (hereafter simply
elements
).
Neighboring elements share their boundaries at an element
interface between them.
Within each element, a spectral expansion can be
performed to represent a field of interest. We also need
to define a prescription for handling boundary corrections
from element interfaces. This family of numerical methods
is broadly called spectral element methods
[85]
. We choose
to adopt the nodal discontinuous Galerkin discretization
[51]
, so our approach is formally referred to as a discon-
tinuous Galerkin spectral element method (DG-SEM),
which is often simply called a discontinuous Galerkin
(DG) method.
Each element is mapped to a reference cube spanning
f
ξ
1
;
ξ
2
;
ξ
3
g
∈
½
−
1
;
1
3
in the reference coordinate system
f
ξ
i
g
. A coordinate map
x
i
ð
ξ
j
Þ
relates the reference coor-
dinates
ξ
j
to physical coordinates
x
i
. A set of collocation
points
f
ξ
1
i
;
ξ
2
j
;
ξ
3
k
g
are chosen to represent the solution
u
ð
ξ
Þ¼
X
i;j;k
u
i;j;k
φ
i;j;k
ð
ξ
Þð
29
Þ
where
u
i;j;k
¼
u
ð
ξ
1
i
;
ξ
2
j
;
ξ
3
k
Þ
is the value of the solution at the
collocation point
ð
ξ
1
i
;
ξ
2
j
;
ξ
3
k
Þ
, and
φ
i;j;k
ð
ξ
Þ
is the nodal basis
function
φ
i;j;k
ð
ξ
1
l
;
ξ
2
m
;
ξ
3
n
Þ¼
1
;
for
i
¼
l; j
¼
m; k
¼
n
0
;
otherwise
:
ð
30
Þ
We use the tensor product basis
φ
i;j;k
ð
ξ
Þ¼
l
i
ð
ξ
1
Þ
l
j
ð
ξ
2
Þ
l
k
ð
ξ
3
Þð
31
Þ
where
l
d
ð
x
Þ
is the 1D Lagrange polynomial interpolating
collocation points along the
d
th axis. We choose to use an
isotropic DG mesh with the same polynomial degree
N
for
each spatial dimension. The resulting nodal expansion of
the solution is
u
ð
ξ
Þ¼
X
N
i
¼
0
X
N
j
¼
0
X
N
k
¼
0
u
i;j;k
l
i
ð
ξ
1
Þ
l
j
ð
ξ
2
Þ
l
k
ð
ξ
3
Þ
:
ð
32
Þ
The solution
(32)
can be also represented in a modal form
u
ð
ξ
Þ¼
X
N
p
¼
0
X
N
q
¼
0
X
N
r
¼
0
c
p;q;r
L
p
ð
ξ
1
Þ
L
q
ð
ξ
2
Þ
L
r
ð
ξ
3
Þ
;
ð
33
Þ
where
L
p
ð
x
Þ
is the Legendre polynomial of degree
p
. See
also Ref.
[86]
for a detailed derivation of formulating the
DG scheme in a curved spacetime.
In this article, we denote a scheme using the
N
th degree
polynomial basis (i.e.
N
þ
1
collocation points) in each
spatial dimension as a DG-
P
N
scheme. For instance, a DG-
P
5
scheme uses
6
3
collocation points in each element and a
solution is approximated as a fifth degree polynomial in
each spatial direction. When the solution is smooth, a DG-
P
N
scheme exhibits
O
ð
L
N
þ
1
Þ
spatial convergence where
L
is the spatial size of an element.
We mainly use a DG-
P
5
scheme, although we present
results for different DG orders where necessary. We use the
Legendre-Gauss-Lobatto collocation points with the mass
lumping approximation
[87]
. For a reduced aliasing error,
an exponential filter is applied to rescale the modal
coefficients
c
p;q;r
in Eq.
(33)
:
c
p;q;r
→
c
p;q;r
Y
n
¼f
p;q;r
g
exp
−
a
n
N
2
b
ð
34
Þ
after every DG time (sub)step. We use
a
¼
36
and
b
¼
50
,
which effectively zeros only the highest mode (
i
¼
N
) and
leaves other modes intact.
3
Filtering out the highest mode
reduces expected spatial converge of a DG-
P
N
scheme
from
O
ð
L
N
þ
1
Þ
to
O
ð
L
N
Þ
.
2
In a MHD model, it is captured as the same amount of
increase in the internal (thermal) energy of the plasma.
3
We note that this is a common practice adopted in spectral
methods for curing aliasing and has marginal effects on capturing
discontinuities, since typically a Gibbs phenomenon near a
discontinuity excites not only the highest mode but multiple
high modes simultaneously.
GENERAL RELATIVISTIC FORCE-FREE ELECTRODYNAMICS
...
PHYS. REV. D
109,
123019 (2024)
123019-5
B. Time integration
Based on the spatial discretization presented in the
previous section, evolution equations can be integrated
over time using the method of lines.
The maximum admissible time step size for a DG-
P
N
scheme is
[65,88]
Δ
t
≤
L
λ
max
ð
2
N
þ
1
Þ
c
D
ð
35
Þ
where
L
is the minimum (Cartesian) edge length of an
element,
λ
max
is the maximum characteristic speed inside
the element,
c
is a stability constant specific to a time
stepper, which is usually of order unity,
4
and
D
is the
number of spatial dimensions.
However, usage of a nontrivial coordinate map
x
ð
ξ
Þ
and
a complex geometry of elements deforms the spatial
distribution of grid points, and an actual upper bound
can differ from Eq.
(35)
. As a practical strategy, we adopt
the following expression
Δ
t
¼
f
ð
Δ
x
Þ
min
λ
max
c
D
ð
36
Þ
for the DG time step size, where
ð
Δ
x
Þ
min
is the minimum
grid spacing between DG collocation points in physical
coordinates and
f
is the CFL factor.
In order to keep the force-free constraint violations as
small as possible during evolution, we aim to use a large
value of the damping coefficient
η
for the driver term in
Eq.
(28)
, possibly up to
η
Δ
t
≳
10
. This implies that the
characteristic timescale of constraint damping
η
−
1
is
smaller than the time step size, which introduces stiffness
in evolution equations and makes explicit time integration
unstable unless an unreasonably small time step is used.
To address the stiffness from rapid constraint damping,
we adopt the implicit-explicit (IMEX) time stepping
technique. In particular, we make use of the IMEX-
SSP3(4,3,3) scheme
[71]
, which is third order in time.
In this IMEX approach, we evolve all quantities explicitly
using a standard 3rd-order Runge-Kutta scheme, and treat
only the stiff part of the source terms
(18)
implicitly.
Specifically, in the evolution of electric fields this requires
us to solve the following nonlinear algebraic equation at all
substeps,
E
i
¼ð
E
i
Þ
−
αη
Δ
t
0
E
j
B
j
B
2
B
i
þ
R
ð
E
2
−
B
2
Þ
B
2
E
i
ð
37
Þ
where
ð
E
i
Þ
are provided values and
Δ
t
0
is an IMEX-
scheme-dependent corrector step size. When
E
2
<B
2
, the
solution to this equation is analytical whereas in general
cases we employ a three-dimensional Newton-Raphson
solver with a specific initial guess.
In addition to the stiff electric current, we also apply the
IMEX time integration to the hyperbolic divergence clean-
ing parts to ensure stability,
ψ
¼
ψ
−
κ
ψ
Δ
t
0
ψ
;
ð
38
Þ
φ
¼
φ
−
κ
φ
Δ
t
0
φ
;
ð
39
Þ
which are linear equations and have exact analytic
inversions.
Because of the simplicity of the implicit equations in our
evolution system, the cost overhead from using an IMEX
scheme is less than 5% of the total runtime. Being able to
use much larger time steps more than compensates for this.
C. The discontinuous Galerkin-finite difference
hybrid method
This section describes our implementation of the DG-FD
hybrid solver for GRFFE equations. We closely follow the
original implementation of
[64]
, which was designed for
GRMHD, with several improvements and adaptations.
1. Overview of the algorithm
Consider an element performing a time step on the DG
grid. After each substep of the time integrator, the candidate
solution is monitored by the
troubled cell indicator
(TCI) to
check if the solution is admissible on the DG grid. If it is
admissible, we continue with the updated solution on the
DG grid. If the candidate solution is inadmissible, the
troubled cell indicator is flagged, we undo the DG substep,
project the DG solution onto the subelement FD grid, then
repeat the substep using the FD solver. Evolution on the FD
grid proceeds in a similar way; after every time step the
solution gets monitored by the troubled cell indicator,
which determines whether the solution needs to stay on
the FD grid or it is admissible on the DG grid. If the
candidate solution looks admissible on the DG grid, the
solution is projected back to the DG grid and the evolution
proceeds using the DG solver.
An optimal number of subelement finite-difference grid
points for a DG-
P
N
scheme is
2
N
þ
1
[65]
. We follow such
a prescription, and an element with
ð
N
þ
1
Þ
D
collocation
points on the DG grid is switched to
ð
2
N
þ
1
Þ
D
FD cells
with a uniform grid spacing
Δ
ξ
i
¼
2
=
ð
2
N
þ
1
Þ
in the
reference coordinates.
At the code initialization phase, all physical quantities
are evaluated on the FD grid to avoid potential spurious
oscillations arising from a spectral representation of the
initial data. Next, each element projects evolved variables
onto the DG grid, then either switches to the DG grid or
stays on the FD grid depending on the decision made by the
troubled cell indicator.
4
For example, the classic 4th-order Runge-Kutta method has
c
≈
1
.
39
[89]
.
KIM, MOST, THROWE, TEUKOLSKY, and DEPPE
PHYS. REV. D
109,
123019 (2024)
123019-6
The projection algorithm of scalar and tensor quantities
between DG and FD grids is described in detail in
[64]
.We
use a general sixth-order accurate interpolation scheme.
Since the scheme is general and does not respect the
physical constraints,
5
repeated applications (i.e., switching
back and forth between DG and FD too frequently) can
introduce spurious errors in the solution. To suppress this
behavior, we design the troubled cell indicator to apply
tighter criteria when switching back from FD to DG grid;
see Sec.
III C 3
.
2. Finite difference solver
Evolution on the finite-difference grid is performed using
a conservative finite-difference scheme
[90,91]
. For an
element using a DG-
P
N
scheme, we divide the reference
coordinate interval
½
−
1
;
1
into
2
N
þ
1
finite-difference
cells and project a solution from the DG grid onto cell-
centered values
f
U
ˆ
{
g
. A flux-balanced law
(15)
is discre-
tized as
dU
ˆ
{
dt
þ
∂
ξ
k
∂
x
j
ˆ
F
j
ˆ
{
þ
1
=
2
−
ˆ
F
j
ˆ
{
−
1
=
2
Δ
ξ
k
¼
S
ð
U
ˆ
{
Þð
40
Þ
where we used hat indices to label FD cells and plain indices
to label spatial directions.
Computation of a numerical flux
ˆ
F
j
ˆ
{
þ
1
=
2
is dimensionally
split, and closely follows that of the
ECHO
scheme
[92]
(see
also Ref.
[93]
for an application to neutron star mergers). At
the left and right sides of the FD cell interface
x
ˆ
{
þ
1
=
2
,
evolved variables are reconstructed using their cell-cen-
tered values
f
U
ˆ
{
g
. In our implementation, densitized
electric current density
̃
J
i
is also reconstructed to compute
fluxes associated with
̃
q
(see also Ref.
[73]
).
Once face-centered values
U
L;R
ˆ
{
þ
1
=
2
are reconstructed, the
interface Riemann flux
F
ˆ
{
þ
1
=
2
is computed using the
Rusanov (local-Lax-Friedrichs) flux formula
[94]
. Since
the principal part of our equations is linear, this solver will
reduce to the exact solution (see Ref.
[77]
).
In order to achieve high-order accuracy, a high-order
derivative corrector is added to the interface Riemann flux
to obtain the final numerical flux:
ˆ
F
ˆ
{
þ
1
=
2
¼
F
ˆ
{
þ
1
=
2
−
G
ð
4
Þ
ˆ
{
þ
1
=
2
:
ð
41
Þ
The original
ECHO
scheme uses the Riemann fluxes from
cell interfaces (e.g.
F
ˆ
{
3
=
2
) for the higher-order correction
term
G
ð
4
Þ
ˆ
{
þ
1
=
2
. Since we do not employ a constrained-trans-
port algorithm requiring a consistent and fixed stencil, we
opt for simpler cell-centered fluxes (e.g.,
F
ˆ
{
1
) for a more
compact stencil and reduced amount of data communica-
tions (see Refs.
[95,96]
).
For the simulations presented in this work, we use the
WENO5-Z reconstruction with the nonlinear weight expo-
nent
q
¼
2
[97]
. The high-order finite-difference corrector
is currently implemented only on Cartesian meshes. We
therefore use it for all of our one-dimensional test prob-
lems, where we assess numerical convergence of the
scheme, and defer to future work its applications in
multi-dimensional contexts. Consistent with previous
assessments, we find it sufficient to use only a fourth-
order accurate derivative correction when combined with
WENO5-Z
[93]
.
3. Troubled cell indicator
In order to decide when to switch between DG and FD
grids, our numerical scheme requires a robust criterion to
identify regions of nonsmoothness. Such an approach
somewhat shares its idea with popular adaptive-mesh-
refinement criteria. These criteria are inherently problem
dependent, and an optimal design of the troubled cell
indicator is at the heart of the DG-FD hybrid method.
Requirements on the indicator include
(i) A relatively low computational cost.
(ii) Early and robust detection of spurious oscillations
developing on the DG grid.
(iii) Being unflagged as soon as the oscillation no longer
exists, so that evolution can be performed by a more
efficient DG solver.
Motivated by the idea of the modal shock indicator
devised by
[98]
, we adopt the oscillation detection criterion
ffiffiffiffiffiffiffiffiffiffiffi
P
i
ˆ
u
2
i
P
i
u
2
i
s
>
ð
N
þ
1
−
M
Þ
−
α
;
ð
42
Þ
where
ˆ
u
is the solution with the lowest
M
modes filtered
out i.e.
ˆ
u
ð
ξ
Þ¼
X
N
p
¼
M
X
N
q
¼
M
X
N
r
¼
M
c
p;q;r
L
p
ð
ξ
1
Þ
L
q
ð
ξ
2
Þ
L
r
ð
ξ
3
Þ
;
ð
43
Þ
and the summations
P
i
in Eq.
(42)
with the nodal values
u
i
;
ˆ
u
i
are performed over all DG grid points. The exponent
α
in the criterion
(42)
controls the sensitivity of the
indicator. Since we filter out the highest mode on the
DG grid, the troubled cell indicator needs to use
M
≥
2
.
We use
M
¼
3
for the troubled cell indicator, effectively
monitoring power from the second and third highest
modes. Empirically we find that
M
≲
bð
N
þ
1
Þ
=
2
c
pro-
vides robust detections of discontinuities without the
indicator being excessively triggered. We use
α
¼
4
.
0
following Refs.
[58,64]
.
To avoid an element switching back and forth between
DG and FD grid in an unnecessarily frequent manner, we
5
For example, the interpolation scheme between the DG and
FD grids does not strictly preserve the force-free conditions or the
divergence (Gauss) constraints.
GENERAL RELATIVISTIC FORCE-FREE ELECTRODYNAMICS
...
PHYS. REV. D
109,
123019 (2024)
123019-7
use
α
0
¼
α
þ
1
when an element is evolving on FD grid.
The tighter bound
α
0
ensures an extra smoothness of
solution when the grid is switched back to DG, preventing
it from switching again to FD within only a few time steps.
Depending on the specific type of an evolved system,
one may consider additional physical admissibility criteria
(e.g. positivity of the mass density in the case of hydro-
dynamics) for the troubled cell indicator. Since the only
physical constraints in our evolution system, the force-
free conditions, are handled by the stiff parallel
electric current, we do not impose any physics-motivated
criteria.
In our implementation of the DG-FD hybrid scheme for
GRFFE, we adopt only one criterion for the troubled cell
indicator: application of the modal sensor
(42)
to the
magnitude of
̃
B
i
. While it looks somewhat oversimplified
that the information of a single scalar quantity is used for
monitoring a system with nine evolution variables
f
̃
E
i
;
̃
B
i
;
̃
ψ
;
̃
φ
;
̃
q
g
, we show in Sec.
IV
that it is capable
of detecting troubled elements in a satisfactory manner.
D. Outer boundary condition
In 3D simulations, the outer boundary of the computa-
tional domain is usually placed far out to avoid spurious
boundary effects leaking into the internal evolution. Still,
in order to suppress potential unphysical noise or reflec-
tions at the outer boundary, we implement a no-incoming
Poynting flux boundary condition as follows. The evolved
variables at the outer boundary
U
out
¼ð
̃
E
i
;
̃
B
i
;
̃
ψ
;
̃
φ
;
̃
q
Þ
out
are prescribed as follows. First, we copy the values of
f
̃
E
i
;
̃
B
i
;
̃
q
g
from the outermost grid points. Then, if the
Poynting flux is pointing inward, we set
ð
̃
E
i
Þ
out
to zero.
Divergence cleaning scalar fields
ð
̃
ψ
Þ
out
and
ð
̃
φ
Þ
out
are always set to zero.
6
On the DG grid,
U
out
is fed as
an external state when computing the boundary correction
terms. On the FD grid, ghost zones are filled with
U
out
during the FD reconstruction step.
IV. RESULTS
In this section, we test and assess our implementation of
the DG-FD hybrid method for evolving GRFFE equations
with a suite of robust code validation problems. We perform
1D tests in Sec.
IVA
, curved spacetime tests with black
holes in Sec.
IV B
, and pulsar magnetosphere tests in
Sec.
IV C
. We also discuss accuracy and efficiency aspects
of the DG-FD hybrid method in Sec.
IV D
.
A. One-dimensional problems
One-dimensional test problems evolve initial data that
only has dependence in the
x
direction. We use a computa-
tional domain consisting of a single element along the
y
and
z
axes, and impose periodic boundary condition on those
directions. Our lowest grid resolution has 32 elements
along the
x
axis, resulting in 192 DG grid points. To
facilitate comparisons with other results available in the
literature, we note that this resolution is equivalent to 352
grid points if all elements are switched to an FD grid. The
number of elements along the
x
axis is increased by a factor
of two to run medium (64 elements) and high (128
elements) resolutions. Dirichlet boundary conditions are
applied at both ends of the
x
axis. We use the CFL factor
0.3 and parallel conductivity
η
¼
10
6
. Simulation setups are
summarized in Table
I
. Initial conditions for 1D test
problems are summarized in Appendix
A
.
1. Fast wave
Originally due to
[40]
, this test problem evolves a pure
fast mode propagating in an electrovacuum. The initial
profile advects to the
þ
x
direction with the wave speed
μ
¼
1
. The analytic solution is
Q
ð
x; t
Þ¼
Q
ð
x
−
t;
0
Þ
for
any physical quantity
Q
.
As shown in Fig.
1
, our scheme shows good convergence
in flat regions with increasing grid resolution. We observe
that the accuracy and numerical convergence of the solution
is substantially lost around two kinks present in the initial
data (corresponding to
x
¼
0
.
5
0
.
1
in Fig.
1
) at which
spatial derivatives of fields are discontinuous.
2. Alfv ́
en wave
The stationary Alfv ́
en wave problem
[24]
has a transition
layer
j
x
j
<
0
.
1
that maintains a strong parallel current, and
TABLE I. Simulation setup for 1D tests in Sec.
IVA
. Grid resolution is increased with
n
¼
0
(Low),
n
¼
1
(Med), and
n
¼
2
(High).
For the FFE breakdown problem, we use
n
¼
3
as a reference solution. Each resolution, if all elements are switched to FD, is equivalent
to
352
×
2
n
finite-difference grid points along the
x
axis.
Domain size
DG grid points
η
CFL factor
Time step size
ð
×
2
−
n
Þ
Fast wave
½
−
0
.
5
;
1
.
5
×
½
−
0
.
1
;
0
.
1
2
ð
192
×
2
n
Þ
×
6
2
10
6
0.3
9
.
22
×
10
−
4
Alfv ́
en wave
½
−
1
.
5
;
1
.
5
×
½
−
0
.
1
;
0
.
1
2
1
.
38
×
10
−
3
FFE breakdown
½
−
0
.
5
;
0
.
5
×
½
−
0
.
1
;
0
.
1
2
4
.
61
×
10
−
4
6
Normally, the level of errors associated with the divergence
cleaning part
ð
ψ
;
φ
Þ
is much smaller than that of the physical
variables
ð
E
i
;B
i
;q
Þ
. Spurious reflections, if any, in the diver-
gence cleaning parts are subdominant to Poynting fluxes trans-
mitting through the outer boundaries.
KIM, MOST, THROWE, TEUKOLSKY, and DEPPE
PHYS. REV. D
109,
123019 (2024)
123019-8
the accuracy of the test results essentially reflects how well
a numerical scheme can maintain the force-free conditions.
We show the result at
t
¼
2
.
0
in Fig.
2
. It needs to be
noted that time derivatives of fields at
t
¼
0
, from the initial
condition
(A2)
, vanish only if the parallel current
J
i
B
i
equals Eq.
(26)
. In our approach, the region
j
x
j
<
0
.
1
initially develops a small transient until the stiff relaxation
term becomes fully active within several time steps and
effectively recovers the same value of
J
l
B
l
. The amplitude
of the initial transient rapidly decreases at higher grid
resolutions. Owing to the higher-order accuracy of the
discontinuous Galerkin discretization, our result shows
good convergence and low amounts of grid dissipation.
3. FFE breakdown
The force-free electrodynamics breakdown problem,
originally designed by
[40]
, demonstrates that a state
initially satisfying the force-free conditions can later
develop into a state violating them.
B
2
−
E
2
decreases
over time toward zero in the transition layer
j
x
j
<
0
.
1
and
the magnetic dominance condition eventually breaks down.
Figure
3
shows numerical results. At
t
≳
0
.
02
, the
rectifier term restoring the
B
2
−
E
2
>
0
condition is
switched on and robustly maintains the magnetic domi-
nance at later times. Since this problem does not have a
closed form solution, we perform an additional higher
resolution run using 256 elements along the
x
axis and use
it as a reference solution to check the convergence. Similar
to the fast wave test, we note the loss of accuracy and
numerical convergence near the kinks present in the
solution.
B. Three-dimensional tests: Black hole magnetospheres
We perform a set of 3D tests in a curved spacetime
using black hole magnetosphere problems. The grid struc-
ture of the computational domain is portrayed in Fig.
4
.
A spherical shell spanning the radius
½
r
in
;r
out
is split into
six cubed-sphere wedges, which are then further refined
into elements. We use an equiangular coordinate map
along angular directions and a logarithmic map along
the radial direction.
FIG. 1. Fast wave at
t
¼
0
.
5
. Top: comparison between the
exact solution and a numerical solution with the lowest grid
resolution. Bottom: error of
E
z
for three different grid resolutions.
FIG. 2. Stationary Alfv ́
en wave at
t
¼
2
.
0
. Same plot descrip-
tion as Fig.
1
.
FIG. 3. FFE breakdown problem. Top: initial data (
t
¼
0
) and a
numerical solution with the lowest grid resolution at
t
¼
0
.
02
and
t
¼
0
.
04
. Bottom: error of
B
z
with respect to the reference
solution.
GENERAL RELATIVISTIC FORCE-FREE ELECTRODYNAMICS
...
PHYS. REV. D
109,
123019 (2024)
123019-9
1. Exact Wald solution
Wald
[99]
found a stationary electrovacuum solution of
Maxwell
’
s equations in the Kerr spacetime. The solution
for the 4-potential is given as
A
b
¼
B
0
2
½ð
∂
φ
Þ
b
þ
2
a
ð
∂
t
Þ
b
;
ð
44
Þ
where
B
0
is the field amplitude,
∂
t
and
∂
φ
are the Killing
vector fields in time and azimuthal directions, and
a
¼
J=M
2
is the dimensionless spin of the Kerr black hole.
The Wald solution with
a
¼
0
satisfies the force-free
conditions outside the horizon. Electric and magnetic fields
in Kerr-Schild coordinates are given by
̃
B
x
¼
̃
B
y
¼
0
;
̃
B
z
¼
B
0
;
̃
E
x
¼
−
2
MB
0
y
r
2
;
̃
E
y
¼
2
MB
0
x
r
2
;
̃
E
z
¼
0
:
ð
45
Þ
We evolve the initial condition
(45)
to
t
¼
5
M
and
measure the L2 error norm
L
2
ð
v
i
Þ
≡
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1
n
X
n
k
¼
1
½ð
v
x
k
Þ
2
þð
v
y
k
Þ
2
þð
v
z
k
Þ
2
s
;
ð
46
Þ
where
v
i
¼
̃
B
i
numerical
−
̃
B
i
exact
and
n
is the number of grid
points. The inner domain boundary is placed at
r
in
¼
1
.
99
M
,
at which no specific boundary condition is imposed. A
Dirichlet boundary condition is imposed at the outer
boundary
r
out
¼
20
M
. Conductivity of the magnetosphere
is turned off by setting
η
¼
0
.
In Table
II
, we show convergence studies for different
orders of DG schemes
N
¼
5
, 7, 9. Measured convergence
of DG-
P
5
and DG-
P
7
schemes is consistent with the order
of DG discretization. A somewhat slower convergence of
the DG-
P
9
scheme can be attributed to other limiting
factors such as the truncation error from time integration or
the sixth-order interpolation from the initial FD grid to DG
grid. In all test cases shown in Table
II
, all elements stayed
on the DG grid throughout the evolution.
2. Vacuum Wald problem
A time-dependent evolution of electromagnetic fields
around a Kerr black hole can be simulated with the initial
magnetic fields given by the Wald solution
(44)
where
electric fields are set to zero at
t
¼
0
. The system reaches a
FIG. 4. A half-cut illustration of the spherical grid used for
black hole tests in Sec.
IV B
. Blue lines show boundaries between
the cubed-sphere wedges, where black lines show boundaries
between each element (in this example, there are 8 elements in
each wedge). The DG-
P
5
mesh consisting of
6
3
Legendre-Gauss-
Lobatto collocation points is shown with gray lines. The total
number of elements in this example is
N
r
×
N
Ω
¼
2
×
24
.
TABLE II. Convergence tests of different DG-
P
N
schemes on the Wald solution (Sec.
IV B 1
). For each level of grid resolution, we
show the number of elements used in radial and angular directions, time step size
Δ
t=M
, L2 error norm of
̃
B
i
at
t
¼
5
M
, and the
measured order of numerical convergence.
Resolution
Elements (
N
r
×
N
Ω
)
Δ
t=M
Error(
̃
B
i
)
Convergence order
DG-
P
5
Low
1
×
6
2
.
90
×
10
−
2
8
.
71
×
10
−
4
Medium
2
×
24
1
.
15
×
10
−
2
2
.
74
×
10
−
5
4.99
High
4
×
96
5
.
76
×
10
−
3
4
.
11
×
10
−
8
4.72
DG-
P
7
Low
1
.
62
×
10
−
2
1
.
23
×
10
−
5
Medium
6
.
29
×
10
−
3
1
.
51
×
10
−
7
6.35
High
3
.
14
×
10
−
3
1
.
29
×
10
−
9
6.87
DG-
P
9
Low
1
.
03
×
10
−
2
2
.
75
×
10
−
7
Medium
3
.
95
×
10
−
3
2
.
03
×
10
−
9
7.08
High
1
.
97
×
10
−
3
8
.
04
×
10
−
12
7.98
KIM, MOST, THROWE, TEUKOLSKY, and DEPPE
PHYS. REV. D
109,
123019 (2024)
123019-10