# THE BONDI PROBLEM REVISITED:
A SPECTRAL DOMAIN DECOMPOSITION CODE

M. A. ALCOFORADO  Departamento de Física Teórica - Instituto de Física A. D.
Tavares, Universidade do Estado do Rio de Janeiro, R. São Francisco Xavier,
524. Rio de Janeiro, RJ, 20550-013, Brazil<EMAIL_ADDRESS>W. O.
BARRETO  Centro de Ciências Naturais e Humanas, Universidade Federal do ABC,
Av. dos Estados 5001, 09210-580 Santo André, São Paulo, Brazil
<EMAIL_ADDRESS>
Centro de Física Fundamental, Universidad de Los Andes, Mérida 5101, Venezuela
H. P. DE OLIVEIRA  Departamento de Física Teórica - Instituto de Física A. D.
Tavares, Universidade do Estado do Rio de Janeiro, R. São Francisco Xavier,
524. Rio de Janeiro, RJ, 20550-013, Brazil
<EMAIL_ADDRESS>
###### Abstract

We present a simple domain decomposition code based on the Galerkin-
Collocation method to integrate the field equations of the Bondi problem. The
algorithm is stable, exhibits exponential convergence when considering the
Bondi formula as an error measure, and is computationally economical. We have
incorporated features of both Galerkin and Collocation methods along with the
establishment of two non-overlapping subdomains. We have further applied the
code to show the decay of the Bondi mass in the nonlinear regime and its
power-law late time decay. Another application is the determination of the
wave-forms at the future null infinity connected with distinct initial data.

##  I Introduction

The recent direct detection of gravitational waves by the LIGO consortium
LIGO_gws  opened a new window to observe the universe. Gravitational waves
have a unique feature of extracting mass from the source besides carrying
relevant information about their origin. The observed signals were identified
as produced by a collision of spinning black holes after comparing them with
those waveforms generated by successful long-term simulations of binary black
holes pretorius ; campanelli ; baker ; mroue ; husa . For this reason,
numerical relativity has become a relevant and mature area of investigation
where we can uncover the consequences of the gravitational strong field
regime.

A crucial step of understanding the emission of gravitational waves by an
isolated source was provided by the seminal paper of Bondi and collaborators
bondi . The decay of the Bondi mass as a consequence of the emission of
gravitational waves in connection with the news function is the central result
that is summarized in the Bondi formula. With Sachs sachs  and Penrose’s works
penrose , Bondi et al.’s contribution constitute the cornerstones of general
relativity’s characteristic formulation, particularly tailored to study
gravitational radiation. The review of Winicour winicour_lrr  presents a self-
contained and complete discussion on the characteristic formulation of the
field equations.

The first numerical code to evolve the Bondi field equations, and later
extended to the more general Bondi-Sachs problem, is the PITT code
papadopoulos . In a series of papers, the code showed up to be stable, second-
order accurate, and fully nonlinear, besides several applications in
situations of physical interest isaacson ; bishop ; babiuc . Another feature
of the PITT code is the calculation of the waveforms at null infinity
bishop_97 ; babiuc_2011_1 ; babiuc_2011_2 .

We have proposed the first code based on the Galerkin and Collocation methods
to evolve the Bondi equations several years ago rodrigues , later Handmer and
Szilagyi implemented a spectral code for the Bondi-Sachs equations handmer .
We highlight two features of the spectral algorithm: the low computational
cost to achieve good accuracy and combined aspects of the Galerkin and
Collocation methods. For the sake of clarifying this last aspect, we provide
below a brief description of the basic structure we have employed in the
present formulation.

Let us consider a function $f(x,t)$ that satisfies a differential equation,
say $f_{,t}-L(f,f_{,x},..)=0$ in a spatial domain $a\leq x\leq b$, where $L$
is a nonlinear function. Spectral methods belong to a general class of the
weighted residual methods (WRM) finlayson  that establish an alternative
strategy to solve any differential equation. The first step is to approximate
the function $f(x,t)$ by a finite series expansion with respect to a set of
analytical functions (for instance, Chebyshev or Legendre polynomials) known
as the trial or basis functions. Accordingly, we have

$f_{N}(x,t)=\sum_{k=0}^{N}\,\hat{f}_{k}(t)\psi_{k}(x),$ (1)

where $\hat{f}_{k}(t)$ are the unknown coefficients or modes, $N$ is the
truncation order that dictates the number of modes, and $\psi_{k}(x)$
represents the basis functions. We choose the modes $\hat{f}_{k}(t)$ such that
an error measure provided by the residual equation
$\mathrm{Res}_{N}(x,t)=f_{N,t}-L(f_{N},f_{N,x},..)$ is forced to be zero in an
average sense, meaning that the weighted integral of the residual equation are
set to zero, or

$\left<\mathrm{Res}_{N}(x,t),\phi_{j}(x)\right>=\int_{a}^{b}\,\mathrm{Res}_{N}(x,t)\phi_{j}(x)w(x)dx=0,$
(2)

for all $j=0,1,..,N$. Here $\phi_{j}(x)$ are the test functions and $w(x)$ is
the associated weight. The choice of the test functions specifies the type of
spectral method we are adopting boyd ; canuto . For instance, if
$\phi_{j}(x)=\psi_{k}(x)$, and choosing the basis functions to satisfy the
boundary conditions, we have the traditional Galerkin method. On the other
hand, if $\phi_{j}(x)=\delta(x-x_{j})$, where $x_{j},\,j=0,1,..,N$ are the
collocation points, we obtain the Collocation method. In this case, we can
infer from Eq. (2) that the residual equation vanishes at each collocation
point, $\mathrm{Res}_{N}(x_{j},t),\,j=0,1,..,N$ (here $w=1$). Consequently,
any spectral method approximates an evolution partial differential equation as
a finite set of ordinary differential equations for the modes $\hat{f}_{k}$ or
the values of $f(x,t)$ at the collocation points, $f_{j}(t)=f_{N}(x_{j},t)$ in
the case of the Collocation method. For an elliptic-type equation, spectral
methods approximate it as a finite set of algebraic equations related to the
modes.

We have combined features of the Galerkin and Collocation methods to adopt
basis functions that satisfy the boundary conditions and assume that the test
functions are Dirac functions. Further, in some cases, we evaluate the
integrals (2) with quadrature formulae, a prescription of the Galerkin method
with numerical integration (G-NI) canuto . We have coined this combination by
the Galerkin-Collocation method and applied it to several problems: the Bondi
problem rodrigues , the spherical collapse of scalar fields and critical
phenomena rodrigues_spherical ; crespo_bh_scalar_field ; crespo_kink ;
crespo_affine ; alcoforado_critical ; alcoforado_cauchy , the evolution of
cylindrical waves celestino_cylindrical ; barreto_cylindrical_DD , and the
initial data for numerical relativity matias_sd ; barreto_DD ; barreto_DD_2 .

Recently, we have incorporated the technique of domain decomposition into the
spectral Galerkin-Collocation scheme (see Refs. crespo_affine ;
alcoforado_critical ; alcoforado_cauchy ; barreto_cylindrical_DD ; barreto_DD
; barreto_DD_2 ). The domain decomposition or multidomain method consists of
dividing the spatial domain into two or more subdomains and establishing
appropriate transmission conditions to connect the solutions in each
subdomain. Domain decomposition is adequate to tackle complicated geometries
or exhibit strong field gradients in some areas of the spatial domain. The use
of this technique is not new in numerical relativity when combined with
spectral methods bona ; pfeiffer ; ansorg ; spec ; lorene ; szilagyi_09 ;
kidder1 ; kidder2 ; hemberger_13 ; sxs_col .

The present paper aims to incorporate the domain decomposition technique in
the Galerkin-Collocation code for the Bondi equations. In Section II, we
present the metric, the field equations of the Bondi problem, and we discuss
the evolution scheme provided by the characteristic scheme. We also show the
boundary conditions imposed on the metric functions to satisfy the
requirements of spacetime regularity and asymptotic flatness and the Bondi
formula. We provide a detailed description of the domain decomposition
algorithm in Section III. We have defined the basis functions in each
subdomain, the transmission conditions, and the procedure to approximate the
hypersurface and evolution equations into sets of finite algebraic and
dynamical systems, respectively. Section IV is devoted to the code validation
with the convergence tests using the Bondi formula.

Some physical applications are presented in Section V. We have exhibited the
power-law Bondi mass decay for late times. We have also considered several
examples of the decay of the Bondi mass for higher initial amplitudes and the
gravitational waveforms generated in each case. Finally, in Section VI we
conclude.

##  II The field equations

The metric established by Bondi, van der Burgh and Metzner bondi  that
describes axisymmetric and asymptotic spacetimes takes the form

$\displaystyle
ds^{2}=-\left(\frac{V}{r}\mathrm{e}^{2\beta}-U^{2}r^{2}\mathrm{e}^{2\gamma}\right)du^{2}-2\mathrm{e}^{2\beta}dudr$
$\displaystyle-2Ur^{2}\mathrm{e}^{2\gamma}dud\theta+r^{2}(\mathrm{e}^{2\gamma}d\theta^{2}+\mathrm{e}^{-2\gamma}\sin^{2}\theta
d\varphi^{2}).$ (3)

Here, $u$ is the retarded time coordinate such that $u=constant$ denotes the
outgoing null cones; the radial coordinate $r$ is chosen demanding that the
surfaces $(u,r)$ have area equal to $4\pi r^{2}$ and the angular coordinates
$(\theta,\varphi)$ are constant alogn the outgoing null geodesics. The metric
functions $\gamma,\beta,U$ and $V$ depend on the coordinates $u,r,\theta$ and
satisfy the vacuum field equations $R_{\mu\nu}=0$. Following Bondi et al bondi
, the equations are organized into two sets: three hypersurface equations and
one evolution equations, given, respectively, by:

$\displaystyle\beta_{,r}$ $\displaystyle=$
$\displaystyle\displaystyle{\frac{1}{2}r(\gamma_{,r})^{2}},$ (4)
$\displaystyle\left(r^{4}\,\mathrm{e}^{2(\gamma-\beta)}\displaystyle{U_{,r}}\right)_{,r}$
$\displaystyle=$
$\displaystyle\displaystyle{2r^{2}\,\left[r^{2}\left(\frac{\beta}{r^{2}}\right)_{,r\theta}-\frac{\left(\sin^{2}\theta\,\gamma\right)_{,r\theta}}{\sin^{2}\theta}+2\gamma_{,r}\gamma_{,\theta}\right]},$
(5) $\displaystyle V_{,r}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}r^{4}\mathrm{e}^{2(\gamma-\beta)}\left(U\right)_{,r}^{2}+\frac{\left(r^{4}\sin\theta\,U\right)_{,r\theta}}{2r^{2}\sin\theta}$
(6)
$\displaystyle+\mathrm{e}^{2(\beta-\gamma)}\bigg{[}1-\frac{(\sin\theta\,\beta_{,\theta})_{,\theta}}{\sin\theta}+\gamma_{,\theta\theta}+3\cot\theta\,\gamma_{,\theta}$
$\displaystyle-(\beta_{,\theta})^{2}-2\gamma_{,\theta}(\gamma_{,\theta}-\beta_{,\theta})\bigg{]},$
$\displaystyle 4r(r\gamma)_{,ur}$ $\displaystyle=$
$\displaystyle\Bigg{\\{}2r\gamma_{,r}V-r^{2}\left[2\gamma_{,\theta}U+\sin\theta\,\left(\frac{U}{\sin\theta}\right)_{,\theta}\right]\Bigg{\\}}_{,r}$
(7)
$\displaystyle-2r^{2}\frac{(\gamma_{,r}U\sin\theta)_{,\theta}}{\sin\theta}+\frac{1}{2}r^{4}\mathrm{e}^{2(\gamma-\beta)}(U_{,r})^{2}$
$\displaystyle+2\mathrm{e}^{2(\beta-\gamma)}\bigg{[}(\beta_{,\theta})^{2}+\sin\theta\left(\frac{\beta_{,\theta}}{\sin\theta}\right)_{,\theta}\bigg{]}.$

The subscripts $u,r$ and $\theta$ denote derivatives with respect to these
coordinates. In the seminal work of Bondi et al.bondi , the characteristic
formulation of General Relativity was introduced and studied for the first
time. As such, the evolution scheme obeys a nice hierarchical structure: once
the initial data $\gamma(u_{0},r,\theta)=\gamma_{0}(r,\theta)$ is fixed, the
hypersurface equations (4) - (6) determine the metric functions $\beta,U$ and
$V$ (modulo integration constants) on the initial null cone $u=u_{0}$. Taking
into account these results, the evolution equation (7) provides $\gamma_{,u}$
at the initial null cone $u=u_{0}$, and consequently allows the determination
of $\gamma$ at the next null surface $u=u_{0}+\delta u$. By repeating the
whole cycle, we ended up with the evolution of spacetime.

By inspecting the field equations, the conditions of regularity of the
spacetime at the origin $r=0$ are

$\displaystyle\gamma=\mathcal{O}(r^{2}),\;\beta=\mathcal{O}(r^{4}),\;U=\mathcal{O}(r),\;V=r+\mathcal{O}(r^{3}),$
(8)

and the conditions of smootheness on the symmetry axis ($\theta=0,\pi$) demand
that

$\displaystyle\bar{\gamma}\equiv\frac{\gamma}{\sin^{2}\theta}\;\;\mathrm{and}\;\;\bar{U}\equiv\frac{U}{\sin\theta},$
(9)

are continuous functions at $\theta=0,\pi$. Taking into consideration Eq. (4)
and the above relation, we have

$\displaystyle\bar{\beta}\equiv\frac{\beta}{\sin^{4}\theta}.$ (10)

Then, when trackling the field equations numerically, we will consider the
metric functions $\bar{\gamma}$, $\bar{U}$, $\bar{\beta}$ and $S$ defined in
Ref. papadopoulos

$\displaystyle S\equiv\frac{V-r}{r^{2}}.$ (11)

In this case, the function $S$ satisfies the same condition of the metric
function $\bar{U}$ with respect to the radial dependence, that is
$S=\mathcal{O}(r)$.

Another relevant piece of information regards the asymptotic form of the
metric functions. We demand the spacetime to be asymptotically flat, and
assuming, for the sake of convenience, the Winicour-Tamburino frame
winicour_tamburino ; isaacson  that consists in reproducing the Minkowski
metric at the center of symmetry, the following conditions must be fulfilled:

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle
K(u,\theta)+\displaystyle{\frac{c(u,\theta)}{r}}+\mathcal{O}(r^{-2}),$ (12)
$\displaystyle\beta$ $\displaystyle=$ $\displaystyle
H(u,\theta)+\mathcal{O}(r^{-2}),$ (13) $\displaystyle U$ $\displaystyle=$
$\displaystyle L(u,\theta)+\mathcal{O}(r^{-1}),$ (14) $\displaystyle V$
$\displaystyle=$
$\displaystyle\frac{(L\sin\theta)_{,\theta}}{\sin\theta}r^{2}+r\mathrm{e}^{2(H-K)}V_{1}(u,\theta)-2\mathrm{e}^{2H}M(u,\theta)$
(15) $\displaystyle+\mathcal{O}(r^{-1}),$

