of 15
Quasistationary hair for binary black hole initial data
in scalar Gauss-Bonnet gravity
Peter James Nee ,
1
,*
Guillermo Lara ,
1
,
Harald P. Pfeiffer ,
1
and Nils L. Vu
2
1
Max Planck Institute for Gravitational Physics (Albert Einstein Institute)
,
Am Mühlenberg 1, D-14476 Potsdam, Germany
2
Theoretical Astrophysics, Walter Burke Institute for Theoretical Physics,
California Institute of Technology
, Pasadena, California 91125, USA
(Received 14 June 2024; accepted 25 November 2024; published 24 January 2025)
Recent efforts to numerically simulate compact objects in alternative theories of gravityhave largelyfocused
on the time-evolution equations. Another critical aspect is the construction of constraint-satisfying initial data
with precise control over the properties of the systems under consideration. Here, we augment the extended
conformal thin sandwich framework to construct quasistationary initial data for black hole systems in scalar
Gauss-Bonnet theory and numerically implement it in the open-source
S
p
ECTRE
code. Despite the resulting
elliptic system being singular at black hole horizons, we demonstrate how to construct numerical solutions that
extend smoothly across the horizon. We obtain quasistationary scalar hair configurations in the test-field limit
for black holes with linear/angular momentum as well as for black hole binaries. For isolated black holes,
we explicitly show that the scalar profile obtained is stationary by evolving the system in time and compare
against previous formulations of scalar Gauss-Bonnet initial data. In the case of the binary, we find that the
scalar hair near the black holes can be markedly altered by the presence of the other black hole. The initial data
constructed hereenable targeted simulations in scalar Gauss-Bonnetsimulationswithreduced initial transients.
DOI:
10.1103/PhysRevD.111.024061
I. INTRODUCTION
Since thefirst gravitational wave (GW) event from a binary
black hole coalescence, GW150914
[1]
, the possibility of
testing our current theories of gravity against observational
GW data in the highly dynamical strong-field regime has
become a reality. To date, while general relativity (GR) has
been found to be consistent with current observations
[2
6]
,
strong field tests for theories beyond GR have not yet been as
thorough. In the context of GWs, this is mostly due to the
substantial effort required to compute the detailed predictions
needed to construct complete waveform models encompass-
ing all stages of compact binary coalescence. Crucially,
accurate modeling of the highly nonlinear late-inspiral and
merger stages relies on the ability to perform large-scale
numerical relativity (NR) simulations
[7]
.
In recent years, there has been growing interest in
extending the techniques of NR to alternative theories of
gravity. Such theories are often motivated by open issues
in gravity and cosmology, e.g., to provide a dynamical
explanation to the observed accelerated expansion of the
Universe or to connect GR to a more fundamental theory
of quantum gravity. For scalar tensor theories with two
propagating tensor modes and one scalar mode
[8
12]
,
interactions between the metric and a dynamical scalar may
lead to significant differences in the phenomenology of
compact binaries. For instance, in scalar Gauss-Bonnet
gravity (sGB), the component black holes (BHs) in the
binary may be endowed with
scalar hair
[13,14]
and
energy may be dissipated through radiation channels in
addition to the two GW polarizations of GR.
As is the case for GR, one can attempt to split the field
equations in alternative theories of gravity into two sets of
partial differential equations: a set of hyperbolic evolution
equations, such as the generalized harmonic equations
in GR; and a set of elliptic constraint equations, such
as the Hamiltonian and momentum constraints in GR.
Nevertheless, the mathematical structure of both sets of
equations differs from GR as the additional interactions
contribute to new terms in the principal part. In the strong
field regime, such new terms may alter the hyperbolic/
elliptic character of the equations and even lead to
equations of mixed type [see, e.g., Refs.
[15,16]
]. In this
respect, numerical relativity efforts have thus far focused on
finding appropriate formulations for the set of evolution
*
Contact author: peter.nee@aei.mpg.de
Contact author: glara@aei.mpg.de
Published by the American Physical Society under the terms of
the
Creative Commons Attribution 4.0 International
license.
Further distribution of this work must maintain attribution to
the author(s) and the published article
s title, journal citation,
and DOI. Open access publication funded by the Max Planck
Society.
PHYSICAL REVIEW D
111,
024061 (2025)
2470-0010
=
2025
=
111(2)
=
024061(15)
024061-1
Published by the American Physical Society
equations that allow for stable numerical evolutions. These
newly developed evolution strategies, which include novel
gauges
[17]
, traditional perturbation theory techniques and
proposals based on viscous hydrodynamics
[18]
, and their
numerical implementation, have already produced a num-
ber of successful merger simulations in alternative gravity
theories [see, e.g., Refs.
[19
36]
; for studies in the test-field
approximation; and the review
[37]
].
In this work, we take a step back and focus on the set of
elliptic constraint equations. Many of the current simulations
for compact binary objects in scalar tensor theories either
start off from initial data constructed for GR or use a
superposition of isolated solutions. While such approaches
are practical and useful for first qualitative explorations, they
are not guaranteed to satisfy the full constraint equations of
the extended theory and will in general not be in quasista-
tionary equilibrium. Indeed, constraint-satisfying solutions
can be obtained after an initial transient stage by employing
standard techniques, e.g., by including constraint-damping
terms or by smoothly turning on the additional interactions.
The cost, however, is a loss in control of the initial physical
parameters (e.g., mass, spin, eccentricity) during the relax-
ation stage (which may migrate to different values), as well
as the additional computational resources spent in simulating
this phase. If our aim is to efficiently obtain accurate
waveforms and to adequately cover the parameter space
for the calibration of waveform models, experience with GR
has shown that constructing constraint-satisfying initial data
in quasistationary equilibrium is important.
In GR, the most common way of formulating the
Hamiltonian and momentum constraints as a set of elliptic
equations is the conformal method, where instead of
solving for geometric quantities directly one performs a
conformal decomposition
[38]
. This is the basis for two of
the most well-known approaches, namely, the conformal
transverse traceless (CTT)
[39]
and the extended conformal
thin sandwich (XCTS) methods
[40,41]
.
For the case of alternative theories, Kovacs
[42]
has
recently examined the mathematical properties of the
elliptic systems arising in weakly coupled four-derivative
scalar tensor theories (a class of theories which includes the
sGB theory investigated here) and provides theorems
regarding the well-posedness of the boundary value prob-
lem using extensions of the CTT and XCTS methods.
On the practical side, several authors have constructed
constraint-satisfying initial data for compact binaries in
theories beyond GR. Considering four-derivative scalar
tensor theory, Ref.
[43]
prescribes an
ad hoc
scalar field
configuration, solving the constraint equations via a modi-
fication of the CTT approach
[44]
, in which the elliptic
equation for the conformal factor is reinterpreted as an
algebraic one for the mean curvature. While the initial data
constructed in this way is constraint-satisfying, since the
scalar hair configuration is not in quasistationary equilib-
rium, it should be expected to lead to significant transients
during the initial stage of evolution. A similar numerical
approach is taken in Refs.
[45,46]
to obtain constraint
satisfying initial data for boson star binaries, where the
constraints are solved for free data specified by the super-
position of isolated boson stars. In the context of Damour-
Esposito-Far`
ese theory
[47]
for neutron star binaries,
Ref.
[48]
solved the constraints for the metric alongside
an additional Poisson equation for the scalar field.
This paper develops and implements a method to
construct constraint-satisfying initial data where the scalar
field is
in equilibrium
. We focus on the decoupling limit
(i.e., the scalar does not back-react onto the metric) of
scalar Gauss-Bonnet gravity in vacuum
S
½
g
ab
;
Ψ

Z
d
4
x
ffiffiffiffiffiffi
g
p

R
2
κ
1
2
a
Ψ
a
Ψ
þ
l
2
f
ð
Ψ
Þ
G

;
ð
1
Þ
where
κ
8
π
G
,
l
denotes the coupling constant,
g
¼
det
ð
g
ab
Þ
is the determinant of the metric
g
ab
, and
Ψ
is the scalar field. To obtain spontaneously scalarized
BHs
[49
51]
, we choose the free function
f
ð
Ψ
Þ
as
[52]
FIG. 1.
Scalar hair around a black hole binary
. Left: opposite sign scalar charges. Right: scalar charges with same sign. This
configuration corresponds to an equal mass, nonspinning binary black hole system in an approximate circular orbit. The sGB coupling
constants are
l
2
η
3
.
34
M
2
,
l
2
ζ
≃−
31
.
1
M
2
, where
M
denotes the mass of either BH, and the BH separation is
D
28
M
.
NEE, LARA, PFEIFFER, and VU
PHYS. REV. D
111,
024061 (2025)
024061-2
f
ð
Ψ
Þ
η
8
Ψ
2
þ
ζ
16
Ψ
4
:
ð
2
Þ
This function couples
Ψ
to the Gauss-Bonnet scalar
G
R
abcd
R
abcd
4
R
ab
R
ab
þ
R
2
;
ð
3
Þ
which is in turn defined in terms of the Riemann tensor
R
abcd
, the Ricci tensor
R
ab
, and the Ricci scalar
R
.
Following Ref.
[42]
, we revisit the conditions for
obtaining quasistationary configurations for the scalar hair
around isolated and binary black holes. We argue that, in
the initial data slice, one must impose a vanishing scalar
momentum
defined in terms of the directional derivative
along an approximate Killing vector of the spacetime
as
opposed to the directional derivative along the normal to
the foliation as in Ref.
[42]
. The adapted coordinates from
the background spacetime given by a solution to the XCTS
equations naturally yield the required approximate Killing
vector. Imposing the appropriate momentum condition on
the scalar equation, we derive a
singular
boundary-value
problem for BH spacetimes.
We demonstrate that this singular boundary-value prob-
lem can be solved without an inner boundary condition
in the spectral elliptic solver of the open-source
S
p
ECTRE
code
[53]
. We thus obtain quasistationary hair for both
single and binary black hole spacetimes, as illustrated in
Fig.
1
. Moreover, for the case of single BHs, we further
evolve the obtained configuration to confirm that the
solution is indeed quasistationary and does not lead to
large transients, and compare against the prescription given
in Ref.
[42]
.
This paper is organized as follows. Section
II
recalls
basic aspects of sGB theory and of the XCTS method. In
Sec.
III
, we revisit different formulations for the scalar
equation and define a scheme that imposes quasiequili-
brium on the scalar hair. We further discuss the singular
boundary value problem and describe our numerical
implementation to solve for single BHs. Section
IV
show-
cases the different single BH systems that we can access
with both the 1D and 3D code, and demonstrates the
stationarity of the solutions obtained via direct evolution.
Section
V
constructs initial data for binary black holes
with scalar hair. We first deal with conceptual issues
regarding scalar configurations on arbitrarily large spatial
domains and then proceed to present our solutions for
quasistationary scalar hair. We summarize and discuss our
results in Sec.
VI
. Throughout this paper, we use geometric
units such that
G
¼
c
¼
1
and
ð
þ þþÞ
signature. Early
alphabet letters
f
a; b; c
g
represent 4-dimensional space-
time indices, while middle alphabet letters
f
i; j; k
g
corre-
spond to 3-dimensional spatial indices.
II. THEORY
Variation of the action of scalar Gauss-Bonnet theory
[Eq.
(1)
] yields a scalar equation
Ψ
¼
l
2
f
0
ð
Ψ
Þ
G
ð
4
Þ
and a tensor equation
R
ab
¼
H
ab
½
g
cd
;
Ψ

;
ð
5
Þ
where
H
ab
½
g
cd
;
Ψ

contains up to second derivatives of
g
ab
and
Ψ
see, e.g., Ref.
[54]
for the full expression. In the
decoupling limit of the theory (i.e., when
Ψ
is considered a
test field), the right-hand-side of Eq.
(5)
vanishes,
H
ab
½
g
cd
;
Ψ

0
.
A. Spontaneous scalarization
Stationary BH solutions of Eqs.
(4)
and
(5)
are often
nonunique. When
f
0
ð
0
Þ¼
0
, as for our choice of
f
ð
Ψ
Þ
,a
GR solution with
Ψ
0
trivially solves Eqs.
(4)
and
(5)
.
However, GR solutions can be energetically disfavored for
a large enough coupling parameter
l
2
η
0
.Thiscanbe
seen
[55]
by expanding around
Ψ
0
to derive an
equation describing the scalar perturbations around the
GR solution,
ð
m
2
Ψ
;
eff
Þ
δ
Ψ
¼
0
;
ð
6
Þ
where
m
2
Ψ
;
eff
l
2
η
G
plays the role of an effective,
spatially varying mass term. If
m
2
Ψ
;
eff
is negative enough,
due to either high curvature or very large spins
[56,57]
,
GR solutions in sGB may become dynamically unstable,
and will
spontaneously scalarize
to yield a second set of
solutions with nonvanishing scalar
hair
[49
51,58]
.
Therefore, BHs in sGB theory are characterized by their
mass, spin, and an additional
scalar charge
parameter
q
,
defined by the asymptotic behavior of the scalar as
Ψ
ð
r
Þ¼
Ψ
þ
qM
2
r
þ
O

1
r
2

;
ð
7
Þ
where
Ψ
is the asymptotic value of the scalar field and
M
is the mass of the BH. Given the
Ψ
Ψ
symmetry of the
theory described by Eqs.
(1)
and
(2)
, any hairy solutions
will have a corresponding equivalent solution related by
Ψ
Ψ
, and which is characterized by a scalar charge of
equal magnitude and opposite sign.
B. The XCTS formulation
In the decoupling limit, the constraint equations arising
from Eq.
(5)
are the usual Hamiltonian and momentum
constraints of GR. To obtain them, we perform a (
3
þ
1
)-
decomposition of the metric,
ds
2
¼
g
ab
dx
a
dx
b
¼
α
2
dt
2
þ
γ
ij
ð
β
i
dt
þ
dx
i
Þð
β
j
dt
þ
dx
j
Þ
;
ð
8
Þ
QUASISTATIONARY HAIR FOR BINARY BLACK HOLE
...
PHYS. REV. D
111,
024061 (2025)
024061-3
where
α
is the lapse,
β
i
¼
γ
ij
β
j
is the shift, and
γ
ij
is
the spatial metric (with inverse
γ
ij
). The constraints in
vacuum read
[7]
ð
3
Þ
R
þ
K
2
K
ij
K
ij
¼
0
;
ð
9a
Þ
D
j
ð
K
ij
γ
ij
K
Þ¼
0
;
ð
9b
Þ
where
ð
3
Þ
R
denotes the Ricci scalar of
γ
ij
and
D
j
is the
3-dimensional covariant derivative compatible with
γ
ij
.
Finally,
K
ab
ð
1
=
2
Þ
L
n
γ
ab
denotes the extrinsic curva-
ture, with trace
K
, where the Lie derivative is taken along
the future-pointing unit normal to the foliation,
n
a
.
We further decompose the spatial metric as
γ
ij
¼
ψ
4
̄
γ
ij
;
ð
10
Þ
where
ψ
>
0
is the conformal factor and
̄
γ
ij
is the
conformal spatial metric, which we are free to specify.
The XCTS formalism
[40,41]
is centered around specifying
certain free data and their
time derivatives
. Specifically,
the conformal metric
̄
γ
ij
and
K
are free data, as well as
t
̄
γ
ij
̄
u
ij
and
t
K
. It is useful to decompose the extrinsic
curvature as
K
ij
¼
1
3
γ
ij
K
þ
ψ
10
̄
A
ij
ð
11
Þ
with
̄
A
ij
¼
1
2
̄
α
½ð
̄
L
β
Þ
ij
̄
u
ij

:
ð
12
Þ
Here,
̄
α
¼
ψ
6
α
is the conformal lapse function, and the
conformal longitudinal operator is defined as
ð
̄
L
β
Þ
ij
¼
2
̄
D
ð
i
β
j
Þ
2
3
̄
γ
ij
̄
D
k
β
k
;
ð
13
Þ
where
̄
D
i
denotes the covariant derivative operator com-
patible with the conformal metric
̄
γ
ij
.
The final XCTS equations are then obtained from
Eqs.
(9)
and from the evolution equation for
K
, and are
given by
[41]
̄
D
2
ψ
1
8
ψ
ð
3
Þ
̄
R
1
12
ψ
5
K
2
þ
1
8
ψ
7
̄
A
ij
̄
A
ij
¼
0
;
ð
14a
Þ
̄
D
i

1
̄
α
½ð
̄
L
β
Þ
ij
̄
u
ij


2
3
ψ
6
̄
D
j
K
¼
0
;
ð
14b
Þ
̄
D
2
ð
αψ
Þ
αψ

7
8
ψ
8
̄
A
ij
̄
A
ij
þ
5
12
ψ
4
K
2
þ
1
8
ð
3
Þ
̄
R

þ
ψ
5
ð
t
K
þ
β
i
̄
D
i
K
Þ¼
0
;
ð
14c
Þ
where
ð
3
Þ
̄
R
is the spatial conformal Ricci scalar. In the XCTS
formalism, the notion of quasistationary equilibrium can be
imposed
[59]
by demanding that the conformal metric and
trace of the extrinsic curvature remain unchanged along
infinitesimally separated spatial slices, i.e.,
̄
u
ij
¼
0
;
ð
15
Þ
t
K
¼
0
:
ð
16
Þ
Combined with appropriate boundary conditions [see
Ref.
[59]
for details], the XCTS system [Eqs.
(14)
] is then
solved for
f
ψ
;
αψ
;
β
i
g
, thus providing not only a solution to
the constraint equations
(9)
, but also a coordinate system
adapted to the symmetry along the approximate Killing
vector
t
a
a
¼ð
α
n
a
þ
β
a
Þ
a
¼
t
:
ð
17
Þ
In Sec.
III
, we will extend this property to the scalar
equation in sGB theory.
III. QUASISTATIONARY SCALAR HAIR
In this section, we revisit the scalar equation Eq.
(4)
and
consider different strategies to include it in the XCTS
scheme. The aim is to obtain solutions for the metric
and
the scalar hair of the BH in a general 3-dimensional space
without symmetry. We further describe our numerical
implementation, which will also be applicable to the more
general case of BH binaries treated in Sec.
V
.
A. Spherical symmetry
We first consider a spherically symmetric BH in horizon-
penetrating Kerr-Schild coordinates
ds
2
¼

1
2
M
r

dt
2
þ
4
M
r
dtdr
þ

1
þ
2
M
r

dr
2
þ
r
2
d
Ω
2
;
ð
18
Þ
with
d
Ω
2
¼
d
θ
2
þ
sin
2
ð
θ
Þ
d
φ
2
. Under the assumption
that the scalar field is time-independent, Eq.
(4)
yields

1
2
M
r

2
r
Ψ
2
ð
M
r
Þ
r
2
r
Ψ
¼
48
M
2
r
6
f
0
ð
Ψ
Þ
;
ð
19
Þ
where
G
¼
48
M
2
=r
6
and where
M
is the mass of the BH.
We will be looking for solutions of Eq.
(19)
with asymp-
totic behavior
(7)
by imposing
1
1
While one can easily place the outer boundary
r
max
at spatial
infinity in the spherically symmetric case, we impose the condition
(20)
to connect with the 3D-implementation in Sec.
IVA
.
NEE, LARA, PFEIFFER, and VU
PHYS. REV. D
111,
024061 (2025)
024061-4
½
r
r
Ψ
þ
Ψ
Ψ

r
¼
0
;
ð
20
Þ
with
Ψ
¼
0
.
This is our first encounter with a
singular
boundary
value problem. Notice that Eq.
(19)
is singular at the BH
horizon
r
h
¼
2
M
, where the factor in front of the highest-
derivative operator vanishes. Despite this observation,
Eq.
(19)
can be easily solved via the shooting method
[60]
. Regularity of
Ψ
at the horizon can be imposed by
expanding
Ψ
as an analytic series around
r
h
. The solutions
satisfying Eq.
(20)
can then be found by numerically
integrating outward starting from
r
h
þ
ε
and performing
a line search in the unknown value
Ψ
j
r
h
at the inner
boundary.
In order to prepare for our later 3D solutions, we shall
take a different approach, and solve Eq.
(19)
on a domain
that extends through the horizon into the interior region.
Utilizing a spectral method, we represent
Ψ
as a series in
Chebychev polynomials
T
i
ð
x
Þ
,
Ψ
ð
x
Þ¼
X
N
i
¼
0
Ψ
ð
i
Þ
T
i
ð
x
Þ
;
ð
21
Þ
where the argument
x
½
1
;
1

is related to radius
r
by the
transformation
x
¼
A=
ð
r
B
Þþ
C
for suitable constants
A
,
B
, and
C
. To cover
r
½
r
min
;r
max

, we set
A
¼
ð
r
min
þ
r
max
þ
2
C
Þ
=
ð
r
max
r
min
Þ
,
B
¼
r
max
Ar
max
þ
C
AC
, and leave
C
as a specifiable constant to adjust
the distribution of resolution throughout the interval.
We choose a spatial grid
f
x
i
g
N
i
¼
0
defined by the nodes
(or zeros) of
T
i
ð
x
Þ
, and compute spatial derivatives of
Ψ
analytically from Eq.
(21)
. Using a Newton-Raphson
scheme, we iteratively solve the scalar equation by expand-
ing
Ψ
Ψ
þ
δ
Ψ
and linearizing Eq.
(19)
. We obtain

1
2
M
r

2
r
δ
Ψ
ð
K
Þ
2
ð
M
r
Þ
r
2
r
δ
Ψ
ð
K
Þ
þ
48
M
2
r
6
f
00
ð
Ψ
ð
K
Þ
Þ
δ
Ψ
ð
K
Þ
¼

1
2
M
r

2
r
Ψ
ð
K
Þ
þ
2
ð
M
r
Þ
r
2
r
Ψ
ð
K
Þ
48
M
2
r
6
f
0
ð
Ψ
ð
K
Þ
Þ
;
ð
22
Þ
where, at a given iteration step
K
, the improved solution is
given by
Ψ
ð
K
þ
1
Þ
¼
Ψ
ð
K
Þ
þ
δ
Ψ
ð
K
Þ
:
ð
23
Þ
For a solution interval crossing the horizon, i.e.,
r
min
<r
h
<r
max
, we impose boundary conditions of the
form
(20)
only at the outer boundary. We do
not
impose
regularity across the entire domain (in particular, at a
singular boundary at
r
¼
r
h
) via boundary conditions, as
it is already built into the spectral expansion
(21)
since all
Chebychev polynomials are regular. We implement this
algorithm in
P
ython
, and for each iteration step
K
we solve
the discretized version of Eq.
(22)
via explicit matrix
inversion using
N
um
P
y
. An exemplary solution of
Eq.
(19)
is shown as the blue line in Fig.
2
, where we
set
r
min
¼
1
.
9
M
and
r
max
¼
10
10
M
.
B. 3D normal formulation
n
=0
To solve for scalar hair in a general 3-dimensional space,
Ref.
[42]
requires the
momentum
Π
n
a
a
Ψ
ð
24
Þ
to vanish everywhere on the initial spatial slice at
t
¼
0
, i.e.,
Π
j
t
¼
0
0
. The scalar equation
(4)
then becomes
i
ð
γ
ij
j
Ψ
Þþ
γ
ij
j
Ψ
ð
i
ln
α
þ
Γ
k
ki
Þ¼
l
2
f
0
ð
Ψ
Þ
G
;
ð
25
Þ
where
Γ
k
ij
is the 3-dimensional spatial Christoffel symbol
with respect to
γ
ij
. Equation
(25)
is both elliptic and
regular everywhere. In Ref.
[42]
, the inner boundary
S
in
is
placed on the apparent horizon of the BH and is supple-
mented with boundary conditions at both inner and outer
boundaries,
ˆ
s
i
i
Ψ
j
S
in
¼
0
;
ð
26
Þ
lim
r
Ψ
¼
Ψ
;
ð
27
Þ
FIG. 2.
Radial profile of the scalar field around a nonrotating
BH
. The curves correspond to different formulations used for
obtaining
Ψ
. We set
l
2
η
=M
2
¼
6
,
l
2
ζ
=M
2
¼
60
. The vertical
gray line indicates the location of the horizon.
QUASISTATIONARY HAIR FOR BINARY BLACK HOLE
...
PHYS. REV. D
111,
024061 (2025)
024061-5
where
ˆ
s
i
is the unit outward normal vector to the BH
horizon(s). For computational domains extending inside
the apparent horizon, we instead impose a constant
Dirichlet boundary condition (i.e.,
Ψ
j
S
in
¼
const.), chosen
such that
ˆ
s
i
i
Ψ
¼
0
on each apparent horizon. On a finite
spatial domain, and assuming an asymptotic decay of
the scalar of the form of Eq.
(7)
, we replace the outer
boundary condition with [c.f. Eq.
(20)
] a Robin type
boundary condition
ð
r
ˆ
s
i
i
Ψ
þ
Ψ
Ψ
Þj
S
out
¼
0
;
ð
28
Þ
where
ˆ
s
i
is now the unit outward normal vector to the
outer spherical boundary
S
out
.Weset
Ψ
0
.Weshall
refer to this as the
n
¼
0
formulation.
1. Caveats of the normal formulation
While the
n
¼
0
formulation provides a readily solvable
elliptic system, the most common use case for the XCTS
formulation is the calculation of quasiequilibrium initial
conditions. Unfortunately, the normal formulation will not
generically lead to stationary spacetimes. Consider, for
example, the case of a Schwarzschild BH in Kerr-Schild
coordinates [Eq.
(18)
]. The timelike Killing vector
ξ
of the
spacetime is
ξ
a
¼
t
a
¼
α
n
a
þ
β
a
:
ð
29
Þ
Assuming that the momentum
Π
is initially zero, the initial
time derivative of the scalar field
Ψ
is
t
c
c
Ψ
¼
α
n
c
c
Ψ
þ
β
i
i
Ψ
¼
β
i
i
Ψ
:
ð
30
Þ
Therefore, whenever
β
i
0
and
i
Ψ
0
,
L
ξ
Ψ
will not
vanish. For our example,
L
ξ
Ψ
¼
β
r
r
Ψ
¼
2
M=
ð
r
þ
2
M
Þ
r
Ψ
0
and the scalar hair obtained will not be
stationary. Indeed, solving the spherically symmetric
version of Eq.
(25)
in our 1D code, we find a profile
different from the
t
¼
0
solution constructed in Sec.
III A
.
This profile is also shown in Fig.
2
.
Finally, we note that the inner boundary condition
[Eq.
(26)
] is inconsistent with stationarity. Indeed, if
Ψ
is regular at the horizon, then it can be expanded as a series
about
r
h
of the form
Ψ
ð
r
Þ¼
X
n
¼
0
Ψ
ð
i
Þ
ð
r
r
h
Þ
n
:
ð
31
Þ
Solving Eq.
(19)
order-by-order perturbatively in
Δ
r
¼
r
r
h
,
we obtain that
r
Ψ
j
r
h
¼
Ψ
ð
1
Þ
¼
3
8
M
3
l
2
ð
η
Ψ
ð
0
Þ
þ
ζ
Ψ
3
ð
0
Þ
Þ
;
ð
32
Þ
which is nonzero in general, contradicting Eq.
(26)
.
C. 3D approximate Killing formulation
t
=0
Motivated by the existence of a symmetry along a Killing
vector, we present a new procedure for extending the XCTS
formulation to sGB gravity, which we refer to as the
t
¼
0
formulation. The main assumption will now be
that the
momentum
with respect to the (approximate)
Killing vector
ξ
a
,givenby
P
L
ξ
Ψ
;
ð
33
Þ
vanishes on the initial slice.
From the previous discussion, for a stationary GR black
hole in coordinates adapted to the symmetry, as well as for
solutions of the XCTS equations, the Killing vector
corresponds to
ξ
¼
t
. By imposing
P
¼
t
Ψ
0
, Eq.
(4)
becomes
i
ð
M
ij
j
Ψ
Þþð
i
ln
α
þ
Γ
k
ki
Þ
M
ij
j
Ψ
¼
l
2
f
0
ð
Ψ
Þ
G
;
ð
34
Þ
where
M
ij
γ
ij
α
2
β
i
β
j
:
ð
35
Þ
Equation
(34)
is the 3D generalization of Eq.
(19)
. In the
spirit of quasiequilibrium, we have also set
t
α
¼
0
and
t
β
i
¼
0
in the derivation of Eq.
(34)
. We note that these
simplifications could be relaxed and their values can be set
according to a desired gauge choice.
The principal part of Eq.
(34)
is
M
ij
i
j
Ψ
. The singu-
larity at
r
h
in the 1D formulation [Eq.
(19)
] now corre-
sponds to the situation where
det
M
j
S
h
¼
0
;
ð
36
Þ
i.e., when (at least) one of the eigenvalues of
M
ij
vanishes
at the apparent horizon
S
h
.
M
is singular on the BH horizon
in general. For a stationary BH in time-independent
coordinates, the time vector on the horizon must be parallel
to the horizon generators as argued in Ref.
[59]
, which
implies that on the horizon
β
i
ˆ
s
i
¼
α
, where
ˆ
s
i
is the
outward-pointing spatial unit normal to the horizon.
Using this equality, it follows that
ˆ
s
i
M
ij
ˆ
s
j
¼
0
.
As for the spherically symmetric example above, our
approach will be to rely on the inherent smoothness of
spectral expansions to single out solutions of Eq.
(34)
that
smoothly pass through the horizon. Regularity at the
horizon reduces the number of possible solutions, and so
we will not impose a boundary condition at the excision
surface in the interior of the horizon. We note that Lau
et al.
[61,62]
encountered the same principal part as Eq.
(34)
in
the context of IMEX evolutions on curved backgrounds.
Reference
[61]
in particular contains an analysis of the
singular boundary value problem. We impose the boundary
NEE, LARA, PFEIFFER, and VU
PHYS. REV. D
111,
024061 (2025)
024061-6