Scalarization of isolated black holes in scalar Gauss-Bonnet theory in the fixing-the-equations approach

Scalarization of isolated black holes in scalar Gauss-Bonnet theory

in the fixing-the-equations approach

Guillermo Lara ,

Harald P. Pfeiffer ,

Nikolas A. Wittek ,

Nils L. Vu ,

Kyle C. Nelli ,

Alexander Carpenter ,

Geoffrey Lovelace ,

Mark A. Scheel ,

and William Throwe

Max Planck Institute for Gravitational Physics

(Albert Einstein Institute), Am Mühlenberg 1, 14476 Potsdam, Germany

Theoretical Astrophysics, Walter Burke Institute for Theoretical Physics,

California Institute of Technology

, Pasadena, California 91125, USA

Nicholas and Lee Begovich Center for Gravitational-Wave Physics and Astronomy,

California State University Fullerton

, Fullerton, California 92831, USA

Cornell Center for Astrophysics and Planetary Science,

Cornell University

, Ithaca, New York 14853, USA

(Received 22 March 2024; accepted 14 June 2024; published 15 July 2024)

One of the most promising avenues to perform numerical evolutions in theories beyond general relativity

is the

fixing-the-equations

approach, a proposal in which new

“

driver

”

equations are added to the evolution

equations in a way that allows for stable numerical evolutions. In this direction, we extend the numerical

relativity code

ECTRE

to evolve a

“

fixed

”

version of scalar Gauss-Bonnet theory in the decoupling limit, a

phenomenologically interesting theory that allows for hairy black hole solutions in vacuum. We focus on

isolated black hole systems both with and without linear and angular momentum, and propose a new driver

equation to improve the recovery of such stationary solutions. We demonstrate the effectiveness of the latter

by numerically evolving black holes that undergo spontaneous scalarization using different driver

equations. Finally, we evaluate the accuracy of the obtained solutions by comparing with the original

unaltered theory.

DOI:

10.1103/PhysRevD.110.024033

I. INTRODUCTION

The increasing availability of gravitational wave (GW)

data, by current (LIGO-Virgo-KAGRA

–

) and future

(LISA

[4]

, Einstein Telescope

[5]

, Cosmic Explorer

[6]

)

ground and space-based interferometers, promises to lead to

the strictest tests yet of Einstein

’

s general theory of relativity

(GR) as the fundamental description of gravitational phe-

nomena. Adding on to the vast array of weak-field tests

[7]

growing suite of model-independent tests on the GW signal

from binary compact coalescences (BCCs)

–

10]

show no

apparent disagreement with GR to date. At the same time,

direct observation by means of very-large-baseline interfer-

ometry show black hole (BH) images consistent with expect-

ations based on Einstein

’

s gravity

[11,12]

Nevertheless, it is conceivable that GR is not the ultimate

description of gravity. On the theoretical side, the quantum

nature of all matter and force fields in the Standard Model,

firmly established in the past half-century, suggests the

existence of a quantum theory of gravity. From this

perspective, GR arises as an

effective field theory

(EFT)

with limited range of applicability and is subject to

corrections relevant at shorter (higher) length (curvature)

scales. On the phenomenological side, modifying gravity at

the largest scales (e.g. by introducing new light degrees of

freedom) could help explain cosmological issues such as

the observed self-accelerated expansion of the Universe (or

dark energy)

[13]

. Among other motivations, numerous

extensions of GR have been devised. Some of the most

important classes of theories include degenerate higher-

order scalar-tensor theories

[14

–

16]

—

a generalization of

the Horndeski class

[17]

–

and EFT expansions of gravity

—

see e.g. Refs.

[18,19]

It is expected that GW signals will carry information

about possible GR corrections. Already, tight constraints on

the speed of tensor modes from GW170917

[20,21]

have

cast widely encompassing constraints on dark energy

models

[22

–

25]

. However, in order to fully exploit the

potential of GWs to cast the strictest bounds on beyond-GR

parameters, model-dependent tests need to be carried out

for (if only a few) alternative theories of gravity. In this

direction, there is a need to construct full waveform models,

encompassing all stages of BCC. For GR, numerical

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

110,

024033 (2024)

2470-0010

2024

110(2)

024033(16)

024033-1

Published by the American Physical Society

simulations of the late-inspiral and merger phases of

compact binary coalescence are essential to construct

and calibrate models describing these highly dynamical

and nonlinear stages

[26

–

29]

. Therefore, a key (and

challenging) step is to extend the techniques of numerical

relativity (NR)

[30]

beyond GR.

As it was the case for GR decades ago, one of the main

difficulties is the mathematical structure of equations of