where $V_{1}(u,\theta)$ is related to the functions $H,K$ and $L$;
$M(u,\theta)$ is the Bondi mass aspect bondi ; papadopoulos . It is worth
mentioning that the above asymptotic expansions do not belong to the standard
Bondi frame, where the functions $K,H,L$ vanish bondi . Regarding the function
$S$ its behavior is similar to the function $U$.

The Bondi formula is the main result presented in Ref. bondi . It relates the
loss of mass of a localized mass distribution due to the gravitational wave
extraction. In the present frame or gauge, the Bondi formula reads as

$\displaystyle\frac{dM_{B}}{du}=-\frac{1}{2}\,\int_{0}^{\pi}\,\frac{\mathrm{e}^{2H}}{\omega}\,\mathcal{N}^{2}\sin\theta\,d\theta,$
(16)

where $M_{B}(u)$ is the Bondi mass, $\mathcal{N}(u,\theta)$ is the news
function and $\omega(u,\theta)$ is a function belonging on the gauge we are
considering. We present in the Appendix the corresponding expressions for the
Bondi mass and the news function.

##  III The domain decomposition Galerkin-Collocation method

###  III.1 Preliminary definitions: new variables and mapps

Before presenting the description of the domain decomposition GC method, it
will be convenient to introduce an additional function
$\bar{Q}=\bar{Q}(u,r,\theta)$ by

$\displaystyle\bar{Q}=\mathrm{e}^{2(\gamma-\beta)}r^{2}\bar{U}_{,r}.$ (17)

As a consequence, the second-order hypersurface equation (5) is split into two
first-order equations: the above equation for the metric function $\bar{U}$
and the resulting equation for the function $\bar{Q}$ after substituting Eq.
(17) into Eq. (5). After a simple inspection, one can shown that the auxiliary
function $\bar{Q}$ satisfies the same boundaries conditions of the metric
function $\bar{\gamma}$.

We divide the spatial domain $\mathcal{D}$ covered by $0\leq
r<\infty,\,0\leq\theta\leq\pi$ into two subdomains denoted by
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$:

$\displaystyle\mathcal{D}_{1}:\,0\leq r\leq r^{(1)},\,0\leq\theta\leq\pi,$
$\displaystyle\mathcal{D}_{2}:\,r^{(1)}\leq r\leq\infty,\,0\leq\theta\leq\pi,$

where $r=r^{(1)}$ is the interface between both subdomains. For the sake of
convenience, we introduce a new variable $x$ by

$x=\cos\theta,$ (18)

meaning that the angular patch is covered by $-1\leq x\leq 1$.

We have followed a similar strategy presented in the spherical case
alcoforado_critical ; alcoforado_cauchy ; namely, we first introduce an
intermediary computational variable $y$ using the algebraic map boyd

$\displaystyle r=L_{0}\frac{(1+y)}{1-y},$ (19)

where $L_{0}$ is the map parameter. Then, we have that $0\leq r<\infty$ is
mapped out into $-1\leq y\leq 1$. And in second place, we introduce the
following linear transformations to define the variables $-1\leq\xi^{(l)}\leq
1$, $l=1,2$ that cover the radial sector of the subdomains:

$\displaystyle
y(\xi^{(l)})=\frac{1}{2}\,\left[\left(y^{(l)}-y^{(l-1)}\right)\xi^{(l)}+y^{(l)}+y^{(l-1)}\right],$
(20)

where $y^{(0)}=-1$ and $y^{(2)}=1$. In terms of the new variables, we have
$\mathcal{D}_{l}:\,-1\leq\xi^{(l)}\leq 1$ and $-1\leq x\leq 1$ with $l=1,2$.
Fig. 1 illustrates the present scheme of domain decomposition. We define the
collocation points $(\xi^{(l)}_{k},x_{j})$ using these variables and then
mapped back to $(r_{k},\theta_{j})$ in the physical domain $\mathcal{D}$.

It will be useful for the definition of the rational Chebyshev polynomials
boyd  the relation between the radial coordinate $r$ and the computational
variables $\xi^{(l)}$ for $l=1,2$. By combining Eqs. (19) and (20), we obtain

$\displaystyle\xi^{(l)}=\frac{a^{(l)}r+b^{(l)}}{(r+L_{0})},$ (21)

where $l=1,2$ and

$\displaystyle a^{(l)}$ $\displaystyle=$
$\displaystyle\frac{2L_{0}+r^{(l)}+r^{(l-1)}}{r^{(l)}-r^{(l-1)}}$ (22)
$\displaystyle b^{(l)}$ $\displaystyle=$
$\displaystyle-\frac{2r^{(l)}r^{{(l-1)}}+L_{0}\left(r^{(l)}+r^{(l-1)}\right)}{r^{(l)}-r^{(l-1)}}.$
(23)

Here $r^{(0)}=0$ and $r^{(2)}$ is located at infinity. Notice that
$r^{(l-1)}\leq r\leq r^{(l)}$ corresponds to $-1\leq\xi^{(l)}\leq 1$ for
$l=1,2$.

Figure 1: We illustrate in the first panel the subdomains
$\mathcal{D}_{1}:\,0\leq r\leq r^{(1)},\,0\leq\theta\leq\pi$ and
$\mathcal{D}_{2}:\,r^{(1)}\leq r<\infty,\,0\leq\theta\leq\pi$. The radial
domain $0\leq r<\infty$ is mapped out into $-1\leq y\leq 1$ using the map
(19), and the angular domain is covered by $-1\leq x\leq 1$ with
$x=\cos\theta$. In the lower panels, we have introduced the computational
variables by $-1\leq\xi_{j}\leq 1$, $j=1,2$, corresponding to the compactified
radial parts of both subdomains, $-1\leq y\leq y^{(1)}$ and $y^{(1)}\leq y\leq
1$, where $y=y^{(1)}$, $-1\leq x\leq 1$ represents the interface.

###  III.2 Spectral approximations

We establish the following spectral approximations of the metric functions
$\bar{\gamma}(u,r,\theta),\,\bar{\beta}(u,r,\theta),\,\bar{U}(u,r,\theta),\,\bar{Q}(u,r,\theta)$
and $S(u,r,\theta)$ in each subdomain:

$\displaystyle\bar{\gamma}^{(l)}(u,r,\theta)$ $\displaystyle=$
$\displaystyle\sum^{N_{\gamma}^{(l)}}_{k=0}\sum^{\tilde{N}_{\gamma}^{(l)}}_{j=0}\,a^{(l)}_{kj}(u)\psi^{(l)}_{k}(r)P_{j}(\cos\theta),$
(24) $\displaystyle\bar{\beta}^{(l)}(u,r,\theta)$ $\displaystyle=$
$\displaystyle\sum^{N_{\beta}^{(l)}}_{k=0}\sum^{\tilde{N}_{\beta}^{(l)}}_{j=0}\,c^{(l)}_{kj}(u)\Phi^{(l)}_{k}(r)P_{j}(\cos\theta),$
(25) $\displaystyle\bar{U}^{(l)}(u,r,\theta)$ $\displaystyle=$
$\displaystyle\sum^{N_{U}^{(l)}}_{k=0}\sum^{\tilde{N}_{U}^{(l)}}_{j=0}\,b^{(l)}_{kj}(u)\chi^{(l)}_{k}(r)P_{j}(\cos\theta),$
(26) $\displaystyle\bar{Q}^{(l)}(u,r,\theta)$ $\displaystyle=$
$\displaystyle\sum^{N_{Q}^{(l)}}_{k=0}\sum^{\tilde{N}_{Q}^{(l)}}_{j=0}\,q^{(l)}_{kj}(u)\psi^{(l)}_{k}(r)P_{j}(\cos\theta),$
(27) $\displaystyle S^{(l)}(u,r,\theta)$ $\displaystyle=$
$\displaystyle\sum^{N_{S}^{(l)}}_{k=0}\sum^{\tilde{N}_{S}^{(l)}}_{j=0}\,s^{(l)}_{kj}(u)\chi^{(l)}_{k}(r)P_{j}(\cos\theta),$
(28)

where $l=1,2$ indicates the specific subdomain;
$\\{a^{(l)}_{kj}(u),\,b^{(l)}_{kj}(u),\,c^{(l)}_{kj}(u),\,q^{(l)}_{kj}(u),\,s^{(l)}_{kj}(u)\\}$
are the modes or the unknown coefficients. The number of modes in each
subdomain depends on the radial and angular truncation orders denoted
generically by $N^{(l)}$ and $\tilde{N}^{(l)}$, respectively. The angular
basis functions common in each subdomain are the Legendre polynomials,
$P_{j}(\cos\theta)$, whereas the radial basis functions
$\\{\psi^{(l)}_{k}(r),\,\Phi^{(l)}_{k}(r),\,\chi^{(l)}_{k}(r)\\}$ are
constructed to satisfy the conditions at $r=0$ and $r=\infty$. This is
possible if we combine appropriately the rational Chebyshev polynomials
defined in each subdomain as shown in the sequence.

###  III.3 Radial basis functions

The rational Chebyshev polynomials are defined in each subdomains as

$\displaystyle
TL_{k}^{(l)}=T_{k}\left(\xi^{(l)}=\frac{a^{(l)}r+b^{(l)}}{(r+L_{0})}\right),$
(29)

where $T_{k}(\xi)$ is the Chebyshev polynomial of kth order; the parameters
$a^{(l)}$ and $b^{(l)}$ are given by the relations (22) and (23),
respectively.

We can now define the radial basis functions for the spectral approximations
(24) - (28). As we have mentioned, the radial basis functions must satisfy the
conditions near the origin and infinity ($\mathcal{J}^{+}$) provided by the
relations (8) and (12) - (15).

We start with the common radial basis for the functions $\bar{U}$ and $S$ in
the first subdomain $\mathcal{D}_{1}$. Since near the origin
$\bar{U}=\mathcal{O}(r)$, it follows that

$\displaystyle\chi_{k}^{(1)}=\frac{1}{2}\,\left(TL_{k+1}^{(1)}(r)+TL_{k}^{(1)}(r)\right).$
(30)

The radial basis functions for $\bar{\gamma}$ and $\bar{Q}$ in the first
subdomains is a simple linear combination of the aboves basis functions:

$\displaystyle\psi_{k}^{(1)}=\frac{1}{4}\,\left[\frac{(2k+1)}{(2k+3)}\,\chi_{k+1}^{(1)}(r)+\chi_{k}^{(1)}(r)\right],$
(31)

that satisfy $\psi_{k}^{(1)}(r)=\mathcal{O}(r^{2})$ near the origin for all
$k$. The radial basis function for $\bar{\beta}$ is an elaborated combination
of the above basis. The basis functions $\Phi_{k}^{(1)}(r)$ are defined by

$\displaystyle\Phi_{k}^{(1)}$
$\displaystyle=\frac{1}{2}\,\frac{(k+1)(2k+1)(2k+3)}{(2k+5)^{2}(k+3)}\,\psi_{k+2}^{(1)}(r)$
(32)
$\displaystyle+2\frac{(k+1)(2k+1)}{(2k+5(2k+3))}\psi_{k+1}^{(1)}(r)+\frac{1}{2}\psi_{k}^{(1)}(r).$

In this case, we have $\Phi_{k}^{(1)}(r)=\mathcal{O}(r^{4})$ for all $k$ as
required by (8).

In the second subdomain, we choose the rational Chebyshev polynomials defined
in $\mathcal{D}_{2}$ as the radial basis functions for all fields. We
summarize these results in Table 1.

Table 1

| $\mathcal{D}_{1}$ | $\mathcal{D}_{2}$
---|---|---
$\bar{\gamma},\,\bar{Q}$ | $\phi_{k}^{(1)}(r)$ |
$\bar{\beta}$ | $\Phi_{k}^{(1)}(r)$ | $TL_{k}^{(2)}(r)$
$\bar{U},\,S$ | $\chi_{k}^{(1)}(r)$ |

Table 1: Radial basis functions for the spectral approximations (24) - (28).

###  III.4 Transmission conditions

The next step is to establish the transmission conditions to guarantee that
all pieces of the spectral approximations (24) - (28) represent the same
corresponding functions, i. e. $\bar{\gamma},\bar{\beta},\bar{U},\bar{Q}$ and
$S$ defined in each subdomain. As in the spherical case alcoforado_critical ;
alcoforado_cauchy , we adopt the patching method canuto  that demands that a
function and all their $d-1$ spatial derivatives must be continuous at the
contiguous subdomains’ interface. Since the interface between both subdomains
is defined by $r-r^{(1)}=0,\,0\leq\theta\leq\pi$, the spatial derivatives
involves only the radial coordinate $r$. After inspecting the field equations
(4) - (7), we come up with the following transmission conditions:

$\displaystyle\bar{\beta}^{(1)}(u,r^{(1)},\theta)$ $\displaystyle=$
$\displaystyle\bar{\beta}^{(2)}(u,r^{(1)},\theta),$ (33)
$\displaystyle\bar{U}^{(1)}(u,r^{(1)},\theta)$ $\displaystyle=$
$\displaystyle\bar{U}^{(2)}(u,r^{(1)},\theta),$ (34)
$\displaystyle\bar{Q}^{(1)}(u,r^{(1)},\theta)$ $\displaystyle=$
$\displaystyle\bar{Q}^{(2)}(u,r^{(1)},\theta),$ (35) $\displaystyle
S^{(1)}(r,r^{(1)},\theta)$ $\displaystyle=$ $\displaystyle
S^{(2)}(r,r^{(1)},\theta),$ (36)
$\displaystyle\bar{\gamma}^{(1)}(u,r^{(1)},\theta)$ $\displaystyle=$
$\displaystyle\bar{\gamma}^{(2)}(u,r^{(1)},\theta),$ (37)
$\displaystyle\left(\frac{\partial\bar{\gamma}^{(1)}}{\partial
r}\right)_{r^{(1)}}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial\bar{\gamma}^{(2)}}{\partial
r}\right)_{r^{(1)}},$ (38)

that are valid for $0\leq\theta\leq\pi$. In the numerical implementation, the
above conditions are to be taken approximately.

###  III.5 Implementing the GC domaind decomposition method

The field equations under consideration are formed by four hypersurface
equations for the functions $\bar{\beta},\bar{Q},\bar{U}$ and $S$ and an wave
equation for $\bar{\gamma}$. In this case, any spectral method will
approximate the hypersurface equations into sets of algebraic equations for
the corresponding modes and a set of ordinary differential equations for the
modes $a_{kj}^{(l)}(u)$,
$k=0,1..,N_{\gamma}^{(l)},j=0,1..,\tilde{N}_{\gamma}^{(l)}$ with $l=1,2$. In
the sequence, we summarize the description of the procedure that follows
closely our previous references alcoforado_critical ; alcoforado_cauchy ;
barreto_DD ; crespo_affine .

