# Two lower bounds for maximum matching

Fedor<EMAIL_ADDRESS>HSE University, Moscow, Russian
Federation

###### Abstract

In this paper, we prove two lower bounds for the maximum matching size in an
arbitrary undirected graph. Despite their simplicity, these results are not
widely known. This article aims to bring pleasure to the reader by giving
short combinatorial proofs of these easy-to-state estimates.

##  1 Main results

We focus on simple undirected graphs with $n$ vertices, $m$ edges, and having
the maximum degree $d$. We investigate the following question: how small can
its maximum matching be? Surprisingly, the answer to this question is hard to
find in the literature. In the paper by P. Erdős and T. Gallai [2], the
authors present a lower bound (Theorem 4.1); however, their proof seems quite
complicated. Another related work is the article by V. Chvátal and D. Hanson
[1] providing a sharp estimate for the maximum number of edges in a graph not
containing a matching of a specific size. Their estimate involves three cases
and is not particularly easy to state. In our work, we present two easy-to-
state estimates and give short combinatorial proofs. Note that these estimates
follow immediately from the much harder Vizing’s theorem [3].

Let $G=(V,E)$ be a simple undirected graph (i.e. without loops and multiple
edges). Denote $n:=|V|$, $m:=|E|$, and $d:=\max_{v\in V}\deg(v)$. A subset of
edges is called a matching if no two edges have a common vertex. Let $k$ be
the maximum matching size. We prove two following lower bounds:

  1. 1.

$k\geq m/n$ (Theorem 1)

  2. 2.

$k\geq 2/3\cdot m/d$ (Theorem 2)

Note that these results strengthen the trivial lower bounds of $m/(2n)$ and
$m/(2d)$ by a constant factor. Also, the constant $2/3$ in the second estimate
is sharp: consider a graph consisting of $r$ disjoint triangles; its maximum
matching has $r=2/3\cdot 3r/2$ edges.

Let $M$ be the set of vertices covered by the maximum matching $\pi_{max}$.
For $u\in M$, denote by $d_{in}(u)$ and $d_{out}(u)$ the number of the
incident edges in $G$ that connect $u$ with vertices inside and outside $M$,
respectively. First, let us present the following observations.

######  Lemma 1.

For each $(u,v)\in\pi_{max}$ we have

  1. 1.

$d_{in}(u)+d_{out}(u)\leq d$,

  2. 2.

$d_{in}(u)\leq 2k$,

  3. 3.

$d_{out}(u)\leq n-2k$,

  4. 4.

either $d_{out}(u)=d_{out}(v)=1$ or at least one of
$\\{d_{out}(u),d_{out}(v)\\}$ is zero.

###### Proof.

The first claim is trivial because $d_{in}(u)+d_{out}(u)=\deg(u)$. The second
and third claims are also clear because the graph has no multiple edges and
thus $d_{in}(u)\leq|M|=2k$ and $d_{out}(u)\leq|V\setminus M|=n-2k$. Let us
prove the last claim.

Consider the sets of vertices $N(u),N(v)\subset V\setminus M$ which are
adjacent to $u$ and $v$, respectively. We have $d_{out}(u)=|N(u)|$ and
$d_{out}(v)=|N(v)|$. If at least one of them is empty, then all is done.
Otherwise, assume that $|N(v)|\geq 2$. Take $u^{\prime}\in N(u)$ and
$v^{\prime}\in N(v)$ such that $u^{\prime}\neq v^{\prime}$. In this case, we
can improve the matching $\pi_{max}$ by removing the edge $(u,v)$ and adding
$(u,u^{\prime})$ and $(v,v^{\prime})$ instead. This contradicts the maximality
of $\pi_{max}$. ∎

######  Lemma 2.

We have
$m=\sum_{(u,v)\in\pi_{max}}\left[\frac{d_{in}(u)+d_{in}(v)}{2}+d_{out}(u)+d_{out}(v)\right]$.

###### Proof.

Note that all edges of $G$ are either inside $M$ or go from $M$ to $V\setminus
M$; there are no edges completely outside $M$ because one could add them to
$\pi_{max}$. Let $m_{in}$ and $m_{out}$ be the number of those edges. By a
simple counting argument,

$\displaystyle
2m_{in}=\sum_{(u,v)\in\pi_{max}}d_{in}(u)+d_{in}(v)\qquad\text{and}\qquad
m_{out}=\sum_{(u,v)\in\pi_{max}}d_{out}(u)+d_{out}(v),$

from which the statement of the lemma follows. ∎

Now we are ready to prove two lower bounds for the size of the maximum
matching.

######  Theorem 1.

For each graph G we have $k\geq m/n$.

###### Proof.

We have the following two cases:

  1. 1.

$n-2k\leq 1$. In this case, $k\geq(n-1)/2\geq m/n$, where the last inequality
follows from the estimate $m\leq n(n-1)/2$.

  2. 2.

$n-2k\geq 2$. In this case, for $(u,v)\in\pi_{max}$, from the last two claims
of Lemma 1 it follows that $d_{out}(u)+d_{out}(v)\leq n-2k$. Substituting this
inequality into Lemma 2 and applying the second claim of Lemma 1, we get

$\displaystyle
m\leq\sum_{(u,v)\in\pi_{max}}\left[\frac{4k}{2}+n-2k\right]=kn,$

proving the desired estimate.

∎

######  Theorem 2.

For each graph G with the maximum degree $d$ we have $k\geq 2/3\cdot m/d$.

###### Proof.

If $d=1$, then $G$ is a matching and thus $k=m\geq\frac{2m}{3}$. Hereafter we
assume $d\geq 2$. For each $(u,v)\in\pi_{max}$ we will show that
$s(u,v):=\frac{d_{in}(u)+d_{in}(v)}{2}+d_{out}(u)+d_{out}(v)\leq\frac{3}{2}d$.
From the last claim of Lemma 1 we have the following two cases:

  1. 1.

$d_{out}(u)=d_{out}(v)=1$. From the first claim of Lemma 1 it follows that
$d_{in}(u),d_{in}(v)\leq d-1$. Thus,

$\displaystyle s(u,v)\leq d+1\leq\frac{3d}{2}.$

  2. 2.

$d_{out}(v)=0$ (wlog). Again, by the first claim of Lemma 1,

$\displaystyle
s(u,v)=\left(\frac{d_{in}(u)}{2}+d_{out}(u)\right)+\frac{d_{in}(v)}{2}\leq(d_{in}(u)+d_{out}(u))+d_{in}(v)/2\leq
d+d/2=\frac{3d}{2}.$

Finally, plugging this inequality into Lemma 2 we have

$\displaystyle m=\sum_{(u,v)\in\pi_{max}}s(u,v)\leq\frac{3kd}{2}\implies k\geq
2/3\cdot m/d.$

∎

##  2 Acknowledgements

This paper was prepared within the framework of the HSE University Basic
Research Program. We thank Benjamin Sudakov, Mikhail Vyalyi, Vladimir
Podolskii, and Stasys Jukna for useful discussions on this topic. We also
thank Martin Milanič for providing the bibliographical references.

## References

  * [1] Chvátal, V., Hanson, D. Degrees and matchings. Journal of Combinatorial Theory, Series B, Volume 20, Issue 2, 1976, Pages 128-138, ISSN 0095-8956, https://doi.org/10.1016/0095-8956(76)90004-6.
  * [2] Erdős, P., Gallai, T. On maximal paths and circuits of graphs. Acta Mathematica Academiae Scientiarum Hungaricae 10, 337–356 (1959), https://doi.org/10.1007/BF02024498.
  * [3] Vizing, V. On an estimate of the chromatic class of a p-graph. Diskret. Analiz. 3 (1964) 25–30 (in Russian).


# Neutron Stars in Scalar-Tensor Theories:
Analytic Scalar Charges and Universal Relations

Kent Yagi  Department of Physics, University of Virginia, P.O. Box 400714,
Charlottesville, VA 22904-4714, USA  Michael Stepniczka  Department of
Physics, University of Virginia, P.O. Box 400714, Charlottesville, VA
22904-4714, USA

###### Abstract

Neutron stars are ideal astrophysical sources to probe general relativity due
to their large compactnesses and strong gravitational fields. For example,
binary pulsar and gravitational wave observations have placed stringent bounds
on certain scalar-tensor theories in which a massless scalar field is coupled
to the metric through matter. A remarkable phenomenon of neutron stars in such
scalar-tensor theories is spontaneous scalarization, where a normalized scalar
charge remains order unity even if the matter-scalar coupling vanishes
asymptotically far from the neutron star. While most works on scalarization of
neutron stars focus on numerical analysis, this paper aims to derive accurate
scalar charges _analytically_. To achieve this, we consider a simple energy
density profile of the Tolman VII form and work in a weak-field expansion. We
solve the modified Tolman-Oppenheimer-Volkoff equations order by order and
apply Padé resummation to account for higher order effects. We find that our
analytic scalar charges in terms of the stellar compactness beautifully model
those computed numerically. We also find a quasi-universal relation between
the scalar charge and stellar binding energy that is insensitive to the
underlying equations of state. Comparison of analytic scalar charges for
Tolman VII and constant density stars mathematically supports this quasi-
universal relation. The analytic results found here provide physically
motivated, ready-to-use accurate expressions for scalar charges.

##  I Introduction

Neutron stars (NSs) are ideal compact astrophysical objects to probe
fundamental physics. Due to their high central density that exceeds the
saturation density of nuclear matter by several-fold, NSs can efficiently test
nuclear physics. For example, recent observations of $2M_{\odot}$ pulsars
Cromartie  _et al._ (2019), x-ray emissions from hot spots on a rotating NS
surface by NICER Riley  _et al._ (2019); Miller  _et al._ (2019), and
gravitational waves from colliding NSs by LIGO and Virgo Abbott  _et al._
(2017) have constrained properties of nuclear/quark matter and certain nuclear
parameters (see e.g. Abbott  _et al._ (2018); Annala  _et al._ (2018); Raithel
_et al._ (2018); Lim and Holt (2018); Bauswein  _et al._ (2017); De  _et al._
(2018); Most  _et al._ (2018); Annala  _et al._ (2019); Malik  _et al._
(2018); Carson  _et al._ (2019a, b); Raithel and Ozel (2019); Chatziioannou
(2020); Raaijmakers  _et al._ (2020); Zimmerman  _et al._ (2020); Jiang  _et
al._ (2020); Dietrich  _et al._ (2020)). Due to their large compactnesses, NSs
are also perfect sources to probe strong-field gravity. Indeed, observations
of binary pulsars Freire  _et al._ (2012); Yagi  _et al._ (2014a, b); Berti
_et al._ (2015); Shao  _et al._ (2017); Archibald  _et al._ (2018); Seymour
and Yagi (2020); Anderson  _et al._ (2019); Gupta  _et al._ (2021) and the
binary neutron star merger GW170817 Abbott  _et al._ (2019); Baker  _et al._
(2017); Sakstein and Jain (2017); Ezquiaga and Zumalacárregui (2017); Zhao
_et al._ (2019); Silva  _et al._ (2020a) have constrained various
modifications to GR.

