# Stability of neutron stars in Horndeski theories with Gauss-Bonnet couplings

Masato Minamitsuji1 and Shinji Tsujikawa2 1Centro de Astrofísica e Gravitação
- CENTRA, Departamento de Física, Instituto Superior Técnico - IST,
Universidade de Lisboa - UL, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal
2Department of Physics, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo
169-8555, Japan

###### Abstract

In Horndeski theories containing a scalar coupling with the Gauss-Bonnet (GB)
curvature invariant $R_{\rm GB}^{2}$, we study the existence and linear
stability of neutron star (NS) solutions on a static and spherically symmetric
background. For a scalar-GB coupling of the form $\alpha\xi(\phi)R_{\rm
GB}^{2}$, where $\xi$ is a function of the scalar field $\phi$, the existence
of linearly stable stars with a nontrivial scalar profile without
instabilities puts an upper bound on the strength of the dimensionless
coupling constant $|\alpha|$. To realize maximum masses of NSs for a linear
(or dilatonic) GB coupling $\alpha_{\rm GB}\phi R_{\rm GB}^{2}$ with typical
nuclear equations of state, we obtain the theoretical upper limit
$\sqrt{|\alpha_{\rm GB}|}<0.7~{}{\rm km}$. This is tighter than those obtained
by the observations of gravitational waves emitted from binaries containing
NSs. We also incorporate cubic-order scalar derivative interactions, quartic
derivative couplings with nonminimal couplings to a Ricci scalar besides the
scalar-GB coupling and show that NS solutions with a nontrivial scalar profile
satisfying all the linear stability conditions are present for certain ranges
of the coupling constants. In regularized 4-dimensional Einstein-GB gravity
obtained from a Kaluza-Klein reduction with an appropriate rescaling of the GB
coupling constant, we find that NSs in this theory suffer from a strong
coupling problem as well as Laplacian instability of even-parity
perturbations. We also study NS solutions with a nontrivial scalar profile in
power-law $F(R_{\rm GB}^{2})$ models, and show that they are pathological in
the interior of stars and plagued by ghost instability together with the
asymptotic strong coupling problem in the exterior of stars.

††preprint: WUCG-22-05

##  I Introduction

After the dawn of gravitational waves (GW) astronomy from a binary system of
black holes (BHs) Abbott  _et al._ (2016), a new observational probe of the
physics on a strong gravitational background has begun. In particular, the
GW170817 event Abbott  _et al._ (2017) allowed us to put constraints on the
mass-radius relation of neutron stars (NSs) from their tidal deformation
before the coalescence Abbott  _et al._ (2018). After a merger of binaries, a
compact object exhibits a damped sinusoidal oscillation with quasi-normal
frequencies Kokkotas and Schmidt (1999); Nollert (1999). Upcoming
observational data of GWs will provide us further detailed information for new
physics in strong gravity regimes in the vicinity of BHs and in the interior
of NSs.

From cosmological observational data, we know that about 95 % of the energy
density of today’s Universe is dominated by mysterious components dubbed dark
energy and dark matter Riess  _et al._ (1998); Perlmutter  _et al._ (1999);
Spergel  _et al._ (2003); Tegmark  _et al._ (2004); Ade  _et al._ (2014).
Since it is still challenging to explain their origins within the framework of
General Relativity (GR) and Standard Model of particle physics, there is a
motivation to introduce new degrees of freedom Copeland  _et al._ (2006);
Silvestri and Trodden (2009); De Felice and Tsujikawa (2010); Clifton  _et
al._ (2012); Joyce  _et al._ (2015); Koyama (2016); Heisenberg (2019); Kase
and Tsujikawa (2019a). A scalar field is one of the simplest candidates, and
widely used to explain physical phenomena relevant to the dark sector of the
Universe.

On a spherically symmetric background, an asymptotically-flat vacuum solution
in GR is uniquely described by the Schwarzschild metric with the mass of a
compact body. The background Schwarzschild geometry can be modified by
introducing a new degree of freedom. For asymptotically-flat BHs, however, the
property of the absence of a nontrivial scalar profile holds for a wide class
of scalar-tensor theories–including a canonical scalar field Hawking (1972a);
Bekenstein (1972), k-essence Graham and Jha (2014), nonminimally coupled
scalar field with Ricci scalar Hawking (1972b); Bekenstein (1995); Sotiriou
and Faraoni (2012); Faraoni (2017), and regular derivative couplings in shift-
symmetric Horndeski theories Hui and Nicolis (2013). An exceptional case is a
scalar coupling with the GB curvature invariant $R_{\rm GB}^{2}$ of the form
$\alpha\xi(\phi)R_{\rm GB}^{2}$, where $\alpha$ is a dimensionless coupling
constant and $\xi$ is a regular function of the scalar field $\phi$ Kanti  _et
al._ (1996); Torii  _et al._ (1997); Kanti  _et al._ (1998); Chen  _et al._
(2007); Guo  _et al._ (2008, 2009); Pani and Cardoso (2009); Sotiriou and Zhou
(2014a, b); Ayzenberg and Yunes (2014); Maselli  _et al._ (2015); Kleihaus
_et al._ (2011, 2016); Doneva and Yazadjiev (2018a); Silva  _et al._ (2018);
Antoniou  _et al._ (2018); Blázquez-Salcedo  _et al._ (2018); Minamitsuji and
Ikeda (2019); Silva  _et al._ (2019); Macedo  _et al._ (2019); Doneva and
Yazadjiev (2022); Minamitsuji  _et al._ (2022a). This scalar-GB coupling
belongs to a subclass of Horndeski theories Horndeski (1974) containing
nonanalytic functions like $\ln|X|$ in the coupling functions $G_{2,3,4,5}$
Kobayashi  _et al._ (2011), where $X$ is a canonical field kinetic term.

For a NS, the presence of matter inside the star can modify the no-hair
property of BHs in several subclasses of Horndeski theories. Nonminimal
couplings with the Ricci scalar $R$ of the form $G_{4}(\phi)R$ allow a
possibility for realizing NS solutions endowed with a nontrivial scalar
profile, especially via the mechanism of spontaneous scalarization Damour and
Esposito-Farese (1993); Harada (1998); Novak (1998); Sotani and Kokkotas
(2004); Cooney  _et al._ (2010); Arapoglu  _et al._ (2011); Orellana  _et al._
(2013); Aparicio Resco  _et al._ (2016); Kase and Tsujikawa (2019b). This
includes Brans-Dicke theories Brans and Dicke (1961) and $f(R)$ gravity
Starobinsky (1980), where the latter corresponds to a particular class of the
former with a scalar potential O’Hanlon (1972); Chiba (2003). These
nonminimally coupled theories belong to a subclass of non-shift-symmetric
Horndeski theories.

In shift-symmetric subclass of Horndeski theories where the field equations of
motion are invariant under the shift $\phi\to\phi+c$, there is a no-hair
argument of stars under several assumptions Lehébel  _et al._ (2017) analogous
to BHs discussed in Ref. Hui and Nicolis (2013). The assumptions are as
follows:

  1. (i)

the scalar field and metrics are regular, static, and spherically symmetric
with an asymptotically-flat spacetime geometry,

  2. (ii)

a canonical kinetic term $X$ is present in the action,

  3. (iii)

the action is analytic with regular coupling functions $G_{2,3,4,5}$.

Under these hypotheses, we end up with a no-hair solution $\phi={\rm
constant}$.

One way of breaking the assumption (i) is to postulate a scalar field of the
form $\phi=qt+\psi(r)$, where $q$ is a nonvanishing constant, $t$ and $r$ are
time and radial coordinates, respectively. Originally, this type of field
configuration was considered to search for BH solutions with a nontrivial
profile of the scalar field in shift-symmetric Horndeski theories Babichev and
Charmousis (2014). The analysis was further extended to relativistic stars in
DHOST theories Babichev  _et al._ (2016); Sakstein  _et al._ (2017); Kobayashi
and Hiramatsu (2018); Ogawa  _et al._ (2020); Ikeda  _et al._ (2021). In this
paper, we do not consider such a time-dependent background scalar field and
focus on the case $q=0$, i.e., the static scalar field.

If we break the assumption (ii), i.e., no canonical kinetic term in the
action, it is known that the quartic nonminimal derivative coupling
$G_{4}\supset\mu_{4}X$ gives rise to a solution endowed with a nontrivial
scalar profile inside the star Cisterna  _et al._ (2015). However, it was
recognized that this solution with a nontrivial scalar profile is plagued by
an angular Laplacian instability of even-parity perturbations for large
multipoles around the surface of star Kase and Tsujikawa (2021, 2022).

Finally, one can break the above assumption (iii) by introducing nonanalytic
coupling functions. In the presence of $X$, the quintic-order coupling
$G_{5}=-4\alpha\ln|X|$, which is equivalent to the linear scalar-GB coupling
$\alpha\phi R_{\rm GB}^{2}$, gives rise to NS solutions endowed with a
nontrivial scalar profile Maselli  _et al._ (2016). This linear coupling can
also accommodate a dilatonic coupling $\xi(\phi)\propto{\rm e}^{\mu\phi}$ in
the limit $|\mu\phi|\ll 1$. For the same scalar-GB coupling, there are also
asymptotically-flat hairy BHs Sotiriou and Zhou (2014a, b) consistent with all
the linear stability conditions against odd- and even-parity perturbations
Minamitsuji  _et al._ (2022b). We note that BH solutions present for other
nonanalytic functions in $G_{2,3,4}$ Babichev  _et al._ (2017) are either
unstable around the horizon or asymptotically non-flat Creminelli  _et al._
(2020); Minamitsuji  _et al._ (2022b). For the linear scalar-GB coupling, it
is not yet clear whether NS solutions with a nontrivial scalar profile satisfy
all the stability conditions against odd- and even-parity perturbations. For
this purpose, we can exploit conditions for the absence of ghost/Laplacian
instabilities recently derived in full Horndeski theories Kase and Tsujikawa
(2022) (see also Refs. Kobayashi  _et al._ (2012, 2014); Kase  _et al._
(2020)). NSs can have a nontrivial scalar profile in more general GB couplings
$\alpha\xi(\phi)R_{\rm GB}^{2}$ with a canonical scalar kinetic term Pani  _et
al._ (2011); Kleihaus  _et al._ (2014); Silva  _et al._ (2018); Doneva and
Yazadjiev (2018b); Blázquez-Salcedo  _et al._ (2016); Olmo  _et al._ (2020).
This belongs to a subclass of non-shift-symmetric Horndeski theories. In this
paper, for such general GB couplings, we will first study the background NS
solutions and their linear stability both analytically and numerically. We
show that the strength of dimensionless coupling constant $|\alpha|$ has an
upper bound to ensure the existence of NSs with a nontrivial scalar profile
free from instabilities around the center of star. Besides
$\alpha\xi(\phi)R_{\rm GB}^{2}$, we also incorporate regular coupling
functions like $G_{3}\supset\mu_{3}X$, $G_{4}\supset\mu_{4}X$, and
$G_{4}\supset\lambda_{4}\phi$ as the representative cases, and study their
effects on the existence and stability of solutions. Provided the coupling
constants are in certain ranges, NS solutions with a nontrivial scalar profile
consistent with all the linear stability conditions are still present in such
combined theories.

Moreover, there are also several other gravitational theories containing the
GB term in the action. In so-called 4-dimensional-Einstein-GB (4DEGB) gravity
Glavan and Lin (2020), the contribution of the GB term in spacetime dimensions
$D$ higher than 4 can be extracted by rescaling the GB coupling constant as
$\alpha\to\alpha/(D-4)$. If we perform a Kaluza-Klein reduction on a flat
internal space whose volume is characterized by the scalar field $\phi$, the
effective 4DEGB theory after the rescaling of $\alpha$ belongs to a subclass
of shift-symmetric Horndeski theories containing the linear GB coupling but
without the canonical scalar kinetic term Lu and Pang (2020); Kobayashi (2020)
(see also Ref. Fernandes  _et al._ (2020); Hennigar  _et al._ (2020) for a
conformal regularization equivalent to the Kaluza-Klein reduction). It is
known that NS solutions with a nontrivial scalar profile are present in the
regularized 4DEGB theory Doneva and Yazadjiev (2021), but we will show that
they are plagued by a strong coupling problem and Laplacian instability of
even-parity perturbations. In theories given by the Lagrangian $R+F(R_{\rm
GB}^{2})$ Nojiri and Odintsov (2005); De Felice and Hindmarsh (2007); Li  _et
al._ (2007); De Felice and Tsujikawa (2009a, b), where $F$ is a positive
power-law function of $R_{\rm GB}^{2}$, we will also show that a ghost
instability and strong coupling at spatial infinity arise for NS solutions
with a nontrivial scalar profile. In these 4DEGB and $F(R_{\rm GB}^{2})$
theories, there are no canonical kinetic terms in the action, which implies
the unhealthy propagation of scalar field perturbations, as it also happens in
derivative coupling theories without the canonical kinetic term Kase and
Tsujikawa (2021, 2022).

This paper is organized as follows. In Sec. II, we present the linear
stability conditions for relativistic stars on the static and spherically
symmetric background. In Sec. III, we derive solutions expanded around the
center of star and at spatial infinity in the presence of scalar-GB couplings
$\alpha\xi(\phi)R_{\rm GB}^{2}$. We give a new theoretical bound on $\alpha$
for the existence of NS solutions with a nontrivial scalar profile free from
instabilities and confirm it numerically for the linear scalar-GB coupling. In
Sec. IV, we implement several regular coupling functions besides the scalar-GB
coupling and explore the parameter space of coupling constants in which there
are NS solutions with a nontrivial scalar profile without ghost or Laplacian
instabilities. In Sec. V, we show the existence of a strong coupling problem
and Laplacian instability for NSs with a nontrivial scalar profile arising in
4DEGB gravity. In Sec. VI, we prove that a non-vanishing scalar-field branch
appearing in power-law $F(R_{\rm GB}^{2})$ models is plagued by ghost and
strong coupling problems at large distances. Sec. VII is devoted to
conclusions.

##  II Background equations and linear stability conditions

We study the existence and stability of NS solutions in Horndeski theories
Horndeski (1974); Deffayet  _et al._ (2011); Kobayashi  _et al._ (2011);
Charmousis  _et al._ (2012), whose action is given by

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\,{\cal L}_{H}+{\cal
S}_{m}(g_{\mu\nu},\Psi_{m})\,,$ (1)

where $g$ is a determinant of the metric tensor $g_{\mu\nu}$, and

$\displaystyle{\cal L}_{H}$ $\displaystyle=$ $\displaystyle
G_{2}(\phi,X)-G_{3}(\phi,X)\square\phi+G_{4}(\phi,X)\,R+G_{4,X}(\phi,X)\left[(\square\phi)^{2}-(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)\right]+G_{5}(\phi,X)G_{\mu\nu}\nabla^{\mu}\nabla^{\nu}\phi$
(2)
$\displaystyle-\frac{1}{6}G_{5,X}(\phi,X)\left[(\square\phi)^{3}-3(\square\phi)\,(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)+2(\nabla^{\mu}\nabla_{\alpha}\phi)(\nabla^{\alpha}\nabla_{\beta}\phi)(\nabla^{\beta}\nabla_{\mu}\phi)\right]\,,$

where the coupling functions $G_{j}$ ($j=2,3,4,5$) depend on the scalar field
$\phi$ and its kinetic term $X=-g^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi/2$,
with the covariant derivative operator $\nabla_{\mu}$. We will use the
notations $\square\phi\equiv\nabla^{\mu}\nabla_{\mu}\phi$ and
$G_{j,\phi}\equiv\partial G_{j}/\partial\phi$, $G_{j,X}\equiv\partial
G_{j}/\partial X$, $G_{j,\phi X}\equiv\partial^{2}G_{j}/(\partial
X\partial\phi)$, and so on. The scalar field $\phi$ and its derivatives are
nonminimally coupled to the Ricci scalar $R$ and Einstein tensor $G_{\mu\nu}$
through the couplings $G_{4}$ and $G_{5}$, respectively. For the matter fields
$\Psi_{m}$ inside NSs, we consider a perfect fluid given by the energy-
momentum tensor

$T_{\mu\nu}=\left(\rho+P\right)u_{\mu}u_{\nu}+Pg_{\mu\nu}\,,$ (3)

where $\rho$ and $P$ are the density and pressure, respectively, and $u_{\mu}$
is the four-velocity of the fluid satisfying the normalization relation
$u_{\mu}u^{\mu}=-1$. Assuming that the perfect fluid is minimally coupled to
gravity, it obeys the continuity equation

$\nabla^{\mu}T_{\mu\nu}=0\,.$ (4)

In terms of the action approach, the perfect fluid can be described by a
Schutz-Sorkin action Schutz and Sorkin (1977); Brown (1993); De Felice  _et
al._ (2010).

###  II.1 Background equations of motion

A static and spherically symmetric background is described by the line element

${\rm d}s^{2}=-f(r){\rm d}t^{2}+h^{-1}(r){\rm d}r^{2}+r^{2}\left({\rm
d}\theta^{2}+\sin^{2}\theta\,{\rm d}\varphi^{2}\right)\,,$ (5)

where $f(r)$ and $h(r)$ are functions of the radial coordinate $r$. On this
background, we consider the scalar field that depends only on the radial
coordinate

$\phi=\phi(r)\,,$ (6)

together with the four-velocity of the fluid $u^{\mu}=[f(r)^{-1/2},0,0,0]$.
Then, the mixed energy-momentum tensor $T^{\mu}{}_{\nu}$ has the following
diagonal components

$T^{\mu}{}_{\nu}={\rm diag}\left[-\rho(r),P(r),P(r),P(r)\right]\,,$ (7)

where $\rho$ and $P$ are functions of $r$ alone. The continuity Eq. (4) gives

$P^{\prime}+\frac{f^{\prime}}{2f}\left(\rho+P\right)=0\,,$ (8)

where a prime represents the derivative with respect to $r$.

The (00), (11), (22) components of gravitational field equations of motion are

$\displaystyle\left(A_{1}+\frac{A_{2}}{r}+\frac{A_{3}}{r^{2}}\right)\phi^{\prime\prime}+\left(\frac{\phi^{\prime}}{2h}A_{1}+\frac{A_{4}}{r}+\frac{A_{5}}{r^{2}}\right)h^{\prime}+A_{6}+\frac{A_{7}}{r}+\frac{A_{8}}{r^{2}}=\rho\,,$
(9)
$\displaystyle-\left(\frac{\phi^{\prime}}{2h}A_{1}+\frac{A_{4}}{r}+\frac{A_{5}}{r^{2}}\right)\frac{hf^{\prime}}{f}+A_{9}-\frac{2\phi^{\prime}}{r}A_{1}-\frac{1}{r^{2}}\left[\frac{\phi^{\prime}}{2h}A_{2}+(h-1)A_{4}\right]=P\,,$
(10)
$\displaystyle\left[\left\\{A_{2}+\frac{(2h-1)\phi^{\prime}A_{3}+2hA_{5}}{h\phi^{\prime}r}\right\\}\frac{f^{\prime}}{4f}+A_{1}+\frac{A_{2}}{2r}\right]\phi^{\prime\prime}+\frac{1}{4f}\left(2hA_{4}-\phi^{\prime}A_{2}+\frac{2hA_{5}-\phi^{\prime}A_{3}}{r}\right)\left(f^{\prime\prime}-\frac{f^{\prime
2}}{2f}\right)$
$\displaystyle+\left[A_{4}+\frac{2h(2h+1)A_{5}-\phi^{\prime}A_{3}}{2h^{2}r}\right]\frac{f^{\prime}h^{\prime}}{4f}+\left(\frac{A_{7}}{4}+\frac{A_{10}}{r}\right)\frac{f^{\prime}}{f}+\left(\frac{\phi^{\prime}}{h}A_{1}+\frac{A_{4}}{r}\right)\frac{h^{\prime}}{2}+A_{6}+\frac{A_{7}}{2r}=-P\,,$
(11)

where the coefficients $A_{1}$-$A_{10}$ are given in Appendix A. The scalar-
field equation of motion is expressed in the form

$\frac{1}{r^{2}}\sqrt{\frac{h}{f}}\left(r^{2}\sqrt{\frac{f}{h}}J^{r}\right)^{\prime}+{\cal
P}_{\phi}=0\,,$ (12)

with

$\displaystyle J^{r}$ $\displaystyle=$ $\displaystyle
h\phi^{\prime}\biggl{[}G_{2,X}-\left(\frac{2}{r}+\frac{f^{\prime}}{2f}\right)h\phi^{\prime}G_{3,X}+2\left(\frac{1-h}{r^{2}}-\frac{hf^{\prime}}{rf}\right)G_{4,X}+2h\phi^{\prime
2}\left(\frac{h}{r^{2}}+\frac{hf^{\prime}}{rf}\right)G_{4,XX}$ (13)
$\displaystyle\qquad-\frac{f^{\prime}}{2r^{2}f}(1-3h)h\phi^{\prime}G_{5,X}-\frac{f^{\prime}h^{3}\phi^{\prime
3}}{2r^{2}f}G_{5,XX}\biggr{]}\,,$ $\displaystyle{\cal P}_{\phi}$
$\displaystyle=$ $\displaystyle
G_{2,\phi}+\lambda_{1}G_{3,\phi}+\lambda_{2}G_{3,\phi\phi}+\lambda_{3}G_{3,\phi
X}+\lambda_{4}G_{4,\phi}+\lambda_{5}G_{4,\phi X}+\lambda_{6}G_{4,\phi\phi
X}+\lambda_{7}G_{4,\phi XX}$ (14)
$\displaystyle+\lambda_{8}G_{5,\phi}+\lambda_{9}G_{5,\phi\phi}+\lambda_{10}G_{5,\phi
X}+\lambda_{11}G_{5,\phi\phi X}+\lambda_{12}G_{5,\phi XX}\,,$