####  III.5.1 Hypersurface equations for $\beta$, $Q$ and $U$.

Let us consider the hypersurface equation for $\bar{\beta}$, Eq. (4). By
substituting the spectral approximations for $\bar{\gamma}$ and $\bar{\beta}$
(Eqs. (24) and (25), respectively), we obtain the corresponding residual
equation

$\displaystyle\mathrm{Res}_{\beta}^{(l)}(u,r,\theta)=\bar{\beta}_{,r}^{(l)}-\displaystyle{\frac{1}{2}r\left(\bar{\gamma}_{,r}^{(l)}\right)^{2}},$
(39)

where $l=1,2$. These equations do not vanish due to the introduced spectral
approximations. Next, to obtain a set of equations for the modes
$c_{mn}^{(l)}(u)$, $m=0,1,..,N_{\beta}^{(l)},n=0,1,..,\tilde{N}_{\beta}^{(l)}$
and $l=1,2$, we follow the prescription of the Collocation method that
consists in vanishing the residual equations (39) at the collocation points
$(r^{(l)}_{k},\theta_{j})$ in each subdomain. Equivalently, the test functions
are deltas of Dirac, $\delta(r-r^{(l)}_{k})\delta(\theta-\theta_{j})$. Then,
we have

$\displaystyle(\bar{\beta}_{,r}^{(l)})_{kj}-\frac{1}{2}r^{(l)}_{k}\,(\bar{\gamma}_{,r}^{(l)})_{kj}^{2}=0,\;l=1,2.$
(40)

The quantities where $(\bar{\beta}_{,r}^{(l)})_{kj}$ and
$(\bar{\gamma}_{,r}^{(l)})_{kj}$ are the values of the derivatives of
$\bar{\beta}^{(l)}$ and $\bar{\gamma}^{(l)}$ at the collocation points:

$\displaystyle(\bar{\beta}_{,r}^{(l)})_{kj}=\sum_{m,n}\,c_{mn}^{(l)}(u)\left(\frac{\partial\Phi_{m}^{(l)}}{\partial
r}\right)_{r^{(l)}_{k}}\,P_{n}(\cos\theta_{j}),$ (41)
$\displaystyle(\bar{\gamma}_{,r}^{(l)})_{kj}=\sum_{m,n}\,a_{mn}^{(l)}(u)\left(\frac{\partial\psi_{m}^{(l)}}{\partial
r}\right)_{r^{(l)}_{k}}\,P_{n}(\cos\theta_{j}).$ (42)

We have to establish the collocation points in each subdomain such that their
number together with the transmission conditions must be equal to

$(N_{\beta}^{(1)}+1)(\tilde{N}_{\beta}^{(1)}+1)+(N_{\beta}^{(2)}+1)(\tilde{N}_{\beta}^{(2)}+1),$

that is the total number of modes $c_{mn}^{(l)}$.

To obtain the approximate equations from the transmission condition (33), we
first need to write the corresponding residual equation

$\displaystyle\mathrm{Res}_{\beta}^{(TC)}(u,\theta)=\bar{\beta}^{(1)}(u,r^{(1)},\theta)-\bar{\beta}^{(2)}(u,r^{(1)},\theta).$
(43)

We can make the above equation to vanish at the angular collocation points
$\theta=\theta_{j}$, or alternatively to impose that
$\mathrm{Res}_{\beta}^{(TC)}(u,\theta)$ vanish in an average sense, meaning
that

$\displaystyle\left<\mathrm{Res}_{\beta}^{(TC)}(u,\theta),P_{j}(\cos\theta)\right>=\int_{-1}^{1}\,\mathrm{Res}_{\beta}^{(TC)}(u,x)P_{j}(x)dx=0,$
(44)

for all $j=0,1,..,\tilde{N}_{\beta}^{(2)}$. Therefore, we have
$\tilde{N}_{\beta}^{(2)}+1$ equations related to the transmission condition
(33) demanding

$(N_{\beta}^{(1)}+1)(\tilde{N}_{\beta}^{(1)}+1)+N_{\beta}^{(2)}(\tilde{N}_{\beta}^{(2)}+1),$

collocation points $(r_{k}^{(l)},x_{j}),\,l=1,2$. Based on the above
considerations, we choose the following collocation points in the first
subdomain:

$\displaystyle\mathcal{D}_{1}:\begin{cases}\xi^{(1)}_{k}=\displaystyle{\cos\left(\frac{k\pi}{N_{\beta}^{(1)}+2}\right)},\;k=1,2,..,N_{\beta}^{(1)}+1,\\\
\\\ x_{j}=-1,\,\mbox{zeros
of}\,\,\displaystyle{\frac{dP_{\tilde{N}_{\beta}^{(1)}}}{dx}},1,\;j=0,1,..,\tilde{N}_{\beta}^{(1)},\end{cases}$
(45)

where $\xi_{0}^{(1)}=1$ and $\xi_{N_{\beta}^{(1)}+2}^{(1)}=-1$ correspond to
the interface $r=r^{(1)}$ and the origin, $r=0$, respectively, that are not
included in the first subdomain. For the second subdomain, we have

$\displaystyle\mathcal{D}_{2}:\begin{cases}\xi^{(2)}_{k}=\displaystyle{\cos\left(\frac{k\pi}{N_{\beta}^{(2)}}\right)},\;k=1,2,..,N_{\beta}^{(2)},\\\
\\\ x_{j}=-1,\,\mbox{zeros
of}\,\,\displaystyle{\frac{dP_{\tilde{N}_{\beta}^{(2)}}}{dx}},1,\;j=0,1,..,\tilde{N}_{\beta}^{(2)}.\end{cases}$
(46)

Here $\xi_{0}^{(2)}=1$ corresponds to the infinity $\mathcal{J}^{+}$ and is
excluded since the residual equation vanish identically;
$\xi^{(2)}_{N_{\beta}^{(2)}}=-1$ is the interface $r=r^{(1)}$.

Now, we can summarize the determination or the update to the modes
$c_{mn}^{(l)}$ from the residual equations (40) together with the transmission
conditions (44). It is convenient to substitute the two indexes that
characterize the values $(\bar{\beta}^{(l)}_{,r})_{kj}$ as well as the modes
$c^{(l)}_{mn}$ by one index, or

$(\bar{\beta}^{(l)}_{,r})_{kj}\,\rightarrow\,(\bar{\beta}_{,r})_{i},\quad
c^{(l)}_{mn}\,\rightarrow\,c_{i}.$

For instance,
$(\bar{\beta}^{(1)}_{,r})_{11}=(\bar{\beta}_{,r})_{1},\,(\bar{\beta}^{(1)}_{,r})_{12}=(\bar{\beta}_{,r})_{2},...$
and so on. The relation between the values $(\bar{\beta}_{,r})_{i}$ and the
modes $c_{i}$ is written in a matrix form

$\displaystyle\begin{pmatrix}(\bar{\beta}_{,r})_{1}\\\
(\bar{\beta}_{,r})_{2}\\\ \vdots\\\ (\bar{\beta}_{,r})_{K_{\beta}}\\\ 0\\\
\vdots\\\ 0\end{pmatrix}=\mathbb{B}\,\begin{pmatrix}c_{1}\\\ c_{2}\\\
\vdots\\\ c_{K_{\beta}}\\\ c_{K_{\beta}+1}\\\ \vdots\\\
c_{N_{\beta}}\end{pmatrix},$ (47)

where the elements of the matrix $\mathbb{B}_{pq}$ are given by

$\left(\frac{d\Phi_{m}^{(l)}}{dr}\right)_{r_{k}}P_{n}(x_{j}),$

with the convention $kj\rightarrow p$ and $mn\rightarrow q$ (including
$l=1,2$). Still,
$K_{\beta}=(N_{\beta}^{(1)}+1)(\tilde{N}_{\beta}^{(1)}+1)+N_{\beta}^{(2)}(\tilde{N}_{\beta}^{(2)}+1)$
is the number of collocation points and
$N_{\beta}=K_{\beta}+\tilde{N}_{\beta}^{(2)}+1$ is the number of modes
$c_{mn}^{(l)},\,l=1,2$. In the above matricial equation, the zeros on the
first column vector correspond the inclusion of the transmission conditions.
Thus, the determination of modes $c_{i}$ is obtained after inverting the
matrix $\mathbb{B}$ with the values $(\bar{\beta}_{,r})_{j}$ determined from
the relations (40).

We have implemented a similar procedure with respect to the equations for
$\bar{Q}$ and $\bar{U}$. The only difference is the number of matrices
relating the values of
$\bar{\beta}_{kj},(\bar{\beta}_{,x})_{kj},(\bar{\beta}_{,xr})_{kj},\bar{\gamma}_{kj},(\bar{\gamma}_{,r})_{kj},(\bar{\gamma}_{,xr})_{kj}$
and $(\bar{\gamma}_{,xr})_{kj}$ (we have dropped the index that indicates the
subdomain for simplicity). These matrices and the hypersurface equations for
$\bar{Q}$ and $\bar{U}$ determine the values $(\bar{Q}_{,r})_{kj}$ and
$(\bar{U}_{,r})_{kj}$ that together with the corresponding transmission
conditions, we can update the modes $q_{kj},\,b_{kj}$ in the same way as we
have implemented previously for the modes $c_{kj}$. For the sake of
simplicity, have considered the same collocation points with coordinates (45)
and (46), in other words we are assuming that
$N_{\beta}^{(l)}=N_{Q}^{(l)}=N_{U}^{(l)}$ and
$\tilde{N}_{\beta}^{(l)}=\tilde{N}_{Q}^{(l)}=\tilde{N}_{U}^{(l)}$, $l=1,2$.

####  III.5.2 Hypersurface equation for $S$ and the evolution equation

In the sequence, we have considered the hypersurface equation for $S$ and the
evolution equation, Eqs. (6) and (7), respectively. The common aspect shared
in treating both equations spectrally is a new set of collocation points
determined by the radial and angular truncation orders of $S$, $N_{S}^{(l)}$
and $\tilde{N}_{S}^{(l)}$, respectively, where, in general, we have
$N_{S}^{(l)}\geq N_{\beta}^{(l)}$ and
$\tilde{N}_{S}^{(l)}\geq\tilde{N}_{\beta}^{(l)}$, $l=1,2$.

The new set of grid points defined in the first and second subdomains is

$\displaystyle\mathcal{D}_{1}:\begin{cases}\xi^{(1)}_{k}=\displaystyle{\cos\left(\frac{k\pi}{N_{S}^{(1)}+2}\right)},\;k=1,2,..,N_{S}^{(1)}+1,\\\
\\\ x_{j}=-1,\,\mbox{zeros
of}\,\,\displaystyle{\frac{dP_{\tilde{N}_{S}^{(1)}}}{dx}},1,\;j=0,1,..,\tilde{N}_{S}^{(1)}.\end{cases}$
(48)
$\displaystyle\mathcal{D}_{2}:\begin{cases}\xi^{(2)}_{k}=\displaystyle{\cos\left(\frac{k\pi}{N_{S}^{(2)}-1}\right)},\;k=0,1,..,N_{S}^{(2)}-1,\\\
\\\ x_{j}=-1,\,\mbox{zeros
of}\,\,\displaystyle{\frac{dP_{\tilde{N}_{S}^{(2)}}}{dx}},1,\;j=0,1,..,\tilde{N}_{S}^{(2)}.\end{cases}$
(49)

As in the first grid, the radial points of $\mathcal{D}_{1}$ do not include
the origin $r=0$ and the interface $r=r^{(1)}$. On the other hand,
$\xi_{0}^{(2)}=1$ corresponds to $r=\infty$ and $\xi_{N_{S}^{(2)}-1}^{(2)}=-1$
corresponds to the interface $r=r^{(1)}$, both radial points belong to the
second subdomain. Therefore, we have a total of

$(N_{S}^{(1)}+1)(\tilde{N}_{S}^{(1)}+1)+N_{S}^{(2)}(\tilde{N}_{S}^{(2)}+1)$

collocation points and $\tilde{N}_{S}^{(2)}+1$ relations arising from the
transmission conditions for the function $S$:

$\displaystyle\left<\mathrm{Res}_{S}^{(TC)}(u,\theta),P_{j}(\cos\theta)\right>=\int_{-1}^{1}\,\mathrm{Res}_{S}^{(TC)}(u,x)P_{j}(x)dx=0,$
(50)

where $j=0,1,..,\tilde{N}_{S}^{(2)}$ and
$\mathrm{Res}_{S}^{(TC)}(u,\theta)=S^{(1)}(u,r^{(1)},\theta)-S^{(2)}(u,r^{(1)},\theta)$.

To update the modes $s_{mn}^{(l)}$, $l=1,2$, we proceed in the same way as
outlined for the modes $c_{mn}^{(l)}$. We first construct the matrices
relating the values of
$\bar{\gamma}_{kj},(\bar{\gamma}_{,x})_{kj},(\bar{\gamma}_{,xx})_{kj},\bar{\beta}_{kj},(\bar{\beta}_{,x})_{kj},(\bar{\beta}_{,xx})_{kj},\bar{U}_{kj},(\bar{U}_{,r})_{kj}$
and $(\bar{U}_{,xr})_{kj}$ at the second set of collocation points. These
values allow the determination of the values $(S_{,r})_{kj}$ at the same
collocation points from the residual equation associated to Eq. (6). Together
with the transmission conditions (50), we can write

$\displaystyle\begin{pmatrix}(S_{,r})_{1}\\\ (S_{,r})_{2}\\\ \vdots\\\
(S_{,r})_{K_{S}}\\\ 0\\\ \vdots\\\
0\end{pmatrix}=\mathbb{S}\,\begin{pmatrix}s_{1}\\\ s_{2}\\\ \vdots\\\
s_{K_{S}}\\\ s_{K_{S}+1}\\\ \vdots\\\ s_{N_{S}}\end{pmatrix},$ (51)

where
$K_{S}=(N_{S}^{(1)}+1)(\tilde{N}_{S}^{(1)}+1)+N_{S}^{(2)}(\tilde{N}_{S}^{(2)}+1)$
and $N_{S}=K_{S}+\tilde{N}_{S}^{(2)}+1$ is the total number of modes
$s_{mn}^{(l)}$ or $s_{i}$ after transforming two indexes into one. The matrix
components $\mathbb{S}_{pq}$ are determined after following similar convention
we have presented in Eq. (47).

The last step is to consider the evolution equation (7). We established the
corresponding residual equation, $\mathrm{Res}_{\gamma}^{(l)}(u,r,x)$, at each
subdomain after substituting the spectral approximations (24) - (28) into Eq.
(7) yielding

