# Neutron stars in massive scalar-Gauss-Bonnet gravity: Spherical structure
and time-independent perturbations

Rui Xu<EMAIL_ADDRESS>Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Yong Gao  Department of Astronomy, School
of Physics, Peking University, Beijing 100871, China  Kavli Institute for
Astronomy and Astrophysics, Peking University, Beijing 100871, China  Lijing
Shao<EMAIL_ADDRESS>Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China

###### Abstract

The class of scalar-tensor theories with the scalar field coupling to the
Gauss-Bonnet invariant has drawn great interest since solutions of spontaneous
scalarization were found for black holes in these theories. We contribute to
the existing literature a detailed study of the spontaneously scalarized
neutron stars (NSs) in a typical theory where the coupling function of the
scalar field takes the quadratic form and the scalar field is massive. The
investigation here includes the spherical solutions of the NSs as well as
their perturbative properties, namely the tidal deformability and the moment
of inertia, treated in a unified and extendable way under the framework of
spherical decomposition. We find that while the mass, the radius, and the
moment of inertia of the spontaneously scalarized NSs show very moderate
deviations from those of the NSs in general relativity (GR), the tidal
deformability exhibits significant differences between the solutions in GR and
the solutions of spontaneous scalarization for certain values of the
parameters in the scalar-Gauss-Bonnet theory. As a result, the celebrated
universal relation between the moment of inertia and the tidal deformability
of neutron stars breaks down. With the mass and the tidal deformability of NSs
attainable in the gravitational waves from binary NS mergers, the radius
measurable using the X-ray satellites, and the moment of inertia accessible
via the high-precision pulsar timing techniques, future multi-messenger
observations can be contrasted with the theoretical results and provide us
necessary information for building up theories beyond GR.

##  I Introduction

One aspect of the elegance in Einstein’s general relativity (GR) is that the
geometric part of the field equations comes from the Einstein-Hilbert action,
which merely consists of the simplest curvature invariant of the spacetime
manifold, namely, the Ricci scalar. This simplicity is sacrificed in many
alternative theories of gravity, where other curvature invariants are taken in
to give modifications to GR, attempting to resolve its fundamental
incompatibility with quantum field theory Stelle (1977), or to address issues
in cosmology including the initial singularity of the Big Bang and the mystery
of dark energy Starobinsky (1980); Copeland  _et al._ (2006); Nojiri and
Odintsov (2011). Among the quadratic curvature invariants, the Gauss-Bonnet
term

$\displaystyle{\mathcal{G}}=-G_{\alpha\beta\gamma\delta}R^{\alpha\beta\gamma\delta}=R_{{\alpha\beta\gamma\delta}}R^{{\alpha\beta\gamma\delta}}-4R_{\alpha\beta}R^{\alpha\beta}+R^{2},$
(1)

is special in the sense that adding it multiplied by a constant to the action
does not alter the field equations. In the above expression, the double dual
of the Riemann curvature tensor is defined as

$\displaystyle
G_{\alpha\beta\gamma\delta}:=\frac{1}{4}E_{\alpha\beta\kappa\lambda}E_{\gamma\delta\rho\sigma}R^{\kappa\lambda\rho\sigma},$
(2)

with $E_{{\alpha\beta\gamma\delta}}$ being the Levi-Civita tensor. Here the
convention adopted for the Levi-Civita symbol
$\epsilon_{{\alpha\beta\gamma\delta}}$ is $\epsilon_{0123}=1$.

This special property of the Gauss-Bonnet term is only true in the four-
dimensional spacetime where the expression
$g_{\mu\nu}E_{{\alpha\beta\gamma\delta}}$, when antisymmetrized over any five
indices, vanishes. Then, the identity

$\displaystyle
g_{\mu\nu}{\mathcal{G}}=-4G_{\mu\beta\gamma\delta}R_{\nu}^{\phantom{\mu}\beta\gamma\delta},$
(3)

can be proved, and the variation of ${\mathcal{G}}$ with respect to the metric
can be calculated to obtain

$\displaystyle\delta\left(\sqrt{-g}\,{\mathcal{G}}\right)=-4\sqrt{-g}\,D_{\gamma}\left(G_{\alpha}^{\phantom{\alpha}\beta\gamma\delta}\delta\Gamma^{\alpha}_{\phantom{\alpha}\beta\delta}\right),$
(4)

where $g$ is the determinant of the metric $g_{\mu\nu}$, and
$\Gamma^{\alpha}_{\phantom{\alpha}\beta\delta}$ are the Christoffel symbols
with $D_{\alpha}$ being the covariant derivative associated with them. The
above variation contributes as a surface term in the action, and hence is
dropped in the equations of motion Lanczos (1938).

The futility of adding ${\mathcal{G}}$ multiplied by a constant to the action
did not stop the interest in the Gauss-Bonnet term. On the contrary, it
stimulated further investigations in dimensions higher than four Cai and Guo
(2004); Konoplya (2005); Cai  _et al._ (2005); Camanho and Edelstein (2010),
motivated by string theory Zwiebach (1985); Boulware and Deser (1985).
Remarkably, a nontrivial theory in the four-dimensional spacetime was recently
established by setting the constant in front of ${\mathcal{G}}$ proportional
to $1/(D-4)$ with $D$ being the dimension of the spacetime Glavan and Lin
(2020). When taking the limit $D\rightarrow 4$, the Gauss-Bonnet term gives a
finite contribution to the Einstein field equations and thus modifies GR.

Another direction to investigate possible modifications to GR generated by the
Gauss-Bonnet term is introducing a scalar field that couples to it Torii and
Maeda (2005); Nojiri  _et al._ (2005); Boulware and Deser (1986); Koivisto and
Mota (2007); Cognola  _et al._ (2007); Leith and Neupane (2007); Akbarieh  _et
al._ (2021). The idea also originates from string theory where a scalar field,
called the dilaton, appears when the extra dimensions are compactified to
recover the metric tensor of the four-dimensional spacetime Scherk and Schwarz
(1974). One well-studied case is when the Gauss-Bonnet term couples to an
exponential function of the scalar field in the action, and inherits the name
dilatonic Einstein-Gauss-Bonnet theory Kanti  _et al._ (1996); Torii  _et al._
(1997); Kanti  _et al._ (1998); Torii and Maeda (1998); Pani and Cardoso
(2009); Kanti  _et al._ (2011); Ayzenberg and Yunes (2014); Cunha  _et al._
(2017); Ayzenberg and Yunes (2018); Gonzalez  _et al._ (2019); Blázquez-
Salcedo  _et al._ (2016); Konoplya  _et al._ (2019); Pierini and Gualtieri
(2021); Yagi  _et al._ (2016); Witek  _et al._ (2019); Wang  _et al._ (2021);
another basic case has the scalar field directly coupled to the Gauss-Bonnet
term so that the field equations are invariant under a constant shift of the
scalar field, and is given the name shift-symmetric Gauss-Bonnet theory
Sotiriou and Zhou (2014); Saffer and Yagi (2021). In these theories, the
variation of the term involving ${\mathcal{G}}$ no long contributes as a
surface term with the help of a nontrivial scalar field, and therefore
modifies the Einstein field equations. Notice that the scalar field equations
in these theories do not admit a constant scalar as the solution. So the
metric solutions in these theories are generally different from that in GR. In
this work, we are interested in another type of theories where a phenomenon
called “spontaneous scalarization” can happen Damour and Esposito-Farèse
(1993, 1996); Harada (1998); Salgado  _et al._ (1998); Ramazanoğlu and
Pretorius (2016); Mendes and Ortiz (2016); Herdeiro  _et al._ (2018); Cunha
_et al._ (2019); Antoniou  _et al._ (2018); Doneva and Yazadjiev (2018a);
Silva  _et al._ (2018); Doneva and Yazadjiev (2018b); Kuan  _et al._ (2021).
This requires the theory to admit the trivial GR solution where the scalar
field is a constant by assigning an even function of the scalar field coupled
to ${\mathcal{G}}$ in the action. The theories in Refs. Antoniou  _et al._
(2018); Doneva and Yazadjiev (2018a); Silva  _et al._ (2018); Doneva and
Yazadjiev (2018b); Kuan  _et al._ (2021); Antoniou  _et al._ (2021); Ventagli
_et al._ (2021), as proved therein, are such theories having solutions of
spontaneous scalarization for both black holes and neutron stars (NSs). Here
we refer them as scalar-Gauss-Bonnet theories of gravity with spontaneous
scalarization, or simply scalar-Gauss-Bonnet theories.

Theories in which compact objects can be scalarized are compelling, because
they make identical predictions as those of GR in the weak-field regime,
guaranteeing their validity as alternative theories of gravity, but are ready
to be distinguished from GR via astrophysical observations of strong-field
systems Damour and Esposito-Farèse (1996); Wex (2014); Shao and Wex (2016);
Shao  _et al._ (2017). These observations include, for example, the shadow
image of the supermassive black hole in the galaxy M87 obtained by the Event
Horizon Telescope (EHT) Akiyama  _et al._ (2019), the orbital decay of binary
pulsars observed by pulsar timing arrays Hulse and Taylor (1975), and the
gravitational waves (GWs) from coalescences of compact binary systems detected
by the LIGO and Virgo detectors Abbott  _et al._ (2019a, 2021). The first
comparison between the M87 shadow image and the numerically calculated black
hole shadows in the typical scalar-Gauss-Bonnet theory proposed in Ref. Doneva
and Yazadjiev (2018b) has been done by Cunha  _et al._ (2019), where a
preliminary constraint on the parameter of the theory is obtained. A recent
work by Danchev  _et al._ (2021) first constrains the same theory using the
orbital data of three binary pulsars under the framework of Bayesian analysis
and MCMC strategy Shao  _et al._ (2017). As for using GW data to test scalar-
Gauss-Bonnet theories, though it seems still early to talk about numerical
bounds, the field is actively developing toward it. Particularly, Shiralilou
_et al._ (2021a, b) recently derives the post-Newtonian solution in scalar-
Gauss-Bonnet theories to the first order beyond the quadrupole approximation,
providing waveforms for the inspiralling stage of compact binary systems.
Numerical simulations that required to build the waveforms for the merger
stage of compact binary systems are unavailable but the investigation by
Kovács and Reall (2020) on the well-posedness of the initial-value problem in
scalar-Gauss-Bonnet theories has set the first step in the direction. The
ringdown stage of a coalescence can be approximated using perturbation theory
in a single-body background spacetime, and the quasinormal modes (QNMs) of
nonrotating black holes, which might serve as the basis for constructing the
GWs emitted by the remnant black hole, have been studied in scalar-Gauss-
Bonnet theories by Bryant  _et al._ (2021) and Staykov  _et al._ (2021). The
QNMs of NSs, though in scalar-tensor theories without the Gauss-Bonnet term,
have been studied by Mendes and Ortiz (2018).

Our work aims at giving a detailed account of the spontaneously scalarized NSs
in a typical scalar-Gauss-Bonnet theory with a massive scalar field, which
complements the existing literature on spontaneously scalarized black holes as
well as the work by Doneva and Yazadjiev (2018b) where NS solutions of
spontaneous scalarization are obtained in a massless theory of such kind. In
this paper, the study includes an analysis of the linearized scalar equation
to find the onsets of spontaneous scalarization, illustrations of the
numerical solutions to the group of nonlinear spherical equations, and a
unified treatment of the slow-rotation and the tidal perturbations using the
formalism of spherical decomposition introduced by Regge and Wheeler (1957)
and developed for NSs in Refs. Thorne and Campolattaro (1967); Price and
Thorne (1969); Thorne (1969). Our numerical solutions of the scalarized NSs
are ready to use for constructing x-ray pulsar pulse profiles via the similar
process demonstrated in Refs. Silva and Yunes (2019); Xu  _et al._ (2020); Hu
_et al._ (2021) when adequate observational data, for example, from the NS
Interior Composition Explorer (NICER) mission, are collected and analyzed to
be compared with the theoretical predictions. The x-ray pulsar observations
are also expected to measure the moment of inertia for NSs Hu  _et al._
(2020); Miao and Li (2021), providing another avenue to constrain the theory
when compared to the numerical results for the moment of inertia calculated
here. The tidal deformability, extracted from GW observations of NS mergers
like GW170817 Abbott  _et al._ (2017), provides a third way to test the
theory, and certain values for the parameters of the theory have been excluded
by a comparison between the measurement from GW170817 and the numerical
results calculated here.

The paper is organized as follows. We start with a general setup for spherical
NSs in scalar-Gauss-Bonnet theories in Sec. II.1, and then numerically solve
the linearized scalar equation in Sec. II.2 and the group of nonlinear
equations in Sec. II.3 for the typical scalar-Gauss-Bonnet theory with a
quadratic coupling function. Section III deals with time-independent
perturbations of the scalarized NSs. It starts with a brief review of the
perturbation theory using spherical decomposition in Sec. III.1. Then
numerical results for the tidal deformability and the moment of inertia of the
scalarized NSs are presented in Sec. III.2 and Sec. III.3. Finally,
conclusions drawn from our study are summarized in Sec. IV. The appendices
provide detailed equations and analyses in our study. Appendix A is about the
expansion of the scalarized solution at the center of the star. Appendix B is
for the expansion of the scalarized solution at infinity. Appendix C displays
the radial equations for the even-parity perturbations with the angular
momentum number $l\geq 2$. Appendix D discusses a numerical issue in
calculating the tidal deformability of the scalarized NSs. Appendix E displays
the radial equations for the odd-parity perturbations.

Throughout the work, equations are written in the geometrized unit system
where $G=c=1$ unless conventional units appear in numerical results for the
intuitive perception of the physical quantities. The sign convention of the
metric is $(-,+,+,+)$.

##  II spherical NSs in scalar-Gauss-Bonnet gravity

###  II.1 General setup

The action for the scalar-Gauss-Bonnet theory of gravitation to be studied has
the general form Antoniou  _et al._ (2018); Doneva and Yazadjiev (2018a);
Silva  _et al._ (2018)

$\displaystyle S$ $\displaystyle=$ $\displaystyle\frac{1}{16\pi}\int
d^{4}x\sqrt{-g}\left(R-g^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi-U(\Phi)+F(\Phi){\mathcal{G}}\right)$
(5) $\displaystyle+S_{m},$

where the scalar field nonminimally couples to the metric through the product
of the coupling function $F$ and the Gauss-Bonnet term ${\mathcal{G}}$. We
assume that the scalar field $\Phi$ is absent in the matter action $S_{m}$.
The field equations

$\displaystyle G_{\mu\nu}$ $\displaystyle=$ $\displaystyle
8\pi\left(T_{m}\right)_{\mu\nu}+\left(T_{\Phi}\right)_{\mu\nu}+\left(T_{\rm
GB}\right)_{\mu\nu},$ $\displaystyle D_{\alpha}D^{\alpha}\Phi$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{dF}{d\Phi}{\mathcal{G}}+\frac{1}{2}\frac{dU}{d\Phi},$
(6)

can be obtained by taking variations with respect to the metric and the scalar
field. In the above equations, $G_{\mu\nu}$ is the Einstein tensor,
$\left(T_{m}\right)_{\mu\nu}$ is the energy-momentum tensor for conventional
matter, $\left(T_{\Phi}\right)_{\mu\nu}$ is the energy-momentum tensor of a
minimally coupled scalar

$\displaystyle\left(T_{\Phi}\right)_{\mu\nu}:=\partial_{\mu}\Phi\partial_{\nu}\Phi-\frac{1}{2}g_{\mu\nu}\left(g^{\alpha\beta}\partial_{\alpha}\Phi\partial_{\beta}\Phi+U\right),$
(7)

and $\left(T_{\rm GB}\right)_{\mu\nu}$ representing the contribution from the
Gauss-Bonnet term is

$\displaystyle\left(T_{\rm
GB}\right)_{\mu\nu}=4G_{\mu\alpha\nu\beta}D^{\alpha}D^{\beta}F,$ (8)

where $G_{\alpha\beta\gamma\delta}$ defined in Eq. (2) is the double dual of
the Riemann tensor while $D_{\alpha}$ is the covariant derivative.

Before we proceed to write down the ordinary differential equations (ODEs)
using the spherical ansatz, let us mention that though the energy-momentum
conservation equation for conventional matter

$\displaystyle D_{\alpha}\left(T_{m}\right)^{\alpha\beta}=0,$ (9)

is merely a consequence of the field equations,111One can prove
$D_{\alpha}\left(T_{\Phi}\right)^{\alpha\beta}+D_{\alpha}\left(T_{\rm
GB}\right)^{\alpha\beta}=0$ with the help of the scalar field equation.
components of it turn out to be handy to manipulate compared to the field
equations themselves.

Now let us take the static spherical ansatz

$\displaystyle
ds^{2}=-e^{2\nu}dt^{2}+e^{2\mu}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right),$ (10)

and the perfect-fluid assumption for the conventional matter

$\displaystyle\left(T_{m}\right)_{\alpha\beta}=(\epsilon+p)u_{\alpha}u_{\beta}+pg_{\alpha\beta},$
(11)

where the metric functions $\mu,\,\nu$, the fluid proper energy density
$\epsilon$, and the fluid proper pressure $p$ are all functions of the radial
coordinate $r$ only, while the 4-velocity field of the star is simply
$u^{\alpha}=(e^{-\nu},0,0,0)$. Then, the $tt$\- and $rr$-components of the
modified Einstein field equations, the scalar field equation, and the
$r$-component of the matter conservation equation conveniently give the group
of ODEs for $\mu,\,\nu,\,\Phi$ and $p$:

