11institutetext: Hector O. Silva 22institutetext: Max Planck Institute for
Gravitational Physics (Albert Einstein Institute),
D-14476 Potsdam, Germany
22email<EMAIL_ADDRESS>
# Neutron stars as extreme gravity probes

Hector O. Silva

###### Abstract

Neutron stars are powerful probes into the extremes of physics. In this
chapter, we will discuss how observations of neutron stars, either in
isolation or in binaries, can be leveraged to test general relativity and
constrain competing theories of gravity.111Note added in the arXiv version:
this text reflects the subject’s status as of January 16th, 2024.

##  1 Introduction

###  1.1 Neutron stars, gravity and relativity

Neutron stars are outstanding laboratories for studying physics at its
extreme, particularly gravity. Let us do a few order-of-magnitude estimates to
understand why this the case. A typical neutron star has a mass of $M_{\rm
NS}=1.4~{}{\rm M}_{\odot}$, where ${\rm M}_{\odot}$ is one solar mass ${\rm
M}_{\odot}\approx 1.99\times 10^{30}$ kg, and a radius $R$ of approximately
$12$ km. We can quantify how relativistic a system is by defining the
compactness $C$, the dimensionless ratio between the system’s characteristic
mass and size

$C=\frac{GM}{Rc^{2}}\,,$ (1)

where $c=299,792,458$ m s-1 is the speed of light in vacuum and $G\approx
6.67\times 10^{-11}$ ${\rm m}^{3}~{}{\rm kg}^{-1}~{}{\rm s}^{-2}$ is Newton’s
gravitational constant. For a neutron star, we have $C_{\textrm{NS}}\approx
0.2$. To give a sense of what this number means, the Sun’s compactness is a
minuscule $C_{\odot}\approx 10^{-6}$, while for a black hole
$C_{\textrm{BH}}=0.5$, where we took $R$ to be hole’s Schwarzschild radius
$R_{\textrm{S,\,{\rm BH}}}=2GM/c^{2}$. The Schwarzschild radius of a neutron
star is approximately,

$R_{\textrm{S,\,NS}}\approx 4.1~{}{\rm km}\,.$ (2)

Comparing this value to their typical $12$ km radius, we find that neutron
stars are at the verge of becoming black holes. To have a measure of how
strong the gravitational field of a neutron star is, we can use the
(Newtonian) escape velocity, $v_{\rm e}$,

$v_{\rm e}/c=\sqrt{\frac{2GM}{Rc^{2}}}=\sqrt{2C}\approx 0.6\,.$ (3)

That is, the velocity required for a test mass, leaving from the surface of a
neutron star, to escape to infinity is more than half the speed of light. For
all these reasons, black holes, neutron stars (and white dwarfs) are
collectively named “compact objects.”

These order-of-magnitude estimates show that we cannot make use the
nonrelativistic gravity of Newton to model a neutron star. Instead, we must
use a _relativistic theory of gravity_ Oppenheimer:1939ne . Einstein’s general
relativity is the standard theory of relativistic gravity. In part, this
privileged position has been achieved by general relativity’s outstanding
success in describing astrophysical phenomena in agreement with observations.
Examples include, in increasing order of how relativistic the system is, the
advance of the periastron of Mercury’s orbit around the Sun, the secular
decrease of the orbital period of binary pulsars and the generation and
propagation of gravitational waves produced by coalescing compact binaries.

###  1.2 Motivations for testing gravity

In light of these (and many other) successes, one can ask in all fairness “
_why do we need to test general relativity further?_ ” One may answer this
question as follows.

First, it is important to scrutinize general relativity against experiments as
it stands as our best theory to describe one of the fundamental interactions
of nature. As emphasized by Will Will:2014kxa , the rigidity of general
relativity, in the sense that it has no tunable parameters, makes any new
experimental test (or any refinement of an existing test) a potentially fatal
one. In the event this ever happens, it would necessarily mean the existence
of new physics and force us to reevaluate our understanding of gravity.

Second, from a theoretical perspective, general relativity is an effective
theory of the gravitational interaction that is found to agree with
observations at the energy (curvature) scales that we have probed so far. It
is generally believed that general relativity is corrected at some higher
energy scale. These corrections can be constructed from a “bottom-up”
approach, by systematically adding to the theory all terms that conform to
some physical principles, or from a “top-down” approach, in which such terms
would be determined from the low-energy limit of a more fundamental theory,
e.g., string theory. We can then wonder whether these corrections have any
observational signatures, for instance, at the energy scales we can probe with
astronomical observations.

Third, there are known examples of theories of gravity which at certain energy
scales are also compatible with observations, whereas at other energy scales
predict deviations of order unity with respect to the predictions of general
relativity. Hence, whenever a new energy scale becomes experimentally
accessible, it is sensible to be open to the possibility that there may be new
gravitational effects that have so far remained “screened” to our
observations, and are only revealed at this new scale.

As we have seen, neutron stars are, second only to black holes, the strongest
sources of gravity we know of. As such, they are natural laboratories to
examine the validity of general relativity in its strong field regime. Quite
remarkably, the earliest work on neutron stars in modified gravity is, to our
knowledge, from 1967 Salmona:1967zz . This is the same year in which Bell
Burnell discovered the first pulsar, PSR B1919$+$21 Hewish:1968bj  and
Shklovsky correctly interpreted the x-ray source Sco-X1 as being an accreting
neutron star Shklovsky:1967 . We can say, half whimsically, that _the study of
neutron stars in modified theories of gravity is as old as the discovery of
neutron stars themselves!_ Quoting Andersson Andersson:2017kru , neutron stars
allows us to hear “whispers from the edge of physics.” We want to hear, if
any, whispers from new physics beyond general relativity.

###  1.3 Goals and the structure of this chapter

In the remainder of this chapter we will discuss some examples of new
phenomena related to neutron stars that take place when we substitute general
relativity with another gravity theory. It is not our intent to provide a
comprehensive review of the topic, a goal that Refs. Berti:2015itd ;
Doneva:2017jop ; Psaltis:2008bb ; Olmo:2019flu  already accomplish. In
addition, despite their importance in testing gravity, binary pulsars will not
be covered here. The reader can find accounts of this subject in Ref.
Wex:2014nva  and in Chapter 3 in this volume by Hu et al. Hu:2023vsq .

The approach we adopt here is instead a pedagogical one, where we will go in
depth in some specific examples. We will assume some familiarity with
Newtonian gravity and general relativity at the level of the book by Poisson
and Will PW:book . For the novice, we hope that this chapter can provide a
gentle introduction to the topic and serve as a starting point for further
reading or, even, inspire original research problems. For the seasoned
practitioner, we hope that this chapter contains one or two new pieces of
information, or that a topic familiar to her is presented under a different
light.

This chapter is organized as follows. In Sec. 2 we will review how we obtain
neutron star solutions in general relativity. This will set the stage for us
to present how we obtain neutron star solutions in scalar-tensor gravity, an
instructive example of a modified gravity theory. We will also glance over the
properties of neutron stars in a few other gravity theories. In Sec. 3 we
discuss the challenges in testing gravity with neutron star observations and
discuss some ideas that have been proposed to circumvent these problems.
Section 4 we review the status of x-ray emission from rotating neutron stars
(pulse profile modeling) in modified gravity. This a promising, but still not
fully explored avenue to perform tests of gravity. In Sec. 5, we review a
small selection of tests of gravity using the observations of gravitational
waves signals produced by compact binaries containing at least one neutron
star. Section 6 summarizes this chapter. We use geometrical units $c=G=1$
unless stated otherwise, and employ the mostly plus metric signature.

##  2 The equilibrium configuration of a nonrotating neutron star

We will start by presenting what modifications to general relativity can cause
on the structure of isolated nonrotating neutron stars. We begin with a
refresher of how we model neutron stars in general relativity, and then move
to extensions thereof.

###  2.1 General relativity

General relativity is described by the action,

$S=\frac{1}{16\pi}\int{\rm d}^{4}x\sqrt{-g}\,R+S_{\rm m}[\Psi_{\rm
m};\,g_{\mu\nu}]\,,$ (4)

where $g_{\mu\nu}$ is the metric (with determinant $g$), $R$ the Ricci scalar,
and $\Psi_{\rm m}$ are the matter fields. Matter fields are coupled
universally to $g_{\mu\nu}$ thus ensuring the weak equivalence principle, that
has been tested to exquisite precision (see, e.g., Ref. Will:2014kxa .) The
field equations of the theory, obtained by varying the action with respect to
the inverse metric $g^{\,\mu\nu}$ (defined such that
$g^{\,\mu\alpha}\,g_{\alpha\nu}=\delta^{\,\mu}{}_{\nu}$) are,

$R_{\mu\nu}-\tfrac{1}{2}R\,g_{\mu\nu}=8\pi T_{\mu\nu}\,,$ (5)

where $T_{\mu\nu}$ is the energy-momentum tensor of the matter fields. This
tensor is expressed as a variation of $S_{\rm m}$ with respect to
$g^{\,\mu\nu}$ as follows

$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta S_{\rm m}}{\delta
g^{\,\mu\nu}}\,.$ (6)

The Bianchi identity,

$\nabla^{\mu}\left(R_{\mu\nu}-\tfrac{1}{2}R\,g_{\mu\nu}\right)=0\,,$ (7)

ensures that $T_{\mu\nu}$ is conserved

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

Now say that we want to use these field equations to model a nonrotating star.
We will need two ingredients for this. The first one is a choice of the matter
content of the star, while the second one is a choice of the spacetime metric
that will be associated to this matter distribution. In isolation, the star
will be static and spherically symmetric, and we can use the following line
element as an ansatz

${\rm d}s^{2}=g_{\mu\nu}\,{\rm d}x^{\,\mu}\,{\rm d}x^{\nu}=-e^{2\nu(r)}\,{\rm
d}t^{2}+e^{2\lambda(r)}\,{\rm d}r^{2}+r^{2}\,{\rm
d}\theta^{2}+r^{2}\sin^{2}\theta\,{\rm d}\varphi^{2}\,.$ (9)

Here, $\nu$ and $\lambda$ are functions of $r$, called the areal radial
coordinate. To understand this nomenclature, we can take a constant-in-time
hypersurface (i.e., ${\rm d}t=0$) and look at a fixed value of $r=R$ (i.e.,
${\rm d}r=0$). We then see from the line element (9) that $r$ measures the
radius of a 2-sphere. For reasons that will become clear in a moment, it is
convenient to write

$\lambda=-\tfrac{1}{2}\log[1-2m(r)/r]\,,$ (10)

where $m$ is called the mass function. Finally, we will call $\nu$ the lapse
function.222The reason for this nomenclature is that when we do a $3+1$
decomposition of the spacetime metric, the $tt$ component is
$\alpha^{2}-\beta_{k}\,\beta^{k}$, where $\alpha$ is the lapse function and
$\beta_{k}$ the shift vector. For $\beta_{k}=0$ (called normal or Eulerian
coordinates), we see that $\nu$ is related to the lapse function $\alpha$
Gourgoulhon:2007ue . As for matter, we will assume that the star can be
described by a perfect fluid. The energy-momentum tensor for a perfect fluid
with pressure $p$ and total energy density $\varepsilon$ is

$T_{\mu\nu}=(\varepsilon+p)u_{\mu}u_{\nu}-p\,g_{\mu\nu}\,,$ (11)

where $u^{\,\mu}$ is the four-velocity of an observer comoving with a fluid
element. Because $u^{\,\mu}$ is a timelike vector, it satisfies the
normalization $u^{\,\mu}u_{\mu}=-1$, which combined with Eq. (9) yields

$u^{\,\mu}=\left\\{-\exp(-\nu),\,0,\,0,\,0\right\\}.$ (12)

The pressure $p$ and energy density $\varepsilon$ are defined with respect to
this observer. The energy density can be decomposed as
$\varepsilon=\rho\,(1+\Pi)\,,$ where $\rho$ is the baryonic mass density and
$\Pi=(\varepsilon-\rho)/\rho$ is the internal energy per unit baryonic mass.

We can now insert Eqs. (9)-(11) in the field equation (5) and energy-momentum
conservation equation (8) to obtain the Tolman-Oppenheimer-Volkoff (TOV)
equations Tolman:1939jz ; Oppenheimer:1939ne :

$\displaystyle\frac{{\rm d}m}{{\rm d}r}$ $\displaystyle=4\pi
r^{2}\varepsilon\,,$ (13a) $\displaystyle\frac{{\rm d}\nu}{{\rm d}r}$
$\displaystyle=\frac{1}{r^{2}}\frac{m+4\pi r^{3}p}{1-2m/r}\,,$ (13b)
$\displaystyle\frac{{\rm d}p}{{\rm d}r}$
$\displaystyle=-(\varepsilon+p)\frac{{\rm d}\nu}{{\rm d}r}\,.$ (13c)

We see that we have four variables to solve for $(m,\nu,\varepsilon,p)$, but
only three equations. To close the system, we need to provide an equation of
state, that is, a relationship between the pressure $p$ and energy density
$\varepsilon$

$p=p(\varepsilon)\,.$ (14)

A simple neutron-star equation of state is a polytropic
$p=\kappa\,\varepsilon^{1+1/n}$, where $n=0.5$ and $\kappa=2\times 10^{5}$ km4
can be used as representative values Kokkotas:2000up .333See Tooper
Tooper:1964 ; Tooper:1966  for a detailed discussion relativistic polytropic
equations and for the general-relativistic version of the Lane-Emden equations
that describe a polytrope: a self-gravititation body obeying a polytropic
equation of state in Newtonian gravity Chandrasekhar:1939isss.book . Equations
of state often come in tabulated form, and are derived from nuclear physics
calculations that employ a variety of methods and physical assumptions. See,
e.g., Refs. Baym:2017whm ; Lattimer:2021emm  for reviews articles, and Refs.
Typel:2013rza ; Oertel:2016bki ; CompOSECoreTeam:2022ddl  for an equation-of-
state database. An equation of state of the form (14) is said to be
barotropic. In general, $p$ could also depend on the temperature $T$, entropy
$s$, etc. However, old neutron stars have temperature well below the Fermi
temperature meaning that in practice they are “cold” and can be assumed to be
at zero temperature.444 In the literature one often finds the terms “soft” and
“stiff” to describe the equation of state. This terminology stems from how
fast the pressure increases as a function of the energy density. The faster
this increase is, the harder it is to compress the fluid and the stiffer the
equation of state is said to be. Hence, stiffer equations of state provide
more pressure against gravitational attraction and can support stars with
larger masses. Conversely, softer equations result in stars with smaller
radii.

The star’s radius $R$ is the location in which the pressure $p$ vanishes. The
star’s gravitational mass then follows from Eq. (13a)

$M=m(R)=\int_{0}^{\,R}4\pi r^{2}\varepsilon\,{\rm d}r\,,$ (15)

which justifies the naming $m(r)$ as the mass function. The upper integration
limit follows from $\varepsilon$ being zero in vacuum; see Eq. (13a).

It is also useful to define the star’s baryonic mass $M_{\rm b}$ as an
integral of the baryonic mass density $\rho$ in the spatial volume occupied by
the star,

$M_{\rm b}=4\pi m_{\rm b}\int_{0}^{\,R}\frac{n_{\rm
b}r^{2}}{\sqrt{1-2m/r}}\,{\rm d}r\,.$ (16)

Here, we wrote $\rho=m_{\rm b}n_{\rm b}$, where $m_{\rm b}\approx 1.66\times
10^{-24}$ g is the atomic mass unit and $n_{\rm b}$ the number density of
baryons. From $M$ and $M_{\rm b}$ we can define the normalized (gravitational)
binding energy $E_{\rm b}$,

$E_{\rm b}=(M_{\rm b}-M)/M\,,$ (17)

which is positive for a bound system as a star. Typical values of $E_{\rm b}$
for neutron stars range between 0.1 and 0.2, approximately.

An important result for us is the Jebsen-Birkhoff theorem. This theorem states
that the _vacuum exterior_ spacetime of any, possibly time dependent,
spherically symmetric mass configuration, is given by the Schwarzschild
metric. Thus, the spacetime outside the star has line element

${\rm d}s^{2}=-\,(1-2M/r)\,{\rm d}t^{2}+(1-2M/r)^{-1}\,{\rm
d}r^{2}+r^{2}\,{\rm d}\theta^{2}+r^{2}\sin^{2}\theta\,{\rm d}\varphi^{2}\,,$
(18)

where $M=m(R)$. If one is interested only in studying the exterior spacetime
of the neutron star, one can also work directly with Eq. (18), using values of
$M$ and $R$ that are typical of neutron stars.

We now have all the ingredients to obtain neutron star solutions in general
relativity. To solve the TOV equations we would in principle start from the
star’s center $r=0$ and proceed outwards until we locate the star’s radius. In
practice, however, we find that the TOV equations are singular at $r=0$.
Instead,we use an approximate, power-series solution of the relevant variables
around $r=0$. For instance, the first terms in the series expansion are:

$\displaystyle m$ $\displaystyle=\tfrac{4}{3}\pi\varepsilon_{\rm c}r^{3}+{\cal
O}(r^{5})\,,$ (19a) $\displaystyle p$ $\displaystyle=p_{\rm
c}-\tfrac{2}{3}\pi(\varepsilon_{\rm c}+p_{\rm c})(\varepsilon_{\rm c}+3p_{\rm
c})r^{2}+{\cal O}(r^{4})\,,$ (19b) $\displaystyle\nu$ $\displaystyle=\nu_{\rm
c}+\tfrac{2}{3}\pi(\varepsilon_{\rm c}+3p_{\rm c})r^{2}+{\cal O}(r^{4})\,,$
(19c)

