Classical and Quantum
Gravity
PAPER • OPEN ACCESS
Simulating binary black hole mergers using
discontinuous Galerkin methods
To cite this article: Geoffrey Lovelace
et al
2025
Class. Quantum Grav.
42 035001
View the
article online
for updates and enhancements.
You may also like
Entanglement entropy in quantum black
holes
Alessio Belfiglio, Orlando Luongo, Stefano
Mancini et al.
-
Observational prospects of self-interacting
scalar superradiance with next-generation
gravitational-wave detectors
Spencer Collaviti, Ling Sun, Marios
Galanis et al.
-
AsterX: a new open-source GPU-
accelerated GRMHD code for dynamical
spacetimes
Jay V Kalinani, Liwei Ji, Lorenzo Ennoggi
et al.
-
This content was downloaded from IP address 24.24.204.12 on 14/01/2025 at 16:35
Classical and Quantum Gravity
Class. Quantum Grav.
42
(2025) 035001 (31pp)
https://doi.org/10.1088/1361-6382/ad9f19
Simulating binary black hole mergers using
discontinuous Galerkin methods
Geoffrey Lovelace
1
,
12
,
∗
, Kyle C Nelli
2
,
12
,
∗
,
Nils Deppe
3
,
4
,
5
, Nils L Vu
2
, William Throwe
5
,
Marceline S Bonilla
1
, Alexander Carpenter
1
,
Lawrence E Kidder
5
, Alexandra Macedo
1
,
Mark A Scheel
2
, Azer Afram
1
, Michael Boyle
5
,
Andrea Ceja
1
,
6
, Matthew Giesler
5
, Sarah Habib
2
,
Ken Z Jones
1
, Prayush Kumar
7
, Guillermo Lara
8
,
Denyz Melchor
1
,
9
, Iago B Mendes
2
,
10
, Keefe Mitman
5
,
Marlo Morales
1
,
11
, Jordan Moxon
2
, Eamonn O’Shea
5
,
Kyle Pannone
1
, Harald P Pfeiffer
8
,
Teresita Ramirez-Aguilar
1
,
6
, Jennifer Sanchez
1
,
6
,
Daniel Tellez
1
, Saul A Teukolsky
2
,
5
and Nikolas A Wittek
8
1
Nicholas and Lee Begovich Center for Gravitational-Wave Physics and
Astronomy, California State University Fullerton, Fullerton, CA 92834, United
States of America
2
Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena,
CA 91125, United States of America
3
Department of Physics, Cornell University, Ithaca, NY 14853, United States of
America
4
Laboratory for Elementary Particle Physics, Cornell University, Ithaca, NY 14853,
United States of America
5
Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca,
NY 14853, United States of America
6
Center for Interdisciplinary Exploration and Research in Astrophysics,
Northwestern University, Evanston, IL 60201, United States of America
7
International Centre for Theoretical Sciences, Tata Institute of Fundamental
Research, Bangalore 560089, India
8
Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
D-14467 Potsdam, Germany
12
The first two authors contributed equally.
∗
Authors to whom any correspondence should be addressed.
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.
© 2025 The Author(s). Published by IOP Publishing Ltd
1
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
9
Department of Physics and Astronomy, UCLA, Mani L. Bhaumik Institute for
Theoretical Physics, Los Angeles, CA 90095, United States of America
10
Department of Physics and Astronomy, Oberlin College, Oberlin, OH 44074,
United States of America
11
Department of Physics & Astronomy, Washington State University, Pullman, WA
99164, United States of America
E-mail:
glovelace@fullerton.edu
and
knelli@caltech.edu
Received 2 October 2024; revised 9 December 2024
Accepted for publication 13 December 2024
Published 7 January 2025
Abstract
Binary black holes are the most abundant source of gravitational-wave obser-
vations. Gravitational-wave observatories in the next decade will require tre-
mendous increases in the accuracy of numerical waveforms modeling binary
black holes, compared to today’s state of the art. One approach to achieving
the required accuracy is using spectral-type methods that scale to many pro-
cessors. Using the
SpECTRE
numerical-relativity (NR) code, we present the
first simulations of a binary black hole inspiral, merger, and ringdown using
discontinuous Galerkin (DG) methods. The efficiency of DG methods allows us
to evolve the binary through
∼
18 orbits at reasonable computational cost. We
then use
SpECTRE
’s Cauchy Characteristic Evolution (CCE) code to extract the
gravitational waves at future null infinity. The open-source nature of
SpECTRE
means this is the first time a spectral-type method for simulating binary black
hole evolutions is available to the entire NR community.
Keywords: discontinuous Galerkin
,
binary black holes
,
numerical relativity
1. Introduction
Binary black holes are the most abundant source of gravitational-wave observations to date [
1
].
Realizing the scientific potential of these observations requires accurate models of the emit-
ted gravitational waves as the black holes inspiral, merge, and ring down to a final, stationary
state. Building these models requires numerical-relativity (NR) simulations of binary black
holes, because analytic approximations (e.g. the post-Newtonian [
2
] approximation) alone
break down near the time of merger.
Since the first breakthrough simulations [
3
–
5
], the NR community has developed codes
capable of evolving two black holes through inspiral, merger, and ringdown (see [
6
,
7
] for
a review). Several groups have used NR codes to build catalogs of gravitational waveforms
for applications to gravitational-wave astronomy [
8
–
15
]. Today’s NR codes are sufficiently
accurate for the observations that LIGO and Virgo are making. However, observatories planned
for the next decade, including the Einstein Telescope [
16
] and Cosmic Explorer [
17
] on Earth
and the Laser Interferometer Space Antenna (LISA) [
18
] in space, will be so sensitive that
they will require NR waveforms with a substantial increase in accuracy [
19
–
21
].
2
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
Spectral-type methods are extremely efficient; this makes them a promising avenue
toward the ultimate goal of achieving the needed accuracy for future gravitational-wave
observatories. In comparison, almost all current NR codes for evolving binary black holes use
finite-difference methods, with numerical errors decreasing as a power law with increasing
resolution. However, recent results from the
AthenaK
code [
22
] show that finite-difference
methods using graphics processing units (GPUs) might be another approach to achieving
the required accuracy. The spectral Einstein code (
SpEC
) [
23
] uses a pseudospectral method
(see [
24
] for a review of these methods) to construct and evolve binary-black-hole initial data.
With pseudospectral methods, errors decrease exponentially with increasing number of grids
points in the computational domain’s elements (‘
p
-refinement’).
SpEC
’s exponential conver-
gence makes it highly efficient, but its performance, and therefore the achievable accuracy, is
limited by aspects of its design. For instance, because it uses computational domains divided
into a small number of high-resolution elements,
SpEC
simulations of binary black holes can-
not scale beyond
O
(
10
2
)
CPU cores.
SpEC
is also a closed-source code, unavailable to most
of the NR community. Other examples of pseudospectral or spectral methods being used for
solving the initial value problem are
Elliptica
[
25
],
FUKA
[
26
], and
bamps
[
27
]. In terms
of evolving spacetimes, the
Nmesh
[
28
] code has been used to successfully simulate single
black holes using discontinuous Galerkin (DG) methods, and the
bamps
[
29
,
30
] code uses
pseudospectral methods to evolve spacetimes with single dynamical black holes with a focus
on critical behavior [
31
–
38
] but has also simulated boson stars [
39
]. Recently, [
40
] performed
a 0.5 orbit grazing collision of two black holes (a similar setup to [
41
]) using a finite volume
grid in the strong field region and a DG method in the wave zone.
We present the first simulations of a binary black hole inspiral, merger, and ringdown using
a DG method [
42
] (see [
43
] for a review of DG). The efficiency of DG methods allows us to
evolve the binary through
∼
18 orbits at reasonable computational cost: DG, being a spectral-
type method, has exponential convergence with
p
-refinement. For context, 18 orbits is slightly
less than the median (20) for binary-black-hole simulations in the SXS catalog [
12
] (which also
uses a spectral-type method) but larger than almost all of the simulations in the RIT and Maya
catalogs [
14
,
15
] (which use finite-difference methods). We chose the length of
SpECTRE
’s
first binary-black-hole simulation largely out of convenience, balancing a desire to demon-
strate
SpECTRE
’s capability with minimizing turnaround time as we tested and fine-tuned our
methods. We expect that using
SpECTRE
to simulate more orbits would be straightforward,
without requiring changes to the code, although extending the length beyond 100 orbits would
likely require implementing in
SpECTRE
similar techniques as those discussed in [
44
], which
presents a 175-orbit
SpEC
binary-black-hole simulation.
Specifically, in this work, we present
SpECTRE
’s [
45
] first simulations of
∼
18 orbits of
inspiral, merger, and ringdown of an equal-mass, non-spinning binary black hole, using DG
methods. We then use
SpECTRE
’s Cauchy Characteristic Evolution (CCE) module [
45
–
47
] to
evolve the gravitational waves to future null infinity. These results demonstrate that DG meth-
ods can provide high-accuracy gravitational waveforms from binary black hole mergers for
application to gravitational-wave data analysis. By implementing our approach in
SpECTRE
,
an open-source NR code, we are also making a spectral-type binary-black-hole evolution code
available to the entire NR community for the first time.
The rest of this paper is organized as follows. In section
2
, we discuss the numerical methods
used in
SpECTRE
’s binary-black-hole simulations. Then, in section
3
, we first test our method’s
stability with simulations of single black holes before presenting results for simulations of
binary black holes with
SpECTRE
. We briefly conclude in section
4
.
3
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
2. Methods
2.1. Equations of motion
We adopt the standard 3
+
1 form of the spacetime metric, (see, e.g. [
48
,
49
]),
d
s
2
=
g
ab
d
x
a
d
x
b
=
−
α
2
d
t
2
+
γ
ij
(
d
x
i
+
β
i
d
t
)(
d
x
j
+
β
j
d
t
)
,
(1)
where
α
is the lapse,
β
i
the shift vector, and
γ
ij
is the spatial metric.We use the Einstein sum-
mation convention, summing over repeated indices. Latin indices from the first part of the
alphabet
a
,
b
,
c
,...
denote spacetime indices ranging from 0 to 3, while Latin indices
i
,
j
,...
are purely spatial, ranging from 1 to 3. We work in units where
c
=
G
=
1.
We evolve the first-order generalized harmonic (FOGH) system, given by [
50
],
∂
t
g
ab
=(
1
+
γ
1
)
β
k
∂
k
g
ab
−
α
Π
ab
−
γ
1
β
i
Φ
iab
+
γ
1
v
k
g
(
∂
k
g
ab
−
Φ
kab
)
,
(2)
∂
t
Φ
iab
=
β
k
∂
k
Φ
iab
−
α∂
i
Π
ab
+
αγ
2
∂
i
g
ab
+
1
2
α
n
c
n
d
Φ
icd
Π
ab
+
αγ
jk
n
c
Φ
ijc
Φ
kab
−
αγ
2
Φ
iab
,
(3)
∂
t
Π
ab
=
β
k
∂
k
Π
ab
−
αγ
ki
∂
k
Φ
iab
+
γ
1
γ
2
β
k
∂
k
g
ab
+
2
α
g
cd
(
γ
ij
Φ
ica
Φ
jdb
−
Π
ca
Π
db
−
g
ef
Γ
ace
Γ
bdf
)
−
2
α
∇
(
a
H
b
)
−
1
2
α
n
c
n
d
Π
cd
Π
ab
−
α
n
c
Π
ci
γ
ij
Φ
jab
+
αγ
0
(
2
δ
c
(
a
n
b
)
−
g
ab
n
c
)
C
c
−
γ
1
γ
2
β
i
Φ
iab
,
(4)
where
g
ab
is the spacetime metric,
Φ
iab
=
∂
i
g
ab
,
Π
ab
=
n
c
∂
c
g
ab
,
n
a
is the unit normal vector to
the spatial slice,
γ
0
damps the 1-index or gauge constraint
C
a
=
H
a
+Γ
a
,
γ
1
controls the linear
degeneracy of the system,
γ
2
damps the 3-index constraint
C
iab
=
∂
i
g
ab
−
Φ
iab
,
Γ
abc
are the
spacetime Christoffel symbols of the first kind,
Γ
a
=
g
bc
Γ
bca
, and
v
k
g
is the grid/mesh velocity
as discussed in section
2.2
.
The gauge source function
H
a
can be any arbitrary function depending only upon the space-
time coordinates
x
a
and
g
ab
, but not derivatives of
g
ab
, since that may spoil the strong hyper-
bolicity of the system [
51
,
52
].
Defining
s
i
to be the unit normal covector to a 2d surface with
s
a
=(
0
,
s
i
)
, and
s
a
=
g
ab
s
b
,
the characteristic fields for the FOGH system are [
50
]
w
g
ab
=
g
ab
,
(5)
w
0
iab
=
(
δ
k
i
−
s
k
s
i
)
Φ
kab
,
(6)
w
±
ab
=Π
ab
±
s
i
Φ
iab
−
γ
2
g
ab
,
(7)
with associated characteristic speeds
λ
w
g
=
−
(
1
+
γ
1
)
β
i
s
i
−
(
1
+
γ
1
)
v
i
g
s
i
,
(8)
λ
w
0
=
−
β
i
s
i
−
v
i
g
s
i
,
(9)
λ
w
±
=
±
α
−
β
i
s
i
−
v
i
g
s
i
,
(10)
where we denote the velocity of the grid/mesh as
v
i
g
(see section
2.2
for details on our moving
mesh method). The evolved variables as a function of the characteristic fields are given by
g
ab
=
w
g
ab
,
(11)
4
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
Π
ab
=
1
2
(
w
+
ab
+
w
−
ab
)
+
γ
2
w
g
ab
,
(12)
Φ
iab
=
1
2
(
w
+
ab
−
w
−
ab
)
s
i
+
w
0
iab
.
(13)
The constraints for the FOGH system are [
50
]
C
a
=
H
a
+Γ
a
,
(14)
C
ia
=
γ
jk
∂
j
Φ
ika
−
1
2
γ
j
a
g
cd
∂
j
Φ
icd
+
n
b
∂
i
Π
ba
−
1
2
n
a
g
cd
∂
i
Π
cd
+
∂
i
H
a
+
1
2
γ
j
a
Φ
jcd
Φ
ief
g
ce
g
df
+
1
2
γ
jk
Φ
jcd
Φ
ike
g
cd
n
e
n
a
−
γ
jk
γ
mn
Φ
jma
Φ
ikn
+
1
2
Φ
icd
Π
be
n
a
(
g
cb
g
de
+
1
2
g
be
n
c
n
d
)
−
Φ
icd
Π
ba
n
c
(
g
bd
+
1
2
n
b
n
d
)
+
1
2
γ
2
(
n
a
g
cd
−
2
δ
c
a
n
d
)
C
icd
,
(15)
C
iab
=
∂
i
g
ab
−
Φ
iab
,
(16)
C
ijab
=
∂
i
Φ
jab
−
∂
j
Φ
iab
,
(17)
and
F
a
=
1
2
γ
i
a
g
bc
∂
i
Π
bc
−
γ
ij
∂
i
Π
ja
−
γ
ij
n
b
∂
i
Φ
jba
+
1
2
n
a
g
bc
γ
ij
∂
i
Φ
jbc
+
n
a
γ
ij
∂
i
H
j
+
γ
i
a
Φ
ijb
γ
jk
Φ
kcd
g
bd
n
c
−
1
2
γ
i
a
Φ
ijb
γ
jk
Φ
kcd
g
cd
n
b
−
γ
i
a
n
b
∂
i
H
b
+
γ
ij
Φ
icd
Φ
jba
g
bc
n
d
−
1
2
n
a
γ
ij
γ
mn
Φ
imc
Φ
njd
g
cd
−
1
4
n
a
γ
ij
Φ
icd
Φ
jbe
g
cb
g
de
+
1
4
n
a
Π
cd
Π
be
g
cb
g
de
−
γ
ij
H
i
Π
ja
−
n
b
γ
ij
Π
bi
Π
ja
−
1
4
γ
i
a
Φ
icd
n
c
n
d
Π
be
g
be
+
1
2
n
a
Π
cd
Π
be
g
ce
n
d
n
b
+
γ
i
a
Φ
icd
Π
be
n
c
n
b
g
de
−
γ
ij
Φ
iba
n
b
Π
je
n
e
−
1
2
γ
ij
Φ
icd
n
c
n
d
Π
ja
−
γ
ij
H
i
Φ
jba
n
b
+
γ
i
a
Φ
icd
H
b
g
bc
n
d
+
γ
2
(
γ
id
C
ida
−
1
2
γ
i
a
g
cd
C
icd
)
+
1
2
n
a
Π
cd
g
cd
H
b
n
b
−
n
a
γ
ij
Φ
ijc
H
d
g
cd
+
1
2
n
a
γ
ij
H
i
Φ
jcd
g
cd
.
(18)
While only the gauge constraint (
14
) and 3-index constraint (
16
) are damped, all constraints
can be monitored to check the accuracy of the numerical simulation. All the constraints can
be combined into a scalar, the constraint energy, given by [
50
]
E
=
δ
ab
[
C
a
C
b
+
(
F
a
F
b
+
C
ia
C
jb
γ
ij
)]
+
δ
ab
δ
cd
(
C
iac
C
jbd
γ
ij
+
C
ikac
C
jlbd
γ
ij
γ
kl
)
.
(19)
In practice we have also found that it is typically only necessary to monitor violations of the
constraints
C
a
and
C
iab
, because they typically grow first and the other violations grow as a
consequence.
5
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
Figure 1.
A diagram of the forward and reverse mappings
x
i
(
ˆ
t
,ξ
ˆ
ı
)
and
ξ
ˆ
ı
(
t
,
x
i
)
, respect-
ively, from logical (right side of the diagram) to inertial coordinates (left side of the
diagram) for an element
Ω
k
.
2.2. DG Method
SpECTRE
uses a DG method for the spatial discretization. We refer readers to [
53
,
54
] and refer-
ences therein for a detailed discussion of the method and its implementation in
SpECTRE
; here
we summarize the necessary results. The FOGH equations are a first-order strongly hyperbolic
system in non-conservative form, which takes the general form
∂
t
u
α
+
B
i
αβ
(
u
α
)
∂
i
u
β
=
S
α
(
u
α
)
,
(20)
where
u
α
=
{
g
ab
,
Φ
iab
,
Π
ab
}
is the state vector of evolved variables,
B
i
αβ
(
u
)
is a matrix that
depends only on
u
α
, and
S
α
(
u
)
are source terms that also only depend on
u
α
. We denote
the logical coordinates of our Legendre–Gauss–Lobatto DG scheme by
{
ˆ
t
,ξ
ˆ
ı
}
=
{
ˆ
t
,ξ,η,ζ
}
and the inertial coordinates as
{
t
=
ˆ
t
,
x
i
(
ˆ
t
,ξ
ˆ
ı
)
}
. We are using a moving mesh as in [
53
,
55
];
therefore, the mapping from logical to inertial coordinates is time dependent. We denote the
determinant of the spatial Jacobian of this map as
J
=
det
(
∂
x
i
∂ξ
ˆ
ı
)
,
(21)
and the grid or mesh velocity by [
53
,
55
]
v
i
g
=
∂
ˆ
t
x
i
.
(22)
SpECTRE
decomposes the domain into a set of non-overlapping hexahedra which are
deformed using the map
x
i
(
ˆ
t
,ξ
ˆ
ı
)
illustrated in figure
1
.
SpECTRE
uses a different number of
grid points in each logical direction, which we denote by
N
̆
ı
in the
ξ
direction,
N
̆
ȷ
in the
η
direction, and
N
̆
k
in the
ζ
direction below.
Since we use a moving mesh, we evolve the system
∂
ˆ
t
u
α
+
[
B
i
αβ
(
u
α
)
−
v
i
g
δ
αβ
]
∂
i
u
β
=
S
(
u
α
)
.
(23)
We denote grid point and modal indices with a breve, i.e.
u
̆
ı
̆
ȷ
̆
k
is the value of
u
at the grid point
( ̆
ı,
̆
ȷ,
̆
k
)
. The semi-discrete equations are given by [
53
,
55
]
6
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
(
∂
ˆ
t
u
α
)
̆
ı
̆
ȷ
̆
k
=(
S
α
)
̆
ı
̆
ȷ
̆
k
−
(
B
i
αβ
−
v
i
g
δ
αβ
)
̆
ı
̆
ȷ
̆
k
(
∂ξ
ˆ
1
∂
x
i
)
̆
ı
̆
ȷ
̆
k
∑
̆
l
D
(
ˆ
1
)
̆
ı
̆
l
(
u
β
)
̆
l
̆
ȷ
̆
k
−
(
∂ξ
ˆ
2
∂
x
i
)
̆
ı
̆
ȷ
̆
k
∑
̆
l
D
(
ˆ
2
)
̆
ȷ
̆
l
(
u
β
)
̆
ı
̆
l
̆
k
+
(
∂ξ
ˆ
3
∂
x
i
)
̆
ı
̆
ȷ
̆
k
∑
̆
l
D
(
ˆ
3
)
̆
k
̆
l
(
u
β
)
̆
ı
̆
ȷ
̆
l
−
δ
N
̆
k
̆
k
w
̆
k
J
̆
ı
̆
ȷ
̆
k
J
√
∂ξ
ˆ
3
∂
x
i
γ
ij
∂ξ
ˆ
3
∂
x
j
D
α
̆
ı
̆
ȷ
N
̆
k
+
J
√
∂ξ
ˆ
3
∂
x
i
γ
ij
∂ξ
ˆ
3
∂
x
j
D
α
̆
ı
̆
ȷ
0
−
δ
N
̆
ȷ
̆
ȷ
w
̆
ȷ
J
̆
ı
̆
ȷ
̆
k
J
√
∂ξ
ˆ
3
∂
x
i
γ
ij
∂ξ
ˆ
3
∂
x
j
D
α
̆
ı
N
̆
ȷ
̆
k
+
J
√
∂ξ
ˆ
3
∂
x
i
γ
ij
∂ξ
ˆ
3
∂
x
j
D
α
̆
ı
0
̆
k
−
δ
N
̆
ı
̆
ı
w
̆
ı
J
̆
ı
̆
ȷ
̆
k
J
√
∂ξ
ˆ
3
∂
x
i
γ
ij
∂ξ
ˆ
3
∂
x
j
D
α
N
̆
ı
̆
ȷ
̆
k
+
J
√
∂ξ
ˆ
3
∂
x
i
γ
ij
∂ξ
ˆ
3
∂
x
j
D
α
0
̆
ȷ
̆
k
,
(24)
where
w
̆
ı
are the Legendre–Gauss–Lobatto integration weights. We use a method of lines
approach to integrate these in time, with the details discussed in section
2.5
below.
For the boundary terms
D
α
, we use an upwind multi-penalty method [
24
,
56
–
58
] given by
D
g
ab
=
̃
λ
ext
w
g
w
ext
,
g
ab
−
̃
λ
int
w
g
w
int
,
g
ab
,
(25)
D
Π
ab
=
1
2
(
̃
λ
ext
w
+
w
ext
,
+
ab
+
̃
λ
ext
w
−
w
ext
,
−
ab
)
+
̃
λ
ext
w
g
γ
2
w
ext
,
g
ab
−
1
2
(
̃
λ
int
w
+
w
int
,
+
ab
+
̃
λ
int
w
−
w
int
,
−
ab
)
−
̃
λ
int
w
g
γ
2
w
int
,
g
ab
,
(26)
D
Φ
iab
=
1
2
(
̃
λ
ext
w
+
w
ext
,
+
ab
−
̃
λ
ext
w
−
w
ext
,
−
ab
)
s
ext
i
+
̃
λ
ext
w
0
w
ext
,
0
iab
−
1
2
(
̃
λ
int
w
+
w
int
,
+
ab
−
̃
λ
int
w
−
w
int
,
−
ab
)
s
int
i
−
̃
λ
int
w
0
w
int
,
0
iab
,
(27)
where the spatial normal vector to the element interface
s
int
i
is pointing out of the DG element
and
̃
λ
=
0 if
λ>
0, otherwise
̃
λ
=
λ
, i.e.
̃
λ
=
λ
Θ(
−
λ
)
. Note that we assume
s
ext
i
and
s
int
i
point
in the same direction. Also note that these boundary flux terms differ from the multi-penalty
approach used in
SpEC
by a factor of 2. That is,
D
SpECTRE
=
2
D
SpEC
,
(28)
ultimately because the lifting terms are different. In
SpEC
and, similarly in
bamps
[
29
], the pen-
alty term is derived from requiring that the total energy be non-increasing, while in
SpECTRE
the terms come from an integration by parts when deriving the semi-discrete DG equations.
2.3. Boundary conditions
At the outer radial boundary, we apply constraint-preserving boundary conditions [
50
,
59
] by
adding terms to the time derivative of the characteristic fields and thus also the time derivatives
of the evolved variables. We use the characteristic fields and speeds defined in section
2.1
.
7
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
We define
d
ˆ
t
w
ˆ
α
as the time derivatives substituted into the transformation equations to the
characteristic fields. That is,
d
ˆ
t
w
g
ab
=
∂
ˆ
t
g
ab
,
(29)
d
ˆ
t
w
0
iab
=
(
δ
k
i
−
s
k
s
i
)
∂
ˆ
t
Φ
kab
,
(30)
d
ˆ
t
w
±
ab
=
∂
ˆ
t
Π
ab
±
s
i
∂
ˆ
t
Φ
iab
−
γ
2
∂
ˆ
t
g
ab
.
(31)
We also define
D
ˆ
t
w
ˆ
α
as the characteristic field transformation of the volume right-hand-side,
i.e.
∂
ˆ
t
u
α
without
any boundary terms. Finally, for brevity we define the projection tensor
P
ab
=
g
ab
+
n
a
n
b
−
s
a
s
b
, the inward directed null vector field
k
a
=(
n
a
−
s
a
)
/
√
2, and the outgoing
null vector field
l
a
=(
n
a
+
s
a
)
/
√
2.
The fields
d
ˆ
t
w
g
ab
and
d
ˆ
t
w
0
iab
are determined solely by the constraint-preserving boundary
condition, while the boundary condition for
d
ˆ
t
w
−
ab
is composed of three parts: the constraint
preserving part, the physical part, and the gauge part. We denote these as
B
C
ab
,
B
P
ab
and
B
G
ab
.
With this, the boundary conditions imposed on the fields are
d
ˆ
t
w
g
ab
=
D
ˆ
t
w
g
ab
+
λ
w
g
s
i
C
iab
,
(32)
d
ˆ
t
w
0
kab
=
D
ˆ
t
w
0
ab
+
λ
w
0
s
i
P
j
k
C
ijab
,
(33)
d
ˆ
t
w
−
ab
=
D
ˆ
t
w
−
ab
+
λ
w
−
[
B
C
ab
+
B
P
ab
+
B
G
ab
]
.
(34)
Transforming to the evolved variables we find that the following terms need to be added in
order to impose the boundary condition,
∂
ˆ
t
g
ab
→
∂
ˆ
t
g
ab
+
λ
w
g
s
i
C
iab
,
(35)
∂
ˆ
t
Π
ab
→
∂
ˆ
t
Π
ab
+
1
2
λ
w
−
[
B
C
ab
+
B
P
ab
+
B
G
ab
]
+
γ
2
λ
w
g
s
i
C
iab
,
(36)
∂
ˆ
t
Φ
iab
→
∂
ˆ
t
Φ
iab
−
s
i
2
λ
w
−
[
B
C
ab
+
B
P
ab
+
B
G
ab
]
+
λ
w
0
s
i
P
j
k
C
ijab
.
(37)
We now need to specify the
B
ab
boundary conditions. The constraint-preserving part is
B
C
ab
=
√
2
(
1
2
P
ab
l
c
+
1
2
l
a
l
b
k
c
−
l
(
a
P
b
)
c
)
(
F
c
−
s
k
C
kc
)
.
(38)
The physical boundary conditions are determined by the propagating parts of the Weyl
curvature tensor. That is,
B
P
ab
=
(
P
a
c
P
b
d
−
1
2
P
ab
P
cd
)
[
C
−
cd
−
γ
2
s
i
C
icd
]
,
(39)
where
C
−
ab
is the inward propagating part of the Weyl tensor, given by
C
±
ab
=
(
P
a
c
P
b
d
−
1
2
P
ab
P
cd
)
(
n
e
∓
s
e
)
(
n
f
∓
s
f
)
C
cedf
.
(40)
For the simulations presented here, we set
C
−
ab
=
0, though Cauchy-Characteristic Matching
(CCM) [
60
] can be used to prescribe a more physically motivated boundary condi-
tion. Recently [
61
] presented an alternative approach to CCM for providing high-order
non-reflecting boundary conditions. Finally, the gauge boundary condition is set using a
8
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
Sommerfeld condition on the components not set by the constraint-preserving and physical
boundary conditions. The projector for the gauge boundary condition is given by
δ
c
a
δ
d
b
−
P
C
ab
cd
−
P
P
ab
cd
=
δ
c
a
δ
d
b
−
1
2
P
ab
P
cd
+
2
l
(
a
P
b
)
(
c
k
d
)
−
l
a
l
b
k
c
k
d
−
P
a
c
P
b
d
+
1
2
P
ab
P
cd
=
δ
c
a
δ
d
b
+
2
l
(
a
P
b
)
(
c
k
d
)
−
l
a
l
b
k
c
k
d
−
P
a
c
P
b
d
.
(41)
The Sommerfeld condition is
B
G
ab
=
1
λ
w
−
(
2
l
(
a
P
b
)
(
c
k
d
)
−
2
k
(
a
l
b
)
k
(
c
l
d
)
−
k
a
k
b
l
c
l
d
)
(
γ
2
−
1
r
)
∂
t
g
cd
.
(42)
When evolving spacetimes with black holes, we excise the interior of the black hole as
is done in
SpEC
[
50
]. At excision boundaries, all information is flowing out of the grid and
into the black hole, so no boundary condition needs to be applied. However, we monitor the
characteristic speeds, (
8
)–(
10
), and terminate the code if any of them point into the compu-
tational domain. We denote the radius of the excision surfaces by
r
exc
. See section
2.7
for a
brief explanation of how we control
r
exc
to avoid any characteristic speed pointing into the
computational domain.
2.4. Spectral filter
We use an exponential filter applied to the spectral coefficient
c
i
in order to eliminate aliasing-
driven instabilities. Specifically, for a 1d spectral expansion
u
(
x
)=
N
∑
̆
ı
=
0
c
̆
ı
P
̆
ı
(
x
)
,
(43)
where
P
̆
ı
(
x
)
are the Legendre polynomials, we use the filter
c
̆
ı
→
c
̆
ı
exp
[
−
a
(
i
N
)
2
b
]
.
(44)
We choose the parameters
a
=
64 and
b
=
210 so that only the highest spectral mode is filtered.
Weapply the filterto all FOGH variables
g
ab
,
Φ
iab
and
Π
ab
. Notethat the filter dropsthe order of
convergence for the FOGH variables from
O
(
N
+
1
)
to
O
(
N
)
on the DG grid, but is necessary
for stability.
2.5. Time integration
We decompose the system using the method of lines and solve the resulting differen-
tial equations using a local adaptive time-stepper based on the Adams–Moulton predictor–
corrector method [
62
].The step size in each element is chosen based on an estimate of the
truncation error of the time step, using the algorithm described in [
63
] section 17.2.1. The spe-
cific values for the absolute and relative tolerances are given in section
3
. As the time-stepping
algorithm is more efficient for aligned steps of the same size, the step size in each element is
rounded down to a value of the form
0
.
1
M
/
2
n
for some non-negative integer
n
. For the highest-
resolution binary-black-hole run in section
3.3
, this results in the most-demanding element
9
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
taking 2
6
–2
8
steps for each step on the least demanding element for most of the inspiral. At
the time of merger, this can increase to as high as 2
11
steps for the most-demanding element
for each step on the least demanding element.
2.6. Gauge condition
We evolve binary black holes (section
3.3
) using the Damped Harmonic gauge condition [
52
,
64
]:
H
a
=[
μ
L
1
log
(
√
γ/α
)+
μ
L
2
log
(
1
/α
)]
n
a
−
μ
S
g
ai
β
i
/α,
(45)
using
μ
L
1
=
A
L
1
e
−
(
r
/σ
r
)
2
[
log
(
√
γ/α
)]
e
L
1
,
(46)
μ
L
2
=
A
L
2
e
−
(
r
/σ
r
)
2
[
log
(
√
γ/α
)]
e
L
2
,
(47)
μ
S
=
A
S
e
−
(
r
/σ
r
)
2
[
log
(
√
γ/α
)]
e
S
,
(48)
where
r
is the coordinate distance from the origin. This condition is designed to drive
√
γ
and
α
to one, while damping out oscillations in the shift. This is because we observe an explosive
growth in
√
γ
and a rapid collapse in
α
as the black holes merge. In practice, this ensures
coordinates remain sufficiently well behaved throughout inspiral, merger, and ringdown. The
amplitudes
A
L
1
,
A
L
2
, and
A
S
and exponents
e
L
1
,
e
L
2
and
e
S
control the amount of damping,
and the spatial decay width
σ
r
ensures that at large distances, the gauge reduces to harmonic
gauge (i.e. to
H
a
=
0). In this paper, we choose
A
L
1
=
A
S
=
1,
A
L
2
=
0,
e
L
1
=
e
L
2
=
e
S
=
2, and
σ
r
=
100
/
√
34
.
54. This choice for
σ
r
ensures that the spatial decay Gaussian falls to 10
−
15
at
a distance
r
=
100
M
from the origin.
For some of the single black-hole evolutions (section
3.1
), we instead choose
H
a
to be
Γ
a
of
the analytic initial data. For other single black-hole evolutions, we evolve in harmonic gauge,
setting
H
a
=
0 everywhere. For the gauge wave evolution (section
3.2
), we set
H
a
(
t
,
x
i
)
to the
value of
Γ
a
(
t
,
x
i
)
of the gauge wave analytic solution.
2.7. Control systems
When evolving the FOGH system, if there are black holes, the physical singularities inside of
the black holes must be excised from the computational domain. To position the excisions with
our moving mesh (described in section
2.2
), we use a feedback control system similar to what
is presented in [
55
,
65
]. As discussed in section
2.3
, the excision surfaces must have all charac-
teristic speeds pointing out of the computational domain, so that no boundary condition must
be imposed. In practice this means that the excision surfaces must remain inside the apparent
horizons, with the caveat that having them too close to the singularity causes instabilities. In
practice the excision surfaces are kept at approximately 95%–99% of the apparent horizons’
radii.
Since we
a priori
do not know the motion or shape of the apparent horizons, we use control
theory to dynamically update the parameters of the moving mesh periodically during the sim-
ulation. The time-dependent coordinate maps of the moving mesh and control signals used to
update them are discussed in [
65
] in sections 4.1–4.3, 4.5, and 5 for the inspiral and section 6
for the ringdown.
10
Class. Quantum Grav.
42
(2025) 035001
G Lovelace
et al
The details of how the control systems are implemented within the context of asynchronous
task-based parallelism along with the local adaptive time stepping described in section
2.5
are
described in [
66
].
3. Results
In this section, we begin by testing
SpECTRE
’s long-term stability and convergence; first with
evolutions of single black holes in different coordinate systems (section
3.1
) and then with
an evolution of a time-dependent gauge wave on a flat spacetime background (section
3.2
).
Finally, we present results from a complete simulation of the inspiral, merger, and ringdown
of two black holes (section
3.3
). The
SpECTRE
input files used for simulations, including gen-
erating the BBH initial data, are provided as ancillary material with the paper.
3.1. Single black hole evolutions
In this section, we use
SpECTRE
to evolve a single, stationary, black hole that, unless otherwise
noted, is non-spinning. We evolve a black hole from analytic initial data corresponding to a
black hole at rest centered at the origin with zero spin. We choose the mass of the black hole to
be
M
=
1 and work in units of
M
. In each evolution we use the following values for the FOGH
constraint damping parameters,
γ
0
=
γ
2
=
A
0
e
−
r
2
/
w
2
0
+
A
1
e
−
r
2
/
w
2
1
,
(49)
γ
1
=
−
1
,
(50)
where
r
is the coordinate distance from the origin,
A
0
=
7
.
0
/
M
,
A
1
=
0
.
1
/
M
,
w
0
=
2
.
5
M
, and
w
1
=
100
.
0
M
. The computational domain of each evolution covers a spherical shell volume
(figure
2
) with inner radius
r
in
=
r
exc
which differs for our different test cases, and outer bound-
ary coordinate radius
r
out
=
1000
M
. We apply boundary conditions as described in section
2.3
.
We use a fourth-order Adams–Moulton predictor–corrector time integrator with absolute and
relative time stepper tolerances of 10
−
8
and 10
−
6
, respectively, unless otherwise stated.
3.1.1. Kerr–Schild coordinates.
We first evolve a single black hole in Kerr–Schild coordin-
ates from Kerr–Schild initial data. The inner radius of the computational domain is
r
in
=
r
exc
=
1
.
9
M
. In this case there are no coordinate dynamics, so a feedback control system is not neces-
sary, though it is enabled in the simulations presented here. The left panel of figure
3
shows the
gauge constraint
C
a
and the 3-index constraint
C
iab
as a function of time for several different
resolutions. We evolve the lowest resolution to time
t
=
10 000
M
to assess long-term stability,
and we evolve the medium and high resolutions to assess convergence. To limit the computa-
tional cost of these tests, we choose to evolve the medium and high resolutions only to time
t
=
2000
M
. All simulations are stable to time
t
=
2000
M
, and the lowest resolution remains
stable to
t
=
10 000
M
. The amount of violation of the gauge constraint
C
a
and the 3-index
constraint
C
iab
is indicative of the overall constraint violation in the simulation. In this evol-
ution, the constraints remain approximately constant, and they decrease exponentially with
increasing
p
-refinement (that is, increasing points per cell per dimension), as expected. We
suspect the transient at
t
=
2000
Mt
=
2000
M
results from constraint violations reflecting off
the outer boundary back to the interior.
The right panel of figure
3
demonstrates long-term stability and convergence for the same
setup but with a black hole of dimensionless spin
χ
≡
S
/
M
2
=
0
.
8 (with
r
in
=
r
exc
=
1
.
57
M
).
11