$\displaystyle\mu^{\prime}$ $\displaystyle=$
$\displaystyle\frac{-2(e^{2\mu}-1)+16\pi
r^{2}e^{2\mu}\epsilon+r^{2}e^{2\mu}U+r^{2}\Phi^{\prime\,2}}{4\left(r+2\left(1-3e^{-2\mu}\right)F^{\prime}\right)}$
$\displaystyle+\frac{8\left(1-e^{-2\mu}\right)F^{\prime\prime}}{4\left(r+2\left(1-3e^{-2\mu}\right)F^{\prime}\right)},$
$\displaystyle\nu^{\prime}$ $\displaystyle=$
$\displaystyle\frac{2(e^{2\mu}-1)+16\pi
r^{2}e^{2\mu}p-r^{2}e^{2\mu}U+r^{2}\Phi^{\prime\,2}}{4\left(r+2(1-3e^{-2\mu})F^{\prime}\right)},$
$\displaystyle\Phi^{\prime\prime}$ $\displaystyle=$
$\displaystyle\left(\mu^{\prime}-\nu^{\prime}-\frac{2}{r}\right)\Phi^{\prime}-\frac{1}{2}e^{2\mu}{\mathcal{G}}\frac{dF}{d\Phi}+\frac{1}{2}e^{2\mu}\frac{dU}{d\Phi},$
$\displaystyle p^{\prime}$ $\displaystyle=$
$\displaystyle-(\epsilon+p)\,\nu^{\prime},$ (12)

where the prime denotes the derivative with respect to $r$, and hence

$\displaystyle
F^{\prime}=\frac{dF}{d\Phi}\Phi^{\prime},~{}~{}U^{\prime}=\frac{dU}{d\Phi}\Phi^{\prime},~{}~{}F^{\prime\prime}=\frac{dF}{d\Phi}\Phi^{\prime\prime}+\frac{d^{2}F}{d\Phi^{2}}\Phi^{\prime\,2}.$
(13)

With the expression for the Gauss-Bonnet term

$\displaystyle{\mathcal{G}}=\frac{8e^{-2\mu}\left[\left(1-3e^{-2\mu}\right)\mu^{\prime}\nu^{\prime}-\left(1-e^{-2\mu}\right)\left(\nu^{\prime\prime}+\nu^{\prime\,2}\right)\right]}{r^{2}},$
(14)

and an equation of state (EOS) that expresses $\epsilon$ in terms of $p$, the
above ODEs can be cast into the standard form of $dy/dx=f(x,y)$ for numerical
integration.

To conduct the numerical integration, we take the scalar potential $U$ as

$\displaystyle
U(\Phi)=\left(\frac{2\pi}{\lambda_{\Phi}}\right)^{2}\Phi^{2}=\frac{1}{\lambdabar_{\Phi}^{2}}\Phi^{2},$
(15)

to represent a massive scalar field whose mass is $m_{\Phi}=h/\lambda_{\Phi}$
with $\lambda_{\Phi}$ being its Compton wavelength and $h$ being the Planck
constant. The coupling function $F$ in our work is assigned to be

$\displaystyle F(\Phi)=\xi\Phi^{2},$ (16)

with $\xi$ being a constant. It matches the coupling function in Ref. Doneva
and Yazadjiev (2018b)

$\displaystyle\lambda^{2}f(\Phi)=\pm\frac{\lambda^{2}}{2\beta}\left(e^{-\beta\Phi^{2}}-1\right),$
(17)

via $\xi=\mp\lambda^{2}/2$ when $\Phi$ is small.

In the rest of this section, we first use the method demonstrated by Xu  _et
al._ (2020) to find the values of, say, the central energy density of the
star, for the onsets of spontaneous scalarization. Then, full nonlinear
solutions of spontaneous scalarization will be found around these critical
values of the central energy density.

Figure 1: $\Phi^{\prime}/\Phi\big{|}_{s}$ with respect to the central energy
density $\epsilon_{0}$. The dashed curves are from the interior solutions
while the solid curves are from the exterior solutions. The value of $a$ is
fixed to $1$ and the values of $b$ are shown in the panels; $a$ and $b$ are
defined in Eqs. (19). The EOS used is AP4. Notice that there are no
intersections for $b=0.9$ and $b=-0.04$. Figure 2: The boundary of the
parameter space for spontaneous scalarization. The solid curves give the
values of $a$ and $b$ for spontaneous scalarization to happen at the highest
central energy density for a given EOS, while the dashed curves demand
spontaneous scalarization to happen at the central energy density
corresponding to the maximum mass of the NS for a given EOS.

###  II.2 Parameter space and onsets of spontaneous scalarization

Figure 3: Mass-radius relation for NSs. The dashed curves are GR results,
while the solid lines and the dots are the $M$-$R$ segments for the scalarized
NSs. The same five EOSs as in Fig. 2 are used: WFF1, SLy4, AP4, MPA1, and
PAL1. Note that there are no scalarization solutions when $a=1,\,b=-0.2$ for
the EOS MPA1 as its critical value $b_{-}$ at $a=1$ shown in Fig. 2 is
slightly less than $-0.2$.  Figure 4: The scalar field at the center of the
star versus the mass of the star for the solutions of spontaneous
scalarization in Fig. 3. To keep the plots clear, only the solutions in the
lowest mass band are plotted when there are more than one interval of
spontaneous scalarization along the $M$-$R$ curve. The same set of colors as
in Fig. 3 is used for the five EOSs.

The set of the ODEs in Eqs. (12) always admits the GR solutions for
$\mu,\,\nu$ and $p$ together with a trivial scalar field $\Phi=0$. The
solutions of spontaneous scalarization where $\Phi$ is nontrivial however have
requirements for both the star parameter which we take as the central energy
density $\epsilon_{0}$, and the theory parameters $m_{\Phi}$ and $\xi$. We
reveal the parameter space of $m_{\Phi}$ and $\xi$ for spontaneous
scalarization as well as the critical values of the central energy density of
the star for given values of $m_{\Phi}$ and $\xi$ by analyzing the linearized
scalar field equation.

The linearized scalar field equation is the scalar field equation in Eqs. (12)
with the metric functions $\mu$ and $\nu$ taking their GR solution. It is
solved in two regions: the interior and the exterior of the star. The interior
scalar field is obtained by directly integrating from the center of the star
where the value of $\Phi$ can be assigned arbitrarily as it only scales the
solution and the value of $\Phi^{\prime}$ is zero to keep
$\Phi^{\prime\prime}$ finite. The exterior solution is computed by integrating
from the surface of the star where the value of $\Phi$ can again be assigned
arbitrarily but the value of $\Phi^{\prime}$ needs to be deliberately selected
to make sure that the asymptotically flat condition $\Phi\rightarrow 0$ at
infinity is satisfied. Specifically speaking, a general exterior solution to
the linearized scalar field equation has the asymptotical expression

$\displaystyle\Phi\rightarrow\Phi_{+}e^{r/\lambdabar_{\Phi}}+\Phi_{-}e^{-r/\lambdabar_{\Phi}},$
(18)

where the functions $\Phi_{+}$ and $\Phi_{-}$ are finite when $r$ goes to
infinity. The physical solution that $\Phi$ goes to zero at infinity excludes
the $\Phi_{+}$ branch.

For calculating the interior solution, practical numerical integrations cannot
begin exactly at the center of the star due to the $\propto 1/r$ factor in the
scalar field equation. They begin at a radius very close to the center, and
hence corrections to the values of $\Phi$ and $\Phi^{\prime}$ at the starting
point are necessary. The corrections can be approximated by the Taylor
expansion of $\Phi$ at the center, whose expansion coefficient at the order of
$r^{2}$ is displayed in Eq. (40) in Appendix A. For calculating the exterior
solution, the shooting strategy combined with the bisection method is used to
find suitable values of $\Phi^{\prime}$ at the surface of the star to make
$\Phi$ vanish at infinity. The explicit expressions for the asymptotic
functions $\Phi_{+}$ and $\Phi_{-}$ are irrelevant to our numerical approach,
but for the completeness of the analysis as well as for illustrating the
asymptotic behavior of the nonlinear solution, they are derived in Appendix B.

Once the interior and the exterior solutions to the linearized scalar field
equation are computed separately, the condition for them to form an entire
solution is that the values of $\Phi^{\prime}/\Phi$ at the surface of the star
(denoted as $\Phi^{\prime}/\Phi|_{s}$) given by the solutions on both sides
match each other. For given metric functions $\mu$ and $\nu$ of a NS in GR,
the match is not guaranteed. In other words, only specific NSs in GR possess
solutions to the linearized scalar field equation. The parameters of these
NSs, for example, their central energy densities or their masses, are critical
in the sense that they mark possible onsets of spontaneous scalarization. To
numerically find the critical values of the central energy density, we compute
$\Phi^{\prime}/\Phi|_{s}$ from the interior and the exterior solutions for
different NSs in GR, and plot it versus the central energy densities of the
NSs. When the curve of $\Phi^{\prime}/\Phi|_{s}$ from the interior solution
intersects the curve of $\Phi^{\prime}/\Phi|_{s}$ from the exterior solution,
we have a match and obtain a critical central energy density.

For the convenience in numerical calculation, we use $l_{0}=10\,{\rm km}$ as
the length unit and parametrize other quantities using it in our work.
Especially, the parametrized versions of $m_{\Phi}$ and $\xi$ are denoted as
$a$ and $b$. In conventional units, $m_{\Phi}$ and $\xi$ in terms of $a$ and
$b$ are

$\displaystyle m_{\Phi}$ $\displaystyle=$ $\displaystyle\frac{ah}{2\pi
l_{0}}=a\times 1.97\times 10^{-11}\,{\rm eV},$ $\displaystyle\xi$
$\displaystyle=$ $\displaystyle b\,l_{0}^{2}=b\times 10^{8}\,{\rm m^{2}}.$
(19)

In Fig. 1, $\Phi^{\prime}/\Phi|_{s}$ from the interior and the exterior
solutions are plotted for changing values of the central energy density of the
NS with $a=1$ and $b$ taking positive values $b=0.9,\,3,\,10$ and negative
values $b=-0.04,\,-0.2,\,-1$. The metric of the NS in GR is solved using the
EOS AP4. Four other EOSs, including WFF1, SLy4, MPA1, and PAL1, are verified
to give similar results. The main observation from the plots is that the curve
of $\Phi^{\prime}/\Phi|_{s}$ from the interior solution and the curve of
$\Phi^{\prime}/\Phi|_{s}$ from the exterior solution start to intersect at the
highest central energy density when $|b|$ is sufficiently large. As $|b|$
increases, the existing intersection moves to left and a new intersection
comes out again at the highest central energy density. The procedure repeats
as long as $|b|$ continuously increases. So the larger $|b|$ is, the more
intersections accumulate.

For given values of $a$, the values of $|b|$ when the first intersection just
appears define the boundary of the theory parameter space for spontaneous
scalarization. The solid curves in Fig. 2 depict the numerical results of the
boundary when different EOSs are used, namely that the value of $b$ needs to
be either greater than $b_{+}$ or less than $b_{-}$ for spontaneous
scalarization to happen when $a$ is given. As the first intersection appears
at the highest central energy density and then moves to left when $|b|$
increases, we can ask what the value of $|b|$ is when the intersection goes
down to the central energy density that corresponds to the maximum mass of the
NS for the EOS used. This value of $|b|$ gives a naive estimate for the
boundary of the theory parameter space if the spontaneous scalarization is to
be stable. The dashed curves in Fig. 2 show the results for the five different
EOSs. For a given value of $a$, if the value of $b$ is between the solid curve
and the dashed curve for an EOS in consideration, then there are spontaneously
scalarized NSs in the theory though they are expected to be unstable due to
their extremely high central energy densities.

###  II.3 Solutions of spontaneously scalarized NSs

Figure 5: The metric components and the scalar field of three representative
scalarized NSs. The dashed and dotted curves show the metric functions of GR
NSs with the same masses for comparison. Notice the difference between the
profiles of the scalar field for $b>0$ and $b<0$. Also note that the scalar
field of the scalarized NS with $2.2\,M_{\odot}$ in the case of $a,b=1,-1$ has
a node, as it belongs to the second scalarization band in the $M$-$R$ diagram.

For values of $m_{\Phi}$ and $\xi$ in the parameter space for spontaneous
scalarization, solutions of spontaneously scalarized NSs can be searched
around the critical central energy densities given by the linearized analysis
described above. Noticing that only $\mu,\,\Phi$ and $p$ are coupled in the
group of the ODEs in Eqs. (12), the search can be conducted through adjusting
the value of $\Phi$ at the center of the star (denoted as $\Phi_{0}$) while
the value of $\mu$ at the center of the star is set to zero and the value of
$p$ at the center is given by the EOS once the central energy density
$\epsilon_{0}$ is set to a value close to one of the critical values obtained
in the linearized analysis. Similarly to the linearized case, the numerical
integrations have to start at a tiny but nonzero radius. So corrections to the
starting values of $\mu,\,\Phi$ and $p$ that can be approximated by their
Taylor expansions at the center of the star are requisite. Appendix A
discusses the expansions and gives the equations from which the expansion
coefficients at the order of $r^{2}$ can be numerically solved.

Differently from the linearized case, now it is unnecessary to integrate the
interior and the exterior regions separately. The integration goes directly
from the center of the star to a large enough radius, where $\Phi$ likely
blows up due to the same $e^{r/\lambdabar_{\Phi}}$ factor unless the central
value $\Phi_{0}$ is properly chosen so that the $e^{r/\lambdabar_{\Phi}}$
branch of the solution vanishes. By applying the shooting strategy to find out
the appropriate values of $\Phi_{0}$ for valid values of $\epsilon_{0}$,
solutions of spontaneous scalarization with $\Phi\rightarrow 0$ at infinity
are attained.

Representative numerical results are plotted on the mass-radius relation
diagram in Fig. 3. The first thing to be noticed is that solutions of
spontaneous scalarization are found around every critical central energy
density given by the linearized analysis for the case of $\xi<0$ while for the
case of $\xi>0$, only the critical central energy density corresponding to the
first intersection in the $\Phi^{\prime}/\Phi|_{s}$ plot (see Fig. 1) is an
onset of spontaneous scalarization. Detailed numerical integrations show that
around the critical central energy densities corresponding to the second and
higher intersections in the $\Phi^{\prime}/\Phi|_{s}$ plot for the case of
$\xi>0$, solutions of spontaneous scalarization fail to develop because of an
unexpected divergence soon after the integrations start. The divergence can
only originate from the high nonlinearity in Eqs. (12).

Table 1: Numerical properties of the five representative scalarized neutron stars plotted in Fig. 5.  Parameters |  ISCO properties |  perturbation properties
---|---|---
$a,b$ | $\Phi_{0}$ |  $\epsilon_{0}/c^{2}$ ($10^{15}\,{\rm g/cm^{3}}$) | $M\,(M_{\odot})$ |  $R$ (km) | Baryonic mass$\,(M_{\odot})$ |  $r$ (km) | $\frac{d\phi}{dt}\,({\rm ms}^{-1})$ | $\lambda/M^{5}$ |  $I$ ($10^{45}\,{\rm g\,cm^{2}}$)
$0,0$ (GR) | 0 | 0.9833 | 1.400 | 11.42 | 1.571 | 12.40 | 9.867 | $2.62\times 10^{2}$ | 1.324
$0.1,-1$ | $1.752\times 10^{-2}$ | 0.9860 | 1.400 | 11.45 | 1.571 | 12.35 | 9.929 | $2.77\times 10^{2}$ | 1.323
$0.1,3$ | $5.237\times 10^{-3}$ | 1.083 | 1.400 | 10.89 | 1.569 | 12.41 | 9.571 | $1.86\times 10^{3}$ | 1.241
$0,0$ (GR) | 0 | 1.532 | 2.000 | 11.00 | 2.392 | 17.72 | 6.907 | $1.64\times 10^{1}$ | 2.152
$1,3$ | $2.741\times 10^{-4}$ | 1.682 | 2.000 | 10.63 | 2.390 | 16.90 | 7.959 | $7.35\times 10^{1}$ | 2.067
$0,0$ (GR) | 0 | 2.308 | 2.200 | 10.31 | 2.706 | 19.49 | 6.279 | 3.62 | 2.309
$1,-0.2$ | $8.116\times 10^{-2}$ | 2.570 | 2.200 | 10.15 | 2.706 | 19.50 | 6.822 | 4.42 | 2.260
$1,-1$ | $1.388\times 10^{-2}$ | 2.319 | 2.200 | 10.32 | 2.706 | 19.69 | 6.717 | 5.75 | 2.315

Secondly, one cannot help noticing the very short ranges on the $M$-$R$
diagram for the solutions of spontaneous scalarization. This is especially
true for the case of $\xi<0$ where for each interval of spontaneous
scalarization the mass of the NS only changes less than $0.1\,M_{\odot}$,
making the solutions actually difficult to find without the hint of the
critical central energy densities from the linearized analysis. Our numerical
results also show that increasing the value of $|\xi|$ does not widen the mass
intervals of spontaneous scalarization. The fact that the occurrence of
spontaneous scalarization is limited to NSs in the narrow mass bands is a
consequence of the quadratic coupling function chosen. Using a different type
of coupling function, for example the exponential one studied by Doneva and
Yazadjiev (2018b), the intervals of spontaneous scalarization can be much
wider.

Thirdly, we notice that all the solutions of spontaneous scalarization
distribute themselves on the left side of the $M$-$R$ curves of GR for the
case of $\xi>0$, implying that these scalarized NSs are smaller than the NSs
in GR with the same mass. The situation for the case of $\xi<0$ is a little
subtle as there can be several scalarization bands along the $M$-$R$ curve of
GR and the scalarization bands are extremely narrow. A zooming-in observation
on them shows that the $M$-$R$ segments of the scalarized NSs can actually
depart the $M$-$R$ curves of GR on their right sides and then turn left and
cross the $M$-$R$ curves of GR when $\xi$ is sufficiently negative. In this
case, a scalarized NS can have a radius smaller than, equal to, or greater
than that of a NS in GR with the same mass, depending on where it sits along
the $M$-$R$ segment. Though it is theoretically interesting to capture the
details, it is unlikely useful in testing the theories with negative $\xi$ in
observations because the deviations from GR are anyhow too small.