motion. These systems of partial differential equations

(PDEs) need to be recast (if possible) in a way suitable

for numerical evolution, such that the system admits a

well-

posed

initial-value problem (IVP)

—

a property that ensures

good behavior of the solutions: uniqueness and continuity

on the space of initial data

[31]

. Failure to satisfy this

property can give rise to unstable solutions and even

seemingly pathological behavior

—

e.g. hyperbolic PDE

systems have been observed to change character to elliptic

in regions outside BH horizons

[32

–

40]

. Nevertheless,

recent developments have fuelled renewed interest and a

steady increase in beyond-GR codes and GW-waveform

examples

—

see Ref.

[41]

for a review. Indeed, some of us

have obtained beyond-GR waveforms in dynamical Chern-

Simons and scalar Gauss-Bonnet (sGB) theory using an

order-reduction scheme to obtain perturbative corrections

to the GW signal

[42

–

45]

—

see also Refs.

[46

–

48]

. Other

studies have adopted theory-specific approaches, in some

cases involving explorations of the beyond-GR parameter

space

[49

–

56]

. More recently, novel modified generalized

harmonic (GH) gauges applicable for a broad class of

weakly coupled Lovelock and Horndeski theories have

been proposed by Kovacs and Reall

[57,58]

. The latter have

been numerically implemented in Refs.

[32,59,60]

, as well

as in similar extensions to the conformal and covariant Z4

(CCZ4) system

[61,62]

, to evolve binary systems in

sGB EFTs.

Here we will explore a fourth avenue, sometimes referred

to in the literature as the

fixing-the-equations

approach

[63]

. In this proposal (inspired by Müller-Israel-Stewart

hydrodynamics), the equations of motion are deformed (or

“

fixed

”

) by the introduction of auxiliary variables subject to

suitable

driver

evolution equations. New timescales appear-

ing in these driver equations control how the fixed solutions

track the solutions of the original system, thus allowing for

the evolution of the scales relevant to the problem, while

“

softening

”

(possibly) problematic high-frequency modes.

In this paper, we report on the first steps to implement the

fixing-the-equations approach in the numerical relativity

code

ECTRE

[64]

, developed by the SXS Collaboration.

While this method can be applied to a wide variety of

theories (see e.g. Refs.

[52,65

–

70]

), for concreteness and

motivated by positive results in spherical symmetry

[68]

we will focus on sGB gravity in vacuum. The action is

;

≡

ffiffiffiffiffiffi

−

∇

Þþ

l

;

where

≡

and

l

are constant couplings,

det

is the determinant of the metric

, the scalar field

is a free function describing the specific

coupling to the Gauss-Bonnet scalar

≡

abcd

−

;

which is defined in terms of combinations of the Riemann

tensor

abcd

, the Ricci tensor

and the Ricci scalar

Scalar Gauss-Bonnet theory is of phenomenological inter-

est mainly because it is a model for black holes endowed

with scalar

hair

[71,72]

. As opposed their GR counterparts

(described by the Kerr metric

[73]

), BHs in sGB theory

evade common no-hair theorems

[74

–

77]

and are charac-

terized by a

scalar charge

in addition to their mass and spin

parameters. BHs can acquire scalar hair through a variety of

scalarization

mechanisms (see Ref.

[78]

for a review) and it

is expected that this hair will significantly impact the GW

signal from binary black hole (BBH) systems. Recently,

model-dependent bounds using the inspiral phase of GW

observations have been placed on shift-symmetric sGB

models

—

giving values of

l

≲

[79,80]

. Addi-

tional theoretical constraints on the form of the coupling

function

l

have been derived from the assumption of

consistency ofsGB theory withawell-behaved(highenergy)

ultraviolet completion

[81]

In our implementation of the fixing-the-equations

approach, we consider the following (related) issues:

(i) Can we accurately describe the stationary solutions of

the theory? (ii) Can we be confident that the fixed equations

are able to reproduce the original physics? We address both

questions empirically in a case study: the scalarization of

isolated BHsin the

decoupling limit

oftheory

(1)

—

wherethe

scalar field backreaction on the metric is neglected. The

reason why this approximation isillustrative istwofold. First,

it will help us test the properties of different driver equations

in a localized sector: the scalar sector. Second, the

“

unmodi-

fied

”

decoupled theory can be recast in strongly hyperbolic

form and thus readily evolved so as to compare different

“

fixing

”

prescriptions against the

“

correct

”

solution. Finally,

as a main result of this paper, we present a new driver

equation that leads to an improved recovery of the stationary

solutions of the theory, with respect to other equations

studied in the literature

–

see Fig.

This paper is organized as follows. In Sec.