Some of the most well-studied modified theories of gravity are scalar-tensor
theories in which one introduces scalar fields (either minimally or non-
minimally coupled to the metric) to the action. A simple scalar-tensor theory
proposed by Damour and Esposito-Farèse has two theoretical parameters
$(\alpha_{0},\beta_{0})$ and allows NSs to scalarize _spontaneously_ Damour
and Esposito-Farese (1992, 1993, 1996). Namely, when $\beta_{0}$ is
sufficiently negative, NSs can have scalar charges of order unity even when
$\alpha_{0}\ll 1$ (or equivalently, when the value of the scalar field at
infinity is small). Such a phenomenon can be understood as NSs undergoing a
tachyonic instability, where the effective mass of the scalar field becomes
imaginary Ramazanoğlu (2017). A catalog of NS scalar charges is provided in
Anderson and Yunes (2019) while a surrogate model has recently been
constructed in Zhao  _et al._ (2019) based on _numerical_ calculations. The
parameters $(\alpha_{0},\beta_{0})$ have been constrained with solar system
experiments, various binary pulsars Freire  _et al._ (2012); Shao  _et al._
(2017); Archibald  _et al._ (2018); Anderson  _et al._ (2019) and the binary
NS merger event GW170817 Zhao  _et al._ (2019) to be $|\alpha_{0}|\lesssim
3\times 10^{-4}$ and $\beta_{0}\gtrsim-4.4$. Future forecasts on probing this
scalar-tensor theory with black-hole/NS binaries have been made with pulsar
Berti  _et al._ (2015) and gravitational-wave Shao  _et al._ (2017); Carson
_et al._ (2020); Zhao  _et al._ (2019) observations.

The goal of this paper is to provide accurate expressions for scalar charges
of NSs in the scalar-tensor theory in Damour and Esposito-Farese (1993) by
solving the field equations _analytically_. Such expressions are complementary
to e.g. the surrogate model Zhao  _et al._ (2019) mentioned earlier. To
achieve this, instead of using realistic equations of state (EoSs) given in
tables, we adopt the Tolman VII model Tolman (1939) that approximates the
energy density profile inside a NS to a quadratic function in a radial
coordinate. Such a simple profile allows one to solve the Tolman-Oppenheimer-
Volkoff (TOV) equations analytically in GR. Unfortunately, due to the
complication of the modified TOV equations in scalar-tensor theories including
the coupling with the scalar field, it is challenging to solve them
analytically. To overcome this, we work in a weak-field approximation and
solve the field equations order by order in the expansion111Similar
calculations were carried out in quadratic gravity Yagi  _et al._ (2016) and
Einstein-Æther theory Gupta  _et al._ (2021) to respectively find analytic
scalar charges and sensitivities, which were then used in tandem with binary
pulsar observations to place bounds on these theories.. To take into account
higher order contributions in the expansion, we apply Padé resummation. We
consider both spontaneous scalarization mentioned earlier and “perturbative”
scalarization where the scalar charge is proportional to $\alpha_{0}$.

Figure 1 summarizes our findings, which compares the analytic scalar charges
with the Tolman VII model in terms of the stellar compactness (the ratio
between the mass and radius) against those with realistic EoSs (WFF1 Wiringa
_et al._ (1988), SLy Douchin and Haensel (2001), MPA1 Müther  _et al._ (1987),
MS1 Mueller and Serot (1996)) computed numerically. Notice that the former can
accurately model the latter, especially for the SLy EoS, for both the
spontaneous and perturbative scalarization cases. We also find that when we
plot the scalar charges as a function of the stellar binding energy (the
difference between the gravitational and baryonic mass), the relation between
these quantities becomes quasi-universal and is insensitive to EoSs with an
error of $\sim 1\%$ (see Pappas  _et al._ (2019) for other universal relations
involving scalar charges in scalar-tensor theories). We estimate analytically
the amount of the quasi-universality by comparing analytic expressions for
scalar charges with the Tolman VII model and with constant density stars and
find that it is quasi-universal with an error of 0.3%. The technique developed
here should easily be applicable to other theories beyond GR to compute the
stellar charges (or sensitivities).

Figure 1:  Comparison of analytic scalar charges with the Tolman VII model
against numerical scalar charges with realistic EoSs as a function of the
compactness. For the former, we show the results for both perturbative (solid)
and spontaneous (dashed) scalarization. We present the results with four
different choices of $\beta_{0}$. $\alpha_{0}$ is fixed to
$\alpha_{0}=10^{-5}$ except for spontaneous scalarization for the Tolman case
which assumes $\alpha_{0}=0$. Observe that the analytic Tolman results
accurately model numerical results, especially with the SLy EoS.

The rest of the paper is organized as follows. In Sec. II, we review the
scalar-tensor theory proposed in Damour and Esposito-Farese (1993), explain
the modified TOV equations, present the relation between mass and radius, and
review spontaneous scalarization following Horbatsch and Burgess (2011). In
Sec. III, we present the formalism for computing scalar charges analytically
in the Tolman VII model by combining weak-field expansions and Padé
resummation. We compare such analytic results against numerical ones in Sec.
IV and also present the quasi-universal relation between the scalar charge and
binding energy. We conclude in Sec. V and give several different directions
for future work. In Appendix A, we repeat the analytic calculation for
constant density stars and derive scalar charges similar to the Tolman VII
case. We use the geometric units of $c=G=1$. The main expressions for scalar
charges are summarized in a Supplemental Mathematica notebook mat .

##  II Scalar-tensor Theories and Neutron Stars

We begin by reviewing scalar-tensor theories and spontaneous scalarization of
NSs.

###  II.1 Theory

Presented here are the action and field equations for scalar-tensor theories
in the Einstein frame. The former is given by Damour and Esposito-Farese
(1992, 1993, 1996)

$S=\int
d^{4}x\frac{\sqrt{-g}}{16\pi}\left(R-2g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi\right)+S_{\mathrm{mat}}\left[\psi,A^{2}(\varphi)g_{\mu\nu}\right]\,,$
(1)

where $R$ is the Ricci scalar for the metric $g_{\mu\nu}$ in this frame, $g$
is its determinant, $\varphi$ is the scalar field and $\psi$ is the matter
field. $A$ is the conformal factor that relates $g_{\mu\nu}$ and the metric
$\tilde{g}_{\mu\nu}$ in the (physical) Jordan frame as
$\tilde{g}_{\mu\nu}=A(\varphi)g_{\mu\nu}$. Notice that we have set $G=1$ for
the bare gravitational constant $G$. Varying the above action with respect to
$g_{\mu\nu}$ and $\varphi$, the field equations are given by Damour and
Esposito-Farese (1993, 1996)

$\displaystyle R_{\mu\nu}$ $\displaystyle=$ $\displaystyle
2\partial_{\mu}\varphi\partial_{\nu}\varphi+8\pi\left(T_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T\right)\,,$
(2) $\displaystyle\square\varphi$ $\displaystyle=$
$\displaystyle-4\pi\alpha(\varphi)T\,,$ (3)

where $T_{\mu\nu}$ is the matter stress-energy tensor and

$\alpha(\varphi)\equiv\frac{\partial\ln A(\varphi)}{\partial\varphi}\,.$ (4)

In this paper, we consider an example scalar-tensor theory first considered by
Damour and Esposito-Farèse Damour and Esposito-Farese (1992, 1993, 1996) with

$A(\varphi)=\exp\left(\frac{\beta_{0}}{2}\varphi^{2}\right)\,.$ (5)

Another parameter of the theory is
$\alpha_{0}\equiv\alpha(\varphi_{0})=\beta_{0}\varphi_{0}$ with $\varphi_{0}$
representing the scalar field at infinity.

###  II.2 Neutron Stars

To construct a NS solution, we begin with the metric ansatz given by

$ds^{2}=-e^{\nu(r)}dt^{2}+e^{\lambda(r)}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\,.$ (6)

For the matter sector, we consider a perfect fluid whose stress-energy tensor
in the Einstein frame is given by Pani and Berti (2014)

$T_{\mu\nu}=A^{4}(\varphi)\left[(\widetilde{\rho}+\widetilde{P})u_{\mu}u_{\nu}+g_{\mu\nu}\widetilde{P}\right]\,,$
(7)

where $u^{\mu}$ is the 4-velocity of the fluid. Notice that $\widetilde{\rho}$
and $\widetilde{P}$ are the energy density and pressure in the Jordan frame
that directly enters in the EoS $\widetilde{P}(\widetilde{\rho})$. Plugging
these into the field equations in Eqs. (2) and (3), one finds Damour and
Esposito-Farese (1993, 1996)222There is also an equation for $\nu$ that we do
not present here since it is unnecessary for deriving scalar charges.

$\displaystyle m^{\prime}$ $\displaystyle=$ $\displaystyle 4\pi
r^{2}A^{4}(\varphi)\widetilde{\rho}+\frac{1}{2}r(r-2m)\varphi^{\prime 2}\,,$
(8) $\displaystyle\widetilde{P}^{\prime}$ $\displaystyle=$
$\displaystyle-(\widetilde{\rho}+\widetilde{P})\left[\frac{m+4\pi
A^{4}\widetilde{P}r^{3}}{r(r-2m)}+\frac{1}{2}r\varphi^{\prime
2}+\alpha(\varphi)\varphi^{\prime}\right]\,,$
$\displaystyle\varphi^{\prime\prime}$ $\displaystyle=$ $\displaystyle
4\pi\frac{r}{r-2m}A^{4}(\varphi)\left[\alpha(\varphi)(\widetilde{\rho}-3\widetilde{P})+r\varphi^{\prime}(\widetilde{\rho}-\widetilde{P})\right]$
(10) $\displaystyle-\frac{2(r-m)}{r(r-2m)}\varphi^{\prime}\,,$

where a prime represents an $r$ derivative and

$e^{-\lambda(r)}=1-\frac{2m(r)}{r}\,.$ (11)

The asymptotic behavior of the scalar field at infinity is given by

$\varphi=\varphi_{0}-\alpha_{A}\frac{M_{A}}{r}+\mathcal{O}\left(\frac{M_{A}^{2}}{r^{2}}\right)\,,$
(12)

where

$\alpha_{A}=\frac{\partial\ln M_{A}}{\partial\varphi_{0}}$ (13)

is the scalar charge while $M_{A}$ is the total gravitational mass that can be
read off from the asymptotic behavior of $m(r)$ at infinity as

$m(r)=M_{A}+\mathcal{O}\left(\frac{M_{A}}{r}\right)\,.$ (14)

We can construct NS solutions as follows. First, we impose the initial
conditions of