Fourthly, a noticeable feature, at least for the case of $\xi>0$, is that the
$M$-$R$ segments of the scalarized NSs have only one end attached to the
$M$-$R$ curves of GR. This suggests that the $M$-$R$ curves bifurcate at the
onsets of spontaneous scalarization, and that the branches of spontaneous
scalarization develop as the scalar field at the center of the star increases
to finite values where the solutions of spontaneous scalarization cease to
exist. Plots of the scalar field at the center of the star versus the mass of
the star, as shown in Fig. 4, help to illustrate the existence and the
evolution of spontaneous scalarization. Starting at the onsets of spontaneous
scalarization where the central scalar $\Phi_{0}$ is technically zero, the
solutions evolve to reach smaller NS masses as $\Phi_{0}$ increases to finite
values and stops. Naturally the value of $\Phi_{0}$ should drop back to zero
after it reaches a maximum. The reason for the abrupt cutoff is that the
Taylor expansion coefficients $\mu_{2},\,\nu_{2}$ and $\Phi_{2}$ solved from
Eqs. (36) start to have imaginary parts when $\Phi_{0}$ and $\epsilon_{0}$
happen to be in certain regions. Noticing that it is because of the quartic
nature of the recurrence equations that they admit complex solutions for
$\mu_{2},\,\nu_{2}$ and $\Phi_{2}$, we can draw the conclusion that the high
nonlinearity in Eqs. (12) again restricts the existence of spontaneous
scalarization. A similar restriction for spontaneously scalarized black holes
in scalar-Gauss-Bonnet theories, expressed in an analytical inequality, can be
found in Refs. Antoniou  _et al._ (2018); Doneva and Yazadjiev (2018a); Silva
_et al._ (2018).

Lastly, a careful comparison of the 8 panels in Fig. 3 tells us that as
$m_{\Phi}$ increases or $|\xi|$ decreases, the intervals of spontaneous
scalarization shift to the higher end along the $M$-$R$ curves of GR, and vice
versa.

From Fig. 3 we see that the two main parameters of the scalarized NSs, namely
their masses and radii, barely deviate from the values of NSs in GR for the
case of $\xi<0$. This makes distinguishing those scalarized NSs from NSs in GR
using observations very difficult. In fact, even for the case of $\xi>0$, the
largest deviation in radius between a scalarized NS and a NS in GR with the
same mass is only about $1\,{\rm km}$ as shown in Fig. 3. This just lies
within the uncertainties of current measurements of NS radius Abbott  _et al._
(2018a); Riley  _et al._ (2019); Miller  _et al._ (2019); Riley  _et al._
(2021); Miller  _et al._ (2021). The degeneracy is even obvious when
individual solutions are compared. In Fig. 5 we plot the metric components of
three representative scalarized NSs calculated using the EOS AP4 to compare
with the metric components of the NSs in GR with the same EOS and masses. To
conclude our study of the spherical solution of spontaneously scalarized NSs,
certain parameters of five representative scalarized NSs are displayed in
Table 1 for reference.

##  III Time-independent perturbations

In this section, we deal with a slowly rotating and tidally deformed NS under
spontaneous scalarization using perturbation theory. The formalism that
decomposes the metric perturbation into spherical harmonics introduced by
Regge and Wheeler (1957) and later developed for NSs by Thorne and
Campolattaro (1967) is adopted. We first give a brief review of the formalism.

###  III.1 Review of the perturbation theory in spherical background

Spherical harmonics form complete sets of orthogonal functions on a 2-sphere,
and therefore are naturally the bases to expand the angular dependence of
perturbations around a spherical background. Depending on the rotational
properties of the quantities on the 2-sphere, scalar, vector, and tensor
spherical harmonics may be used. Following the definitions of Thorne (1980),
the transverse vector spherical harmonics and the traceless transverse tensor
spherical harmonics with rank-2 used here are constructed from the usual
scalar spherical harmonics $Y_{lm}$ via,

$\displaystyle\bm{Y}^{E}_{lm}$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{l(l+1)}}r\bm{\nabla}Y_{lm}=\frac{r}{\sqrt{l(l+1)}}\left({\bm{d\theta}}\,\partial_{\theta}+{\bm{d\phi}}\,\partial_{\phi}\right)Y_{lm},$
$\displaystyle\bm{Y}^{B}_{lm}$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{l(l+1)}}\bm{r}\times\bm{\nabla}Y_{lm}=\frac{r}{\sqrt{l(l+1)}}\left(-\frac{{\bm{d\theta}}}{\sin\theta}\,\partial_{\phi}+{\bm{d\phi}}\,\sin\theta\,\partial_{\theta}\right)Y_{lm},$
$\displaystyle\bm{T}^{E2}_{lm}$ $\displaystyle:=$
$\displaystyle\sqrt{\frac{2}{(l-1)(l+2)}}\left[r\bm{\nabla}\bm{Y}^{E}_{lm}\right]^{\rm
ST}$ $\displaystyle=$ $\displaystyle
N_{l}\,r^{2}\left[{\bm{d\theta}}\otimes{\bm{d\theta}}\,\partial_{\theta}^{2}+{\bm{d\phi}}\otimes{\bm{d\phi}}\left(\sin\theta\cos\theta\,\partial_{\theta}+\partial_{\phi}^{2}\right)+\left({\bm{d\theta}}\otimes{\bm{d\phi}}+{\bm{d\phi}}\otimes{\bm{d\theta}}\right)\partial_{\phi}\left(\partial_{\theta}-\cot\theta\right)\right]^{\rm
T}Y_{lm},$ $\displaystyle\bm{T}^{B2}_{lm}$ $\displaystyle:=$
$\displaystyle\sqrt{\frac{2}{(l-1)(l+2)}}\left[r\bm{\nabla}\bm{Y}^{B}_{lm}\right]^{\rm
ST}$ (20) $\displaystyle=$ $\displaystyle
N_{l}\,r^{2}\left[\left(-\frac{{\bm{d\theta}}\otimes{\bm{d\theta}}}{\sin\theta}+{\bm{d\phi}}\otimes{\bm{d\phi}}\,\sin\theta\right)\partial_{\phi}\left(\partial_{\theta}-\cot\theta\right)+{\bm{d\theta}}\otimes{\bm{d\phi}}\,\sin\theta\,\partial_{\theta}^{2}-{\bm{d\phi}}\otimes{\bm{d\theta}}\left(\cos\theta\,\partial_{\theta}+\frac{\partial_{\phi}^{2}}{\sin\theta}\right)\right]^{\rm
S}Y_{lm},$

where $\bm{\nabla}$ is the gradient in three-dimensional Euclidean space,
$N_{l}:=\sqrt{2/\left[(l-1)l(l+1)(l+2)\right]}$, and
$\\{{\bm{dr}},\,{\bm{d\theta}},\,{\bm{d\phi}}\\}$ is used as the basis for
covariant components of a vector in the spherical coordinates
$(r,\theta,\phi)$. The superscript “S” means symmetrizing the tensor and the
superscript “T” means removing the trace. By dropping the factor
$N_{l}\,r^{2}$, the covariant components of $\bm{T}^{E2}_{lm}$ without
removing the trace and the covariant components of $\bm{T}^{B2}_{lm}$ can be
recognized as $\psi_{l{\phantom{m}}{\mu\nu}}^{{\phantom{l}}m}$ and
$\chi_{l{\phantom{m}}{\mu\nu}}^{{\phantom{l}}m}$ used by Regge and Wheeler
(1957).

The perturbative quantities to be expanded using the spherical harmonics are
the metric perturbation $\delta g_{{\mu\nu}}$, the fluid displacement
$\xi_{\mu}$ and the scalar perturbation $\delta\Phi$. On the 2-sphere, the
metric perturbation can be decomposed as four scalars:

$\displaystyle\delta g_{tt},\ \delta g_{tr},\ \delta g_{rr},\ \delta
g_{\theta\theta}+\frac{1}{\sin^{2}\theta}\delta g_{\phi\phi},$

two vectors:

$\displaystyle\delta g_{t\theta}\,{\bm{d\theta}}+\delta
g_{r\phi}\,{\bm{d\phi}},\ \delta g_{r\theta}\,{\bm{d\theta}}+\delta
g_{r\phi}\,{\bm{d\phi}},$

and a traceless tensor:

$\displaystyle\Big{[}\delta
g_{\theta\theta}\,{\bm{d\theta}}\otimes{\bm{d\theta}}+\delta
g_{\theta\phi}\left({\bm{d\theta}}\otimes{\bm{d\phi}}+{\bm{d\phi}}\otimes{\bm{d\theta}}\right)+\delta
g_{\phi\phi}\,{\bm{d\phi}}\otimes{\bm{d\phi}}\Big{]}^{\rm T},$

while the fluid displacement can be decomposed as two scalars:

$\displaystyle\xi_{t},\ \xi_{r},$

and a vector:

$\displaystyle\xi_{\theta}\,{\bm{d\theta}}+\xi_{\phi}{\bm{d\phi}}.$

Furthermore, the expansions of all the 11 quantities can be conveniently
divided into two parts according to the parity of the spherical harmonics: the
so-called even-parity ones, $\\{Y_{lm},\bm{Y}^{E}_{lm},\bm{T}^{E2}_{lm}\\}$,
that pick up a sign of $(-1)^{l}$ under a parity transformation, and the so-
called odd-parity ones, $\\{\bm{Y}^{B}_{lm},\bm{T}^{B2}_{lm}\\}$, that pick up
a sign of $(-1)^{l+1}$ under the same transformation. Because the perturbative
field equations directly inherit the parity properties of $\delta g_{\mu\nu}$
and $\xi_{\mu}$, perturbations with different types of parity do not mix. The
same conclusion is also true for the angular momentum numbers $l$ and $m$,
allowing us to consider perturbations with specific values of $l$ and $m$
rather than a sum over all the values of them.

With the separation of even and odd parities, the 11 quantities become 22. But
since there is no scalar spherical harmonics of odd parity, 7 of them
trivially vanish. Moreover, the freedom to choose coordinates at the order of
the perturbation can help us eliminate another 4 of them. Using the Regge-
Wheeler gauge Regge and Wheeler (1957), the four quantities eliminated are all
in the metric perturbation: the two vectors with even parity and the traceless
tensor with even parity and with odd parity. This finally leaves us $11$
nonvanishing quantities to express in terms of the spherical harmonics of a
specific parity as well as specific $l$ and $m$, namely, for even parity

$\displaystyle\delta g_{tt}=e^{2\nu}H_{0}Y_{lm},\quad\delta g_{tr}$
$\displaystyle=$ $\displaystyle H_{1}Y_{lm},\quad\delta
g_{rr}=e^{2\mu}H_{2}Y_{lm}$ $\displaystyle\delta
g_{\theta\theta}+\frac{1}{\sin^{2}\theta}\delta g_{\phi\phi}$ $\displaystyle=$
$\displaystyle 2r^{2}KY_{lm},$ $\displaystyle\xi_{t}=W_{0}Y_{lm},\quad\xi_{r}$
$\displaystyle=$ $\displaystyle W_{1}Y_{lm}$
$\displaystyle\xi_{\theta}\,{\bm{d\theta}}+\xi_{\phi}\,{\bm{d\phi}}$
$\displaystyle=$ $\displaystyle V\bm{Y}^{E}_{lm}\frac{\sqrt{l(l+1)}}{r},$
$\displaystyle\delta\Phi$ $\displaystyle=$ $\displaystyle XY_{lm},$ (21)

and for odd parity

$\displaystyle\delta g_{t\theta}\,{\bm{d\theta}}+\delta
g_{t\phi}\,{\bm{d\phi}}$ $\displaystyle=$ $\displaystyle
h_{0}\bm{Y}^{B}_{lm}\frac{\sqrt{l(l+1)}}{r},$ $\displaystyle\delta
g_{r\theta}\,{\bm{d\theta}}+\delta g_{r\phi}\,{\bm{d\theta}}$ $\displaystyle=$
$\displaystyle h_{1}\bm{Y}^{B}_{lm}\frac{\sqrt{l(l+1)}}{r},$
$\displaystyle\xi_{\theta}\,{\bm{d\theta}}+\xi_{\phi}\,{\bm{d\phi}}$
$\displaystyle=$ $\displaystyle h_{2}\bm{Y}^{B}_{lm}\frac{\sqrt{l(l+1)}}{r},$
(22)

where $H_{0},\,H_{1},\,H_{2},\,K,\,W_{0},\,W_{1},\,V,\,X$ and
$h_{0},\,h_{1},\,h_{2}$ are unknown functions of $t$ and $r$ to be solved from
the perturbative field equations.

To obtain differential equations for the 11 perturbative functions using the
field equations (6), the fluid displacement $\xi_{\mu}$ needs to be related to
the perturbation of the 4-velocity, the perturbation of the energy density,
and the perturbation of the pressure. The relations are given, for example, by
Gittins  _et al._ (2020), as

$\displaystyle\Delta u^{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}u^{\mu}u^{\alpha}u^{\beta}\Delta g_{\alpha\beta},$
$\displaystyle\frac{\Delta n}{n}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(g^{\alpha\beta}+u^{\alpha}u^{\beta}\right)\Delta
g_{\alpha\beta},$ (23)

for the Lagrangian perturbation of the 4-velocity $\Delta u^{\mu}$ and the
Lagrangian perturbation of the particle number density $\Delta n$ in terms of
the Lagrangian perturbation of the metric $\Delta g_{\alpha\beta}$. The fluid
displacement comes into Eqs. (23) when we change the Lagrangian perturbations
to the Eulerian perturbations via the definition
$\Delta:=\delta+{\mathcal{L}}_{\xi}$, where ${\mathcal{L}}_{\xi}$ is the Lie
derivative along the fluid displacement. Therefore, the Eulerian perturbation
of the 4-velocity to be used in the perturbative field equations is

$\displaystyle\delta u^{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}u^{\mu}u^{\alpha}u^{\beta}\left(\delta
g_{\alpha\beta}+D_{\alpha}\xi_{\beta}+D_{\beta}\xi_{\alpha}\right)$ (24)
$\displaystyle-\left(\xi^{\lambda}D_{\lambda}u^{\mu}-u^{\lambda}D_{\lambda}\xi^{\mu}\right).$

The Eulerian perturbations of the energy density and the pressure can be
calculated from their Lagrangian perturbations

$\displaystyle\Delta\epsilon=\left(\epsilon+p\right)\frac{\Delta
n}{n},\quad\Delta p=\gamma p\frac{\Delta n}{n},$ (25)

where the first equation is a result of the first law of thermodynamics under
the assumption that the fluid elements are adiabatic, and the second equation
can be regarded as the definition of the adiabatic index $\gamma$. Given an
EOS relating $\epsilon$ and $p$, the adiabatic index $\gamma$ can be obtained
via

$\displaystyle\gamma=\frac{\epsilon+p}{p}\frac{dp}{d\epsilon}.$ (26)

Now substituting $\delta g_{\alpha\beta}$ and $\xi_{\mu}$ from Eqs. (21) and
(22), as well as the zeroth-order quantities, into Eqs. (23), (24) and (25),
we can obtain the spherical decompositions for the perturbations of the
4-velocity, the energy density, and the pressure. The perturbative field
equations then can be obtained. Once the angular dependence is detached, they
give equations for $H_{0},\,H_{1},\,H_{2},\,K,\,W_{0},\,W_{1},\,V,\,X$ and
$h_{0},\,h_{1},\,h_{2}$. In Appendices C and E, these equations are displayed
with the time dependence set to $e^{i\sigma t}$, where the constant $\sigma$
stands for the frequency of a normal or quasinormal mode of the star.

Figure 6: The dimensionless tidal deformability $\lambda/M^{5}$ versus the
mass of NSs. The dashed curves are GR results, while the solid curves are
results of the scalarized NSs in the lowest mass bands in Fig. 3. Please
notice the discontinuities in the panels of $a=1,b=3,10$. Details are
discussed in the text.

###  III.2 Tidal deformation: static even-parity $l=2$

Our first application of the perturbation theory in spherical background is
calculating the perturbation of the metric under a static external tidal
field. The equations and the numerical approach to solve them are summarized
in Appendix C. The quantity to be calculated is the tidal deformability
$\lambda$, defined via Hinderer (2008)

$\displaystyle Q_{2m}=-\lambda{\mathcal{E}}_{2m},$ (27)

where $Q_{2m}$ with $m=0,\pm 1,\pm 2$ are the induced quadrupole moments of an
originally spherical star and ${\mathcal{E}}_{2m}$ are the spherical
components of the static external tidal field. Both $Q_{2m}$ and
${\mathcal{E}}_{2m}$ are defined in the asymptotic expansion of the metric
perturbation $\delta g_{tt}$ of even parity and angular momentum number $l=2$,
namely

$\displaystyle\delta g_{tt}$ $\displaystyle=$ $\displaystyle
e^{2\nu}H_{0}\sum\limits_{m=-2}^{2}A_{m}Y_{2m}$ (28)
$\displaystyle\rightarrow$
$\displaystyle\frac{3}{r^{3}}\sum\limits_{m=-2}^{2}Q_{2m}Y_{2m}-r^{2}\sum\limits_{m=-2}^{2}{\mathcal{E}}_{2m}Y_{2m},$

where $A_{m}$ are $m$-related constants to give angular dependence that
matches the second line. Because the perturbative function $H_{0}$ is a
solution to Eqs. (54) that are independent of $m$, its leading asymptotic
terms