,wegivea

brief theoretical background of sGB theory and the fixing-

the-equations approach. We describe our numerical imple-

mentation and diagnostics in Sec.

III

.InSec.

, we illustrate

with different examples of isolated BH scalarization how

stationary solutions are recovered by the fixed systems and

how the new driver equation improves their accuracy using a

variety of diagnostics. We summarize our results and com-

ment on the future steps of this program in Sec.

. Further

material is presented in the Appendices. In Appendix

the

driver equations of the main text are compared with a third

GUILLERMO LARA

et al.

PHYS. REV. D

110,

024033 (2024)

024033-2

wavelike driver. Finally, in Appendix

, convergence plots

are presented. Throughout this paper we set

and

work in

−

þþþ

signature. Early alphabet letters

a; b; c;

...

represent 4-dimensional spacetime indices,

whereas middle letters

i; j; k;

...

represent 3-dimensional

spatial indices.

II. THEORY

Taking the variation of action

(1)

with respect to the

scalar field

, we obtain the scalar equation

□

−

l

Variation with respect to the metric yields

l

;

where

≡

−

is the Einstein tensor,

≡∇

∇

−

∇

is the canonical stress-energy tensor of the scalar, and

[68]

≡

−

acbd

∇

;

with

abcd

≡

abcd

−

In the

decoupling limit

of the theory, the backreaction of

the scalar field on Eq.

(4)

is neglected. The resulting system

is equivalent to a test scalar field, described by Eq.

(3)