where $\lambda_{1}$–$\lambda_{12}$ are presented in Appendix A. This equation
also follows by combining Eqs. (9)-(11). Note that, in shift-symmetric
Horndeski theories where the coupling functions $G_{j}$ contain the $X$
dependence alone, we have ${\cal P}_{\phi}=0$. In this case, Eq. (12) gives
the solution $J^{r}=(Q/r^{2})\sqrt{h/f}$, where $Q$ is a constant.

###  II.2 Linear stability conditions

To study the linear stability of NS solutions, we consider metric
perturbations $h_{\mu\nu}$ on top of the background (5) besides perturbations
of the scalar field and perfect fluid. Expanding perturbations on the
background (5) in terms of the spherical harmonics of the unit two-sphere
$Y_{lm}(\theta,\varphi)$, one can decompose them into the two different
sectors depending on the parity under the rotation along two-dimensional
sphere Regge and Wheeler (1957); Zerilli (1970). The odd- and even-parity
perturbations have the parities $(-1)^{l+1}$ and $(-1)^{l}$, respectively. Any
scalar perturbation has the even mode alone, whereas vector and tensor
perturbations contain both odd and even modes. The decomposition of
perturbations of metrics, scalar field, and perfect fluids into the odd- and
even-parity modes was addressed in Ref. De Felice  _et al._ (2011); Motohashi
and Suyama (2011); Kobayashi  _et al._ (2012, 2014); Kase  _et al._ (2020);
Kase and Tsujikawa (2022). In the presence of perfect fluids the stability
conditions against odd- and even-parity perturbations were already derived in
Ref. Kase and Tsujikawa (2022), so we briefly summarize them in the following.

In the odd-parity sector, there is a dynamical perturbation $\chi$ arising
from the gravity sector besides a nondynamical perturbation $\delta j$ related
to the $(\theta,\varphi)$ components of fluid four velocity Kase and Tsujikawa
(2022). In the limit of large frequencies and multipoles $l$, the no-ghost
condition for the dynamical field $\chi$ translates to

${\cal G}\equiv 2G_{4}+2h\phi^{\prime 2}G_{4,X}-h\phi^{\prime
2}\left(G_{5,\phi}+\frac{f^{\prime}h\phi^{\prime}G_{5,X}}{2f}\right)>0\,.$
(15)

Under this condition, the Laplacian instability along the radial and angular
directions can be avoided for

$\displaystyle{\cal H}$ $\displaystyle\equiv$ $\displaystyle
2G_{4}+2h\phi^{\prime 2}G_{4,X}-h\phi^{\prime
2}G_{5,\phi}-\frac{h^{2}\phi^{\prime 3}G_{5,X}}{r}>0\,,$ (16)
$\displaystyle{\cal F}$ $\displaystyle\equiv$ $\displaystyle
2G_{4}+h\phi^{\prime 2}G_{5,\phi}-h\phi^{\prime
2}\left(\frac{1}{2}h^{\prime}\phi^{\prime}+h\phi^{\prime\prime}\right)G_{5,X}>0\,,$
(17)

under which the squared propagation speeds $c_{r}^{2}={\cal G}/{\cal F}$ and
$c_{\Omega}^{2}={\cal G}/{\cal H}$ are positive. The expressions (15)-(17)
coincide with those originally derived in Ref. Kobayashi  _et al._ (2012) in
the absence of the perfect fluid.

In the even-parity sector, there are three dynamical perturbations: matter
perturbation $\delta\rho$, gravitational perturbations $\psi$, and scalar-
field perturbation $\delta\phi$. We focus on the linear stability conditions
of high radial and angular momentum modes. With the condition (16), there are
no ghosts for even-parity perturbations if

$\displaystyle\rho+P>0\,,$ (18) $\displaystyle{\cal K}\equiv(2{\cal
P}_{1}-{\cal F})h\mu^{2}-2{\cal H}^{2}r^{4}(\rho+P)>0\,,$ (19)

where

${\cal P}_{1}\equiv\frac{h\mu}{2fr^{2}{\cal H}^{2}}\left(\frac{fr^{4}{\cal
H}^{4}}{\mu^{2}h}\right)^{\prime}\,,\qquad\mu\equiv\frac{2(\phi^{\prime}a_{1}+r\sqrt{fh}{\cal
H})}{\sqrt{fh}}\,.$ (20)

The definition of $a_{1}$ is given in Appendix B.

The radial propagation speed squared of $\psi$ is given by $c_{r2}^{2}={\cal
G}/{\cal F}$, which is the same as that of $\chi$. Thus the gravitational
perturbations in the odd- and even-parity sectors propagate in the same manner
along the radial direction. The radial Laplacian instabilities of $\delta\rho$
and $\delta\phi$ can be avoided for

$\displaystyle c_{m}^{2}\equiv\frac{n\rho_{,nn}}{\rho_{,n}}>0\,,$ (21)
$\displaystyle c_{r3}^{2}\equiv\frac{2\phi^{\prime}[4r^{2}(fh)^{3/2}{\cal
H}c_{4}(2\phi^{\prime}a_{1}+r\sqrt{fh}\,{\cal
H})-2a_{1}^{2}f^{3/2}\sqrt{h}\phi^{\prime}{\cal
G}+(a_{1}f^{\prime}+2c_{2}f)r^{2}fh{\cal H}^{2}]}{f^{5/2}\sqrt{h}\,{\cal
K}}>0\,,$ (22)

where $c_{2}$ and $c_{4}$ are given in Appendix B. The conditions (18) and
(21) correspond to those in the perfect fluid sector. We will consider the
fluid equation of state (EOS) satisfying these inequalities. We note that
$c_{r3}^{2}$ corresponds to the propagation speed squared of scalar-field
perturbation $\delta\phi$.

Along the angular direction, the perfect fluid in the even-parity sector has
the propagation speed squared same as $c_{m}^{2}=n\rho_{,nn}/\rho_{,n}$. The
angular Laplacian instabilities for $\psi$ and $\delta\phi$ are absent if

$c_{\Omega\pm}^{2}=-B_{1}\pm\sqrt{B_{1}^{2}-B_{2}}>0\,,$ (23)

where

$\displaystyle B_{1}\equiv\frac{r^{3}\sqrt{fh}{\cal
H}[4h\beta_{0}\beta_{1}+\beta_{2}-4\phi^{\prime}a_{1}\beta_{3}+r\sqrt{f^{3}h}{\cal
G}{\cal H}(\rho+P)]-2fh{\cal G}[r\sqrt{fh}(2{\cal P}_{1}-{\cal F}){\cal
H}(\beta_{0}+\phi^{\prime}a_{1})+2\phi^{\prime 2}a_{1}^{2}{\cal
P}_{1}]}{4fh(2{\cal P}_{1}-{\cal F}){\cal H}\beta_{0}^{2}-2r^{4}f^{2}h{\cal
H}^{3}(\rho+P)},$ (24) $\displaystyle
B_{2}\equiv-2r^{2}{\frac{r^{2}h\beta_{1}(2fh{\cal F}{\cal
G}\beta_{0}+r^{2}\beta_{2})-{r}^{4}\beta_{2}\beta_{3}-fh{\cal F}{\cal
G}(\phi^{\prime}fh{\cal F}{\cal G}a_{1}+2r^{3}\sqrt{fh}{\cal
H}\beta_{3})}{2fh\phi^{\prime}a_{1}(2{\cal P}_{1}-{\cal F}){\cal
F}\beta_{0}^{2}-r^{4}f^{2}h\phi^{\prime}a_{1}{\cal F}{\cal
H}^{2}(\rho+P)}}\,.$ (25)

The explicit forms of $\beta_{0},\beta_{1},\beta_{2},\beta_{3}$ are presented
in Appendix B. The stability conditions (23) are satisfied if

$B_{1}^{2}\geq B_{2}>0\,,\quad{\rm and}\quad B_{1}<0\,.$ (26)

In summary, we require that NS solutions with a nontrivial scalar profile
should satisfy the inequalities (15), (16), (17), (19), (22), and (23) besides
the fluid stability conditions (18) and (21).

##  III Scalar-Gauss-Bonnet couplings

Let us first consider the Einstein-scalar-GB theory given by the action

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm Pl}^{2}}{2}R+\eta
X+\alpha\xi(\phi)R_{\rm GB}^{2}\right]\,,$ (27)

where $\eta$ is a constant, $M_{\rm Pl}$ is the reduced Planck mass, $\alpha$
is a dimensionless coupling, $\xi$ is a function of $\phi$, and

$R_{\rm GB}^{2}\equiv
R^{2}-4R_{\alpha\beta}R^{\alpha\beta}+R_{\alpha\beta\mu\nu}R^{\alpha\beta\mu\nu}\,,$
(28)

with $R_{\alpha\beta}$ and $R_{\alpha\beta\mu\nu}$ being the Ricci and Riemann
tensors respectively. The action (27) belongs to a subclass of Horndeski
theories with the coupling functions

$\displaystyle G_{2}=\eta X+8\alpha\xi^{(4)}(\phi)X^{2}(3-\ln|X|)\,,\qquad
G_{3}=4\alpha\xi^{(3)}(\phi)X(7-3\ln|X|)\,,$ $\displaystyle G_{4}=\frac{M_{\rm
Pl}^{2}}{2}+4\alpha\xi^{(2)}(\phi)X(2-\ln|X|)\,,\qquad
G_{5}=-4\alpha\xi^{(1)}(\phi)\ln|X|\,,$ (29)

where $\xi^{(n)}(\phi)\equiv{\rm d}^{n}\xi(\phi)/{\rm d}\phi^{n}$.

The background Eqs. (9), (10), and (12) reduce, respectively, to

$\displaystyle h^{\prime}=\frac{2(1-h)[M_{\rm Pl}^{2}-8\alpha
h(\phi^{\prime\prime}\xi_{,\phi}+\phi^{\prime
2}\xi_{,\phi\phi})]-r^{2}(2\rho+\eta h\phi^{\prime 2})}{2M_{\rm
Pl}^{2}r+8\alpha(1-3h)\phi^{\prime}\xi_{,\phi}}\,,$ (30) $\displaystyle
f^{\prime}=\frac{f[2(1-h)M_{\rm Pl}^{2}+r^{2}(2P+\eta h\phi^{\prime
2})]}{h[2M_{\rm Pl}^{2}r+8\alpha(1-3h)\phi^{\prime}\xi_{,\phi}]}\,,$ (31)
$\displaystyle\eta
h\phi^{\prime\prime}+\frac{\eta(f^{\prime}hr+fh^{\prime}r+4fh)}{2rf}\phi^{\prime}-\frac{2\alpha\xi_{,\phi}[h(1-h)(2ff^{\prime\prime}-f^{\prime
2})+ff^{\prime}h^{\prime}(1-3h)]}{r^{2}f^{2}}=0\,.$ (32)

The linear stability conditions (15)-(17) in the odd-parity sector translate
to

$\displaystyle{\cal G}$ $\displaystyle=$ $\displaystyle M_{\rm
Pl}^{2}-\frac{4\alpha\xi_{,\phi}\phi^{\prime}f^{\prime}h}{f}>0\,,$ (33)
$\displaystyle{\cal H}$ $\displaystyle=$ $\displaystyle M_{\rm
Pl}^{2}-\frac{8\alpha\xi_{,\phi}\phi^{\prime}h}{r}>0\,,$ (34)
$\displaystyle{\cal F}$ $\displaystyle=$ $\displaystyle M_{\rm
Pl}^{2}-4\alpha\left(2h\xi_{,\phi\phi}\phi^{\prime
2}+2h\xi_{,\phi}\phi^{\prime\prime}+h^{\prime}\xi_{,\phi}\phi^{\prime}\right)>0\,.$
(35)

In the limit of a small GB coupling $\alpha\to 0$, all of ${\cal G}$, ${\cal
H}$, and ${\cal F}$ approach $M_{\rm Pl}^{2}$, so the stability against odd-
parity perturbations is ensured. As we will see later in Sec. III.1, the
leading-order term of $\phi^{\prime}$ around $r=0$ is proportional to $r$, so
${\cal H}$ does not diverge at the center of star.

On using Eqs. (30) and (31) to eliminate $\rho$ and $P$, the no-ghost
condition (19) in the even-parity sector yields

$\displaystyle{\cal K}$ $\displaystyle=$ $\displaystyle 2\eta h\phi^{\prime
2}r^{2}(M_{\rm Pl}^{2}r-8\alpha h\xi_{,\phi}\phi^{\prime})^{2}$ (36)
$\displaystyle+64\alpha^{2}h(h-1)\phi^{\prime 2}\xi_{,\phi}^{2}[M_{\rm
Pl}^{2}(2rh^{\prime}-h+1)-4\alpha\xi_{,\phi}\phi^{\prime}h^{\prime}(1+3h)+8\alpha
h(h-1)(\xi_{,\phi}\phi^{\prime\prime}+\xi_{,\phi\phi}\phi^{\prime 2})]>0\,.$

In the limit $\alpha\to 0$, we have ${\cal K}\to 2\eta h\phi^{\prime 2}M_{\rm
Pl}^{4}r^{4}$ and hence the ghost can be avoided for

$\eta>0\,.$ (37)

In the following, we will focus on the case in which $\eta$ is positive.

The condition (22) for the absence of the Laplacian instability along the
radial direction reduces to

$c_{r3}^{2}=\frac{2h\phi^{\prime 2}[\eta fr^{2}(M_{\rm Pl}^{2}r-8\alpha
h\xi_{,\phi}\phi^{\prime})^{2}-32\alpha^{2}(h-1)\xi_{,\phi}^{2}\\{M_{\rm
Pl}^{2}[f(h-1)-2rf^{\prime}h]+4\alpha\xi_{,\phi}\phi^{\prime}f^{\prime}h(1+3h)\\}]}{f{\cal
K}}>0\,.$ (38)

Expanding $c_{r3}^{2}$ around $\alpha=0$, we obtain

$c_{r3}^{2}=1-\frac{64\xi_{,\phi}^{2}(1-h)(f^{\prime}h-fh^{\prime})}{\eta
fM_{\rm Pl}^{2}r^{3}}\alpha^{2}+{\cal O}(\alpha^{3})\,,$ (39)

and hence $c_{r3}^{2}\to 1$ as $\alpha\to 0$.

The squared angular propagation speeds of even-parity perturbations are
complicated, but expanding $c_{\Omega,\pm}^{2}$ around $\alpha=0$ leads to

$\displaystyle
c_{\Omega,\pm}^{2}=1+\frac{\xi_{,\phi}[2\sqrt{\eta}fh(2f-rf^{\prime})r\phi^{\prime}\pm\sqrt{2{\cal
B}}]}{\sqrt{\eta}f^{2}M_{\rm Pl}^{2}r^{2}}|\alpha|+{\cal O}(\alpha^{2})\,,$
(40)

where

${\cal B}\equiv M_{\rm Pl}^{2}[f^{\prime
2}hr^{2}+2f^{2}(rh^{\prime}-2h+2)-rf(2f^{\prime\prime}hr+rf^{\prime}h^{\prime}-2f^{\prime}h)]^{2}+2\eta
f^{2}h^{2}(rf^{\prime}-2f)^{2}r^{2}\phi^{\prime 2}\,.$ (41)

Under the no-ghost condition $\eta>0$, ${\cal B}$ is positive and hence
$c_{\Omega,\pm}^{2}$ are the real values. In the limit $\alpha\to 0$,
$c_{\Omega,\pm}^{2}$ approach 1.

From the above discussion, all the linear stability conditions should be
consistently satisfied for $\eta>0$ and $|\alpha|\ll 1$. We note, however,
that the quantities like ${\cal H}$, $c_{r3}^{2}$, $c_{\Omega,\pm}^{2}$
contain positive power-law terms of $r$ in the denominators. To show the
finiteness of these quantities at $r=0$, we derive the solutions to $f$, $h$,
$\phi$ expanded around $r=0$ in Sec. III.1. For the existence of hairy stars
free from instabilities, we then put a limit on the coupling constant
$\alpha$.

###  III.1 Solutions expanded around $r=0$ and their stability

Around the center of star, we impose the regular boundary conditions
$f(0)=f_{c}$, $h(0)=1$, $\phi(0)=\phi_{c}$, $\rho(0)=\rho_{c}$, $P(0)=P_{c}$
and
$f^{\prime}(0)=h^{\prime}(0)=\phi^{\prime}(0)=\rho^{\prime}(0)=P^{\prime}(0)=0$.
Then, around $r=0$, the scalar field is expanded as

$\phi=\phi_{c}+\phi_{2}r^{2}+{\cal O}(r^{3})\,,$ (42)

likewise for $f$, $h$, $\rho$, and $P$. We can also expand the coupling
function $\xi(\phi)$ (and its $\phi$ derivatives), as

$\xi(\phi)=\xi(\phi_{c})+\sum_{n\geq
1}\xi^{(n)}(\phi_{c})\frac{(\phi-\phi_{c})^{n}}{n!}\,.$ (43)

On using the background Eqs. (30)-(32) with Eq. (8), the quantity $\phi_{2}$
in Eq. (42) obeys the following algebraic equation

$Y(x)\equiv\frac{9x(1-16x)^{3}}{2[1+3w_{c}(1-16x)]}=\kappa\alpha^{2}\,,$ (44)

where

$x\equiv\frac{\alpha\xi_{,\phi}(\phi_{c})}{M_{\rm
Pl}^{2}}\phi_{2}\,,\qquad\kappa\equiv\frac{\rho_{c}^{2}\xi_{,\phi}^{2}(\phi_{c})}{\eta
M_{\rm Pl}^{6}}\,,\qquad w_{c}=\frac{P_{c}}{\rho_{c}}\,.$ (45)

We assume that the EOS parameter $w_{c}$ is in the range $w_{c}>0$. In the
small-coupling limit $|\alpha|\ll 1$, we have
$9x/[2(1+3w_{c})]\simeq\kappa\alpha^{2}$ and hence

$\phi_{2}=\frac{2\xi_{,\phi}(\phi_{c})\rho_{c}^{2}(1+3w_{c})}{9\eta M_{\rm
Pl}^{4}}\alpha+{\cal O}(\alpha^{3})\,.$ (46)

Substituting the leading-order solution of Eq. (46) into the definition of $x$
in Eq. (45), we have

$x>0\,,$ (47)

under the condition $\eta>0$. By the end of this section, we will not exploit
the small $\alpha$ expansion to discuss the linear stability of NSs.

The solutions to $f$, $h$, $P$ expanded around $r=0$ are given, respectively,
by

$\displaystyle f$ $\displaystyle=$ $\displaystyle f_{c}+\frac{f_{c}[M_{\rm
Pl}^{2}\rho_{c}(1+3w_{c})-48\alpha\xi_{,\phi}(\phi_{c})\rho_{c}w_{c}\phi_{2}]}{6[M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]^{2}}r^{2}+{\cal O}(r^{3})\,,$
(48) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{\rho_{c}}{3[M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]}r^{2}+{\cal O}(r^{3})\,,$ (49)
$\displaystyle P$ $\displaystyle=$ $\displaystyle
P_{c}-\frac{\rho_{c}(1+w_{c})[M_{\rm
Pl}^{2}\rho_{c}(1+3w_{c})-48\alpha\xi_{,\phi}(\phi_{c})\rho_{c}w_{c}\phi_{2}]}{12[M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]^{2}}r^{2}+{\cal O}(r^{3})\,.$
(50)

The scalar field is of the form (42) with $\phi_{2}$ satisfying the relation
(44). On using these solutions, the quantities ${\cal G}$, ${\cal H}$, and
${\cal F}$ at $r=0$ reduce to

${\cal G}(r=0)=M_{\rm Pl}^{2}\,,\qquad{\cal H}(r=0)={\cal F}(r=0)=M_{\rm
Pl}^{2}\left(1-16x\right)\,,$ (51)

where ${\cal H}$ is finite at $r=0$ due to the property that
$\phi^{\prime}=2\phi_{2}r$ at leading order. The squared propagation speeds of
odd-parity perturbations along the radial and angular directions are given by

$c_{r}^{2}(r=0)=c_{\Omega}^{2}(r=0)=\frac{1}{1-16x}\,.$ (52)

Then, there are neither ghost nor Laplacian instabilities in the odd-parity
sector if

$0<x<\frac{1}{16}\,,$ (53)

where we have also taken into account the condition (47). Both
$c_{r}^{2}(r=0)$ and $c_{\Omega}^{2}(r=0)$ are larger than 1. In the limit
that $\alpha\to 0$, we have $x\to 0$ and hence
$c_{r}^{2}(r=0)=c_{\Omega}^{2}(r=0)\to 1$ as expected. In the subject of NSs,
“superluminality” is sometimes argued as a sign of “acausality.” We emphasize
that the propagation speed of scalar-field perturbations just fixes a causal
boundary of the scalar field at each position, and the superluminal speeds do
not mean acausality.

In the even-parity sector, the leading-order term of ${\cal K}$ expanded
around $r=0$ is proportional to $r^{6}$, i.e.,

${\cal K}=\frac{16\rho_{c}^{2}M_{\rm
Pl}^{2}x[1+48x+3w_{c}(1-16x)]}{9(1-16x)}r^{6}+{\cal O}(r^{7})\,,$ (54)