where $\nu_{\rm c}$ is a free constant. The integration of the TOV equations
(13) then proceed as follows:

Box 1: TOV integration in general relativity

  1. 1.

Choose a value of the energy density $\varepsilon_{{\rm c}}$ at the star’s
center $r=r_{{\rm c}}$. The central pressure $p_{{\rm c}}=p(\varepsilon_{{\rm
c}})$ is fixed by the equation of state (14).

  2. 2.

Integrate the TOV equations (13) outwards starting from $r=r_{\rm c}$, with
initial conditions for $m$, $p$ and $\nu$ given by Eqs. (19) evaluated at
$r=r_{\rm c}$.

  3. 3.

Proceed with the integration until the pressure $p$ vanishes. This fixes the
radius $R$ and the mass $M=m(R)$.

  4. 4.

Fix the constant $\nu_{\rm c}$, such that $\nu$ matches smoothly
$\exp(2\nu(R))=1-2M/R$ at $r=R$.

If we repeat these steps for a range of values of $\varepsilon_{\rm c}$ we
obtain a one-parameter sequence of solutions. These solutions can be plotted
in a $MR$-plane, producing what is known as the mass-radius relation.

In Fig. 1, we show several of such curves for a selection of equations of
state that support neutron stars with masses $M\geqslant 2~{}{\rm M}_{\odot}$,
as required by the precise inferred masses of pulsars J0348$+$0432, $M=2.01\pm
0.04~{}{\rm M}_{\odot}$ Antoniadis:2013pzd , and J0740$+$6620, $M=2.08\pm
0.07~{}{\rm M}_{\odot}$ NANOGrav:2019jur ; Fonseca:2021wxt , both at $1\sigma$
level. The equations of state we used form a subset of those referred to in
Ref. Silva:2020acr . The large variability in the predicted masses and radii
is a consequence of our current poor understanding of the equation of state at
high densities.

Figure 1: Neutron star mass-radius relations obtained by integrating the TOV
equation for a large sample of equations of state.

###  2.2 Bergmann-Wagoner scalar-tensor gravity

#### Motivation

In general, the steps to obtain neutron star solutions taken previously can be
repeated in any other gravity theory with some modifications. As a simple, but
nontrivial example, we will consider the Bergmann-Wagoner scalar-tensor
gravity Bergmann:1968ve ; Wagoner:1970vr . This is arguably the simplest
extension to general relativity, and it is, to our knowledge, the earliest
beyond-general-relativity theory in which neutron stars were studied
Salmona:1967zz .

#### Theoretical minimum

The action of this theory is given by,

$S=\frac{1}{16\pi G_{\ast}}\int{\rm
d}^{4}x\sqrt{-g_{\ast}}\,\left[R_{\ast}-2(\partial\varphi)^{2}-4V(\varphi)\right]+S_{\rm
m}\left[\Psi_{\rm m};\,A^{2}(\varphi)\,g_{\ast\mu\nu}\right]\,,$ (20)

where $\varphi$ is scalar field with canonical kinetic term
$(\partial\varphi)^{2}=\partial_{\alpha}\varphi\,\partial^{\alpha}\varphi$,
and self-interaction potential $V$. The spacetime metric is $g_{\ast\mu\nu}$
with respect to which we define the Ricci scalar $R_{\ast}$, as usual. Lastly,
$A$ is a nonvanishing a function of $\varphi$.

Why are we introducing asterisks in the metric and Ricci scalar here? And why
did we reinstate the gravitational constant $G_{\ast}$? [Compare against Eq.
(4).] The reason is that in scalar-tensor theories we have the freedom to
perform field redefinitions and metric rescalings to write the action in
different forms, known as “frames.”

Equation (20) describes the theory in the _Einstein frame_ characterized by
the absence of a coupling between $\varphi$ and the Ricci scalar $R_{\ast}$.
The main difference between Eq. (20) and the action describing a scalar field
minimally coupled to gravity is in the universal coupling between all matter
fields $\Psi_{\rm m}$ with the metric

$\tilde{g}_{\mu\nu}=A^{2}(\varphi)\,g_{\ast\mu\nu}\,.$ (21)

Hence, the theory obeys the weak equivalence principle, and massive test
particles follow geodesics in $\tilde{g}_{\mu\nu}$. Massless particles, such
as photons, follow null geodesics which are invariant under the transformation
(21), sometimes called a conformal transformation. We call
$\tilde{g}_{\mu\nu}$ the _Jordan frame_ metric. The reason for this
nomenclature is that we can rewrite Eq. (20) replacing $g_{\ast}$ and
$R_{\ast}$ by their counterparts calculated with respect to
$\tilde{g}_{\mu\nu}$. We use tildes and asterisks to denote quantities in the
Jordan and Einstein frames respectively. The advantage of working in the
Einstein frame is that the field equations are considerably simpler.

In the following, we will assume that the potential $V$ is zero, i.e., the
scalar field is massless and non-self-interacting. Performing the variation of
the action (20) with respect to $g^{\,\mu\nu}_{\ast}$ and $\varphi$ gives us
the field equations,

$\displaystyle
R_{\ast\mu\nu}-\tfrac{1}{2}R_{\ast}\,g_{\ast\mu\nu}=2\partial_{\mu}\varphi\,\partial_{\nu}\varphi-
g_{\ast\mu\nu}(\partial\varphi)^{2}+8\pi T_{\ast}\,,$ (22a)
$\displaystyle\Box_{\ast}\varphi=-4\pi\alpha T_{\ast}\,,$ (22b)

where $\Box_{\ast}=g_{\ast\mu\nu}\nabla_{\ast}^{\mu}\nabla_{\ast}^{\nu}$ is
the d’Alembertian operator, $T_{\ast}=T_{\ast}^{\mu\nu}g_{\ast\mu\nu}$ is the
trace of the energy-momentum tensor, which is defined in analogy with Eq. (6),

$T_{\ast\mu\nu}=-\frac{2}{\sqrt{-g_{\ast}}}\frac{\delta S_{\rm m}}{\delta
g_{\ast}^{\,\mu\nu}}\,,$ (23)

and

$\alpha(\varphi)=\frac{{\rm d}}{{\rm d}\varphi}\log A(\varphi)\,.$ (24)

We call $\alpha$ the scalar-field-matter coupling as it controls how strongly
matter sources $\varphi$ in Eq. (22b). We have also set $G_{\ast}=1$.

We define the Jordan frame energy-momentum tensor $\tilde{T}^{\mu\nu}$ in
analogy with Eq. (6),

$\tilde{T}_{\mu\nu}=-\frac{2}{\sqrt{-\tilde{g}}}\frac{\delta S_{\rm
m}}{\delta\tilde{g}^{\,\mu\nu}}\,.$ (25)

From their definitions, we see that the Einstein and Jordan frames energy-
momentum tensors satisfy the following relations

$T_{\ast}^{\mu\nu}=A^{6}(\varphi)\tilde{T}^{\mu\nu},\quad
T_{\ast\mu\nu}=A^{2}(\varphi)\tilde{T}_{\mu\nu},\quad\textrm{and}\quad
T_{\ast}=A^{4}(\varphi)\tilde{T}\,,$ (26)

and their divergences are

$\tilde{\nabla}^{\mu}\tilde{T}_{\mu\nu}=0\,,\quad\textrm{and}\quad\nabla_{\ast}^{\mu}T_{\ast\mu\nu}=\alpha(\varphi)\,T_{\ast}\,\nabla_{\ast\nu}\varphi\,.$
(27)

As already mentioned, this shows that test particles follow geodesics in the
Jordan, but not in the Einstein frame.

The physically measured pressure and densities of a fluid are those in the
Jordan frame to which all matter (including our measurement devices) are
universally coupled to. In analogy with Eq. (11), we write for a perfect
fluid:

$T_{\ast\mu\nu}=(\varepsilon_{\ast}+p_{\ast})u_{\ast\mu}u_{\ast\nu}-p_{\ast}\,g_{\ast\mu\nu}\,.$
(28)

Taking the trace, we obtain,

$T_{\ast}=T_{\ast\mu\nu}\,g_{\ast}^{\mu\nu}=3p_{\ast}-\varepsilon_{\ast}=A^{4}(\varphi)\,(3\tilde{p}-\tilde{\varepsilon})\,,$
(29)

which is nonzero in general.555Note that $T_{\ast}=\tilde{T}=0$ if
$\tilde{p}=\tilde{\varepsilon}/3$. This corresponds to conformal fluids and
they are interesting in the context of the gravity-fluid correspondence
Rangamani:2009xk . Conformal fluids have a sound speed squared $v_{\rm
s}^{2}={{\rm d}p/{\rm d}\varepsilon}=1/3$ that was thought to be an upper
bound on the sound speed of any medium. This bound has been challenged by the
observation of approximately $2~{}{\rm M}_{\odot}$ neutron stars. See, for
instance, Refs. Alsing:2017bbc ; Bedaque:2014sqa  for a discussion. Thus,
stars can act as scalar field sources; see Eq. (22b).

The limit of general relativity (plus a minimally coupled scalar field) is
obtained when we take $\alpha$ to be zero. In this limit, the Jordan- and
Einstein-frame metrics coincide and we can safely drop their identifiers.

#### The equations of hydrostatic equilibrium

The generalization of the TOV equations can be obtained in the same way as in
general relativity. We assume the line element:

${\rm d}s_{\ast}^{2}=g_{\ast\mu\nu}\,{\rm d}x^{\,\mu}\,{\rm
d}x^{\nu}=-e^{2\nu(r)}\,{\rm d}t^{2}+e^{2\lambda(r)}\,{\rm d}r^{2}+r^{2}\,{\rm
d}\theta^{2}+r^{2}\sin^{2}\theta\,{\rm d}\varphi^{2}\,,$ (30)

and a Jordan-frame perfect fluid described by

$\tilde{T}_{\mu\nu}=(\tilde{\varepsilon}+\tilde{p})\tilde{u}_{\mu}\tilde{u}_{\nu}-\tilde{p}\,\tilde{g}_{\mu\nu}\,,$
(31)

where $\tilde{u}_{\mu}\tilde{u}^{\,\mu}=-1$ and hence
$\tilde{u}^{\,\mu}=A^{-1}(\varphi)\,\\{-\exp(-\nu),\,0,\,0,\,0\\}$.

Using Eqs. (30) and (11), we find from the field equations (22) the TOV
equations in scalar-tensor gravity

$\displaystyle\frac{{\rm d}m}{{\rm d}r}$ $\displaystyle=4\pi
r^{2}A^{4}\tilde{\varepsilon}+\frac{1}{2}r(r-2m)\psi^{2},$ (32a)
$\displaystyle\frac{{\rm d}\nu}{{\rm d}r}$ $\displaystyle=4\pi
r^{2}A^{4}\frac{\tilde{p}}{r-2m}+\frac{1}{2}r\psi^{2}+\frac{m}{r(r-2m)},$
(32b) $\displaystyle\frac{{\rm d}\tilde{p}}{{\rm d}r}$
$\displaystyle=-(\tilde{p}+\tilde{\varepsilon})\left(\frac{{\rm d}\nu}{{\rm
d}r}+\alpha\psi\right),$ (32c) $\displaystyle\frac{{\rm d}\psi}{{\rm d}r}$
$\displaystyle=4\pi\frac{rA^{4}}{r-2m}[\alpha(\tilde{\varepsilon}-3\tilde{p})+r\psi(\tilde{\varepsilon}-\tilde{p})]-\frac{2(r-m)}{r(r-2m)}\psi,$
(32d)

where we defined $\psi={\rm d}\varphi/{\rm d}r$, and Eq. (32d) follows from
(22b).

The main difference with respect to general relativity is the presence of the
scalar field, and this will force us to introduce some modifications to the
integration procedure to obtain neutron star solutions and determine their
properties. To understand why this is the case, we evaluate Eqs. (32a) and
(32d) outside the star by setting $\tilde{\varepsilon}=\tilde{p}=0$:

$\displaystyle\frac{{\rm d}m}{{\rm d}r}$
$\displaystyle=\frac{1}{2}r(r-2m)\psi^{2},$ (33a) $\displaystyle\frac{{\rm
d}\nu}{{\rm d}r}$ $\displaystyle=\frac{1}{2}r\psi^{2}+\frac{m}{r(r-2M)},$
(33b) $\displaystyle\frac{{\rm d}\psi}{{\rm d}r}$
$\displaystyle=-\frac{2(r-m)}{r(r-2m)}\psi.$ (33c)

We see that unless the scalar field is constant, $\psi=0$, the mass function
$m$ has a nonzero gradient in vacuum and the lapse behaves differently with
respect to general relativity. These facts have two important consequences.
The first is that the mass function evaluated at $r=R$ does not determine the
star’s total mass anymore. The second is that the spacetime in the star’s
exterior is no longer Schwarzschild.

To determine the star’s mass, we can expand Eqs. (33) in a power series in
$1/r$ and take the far-field limit $R/r\ll 1$. The result is

$\displaystyle m$ $\displaystyle=M-q^{2}/(2r)-q^{2}M/(2r^{2})+O(r^{-3})\,,$
(34a) $\displaystyle e^{2\nu}$ $\displaystyle=1-2M/r+O(r^{-3})\,,$ (34b)
$\displaystyle\varphi$ $\displaystyle=\varphi_{0}+q/r+qM/r^{2}+O(r^{-3})\,,$
(34c)

where we have identified the $1/r$ coefficient in the lapse’s expansion as the
mass $M$, and $\varphi_{0}$ and $q$ are constants. The constant $\varphi_{0}$
can be interpreted as the ambient scalar field in which the star is embedded,
while $q$, in analogy with electrostatics, is called the (monopole) scalar
charge. This charge is not a conserved charge as in electromagnetism. For a
constant scalar field the charge vanishes and then $\varphi=\varphi_{0}$,
$m=M$ and $\exp(2\nu)=1-2M/r$. Hence, we recover the Schwarzschild exterior
solution.

Equations (34) suggest that the vacuum spacetime of a neutron star in scalar-
tensor gravity cannot be obtained in closed form, contrary to general
relativity. It turns out that this is possible if we adopt a different
coordinate system, named after Just Buchdahl:1959nk ; Coquereaux:1990qs ;
Damour:1992we ; Just:1959 . The line element is:

${\rm d}s_{\ast}^{2}=-f^{b/a}\,{\rm d}t^{2}+f^{-b/a}\,{\rm
d}\varrho^{2}+\varrho^{2}f^{1-b/a}\,{\rm
d}\theta^{2}+\varrho^{2}f^{1-b/a}\,\sin^{2}\theta\,{\rm d}\phi^{2}\,,$ (35)

where $b=2M$ and $f=1-a/\varrho$. The constant $a$ has units of length and is
related to the star’s mass $M$ and the scalar charge $q$, as we will see
shortly. The Just and Schwarzschild radial coordinates are related as

$r=\varrho\,(1-a/\varrho)^{(1-b/a)/2}\,,$ (36)

which cannot be inverted analytically. However, $r\simeq\varrho$ when
$a/\varrho\ll 1$, for instance, far away from the star.

Using Eq. (35), we can solve analytically Eq. (22b) in the exterior region of
a star:

$\varphi=\varphi_{0}+(q/a)\log(1-a/\varrho)\,.$ (37)

In the large $\varrho$-limit,
$\varphi\simeq\varphi_{0}+q/\varrho\simeq\varphi_{0}+q/r$, confirming that $q$
is the scalar charge. The constants $a$, $b$ and $q$ are not independent, and
satisfy the constraint

$a^{2}-b^{2}-(2q)^{2}=0\,,$ (38)

that follows from the $\theta\theta$ component of Eq. (22a) in vacuum. From
Eq. (38), we can express the ratio $a/b$ appearing in Eq. (35) as

$a/b=\sqrt{1+Q^{2}}\,,$ (39)

where we introduced the dimensionless scalar-charge-to-mass ratio

$Q=q/M\,.$ (40)

Therefore, the exterior metric of a neutron star in scalar-tensor theory is
fully determined by $M$ and $Q$. The constant $Q$ controls how much the
spacetime deviates from the Schwarzschild spacetime, which is recovered for
$Q=0$.

The benefit of having this analytical solution is that we can now integrate
the TOV equations (32) only within the star, and determine $M$,
$\varphi_{\infty}$ and $q$ by imposing that the variables match smoothly Eqs.
(35) and (37) at the star’s radius. This makes the problem of finding neutron
star solutions similar to that in general relativity. This calculation was
first done in Ref. Damour:1993hw  (see also Ref. Creci:2023cfx ) and yields:

$\displaystyle q$ $\displaystyle=-M\psi_{\rm s}/\nu^{\prime}_{\rm s}\,,$ (41a)
$\displaystyle\varphi_{0}$ $\displaystyle=\varphi_{\rm s}-(q/M)\nu_{\rm
s}=\varphi_{{\rm s}}+\frac{\psi_{\rm s}}{(\nu^{\prime 2}_{\rm s}+\psi_{\rm
s}^{2})^{1/2}}{\rm arctanh}\left[\frac{(\nu^{\prime 2}_{\rm s}+\psi_{\rm
s}^{2})^{1/2}}{\nu^{\prime}_{\rm s}+1/R}\right]\,,$ (41b) $\displaystyle M$
$\displaystyle={R^{2}\nu^{\prime}_{\rm s}}\left(1-\frac{2m_{\rm
s}}{R}\right)^{1/2}\exp\left\\{-\frac{\nu^{\prime}_{\rm s}}{(\nu^{\prime
2}_{\rm s}+\psi_{\rm s}^{2})^{1/2}}{\rm arctanh}\left[\frac{(\nu^{\prime
2}_{\rm s}+\psi_{\rm s}^{2})^{1/2}}{\nu^{\prime}_{\rm
s}+1/R}\right]\right\\}\,,$ (41c)