$\displaystyle
H_{0}\rightarrow\frac{Q_{2m}}{A_{m}}\frac{3}{r^{3}}-\frac{{\mathcal{E}}_{2m}}{A_{m}}r^{2},$
(29)

read off from Eq. (28), must be $m$-independent. Therefore, the tidal
deformability $\lambda=-Q_{2m}/{\mathcal{E}}_{2m}$ does not depend on $m$.

To calculate $H_{0}$ and extract its leading asymptotic expansion coefficients
as shown in Eq. (29), we set $\sigma=0$ in Eqs. (54) so they eventually
simplify to three coupled ODEs for $H_{0},\,K$ and $X$ that can be numerically
integrated from the center of the star. As usual, practical integrations start
at a tiny radius rather than exactly at the center of the star to avoid
numerical zeros in denominators. The values of $H_{0},\,K,\,X$ and
$X^{\prime}$ at the tiny radius where the integrations start are then
necessary inputs in the integrator. They are given by the Taylor expansions in
Eqs. (55) with two of the expansion coefficients, say $H_{00}$ and $X_{0}$ as
used in Eqs. (56), free to choose for the moment. The two coefficients are
actually not free as soon as the physical asymptotic behavior is required for
the solution. Similarly to the spherical solution, the asymptotic expression
of a general solution to the perturbative equations contains both the factors
of $e^{-r/\lambdabar_{\Phi}}$ and the factors of $e^{r/\lambdabar_{\Phi}}$,
while the contributions from the latter is supposed to be eliminated in the
physical solution by adjusting the inputs of the integration through $H_{00}$
and $X_{0}$. The problem of generating the physical solution is much
simplified since the perturbative equations are homogeneous in the
perturbative functions. Following Pani and Berti (2014), the physical solution
can be found by linearly combining any two independent solutions to make $X$
vanish at infinity instead of employing the shooting strategy that aims at
finding suitable values of $H_{00}$ and $X_{0}$.

In our numerical approach, two independent solutions are found with
$(H_{00},\,X_{0})=(1,\,0)$ and $(H_{00},\,X_{0})=(1,1)$, and then after the
physical solution is constructed from them, two methods for computing the
tidal deformability are available. The first one is extracting the
coefficients $Q_{2m}/A_{m}$ and ${\mathcal{E}}_{2m}/A_{m}$ shown in Eq. (29)
by directly fitting the numerical solution of $H_{0}$ in the asymptotic region
and then computing $\lambda$. The second one is using the formula developed by
Hinderer (2008) with all the quantities defined on the surface of the star
therein substituted by the corresponding values at a sufficiently large $r$
where the terms containing the factors of $e^{-r/\lambdabar_{\Phi}}$ can be
safely dropped so the expansion of $H_{0}$ is exactly the same as that in GR
for the formula to apply Hu  _et al._ (2021). The two methods are no doubt
theoretically equivalent but the second one turns out to be too sensitive to
numerical errors when $\lambdabar_{\Phi}$ is large. The first method is
adopted in this work. The issue of the second method with large
$\lambdabar_{\Phi}$ is discussed in detail in Appendix D.

Figure 6 shows the results of the tidal deformability for the NSs in Fig. 3.
For $\xi<0$, the deviations are anyhow tiny so we only plotted the results for
the scalarized NSs in the lowest mass bands. The results for $\xi>0$ however
are quite appealing. First and foremost, the scalarized NSs can have the
dimensionless tidal deformability $\lambda/M^{5}$ an order of magnitude larger
than that of the NSs in GR. As a matter of fact, it is counterintuitive that
the tidal deformability of the scalarized NSs is larger than that of the NSs
in GR at all since according to the $M$-$R$ diagram in Fig. 3, these
scalarized NSs are more compact than their GR counterparts. Secondly, the
plots for $a=1$ show jumps where the sign of the tidal deformability changes,
suggesting that the tidal deformability of the scalarized NSs approaches
infinities near the jumps. By investigating the numerically fitted expansion
coefficients of $H_{00}$, we find that the coefficient in front of $r^{2}$
decreases to zero for the special scalarized NSs that have masses at the
discontinuities in the plots for $a=1,\,b=3,10$ in Fig. 6. For scalarized NSs
near the special masses, the sign of the $r^{2}$-coefficient in the asymptotic
expansion of $H_{00}$ changes. This feature neither shows up in GR nor in the
nonminimal scalar-tensor theories studied before Hu  _et al._ (2021). It is a
qualitative difference between the scalar-Gauss-Bonnet gravity studied here
and GR.

The large deviations in the tidal deformability of the scalarized NSs from
their GR counterparts are no doubt most helpful in testing the scalar-Gauss-
Bonnet gravity in astrophysical observations. Simply taking the result from
the observation of GW170817: $\lambda/M^{5}\lesssim 720$ for $M\approx
1.4M_{\odot}$ Abbott  _et al._ (2019b), we can already exclude the point
$(a,\,b)=(0.1,\,3)$ in the parameter space for the two popular EOSs AP4 and
SLy4. To rule out a region in the parameter space of $(a,\,b)$, numerical
results of the tidal deformability with $a$ and $b$ taking values near the
point $(0.1,3)$ are necessary to compare with the observational constraint.
The implementation of the comparison algorithm into numerical code is beyond
the scope of the present work but certainly worth further researching as more
and more constraints on the tidal deformability of NSs are expected from
future GW observations Abbott  _et al._ (2018b).

###  III.3 Slow rotation: static odd-parity $l=1$

Figure 7: The moment of inertia versus the mass of NSs. The dashed curves are
GR results, while the solid segments are results of the scalarized NSs in the
lowest mass bands in Fig. 3.

The second application of the perturbation theory in spherical background is
calculating the moment of inertia of slowly rotating scalarized NSs. The
problem boils down to solving the first equation in Eqs. (67) with $h_{1}=0$
and $\sigma\rightarrow 0$ while the product $\sigma h_{2}$ keeps finite and
represents the radial part of the angular velocity of the star via the match
(73). Following Damour and Esposito-Farèse (1996), we assume the angular
velocity $\Omega$ to be constant so that the moment of inertia of the star can
be defined as

$\displaystyle I:=\frac{J}{\Omega},$ (30)

where the angular momentum $J$ is extracted from the leading asymptotic term
of $g_{t\phi}$. The match (73) indicates $l=1$ for a constant $\Omega$, and
$i\sigma h_{2}$ is to be substituted by $-\sqrt{3/(4\pi)}\,\Omega r^{2}$ in
the first equation in Eqs. (67), which becomes a second-order ODE with a
source term for $h_{0}$.

The numerical method to solve $h_{0}$ is integrating from the center of the
star where $h_{0}$ has the expansion

$\displaystyle h_{0}=h_{00}r^{2}+h_{02}r^{4}+\cdots\,,$ (31)

when $l=1$. The leading coefficient $h_{00}$ can be set to 1 and the next-to-
leading coefficient $h_{02}$ can be calculated using Eq. (76) for given values
of $\Omega$. For an arbitrary value of $\Omega$, the solution of $h_{0}$
includes an asymptotic term proportional to $r^{2}$ as well as an asymptotic
term proportional to $1/r$. The former needs to be eliminated by elaborately
choosing a proper value of $\Omega$ so that the metric component $g_{t\phi}$
has the desired asymptotic behavior

$\displaystyle g_{t\phi}=h_{0}\sin\theta\,\partial_{\theta}Y_{10}$
$\displaystyle=$ $\displaystyle-\sqrt{\frac{3}{4\pi}}\sin^{2}\theta\,h_{0}$
(32) $\displaystyle\rightarrow$ $\displaystyle-\frac{2J}{r}\sin^{2}\theta,$

from which the angular momentum $J$ can be read off.

Before showing the numerical results, a comment on the equivalence between the
approach used here and the approach used by Damour and Esposito-Farèse (1996)
might be worthwhile. The metric ansatz
$g_{t\phi}=\left(\omega-\Omega\right)r^{2}\sin^{2}\theta$ can be compared with
Eq. (32) to find $h_{0}=-\sqrt{4\pi/3}\left(\omega-\Omega\right)r^{2}$.
Therefore solving $h_{0}$ is equivalent to solving $\omega$ only that the
constant angular velocity of the star $\Omega$ appears as an input when
solving $h_{0}$, while it becomes an output together with the angular momentum
$J$ in the asymptotic expansion of $\omega$,

$\displaystyle\omega\rightarrow\Omega-\frac{2J}{r^{3}}+O\left(r^{-4}\right).$
(33)

We have to admit that though we derived the equation for $h_{0}$, it is
converted to the equation for $\omega$ to solve numerically. Because there is
no adjustable input in solving $\omega$, it is practically more convenient
than solving $h_{0}$.

Figure 8: Plots for $I/M^{3}$ versus $\lambda/M^{5}$. The dashed curve is the
numerically fitted universal relation in GR Yagi and Yunes (2017), while the
solid curves are results of the scalarized NSs in the lowest mass bands in
Fig. 3. The panels of $a=1,\,b=3,10$ only plot values of $\lambda$ on the
right side of the discontinuities shown in Fig. 6.

In Fig. 7, the moment of inertia of the scalarized NSs in the lowest mass
bands in Fig. 3 is compared with that of NSs in GR. We see that the scalarized
ones all have slightly smaller values of the moment of inertia compared with
NSs in GR when the masses are the same, intuitively agreeing with the fact
that the radii of most of them are smaller than those of their GR
counterparts.222However, it also raises a puzzle for the remaining of them
that have larger radii than their GR counterparts.

With the tidal deformability and the moment of inertia calculated, we are
ready to check if the universal relation between the dimensionless quantities
$\bar{\lambda}:=\lambda/M^{5}$ and $\bar{I}:=I/M^{3}$ discovered by Yagi and
Yunes (2017) in GR applies to the scalar-Gauss-Bonnet theory studied here.
Since Fig. 6 shows noticeable difference between the tidal deformability of
the scalarized NSs and the tidal deformability of the NSs in GR for $\xi>0$
while the difference in the moment of inertia shown in Fig. 7 is very
restrained, we expect that at least for $\xi>0$, the scalarized NSs have
different $\bar{I}$-$\bar{\lambda}$ relations rather than the same universal
relation as in GR. The plots in Fig. 8 verify the expectation. The
$\bar{I}$-$\bar{\lambda}$ segments corresponding to the scalarized NSs for
$\xi>0$ visibly deviate from the universal relation in GR. Moreover, the
$\bar{I}$-$\bar{\lambda}$ segments for different EOSs do not overlap, so they
fail to form a universal relation that is independent of EOSs. This could have
implication for observations. Universal relations involving quantities derived
from higher-order calcualtion in the spin, e.g. the quadrupole moment and the
shape parameter (cf. Ref. Gao  _et al._ (2022)), are beyond the scope of this
study.

##  IV Summary

In this work which complements the study of spontaneous scalarization of NSs
in the massless scalar-Gauss-Bonnet theory of gravitation by Doneva and
Yazadjiev (2018b), the same phenomenon is investigated in detail for the
massive scalar-Gauss-Bonnet theory. A linearized analysis on the scalar field
equation is conducted to discover the parameter space as well as the onsets of
spontaneous scalarization first. Then, numerical solutions of spontaneous
scalarization are computed and their properties are discussed. To gain more
insights into the scalarized NSs, their static perturbative properties are
investigated. These properties include the tidal deformability and the moment
of inertia, both of which we treat under the unified framework of spherical
decomposition.

To summarize the results, we find that our solutions share similar qualitative
characteristics with the solutions in the massless theory. First, solutions of
spontaneous scalarization exist for both positive and negative coupling
constants. Second, the scalar field at the center of the star emerges from
zero at the onsets of spontaneous scalarization, and develops to finite values
while the mass of the star changes (see Fig. 4). But it does not drop back to
zero again, namely that the onsets of spontaneous scalarization obtained by
the linearized analysis are one-sided. Solutions of spontaneous scalarization
cease to exist on the other end due to the nonlinearity of the leading-order
recurrence equations (36) at the center of the star. Equations (36) only have
real solutions for a certain range of the central parameters
$\epsilon_{0},p_{0}$ and $\Phi_{0}$, restricting the possible combinations of
the star mass $M$ and the central scalar field $\Phi_{0}$ for solutions of
spontaneous scalarization. Third, for negative coupling constants, spontaneous
scalarization can happen in more than one interval along the $M$-$R$ curve
(see Fig. 3). The scalarization interval corresponding to the lowest mass band
has scalar solutions possessing zero node. The scalar solutions in the second
lowest mass band possess one node, and so on.

More importantly, the quantitative details of our results show differences
from those in Ref. Doneva and Yazadjiev (2018b) since we use a different type
of coupling function and consider a massive scalar field. The coupling
function used here is quadratic in $\Phi$, which is intuitively weaker than
the exponential coupling function used by Doneva and Yazadjiev (2018b).
Therefore, the ranges of the star mass for spontaneous scalarization in our
theory are narrow compared to the results in Ref. Doneva and Yazadjiev
(2018b). We notice that increasing the absolute value of the coupling constant
$\xi$ in our theory does not lead to wider bands of the star mass for
spontaneous scalarization, indicating that the widths of the mass bands depend
on the type of coupling function rather than the magnitude of the coupling
constant. Because of the smallness of the spontaneous scalarization bands, we
have to zoom in on their $M$-$R$ diagrams to observe so that we find that the
$M$-$R$ segments of the scalarized NSs can cross the $M$-$R$ curves in GR when
the coupling constant $\xi$ is sufficiently negative, a fact that might be
easily ignored if the $M$-$R$ segments of the scalarized NSs spread out
longer.

The mass of the scalar field, according to Fig. 3, drives the spontaneously
scalarized NSs to the higher end along the $M$-$R$ curves in GR. Apart from
that, its effect is mainly the Yukawa suppression factor
$e^{-r/\lambdabar_{\Phi}}$ that once dropped, makes the asymptotic behavior of
the scalarized solutions identical to that of the solutions in GR. In case the
exponential terms are desired, for example, when the reduced Compton
wavelength of the scalar field $\lambdabar_{\Phi}$ is much larger than the
radius of the region that one is interested in, we employ the idea of double
expansion, namely expanding the asymptotic solution in factors of $M/r$ and
factors of $e^{-r/\lambdabar_{\Phi}}$, to come up with the asymptotic
expressions (46) for the first time to our knowledge. Unexpectedly, we find
that Eqs. (46) do not have limits when $\lambdabar_{\Phi}\rightarrow\infty$,
invalidating our attempt of matching them to the asymptotic expansions (47) in
the massless scalar-tensor theory.

The scalarized NSs found in our work show very moderate deviations from NSs in
GR on the $M$-$R$ diagram, making the theory very easy to escape constraints
from astrophysical observations of the mass and radius of NSs, especially with
the uncertainty from various EOSs. Therefore, we further calculated their
time-independent perturbative properties, including the tidal deformability
and the moment of inertia. While the moment of inertia of the scalarized NSs
only has insignificant differences from that of NSs in GR, remarkable
deviations in the tidal deformability of the scalarized NSs from their GR
counterparts are found for the case of a positive coupling constant. Not only
can the magnitude of the tidal deformability of these scalarized NSs be
significantly higher than that of the NSs in GR, but also the sign of the
tidal deformability changes to negative where a discontinuity pops up at a
special NS mass when the mass of the scalar field in the theory is large (see
the panels of $a=1,\,b=3,10$ in Fig. 6). The discontinuity is caused by the
vanishing of the $r^{2}$-coefficient in the asymptotic expansion of the metric
perturbation $\delta g_{tt}$, a feature that does not occur in GR.

Distinguishing GR and the scalar-Gauss-Bonnet theory with a positive coupling
constant becomes feasible with the discrepancies in the tidal deformability of
the NSs. For the scalar-Gauss-Bonnet theory with a negative coupling constant,
the quasinormal modes of the scalarized NSs might shed light on the way to
distinguish them from NSs in GR using observations, which will be studied in a
follow-up work. In addition, as the moment of inertia barely shows differences
between the scalarized NSs and the NSs in GR, the renowned universal relation
between the dimensionless tidal deformability $\lambda/M^{5}$ and the
dimensionless moment of inertia $I/M^{3}$ discovered by Yagi and Yunes (2017)
in GR no longer holds for the scalarized NSs (see Fig. 8). The tidal
deformability of NSs has begun to be accessible in the gravitational wave data
by the LIGO-Virgo-KAGRA Collaboration Abbott  _et al._ (2017, 2019b, 2018b),
and the moment of inertia of binary pulsars is expected to be extractable via
pulsar timing in the coming era of the Square Kilometre Array (SKA) Hu  _et
al._ (2020); Shao  _et al._ (2015); Weltman  _et al._ (2020). Once they are
measured for a sufficient group of NSs with necessary precision, we might be
able to tell whether GR should be regarded as an approximation of a
fundamentally motivated scalar-Gauss-Bonnet theory Zwiebach (1985); Boulware
and Deser (1985); Scherk and Schwarz (1974), or at least, to exclude an
appreciable region in the parameter space of it.

###### Acknowledgements.

We are grateful to the anonymous referee for the useful suggestions and
comments, and Zexin Hu for helpful discussions. This work was supported by the
National SKA Program of China (2020SKA0120300), the National Natural Science
Foundation of China (11975027, 11991053, 11721303), the Young Elite Scientists
Sponsorship Program by the China Association for Science and Technology
(2018QNRC001), the Max Planck Partner Group Program funded by the Max Planck
Society, and the High-performance Computing Platform of Peking University. RX
is supported by the Boya Postdoctoral Fellowship at Peking University.

##  Appendix A Expansion of the spherical solution at the center of the star