where we used Eq. (44) to eliminate $\kappa\alpha^{2}$. Under the condition
(53), the coefficient of $r^{6}$ in Eq. (54) is positive. The squared radial
and angular propagation speeds of $\delta\phi$ at $r=0$ are given by

$c_{r3}^{2}(r=0)=c_{\Omega-}^{2}(r=0)=\frac{(1+3w_{c})(1-64x)+2304w_{c}x^{2}}{(1-16x)[1+48x+3w_{c}(1-16x)]}\,,$
(55)

whereas $c_{\Omega+}^{2}(r=0)$ is equivalent to Eq. (52). Under the condition
(53), the absence of Laplacian instability of $\delta\phi$ requires a
positivity of the numerator of Eq. (55), so that

$0<x<x_{m}\equiv\frac{2(1+3w_{c})-\sqrt{(1+3w_{c})(4+3w_{c})}}{144w_{c}}\,.$
(56)

For $0<w_{c}<\infty$ we have $1/64<x_{m}<1/48$, so $x_{m}$ is smaller than
$1/16$. The function $Y(x)$ in Eq. (44) has a maximum value $Y(x_{m})$ at
$x=x_{m}$. Provided that $0<\kappa\alpha^{2}<Y(x_{m})$, there are solutions to
Eq. (44). This gives an upper bound on $|\alpha|$, as

$|\alpha|<\frac{\sqrt{\eta}M_{\rm
Pl}^{3}}{\rho_{c}\xi_{,\phi}(\phi_{c})}\frac{\sqrt{6}(2\sqrt{1+3w_{c}}-\sqrt{4+3w_{c}})^{1/2}[\sqrt{(1+3w_{c})(4+3w_{c})}-2+3w_{c}]^{3/2}}{216w_{c}^{2}(\sqrt{1+3w_{c}}+\sqrt{4+3w_{c}})^{1/2}}\,.$
(57)

Among the two solutions of $\kappa\alpha^{2}=Y(x)$, one of them ($x=x_{1}$)
exists in $0<x_{1}<x_{m}$, while the other ($x=x_{2}$) is in the region
$x_{m}<x_{2}<1/16$. The former is in the region (56) satisfying the condition
of Laplacian stability. In summary, as long as $\alpha$ is in the range (57),
there is a solution $x=x_{1}$ consistent will all the linear stability
conditions at $r=0$. The criterion (57) is also valid for nonrelativistic
stars with $w_{c}\ll 1$. In this case, the scalar-GB coupling is constrained
to be

$|\alpha|<\frac{9\sqrt{6\eta}M_{\rm
Pl}^{3}}{128\rho_{c}\xi_{,\phi}(\phi_{c})}\,,\quad{\rm for}\quad w_{c}\to
0\,,$ (58)

which can be applied to nonrelativistic objects such as Sun and Earth.

Figure 1: Maximum value of $|\alpha|$ versus $w_{c}$ constrained by the bound
(61) for the linear GB coupling ($n=1$) with $\eta=1$. The solid and dashed
lines correspond to the upper limits $|\alpha_{\rm max}|$ for SLy and FPS
EOSs, respectively, with the central density in the range
$10^{-2}\rho_{0}\leq\rho_{c}\leq 20\rho_{0}$.

Let us consider the power-law scalar-GB coupling

$\xi(\phi)=M_{\rm Pl}^{2-n}r_{0}^{2}\phi^{n}\,,$ (59)

where

$r_{0}=\sqrt{\frac{8\pi M_{\rm Pl}^{2}}{\rho_{0}}}=89.664~{}{\rm
km}\,,\qquad\rho_{0}=m_{n}n_{0}=1.6749\times 10^{14}~{}{\rm g}\cdot{\rm
cm}^{-3}\,.$ (60)

Here, $m_{n}=1.6749\times 10^{-24}$ g is the neutron mass and
$n_{0}=0.1~{}(\rm fm)^{-3}$ is the typical density of NSs. For this coupling,
the bound (57) yields

$|\alpha|<|\alpha_{\rm max}|\equiv\frac{\sqrt{\eta}}{8\pi
n}\frac{\rho_{0}}{\rho_{c}}\left(\frac{\phi_{c}}{M_{\rm
Pl}}\right)^{1-n}\frac{\sqrt{6}(2\sqrt{1+3w_{c}}-\sqrt{4+3w_{c}})^{1/2}[\sqrt{(1+3w_{c})(4+3w_{c})}-2+3w_{c}]^{3/2}}{216w_{c}^{2}(\sqrt{1+3w_{c}}+\sqrt{4+3w_{c}})^{1/2}}\,.$
(61)

For the linear GB coupling ($n=1$), the central density $\rho_{c}$ and EOS
parameter $w_{c}$ determine the upper limit of $|\alpha|$, without having the
dependence of $\phi_{c}$. In Fig. 1, we plot $|\alpha_{\rm max}|$ as a
function of $w_{c}$ for SLy (solid) and FPS (dashed) EOSs Haensel and Potekhin
(2004) in the range $10^{-2}\rho_{0}\leq\rho_{c}\leq 20\rho_{0}$. As
$\rho_{c}$ increases, $w_{c}$ grows from nonrelativistic values of order
$10^{-3}$ to relativistic values of order $10^{-1}$. We find that
$|\alpha_{\rm max}|$ is a decreasing function of $w_{c}$. For SLy EOS, the
maximum mass of NS is reached around the central density $\rho_{c}\approx
15\rho_{0}$ with the EOS parameter $w_{c}\approx 0.5$. To realize such a
maximum mass of NS, the scalar-GB coupling is constrained to be

$|\alpha|<3\times 10^{-4}\,.$ (62)

For FPS EOS, we also obtain a similar upper bound. The existence of NS
solutions with the central density $\rho_{c}\gtrsim 15\rho_{0}$ gives the
value of $|\alpha_{\rm max}|$ even smaller than $3\times 10^{-4}$.

The dilatonic GB coupling $\alpha_{\rm GB}$ in the regime $\phi/M_{\rm Pl}\ll
1$, which was discussed in Refs. Saffer and Yagi (2021); Lyu  _et al._ (2022)
with the unit $M_{\rm pl}=1/\sqrt{8\pi}$, is related to our linear scalar-GB
coupling $\alpha$ by $\alpha_{\rm GB}=\alpha M_{\rm Pl}r_{0}^{2}=\alpha
r_{0}^{2}/\sqrt{8\pi}$. Then, the theoretical bound (62) translates to

$\sqrt{|\alpha_{\rm GB}|}<0.7~{}{\rm km}\,.$ (63)

This is tighter than the typical observational bounds $\sqrt{\alpha_{\rm
GB}}<{\cal O}(1)~{}{\rm km}$ derived from the GW measurements of NS-NS, BH-NS,
and BH-BH binaries, see Table I of Ref. Lyu  _et al._ (2022). Thus, the
existence of NS solutions consistent with the linear stability conditions
gives a new theoretical upper bound on the dilatonic GB coupling.

###  III.2 Solutions expanded at spatial infinity and their stability

The surface of star is defined by the radius $r_{s}$ at which the fluid
pressure $P$ vanishes. Outside the star ($r>r_{s}$), we have $\rho=0=P$ in the
background Eqs. (9)-(12). Imposing the asymptotic flatness at spatial
infinity, we can expand $f$, $h$, and $\phi$ in the forms
$f=1+\sum_{i=1}\hat{f}_{i}/r^{i}$, $h=1+\sum_{i=1}\hat{h}_{i}/r^{i}$, and
$\phi=\phi_{0}+\sum_{i=1}\hat{\phi}_{i}/r^{i}$. We also use the expansions of
$\xi(\phi)$ and its $\phi$ derivatives analogous to Eq. (43), with the
replacement $\phi_{c}\to\phi_{0}$. The large-distance solutions consistent
with the background Eqs. (30)-(32) are given by

$\displaystyle f$ $\displaystyle=$ $\displaystyle 1-\frac{2M}{r}+\frac{\eta
M\hat{\phi}_{1}^{2}}{6M_{\rm Pl}^{2}r^{3}}+\frac{M\hat{\phi}_{1}[\eta
M\hat{\phi}_{1}+24\alpha\xi_{,\phi}(\phi_{0})]}{3M_{\rm Pl}^{2}r^{4}}+{\cal
O}(r^{-5})\,,$ (64) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{2M}{r}+\frac{\eta\hat{\phi}_{1}^{2}}{2M_{\rm Pl}^{2}r^{2}}+\frac{\eta
M\hat{\phi}_{1}^{2}}{2M_{\rm Pl}^{2}r^{3}}+\frac{2M\hat{\phi}_{1}[\eta
M\hat{\phi}_{1}+24\alpha\xi_{,\phi}(\phi_{0})]}{3M_{\rm Pl}^{2}r^{4}}+{\cal
O}(r^{-5})\,,$ (65) $\displaystyle\phi$ $\displaystyle=$
$\displaystyle\phi_{0}+\frac{\hat{\phi}_{1}}{r}+\frac{M\hat{\phi}_{1}}{r^{2}}+\frac{(16M^{2}M_{\rm
Pl}^{2}-\eta\hat{\phi}_{1}^{2})\hat{\phi}_{1}}{12M_{\rm
Pl}^{2}r^{3}}+\frac{M[6\eta M^{2}M_{\rm
Pl}^{2}\hat{\phi}_{1}-\eta^{2}\hat{\phi}_{1}^{3}-12\alpha MM_{\rm
Pl}^{2}\xi_{,\phi}(\phi_{0})]}{3\eta M_{\rm Pl}^{2}r^{4}}+{\cal O}(r^{-5})\,,$
(66)

where we set $\hat{f}_{1}=-2M$. Then the quantities (15)-(17) can be estimated
as ${\cal G}=M_{\rm Pl}^{2}+8\alpha
M\xi_{,\phi}(\phi_{0})\hat{\phi}_{1}r^{-4}+{\cal O}(r^{-5})$, ${\cal H}=M_{\rm
Pl}^{2}+8\alpha\xi_{,\phi}(\phi_{0})\hat{\phi}_{1}r^{-3}+{\cal O}(r^{-4})$,
and ${\cal F}=M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{0})\hat{\phi}_{1}r^{-3}+{\cal O}(r^{-4})$,
so the squared radial and angular propagation speeds reduce, respectively, to

$c_{r}^{2}=1+\frac{16\alpha\xi_{,\phi}(\phi_{0})\hat{\phi}_{1}}{M_{\rm
Pl}^{2}r^{3}}+{\cal O}(r^{-4})\,,\qquad
c_{\Omega}^{2}=1-\frac{8\alpha\xi_{,\phi}(\phi_{0})\hat{\phi}_{1}}{M_{\rm
Pl}^{2}r^{3}}+{\cal O}(r^{-4})\,.$ (67)

As $r\to\infty$, both $c_{r}^{2}$ and $c_{\Omega}^{2}$ approach 1.

In the even-parity sector, the quantity ${\cal K}$ is expressed as

${\cal K}=2\eta M_{\rm Pl}^{4}\hat{\phi}_{1}^{2}+\frac{4\eta M_{\rm
Pl}^{4}M\hat{\phi}_{1}^{2}}{r}+{\cal O}(r^{-2})\,,$ (68)

whose positivity is ensured for $\eta>0$. The radial propagation speed squared
(38) yields

$c_{r3}^{2}=1-\frac{128\alpha^{2}M\xi_{,\phi}^{2}(\phi_{0})\hat{\phi}_{1}^{2}}{M_{\rm
Pl}^{4}r^{7}}+{\cal O}(r^{-8})\,,$ (69)

which quickly approaches 1 at large distances even compared to $c_{r}^{2}$ and
$c_{\Omega}^{2}$. For the angular propagation, we obtain

$c_{\Omega\pm}^{2}=1\pm\frac{4\xi_{,\phi}(\phi_{0})[(\eta\hat{\phi}_{1}^{2}+72M^{2}M_{\rm
Pl}^{2})^{1/2}\mp\sqrt{\eta}\hat{\phi}_{1}]}{\sqrt{\eta}M_{\rm
Pl}^{2}r^{3}}|\alpha|+{\cal O}(r^{-4})\,.$ (70)

For $r\gg r_{s}$ all the squared propagation speeds given above rapidly
approach 1, so the sign of $\alpha$ does not matter for the discussion of
Laplacian instability. Provided that $\eta>0$, there are neither ghost nor
Laplacian instabilities at spatial infinity.

###  III.3 Numerical solutions and stability conditions

The discussions in Secs. III.1 and III.2 show that, under the bound (57) with
$\eta>0$, there are NS solutions with a nontrivial scalar profile consistent
with all the linear stability conditions around $r=0$ and $r\to\infty$.
However, they do not necessarily guarantee the linear stability of solutions
at intermediate distances, so we will numerically study whether neither ghost
nor Laplacian instabilities appear at any radius $r$.

For concreteness, we study the linear scalar-GB coupling given by the $n=1$
case of Eq. (59). Since this corresponds to $G_{2}=\eta X$, $G_{3}=0$,
$G_{4}=M_{\rm Pl}^{2}/2$, and $G_{5}=-4\alpha M_{\rm Pl}r_{0}^{2}\ln|X|$, it
belongs to a subclass of shift-symmetric Horndeski theories. The same scalar-
GB coupling can also accommodate the dilatonic coupling $\xi(\phi)\propto{\rm
e}^{\mu\phi}$ in the limit $|\mu\phi|\ll 1$. Since ${\cal P}_{\phi}=0$ in the
scalar-field Eq. (12), we have $r^{2}\sqrt{f/h}\,J^{r}=Q={\rm constant}$ and
hence

$\frac{h}{f}\left[\eta fr^{2}\phi^{\prime}+4\alpha M_{\rm
Pl}r_{0}^{2}f^{\prime}(h-1)\right]=Q\sqrt{\frac{f}{h}}\,.$ (71)

To satisfy the boundary conditions of $f$, $h$, and $\phi^{\prime}$ at $r=0$,
we require that $Q=0$. Then, the field derivative can be expressed as

$\phi^{\prime}(r)=-\frac{4\alpha M_{\rm Pl}r_{0}^{2}f^{\prime}(h-1)}{\eta
fr^{2}}\,.$ (72)

Substituting Eqs. (48) and (49) into Eq. (72) around $r=0$, we have
$\phi^{\prime}(r)\propto r$ as consistent with Eq. (42). At spatial infinity,
using the expanded solutions (64) and (65) in Eq. (72) leads to
$\phi^{\prime}(r)\propto r^{-5}$. Since the integration constant $Q$
corresponds to the scalar charge, the choice of $Q=0$ means that
$\hat{\phi}_{1}=0$ in the expansion of Eq. (66) and hence
$\phi=\phi_{0}-4\alpha M^{2}\xi_{,\phi}(\phi_{0})/(\eta r^{4})$. In Eq. (68)
the leading-order term of ${\cal K}$ vanishes, but it is replaced by ${\cal
K}=512\alpha^{2}M^{4}M_{\rm Pl}^{4}\xi_{,\phi}(\phi_{0})^{2}/(\eta
r^{6})+{\cal O}(r^{-7})$. Then, the no-ghost condition is satisfied for
$\eta>0$.

To perform the numerical integration, we introduce the following variables

$s=\ln\frac{r}{r_{0}}\,,\qquad{\cal M}(r)=4\pi rM_{\rm Pl}^{2}(1-h)\,,\qquad
m(r)=\frac{3{\cal M}(r)}{4\pi r_{0}^{3}\rho_{0}}\,,\qquad
y=\frac{\rho}{\rho_{0}}\,,\qquad z=\frac{P}{\rho_{0}}\,.$ (73)

The ADM mass of star can be computed as

$M\equiv{\cal M}(r\to\infty)=2.5435\times 10^{2}\,m_{\rm\infty}\,M_{\odot}\,,$
(74)

where $m_{\infty}\equiv m(r\to\infty)$ and $M_{\odot}=1.9884\times 10^{33}$ g
is the solar mass.

For the NS EOS, we exploit the analytic representation of SLy EOS parametrized
by

$\xi=\log_{10}(\rho/{\rm g\cdot cm}^{-3})\,,\qquad\zeta=\log_{10}(P/{\rm
dyn\cdot cm}^{-2})\,,$ (75)

where the relation between $\xi$ and $\zeta$ is given in Ref. Haensel and
Potekhin (2004). In terms of $y$ and $z$, we can express $\xi$ and $\zeta$ as

$\xi=\alpha_{1}+\alpha_{2}\ln y\,,\qquad\zeta=\alpha_{3}+\alpha_{2}\ln z\,,$
(76)

where $\alpha_{1}=\ln(\rho_{0}/{\rm g\cdot cm^{-3}})/\ln 10$, $\alpha_{2}=(\ln
10)^{-1}$, and $\alpha_{3}=\ln(\rho_{0}/{\rm dyn\cdot cm^{-2}})/\ln 10$. Then,
the EOS translates to

$z=\exp\left[\frac{\zeta(\xi)-\alpha_{3}}{\alpha_{2}}\right]\,,$ (77)

so that $z$ is known by the analytic representation of $\zeta(\xi)$. The
continuity Eq. (8) gives the differential equation for $y$, as

$\frac{{\rm d}y}{{\rm d}s}=-\frac{y(y+z)}{2z}\left(\frac{{\rm d}\zeta}{{\rm
d}\xi}\right)^{-1}\frac{1}{f}\frac{{\rm d}f}{{\rm d}s}\,.$ (78)

We derive the differential equations of $f$ and $h$ by solving Eqs. (9), (11),
and (12) for $f^{\prime\prime}$, $h^{\prime}$, and $\phi^{\prime\prime}$. We
replace the first-order field derivative $\phi^{\prime}$ in the differential
equations of $f$ and $h$ by using Eq. (72).

Figure 2:  (Left) Mass function ${\cal M}$ (normalized by the solar mass
$M_{\odot}$), field derivative $\phi^{\prime}$ (normalized by $M_{\rm
Pl}/r_{0}$), fluid pressure $P$ (normalized by $\rho_{c}$) versus $r/r_{0}$
for the linear scalar-GB coupling $\xi(\phi)=M_{\rm Pl}r_{0}^{2}\phi$ with
$\alpha=2\times 10^{-4}$ and $\eta=1$. For the perfect fluid, we choose SLy
EOS with the central density $\rho_{c}=15\rho_{0}$. (Right) $c_{r}^{2}$,
$c_{\Omega}^{2}$, $c_{r3}^{2}$, $c_{\Omega+}^{2}$, and $c_{\Omega-}^{2}$
versus $r/r_{0}$ for the same model parameters as those used in the left.

In the left panel of Fig. 2, we plot ${\cal M}$, $\phi^{\prime}$, and $P$ as a
function of $r/r_{0}$ for $\alpha=2\times 10^{-4}$, $\eta=1$, and
$\rho_{c}=15\rho_{0}$. We choose the boundary conditions at $r/r_{0}=10^{-5}$
to be consistent with Eqs. (42) and (48)-(50). As we estimated in Sec. III.1,
the field derivative increases as $\phi^{\prime}(r)\propto r$ around $r=0$,
with ${\cal M}(r)\propto r^{3}$ and $P(r)$ decreasing according to Eq. (50).
The radius at which $P$ vanishes is $r_{s}=0.114r_{0}\simeq 10.22$ km, which
corresponds to the surface of star. The field derivative starts to decrease
around $r=r_{s}$ and it has the dependence $\phi^{\prime}(r)\propto r^{-5}$ at
large distances. The growth of mass function ${\cal M}(r)$ saturates around
the surface of star and it quickly approaches the ADM mass $M$. In the
numerical simulation of Fig. 2 we have $M=2.003M_{\odot}$, which is smaller
than the corresponding mass $2.044M_{\odot}$ in GR ($\alpha=0$). This
suppressed ADM mass in comparison to GR is consistent with the result obtained
for the dilatonic scalar-GB coupling $\xi(\phi)\propto{\rm e}^{\mu\phi}$ in
Ref. Pani  _et al._ (2011).

In the right panel of Fig. 2, we show $c_{r}^{2}$, $c_{\Omega}^{2}$,
$c_{r3}^{2}$, $c_{\Omega+}^{2}$, and $c_{\Omega-}^{2}$ versus $r/r_{0}$ for
the same model parameters as those used in the left. As estimated from Eq.
(52), the radial and angular propagation speeds in the odd-parity sector are
superluminal in the regime $r\ll r_{s}$. In the even-parity sector the
evolution of $c_{\Omega+}^{2}$ around $r=0$ is similar to $c_{r}^{2}$ and
$c_{\Omega}^{2}$, whereas $c_{r3}^{2}$ and $c_{\Omega-}^{2}$ are in the
subluminal range as estimated by Eq. (39). Around $r=r_{s}$, $c_{\Omega+}^{2}$
and $c_{\Omega-}^{2}$ exhibit temporal increase and decrease, respectively,
but they remain to be finite positive values without Laplacian instabilities.
Outside the star, all the squared propagation speeds shown in Fig. 2 quickly
approach 1, as consistent with the discussion in Sec. III.2. Numerically, we
have also confirmed that the no-ghost conditions of odd- and even-parity
perturbations are satisfied at any distance $r$.