where a prime indicates a derivative with respect to $r$ and the subscript “s”
means that a variable is evaluated at the star’s surface $r=R$.

We now have all the ingredients to obtain neutron star solutions in scalar-
tensor gravity. To solve the TOV equations (32) we proceed as follows: Box 2:
TOV integration in scalar-tensor gravity given a conformal factor $A(\varphi)$

  1. 1.

Choose a value of the energy density $\tilde{\varepsilon}_{{\rm c}}$ at the
star’s center $r=r_{{\rm c}}$. The central pressure $\tilde{p}_{{\rm
c}}=\tilde{p}(\tilde{\varepsilon}_{{\rm c}})$ is fixed by the equation of
state (14).

  2. 2.

Integrate the TOV equations (32) outwards starting from $r=r_{\rm c}$, with
initial conditions $\varphi_{\rm c}=\varphi(r_{\rm c})$, $\psi_{\rm c}=0$,
$\nu_{\rm c}=\nu(r_{\rm c})$ (a freely chosen constant), and $m_{\rm c}=0$ for
the scalar field, its derivative, lapse and mass functions.

  3. 3.

Proceed with the integration until the pressure $p$ vanishes. This fixes the
radius $R$. From the values of $\nu^{\prime}_{\rm s}$, $\psi_{\rm s}$, and
$m_{\rm s}$ determine $q$, $\varphi_{0}$ and $M$ using Eqs. (41).

  4. 4.

Fix the constant $\nu_{\rm c}$, such that $\nu$ matches smoothly $\nu_{\rm s}$
defined in the second equality in Eq. (41b).

Often one requires that $\varphi_{0}$ has a specific value fixed e.g., either
by observational constraints or by the scalar field’s cosmological evolution.
In this case, $\varphi_{\rm c}$ _is not a free constant_. Instead, we must
apply a root-finding algorithm to determine $\varphi_{\rm c}$ such that the
outcome of steps 1 to 4 is $\varphi_{0}$ with the targeted value.

We stress that one can also numerically integrate the TOV equations (33) in
vacuum up to a large $r=r_{\rm max}$. We then fix the three constants
$\varphi_{0}$, $q$ and $M$, from the knowledge of $m$, $\varphi$ and $\psi$
using the three equations (34) evaluated at $r_{\rm max}$. The advantage of
this approach is that it generalizes to other theories of gravity in which an
exact exterior spacetime solution does not exist. Of course, one has to first
obtain the equivalent of Eqs. (34) in the theory of interest. This method
gives the same results as the approach presented in Box 2.

Finally, we emphasize that $q$, $M$ and $R$ are computed in the Einstein
frame. From ${\rm d}\tilde{s}^{2}=A^{2}(\varphi)\,{\rm d}s_{\ast}^{2}$, it is
easy to see that the Jordan-frame radius $\tilde{R}$ is:

$\tilde{R}=A^{2}(\varphi_{\rm s})\,R\,,$ (42)

whereas a translation from Einstein to Jordan frame of Eqs. (34) shows that
the Jordan-frame gravitational mass $\tilde{M}$ and scalar charge $\tilde{q}$
are:

$\tilde{M}=M\,\frac{1}{A_{0}}\,\frac{1-Q\,\alpha_{0}}{1+\alpha^{2}_{0}}\,,\quad\tilde{q}=q\,\frac{\alpha_{0}}{A^{3}_{0}}\frac{2}{1+\alpha^{2}_{0}}\,,$
(43)

respectively, and $A_{0}=A(\varphi_{0})$ and $\alpha_{0}=\alpha(\varphi_{0})$.
Details of this calculation can be found in Refs. Creci:2023cfx ;
Damour:1992we ; Pani:2014jra . In scalar-tensor theories, it is the Einstein-
frame mass that has desired physical properties. See Ref. Doneva:2022ewd ,
Sec. III.A and references therein.

#### Spontaneous scalarization

The attentive reader may at this point start becoming anxious about the fact
that the function $A(\varphi)$ has been left undetermined so far. What
functional form do we chose for it? Let us Taylor-expand $\alpha(\varphi)$
[recall Eq. (24)] as follows

$\alpha(\varphi)=\alpha_{0}+\beta_{0}(\varphi-\varphi_{0})+\dots\,,\quad\beta_{0}=({\rm
d}\alpha/{\rm d}\varphi)_{\varphi=\varphi_{0}}\,,$ (44)

which we truncate at linear order in $\varphi-\varphi_{0}$. One can show that
the deviations from general relativity in the post-Newtonian approximation are
proportional to the square of $\alpha_{0}$. This term is related to the
constant $\omega_{\rm BD}$ in Brans-Dicke theory Brans:1961sx , and is
constrained by Solar System test to be smaller than $3.5\cdot 10^{-3}$. For
this reason, it is tempting to imagine that higher-order powers in this series
will contribute very little to the theory’s prediction. This is not the case,
however, as first observed by Damour and Esposito-Farèse Damour:1993hw ;
Damour:1996ke .

Consider for the moment that the value of the scalar field is so small, that
we can treat it as a small perturbation $\delta\varphi$ around a general-
relativistic background. Let us write $\varphi=\varphi_{0}+\delta\varphi$, set
without loss of generality $\varphi_{0}=0$, and $\alpha_{0}=0$. By linearizing
Eq. (22b) in $\delta\varphi$, we find

$\Box_{\ast}\delta\varphi=-4\pi\beta_{0}\tilde{T}\delta\varphi\,,$ (45)

where the conformal factor $A\approx 1$ in our approximation and hence
$\tilde{T}\simeq T_{\ast}=T$. This equation has the same form of a massive
Klein-Gordon equation,

$\left[\Box_{\ast}-\mu^{2}_{\rm
eff}(r)\right]\delta\varphi=0\,,\quad\mu^{2}_{\rm eff}=-4\pi\beta_{0}T\,,$
(46)

where $\mu_{\rm eff}$ is a position-dependent effective mass for the scalar
field.

There are two important observations we can make from this equation. First,
$\delta\varphi=0$ is a solution. Second, as long as $\beta_{0}T>0$, the
perturbations can have an imaginary effective mass. This means that some modes
of $\delta\varphi$ can become tachyonic unstable, causing $\delta\varphi$ to
exponentially grow in time; see Ref. Doneva:2022ewd , Sec. II, for a review.
In general, the trace $T$ is negative, meaning that we need $\beta_{0}$
negative for this to happen.

To find the threshold for this instability, we can first exploit the spherical
symmetry of the background spacetime and do a harmonic mode decomposition as
follows:

$\delta\varphi(t,r,\theta,\varphi)=\sum_{l,\,m}\frac{\psi(r)}{r}Y_{lm}(\theta,\varphi)\,{\rm
e}^{{\rm i}\sigma t}\,,$ (47)

where $Y_{lm}$ are the spherical harmonics and $\sigma$ is a frequency.
Substituting Eq. (47) in Eq. (46) and performing separation of variables, we
obtain an ordinary differential equation for $\psi$. This equation can be
brought to the form of a Schrödinger-like equation

$\left[\frac{{\rm d}^{2}}{{\rm d}x^{2}}+\sigma^{2}-V_{\rm
eff}(x)\right]\,\psi=0\,,$ (48)

by introducing a new radial coordinate coordinate $x$, related to $r$ as

$\frac{{\rm d}x}{{\rm d}r}=(1-2m/r)^{-1/2}\,{\rm e}^{-\nu}\,.$ (49)

Here, $\nu$ and $m$ are the lapse and mass functions of the background
spacetime, respectively. In Eq. (48), we identified an effective potential
$V_{\rm eff}$ for $\psi$, namely,

$V_{\rm eff}={\rm
e}^{2\nu}\left[\frac{l(l+1)}{r^{2}}+\frac{2m}{r^{3}}+4\pi(p-\varepsilon)+\mu^{2}_{\rm
eff}\right]\,.$ (50)

We can now identify the threshold for the instability by searching for
$\sigma=0$ (zero mode or bound state) solutions of Eq. (48). Bound state
solutions are characterized by being regular at the center of the star ($r=0$)
and for approaching a constant (in our case chosen to be zero) at spatial
infinity.666The black-hole perturbation theory enthusiast will note that in
vacuum, the effective potential (50) reduces to that of massless scalar
perturbations in a Schwarzschild spacetime, as expected. See, for instance,
Chapter 9, Sec 9.2.2, in this volume Franchini:2023eda . This constitutes a
boundary-value problem, where given a spacetime background and a value of $l$,
we try to find a value of $\beta_{0}$ for which a solution $\psi$ that
satisfies these boundary conditions exists. In practice, we can solve this
boundary-value problem using a shooting method.

The upshot of this exercise (see Ref. Silva:2014fca , Sec. 4.3., for details)
is that bound states (for $l=0$) exist when $\beta_{0}\lesssim-4.35$, quite
independently of the equation of state. These bound states mark the
bifurcation point from the general-relativistic sequence of solutions of a new
branch of neutron stars endowed with a nontrivial scalar field. This is known
as spontaneous scalarization.777 The terminology is inspired by the phenomenon
of spontaneous magnetization of ferromagnets Damour:1996ke . See also Ref.
Harada:1998ge  for an analysis of scalarization within catastrophe theory.
Figure 2 shows the mass-radius relation of neutron stars in this theory. These
curves were obtained by solving the TOV equations (32) with the
$\alpha=\beta_{0}\varphi$ and asymptotic boundary condition $\varphi_{0}=0$,
using a shooting method as outlined in Box 2.

Figure 2:  Neutron star scalarization. We show a sequence of neutron stars
solutions modeled using, as an example, the APR equation of state Akmal:1998cf
, for $\alpha_{0}=10^{-5}$ and $\beta_{0}\in[-6,0]$ in intervals of
$\Delta\beta_{0}=0.1$. The left panel shows the mass-radius relations, while
the right panel shows the charge-to-mass ratio $Q$ as a function to central
density. When $\beta_{0}\lesssim-4.5$, spontaneous scalarization takes place,
causing the stars to have a nonzero $Q$ and masses and radii that can deviate
considerably from general relativity.

The fact that the theory admits multiple branches of solutions begs the
question: which of them is favored? To analyze this question we can look at
the binding energy $E_{\rm b}$, defined in Eq. (17). One finds that stars in
the scalarized branch have larger binding energy and are energetically
favorable over the nonscalarized solutions Damour:1993hw ; Harada:1997mr ;
Horbatsch:2010hj . This has been confirmed in numerical simulations by Novak
Novak:1998rk , who studied the dynamical evolution of neutron stars from the
nonscalarized to the scalarized branch.

The presence of the scalar charge causes the neutron star, when placed in a
binary, to source scalar dipole radiation. This causes the orbital motion of
the binary to evolve differently if were to take into account gravitational
radiation only. Binary pulsar observations have severely constrained the
existence of scalar dipole radiation and have essentially ruled out
spontaneous scalarization in the model proposed by Damour and Esposito-Farèse
Kramer:2021jcw ; Zhao:2022vig

To close our brief discussion on scalarization, let us revisit two statements
we made here. First, we set $\alpha_{0}=0$. If we had not done this, it is
easy to see from Eq. (22b) that $\varphi=0$ is no longer a solution of the
theory in the presence of matter. This means neutron stars (and any other
material object) in Brans-Dicke theory always have some scalar field.
Ultimately, this is what makes it possible to constraint this theory with
Solar System observations. Scalarization, in contrast, requires strongly self-
gravitating objects. Second, we said that in general $\tilde{T}<0$ and for
this reason scalarization happens when $\beta_{0}<0$. It turns out that for
some equations of state, $\tilde{T}$ can become positive at sufficiently large
densities. Hence, scalarization can also happen for positive values of
$\beta_{0}$ in these cases. See Refs. Mendes:2014ufa ; Mendes:2016fby ;
Palenzuela:2015ima ; Podkowka:2018gib  for details.

#### The black-hole limit

Let us close our discussion by contrasting neutron stars with black holes in
scalar-tensor gravity. Equations (35) and (37) remain valid in a black-hole
spacetime. If Eq. (35) describes a black hole, the event horizon
$\varrho_{{\rm h}}$ will be located at $\varrho_{{\rm h}}=a$. However, from
Eq. (37), the scalar field $\varphi$, and all its derivatives, become singular
at the horizon, unless $q=0$. As we saw earlier, for $q=0$ we obtain the
Schwarzschild metric plus a trivial scalar field.

This result is nothing but a statement of the no-hair theorems in scalar-
tensor gravity. It teaches us an important lesson: _to test gravity, sometimes
it is better to rely on neutron stars rather than on black holes._ See Ref.
Herdeiro:2015waa  for a review on no-hair theorems.

###  2.3 Remarks about other gravity theories

We saw how to construct nonrotating neutron star solutions in general
relativity and in scalar-tensor gravity. In the latter, we found that neutron
stars can be significantly different with respect to their general-
relativistic counterparts through an effect known as spontaneous
scalarization. What can happen in other gravity theories? Before moving on,
let us highlight the salient properties of neutron stars in a few other
gravity theories.

#### Massive scalar-tensor gravity

In this theory, one chooses the potential appearing in Eq. (20) to have the
form $V(\varphi)=m_{\varphi}^{2}\varphi^{2}/(2\hbar^{2})$, giving the scalar
field a mass $m_{\varphi}$. This simple modification is sufficient to make the
Damour-Esposito-Farèse model consistent with binary pulsar observations. The
reason can be understood qualitatively as follows. The inclusion of the mass
term modifies the spatial behavior of the scalar field outside the star
causing it to decay exponentially as
$\varphi\sim(1/r)\,\exp[-(m_{\varphi}/\hbar)r]$, in the far-field limit. This
strongly suppresses the scalar interaction between the binary components when
the orbital separation $R_{\rm orb}$ is larger than the Compton wavelength

$\lambda_{\varphi}=2\pi\hbar/(m_{\varphi}c)\,,$ (51)

associated to the scalar field. In this case, the emission of scalar dipole
radiation is suppressed Alsing:2011er . Then, to evade current binary pulsar
constraints on the emission of scalar dipole radiation one can chose
$m_{\varphi}\gg 10^{-16}~{}{\rm eV}$, and $3\lesssim-\beta\lesssim 10^{3}$ if
one wants neutron stars, but not white dwarfs, to scalarize
Ramazanoglu:2016kul . See also the related “asymmetron” theory proposed in
Refs. Chen:2015zmx ; Morisaki:2017nit .

Regarding scalarization, the main difference with respect to the massless
case, is that the scalar field mass raises the threshold for scalarization,
for a fixed value of $\beta_{0}$. Qualitatively, we can imagine that by giving
the scalar field a mass we would require more energy to excite its
perturbations, hence making scalarization harder to happen.

Nonetheless, this is compensated by the wider range of value of $\beta_{0}$
compatible with observations that one can work with. A detailed investigation
of neutron star solutions was done in Ref. Rosca-Mead:2020bzt . This work also
discusses numerical methods to solve the TOV equations using a relaxation
scheme.

#### Einstein-scalar-Gauss-Bonnet gravity

This theory is described by the action

$S=\frac{1}{16\pi}\int{\rm
d}^{4}x\sqrt{-g}\,\left[R-2(\partial\varphi)^{2}+\ell^{2}_{\rm
GB}\,f(\varphi)\,R^{2}_{\rm GB}\right]+S_{\rm m}[\Psi_{\rm m};g_{\mu\nu}]\,,$
(52)

where $\ell_{\rm GB}$ is a lengthscale. The theory has as main feature the
coupling between the scalar field $\varphi$ with the higher-order curvature
scalar, known as Gauss-Bonnet invariant:

$R^{2}_{\rm
GB}=R^{2}-4R_{\alpha\beta}R^{\alpha\beta}+R_{\alpha\beta\gamma\delta}R^{\alpha\beta\gamma\delta}\,.$
(53)

In four-dimensional spacetimes, $R_{\rm GB}^{2}$ is total derivative and
therefore cannot contribute to the field equations. However, the coupling of
$R_{\rm GB}^{2}$ to $\varphi$, causes the former to contribute nontrivially to
the theory’s field equations.

The theory described by the action (52) has a few interesting properties that
makes it an attractive toy model to explore beyond-general-relativity physics.
Firstly, the field equations of the theory are at most second-order in time
and therefore the theory is free from negative-energy propagating degrees of
freedom known as ghosts. (This is known as the “Ostrogradski instability”;
see, e.g., Refs. Kobayashi:2019hrl ; Woodard:2015zca .) In fact, Einstein-
scalar-Gauss-Bonnet gravity, is a subclass of the broader class of Horndeski
theories of gravity Horndeski:1974wa , the most general scalar-tensor theory
in four dimensions with this property; see Ref. Kobayashi:2019hrl  for a
review. Secondly, the theory introduces a single free parameter, the
lengthscale $\ell_{\rm GB}$. Thirdly, the coupling between $\varphi$ and
$R_{\rm GB}^{2}$ is present in low-energy limits of string theory; see e.g.,
Refs. Cano:2021rey ; Metsaev:1987zx .