$m(0)=0\,,\quad\widetilde{P}(0)=\widetilde{P}_{c}\,,\quad\varphi(0)=\varphi_{c}\,,\quad\varphi^{\prime}(0)=0\,,$
(15)

where $\widetilde{P}_{c}$ is the central pressure while $\varphi_{c}$ is the
central value for the scalar field. Under these initial conditions together
with a choice of $(\alpha_{0},\beta_{0})$ and an EoS, we solve Eqs. (8)–(10)
in the interior region numerically and the stellar radius $R$ is determined by
$\widetilde{P}(R)=0$. We then solve the equations in the exterior region with
$\widetilde{P}=\widetilde{\rho}=0$. Finally, we read off $M_{A}$,
$\varphi_{0}$ and $\alpha_{A}$ by comparing the asymptotic behavior of the
numerically-computed $\varphi$ and $m$ with Eqs. (12) and (14). Figure 2
presents the mass-radius relation of NSs in GR and the scalar-tensor theory
for four representative EoSs (WFF1 Wiringa  _et al._ (1988), SLy Douchin and
Haensel (2001), MPA1 Müther  _et al._ (1987), MS1 Mueller and Serot (1996))
with different stiffness333These numerical results are computed with a
Mathematica notebook developed in Pani and Berti (2014).. “Humps” in the
relation for the scalar-tensor theory correspond to NSs with spontaneous
scalarization.

Figure 2:  Mass-radius relation for NSs with four representative EoSs in GR
($\alpha_{0}=\beta_{0}=0$) and the scalar-tensor theory ($\alpha_{0}=10^{-5}$
and $\beta_{0}=-4.4$, $-4.5$). The latter deviates from the former due to
spontaneous scalarization. Two cyan dots correspond to the NSs of
$1.8M_{\odot}$ with SLy EoS whose energy density profiles are shown in Fig. 4.

###  II.3 Spontaneous Scalarization

An analytic attempt of computing scalar charges with spontaneous scalarization
was taken in Horbatsch and Burgess (2011). First, when $\alpha$ at the center
of a star, $\alpha_{c}$, is small, one can expand both $\alpha_{0}$ and
$\alpha_{A}$ in terms of $\alpha_{c}$ as444The scalar field (and its
derivative) enter in even powers in Eqs. (8) and (II.2), and in odd powers in
Eq. (10). This means that $\alpha_{0}$ enters in odd powers in $\alpha_{c}$
and $\alpha_{A}$. Inverting the former order by order (and substituting it to
the latter), one finds that $\alpha_{0}$ and $\alpha_{A}$ enter in odd powers
in $\alpha_{c}$.

$\displaystyle\alpha_{0}$ $\displaystyle=$ $\displaystyle
d_{1}\alpha_{c}+d_{2}\alpha_{c}^{3}+\mathcal{O}(\alpha_{c}^{5})\,,$ (16)
$\displaystyle\alpha_{A}$ $\displaystyle=$ $\displaystyle
e_{1}\alpha_{c}+e_{2}\alpha_{c}^{3}+\mathcal{O}(\alpha_{c}^{5})\,.$ (17)

Next, one can solve Eq. (16) for $\alpha_{c}$ as

$\alpha_{c}=\omega C_{+}+\bar{\omega}C_{-}\,,$ (18)

with

$\omega=1,\quad-e^{\pm i\pi/3}\,,$ (19)

and

$C_{\pm}=\left(\frac{\alpha_{0}}{2d_{2}}\pm\sqrt{D}\right)^{1/3}\,,\quad
D=\left(\frac{\alpha_{0}}{2d_{2}}\right)^{2}+\left(\frac{d_{1}}{3d_{2}}\right)^{3}\,.$
(20)

From Eqs. (17) and (18), the scalar charge is given by

$\alpha_{A}=(\omega C_{+}+\bar{\omega}C_{-})e_{1}+(\omega
C_{+}+\bar{\omega}C_{-})^{3}e_{2}+\mathcal{O}(\alpha_{c}^{5})\,.$ (21)

For example, when $\alpha_{0}=0$555Given current stringent bounds on
$\alpha_{0}$, $\alpha_{0}\approx 0$ is a valid approximation when discussing
spontaneous scalarization., $\alpha_{c}$ can take non-vanishing values as
$\alpha_{c}=\pm\sqrt{-d_{1}/d_{2}}$ and

$\alpha_{A}=\pm\left[\left(-\frac{d_{1}}{d_{2}}\right)^{1/2}e_{1}+\left(-\frac{d_{1}}{d_{2}}\right)^{3/2}e_{2}\right]\,.$
(22)

Thus, although $\alpha_{0}=0$, the scalar charge becomes non-vanishing when
$d_{1}/d_{2}<0$.

For constant density stars, $d_{1}$ and $e_{1}$ are given in a closed analytic
form as Horbatsch and Burgess (2011)

$\displaystyle d_{1}$ $\displaystyle=$
$\displaystyle\operatorname{HeunG}\left(\widetilde{a},\widetilde{q};\widetilde{\alpha}_{-},\widetilde{\alpha}_{+},\frac{3}{2},\frac{3}{2};Z\right)$
$\displaystyle-\left(\widetilde{Z}\log\widetilde{Z}\right)\,\operatorname{HeunG}^{\prime}\left(\widetilde{a},\widetilde{q};\widetilde{\alpha}_{-},\widetilde{\alpha}_{+},\frac{3}{2},\frac{3}{2};Z\right)\,,$
$\displaystyle e_{1}$ $\displaystyle=$
$\displaystyle\beta_{0}\widetilde{Z}\operatorname{HeunG}^{\prime}\left(\widetilde{a},\widetilde{q};\widetilde{\alpha}_{-},\widetilde{\alpha}_{+},\frac{3}{2},\frac{3}{2};Z\right)\,,$
(24)

with

$\displaystyle\widetilde{a}$ $\displaystyle=$
$\displaystyle-\frac{1}{1+3\widetilde{p}_{c}}\,,$ (25)
$\displaystyle\widetilde{q}$ $\displaystyle=$
$\displaystyle\frac{3\beta_{0}}{2}\left(\frac{3\widetilde{p}_{c}-1}{3\widetilde{p}_{c}+1}\right)\,,$
(26) $\displaystyle\widetilde{\alpha}_{\pm}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\left(1\pm\sqrt{1-\frac{8\beta_{0}}{3}}\right)\,,$
(27) $\displaystyle Z$ $\displaystyle=$
$\displaystyle\frac{\widetilde{p}_{c}}{1+3\widetilde{p}_{c}}\,,$ (28)
$\displaystyle\widetilde{Z}$ $\displaystyle=$
$\displaystyle\frac{1+\widetilde{p}_{c}}{1+3\widetilde{p}_{c}}\,,$ (29)

and $\widetilde{p}_{c}=\widetilde{P}_{c}/\widetilde{\rho}_{c}$ with
$\widetilde{P}_{c}$ and $\widetilde{\rho}_{c}$ being the central pressure and
energy density. HeunG is the general Heun function (a generalization of the
hypergeometric function) while its prime refers to a derivative with respect
to $Z$.

Figure 3 shows $d_{1}$ for constant density stars as a function of the stellar
compactness $C_{A}$ defined by

$C_{A}=\frac{M_{A}}{R}\,.$ (30)

One can check numerically that spontaneous scalarization occurs when $d_{1}<0$
Horbatsch and Burgess (2011). This means that $d_{2}>0$ in the relevant
parameter region, which we checked with the analytic expression for $d_{2}$ in
Appendix A. Notice that $\beta_{0}=-4.33$ 666See e.g. Harada (1997) for a
related, analytic estimate of the critical $\beta_{0}$ with constant density
stars. The critical value may vary slightly for neutron stars with other
equations of state. is the critical value below which the spontaneous
scalarization is realized for a certain range of $C_{A}$. In Horbatsch and
Burgess (2011), closed analytic expressions are not found for $d_{2}$ and
$e_{2}$.

Figure 3:  The linear coefficient $d_{1}$ in Eq. (16) as a function of the
stellar compactness $C_{A}$ for constant density stars with $\alpha_{0}=0$ and
various choices of $\beta_{0}$. Spontaneous scalarization is realized when
$d_{1}<0$.

##  III Analytic Scalar Charges

We now explain how we can use the weak-field expansion and Padé resummation to
analytically calculate NS scalar charges with the Tolman VII model.

###  III.1 Formalism

To construct approximate, analytic stellar solutions in the scalar-tensor
theory, we work under a weak-field expansion. Namely, we decompose each
unknown function $f(r)$ as

$f(r)=\sum_{k=0}f_{k}(r)\epsilon^{k}\,,$ (31)

where $\epsilon$ is a book-keeping parameter that counts the order of the GR
compactness $C_{0}=M_{0}/R$ where $M_{0}$ is the GR mass. We have

$m_{0}(r)=\widetilde{\rho}_{0}(r)=\widetilde{p}_{0}(r)=\widetilde{p}_{1}(r)=0\,,$
(32)

and thus to leading order, $m(r)=\mathcal{O}(\epsilon)$,
$\widetilde{\rho}(r)=\mathcal{O}(\epsilon)$ and
$\widetilde{p}(r)=\mathcal{O}(\epsilon^{2})$. Notice that
$\widetilde{\rho}\gg\widetilde{p}$ in the weak-field limit. For the scalar
field, we find it convenient to introduce

$\bar{\varphi}(r)\equiv\frac{\varphi(r)}{\alpha_{0}}\,,$ (33)

so that $\bar{\varphi}(r)=\mathcal{O}(\alpha_{0}^{0})$ when $\alpha_{0}\ll 1$.
We decompose $\bar{\varphi}(r)$ as in Eq. (31) and we have

$\bar{\varphi}_{0}(r)=\frac{1}{\beta_{0}}\,.$ (34)

To see the effect of spontaneous scalarization mentioned in Sec. II.3, we need
to derive the scalar charge $\alpha_{A}$ valid up to
$\mathcal{O}(\alpha_{0}^{3})$ (which, in turn, means that it is valid to
$\mathcal{O}(\alpha_{c}^{3})$ as in Eq. (16)). Namely, we need to find
solutions up to to $\mathcal{O}(\alpha_{0}^{2})$ higher than the leading. For
this, we seek for a solution $\bar{\varphi}$ to $\mathcal{O}(\alpha_{0}^{2})$.
For $m(r)$, since $m_{k}=\mathcal{O}(\alpha_{0}^{2})$ when $k\geq 2$ (namely,
$m(r)=m_{1}(r)$ in GR), we need to find a solution for $m_{k}(r)$ valid to
$\mathcal{O}(\alpha_{0}^{4})$ such that $\alpha_{A}$ is valid up to
$\mathcal{O}(\alpha_{0}^{3})$. Since $\widetilde{p}$ does not enter in the
$m^{\prime}(r)$ equation in Eq. (8), we only need a solution for
$\widetilde{p}$ to $\mathcal{O}(\alpha_{0}^{2})$. Next, we substitute the
weak-field expansion of each unknown function in Eq. (31) to the differential
equations in Eqs. (8)–(10) and solve order by order in $\epsilon$. Below, we
present the equations and solutions mainly to leading order in $\epsilon$ for
an illustration purpose. In the actual calculation, we derived higher order
contributions and applied a Padé resummation that we explain later.