The numerical simulation of Fig. 2 corresponds to $\rho_{c}=15\rho_{0}$,
$w_{c}=0.48$, $\eta=1$, and $n=1$, so the bound (61) gives
$|\alpha|<3.19\times 10^{-4}$. For the coupling $\alpha$ in this range, we
numerically confirmed that all the linear stability conditions are satisfied
at any distance $r$. Since $c_{r3}^{2}$ and $c_{\Omega-}^{2}$ are smallest at
the center of NS, the upper limit of $|\alpha|$ is determined by their values
at $r=0$ as we performed in Sec. III.1. When $\alpha<0$, the background
solution and its linear stability are similar to those for the corresponding
positive value $|\alpha|$. If $|\alpha|$ exceeds the upper limit $|\alpha_{\rm
max}|$, there are Laplacian instabilities associated with negative values of
$c_{r3}^{2}$ and $c_{\Omega-}^{2}$ at $r=0$. Thus, the theoretical bound (61)
is sufficiently accurate for the estimation of maximum allowed values of
$|\alpha|$. We have also performed numerical simulations for the scalar-GB
couplings (59) with different powers $n$ and confirmed that, for $\alpha$ in
the range (61) with $\eta>0$, there are hairy NS solutions consistent with all
the linear stability conditions.

##  IV Scalar-Gauss-Bonnet and regular couplings

In this section, we study the existence and the linear stability of NS
solutions with a nontrivial scalar profile in the presence of several regular
couplings besides the scalar-GB coupling $\alpha\xi(\phi)R_{\rm GB}^{2}$. We
also take into account the Einstein-Hilbert term $M_{\rm Pl}^{2}R/2$ and the
canonical kinetic term $\eta X$ (with $\eta>0$) in the action. We exploit SLy
EOS for the numerical analysis in this section.

###  IV.1 Cubic Galileon and scalar-GB couplings

The cubic Galileon corresponds to the coupling function
$G_{3}\supset\mu_{3}X$, where $\mu_{3}$ is a constant. Let us consider
theories given by the action

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm Pl}^{2}}{2}R+\eta
X+\alpha\xi(\phi)R_{\rm GB}^{2}+\mu_{3}X\square\phi\right]\,.$ (79)

We first derive the solutions expanded around $r=0$ and then discuss the
linear stability of them. In doing so, we use the approximation $|\alpha|\ll
1$ and pick up terms up to the order of $\alpha^{2}$. Then, the coefficient
$\phi_{2}$ in the field expansion (42) is

$\phi_{2}=\frac{3\eta M_{\rm
Pl}^{6}-32\rho_{c}^{2}(1+2w_{c})\xi_{,\phi}^{2}(\phi_{c})\alpha^{2}}{24\mu_{3}M_{\rm
Pl}^{6}}\left[1-\sqrt{1-\frac{32M_{\rm
Pl}^{8}\rho_{c}^{2}(1+3w_{c})\xi_{,\phi}(\phi_{c})\mu_{3}\alpha}{[3\eta M_{\rm
Pl}^{6}-32\rho_{c}^{2}(1+2w_{c})\xi_{,\phi}^{2}(\phi_{c})\alpha^{2}]^{2}}}\right]+{\cal
O}(\alpha^{3})\,,$ (80)

where we have chosen the branch recovering Eq. (46) in the limit $\mu_{3}\to
0$. We do not necessarily assume that the cubic Galileon coupling is of the
same order as the scalar-GB coupling. The consistency of the small $\alpha$
expansion requires that

$\xi_{,\phi}(\phi_{c})|\mu_{3}\alpha|\ll\frac{\eta^{2}M_{\rm
Pl}^{4}}{\rho_{c}^{2}(1+3w_{c})}\,,$ (81)

where we also assumed $\alpha^{2}\rho_{c}^{2}\xi_{,\phi}^{2}(\phi_{c})\ll\eta
M_{\rm Pl}^{6}$. Up to the order of $r^{2}$, the metric components and fluid
pressure are the same forms as Eqs. (48)-(50) with
$\phi^{\prime}(r)=2\phi_{2}r$. In the limit that $\alpha\to 0$, we have
$\phi_{2}=0$ even for $\mu_{3}\neq 0$. This shows that the cubic-order
coupling alone does not give rise to NS solutions with $\phi^{\prime}(r)\neq
0$. This is consistent with the no-hair argument of Ref. Lehébel  _et al._
(2017) for regular couplings in shift-symmetric Horndeski theories. In other
words, the scalar-GB coupling is needed for the realization of NSs with a
nontrivial scalar profile.

On using the above background solutions around $r=0$, we find that ${\cal G}$,
${\cal H}$, and ${\cal F}$ are of the same form as Eq. (51), where
$x=\alpha\xi_{,\phi}(\phi_{c})\phi_{2}/M_{\rm Pl}^{2}$ and $\phi_{2}$
satisfying Eq. (80). Then, the absence of ghost/Laplacian instabilities in the
odd-parity sector requires that $x<1/16$. The squared propagation speeds
$c_{r}^{2}$, $c_{\Omega}^{2}$ and $c_{\Omega+}^{2}$ at $r=0$ are equivalent to
$1/(1-16x)$. The other two squared propagation speeds at $r=0$, which are
relevant to the radial and angular propagations of $\delta\phi$, are modified
to

$c_{r3}^{2}(r=0)=c_{\Omega-}^{2}(r=0)=\frac{\eta-8\mu_{3}\phi_{2}}{\eta-12\mu_{3}\phi_{2}}-\frac{64\rho_{c}^{2}\xi_{,\phi}^{2}(\phi_{c})[\eta(1+w_{c})-2\mu_{3}\phi_{2}(5+6w_{c})]}{3M_{\rm
Pl}^{6}(\eta-12\mu_{3}\phi_{2})^{2}}\alpha^{2}+{\cal O}(\alpha^{4})\,.$ (82)

Expanding Eq. (80) with respect to $\alpha$ and using the leading-order
solution $\phi_{2}=2\rho_{c}^{2}(1+3w_{c})\xi_{,\phi}(\phi_{c})\alpha/(9\eta
M_{\rm Pl}^{4})$, the inequality (81) translates to $\mu_{3}\phi_{2}\ll
2\eta/9$. In the limit that $|\mu_{3}\alpha|\ll 1$, the first term on the
right hand-side of Eq. (82) approaches 1. On the other hand, even if
$|\alpha|\ll 1$, the product $|\mu_{3}\phi_{2}|$ can be as large as the order
$0.1\eta$ for a large coupling $|\mu_{3}|$. To avoid that the leading-order
term of Eq. (82) becomes negative, we require the condition
$\eta-12\mu_{3}\phi_{2}>0$. On using the solution (80) under the approximation
$\alpha^{2}\rho_{c}^{2}\xi_{,\phi}^{2}(\phi_{c})\ll\eta M_{\rm Pl}^{6}$, this
condition translates to

$\xi_{,\phi}(\phi_{c})\mu_{3}\alpha<\frac{\eta^{2}M_{\rm
Pl}^{4}}{4\rho_{c}^{2}(1+3w_{c})}\,.$ (83)

Around $r=0$, the leading-order contribution to (19) is proportional to
$r^{6}$, such that

${\cal K}=\frac{8M_{\rm Pl}^{2}[3M_{\rm
Pl}^{2}(1-16x)^{3}\phi_{2}^{2}(\eta-12\mu_{3}\phi_{2})+32\rho_{c}^{2}x^{2}]}{3(1-16x)}r^{6}+{\cal
O}(r^{7})\,.$ (84)

Here, we have not used the expansion with respect to $\alpha$. Provided that
$x<1/16$ and $\eta-12\mu_{3}\phi_{2}>0$, the leading-order term of ${\cal K}$
is positive. From the above discussion, the linear stability of NSs around
$r=0$ is ensured under the condition (83) with $x<1/16$.

At large distances, the solutions to $f$ and $h$ are the same forms as Eqs.
(64) and (65) up to the order of $r^{-4}$. The solution to the scalar field is
modified to

$\phi=\phi_{0}+\frac{\hat{\phi}_{1}}{r}+\frac{M\hat{\phi}_{1}}{r^{2}}+\frac{(16M^{2}M_{\rm
Pl}^{2}-\eta\hat{\phi}_{1}^{2})\hat{\phi}_{1}}{12M_{\rm
Pl}^{2}r^{3}}+\frac{3M_{\rm Pl}^{2}[4\eta
M^{3}\hat{\phi}_{1}-8M^{2}\xi_{,\phi}(\phi_{0})\alpha-\mu_{3}\hat{\phi}_{1}^{2}]-2\eta^{2}M\hat{\phi}_{1}^{3}}{6\eta
M_{\rm Pl}^{2}r^{4}}+{\cal O}(r^{-5})\,,$ (85)

in which the cubic Galileon coupling appears at the order of $r^{-4}$. On
using these large-distance solutions, it follows that $c_{r}^{2}$,
$c_{\Omega}^{2}$, and ${\cal K}$ are the same as those given in Eqs. (67) and
(68) up to the next-to-leading order. On the other hand, the other squared
propagation speeds are modified to

$\displaystyle c_{r3}^{2}$ $\displaystyle=$ $\displaystyle
1+\frac{4\hat{\phi}_{1}\mu_{3}}{\eta r^{3}}+{\cal O}(r^{-4})\,,$ (86)
$\displaystyle c_{\Omega\pm}^{2}$ $\displaystyle=$ $\displaystyle
1\pm\frac{\sqrt{16\alpha^{2}\eta\xi_{,\phi}^{2}(\phi_{0})(72M^{2}M_{\rm
Pl}^{2}+\eta\hat{\phi}_{1}^{2})+\mu_{3}M_{\rm
Pl}^{2}\hat{\phi}_{1}^{2}[\mu_{3}M_{\rm
Pl}^{2}-8\alpha\eta\xi_{,\phi}(\phi_{0})]}\mp\hat{\phi}_{1}[\mu_{3}M_{\rm
Pl}^{2}+4\alpha\eta\xi_{,\phi}(\phi_{0})]}{\eta M_{\rm Pl}^{2}r^{3}}$ (87)
$\displaystyle+{\cal O}(r^{-4})\,,$

both of which approach 1 in the limit $r\to\infty$.

To study the linear stability of NS solutions with a nontrivial scalar profile
discussed above, we perform the numerical integration for the linear scalar-GB
coupling $\xi(\phi)=M_{\rm Pl}r_{0}^{2}\phi$ in the presence of cubic
Galileons. Since this belongs to a subclass of shift-symmetric Horndeski
theories, we have $J^{r}=(Q/r^{2})\sqrt{f/h}$ from Eq. (12). The regularity at
$r=0$ leads to $Q=0$ and hence

$\phi^{\prime}(r)=\frac{\eta
fr}{\mu_{3}h(rf^{\prime}+4f)}\left[1-\sqrt{1+\frac{8\tilde{\mu}_{3}\alpha
r_{0}^{4}f^{\prime}h(rf^{\prime}+4f)(h-1)}{\eta^{2}r^{3}f^{2}}}\right]\,,$
(88)

where $\tilde{\mu}_{3}$ is a dimensionless coupling defined by

$\tilde{\mu}_{3}\equiv\mu_{3}\frac{M_{\rm Pl}}{r_{0}^{2}}\,.$ (89)

In the limit $\tilde{\mu}_{3}\to 0$, the branch (88) smoothly approaches (72).
In the following, we consider the positive scalar-GB coupling ($\alpha>0$) and
study the effect of cubic Galileons on the linear stability of NS.

From Eq. (80), the positive Galileon coupling $\tilde{\mu}_{3}>0$ leads to the
enhancement of $\phi^{\prime}(r)$ around $r=0$ in comparison to the case
$\tilde{\mu}_{3}=0$. On the other hand, Eq. (83) gives the following upper
bound

$\tilde{\mu}_{3}<\tilde{\mu}_{3}^{{\rm
max}}\equiv\frac{\eta^{2}}{256\pi^{2}\alpha}\left(\frac{\rho_{0}}{\rho_{c}}\right)^{2}\frac{1}{1+3w_{c}}\,.$
(90)

Due to this limit, the enhancement of $\phi^{\prime}(r)$ induced by the
positive $\tilde{\mu}_{3}$ is not so significant. With the model parameters
$\alpha=2.0\times 10^{-4}$, $\rho_{c}=15\rho_{0}$, $w_{c}=0.48$, and $\eta=1$,
for example, we have $\tilde{\mu}_{3}^{{\rm max}}=3.6\times 10^{-3}$ from Eq.
(90). In this case, we numerically find that the Galileon coupling needs to be
in the range $\tilde{\mu}_{3}<3.0\times 10^{-3}$ to satisfy all the linear
stability conditions at any distance $r$. Thus, the condition (90) gives a
good approximate criterion for the existence of NSs with a nontrivial scalar
profile without instabilities. For $\tilde{\mu}_{3}=2.5\times 10^{-3}$, the
ADM mass and radius of NS are $M=1.998M_{\odot}$ and $r_{s}=10.22$ km,
respectively. They are similar to the values $M=2.003M_{\odot}$ and
$r_{s}=10.22$ km derived for $\tilde{\mu}_{3}=0$. For increasing
$\tilde{\mu}_{3}$, the squared propagation speeds $c_{r3}^{2}$ and
$c_{\Omega-}^{2}$ tend to be larger than those for $\tilde{\mu}_{3}$ because
of the approach of the term $\eta-12\mu_{3}\phi_{2}$ to $+0$ in Eq. (82). We
note that $c_{r}^{2}$, $c_{\Omega}^{2}$, and $c_{\Omega}^{2}$ are not strongly
affected by a positive coupling $\tilde{\mu}_{3}$.

Figure 3:  (Left) The solid lines show the mass function ${\cal M}$
(normalized by the solar mass $M_{\odot}$) and field derivative
$\phi^{\prime}$ (normalized by $M_{\rm Pl}/r_{0}$) versus $r/r_{0}$ for the
negative cubic Galileon coupling $\tilde{\mu}_{3}=-0.1$ with $\alpha=2\times
10^{-4}$, $\eta=1$, and $\rho_{c}=15\rho_{0}$. The dotted lines correspond to
${\cal M}$ and $\phi^{\prime}$ for $\tilde{\mu}_{3}=0$, $\alpha=2\times
10^{-4}$, $\eta=1$, and $\rho_{c}=15\rho_{0}$, i.e., the case plotted in Fig.
2. (Right) $c_{r}^{2}$, $c_{\Omega}^{2}$, $c_{r3}^{2}$, $c_{\Omega+}^{2}$, and
$c_{\Omega-}^{2}$ versus $r/r_{0}$ for $\tilde{\mu}_{3}=-0.1$, $\alpha=2\times
10^{-4}$, $\eta=1$, and $\rho_{c}=15\rho_{0}$.

The negative value of $\tilde{\mu}_{3}$ is not bounded by Eq. (90). In the
left panel of Fig. 3, we plot ${\cal M}$ and $\phi^{\prime}$ versus $r/r_{0}$
for $\tilde{\mu}_{3}=-0.1$ as solid lines. In comparison to the case
$\tilde{\mu}_{3}=0$, the field derivative is suppressed by the negative
Galileon coupling. For $\tilde{\mu}_{3}=-0.1$ the ADM mass is
$M=2.026M_{\odot}$, which is larger than $M=2.003M_{\odot}$ derived for
$\tilde{\mu}_{3}=0$. This increase of $M$ is attributed to the suppression of
$\phi^{\prime}$ induced by the negative $\tilde{\mu}_{3}$. In the right panel
of Fig. 3, we show the five squared propagation speeds versus $r/r_{0}$ for
$\tilde{\mu}_{3}=-0.1$. They are different from 1 inside the star, but all of
them are larger than 0. In comparison to the case $\tilde{\mu}_{3}=0$ shown in
Fig. 2, there are more rapid temporal increase of $c_{r3}^{2}$ and decrease of
$c_{\Omega-}^{2}$ around the surface of star. For decreasing
$\tilde{\mu}_{3}$, we find that these transient variations of $c_{r3}^{2}$ and
$c_{\Omega-}^{2}$ tend to be more significant. Numerically, we find that there
is a lower limit of $\tilde{\mu}_{3}$ in general to avoid that
$c_{\Omega-}^{2}$ becomes negative around $r=r_{s}$. With the model parameters
$\alpha=2\times 10^{-4}$, $\eta=1$, and $\rho_{c}=15\rho_{0}$, for example,
this instability arises for $\tilde{\mu}_{3}<\tilde{\mu}_{3}^{{\rm
min}}=-0.22$.

The above argument shows that, for the coupling $\tilde{\mu}_{3}$ between
$\tilde{\mu}_{3}^{{\rm min}}$ and $\tilde{\mu}_{3}^{{\rm max}}$, there are NS
solutions with a nontrivial scalar profile free from ghost/Laplacian
instabilities. As $\tilde{\mu}_{3}$ approaches $\tilde{\mu}_{3}^{{\rm min}}$
or $\tilde{\mu}_{3}^{{\rm max}}$, the squared propagation speeds $c_{r3}^{2}$
and $c_{\Omega-}^{2}$ of $\delta\phi$ exhibit difference from those for
$\tilde{\mu}_{3}=0$.

###  IV.2 Quartic derivative and GB couplings

We proceed to a theory with the linear quartic derivative coupling
$G_{4}\supset\mu_{4}X$ besides the scalar-GB coupling $\alpha\xi(\phi)R_{\rm
GB}^{2}$. We note that the quartic coupling $\mu_{4}X$ in $G_{4}$ is
equivalent to the quintic coupling $-\mu_{4}\phi$ in $G_{5}$ Kobayashi  _et
al._ (2011). This theory is given by the action

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm Pl}^{2}}{2}R+\eta
X+\alpha\xi(\phi)R_{\rm
GB}^{2}+\mu_{4}XR+\mu_{4}\left\\{(\square\phi)^{2}-(\nabla_{\mu}\nabla_{\nu}\phi)(\nabla^{\mu}\nabla^{\nu}\phi)\right\\}\right]\,.$
(91)

The quartic derivative coupling alone does not give rise to asymptotically-
flat NS solutions Lehébel  _et al._ (2017), but the presence of scalar-GB
couplings modifies this no-hair property.

Then, the solutions expanded around $r=0$ are given by

$\displaystyle f$ $\displaystyle=$ $\displaystyle f_{c}+\frac{f_{c}[M_{\rm
Pl}^{2}\rho_{c}(1+3w_{c})-48\alpha\xi_{,\phi}(\phi_{c})\phi_{2}(w_{c}\rho_{c}-8\mu_{4}\phi_{2}^{2})]}{6[M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]^{2}}r^{2}+{\cal O}(r^{3})\,,$
(92) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{\rho_{c}+24\mu_{4}\phi_{2}^{2}}{3[M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]}r^{2}+{\cal O}(r^{3})\,,$ (93)
$\displaystyle P$ $\displaystyle=$ $\displaystyle
P_{c}-\frac{\rho_{c}(1+w_{c})[M_{\rm
Pl}^{2}\rho_{c}(1+3w_{c})-48\alpha\xi_{,\phi}(\phi_{c})\phi_{2}(w_{c}\rho_{c}-8\mu_{4}\phi_{2}^{2})]}{12[M_{\rm
Pl}^{2}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]^{2}}r^{2}+{\cal O}(r^{3})\,,$
(94)

where $\phi_{2}$ is the coefficient appearing in the scalar-field expansion
$\phi=\phi_{0}+\phi_{2}r^{2}+\cdots$. Since the leading-order contribution to
the scalar-field equation of motion around $r=0$ is highly nonlinear in
$\phi_{2}$, it is impossible to solve it for any value $\alpha$. So, we
restrict ourselves to the case $|\alpha|\ll 1$. For $|\alpha|\ll 1$,
$\phi_{2}$ satisfies the cubic-order equation

$\frac{9M_{\rm Pl}^{2}(\eta M_{\rm
Pl}^{2}+16\mu_{4}^{2}\phi_{2}^{2}-2w_{c}\rho_{c}\mu_{4})}{2\xi_{,\phi}(\phi_{c})\rho_{c}(1+3w_{c})(\rho_{c}+72\mu_{4}\phi_{2}^{2})}\phi_{2}=\alpha\,,$
(95)

at linear order in $\alpha$. In the limit $\alpha\to 0$ with a small
derivative coupling constant $\mu_{4}$, the only solution to Eq. (95) is
$\phi_{2}=0$. Hence the nonvanishing scalar-GB coupling $\alpha$ is required
to have the solution with a nontrivial scalar profile with $\phi_{2}\neq 0$.
On using these expanded solutions, the quantities associated with the
stability of odd-parity perturbations at $r=0$ are of the same form as Eq.
(51), where $x=\alpha\xi_{,\phi}(\phi_{c})\phi_{2}/M_{\rm Pl}^{2}$ and
$\phi_{2}$ satisfies the relation (95). The stability in the odd-parity sector
is ensured for $x<1/16$.

In the following, we will exploit the expansion with respect to the small
coupling $\alpha$. Around $r=0$, the no-ghost parameter of even-parity
perturbations has the dependence

${\cal K}=\left[8\phi_{2}^{2}M_{\rm Pl}^{2}\left(\eta M_{\rm
Pl}^{2}+144\phi_{2}^{2}\mu_{4}^{2}+2\rho_{c}\mu_{4}\right)-256\phi_{2}^{3}\xi_{,\phi}(\phi_{c})\left(\eta
M_{\rm Pl}^{2}+24\phi_{2}^{2}\mu_{4}^{2}-\mu_{4}\rho_{c}\right)\alpha+{\cal
O}(\alpha^{2})\right]r^{6}+{\cal O}(r^{7})\,.$ (96)