The action (52) yields the following field equation for the scalar field

$\Box\varphi=-\tfrac{1}{4}\ell_{\rm GB}^{2}\,f_{,\varphi}(\varphi)\,R_{\rm
GB}^{2}\,,$ (54)

where $f_{,\varphi}(\varphi)={\rm d}f/{\rm d}\varphi$. We can identify two
broad subclasses of Einstein-scalar-Gauss-Bonnet gravity according to the form
of the coupling function $f(\varphi)$. One subclass, call it Type I, has
$f_{,\varphi}(\varphi)\neq 0$ for all $\varphi$. Examples, include the
dilatonic $f\sim\exp(\varphi)$ and the shift-symmetric $f\sim\varphi$
couplings.888The name dilatonic arises from the context of string theory,
where $\varphi$ is called a “dilaton field.” The name “shift symmetric” is due
to the fact that for this coupling function the action is invariant under a
shift $\varphi\to\varphi+{\rm const.}$. Another subclass, say Type II, admits
$f_{,\varphi}(\varphi_{0})=0$ for some constant $\varphi_{0}$. Examples
include the quadratic $f\sim\varphi^{2}$ and Gaussian
$f\sim\exp(-\varphi^{2})$ couplings.

Let us first briefly discuss neutron star solutions in Type I theories. Yagi
et al. Yagi:2015oca  showed that whether or not neutron stars have monopole
scalar charge in these theories depends on the functional form of $f$. In
particular, they showed that neutron stars in the shift-symmetric theory _do
not_ posses monopole scalar charge, whereas in the dilatonic case they do.
This is interesting for it identifies a subclass of theories in which scalar
dipole radiation emission would vanish in binary pulsars. In contrast, black
holes in the shift-symmetric theory _do have_ monopole scalar charge
Yunes:2009ch . More generally, the interaction term $\varphi\,R_{\rm GB}^{2}$
is necessary for black holes to have monopole scalar charge in any shift-
symmetric Horndeski theory Sotiriou:2013qea ; Sotiriou:2014pfa ;
Saravani:2019xwx . This is means that to test this theory, a binary with at
least one black hole is desirable over a double neutron star system. Compare
with our discussion in Sec. 2.2.

We now turn to the Type II theories. Their most important feature is that they
admit spontaneous scalarization Silva:2017uqg ; Doneva:2017jop . We can
understand why this happens by analogy with the equations we met in Sec. 2.2.
More precisely, Eq. (45) can be compared against Eq. (54), with the latter
particularized for the coupling function $f=\eta~{}\varphi^{2}/2$ and likewise
linearized around a neutron-star background in general relativity,

$\Box\delta\varphi=-\tfrac{1}{4}\eta~{}\ell_{\rm GB}^{2}R_{\rm
GB}^{2}\delta\varphi\,,$ (55)

where $\eta$ is a dimensionless constant. We can then make the identifications

$\eta\sim 16\pi\beta_{0},\quad\textrm{and}\quad\ell_{\rm GB}^{2}~{}R_{\rm
GB}^{2}\sim\tilde{T}\,,$ (56)

establishing the analogy between the two theories. Despite the complexity of
Eq. (53), the Gauss-Bonnet invariant for a neutron-star spacetime in general
relativity is simple Silva:2017uqg

$R_{\rm GB}^{2}=\frac{48m^{2}}{r^{6}}-\frac{128\pi(m+2\pi
r^{3}p)\varepsilon}{r^{3}}\,.$ (57)

We can use the near-center expansion formulas (19), to conclude that $R_{\rm
GB}^{2}$ is negative at the star’s core. We also see that $R_{\rm GB}^{2}$
matches smoothly the Schwarzschild result $R_{\rm GB}^{2}=48M^{2}/r^{6}$ at
the star’s surface $r=R$, where $m(r)=M$, as we should expect from the Jebsen-
Birkhoff theorem. The fact that $R^{2}_{\rm GB}$ is not positive-definite
suggests that scalarization can happen for both signs of the coupling constant
$\eta$ Silva:2017uqg . When $\eta$ is negative, scalarization is triggered
deep within the star, while for $\eta$ positive, scalarization is triggered
near the star’s surface. Scalarized neutron stars exhibit the familiar
bifurcation off the nonscalarized sequence of solutions Doneva:2017duq . Stars
in the scalarized branch have smaller mass (at fixed central density) with
respect to nonscalarized solution, and are energetically favored over the
latter. Finally, the scalarized solutions have a monopole scalar charge
Doneva:2017duq .

The most striking difference between scalarization in the Damour-Esposito-
Farèse model and in Type II Einstein-scalar-Gauss-Bonnet gravity is that in
the latter scalarization also occurs for black holes Doneva:2017jop ;
Silva:2017uqg . This follows from $R^{2}_{\rm GB}$ being nonzero for vacuum
black-hole spacetimes whereas $\tilde{T}$ is zero. We refer the reader to Ref.
Doneva:2022ewd , Sec. IV, for a review of black-hole scalarization.

#### Dynamical Chern-Simons gravity

This theory is described by the action Jackiw:2003pm

$S=\frac{1}{16\pi}\int{\rm
d}^{4}x\sqrt{-g}\,\left[R-2(\partial\vartheta)^{2}+\ell^{2}_{\rm
CS}\,\vartheta\,{}^{\ast}RR\right]+S_{\rm m}[\Psi_{\rm m};g_{\mu\nu}]\,,$ (58)

where $\ell_{\rm CS}$ is a lengthscale. The theory has as main feature the
coupling between the pseudoscalar field $\vartheta$ with the higher-order
curvature scalar, known as the Pontryagin density

${}^{\ast}RR=R_{\nu\mu\rho\sigma}\,{}^{\ast}R^{\,\mu\nu\rho\sigma}\,,$ (59)

where
${}^{\ast}R^{\mu\nu\rho\sigma}=\varepsilon^{\,\rho\sigma\alpha\beta}R^{\mu\nu}{}_{\alpha\beta}/2$
is the dual of the Riemann tensor and $\varepsilon^{\,\rho\sigma\alpha\beta}$
is the totally antisymmetric Levi-Civita tensor; see Ref. Alexander:2009tp
for a review. As with the case of the Gauss-Bonnet invariant, ${}^{\ast}RR$ is
a total derivative in four dimensions and it cannot contribute to the field
equations unless we couple it to $\vartheta$. Contrary to Einstein-scalar-
Gauss-Bonnet gravity, dynamical Chern-Simons theory has field equations which
are third-order in time. For this reason, the theory should be seen as an
effective theory, with $\ell_{\rm CS}/\ell\ll 1$, where $\ell$ is some
characteristic lengthscale associate to the problem; see Ref. Delsate:2014hba
. For neutron stars, one imposes, for example, that $\ell_{\rm CS}/M\ll 1$.

An interesting aspect of this theory, is that the Pontryagin density is zero
for any spherically symmetric spacetime Jackiw:2003pm . For this reason,
nonrotating neutron stars, which is what we have been considering so far, _are
identical to those of general relativity._ Deviations from general relativity
only appear when the neutron-star spacetime is no longer spherically
symmetric. The simplest example of such a spacetime is when one gives a small
degree of rotation to star. In this case, the line element (18) is augmented
by a mixed ${\rm d}t\,{\rm d}\phi$ term as follows:

${\rm d}s^{2}=-e^{2\nu(r)}\,{\rm d}t^{2}+e^{2\lambda(r)}\,{\rm
d}r^{2}+r^{2}\,{\rm d}\theta^{2}+r^{2}\sin^{2}\theta\,{\rm
d}\varphi^{2}+2\omega(r,\theta)r^{2}\sin^{2}\theta\,{\rm d}t\,{\rm d}\phi\,,$
(60)

where $\omega$ is a function that represents the angular velocity acquired by
a particle falling from infinity as measured by a static asymptotic observer
(“frame dragging effect”) Hartle:1967he . The scalar field sourced by the
slowly-rotating star is not spherical symmetric, and has instead a dipolar
structure Yunes:2009ch . Asymptotically, the scalar field decays as

$\vartheta\simeq\mu\cos\theta/r^{2}+{\cal O}(r^{-3})\,,$ (61)

where $\mu$ is a scalar dipole charge, whose magnitude is proportional to
$\ell_{\rm CS}/M$ and the angular momentum of the star. This is different from
the other theories we encountered so far. See Ref. Yagi:2013mbt  for details.

An important consequence of these results is that dynamical Chern-Simons
gravity is poorly constrained by binary pulsar observations. This is because
the binary components are slowly rotating999 Of course, pulsars are fast
rotating in comparison to any other astronomical object we know of. Here we
mean “slow” with respect to to the mass-shedding frequency $f_{\rm
ms}=(2\pi)^{-1}~{}\sqrt{GM/R_{\rm eq}}$, at which the angular velocity of the
star equals that of a particle in Keplerian orbit at the star’s equator. For a
representative millisecond pulsar with $M=1.4~{}{\rm M}_{\odot}$, $R_{\rm
eq}=12$ km and spinning at $f=700$ Hz, we have $f/f_{\rm ms}\approx 0.4$. This
justifies the use of a slow-rotation expansion to model rotating neutron
stars, as pioneered by Hartle Hartle:1967he . Equation (60) is the leading
order in this expansion. See also Refs. Hartle:1968si ; Chandrasekhar:1974
and Benhar:2005gi  for details on the expansion to second and third orders in
spin, respectively. The subject of _rapidly_ rotating neutron stars is
presented in the monographs Friedman:2013xza ; Gourgoulhon:2010ju . and the
dipole scalar charge is proportional to the spin angular momentum through the
quadratic-in-curvature terms Yagi:2013mbt . Nevertheless, the strongest
constraint to date on this theory comes from neutron star observations
Silva:2020acr . We will explain how this is possible in Sec. 3.2.

##  3 Tests of gravity with isolated neutron stars

Now that we have seen in a few examples how neutron stars are different in
theories beyond general relativity, we will proceed to discuss the
observational implications of these results. Some, such the emission of scalar
radiation in a binary, were already mentioned. First, let us clarify what we
mean by the statement of “testing a theory of gravity with neutron star
observations”.

###  3.1 Digression: what do we mean by testing gravity?

In general, observations provides us with some data $d$ and we use a model (or
hypothesis) ${\cal H}$ to attempt to explain this data. The model ${\cal H}$
would usually assume general relativity combined with other physics (say
hydromagnetics) that are relevant to produce an accurate description of the
observed data $d$. Any statistical inference from $d$ using ${\cal H}$ (say,
an estimate of a neutron star mass) will come with an associated error, with a
statistical and a systematic components. Both error components are always
present. The former because our observations are subject to noise sources and
the latter because our models are not a perfect description of nature.
Qualitatively, the statistical error gives the width of the posterior
distribution of a fiducial inferred parameter, while the systematic error
tells how much the peak of the posterior distribution (i.e., the most likely
value of the parameter) disagrees from the true value of the parameter. We are
fine as long as our systematic errors are smaller than our statistical errors.
That is because our knowledge about the value of this fiducial parameter will
still overlap with its true value (or accepted value) to some confidence
level. However, this will no longer be the case if the systematic errors
become dominant. If this happens, our inferences become biased.

It is often the case that we make statistical inferences assuming the validity
of general relativity. We then fiddle around with the free parameters of some
alternative theory and identify how large they can be before getting in
tension with the error in our parameter inference. This gives an upper bound
on the parameters of the modified theory. Many of the results that will quote
should be interpreted in this way.

If nature happens to not be described by general relativity, this would be a
source of systematic error in any analysis that uses a model that assumes
general relativity to be true. From our previous discussion, this may lead to
a fundamental bias (using the language of Ref. Yunes:2009ke ) in our parameter
inference. This nomenclature is used to distinguish this type of bias from
other biases originating, for instance, from some simplified or overlooked
physics that went into ${\cal H}$.

In reality, of course, when we collect data with our telescopes we have no
idea whether the photons hitting the detector were produced in a process in
which general relativity or an alternative holds. A solution to this problem
is to come up with a model that is to some extend theory agnostic. We say “to
some extend” because one still relies on some physical assumptions (i.e.,
symmetry requirements or the hypothesis that we deviate a little from general-
relativistic predictions) to develop such model. This is the approach adopted
in the parametrized post-Newtonian (to test gravity at first post-Newtonian
order with Solar System tests) Will:1993hxu , the parametrized post-Keplerian
Deruelle:1986 ; Damour:1992  (to test gravity with binary pulsars), and the
parametrized post-Einsteinian formalism Yunes:2009ke ; Chatziioannou:2012rf
and parametrized effective-one-body waveform models Brito:2018rfr ;
Ghosh:2021mrv ; Maggio:2022hre  (to test gravity with gravitational waves), to
name a few. All these theory-agnostic models (and others that we will see in
Secs. 3.2 and 3.2) have free parameters which have specific values in general
relativity. We can then include them into our parameter inference. If their
posterior distributions agree with their predicted values in general
relativity we have confirmed Einstein’s theory, if not, it signifies the
presence of some systematical effect caused either by the true value of the
parameter not being the one predicted by general relativity (exciting) or by
some extra physics missed out by our formalism (often boring.) It often
happens that these constraints on theory-agnostic parameters can be mapped in
theory-specific constraints. One should also have in mind that sometimes
missing nongravitational physics may emulate beyond-general-relativity
effects, and potentially lead to false violations of general relativity.

Finally, one could develop also two models, ${\cal H}_{1}$ and ${\cal H}_{2}$,
which include the same physics and differing only on the theory of gravity
assumed. We will see one such example in Sec. 4. We could then perform
Bayesian model selection between the two models and see which one is favored
by the data. Theory-agnostic and theory-specific approaches are complementary
to one another and pursuing both is a healthy approach in testing gravity.

It is important to keep in mind that sometimes scientific progress allows us
to take big leaps across different regimes of a theory. This was certainly the
case with the first detection of gravitational waves produced by the
coalescence of a black-hole binary LIGOScientific:2016aoc . Until 2015, we had
never explored with observations the regime of gravity in which two strongly
self-gravitating bodies move at very high speeds ($v/c\approx 0.5$ at the late
inspiral) and the gravitational field is strong and rapidly varying in time. A
similar leap happened with the discovery of binary pulsars, which opened a new
frontier in tests of gravity and expanded our testing ground beyond the Solar
System. When faced with these situations, it makes sense to be more cautious
about the hypothesis used in modeling our new exciting data, the theory of
gravity included.

###  3.2 Difficulties in testing theories of gravity with neutron star
observations

Figure 3: Illustration of the gravity-theory degeneracy. We show the mass-
radius relation for neutron stars described by the APR equation of state
Akmal:1998cf , for a few theories of gravity, including some not covered in
Sec. 3 and values of their free constants, if any. The theories are: general
relativity (“GR,” solid line), dynamical Chern-Simon gravity (“dCS,”
coincident solid line with general relativity), Eddington-inspired Born-Infeld
(“EiBI,” dashed line, $\kappa=\pm\,5\cdot 10^{-3}$) Pani:2011mg , Einstein-
scalar-Gauss-Bonnet (“ESGB,” dilatonic coupling, dot-dashed line,
$\alpha=20~{}{\rm M}_{\odot}$) Pani:2011xm , and Einstein-Aether (“EA,” dotted
line, $c_{14}=0.3$) Yagi:2013ava . All curves terminate at the maximum-mass
neutron star. We see that even if knew the correct equation of state, it would
be hard to distinguish these theories with neutron star masses and radii
observations alone.

In principle, precise measurements of masses and radii of neutron stars can be
used to constrain modifications to general relativity. Unfortunately, things
are not that simple. See recent efforts in Refs. Demirboga:2023ktt ;
Tuna:2022qqr . As we mentioned already, there are still large uncertainties on
the equation of state of neutron stars, particularly at densities a few times
larger than the nuclear saturation density $\varepsilon_{\rm n}\simeq
2.67\times 10^{14}$ g cm-3, roughly above which the neutron-star equation of
state starts becoming poorly known.

To makes matters worse, even if we knew the correct equation of state,
different gravity theories predict similar deformations of the general-
relativistic mass-radius relation. The only exception we are aware of are
theories which allow for spontaneous scalarization. Figure 3 illustrates this
discussion. In summary, we must face:

  1. 1.

_The equation-of-state-gravity degeneracy:_ the predictions of the bulk
properties (mass, radius, moment of inertia, etc.) of a neutron star in a
theory $\mathscr{T}$ belonging to the theory space $\\{\mathscr{T}_{i}\\}$ are
degenerate with general-relativistic predictions due to uncertainties in the
equation of state.

  2. 2.

_The gravity-theory degeneracy:_ even if we knew the equation of state
exactly, the prediction of the bulk properties of neutron stars for the
theories in the set $\\{\mathscr{T}_{i}\\}$ are degenerate with one another.

We will discuss ways to address these two problems separately, although, of
course, both could be analyzed simultaneously.

#### The equation-of-state-gravity degeneracy, quasiuniversal relations and a
parametrized I-Love test of gravity

One way to break this degeneracy is by making use of relations among neutron
star observables which depend weakly on the equation of state. We call these
“quasiuniversal relations” and a number of them are known to exist in general
relativity, as reviewed in Ref. Yagi:2016bkt . For example, in 1989, Lattimer
and Yahil Lattimer:1989zz  proposed that the neutron star binding energy (17)
is given by