###  III.2 Interior Solutions

Figure 4:  Energy density profiles (normalized by the central value) for NSs
with $1.8M_{\odot}$ in GR and in the scalar-tensor theory with
$\alpha_{0}=10^{-5}$ and $\beta_{0}=-4.5$ (corresponding to cyan dots in Fig.
2). We use the SLy EoS. For comparison, we also show the profile for the
Tolman VII model given in Eq. (35). Observe that such a model can accurately
describe the profile for realistic NSs.

Since realistic EoSs are given in tables, one can only construct stellar
solutions numerically. To overcome this, we use the Tolman VII model Tolman
(1939); Jiang and Yagi (2019) that is known to accurately model realistic NSs
in GR.

The density profile is given by a simple quadratic form as

$\widetilde{\rho}_{1}=\frac{15M_{0}}{8\pi
R^{3}}\left(1-\frac{r^{2}}{R^{2}}\right)\,,\quad\widetilde{\rho}_{i}(r)=0\quad(i\geq
2)\,.$ (35)

Thus, these stars are parameterized by $M_{0}$ and $R$, where the former is
the stellar mass in GR. To check the validity of the above profile, we present
in Fig. 4 the normalized energy density profile for NSs with $1.8M_{\odot}$
and the SLy EoS found numerically for both GR and scalar-tensor theory. For
the latter, the star is spontaneously scalarized (see Fig. 2). Notice that the
energy density profiles in both theories are almost identical and can be
approximated by Eq. (35) shown by the black solid curve. This justifies the
use of the Tolman VII model even for NSs in scalar-tensor theory.

We next derive the field equations and solutions in the interior region. The
leading differential equations are given by

$\displaystyle m_{1}^{\prime}(r)$ $\displaystyle=$
$\displaystyle\frac{15}{2}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)\frac{M_{0}}{R^{5}}r^{2}(R^{2}-r^{2})+\mathcal{O}(\alpha_{0}^{6})\,,$
$\displaystyle\bar{\varphi}_{1}^{\prime\prime}(r)$ $\displaystyle=$
$\displaystyle-\frac{2\bar{\varphi}_{1}^{\prime}}{r}+\frac{15}{2}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)\frac{M_{0}}{R^{5}}(R^{2}-r^{2})+\mathcal{O}(\alpha_{0}^{4})\,,$
$\displaystyle\widetilde{P}_{2}^{\prime}(r)$ $\displaystyle=$
$\displaystyle-\frac{15M_{0}}{8\pi
R^{5}}(R^{2}-r^{2})\left(\frac{m_{1}}{r^{2}}+\alpha_{0}^{2}\bar{\varphi}_{1}^{\prime}\right)\,.$
(38)

Imposing the boundary conditions

$m_{k}(0)=0\,,\quad\bar{\varphi}_{k}(0)=\bar{\varphi}_{1c}\,,\quad\widetilde{P}_{k}(R)=0\,,$
(39)

we find

$\displaystyle m_{1}^{\mathrm{(int)}}(r)$ $\displaystyle=$
$\displaystyle\frac{5}{2}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)\frac{M_{0}}{R^{3}}r^{3}\left(1-\frac{3}{5}\frac{r^{2}}{R^{2}}\right)$
(40) $\displaystyle+\mathcal{O}(\alpha_{0}^{6})\,,$
$\displaystyle\bar{\varphi}_{1}^{\mathrm{(int)}}(r)$ $\displaystyle=$
$\displaystyle\bar{\varphi}_{1c}+\frac{5}{4}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)\frac{M_{0}}{R^{3}}r^{2}\left(1-\frac{3}{10}\frac{r^{2}}{R^{2}}\right)$
(41) $\displaystyle+\mathcal{O}(\alpha_{0}^{4})\,,$
$\displaystyle\widetilde{P}_{2}^{\mathrm{(int)}}(r)$ $\displaystyle=$
$\displaystyle\frac{15}{16\pi}\left(1+\frac{4+\beta_{0}}{\beta_{0}}\alpha_{0}^{2}\right)\frac{M_{0}^{2}}{R^{4}}$
(42)
$\displaystyle\times\left(1-\frac{r^{2}}{R^{2}}\right)^{2}\left(1-\frac{1}{2}\frac{r^{2}}{R^{2}}\right)+\mathcal{O}(\alpha_{0}^{4})\,.$

###  III.3 Exterior Solutions and Perturbative Scalar Charges

Next, we study the exterior region. By setting
$\widetilde{P}=\widetilde{\rho}=0$ in Eqs. (8)–(10), we find

$m_{1}^{\prime}(r)=0\,,\quad\bar{\varphi}_{1}^{\prime\prime}(r)=-\frac{2\bar{\varphi}_{1}^{\prime}}{r}\,.$
(43)

Imposing regularity at infinity, we can solve the above equations to find the
exterior solution as777An analytic solution for a NS exterior in scalar-tensor
theories without the weak-field expansion has been found in a different
coordinate system Damour and Esposito-Farese (1992, 1996). We found that
working in the same coordinate system as the interior case is easier due to
some nonlinearity in the coordinate transformation.

$m_{1}^{\mathrm{(ext)}}(r)=M_{1}\,,\quad\bar{\varphi}_{1}^{\mathrm{(ext)}}(r)=\frac{\bar{\omega}_{1}}{r}\,,$
(44)

where $M_{1}$ and $\bar{\omega}_{1}$ are integration constants. We can
determine the integration constants both in the interior and exterior
solutions by imposing the boundary conditions at the surface:

$\displaystyle m_{k}^{\mathrm{(int)}}(R)$ $\displaystyle=$ $\displaystyle
m_{k}^{\mathrm{(ext)}}(R)\,,$ (45)
$\displaystyle\bar{\varphi}_{k}^{\mathrm{(int)}}(R)$ $\displaystyle=$
$\displaystyle\bar{\varphi}_{k}^{\mathrm{(ext)}}(R)\,,$ (46)
$\displaystyle\bar{\varphi}_{k}^{\prime\mathrm{(int)}}(R)$ $\displaystyle=$
$\displaystyle\bar{\varphi}_{k}^{\prime\mathrm{(ext)}}(R)\,.$ (47)

Similar to Eq. (44), we can introduce integration constants $M_{k}$ and
$\bar{\omega}_{k}$ at order $\epsilon^{k}$ as

$\displaystyle m_{k}^{\mathrm{(ext)}}(r)$ $\displaystyle=$ $\displaystyle
M_{k}\left[1+\mathcal{O}\left(\frac{M_{0}}{r}\right)\right]\,,$ (48)
$\displaystyle\bar{\varphi}_{k}^{\mathrm{(ext)}}(r)$ $\displaystyle=$
$\displaystyle\frac{\bar{\omega}_{k}}{r}+\mathcal{O}\left(\frac{M_{0}^{2}}{r^{2}}\right)\,.$
(49)

Once the integration constants are determined, we can compute the perturbative
scalar charge from Eq. (12) as

$\alpha_{A}=-\alpha_{0}\frac{\sum_{k=1}^{N}\bar{\omega}_{k}\epsilon^{k}}{\sum_{k=1}^{N}M_{k}\epsilon^{k}}\,.$
(50)

In this paper, we computed up to $N=10$. Equation (50) allows us to express
$\alpha_{A}$ in a series of the GR compactness $C_{0}$. It would be more
useful to express the scalar charges in terms of the physical compactness in
the scalar-tensor theory defined by

$C_{A}=\frac{M_{A}}{R}\,,\quad M_{A}=\sum_{k=1}M_{k}\epsilon^{k}\,.$ (51)

We can solve order by order to find $C_{0}$ in terms of $C_{A}$ and plug this
into Eq. (50). To have the results consistent up to the order analyzed, we
expand Eq. (50) about $\alpha_{0}=0$ and $C_{A}=0$ and keep only to
$\mathcal{O}(\alpha_{0}^{3},C_{A}^{N})$. To have the series converge, we then
construct a Padé approximant of order $N/2$ (when $N$ is an even number) in
$C_{A}$.

At leading order, the integration constants are derived as

$\displaystyle\bar{\varphi}_{1c}$ $\displaystyle=$
$\displaystyle-\frac{15}{8}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)\frac{M_{0}}{R}+\mathcal{O}(\alpha_{0}^{4})\,.$
(52) $\displaystyle M_{1}$ $\displaystyle=$
$\displaystyle\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)M_{0}+\mathcal{O}(\alpha_{0}^{6})\,,$
(53) $\displaystyle\bar{\omega}_{1}$ $\displaystyle=$
$\displaystyle-\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)M_{0}+\mathcal{O}(\alpha_{0}^{4})\,.$
(54)

Using Eq. (50), the scalar charge to this order is given by

$\displaystyle\alpha_{A}$ $\displaystyle=$
$\displaystyle\alpha_{0}\frac{\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)}{\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)}$
(55) $\displaystyle=$
$\displaystyle\alpha_{0}+\mathcal{O}(\alpha_{0}^{5},C_{A})\,.$

We have derived $\alpha_{A}$ valid to $\mathcal{O}(C_{A}^{10})$. The first few
terms are given by

$\displaystyle\alpha_{A}$ $\displaystyle=$
$\displaystyle\alpha_{0}-\frac{10}{7}\alpha_{0}\left(\beta_{0}+1\right)C_{A}+\frac{5\alpha_{0}}{3003}$
$\displaystyle\times\left(1253\beta_{0}^{2}+1514\beta_{0}-1126\right)C_{A}^{2}+\mathcal{O}(\alpha_{0}^{3},C_{A}^{3})\,.$

We then construct a Padé approximant to 5th order in $C_{A}$, which we provide
in a Supplemental Mathematica notebook mat .

###  III.4 Spontaneous Scalarization

So far, we have constructed $\alpha_{A}$ in a series of $\alpha_{0}$ and
$C_{A}$ in a form

$\alpha_{A}=\bar{e}_{1}(C_{A})\alpha_{0}+\bar{e}_{3}(C_{A})\alpha_{0}^{3}+\mathcal{O}(\alpha_{0}^{5})\,,$
(57)

where $\bar{e}_{1}$ and $\bar{e}_{3}$ are functions of $C_{A}$ valid to
$\mathcal{O}(C_{A}^{N})$. We also have

$\displaystyle\alpha_{c}$ $\displaystyle=$
$\displaystyle\alpha_{0}\beta_{0}\sum_{k=1}^{N}\bar{\varphi}_{kc}\epsilon^{k}$
(58) $\displaystyle=$
$\displaystyle\bar{d}_{1}(C_{A})\alpha_{0}+\bar{d}_{3}(C_{A})\alpha_{0}^{3}+\mathcal{O}(\alpha_{0}^{5})\,.$