$\displaystyle\mathrm{Res}_{\gamma}^{(l)}(u,r,x)$ $\displaystyle=$
$\displaystyle(r\bar{\gamma}_{,u})_{,r}$ (52)
$\displaystyle-\frac{1}{4r}\,\bigg{[}2r\bar{\gamma}_{,r}(r+r^{2}S)-2(1-x^{2})r^{2}\bar{\gamma}_{,x}+r^{2}\,(4x\bar{\gamma}\bar{U}+\bar{U}_{,x})\bigg{]}_{,r}$
$\displaystyle+2xr\bar{\gamma}_{,r}\bar{U}-\frac{1}{2}(1-x^{2})r(\bar{\gamma}_{,r}\bar{U})_{,x}-\frac{r^{3}}{8}\mathrm{e}^{2(\gamma-\beta)}\bar{U}_{,r}^{2}$
$\displaystyle-\frac{1}{2r}\mathrm{e}^{2(\beta-\gamma)}\bigg{\\{}(1-x^{2})^{2}\big{[}4x\bar{\beta}-(1-x^{2})\bar{\beta}_{,x}\big{]}^{2}$
$\displaystyle+\left[(1-x^{2})^{2}\bar{\beta}_{,x}-4x(1-x^{2})\bar{\beta}\right]_{,x}\bigg{\\}},$

where we have divided Eq. (7) by $4r\sin^{2}\theta$ and introduced explicitely
the variable $x$. For simplicity, we have dropped the index that indicates the
subdomain on the RHS of the above equation.

At this point, we can follow the procedure adopted for all residual equations,
i. e. imposing the vanishing of the residual equation (52) at a convenient set
of collocation points. Instead, we propose to extend the strategy of Ref.
rodrigues  in the case of a single domain. The procedure consists in to choose
the modes $a_{mn,u}^{(l)}$ such that the residual equation (52) is forced to
vanish in an average sense according to the Residual Weighted Methods
finlayson . It means that the test functions are no longer Dirac functions but
the same as the basis functions for $\bar{\gamma}^{(l)}$ (cf. Eq. (2)). Then,
we evaluate the inner products using quadrature formulae as indicated by the
G-NI method taking advantage of the second set of collocation points (48) and
(49) as

$\displaystyle\displaystyle{\left<\mathrm{Res}_{\gamma}(u,r(\xi),x),\psi_{m}(y)P_{n}(x)\right>^{(l)}}$
$\displaystyle\approx\sum_{k,j}^{N_{S}^{(l)},\tilde{N}_{S}^{(l)}}\,(\mathrm{Res}_{\gamma})^{(l)}_{kj}\psi_{m}^{(l)}(\xi_{k})P_{n}(x_{j})w_{k}^{(l)}v_{j}^{(l)},$
(53)

where $l=1,2$, $m=0,1,..,N_{\gamma}^{(l)}-1$,
$n=0,1,..,\tilde{N}_{\gamma}^{(l)}$, and $w_{k}^{(l)},v_{j}^{(l)}$ are the
weights associated to the radial and angular test functions, respectively. It
is worth mentioning that we have
$N_{\gamma}^{(1)}(\tilde{N}_{\gamma}^{(1)}+1)+N_{\gamma}^{(2)}(\tilde{N}_{\gamma}^{(2)}+1)$
equations relating the modes $a_{mn,u}^{(l)}$ with the values of the several
quantities present in the residual evolution equation. The remaining
$(\tilde{N}_{\gamma}^{(1)}+1)+(\tilde{N}_{\gamma}^{(2)}+1)$ relations are from
the transmission conditions:

$\displaystyle\int_{-1}^{1}\,\left(\bar{\gamma}_{,u}^{(1)}(u,r^{(1)},x)-\bar{\gamma}_{,u}^{(2)}(u,r^{(1)},x)\right)\,P_{j}(x)dx=0,$
(54)

where $j=0,1,..,\tilde{N}_{\gamma}^{(2)}$ and

$\displaystyle\int_{-1}^{1}\,\left(\bar{\gamma}_{,ur}^{(1)}(u,r^{(1)},x)-\bar{\gamma}_{,ur}^{(2)}(u,r^{(1)},x)\right)\,P_{j}(x)dx=0,$
(55)

where $j=0,1,..,\tilde{N}_{\gamma}^{(1)}$. The set of equations (53), (54) and
(55) can be solved to obtain the set of dynamical equations for the modes
$a_{mn}^{(l)}(u)\,\rightarrow\,a_{k}(u)$. Schemactically, we have

$\displaystyle\frac{da_{k}}{du}=F_{k}\left(\bar{\gamma}_{ij},(\bar{\gamma}_{,x})_{ij},(\bar{\gamma}_{,xx})_{ij},..\right),$
(56)

where
$k=1,2,..,(N_{\gamma}^{(1)}+1)(\tilde{N}_{\gamma}^{(1)}+1)+(N_{\gamma}^{(2)}+1)(\tilde{N}_{\gamma}^{(2)}+1)$,
and $\bar{\gamma}_{ij},(\bar{\gamma}_{,x})_{ij},(\bar{\gamma}_{,xx})_{ij},..$
are the values of these quantities calculated in the second grid (48) and
(49).

##  IV Code validation

In all numerical experiments, we take advantage of the field equation’s
symmetry concerning the angular dependence of the functions
$\bar{\gamma},\bar{\beta},\bar{U}$ and $S$. Let us consider that
$\bar{\gamma}(u,r,x)$ is an even function of $x$. Eq. (4) always demands that
$\bar{\beta}(u,r,x)$ is an even function of $x$. After inspecting Eq. (5), it
follows that $\bar{U}(u,r,x)$, and consequently $\bar{Q}(u,r,x)$, must be an
odd function of $x$. Inserting these informations into Eq. (6), one can show
that $S(u,r,x)$ is also an even function of $x$. We can define parity in $x$
in the spectral approximations (24) - (28) by setting the angular basis
functions to $P_{2j}(x)$ and $P_{2j+1}(x)$, respectively, for the even and odd
functions of $x$. A direct consequence, we increase the angular resolution by
placing the collocation points $x_{j}$ in the interval $0\leq x\leq 1$.

Figure 2: Exponential decay of maximum deviation $C_{\mathrm{max}}$ for the
initial data (58) with $A_{0}=30$ (squares) and $50$ (circles).
$N_{\gamma}^{(1)}$ is the truncation order in the first domain taken as a
reference. Here $y^{(1)}=0$ (cf. Eq. (19)) and $L_{0}=0.5$ for $A_{0}=30$, and
$L_{0}=0.25$ for $A_{0}=50$. The interface is located at $r^{(1)}=0.5$ and
$r^{(1)}=0.25$, respectively.

Figure 3: Exponential decay of maximum deviation $C_{\mathrm{max}}$ for the
initial data (58) with $A_{0}=30,40$ and $50$ represented by black circles,
blue squares, and red diamonds, respectively . $N_{\gamma}^{(1)}$ is the
truncation order in the first domain taken as a reference. Here $L_{0}=1.0$
(upper panel) and $L_{0}=0.5$ (lower panel). We have set $y^{(1)}=-0.5$ (cf.
Eq. (19)) yielding the interface located at $r^{(1)}=L_{0}/3$.

We have additionally validated the domain decomposition code with the
verification of the Bondi formula. Following Gomez et al. papadopoulos , we
introduce the function $C(u)$ that measures the deviation from the Bondi
formula or the global energy conservation:

$\displaystyle
C(u)=\frac{1}{M_{B}(u_{0})}\Bigg{[}M_{B}(u)-M_{B}(u_{0})+\frac{1}{2}\times$
$\displaystyle+\int_{u_{0}}^{u}\int_{-1}^{1}\left(\frac{\mathrm{e}^{2H}}{\omega}\mathcal{N}(u,x)^{2}dx\right)\,du\Bigg{]}\times
100,$ (57)

Figure 4: Exponential decay of maximum deviation $C_{\mathrm{max}}$ for the
initial data (59) with $r_{0}=3.0$, $\sigma=1$, and $A_{0}=0.1$.
$N_{\gamma}^{(1)}$ is the truncation order in the first domain taken as a
reference. Here $L_{0}=5$ and $y^{(1)}=0$ (cf. Eq. (19)) yielding the
interface located at $r^{(1)}=L_{0}=5$.

where $M_{B}(u_{0})$ is the initial Bondi mass. We have presented in Appendix
the corresponding expressions for the Bondi mass, the news function,
$\mathcal{N}(u,x)$, and $\omega(u,x)$. For this test, we have considered the
following initial data

$\displaystyle\bar{\gamma}_{0}(r,x)$ $\displaystyle=$ $\displaystyle
216A_{0}\frac{r^{6}(1-x^{2})^{2}}{(1+2r)^{9}},$ (58)
$\displaystyle\bar{\gamma}_{0}(r,x)$ $\displaystyle=$
$\displaystyle\frac{6A_{0}r^{2}(1+r)}{(1+2r)^{3}}(1-x^{2})\mathrm{e}^{-(r-r_{0})^{2}/\sigma^{2}},$
(59)

where $A_{0}$ plays the role of the initial amplitude. The first initial data
is a pulse of gravitational wave more concentrated near the origin, whereas
the second represents a gravitational wave Gaussian-like packet centered at
$r_{0}$ with width $\sigma$.

To verify the convergence of the error associated with the Bondi formula,
$C(u)$, we have evolved the field equations with increasing resolutions. We
opted for a simple selection of the same truncation orders or the same number
of collocation points in each subdomain. The radial and angular resolutions
for the function $\bar{\gamma}$ in each subdomain satisfy are related by

$N_{\gamma}=2\tilde{N}_{\gamma}.$

The radial and angular truncation orders of $\bar{\beta},\bar{U},\bar{Q}$
field functions are $N_{\beta}=N_{U}=N_{Q}=N_{\gamma}+2$ and
$\tilde{N}_{\beta}=\tilde{N}_{U}=\tilde{N}_{Q}=\tilde{N}_{\gamma}+2$,
respectively. For the metric function $S$, we have $N_{S}=N_{U}+2$ and
$\tilde{N}_{S}=\tilde{N}_{U}+1$.

One of the considerable technical difficulties resided in calculating the
Bondi mass $M_{B}(u)$. This task involves the determination of the mass aspect
$M(u,\theta)$ (cf. Eq. (15)) along with the gauge terms $H,K,c$ and $L$ that
appears in the asymptotic expansions of the metric functions $\beta,\gamma$,
and $U$ as depicted by Eqs. (12) - (14). In the present numerical scheme,
these quantities are obtained directly from the spectral approximations (24) -
(28) since all basis functions satisfy the asymptotic conditions (12) - (15).
Further, we need the same quantities list above to calculate the news function
$\mathcal{N}(u,x)$. Therefore, we adapted in the code the calculations of the
Bondi mass $M_{B}(u)$ and the news functions allowing us to obtain the
deviation $C(u)$ given by Eq. (57).

In the Appendix, we indicated the approximations used to calculate the Bondi
mass and the integrals with the news functions found in the definition of
$C(u)$ (cf. Eq. (57)). For the initial data (58), we have tested distinct
values of the initial amplitude $A_{0}$, the map parameter, the location of
interfaces. The exponential decay of the maximum deviation was observed in all
cases. In Fig. 2, the red (squares) lines correspond to the initial amplitude
$A_{0}=30$, $y^{(1)}=0$, and map parameter $L_{0}=0.5$ with the interface at
$r^{(1)}=0.5$. The saturation of $C_{\mathrm{max}}$ occurs for $N_{\gamma}\geq
32$ and is about $10^{-6}\%$. In the same panel, the blue line (circles)
corresponds to the following parameters: $A_{0}=50$, $y^{(1)}=0$, and
$L_{0}=0.25$ producing an interface at $r^{(1)}=0.25$. Even for
$N^{(1)}_{\gamma}=50$ the saturation is not achieved, but the exponential
decay is satisfactory. The difference in the interface locations in both cases
was necessary to concentrating more collocation points near the origin for the
higher amplitude. In general, $C_{\mathrm{max}}$ occurs when the Bondi mass
varies more rapidly, therefore being dependend on the initial amplitude. For
$A_{0}=30,50$ (Fig. 2) we found $C_{\mathrm{max}}\approx 0.5$. Although not
presented here, the results with domain decomposition are better than those
obtained for a single domain with the same number of collocation points.

We additionally explored the error decay after fixing $y^{(1)}=-0.5$ such that
the interface is now located at $r^{(1)}=L_{0}/3$. In Fig. 3 we present the
exponential decay of $C_{\mathrm{max}}$ for the initial data (58) with
amplitudes of $A_{0}=30$ (black circles), $A_{0}=40$ (blue squares) and
$A_{0}=50$ (red diamonds). The results are better for $L_{0}=1.0$ (upper
panel) than for $L_{0}=0.5$ (lower panel).

We have set the initial amplitude $A_{0}=0.1$ and $\sigma=1$ for the Gaussian-
type initial pulse of gravitational waves located at $r_{0}=3$. Although the
amplitude seems small, it is enough to excite the nonlinearities of the field
equations. Fig. 4 shows the log-plot of the maximum deviation of the Bondi
formula, $C_{\mathrm{max}}$. It becomes clear the exponential decay until
saturation of about $10^{-5}\%$ is achieved, in other words, a deviation of
one part in $10^{7}$. We have integrated the approximate field equations with
a fourth-order Runge-Kutta integrator. The map parameter is $L_{0}=5.0$ and
$y^{(1)}=0$ yielding the interface located at $r^{(1)}=L_{0}=5.0$.

We close this section by providing some information about the computational
time-consuming in the numerical experiments. As mentioned, we have performed
the numerical simulations using a fourth-order Runge-Kutta integrator with a
fixed step size.

The initial profiles (58) and (59) demand distinct choices for the map
parameter and, together with the resolution or truncation orders, influence
the step size selection. The initial data (58) represents a gravitational wave
packet more concentrated near the origin, which requires a smaller map
parameter. On the other hand, a relatively more significant map parameter is
more appropriate for a Gaussian-like wave packet (59) located at $r_{0}=3$.

We indicate the value of $N_{\gamma}$ to characterize the other truncation
orders and establish the time needed to integrate one unit of the central
time, i. e. $\Delta u=1$. All numerical experiments were performed on a
notebook Linux Ubuntu with processor $i7\,9750H$ and memory RAM of 32 GB. With
the initial data (58), interface located at $y^{(1)}=0$, map parameter
$L_{0}=0.5$, we have obtained: $57.28$ minutes for $N_{\gamma}=28$, $96.81$
minutes for $N_{\gamma}=32$, and $4.59$ hours for $N_{\gamma}=36$ with the
step sizes of $1.0\times 10^{-5},1.0\times 10^{-5}$, and $0.5\times 10^{-5}$,
respectively. There is a caveat: in general the deviation reaches its maximum
value for $u<1.0$. By choosing now $y^{(1)}=-0.5$ rendering $r^{(1)}=L_{0}/3$
and $L_{0}=0.5$, we have obtained the following time costs: $119.2$ minutes
for $N_{\gamma}=28$, $473.7$ minutes for $N_{\gamma}=32$, and $11$ hours for
$N_{\gamma}=36$ with the step sizes of $5.0\times 10^{-6},2.0\times 10^{-6}$,
and $2.0\times 10^{-6}$, respectively.