$E_{\rm b}\approx 1.5\cdot 10^{51}\,(M/{\rm M}_{\odot})^{2}~{}{\rm
ergs}=0.084~{}{\rm M}_{\odot}\,(M/{\rm M}_{\odot})^{2}\,,$ (62)

which is valid to about 20%. The “quasiuniversal” character of Eq. (62) over a
wide range of neutron star models can be improved by including the dependence
on their radii. This can be accomplished using the compactness $C$. In
particular, Lattimer and Prakash Lattimer:2000nx  suggested the relation

$E_{\rm b}/M=\frac{3}{5}\,\frac{C}{1-C/2}\,.$ (63)

In the same work, Lattimer and Prakash, also obtained formulas relating the
dimensionless moment of inertia $I/(MR^{2})$ as a function of the compactness
$C$ by working with some exact solutions of the TOV equations in general
relativity. For instance, they derived from what is known as the Tolman VII
solution Tolman:1939jz  the relation,

$I/(MR^{2})=\tfrac{2}{7}\,(1-1.1~{}C-0.6~{}C^{2})^{-1}\,,$ (64)

which yields a fair approximation for $I/(MR^{2})$ when $C\gtrsim 0.1$.101010
The Tolman VII solution is an analytical solution of the TOV equations that
models fairly well neutron stars. Recent works have generalized the Tolman VII
solution with the inclusion of new parameters which are then fitted against
realistic equation-of-state neutron star models, see Refs. Jiang:2019vmf ;
Posada:2022lij . This solution was used by Ref. Yagi:2021loe  (following Ref.
Horbatsch:2010hj ) to derive, analytically, scalar charges in the Damour-
Esposito-Farèse scalarization model discussed in Sec. 2.2.  See also the
earlier work of Ravenhall and Pethick Ravenhall:1994ApJ . Other early
quasiuniversal relations relate the oscillation and damping timescales of
pulsation modes of neutron stars to their masses and radii Andersson:1997rn .
The motivation for this work is the possibility of doing neutron star
seismology. That is, infer the mass and radius of a neutron star, and
consequently learn about its internal structure, by measuring of one or more
of its pulsation modes. See Ref. Andersson:2021qdq  for a review.

It was later observed that the moment of inertia $I$, spin-induced quadrupole
moment $Q$, and tidal Love number $\lambda$ are remarkably equation-of-state
independent when made dimensionless by appropriate powers of the mass $M$ and
of the dimensionless spin parameter $\chi=2\pi f\bar{I}$, namely

$\bar{I}=I/M^{3},\quad\bar{Q}=-Q/(\chi^{2}M^{3}),\quad\bar{\lambda}=\lambda/M^{5}.$
(65)

The functions $y(x)$ among any two of these three quantities are known by the
playful name “I-Love-Q” relations Yagi:2013bca ; Yagi:2013awa . For example,
in its original version, $\bar{I}$ and $\bar{\lambda}$ are related as:

$\ln\bar{I}=1.47+0.0817\ln\bar{\lambda}+0.0149\ln\bar{\lambda}^{2}+2.87\cdot
10^{-4}\ln\bar{\lambda}^{3}-3.64\cdot 10^{-5}\ln\bar{\lambda}^{4}\,.$ (66)

Similar relations exist between each of these three quantities and the
compactness $C$, but to a smaller degree of equation-of-state independence.
See, e.g., Refs. Breu:2016ufb ; Maselli:2013mva .

Figure 4: The equation-of-state independence of the moment of inertia. In the
left panel we show the dimensionless moment of inertia $\bar{I}=I/M^{3}$ as a
function of the mass $M$ for several equations of state. In the right panel we
show $\bar{I}$ now as a function of the compactness $C$. Observe how the
variability of $\bar{I}$ on the equation of state is dramatically reduced when
plotting it as a function of $C$. The red dashed line represents the fit
proposed in Ref. Breu:2016ufb , valid in the slow-rotation approximation.

We illustrate our discussion in Figs. 4 and 5. Figure 4 shows what we mean by
equation-of-state independence. In the left panel we show $\bar{I}$ as a
function of the mass $M$, where each curve corresponds to a given equation of
state; the same ones used in Fig. 1. In the right panel we plot the same
curves but now as a functions of the compactness $C$. The reduction of the
variability with respect to the equation of state is clearly visible and
justifies the terminology “quasiuniversal.” The dashed line corresponds to the
“I-C” relation obtained in Ref. Breu:2016ufb , where the moment of inertia was
calculated using the slow-rotation expansion. The maximum absolute relative
errors, $\delta=|1-y_{\rm fit}/y_{\rm data}|$, between data and the fit is
about $\sim 10\%$.

Figure 5: The equation-of-state independence between dimensionless moment of
inertia and tidal Love number. In the top panel we show $\bar{I}$ and
$\bar{\lambda}$ calculated for for several equations of state (barely visible
bundle of solid lines) and the fit (66). The absolute relative error between
the data and fit below $\lesssim 1\%$.

In Fig. 5 we show the “I-Love” relation between $\bar{I}$ and $\bar{\lambda}$.
Although not apparent, the top panel contains the same number of curves as in
Fig. 4, but now the variability between different equations of state is hardly
visible. The dashed line corresponds to Eq. (66). In the bottom panel we show
$\delta$ between $\bar{I}_{\rm data}$ as predicted by each equation of state
with respect to the fit (66). We see that the I-Love relation holds to within
$\sim 1\%$ accros all stellar models and equations of state.

An alternative to Eq. (66) can be worked out as follows. In Newtonian gravity,
we know that the moment of inertia $I\propto MR^{2}$, and hence
$I/M^{3}\propto C^{-2}$, while the Love number $\lambda\propto R^{5}$ and
hence $\lambda/M^{5}\propto C^{-5}$. Then, we can eliminate the compactness
and write $\bar{I}=k\,\bar{\lambda}^{2/5}$. Here, $k\approx 0.52$ is a
constant that depends weakly on the equation of state Yagi:2013awa . This
calculation can be extended systematically in powers of $C$, i.e., in a post-
Minkowskian expansion Chan:2014tva ; Jiang:2020uvb , and yields a powers
series in $C$ for both and $\bar{I}$ and $\bar{\lambda}$. These two equations
can be combined to result in an equation for $\bar{I}(\bar{\lambda})$. The
simple relation

$\bar{I}=\bar{\lambda}^{2/5}\,(c_{0}+c_{1}\bar{\lambda}^{-1/5}+c_{2}\bar{\lambda}^{-2/5})\,,$
(67)

where $c_{0}=0.584$, $c_{1}=0.980$, $c_{2}=2.965$, achieves the same degree of
quasiuniversality as Eq. (66). The prefactor $\bar{\lambda}^{2/5}$ is the
Newtonian result, and each power of $\bar{\lambda}^{-1/5}$ inside the
parenthesis represents a post-Minkowskian correction. The reason is that, as
we saw previously, $C\propto\bar{\lambda}^{-1/5}$, as long as $C$ is small
Silva:2020acr .

Remarkably, quasiuniversal relations also hold in modified theories of
gravity, see Refs. Doneva:2017jop  and Yagi:2016bkt , Sec. 6, for reviews.
However, this does not mean that they have to be described by the same
formulas as in general relativity. This suggests that one can try to
parametrize deviations from the quasiuniversal relations in general relativity
in a theory independent way. This was done in Ref. Silva:2020acr , using the
post-Minkowskian formula (67) as baseline. The parametrized I-Love relation is

$\bar{I}=\bar{\lambda}^{2/5}\,(c_{0}+c_{1}\bar{\lambda}^{-1/5}+c_{2}\bar{\lambda}^{-2/5})+\beta\bar{\lambda}^{-b/5}\,.$
(68)

Here, $\beta$ is a positive real number and $b$ is an integer. These
parameters encode the magnitude and the post-Minkowskian order in which the
modified gravity corrections affect the general-relativistic I-Love relation
(67), respectively. For example, a correction at the Newtonian level
corresponds to $b=-2$.

How can the I-Love relation be used to perform tests of gravity? The first
test we can perform is a null test of general relativity. Assume, for example,
that we have two independent measurements, in general relativity, of
$\bar{I}_{1.4}$ and $\bar{\lambda}_{1.4}$ for a neutron with canonical mass
$M=1.4~{}{\rm M}_{\odot}$, and their respective errors Yagi:2013awa . These
measurements delimit an “error box” in the $\bar{I}$-$\bar{\lambda}$ plane.
For the I-Love relation to be consistent with these measurements it has to
pass through this error box. Another test we can envisage makes use of the
parametrized I-Love relation (68). For a given value of $b$, we search for the
critical value $\beta_{\rm crit}$ of $\beta$ above which the resultant I-Love
curve would be in tension with the earlier mentioned error box Silva:2020acr .

Figure 6: The I-Love test of gravity. The vertical and horizontal lines
delimit the 90% confidence region of $\bar{\lambda}_{1.4}$ and
$\bar{I}_{1.4}$, while the circle marks the median $(190,\,14.6)$. The solid
black line corresponds to Eq. (67) and is consistent with observation at 90%
confidence. The various magenta curves represent the parametrized I-Love
relation (68). The labels correspond to the pair $(b,\,\beta_{\rm crit})$. For
each exponent $b$, the quantity $\beta_{\rm crit}$ represent the critical
value of $\beta$ for which Eq. (68) would be in tension with the measurements
of $\bar{I}_{1.4}$ and $\bar{\lambda}_{1.4}$ at 90% confidence. Figure taken
from Ref. Silva:2020acr .

Figure 6 shows an application of this idea. We use thevalue inferred for the
moment of inertia $\bar{I}_{1.4}=14.6^{+4.5}_{-3.3}$ from x-ray observations
of the pulsar J0737-3039A Silva:2020acr , and of the tidal deformability
$\bar{\lambda}_{1.4}=190^{+390}_{-120}$ from the binary-neutron-star
gravitational-wave signal GW170817 LIGOScientific:2018cki . The values quoted
are at 90% confidence. The tidal deformability is obtained conditional on the
assumptions that the binary components were described by the same equation of
state and were slowly spinning. The circle marks $(14.6,190)$ and the gray box
the respective 90% confidence region. The solid black line is Eq. (67), which
is consistent with the two measurements. Hence, general relativity has passed
this null test. The different magenta lines correspond to post-Minkowskian
exponents $b\in[-2,~{}5]$ and their respective values of $\beta_{\rm crit}$;
both numbers are indicated in parenthesis $(b,~{}\beta_{\rm crit})$ in each
curve. We remark that the constraints for $\beta\leqslant 0$ depends on the
assumptions used to infer $\bar{\lambda}_{1.4}$. In particular, if the tidal
deformabilities were treated as independent free parameters in the waveform
model LIGOScientific:2017vwq , then the $\bar{\lambda}_{1.4}$ posterior would
not have a lower limit. This would make all parametrized I-Love relations with
$b\leqslant 0$ be consistent with the observations.

We can wonder whether this theory-independent constraints can be translated
into constraints on any specific theory. As an example, Ref. Silva:2020acr
considered dynamical Chern-Simons gravity, a theory we met back in Sec. 2.3.
The I-Love relation in this theory is different from that of general
relativity and comparison against a large catalog of neutron stars shows that

$b_{\rm CS}=4\,,\quad\beta_{\rm CS}=6.15\times 10^{-2}\,\bar{\xi}\,,$ (69)

where $\bar{\xi}=4\ell_{\rm CS}^{4}/M^{4}$, where $\ell_{\rm CS}$ is the
lengthscale associated with the theory; cf. Eq. (58). For $b=4$, we see from
Fig. 6 that $\beta_{\rm crit}=8.84\times 10^{2}$, which can be translated into
a bound on $\ell_{\rm CS}$, namely,

$\ell_{\rm CS}\leqslant 16~{}{\rm km}\,,$ (70)

at 90% confidence, if the theory is to be consistent with the observational
bounds on $\bar{I}_{1.4}$ and $\bar{\lambda}_{1.4}$.111111This corresponds to
$\alpha_{\rm CS}^{1/2}\leqslant 8.5$ km in the notation used in Ref.
Silva:2020acr . The two quantities are related by $\ell_{\rm
CS}=(4\pi)^{1/4}\,\alpha_{\rm CS}^{1/2}$. Reference Silva:2020acr  also showed
that this value of $\ell_{\rm CS}$, combined with the estimates on the radius
and compactness of a $1.4~{}{\rm M}_{\odot}$ neutron star derived there,
satisfies the small-coupling approximation as required by the mathematical
consistency of the theory. This is the most stringent constrain on this theory
to date and is seven orders of magnitude stronger than previous Solar System
constraints Ali-Haimoud:2011zme ; Nakamura:2018yaw . We see that
quasiuniversal relations can be a powerful way to test gravity, avoiding
uncertainties regarding the equation of state.

#### The gravity-theory degeneracy and parametrized TOV equations

The observation that different theories of gravity, given an equation of
state, yield mass-radius relations which are in general smooth deformations of
the general-relativity curve suggests that one construct a theory-agnostic
parametrized system of TOV equations.

This was the main motivation for the post-TOV formalism developed in Refs.
Glampedakis:2015sua ; Glampedakis:2016pes ; Silva:2019wfn . The formalism is
valid for static spherical stars and is inspired by the parametrized post-
Newtonian theory, extended to second PN order by adding suitable correction
terms to the fully relativistic TOV equations. In the post-TOV formalism, the
stellar structure equations are given by

$\displaystyle\frac{{\rm d}p}{{\rm d}r}$ $\displaystyle=\left(\frac{{\rm
d}p}{{\rm d}r}\right)_{\rm GR}-\frac{\rho
m}{r^{2}}(\mathscr{P}_{1}+\mathscr{P}_{2})\,,$ (71a) $\displaystyle\frac{{\rm
d}m}{{\rm d}r}$ $\displaystyle=\left(\frac{{\rm d}m}{{\rm d}r}\right)_{\rm
GR}+4\pi r^{2}\rho(\mathscr{M}_{1}+\mathscr{M}_{2})\,,$ (71b)

where $({\rm d}p/{\rm d}r)_{\rm GR}$ and $({\rm d}m/{\rm d}r)_{\rm GR}$ are
given by Eqs. (13) and,

$\displaystyle\mathscr{P}_{1}$
$\displaystyle=\delta_{1}\frac{m}{r}+4\pi\delta_{2}\frac{r^{3}p}{m}\,,$ (72a)
$\displaystyle\mathscr{M}_{1}$
$\displaystyle=\delta_{3}\frac{m}{r}+\delta_{4}\Pi\,,$ (72b)
$\displaystyle\mathscr{P}_{2}$
$\displaystyle=\pi_{1}\frac{m^{3}}{r^{5}\rho}+\pi_{2}\frac{m^{2}}{r^{2}}+\pi_{3}r^{2}p+\pi_{4}\frac{\Pi
p}{\rho}\,,$ (72c) $\displaystyle\mathscr{M}_{2}$
$\displaystyle=\mu_{1}\frac{m^{3}}{r^{5}\rho}+\mu_{2}\frac{m^{2}}{r^{2}}+\mu_{3}r^{2}p+\mu_{4}\frac{\Pi
p}{\rho}+\mu_{5}\Pi^{3}\frac{r}{m}\,.$ (72d)

In these equations, $r$ is the circumferential radius, $m$ is the mass
function, $p$ is the fluid pressure, $\rho$ is the baryonic rest mass density,
$\varepsilon$ is the total energy density, and $\Pi=(\varepsilon-\rho)/\rho$
is the internal energy per unit baryonic mass. The parameters $\delta_{i}$,
$\pi_{i}$ and $\mu_{i}$ are the phenomenological post-TOV parameters. The TOV
equations are recovered when all the post-TOV parameters vanish. The
dimensionless terms $\mathscr{P}_{1}$, $\mathscr{M}_{1}$ and
$\mathscr{P}_{2}$, $\mathscr{M}_{2}$ describe the deviations from the TOV
equations at 1PN and 2PN orders respectively. In fact, the 1PN post-TOV
parameters $\delta_{i}$ can all be expressed in terms of the more familiar
parametrized post-Newtonian parameters which are strongly constrained by
observations.121212Stellar solutions at 1PN order in the parametrized post-
Newtonian formalism were obtained in Refs. Ciufolini:1981 ; Ciufolini:1983 ;
Shapiro:1976 ; WagonerMalone . Neutron stars were also modelled in the post-
Newtonian approximation in Ref. Shinkai:1998mg ; Gupta:2000zb , and more
recently in Ref. Andersson:2022cax . In practice, one can then neglect
$\mathscr{P}_{1}$ and $\mathscr{M}_{1}$ from Eqs. (71). In contrast, $\pi_{i}$
and $\mu_{i}$ are unconstrained and $\mathscr{P}_{2}$ and $\mathscr{M}_{2}$
should be viewed as describing the dominant departures from general
relativity.

The post-TOV formalism is also in a sense “complete” in that, when neglecting
the 1PN terms, the post-TOV equations can be derived by the covariant
conservation of a perfect fluid energy-momentum tensor. That is,

$\nabla_{\nu}T^{\mu\nu}=0\,,\quad T^{\mu\nu}=(\varepsilon_{\rm
eff}+p)u^{\mu}u^{\nu}+pg^{\mu\nu}\,,$ (73)

where $\varepsilon_{\rm eff}$ is an effective, gravity-modified energy
density,

$\varepsilon_{\rm eff}=\varepsilon+\rho\mathscr{M}_{2}\,,$ (74)

and the covariant derivative is compatible with the effective post-TOV metric