We can invert Eq. (58) order by order in $\alpha_{c}$ to find
$\alpha_{0}(\alpha_{c})$ in a form given by Eq. (16). Then, we substitute this
to Eq. (57) and expand about $\alpha_{c}=0$ to find $\alpha_{A}(\alpha_{c})$
in a form in Eq. (17). Following Sec. II.3, we can then find $\alpha_{A}$ for
NSs under spontaneous scalarization.

In fact, it turns out that the contribution for the second term in Eq. (22) is
negligible and we focus on deriving $d_{1}$, $d_{2}$, and $e_{1}$. We first
find these coefficients in terms of a series expansion of $C_{A}$ valid to
$\mathcal{O}(C_{A}^{10})$. The first few terms of these functions are given by

$\displaystyle d_{1}$ $\displaystyle=$ $\displaystyle
1+\frac{15\beta_{0}}{8}C_{A}+\frac{5\beta_{0}}{128}(17\beta_{0}-14)C_{A}^{2}+\frac{\beta_{0}}{7168}\left(695\beta_{0}^{2}-2294\beta_{0}-27720\right)C_{A}^{3}+\mathcal{O}(C_{A}^{4})\,,$
(59) $\displaystyle d_{2}$ $\displaystyle=$
$\displaystyle-\frac{1955\beta_{0}}{448}C_{A}^{2}+\frac{(1630328-3417469\beta_{0})\beta_{0}}{512512}C_{A}^{3}+\mathcal{O}(C_{A}^{4})\,,$
(60) $\displaystyle e_{1}$ $\displaystyle=$ $\displaystyle
1+\frac{5}{56}(5\beta_{0}-16)C_{A}+\frac{5}{384384}\left(5515\beta_{0}^{2}-54170\beta_{0}-144128\right)C_{A}^{2}$
(61)
$\displaystyle+\frac{1}{156828672}\left(941985\beta_{0}^{3}-19472682\beta_{0}^{2}-258360664\beta_{0}-441114624\right)C_{A}^{3}+\mathcal{O}(C_{A}^{4})\,.$

Next, we construct the Padé approximants on these. We first tried 5th order
Padé approximants but found that there are some unphysical divergences (Fig.
5). Instead, we use 4th order Padé approximants. We substituted these into the
following expression to find scalar charges for NSs under spontaneous
scalarization:

$\alpha_{A}=\pm\left(-\frac{d_{1}}{d_{2}}\right)^{1/2}e_{1}\,.$ (62)

The final expression is available in the Supplemental Mathematica notebook mat
.

##  IV Comparison with Numerical Results

Let us now compare the analytic scalar charge expression for the Tolman VII
model with scalar charges for realistic NSs found numerically. Scalar charges
can be both positive and negative but we focus on the former.

###  IV.1 Scalar Charges vs Compactness

Figure 5:  Convergence of scalar charge calculations with the Tolman model for
$(\alpha_{0},\beta_{0})=(10^{-5},-4.2)$ (top) and
$(\alpha_{0},\beta_{0})=(0,-4.5)$ (bottom) under various approximations. We
compare the Taylor series results and various Padé resummation results against
numerical results for the SLy EoS. Observe that Padé resummation is crucial to
accurately model the numerical result. The expressions for the black solid
curves are the ones that we will use in the remaining part of this paper.

First, let us check the convergence of the analytic result for scalar charges
in terms of the order of the Padé resummation. Figure 5 compares the analytic
scalar charges at various orders of the Padé approximants. We show the results
for the perturbative scalarization case
$(\alpha_{0},\beta_{0})=(10^{-5},-4.2)$ and the spontaneous scalarization case
$(\alpha_{0},\beta_{0})=(0,-4.5)$. We also present Taylor series (without the
Padé resummation) up to $\mathcal{O}(C_{A}^{10})$ for perturbative
scalarization and to $\mathcal{O}(C_{A}^{8})$ for spontaneous scalarization,
as well as numerical results computed with the SLy EoS. Observe that for the
perturbative scalarization case (top panel), the Padé resummation at 3rd order
is already a good approximation and is converging fast to the numerical
result888The 3rd and 5th order Padé results for the _perturbative_ scalar
charge expressions have a noticeable difference for lower $\beta_{0}$ where
the expressions diverge due to the onset of spontaneous scalarization (not
shown in Fig. 5) and thus we use the 5th order result in the remainder of this
section.. On the other hand, the Taylor series result becomes less accurate
when the compactness is relatively high. For the spontaneous scalarization
case (bottom panel), there is little difference between the 4th and 5th order
Padé results, though the latter has some unphysical divergence. We therefore
use the former to avoid this divergence. Additionally, the former has the
advantages of a simpler functional expression and better agreement with the
numerical SLy data. All in all, these findings suggest that it is crucial to
perform the resummation to find an accurate modeling.

Having understood the convergence, let us now carry out the comparison in more
detail for various values of $\beta_{0}$ and EoSs. Figure 1 shows such a
comparison. Observe that the analytic scalar charges for the Tolman VII case
beautifully captures the numerical results for realistic NSs, especially for
the SLy EoS. Notice that the agreement between the analytic and numerical
results is good even for stars with spontaneous scalarization. This shows that
the analytic result serves as an accurate, ready-to-use expression for the NS
scalar charge in this scalar-tensor theory.

###  IV.2 Scalar Charges vs Binding Energy

We next look at the relation between the scalar charge and the stellar binding
energy, where the latter is defined by

$\displaystyle\Omega_{A}$ $\displaystyle=$ $\displaystyle-\frac{1}{2}\int
d^{3}x\widetilde{\rho}(r)\int
d^{3}x^{\prime}\frac{\widetilde{\rho}\left(r^{\prime}\right)}{\left|\bm{x}-\bm{x}^{\prime}\right|}\,,$
(63) $\displaystyle=$
$\displaystyle-16\pi^{2}\int_{0}^{R}drr^{2}\widetilde{\rho}(r)\left(\int_{r}^{R}dr^{\prime}r^{\prime}\widetilde{\rho}(r^{\prime})\right)\,,$
(64)

where the second equality is valid only for spherically symmetric systems.
Physically, this quantity measures the difference between the gravitational
and baryonic mass of a star and can be used to e.g. probe certain formation
scenario of NSs Newton  _et al._ (2018). For the Tolman VII model, the
relation between the binding energy and compactness is given by

$\frac{\Omega_{A}}{M_{A}}=-\frac{5}{7}C_{A}+\mathcal{O}(\alpha_{0}^{2})\,.$
(65)

Such a GR relation is sufficient to find the scalar charge expression in terms
of the binding energy for perturbative scalar charges since the latter is
already proportional to $\alpha_{0}$. We can invert this relation and
substitute it to the $\alpha_{A}(C_{A})$ expressions.

Figure 6:  Similar to Fig. 1 but as a function of the binding energy. Notice
that the results are quasi-universal and insensitive to the EoSs for fixed
$\beta_{0}$.

Figure 6 presents the scalar charge as a function of the binding energy
(normalized by the stellar mass) for the same choices of $\alpha_{0}$ and
$\beta_{0}$ and EoSs as in Fig. 1. Observe that numerical results are quasi-
universal among different EoSs for a fixed $\beta_{0}$. Observe also that the
analytic result accurately describes the numerical ones for perturbative
scalarization. By taking the fractional difference between the numerical and
analytic results for $\beta_{0}=-4.2$ or $\beta_{0}=-4.3$, we found that the
quasi-universality holds to $\sim 1\%$ for $\Omega_{A}/M_{A}>-0.2$.

Analytic expressions help us to study the quasi-universal relation in more
detail. When we expand the perturbative scalar charge expression for the
Tolman VII model for small binding energy, we find

$\displaystyle\alpha_{A}^{(\mathrm{Tol})}$ $\displaystyle=$
$\displaystyle\alpha_{0}+2\alpha_{0}\left(\beta_{0}+1\right)\frac{\Omega_{A}}{M_{A}}$
(66)
$\displaystyle+\frac{7\alpha_{0}}{2145}\left(1253\beta_{0}^{2}+1514\beta_{0}-1126\right)\frac{\Omega_{A}^{2}}{M_{A}^{2}}$
$\displaystyle+\mathcal{O}\left(\alpha_{0}^{3},\frac{\Omega_{A}^{3}}{M_{A}^{3}}\right)\,.$

In Appendix A, we derive analytic scalar charges for constant density stars.
We show that although $\alpha_{A}(C_{A})$ is quite different from those for
Tolman VII and realistic EoS cases, $\alpha_{A}(\Omega_{A}/M_{A})$ is similar
to the latter two cases, supporting the quasi-universal relation. For the
perturbative scalar charge, the Taylor-series expansion in small binding
energy is given by

$\displaystyle\alpha_{A}^{(\mathrm{CD})}$ $\displaystyle=$
$\displaystyle\alpha_{0}+2\alpha_{0}\left(\beta_{0}+1\right)\frac{\Omega_{A}}{M_{A}}$
(67)
$\displaystyle+\frac{5\alpha_{0}}{21}\left(17\beta_{0}^{2}+20\beta_{0}-16\right)\frac{\Omega_{A}^{2}}{M_{A}^{2}}$
$\displaystyle+\mathcal{O}\left(\alpha_{0}^{3},\frac{\Omega_{A}^{3}}{M_{A}^{3}}\right)\,.$

Comparing this with Eq. (66), we see that the two expressions are identical up
to $\mathcal{O}(\Omega_{A}/M_{A})$ and the difference only appears at
$\mathcal{O}(\Omega_{A}^{2}/M_{A}^{2})$. When taking the fractional difference
between $\alpha_{A}^{(\mathrm{Tol})}$ and $\alpha_{A}^{(\mathrm{CD})}$, we
find

$\displaystyle\frac{\alpha_{A}^{(\mathrm{Tol})}-\alpha_{A}^{(\mathrm{CD})}}{\alpha_{0}}$
$\displaystyle=$
$\displaystyle\frac{2\left(311\beta_{0}^{2}+1343\beta_{0}+1013\right)}{15015}\frac{\Omega_{A}^{2}}{M_{A}^{2}}$
(68)
$\displaystyle+\mathcal{O}\left(\alpha_{0}^{3},\frac{\Omega_{A}^{3}}{M_{A}^{3}}\right)\,.$

Figure 7 presents this fractional difference for $\beta_{0}=-4.2$ and $-4.3$.
Observe that the Tolman and constant density cases agree within an error of
$\mathcal{O}(10^{-3})$. This provides strong analytic support for the quasi-
universality of the relation between $\alpha_{A}$ and $\Omega_{A}/M_{A}$.

Figure 7:  Fractional difference in the scalar charge between Tolman VII and
constant density cases in terms of the binding energy (Eq. (68)) for two
choices of $\beta_{0}$. Notice that the two agree within an error of
$\mathcal{O}(10^{-3})$, indicating the quasi-universality of the relation
between the scalar charge and binding energy.