For the initial data (59) with $r_{0}=3$, we have set $L_{0}=5$ and the
interface $y^{(1)}=0$. The time-costs are the following: $5.7$ minutes for
$N_{\gamma}=28$, $19.42$ minutes for $N_{\gamma}=32$, $27.45$ minutes for
$N_{\gamma}=36$, and $100$ minutes for $N_{\gamma}=40$. The step sizes are:
$1.0\times 10^{-4}$, $5.0\times 10^{-5}$, $5.0\times 10^{-5}$, and $2.0\times
10^{-5}$, respectively.

Figure 5: Quadratic dependence of the Bondi mass with the initial amplitude
for $u=0,\,1,\,2$ (top plot). For each time the Bondi mass is normalized with
the Bondi mass for the smallest considered amplitude ($A_{0}=10$). Bondi mass
normalized by the initial Bondi mass for each amplitude as a function of time
(bottom plot). We have used $N_{\gamma}^{(1)}=30$ for $A_{0}=10,\,30$ and
$N_{\gamma}^{(1)}=40$ for other amplitudes. Figure 6: Power-law late time
decay of the Bondi masses for $A_{0}=10$ for the initial data given by (58).
Here $M_{B}(u)\sim u^{-6.59}$.

##  V The Bondi mass decay and the gravitational wave patterns

As a physical application, we generate the gravitational waveforms emitted by
an initially ingoing pulse towards the origin. It is well established that
asymptotic flat spacetimes admitting gravitational waves have a nonvanishing
term of order $r^{-1}$ in the Weyl tensor when viewed far from the source
newman_penrose . It means, according to the Peeling Theorem sachs , that the
curvature tensor has the same algebraic structure as a plane wave.
Consequently, we are interested in the term of the Weyl scalar $\Psi_{4}$ that
falls as $r^{-1}$ at the future null infinity.

In the original Bondi frame bondi  characterized by $H=L=K=0$, one can show,
after choosing an appropriate null tetrad basis (see Appendix), that

$\Psi_{4}(u,r,\theta)\simeq-\frac{c(u,\theta)_{,uu}}{r},$ (60)

where $c(u,\theta)=\lim_{r\rightarrow\infty}\,(r\gamma)$ (cf. Eq. (12)). In
the frame we are adopting, the Weyl scalar has a similar structure:

$\displaystyle\Psi_{4}(u,r,\theta)=\frac{\Psi^{(1)}_{4}(u,\theta)}{r}+\mathcal{O}(r^{-2}),$
(61)

with the term $\Psi^{(1)}_{4}(u,\theta)$ describing the waveforms perceive by
a distant observer. The full expression of $\Psi^{(1)}_{4}(u,\theta)$ is quite
large bishop_rezolla , but, for the sake of convenience, we exhibit this
function evaluated at $\theta=\pi/2$ or $x=0$. It reads

$\displaystyle\Psi^{(1)}_{4}$
$\displaystyle(u,x=0)=\bigg{\\{}\mathrm{e}^{2(\bar{H}-\bar{K})}\bigg{[}\left(\bar{H}_{,xx}-\frac{1}{2}\bar{K}_{,xx}+\bar{K}-4\bar{H}+1\right)\times$
(62)
$\displaystyle\bar{L}_{,x}+\left(\bar{K}_{,xx}+2\bar{H}_{,xx}-2\bar{K}-8\bar{H}+1\right)\,\bar{K}_{,u}-\frac{1}{4}\bar{L}_{,x,x,x}$
$\displaystyle+4\bar{H}_{,u}-\bar{H}_{,uxx}\bigg{]}-\bar{H}_{,u}\bar{L}_{,x}\bar{c}-\frac{1}{2}\bar{L}_{,x}^{2}\bar{c}+\frac{3}{2}\bar{L}_{,x}\bar{c}_{,u}$
$\displaystyle+\frac{1}{2}\bar{L}_{,ux}\bar{c}+2\bar{H}_{,u}\bar{c}_{,u}-\bar{c}_{,uu}\bigg{\\}}_{x=0},$

where we have taken into account the symmetries of the metric functions with
respect to $x=0$, $H=(1-x^{2})^{2}\bar{H}$, $K=(1-x^{2})\bar{K}$,
$L=\sqrt{1-x^{2}}\bar{L}$ and $c=(1-x^{2})\bar{c}$. As expected, we recover
the original expression (60) of $\Psi_{4}$ in the Bondi frame.

Figure 7: Upper panel: plot of $\Psi_{4}^{(1)}(u)$ together with the inset
showing the corresponding logplot. The amplitude is $A_{0}=10$. Lower panel:
logplot of $\Psi_{4}^{(1)}(u)$ for $A_{0}=10$ (dashed line) and $A_{0}=30$
continuous line.

In the sequence, we present the wave patterns generated by different initial
incoming pulses described by $\bar{\gamma}(u_{0},r,x)=\bar{\gamma}_{0}(r,x)$.
We first considered the initial data Eq. (58) previously used for the code
validation. In that case, the pulse is close to the symmetry axis, meaning
that it bounces after hitting the symmetry axis in a short time interval. As a
consequence, gravitational waves carry away most of the mass.

Figure 5 (top plot) displays the expected quadratic dependence of the Bondi
mass with the initial amplitude for $u=0,1,2$. For each time the Bondi mass is
normalized with the Bondi mass for the smallest considered amplitude
($A_{0}=10$). Clearly, with evolution the initial linear dependence (in the
log-log plot of the Bondi mass) deviates. This is a measure of the transient
non-linear effects. In the sequence, in the bottom plot we display the Bondi
mass normalized by the initial Bondi mass for each amplitude (whose dependence
with the initial amplitude is shown in the top plot) as a function of time.
Certainly, in a different way, the bottom plot displays the non-linear effects
for different initial amplitudes and under the action of time. The Bondi mass
decay corresponding to the initial amplitudes $A_{0}=10,\,30,\,60,\,70$ and
$75$. In all cases, the gravitational waves extract about more than $50\%$ of
the initial mass when $u\approx 1.0$. The relative amount of mass loss is
greater for smaller initial amplitudes; more precisely, we have found that for
$A_{0}=10,30,60,70,75$, about $80.8\%,78.3\%,69.4\%,67.15\%,64.28\%$
respectively, of the initial mass is radiated away at $u\approx 1.0$. For
higher amplitudes, $A_{0}=60,70,75$, there is a fast initial mass decay as
soon as the pulse moves towards the symmetry axis. It seems to us that the
behavior of $\dot{M}_{B}\neq 0$ at $u=0$ is related with the backscattering
and consequently with the initial flux of energy at infinity; the
backscattering is stronger for the higher initial amplitude. One can think it
is a non linear effect, but a more careful study is required because a
previous result indicates that could be present even in the linear regime. In
this sense, we note that the same behavior appears in the setting studied in
Ref. b14  for the Einstein-Klein-Gordon system, using ingoing null cones. In
this work, which assumes the same symmetry, $\dot{M}_{B}\neq 0$ at $u=0$ is
stronger with the increasing of the angular structure (in the linear regime).

We have also observed (see Fig. 6) the late time power-law decay of the Bondi
mass, $M_{B}(u)\sim u^{-6.59}$. We have set $A_{0}=10$, but the same law is
obtained for greater initial Bondi masses.

(a)

(b)

Figure 8: (a) Log-linear plots of $\Psi_{4}^{(1)}(u)$ for $A_{0}=60,70$
represented by dashed and dotted lines, respectively. (b) Zoom of the region
$0\leq u\leq 0.8$ to visualize the appearence of small oscillations. The
continuous line corresponds to $A_{0}=75$.

We expect that the nuances of the mass decay leave imprints in the observed
pattern of the gravitational waves described by Eq. (62). Fig. 7 (upper panel)
shows the plots of $\Psi_{4}^{(1)}(u)$ corresponding to $A_{0}=10$ together
with its log-plot in the inset. There is a rapid variation $\Psi_{4}^{(1)}(u)$
in the interval $0\leq u\leq 1$ and an oscillation. We can better visualize
the effect of increasing the initial amplitude from $A_{0}=10$ to $A_{0}=30$
on the wave pattern with the log-plot of $\Psi_{4}^{(1)}(u)$ as shown in Fig.
7 (lower panel). We observe that the original oscillation takes place in a
shorter time interval, meaning increasing the frequency.

Furthermore, a more significant change occurs in the highly nonlinear regime
in which $A_{0}=60,70$ and $75$. The Bondi mass decay is similar in all cases,
but increasing the initial amplitude produces more oscillations in shorter
time intervals, as shown in Fig. 8. The dashed and dotted lines in Fig. 8(a)
correspond to $A0=60$ and $70$, respectively. In the zoom displayed in Fig.
8(b), with the plots of $\Psi_{4}^{(1)}(u)$ for $A_{0}=60,70$ and $75$
(continuous line).

(a)

(b)

Figure 9: (a) Decay of the Bondi mass for the initial data (59) with
$A_{0}=0.05$ (dashed-dotted line), $A_{0}=0.10$ (dashed line) and $A_{0}=0.12$
(continuous line). (b) Following the same convention, we show the behavior of
$\Psi_{4}^{(1)}(u)$ for $A_{0}=0.05,0.10$ and $A_{0}=0.12$. Figure 10: Power-
law late time decay of the Bondi mass for $A_{0}=0.3$ for the initial data
given by (59). Here $M_{B}(u)\sim u^{-9.83}$.

We have now considered the initial data (59), but with $r_{0}=4$ and
$\sigma=1$. The simulations have taken a long time necessary for the pulse to
hit the symmetry axis and bounce away. In Fig. 9(a), we present the decay of
the Bondi mass for $A_{0}=0.05,0.10$ and $A_{0}=0.12$, where the nonlinear
regime progressively enters into action. Notice the substantial increase in
the initial Bondi mass when the amplitude $A_{0}$ changes. From $A_{0}=0.05$
to $A_{0}=0.10$, the initial Bondi mass increases almost four times, and from
$A_{0}=0.10$ to $A_{0}=0.12$, the initial Bondi mass increases approximately
$50\%$. Since the initial pulse is not close to the symmetry axis, there is no
considerable mass loss in the first moments; only starting at $u\approx 6$
there is a huge mass extraction by gravitational radiation, mainly for
$A_{0}=0.10$ and $A_{0}=0.12$. The resulting patterns of $\Psi_{4}^{(1)}(u)$
depicted in Fig. 9(b) are sensitive to the decay of the Bondi mass
corresponding to the initial data under consideration. The considerable
variation of the Bondi mass left its imprint producing the oscillatory pattern
$\Psi_{4}^{(1)}(u)$. Finally, in Fig. 10, we present the power-law behavior of
the Bondi mass at late times, where we have found $u\sim u^{-9.83}$ which is
independent of the initial amplitude.

##  VI Final comments

We have implemented a simple domain decomposition spectral code based on the
Galerkin-Collocation method applied to the field equations of the celebrated
Bondi problem. The present task is a part of a program of implementing domain
decomposition in problems of interest in numerical relativity
alcoforado_critical ; alcoforado_cauchy ; barreto_DD ; barreto_cylindrical_DD
. We have presented with detail the following topics:

  1. (i)

The new computational scheme of establishing the non-overlapping subdomains
after the compactification of the spatial domains (cf. Fig. 1).

  2. (ii)

The definition of radial and angular basis functions for the metric functions
$\gamma,\,\beta,\,U$, and $V$.

  3. (iii)

The use of two sets of collocation points in each subdomain connected with the
functions $\beta,\,U$ and $V$.

  4. (iv)

The incorporation of the transmission conditions into the field equations. In
particular, we have considered the most simple form (see Refs.
alcoforado_critical ; alcoforado_cauchy  ) for the evolution equation, since
for hyperbolical problems, there is no unique form for the transmission
conditions.

  5. (v)

The approximations of the hypersurface equations as systems of algebraic
equations and the evolution equation as a system of ordinary differential
equations for the modes associated with the metric function $\gamma$.

The combination of features of both Galerkin and Collocation methods rendered
a stable and exponentially convergent code which is attested with the
verification of the Bondi formula. We have proceeded with the decay of the
Bondi mass for high amplitude initial data, and we have exhibited its late
time power-law decay. Since the decay of the Bondi mass is connected with the
gravitational wave extraction, we have shown the waveforms at the future null
infinity. We obtain the wave-forms from the Weyl scalar $\Psi_{4}$, which we
have calculated in the present gauge, that is, the one characterized by non-
vanishing asymptotic quantities $H,\,K$, and $L$.

We indicate some directions of the present work that are under investigation.
The first is to study the collapse of gravitational waves more thoroughly with
or without a matter content, such as a massless scalar field. In this case,
examining the critical behavior in axisymmetric collapse is of great interest.
Several works on the collapse of Brill waves have touched this issue without a
complete understanding choptuik_03 ; hilditch_13 ; ledvinka_21 . It is
possible to include more subdomains to increase the code resolution for a
better description of the implosion of gravitational waves.

## Acknowledgments

W. B. acknowledges the financial support by Universidade do Estado do Rio de
Janeiro (UERJ) and Fundação Carlos Chagas Filho de Amparo a Pesquisa do Estado
do Rio de Janeiro (FAPERJ); also thanks the hospitality of the Departamento de
Física Teórica at Instituto de Física (UERJ). H. P. O. acknowledges the
financial support of the Brazilian Agency CNPq and the support of FAPERJ. This
study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de
Nível Superior - Brasil (CAPES) - Finance Code 001.

##  Appendix A Bondi mass and the news function

The Bondi mass in the present gauge isaacson  is given by

$\displaystyle M_{B}(u)$
$\displaystyle=\int_{-1}^{1}\,\frac{\mathrm{e}^{-2K}}{\omega}\bigg{[}\frac{\mathrm{e}^{2K}}{2}M+\frac{1}{4}(1-x^{2})c^{\prime\prime}-xc^{\prime}-\frac{1}{2}c\times$
(63)
$\displaystyle(1-x^{2})c^{\prime}(H^{\prime}+K^{\prime})-(1-x^{2})c\bigg{(}H^{\prime
2}-2H^{\prime}K^{\prime}-$ $\displaystyle-K^{\prime
2}+\frac{1}{2}H^{\prime\prime}+\frac{1}{2}K^{\prime\prime}\bigg{)}+2xc(H^{\prime}+K^{\prime})\bigg{]}dx,$

where

$\omega(u,x)=\frac{2\,\mathrm{e}^{2K}}{(1+x)\mathrm{e}^{\Delta}+(1-x)\mathrm{e}^{-\Delta}},$

and

$\Delta(u,x)=\int_{0}^{x}\,\frac{(\mathrm{e}^{2K}-1)}{1-x^{2}}\,dx.$

Here prime means derivative with respect to $x=\cos\theta$. The news function
can be written as

$\displaystyle\mathcal{N}$
$\displaystyle=\mathrm{e}^{2H}c_{,u}-\frac{\mathrm{e}^{-2H}}{2c}(c^{2}L)^{\prime}+\frac{1}{2}(1-x^{2})\mathrm{e}^{-2K}+$
(64)
$\displaystyle+\frac{1}{2}(1-x^{2})\mathrm{e}^{-2K}\bigg{(}2H^{\prime\prime}+4H^{\prime
2}+\frac{\omega^{\prime\prime}}{\omega}-2\frac{\omega^{\prime
2}}{\omega^{2}}\bigg{)}.$