$g_{\mu\nu}=\textrm{diag}[\exp(2\nu),\,(1-2m/r)^{-1},\,r^{2},\,r^{2}\sin^{2}\theta]\,,$
(75)

where the lapse $\nu$ obeys,

$\frac{{\rm d}\nu}{{\rm
d}r}=\frac{1}{r^{2}}\left[(1-\mathscr{M}_{2})\frac{m+4\pi
r^{3}p}{1-2m/r}+m\mathscr{P}_{2}\right]\,.$ (76)

In the post-TOV formalism, the star’s mass and radius are determined in the
same way as in general relativity. One limitation of the post-TOV formalism is
that it cannot, by construction, capture nonperturbative effects, such as
spontaneous scalarization that we saw in Sec. 2.2. Some astrophysical
applications of the post-TOV formalism were worked out in Refs.
Glampedakis:2016pes ; Silva:2019wfn .

The post-TOV formalism is not the only parametrization of the TOV equations.
An ad-hoc alternative was proposed in Ref. Velten:2016bdk , building on Ref.
Schwab:2008ce . To understand this parametrization, it is convenient to first
rewrite the pressure gradient TOV equation (13) in the equivalent factorized
form

$\displaystyle\frac{{\rm d}p}{{\rm d}r}=-(\varepsilon+p)\,\frac{{\rm
d}\nu}{{\rm
d}r}=-\frac{Gm\varepsilon}{r^{2}}\,\left(1+\vphantom{\frac{2Gm\gamma}{r}}\frac{p}{\varepsilon}\right)\,\left(1-\frac{2Gm}{r}\right)^{-1}\,\left(1+\frac{4\pi
r^{3}p}{m}\right)\,.$ (77)

Each of the factors above has a specific origin. In the Newtonian limit, only
the first term survives, with the energy density $\varepsilon$ becoming the
mass density $\rho$, i.e.,

$\frac{{\rm d}p}{{\rm d}r}=-\rho\frac{Gm}{r^{2}}=-\rho\frac{{\rm d}\Phi}{{\rm
d}r}\,,$ (78)

where in the second equality, $\Phi$ is the Newtonian gravitational potential,
determined from Poisson’s equation. All other terms in Eq. (13) are
relativistic in origin. The first term, $(1+p/\varepsilon)$, acts as passive
(inertial) gravitational mass (normalized by the energy density). This can be
seen by comparing the first equality in Eq. (77) with the second equality in
Eq. (78). The second term, $(1+4\pi r^{3}p/m)$, acts as the active
gravitational mass density (normalized by the mass function). This can be seen
by identifying the factor $-(Gm/r^{2})$ in Eq. (78) with $-(Gm/r^{2})(1+4\pi
r^{3}p/m)$ in Eq. (77). At last, the third term, $1/(1-2Gm/r)$, is a
“geometric form factor” associated to spacetime curvature.

The idea of Refs. Schwab:2008ce ; Velten:2016bdk  is to introduce
dimensionless parameters to each of these factors. In practice, the
parametrized TOV equations are

$\displaystyle\frac{{\rm d}m}{{\rm d}r}$ $\displaystyle=4\pi
r^{2}\varepsilon\,\left(1+\frac{\sigma p}{\varepsilon}\right)\,,$ (79a)
$\displaystyle\frac{{\rm d}p}{{\rm d}r}$
$\displaystyle=-\frac{G(1+\alpha)m\varepsilon}{r^{2}}\,\left(1-\frac{2Gm\gamma}{r}\right)^{-1}\,\left(1+\vphantom{\frac{2Gm\gamma}{r}}\frac{\beta
p}{\varepsilon}\right)\,\left(1+\frac{4\pi\chi r^{3}p}{m}\right)\,,$ (79b)

where $\alpha$, $\beta$, $\gamma$, $\chi$ and $\sigma$ are free constants. The
TOV equations (13) are recovered for $\alpha=\sigma=0$ and
$\beta=\gamma=\chi=1$. From our previous discussion, we see that $\beta$,
$\gamma$ and $\chi$ affect the passive gravitational mass, geometric form
factor and active gravitational mass, respectively. The term $\alpha$
introduces an “effective gravitational constant”, $G_{\rm eff}=G(1+\alpha)$,
whereas $\sigma$ introduces a pressure contribution to the ${\rm d}m/{\rm d}r$
equation which does not exist in general relativity.

We note that neither parametrized TOV systems of equations (71) nor (79) has
been generalized to rotating or tidally deformed stars.

#### Extreme neutron star properties

Another way to circumvent the equation-of-state-gravity-theory degeneracy is
to look at _extremal_ properties of neutron stars. The idea is simple: with
few, physically sensible assumptions about the equation of state, one uses the
TOV equations to determine, for instance, what is the maximum mass $M_{\rm
max}$ a neutron star can have.

The first work in this direction was made by Rhoades and Ruffini in 1974
Rhoades:1974fn . Their motivation was to obtain a general upper bound on the
neutron star mass, which could then be used to distinguish them from black
holes in observations of x-ray sources in binary systems.

Specifically, Rhoades and Ruffini assumed that the star consists of two
layers:

  1. 1.

an outer layer, whose contribution to the star’s total mass $M$ is determined
according to a known equation of state;

  2. 2.

an inner core, whose contribution to $M$ must be extremized according to a
variational principle.

The two layers are matched at a density $\varepsilon_{0}$ near the nuclear
saturation density $\varepsilon_{\rm n}$. The variational calculation is then
performed under the assumptions:

  1. 1.

The star is described by the TOV equations (13).

  2. 2.

The energy density is positive, and related to the pressure by a barotropic
equation of state (14).

  3. 3.

That ${\rm d}p/{\rm d}\varepsilon>0$. Thus the speed of sound $v^{2}_{\rm
s}={\rm d}p/{\rm d}\varepsilon$ is positive, and matter is stable against
microscopic collapse.

  4. 4.

That $v_{\rm s}\leqslant 1$, and the equation of state is causal.

Items 2 to 4 are minimal physics assumptions on the unknown equation of state
at the inner core. Rhoades and Ruffini showed that the equation of state at
the core that maximizes $M$ is

$p=p_{0}+(\varepsilon-\varepsilon_{0}),\quad\textrm{for}\quad\varepsilon\geqslant\varepsilon_{0}$
(80)

where $p_{0}$ and $\varepsilon_{0}$ are the pressure and energy density at the
transition density of the two stellar layers. They found that,

$M_{\rm max}\simeq 3.2\,{\rm M}_{\odot}\,.$ (81)

This result depends on the value of $\varepsilon_{0}$ and the equation of
state assumed for the outer layer. In the late 90s, Kalogera and Baym
Kalogera:1996ci  revisited this upper bound obtaining $M_{\rm max}\simeq
2.9~{}{\rm M}_{\odot}$, assuming that $\varepsilon=2\,\varepsilon_{\rm n}$ and
that the outer layer is described by the equation of state of Wiringa et al.
Wiringa:1988tp .

Hartle and collaborators pushed this program further in the 1970’s,
culminating in Ref. Hartle:1978 . What is particularly interesting for us, is
that Ref. Hartle:1978 , Sec. 6.3, shows a handful of extensions of this
variational calculation to other theories of gravity, most of which have now
been severely constrained by observations.131313For the aficionado, the
examples are: Brans-Dicke theory, Rosen’s bimetric theory and Ni’s theory. The
reader can find summaries about these theories in Ref. Will:1993hxu , Chapter
5.

We encourage this program to be pursued further and applied to the “modern”
gravity theories of today. Observations of massive neutron stars consistent
with general relativity can be used to place minimally-equation-of-state-
dependent constraints on other theories. To do so, one has to derive the
prediction of $M_{\rm max}$ following Ref. Rhoades:1974fn  in a specific
theory.141414We note that Refs. Pani:2011xm ; Sotani:2017pfj  used the
equation of state (80) to compute the maximum neutron-star mass in modified
gravity theories. However, neither demonstrated through a calculation as in
Rhoades and Ruffini Rhoades:1974fn  that this _is_ the equation of state that
maximizes the neutron star mass in the context of theories they studied.
Next, one looks for the regions of the theory’s parameter space in which
neutron stars with $M_{\rm max}\geqslant M_{\rm obs}$ do not exist. This would
identify the regions of the theory’s parameter space that would be in tension
with observations. If we are ever to find a neutron star with a mass $M\gtrsim
3.2~{}{\rm M}_{\odot}$, it would be hard to explain this observation within
general relativity. It would be evidence that one (or more) of the assumptions
in Ref. Rhoades:1974fn  is wrong; including that the TOV equations of general
relativity hold.

##  4 Neutron-star pulse profiles

Much of our understanding of the extreme physics of neutron stars comes from
observations across the electromagnetic spectrum of astrophysical processes
taking place at the surface or at the vicinity of these stars. See, for
instance, Refs. Nattila:2022evn ; Ozel:2012wu ; Psaltis:2008bb  for reviews.
In this section, we focus on one promising way of testing theories of gravity
in the x-ray spectrum by observations of neutron-star pulse profiles, also
known as lightcurves. These observations in the soft x-ray band made by the
Neutron Star Interior Composition Explorer (NICER) telescope are what enabled
the moment of inertia inference quoted in Sec. 3.2, and used in the
parametrized I-Love test described there.

The general picture is the following. By some physical process a region on the
surface of a rotating neutron stars gets heated up (relative to the rest of
the star), forming what is called a hotspot. In general, the radiation
produced by the hotspot will appear to a distant observer as a periodic flux,
modulated by the star’s rotation frequency. Photons emitted by the hotspot are
subject to a number of relativistic effects as they travel through the
spacetime. These effects, for example, gravitational time delay and
gravitational deflection leave an imprint of the properties of the neutron-
star spacetime onto the observed x-ray flux.

We saw in Sec. 2.1 that in general relativity the exterior spacetime is
described by the Schwarzschild metric. It is then conceivable that with
sufficient sensitivity (to detect as many photon as possible) and time-
resolution (to resolve in time as well as possible the incoming flux) one may
be able to infer the mass and radius of the radiating star. In a gross
oversimplification, this is what NICER has achieved with observations of the
isolated rotation-powered millisecond pulsar PSR J0030+0451 Miller:2019cac ;
Riley:2019yda  and the massive pulsar PSR J0740+6620, the latter combining
data from XMN-Newton Miller:2021qha ; Riley:2021pdl . The results of the
former are particularly spectacular since PSR J0030+0451 is an isolated
pulsar. Two independent analysis of this system suggest that it has mass
$M=1.44^{+0.15}_{-0.14}~{}{\rm M}_{\odot}$ and (equatorial) radius
$R=13.02^{+1.24}_{-1.06}$ km as inferred by Miller:2021qha  and
$M=1.34^{+0.15}_{-0.16}~{}{\rm M}_{\odot}$ and $R=12.71^{+1.14}_{-1.19}$ km as
inferred by Riley:2019yda

Figure 7: The bolometric flux produced by a single hotpot during one
revolution of neutron star for two values of the scalar-charge-to-mass ration
$Q$ according to the J+D model. In both cases, the star is assumed to have the
same mass, $1.4~{}{\rm M}_{\odot}$, compactness, $C=0.25$ and rotational
frequency of $600$ Hz. All other parameters of the model, such as the hotspot
location and its temperature, are also fixed. The only difference is the value
of $Q$. Figure taken from Ref. Silva:2017uqg .

Given that pulse profiles are sensitive to the neutron-star spacetime, it is
natural to ask: can they be used to test general relativity and other
alternatives? To answer this question one must first develop a pulse profile
model in any given theory of gravity. In general relativity, the earliest work
is this direction dates back to 1983, when Pechenick et al. Pechenick:1983apj
computed the flux of circular antipodal polar caps for a slowly-rotating
neutron star including the effect of gravitational light bending, but
neglecting frame-dragging effects.151515We mention in passing that the first
visualization of a neutron star surface was done by H.-P. Nollert et al.
Nollert:1989 . Light deflection causes the visible fraction of the star’s
surface to be larger than it would be if the star was embedded in flat
spacetime. This makes the star to appear larger than it would be otherwise.
See also, e.g., Ref. Nattila:2022evn , Fig. 1, for a similar visualization
including the effects of rapid rotation obtained using the ray-tracing code of
Ref. Pihajoki:2018ihj . Later works introduced special-relativistic effects of
Doppler boost and aberration caused by the rapid motion of the hotspot, and
the general-relativistic contribution due to gravitational-time delay in the
photon’s arrival time at the observer. This culminated in what is known as the
“Schwarzschild plus Doppler” (S+D) pulse profile model of Miller and Lamb
Miller:1997if  (see also Refs. Poutanen:2003yd ; Poutanen:2006hw ). At high
rotation frequencies ($f\gtrsim 300$ Hz), the dominant effect on the pulse
profile comes not from frame dragging, but from the spin-induced quadrupolar
deformation of the star Braje:2000qb ; Cadeau:2004gm ; Cadeau:2006dc ;
Morsink:2007tv ; Psaltis:2013zja ; Nattila:2017hdb . One can then assume that
photons are emitted from an oblate surface and then propagate in an ambient
Schwarzschild background.161616The OS model also uses an analytical
description of how the star’s surface deforms as a function of its mass,
equatorial radius and rotational frequency. See Refs. AlGendy:2014eua ;
Morsink:2007tv ; Silva:2020acr  for such formulas. This “Oblate-Schwarzschild”
(OS) approximation Cadeau:2006dc ; Morsink:2007tv  is the model used by NICER
Bogdanov:2019qjb .171717In principle, one could do ray-tracing on the
numerical spacetime of a rotating neutron star. In practice, due to the large
parameter space of the problem, it is computationally prohibitive to use this
approach for parameter estimation using Bayesian inference. The alternative is
to use the approximate models described here.

Figure 8: Projected constraints on $\alpha_{0}$ and $\beta_{0}$ parameters of
the scalar-tensor gravity model of Refs. Damour:1993hw ; Damour:1996ke  using
pulse profile observations. In this plane, general relativity is located at
$|\alpha_{0}|=0$ and $\beta_{0}=0$. The shaded area is the region where
spontaneous scalarization happens, and depend onle weakly on the equation of
state. The region to the left of the nearly vertical solid and dashed lines
are excluded to 68% for a star of mass $M=1.93~{}{\rm M}_{\odot}$, described
either by the SLy4 Douchin:2001sv  or APR Akmal:1998cf  equation of state,
spinning at 200 or 400 Hz, respectively. For comparison, the figure also shows
68% credible regions constrained by the Cassini mission Bertotti:2003rm
(dotted horizon line), the binary pulsars J1738-0333, J0348-0432
Anderson:2019eay , and the hierarchical triple system PSR J0447-1715
Archibald:2018oxs . This figure is taken from Ref. Silva:2019leq . Since then,
stronger binary-pulsar constraints were obtained in Refs. Kramer:2021jcw ;
Zhao:2022vig . These constraint strongly constrain the neutron-star
scalarization, at least in the context of the original Damour-Esposito-Farèse
model Damour:1993hw ; Damour:1996ke .

Da̧browski and Osarczuk Dabrowski:1995  and Sotani and Miyamoto Sotani:2017bho
; Sotani:2017rrt  explored how lifting the assumption of the spacetime being
Schwarzschild modifies the pulse profile. The first model including the
minimal level of astrophysical realism necessary for pulse-profile analysis in
scalar-tensor gravity was developed by Silva and Yunes Silva:2018yxz . They
took into account all the physical effects included in the S+D model, assuming
that the exterior spacetime is described by the Just metric; see Sec. 2.2.
This “Just+Doppler” (J+D) model has the S+D model as the particular limit in
which the scalar-charge-to-mass ratio $Q$ (40) becomes zero. Figure 7 shows an
example pulse profile computed within the J+D approximation.

Reference Silva:2019leq  then used the J+D model to make a first forecast of
what constrains one could expect to place on spontaneous scalarization through
NICER-like observations consistent with general relativity. The main
conclusion of this work is that, in principle, these constraints can be
competitive with those obtained from binary-pulsar observations, as summarized
in Fig. 8.

Despite this promising result it is fair to say that the subject has not
matured enough to give a definite answer on how useful pulse profile
observations are to test theories of gravity. For example, the analysis of
Ref. Silva:2019leq  assumed for simplicity a single pointlike hotspot and
adopted a restricted Bayesian calculation. It would be interesting to perform
a more complete analysis, similar to that of Refs. Lo:2013ava ; Miller:2014mca
, in the J+D model to explore relative contributions (and identify potential
degeneracies) between hotspot properties and modifications to general
relativity. Moreover, none of the beyond-general-relativity pulse profile
studies mentioned previously, as well as posterior works in massive scalar-
tensor gravity Hu:2021tyw ; Xu:2020vbs , considered stellar oblateness.
Eventually, of course, it would be important to use these pulse-profile models
in a Bayesian inference analysis of the data collected by NICER
Bogdanov:2019ixe . It is quite clear that this subject remains a fresh, and
likely fertile, ground for work.

##  5 Compact binaries and gravitational waves

The first detection of gravitational waves produced by the coalescence of two
black holes on September 14, 2015, by the LIGO interferometers in Hanford and
Livingstone, USA, marked the beginning of the field of gravitational-wave
astronomy LIGOScientific:2016aoc . Since then, the detectors of the LIGO-
Virgo-Kagra collaboration have detected more than 90 gravitational-wave events
produced by binaries with a black hole or a neutron star as components. The
events have been used to perform a multitude of tests of gravity, as
exemplified in Refs. LIGOScientific:2016lio ; LIGOScientific:2018dkp ;
LIGOScientific:2019fpa ; LIGOScientific:2020tif ; LIGOScientific:2021sio ;
Yunes:2016jcc . So far, all tests have shown consistency between the
predictions of general relativity with observations. In this section, we
review the implications to tests of gravity of three particular gravitational-
wave events, where at least one binary component was inferred to be a neutron
star.