For spontaneous scalarization, although the analytic expressions are valid as
an order-of-magnitude estimate, the agreement with numerical results are not
as good as the perturbative scalarization case. This is partially because we
have used the relation between the binding energy and compactness in GR in Eq.
(65). Though, due to multi-valued scalar charges for some fixed binding energy
with small $\beta_{0}$ (e.g. $\beta_{0}=-4.5$), it would be difficult to find
a closed, analytic expression for $\alpha_{A}(\Omega_{A}/M_{A})$.

##  V Conclusion

We derived analytically scalar charges for NSs in a scalar-tensor theory
proposed in Damour and Esposito-Farese (1993). This was achieved by
considering the Tolman VII energy density profile. We first worked within the
weak-field approximation and then resummed the series through Padé
approximants, deriving scalar charge expressions for both perturbative and
spontaneous scalarization. We found that the analytic scalar charges are in
excellent agreement with those computed numerically from realistic EoSs,
especially for the SLy EoS. We also found a quasi-universal relation between
the scalar charge and binding energy. The analytic expressions derived here
allow one to mathematically support the quasi-universality by comparing the
Tolman VII result with the constant density one. The analytic result provides
an accurate, ready-to-use, and physically-motivated expression for scalar
charges.

A similar quasi-universality between the stellar sensitivities and the binding
energy were found recently in Einstein-Æther theory Gupta  _et al._ (2021).
The relation in the weak-field limit was first derived in Foster (2007). This
was based on the result for weakly-gravitating stars in Foster (2006), where
the dipole moment for the vector field (which depends on the sensitivities)
was derived within the parameterized post-Newtonian framework and all the
internal structure dependence was found to be encoded in the binding energy.
It would be interesting to study if a similar analysis can be carried out in
scalar-tensor theories to explain the quasi-universal relation found here.

Various avenues exist for other possible future work. For example, one obvious
extension is to apply the analysis presented here to other scalar-tensor
theories with spontaneous scalarization, such as the ones with the mass
potential Chen  _et al._ (2015); Ramazanoğlu and Pretorius (2016); Yazadjiev
_et al._ (2016); Doneva and Yazadjiev (2016) or quartic interaction Staykov
_et al._ (2018), scalar-Gauss-Bonnet gravity Doneva and Yazadjiev (2018a);
Silva  _et al._ (2018); Doneva and Yazadjiev (2018b); Silva  _et al._ (2019);
Minamitsuji and Ikeda (2019); Macedo  _et al._ (2019); Dima  _et al._ (2020);
Doneva  _et al._ (2020); Herdeiro  _et al._ (2021); Silva  _et al._ (2020b),
or Horndeski theories Andreou  _et al._ (2019). One may also apply the
calculation to spontaneous vectorization Ramazanoğlu (2017); Annulli  _et al._
(2019); Kase  _et al._ (2020), tensorization Ramazanoğlu (2019), or
spinorization Ramazanoğlu (2018). Another possibility might be to consider an
analytic representation of the dynamical/induced scalarization Barausse  _et
al._ (2013); Palenzuela  _et al._ (2014); Shibata  _et al._ (2014); Sampson
_et al._ (2014); Sennett and Buonanno (2016); Khalil  _et al._ (2019) in
compact binary mergers. Finally, the parameter region of the coupling
parameter $\beta_{0}$ that gives rise to spontaneous scalarization in the
scalar-tensor theories studied here has been shown to be inconsistent with
solar system experiments if one includes cosmological evolution of the scalar
field Damour and Nordtvedt (1993a, b); Sampson  _et al._ (2014); Anderson  _et
al._ (2016); Anderson and Yunes (2017). It would be interesting to consider
analytic scalar charges in scalar-tensor theories that has a consistent
cosmological evolution of the scalar field Mendes (2015); Mendes and Ortiz
(2016); Antoniou  _et al._ (2021).

###### Acknowledgements.

K.Y. acknowledges support from NSF Grant PHY-1806776, NASA Grant
80NSSC20K0523, a Sloan Foundation Research Fellowship and the Owens Family
Foundation. K.Y. would like to also acknowledge support by the COST Action
GWverse CA16104 and JSPS KAKENHI Grants No. JP17H06358.

##  Appendix A Scalar Charges for Constant Density Stars

In this appendix, we repeat the calculations in the main text to find analytic
expressions for scalar charges for constant density stars. For such stars, we
give the density profile as

$\widetilde{\rho}_{1}=\frac{3M_{0}}{4\pi
R^{3}}\,,\quad\widetilde{\rho}_{i}(r)=0\quad(i\geq 2)\,.$ (69)

Thus, these stars are parameterized by $M_{0}$ and $R$. The leading
differential equations are given by

$\displaystyle m_{1}^{\prime}(r)$ $\displaystyle=$ $\displaystyle
3\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)\frac{M_{0}}{R^{3}}r^{2}+\mathcal{O}(\alpha_{0}^{6})\,,$
(70) $\displaystyle\bar{\varphi}_{1}^{\prime\prime}(r)$ $\displaystyle=$
$\displaystyle-\frac{2\bar{\varphi}_{1}^{\prime}}{r}+3\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)\frac{M_{0}}{R^{3}}+\mathcal{O}(\alpha_{0}^{4})\,,$
(71) $\displaystyle\widetilde{P}_{2}^{\prime}(r)$ $\displaystyle=$
$\displaystyle-\frac{3M_{0}}{4\pi
R^{3}}\left(\frac{m_{1}}{r^{2}}+\alpha_{0}^{2}\bar{\varphi}_{1}^{\prime}\right)\,.$
(72)

The last equation is valid to full order in $\alpha_{0}$. Imposing the
boundary condition as in Eq. (39), one can solve the above differential
equations in the interior region to find

$\displaystyle m_{1}^{\mathrm{(int)}}(r)$ $\displaystyle=$
$\displaystyle\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)\frac{M_{0}}{R^{3}}r^{3}+\mathcal{O}(\alpha_{0}^{6})\,,$
(73) $\displaystyle\bar{\varphi}_{1}^{\mathrm{(int)}}(r)$ $\displaystyle=$
$\displaystyle\bar{\varphi}_{1c}+\frac{1}{2}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)\frac{M_{0}}{R^{3}}r^{2}+\mathcal{O}(\alpha_{0}^{4})\,,$
(74) $\displaystyle\widetilde{P}_{2}^{\mathrm{(int)}}(r)$ $\displaystyle=$
$\displaystyle\frac{3}{8\pi}\left(1+\frac{4+\beta_{0}}{\beta_{0}}\alpha_{0}^{2}\right)\frac{M_{0}^{2}}{R^{4}}\left(1-\frac{r^{2}}{R^{2}}\right)+\mathcal{O}(\alpha_{0}^{4})\,.$

Next, we study the exterior solution and the scalar charge. At leading order
in the weak-field expansion, the integration constants are determined as

$\displaystyle\bar{\varphi}_{1c}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)\frac{M_{0}}{R}+\mathcal{O}(\alpha_{0}^{4})\,.$
(76) $\displaystyle M_{1}$ $\displaystyle=$
$\displaystyle\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}+8\frac{\alpha_{0}^{4}}{\beta_{0}^{2}}\right)M_{0}+\mathcal{O}(\alpha_{0}^{6})\,,$
(77) $\displaystyle\bar{\omega}_{1}$ $\displaystyle=$
$\displaystyle-\left(1+4\frac{\alpha_{0}^{2}}{\beta_{0}}\right)M_{0}+\mathcal{O}(\alpha_{0}^{4})\,.$
(78)

Notice that $M_{1}$ and $\bar{\omega}_{1}$ are the same as the Tolman VII case
as in Eqs. (53) and (54). Thus, the scalar charge is also same as in Eq. (55).
Similar to the Tolman case, we then construct a Padé approximant to 5th order
in $C_{A}$ that we provide in the Supplemental Mathematica notebook mat .

Let us now find the expression for spontaneous scalarization. Although $d_{1}$
and $e_{1}$ are fully given in Eqs. (II.3) and (24), we derived $d_{1}$ and
$e_{1}$ in a 4th order Padé approximant form to make the scalar charge
expression similar to that for the Tolman case discussed in the main text. We
have carried out a similar analysis for $d_{2}$ whose full expression has not
been found yet. The first few terms are given by

$\displaystyle d_{2}=$
$\displaystyle-\frac{21\beta_{0}}{8}C_{A}^{2}-\frac{3}{560}\beta_{0}(547\beta_{0}-278)C_{A}^{3}+\mathcal{O}(C_{A}^{4})\,.$
(79)

We have derived $d_{2}$ to $\mathcal{O}(C_{A}^{8})$ and constructed a 4th
order Padé approximant. We can substitute these Padé resummed forms for
$d_{1}$, $d_{2}$, and $e_{1}$ to Eq. (62) to find the scalar charge for
constant density stars under spontaneous scalarization when $\alpha_{0}\ll 1$.
We provide the final expression in the Supplemental Mathematica notebook mat .

Figure 8:  Similar to Fig. 1 but including analytic results for constant
density stars. Notice that scalar charges for such stars are quite different
from those for realistic NSs.

Figure 8 compares the analytic scalar charges for constant density stars (as a
function of the compactness) with those for Tolman VII model and numerical
charges with two representative EoSs. Notice that the spontaneous
scalarization happens for larger compactnesses compared to the Tolman VII
model and realistic NSs. This shows that the former is not an accurate model
of the latter.

One can further convert the scalar charge expression in terms of the binding
energy. For constant density stars, the leading relation is
$\Omega_{A}/M_{A}=-(3/5)C_{A}+\mathcal{O}(\alpha_{0}^{2})$.

Figure 9:  Similar to Fig. 6 but including analytic results for constant
density stars. Unlike in Fig. 8, scalar charges for constant density stars are
now similar to realistic NSs.

Inverting this and substituting it into the analytic expression for the scalar
charge in terms of the compactness, one finds the scalar charge as a function
of the binding energy, which is shown in Fig. 9. We compare the constant
density result with the Tolman VII one and numerical results with realistic
EoSs. Observe that the scalar charges for constant density stars are very
similar to other results, supporting the quasi-universality of the relation
between the scalar charge and binding energy.