We have calculated the integral (A1) using Gauss-Legendre quadrature formulae
schematically as

$M_{B}(u)\approx\sum_{k=0}^{N_{q}}\,\left(...\right)_{k}wq_{k},$ (65)

where $\left(...\right)_{k}$ means the integrand evaluated at the quadrature
collocation points, $N_{q}$ indicates the number of quadrature points, and
$wq_{k}$ represents the weights. Similarly, we calculate the integral
involving the news function (RHS of Eq. (16)), and the integration with
respect to $u$ (cf. Eq. 57) is performed using the Newton-Cotes formula.

##  Appendix B Tetrad basis

We have used the following tetrad basis for the calculation of the Weyl scalar
$\Psi_{4}$ nerozzi :

$\displaystyle l^{\mu}$ $\displaystyle=$
$\displaystyle\left(0,-\mathrm{e}^{-2\beta},0,0\right)$ (66) $\displaystyle
n^{\mu}$ $\displaystyle=$
$\displaystyle\left(-1,\frac{V}{2r}-\frac{1}{2}U^{2}r^{2}\mathrm{e}^{-2(\gamma-\beta)},0,0\right)$
(67) $\displaystyle m_{\mu}$ $\displaystyle=$
$\displaystyle\left(0,-\frac{1}{\sqrt{2}}U\mathrm{e}^{\gamma-2\beta},\frac{1}{\sqrt{2}r}\mathrm{e}^{-\gamma},\frac{i}{\sqrt{2}r\sin\theta}\mathrm{e}^{\gamma}\right)$
(68)

## References

  * (1) B. P. Abbott et al. (LIGO Scientific Collaboration and V. Collaboration), Phys. Rev. Lett. 116, 061102 (2016).
  * (2) F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005).
  * (3) M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006).
  * (4) J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Centrella, and B. J. Kelly, Phys. Rev. Lett. 99, 181101 (2007).
  * (5) A. H. Mroué et al., Phys. Rev. Lett. 111, 241104 (2013).
  * (6) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. J. Forteza, and A. Bohé, Phys. Rev. D 93, 044006 (2016).
  * (7) H. Bondi, M. van der Burgh and A. Metzner, Proc. Roy. Soc. (London) A 269, 21 (1962).
  * (8) R. Sachs, Proc. Roy. Soc. (London) A 270, 103 (1962); ibid A 265, 463 (1962).
  * (9) R. Penrose, Phys. Rev. Lett. 10, 66 (1963).
  * (10) J. Winicour, Characteristic Evolution and Matching, Living Rev. Rel. 15, 2 (2012).
  * (11) R. Goméz, P. Papadopoulos and J. Winicour, J. Math. Phys. 35, 4184 (1994).
  * (12) R. A. Isaacson, J. S. Welling and J. Winicour, J. Math. Phys. 24, 1824 (1983).
  * (13) N. T. Bishop, R. Gómez, L. Lehner, M. Maharaj and J. Winicour, Phys. Rev. D 56 6298 (1997).
  * (14) M. Babiuc, B. Szilágyi, I. Hawke, and Y. Zlochower, Class. Quantum Grav. 22, 5089 (2005).
  * (15) R. Gómez, R. L. Marsa and J. Winicour, Phys. Rev. D 56 6310 (1997).
  * (16) M. Babiuc, J. Winicour and Y. Zlochower, Class. Quantum Grav. 28, 134006 (2011).
  * (17) M. Babiuc, B. Szilágyi, J. Winicour and Y. Zlochower, Phys. Rev. D 84, 044057 (2011).
  * (18) H. P. de Oliveira and E. L. Rodrigues, Class. Quantum Grav. 28, 235011 (2011).
  * (19) C. J. Handmer and B. Szilagyi, Class. Quantum Grav. 32, 025008 (2015).
  * (20) B. A. Finlayson, The method of weighted residuals and variational principles, Academic Press (1972).
  * (21) J. P. Boyd, Chebyshev and Fourier Spectral Methods (Dover Publications, New York, 2001).
  * (22) C. Canuto, A. Quarteroni, M. Y. Hussaini and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag (1988).
  * (23) H. P. de Oliveira, E. L. Rodrigues, and J. E. F. Skea, Phys. Rev. D 82, 104023 (2010).
  * (24) J. A. Crespo and H. P. de Oliveira, Phys. Rev. D 92, 064004 (2015).
  * (25) W. Barreto, J. A. Crespo, H. P. de Oliveira, E. L. Rodrigues, and B. Rodriguez-Mueller, Phys. Rev. D 93, 064042 (2016).
  * (26) J. A. Crespo, H. P. de Oliveira, and J. Winicour, Phys. Rev. D 100, 104017 (2019).
  * (27) M. A. Alcoforado, W.O. Barreto, and H. P. de Oliveira, Gen. Relativ. Gravit. 53, 42 (2021).
  * (28) M. A. Alcoforado, R. F. Aranha, W. O. Barreto, and H. P. de Oliveira, Multidomain Galerkin-Collocation method: spherical collapse of scalar fields II, arXiv:2012.01302v1.
  * (29) J. Celestino, H. P. de Oliveira and E. L. Rodrigues, Phys. Rev. D 93, 104018 (2016).
  * (30) W. O. Barreto, J. A. Crespo, H. P. de Oliveira, and E. L. Rodrigues, Class. Quantum Grav. 36, 215011 (2019).
  * (31) P. C. M. Clemente and H. P. de Oliveira, Phys. Rev. D 96, 024035 (2017).
  * (32) W. Barreto, P. C. M. Clemente, H. P. de Oliveira, and B. Rodriguez-Mueller, Phys. Rev. D 97, 104035 (2018)
  * (33) W. Barreto, P. C. M. Clemente, H. P. de Oliveira, and B. Rodriguez-Mueller, Gen. Relativ. Grav. D 50, 71 (2018).
  * (34) S. Bonazzola, E. Gourgoulhon, M. Salgado and J. A. Marck, Astron. Astrophys. 278, 421 (1993).
  * (35) H. P. Pfeiffer, Initial data for black hole evolutions, Ph.D. thesis (2003), arXiv:gr-qc/0510016.
  * (36) M. Ansorg, Class. Quantum Grav. 24 S1 (2007).
  * (37) Spectral Einstein Code, https://www.black-holes.org/code/SpEC.html
  * (38) LORENE (Langage Objet pour la Relativité Numérique), http://www.lorene.obspm.fr
  * (39) B. Szilágyi, L. Lindblom and M. A. Scheel, Phys. Rev. D 80, 124010 (2009).
  * (40) L. E. Kidder, S. E. Field, F Foucart, E. Schnetter, S. A. Teukolsky, A. Bohn, N. Deppe, P. Diener, F. Hébert, J. Lippuner, J. Miller, C. D. Ott, M. A. Scheel, and T. Vicent, J. Comp. Phys. 335, 84 (2017).
  * (41) F. Hébert, L. E. Kidder and S. A. Teukolsky, Phys. Rev. D 98, 0444041 (2018).
  * (42) D. A. Hemberger, M. A. Scheel, L. E. Kidder, B. Szilágyi, G. Lovelace, N. W. Taylor, S. A. Teukolsky, Class. Quantum Grav. 30, 115001 (2013).
  * (43) M. Boyle et al., The SXS Collaboration catalog of binary black hole simulations, arXiv: 1904.04831 (2019).
  * (44) L. A. Tamburino and J. Wincour, Phys. Rev. 150, 1039 (1966).
  * (45) R. Isaacson, J. Welling and J. Winicour, J. Math. Phys. 24, 1824 (1983).
  * (46) E. T. Newman and R. Penrose, J. Math. Phys. 3, 566 (1962); J. Math. Phys. 4, 998 (1963).
  * (47) N. T. Bishop and L. Rezolla, Extraction of gravitional waves in numerical relativity, Living Rev. Rel. 19, 2 (2016).
  * (48) W. Barreto, Phys. Rev. D 90, 024055 (2014).
  * (49) M. W. Choptuik, E. W. Hirschmann, S. L. Liebling, and F. Pretorius, Phys. Rev. D 68, 044007 (2003).
  * (50) D. Hilditch, T. W. Baumgarte, A. Weyhausen, T. Dietrich, B. Brügmann, P. J. Montero, and E. Müller, Phys. Rev. D 88, 13009 (2013).
  * (51) T. Ldevinka and A. Khirnov, Phys. Rev. lett. 127, 011104 (2021).
  * (52) A. Nerozzi, M. Bruni, V. Re, and L. M. Burko, Phys. Rev. D 73, 044020 (2006).