To avoid ghosts for small $\alpha$, the leading-order term
$8\phi_{2}^{2}M_{\rm Pl}^{2}(\eta M_{\rm
Pl}^{2}+144\phi_{2}^{2}\mu_{4}^{2}+2\rho_{c}\mu_{4})r^{6}$ in ${\cal K}$ needs
to be positive. This amounts to the condition

$\mu_{4}>-\frac{\eta M_{\rm Pl}^{2}+144\phi_{2}^{2}\mu_{4}^{2}}{2\rho_{c}}\,,$
(97)

which indicates the existence of a lower bound on negative values of
$\mu_{4}$. We note that $\phi_{2}$ depends on $\mu_{4}$ through Eq. (95). At
$r=0$, $c_{\Omega+}^{2}$ as well as $c_{r}^{2}$ and $c_{\Omega}^{2}$ are of
the same forms as Eq. (52), with $\phi_{2}$ satisfying Eq. (95). The other two
squared propagation speeds at $r=0$ are

$c_{r3}^{2}(r=0)=c_{\Omega-}^{2}(r=0)=\frac{\eta M_{\rm
Pl}^{2}+48\phi_{2}^{2}\mu_{4}^{2}-2w_{c}\rho_{c}\mu_{4}}{\eta M_{\rm
Pl}^{2}+144\phi_{2}^{2}\mu_{4}^{2}+2\rho_{c}\mu_{4}}+{\cal O}(\alpha)\,.$ (98)

Under the condition (97), the absence of Laplacian instability requires that

$\mu_{4}<\frac{\eta M_{\rm
Pl}^{2}+48\phi_{2}^{2}\mu_{4}^{2}}{2w_{c}\rho_{c}}\,,$ (99)

which indicates the existence of an upper bound on positive values of
$\mu_{4}$.

The solutions to $f$ and $\phi$ expanded at large distances are of the same
forms as Eqs. (64) and (66), respectively, while the solution to $h$ is

$h=1-\frac{2M}{r}+\frac{\eta\hat{\phi}_{1}^{2}}{2M_{\rm
Pl}^{2}r^{2}}+\frac{\eta M\hat{\phi}_{1}^{2}}{2M_{\rm
Pl}^{2}r^{3}}+\frac{2\hat{\phi}_{1}[(M^{2}\eta-3\mu_{4})\hat{\phi}_{1}+24M\alpha\xi_{,\phi}(\phi_{0})]}{3M_{\rm
Pl}^{2}r^{4}}+{\cal O}(r^{-5})\,.$ (100)

Up to next-to-leading order, the quantities $c_{r}^{2}$, $c_{\Omega}^{2}$,
${\cal K}$, and $c_{\Omega\pm}^{2}$ are the same as those in Eqs. (67), (68),
and (70), while the third radial propagation speed squared is modified to

$c_{r3}^{2}=1-\frac{2\mu_{4}\hat{\phi}_{1}^{2}}{M_{\rm Pl}^{2}r^{4}}+{\cal
O}(r^{-5})\,.$ (101)

Provided that $\eta>0$, the linear stability conditions are consistently
satisfied at spatial infinity.

For concreteness, we consider the linear scalar-GB coupling $\xi(\phi)=M_{\rm
Pl}r_{0}^{2}\phi$ besides the quartic derivative coupling. Since this theory
falls in a subclass of shift-symmetric Horndeski theories, we have
$J^{r}=(Q/r^{2})\sqrt{f/h}$. The constant $Q$ must be 0 to satisfy the
boundary conditions at $r=0$. Then, it follows that

$\phi^{\prime}(r)=-\frac{4\alpha f^{\prime}(h-1)M_{\rm Pl}r_{0}^{2}}{\eta
r^{2}f-2\tilde{\mu}_{4}r_{0}^{2}[rf^{\prime}h+f(h-1)]}\,,$ (102)

where

$\tilde{\mu}_{4}\equiv\frac{\mu_{4}}{r_{0}^{2}}\,.$ (103)

Substituting the expanded solutions (92) and (93) into Eq. (102), we find that
the field derivative behaves as $\phi^{\prime}(r)\propto r$ around $r=0$.
Substitution of Eqs. (64) and (100) into Eq. (102) gives
$\phi^{\prime}(r)=16\alpha M^{2}M_{\rm Pl}r_{0}^{2}/(\eta r^{5})$ at spatial
infinity and hence $\hat{\phi}_{1}=0$ in Eq. (66). Note that the condition for
the regularity at the center, $Q=0$, makes the leading scalar charge vanish at
spatial infinity (${\hat{\phi}}_{1}=0$) in the expansion of Eq. (66).

Figure 4:  (Left) Mass function ${\cal M}$ (normalized by the solar mass
$M_{\odot}$) and field derivative $\phi^{\prime}$ (normalized by $M_{\rm
Pl}/r_{0}$) versus $r/r_{0}$ for the quartic derivative coupling
$\tilde{\mu}_{4}=9\times 10^{-4}$, $\alpha=2\times 10^{-4}$, $\eta=1$, and
$\rho_{c}=15\rho_{0}$ (solid lines). The dotted lines show ${\cal M}$ and
$\phi^{\prime}$ for $\tilde{\mu}_{4}=0$, $\alpha=2\times 10^{-4}$, $\eta=1$,
and $\rho_{c}=15\rho_{0}$. (Right) $c_{r}^{2}$, $c_{\Omega}^{2}$,
$c_{r3}^{2}$, $c_{\Omega+}^{2}$, and $c_{\Omega-}^{2}$ versus $r/r_{0}$ for
$\tilde{\mu}_{4}=9\times 10^{-4}$, $\alpha=2\times 10^{-4}$, $\eta=1$, and
$\rho_{c}=15\rho_{0}$.

For given values of $\alpha$, $\rho_{c}$, $w_{c}$, and $\eta$, $\phi_{2}$ is
known from Eq. (95) as a function of $\tilde{\mu}_{4}$. Then, the minimum
value of $\tilde{\mu}_{4}$ can be found by Eq. (97). When $\alpha=2\times
10^{-4}$, $\rho_{c}=15\rho_{0}$, $w_{c}=0.48$, and $\eta=1$, we have
$\tilde{\mu}_{4}>-1.36\times 10^{-3}$. Numerically, we find that the values of
$c_{r3}^{2}$ and $c_{\Omega-}^{2}$ around $r=0$ become negative for
$\tilde{\mu}_{4}<-1.30\times 10^{-4}$. Hence the condition (97) gives a good
approximate criterion for the existence of NS solutions with a nontrivial
scalar profile consistent with linear stability conditions. At the background
level the negative coupling $\tilde{\mu}_{4}$ leads to tiny suppression of
$\phi^{\prime}(r)$, so the ADM mass of NS is only slightly increased. For
$\tilde{\mu}_{4}=-1\times 10^{-3}$ and $\alpha=2\times 10^{-4}$, we obtain
$M=2.008M_{\odot}$, which is close to the value $M=2.003M_{\odot}$ derived for
$\tilde{\mu}_{4}=0$ and $\alpha=2\times 10^{-4}$. As $\mu_{4}$ approaches the
lower bound (97), the leading-order contributions to $c_{r3}^{2}(r=0)$ and
$c_{\Omega-}^{2}(r=0)$ become highly superluminal. Outside the star, we
numerically confirm that all of the propagation speeds quickly approach 1.

The positive quartic coupling $\mu_{4}$ is constrained to be in the range
satisfying (99). With the model parameters $\alpha=2\times 10^{-4}$,
$\rho_{c}=15\rho_{0}$, $w_{c}=0.48$, and $\eta=1$, the condition (99) is
satisfied for any positive $\mu_{4}$. However, we need to caution that the
next-to-leading order correction to Eq. (98) gives rise to a negative term of
order $-0.1$. For $\tilde{\mu}_{4}>1\times 10^{-3}$, we numerically find that
$c_{r3}^{2}$ and $c_{\Omega-}^{2}$ become negative around the center of star.
Hence there is actually an upper bound of $\tilde{\mu}_{4}$ to avoid the
Laplacian instability of even-parity perturbations. The plots in Fig. 4
correspond to the coupling $\tilde{\mu}_{4}=9\times 10^{-4}$, in which case
$c_{r3}^{2}$ and $c_{\Omega-}^{2}$ are as close as 0.1 in the central region
of star. Unlike the case $\tilde{\mu}_{4}=0$ shown in Fig. 2,
$c_{\Omega-}^{2}$ grows smoothly as a function of $r$ toward the asymptotic
value 1. The behavior of other squared propagation speeds $c_{r}^{2}$,
$c_{\Omega}^{2}$, and $c_{\Omega+}^{2}$ is not much different from those seen
in Fig. 2. As we observe in the left panel of Fig. 4, the field derivative
inside the star is slightly enhanced in comparison to the case
$\tilde{\mu}_{4}=0$. The resulting ADM mass $M=1.995M_{\odot}$ is a bit
smaller than the value $M=2.003M_{\odot}$ obtained for $\tilde{\mu}_{4}=0$.

In summary, there are NS solutions with a nontrivial scalar profile free from
ghost/Laplacian instabilities in certain ranges of the coupling
$\tilde{\mu}_{4}$. For $\tilde{\mu}_{4}$ close to its lower limit,
$c_{r3}^{2}$ and $c_{\Omega-}^{2}$ are highly superluminal deep inside the
star. For $\tilde{\mu}_{4}$ close to its upper limit, $c_{r3}^{2}$ and
$c_{\Omega-}^{2}$ approach $+0$ around $r=0$.

###  IV.3 Nonminimal Ricci scalar and GB couplings

At the end of this section, we study the existence and stability of NS
solutions in theories with nonminimal Ricci scalar and scalar-GB couplings. We
incorporate a linear nonminimal coupling of the form $\lambda_{4}\phi R$ in
$G_{4}$. This also accommodates the dilatonic coupling ${\rm
e}^{\lambda_{4}\phi}R$ in the regime $|\lambda_{4}\phi|\ll 1$. The action in
such theories is given by

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm
Pl}^{2}}{2}R+\lambda_{4}\phi R+\eta X+\alpha\xi(\phi)R_{\rm GB}^{2}\right]\,.$
(104)

Provided that the scalar-GB coupling is suppressed on a weak gravitational
background, the dimensionless coupling constant
$\tilde{\lambda}_{4}=\lambda_{4}/M_{\rm Pl}$ is constrained to be
$\tilde{\lambda}_{4}\leq 2.5\times 10^{-3}$ from Solar System experiments
Tsujikawa  _et al._ (2008); De Felice and Tsujikawa (2010). On the strong
gravitational background, we will study the effect of nonminimal coupling on
the existence and the linear stability of NSs.

In the absence of the scalar-GB coupling, the squared propagation speeds are
Kase and Tsujikawa (2022)

$c_{r}^{2}=c_{\Omega}^{2}=c_{r3}^{2}=c_{\Omega+}^{2}=c_{\Omega-}^{2}=1\,,\qquad{\rm
for}\qquad\alpha=0\,,$ (105)

at any distance $r$. These values are the same as those in GR, but the
nonminimal coupling can give rise to NS solutions with a nontrivial scalar
profile with $\phi^{\prime}(r)\neq 0$ even in the absence of the scalar-GB
coupling.

Let us consider the case in which both nonminimal and scalar-GB couplings are
present. The solutions expanded around $r=0$ are given by

$\displaystyle f$ $\displaystyle=$ $\displaystyle f_{c}+\frac{f_{c}\\{(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c})[\rho_{c}(1+3w_{c})-12\lambda_{4}\phi_{2}]-48\alpha\xi_{,\phi}(\phi_{c})\phi_{2}(w_{c}\rho_{c}-8\lambda_{4}\phi_{2})\\}}{6[M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]^{2}}r^{2}+{\cal
O}(r^{3})\,,$ (106) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{\rho_{c}+12\lambda_{4}\phi_{2}}{3[M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]}r^{2}+{\cal
O}(r^{3})\,,$ (107) $\displaystyle P$ $\displaystyle=$ $\displaystyle
P_{c}-\frac{\rho_{c}(1+w_{c})\\{(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c})[\rho_{c}(1+3w_{c})-12\lambda_{4}\phi_{2}]-48\alpha\xi_{,\phi}(\phi_{c})\phi_{2}(w_{c}\rho_{c}-8\lambda_{4}\phi_{2})\\}}{12[M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}]^{2}}r^{2}+{\cal
O}(r^{3})\,,$ (108)

where $\phi_{2}$ appears as a coefficient in Eq. (42). Performing the
expansion with respect to a small scalar-GB coupling $\alpha$, $\phi_{2}$
obeys the second-order algebraic equation

$\frac{3(M_{\rm Pl}^{2}+2\lambda_{4}\phi_{c})\\{6\eta M_{\rm
Pl}^{2}\phi_{2}+\lambda_{4}[\rho_{c}(1-3w_{c})+12\eta\phi_{c}\phi_{2}]+36\lambda_{4}^{2}\phi_{2}\\}}{4\xi_{,\phi}(\phi_{c})[\rho_{c}^{2}(1+3w_{c})+72\lambda_{4}\phi_{2}(w_{c}\rho_{c}-6\lambda_{4}\phi_{2})]}=\alpha\,,$
(109)

up to linear order in $\alpha$. Choosing the regular branch in the limit
$\alpha\to 0$, we obtain

$\phi_{2}(\alpha=0)=\frac{\lambda_{4}\rho_{c}(3w_{c}-1)}{6[\eta\left(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}\right)+6\lambda_{4}^{2}]}\,.$ (110)

Thus the nonminimal coupling $\lambda_{4}\neq 0$ alone leads to a nonvanishing
value $\phi_{2}(\alpha=0)\neq 0$ around the center of body.

At $r=0$, the absence of ghost/Laplacian instabilities for odd-parity
perturbations requires that

$\displaystyle{\cal G}(r=0)=M_{\rm Pl}^{2}+2\lambda_{4}\phi_{c}>0\,,$ (111)
$\displaystyle{\cal H}(r=0)={\cal F}(r=0)=M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}-16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}>0\,,$ (112)

with the associated squared propagation speeds

$c_{r}^{2}(r=0)=c_{\Omega}^{2}(r=0)=\left[1-\frac{16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}}{M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}}\right]^{-1}\,.$ (113)

Around $r=0$, the no-ghost parameter ${\cal K}$ is in proportion to $r^{6}$
and hence

${\cal K}/r^{6}=8\phi_{2}^{2}\left(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}\right)\left[\eta\left(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}\right)+6\lambda_{4}^{2}\right]+{\cal
O}(\alpha)\,,$ (114)

where we performed the expansion with respect to small $\alpha$. Under the
condition (111), the absence of ghosts in the limit $\alpha\to 0$ requires
that

$\eta\left(M_{\rm Pl}^{2}+2\lambda_{4}\phi_{c}\right)+6\lambda_{4}^{2}>0\,.$
(115)

Provided that $\eta>0$ and $\lambda_{4}\phi_{c}>0$, the conditions (111),
(112), and (115) are automatically satisfied for $\alpha\to 0$. At $r=0$, we
also obtain

$\displaystyle c_{r3}^{2}(r=0)$ $\displaystyle=$ $\displaystyle
c_{\Omega-}^{2}(r=0)=1-\frac{16\lambda_{4}\xi_{,\phi}(\phi_{c})(\rho_{c}+P_{c}+6\lambda_{4}\phi_{2})}{(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c})[\eta\left(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}\right)+6\lambda_{4}^{2}]}\alpha+{\cal
O}(\alpha^{2})\,,$ (116) $\displaystyle c_{\Omega+}^{2}(r=0)$ $\displaystyle=$
$\displaystyle 1+\frac{16\alpha\xi_{,\phi}(\phi_{c})\phi_{2}}{M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{c}}+{\cal O}(\alpha^{2})\,.$ (117)

From Eqs. (113) and (117) we have
$c_{r}^{2}(r=0)=c_{\Omega}^{2}(r=0)=c_{\Omega+}^{2}(r=0)$ up to linear order
in $\alpha$. In the limit $\alpha\to 0$, all the above squared propagation
speeds approach 1.

The solutions expanded at spatial infinity are given by

$\displaystyle f$ $\displaystyle=$ $\displaystyle
1-\frac{2M}{r}+\frac{\lambda_{4}\hat{\phi}_{1}(4\eta MM_{\rm Pl}^{2}+8\eta
M\lambda_{4}\hat{\phi}_{0}+\eta\lambda_{4}\hat{\phi}_{1}+24M\lambda_{4}^{2})}{(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})[\eta(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})+6\lambda_{4}^{2}]r^{2}}+{\cal
O}(r^{-3})\,,$ (118) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{2(MM_{\rm
Pl}^{2}+2M\lambda_{4}\hat{\phi}_{0}-2\lambda_{4}\hat{\phi}_{1})}{(M_{\rm
Pl}^{2}+2\lambda_{4}\phi_{0})r}+\frac{{\hat{\phi}}_{1}[\eta^{2}M_{\rm
Pl}^{2}\hat{\phi}_{1}+2\eta\lambda_{4}(\eta\hat{\phi}_{0}\hat{\phi}_{1}+2MM_{\rm
Pl}^{2})+2\eta\lambda_{4}^{2}(4M\hat{\phi}_{0}+\hat{\phi}_{1})+24M\lambda_{4}^{3}]}{2(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})[\eta(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})+6\lambda_{4}^{2}]r^{2}}$ (119)
$\displaystyle+{\cal O}(r^{-3})\,,$ $\displaystyle\phi$ $\displaystyle=$
$\displaystyle\hat{\phi}_{0}+\frac{\hat{\phi}_{1}}{r}+\frac{\hat{\phi}_{1}[2\eta
MM_{\rm Pl}^{4}+\eta\lambda_{4}M_{\rm
Pl}^{2}(8M\hat{\phi}_{0}-3\hat{\phi}_{1})+\lambda_{4}^{2}\\{2\eta\hat{\phi}_{0}(4M\hat{\phi}_{0}-3\hat{\phi}_{1})+12MM_{\rm
Pl}^{2}\\}+12\lambda_{4}^{3}(2M\hat{\phi}_{0}-\hat{\phi}_{1})]}{2(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})[\eta(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})+6\lambda_{4}^{2}]r^{2}}$ (120)
$\displaystyle+{\cal O}(r^{-3})\,,$

where $M$, $\hat{\phi}_{0}$, $\hat{\phi}_{1}$ are integration constants. The
coupling $\alpha$ does not appear in Eqs. (118)-(120) up to the order of
$r^{-2}$, so the nonminimal coupling $\lambda_{4}$ provides larger
contributions to $f$, $h$, and $\phi$ far outside the star in comparison to
the scalar-GB term. On using these large-distance solutions, the dominant
contributions to ${\cal G}$, ${\cal H}$, and ${\cal F}$ are $M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0}$, so the linear stability conditions of
odd-parity perturbations are satisfied if $M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0}>0$. As for $c_{r}^{2}$ and
$c_{\Omega}^{2}$, the term $M_{\rm Pl}^{2}r^{3}$ appearing in the denominators
of Eq. (67) is modified to $(M_{\rm Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})r^{3}$.
In the even-parity sector, we have

$\displaystyle{\cal K}$ $\displaystyle=$ $\displaystyle
2\hat{\phi}_{1}^{2}\left(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0}\right)\left[\eta(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})+6\lambda_{4}^{2}\right]+{\cal
O}(r^{-1})\,,$ (121) $\displaystyle c_{r3}^{2}$ $\displaystyle=$
$\displaystyle
1-\frac{96\alpha\xi_{,\phi}(\phi_{c})\lambda_{4}^{2}\hat{\phi}_{1}}{(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})[\eta(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})+6\lambda_{4}^{2}]r^{3}}+{\cal
O}(r^{-4})\,.$ (122)

Under the condition $M_{\rm Pl}^{2}+2\lambda_{4}\hat{\phi}_{0}>0$, the ghosts
are absent for $\eta(M_{\rm
Pl}^{2}+2\lambda_{4}\hat{\phi}_{0})+6\lambda_{4}^{2}>0$. In comparison to Eq.
(69), the nonminimal coupling gives rise to a term proportional to $r^{-3}$ in
$c_{r3}^{2}-1$. The expressions of squared angular propagation speeds
$c_{\Omega\pm}^{2}$ are complicated, but they have the dependence
$|c_{\Omega\pm}^{2}-1|\propto|\alpha|/r^{3}$ under the small $\alpha$
expansion.

Figure 5:  (Left) Mass function ${\cal M}$ (normalized by the solar mass
$M_{\odot}$) and field derivative $\phi^{\prime}$ (normalized by $M_{\rm
Pl}/r_{0}$) versus $r/r_{0}$ for $\tilde{\lambda}_{4}=-2.5\times 10^{-3}$,
$\phi_{c}=M_{\rm Pl}$, $\alpha=2.5\times 10^{-4}$, $\eta=1$, and
$\rho_{c}=10\rho_{0}$ (solid lines). The dashed lines correspond to the case
$\alpha=0$, while the other model parameters are unchanged. (Right)
$c_{r}^{2}$, $c_{\Omega}^{2}$, $c_{r3}^{2}$, $c_{\Omega+}^{2}$, and
$c_{\Omega-}^{2}$ versus $r/r_{0}$ for the same model parameters as those used
for plotting solid lines in the left panel (i.e., both nonminimal and scalar-
GB couplings are present).