###  5.1 GW170817: the first binary neutron star event

In August 17, 2017, the LIGO and Virgo detectors observed the gravitational-
wave signal produced by the inspiral of two neutron stars. The signal stayed
approximately 100 s in the sensitivity band of the detectors and had a signal-
to-noise ratio of 32.4. The data from the three detectors allowed a precise
position in sky localization which together with an almost simultaneous (only
$\approx 1.7$ s time delay) detection of $\gamma$-rays by the Fermi Gamma-ray
Burst Monitor enabled a massive effort to search for an electromagnetic
counterpart to this event, confirmed by several observatories. These
multimessenger observations gave support to the interpretation of GW170817
being a neutron-star binary.

An important outcome of this event is that the almost simultaneous
observations of gravitational-wave and electromagnetic signals that traveled a
distance of approximately 130 million light years, placed the strongest
constrain to date on the propagation speed $c_{g}$ of gravitational waves

$-3\cdot 10^{-15}\leqslant c_{g}/c-1\leqslant 7\cdot 10^{-16}.$ (82)

This result had a profound impact to extension to general relativity that are
studied as possible explanations for dark energy or as alternatives to the
standard cosmological model, $\Lambda$CDM.181818This is often not the case for
the theories discussed in Sec. 3. These theories generically predict that
gravitational waves propagate with speed $c_{g}\neq c$ over a cosmological
background. Hence, the strong bound (82) severely constrains such theories.
See, e.g., Refs. Baker:2017hug ; Creminelli:2017sry ; Ezquiaga:2017ekz ;
Sakstein:2017xjx .

Let us mention two additional consequences of this event for tests of gravity.
First, the agreement between the gravitational-wave signal with the
predictions of general relativity placed a bound on the emission of scalar
dipole radiation from the system. However, this bound is weaker in comparison
with those obtained with binary pulsar observations LIGOScientific:2018dkp ;
Kramer:2021jcw . Second, GW170817 enabled the previously mentioned estimate of
the tidal deformability of a $1.4~{}{\rm M}_{\odot}$ neutron star,
$\bar{\lambda}_{1.4}=190^{+390}_{-120}$ LIGOScientific:2018cki  that was used
to perform the parametrized I-Love test of gravity discussed in Sec. 3.2. See
Ref. Chatziioannou:2020pqz  for a review on the role of tidal deformabality in
gravitational-wave observations and its implications to constraints on the
neutron-star equation of state.

###  5.2 GW200105 and GW200115: the neutron-star-black-hole binaries

In the interval of a few days in January, 2020, the LIGO and Virgo detectors
observed the first neutron-star-black-hole binaries LIGOScientific:2021qlt .
These two events, GW200105 and GW200115, are interesting for tests of gravity
for the following reasons. The contribution to the gravitational-wave phase
due to scalar dipole radiation is proportional to the square of the difference
of the scalar charges of the two objects. As we saw back in Sec. 2.3, in
shift-symmetric Einstein-scalar-Gauss-Bonnet gravity neutron stars have no
monopolar scalar charge, whereas black holes do. This means that the emission
of scalar dipole radiation is controlled only by the charge of the black hole,
that in turn depends on $\ell_{\rm GB}$. In addition, the relatively low mass
of the black hole components in these two events ($8.9^{+1.2}_{-1.5}~{}{\rm
M}_{\odot}$ and $5.7^{+1.8}_{-0.3}~{}{\rm M}_{\odot}$ for GW200105 and
GW200115, respectively) magnifies the effects of higher-curvature corrections,
enhancing the magnitude of the scalar charge of the black holes. Consequently,
the consistency of the data with general relativity translates into a strong
bound on this theory. Lyu et al. Lyu:2022gdr  exploited these facts to place
the bound $\ell_{\rm GB}\lesssim 5$ km using GW200115.191919This corresponds
to $\alpha_{\rm GB}^{1/2}\lesssim 1.33$ km in the notation used in Ref.
Lyu:2022gdr . The two quantities are related by $\ell_{\rm
GB}=(64\pi)^{1/4}~{}\alpha_{\rm GB}^{1/2}$. This bound is compatible with
other astrophysical constraints on this theory Yagi:2012gp ; Saffer:2021gak
and improves the bound $\ell_{\rm GB}\lesssim 6$ km obtained previously in
Perkins et al Perkins:2021mhb , who used binary black hole events that had
been detected by then. See Refs. Julie:2019sab ; Yagi:2011xp ;
Shiralilou:2020gah ; Shiralilou:2021mfl  for works in post-Newtonian theory
that enable these constraints to be obtained.

##  6 Summary

We presented an overview of how modifications to general relativity impact the
structure of neutron stars and discussed some of its astrophysical
implications. We discussed the difficulties in testing these modified theories
of gravity using observations of neutron stars, and how these difficulties
originate from the large uncertainties in our knowledge about the neutron star
equation of state at supranucler densities. We saw how quasiuniversal
relations can be helpful to mitigate this problem, and reviewed some examples
of “theory-agnostic” approaches to describe the structure of neutron stars.
With the hope of motivating further work, we also reviewed the current status
of pulse profile modeling in modified gravity theories, and identified some
directions for future research. At last, we presented a small selection of
constraints on modified theories of gravity that were obtained using
gravitational wave observations of compact binaries that contain at least one
neutron star.

It is not an exaggeration to say that the observations made by NICER and the
detection of gravitational-wave signals by the LIGO-Virgo-Kagra collaboration
have opened a new era in the study of neutron stars. Will these objects
provide us the first whispers of new physics beyond general relativity? This
is for time to answer. However, one thing is certain: neutron stars will
continue to play a fundamental role in our understanding of gravity at its
extreme.

###### Acknowledgements.

We thank Cosimo Bambi and Alejandro Cárdenas-Avendaño for the invitation to
write this chapter, and to participate in the “Recent progress on gravity
tests” workshop in 2022. We are also indebted to our collaborators with whom
we have studied various aspects of tests of gravity through the neutron-star
perspective over the years. We thank Helvi Witek for the encouragement in the
final stages of writing. We acknowledge funding from the Deutsche
Forschungsgemeinschaft (DFG) – Project No. 386119226