* Landulfo  _et al._ [2012] Landulfo, A. G. S., W. C. C. Lima, G. E. A. Matsas, and D. A. T. Vanzella (2012), Phys. Rev. D 86, 104025, arXiv:1204.3654 [gr-qc] .EOS
  * Landulfo  _et al._ [2015] Landulfo, A. G. S., W. C. C. Lima, G. E. A. Matsas, and D. A. T. Vanzella (2015), Phys. Rev. D 91 (2), 024011, arXiv:1410.2274 [gr-qc] .EOS
  * Lang [2014] Lang, R. N. (2014), Phys. Rev. D 89 (8), 084014, arXiv:1310.3320 [gr-qc] .EOS
  * Lattimer and Prakash [2016] Lattimer, J. M., and M. Prakash (2016), Phys. Rept. 621, 127, arXiv:1512.07820 [astro-ph.SR] .EOS
  * Lattimer and Schutz [2005] Lattimer, J. M., and B. F. Schutz (2005), Astrophys. J. 629, 979, arXiv:astro-ph/0411470 .EOS
  * Lee [1974] Lee, D. L. (1974), Phys. Rev. D 10, 2374.EOS
  * Liebling and Palenzuela [2012] Liebling, S. L., and C. Palenzuela (2012), Living Rev. Rel. 15, 6, arXiv:1202.5809 [gr-qc] .EOS
  * Lima  _et al._ [2010] Lima, W. C. C., G. E. A. Matsas, and D. A. T. Vanzella (2010), Phys. Rev. Lett. 105, 151102, arXiv:1009.1771 [gr-qc] .EOS
  * Lima  _et al._ [2013] Lima, W. C. C., R. F. P. Mendes, G. E. A. Matsas, and D. A. T. Vanzella (2013), Phys. Rev. D 87 (10), 104039, arXiv:1304.0582 [gr-qc] .EOS
  * Lima and Vanzella [2010] Lima, W. C. C., and D. A. T. Vanzella (2010), Phys. Rev. Lett. 104, 161102, arXiv:1003.3421 [gr-qc] .EOS
  * Lindblom [2010] Lindblom, L. (2010), Phys. Rev. D 82, 103011, arXiv:1009.0738 [astro-ph.HE] .EOS
  * Macedo  _et al._ [2019] Macedo, C. F. B., J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou (2019), Phys. Rev. D 99 (10), 104041, arXiv:1903.06784 [gr-qc] .EOS
  * Maione  _et al._ [2016] Maione, F., R. De Pietri, A. Feo, and F. Löffler (2016), Class. Quant. Grav. 33 (17), 175009, arXiv:1605.03424 [gr-qc] .EOS
  * Maselli  _et al._ [2020] Maselli, A., N. Franchini, L. Gualtieri, and T. P. Sotiriou (2020), Phys. Rev. Lett. 125 (14), 141101, arXiv:2004.11895 [gr-qc] .EOS
  * Maselli  _et al._ [2022] Maselli, A., N. Franchini, L. Gualtieri, T. P. Sotiriou, S. Barsanti, and P. Pani (2022), Nature Astron. 6 (4), 464, arXiv:2106.11325 [gr-qc] .EOS
  * Mayo and Bekenstein [1996] Mayo, A. E., and J. D. Bekenstein (1996), Phys. Rev. D 54, 5059, arXiv:gr-qc/9602057 .EOS
  * McDermott  _et al._ [1983] McDermott, P. N., H. M. van Horn, and J. F. Scholl (1983), ApJ 268, 837.EOS
  * Mendes [2015] Mendes, R. F. (2015), Phys. Rev. D 91 (6), 064024, arXiv:1412.6789 [gr-qc] .EOS
  * Mendes and Ortiz [2016] Mendes, R. F., and N. Ortiz (2016), Phys. Rev. D 93 (12), 124035, arXiv:1604.04175 [gr-qc] .EOS
  * Mendes and Ortiz [2018] Mendes, R. F., and N. Ortiz (2018), Phys. Rev. Lett. 120 (20), 201104, arXiv:1802.07847 [gr-qc] .EOS
  * Mendes and Ottoni [2019] Mendes, R. F., and T. Ottoni (2019), Phys. Rev. D 99 (12), 124003, arXiv:1903.11638 [gr-qc] .EOS
  * Mendes  _et al._ [2014a] Mendes, R. F. P., G. E. A. Matsas, and D. A. T. Vanzella (2014a), Phys. Rev. D 90 (4), 044053, arXiv:1407.6405 [gr-qc] .EOS
  * Mendes  _et al._ [2014b] Mendes, R. F. P., G. E. A. Matsas, and D. A. T. Vanzella (2014b), Phys. Rev. D 89 (4), 047503, arXiv:1310.2185 [gr-qc] .EOS
  * Miller  _et al._ [2019] Miller, M. C., _et al._ (2019), Astrophys. J. Lett. 887 (1), L24, arXiv:1912.05705 [astro-ph.HE] .EOS
  * Minamitsuji [2020a] Minamitsuji, M. (2020a), Phys. Rev. D 101 (10), 104044, arXiv:2003.11885 [gr-qc] .EOS
  * Minamitsuji [2020b] Minamitsuji, M. (2020b), Phys. Rev. D 102 (4), 044048, arXiv:2008.12758 [gr-qc] .EOS
  * Minamitsuji [2021] Minamitsuji, M. (2021), Phys. Rev. D 103 (8), 084002, arXiv:2104.03660 [gr-qc] .EOS
  * Minamitsuji and Ikeda [2019a] Minamitsuji, M., and T. Ikeda (2019a), Phys. Rev. D 99 (4), 044017, arXiv:1812.03551 [gr-qc] .EOS
  * Minamitsuji and Ikeda [2019b] Minamitsuji, M., and T. Ikeda (2019b), Phys. Rev. D 99 (10), 104069, arXiv:1904.06572 [gr-qc] .EOS
  * Minamitsuji and Silva [2016] Minamitsuji, M., and H. O. Silva (2016), Phys. Rev. D 93 (12), 124041, arXiv:1604.07742 [gr-qc] .EOS
  * Minamitsuji and Tsujikawa [2022] Minamitsuji, M., and S. Tsujikawa (2022), “Symmetry restoration in the vicinity of neutron stars with a nonminimal coupling,” arXiv:2208.08107 [gr-qc] .EOS
  * Mirshekari and Will [2013] Mirshekari, S., and C. M. Will (2013), Phys. Rev. D 87 (8), 084070, arXiv:1301.4680 [gr-qc] .EOS
  * Morisaki and Suyama [2017] Morisaki, S., and T. Suyama (2017), Phys. Rev. D 96 (8), 084026, arXiv:1707.02809 [gr-qc] .EOS
  * Motohashi and Minamitsuji [2018] Motohashi, H., and M. Minamitsuji (2018), Phys. Lett. B 781, 728, arXiv:1804.01731 [gr-qc] .EOS
  * Motohashi and Suyama [2012] Motohashi, H., and T. Suyama (2012), Phys. Rev. D 85, 044054, arXiv:1110.6241 [gr-qc] .EOS
  * Mou and Zhang [2022] Mou, Z.-G., and H.-Y. Zhang (2022), Phys. Rev. Lett. 129 (15), 151101, arXiv:2204.11324 [hep-th] .EOS
  * Müther  _et al._ [1987] Müther, H., M. Prakash, and T. L. Ainsworth (1987), Phys. Lett. B 199, 469.EOS
  * Myung and Zou [2018] Myung, Y. S., and D.-C. Zou (2018), Phys. Rev. D 98 (2), 024030, arXiv:1805.05023 [gr-qc] .EOS
  * Myung and Zou [2019a] Myung, Y. S., and D.-C. Zou (2019a), Int. J. Mod. Phys. D 28 (09), 1950114, arXiv:1903.08312 [gr-qc] .EOS
  * Myung and Zou [2019b] Myung, Y. S., and D.-C. Zou (2019b), Eur. Phys. J. C 79 (3), 273, arXiv:1808.02609 [gr-qc] .EOS
  * Myung and Zou [2019c] Myung, Y. S., and D.-C. Zou (2019c), Eur. Phys. J. C 79 (8), 641, arXiv:1904.09864 [gr-qc] .EOS
  * Myung and Zou [2021] Myung, Y. S., and D.-C. Zou (2021), Phys. Lett. B 814, 136081, arXiv:2012.02375 [gr-qc] .EOS
  * Nelmes and Piette [2012] Nelmes, S. G., and B. M. A. G. Piette (2012), Phys. Rev. D 85, 123004, arXiv:1204.0910 [astro-ph.SR] .EOS
  * Niu  _et al._ [2021] Niu, R., X. Zhang, B. Wang, and W. Zhao (2021), Astrophys. J. 921 (2), 149, arXiv:2105.13644 [gr-qc] .EOS
  * Novak [1998a] Novak, J. (1998a), Phys. Rev. D 58, 064019, arXiv:gr-qc/9806022 .EOS
  * Novak [1998b] Novak, J. (1998b), Phys. Rev. D 57, 4789, arXiv:gr-qc/9707041 .EOS
  * Novak and Ibanez [2000] Novak, J., and J. M. Ibanez (2000), Astrophys. J. 533, 392, arXiv:astro-ph/9911298 .EOS
  * Ofengeim [2020] Ofengeim, D. D. (2020), Phys. Rev. D 101 (10), 103029, arXiv:2005.03549 [astro-ph.HE] .EOS
  * Oliveira and Pombo [2021] Oliveira, J. a. M. S., and A. M. Pombo (2021), Phys. Rev. D 103 (4), 044004, arXiv:2012.07869 [gr-qc] .EOS
  * Ozel  _et al._ [2016] Ozel, F., D. Psaltis, T. Guver, G. Baym, C. Heinke, and S. Guillot (2016), Astrophys. J. 820 (1), 28, arXiv:1505.05155 [astro-ph.HE] .EOS
  * Palenzuela  _et al._ [2014] Palenzuela, C., E. Barausse, M. Ponce, and L. Lehner (2014), Phys. Rev. D 89 (4), 044024, arXiv:1310.4481 [gr-qc] .EOS
  * Palenzuela and Liebling [2016] Palenzuela, C., and S. L. Liebling (2016), Phys. Rev. D 93 (4), 044009, arXiv:1510.03471 [gr-qc] .EOS
  * Pani and Berti [2014] Pani, P., and E. Berti (2014), Phys. Rev. D 90 (2), 024025, arXiv:1405.4547 [gr-qc] .EOS
  * Pani  _et al._ [2011a] Pani, P., E. Berti, V. Cardoso, and J. Read (2011a), Phys. Rev. D 84, 104035, arXiv:1109.0928 [gr-qc] .EOS
  * Pani  _et al._ [2011b] Pani, P., V. Cardoso, E. Berti, J. Read, and M. Salgado (2011b), Phys. Rev. D 83, 081501, arXiv:1012.1343 [gr-qc] .EOS
  * Pappas and Apostolatos [2014] Pappas, G., and T. A. Apostolatos (2014), Phys. Rev. Lett. 112, 121101, arXiv:1311.5508 [gr-qc] .EOS
  * Pappas  _et al._ [2019] Pappas, G., D. D. Doneva, T. P. Sotiriou, S. S. Yazadjiev, and K. D. Kokkotas (2019), Phys. Rev. D 99 (10), 104014, arXiv:1812.01117 [gr-qc] .EOS
  * Pappas and Sotiriou [2015] Pappas, G., and T. P. Sotiriou (2015), Phys. Rev. D 91 (4), 044011, arXiv:1412.3494 [gr-qc] .EOS
  * Paschalidis and Stergioulas [2017] Paschalidis, V., and N. Stergioulas (2017), Living Rev. Rel. 20 (1), 7, arXiv:1612.03050 [astro-ph.HE] .EOS
  * Podkowka  _et al._ [2018] Podkowka, D. M., R. F. P. Mendes, and E. Poisson (2018), Phys. Rev. D 98 (6), 064057, arXiv:1807.01565 [gr-qc] .EOS
  * Ponce  _et al._ [2015] Ponce, M., C. Palenzuela, E. Barausse, and L. Lehner (2015), Phys. Rev. D 91 (8), 084038, arXiv:1410.0638 [gr-qc] .EOS
  * Popchev [2015] Popchev, D. (2015),  _Bifurcation of neutron star solutions in scalar-tensor theories of gravity_ , Master’s thesis (University of Sofia).EOS
  * Popchev  _et al._ [2019] Popchev, D., K. V. Staykov, D. D. Doneva, and S. S. Yazadjiev (2019), Eur. Phys. J. C 79 (2), 178, arXiv:1812.00347 [gr-qc] .EOS
  * Pretorius [2005] Pretorius, F. (2005), Phys. Rev. Lett. 95, 121101, arXiv:gr-qc/0507014 .EOS
  * Proca [1936] Proca, A. (1936), J. Phys. Radium 7, 347.EOS
  * Ramazanoğlu [2017] Ramazanoğlu, F. M. (2017), Phys. Rev. D96 (6), 064009, arXiv:1706.01056 [gr-qc] .EOS
  * Ramazanoğlu [2018a] Ramazanoğlu, F. M. (2018a), Phys. Rev. D97 (2), 024008, arXiv:1710.00863 [gr-qc] .EOS
  * Ramazanoğlu [2018b] Ramazanoğlu, F. M. (2018b), Phys. Rev. D98 (4), 044013, arXiv:1804.03158 [gr-qc] .EOS
  * Ramazanoğlu [2018c] Ramazanoğlu, F. M. (2018c), Phys. Rev. D98 (4), 044011, arXiv:1804.00594 [gr-qc] .EOS
  * Ramazanoğlu [2019a] Ramazanoğlu, F. M. (2019a), Phys. Rev. D 99 (4), 044003, arXiv:1901.00194 [gr-qc] .EOS
  * Ramazanoğlu [2019b] Ramazanoğlu, F. M. (2019b), Phys. Rev. D99 (8), 084015, arXiv:1901.10009 [gr-qc] .EOS
  * Ramazanoğlu [2019c] Ramazanoğlu, F. M. (2019c), Turk. J. Phys. 43 (6), 586.EOS
  * Ramazanoğlu and Pretorius [2016] Ramazanoğlu, F. M., and F. Pretorius (2016), Phys. Rev. D 93 (6), 064005, arXiv:1601.07475 [gr-qc] .EOS
  * Ramazanoğlu and Ünlütürk [2019] Ramazanoğlu, F. M., and K. I. Ünlütürk (2019), Phys. Rev. D 100 (8), 084026, arXiv:1910.02801 [gr-qc] .EOS
  * Regge and Wheeler [1957] Regge, T., and J. A. Wheeler (1957), Phys. Rev. 108, 1063.EOS
  * Rezzolla and Takami [2016] Rezzolla, L., and K. Takami (2016), Phys. Rev. D 93 (12), 124051, arXiv:1604.00246 [gr-qc] .EOS
  * de Rham [2014] de Rham, C. (2014), Living Rev. Rel. 17, 7, arXiv:1401.4173 [hep-th] .EOS
  * de Rham  _et al._ [2011] de Rham, C., G. Gabadadze, and A. J. Tolley (2011), Phys. Rev. Lett. 106, 231101, arXiv:1011.1232 [hep-th] .EOS
  * Rhoades and Ruffini [1974] Rhoades, C. E., Jr., and R. Ruffini (1974), Phys. Rev. Lett. 32, 324.EOS
  * Ribeiro and Vanzella [2020] Ribeiro, C. C. H., and D. A. T. Vanzella (2020), Phys. Rev. Res. 2 (1), 013281, arXiv:1912.01971 [gr-qc] .EOS
  * Riley  _et al._ [2019] Riley, T. E., _et al._ (2019), Astrophys. J. Lett. 887 (1), L21, arXiv:1912.05702 [astro-ph.HE] .EOS
  * Ripley [2022] Ripley, J. L. (2022), “Numerical relativity for Horndeski gravity,” arXiv:2207.13074 [gr-qc] .EOS
  * Ripley and Pretorius [2020] Ripley, J. L., and F. Pretorius (2020), Class. Quant. Grav. 37 (15), 155003, arXiv:2005.05417 [gr-qc] .EOS
  * Rosca-Mead  _et al._ [2019] Rosca-Mead, R., C. J. Moore, M. Agathos, and U. Sperhake (2019), Class. Quant. Grav. 36 (13), 134003, arXiv:1903.09704 [gr-qc] .EOS
  * Rosca-Mead  _et al._ [2020a] Rosca-Mead, R., C. J. Moore, U. Sperhake, M. Agathos, and D. Gerosa (2020a), Symmetry 12 (9), 1384, arXiv:2007.14429 [gr-qc] .EOS
  * Rosca-Mead  _et al._ [2020b] Rosca-Mead, R., U. Sperhake, C. J. Moore, M. Agathos, D. Gerosa, and C. D. Ott (2020b), Phys. Rev. D 102 (4), 044010, arXiv:2005.09728 [gr-qc] .EOS
  * Ruegg and Ruiz-Altaba [2004] Ruegg, H., and M. Ruiz-Altaba (2004), Int. J. Mod. Phys. A 19, 3265, arXiv:hep-th/0304245 .EOS
  * Ruiz  _et al._ [2012] Ruiz, M., J. C. Degollado, M. Alcubierre, D. Nunez, and M. Salgado (2012), Phys. Rev. D 86, 104044, arXiv:1207.6142 [gr-qc] .EOS
  * Saffer  _et al._ [2019] Saffer, A., H. O. Silva, and N. Yunes (2019), Phys. Rev. D 100 (4), 044030, arXiv:1903.07779 [gr-qc] .EOS
  * Salgado [2006] Salgado, M. (2006), Class. Quant. Grav. 23, 4719, arXiv:gr-qc/0509001 .EOS
  * Salgado  _et al._ [2008] Salgado, M., D. Martinez-del Rio, M. Alcubierre, and D. Nunez (2008), Phys. Rev. D 77, 104010, arXiv:0801.2372 [gr-qc] .EOS
  * Salgado  _et al._ [1998] Salgado, M., D. Sudarsky, and U. Nucamendi (1998), Phys. Rev. D 58, 124003, arXiv:gr-qc/9806070 .EOS
  * Salgado  _et al._ [2004] Salgado, M., D. Sudarsky, and U. Nucamendi (2004), Phys. Rev. D 70, 084027, arXiv:gr-qc/0402126 .EOS
  * Sampson  _et al._ [2014] Sampson, L., N. Yunes, N. Cornish, M. Ponce, E. Barausse, A. Klein, C. Palenzuela, and L. Lehner (2014), Phys. Rev. D 90 (12), 124091, arXiv:1407.7038 [gr-qc] .EOS
  * Sanchis-Gual  _et al._ [2016a] Sanchis-Gual, N., J. C. Degollado, C. Herdeiro, J. A. Font, and P. J. Montero (2016a), Phys. Rev. D 94 (4), 044061, arXiv:1607.06304 [gr-qc] .EOS
  * Sanchis-Gual  _et al._ [2016b] Sanchis-Gual, N., J. C. Degollado, P. J. Montero, J. A. Font, and C. Herdeiro (2016b), Phys. Rev. Lett. 116 (14), 141101, arXiv:1512.05358 [gr-qc] .EOS
  * Santiago  _et al._ [2016] Santiago, J., A. G. S. Landulfo, W. C. C. Lima, G. E. A. Matsas, R. F. P. Mendes, and D. A. T. Vanzella (2016), Phys. Rev. D 93 (2), 024043, arXiv:1512.02120 [gr-qc] .EOS
  * Sathyaprakash  _et al._ [2019] Sathyaprakash, B. S., _et al._ (2019), “Extreme Gravity and Fundamental Physics,” arXiv:1903.09221 [astro-ph.HE] .EOS
  * Scheel  _et al._ [1995a] Scheel, M. A., S. L. Shapiro, and S. A. Teukolsky (1995a), Phys. Rev. D 51, 4208, arXiv:gr-qc/9411025 .EOS
  * Scheel  _et al._ [1995b] Scheel, M. A., S. L. Shapiro, and S. A. Teukolsky (1995b), Phys. Rev. D 51, 4236, arXiv:gr-qc/9411026 .EOS
  * Sen and Banerjee [2001] Sen, S., and N. Banerjee (2001), Pramana 56, 487, arXiv:gr-qc/9809064 .EOS
  * Sennett and Buonanno [2016] Sennett, N., and A. Buonanno (2016), Phys. Rev. D 93 (12), 124004, arXiv:1603.03300 [gr-qc] .EOS
  * Sennett  _et al._ [2017] Sennett, N., L. Shao, and J. Steinhoff (2017), Phys. Rev. D 96 (8), 084019, arXiv:1708.08285 [gr-qc] .EOS
  * Shao  _et al._ [2017] Shao, L., N. Sennett, A. Buonanno, M. Kramer, and N. Wex (2017), Phys. Rev. X 7 (4), 041025, arXiv:1704.07561 [gr-qc] .EOS
  * Shibata  _et al._ [2014] Shibata, M., K. Taniguchi, H. Okawa, and A. Buonanno (2014), Phys. Rev. D 89 (8), 084005, arXiv:1310.0627 [gr-qc] .EOS
  * Shiralilou  _et al._ [2021] Shiralilou, B., T. Hinderer, S. Nissanke, N. Ortiz, and H. Witek (2021), Phys. Rev. D 103 (12), L121503, arXiv:2012.09162 [gr-qc] .EOS
  * Shiralilou  _et al._ [2022] Shiralilou, B., T. Hinderer, S. M. Nissanke, N. Ortiz, and H. Witek (2022), Class. Quant. Grav. 39 (3), 035002, arXiv:2105.13972 [gr-qc] .EOS
  * Silva  _et al._ [2022] Silva, H. O., A. Coates, F. M. Ramazanoğlu, and T. P. Sotiriou (2022), Phys. Rev. D 105 (2), 024046, arXiv:2110.04594 [gr-qc] .EOS
  * Silva  _et al._ [2021a] Silva, H. O., A. M. Holgado, A. Cárdenas-Avendaño, and N. Yunes (2021a), Phys. Rev. Lett. 126 (18), 181101, arXiv:2004.01253 [gr-qc] .EOS
  * Silva  _et al._ [2015] Silva, H. O., C. F. Macedo, E. Berti, and L. C. Crispino (2015), Class. Quant. Grav. 32, 145008, arXiv:1411.6286 [gr-qc] .EOS
  * Silva  _et al._ [2019] Silva, H. O., C. F. B. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein, and E. Berti (2019), Phys. Rev. D 99 (6), 064011, arXiv:1812.05590 [gr-qc] .EOS
  * Silva  _et al._ [2018] Silva, H. O., J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti (2018), Phys. Rev. Lett. 120 (13), 131104, arXiv:1711.02080 [gr-qc] .EOS
  * Silva  _et al._ [2014] Silva, H. O., H. Sotani, E. Berti, and M. Horbatsch (2014), Phys. Rev. D 90 (12), 124044, arXiv:1410.2511 [gr-qc] .EOS
  * Silva  _et al._ [2021b] Silva, H. O., H. Witek, M. Elley, and N. Yunes (2021b), Phys. Rev. Lett. 127 (3), 031101, arXiv:2012.10436 [gr-qc] .EOS
  * Silva and Yunes [2019a] Silva, H. O., and N. Yunes (2019a), Class. Quant. Grav. 36 (17), 17LT01, arXiv:1902.10269 [gr-qc] .EOS
  * Silva and Yunes [2019b] Silva, H. O., and N. Yunes (2019b), Phys. Rev. D 99 (4), 044034, arXiv:1808.04391 [gr-qc] .EOS
  * Soldateschi  _et al._ [2020] Soldateschi, J., N. Bucciantini, and L. Del Zanna (2020), Astron. Astrophys. 640, A44, arXiv:2005.12758 [astro-ph.HE] .EOS
  * Soldateschi  _et al._ [2021] Soldateschi, J., N. Bucciantini, and L. Del Zanna (2021), Astron. Astrophys. 645, A39, arXiv:2010.14833 [astro-ph.HE] .EOS
  * Sotani [2012] Sotani, H. (2012), Phys. Rev. D 86, 124036, arXiv:1211.6986 [astro-ph.HE] .EOS
  * Sotani [2014] Sotani, H. (2014), Phys. Rev. D 89 (6), 064031, arXiv:1402.5699 [astro-ph.HE] .EOS
  * Sotani [2017] Sotani, H. (2017), Phys. Rev. D 96 (10), 104010, arXiv:1710.10596 [astro-ph.HE] .EOS
  * Sotani and Kokkotas [2004] Sotani, H., and K. D. Kokkotas (2004), Phys. Rev. D 70, 084026, arXiv:gr-qc/0409066 .EOS
  * Sotani and Kokkotas [2005] Sotani, H., and K. D. Kokkotas (2005), Phys. Rev. D 71, 124038, arXiv:gr-qc/0506060 .EOS
  * Sotani and Kokkotas [2017] Sotani, H., and K. D. Kokkotas (2017), Phys. Rev. D 95 (4), 044032, arXiv:1702.00874 [gr-qc] .EOS
  * Sotiriou and Faraoni [2012] Sotiriou, T. P., and V. Faraoni (2012), Phys. Rev. Lett. 108, 081103, arXiv:1109.6324 [gr-qc] .EOS
  * Sotiriou and Zhou [2014a] Sotiriou, T. P., and S.-Y. Zhou (2014a), Phys. Rev. Lett. 112, 251102, arXiv:1312.3622 [gr-qc] .EOS
  * Sotiriou and Zhou [2014b] Sotiriou, T. P., and S.-Y. Zhou (2014b), Phys. Rev. D 90, 124063, arXiv:1408.1698 [gr-qc] .EOS
  * Sperhake  _et al._ [2017] Sperhake, U., C. J. Moore, R. Rosca, M. Agathos, D. Gerosa, and C. D. Ott (2017), Phys. Rev. Lett. 119 (20), 201103, arXiv:1708.03651 [gr-qc] .EOS
  * Staykov and Doneva [2022] Staykov, K. V., and D. D. Doneva (2022), “Multi-scalar Gauss-Bonnet gravity: scalarized black holes beyond spontaneous scalarization,” arXiv:2209.01038 [gr-qc] .EOS
  * Staykov  _et al._ [2019] Staykov, K. V., D. D. Doneva, and S. S. Yazadjiev (2019), Astrophys. Space Sci. 364 (10), 178, arXiv:1902.09208 [gr-qc] .EOS
  * Staykov  _et al._ [2015] Staykov, K. V., D. D. Doneva, S. S. Yazadjiev, and K. D. Kokkotas (2015), Phys. Rev. D 92 (4), 043009, arXiv:1503.04711 [gr-qc] .EOS
  * Staykov  _et al._ [2018] Staykov, K. V., D. Popchev, D. D. Doneva, and S. S. Yazadjiev (2018), Eur. Phys. J. C 78 (7), 586, arXiv:1805.07818 [gr-qc] .EOS
  * Stefanov  _et al._ [2007a] Stefanov, I. Z., S. S. Yazadjiev, and M. D. Todorov (2007a), Phys. Rev. D 75, 084036, arXiv:0704.3784 [gr-qc] .EOS
  * Stefanov  _et al._ [2007b] Stefanov, I. Z., S. S. Yazadjiev, and M. D. Todorov (2007b), Mod. Phys. Lett. A 22, 1217, arXiv:0708.3203 [gr-qc] .EOS
  * Stefanov  _et al._ [2008] Stefanov, I. Z., S. S. Yazadjiev, and M. D. Todorov (2008), Mod. Phys. Lett. A 23, 2915, arXiv:0708.4141 [gr-qc] .EOS
  * Strohmayer and Watts [2005] Strohmayer, T. E., and A. L. Watts (2005), Astrophys. J. Lett. 632, L111, arXiv:astro-ph/0508206 .EOS
  * Strohmayer and Watts [2006] Strohmayer, T. E., and A. L. Watts (2006), Astrophys. J. 653, 593, arXiv:astro-ph/0608463 .EOS
  * Takami  _et al._ [2014] Takami, K., L. Rezzolla, and L. Baiotti (2014), Phys. Rev. Lett. 113 (9), 091104, arXiv:1403.5672 [gr-qc] .EOS
  * Takami  _et al._ [2015] Takami, K., L. Rezzolla, and L. Baiotti (2015), Phys. Rev. D 91 (6), 064001, arXiv:1412.3240 [gr-qc] .EOS
  * Taniguchi and Gourgoulhon [2003] Taniguchi, K., and E. Gourgoulhon (2003), Phys. Rev. D 68, 124025, arXiv:gr-qc/0309045 .EOS
  * Taniguchi and Shibata [2010] Taniguchi, K., and M. Shibata (2010), Astrophys. J. Suppl. 188, 187, arXiv:1005.0958 [astro-ph.SR] .EOS
  * Taniguchi  _et al._ [2015] Taniguchi, K., M. Shibata, and A. Buonanno (2015), Phys. Rev. D 91 (2), 024033, arXiv:1410.0738 [gr-qc] .EOS
  * Tasinato [2014] Tasinato, G. (2014), JHEP 04, 067, arXiv:1402.6450 [hep-th] .EOS
  * Taylor and Weisberg [1982] Taylor, J. H., and J. M. Weisberg (1982), Astrophys. J. 253, 908.EOS
  * Taylor  _et al._ [1993] Taylor, J. N., A. Wolszczan, and T. Damour (1993), Nature 355, 132.EOS
  * Thorne and Campolattaro [1967] Thorne, K. S., and A. Campolattaro (1967), Astrophysical Journal , 591.EOS
  * Tuna  _et al._ [2022] Tuna, S., K. I. Ünlütürk, and F. M. Ramazanoğlu (2022), Phys. Rev. D 105 (12), 124070, arXiv:2204.02138 [gr-qc] .EOS
  * Ventagli  _et al._ [2021] Ventagli, G., G. Antoniou, A. Lehébel, and T. P. Sotiriou (2021), Phys. Rev. D 104 (12), 124078, arXiv:2111.03644 [gr-qc] .EOS
  * Ventagli  _et al._ [2020] Ventagli, G., A. Lehébel, and T. P. Sotiriou (2020), Phys. Rev. D 102 (2), 024050, arXiv:2006.01153 [gr-qc] .EOS
  * Voisin  _et al._ [2020] Voisin, G., I. Cognard, P. C. C. Freire, N. Wex, L. Guillemot, G. Desvignes, M. Kramer, and G. Theureau (2020), Astron. Astrophys. 638, A24, arXiv:2005.01388 [gr-qc] .EOS
  * Wald [1993] Wald, R. M. (1993), Phys. Rev. D 48 (8), R3427, arXiv:gr-qc/9307038 .EOS
  * Watts  _et al._ [2016] Watts, A. L., _et al._ (2016), Rev. Mod. Phys. 88 (2), 021001, arXiv:1602.01081 [astro-ph.HE] .EOS
  * Weih  _et al._ [2018] Weih, L. R., E. R. Most, and L. Rezzolla (2018), Mon. Not. Roy. Astron. Soc. 473 (1), L126, arXiv:1709.06058 [gr-qc] .EOS
  * Wex [2014] Wex, N. (2014), “Testing Relativistic Gravity with Radio Pulsars,” arXiv:1402.5594 [gr-qc] .EOS
  * Whinnett [2000] Whinnett, A. (2000), Phys. Rev. D 61, 124014, arXiv:gr-qc/9911052 .EOS
  * Whinnett [1999] Whinnett, A. W. (1999), Class. Quant. Grav. 16, 2797.EOS
  * Whinnett and Torres [2004] Whinnett, A. W., and D. F. Torres (2004), Astrophys. J. Lett. 603, L133, arXiv:astro-ph/0401521 .EOS
  * Will [2014] Will, C. M. (2014), Living Rev. Rel. 17, 4, arXiv:1403.7377 [gr-qc] .EOS
  * Will [2018] Will, C. M. (2018),  _Theory and Experiment in Gravitational Physics_ (Cambridge University Press).EOS
  * Witek  _et al._ [2019] Witek, H., L. Gualtieri, P. Pani, and T. P. Sotiriou (2019), Phys. Rev. D 99 (6), 064035, arXiv:1810.05177 [gr-qc] .EOS
  * Wong  _et al._ [2022] Wong, L. K., C. A. R. Herdeiro, and E. Radu (2022), Phys. Rev. D 106 (2), 024008, arXiv:2204.09038 [gr-qc] .EOS
  * Xu  _et al._ [2020] Xu, R., Y. Gao, and L. Shao (2020), Phys. Rev. D 102 (6), 064057, arXiv:2007.10080 [gr-qc] .EOS
  * Yagi  _et al._ [2014] Yagi, K., K. Kyutoku, G. Pappas, N. Yunes, and T. A. Apostolatos (2014), Phys. Rev. D 89 (12), 124013, arXiv:1403.6243 [gr-qc] .EOS
  * Yagi  _et al._ [2012] Yagi, K., L. C. Stein, N. Yunes, and T. Tanaka (2012), Phys. Rev. D 85, 064022, [Erratum: Phys.Rev.D 93, 029902 (2016)], arXiv:1110.5950 [gr-qc] .EOS
  * Yagi and Stepniczka [2021] Yagi, K., and M. Stepniczka (2021), Phys. Rev. D 104 (4), 044017, arXiv:2105.01614 [gr-qc] .EOS
  * Yagi and Yunes [2013a] Yagi, K., and N. Yunes (2013a), Science 341, 365, arXiv:1302.4499 [gr-qc] .EOS
  * Yagi and Yunes [2013b] Yagi, K., and N. Yunes (2013b), Phys. Rev. D 88 (2), 023009, arXiv:1303.1528 [gr-qc] .EOS
  * Yagi and Yunes [2017] Yagi, K., and N. Yunes (2017), Phys. Rept. 681, 1, arXiv:1608.02582 [gr-qc] .EOS
  * Yamashita  _et al._ [2014] Yamashita, Y., A. De Felice, and T. Tanaka (2014), Int. J. Mod. Phys. D 23, 1443003, arXiv:1408.0487 [hep-th] .EOS
  * Yazadjiev [1999] Yazadjiev, S. (1999), Class. Quant. Grav. 16, L63, arXiv:gr-qc/9906038 .EOS
  * Yazadjiev  _et al._ [2017] Yazadjiev, S. S., D. D. Doneva, and K. D. Kokkotas (2017), Phys. Rev. D 96 (6), 064002, arXiv:1705.06984 [gr-qc] .EOS
  * Yazadjiev  _et al._ [2016] Yazadjiev, S. S., D. D. Doneva, and D. Popchev (2016), Phys. Rev. D 93 (8), 084038, arXiv:1602.04766 [gr-qc] .EOS
  * Yunes and Pretorius [2009] Yunes, N., and F. Pretorius (2009), Phys. Rev. D 79, 084043, arXiv:0902.4669 [gr-qc] .EOS
  * Yunes and Siemens [2013] Yunes, N., and X. Siemens (2013), Living Rev. Rel. 16, 9, arXiv:1304.3473 [gr-qc] .EOS
  * Yunes and Stein [2011] Yunes, N., and L. C. Stein (2011), Phys. Rev. D 83, 104002, arXiv:1101.2921 [gr-qc] .EOS
  * Zenginoglu [2008] Zenginoglu, A. (2008), Class. Quant. Grav. 25, 145002, arXiv:0712.4333 [gr-qc] .EOS
  * Zhang  _et al._ [2020] Zhang, S.-J., B. Wang, A. Wang, and J. F. Saavedra (2020), Phys. Rev. D 102 (12), 124056, arXiv:2010.05092 [gr-qc] .EOS
  * Zhao  _et al._ [2022] Zhao, J., P. C. C. Freire, M. Kramer, L. Shao, and N. Wex (2022), Class. Quant. Grav. 39 (11), 11LT01, arXiv:2201.03771 [astro-ph.HE] .EOS
  * Zumalacárregui and García-Bellido [2014] Zumalacárregui, M., and J. García-Bellido (2014), Phys. Rev. D 89, 064046, arXiv:1308.4685 [gr-qc] .EOS