To investigate the behavior of the solutions to the set of ODEs in Eqs. (12)
at the center of the star, the Taylor expansions

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}\mu_{n}r^{n},\quad\nu=\sum\limits_{n=0}^{\infty}\nu_{n}r^{n},\quad\Phi=\sum\limits_{n=0}^{\infty}\Phi_{n}r^{n},$
$\displaystyle p$ $\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}p_{n}r^{n},\quad\epsilon=\sum\limits_{n=0}^{\infty}\epsilon_{n}r^{n},$
(34)

are substituted into Eqs. (12) to find the recurrence relations for the
expansion coefficients. We are only interested in the physical solutions that
are smooth at $r=0$, namely those with

$\displaystyle\mu_{1}=\nu_{1}=\Phi_{1}=p_{1}=\epsilon_{1}=0.$ (35)

Then, the leading-order recurrence relations fix $\mu_{0}$ to zero while
leaving $p_{0}$ and $\Phi_{0}$ free, from which nonvanishing
$\mu_{n},\,\nu_{n},\,\Phi_{n}$ and $p_{n}$ for even $n$’s are generated by the
higher-order recurrence relations. Note that the coefficients $\epsilon_{n}$
are related to $p_{n}$ via the EOS.

In particular, the coefficients in front of $r^{2}$ are given by

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\mu_{2}-8\frac{dF}{d\Phi}\mu_{2}\Phi_{2}-\frac{4\pi}{3}\epsilon_{0}-\frac{1}{12}U,$
$\displaystyle 0$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\mu_{2}+\nu_{2}-8\frac{dF}{d\Phi}\nu_{2}\Phi_{2}-2\pi
p_{0}+\frac{1}{8}U,$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle-8\frac{dF}{d\Phi}\mu_{2}\nu_{2}+\Phi_{2}-\frac{1}{12}\frac{dU}{d\Phi},$
$\displaystyle 0$ $\displaystyle=$
$\displaystyle(\epsilon_{0}+p_{0})\nu_{2}+p_{2},$ (36)

where $dF/d\Phi,\,U$ and $dU/d\Phi$ are understood as their values at
$\Phi=\Phi_{0}$. For a given pair of $p_{0}$ and $\Phi_{0}$,
$\mu_{2},\,\nu_{2},\,p_{2}$ and $\Phi_{2}$ can be solved. However, noticing
that the first three equations in (36) are nonlinear, real solutions might not
exist. In fact, eliminating two variables among
$\\{\mu_{2},\,\nu_{2},\,\Phi_{2}\\}$ leads to a quartic equation for the third
one. For example, by eliminating $\nu_{2}$ and $\Phi_{2}$, we find

$\displaystyle\mu_{2}^{4}+a_{3}\,\mu_{2}^{3}+a_{1}\,\mu_{2}+a_{0}=0,$ (37)

with

$\displaystyle a_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{4608}\left(16\pi\epsilon_{0}+U\right)^{2}\left(\frac{dF}{d\Phi}\right)^{-2},$
$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{1152}\left(16\pi\epsilon_{0}+U\right)\left(2\frac{dF}{d\Phi}\frac{dU}{d\Phi}-3\right)\left(\frac{dF}{d\Phi}\right)^{-2},$
$\displaystyle a_{3}$ $\displaystyle=$ $\displaystyle 4\pi
p_{0}-\frac{1}{4}U.$ (38)

Discussing the analytical conditions on $p_{0}$ and $\Phi_{0}$ for Eq. (37) to
have real solutions would be too involved. Equation (37) is practically solved
numerically. We just point out that these conditions diminish the parameter
space for spontaneous scalarization, causing abrupt stops in $\Phi_{0}$ as
$p_{0}$ varies as shown in Fig. 4.

We would also like to point out that the center behavior of the linearized
scalar field can be obtained from Eqs. (36) as a byproduct. Dropping all the
$\Phi$-related quantities in the first and the second equations in (36)
recovers $\mu_{2}$ and $\nu_{2}$ in GR, namely

$\displaystyle\mu_{2}=\frac{4\pi}{3}\epsilon_{0},\quad\nu_{2}=\frac{2\pi}{3}\left(\epsilon_{0}+3p_{0}\right).$
(39)

Then, the third equation in (36) gives $\Phi_{2}$ in the linearized case

$\displaystyle\Phi_{2}=\left[\frac{128\pi^{2}\xi}{9}\epsilon_{0}\left(\epsilon_{0}+3p_{0}\right)+\frac{1}{6\lambdabar_{\Phi}^{2}}\right]\Phi_{0},$
(40)

where $dF/d\Phi=2\xi\Phi_{0}$ and $dU/d\Phi=2\Phi_{0}/\lambdabar_{\Phi}^{2}$
have been used. As the recurrence relation is linear, the linearized scalar
field inside the star is fixed up to a scale factor for a given NS
configuration in GR.

##  Appendix B Asymptotic expansion of the spherical solution

It is illustrative to first study the asymptotic expansion of the linearized
scalar field by only considering the third equation in Eqs. (12) with $\mu$
and $\nu$ given by the Schwarzschild metric,

$\displaystyle\mu=-\nu=-\frac{1}{2}\ln{\left(1-\frac{2M}{r}\right)},$ (41)

where $M$ is the Arnowitt-Deser-Misner mass of the star. The ansatz

$\displaystyle\Phi=e^{\pm r/\lambdabar_{\Phi}}\Phi_{\pm}=e^{\pm
r/\lambdabar_{\Phi}}\,r^{s_{\pm}}\sum_{n=1}^{\infty}\left(\Phi_{\pm}\right)_{n}\frac{1}{r^{n}},$
(42)

is to be used to find series solutions for $\Phi$, where the constants
$s_{\pm}$ and the expansion coefficients $\left(\Phi_{\pm}\right)_{n}$ are
expected to be determined once the ansatz is substituted into the linearized
scalar equation together with the series expansions of $\mu$ and $\nu$ in Eq.
(41). Indeed, the leading-order expansion of the linearized scalar equation
provides

$\displaystyle\left(s_{\pm}\mp\frac{M}{\lambdabar_{\Phi}}\right)(\Phi_{\pm})_{1}=0,$
(43)

implying $s_{\pm}=\pm M/\lambdabar_{\Phi}$, while higher-order expansion terms
give equations generating $(\Phi_{\pm})_{n}$ from $(\Phi_{\pm})_{1}$:

$\displaystyle(\Phi_{\pm})_{2}$ $\displaystyle=$
$\displaystyle\frac{M}{2}\left(1-3s_{\pm}\right)(\Phi_{\pm})_{1},$
$\displaystyle(\Phi_{\pm})_{3}$ $\displaystyle=$
$\displaystyle\frac{M^{2}}{8}\frac{9s_{\pm}^{3}-16s_{\pm}^{2}+5s_{\pm}-2}{s_{\pm}}(\Phi_{\pm})_{1},$
(44) $\displaystyle\cdots$ $\displaystyle\,.$

A general solution for $\Phi$ is a linear combination of
$e^{r/\lambdabar_{\Phi}}\Phi_{+}$ and $e^{-r/\lambdabar_{\Phi}}\Phi_{-}$. But
physical solutions require $\Phi$ to be suppressed at infinity, meaning that
$\Phi_{+}$ has to be zero. This leaves us one free coefficient
$(\Phi_{-})_{1}$ in the series solution for $\Phi$ as $(\Phi_{+})_{1}$ must be
set to zero. Since the linearized scalar equation is homogeneous, the free
coefficient $(\Phi_{-})_{1}$ is merely a scale factor apart from which the
exterior solution for $\Phi$ is completely fixed for a star with a given value
of mass $M$. A suitable value of the NS mass at which the exterior solution of
$\Phi$ connects smoothly to the interior solution of $\Phi$ marks a possible
onset of spontaneous scalarization.

Now the series expansion in Eq. (42) is ready to be generalized to the
solution for the full set of ODEs in Eqs. (12). The nonlinearity inevitably
leads to terms containing $e^{\pm 2r/\lambdabar_{\Phi}},\,e^{\pm
3r/\lambdabar_{\Phi}},\,e^{\pm 4r/\lambdabar_{\Phi}},\,\cdots$, so the
physical solution can be formally written as

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle\mu^{(0)}+e^{-r/\lambdabar_{\Phi}}r^{s_{-}}\mu^{(1)}+e^{-2r/\lambdabar_{\Phi}}r^{2s_{-}}\mu^{(2)}+\cdots\,,$
$\displaystyle\nu$ $\displaystyle=$
$\displaystyle\nu^{(0)}+e^{-r/\lambdabar_{\Phi}}r^{s_{-}}\nu^{(1)}+e^{-2r/\lambdabar_{\Phi}}r^{2s_{-}}\nu^{(2)}+\cdots\,,$
$\displaystyle\Phi$ $\displaystyle=$ $\displaystyle
e^{-r/\lambdabar_{\Phi}}r^{s_{-}}\Phi^{(1)}+e^{-2r/\lambdabar_{\Phi}}r^{2s_{-}}\Phi^{(2)}+\cdots\,,$
(45)

where $\mu^{(i)},\,\nu^{(i)}$ and $\Phi^{(i)}$ have Taylor expansions in terms
of $1/r$. When the ansatz (45) is substituted into Eqs. (12), terms containing
the same exponential factor are collected to generate equations by the orders
of $e^{-r/\lambdabar_{\Phi}}$. Starting with the zeroth order of
$e^{-r/\lambdabar_{\Phi}}$, we acquire the equations for $\mu^{(0)}$ and
$\nu^{(0)}$ that are identical to the GR equations for $\mu$ and $\nu$ while
the scalar equation trivially vanishes. Next, at the first order of
$e^{-r/\lambdabar_{\Phi}}$, the equations for $\mu^{(1)},\,\nu^{(1)}$ and
$\Phi^{(1)}$ come out. The equation for $\Phi^{(1)}$ is identical to the
linearized scalar equation. Setting up their Taylor expansions in terms of
$1/r$ and substituting into their equations, we find recurrence equations that
fix the expansion coefficients of $\mu^{(1)}$ and $\nu^{(1)}$ to zero. Similar
procedures apply to orders at
$e^{-2r/\lambdabar_{\Phi}},\,e^{-3r/\lambdabar_{\Phi}},\,\cdots$, so the
Taylor expansions for higher-order $\mu^{(i)},\nu^{(i)}$ and $\Phi^{(i)}$ can
be found.

We present the first few terms of the result:

$\displaystyle g_{tt}$ $\displaystyle=$
$\displaystyle-1+\frac{2M}{r}+\left[\frac{(\Phi_{-})_{1}}{M}\right]^{2}e^{-2r/\lambdabar_{\Phi}}r^{2s_{-}}\left(\frac{1}{4}\frac{M^{2}}{r^{2}}+\cdots\right)+\cdots\,,$
$\displaystyle g_{rr}$ $\displaystyle=$ $\displaystyle
1+\frac{2M}{r}+\frac{4M^{2}}{r^{2}}+\frac{8M^{3}}{r^{3}}+\cdots+\left[\frac{(\Phi_{-})_{1}}{M}\right]^{2}e^{-2r/\lambdabar_{\Phi}}r^{2s_{-}}\left(\frac{s_{-}}{2}\frac{M}{r}-\frac{3s_{-}^{2}-4s_{-}+1}{2}\frac{M^{2}}{r^{2}}+\cdots\right)+\cdots\,,$
$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\frac{(\Phi_{-})_{1}}{M}e^{-r/\lambdabar_{\Phi}}r^{s_{-}}\left(\frac{M}{r}+\frac{1-3s_{-}}{2}\frac{M^{2}}{r^{2}}+\frac{9s_{-}^{3}-16s_{-}^{2}+5s_{-}-2}{8s_{-}}\frac{M^{3}}{r^{3}}+\cdots\right)+\left[\frac{(\Phi_{-})_{1}}{M}\right]^{3}e^{-3r/\lambdabar_{\Phi}}r^{3s_{-}}\left(\frac{s_{-}}{8}\frac{M^{2}}{r^{2}}+\cdots\right)+\cdots\,.$
(46)

Here $M$ and $(\Phi_{-})_{1}$ are the only two free parameters describing the
asymptotic behavior of the solution. The constant $s_{-}=-M/\lambdabar_{\Phi}$
inherits the notation from the linearized case. The first few terms shown in
Eqs. (46) provide a decent approximation for the solution when $r\gg M$ and
$r\gg\lambdabar_{\Phi}$. By assuming that $r$ is much greater than the larger
one of $M$ and $\lambdabar_{\Phi}$, the asymptotic behaviors of $g_{tt}$ and
$g_{rr}$ are basically the same as those in GR since the modifications from
the scalar field are proportional to the exponential factors and therefore are
dropped. In astrophysical applications, the modifications from the scalar
field at the far-field region however might not be dropped when
$\lambdabar_{\Phi}\gg M$ so the suppression from the exponential factors is
inefficient in the region $M\ll r\ll\lambdabar_{\Phi}$. Noticing that the
leading term with the exponential factor $e^{-2r/\lambdabar_{\Phi}}$ in
$g_{tt}$ roughly has the same magnitude as the leading term of $\Phi^{2}$, a
criterion to determine whether the scalar field affects the far-field metric
significantly is checking if $\Phi$ decreases fast enough so that in the
problematic region $M\ll r\ll\lambdabar_{\Phi}$ its value is much less than
$\sqrt{M/r}$ or not.

At the end of the discussion on the asymptotic behavior of the scalarized
solution to Eqs. (12), we especially would like to point out that by setting
$U=0$ in Eqs. (12) and employing Taylor expansions for $\mu,\,\nu$ and $\Phi$
in terms of $1/r$, their asymptotic expansions in the massless-scalar theory
can be obtained. They turn out to be

$\displaystyle g_{tt}$ $\displaystyle=$
$\displaystyle-1+\frac{2M}{r}-\frac{Q_{s}^{2}}{6}\frac{M^{3}}{r^{3}}+\cdots\,,$
$\displaystyle g_{rr}$ $\displaystyle=$ $\displaystyle
1+\frac{2M}{r}+\left(4-Q_{s}^{2}\right)\frac{M^{2}}{r^{2}}+\left(8-\frac{5Q_{s}^{2}}{2}\right)\frac{M^{3}}{r^{3}}+\cdots\,,$
$\displaystyle\Phi$ $\displaystyle=$ $\displaystyle
Q_{s}\left(\frac{M}{r}+\frac{M^{2}}{r^{2}}+\left(\frac{4}{3}-\frac{Q_{s}^{2}}{12}\right)\frac{M^{3}}{r^{3}}+\cdots\right)\,,$
(47)

where $M$ and $Q_{s}$ are the two free parameters. It is surprising to see
that Eqs. (47) are not the limits of Eqs. (46) at
$\lambdabar_{\Phi}\rightarrow\infty$. In fact, Eqs. (46) do not have limits at
$\lambdabar_{\Phi}\rightarrow\infty$, because these infinite series contain
increasing powers of $s_{-}$ in the denominators that unfortunately give rise
to infinities unable to be regularized by redefining the free parameter
$(\Phi_{-})_{1}$ as $\lambdabar_{\Phi}\rightarrow\infty$. The discontinuity
between the solutions of the massive scalar and of the massless scalar
resembles the famous van Dam-Veltman-Zakharov discontinuity regarding the mass
of graviton van Dam and Veltman (1970); Zakharov (1970). Apart from that, we
have to admit that a plausible physical explanation for this puzzle as well as
its rigorous mathematical resolution, is currently missing.

##  Appendix C Equations for even-parity perturbations with $l\geq 2$

In this appendix, we present the equations for the radial functions
$H_{0},\,H_{1},\,H_{2},\,K,\,W_{1},\,V$ and $X$ with the time dependence
assumed to be $e^{i\sigma t}$ for the angular momentum number $l\geq 2$. When
$l=0$ or $l=1$, certain components of the Regge-Wheeler gauge become trivial,
suggesting that the freedom of choosing coordinates is incompletely utilized.
Further coordinate transformations can eliminate one (for $l=1$) or two (for
$l=0$) functions among $H_{0},\,H_{1},\,H_{2},\,K,\,W_{1}$ and $V$, leading to
different sets of equations from Eqs. (54) displayed below (e.g., see Refs.
Chandrasekhar (1964); Campolattaro and Thorne (1970) for the GR case). Here we
only focus on the case of $l\geq 2$ where the Regge-Wheeler gauge completely
fixes the coordinates at the order of the perturbation.

First, the Eulerian perturbation of the 4-velocity is found to be