which evolves on a known GR background (see e.g.

[47,48]

□

−

l

A. Spontaneous scalarization

In order to describe BHs that undergo

spontaneous

scalarization

, we will consider the model of Ref.

[82]

≡

;

where

and

are dimensionless parameters. As can be

readily seen from Eq.

(9)

and action

(1)

, the resulting theory

symmetric and solutions related by

→

−

are

equivalent.

A trivial solution to Eqs.

(3)

and

(4)

is a Kerr metric with

≡

. Thus, the solutions of GR are also solutions for

model

(9)

of sGB theory. Nevertheless, a

tachyonic

instability in the scalar sector may render GR solutions

linearly unstable

[78]

. This can be readily seen in the

decoupling limit by writing the equation for a linear

perturbation

of the scalar, which reads

□

−

;

eff

where

;

eff

≡

−

l

−

l

. Depending

on the combined signs of the coupling

and the Gauss-

Bonnet scalar

, unstable regions where

;

eff

may

exist in vicinity of the BH. For a Schwarzschild BH,

, where

is the mass of the BH, and

thus

≳

suffices to create an instability region leading to

the growth of the scalar

[83]

. More generally, when the BH

is spinning,

may change sign around the poles. For

≳

a similar situation to the spinless case occurs near the

equator, whereas for

≲

, and large enough dimensionless

spin parameter

≳

spin-induced

instability can occur

in the region near the poles

[84]

As this instability triggers an exponential growth of the

scalar in the vicinity of the BH, the nonlinear terms in

will eventually become as large as the linear term.

Spontaneous scalarization for BHs (originally discovered

for neutron stars in Damour-Esposito-Far`

ese theory

[85]

)in

sGB theory occurs when nonlinearities are such that at a

later stage the instability is quenched and the system can

settle to a stable configuration with nontrivial scalar hair

[83,86,87]

—

as illustrated in Fig.

; see also Ref.

[78]

for a

review. For the model

(9)

the transition to the final

configuration will depend on the parameter

since the

nonlinear terms associated to it will be dominant for large

. A back-of-the-envelope calculation suggests that the

end of the tachyonic instability (where the effective mass

vanishes) occurs when the scalar reaches a final (maxi-

mum) scalar amplitude of order

∼

ffiffiffiffiffiffiffiffiffiffiffi

−

, for

with

opposite sign to

. For stability, the value of

is further

required to satisfy

≲

−

[82,88]

. The parameter

space for such hairy black hole solutions can be charac-

terized in terms of the black hole mass

and correspond-

ing

scalar charge

—

see e.g. Ref.

[89]

for the parameter

space in a similar exponential model.

Herewe define the scalar charge (per unit mass squared) as

≡

−

lim

→

∞

ˆ

∂

;

where the integral is taken over a sphere of radius

,and

ˆ

is the outward normal to the sphere. This definition of the

scalar charge coincides with the usual definition through

the scalar falloff at infinity, i.e.

→

∞

Þ¼

∞

−

—

see e.g. Ref.

[89]

. Notice also that,

for a generic model

, the scalar charge does not

necessarily correspond to a conserved Noether charge.

As we will see in Sec.

III

, our basic set up will consist on

perturbing GR solutions to obtain such hairy BHs.

SCALARIZATION OF ISOLATED BLACK HOLES IN SCALAR

...

PHYS. REV. D

110,

024033 (2024)

024033-3

B. Fixing the equations of scalar Gauss-Bonnet gravity

Before using the methods of NR to obtain numerical

predictions, further mathematical treatment is often needed

to bring the evolution equations into the form of a strongly

hyperbolic PDE system

—

which under appropriate initial

conditions can then allow for local well-posedness of the

IVP. For GR, this procedure often makes use of auxiliary

variables, constraint equations and gauge freedom

—

some

widely used families of formulations include the Z4 and

generalized harmonic systems.

For the equations of full sGB theory [Eqs.

(3)

–

(7)

]

rewriting these equations in such a way is no trivial task.

All of the existing approaches (either perturbative

[42]

with modified gauges

[57,58,61,62]

) are valid only at

weak

coupling

l

, where the corrections with respect to GR

are assumed to remain suitably small.

In the fixing-the-equations approach

[63]

, we replace the

terms that alter the principal part of the system (the terms

with the highest order of derivatives) with new variables

subject to

ad hoc driver

equations. The aim is to have the

system react in a time scale

such that it will have the effect

of damping (high-frequency) modes

with

≳

which

are suspected of leading to instabilities during the evolu-

tion. As of yet, there is no unique prescription for the terms

to replace nor on the precise form of that the driver

equations need to have. In Sec.

II D

we will use precisely

this freedom to devise a driver equation adapted to the

problem of recovering stationary solutions.

For the equations of sGB [Eqs.

(3)

and

(4)

], we replace

the right-hand side of the equations with new variables

to obtain

[68]

□

In this way, the principal part of Eq.

(12)

is reduced to that

of the simpler Klein-Gordon and Einstein equations, for

which the standard methods to rewrite them as strongly

hyperbolic systems apply.

To close the system, additional evolution equations that

ensure that the auxiliary variables approach their corre-

sponding physical source terms on a specified time scale

must be prescribed for

. The simplest example for

such

driver

equations is

[65]

∂

−

;

∂

−

;

where

≡

−

l

;

≡

l

;

FIG. 1. Scalarized spinning BH. Curvature-induced scalariza-

tion for a spinning BH with spin

in the

direction (top)

evolved in the decoupling limit [Eq.

(8)

] for model

(9)

. We have

set

l

;

−

. In color, the absolute value of

the scalar field

t=M

1000

. See Sec.

for more details.

FIG. 2. Scalarized spinning BH (equatorial profile). Top: scalar

profile (solid black) for the BH of Fig.

along the equatorial

direction at

t=M

800

. In colored lines, various approximations

with the fixing-the-equations approach. Dashed lines correspond

to driver I [Eq.

(16)

] whereas solid lines correspond to our

new

driver II [Eq.

(19)

]. We fix

. The scalar profiles are

plotted starting from the excision sphere, slightly inside the

apparent horizon of the BH. Bottom: relative error in the scalar

profile with respect to the correct solution given by the decou-

pling limit [Eq.

(8)

GUILLERMO LARA

et al.

PHYS. REV. D

110,

024033 (2024)

024033-4

and

is a relaxation timescale (with units of mass)

controlling how fast the new variables approach the original

source terms. (Notice the parallel of Eq.

(13)

to exponential

decay laws or to constraint damping terms in NR.) In some

sense, if the action

(1)

is thought as an effective field theory,

the damping of high frequencies can be interpreted as a way

“

soften

”

the truncation of the theory in a way that allows

one to obtain quantitative predictions.

Finally, to compute or replace the second-order time-

derivative terms (

∂

∂

) remaining in

;

additional techniques may be required, such as

order

reduction

[65]

, where the leading order equations of motion

(in

l

) would be used to reduce the order of second-order

time-derivative terms.

Although Eq.

(13)

may be the simplest form of driver

equations, in practice, they may not be the most convenient

in an actual numerical implementation since they result in

standing modes for

, instead of propagating

modes. Indeed, in Ref.

[68]

, numerical relativity simula-

tions of the fully coupled system

(12)

have been performed

in spherical symmetry using a driver equation with wave-

like properties. We will consider other alternatives for these

driver equations in Sec.

II D

C. Case study: Fixing the decoupling limit

of scalar Gauss Bonnet

As remarked in Sec.

, in this paper, we take our first

steps towards an implementation in full 3-dimensional

space. Since such an undertaking is far more complex

than that of Ref.

[68]

, for simplicity, we begin by studying

the scalar dynamics of system

(12)

, i.e. we will

“

fix

”

the

decoupling limit of scalar Gauss-Bonnet theory. To be more

precise, we will evolve

□

;

where the auxiliary variable

is subject to a driver equation

to be specified below.

In this particular case, the original decoupling limit

equations [Eq.

(3)

and

] can be cast as a strongly

hyperbolic system and can be readily evolved without

fixing. Therefore, this set up will be useful as an exercise to

test the numerical implementation of the system as well as

to evaluate whether it can accurately describe the dynamics

of the scalar and the final stationary solutions. Indeed, we

can directly compare with the solutions of the unmodified

theory, which are, in some sense the

“

exact

”

solution which

one would like to recover with the fixing-the-equations

approach.

D. Advanced drivers and recovery

of stationary solutions

In this section we consider more complex choices for

driver equations for scalar variables. We will begin by

considering a driver equation first presented in Ref.

[69]

(with some slight modifications

)

∂

−

∂

∂

−

∂

−

;

where

is the lapse,

is the shift, and

is the

spatial metric (with inverse

) in the (

) decompo-

sition of the spacetime metric

−

Þð

;

and also where

≡

−

l

and

;

are dimen-

sionful non-negative parameters (with dimensions

mass

and [mass], respectively) that will control the relaxation

time scale(s) of

towards the desired value given by

.In

the following, we will refer to Eq.

(16)

driver I

.Justas

Eq.

(13)

is reminiscent of exponential decay laws, by

looking at the time-derivative terms, Eq.

(16)

is reminis-

cent of a forced damped harmonic oscillator,

̈

Þþ

̇

Þþ

Þ¼

,where

≡

ffiffiffi

is the oscil-

lator frequency,

≡

ffiffiffi

is the damping parameter,

and the source

enters as a forcing term

≡

Moreover, it is clear that the advective operators

∂

−

∂

entering driver I are such that the propagating modes

will move in the opposite direction of

. For BH space-

times (depending on the specific gauge) these modes can

become ingoing (into the hole) near apparent horizons,

thus allowing inaccuracies in the tracking of the source

term to propagate outside of the region of interest.

While driver I has many desirable properties, one issue

associated with it (already pointed out in Ref.

[69]

), as well

as with related wavelike drivers of the form (such as the one

of Ref.

[68]

)

−

□

∂

−

∂

−

;

is that stationary solutions are not recovered exactly. To see

this, consider a stationary situation, with respect to the

coordinate

associated with an approximate Killing vector

∂

. Assuming that all fields have become time inde-

pendent

t; x

Þ¼

, we can set to zero all the time

derivatives in driver I. The resulting equation,

∂

−

τβ

∂

−

, still contains spatial

derivative operators. If

is spatially dependent (as for

any non-trivial solution),

cannot be a solution, and

driver I will not reproduce the correct answer. Nevertheless,

one would expect that these errors would decrease as

;

are decreased

—

or as the relaxation time scales associated

The precise difference between the driver equation here and

that in Ref.

[69]

lies in the terms with the lowest derivatives (e.g.

terms proportional to first derivatives of the lapse and shift), as

well as multiplicative factors of the lapse.

SCALARIZATION OF ISOLATED BLACK HOLES IN SCALAR

...

PHYS. REV. D

110,

024033 (2024)

024033-5

to these parameters are decreased. In Sec.

we will show

with examples that this indeed seems to be the case.

To address this issue, we propose a

new

driver given by

∂

−

∂

Þð

∂

∂

∂

∂

−

In the following, we will refer to this equation as

driver II

For a BH in uniform motion, we identify

as the velocity

of the BH as obtained from the trajectory of the center of

the apparent horizon. The

new

driver equation has now the

desired limit to recover stationary solutions, which can be

seen as follows: proceeding as before, for a BH with zero

linear momentum (

), setting the time derivatives to

zero gives indeed

. Whereas for a BH with linear

momentum, the same should hold since one can identify

≡

∂

∂

as the

convective

time derivative. In Sec.

we will show that indeed Eq.

(19)

gives an improved

recovery over Eq.

(16)

III. METHODOLOGY

A. First-order reduction

All of the component subsystems of Eqs.

(8)

and

(15)

are

decoupled and can be written individually in first-order

symmetric hyperbolic form

∂

∂

S

;

where

are first-order variables,

is a square

symmetric matrix that may depend on

(but not its

derivatives), and

S

is a source term.

The principal part of the scalar equation for

is that of

Klein-Gordon equation. Therefore, we introduce first-order

variables

≡

;

≡

∂

;

≡

−

∂

;

where

≡

−

αδ

is the unit normal to the hypersurface,

and rewrite the scalar equation as

∂

−

∂

≃

;

∂

−

∂

αγ

∂

−

∂

≃

;

∂

−

∂

∂

−

∂

≃

Here,

;

are constraint damping parameters and the

≃

notation indicates that only the principal part terms are

displayed

—

see Ref.

[90]

for the full expressions.

The Einstein field equations in the generalized harmonic

gauge have a principal part of the form

∂

∂

≃

, and

can therefore be written in a first-order form analogous to

that of the Klein-Gordon equation. We refer to Lindblom

et al.

[91]

for the full first-order equations employed here,

where the first-order variables are given by

≡

;

iab

≡

∂

;

≡

−

∂

The first-order form for the driver sector depends on the

particular driver used. For Eq.

(16)

it reads

∂

−

∂

−

;

∂

−

∂

−

;

whereas for Eq.

(19)

,itisgivenby

∂

∂

−

;

∂

−

∂

−

Both systems depend on the first-order variables

≡

;

, where

is defined by the first equation

of each system. Notice that we do not need to introduce a

first-order variable

for

∂

Finally, the Gauss-Bonnet scalar (in vacuum

≡

)

appearing in the sources of the scalar equations for

;

is computed in practice as

abcd

−

;

where

≡

aibj

and

≡

−

aibj

are the

electric and magnetic parts (respectively) of the Weyl

tensor

abcd

, with left dual

abcd

≡

abef

–

their full

expressions in terms of (

)-decomposition quantities

are given in Eqs. (B7) and (B8) of Ref.

[42]

B. Evolution code

We implement systems

(8)

and

(15)

ECTRE

[64]

,an

open-source code based on task-based parallelism.

The full system is discretized following a discontinuous

Galerkin (DG) scheme employing a numerical upwind flux.

Evolution in time is carried out by means of a fourth-order

Adams-Bashforth time stepper with local adaptive time-

stepping

[92]

, and we apply a weak exponential filter on all

evolved fields at each time step. The spatial domain is

illustrated in Fig.

and consists of a series of concentric

spherical shells. The outer boundary is located at

R=M

500

, while the inner boundary conforms to the

shape of the apparent horizon. Each shell is further

decomposed into elements isomorphic to unit cubes,

endowed with a tensor product of Legendre polynomials.

The implementation of the scalar sector with a Klein-

Gordon equation has been recently been described in detail

in Ref.

[90]

. In particular, we impose constraint-preserving

spherical-radiation boundary conditions

[93]

on the scalar

field

. The constraint damping parameters

;

are

modulated by Gaussian profiles

≡

exp

−

c;J

;

GUILLERMO LARA

et al.

PHYS. REV. D

110,

024033 (2024)

024033-6

centered around the black hole

c;J

. We set

≡

and specify

through the parameters

−

For the metric sector, the full details of the numerical

implementation will bepresentedin a specialized publication

elsewhere

[94]

, and we only broadly summarize the most

relevant aspects here. The evolution system for the metric

variables in the GH gauge is given in Ref.

[91]

and is

implemented in a dual-frame formulation

[95]

.The

constraint-damping parameters

;

follow a spatial

distribution ofthe form

(27)

. We choose

≡

−

and

;

w=M

−

for

and

. The gauge functions

in the GH system are evolved in damped harmonic (DH)

gauge

[96

–

98]

. A portion of the domain inside the apparent

horizonisexcisedandnoboundaryconditionsareimposedat

the excision sphere

—

we monitor, however, that at that

boundary all characteristic fields have velocities flowing

out of the computational domain, so that no boundary

condition is required for the evolution to remain well posed.

Finally, constraint-preserving Bjorhus boundary conditions

[91,99]

are imposed at the outer boundary.

For the scalar driver part, the implementation is similar to

that of the scalar field. We monitor that the characteristic

velocities are flowing out of the computational domain at

the excision sphere and impose zero Dirichlet boundary

conditions (i.e.

out

≡

) at the outer boundary of the

spatial domain. Further, we allow for a Gaussian spatial

dependence of the form

(27)

for the parameters

;

introduced in the driver equation for

. We choose

w=M

150

for

w=M

150

−

for

and vary the amplitude

for both parameters. For

simplicity, when we quote values of

;

, we refer to their

values at the center of the Gaussian profile.

C. Initial data

We choose as initial data a

hairless

Kerr black hole

subject to an initial scalar perturbation on

. As remarked

in Sec.

, this solution is unstable and is expected to

migrate during the evolution to a stable

hairy

black hole.

For the metric variables, we choose Kerr-Schild initial

data (with

), in Cartesian coordinates

¼ð

x; y; z

Here

diag

−

;

is the Minkowski metric, and

the scalar function

and one-form

(which satisfies

∂

∇

) are given by

≡

;

≡

;

−

;

where the spin direction is along the

direction,

implicitly defined through

Þþð

, and

is the dimensionless spin parameter

of the BH.

For the case of a BH with linear momentum, we obtain a

boosted Kerr-Schild solution by applying the appropriate

Lorentz boost to the coordinates

and

To induce the scalarization of the BH, we prescribe a

scalar perturbation of the form

[100]

Þj

≡

;

Þj

≡

−

;

where

;

are spherical harmonics and

;w;r

are

the amplitude, width and radius of the scalar profile.

For the quartic model

(9)

, this initial data for

implies

that initially the scalar source vanishes, i.e.

≡

−

l

. Therefore, we prescribe zero initial

conditions for the scalar driver, i.e.

≡

D. Comparison strategy and diagnostics

For each of the examples in Sec.

we evolve three

systems of equations with the same initial data. These are

(1) The (unmodified) decoupling limit equations

(8)

We label these solutions

dec

(2) The fixed decoupling limit equations

(15)

with

evolved using driver I [Eq.

(16)

]. We label these

solutions

fix

;

(3) The fixed decoupling limit equations

(15)

with

driver II [Eq.

(19)

]. We label these solutions

fix

;

FIG. 3. Spatial domain and apparent horizon. Domain decom-

position for a black hole with (dimensionless) spin

χ

ˆ

s

where

ˆ

¼ð

ffiffiffi

;

ffiffiffi

, in the

plane depicted above. The

black lines correspond to discontinuous Galerkin element boun-

daries. The solid black region is delimited by the apparent horizon

of the BH.

SCALARIZATION OF ISOLATED BLACK HOLES IN SCALAR

...

PHYS. REV. D

110,

024033 (2024)

024033-7

We compare the scalar field volume configurations

across the different systems (decoupling limit against the

fixed versions of the decoupling limit) by computing a

relative

accuracy estimates

of the form

≡

dec

−

fix

dec

kþ

;

where

is the desired quantity to compare,

is a

suitable norm and

is a small number.

The accuracy estimate for the scalar field profile along a

direction, for instance, in the

direction, is given by

Þð

≡

dec

−

fix

Þj

dec

Þj þ

When we compare the across the whole spatial domain we

can define (see Refs.

[65,67]

for similar measures)

≡

dec

−

fix

dec

;

is the

norm with respect to the grid points

≡

t; x

;

with

being the number of grid points.

We monitor the accuracy of the auxiliary field tracking

the original value of the source term with (see Refs.

[67,69]

for similar measures)

≡

−

Finally, we compute an accuracy estimate error for the

scalar waveform, which is an example of an observable

quantity. We define the accuracy estimate for the rms

waveform

ffiffiffiffiffiffiffiffiffiffi

;

≡

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

dec

−

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

fix

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

dec

;

where

ffiffiffiffiffiffiffiffiffiffi

is extracted at 6 equally spaced spheres in the

range

R=M

∈

;

100

and for which the sphere average of

the square of the scalar amplitude is computed as

≡

1. Constraints

We keep track of the first-order constraints for the scalar

field to monitor the evolution. These are given by

≡

∂

−

;

≡

∂

−

∂

For the metric we keep track of the following constraints:

≡

;

≃

∂

;

iab

≡

∂

−

iab

;

≃

∂

;

where

and

are given in full in Eqs. (43) and (44) of

Ref.

[91]

. Here,

≡

is a contraction of the

4-dimensional Christoffel symbol

. Given the symmet-

ric hyperbolic nature of the metric evolution equations, a

symmetrizer can be constructed, and thus these constraints

can be further condensed in a constraint energy

—

given

in Eq. (53) of Ref.

[91]

IV. RESULTS

In the following, we illustrate with various examples the

ability of our implementation of the fixing-the-equations

approach to reproduce the scalar dynamics of the decou-

pling limit of sGB theory.

A. Spinning black holes

We start with an initially hairless Kerr BH with dimen-

sionless spin parameter

and with zero linear

momentum. We set the coupling parameters to

;

−

andinduce spontaneous scalarization by perturbing the

BH with an initial scalar profile [Eq.

(31)

] with parameters

−

;w=M

;

l; m

Þ¼ð

;

Þg

.InFig.

we show the final stationary hairy solution in the decoupling

limit (the correct solution we are seeking to recover) against

which we will compare the results obtained with different

fixing-the-equations schemes.We then evolve the sameset of

initial data in the fixing-the-equations approach for both

driver I and for our new driver II to obtain stationary

solutions. Since the strongest dependence is on the

parameter, for simplicity, we fix

and vary only

. We choose values of

≤

such that the driver

equations are in the underdamped regime

in the

parallel made in Sec.

II D

with the damped harmonic

oscillator.

The most straightforward comparison of the solutions is

that of the scalar profile along a radial direction at late

times. In the top panel of Fig.

we plot the scalar field

along the

axis (in the equatorial plane) for the exact

solution (solid black) and for the different

fixing-the-

equations

approximations (colored lines). In the bottom

panel, we plot the accuracy estimate for the radial profile

with respect to the exact solution. We observe that as

decreased, the solutions obtained with both driver I and II

approach the desired solution. For instance, for the largest

choice of parameters,

, the solution with driver I

(dashed red) completely fails to reproduce the scalarized

solution, giving a vanishing scalar profile everywhere.

GUILLERMO LARA

et al.

PHYS. REV. D

110,

024033 (2024)

024033-8

As sigma is decreased to

−

(dashed green), the

solution obtained with this driver is qualitatively similar in

shape to the exact solution, but underestimated in magni-

tude by about 20%. For smaller values of

, the agreement

with the exact solution is close enough that it is no longer

apparent in the top panel, and they yield an accuracy

estimate of

≲

. This is to be expected since

ffiffiffi

is a

timescale (length scale) that controls how closely the driver

variable

is tracking the intended source term

Nevertheless, it is apparent that driver II (solid lines)

reproduces the correct solution far more accurately than

driver I. Even for the largest parameters

, the

agreement with the stationary solution is better than the

best solution obtained with driver I.

Having focused on late times in the previous paragraph,

we now turn to analyzing the tracking in time throughout

all of the scalarization event. In our implementation of the

fixing-the-equations approach,

acts as the source term of

the scalar field

. To correctly capture the physics of the

system, it needs to accurately track the original value of the

source term

≡

−

l

. In order to evaluate how well

is tracking

, in Fig.

we compute the accuracy estimate

for the scalar source term [Eq.

(37)

] as a function of time.

Since

is initially vanishing, the accuracy estimate is

largest at the beginning of the simulation. Moreover, the

source term is proportional to

, which scales as

Therefore, it remains significantly small during the imme-

diate evolution. We ignore the very early peak since we

consider the relative accuracy estimate to not be meaningful

there. A second peak at

t=M

≃

100

corresponds to the time

when the tachyonic instability in the scalar is active and the

scalar hair is growing at the fastest rate. At late times, as the

scalar settles into a stationary situation, the estimate tends

to an (almost) constant value. As before, the accuracy

estimate decreases as

is decreased. Whereas driver I is

limited at late times to approach

exactly by nonvanishing

spatial derivative terms, our new driver II was designed to

avoid this issue

—

recall the discussion in Sec

II D

Therefore, the improved behavior of our new driver II

(in solid lines) is apparent throughout the entire evolution:

the accuracy estimate of

is generally smaller, decreases

faster and achieves much lower values at late times. We

have observed similar behavior when we decrease the value

with respect to

Since

and

are related through a differential equation,

the tracking the source term is an indirect comparison of the

dynamics of the scalar field. A more direct comparison

however is to compare the full 3-dimensional scalar

configuration of the fixing-the-equations approximations

against that of the correct solution. In Fig.

, we plot the

accuracy estimate in time of the spatial scalar profiles with

respect to the reference, which confirms that the tracking

behavior of the source term

observed before directly

translates into a good tracking of the dynamics of

.As

opposed to the tracking estimate in Fig.

, this estimate is

meaningful throughout the simulation.

Finally, we turn to comparing an observable quantity: the

scalar waveform. In the top panel of Fig.

, we plot the rms

scalar waveform extracted at

R=M

100

as well as the

waveform in the exact solution. For completeness, we also

FIG. 4.

dependence for the accuracy estimate of the scalar

source tracking (spinning). Accuracy estimate for the scalar

source (dashed lines) [Eq.

(37)

] for driver I [Eq.

(16)

]asa

function of time and for different values of sigma and where we

fix

. Also in the plot, the same accuracy estimate (solid

lines) for the driver II [Eq.

(19)

FIG. 5.

dependence for the accuracy estimate in scalar field

slices (spinning). Accuracy estimate for the scalar field 3-

dimensional slices (dashed lines) [Eq.

(35)

] for driver I

[Eq.

(16)

] as a function of time and for different values of sigma

and where we fix

. Also in the plot, the same accuracy

estimate (solid lines) for the driver II [Eq.

(19)

SCALARIZATION OF ISOLATED BLACK HOLES IN SCALAR

...

PHYS. REV. D

110,

024033 (2024)

024033-9