## References

  * (1) B. P. Abbott _et al._ , Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116 (2016) 061102.  doi:10.1103/PhysRevLett.116.061102.  arXiv:1602.03837.
  * (2) B. P. Abbott _et al._ , Tests of general relativity with GW150914, Phys. Rev. Lett. 116 (2016) 221101.  doi:10.1103/PhysRevLett.116.221101.  [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841.
  * (3) B. P. Abbott _et al._ , GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119 (2017) 161101.  doi:10.1103/PhysRevLett.119.161101.  arXiv:1710.05832.
  * (4) B. P. Abbott _et al._ , GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121 (2018) 161101.  doi:10.1103/PhysRevLett.121.161101.  arXiv:1805.11581.
  * (5) B. P. Abbott _et al._ , Tests of General Relativity with GW170817, Phys. Rev. Lett. 123 (2019) 011102.  doi:10.1103/PhysRevLett.123.011102.  arXiv:1811.00364.
  * (6) B. P. Abbott _et al._ , Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1, Phys. Rev. D 100 (2019) 104036.  doi:10.1103/PhysRevD.100.104036.  arXiv:1903.04467.
  * (7) R. Abbott _et al._ , Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences, Astrophys. J. Lett. 915 (2021) L5.  doi:10.3847/2041-8213/ac082e.  arXiv:2106.15163.
  * (8) R. Abbott _et al._ , Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103 (2021) 122002.  doi:10.1103/PhysRevD.103.122002.  arXiv:2010.14529.
  * (9) R. Abbott _et al._ Tests of General Relativity with GWTC-3 (2021).  arXiv:2112.06861.
  * (10) A. Akmal, V. R. Pandharipande and D. G. Ravenhall, The Equation of state of nucleon matter and neutron star structure, Phys. Rev. C 58 (1998) 1804–1828.  doi:10.1103/PhysRevC.58.1804.  arXiv:nucl-th/9804027.
  * (11) S. Alexander and N. Yunes, Chern-Simons Modified General Relativity, Phys. Rept. 480 (2009) 1–55.  doi:10.1016/j.physrep.2009.07.002.  arXiv:0907.2562.
  * (12) M. AlGendy and S. M. Morsink, Universality of the Acceleration Due to Gravity on the Surface of a Rapidly Rotating Neutron Star, Astrophys. J. 791 (2014) 78.  doi:10.1088/0004-637X/791/2/78.  arXiv:1404.0609.
  * (13) Y. Ali-Haimoud and Y. Chen, Slowly-rotating stars and black holes in dynamical Chern-Simons gravity, Phys. Rev. D 84 (2011) 124033.  doi:10.1103/PhysRevD.84.124033.  arXiv:1110.5329.
  * (14) J. Alsing, E. Berti, C. M. Will _et al._ , Gravitational radiation from compact binary systems in the massive Brans-Dicke theory of gravity, Phys. Rev. D 85 (2012) 064041.  doi:10.1103/PhysRevD.85.064041.  arXiv:1112.4903.
  * (15) J. Alsing, H. O. Silva and E. Berti, Evidence for a maximum mass cut-off in the neutron star mass distribution and constraints on the equation of state, Mon. Not. Roy. Astron. Soc. 478 (2018) 1377–1391.  doi:10.1093/mnras/sty1065.  arXiv:1709.07889.
  * (16) D. Anderson, P. Freire and N. Yunes, Binary pulsar constraints on massless scalar–tensor theories using Bayesian statistics, Class. Quant. Grav. 36 (2019) 225009.  doi:10.1088/1361-6382/ab3a1c.  arXiv:1901.00938.
  * (17) N. Andersson, Whispers from the edge of physics, J. Astrophys. Astron. 38 (2017) 58.  doi:10.1007/s12036-017-9463-8.  arXiv:1709.07215.
  * (18) N. Andersson, A gravitational-wave perspective on neutron-star seismology, Universe 7 (2021) 97.  doi:10.3390/universe7040097.  arXiv:2103.10223.
  * (19) N. Andersson, F. Gittins, S. Yin _et al._ , Building post-Newtonian neutron stars, Class. Quant. Grav. 40 (2023) 025016.  doi:10.1088/1361-6382/acace5.  arXiv:2209.05871.
  * (20) N. Andersson and K. D. Kokkotas, Towards gravitational wave asteroseismology, Mon. Not. Roy. Astron. Soc. 299 (1998) 1059–1068.  doi:10.1046/j.1365-8711.1998.01840.x.  arXiv:gr-qc/9711088.
  * (21) J. Antoniadis _et al._ , A Massive Pulsar in a Compact Relativistic Binary, Science 340 (2013) 6131.  doi:10.1126/science.1233232.  arXiv:1304.6875.
  * (22) A. M. Archibald, N. V. Gusinskaia, J. W. T. Hessels _et al._ , Universality of free fall from the orbital motion of a pulsar in a stellar triple system, Nature 559 (2018) 73–76.  doi:10.1038/s41586-018-0265-1.  arXiv:1807.02059.
  * (23) T. Baker, E. Bellini, P. G. Ferreira _et al._ , Strong constraints on cosmological gravity from GW170817 and GRB 170817A, Phys. Rev. Lett. 119 (2017) 251301.  doi:10.1103/PhysRevLett.119.251301.  arXiv:1710.06394.
  * (24) G. Baym, T. Hatsuda, T. Kojo _et al._ , From hadrons to quarks in neutron stars: a review, Rept. Prog. Phys. 81 (2018) 056902.  doi:10.1088/1361-6633/aaae14.  arXiv:1707.04966.
  * (25) P. Bedaque and A. W. Steiner, Sound velocity bound and neutron stars, Phys. Rev. Lett. 114 (2015) 031103.  doi:10.1103/PhysRevLett.114.031103.  arXiv:1408.5116.
  * (26) O. Benhar, V. Ferrari, L. Gualtieri _et al._ , Perturbative approach to the structure of rapidly rotating neutron stars, Phys. Rev. D 72 (2005) 044028.  doi:10.1103/PhysRevD.72.044028.  arXiv:gr-qc/0504068.
  * (27) P. G. Bergmann, Comments on the scalar tensor theory, Int. J. Theor. Phys. 1 (1968) 25–36.  doi:10.1007/BF00668828.
  * (28) E. Berti _et al._ , Testing General Relativity with Present and Future Astrophysical Observations, Class. Quant. Grav. 32 (2015) 243001.  doi:10.1088/0264-9381/32/24/243001.  arXiv:1501.07274.
  * (29) B. Bertotti, L. Iess and P. Tortora, A test of general relativity using radio links with the Cassini spacecraft, Nature 425 (2003) 374–376.  doi:10.1038/nature01997.
  * (30) S. Bogdanov _et al._ , Constraining the Neutron Star Mass–Radius Relation and Dense Matter Equation of State with $NICER$. I. The Millisecond Pulsar X-Ray Data Set, Astrophys. J. Lett. 887 (2019) L25.  doi:10.3847/2041-8213/ab53eb.  arXiv:1912.05706.
  * (31) S. Bogdanov _et al._ , Constraining the Neutron Star Mass–Radius Relation and Dense Matter Equation of State with $NICER$. II. Emission from Hot Spots on a Rapidly Rotating Neutron Star, Astrophys. J. Lett. 887 (2019) L26.  doi:10.3847/2041-8213/ab5968.  arXiv:1912.05707.
  * (32) T. M. Braje, R. W. Romani and K. P. Rauch, Light curves of rapidly rotating neutron stars, Astrophys. J. 531 (2000) 447.  doi:10.1086/308448.  arXiv:astro-ph/0004411.
  * (33) C. Brans and R. H. Dicke, Mach’s principle and a relativistic theory of gravitation, Phys. Rev. 124 (1961) 925–935.  doi:10.1103/PhysRev.124.925.
  * (34) C. Breu and L. Rezzolla, Maximum mass, moment of inertia and compactness of relativistic stars, Mon. Not. Roy. Astron. Soc. 459 (2016) 646–656.  doi:10.1093/mnras/stw575.  arXiv:1601.06083.
  * (35) R. Brito, A. Buonanno and V. Raymond, Black-hole Spectroscopy by Making Full Use of Gravitational-Wave Modeling, Phys. Rev. D 98 (2018) 084038.  doi:10.1103/PhysRevD.98.084038.  arXiv:1805.00293.
  * (36) H. A. Buchdahl, Reciprocal Static Metrics and Scalar Fields in the General Theory of Relativity, Phys. Rev. 115 (1959) 1325–1328.  doi:10.1103/PhysRev.115.1325.
  * (37) C. Cadeau, D. A. Leahy and S. M. Morsink, Pulse shapes from rapidly-rotating neutron stars: Equatorial photon orbits, Astrophys. J. 618 (2005) 451–462.  doi:10.1086/425857.  arXiv:astro-ph/0409261.
  * (38) C. Cadeau, S. M. Morsink, D. Leahy _et al._ , Light Curves for Rapidly-Rotating Neutron Stars, Astrophys. J. 654 (2007) 458–469.  doi:10.1086/509103.  arXiv:astro-ph/0609325.
  * (39) P. A. Cano and A. Ruipérez, String gravity in D=4, Phys. Rev. D 105 (2022) 044022.  doi:10.1103/PhysRevD.105.044022.  arXiv:2111.04750.
  * (40) T. K. Chan, A. P. O. Chan and P. T. Leung, I-Love relations for incompressible stars and realistic stars, Phys. Rev. D 91 (2015) 044017.  doi:10.1103/PhysRevD.91.044017.  arXiv:1411.7141.
  * (41) S. Chandrasekhar.  An introduction to the study of stellar structure (1939).
  * (42) S. Chandrasekhar and J. C. Miller, On slowly rotating homogeneous masses in general relativity, Mon. Not. Roy. Astron. Soc. 167 (1974) 63–80.  doi:10.1093/mnras/167.1.63.
  * (43) K. Chatziioannou, Neutron star tidal deformability and equation of state constraints, Gen. Rel. Grav. 52 (2020) 109.  doi:10.1007/s10714-020-02754-3.  arXiv:2006.03168.
  * (44) K. Chatziioannou, N. Yunes and N. Cornish, Model-Independent Test of General Relativity: An Extended post-Einsteinian Framework with Complete Polarization Content, Phys. Rev. D 86 (2012) 022004.  doi:10.1103/PhysRevD.86.022004.  [Erratum: Phys.Rev.D 95, 129901 (2017)], arXiv:1204.2585.
  * (45) P. Chen, T. Suyama and J. Yokoyama, Spontaneous scalarization: asymmetron as dark matter, Phys. Rev. D 92 (2015) 124016.  doi:10.1103/PhysRevD.92.124016.  arXiv:1508.01384.
  * (46) I. Ciufolini and R. Ruffini, On the value of the masses of neutron stars in the parameterized post-Newtonian formalism, Astron. Astrophys. 97 (1981) L12–L14.
  * (47) I. Ciufolini and R. Ruffini, Equilibrium configurations of neutron stars and the parametrized post-Newtonian metric theories of gravitation, Astrophys. J. 275 (1983) 867–877.  doi:10.1086/161580.
  * (48) R. Coquereaux and G. Esposito-Farèse, The Theory of Kaluza-Klein-Jordan-Thiry revisited, Ann. Inst. H. Poincare Phys. Theor. 52 (1990) 113–150.
  * (49) G. Creci, T. Hinderer and J. Steinhoff, Tidal properties of neutron stars in scalar-tensor theories of gravity, Phys. Rev. D 108 (2023) 124073.  doi:10.1103/PhysRevD.108.124073.  arXiv:2308.11323.
  * (50) P. Creminelli and F. Vernizzi, Dark Energy after GW170817 and GRB170817A, Phys. Rev. Lett. 119 (2017) 251302.  doi:10.1103/PhysRevLett.119.251302.  arXiv:1710.05877.
  * (51) H. T. Cromartie _et al._ , Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar, Nature Astron. 4 (2019) 72–76.  doi:10.1038/s41550-019-0880-2.  arXiv:1904.06759.
  * (52) T. Damour and N. Deruelle, General relativistic celestial mechanics of binary systems. II. The post-Newtonian timing formula., Annales de L’Institut Henri Poincare Section (A) Physique Theorique 44 (1986) 263–292.
  * (53) T. Damour and G. Esposito-Farèse, Tensor multiscalar theories of gravitation, Class. Quant. Grav. 9 (1992) 2093–2176.  doi:10.1088/0264-9381/9/9/015.
  * (54) T. Damour and G. Esposito-Farèse, Nonperturbative strong field effects in tensor - scalar theories of gravitation, Phys. Rev. Lett. 70 (1993) 2220–2223.  doi:10.1103/PhysRevLett.70.2220.
  * (55) T. Damour and G. Esposito-Farèse, Tensor-scalar gravity and binary pulsar experiments, Phys. Rev. D 54 (1996) 1474–1491.  doi:10.1103/PhysRevD.54.1474.  arXiv:gr-qc/9602056.
  * (56) T. Damour and J. H. Taylor, Strong field tests of relativistic gravity and binary pulsars, Phys. Rev. D 45 (1992) 1840–1868.  doi:10.1103/PhysRevD.45.1840.
  * (57) M. P. Da̧browski and J. Osarczuk, Light curves of relativistic charged neutron star, Astrophysics and Space Science 229 (1995) 139–155.  doi:10.1007/BF00658572.
  * (58) T. Delsate, D. Hilditch and H. Witek, Initial value formulation of dynamical Chern-Simons gravity, Phys. Rev. D 91 (2015) 024027.  doi:10.1103/PhysRevD.91.024027.  arXiv:1407.6727.
  * (59) E. S. Demirboğa, Y. E. Şahin and F. M. Ramazanoğlu, Subtleties in constraining gravity theories with mass-radius data, Phys. Rev. D 108 (2023) 024028.  doi:10.1103/PhysRevD.108.024028.  arXiv:2303.01910.
  * (60) D. D. Doneva and G. Pappas, Universal Relations and Alternative Gravity Theories, Astrophys. Space Sci. Libr. 457 (2018) 737–806.  arXiv:1709.08046.
  * (61) D. D. Doneva, F. M. Ramazanoğlu, H. O. Silva _et al._ Scalarization (2022).  arXiv:2211.01766.
  * (62) D. D. Doneva and S. S. Yazadjiev, Neutron star solutions with curvature induced scalarization in the extended Gauss-Bonnet scalar-tensor theories, JCAP 04 (2018) 011.  doi:10.1088/1475-7516/2018/04/011.  arXiv:1712.03715.
  * (63) F. Douchin and P. Haensel, A unified equation of state of dense matter and neutron star structure, Astron. Astrophys. 380 (2001) 151.  doi:10.1051/0004-6361:20011402.  arXiv:astro-ph/0111092.
  * (64) J. M. Ezquiaga and M. Zumalacárregui, Dark Energy After GW170817: Dead Ends and the Road Ahead, Phys. Rev. Lett. 119 (2017) 251304.  doi:10.1103/PhysRevLett.119.251304.  arXiv:1710.05901.
  * (65) E. Fonseca _et al._ , Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620, Astrophys. J. Lett. 915 (2021) L12.  doi:10.3847/2041-8213/ac03b8.  arXiv:2104.00880.
  * (66) N. Franchini and S. H. Völkel.  Testing General Relativity with Black Hole Quasi-Normal Modes (2023).  arXiv:2305.01696.
  * (67) J. L. Friedman and N. Stergioulas.  Rotating Relativistic Stars.  Cambridge Monographs on Mathematical Physics. Cambridge University Press (2013).  ISBN 978-1-107-30217-4, 978-0-521-87254-6.  doi:10.1017/CBO9780511977596.
  * (68) A. Ghosh, R. Brito and A. Buonanno, Constraints on quasinormal-mode frequencies with LIGO-Virgo binary–black-hole observations, Phys. Rev. D 103 (2021) 124041.  doi:10.1103/PhysRevD.103.124041.  arXiv:2104.01906.
  * (69) K. Glampedakis, G. Pappas, H. O. Silva _et al._ , Post-Tolman-Oppenheimer-Volkoff formalism for relativistic stars, Phys. Rev. D 92 (2015) 024056.  doi:10.1103/PhysRevD.92.024056.  arXiv:1504.02455.
  * (70) K. Glampedakis, G. Pappas, H. O. Silva _et al._ , Astrophysical applications of the post-Tolman-Oppenheimer-Volkoff formalism, Phys. Rev. D 94 (2016) 044030.  doi:10.1103/PhysRevD.94.044030.  arXiv:1606.05106.
  * (71) E. Gourgoulhon.  3+1 formalism and bases of numerical relativity (2007).  arXiv:gr-qc/0703035.
  * (72) E. Gourgoulhon.  An Introduction to the theory of rotating relativistic stars (2010).  arXiv:1003.5015.
  * (73) A. Gupta, A. Gopakumar, B. R. Iyer _et al._ , Padé approximants for truncated postNewtonian neutron star models, Phys. Rev. D 62 (2000) 044038\.  doi:10.1103/PhysRevD.62.044038.  arXiv:gr-qc/0002094.
  * (74) T. Harada, Stability analysis of spherically symmetric star in scalar-tensor theories of gravity, Prog. Theor. Phys. 98 (1997) 359–379.  doi:10.1143/PTP.98.359.  arXiv:gr-qc/9706014.
  * (75) T. Harada, Neutron stars in scalar tensor theories of gravity and catastrophe theory, Phys. Rev. D 57 (1998) 4802–4811.  doi:10.1103/PhysRevD.57.4802.  arXiv:gr-qc/9801049.
  * (76) J. B. Hartle, Slowly rotating relativistic stars. 1. Equations of structure, Astrophys. J. 150 (1967) 1005–1029.  doi:10.1086/149400.
  * (77) J. B. Hartle, Bounds on the mass and moment of inertia of nonrotating neutron stars., Phys. Rept. 46 (1978) 201–247.  doi:10.1016/0370-1573(78)90140-0.
  * (78) J. B. Hartle and K. S. Thorne, Slowly Rotating Relativistic Stars. II. Models for Neutron Stars and Supermassive Stars, Astrophys. J. 153 (1968) 807\.  doi:10.1086/149707.
  * (79) C. A. R. Herdeiro and E. Radu, Asymptotically flat black holes with scalar hair: a review, Int. J. Mod. Phys. D 24 (2015) 1542014.  doi:10.1142/S0218271815420146.  arXiv:1504.08209.
  * (80) A. Hewish, S. J. Bell, J. D. H. Pilkington _et al._ , Observation of a rapidly pulsating radio source, Nature 217 (1968) 709–713.  doi:10.1038/217709a0.
  * (81) M. W. Horbatsch and C. P. Burgess, Semi-Analytic Stellar Structure in Scalar-Tensor Gravity, JCAP 08 (2011) 027.  doi:10.1088/1475-7516/2011/08/027.  arXiv:1006.4411.
  * (82) G. W. Horndeski, Second-order scalar-tensor field equations in a four-dimensional space, Int. J. Theor. Phys. 10 (1974) 363–384.  doi:10.1007/BF01807638.
  * (83) Z. Hu, Y. Gao, R. Xu _et al._ , Scalarized neutron stars in massive scalar-tensor gravity: X-ray pulsars and tidal deformability, Phys. Rev. D 104 (2021) 104014.  doi:10.1103/PhysRevD.104.104014.  arXiv:2109.13453.
  * (84) Z. Hu, X. Miao and L. Shao.  Tests of Classical Gravity with Radio Pulsars (2023).  arXiv:2303.17185.
  * (85) R. Jackiw and S. Y. Pi, Chern-Simons modification of general relativity, Phys. Rev. D 68 (2003) 104012.  doi:10.1103/PhysRevD.68.104012.  arXiv:gr-qc/0308071.
  * (86) N. Jiang and K. Yagi, Improved Analytic Modeling of Neutron Star Interiors, Phys. Rev. D 99 (2019) 124029.  doi:10.1103/PhysRevD.99.124029.  arXiv:1904.05954.
  * (87) N. Jiang and K. Yagi, Analytic I-Love-C relations for realistic neutron stars, Phys. Rev. D 101 (2020) 124006.  doi:10.1103/PhysRevD.101.124006.  arXiv:2003.10498.
  * (88) F.-L. Julié and E. Berti, Post-Newtonian dynamics and black hole thermodynamics in Einstein-scalar-Gauss-Bonnet gravity, Phys. Rev. D 100 (2019) 104061.  doi:10.1103/PhysRevD.100.104061.  arXiv:1909.05258.
  * (89) K. Just, Notizen: The Motion of Mercury according to the Theory of Thiry and Lichnerowicz, Zeitschrift Naturforschung Teil A 14 (1959) 751–751.  doi:10.1515/zna-1959-0810.
  * (90) V. Kalogera and G. Baym, The maximum mass of a neutron star, Astrophys. J. Lett. 470 (1996) L61–L64.  doi:10.1086/310296.  arXiv:astro-ph/9608059.
  * (91) T. Kobayashi, Horndeski theory and beyond: a review, Rept. Prog. Phys. 82 (2019) 086901.  doi:10.1088/1361-6633/ab2429.  arXiv:1901.07183.
  * (92) K. D. Kokkotas and J. Ruoff, Radial oscillations of relativistic stars, Astron. Astrophys. 366 (2001) 565.  doi:10.1051/0004-6361:20000216.  arXiv:gr-qc/0011093.
  * (93) M. Kramer _et al._ , Strong-Field Gravity Tests with the Double Pulsar, Phys. Rev. X 11 (2021) 041050.  doi:10.1103/PhysRevX.11.041050.  arXiv:2112.06795.
  * (94) J. M. Lattimer, Neutron Stars and the Nuclear Matter Equation of State, Ann. Rev. Nucl. Part. Sci. 71 (2021) 433–464.  doi:10.1146/annurev-nucl-102419-124827.
  * (95) J. M. Lattimer and M. Prakash, Neutron star structure and the equation of state, Astrophys. J. 550 (2001) 426.  doi:10.1086/319702.  arXiv:astro-ph/0002232.
  * (96) J. M. Lattimer and A. Yahil, Analysis of the neutrino events from supernova 1987A, Astrophys. J. 340 (1989) 426–434.  doi:10.1086/167404.
  * (97) K. H. Lo, M. Coleman Miller, S. Bhattacharyya _et al._ , Determining neutron star masses and radii using energy-resolved waveforms of X-ray burst oscillations, Astrophys. J. 776 (2013) 19.  doi:10.1088/0004-637X/776/1/19.  [Erratum: Astrophys.J. 854, 187 (2018)], arXiv:1801.08031.
  * (98) Z. Lyu, N. Jiang and K. Yagi, Constraints on Einstein-dilation-Gauss-Bonnet gravity from black hole-neutron star gravitational wave events, Phys. Rev. D 105 (2022) 064001.  doi:10.1103/PhysRevD.105.064001.  [Erratum: Phys.Rev.D 106, 069901 (2022), Erratum: Phys.Rev.D 106, 069901 (2022)], arXiv:2201.02543.
  * (99) E. Maggio, H. O. Silva, A. Buonanno _et al._ , Tests of general relativity in the nonlinear regime: A parametrized plunge-merger-ringdown gravitational waveform model, Phys. Rev. D 108 (2023) 024043.  doi:10.1103/PhysRevD.108.024043.  arXiv:2212.09655.
  * (100) A. Maselli, V. Cardoso, V. Ferrari _et al._ , Equation-of-state-independent relations in neutron stars, Phys. Rev. D 88 (2013) 023007.  doi:10.1103/PhysRevD.88.023007.  arXiv:1304.2052.
  * (101) R. F. P. Mendes, Possibility of setting a new constraint to scalar-tensor theories, Phys. Rev. D 91 (2015) 064024.  doi:10.1103/PhysRevD.91.064024.  arXiv:1412.6789.
  * (102) R. F. P. Mendes and N. Ortiz, Highly compact neutron stars in scalar-tensor theories of gravity: Spontaneous scalarization versus gravitational collapse, Phys. Rev. D 93 (2016) 124035.  doi:10.1103/PhysRevD.93.124035.  arXiv:1604.04175.
  * (103) R. R. Metsaev and A. A. Tseytlin, Order alpha-prime (Two Loop) Equivalence of the String Equations of Motion and the Sigma Model Weyl Invariance Conditions: Dependence on the Dilaton and the Antisymmetric Tensor, Nucl. Phys. B 293 (1987) 385–419.  doi:10.1016/0550-3213(87)90077-0.
  * (104) M. C. Miller and F. K. Lamb, Bounds on the compactness of neutron stars from brightness oscillations during x-ray bursts, Astrophys. J. Lett. 499 (1998) L37.  doi:10.1086/311335.  arXiv:astro-ph/9711325.
  * (105) M. C. Miller and F. K. Lamb, Determining Neutron Star Properties by Fitting Oblate-star Waveform Models to X-ray Burst Oscillations, Astrophys. J. 808 (2015) 31.  doi:10.1088/0004-637X/808/1/31.  arXiv:1407.2579.
  * (106) 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 (2019) L24.  doi:10.3847/2041-8213/ab50c5.  arXiv:1912.05705.
  * (107) M. C. Miller _et al._ , The Radius of PSR J0740+6620 from NICER and XMM-Newton Data, Astrophys. J. Lett. 918 (2021) L28.  doi:10.3847/2041-8213/ac089b.  arXiv:2105.06979.
  * (108) S. Morisaki and T. Suyama, Spontaneous scalarization with an extremely massive field and heavy neutron stars, Phys. Rev. D 96 (2017) 084026.  doi:10.1103/PhysRevD.96.084026.  arXiv:1707.02809.
  * (109) S. M. Morsink, D. A. Leahy, C. Cadeau _et al._ , The Oblate Schwarzschild Approximation for Light Curves of Rapidly Rotating Neutron Stars, Astrophys. J. 663 (2007) 1244–1251.  doi:10.1086/518648.  arXiv:astro-ph/0703123.
  * (110) Y. Nakamura, D. Kikuchi, K. Yamada _et al._ , Weakly-gravitating objects in dynamical Chern–Simons gravity and constraints with gravity probe B, Class. Quant. Grav. 36 (2019) 105006.  doi:10.1088/1361-6382/ab04c5.  arXiv:1810.13313.
  * (111) J. Nättilä and J. J. E. Kajava.  Fundamental physics with neutron stars (2022).  arXiv:2211.15721.
  * (112) J. Nättilä and P. Pihajoki, Radiation from rapidly rotating oblate neutron stars, Astron. Astrophys. 615 (2018) A50.  doi:10.1051/0004-6361/201630261.  arXiv:1709.07292.
  * (113) H. P. Nollert, H. Ruder, H. Herold _et al._ , The relativistic ’looks’ of a neutron star, Astron. Astrophys. 208 (1989) 153–156.
  * (114) J. Novak, Neutron star transition to strong scalar field state in tensor scalar gravity, Phys. Rev. D 58 (1998) 064019.  doi:10.1103/PhysRevD.58.064019.  arXiv:gr-qc/9806022.
  * (115) M. Oertel, M. Hempel, T. Klähn _et al._ , Equations of state for supernovae and compact stars, Rev. Mod. Phys. 89 (2017) 015007.  doi:10.1103/RevModPhys.89.015007.  arXiv:1610.03361.
  * (116) G. J. Olmo, D. Rubiera-Garcia and A. Wojnar, Stellar structure models in modified theories of gravity: Lessons and challenges, Phys. Rept. 876 (2020) 1–75.  doi:10.1016/j.physrep.2020.07.001.  arXiv:1912.05202.
  * (117) J. R. Oppenheimer and G. M. Volkoff, On massive neutron cores, Phys. Rev. 55 (1939) 374–381.  doi:10.1103/PhysRev.55.374.
  * (118) F. Ozel, Surface Emission from Neutron Stars and Implications for the Physics of their Interiors, Rept. Prog. Phys. 76 (2013) 016901.  doi:10.1088/0034-4885/76/1/016901.  arXiv:1210.0916.
  * (119) C. Palenzuela and S. L. Liebling, Constraining scalar-tensor theories of gravity from the most massive neutron stars, Phys. Rev. D 93 (2016) 044009\.  doi:10.1103/PhysRevD.93.044009.  arXiv:1510.03471.
  * (120) P. Pani and E. Berti, Slowly rotating neutron stars in scalar-tensor theories, Phys. Rev. D 90 (2014) 024025.  doi:10.1103/PhysRevD.90.024025.  arXiv:1405.4547.