In the left panel of Fig. 5, we plot $\phi^{\prime}$ versus $r/r_{0}$ for
$\alpha=0$, $\tilde{\lambda}_{4}=\lambda_{4}/M_{\rm Pl}=-2.5\times 10^{-3}$,
$\phi_{c}=M_{\rm Pl}$, $\eta=1$, and $\rho_{c}=10\rho_{0}$ as a dashed line.
The field derivative increases around $r=0$ according to the relation
$\phi^{\prime}\simeq 2\phi_{2}r$, where $\phi_{2}$ is given by Eq. (110). In
Fig. 5, we observe that $\phi^{\prime}$ starts to decrease around the surface
of star and it joins the large-distance solution
$\phi^{\prime}\simeq-\hat{\phi}_{1}/r^{2}$ with $\hat{\phi}_{1}<0$. Note that
breaking the shift symmetry leads to the leading scalar charge
$\hat{\phi}_{1}\neq 0$ even if the regular boundary conditions are imposed at
the center. This hairy NS solution satisfies all the linear stability
conditions, with the squared propagation speeds given by Eq. (105). In this
case, the ADM mass of NS is $M=1.908M_{\odot}$ with the radius $r_{s}=10.97$
km. They are almost similar to the values $M=1.912M_{\odot}$ and $r_{s}=11.00$
km in GR with the same central density $\rho_{c}=10\rho_{0}$. Increasing the
value of $|\lambda_{4}\phi_{c}|$ further, it is possible to realize NSs with a
nontrivial scalar profile whose mass and radius exhibit notable difference
from those in GR. If we take into account Solar System constraints, however,
the coupling $\tilde{\lambda}_{4}$ should be less than the order $10^{-3}$.
Hence we do not consider the case in which the product
$|\tilde{\lambda}_{4}\phi_{c}|$ exceeds the order of $10^{-3}$.

If the scalar-GB coupling is present besides the nonminimal coupling, it is
possible to realize NS solutions with a nontrivial scalar profile as well. The
solid curves in the left panel of Fig. 5 correspond to the radial dependence
of $\phi^{\prime}$ and ${\cal M}$ for $\alpha=2.5\times 10^{-4}$,
$\tilde{\lambda}_{4}=-2.5\times 10^{-3}$, and $\phi_{c}=M_{\rm Pl}$. Inside
the NS, the field derivative is significantly enhanced in comparison to the
case $\alpha=0$. Around the surface of star $\phi^{\prime}$ starts to decrease
rapidly, but it enters the region with the radial dependence
$\phi^{\prime}\simeq-\hat{\phi}_{1}/r^{2}$ for $r\gtrsim r_{0}$. As we already
mentioned, this latter property is attributed to the fact that the
contribution to $\phi^{\prime}$ from the nonminimal coupling dominates over
that from the scalar-GB coupling at large distances. For $\alpha=2.5\times
10^{-4}$ the mass and radius of NS are found to be $M=1.873M_{\odot}$ and
$r_{s}=10.93$ km, both of which are slightly smaller than those for $\alpha=0$
mentioned above. In the right panel of Fig. 5, we show the five squared
propagation speeds versus $r/r_{0}$ for $\alpha=2.5\times 10^{-4}$,
$\tilde{\lambda}_{4}=-2.5\times 10^{-3}$, and $\phi_{c}=M_{\rm Pl}$. In
comparison to the case $\alpha=0$ where all the propagation speeds are
equivalent to 1, they are different from 1 deep inside the NS and approach 1
outside the star. For the model parameters used in Fig. 5, there are neither
ghost nor Laplacian instabilities for NS solutions with a nontrivial scalar
profile.

Provided that $|\tilde{\lambda}_{4}\phi_{c}|\lesssim 10^{-3}$, the scalar-GB
coupling $|\alpha|$ larger than the order $10^{-5}$ gives the dominant
contribution to $\phi^{\prime}$ inside the star. Then, for $|\alpha|\gtrsim
10^{-5}$, the background NS solution and its linear stability are not much
different from those for the scalar-GB coupling alone discussed in Sec. III.

##  V Regularized 4D-Einstein-Gauss-Bonnet gravity

In this section, we study the linear stability of NS solutions with a
nontrivial scalar profile in so-called “4DEGB gravity” arising from the
reduction of higher-dimensional GB theory to 4 dimensions. If we consider the
GB term ${\cal R}_{\rm GB}^{2}$ in spacetime dimensions $D$ higher than 4, the
field equations of motion following from the Lagrangian
$L=\sqrt{-g}\,\hat{\alpha}_{\rm GB}{\cal R}_{\rm GB}^{2}$ vanish in 4
dimensions Lovelock (1971). However, rescaling the GB coupling constant as
$\hat{\alpha}_{\rm GB}\to\alpha_{\rm GB}/(D-4)$ allows a possibility for
extracting contributions of the higher-dimensional GB term Glavan and Lin
(2020). Under such a rescaling, in the limit $D\to 4$, it is possible to
construct a regularized 4-dimensional theory by adding a counter-term to
eliminate divergent parts of the theory Fernandes  _et al._ (2020); Hennigar
_et al._ (2020). The other equivalent procedure is to perform a Kaluza-Klein
reduction of $D$-dimensional Einstein-GB gravity on a $(D-4)$-dimensional
maximally symmetric space with a vanishing spatial curvature Lu and Pang
(2020); Kobayashi (2020). The size of such a maximally symmetric space is
characterized by a scalar field $\phi$. The 4-dimensional action obtained from
the Kaluza-Klein reduction of $D$-dimensional Einstein-GB theory belongs to a
subclass of shift-symmetric Horndeski theories given by the coupling functions

$G_{2}=8\alpha_{\rm GB}X^{2}\,,\qquad G_{3}=8\alpha_{\rm GB}X\,,\qquad
G_{4}=1+4\alpha_{\rm GB}X\,,\qquad G_{5}=4\alpha_{\rm GB}\ln|X|\,.$ (123)

Notice that the standard kinetic term $X$ is absent in $G_{2}$. Throughout
this section, we use the unit $M_{\rm Pl}^{2}/2=1$.

From the scalar-field Eq. (12), we have $r^{2}\sqrt{f/h}\,J^{r}=Q={\rm
constant}$ and hence

$4\sqrt{\frac{h}{f}}\left(f^{\prime}+2\phi^{\prime}f\right)\left[1-h(1+r\phi^{\prime})^{2}\right]\alpha_{\rm
GB}=Q\,.$ (124)

Due to the regularity conditions $f^{\prime}(0)=0$ and $\phi^{\prime}(0)=0$ at
the center of star, it follows that $Q=0$. To satisfy Eq. (124) with $Q=0$ at
any radius $r$, we require that $1-h(1+r\phi^{\prime})^{2}=0$. The branch
where $\phi(r)$ decreases at spatial infinity (satisfying the asymptotic
flatness $h\to 1$ as $r\to\infty$) is given by

$\phi^{\prime}=\frac{1}{r}\left(\frac{1}{\sqrt{h}}-1\right)\,.$ (125)

Taking the $r$ derivative of this equation and substituting
$\phi^{\prime\prime}$ and Eq. (125) into Eqs. (9) and (10), we obtain

$\displaystyle h^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{2(h-1)[r^{2}+(h-1)\alpha_{\rm GB}]+\rho
r^{4}}{2r[r^{2}-2(h-1)\alpha_{\rm GB}]}\,,$ (126) $\displaystyle f^{\prime}$
$\displaystyle=$ $\displaystyle-f\frac{2(h-1)[r^{2}+(h-1)\alpha_{\rm
GB}]-Pr^{4}}{2hr[r^{2}-2(h-1)\alpha_{\rm GB}]}\,.$ (127)

Outside the star ($\rho=0=P$), there is the following analytic solution Lu and
Pang (2020); Fernandes  _et al._ (2021)

$h=f=1+\frac{r^{2}}{2\alpha_{\rm GB}}\left[1-\sqrt{1+\frac{8\alpha_{\rm
GB}M}{r^{3}}}\right]\,,$ (128)

where $M$ is an integration constant. At spatial infinity, this has the
asymptotic behavior $h=f=1-2M/r+{\cal O}(r^{-4})$ and
$\phi^{\prime}=M/r^{2}+{\cal O}(r^{-3})$.

Inside the NS, the integrated solutions to $h$ and $f$ depend on the fluid
EOS. For constant density $\rho$, we have the following exact solution Doneva
and Yazadjiev (2021)

$\displaystyle h$ $\displaystyle=$ $\displaystyle 1-\zeta r^{2}\,,$ (129)
$\displaystyle f$ $\displaystyle=$ $\displaystyle\frac{(1-\alpha_{\rm
GB}\zeta)^{2}}{4(1+2\alpha_{\rm GB}\zeta)^{2}}\left[\frac{3(1+\alpha_{\rm
GB}\zeta)}{1-\alpha_{\rm GB}\zeta}\sqrt{1-\zeta r_{s}^{2}}-\sqrt{1-\zeta
r^{2}}\right]^{2}\,,$ (130) $\displaystyle P$ $\displaystyle=$
$\displaystyle\rho\frac{(1-\alpha_{\rm GB}\zeta)[\sqrt{1-\zeta
r^{2}}-\sqrt{1-\zeta r_{s}^{2}}]}{3(1+\alpha_{\rm GB}\zeta)\sqrt{1-\zeta
r_{s}^{2}}-(1-\alpha_{\rm GB}\zeta)\sqrt{1-\zeta
r^{2}}}=\rho\left(\sqrt{\frac{1-\zeta r_{s}^{2}}{f}}-1\right)\,,$ (131)

where

$\zeta\equiv\frac{1}{2\alpha_{\rm GB}}\left(\sqrt{1+\frac{2}{3}\rho\alpha_{\rm
GB}}-1\right)\,.$ (132)

At the surface of star ($r=r_{s}$), the fluid pressure (131) vanishes.
Matching $h$ and $f$ with Eq. (128) at $r=r_{s}$, there is the relation
$M=\rho r_{s}^{3}/12$. In the following, we will study the linear stability of
NSs without assuming their EOSs.

The quantities associated with the linear stability against odd-parity
perturbations are given by

$\displaystyle{\cal G}$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{h}\,r^{4}-2r^{2}(\sqrt{h}-1)(Pr^{2}+2-2h+4\sqrt{h})\alpha_{\rm
GB}-4(\sqrt{h}+1)(\sqrt{h}-1)^{4}\alpha_{\rm
GB}^{2}}{r^{2}\sqrt{h}\,[r^{2}+2(1-h)\alpha_{\rm GB}]}\,,$ (133)
$\displaystyle{\cal F}$ $\displaystyle=$
$\displaystyle\frac{2\\{r^{4}+r^{2}[\rho r^{2}+2(h-1)]\alpha_{\rm
GB}-2(h-1)^{2}\alpha_{\rm GB}^{2}\\}}{r^{2}[r^{2}+2(1-h)\alpha_{\rm GB}]}\,,$
(134) $\displaystyle{\cal H}$ $\displaystyle=$ $\displaystyle
2+\frac{4(1-h)}{r^{2}}\alpha_{\rm GB}\,.$ (135)

In the limit that $\alpha_{\rm GB}\to 0$ we have ${\cal G}={\cal F}={\cal
H}=2$, so the linear stability against odd-parity perturbations is ensured for
the small GB coupling. On using the background Eqs. (126) and (127), it
follows that

${\cal K}=0\,,$ (136)

at any radius $r>0$. Inside the star, the term $(2{\cal P}_{1}-{\cal
F})h\mu^{2}$ exactly cancels the contribution $-2{\cal H}^{2}r^{4}(\rho+P)$,
while, outside the star, ${\cal K}=(2{\cal P}_{1}-{\cal F})h\mu^{2}=0$. The
fact that ${\cal K}$ vanishes everywhere, which is mostly attributed to the
absence of a standard kinetic term in $G_{2}$, is the signal of a strong
coupling problem. Indeed, the squared radial propagation speed associated with
the stability of scalar-field perturbation in the even-parity sector yields

$c_{r3}^{2}\to\infty\,.$ (137)

This divergent property of $c_{r3}^{2}$ arises from the existence of term
${\cal K}$ in the denominator of Eq. (22). We note that the product ${\cal
K}c_{r3}^{2}$ is finite. Exploiting the exact solution (128) outside the star
and performing the expansion with respect to $1/r$ at spatial infinity, we
obtain the dependence

${\cal K}c_{r3}^{2}=-\frac{256M^{3}\alpha_{\rm GB}}{r^{3}}+{\cal
O}(r^{-4})\,.$ (138)

For $\alpha_{\rm GB}<0$, the leading-order contribution to ${\cal
K}c_{r3}^{2}$ is positive at large distances.

The quantities $B_{1}$ and $B_{2}$ also diverge, so this leads to the
divergence of angular propagation speeds in the even-parity sector. Expanding
the product ${\cal K}B_{2}$ at spatial infinity, it follows that

${\cal K}B_{2}=\frac{128M^{3}\alpha_{\rm GB}}{r^{3}}+{\cal O}(r^{-4})\,.$
(139)

Dividing Eq. (138) by Eq. (139), we obtain

$\frac{c_{r3}^{2}}{B_{2}}=-2+{\cal O}(r^{-1})\,.$ (140)

Since the signs of $c_{r3}^{2}$ and $B_{2}$ are different from each other,
either of the linear stability conditions (22) or (26) is violated at large
distances. Thus, the NS solutions in 4DEGB gravity not only suffer from the
strong coupling problem but also the Laplacian instability of even-parity
perturbations. These problems also persist for hairy BHs given by the line
element (128) present in 4DEGB theory. In the case of BHs, there is also the
instability of even-parity perturbations in the vicinity of the event horizon
Tsujikawa (2022).

##  VI $F(R_{\rm GB}^{2})$ gravity

Finally, we study a modified GB gravity in which the 4-dimensional action
contains an arbitrary function $F$ of the GB term $R_{\rm GB}^{2}$ besides the
Einstein-Hilbert term. This theory is given by the action

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm Pl}^{2}}{2}R+F(R_{\rm
GB}^{2})\right]\,,$ (141)

which is equivalent to Kobayashi  _et al._ (2011)

${\cal S}=\int{\rm d}^{4}x\sqrt{-g}\left[\frac{M_{\rm
Pl}^{2}}{2}R+\xi(\phi)R_{\rm GB}^{2}-V(\phi)\right]\,,$ (142)

where

$\varphi\equiv R_{\rm GB}^{2}\,,\qquad\phi\equiv M_{\rm
Pl}r_{0}^{4}\varphi\,,\qquad\xi(\phi)\equiv F_{,\varphi}\,,\qquad
V(\phi)\equiv\varphi\xi-F\,.$ (143)

The action (142) belongs to a subclass of Horndeski theories given by the
coupling functions

$\displaystyle G_{2}=-V(\phi)+8\xi^{(4)}(\phi)X^{2}(3-\ln|X|)\,,\qquad
G_{3}=4\xi^{(3)}(\phi)X(7-3\ln|X|)\,,$ $\displaystyle G_{4}=\frac{M_{\rm
Pl}^{2}}{2}+4\xi^{(2)}(\phi)X(2-\ln|X|)\,,\qquad
G_{5}=-4\xi^{(1)}(\phi)\ln|X|\,.$ (144)

Hence there is no standard scalar kinetic term in $F(R_{\rm GB}^{2})$ gravity.

Let us focus on the power-law $F(R_{\rm GB}^{2})$ models given by

$F(R_{\rm GB}^{2})=\beta\left(R_{\rm GB}\right)^{n}\,,$ (145)

where $\beta$ and $n$ are constants. We consider the positive integers $n$ in
the range $n\geq 2$. Then, the GB coupling function $\xi(\phi)$ and scalar
potential $V(\phi)$ in the action (142) yield

$\xi(\phi)=\alpha n\frac{r_{0}^{2}M_{\rm Pl}^{2}}{2}\left(\frac{\phi}{M_{\rm
Pl}}\right)^{n-1}\,,\qquad V(\phi)=\alpha(n-1)\frac{M_{\rm
Pl}^{2}}{2r_{0}^{2}}\left(\frac{\phi}{M_{\rm Pl}}\right)^{n}\,,$ (146)

where $\alpha\equiv 2r_{0}^{2-4n}M_{\rm Pl}^{-2}\beta$ is a dimensionless
coupling.

From the scalar-field Eq. (12), we obtain

$\alpha\phi^{n-2}\left\\{r^{2}f^{2}\phi-2\left[ff^{\prime}h^{\prime}(3h-1)+f^{\prime
2}h(1-h)+2f^{\prime\prime}fh(h-1)\right]M_{\rm Pl}r_{0}^{4}\right\\}=0\,.$
(147)

For $n=2$, there is only one branch characterized by

$\phi=\frac{2\left[ff^{\prime}h^{\prime}(3h-1)+f^{\prime
2}h(1-h)+2f^{\prime\prime}fh(h-1)\right]M_{\rm Pl}r_{0}^{4}}{r^{2}f^{2}}\,.$
(148)

For $n\geq 3$, we also have the no-hair branch $\phi=0$ besides (148).

Let us first consider the power $n=2$. Around the center of star, the
solutions consistent with the boundary conditions at $r=0$ are

$\displaystyle f$ $\displaystyle=$ $\displaystyle f_{c}+\frac{f_{c}[M_{\rm
Pl}r_{0}^{2}\rho_{c}(1+3w_{c})-M_{\rm
Pl}\phi_{c}^{2}\alpha-24(2w_{c}\rho_{c}r_{0}^{2}-\phi_{c}^{2}\alpha)r_{0}^{2}\phi_{2}\alpha]}{6M_{\rm
Pl}r_{0}^{2}(M_{\rm Pl}-16r_{0}^{2}\phi_{2}\alpha)^{2}}r^{2}+{\cal
O}(r^{3})\,,$ (149) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{2\rho_{c}r_{0}^{2}+\alpha\phi_{c}^{2}}{6M_{\rm Pl}r_{0}^{2}(M_{\rm
Pl}-16\alpha r_{0}^{2}\phi_{2})}r^{2}+{\cal O}(r^{3})\,,$ (150)
$\displaystyle\phi$ $\displaystyle=$
$\displaystyle\phi_{c}+\phi_{2}r^{2}+{\cal O}(r^{3})\,,$ (151) $\displaystyle
P$ $\displaystyle=$ $\displaystyle P_{c}-\frac{\rho_{c}(1+w_{c})\\{M_{\rm
Pl}r_{0}^{2}\rho_{c}(1+3w_{c})-M_{\rm
Pl}\phi_{c}^{2}\alpha-24(2w_{c}\rho_{c}r_{0}^{2}-\phi_{c}^{2}\alpha)r_{0}^{2}\phi_{2}\alpha\\}}{12M_{\rm
Pl}r_{0}^{2}(M_{\rm Pl}-16r_{0}^{2}\phi_{2}\alpha)^{2}}r^{2}+{\cal
O}(r^{3})\,.$ (152)

The equation to determine $\phi_{2}$ is given by

$-2r_{0}^{2}\rho_{c}\left(1+3w_{c}\right)+2\alpha\phi_{c}^{2}+\frac{48r_{0}^{2}\alpha\phi_{2}\left(2r_{0}^{2}w_{c}\rho_{c}-\alpha\phi_{c}^{2}\right)}{M_{\rm
Pl}}-\frac{3\left(M_{\rm
Pl}-16r_{0}^{2}\alpha\phi_{2}\right)^{3}\phi_{c}}{2r_{0}^{2}\rho_{c}+\alpha\phi_{c}^{2}}=0\,.$
(153)

Since Eq. (153) is a cubic-order equation for $\phi_{2}$, there are three
solutions. The two of them are imaginary solutions, and the remaining one is
real. In the limit $\alpha\to 0$, the real solution is approximately given by

$\displaystyle\phi_{2}$ $\displaystyle=$ $\displaystyle\frac{1}{48M_{\rm
Pl}r_{0}^{2}\alpha\phi_{c}}\left[3M_{\rm Pl}^{2}\phi_{c}-2\times 6^{2/3}M_{\rm
Pl}w_{c}\phi_{c}\left(\frac{r_{0}^{8}\rho_{c}^{4}}{3M_{\rm
Pl}^{3}\phi_{c}^{2}+\sqrt{3M_{\rm
Pl}^{3}\phi_{c}^{3}\left(16r_{0}^{4}w_{c}^{3}\rho_{c}^{2}+3M_{\rm
Pl}^{3}\phi_{c}\right)}}\right)^{1/3}\right.$ (154)
$\displaystyle\left.+6^{1/3}r_{0}\left\\{r_{0}\rho_{c}^{2}\left(3M_{\rm
Pl}^{3}\phi_{c}^{2}+\sqrt{3M_{\rm
Pl}^{3}\phi_{c}^{3}\left(16r_{0}^{4}w_{c}^{3}\rho_{c}^{2}+3M_{\rm
Pl}^{3}\phi_{c}\right)}\right)\right\\}^{1/3}\right]+{\cal O}(\alpha^{0}).$

Thus, the scalar field is divergent as $\alpha\to 0$. In general, in the limit
$\alpha\to 0$, the same type of divergence can be also observed for $n\geq 3$.
This indicates a pathology of the interior solutions in the power-law
$f(R_{\rm GB}^{2})$ model.

The solutions expanded far outside the star are

$\displaystyle f$ $\displaystyle=$ $\displaystyle
1-\frac{2M}{r}+\frac{1024M^{3}r_{0}^{6}\alpha}{r^{9}}+{\cal O}(r^{-{10}})\,,$
(155) $\displaystyle h$ $\displaystyle=$ $\displaystyle
1-\frac{2M}{r}+\frac{4608M^{3}r_{0}^{6}\alpha}{r^{9}}+{\cal O}(r^{-{10}})\,,$
(156) $\displaystyle\phi$ $\displaystyle=$ $\displaystyle\frac{48M^{2}M_{\rm
Pl}r_{0}^{4}}{r^{6}}-\frac{1216512M^{4}M_{\rm
Pl}r_{0}^{10}\alpha}{r^{14}}+{\cal O}(r^{-15})\,.$ (157)