$\displaystyle\delta u^{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{-\nu}H_{0}Y_{lm}\,e^{i\sigma t},$
$\displaystyle\delta u^{r}$ $\displaystyle=$ $\displaystyle i\sigma
e^{-\nu-2\mu}W_{1}Y_{lm}\,e^{i\sigma t},$ $\displaystyle\delta u^{\theta}$
$\displaystyle=$ $\displaystyle i\sigma
e^{-\nu}V\frac{\partial_{\theta}Y_{lm}}{r^{2}}\,e^{i\sigma t},$
$\displaystyle\delta u^{\phi}$ $\displaystyle=$ $\displaystyle i\sigma
e^{-\nu}V\frac{\partial_{\phi}Y_{lm}}{r^{2}\sin^{2}\theta}\,e^{i\sigma t},$
(48)

by substituting Eqs. (21) into Eqs. (24). The covariant components

$\displaystyle\delta u_{t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{\nu}H_{0}Y_{lm}\,e^{i\sigma t},$
$\displaystyle\delta u_{r}$ $\displaystyle=$ $\displaystyle
e^{-\nu}\left(H_{1}+i\sigma W_{1}\right)Y_{lm}\,e^{i\sigma t},$
$\displaystyle\delta u_{\theta}$ $\displaystyle=$ $\displaystyle i\sigma
e^{-\nu}V\partial_{\theta}Y_{lm}\,e^{i\sigma t},$ $\displaystyle\delta
u_{\phi}$ $\displaystyle=$ $\displaystyle i\sigma
e^{-\nu}V\partial_{\phi}Y_{lm}\,e^{i\sigma t},$ (49)

are then calculated using

$\displaystyle\delta u_{\alpha}=g^{(0)}_{\alpha\beta}\delta u^{\beta}+\delta
g_{\alpha\beta}u^{(0)\beta},$ (50)

where $g^{(0)}_{\alpha\beta}={\rm
diag}\Big{\\{}-e^{2\nu},\,e^{2\mu},\,r^{2},\,r^{2}\sin^{2}\theta\Big{\\}}$ and
$u^{(0)\beta}=\left(e^{-\nu},0,0,0\right)$ are the zeroth-order quantities.

Next, the Lagrangian perturbation of the particle number density is found to
be

$\displaystyle\frac{\Delta
n}{n}=\left(\frac{1}{2}H_{2}+K-\frac{l(l+1)}{r^{2}}V-\frac{e^{-\mu}}{r^{2}}W^{\prime}\right)Y_{lm}e^{i\sigma
t},$ (51)

where $W:=-r^{2}e^{-\mu}W_{1}$ is the fluid displacement function used by
Thorne and Campolattaro (1967) and the prime denotes the derivative with
respect to $r$, by substituting Eqs. (21) into the second equation in Eqs.
(23). Using Eqs. (25), the Eulerian perturbations of the energy density and
the pressure are then

$\displaystyle\delta\epsilon$ $\displaystyle=$
$\displaystyle(\epsilon+p)\frac{\Delta
n}{n}-e^{-2\mu}W_{1}\epsilon^{\prime}Y_{lm}e^{i\sigma t},$
$\displaystyle\delta p$ $\displaystyle=$ $\displaystyle\gamma p\frac{\Delta
n}{n}-e^{-2\mu}W_{1}p^{\prime}Y_{lm}e^{i\sigma t}.$ (52)

Finally, the perturbative field equations can be calculated using the
expressions of $\delta g_{\alpha\beta},\,\delta u_{\alpha},\delta\epsilon,\
\delta p$ and $\delta\Phi$. For the modified Einstein equations, the
perturbation

$\displaystyle\delta E_{\mu\nu}:=\delta G_{\mu\nu}-8\pi(\delta
T_{m})_{\mu\nu}-(\delta T_{\Phi})_{\mu\nu}-(\delta T_{\rm GB})_{\mu\nu},$ (53)

can be grouped into 4 scalars, 2 vectors, and 1 traceless tensor, providing 7
equations. Together with the perturbation of the scalar field equation, there
are 8 equations in total for the 7 variables
$H_{0},\,H_{1},\,H_{2},\,K,\,W_{1},\,V$ and $X$, indicating that one equation
depends on the others. Please note that $W_{0}$ never shows up in the field
equations. This is expected as $\xi_{t}$ is unnecessary in describing the
motion of the fluid.

The following is a set of 7 equations from which
$H_{0},\,H_{1},\,H_{2},\,K,\,W_{1},\,V$ and $X$ can be solved:

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{1}\left[1-2r\mu^{\prime}+\frac{(l+2)(l-1)}{2}e^{2\mu}+8\pi
r^{2}e^{2\mu}\epsilon+\frac{1}{2}r^{2}e^{2\mu}U+\frac{1}{2}r^{2}\Phi^{\prime\,2}+4\left(1-e^{-2\mu}\right)F^{\prime\prime}-\frac{2l(l+1)}{r}F^{\prime}-4\left(1-3e^{-2\mu}\right)\mu^{\prime}F^{\prime}\right]$
$\displaystyle+i\sigma
H_{2}\left[re^{2\mu}+2\left(e^{2\mu}-3\right)F^{\prime}\right]-i\sigma
K^{\prime}(r^{2}e^{2\mu}-4rF^{\prime})+i\sigma
K\left(r\nu^{\prime}-1\right)\left(re^{2\mu}-4F^{\prime}\right)+8\pi i\sigma
W_{1}r^{2}e^{2\mu}(\epsilon+p)$ $\displaystyle+4i\sigma
X^{\prime}(1-e^{2\mu})\frac{dF}{d\Phi}+i\sigma
X\left[4(1-e^{2\mu})\Phi^{\prime}\frac{d^{2}F}{d\Phi^{2}}-4(1-e^{2\mu})\nu^{\prime}\frac{dF}{d\Phi}-r^{2}e^{2\mu}\Phi^{\prime}\right],$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{1}^{\prime}\left(re^{2\mu}-4F^{\prime}\right)+H_{1}\left[re^{2\mu}(\nu^{\prime}-\mu^{\prime})+4\left(3F^{\prime}\mu^{\prime}-F^{\prime}\nu^{\prime}+\frac{1}{r}F^{\prime}-F^{\prime\prime}\right)\right]-i\sigma
H_{2}e^{2\mu}\left(re^{2\mu}-4F^{\prime}\right)$ $\displaystyle-i\sigma
Kre^{2\mu}\left(e^{2\mu}-4F^{\prime\prime}+4F^{\prime}\mu^{\prime}\right)+16\pi
i\sigma Vre^{4\mu}(\epsilon+p)-8i\sigma Xe^{2\mu}\mu^{\prime}\frac{dF}{d\Phi}$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{0}^{\prime}\left[r+2\left(1-3e^{-2\mu}\right)F^{\prime}\right]-\frac{l(l+1)}{2}H_{0}\frac{1}{r}\left(re^{2\mu}-4F^{\prime}\right)+H_{2}\left[e^{2\mu}+4\pi
r^{2}e^{2\mu}(2-\gamma)p-\frac{1}{2}r^{2}e^{2\mu}U-12e^{-2\mu}\nu^{\prime}F^{\prime}\right]$
$\displaystyle-K^{\prime}r\left(1+r\nu^{\prime}-12e^{-2\mu}\nu^{\prime}F^{\prime}\right)-K\left[\sigma^{2}r^{2}e^{2(\mu-\nu)}-\frac{(l+2)(l-1)}{2}e^{2\mu}+4\pi
r^{2}e^{2\mu}\gamma
p-4\sigma^{2}re^{-2\nu}F^{\prime}+2(l+2)(l-1)\nu^{\prime}F^{\prime}\right]$
$\displaystyle-2i\sigma
H_{1}e^{-2\nu}\left[r+2\left(1-3e^{-2\mu}\right)F^{\prime}\right]-4\pi
W_{1}^{\prime}r^{2}\gamma p-4\pi W_{1}r\left(rp^{\prime}+2\gamma
p-r\mu^{\prime}\gamma p\right)+4\pi l(l+1)Ve^{2\mu}\gamma p$
$\displaystyle+X^{\prime}\left[r^{2}\Phi^{\prime}+4\left(1-3e^{-2\mu}\right)\nu^{\prime}\frac{dF}{d\Phi}\right]+X\left[4\sigma^{2}e^{-2\nu}\left(1-e^{2\mu}\right)\frac{dF}{d\Phi}-\frac{4l(l+1)}{r}\nu^{\prime}\frac{dF}{d\Phi}-4\left(1-3e^{-2\mu}\right)\nu^{\prime}\Phi^{\prime}\frac{d^{2}F}{d\Phi^{2}}-\frac{1}{2}\frac{dU}{d\Phi}r^{2}e^{2\mu}\right]$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{0}^{\prime}\left(re^{2\mu}-4F^{\prime}\right)+H_{0}\frac{1}{r}\left(r\nu^{\prime}-1\right)\left(re^{2\mu}-4F^{\prime}\right)-i\sigma
H_{1}e^{-2\nu}\left(re^{2\mu}-4F^{\prime}\right)+H_{2}\left[(1+r\nu^{\prime})e^{2\mu}-12\nu^{\prime}F^{\prime}\right]$
$\displaystyle-K^{\prime}\left(re^{2\mu}-4r\nu^{\prime}F^{\prime}\right)+8X^{\prime}\nu^{\prime}\frac{dF}{d\Phi}+2X\left(4\nu^{\prime}\Phi^{\prime}\frac{d^{2}F}{d\Phi^{2}}-4\frac{\nu^{\prime}}{r}\frac{dF}{d\Phi}-re^{2\mu}\Phi^{\prime}\right)$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{0}\left(e^{2\mu}-4F^{\prime\prime}+4\mu^{\prime}F^{\prime}\right)-H_{2}\left(e^{2\mu}-4\nu^{\prime}F^{\prime}\right)+8X\left(\nu^{\prime\prime}+\nu^{\prime\,2}-\mu^{\prime}\nu^{\prime}\right)\frac{dF}{d\Phi},$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
H_{0}r^{2}(\epsilon+p)+H_{2}r^{2}\gamma p+2Kr^{2}\gamma
p+2W_{1}^{\prime}r^{2}e^{-2\mu}\gamma
p+2W_{1}re^{-2\mu}\left(rp^{\prime}+2\gamma p-r\mu^{\prime}\gamma
p\right)+2V\left[\sigma^{2}r^{2}e^{-2\nu}(\epsilon+p)-l(l+1)\gamma p\right],$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle
2H_{0}^{\prime\prime}\left(1-e^{-2\mu}\right)\frac{dF}{d\Phi}+H_{0}^{\prime}\left[4\left(1-e^{-2\mu}\right)\nu^{\prime}\frac{dF}{d\Phi}-2\left(1-3e^{-2\mu}\right)\mu^{\prime}\frac{dF}{d\Phi}-\frac{1}{2}r^{2}\Phi^{\prime}\right]-2l(l+1)H_{0}\frac{\mu^{\prime}}{r}\frac{dF}{d\Phi}$
(54) $\displaystyle-4i\sigma
H_{1}^{\prime}e^{-2\nu}\left(1-e^{-2\mu}\right)\frac{dF}{d\Phi}+i\sigma
H_{1}e^{-2\nu}\left(r^{2}\Phi^{\prime}+4\left(1-3e^{-2\mu}\right)\mu^{\prime}\frac{dF}{d\Phi}\right)+H_{2}^{\prime}\left[2\left(1-3e^{-2\mu}\right)\nu^{\prime}\frac{dF}{d\Phi}-\frac{1}{2}r^{2}\Phi^{\prime}\right]$
$\displaystyle+H_{2}\left[4\left(1-2e^{-2\mu}\right)\nu^{\prime\prime}+4\left(1-2e^{-2\mu}\right)\nu^{\prime\,2}-4\left(1-6e^{-2\mu}\right)\mu^{\prime}\nu^{\prime}-\frac{2l(l+1)}{r}\nu^{\prime}-2\sigma^{2}\left(e^{2\mu}-1\right)\right]\frac{dF}{d\Phi}$
$\displaystyle+H_{2}\left[r\left(r\mu^{\prime}-r\nu^{\prime}-2\right)\Phi^{\prime}-r^{2}\Phi^{\prime\prime}\right]+4K^{\prime\prime}re^{-2\mu}\nu^{\prime}\frac{dF}{d\Phi}+K^{\prime}\left[r^{2}\Phi^{\prime}+4e^{-2\mu}\left(r\nu^{\prime\prime}+r\nu^{\prime\,2}+2\nu^{\prime}-3r\mu^{\prime}\nu^{\prime}\right)\frac{dF}{d\Phi}\right]$
$\displaystyle-2K\left[2\sigma^{2}re^{-2\nu}\mu^{\prime}+(l+2)(l-1)\left(\nu^{\prime\prime}+\nu^{\prime\,2}-\mu^{\prime}\nu^{\prime}\right)\right]\frac{dF}{d\Phi}+X^{\prime\prime}r^{2}+X^{\prime}r\left(r\nu^{\prime}-r\mu^{\prime}+2\right)$
$\displaystyle+X\left[\sigma^{2}r^{2}e^{2(\mu-\nu)}-l(l+1)e^{2\mu}+4\left(1-3e^{-2\mu}\right)\mu^{\prime}\nu^{\prime}\frac{d^{2}F}{d\Phi^{2}}-4\left(1-e^{-2\mu}\right)(\nu^{\prime\prime}+\nu^{\prime\,2})\frac{d^{2}F}{d\Phi^{2}}-\frac{1}{2}r^{2}e^{2\mu}\frac{d^{2}U}{d\Phi^{2}}\right].$

The first 5 equations are extracted from $\delta E_{tr},\,\delta
E_{t\theta},\,\delta E_{rr},\,\delta E_{r\theta}$ and $\delta E_{\theta\phi}$.
The penultimate equation is extracted from the $\theta$-component of the
conservation equation (9), rather than from $\delta E_{tt}$ or $\delta
E_{\theta\theta}+\delta E_{\phi\phi}/\sin^{2}\theta$ as they are even more
lengthy. The last equation is from the perturbative scalar field equation. We
notice that all the equations are independent of the angular momentum number
$m$ while some of them depend on $l$. The condition $l\geq 2$ should be
addressed again here, because $l\geq 1$ has to be assumed to extract the
equations from the vector components $\delta E_{t\theta}$ and $\delta
E_{r\theta}$, and $l\geq 2$ has to be assumed to extract the equation from the
tensor component $\delta E_{\theta\phi}$.

Among the 7 perturbative functions, $H_{1},\,H_{2}$ and $V$ can be solved in
terms of other functions using the first, the fifth, and the sixth equations
in Eqs. (54). So the system can be reduced to 4 ODEs for $H_{0},\,K,\,W_{1}$
and $X$, with the highest derivatives involved being
$H_{0}^{\prime},\,K^{\prime},\,W_{1}^{\prime\prime}$ and $X^{\prime\prime}$.
The boundary conditions imposed to get physical solutions are analogous to
those in GR. At the center of the star, the functions have Taylor expansions
taking the following form

$\displaystyle H_{0}$ $\displaystyle=$ $\displaystyle
H_{00}r^{l}+H_{02}r^{l+2}+\cdots\,,$ $\displaystyle H_{1}$ $\displaystyle=$
$\displaystyle H_{10}r^{l+1}+H_{12}r^{l+3}\cdots\,,$ $\displaystyle H_{2}$
$\displaystyle=$ $\displaystyle H_{20}r^{l}+H_{22}r^{l+2}\cdots\,,$
$\displaystyle K$ $\displaystyle=$ $\displaystyle
K_{0}r^{l}+K_{2}r^{l+2}\cdots\,,$ $\displaystyle W_{1}$ $\displaystyle=$
$\displaystyle W_{10}r^{l-1}+W_{12}r^{l+1}\cdots\,,$ $\displaystyle V$
$\displaystyle=$ $\displaystyle V_{0}r^{l}+V_{2}r^{l+2}\cdots\,,$
$\displaystyle X$ $\displaystyle=$ $\displaystyle
X_{0}r^{l}+X_{2}r^{l+2}\cdots\,,$ (55)

where only three of the coefficients are free. For example, using
$H_{00},\,W_{10}$ and $X_{0}$ as the free coefficients, we can find

$\displaystyle H_{10}$ $\displaystyle=$
$\displaystyle\frac{2i\sigma}{l+1}\left[H_{00}\left(1-8\Phi_{2}\frac{dF}{d\Phi}\right)+16X_{0}\nu_{2}\frac{dF}{d\Phi}\right]$
$\displaystyle+\frac{16i\sigma}{l+1}\left[\frac{X_{0}\mu_{2}\frac{dF}{d\Phi}}{1-8\Phi_{2}\frac{dF}{d\Phi}}-\frac{\pi
W_{10}\left(\epsilon_{0}+p_{0}\right)}{l\left(1-8\Phi_{2}\frac{dF}{d\Phi}\right)}\right]\,,$
$\displaystyle H_{20}$ $\displaystyle=$ $\displaystyle
K_{0}=H_{00}\left(1-8\Phi_{2}\frac{dF}{d\Phi}\right)+16X_{0}\nu_{2}\frac{dF}{d\Phi}\,,$
$\displaystyle V_{0}$ $\displaystyle=$ $\displaystyle\frac{W_{10}}{l},$ (56)

where $\mu_{2},\,\nu_{2}$ and $\Phi_{2}$ are solutions from Eqs. (36) and
$dF/d\Phi$ is understood to take its value at the center of the star, as well
as higher-order recurrence equations to generate the higher-order
coefficients. At infinity, the solutions for the perturbative functions
feature the factors $e^{\pm i\sigma r}$ while those with the factor
$e^{i\sigma r}$ need to be discarded to fulfill the boundary condition of no
incoming waves taking into consideration of the assumed time-dependence
$e^{i\sigma t}$. An extra boundary condition at infinity compared to GR comes
from the fact that the factors $e^{\pm r/\lambdabar_{\Phi}},\,e^{\pm
2r/\lambdabar_{\Phi}},\,\cdots$, appear in the solutions for the perturbative
functions, where $\lambdabar_{\Phi}$ is the reduced Compton wavelength of the
scalar field, so those with the factors
$e^{r/\lambdabar_{\Phi}},\,e^{2r/\lambdabar_{\Phi}},\,\cdots$, need to be
discarded for the perturbation to be consistent with the spherical solution.
Last but not least, at the surface of the star every physical solution has a
vanishing Lagrangian perturbation of the particle number density, leaving a
condition on the combination of the functions $H_{2},\,K,\,V$ and $W_{1}$ as
shown in Eq. (51).

For the numerical approach that integrates from the center of the star to
solve Eqs. (54), the use of the Taylor expansions (55) automatically meets the
boundary condition at the center. To satisfy the boundary conditions at the
surface of the star and at infinity, the free parameters
$H_{00},\,W_{10},\,X_{0}$ as well as the complex frequency $\sigma$ need to be
adjusted to suitable values. Noticing that Eqs. (54) are homogeneous for the
set of the perturbative functions, it is the ratios between
$H_{00},\,W_{10},\,X_{0}$ that matter. After taking $H_{00}=1$, the shooting
strategy can be implemented for finding suitable values of $W_{10}$ and
$X_{0}$, aiming at vanishing $\Delta n$ at the surface of the star and the
divergent part in $\delta\Phi$ at infinity caused by the factors of
$e^{r/\lambdabar_{\Phi}}$. The suitable values of $\sigma$, namely the
quasinormal frequencies, are to be found by shooting for eliminating the
$e^{i\sigma r}$ part of the solutions in principle, but the process of trial
and error is much elaborate due to the divergence in the solutions at infinity
caused by the imaginary part of $\sigma$ Lindblom and Detweiler (1983);
Detweiler and Lindblom (1985); Chandrasekhar and Ferrari (1991); Andersson
_et al._ (1993).

For the static case studied in this work, the solutions are much simplified by
setting $\sigma=0$. The first equation in Eqs. (54) gives $H_{1}=0$, and
remarkably, $W_{1}$ and $W_{1}^{\prime}$ can be eliminated altogether with $V$
using the sixth equation, leaving us two first-order ODEs for $H_{0}$ and $K$,
and a second-order ODE for $X$. We only need to find suitable values of one
parameter, namely $X_{0}$, using the shooting strategy. In fact, following
Pani and Berti (2014), due to the linearity of Eqs. (54) in terms of the
perturbative functions, instead of applying the shooting strategy two
independent solutions with $H_{00}=1,\,X_{0}=0$ and $H_{00}=1,\,X_{0}=1$ are
computed first, and the physical solution with $X$ vanishing at infinity can
be conveniently constructed as a linear combination of the two computed
solutions.

##  Appendix D Calculating the tidal deformability using the formula of
Hinderer

Using the asymptotic expressions (46) for the spherical solution we can obtain
the following asymptotic expansion for $H_{0}$ when $\sigma=0$ and $l=2$:

$\displaystyle H_{0}$ $\displaystyle\rightarrow$ $\displaystyle
a_{P}\left(-\frac{3r^{2}}{M^{2}}+\frac{6r}{M}\right)+a_{Q}\left(-\frac{8}{5}\frac{M^{3}}{r^{3}}-\frac{24}{5}\frac{M^{4}}{r^{4}}+\cdots\right)+\left[\frac{\left(\Phi_{-}\right)_{1}}{M}\right]^{2}e^{-2r/\lambdabar_{\Phi}}r^{2s_{-}}\left[a_{P}\left(-\frac{s_{-}}{4}\frac{r}{M}+\cdots\right)+a_{Q}\left(\frac{s_{-}}{5}\frac{M^{4}}{r^{4}}+\cdots\right)\right]+\cdots,$
(57)

where $\left(\Phi_{-}\right)_{1}$ is the coefficient in the spherical solution
for $\Phi$ and $s_{-}=-M/\lambdabar_{\Phi}$. To match the leading terms as
shown in Eq. (29), the coefficients $a_{P}$ and $a_{Q}$ are

$\displaystyle a_{P}=\frac{M^{2}}{3}\frac{{\mathcal{E}}_{2m}}{A_{m}},\quad
a_{Q}=-\frac{15}{8M^{3}}\frac{Q_{2m}}{A_{m}},$ (58)

so the tidal deformability can be calculated using $a_{P}$ and $a_{Q}$ as

$\displaystyle\lambda=\frac{8}{45}\frac{a_{Q}}{a_{P}}M^{5}.$ (59)

Similarly to the spherical solution, the expansion for $H_{0}$ is simply the
same as that in GR after dropping the terms involving the factors of
$e^{-r/\lambdabar_{\Phi}}$. Since for static perturbations, $H_{0}$ in GR
outside the star has the analytical solution

$\displaystyle\left(H_{0}\right)_{\rm
GR}=a_{P}P^{2}_{l}\left(\frac{r}{M}-1\right)+a_{Q}{\rm
Re}\left[Q^{2}_{l}\left(\frac{r}{M}-1\right)\right],$ (60)

with $P^{2}_{l}$ and $Q^{2}_{l}$ being the associated Legendre functions, the
coefficients $a_{P}$ and $a_{Q}$ can be calculated once the values of $H_{0}$
and its derivative are given at a sufficiently large radius $r_{\infty}$,
namely

$\displaystyle a_{P}$ $\displaystyle=$ $\displaystyle\frac{H_{0}^{\prime}{\rm
Re}\left[Q^{2}_{l}\right]-H_{0}{\rm
Re}\left[Q^{2}_{l}\right]^{\prime}}{\left(P^{2}_{l}\right)^{\prime}{\rm
Re}\left[Q^{2}_{l}\right]-P^{2}_{l}{\rm
Re}\left[Q^{2}_{l}\right]^{\prime}}\Bigg{|}_{r=r_{\infty}},$ $\displaystyle
a_{Q}$ $\displaystyle=$
$\displaystyle-\frac{H_{0}^{\prime}P^{2}_{l}-H_{0}\left(P^{2}_{l}\right)^{\prime}}{\left(P^{2}_{l}\right)^{\prime}{\rm
Re}\left[Q^{2}_{l}\right]-P^{2}_{l}{\rm
Re}\left[Q^{2}_{l}\right]^{\prime}}\Bigg{|}_{r=r_{\infty}},$ (61)

where the prime denotes the derivative with respect to $r$ and the argument of
the associated Legendre functions is $r/M-1$. Taking $l=2$ and defining
$y:=rH_{0}^{\prime}/H_{0}\big{|}_{r=r_{\infty}}$, $C:=M/r_{\infty}$, the tidal
deformability given by substituting Eqs. (61) into Eq. (59) has exactly the
same form as the formula of Hinderer (2008)

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\frac{16}{15}M^{5}(1-2C)^{2}\left[2Cy-2C+2-y\right]$ (62)
$\displaystyle\times\Big{[}8C^{5}(y+1)+4C^{4}(3y-2)+4C^{3}(13-11y)$
$\displaystyle\hskip 5.0pt+6C^{2}(5y-8)+2C(6-3y)$ $\displaystyle\hskip
5.0pt+3(1-2C)^{2}(2Cy-2C+2-y)\ln\left(1-2C\right)\Big{]}^{-1}.$

Figure 9: $\lambda/M^{5}$ versus $y$ for different values of $C$.

Now supposing $\lambdabar_{\Phi}\gg M$, we take
$r_{\infty}\approx\lambdabar_{\Phi}$ for the exponential term in Eq. (57) to
drop. Then the issue of $C\approx M/\lambdabar_{\Phi}$ being too small is
exposed in using the formula (62). Simply plotting $\lambda$ versus $y$ for
several values of $C$ in Fig. 9 reveals the rather violent change of $\lambda$
when $y$ is near its asymptotic value $2$ and $C\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}0.01$. Remarkably,
the error in $\lambda$ caused by an inaccurate $y$ with an error $\delta y$
can be estimated by

$\displaystyle\delta\lambda\approx\left|\frac{\partial\lambda}{\partial
y}\right|\delta y\approx\frac{M^{5}}{15C^{5}}\delta y,$ (63)

where $y\approx 2$ has been used. The error in $y$ comes from the error in
computing the solution for $H_{0}$, which can be estimated by the error of the
integrator. For the efficiency of our code, the relative error of the
integrator has been set to $10^{-7}$. Using $\delta y=y\times 10^{-7}\approx
2\times 10^{-7}$, the absolute error in $\lambda$ is, for instance,
$\delta\lambda\approx 0.04\times M^{5}$ when $M/\lambdabar_{\Phi}=0.05$ and
$\delta\lambda\approx 133\times M^{5}$ when $M/\lambdabar_{\Phi}=0.01$. The
mass of NSs is at the order of kilometer in geometrized units. So when
$\lambdabar_{\Phi}\sim 100\,{\rm km}$ ($a=0.1$), the absolute error for the
dimensionless tidal deformability is as high as $10^{2}$, which is
unacceptable if $\lambda/M^{5}$ itself is at the same order of magnitude or
even less.

##  Appendix E Equations for odd-parity perturbations

In this appendix, we present the equations for the radial functions
$h_{0},\,h_{1}$ and $h_{2}$ with the time dependence assumed to be $e^{i\sigma
t}$. Here the angular momentum number $l$ takes positive integers as the odd-
parity spherical harmonics start from $l=1$.

Substituting Eqs. (22) into Eq. (24), the Eulerian perturbation of the
4-velocity turns out to be

$\displaystyle\delta u^{t}$ $\displaystyle=$ $\displaystyle\delta u^{r}=0,$
$\displaystyle\delta u^{\theta}$ $\displaystyle=$ $\displaystyle-i\sigma
e^{-\nu}h_{2}\frac{\partial_{\phi}Y_{lm}}{r^{2}\sin\theta}\,e^{i\sigma t},$
$\displaystyle\delta u^{\phi}$ $\displaystyle=$ $\displaystyle i\sigma
e^{-\nu}h_{2}\frac{\partial_{\theta}Y_{lm}}{r^{2}\sin\theta}\,e^{i\sigma t}.$
(64)

The covariant components are

$\displaystyle\delta u_{t}$ $\displaystyle=$ $\displaystyle\delta u_{r}=0,$
$\displaystyle\delta u_{\theta}$ $\displaystyle=$
$\displaystyle-e^{-\nu}(h_{0}+i\sigma
h_{2})\frac{\partial_{\phi}Y_{lm}}{\sin\theta}\,e^{i\sigma t},$
$\displaystyle\delta u_{\phi}$ $\displaystyle=$ $\displaystyle
e^{-\nu}(h_{0}+i\sigma h_{2})\sin\theta\,\partial_{\theta}Y_{lm}\,e^{i\sigma
t}.$ (65)

Direct calculation using the second equation in Eqs. (23) verifies that the
Lagrangian perturbation of the particle number density vanishes, and hence

$\displaystyle\delta\epsilon=\delta p=0.$ (66)

The perturbative Einstein equations (53) provide 3 equations for
$h_{0},\,h_{1}$ and $h_{2}$ when $l\geq 2$ but only 2 equations when $l=1$,
consistent with the fact that the Regge-Wheeler gauge is trivial for the case
of odd-parity with $l=1$. Using the $\theta$-component of the conservation
equation instead of the tensor component $\delta E_{\theta\phi}$ that is
invalid when $l=1$, we can write down the following three equations:

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
h_{0}^{\prime\prime}\left(r^{2}e^{2\mu}-4rF^{\prime}\right)-h_{0}^{\prime}\left[r^{2}e^{2\mu}(\mu^{\prime}+\nu^{\prime})+4rF^{\prime\prime}-4\left(1+3r\mu^{\prime}+r\nu^{\prime}\right)F^{\prime}\right]$
$\displaystyle-
h_{0}\left[2e^{2\mu}(1-2r\mu^{\prime}+r^{2}\nu^{\prime\prime}+r^{2}\nu^{\prime\,2}-r^{2}\mu^{\prime}\nu^{\prime})+(l+2)(l-1)e^{4\mu}+16\pi
r^{2}e^{4\mu}\epsilon\right]$ $\displaystyle-
h_{0}\left[r^{2}e^{2\mu}\Phi^{\prime\,2}+r^{2}e^{4\mu}U+4(l+2)(l-1)e^{2\mu}(\mu^{\prime}F^{\prime}-F^{\prime\prime})+8\left(\nu^{\prime}+3\mu^{\prime}-r\nu^{\prime\prime}-r\nu^{\prime\,2}+3r\mu^{\prime}\nu^{\prime}\right)F^{\prime}-8(1+r\nu^{\prime})F^{\prime\prime}\right]$
$\displaystyle-i\sigma
h_{1}^{\prime}r\left(re^{2\mu}-4F^{\prime}\right)+i\sigma
h_{1}\left[r^{2}e^{2\mu}(r\mu^{\prime}+r\nu^{\prime}-2)+4\left(1-3r\mu^{\prime}-r\nu^{\prime}\right)F^{\prime}+4rF^{\prime\prime}\right]-16\pi
i\sigma h_{2}r^{2}e^{4\mu}(\epsilon+p),$ $\displaystyle 0$ $\displaystyle=$
$\displaystyle i\sigma
h_{0}^{\prime}\left(r^{2}e^{2\mu}-4rF^{\prime}\right)-2i\sigma
h_{0}\left(re^{2\mu}-4F^{\prime}\right)$
$\displaystyle+h_{1}\left[\sigma^{2}r^{2}e^{2\mu}-(l+2)(l-1)e^{2(\mu+\nu)}-2re^{2\nu}(\nu^{\prime}-\mu^{\prime}+r\nu^{\prime\prime}+r\nu^{\prime\,2}-r\mu^{\prime}\nu^{\prime})+16\pi
r^{2}e^{2(\mu+\nu)}p-r^{2}e^{2\nu}\Phi^{\prime\,2}-r^{2}e^{2(\mu+\nu)}U\right]$
$\displaystyle+4h_{1}\left[-\sigma^{2}r\Phi^{\prime}\frac{dF}{d\Phi}+(l+2)(l-1)e^{2\nu}\nu^{\prime}\Phi^{\prime}\frac{dF}{d\Phi}+2re^{2(\nu-\mu)}(\nu^{\prime\prime}+\nu^{\prime\,2}-3\mu^{\prime}\nu^{\prime})\Phi^{\prime}\frac{dF}{d\Phi}+2re^{2(\nu-\mu)}\nu^{\prime}F^{\prime\prime}\right],$
$\displaystyle 0$ $\displaystyle=$ $\displaystyle\sigma^{2}h_{2}-i\sigma
h_{0},$ (67)

where the first and the second equations come from $\delta E_{t\theta}$ and
$\delta E_{r\theta}$, and the third equation is from the $\theta$-component of
the conservation equation. They are valid for $l$ taking all positive integers
only that one of the equations becomes dependent on the other two when $l=1$.

For $\sigma\neq 0$, a coordinate transformation can be performed at the order
of the perturbation to eliminate one of the perturbative functions in the case
of $l=1$. Choosing it to be $h_{1}=0$, then the second equation in Eqs. (67)
immediately gives $h_{2}\propto r^{2}$, which does not contain the oscillating
factors $e^{\pm i\sigma r}$, implying no odd-parity dipole gravitational wave
in the scalar-tensor theory. When $l\geq 2$, the Regge-Wheeler gauge
completely fixes the coordinates and the three independent equations in Eqs.
(67) are sufficient to solve the three variables $h_{0},\,h_{1}$ and $h_{2}$.
The boundary conditions for physical solutions include the finite behavior at
the center of the star described by the expansions

$\displaystyle h_{0}$ $\displaystyle=$ $\displaystyle
h_{00}r^{l+1}+h_{02}r^{l+3}+\cdots,$ $\displaystyle h_{1}$ $\displaystyle=$
$\displaystyle h_{10}r^{l+2}+h_{12}r^{l+4}+\cdots,$ (68)

as well as the requirement of no incoming waves at infinity. Note that there
is only one free coefficient in the expansions (68). Specifically speaking,
$h_{00}$ and $h_{10}$ are related via

$\displaystyle(l+2)h_{10}=i\sigma
e^{-2\nu}\left(1-8\Phi_{2}\frac{dF}{d\Phi}\right)h_{00},$ (69)

where $\Phi_{2}$ is solved from Eqs. (36), and $\nu$ and $dF/d\Phi$ are
understood to take their values at the center of the star, while the higher-
order coefficients are generated from $h_{00}$ and $h_{10}$ by recurrence
equations. Since Eqs. (67) are homogeneous in the perturbative functions, the
free coefficient in the expansions (68) is merely a scale factor. Therefore
the only parameter that requires suitable values for physical solutions is the
frequency $\sigma$, which should be selected to eliminate the incoming wave at
infinity.

For $\sigma=0$, the case of $l=1$ is not much different from the case of
$l\geq 2$, because the perturbative Einstein equations only provide two
independent equations even when $l\geq 2$ as long as $\sigma=0$. This is
conveniently revealed by the fact that the third equation in Eqs. (67)
trivially vanishes. The problem is further simplified as the second equation
in Eqs. (67) gives $h_{1}=0$. In conclusion, the physical solution to the
first equation in Eqs. (67) turns out to describe the only nonvanishing off-
diagonal metric component

$\displaystyle g_{t\phi}=h_{0}\sin\theta\,\partial_{\theta}Y_{l0},$ (70)

of a slowly rotating star. In fact, writing the 4-velocity of a slowly
rotating star as

$\displaystyle u^{\mu}=\frac{1}{\sqrt{-g_{tt}-2\Omega
g_{t\phi}-\Omega^{2}g_{\phi\phi}}}\left(1,0,0,\Omega\right),$ (71)

and expanding it up to the linear order of the angular velocity $\Omega$, the
perturbation of the 4-velocity is found to be

$\displaystyle\delta u^{\alpha}=\left(0,0,0,e^{-\nu}\Omega\right),$ (72)

which matches Eqs. (64) by taking $m=0$ and

$\displaystyle i\sigma
h_{2}\frac{\partial_{\theta}Y_{l0}}{r^{2}\sin\theta}\rightarrow\Omega,$ (73)

as $\sigma\rightarrow 0$. To solve $h_{0}$, the angular velocity $\Omega$
needs to be specified. But it cannot be specified arbitrarily due to the match
(73). In general it takes the form

$\displaystyle\Omega=\tilde{\Omega}\frac{\partial_{\theta}Y_{l0}}{\sin\theta},$
(74)

where $\tilde{\Omega}$ is a given function of $r$ and has the expansion

$\displaystyle\tilde{\Omega}=\tilde{\Omega}_{0}r^{l-1}+\tilde{\Omega}_{2}r^{l+1}+\tilde{\Omega}_{4}r^{l+3}+\cdots,$
(75)

at the center of the star. The solution of $h_{0}$ then still has an expansion
in the form of the first equation in Eqs. (68) at the center of the star, now
with the coefficients
$\tilde{\Omega}_{0},\,\tilde{\Omega}_{2},\,\tilde{\Omega}_{4},\,\cdots$,
entering into the recurrence equations for $h_{0n}$. For example, $h_{02}$ is
related to $h_{00}$ and $\tilde{\Omega}_{0}$ via

$\displaystyle h_{02}$ $\displaystyle=$
$\displaystyle\frac{\left[4\pi\epsilon_{0}(l+1)+2\pi
p_{0}(l+3)+\left(l+\frac{1}{2}\right)(l-1)\mu_{2}\right]h_{00}}{(2l+3)(1-8\Phi_{2}\frac{dF}{d\Phi})}$
(76)
$\displaystyle-\frac{(l-1)\left[16(l+1)\Phi_{4}\frac{dF}{d\Phi}+8(l+1)\Phi_{2}^{2}\frac{d^{2}F}{d\Phi^{2}}-\frac{1}{8}U\right]h_{00}}{(2l+3)(1-8\Phi_{2}\frac{dF}{d\Phi})}$
$\displaystyle+\frac{8\pi(\epsilon_{0}+p_{0})\tilde{\Omega}_{0}}{(2l+3)(1-8\Phi_{2}\frac{dF}{d\Phi})},$

where $p_{0},\,\epsilon_{0},\,\mu_{2},\,\Phi_{2}$ and $\Phi_{4}$ are the
expansion coefficients defined in Eqs. (34), and
$\frac{dF}{d\Phi},\,\frac{d^{2}F}{d\Phi^{2}}$ and $U$ are understood to take
their values at the center of the star. Differently from the case of
$\sigma\neq 0$, now there are two free parameters, namely $h_{00}$ and
$\tilde{\Omega}_{0}$, and their ratio determines the solution of $h_{0}$,
which in general consists of a part whose leading asymptotic term is
proportional to $r^{l+1}$ and a part whose leading asymptotic term is
proportional to $1/r^{l}$. For $g_{t\phi}$ to vanish at infinity, a suitable
value of $h_{00}/\tilde{\Omega}_{0}$ is demanded so that the $r^{l+1}$ part
does not show up.

## References

  * Stelle (1977) K. S. Stelle, Phys. Rev. D 16, 953 (1977).
  * Starobinsky (1980) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * Copeland  _et al._ (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Nojiri and Odintsov (2011) S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59 (2011), arXiv:1011.0544 [gr-qc] .
  * Lanczos (1938) C. Lanczos, Annals Math. 39, 842 (1938).
  * Cai and Guo (2004) R.-G. Cai and Q. Guo, Phys. Rev. D 69, 104025 (2004), arXiv:hep-th/0311020 .
  * Konoplya (2005) R. Konoplya, Phys. Rev. D 71, 024038 (2005), arXiv:hep-th/0410057 .
  * Cai  _et al._ (2005) R.-G. Cai, H.-S. Zhang, and A. Wang, Commun. Theor. Phys. 44, 948 (2005), arXiv:hep-th/0505186 .
  * Camanho and Edelstein (2010) X. O. Camanho and J. D. Edelstein, JHEP 04, 007 (2010), arXiv:0911.3160 [hep-th] .
  * Zwiebach (1985) B. Zwiebach, Phys. Lett. B 156, 315 (1985).
  * Boulware and Deser (1985) D. G. Boulware and S. Deser, Phys. Rev. Lett. 55, 2656 (1985).
  * Glavan and Lin (2020) D. Glavan and C. Lin, Phys. Rev. Lett. 124, 081301 (2020), arXiv:1905.03601 [gr-qc] .
  * Torii and Maeda (2005) T. Torii and H. Maeda, Phys. Rev. D 71, 124002 (2005), arXiv:hep-th/0504127 .
  * Nojiri  _et al._ (2005) S. Nojiri, S. D. Odintsov, and M. Sasaki, Phys. Rev. D 71, 123509 (2005), arXiv:hep-th/0504052 .
  * Boulware and Deser (1986) D. G. Boulware and S. Deser, Phys. Lett. B 175, 409 (1986).
  * Koivisto and Mota (2007) T. Koivisto and D. F. Mota, Phys. Rev. D 75, 023518 (2007), arXiv:hep-th/0609155 .
  * Cognola  _et al._ (2007) G. Cognola, E. Elizalde, S. Nojiri, S. Odintsov, and S. Zerbini, Phys. Rev. D 75, 086002 (2007), arXiv:hep-th/0611198 .
  * Leith and Neupane (2007) B. M. Leith and I. P. Neupane, JCAP 05, 019 (2007), arXiv:hep-th/0702002 .
  * Akbarieh  _et al._ (2021) A. R. Akbarieh, S. Kazempour, and L. Shao, Phys. Rev. D 103, 123518 (2021), arXiv:2105.03744 [gr-qc] .
  * Scherk and Schwarz (1974) J. Scherk and J. H. Schwarz, Nuclear Physics B 81, 118 (1974).
  * Kanti  _et al._ (1996) P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, Phys. Rev. D 54, 5049 (1996), arXiv:hep-th/9511071 .
  * Torii  _et al._ (1997) T. Torii, H. Yajima, and K.-i. Maeda, Phys. Rev. D 55, 739 (1997), arXiv:gr-qc/9606034 .
  * Kanti  _et al._ (1998) P. Kanti, N. E. Mavromatos, J. Rizos, K. Tamvakis, and E. Winstanley, Phys. Rev. D 57, 6255 (1998), arXiv:hep-th/9703192 .
  * Torii and Maeda (1998) T. Torii and K.-i. Maeda, Phys. Rev. D 58, 084004 (1998).
  * Pani and Cardoso (2009) P. Pani and V. Cardoso, Phys. Rev. D 79, 084031 (2009), arXiv:0902.1569 [gr-qc] .
  * Kanti  _et al._ (2011) P. Kanti, B. Kleihaus, and J. Kunz, Phys. Rev. Lett. 107, 271101 (2011), arXiv:1108.3003 [gr-qc] .
  * Ayzenberg and Yunes (2014) D. Ayzenberg and N. Yunes, Phys. Rev. D 90, 044066 (2014), [Erratum: Phys.Rev.D 91, 069905 (2015)], arXiv:1405.2133 [gr-qc] .
  * Cunha  _et al._ (2017) P. V. P. Cunha, C. A. R. Herdeiro, B. Kleihaus, J. Kunz, and E. Radu, Phys. Lett. B 768, 373 (2017), arXiv:1701.00079 [gr-qc] .
  * Ayzenberg and Yunes (2018) D. Ayzenberg and N. Yunes, Class. Quant. Grav. 35, 235002 (2018), arXiv:1807.08422 [gr-qc] .
  * Gonzalez  _et al._ (2019) G. A. Gonzalez, B. Kleihaus, J. Kunz, and S. Mojica, Phys. Rev. D 99, 024041 (2019), arXiv:1812.02686 [gr-qc] .
  * Blázquez-Salcedo  _et al._ (2016) J. L. Blázquez-Salcedo, C. F. B. Macedo, V. Cardoso, V. Ferrari, L. Gualtieri, F. S. Khoo, J. Kunz, and P. Pani, Phys. Rev. D 94, 104024 (2016), arXiv:1609.01286 [gr-qc] .
  * Konoplya  _et al._ (2019) R. A. Konoplya, A. F. Zinhailo, and Z. Stuchlík, Phys. Rev. D 99, 124042 (2019), arXiv:1903.03483 [gr-qc] .
  * Pierini and Gualtieri (2021) L. Pierini and L. Gualtieri, Phys. Rev. D 103, 124017 (2021), arXiv:2103.09870 [gr-qc] .
  * Yagi  _et al._ (2016) K. Yagi, L. C. Stein, and N. Yunes, Phys. Rev. D 93, 024010 (2016), arXiv:1510.02152 [gr-qc] .
  * Witek  _et al._ (2019) H. Witek, L. Gualtieri, P. Pani, and T. P. Sotiriou, Phys. Rev. D 99, 064035 (2019), arXiv:1810.05177 [gr-qc] .
  * Wang  _et al._ (2021) H.-T. Wang, S.-P. Tang, P.-C. Li, M.-Z. Han, and Y.-Z. Fan, Phys. Rev. D 104, 024015 (2021), arXiv:2104.07590 [gr-qc] .
  * Sotiriou and Zhou (2014) T. P. Sotiriou and S.-Y. Zhou, Phys. Rev. Lett. 112, 251102 (2014), arXiv:1312.3622 [gr-qc] .
  * Saffer and Yagi (2021) A. Saffer and K. Yagi, (2021), arXiv:2110.02997 [gr-qc] .
  * Damour and Esposito-Farèse (1993) T. Damour and G. Esposito-Farèse, Phys. Rev. Lett. 70, 2220 (1993).
  * Damour and Esposito-Farèse (1996) T. Damour and G. Esposito-Farèse, Phys. Rev. D 54, 1474 (1996), arXiv:gr-qc/9602056 [gr-qc] .
  * Harada (1998) T. Harada, Phys. Rev. D 57, 4802 (1998), arXiv:gr-qc/9801049 .
  * Salgado  _et al._ (1998) M. Salgado, D. Sudarsky, and U. Nucamendi, Phys. Rev. D 58, 124003 (1998), arXiv:gr-qc/9806070 .
  * Ramazanoğlu and Pretorius (2016) F. M. Ramazanoğlu and F. Pretorius, Phys. Rev. D 93, 064005 (2016), arXiv:1601.07475 [gr-qc] .
  * Mendes and Ortiz (2016) R. F. P. Mendes and N. Ortiz, Phys. Rev. D 93, 124035 (2016), arXiv:1604.04175 [gr-qc] .
  * Herdeiro  _et al._ (2018) C. A. R. Herdeiro, E. Radu, N. Sanchis-Gual, and J. A. Font, Phys. Rev. Lett. 121, 101102 (2018), arXiv:1806.05190 [gr-qc] .
  * Cunha  _et al._ (2019) P. V. P. Cunha, C. A. R. Herdeiro, and E. Radu, Phys. Rev. Lett. 123, 011101 (2019), arXiv:1904.09997 [gr-qc] .
  * Antoniou  _et al._ (2018) G. Antoniou, A. Bakopoulos, and P. Kanti, Phys. Rev. Lett. 120, 131102 (2018), arXiv:1711.03390 [hep-th] .
  * Doneva and Yazadjiev (2018a) D. D. Doneva and S. S. Yazadjiev, Phys. Rev. Lett. 120, 131103 (2018a), arXiv:1711.01187 [gr-qc] .
  * Silva  _et al._ (2018) H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou, and E. Berti, Phys. Rev. Lett. 120, 131104 (2018), arXiv:1711.02080 [gr-qc] .
  * Doneva and Yazadjiev (2018b) D. D. Doneva and S. S. Yazadjiev, JCAP 04, 011 (2018b), arXiv:1712.03715 [gr-qc] .
  * Kuan  _et al._ (2021) H.-J. Kuan, D. D. Doneva, and S. S. Yazadjiev, Phys. Rev. Lett. 127, 161103 (2021), arXiv:2103.11999 [gr-qc] .
  * Antoniou  _et al._ (2021) G. Antoniou, A. Lehébel, G. Ventagli, and T. P. Sotiriou, Phys. Rev. D 104, 044002 (2021), arXiv:2105.04479 [gr-qc] .
  * Ventagli  _et al._ (2021) G. Ventagli, G. Antoniou, A. Lehébel, and T. P. Sotiriou, (2021), arXiv:2111.03644 [gr-qc] .
  * Wex (2014) N. Wex, in  _Frontiers in Relativistic Celestial Mechanics: Applications and Experiments_ , Vol. 2, edited by S. M. Kopeikin (Walter de Gruyter GmbH, Berlin/Boston, 2014) p. 39, arXiv:1402.5594 [gr-qc] .
  * Shao and Wex (2016) L. Shao and N. Wex, Sci. China Phys. Mech. Astron. 59, 699501 (2016), arXiv:1604.03662 [gr-qc] .
  * Shao  _et al._ (2017) L. Shao, N. Sennett, A. Buonanno, M. Kramer, and N. Wex, Phys. Rev. X 7, 041025 (2017), arXiv:1704.07561 [gr-qc] .
  * Akiyama  _et al._ (2019) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Hulse and Taylor (1975) R. A. Hulse and J. H. Taylor, Astrophys. J. Lett. 195, L51 (1975).
  * Abbott  _et al._ (2019a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019a), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021), arXiv:2010.14527 [gr-qc] .
  * Danchev  _et al._ (2021) V. I. Danchev, D. D. Doneva, and S. S. Yazadjiev, (2021), arXiv:2112.03869 [gr-qc] .
  * Shiralilou  _et al._ (2021a) B. Shiralilou, T. Hinderer, S. Nissanke, N. Ortiz, and H. Witek, Phys. Rev. D 103, L121503 (2021a), arXiv:2012.09162 [gr-qc] .
  * Shiralilou  _et al._ (2021b) B. Shiralilou, T. Hinderer, S. Nissanke, N. Ortiz, and H. Witek, (2021b), arXiv:2105.13972 [gr-qc] .
  * Kovács and Reall (2020) A. D. Kovács and H. S. Reall, Phys. Rev. Lett. 124, 221101 (2020), arXiv:2003.04327 [gr-qc] .
  * Bryant  _et al._ (2021) A. Bryant, H. O. Silva, K. Yagi, and K. Glampedakis, Phys. Rev. D 104, 044051 (2021), arXiv:2106.09657 [gr-qc] .
  * Staykov  _et al._ (2021) K. V. Staykov, J. L. Blázquez-Salcedo, D. D. Doneva, J. Kunz, P. Nedkova, and S. S. Yazadjiev, (2021), arXiv:2112.00703 [gr-qc] .
  * Mendes and Ortiz (2018) R. F. P. Mendes and N. Ortiz, Phys. Rev. Lett. 120, 201104 (2018), arXiv:1802.07847 [gr-qc] .
  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957).
  * Thorne and Campolattaro (1967) K. S. Thorne and A. Campolattaro, Astrophys. J. 149, 591 (1967).
  * Price and Thorne (1969) R. Price and K. S. Thorne, Astrophys. J. 155, 163 (1969).
  * Thorne (1969) K. S. Thorne, Astrophys. J. 158, 1 (1969).
  * Silva and Yunes (2019) H. O. Silva and N. Yunes, Phys. Rev. D 99, 044034 (2019), arXiv:1808.04391 [gr-qc] .
  * Xu  _et al._ (2020) R. Xu, Y. Gao, and L. Shao, Phys. Rev. D 102, 064057 (2020), arXiv:2007.10080 [gr-qc] .
  * Hu  _et al._ (2021) Z. Hu, Y. Gao, R. Xu, and L. Shao, Phys. Rev. D 104, 104014 (2021), arXiv:2109.13453 [gr-qc] .
  * Hu  _et al._ (2020) H. Hu, M. Kramer, N. Wex, D. J. Champion, and M. S. Kehl, Mon. Not. Roy. Astron. Soc. 497, 3118 (2020), arXiv:2007.07725 [astro-ph.SR] .
  * Miao and Li (2021) Z. Miao and A. Li, (2021), arXiv:2107.07979 [astro-ph.HE] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018a), arXiv:1805.11581 [gr-qc] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Thorne (1980) K. S. Thorne, Rev. Mod. Phys. 52, 299 (1980).
  * Gittins  _et al._ (2020) F. Gittins, N. Andersson, and J. P. Pereira, Phys. Rev. D 101, 103025 (2020), arXiv:2003.05449 [astro-ph.HE] .
  * Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Pani and Berti (2014) P. Pani and E. Berti, Phys. Rev. D 90, 024025 (2014), arXiv:1405.4547 [gr-qc] .
  * Abbott  _et al._ (2019b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019b), arXiv:1805.11579 [gr-qc] .
  * Abbott  _et al._ (2018b) B. P. Abbott _et al._ (KAGRA, LIGO Scientific, Virgo, VIRGO), Living Rev. Rel. 21, 3 (2018b), arXiv:1304.0670 [gr-qc] .
  * Yagi and Yunes (2017) K. Yagi and N. Yunes, Phys. Rept. 681, 1 (2017), arXiv:1608.02582 [gr-qc] .
  * Gao  _et al._ (2022) Y. Gao, X.-Y. Lai, L. Shao, and R.-X. Xu, Mon. Not. Roy. Astron. Soc. 509, 2758 (2022), arXiv:2109.13234 [gr-qc] .
  * Shao  _et al._ (2015) L. Shao _et al._ , PoS AASKA14, 042 (2015), arXiv:1501.00058 [astro-ph.HE] .
  * Weltman  _et al._ (2020) A. Weltman _et al._ , Publ. Astron. Soc. Austral. 37, e002 (2020), arXiv:1810.02680 [astro-ph.CO] .
  * van Dam and Veltman (1970) H. van Dam and M. J. G. Veltman, Nucl. Phys. B 22, 397 (1970).
  * Zakharov (1970) V. I. Zakharov, JETP Lett. 12, 312 (1970).
  * Chandrasekhar (1964) S. Chandrasekhar, Astrophys. J. 140, 417 (1964), [Erratum: Astrophys.J. 140, 1342 (1964)].
  * Campolattaro and Thorne (1970) A. Campolattaro and K. S. Thorne, Astrophys. J. 159, 847 (1970).
  * Lindblom and Detweiler (1983) L. Lindblom and S. L. Detweiler, Astrophys. J. Suppl. 53, 73 (1983).
  * Detweiler and Lindblom (1985) S. L. Detweiler and L. Lindblom, Astrophys. J. 292, 12 (1985).
  * Chandrasekhar and Ferrari (1991) S. Chandrasekhar and V. Ferrari, Proc. Roy. Soc. Lond. A 432, 247 (1991).
  * Andersson  _et al._ (1993) N. Andersson, M. E. Araujo, and B. F. Schutz, Classical and Quantum Gravity 10, 735 (1993).