## References

  * Cromartie  _et al._ (2019) H. T. Cromartie _et al._ (NANOGrav), “Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar,” Nature Astron. 4, 72–76 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , “A $NICER$ View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation,” Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , “PSR J0030+0451 Mass and Radius from $NICER$ Data and Implications for the Properties of Neutron Star Matter,” Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (Virgo, LIGO Scientific), “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “GW170817: Measurements of neutron star radii and equation of state,” Phys. Rev. Lett. 121, 161101 (2018).
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, “Gravitational-wave constraints on the neutron-star-matter Equation of State,” Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Raithel  _et al._ (2018) C. Raithel, F. Ozel, and D. Psaltis, “Tidal deformability from GW170817 as a direct probe of the neutron star radius,” Astrophys. J. 857, L23 (2018), arXiv:1803.07687 [astro-ph.HE] .
  * Lim and Holt (2018) Y. Lim and J. W. Holt, “Neutron star tidal deformabilities constrained by nuclear theory and experiment,” Phys. Rev. Lett. 121, 062701 (2018), arXiv:1803.02803 [nucl-th] .
  * Bauswein  _et al._ (2017) A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas, “Neutron-star radius constraints from GW170817 and future detections,” Astrophys. J. 850, L34 (2017), arXiv:1710.06843 [astro-ph.HE] .
  * De  _et al._ (2018) S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, “Tidal Deformabilities and Radii of Neutron Stars from the Observation of GW170817,” Phys. Rev. Lett. 121, 091102 (2018), [Erratum: Phys. Rev. Lett.121,no.25,259902(2018)], arXiv:1804.08583 [astro-ph.HE] .
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, “New constraints on radii and tidal deformabilities of neutron stars from GW170817,” Phys. Rev. Lett. 120, 261103 (2018), arXiv:1803.00549 [gr-qc] .
  * Annala  _et al._ (2019) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, “Quark-matter cores in neutron stars,” (2019), arXiv:1903.09121 [astro-ph.HE] .
  * Malik  _et al._ (2018) T. Malik, N. Alam, M. Fortin, C. Providência, B. K. Agrawal, T. K. Jha, B. Kumar, and S. K. Patra, “GW170817: constraining the nuclear matter equation of state from the neutron star tidal deformability,” Phys. Rev. C98, 035804 (2018), arXiv:1805.11963 [nucl-th] .
  * Carson  _et al._ (2019a) Z. Carson, A. W. Steiner, and K. Yagi, “Constraining nuclear matter parameters with gw170817,” Phys. Rev. D 99, 043010 (2019a).
  * Carson  _et al._ (2019b) Z. Carson, A. W. Steiner, and K. Yagi, “Future Prospects for Constraining Nuclear Matter Parameters with Gravitational Waves,” Phys. Rev. D100, 023012 (2019b), arXiv:1906.05978 [gr-qc] .
  * Raithel and Ozel (2019) C. A. Raithel and F. Ozel, “Measurement of the nuclear symmetry energy parameters from gravitational wave events,”  (2019), 10.3847/1538-4357/ab48e6, arXiv:1908.00018 [astro-ph.HE] .
  * Chatziioannou (2020) K. Chatziioannou, “Neutron star tidal deformability and equation of state constraints,” Gen. Rel. Grav. 52, 109 (2020), arXiv:2006.03168 [gr-qc] .
  * Raaijmakers  _et al._ (2020) G. Raaijmakers _et al._ , “Constraining the dense matter equation of state with joint analysis of NICER and LIGO/Virgo measurements,” Astrophys. J. Lett. 893, L21 (2020), arXiv:1912.11031 [astro-ph.HE] .
  * Zimmerman  _et al._ (2020) J. Zimmerman, Z. Carson, K. Schumacher, A. W. Steiner, and K. Yagi, “Measuring Nuclear Matter Parameters with NICER and LIGO/Virgo,” (2020), arXiv:2002.03210 [astro-ph.HE] .
  * Jiang  _et al._ (2020) J.-L. Jiang, S.-P. Tang, Y.-Z. Wang, Y.-Z. Fan, and D.-M. Wei, “PSR J0030+0451, GW170817 and the nuclear data: joint constraints on equation of state and bulk properties of neutron stars,” Astrophys. J. 892, 1 (2020), arXiv:1912.07467 [astro-ph.HE] .
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, “Multimessenger constraints on the neutron-star equation of state and the Hubble constant,” Science 370, 1450–1453 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Freire  _et al._ (2012) P. C. Freire _et al._ , “The relativistic pulsar-white dwarf binary PSR J1738+0333 II. The most stringent test of scalar-tensor gravity,” Mon.Not.Roy.Astron.Soc. 423, 3328 (2012), arXiv:1205.1450 [astro-ph.GA] .
  * Yagi  _et al._ (2014a) K. Yagi, D. Blas, N. Yunes, and E. Barausse, “Strong Binary Pulsar Constraints on Lorentz Violation in Gravity,” Phys.Rev.Lett. 112, 161101 (2014a), arXiv:1307.6219 [gr-qc] .
  * Yagi  _et al._ (2014b) K. Yagi, D. Blas, E. Barausse, and N. Yunes, “Constraints on Einstein-Æther theory and Horava gravity from binary pulsar observations,” Phys.Rev. D89, 084067 (2014b), arXiv:1311.7144 [gr-qc] .
  * Berti  _et al._ (2015) E. Berti _et al._ , “Testing General Relativity with Present and Future Astrophysical Observations,” Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Shao  _et al._ (2017) L. Shao, N. Sennett, A. Buonanno, M. Kramer, and N. Wex, “Constraining nonperturbative strong-field effects in scalar-tensor gravity by combining pulsar timing and laser-interferometer gravitational-wave detectors,” Phys. Rev. X7, 041025 (2017), arXiv:1704.07561 [gr-qc] .
  * Archibald  _et al._ (2018) A. M. Archibald, N. V. Gusinskaia, J. W. T. Hessels, A. T. Deller, D. L. Kaplan, D. R. Lorimer, R. S. Lynch, S. M. Ransom, and I. H. Stairs, “Universality of free fall from the orbital motion of a pulsar in a stellar triple system,” Nature 559, 73–76 (2018), arXiv:1807.02059 [astro-ph.HE] .
  * Seymour and Yagi (2020) B. C. Seymour and K. Yagi, “Probing Massive Scalar Fields from a Pulsar in a Stellar Triple System,” Class. Quant. Grav. 37, 145008 (2020), arXiv:1908.03353 [gr-qc] .
  * Anderson  _et al._ (2019) D. Anderson, P. Freire, and N. Yunes, “Binary pulsar constraints on massless scalar–tensor theories using Bayesian statistics,” Class. Quant. Grav. 36, 225009 (2019), arXiv:1901.00938 [gr-qc] .
  * Gupta  _et al._ (2021) T. Gupta, M. Herrero-Valea, D. Blas, E. Barausse, N. Cornish, K. Yagi, and N. Yunes, “Updated Binary Pulsar Constraints on Einstein-æther Theory in Light of Gravitational Wave Constraints on the Speed of Gravity,” (2021), arXiv:2104.04596 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of General Relativity with GW170817,” Phys. Rev. Lett. 123, 011102 (2019), arXiv:1811.00364 [gr-qc] .
  * Baker  _et al._ (2017) T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller, and I. Sawicki, “Strong constraints on cosmological gravity from GW170817 and GRB 170817A,” Phys. Rev. Lett. 119, 251301 (2017), arXiv:1710.06394 [astro-ph.CO] .
  * Sakstein and Jain (2017) J. Sakstein and B. Jain, “Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories,” Phys. Rev. Lett. 119, 251303 (2017), arXiv:1710.05893 [astro-ph.CO] .
  * Ezquiaga and Zumalacárregui (2017) J. M. Ezquiaga and M. Zumalacárregui, “Dark Energy After GW170817: Dead Ends and the Road Ahead,” Phys. Rev. Lett. 119, 251304 (2017), arXiv:1710.05901 [astro-ph.CO] .
  * Zhao  _et al._ (2019) J. Zhao, L. Shao, Z. Cao, and B.-Q. Ma, “Reduced-order surrogate models for scalar-tensor gravity in the strong field regime and applications to binary pulsars and GW170817,” Phys. Rev. D100, 064034 (2019), arXiv:1907.00780 [gr-qc] .
  * Silva  _et al._ (2020a) H. O. Silva, A. M. Holgado, A. Cárdenas-Avendaño, and N. Yunes, “Astrophysical and theoretical physics implications from multimessenger neutron star observations,” (2020a), arXiv:2004.01253 [gr-qc] .
  * Damour and Esposito-Farese (1992) T. Damour and G. Esposito-Farese, “Tensor multiscalar theories of gravitation,” Class.Quant.Grav. 9, 2093–2176 (1992).
  * Damour and Esposito-Farese (1993) T. Damour and G. Esposito-Farese, “Nonperturbative strong field effects in tensor - scalar theories of gravitation,” Phys.Rev.Lett. 70, 2220–2223 (1993).
  * Damour and Esposito-Farese (1996) T. Damour and G. Esposito-Farese, “Tensor \- scalar gravity and binary pulsar experiments,” Phys.Rev. D54, 1474–1491 (1996), arXiv:gr-qc/9602056 [gr-qc] .
  * Ramazanoğlu (2017) F. M. Ramazanoğlu, “Spontaneous growth of vector fields in gravity,” Phys. Rev. D96, 064009 (2017), arXiv:1706.01056 [gr-qc] .
  * Anderson and Yunes (2019) D. Anderson and N. Yunes, “Scalar charges and scaling relations in massless scalar–tensor theories,” Class. Quant. Grav. 36, 165003 (2019), arXiv:1901.00937 [gr-qc] .
  * Carson  _et al._ (2020) Z. Carson, B. C. Seymour, and K. Yagi, “Future Prospects for Probing Scalar-Tensor Theories with Gravitational Waves from Mixed Binaries,” Class. Quant. Grav. 37, 065008 (2020), arXiv:1907.03897 [gr-qc] .
  * Tolman (1939) R. C. Tolman, “Static solutions of Einstein’s field equations for spheres of fluid,” Phys.Rev. 55, 364–373 (1939).
  * Yagi  _et al._ (2016) K. Yagi, L. C. Stein, and N. Yunes, “Challenging the Presence of Scalar Charge and Dipolar Radiation in Binary Pulsars,” Phys. Rev. D93, 024010 (2016), arXiv:1510.02152 [gr-qc] .
  * Wiringa  _et al._ (1988) R. B. Wiringa, V. Fiks, and A. Fabrocini, “Equation of state for dense nucleon matter,” Phys.Rev. C38, 1010–1037 (1988).
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, “A unified equation of state of dense matter and neutron star structure,” Astron. Astrophys. 380, 151–167 (2001).
  * Müther  _et al._ (1987) H. Müther, M. Prakash, and T. L. Ainsworth, “The nuclear symmetry energy in relativistic Brueckner-Hartree-Fock calculations,” Physics Letters B 199, 469–474 (1987).
  * Mueller and Serot (1996) H. Mueller and B. D. Serot, “Relativistic mean field theory and the high density nuclear equation of state,” Nucl.Phys. A606, 508–537 (1996), arXiv:nucl-th/9603037 [nucl-th] .
  * Pappas  _et al._ (2019) G. Pappas, D. D. Doneva, T. P. Sotiriou, S. S. Yazadjiev, and K. D. Kokkotas, “Multipole moments and universal relations for scalarized neutron stars,” Phys. Rev. D 99, 104014 (2019), arXiv:1812.01117 [gr-qc] .
  * Horbatsch and Burgess (2011) M. W. Horbatsch and C. P. Burgess, “Semi-Analytic Stellar Structure in Scalar-Tensor Gravity,” JCAP 1108, 027 (2011), arXiv:1006.4411 [gr-qc] .
  * (51) https://github.com/KentYagi/analytic-scalar-charges.git.
  * Pani and Berti (2014) P. Pani and E. Berti, “I-Love-Q, Spontaneously: Slowly Rotating Neutron Stars in Scalar-Tensor Theories,” Phys.Rev. D90, 024025 (2014), arXiv:1405.4547 [gr-qc] .
  * Harada (1997) T. Harada, “Stability analysis of spherically symmetric star in scalar - tensor theories of gravity,” Prog. Theor. Phys. 98, 359–379 (1997), arXiv:gr-qc/9706014 .
  * Jiang and Yagi (2019) N. Jiang and K. Yagi, “Improved Analytic Modeling of Neutron Star Interiors,” Phys. Rev. D99, 124029 (2019), arXiv:1904.05954 [gr-qc] .
  * Newton  _et al._ (2018) W. G. Newton, A. W. Steiner, and K. Yagi, “Testing the formation scenarios of binary neutron star systems with measurements of the neutron star moment of inertia,” Astrophys. J. 856, 19 (2018), arXiv:1611.09399 [astro-ph.SR] .
  * Foster (2007) B. Z. Foster, “Strong field effects on binary systems in Einstein-aether theory,” Phys. Rev. D76, 084033 (2007), arXiv:0706.0704 [gr-qc] .
  * Foster (2006) B. Z. Foster, “Radiation damping in Einstein-aether theory,” Phys. Rev. D73, 104012 (2006), [Erratum: Phys. Rev.D75,129904(2007)], arXiv:gr-qc/0602004 [gr-qc] .
  * Chen  _et al._ (2015) P. Chen, T. Suyama, and J. Yokoyama, “Spontaneous scalarization: asymmetron as dark matter,” Phys. Rev. D92, 124016 (2015), arXiv:1508.01384 [gr-qc] .
  * Ramazanoğlu and Pretorius (2016) F. M. Ramazanoğlu and F. Pretorius, “Spontaneous Scalarization with Massive Fields,” Phys. Rev. D93, 064005 (2016), arXiv:1601.07475 [gr-qc] .
  * Yazadjiev  _et al._ (2016) S. S. Yazadjiev, D. D. Doneva, and D. Popchev, “Slowly rotating neutron stars in scalar-tensor theories with a massive scalar field,” Phys. Rev. D93, 084038 (2016), arXiv:1602.04766 [gr-qc] .
  * Doneva and Yazadjiev (2016) D. D. Doneva and S. S. Yazadjiev, “Rapidly rotating neutron stars with a massive scalar field—structure and universal relations,” JCAP 1611, 019 (2016), arXiv:1607.03299 [gr-qc] .
  * Staykov  _et al._ (2018) K. V. Staykov, D. Popchev, D. D. Doneva, and S. S. Yazadjiev, “Static and slowly rotating neutron stars in scalar–tensor theory with self-interacting massive scalar field,” Eur. Phys. J. C78, 586 (2018), arXiv:1805.07818 [gr-qc] .
  * Doneva and Yazadjiev (2018a) D. D. Doneva and S. S. Yazadjiev, “New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories,” 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, “Spontaneous scalarization of black holes and compact stars from a Gauss-Bonnet coupling,” Phys. Rev. Lett. 120, 131104 (2018), arXiv:1711.02080 [gr-qc] .
  * Doneva and Yazadjiev (2018b) D. D. Doneva and S. S. Yazadjiev, “Neutron star solutions with curvature induced scalarization in the extended Gauss-Bonnet scalar-tensor theories,” JCAP 1804, 011 (2018b), arXiv:1712.03715 [gr-qc] .
  * Silva  _et al._ (2019) H. O. Silva, C. F. B. Macedo, T. P. Sotiriou, L. Gualtieri, J. Sakstein, and E. Berti, “Stability of scalarized black hole solutions in scalar-Gauss-Bonnet gravity,” Phys. Rev. D99, 064011 (2019), arXiv:1812.05590 [gr-qc] .
  * Minamitsuji and Ikeda (2019) M. Minamitsuji and T. Ikeda, “Scalarized black holes in the presence of the coupling to Gauss-Bonnet gravity,” Phys. Rev. D99, 044017 (2019), arXiv:1812.03551 [gr-qc] .
  * Macedo  _et al._ (2019) C. F. B. Macedo, J. Sakstein, E. Berti, L. Gualtieri, H. O. Silva, and T. P. Sotiriou, “Self-interactions and Spontaneous Black Hole Scalarization,” Phys. Rev. D99, 104041 (2019), arXiv:1903.06784 [gr-qc] .
  * Dima  _et al._ (2020) A. Dima, E. Barausse, N. Franchini, and T. P. Sotiriou, “Spin-induced black hole spontaneous scalarization,” Phys. Rev. Lett. 125, 231101 (2020), arXiv:2006.03095 [gr-qc] .
  * Doneva  _et al._ (2020) D. D. Doneva, L. G. Collodel, C. J. Krüger, and S. S. Yazadjiev, “Black hole scalarization induced by the spin: 2+1 time evolution,” Phys. Rev. D 102, 104027 (2020), arXiv:2008.07391 [gr-qc] .
  * Herdeiro  _et al._ (2021) C. A. R. Herdeiro, E. Radu, H. O. Silva, T. P. Sotiriou, and N. Yunes, “Spin-induced scalarized black holes,” Phys. Rev. Lett. 126, 011103 (2021), arXiv:2009.03904 [gr-qc] .
  * Silva  _et al._ (2020b) H. O. Silva, H. Witek, M. Elley, and N. Yunes, “Dynamical scalarization and descalarization in binary black hole mergers,” (2020b), arXiv:2012.10436 [gr-qc] .
  * Andreou  _et al._ (2019) N. Andreou, N. Franchini, G. Ventagli, and T. P. Sotiriou, “Spontaneous scalarization in generalised scalar-tensor theory,” Phys. Rev. D99, 124022 (2019), [Erratum: Phys. Rev.D101,no.10,109903(2020)], arXiv:1904.06365 [gr-qc] .
  * Annulli  _et al._ (2019) L. Annulli, V. Cardoso, and L. Gualtieri, “Electromagnetism and hidden vector fields in modified gravity theories: spontaneous and induced vectorization,” Phys. Rev. D99, 044038 (2019), arXiv:1901.02461 [gr-qc] .
  * Kase  _et al._ (2020) R. Kase, M. Minamitsuji, and S. Tsujikawa, “Neutron stars with a generalized Proca hair and spontaneous vectorization,” Phys. Rev. D102, 024067 (2020), arXiv:2001.10701 [gr-qc] .
  * Ramazanoğlu (2019) F. M. Ramazanoğlu, “Spontaneous tensorization from curvature coupling and beyond,” Phys. Rev. D99, 084015 (2019), arXiv:1901.10009 [gr-qc] .
  * Ramazanoğlu (2018) F. M. Ramazanoğlu, “Spontaneous growth of spinor fields in gravity,” Phys. Rev. D98, 044011 (2018), [erratum: Phys. Rev.D100,no.2,029903(2019)], arXiv:1804.00594 [gr-qc] .
  * Barausse  _et al._ (2013) E. Barausse, C. Palenzuela, M. Ponce, and L. Lehner, “Neutron-star mergers in scalar-tensor theories of gravity,” Phys. Rev. D87, 081506 (2013), arXiv:1212.5053 [gr-qc] .
  * Palenzuela  _et al._ (2014) C. Palenzuela, E. Barausse, M. Ponce, and L. Lehner, “Dynamical scalarization of neutron stars in scalar-tensor gravity theories,” Phys. Rev. D89, 044024 (2014), arXiv:1310.4481 [gr-qc] .
  * Shibata  _et al._ (2014) M. Shibata, K. Taniguchi, H. Okawa, and A. Buonanno, “Coalescence of binary neutron stars in a scalar-tensor theory of gravity,” Phys. Rev. D89, 084005 (2014), arXiv:1310.0627 [gr-qc] .
  * Sampson  _et al._ (2014) L. Sampson, N. Yunes, N. Cornish, M. Ponce, E. Barausse, A. Klein, C. Palenzuela, and L. Lehner, “Projected Constraints on Scalarization with Gravitational Waves from Neutron Star Binaries,” Phys. Rev. D90, 124091 (2014), arXiv:1407.7038 [gr-qc] .
  * Sennett and Buonanno (2016) N. Sennett and A. Buonanno, “Modeling dynamical scalarization with a resummed post-Newtonian expansion,” Phys. Rev. D93, 124004 (2016), arXiv:1603.03300 [gr-qc] .
  * Khalil  _et al._ (2019) M. Khalil, N. Sennett, J. Steinhoff, and A. Buonanno, “Theory-agnostic framework for dynamical scalarization of compact binaries,” Phys. Rev. D100, 124013 (2019), arXiv:1906.08161 [gr-qc] .
  * Damour and Nordtvedt (1993a) T. Damour and K. Nordtvedt, “General relativity as a cosmological attractor of tensor scalar theories,” Phys. Rev. Lett. 70, 2217–2219 (1993a).
  * Damour and Nordtvedt (1993b) T. Damour and K. Nordtvedt, “Tensor - scalar cosmological models and their relaxation toward general relativity,” Phys. Rev. D48, 3436–3450 (1993b).
  * Anderson  _et al._ (2016) D. Anderson, N. Yunes, and E. Barausse, “Effect of cosmological evolution on Solar System constraints and on the scalarization of neutron stars in massless scalar-tensor theories,” Phys. Rev. D94, 104064 (2016), arXiv:1607.08888 [gr-qc] .
  * Anderson and Yunes (2017) D. Anderson and N. Yunes, “Solar System constraints on massless scalar-tensor gravity with positive coupling constant upon cosmological evolution of the scalar field,” Phys. Rev. D96, 064037 (2017), arXiv:1705.06351 [gr-qc] .
  * Mendes (2015) R. F. P. Mendes, “Possibility of setting a new constraint to scalar-tensor theories,” Phys. Rev. D91, 064024 (2015), arXiv:1412.6789 [gr-qc] .
  * Mendes and Ortiz (2016) R. F. P. Mendes and N. Ortiz, “Highly compact neutron stars in scalar-tensor theories of gravity: Spontaneous scalarization versus gravitational collapse,” Phys. Rev. D93, 124035 (2016), arXiv:1604.04175 [gr-qc] .
  * Antoniou  _et al._ (2021) G. Antoniou, L. Bordin, and T. P. Sotiriou, “Compact object scalarization with general relativity as a cosmic attractor,” Phys. Rev. D 103, 024012 (2021), arXiv:2004.14985 [gr-qc] .