The Schwarzschild metrics receive corrections from the coupling $\alpha$ at
the order of $r^{-9}$. If there are NSs with a nontrivial scalar profile, the
interior solutions (149)-(151) should be joined with the exterior solutions
(155)-(157). On using the large-distance solutions (155)-(157), we obtain

${\cal K}=-\frac{63700992M^{6}M_{\rm
Pl}^{6}r_{0}^{12}\alpha^{2}}{r^{16}}+{\cal O}(r^{-17})\,.$ (158)

The leading-order term of ${\cal K}$ is negative, and hence there is a ghost
instability issue at large distances. Moreover, since ${\cal K}$ is suppressed
by a high power ${\cal O}(r^{-16})$, it quickly approaches 0 for increasing
$r$. Then the above solution also has a strong coupling problem in the
asymptotic region. Note that in Ref. Minamitsuji  _et al._ (2022a) the
definition of ${\cal K}$ is $2{\cal P}_{1}-{\cal F}$, in which case ${\cal
K}\propto r^{-18}$ as consistent with Eq. (158). In Ref. Minamitsuji  _et al._
(2022a) a small $\alpha$ expansion was used for deriving the background BH
solution, which means that our result (158) is more general.

For $n\geq 3$, there is also the branch (148) of a nonvanishing scalar field.
For this branch, the solutions in the vicinity of $r=0$ are similar to Eqs.
(149)-(152) with some modifications of coefficients. At large distances, the
leading-order solutions are Schwarzschild metric components as in Eqs.
(155)-(157), in which case we obtain

${\cal K}=-\frac{3^{2n+1}\cdot 256^{n}M^{2(2n-1)}M_{\rm
Pl}^{6}n^{2}(n-1)^{2}r_{0}^{4(2n-1)}\alpha^{2}}{r^{4(3n-2)}}+{\cal
O}(r^{7-12n})\,.$ (159)

Thus the leading-order term of ${\cal K}$ is negative, with a rapid decrease
of ${\cal K}$ toward 0 at large distances. To compute the quantity ${\cal K}$
above, we have not used the expansion with respect to a small coupling
$\alpha$. Hence the ghost instability and asymptotic strong coupling problem
of NS solutions with a nontrivial scalar profile given by the branch (148) are
generally present for an arbitrary nonvanishing coupling $\alpha$.

##  VII Conclusions

In gravitational theories with a coupling to GB curvature invariant $R_{\rm
GB}^{2}$, we studied the existence and stability of NS solutions with a
nontrivial profile of the scalar field on a static and spherically symmetric
background. For this purpose, we exploited conditions for avoiding
ghost/Laplacian instabilities of odd- and even-parity perturbations with high
radial and angular momentum modes Kase and Tsujikawa (2022). These linear
stability conditions in full Horndeski theories, which are summarized in Sec.
II, can be applied not only to BHs, i.e., the vacuum case, but also to NSs,
i.e., the case with matter fluids.

The scalar-GB coupling $\alpha\xi(\phi)R_{\rm GB}^{2}$ gives rise to NSs
endowed with a nontrivial scalar profile both inside and outside the star. In
Sec. III, we derived the approximated background solutions by using the
expansion both around $r=0$ and at spatial infinity. We then studied the
linear stability of them. The no-ghost condition requires that the theory has
to contain a positive scalar kinetic term of the scalar field., i.e.,
$\eta>0$. For the existence of NS solutions with a nontrivial scalar profile
free from instabilities at $r=0$, we derived the upper bound (57) on the
dimensionless coupling constant $\alpha$. For the linear scalar-GB coupling
$\alpha_{\rm GB}\phi R_{\rm GB}^{2}$, this bound translates to
$\sqrt{|\alpha_{\rm GB}|}<0.7~{}{\rm km}$ to realize maximum masses of NSs for
typical EOSs. This theoretical limit is stronger than those observationally
constrained from the GW measurements emitted from binaries containing NSs.
Numerically, we confirmed that the bound (57) is sufficiently accurate for
ensuring the existence of NSs with a nontrivial scalar profile consistent with
all the linear stability conditions. As we see in Fig. 2, the squared
propagation speeds of odd- and even-parity perturbations deviate from 1 deep
inside the NS, but all of them approach 1 outside the star. Note that the
propagation speeds different from the speed of light do not imply any
acausality and instead determine the causal boundaries locally at each
position.

In Sec. IV, we incorporated several regular Horndeski couplings besides the
scalar-GB coupling $\alpha\xi(\phi)R_{\rm GB}^{2}$. In the presence of a
canonical kinetic term $\eta X$, the cubic-order derivative interaction of the
scalar field $G_{3}\supset\mu_{3}X$ or the quartic-order derivative coupling
$G_{4}\supset\mu_{4}X$ alone do not give rise to asymptotically-flat NS
solutions with a nontrivial scalar profile Lehébel  _et al._ (2017). We
confirmed that the presence of the GB coupling is mandatory to obtain the NS
solutions with a nontrivial profile of the scalar field.

In the case where the scalar-GB coupling is present, we showed the existence
of NS solutions with a nontrivial scalar profile free from ghost/Laplacian
instabilities for certain ranges of the coupling constants $\mu_{3}$ and
$\mu_{4}$. As we see in Figs. 3 and 4, these new couplings do not lead to
significant changes to the ADM mass of NSs, but the squared propagation speeds
$c_{r3}^{2}$ and $c_{\Omega-}^{2}$ inside the star are subject to
modifications. We also showed that NS solutions with a nontrivial scalar
profile exist for a linear nonminimal coupling $G_{4}\supset\lambda_{4}\phi$
alone. Under local gravity constraints on the coupling constant $\lambda_{4}$,
however, the scalar-GB coupling gives dominant contributions to the background
scalar-field profile as well as the linear stability of NSs for $|\alpha|$
close to its upper limit.

In Sec. V, we addressed the linear stability of NS solutions with a nontrivial
scalar profile in regularized 4DEGB gravity. After the Kaluza-Klein reduction
of $D$-dimensional Einstein-GB theory on a flat internal space, the resulting
4-dimensional action belongs to a subclass of shift-symmetric Horndeski
theories. For the solution with a nontrivial scalar profile (125), we showed
that ${\cal K}=0$ and $c_{r3}^{2}\to\infty$ at any radius $r$, and hence it is
plagued by a strong coupling problem. Moreover the leading-order contribution
to $c_{r3}^{2}/B_{2}$ is $-2$ at large distances, so there is also the
Laplacian instability for even-parity perturbations. Along with the fact that
the BH solutions in regularized 4DEGB gravity are also unstable Tsujikawa
(2022), there are no stable compact objects endowed with a nontrivial scalar
profile in this theory.

In Sec. VI, we considered the power-law $F(R_{\rm GB}^{2})$ models with the
Einstein-Hilbert term, which are equivalent to the action (142) with
$\xi(\phi)$ and $V(\phi)$ given by Eq. (146). For $n\geq 2$, there is a
nonvanishing scalar-field branch characterized by Eq. (148). If NSs with a
nontrivial scalar profile are present, the interior solutions should join the
large-distance exterior solutions (155)-(157). On using the latter, we find
that the leading-order term of ${\cal K}$ is negative. For increasing $r$,
${\cal K}$ rapidly decreases toward $-0$ with a large negative power-law
dependence on $r$. Hence the hairy branch (148) is excluded by the problems of
ghost instability and strong coupling at spatial infinity. This is analogous
to what was found for BHs in the same theory Minamitsuji  _et al._ (2022a). In
both regularized 4DEGB theory and $F(R_{\rm GB}^{2})$ gravity, the
instabilities of NSs with a nontrivial scalar profile arise from the unhealthy
propagation of $\delta\phi$ associated with the absence of a canonical kinetic
term. Moreover, we found that in power-law $F(R_{\rm GB}^{2})$ models, the
interior solution of the scalar field suffers from the divergence in the zero-
coupling limit, which indicates an intrinsic pathology.

In summary, we have shown that NSs with a nontrivial scalar profile consistent
with the linear stability conditions are present for the scalar-GB coupling
$\alpha\xi(\phi)R_{\rm GB}^{2}$ besides regular Horndeski coupling functions.
As we studied in Sec. IV.3, nonminimal couplings with the Ricci scalar
$G_{4}(\phi)R$ also give rise to NS solutions with a nontrivial scalar profile
even without the scalar-GB coupling. The latter includes NSs with spontaneous
scalarization, which can occur for a large nonminimal coupling constant of
order unity Damour and Esposito-Farese (1993); Harada (1998); Novak (1998). It
will be of interest to extend our linear stability analysis to such large
nonminimal coupling regimes by taking the scalar-GB coupling into account.

## Acknowledgements

MM was supported by the Portuguese national fund through the Fundação para a
Ciência e a Tecnologia (FCT) in the scope of the framework of the Decree-Law
57/2016 of August 29, changed by Law 57/2017 of July 19, and the Centro de
Astrofísica e Gravitação (CENTRA) through the Project No. UIDB/00099/2020. MM
also would like to thank Yukawa Institute for Theoretical Physics (under the
Visitors Program of FY2022) and Department of Physics of Waseda University for
their hospitality. ST was supported by the Grant-in-Aid for Scientific
Research Fund of the JSPS Nos. 19K03854 and 22K03642.

##  Appendix A Coefficients in the background equations

The coefficients in Eqs. (9)-(11) are given by

$\displaystyle A_{1}=-h^{2}(G_{3,X}-2G_{4,\phi X})\phi^{\prime
2}-2G_{4,\phi}h\,,\qquad A_{2}=2h^{3}(2G_{4,XX}-G_{5,\phi X})\phi^{\prime
3}-4h^{2}(G_{4,X}-G_{5,\phi})\phi^{\prime}\,,$ $\displaystyle
A_{3}=-h^{4}G_{5,XX}\phi^{\prime 4}+h^{2}G_{5,X}(3h-1)\phi^{\prime 2}\,,\qquad
A_{4}=h^{2}(2G_{4,XX}-G_{5,\phi X})\phi^{\prime
4}+h(3G_{5,\phi}-4G_{4,X})\phi^{\prime 2}-2G_{4}\,,$ $\displaystyle
A_{5}=-\frac{1}{2}\left[G_{5,XX}h^{3}{\phi^{\prime}}^{5}-hG_{5,X}(5h-1)\phi^{\prime
3}\right]\,,\qquad A_{6}=h(G_{3,\phi}-2G_{4,\phi\phi})\phi^{\prime
2}+G_{2}\,,$ $\displaystyle A_{7}=-2h^{2}(2G_{4,\phi
X}-G_{5,\phi\phi})\phi^{\prime 3}-4G_{4,\phi}h\phi^{\prime}\,,\qquad
A_{8}=G_{5,\phi X}h^{3}\phi^{\prime
4}-h(2G_{4,X}h-G_{5,\phi}h-G_{5,\phi})\phi^{\prime 2}-2G_{4}(h-1)\,,$
$\displaystyle A_{9}=-h(G_{2,X}-G_{3,\phi})\phi^{\prime 2}-G_{2}\,,\qquad
A_{10}=\frac{1}{2}G_{5,\phi X}h^{3}\phi^{\prime
4}-\frac{1}{2}h^{2}(2G_{4,X}-G_{5,\phi})\phi^{\prime 2}-G_{4}h\,.$ (160)

The coefficients in Eq. (14) are

$\displaystyle\lambda_{1}=-\left(h^{\prime}+\frac{4h}{r}+\frac{f^{\prime}h}{f}\right)\phi^{\prime}-2h\phi^{\prime\prime}\,,\qquad\lambda_{2}=-h\phi^{\prime
2}\,,\qquad\lambda_{3}=\frac{1}{2}h\phi^{\prime
2}\left(h^{\prime}\phi^{\prime}+2h\phi^{\prime\prime}\right)\,,$
$\displaystyle\lambda_{4}=\frac{2}{r^{2}}(1-h-rh^{\prime})+\frac{hf^{\prime
2}}{2f^{2}}-\frac{r(2f^{\prime\prime}h+f^{\prime}h^{\prime})+4f^{\prime}h}{2fr}\,,$
$\displaystyle\lambda_{5}=h\phi^{\prime}\left[\left(\frac{8h^{\prime}}{r}+\frac{6h}{r^{2}}-\frac{f^{\prime
2}h}{2f^{2}}+\frac{(f^{\prime\prime}r+6f^{\prime})h+2rf^{\prime}h^{\prime}}{fr}\right)\phi^{\prime}+3h\left(\frac{f^{\prime}}{f}+\frac{4}{r}\right)\phi^{\prime\prime}\right]\,,$
$\displaystyle\lambda_{6}=h^{2}\phi^{\prime
3}\left(\frac{f^{\prime}}{f}+\frac{4}{r}\right)\,,\qquad\lambda_{7}=-\frac{1}{2}h^{2}\phi^{\prime
3}\left(\frac{f^{\prime}}{f}+\frac{4}{r}\right)\left(h^{\prime}\phi^{\prime}+2h\phi^{\prime\prime}\right)\,,$
$\displaystyle\lambda_{8}=\frac{1}{r^{2}}\left[h^{\prime}(3h-1)\phi^{\prime}+2h(h-1)\phi^{\prime\prime}\right]-\frac{f^{\prime
2}h^{2}\phi^{\prime}}{f^{2}r}+\frac{1}{fr^{2}}\left[(2f^{\prime\prime}r+3f^{\prime})h^{2}\phi^{\prime}+f^{\prime}h(3rh^{\prime}-1)\phi^{\prime}+2f^{\prime}h^{2}r\phi^{\prime\prime}\right]\,,$
$\displaystyle\lambda_{9}=\frac{h\phi^{\prime
2}}{fr^{2}}\left[f(h-1)+f^{\prime}hr\right]\,,$
$\displaystyle\lambda_{10}=-\frac{h\phi^{\prime
2}}{2r^{2}}\left[10h^{2}\phi^{\prime\prime}+h(7h^{\prime}\phi^{\prime}-2\phi^{\prime\prime})-h^{\prime}\phi^{\prime}\right]+\frac{f^{\prime
2}h^{3}\phi^{\prime 3}}{2f^{2}r}-\frac{h^{2}\phi^{\prime
2}}{2fr^{2}}\left[(2f^{\prime\prime}r+4f^{\prime})h\phi^{\prime}+10f^{\prime}hr\phi^{\prime\prime}+7f^{\prime}h^{\prime}r\phi^{\prime}\right],$
$\displaystyle\lambda_{11}=-\frac{h^{3}\phi^{\prime
4}}{fr^{2}}(rf^{\prime}+f)\,,\qquad\lambda_{12}=\frac{h^{3}\phi^{\prime
4}}{2fr^{2}}(rf^{\prime}+f)(h^{\prime}\phi^{\prime}+2h\phi^{\prime\prime})\,.$
(161)

##  Appendix B Coefficients in the perturbation equations

The quantities appearing in the linear stability conditions (19), (22) and
(23) are

$\displaystyle a_{1}$ $\displaystyle=$
$\displaystyle\sqrt{fh}\left[\left\\{G_{4,\phi}+\frac{1}{2}h(G_{3,X}-2G_{4,\phi
X})\phi^{\prime
2}\right\\}r^{2}+2h\phi^{\prime}\left\\{G_{4,X}-G_{5,\phi}-\frac{1}{2}h(2G_{4,XX}-G_{5,\phi
X})\phi^{\prime 2}\right\\}r\right.$
$\displaystyle\left.+\frac{1}{2}G_{5,XX}h^{3}\phi^{\prime
4}-\frac{1}{2}G_{5,X}h(3h-1)\phi^{\prime 2}\right]\,,$ $\displaystyle c_{2}$
$\displaystyle=$
$\displaystyle\sqrt{fh}\left[\left\\{\frac{1}{2f}\left(-\frac{1}{2}h(3G_{3,X}-8G_{4,\phi
X})\phi^{\prime 2}+\frac{1}{2}h^{2}(G_{3,XX}-2G_{4,\phi XX})\phi^{\prime
4}-G_{4,\phi}\right)r^{2}\right.\right.$
$\displaystyle\left.\left.-{\frac{h\phi^{\prime}}{f}}\left(\frac{1}{2}{h^{2}(2G_{4,XXX}-G_{5,\phi
XX})\phi^{\prime 4}}-\frac{1}{2}{h(12G_{4,XX}-7G_{5,\phi X})\phi^{\prime
2}}+3(G_{4,X}-G_{5,\phi})\right)r\right.\right.$
$\displaystyle\left.\left.+\frac{h\phi^{\prime
2}}{4f}\left(G_{5,XXX}h^{3}\phi^{\prime 4}-G_{5,XX}h(10h-1)\phi^{\prime
2}+3G_{5,X}(5h-1)\right)\right\\}f^{\prime}\right.$
$\displaystyle\left.+\phi^{\prime}\left\\{\frac{1}{2}G_{2,X}-G_{3,\phi}-\frac{1}{2}h(G_{2,XX}-G_{3,\phi
X})\phi^{\prime 2}\right\\}r^{2}\right.$
$\displaystyle\left.+2\left\\{-\frac{1}{2}h(3G_{3,X}-8G_{4,\phi
X})\phi^{\prime 2}+\frac{1}{2}h^{2}(G_{3,XX}-2G_{4,\phi XX})\phi^{\prime
4}-G_{4,\phi}\right\\}r\right.$
$\displaystyle\left.-\frac{1}{2}h^{3}(2G_{4,XXX}-G_{5,\phi XX})\phi^{\prime
5}+\frac{1}{2}h\left\\{2\left(6h-1\right)G_{4,XX}+\left(1-7h\right)G_{5,\phi
X}\right\\}\phi^{\prime 3}-(3h-1)(G_{4,X}-G_{5,\phi})\phi^{\prime}\right]\,,$
$\displaystyle c_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\frac{\sqrt{f}}{\sqrt{h}}\left[{\frac{h\phi^{\prime}}{f}\left\\{2G_{4,X}-2G_{5,\phi}-h(2G_{4,XX}-G_{5,\phi
X})\phi^{\prime 2}-{\frac{h\phi^{\prime}(3G_{5,X}-G_{5,XX}\phi^{\prime
2}h)}{r}}\right\\}}f^{\prime}\right.$ (162)
$\displaystyle\left.+4G_{4,\phi}+2h(G_{3,X}-2G_{4,\phi X})\phi^{\prime
2}+{\frac{4h(G_{4,X}-G_{5,\phi})\phi^{\prime}-2h^{2}(2G_{4,XX}-G_{5,\phi
X})\phi^{\prime 3}}{r}}\right]\,,$

and

$\displaystyle\beta_{0}$ $\displaystyle=$
$\displaystyle\phi^{\prime}a_{1}+r\sqrt{fh}{\cal H}\,,$ (163)
$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\phi^{\prime 2}\sqrt{fh}{\cal
H}e_{4}-\phi^{\prime}\left(\sqrt{fh}{\cal
H}\right)^{\prime}c_{4}+\frac{\sqrt{fh}}{2}\left[\left({\frac{f^{\prime}}{f}}+{\frac{h^{\prime}}{h}}-\frac{2}{r}\right){\cal
H}+{\frac{2{\cal F}}{r}}\right]\phi^{\prime}c_{4}+{\frac{f{\cal F}{\cal
G}}{2r^{2}}}\,,$ (164) $\displaystyle\beta_{2}$ $\displaystyle=$
$\displaystyle\left[\frac{\sqrt{fh}{\cal F}}{r^{2}}\left(2hr\phi^{\prime
2}c_{4}+\frac{r\phi^{\prime}f^{\prime}\sqrt{h}}{2\sqrt{f}}{\cal
H}-\phi^{\prime}\sqrt{fh}{\cal G}\right)-\frac{\phi^{\prime}fh{\cal G}{\cal
H}}{r}\left(\frac{{\cal G}^{\prime}}{{\cal G}}-\frac{{\cal H}^{\prime}}{{\cal
H}}+\frac{f^{\prime}}{2f}-\frac{1}{r}\right)\right]a_{1}-\frac{2}{r}(fh)^{3/2}{\cal
F}{\cal G}{\cal H}\,,\qquad\,\,$ (165) $\displaystyle\beta_{3}$
$\displaystyle=$ $\displaystyle\frac{\sqrt{fh}{\cal
H}}{2}\phi^{\prime}\left(hc_{4}^{\prime}+\frac{1}{2}h^{\prime}c_{4}-\frac{d_{3}}{2}\right)-\frac{\sqrt{fh}}{2}\left(\frac{\cal
H}{r}+{\cal H}^{\prime}\right)\left(2h\phi^{\prime}c_{4}+\frac{\sqrt{fh}{\cal
G}}{2r}+\frac{f^{\prime}\sqrt{h}{\cal H}}{4\sqrt{f}}\right)$ (166)
$\displaystyle+{\frac{\sqrt{fh}{\cal
F}}{4r}\left(2h\phi^{\prime}c_{4}+\frac{3\sqrt{fh}{\cal
G}}{r}+\frac{f^{\prime}\sqrt{h}{\cal H}}{2\sqrt{f}}\right)}\,,$

with

$\displaystyle e_{4}$ $\displaystyle=$
$\displaystyle{\frac{1}{\phi^{\prime}}}c_{4}^{\prime}-{\frac{f^{\prime}}{4fh\phi^{\prime
2}}}\left(\sqrt{fh}{\cal H}\right)^{\prime}-{\frac{\sqrt{f}}{2\phi^{\prime
2}\sqrt{h}r}}{\cal
G}^{\prime}+{\frac{1}{h\phi^{\prime}r^{2}}\left({\frac{\phi^{\prime\prime}}{\phi^{\prime}}}+\frac{1}{2}{\frac{h^{\prime}}{h}}\right)}a_{1}$
(167) $\displaystyle+{\frac{\sqrt{f}}{8\sqrt{h}\phi^{\prime
2}}\left[{\frac{(f^{\prime}r-6f)f^{\prime}}{f^{2}r}}+\frac{h^{\prime}(f^{\prime}r+4f)}{fhr}-{\frac{4f(2\phi^{\prime\prime}h+h^{\prime}\phi^{\prime})}{\phi^{\prime}h^{2}r(f^{\prime}r-2f)}}\right]}{\cal
H}+{\frac{h^{\prime}}{2h\phi^{\prime}}}c_{4}-\frac{f^{\prime}r-2f}{4\sqrt{fh}r\phi^{\prime}}\frac{\partial{\cal
H}}{\partial\phi}$
$\displaystyle+{\frac{f^{\prime}hr-f}{2r^{2}\sqrt{f}{h}^{3/2}\phi^{\prime
2}}}{\cal F}+{\frac{\sqrt{f}}{2r\phi^{\prime
2}{h}^{3/2}}\left[{\frac{f(2\phi^{\prime\prime}h+h^{\prime}\phi^{\prime})}{h\phi^{\prime}(f^{\prime}r-2f)}}+{\frac{2f-f^{\prime}hr}{2fr}}\right]}{\cal
G}+\frac{\sqrt{f}(\rho+P)}{2h^{3/2}\phi^{\prime 2}}\,,$ $\displaystyle d_{3}$
$\displaystyle=$
$\displaystyle-{\frac{1}{r^{2}}\left({\frac{2\phi^{\prime\prime}}{\phi^{\prime}}}+{\frac{h^{\prime}}{h}}\right)}a_{1}+{\frac{f^{3/2}h^{1/2}}{(f^{\prime}r-2f)\phi^{\prime}}\left({\frac{2\phi^{\prime\prime}}{h\phi^{\prime}r}}+{\frac{{f^{\prime}}^{2}}{f^{2}}}-{\frac{f^{\prime}h^{\prime}}{fh}}-{\frac{2f^{\prime}}{fr}}+{\frac{2h^{\prime}}{hr}}+{\frac{h^{\prime}}{h^{2}r}}\right)}{\cal
H}$ (168)
$\displaystyle+\frac{f^{\prime}r-2f}{2r}\sqrt{\frac{h}{f}}\frac{\partial{\cal
H}}{\partial\phi}+{\frac{\sqrt{f}}{\phi^{\prime}\sqrt{h}r^{2}}}{\cal
F}-{\frac{{f}^{3/2}}{\sqrt{h}(f^{\prime}r-2f)\phi^{\prime}}\left({\frac{f^{\prime}}{fr}}+{\frac{2\phi^{\prime\prime}}{\phi^{\prime}r}}+{\frac{h^{\prime}}{hr}}-\frac{2}{r^{2}}\right)}{\cal
G}-\frac{\sqrt{f}(\rho+P)}{\phi^{\prime}\sqrt{h}}\,.$

## References

  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. Lett._ 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ ,  _Phys. Rev. Lett._ 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo),  _Phys. Rev. Lett._ 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt,  _Living Rev. Rel._ 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Nollert (1999) H.-P. Nollert,  _Class. Quant. Grav._ 16, R159 (1999).
  * Riess  _et al._ (1998) A. G. Riess _et al._ (Supernova Search Team),  _Astron. J._ 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ (1999) S. Perlmutter _et al._ (Supernova Cosmology Project),  _Astrophys. J._ 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Spergel  _et al._ (2003) D. N. Spergel _et al._ (WMAP),  _Astrophys. J. Suppl._ 148, 175 (2003), arXiv:astro-ph/0302209 .
  * Tegmark  _et al._ (2004) M. Tegmark _et al._ (SDSS),  _Phys. Rev. D_ 69, 103501 (2004), arXiv:astro-ph/0310723 .
  * Ade  _et al._ (2014) P. A. R. Ade _et al._ (Planck),  _Astron. Astrophys._ 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Copeland  _et al._ (2006) E. J. Copeland, M. Sami, and S. Tsujikawa,  _Int. J. Mod. Phys. D_ 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Silvestri and Trodden (2009) A. Silvestri and M. Trodden,  _Rept. Prog. Phys._ 72, 096901 (2009), arXiv:0904.0024 [astro-ph.CO] .
  * De Felice and Tsujikawa (2010) A. De Felice and S. Tsujikawa,  _Living Rev. Rel._ 13, 3 (2010), arXiv:1002.4928 [gr-qc] .
  * Clifton  _et al._ (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis,  _Phys. Rept._ 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
  * Joyce  _et al._ (2015) A. Joyce, B. Jain, J. Khoury, and M. Trodden,  _Phys. Rept._ 568, 1 (2015), arXiv:1407.0059 [astro-ph.CO] .
  * Koyama (2016) K. Koyama,  _Rept. Prog. Phys._ 79, 046902 (2016), arXiv:1504.04623 [astro-ph.CO] .
  * Heisenberg (2019) L. Heisenberg,  _Phys. Rept._ 796, 1 (2019), arXiv:1807.01725 [gr-qc] .
  * Kase and Tsujikawa (2019a) R. Kase and S. Tsujikawa,  _Int. J. Mod. Phys. D_ 28, 1942005 (2019a), arXiv:1809.08735 [gr-qc] .
  * Hawking (1972a) S. W. Hawking,  _Commun. Math. Phys._ 25, 152 (1972a).
  * Bekenstein (1972) J. D. Bekenstein,  _Phys. Rev. Lett._ 28, 452 (1972).
  * Graham and Jha (2014) A. A. H. Graham and R. Jha,  _Phys. Rev. D_ 89, 084056 (2014), [Erratum: Phys. Rev. D 92, 069901 (2015)], arXiv:1401.8203 [gr-qc] .
  * Hawking (1972b) S. W. Hawking,  _Commun. Math. Phys._ 25, 167 (1972b).
  * Bekenstein (1995) J. D. Bekenstein,  _Phys. Rev. D_ 51, R6608 (1995).
  * Sotiriou and Faraoni (2012) T. P. Sotiriou and V. Faraoni,  _Phys. Rev. Lett._ 108, 081103 (2012), arXiv:1109.6324 [gr-qc] .
  * Faraoni (2017) V. Faraoni,  _Phys. Rev. D_ 95, 124013 (2017), arXiv:1705.07134 [gr-qc] .
  * Hui and Nicolis (2013) L. Hui and A. Nicolis,  _Phys. Rev. Lett._ 110, 241104 (2013), arXiv:1202.1296 [hep-th] .
  * 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 .
  * Chen  _et al._ (2007) C.-M. Chen, D. V. Gal’tsov, and D. G. Orlov,  _Phys. Rev. D_ 75, 084030 (2007), arXiv:hep-th/0701004 .
  * Guo  _et al._ (2008) Z.-K. Guo, N. Ohta, and T. Torii,  _Prog. Theor. Phys._ 120, 581 (2008), arXiv:0806.2481 [gr-qc] .
  * Guo  _et al._ (2009) Z.-K. Guo, N. Ohta, and T. Torii,  _Prog. Theor. Phys._ 121, 253 (2009), arXiv:0811.3068 [gr-qc] .
  * Pani and Cardoso (2009) P. Pani and V. Cardoso,  _Phys. Rev. D_ 79, 084031 (2009), arXiv:0902.1569 [gr-qc] .
  * Sotiriou and Zhou (2014a) T. P. Sotiriou and S.-Y. Zhou,  _Phys. Rev. Lett._ 112, 251102 (2014a), arXiv:1312.3622 [gr-qc] .
  * Sotiriou and Zhou (2014b) T. P. Sotiriou and S.-Y. Zhou,  _Phys. Rev. D_ 90, 124063 (2014b), arXiv:1408.1698 [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] .
  * Maselli  _et al._ (2015) A. Maselli, P. Pani, L. Gualtieri, and V. Ferrari,  _Phys. Rev. D_ 92, 083014 (2015), arXiv:1507.00680 [gr-qc] .
  * Kleihaus  _et al._ (2011) B. Kleihaus, J. Kunz, and E. Radu,  _Phys. Rev. Lett._ 106, 151104 (2011), arXiv:1101.2868 [gr-qc] .
  * Kleihaus  _et al._ (2016) B. Kleihaus, J. Kunz, S. Mojica, and E. Radu,  _Phys. Rev. D_ 93, 044047 (2016), arXiv:1511.05513 [gr-qc] .
  * 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] .
  * Antoniou  _et al._ (2018) G. Antoniou, A. Bakopoulos, and P. Kanti,  _Phys. Rev. Lett._ 120, 131102 (2018), arXiv:1711.03390 [hep-th] .
  * Blázquez-Salcedo  _et al._ (2018) J. L. Blázquez-Salcedo, D. D. Doneva, J. Kunz, and S. S. Yazadjiev,  _Phys. Rev. D_ 98, 084011 (2018), arXiv:1805.05755 [gr-qc] .
  * Minamitsuji and Ikeda (2019) M. Minamitsuji and T. Ikeda,  _Phys. Rev. D_ 99, 044017 (2019), arXiv:1812.03551 [gr-qc] .
  * Silva  _et al._ (2019) H. O. Silva, C. F. B. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein, and E. Berti,  _Phys. Rev. D_ 99, 064011 (2019), arXiv:1812.05590 [gr-qc] .
  * Macedo  _et al._ (2019) C. F. B. Macedo, J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou,  _Phys. Rev. D_ 99, 104041 (2019), arXiv:1903.06784 [gr-qc] .
  * Doneva and Yazadjiev (2022) D. D. Doneva and S. S. Yazadjiev,  _Phys. Rev. D_ 105, L041502 (2022), arXiv:2107.01738 [gr-qc] .
  * Minamitsuji  _et al._ (2022a) M. Minamitsuji, K. Takahashi, and S. Tsujikawa, arXiv:2204.13837 [gr-qc] .
  * Horndeski (1974) G. W. Horndeski,  _Int. J. Theor. Phys._ 10, 363 (1974).
  * Kobayashi  _et al._ (2011) T. Kobayashi, M. Yamaguchi, and J. Yokoyama,  _Prog. Theor. Phys._ 126, 511 (2011), arXiv:1105.5723 [hep-th] .
  * Damour and Esposito-Farese (1993) T. Damour and G. Esposito-Farese,  _Phys. Rev. Lett._ 70, 2220 (1993).
  * Harada (1998) T. Harada,  _Phys. Rev. D_ 57, 4802 (1998), arXiv:gr-qc/9801049 .
  * Novak (1998) J. Novak,  _Phys. Rev. D_ 58, 064019 (1998), arXiv:gr-qc/9806022 .
  * Sotani and Kokkotas (2004) H. Sotani and K. D. Kokkotas,  _Phys. Rev. D_ 70, 084026 (2004), arXiv:gr-qc/0409066 .
  * Cooney  _et al._ (2010) A. Cooney, S. DeDeo, and D. Psaltis,  _Phys. Rev. D_ 82, 064033 (2010), arXiv:0910.5480 [astro-ph.HE] .
  * Arapoglu  _et al._ (2011) A. S. Arapoglu, C. Deliduman, and K. Y. Eksi,  _JCAP_ 07, 020 (2011), arXiv:1003.3179 [gr-qc] .
  * Orellana  _et al._ (2013) M. Orellana, F. Garcia, F. A. Teppa Pannia, and G. E. Romero,  _Gen. Rel. Grav._ 45, 771 (2013), arXiv:1301.5189 [astro-ph.CO] .
  * Aparicio Resco  _et al._ (2016) M. Aparicio Resco, A. de la Cruz-Dombriz, F. J. Llanes Estrada, and V. Zapatero Castrillo,  _Phys. Dark Univ._ 13, 147 (2016), arXiv:1602.03880 [gr-qc] .
  * Kase and Tsujikawa (2019b) R. Kase and S. Tsujikawa,  _JCAP_ 09, 054 (2019b), arXiv:1906.08954 [gr-qc] .
  * Brans and Dicke (1961) C. Brans and R. H. Dicke,  _Phys. Rev._ 124, 925 (1961).
  * Starobinsky (1980) A. A. Starobinsky,  _Phys. Lett. B_ 91, 99 (1980).
  * O’Hanlon (1972) J. O’Hanlon,  _Phys. Rev. Lett._ 29, 137 (1972).
  * Chiba (2003) T. Chiba,  _Phys. Lett. B_ 575, 1 (2003), arXiv:astro-ph/0307338 .
  * Lehébel  _et al._ (2017) A. Lehébel, E. Babichev, and C. Charmousis,  _JCAP_ 07, 037 (2017), arXiv:1706.04989 [gr-qc] .
  * Babichev and Charmousis (2014) E. Babichev and C. Charmousis,  _JHEP_ 08, 106 (2014), arXiv:1312.3204 [gr-qc] .
  * Babichev  _et al._ (2016) E. Babichev, K. Koyama, D. Langlois, R. Saito, and J. Sakstein,  _Class. Quant. Grav._ 33, 235014 (2016), arXiv:1606.06627 [gr-qc] .
  * Sakstein  _et al._ (2017) J. Sakstein, E. Babichev, K. Koyama, D. Langlois, and R. Saito,  _Phys. Rev. D_ 95, 064013 (2017), arXiv:1612.04263 [gr-qc] .
  * Kobayashi and Hiramatsu (2018) T. Kobayashi and T. Hiramatsu,  _Phys. Rev. D_ 97, 104012 (2018), arXiv:1803.10510 [gr-qc] .
  * Ogawa  _et al._ (2020) H. Ogawa, T. Kobayashi, and K. Koyama,  _Phys. Rev. D_ 101, 024026 (2020), arXiv:1911.01669 [gr-qc] .
  * Ikeda  _et al._ (2021) T. Ikeda, A. Iyonaga, and T. Kobayashi,  _Phys. Rev. D_ 104, 104009 (2021), arXiv:2107.13804 [gr-qc] .
  * Cisterna  _et al._ (2015) A. Cisterna, T. Delsate, and M. Rinaldi,  _Phys. Rev. D_ 92, 044050 (2015), arXiv:1504.05189 [gr-qc] .
  * Kase and Tsujikawa (2021) R. Kase and S. Tsujikawa,  _JCAP_ 01, 008 (2021), arXiv:2008.13350 [gr-qc] .
  * Kase and Tsujikawa (2022) R. Kase and S. Tsujikawa,  _Phys. Rev. D_ 105, 024059 (2022), arXiv:2110.12728 [gr-qc] .
  * Maselli  _et al._ (2016) A. Maselli, H. O. Silva, M. Minamitsuji, and E. Berti,  _Phys. Rev. D_ 93, 124056 (2016), arXiv:1603.04876 [gr-qc] .
  * Minamitsuji  _et al._ (2022b) M. Minamitsuji, K. Takahashi, and S. Tsujikawa,  _Phys. Rev. D_ 105, 104001 (2022b), arXiv:2201.09687 [gr-qc] .
  * Babichev  _et al._ (2017) E. Babichev, C. Charmousis, and A. Lehébel,  _JCAP_ 04, 027 (2017), arXiv:1702.01938 [gr-qc] .
  * Creminelli  _et al._ (2020) P. Creminelli, N. Loayza, F. Serra, E. Trincherini, and L. G. Trombetta,  _JHEP_ 08, 045 (2020), arXiv:2004.02893 [hep-th] .
  * Kobayashi  _et al._ (2012) T. Kobayashi, H. Motohashi, and T. Suyama,  _Phys. Rev. D_ 85, 084025 (2012), [Erratum: Phys. Rev. D 96, 109903 (2017)], arXiv:1202.4893 [gr-qc] .
  * Kobayashi  _et al._ (2014) T. Kobayashi, H. Motohashi, and T. Suyama,  _Phys. Rev. D_ 89, 084042 (2014), arXiv:1402.6740 [gr-qc] .
  * Kase  _et al._ (2020) R. Kase, R. Kimura, S. Sato, and S. Tsujikawa,  _Phys. Rev. D_ 102, 084037 (2020), arXiv:2007.09864 [gr-qc] .
  * Pani  _et al._ (2011) P. Pani, E. Berti, V. Cardoso, and J. Read,  _Phys. Rev. D_ 84, 104035 (2011), arXiv:1109.0928 [gr-qc] .
  * Kleihaus  _et al._ (2014) B. Kleihaus, J. Kunz, and S. Mojica,  _Phys. Rev. D_ 90, 061501 (2014), arXiv:1407.6884 [gr-qc] .
  * Doneva and Yazadjiev (2018b) D. D. Doneva and S. S. Yazadjiev,  _JCAP_ 04, 011 (2018b), arXiv:1712.03715 [gr-qc] .
  * Blázquez-Salcedo  _et al._ (2016) J. L. Blázquez-Salcedo, L. M. González-Romero, J. Kunz, S. Mojica, and F. Navarro-Lérida,  _Phys. Rev. D_ 93, 024052 (2016), arXiv:1511.03960 [gr-qc] .
  * Olmo  _et al._ (2020) G. J. Olmo, D. Rubiera-Garcia, and A. Wojnar,  _Phys. Rept._ 876, 1 (2020), arXiv:1912.05202 [gr-qc] .
  * Glavan and Lin (2020) D. Glavan and C. Lin,  _Phys. Rev. Lett._ 124, 081301 (2020), arXiv:1905.03601 [gr-qc] .
  * Lu and Pang (2020) H. Lu and Y. Pang,  _Phys. Lett. B_ 809, 135717 (2020), arXiv:2003.11552 [gr-qc] .
  * Kobayashi (2020) T. Kobayashi,  _JCAP_ 07, 013 (2020), arXiv:2003.12771 [gr-qc] .
  * Fernandes  _et al._ (2020) P. G. S. Fernandes, P. Carrilho, T. Clifton, and D. J. Mulryne,  _Phys. Rev. D_ 102, 024025 (2020), arXiv:2004.08362 [gr-qc] .
  * Hennigar  _et al._ (2020) R. A. Hennigar, D. Kubizňák, R. B. Mann, and C. Pollack,  _JHEP_ 07, 027 (2020), arXiv:2004.09472 [gr-qc] .
  * Doneva and Yazadjiev (2021) D. D. Doneva and S. S. Yazadjiev,  _JCAP_ 05, 024 (2021), arXiv:2003.10284 [gr-qc] .
  * Nojiri and Odintsov (2005) S. Nojiri and S. D. Odintsov,  _Phys. Lett. B_ 631, 1 (2005), arXiv:hep-th/0508049 .
  * De Felice and Hindmarsh (2007) A. De Felice and M. Hindmarsh,  _JCAP_ 06, 028 (2007), arXiv:0705.3375 [astro-ph] .
  * Li  _et al._ (2007) B. Li, J. D. Barrow, and D. F. Mota,  _Phys. Rev. D_ 76, 044027 (2007), arXiv:0705.3795 [gr-qc] .
  * De Felice and Tsujikawa (2009a) A. De Felice and S. Tsujikawa,  _Phys. Lett. B_ 675, 1 (2009a), arXiv:0810.5712 [hep-th] .
  * De Felice and Tsujikawa (2009b) A. De Felice and S. Tsujikawa,  _Phys. Rev. D_ 80, 063516 (2009b), arXiv:0907.1830 [hep-th] .
  * Deffayet  _et al._ (2011) C. Deffayet, X. Gao, D. A. Steer, and G. Zahariade,  _Phys. Rev. D_ 84, 064039 (2011), arXiv:1103.3260 [hep-th] .
  * Charmousis  _et al._ (2012) C. Charmousis, E. J. Copeland, A. Padilla, and P. M. Saffin,  _Phys. Rev. Lett._ 108, 051101 (2012), arXiv:1106.2000 [hep-th] .
  * Schutz and Sorkin (1977) B. F. Schutz and R. Sorkin,  _Annals Phys._ 107, 1 (1977).
  * Brown (1993) J. D. Brown,  _Class. Quant. Grav._ 10, 1579 (1993), arXiv:gr-qc/9304026 .
  * De Felice  _et al._ (2010) A. De Felice, J.-M. Gerard, and T. Suyama,  _Phys. Rev. D_ 81, 063527 (2010), arXiv:0908.3439 [gr-qc] .
  * Regge and Wheeler (1957) T. Regge and J. A. Wheeler,  _Phys. Rev._ 108, 1063 (1957).
  * Zerilli (1970) F. J. Zerilli,  _Phys. Rev. Lett._ 24, 737 (1970).
  * De Felice  _et al._ (2011) A. De Felice, T. Suyama, and T. Tanaka,  _Phys. Rev. D_ 83, 104035 (2011), arXiv:1102.1521 [gr-qc] .
  * Motohashi and Suyama (2011) H. Motohashi and T. Suyama,  _Phys. Rev. D_ 84, 084041 (2011), arXiv:1107.3705 [gr-qc] .
  * Haensel and Potekhin (2004) P. Haensel and A. Y. Potekhin,  _Astron. Astrophys._ 428, 191 (2004), arXiv:astro-ph/0408324 .
  * Saffer and Yagi (2021) A. Saffer and K. Yagi,  _Phys. Rev. D_ 104, 124052 (2021), arXiv:2110.02997 [gr-qc] .