I-Love-Q relations: from compact stars to 
black holes
Authors: Kent Yagi & Nicolás Yunes
This is an author-created, un-copyedited version of an article published in Classical and 
Quantum Gravity. IOP Publishing Ltd is not responsible for any errors or omissions in this 
version of the manuscript or any version derived from it. The Version of Record is available 
online at https://dx.doi.org/10.1088/0264-9381/33/9/095005.
Yagi, Kent , and Nicolas Yunes. "I-Love-Q relations: from compact stars to black holes." Classical 
and Quantum Gravity 33 (April 2016): 095005. DOI: 10.1088/0264-9381/33/9/095005
Made available through Montana State University’s ScholarWorks 
scholarworks.montana.edu 
I-Love-Q relations: from compact stars to
black holes
Kent Yagi1,2 and Nicolás Yunes2
1 Department of Physics, Princeton University, Princeton, NJ 08544, USA
2 eXtreme Gravity Institute, Department of Physics, Montana State University,
Bozeman, MT 59717, USA
Abstract
The relations between most observables associated with a compact star, such
as the mass and radius of a neutron star or a quark star, typically depend
strongly on their unknown internal structure. The recently discovered I-Love-
Q relations (between the moment of inertia, the tidal deformability and the
quadrupole moment) are however approximately insensitive to this structure.
These relations become exact for stationary black holes (BHs) in General
Relativity as shown by the no-hair theorems, mainly because BHs are vacuum
solutions with event horizons. In this paper, we take the first steps toward
studying how the approximate I-Love-Q relations become exact in the limit as
compact stars become BHs. To do so, we consider a toy model for compact
stars, i.e.incompressible stars with anisotropic pressure, which allows us to
model an equilibrium sequence of stars with ever increasing compactness that
approaches the BH limit arbitrarily closely. We numerically construct such a
sequence in the slow-rotation and in the small-tide approximations by
extending the Hartle–Thorne formalism, and then extract the I-Love-Q trio
from the asymptotic behavior of the metric tensor at spatial infinity. We find
that the I-Love-Q relations approach the BH limit in a nontrivial way, with the
quadrupole moment and the tidal deformability changing sign as the com-
pactness and the amount of anisotropy are increased. Through a generalization
of Maclaurin spheroids to anisotropic stars, we show that the multipole
moments also change sign in the Newtonian limit as the amount of anisotropy
is increased because the star becomes prolate. We also prove analytically that
the stellar moment of inertia reaches the BH limit as the compactness reaches a
critical BH value in the strongly anisotropic limit. Modeling the BH limit
through a sequence of anisotropic stars, however, can fail when considering
other theories of gravity. We calculate the scalar dipole charge and the
moment of inertia in a particular parity-violating modified theory and find that
these quantities do not tend to their BH counterparts as the anisotropic stellar
sequence approaches the BH limit.
Keywords: neutron stars, black holes, modified theories of gravity, relativity
(Some figures may appear in colour only in the online journal)
1. Introduction
A plethora of compact stars with masses between 1 M and 2 M and with radii of
approximately 12 km have been discovered through a variety of astrophysical observations
[1–4]. The limited accuracy of these observations, coupled to degeneracies in the obser-
vables with respect to different models for the nuclear physics at supranuclear densities
encoded in the equation of state (EoS), have prevented observations from elucidating the
internal structure of compact objects. For example, x-ray observations do not typically
allow us to confidently state whether the compact objects observed are standard neutron
stars [5–8], or hybrid stars with quark-gluon plasma cores [9, 10], or perhaps even strange
quark stars [11]. Future observations of compact objects could shed some light on this
problem, as the accuracy of the observations increases and more observables are
obtained [12, 13].
The extraction of information from these future observations is aided by the use of
approximately universal relations, i.e.relations between certain observables that are
roughly insensitive to the EoS [13–17]. For example, the moment of inertia I, the tidal
deformability l2 (or tidal Love number) and the (rotation-induced) quadrupole moment Q
satisfy relations (the so-called I-Love-Q relations) that are EoS insensitive to a few% level
[14, 15]. Such relations are useful to analytically break degeneracies in the models used to
extract information from X-ray and gravitational-wave observations of compact objects.
This information, in turn, allows us to better probe nuclear physics [13] and gravitational
physics [14, 15].
Similar universal relations exist among the multipole moments of compact stars [18–22],
i.e.the coefficients of a multipolar expansion of the gravitational field far from the compact
object. These no-hair like relations resemble the well-known, black hole (BH) no-hair rela-
tions of general relativity (GR) [23–29]. The latter state that all multipole moments of an
uncharged, stationary BH in GR can be prescribed only in terms of the first two (the BH mass
and spin). The no-hair like relations of compact stars differ from the BH ones in that the
former require knowledge of the first three stellar multipole moments to prescribe all higher
moments in a manner that is roughly insensitive to the underlying EoS.
But how are the approximate I-Love-Q and no-hair like relations for compact stars
related to those that hold for BHs exactly? One way to address this question is to carry out
simulations of compact stars that gravitationally collapse into BHs, extract the I-Love-Q and
multipole moments and study how the relations evolve dynamically. However, not only are
such simulations computationally expensive, but the machinery employed in the past would
no longer be useful, as it is valid only for stationary spacetimes, i.e.the Geroch–Hansen
multipole moments [30, 31] used for example in [14, 15, 18, 20] are not well-defined for non-
stationary spacetimes. One would have to employ a dynamical generalization of these
moments and develop a procedure to extract them from dynamical simulations.
A simpler way to gain some insight is to consider how the universal relations evolve in a
sequence of equilibrium stellar configurations4 of ever increasing compactness that approa-
ches the compactness of BHs arbitrarily closely. Such a sequence, however, cannot be
constructed from neutron star solutions with isotropic pressure, as used in the original I-Love-
Q [14, 15] and no-hair like relations [18, 20]; such stars have a maximum stellar compactness
(i.e.the ratio between the stellar mass and radius) that is well below the BH limit. An
alternative approach is to consider a sequence of anisotropic stars5 (see e.g. [33] for a review
of anisotropic stars), which, for example, in the Bowers and Liang (BL) model [34] can reach
BH compactnesses for incompressible stars in the strongly anisotropic limit.
Following this logic, in [35] we studied how the no-hair like relations for compact stars
approach the BH limit. We first showed that the stellar shape transitions from prolate to oblate
as the compactness is increased. We then showed that the multipole moments approach the
BH limit with a power-law scaling and that the no-hair like relations also approach the BH
limit in a very nontrivial way. In this paper we extend these investigations in a variety of ways
and clarify several points that were left out of the initial analysis.
First, in this paper we consider both slowly-rotating stars and tidally-deformed stars,
which allow us to study how the I-Love-Q relations approach the BH relations in the BH
limit. In [36], we constructed tidally-deformed or slowly-rotating, anisotropic compact stars to
third order in spin for various realistic EoSs. We here follow [36] but focus on incompressible
stars, as this allows us to construct an equilibrium sequence of anisotropic stars that
approaches the BH limit arbitrarily closely.
Second, we extend the analysis of [35] by carrying out analytic calculations in various
limits: (i) the weak-field limit, (ii) beyond the weak-field limit, (iii) the strong-field limit and
(iv) the strongly-anisotropic limit. In the first limit, we expand all equations in small com-
pactness and retain only the leading terms in the expansion. This leads to anisotropic stars
modeled as incompressible spheroids with arbitrary rotation that reduce to Maclaurin
spheroids [37–39] in the isotropic limit. When going beyond the weak-field limit, we retain
subleading terms in the small compactness expansion, which is equivalent to a post-Min-
kowskian (PM) expansion; we extend the work of [40] for isotropic stars to the anisotropic
case and derive the moment of inertia and tidal deformability. In the third limit, we expand all
equations about the maximum compactness allowed for incompressible stars, extending the
analysis of [41] to anisotropic stars and deriving the tidal deformability for specific choices of
the anisotropy parameter. In the fourth limit, we expand all equations about the maximum
anisotropy allowed by the BL model, analytically deriving the moment of inertia for
incompressible stars as a function of the compactness.
Third, we study whether an equilibrium sequence of anisotropic compact stars can be
used to study the BH limit of stellar observables in theories other than GR. As an example, we
work in dynamical Chern–Simons (dCS) gravity [42–44], a parity violating modified theory
of gravity that is motivated from string theory [45], loop quantum gravity [46–48] and
effective theories of inflation [49]. We treat this modified theory as an effective field theory
and assume that the GR deformation is small. Such a treatment ensures the well-posedness of
the initial value problem [50]. Slowly-rotating, anisotropic compact stars to linear order in
spin in dCS gravity were constructed in [36] using realistic EoSs within the anisotropy model
4 Another approach is to consider ‘BH mimickers’ whose compactness can reach that of BHs, such as the gravastars
considered in [32]. The latter, however, are very different from neutron stars or quark stars.
5 Anisotropic stars are here only used as a toy model to study an equilibrium sequence of compact stars that can
reach the BH limit, and not as a realistic model for compact stars.
proposed by Horvat et al [51]. We now extend the treatment in [36] to the BL anisotropy
model and focus on the incompressible case.
1.1. Executive summary
Let us now present a brief summary of our results. We find that the I-Love-Q relations for
strongly anisotropic stars in GR indeed approach the BH limit as one increases the com-
pactness. Figure 1 shows evidence for this by presenting the I-Love-Q relations for incom-
pressible stars with a variety of anisotropy parameterslBL in the BL model [34]. The isotropic
case is recovered when l = 0BL , while l p= -2BL corresponds to the strongly anisotropic
limit. The BH limit (l = 02,BH¯ ) corresponds to =I 4BH¯ and =Q 1BH¯ , shown with dashed
horizontal lines. We confirm the validity of our numerical results by comparing them to an
analytic calculation of the I-Love relations in the PM approximation (solid curves in the top
panel of figure 1). Observe that the relations approach the BH limit as the compactness is
increased (shown with arrows) in a way that depends quite strongly on lBL, with l2¯ and Q¯
changing sign as the BH limit is approached.
We also find that the approach of the I-Love-Q relations to the BH limit appears to be
continuous, as shown in figure 1. That is, we find no evidence of the discontinuity
Figure 1. Relations between the dimensionless moment of inertia *ºI I M3¯ and the
dimensionless tidal deformability *l lº M2 2 5¯ (top), and between the dimensionless
quadrupole moment *cº -Q Q M3 2¯ ( ) and l2¯ (bottom) for an equilibrium sequence of
anisotropic, incompressible stars with varying compactness (the arrows indicate
increasing compactness), given some anisotropy parameter lBL, to leading-order in
slow rotation and tidal deformation. *M is the stellar mass, *c º J M 2 is the
dimensionless spin parameter, with J the magnitude of the stellar spin angular
momentum. Isotropic stars correspond to l = 0BL , while strongly-anisotropic stars
correspond to l p= -2BL . Our numerical results are validated by analytic PM
calculations (solid curves). The dashed horizontal lines correspond to the BH values of
I¯ and Q¯, while l = 02,BH¯ is the BH value for the dimensionless tidal deformability.
Observe that the I-Love-Q relations of anisotropic stars approach the BH limit
continuously.
hypothesized in [52], based on a weak-field calculation of the quadrupole moment of strongly
anisotropic, incompressible stars. We in fact prove analytically that the moment of inertia of a
strongly anisotropic, incompressible compact star reaches the BH limit continuously as the
compactness is increased. We do so by constructing slowly-rotating, anisotropic incom-
pressible stars to linear order in spin in the strongly anisotropic limit (l p= -2BL ) and
analytically deriving I¯ as a function of the compactness C in terms of hypergeometric
functions. Taylor expanding I¯ about  c= +C 1 2BH 2( ), with χ the dimensionless spin
parameter, we find that  c= + -I C I C C ,BH BH 2¯ ( ) ¯ ( ).
The quadrupole moment changes sign as it approaches the BH limit, as shown in figure 1,
but is this the case for all multipole moments? We find that this is not the case by constructing
incompressible spheroids with anisotropic pressure and arbitrary rotation in the weak-field
limit. We derive a necessary condition on the anisotropy model such that spheroidal con-
figurations are realized and find that the BL model satisfies such a condition. We then
calculate the ℓth mass and current multipole moments, Mℓ and Sℓ , in the slow-rotation limit
within the BL model and find that +M ℓ2 2 and +S ℓ2 3 are both proportional to
p l+ +1 4 5 ℓBL 1( ) . This means that the sign of only (M2, M6, M10...) and (S3, S7, S11...) is
opposite to that of the isotropic case when l p< -4 5BL , which is consistent with the results
of [52] for the quadrupole moment M2. In particular, the sign of (M4, M8, M12...) and (S5, S9,
S13...) is the same as that of the sign of the isotropic case even in the strongly-anisotropic
limit.
Although the I-Love-Q relations for compact stars approach the BH limit as one increases
the compactness in GR, we find that this is not always the case in other theories of gravity
Figure 2. Dimensionless scalar dipole charge m¯ (equation (72)) (top) and the dCS
correction to the dimensionless moment of inertia dI¯ (normalized by the dimensionless
dCS coupling constant and the GR value of I¯ ) (equation (73)) (bottom) for a sequence
of anisotropic, incompressible stars labeled by stellar compactness C and anisotropy
parameter lBL. Corresponding BH values are shown by black crosses. Our numerical
results are validated by analytic PM calculations (solid curves) (equation (B.12)) in the
top panel. Observe that, unlike in the GR case, m¯ and dI¯ do not approach the BH limit
as one decreases lBL and increases C.
when the limit is modeled through an equilibrium sequence of anisotropic stars. Figure 2
presents evidence for this by showing the scalar dipole charge and the correction to the
dimensionless moment of inertia in dCS gravity as a function of the stellar compactness. Once
more, we validate our numerical results by comparing them to analytic PM relations for the
scalar dipole charge. Observe that unlike in the GR case, these quantities do not approach the
dCS BH limit (shown with black crosses) as one increases the compactness. This result
suggests that modeling the BH limit through strongly anisotropic stars is not appropriate in
certain modified theories of gravity.
The remainder of this paper presents the details of the calculations that led to the results
summarized above. In section 2, we explain the formalism that we use to construct slowly-
rotating and tidally-deformed anisotropic stars. We also describe the BL anisotropy model
and show how the maximum stellar compactness for a non-rotating configuration approaches
the BH one in the strongly anisotropic limit. In section 3, we present analytic calculations of
the stellar moment of inertia, tidal deformability and multipole moments in certain limits. In
section 4, we present numerical results that show how the I-Love-Q relations approach the
BH limit in GR. We also show that the scalar dipole charge and the correction to the moment
of inertia in dCS gravity do not approach the BH limit. Finally, in section 5, we give a short
summary and discuss various avenues for future work. We use the geometric units of
= =c G1 throughout this paper.
2. Formalism and anisotropy model
In this section, we first explain the formalism we use to construct slowly-rotating or tidally-
deformed compact stars with anisotropic pressure and extract the stellar multipole moments
and tidal deformability. We then explain the specific anisotropic model that we will use
throughout the paper. We present the spherically-symmetric background solution and
describe the maximum compactness such a solution can possess for polytropic EoSs of the
form r= +p K n1 1 . Here p and ρ are the stellar radial pressure and energy density, while K
and n are constants. Henceforth, the stellar compactness is defined by * *ºC M R , where *M
and *R are the stellar mass and radius for a non-rotating configuration respectively.
2.1. Formalism
Let us first explain how one can construct slowly-rotating compact stars with anisotropic
pressure, by following [36, 53–55] and extending the Hartle–Thorne approach [56, 57] to
third order in spin. Let us assume the spacetime is stationary and axisymmetric, such that the
metric can be written as
⎡
⎣⎢
⎤
⎦⎥

  

q q
q q q f w q q
=- + + + -
+ + + - W - +
+
n ls h r t m r
r M r
r
r k r r w r t
d e 1 2 , d e 1
2 ,
2
d
1 2 , d sin d , , d
, 1
r r2 2 2
2
2
2 2 2 2 2 2
4
[ ( )] ( )
( )
[ ( )]( { [ ( ) ( )] } )
( ) ( )
( ) ( )
where ν and λ are functions of the radial coordinate r only, while ω, h, k, m and w are
functions of both r and θ. The quantity ò is a book-keeping parameter that labels the order of
an expression in *WM( ), where Ω is the spin angular velocity. The surface is defined as the
location where the radial pressure vanishes. We transform the radial coordinate via
 q x q= + +r R R R, , , 22 4( ) ( ) ( ) ( )
so that the spin perturbation to the radial pressure and density vanish throughout the star
[56, 57]. The enclosed mass function M(r) is defined via
º -l- M r
r
e 1
2
, 3r
( ) ( )( )
and thus, *M is the value ofM(r) evaluated at the stellar surface *R . We decompose ω, h, k, m,
ξ and w in Legendre polynomials [36].
The stress–energy tensor for matter with anisotropic pressure can be written as
[36, 58, 59]
r= + + Pmn m n m n mnT u u p k k q , 4( )
where q is the tangential pressure and mu is the fluid four-velocity, given by
= Wmu u u, 0, 0,0 0( ), with u0 determined through the normalization condition
= -m mu u 1. mk is a unit radial vector that is spacelike ( =m mk k 1) and orthogonal to the
four-velocity ( =m mk u 0) of the fluid, while P º + -mn mn m n m ng u u k k is a projection
operator onto a two-surface orthogonal to mu and mk . We introduce the anisotropy parameter
s º -p q [58, 59] with s = 0 corresponding to isotropic matter. Following the treatment of
metric perturbations, we expand σ in the slow-rotation approximation and decompose each
term in Legendre polynomials as
 s q s s s q= + + +R R R R P, cos . 500 2 02 22 2 4( ) ( ) { ( ) ( ) ( )} ( ) ( )( ) ( ) ( )
Notice that the superscript (subscript) in sℓn( ) corresponds to the order of the spin (Legendre)
decomposition. The function s00( ) needs to be specified a priori, and it in fact defines the
anisotropy model. The function s22( ) is determined consistently by solving the perturbed
Einstein equations, once s00( ) is chosen [36]. The function s02( ) is irrelevant in this paper as it
only affects the stellar mass at subleading order in a small spin expansion.
We construct slowly-rotating compact star solutions with anisotropic pressure as follows.
First, we substitute the metric ansatz and the matter stress–energy tensor mentioned above
into the Einstein equations. We then expand in small spin (or equivalently in ò) and solve the
perturbed Einstein equations order by order in ò. In the interior region, we solve the equations
numerically with a regularity condition at the center. In the exterior region, we solve the
equations analytically with an asymptotic flatness condition at spatial infinity. We finally
match the two solutions at the stellar surface to determine any integration constants. The latter
determine the moment of inertia I, the quadrupole moment Q and the octupole moment S3 of
the exterior solution at linear, quadratic and third order in spin respectively.
In this paper, we also construct non-rotating but tidally-deformed compact stars to extract
the stellar tidal deformability [41, 60]. We are here particularly interested in the quadrupolar,
electric-type tidal deformability, l2, which is defined as the ratio of the tidally-induced
quadrupole moment and the external tidal field strength. We follow [41, 60, 61] and treat tidal
deformations as small perturbations of an isolated compact star solution. Such a tidally-
deformed compact star can be constructed similarly to how we construct slowly-rotating
solution, except that we set wW = = =w 0, as we are only interested in electric-type, even-
parity perturbations.
For convenience, we work with the following dimensionless quantities throughout:
* * * *c c
l lº º - º - ºI I
M
Q
Q
M
S
S
M M
, , , . 6
3 3 2 3
3
4 3 2
2
5
¯ ¯ ¯ ¯ ( )
Here, the dimensionless spin parameter χ is defined through the magnitude of the spin
angular momentum J by *c º J M
2, with J only kept to ( ). The BH value of each
dimensionless quantity above is =I 4BH¯ [62], =Q 1BH¯ [31], =S 13,BH¯ [31] and l = 02,BH¯
[29, 41, 60, 63, 64]. We here choose to work with the above choice of normalization
introduced in [14, 15, 18–20], but clearly this choice is not unique. In fact, one can choose
other normalizations, for example involving the stellar compactness, which may improve the
universality in the I-Love-Q relations and no-hair like relations for compact stars among
stellar multipole moments [22]. Other choices of normalization, nonetheless, will not affect
the conclusions we arrive at in this paper.
2.2. Anisotropy model
Let us now describe the specific anisotropy model that we use in this paper. Following BL
[34], we choose
⎜ ⎟⎛⎝
⎞
⎠s
l r r= + + -
-
R p p
M
R
R
3
3 1
2
. 70
0 BL
1
2( ) ( )( ) ( )( )
Here, lBL is a constant parameter that characterizes the amount of anisotropy. Isotropic
pressure corresponds to l = 0BL , since then both s00( ) and s22( ) vanish. The particular form of
s00( ) in equation (7) was proposed such that the Tolman–Oppenheimer–Volkoff equation
could be solved analytically for a spherically-symmetric, incompressible (polytropic index
n = 0, i.e. r = const.) anisotropic star to yield [34]
*
=M R C
R
R , 8
2
3( ) ( )
*
r p=R
C
R
3
4
, 9
2
( ) ( )
*
*
*
p= -
- - -
- - -
g g
g gp R
C
R
C CR R
C CR R
3
4
1 2 1 2
3 1 2 1 2
, 10
2
2 2
2 2
( ) ( ) ( )
( ) ( )
( )
⎡
⎣⎢
⎤
⎦⎥*n g=
- - -g g
R
C CR R1
ln
3 1 2 1 2
2
, 11
2 2
( ) ( ) ( ) ( )
where γ is defined by
⎜ ⎟⎛⎝
⎞
⎠g
l
pº +
1
2
1
2
. 12BL ( )
How do the maximum compactness Cmax of anisotropic stars differ from the isotropic
case? One can find Cmax for anisotropic stars by finding the value of C for which the central
radial pressure =p R 0( ) diverges [34]:
= - g-C 1
2
1 3 . 13max 1( ) ( )
Clearly, the solution in equations (10) and(11) is only well-defined for l p-2BL , such that
C 0max . Such a condition on lBL also ensures that p 0 and the solution does not diverge
when C Cmax . The red solid curve in figure 3 presents Cmax as a function of lBL for
incompressible stars (equation (13)). Observe that in the isotropic incompressible case,
= »C 4 9 0.444 ...max , while Cmax approaches 1/2, the compactness of a non-rotating BH,
in the l p -2BL limit. This is exactly why we consider anisotropic stars in this paper, as
they allow us to construct a sequence of equilibrium stars that approaches the BH limit
arbitrarily closely. Henceforth, ºC 1 2BH , the compactness of a non-rotating BH, since we
work to leading order in the slow rotation approximation and  c= +C 1 2Kerr 2( ). For
reference, figure 3 also shows the maximum compactness for anisotropic stars with an n=1
polytropic EoS (blue dashed curve), which is always smaller than that of incompressible stars.
Although the maximum compactness of non-rotating, anisotropic stars can reach the
compactness of non-rotating BHs in the strongly anisotropic limit (l p -2BL ), the causal
structure inside such a star is quite different from that of a BH. The left, middle and right
Figure 3. Maximum compactness realized for the n=0 (solid) and n=1 (dashed)
polytopes as a function of the anisotropy parameter lBL. The horizontal dashed line
corresponds to the compactness of a non-rotating BH.
Figure 4. Causal structure of a non-rotating compact anisotropic star with C = 0.3 (left)
and C = 0.5 (middle) with l p= -2BL and a non-rotating BH (right). Blue and red
curves correspond to ingoing and outgoing null geodesics respectively. The opening
angle between these curves at each crossing point shows that of a light cone at each
point. Observe that the surface of an anisotropic star with C = 0.5 is a trapped surface
and radiation inside the star cannot escape to outside. Observe also how the causal
structure of the interior region for such a star is different from that of a BH.
panels of figure 4 show the causal structure of non-rotating anisotropic compact stars with
C = 0.3 and C = 0.5, and that of a BH respectively. To construct these panels, we introduce a
new (retarded) time coordinate = -T v R [65], where *= +v t r is a null coordinate with r*
the tortoise coordinate in the exterior and interior regions, given by equations (A.1) and(A.2)
respectively. The ingoing null geodesics (blue lines in the figure) are given by =v const.,
while the outgoing null geodesics (red curves in the figure) are given by *- =t r const. The
opening angle between blue and red curves at each point represents that of the light cone,
while the stellar surface or the event horizon are denoted by a black dashed vertical line.
Observe that a photon emitted inside a star with C = 0.3 can escape out to spatial infinity,
while that from a star with C = 0.5 cannot. This is because the surface for the latter acts as a
trapped surface, just like the event horizon of a BH. However, notice that the causal structure
in the interior region between an anisotropic compact star with C = 0.5 and a BH is different.
In particular, a photon emitted inside an anisotropic star with C = 0.5 stays at a constant
radius R, while the one emitted inside a BH eventually falls into singularity.
3. Analytic calculations
Before diving into a full numerical analysis, let us first present some analytic calculations of
the multipole moments, of I¯ and of l2¯ for incompressible, anisotropic stars in certain limits.
Later on, in section 4, we will use these analytic calculations to verify the validity of our
numerical analysis. In section 3.1, we calculate multipole moments with arbitrary ℓ in the
weak-field (or so-called ‘Newtonian’) limit by constructing a spheroid with arbitrary rotation
that reduces to Maclaurin spheroids [37–39] in the isotropic limit. In section 3.2, we calculate
I¯ and l2¯ using a PM analysis, which is valid beyond the weak-field limit. In section 3.3, we
derive l2¯ in the strong-field, or maximum-compactness limit, for specific choices of lBL,
while in section 3.4 we calculate I¯ as a function of C in the strongly anisotropic limit,
l p= -2BL¯ . In the latter, we prove that I¯ approaches the moment of inertia of a non-spinning
BH in the limit as the compactness goes to 1/2.
3.1. Weak-field limit
We here derive the multipole moments of anisotropic, incompressible stars in the weak-field
limit. We begin by constructing an anisotropic stellar solution that is spheroidal and valid to
arbitrary order in rotation. In the isotropic case, such a solution reduces to Maclaurin
spheroids [37–39]. We then use the anisotropic spheroidal solution to find the multipole
moments of the star. We conclude this subsection by providing a phenomenological expla-
nation for why rotating strongly-anisotropic stars are prolate in the Newtonian limit, and
whether such anisotropic stars are stable to perturbations in the amount of anisotropy.
3.1.1. Maclaurin-like spheroids. Let us first prove that σ is purely a function of r in
axisymmetry, working in spherical coordinates q fr, ,( ). The r and θ components of the
hydrostatic equilibrium equation are given by [52]
r r
s¶
¶ +
¶F
¶ + =
p
r r r
1 1 2
0, 14( )
r q q r
s
q
¶
¶ +
F
¶ -
¶
¶ =
p1 d 1
0, 15( )
where F = F + FG c with FG and Fc representing the gravitational and centrifugal potentials
respectively. We decompose s= FB p, ,( ) using Legendre polynomials Pℓ as
å q=B B r P cos . 16
ℓ
ℓ ℓ( ) ( ) ( )
Substituting this into equations (14) and(15), one finds
r r
s+ F + =p
r r r
ℓ
1 d
d
d
d
1 2
0 0 , 17ℓ ℓ ℓ ( ) ( )
r r s+ F - = >p ℓ
1 1
0 0 , 18ℓ ℓ ℓ ( ) ( )
while equation (15) is automatically satisfied when =ℓ 0. Taking a derivative of equation (18)
with respect to r and combining this with equation (17), one finds
s s+ = >
r r
ℓ
d
d
2
0 0 . 19ℓ ℓ ( ) ( )
Imposing regularity at the stellar center, one finds s = 0ℓ for >ℓ 0. This is consistent with
[52], in which the authors showed that s = 02 at quadratic order in spin. Since the only non-
vanishing Legendre mode of σ is the =ℓ 0 mode, σ cannot depend on any angular
coordinates for stationary and axisymmetric, incompressible and anisotropic stars in the
weak-field limit.
Let us now derive a necessary condition for s r0 ( ) such that the star is a spheroid. This
condition comes from the equations of structure that determine the gravitational potential, FG,
and the radial pressure p. The Poisson equation determines the former, which is not modified
from its form in the isotropic case:
⎛
⎝⎜
⎞
⎠⎟ pr
¶
¶ +
¶
¶ +
¶
¶ F =x y z 4 , 20
2
2
2
2
2
2 G
( )
working in Cartesian coordinates x y z, ,( ), with the z-axis identified with the axis of rotation.
The hydrostatic equilibrium equation determines p x y z, ,( ), and with s s= r0 ( ) its
components in the x and z directions are
r r
s¶
¶ = -
¶F
¶ - + W
p
x x r
n x
1 1 2
, 21xG 0 2 ( )
r r
s¶
¶ = -
¶F
¶ -
p
z z r
n
1 1 2
, 22zG 0 ( )
where = = + +xr x y z2 2 2∣ ∣ and =n x ri i . The second term on the right-hand side in
equations (21) and(22) correspond to an extra force induced by pressure anisotropy.
We now solve the above set of differential equations to find a condition on s r0 ( ).
Equation (20) can be easily solved using Green’s function as in the isotropic case, and the
solution for a spheroid is given by [38, 39]
prF = - - + -x y z A a A x y A z, , . 23G 0 12 1 2 2 3 2( ) [ ( ) ] ( )
Here, A0, A1 and A3 are given in terms of the stellar eccentricity º -e a a12 32 12( ) by
= -A e e
e
2 1 arcsin
, 240
2 ( ) ( )
= - - -A e e
e
e
e
1 arcsin 1
, 251
2
3
2
2
( ) ( )
= - -A
e
e e
e
2
2
1 arcsin
, 263 2
2
3
( ) ( )
where the stellar radius on the x- (or y)- and z-axes are denoted by a1 and a3 respectively.
Note that < <e0 12 when a star is oblate, while <e 02 when a star is prolate, in which case
a1 (a3) do not correspond to the semi-major (semi-minor) radius of the star. Imposing that the
solution is spheroidal implies that the radial pressure p must be given by [38]
⎛
⎝⎜
⎞
⎠⎟= -
+ -p x y z p x y
a
z
a
, , 1 , 27c
2 2
1
2
2
3
2
( ) ( )
where pc is the central radial pressure. Substituting equations (23) and(27) into (22) and(21),
one finds that s0 must satisfy
⎛
⎝⎜
⎞
⎠⎟s p r= -
p
a
A r , 280
c
3
2 3
2 2 ( )
while W2 must satisfy
pr
rW =
- + -A A a p a a p
a a
2 2 2
. 292
2
1 3 3
2
c 1
2
3
2
c
1
2
3
2
[ ( ) ] ( )
Equation (28) shows that σ needs to be proportional to r2 to realize a spheroid, since recall
that for an incompressible fluid r = const.
3.1.2. Spheroidal shape and multipole moments. Consider now a specific choice of s0 that is
consistent with the BL model and that satisfies equation (28):
s l r= r
3
. 300
BL 2 2 ( )
Note that if one transforms equation (30) to the R coordinate with equations (2), (30) agrees
with equation (7) in the Newtonian limit at 0th order in spin. The radial part of the =ℓ 2
contribution at second-order in spin agrees with s R22 ( )( ) in equation (5) in the Newtonian
limit. Combining equations (30) with (28), one finds
⎜ ⎟⎛⎝
⎞
⎠p
l r= +p A a
3
. 31c 3
BL 2
3
2 ( )
Substituting this into equation (29), one finds the relation between Ω and e2 as
p p
lW = W
-
- - +
e
A e A e
3
2 1
1
3
, 322 K
2
2
1
2
3
BL 2{ }[ ( ) ] ( )
where
* p rW º = -M
a
e
4
3
1 33K
2
1
3
2 ( )
corresponds to the (squared) Keplerian angular velocity at the equatorial surface.
The left panel of figure 5 shows W W2 K2 as a function of e2¯ in equation (32) for various
values of lBL. Here, e2¯ is defined as
⎧
⎨⎪
⎩⎪

=
- =
- = <-
e
e a a
a a
1 ,
1 ,
34
a
a
a
a
e
e
2
2
1 3
1 1 3
3
2
1
2
1
2
3
2
2
2
¯
( )
( )
( )
such that the star is prolate (oblate) for - < <e1 02¯ ( < <e0 12¯ ). The right panel of
figure 5 shows regions in the lBL–e2¯ plane in which W2 is positive (i.e. the stellar solution
exists). Observe that when l 0BL the stellar shape can only be oblate, while when
l p-0.85BL it can only be prolate. Observe also that the shape can be either oblate or
prolate when p l- <0.85 0BL , but one of the branches is probably unstable.
We are now in a position to derive the multipole moments for such an anisotropic
spheroidal solution, following the analysis in [19]. The mass and current multipole moments
for a Newtonian incompressible star with arbitrary rotation are given by
= -+ + W+
+ + +M
ℓ ℓ
e a
1 3
2 3 2 5
, 35ℓ
ℓ
ℓ ℓ
2 2
1
K
2 2 2
1
2 5( )
( )( )
( )
= -+ + W W+
+S
ℓ ℓ
e a
1 6
2 3 2 5
. 36ℓ
ℓ
ℓ ℓ
2 1 K
2 2
1
2 5( )
( )( )
( )
One could numerically invert equation (32) and insert the solution into the equations above in
order to express the multipole moments entirely in terms of the angular velocity.
Let us instead consider the slow-rotation approximation, as this allows us to invert
equation (32) analytically. One can perturbatively solve equation (32) for e2 within the slow-
rotation approximation to find
⎛
⎝⎜
⎞
⎠⎟
⎡
⎣
⎢⎢
⎛
⎝⎜
⎞
⎠⎟
⎤
⎦
⎥⎥
p
p l= +
W
W +
W
We
10
4 5
. 372
BL K
2
K
4
( )
Notice that the star is oblate ( >e 02 ) when l p> -4 5BL , while it is prolate ( <e 02 ) when
l p< -4 5BL . Notice, however, that there is a range oflBL in which one cannot construct an
equilibrium configuration within the small-rotation approximation, because then >e 12 ,
Figure 5. (Left) W W2 K2 (equation (32)) as a function of e2¯ (equation (34)) for various
lBL. The solid and dashed curves correspond to the oblate and prolate configurations
respectively. The dotted–dashed curve can have both oblate and prolate configurations.
(Right) Regions in the lBL–e2¯ plane where W 02 and the stellar solution exists
(shaded). Observe that when l > 0BL ( l p0.85BL ) only oblate (prolate) configura-
tions exist. Both prolate and oblate branches exist when l pÎ -0.85 , 0BL ( ).
which would force one of the semi-axes to be imaginary. This occurs in the range
p l l- < < =4 5 eBL BL 12( ), where we have defined
⎛
⎝⎜
⎞
⎠⎟l
p= - - WW
= 4
5
1
5
2
. 38eBL
1
2
K
2
2 ( )( )
Notice also that when -e 12 (l p+ 4 5 1BL∣ ∣ and l p< -4 5BL ), the star remains
prolate and the semi-axes can remain real.
Let us now evaluate the multipole moments in the slow-rotation approximation.
Substituting equations (37) into (35) and(36) and eliminating ρ using WK2 in equation (33),
one finds
⎛
⎝⎜
⎞
⎠⎟
⎡
⎣
⎢⎢
⎛
⎝⎜
⎞
⎠⎟
⎤
⎦
⎥⎥
p
p l=
-
+ + +
W
W W +
W
W+
+ +
+
+
+
+
M
ℓ ℓ
a
1 3 10
2 3 2 5 4 5
,
39
ℓ
ℓ ℓ
ℓ
ℓ
K
ℓ
ℓ
2 2
1 1
BL
1
K
2 2
2
1
2 5
K
2 4( ) ( )
( )( )( )
( )
⎛
⎝⎜
⎞
⎠⎟
⎡
⎣
⎢⎢
⎛
⎝⎜
⎞
⎠⎟
⎤
⎦
⎥⎥
p
p l=
-
+ + +
W
W W +
W
W+
+
+
+
S
ℓ ℓ
a
1 6 10
2 3 2 5 4 5
. 40ℓ
ℓ ℓ
ℓ
ℓ
ℓ
ℓ
2 1
BL K
2 1
K
3
1
2 5
K
2 3( ) ( )
( )( )( )
( )
Since p lµ ++ - +M 4 5ℓ ℓ2 2 BL 1( ) ( ), the sign of +M ℓ2 2 is opposite to that in the isotropic case
for ℓ even, i.e.the moments M2, M6, M10..., or simply +M ℓ4 2 with integer ℓ, flip sign when
l p< -4 5BL . Similarly, since p lµ ++ -S 4 5ℓ ℓ2 1 BL( ) , the sign of +S ℓ2 1 is opposite to that
in the isotropic case when ℓ is odd, i.e.the moments S3, S7, S11..., or simply +S ℓ4 3 with integer
ℓ, flip sign when l p< -4 5BL . In the quadrupole case, p lµ + -M 4 52 BL 1( ) , which is
consistent with [52].
Before proceeding, notice that the multipole moments in equations (39) and(40) diverge
atl p= -4 5BL . This divergence originates from equation (37) and is an artifact of the slow-
rotation approximation. Such an approximation breaks down near the divergence, and such a
feature is absent in the multipole moments valid for arbitrary rotation in equations (35)
and(36). For example, the left panel of figure 5 shows that e2 (or e2¯ ) is finite for a given
W W2 K2 at l p= -4 5BL , and hence, the multipole moments in equations (35) and(36) are
also finite at l p= -4 5BL . We will see how different manifestations of this divergence
contaminate our calculations later on; fortunately, it does not affect the behavior of the
multipole moments near the BH limit, since the latter requires we take the l p -2BL limit,
which is far from the divergent value of lBL.
3.1.3. From oblate to prolate anisotropic stars. Let us now try to develop a better
understanding of why the stellar shape changes from oblate to prolate as one decreaseslBL. In
particular, the goal of this subsection is to understand the behavior of equation (37) from a
force balance argument on a fluid element inside an isolated star.
We begin by defining the sum of the pressure gradient and potential gradient of a fluid
element (normalized by WK2 for later convenience) acting along the x-axis as
⎛
⎝⎜
⎞
⎠⎟rº -W
¶
¶ +
¶F
¶F
p
x x
1 1
, 41x
K
2
G ( )
and decompose it within the small eccentricity approximation ( e 12∣ ∣ ) as
å d=
=
F F , 42x
k
x
k k
0
2 ( )( )
where δ is a book keeping parameter that counts orders in e.
Let us now look at force balance along the x-axis at second-order in rotation. At this
order, the second term on the right hand side of equation (21) vanishes (as it is spin
independent), and hence the force balance equation is given by
+ + WW =lF F x 0, 43x x,0
1
,
1
2
K
2BL
( )( ) ( )
where Fx,0
1( ) corresponds to Fx
1( ) for an isotropic star while lFx,
1
BL
( ) represents its anisotropic
correction for a fixed e2. From the exact solution of p and FG given in equations (23)
and(27), one finds
l
p= - = -lF x e F x e
2
5
,
2
. 44x x,0
1 2
,
1 BL 2
BL
( )( ) ( )
Notice that the direction of the above forces changes depending on the sign of e2,
i.e.depending on whether the star is oblate or prolate. Notice also that equation (43) states
that the centrifugal force, which acts always away from the stellar center, needs to balance the
sum of Fx,0
1( ) and lFx,
1
BL
( ) , which from equation (44) is given by
p l
p+ = -
+
lF F xe
4 5
10
. 45x x,0
1
,
1 BL 2
BL
( )( ) ( )
With these expressions at hand, let us now discuss how force balance and the stellar
shape change as one varies lBL, with a fixed W WK. Consider first the isotropic case. Setting
l = 0BL in equation (44), the anisotropic correction to the pressure gradient force vanishes.
Since Fx,0
1( ) needs to balance the centrifugal force, from equations (43) and(44), one finds that
>e 02 , which implies the star must be oblate.
Let us now imagine we were to add a small, negative amount of anisotropy (l  1BL∣ ∣
and l < 0BL ) to an isotropic configuration. The additional anisotropy force, lFx,1 BL( ) , no longer
vanishes, and in fact, it must be centrifugal, as depicted on the top of figure 6, since >e 02
for the background, isotropic configuration. In order for balance to be restored, the star must
compensate for this additional force, and the only way to do so is for Fx,0
1( ) to increase, which
implies e2 must also increase and the star becomes more oblate.
Imagine now that we increased the magnitude of lBL further, always with l < 0BL . As
mentioned above, e2 is forced to increase as lBL becomes more negative to maintain
equilibrium, forcing the star to become more and more oblate. Eventually, lBL reaches the
critical value l =eBL 12( ) of equation (38) where =e 12 and the star would seem to become
infinitely oblate. Notice, however, that the small eccentricity approximation of equation (42)
becomes invalid near this critical point. If one further decreases lBL, one enters into a
forbidden and unphysical region p l l- < < =4 5 eBL BL 12( )( ) , where e2 becomes larger than
unity, which is an artifact of the slow-rotation approximation.
Consider now decreasing lBL even further, such that p l p- < < -2 4 5BL .
Equation (45) shows that the total pressure and potential gradient force cannot balance the
centrifugal force when l p< -4 5BL unless e2 flips sign, i.e.unless the star becomes prolate
as shown in the bottom panel of figure 6. Moreover, this equation also shows that as lBL
decreases the magnitude of the total pressure and potential gradient force increases for a fixed
e2 (and x). This implies that e2∣ ∣ must also decrease to keep the sum of forces balanced against
the centrifugal force, making the star less prolate aslBL decreases. Thus, the change in sign of
the pressure and potential gradient force and its anisotropy correction is responsible for the
change in stellar shape.
Let us conclude with a brief discussion of the stability of strongly-anisotropic, rotating
stars by investigating how the fluid elements shift and whether the stellar configuration
approaches an equilibrium configuration as one varies lBL. Notice that since equation (37)
only determines the ratio between a1 and a3, for simplicity, we will consider the =a const3
case. This choice allows us to determine how e2 varies purely from the change in a1.
Let us first look at the stability of an oblate, rotating anisotropic star by considering
adding a finite amount of negative anisotropy to an isotropic configuration (l = 0BL ). The
additional force due to anisotropy shifts the fluid elements in the centrifugal direction (see the
top of figure 6). This increases a1, and thus increasing e
2 and making the star more oblate. As
we argued before, the equilibrium configuration becomes more oblate as the amount of
anisotropy is increased. Thus, the anisotropy force pushes the star toward the equilibrium
configuration, making them stable to perturbations in lBL.
Let us now study the stability of a prolate, rotating star by considering adding a small
amount of negative anisotropy to an equilibrium, prolate configuration with some value of
Figure 6. A schematic picture that shows a free body diagram for a fluid element in an
anisotropic star at second order in spin on the equatorial plane. W W xK 2( ) is the
centrifugal force, while Fx,0
1( ) is the sum of the pressure gradient and gravitational force
with isotropic pressure and lFx,
1
BL
( ) corresponds to its anisotropic correction. The top
diagram corresponds to the case with l l> =eBL BL 12( ) and < <e0 12 . Since lFx,1 BL( ) acts
outward, the oblateness of the stellar shape increases as one decreases lBL. Once the
anisotropy parameter reaches the critical value l =eBL 12( ) (defined in equation (38)), the
star becomes infinitely oblate. One cannot construct an equilibrium configuration with
l l p< < -= 4 5eBL 1 BL2( ) within the slow-rotation approximation. If one further
decreases lBL, an equilibrium oblate configuration does not exist anymore and the
stellar shape changes to prolate, whose force balance is shown in the bottom diagram
for p l p- < < -2 4 5BL and <e 02 .
lBL in p l p- < < -2 4 5BL . The force due to additional anisotropy now shifts the fluid
elements in the centripetal direction for a fixed e2 (see the top of figure 6). This decreases a1,
and thus e2∣ ∣ increases, making the star even more prolate. However, we argued before that the
equilibrium configuration of prolate stars becomes less prolate as anisotropy is increased.
Therefore, the anisotropy force pushes the star away from the prolate equilibrium
configuration. This suggests that prolate anisotropic stars are unstable to perturbations
in lBL6.
3.2. Beyond the weak-field limit
We here derive an analytic expression for the I-Love relation by perturbing about the weak-
field limit in a PM approximation, i.e.in an expansion about small compactness. We extend
[40] for isotropic stars to anisotropic stars, starting from the exact solution to the Einstein
equations for incompressible, non-rotating, anisotropic stars with the BL model given in
equations (8)–(11) [34] as our background.
Let us begin by finding the metric perturbation at linear order in spin. To achieve this, we
impose the following ansatz for the metric perturbation and the moment of inertia:
å=
=
A A C , 46
k
k k
0
( )( )
where A is either w1 or I, with w R1( ) the l=1 mode of the ft,( ) component of the metric
perturbation w qR,( ) in a Legendre decomposition. We substitute this ansatz for w1, together
with equations (8)–(11), into the Einstein equations and perturb about C=0. We then solve
the perturbed Einstein equations order by order in C in the interior region, with regularity
imposed at the center. We also substitute the above ansatz for I in the exterior solution and
perturb about C=0. Finally, we match the perturbed interior and exterior solutions order by
order in C at the stellar surface to calculate I within the PM approximation.
Through this procedure, we find that I¯ is given by
⎜ ⎟
⎡
⎣⎢
⎛
⎝
⎞
⎠
⎤
⎦⎥åå
l
p= + += =
I
C
c C C
2
5
1
1 , 47
i j
i
ij
I
j
i
2
1
6
1
BL 7¯ ( ) ( )( ¯)
where the coefficients cij
I( ¯) are explicitly given in table B1. Notice that the leading term in C
does not depend on lBL, which shows that the I¯ –C relation is unaffected by anisotropy in the
weak-field limit.
One can derive the PM expression for l2 in a similar manner, namely by perturbing the
metric and the Einstein equations, solving the perturbed equations order by order in C, and
carrying out a matching calculation at the stellar surface. Doing so, we find
⎜ ⎟ ⎜ ⎟
⎪ ⎪
⎪ ⎪
⎧
⎨
⎩
⎛
⎝
⎞
⎠
⎡
⎣⎢
⎛
⎝
⎞
⎠
⎤
⎦⎥
⎫
⎬
⎭
åål pp l
l
p
l
p= + + + +
l
= =
+ -
C
c
C
C
2
4 5
1
1 1
5
4 4
, 48
i j
i
ij
j i
2
BL
5
1
6
1
2
BL BL
1
72¯ ( ) ( )( ¯ )
where the coefficients lcij 2(
¯ ) are given in table B2.
We can now invert equation (48) perturbatively (in small C or large l2¯ ) to obtain C as a
function of l2¯ . Inserting such an expression into equation (47), we find the I-Love relation:
6 Reference [52] discusses the instability of prolate anisotropic stars in the context of tidal deformations, though
further study is needed to understand how that instability relates to the ones discussed here (if at all).
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎪
⎪
⎡
⎣⎢
⎛
⎝
⎞
⎠
⎤
⎦⎥
⎧
⎨
⎩
⎛
⎝
⎞
⎠
⎡
⎣⎢
⎛
⎝
⎞
⎠
⎤
⎦⎥

åålp l
l
p
l
p l
l
= + + +
+
l
= =
+ -
-
I c
2
5
2 1
5
4
1
1
2
1
5
4
1
. 49
i j
i
ij
I
j i
BL
2
2 5
1
6
1
2
BL
11
BL
6
2
5
2
7 5
2¯ ¯
¯
( ¯ )} ( )
( ¯ ¯ )
Here, the coefficients lcijI 2(
¯ ¯ ) are given in table B3. Equations (47)–(49) with l = 0BL agree
with those found in [40] in the isotropic limit.
3.3. Strong-field limit
Let us now calculate the tidal deformability l2¯ for incompressible, anisotropic stars in the
strong-field limit for specific choices of lBL. The strong-field limit here refers to the max-
imum-compactness limit (not to be confused with the BH limit), i.e.the limit in which
C Cmax . Following the analytic strong-field analysis for isotropic stars of [41], we intro-
duce the new radial coordinate
⎛
⎝⎜
⎞
⎠⎟*
º -
g
x
C R
R
1 2 , 50
2
2
( )
where γ is given by equation (12). Notice that the stellar center corresponds to x=1, while
the stellar surface corresponds to =x 1 3 for a maximum compactness configuration (see
equation (13)).
With this coordinate choice, the radial pressure p and the metric function ν in
equations (10) and(11) simplify to
*
p= -
- -
- -
g
gp x
C C x
R C x
3
4
1 2
3 1 2
, 51
2
( ) [( ) ]
[ ( ) ]
( )
n g=
- - -g
x
C xln 3 1 2 ln 2
. 52( ) [ ( ) ] ( ) ( )
Moreover, this coordinate choice allows us to decouple the differential equation for the metric
perturbation h R2 ( ) (that corresponds to the =ℓ 2 mode of qh R,( ) in equation (1)) from other
metric perturbations, leading to a homogeneous, third-order differential equation.
Before solving for the tidal deformability of anisotropic stars, let us review how this is
done analytically in the isotropic case. The third-order differential equation for h2 can be
integrated once, reducing it to a second-order equation. Reference [41] solved this equation
for a maximum compactness configuration ( =C 4 9max ) both in the interior and in the
exterior regions, imposing regularity at the center. The authors then matched the interior
solution to the exterior solution at the surface, with a certain jump condition due to the
discontinuous density at the surface. Doing so, they found the dimensionless tidal deform-
ability l lc2 BLmax¯ ( ) to be
l = - 0
72
5 308 81 ln 3
0.0658, 53c2
max¯ ( )
( )
( )
for an isotropic, incompressible star of maximum compactness.
One can carry out a similar analysis completely analytically for certain choices of lBL.
For example, when l p= 2BL , the third-order differential equation for h2 for a maximum
compactness configuration is given by
- + - -
+ - + - - =
x x
h
x
x x
h
x
x x
h
x
x h
4 1
d
d
6 1 1 4
d
d
19 46 27
d
d
5 3 0, 54
x x
x
x
3
3
2
3
2
2
2
2
2 2
2
( ) ( ) ( )
( ) ( ) ( )
where ºh x h R xx2 2( ) [ ( )]. Imposing regularity at the center (x = 1), the solution to the above
equation is
= - - -h x
C x
x
2 exp 5 arctanh
5 1
, 55x h2 ( )
( )
( )
( )
where Ch is an integration constant, from which one can find
*
* *
*
º ¢ = -y R h R R
h R
15 1. 562
2
( ) ( )
( )
( )
We find that one does not need to worry about the density discontinuity at the surface when
l p= 2BL , because the solution is not discontinuous at the surface for that value of lBL.
Hence, one can use equation (23) in [61] with y given by equation (56) to find the tidal Love
number k2. The latter can be turned into the tidal deformability l2 as explained in [61] to find
l p = + - 2
48
180 8 15 135 ln 3
0.766. 57c2
max¯ ( ) ( )
Notice that the tidal deformability is much larger in the anisotropic case with l p= 2BL than
in the isotropic one.
3.4. Strongly anisotropic limit
Let us now analytically investigate the strongly anisotropic limit, l p -2BL , for incom-
pressible anisotropic compact stars of arbitrary compactness. The goal of this subsection is to
prove that =I 4¯ in such a limit as C 1 2, which agrees with the expected result for BHs.
To achieve this goal, we analytically construct slowly-rotating anisotropic incompressible
compact stars to linear order in spin.
Let us first discuss the (background) 0th-order in spin solution. From equation (10), one
finds that the radial pressure vanishes in the strongly anisotropic limit [34, 52]. Taking the
limit of l p -2BL in equation (11), one further finds that
⎛
⎝⎜
⎞
⎠⎟*
n = - - -R C C R
R
3
2
ln 1 2
1
2
ln 1 2 , 58int
2
2
( ) ( ) ( )( )
where the superscript (int) reminds us that this is the solution in the interior region.
Let us now find the interior solution at linear order in spin. Substituting the above
background solution into the differential equation for w1 (the =ℓ 1 mode of w qR,( ) in
equation (1) of [36]), one finds
*
*
*
*
w w w+ -- +
-
- =R
CR R
R CR R R
C
CR R
CR R
d
d
11 4
2
d
d
6
3 2
2
0. 59
2
1
int
2
2 2
2 2
1
int 2 2
2 2 2 1
int
( ) ( )
( )
( ) ( )
( )
Imposing regularity at the center, we find an analytic solution for w1 in the interior region:
⎛
⎝⎜
⎞
⎠⎟* *
w = -wR C R CR F C R
R
2
3
2
,
9
4
;
5
2
; 2 , 601
int 2 2 3 4
2 1
2
21
( ) ( ) ( )( )
where wC 1 is an integration constant and F , ; ;2 1(· · · ·) represents hypergeometric functions.
Matching the above interior solution to the exterior solution, w = W -R I R1 21ext 3( ) ( ),
at the stellar surface with the conditions * *w w=R R1int 1ext( ) ( )( ) ( ) and * *w w¢ = ¢R R1int 1ext( ) ( )( ) ( ) ,
one finds [35]
å
å=
=
=
I C
C
a C C
b C C
1
2
, 61i
i
i
j i
j2
0
1
0
2
¯ ( )
( )
( )
( )
where the coefficients of the numerator are
⎜ ⎟ ⎜ ⎟⎛⎝
⎞
⎠
⎛
⎝
⎞
⎠= - +a C F C F C9
5
2
,
13
4
;
7
2
; 2 5
3
2
,
9
4
;
5
2
; 2 , 620 2 1 2 1( ) ( )
⎜ ⎟⎛⎝
⎞
⎠=a C F C18
5
2
,
13
4
;
7
2
; 2 , 631 2 1( ) ( )
while the coefficients of the denominator are
⎜ ⎟⎛⎝
⎞
⎠= -b C F C5
3
2
,
9
4
;
5
2
; 2 , 640 2 1( ) ( )
⎜ ⎟ ⎜ ⎟⎛⎝
⎞
⎠
⎛
⎝
⎞
⎠= - +b C F C F C9
5
2
,
13
4
;
7
2
; 2 15
3
2
,
9
4
;
5
2
; 2 , 651 2 1 2 1( ) ( )
⎜ ⎟⎛⎝
⎞
⎠=b C F C18
5
2
,
13
4
;
7
2
; 2 . 662 2 1( ) ( )
We can now Taylor expand the above expression about =C CBH. For a slowly-rotating BH,
 c= +C 1 2BH 2( ), where χ is the dimensionless spin parameter. With this, one then finds
that
 c= - - + - + -I C C C C C C C4 40 224 , . 67BH BH 2 BH 3 2¯ ( ) ( ) ( ) [( ) ] ( )
In particular, in the C CBH limit, one finds that =I 1 2 4¯ ( ) , which agrees with the BH
result. Therefore, the above analysis analytically proves that the moment of inertia for
strongly anisotropic compact stars approaches the BH limit exactly as the compactness goes
to the compactness of a non-spinning BH.
4. Numerical results
In this section, we investigate numerically how I¯ , l2¯ , Q¯ and S3¯ approach the BH limit as one
increases the stellar compactness. We first consider these quantities for anisotropic incom-
pressible stars in GR, and then consider stars in dCS gravity, exploring both the scalar dipole
charge and the dCS correction to the moment of inertia as a function of C. The results of the
previous section are used here to validate our numerical calculations. This section thus shows
explicitly how the I-Love-Q relations and the multipole moments approach the BH limit using
an equilibrium sequence of strongly anisotropic stars of ever increasing compactness.
4.1. Approaching the BH limit in GR
Let us first investigate how Q¯ and l2¯ depend on the compactness within the range that is
realistic for neutron stars and quark stars. The left panel of figure 7 presents these quantities as
a function of compactness within  C0.1 0.3. Observe that Q¯ and l2¯ are always positive
for incompressible stars with isotropic pressure or with l p> -0.8BL , where the latter
corresponds to the critical value in the weak-field limit [52]. On the other hand, when
l p= -BL , Q¯ and l2¯ are at first negative when C is small, and then they diverge at around
C = 0.26. For even more anisotropic configurations, such as l p= -1.5BL , Q¯ and l2¯ are
always negative when C 0.3 and they do not diverge.
Let us find all the values of compactness for which l2¯ diverges. We can do so analy-
tically through the (3,3)-Padé resummation in equation (B.2) of the Taylor series of l2¯ in
equation (48). The bottom, left panel of figure 7 plots this Padé resummation7, which agrees
very well with the numerical results, validating the latter. We can also find the location of the
divergencies numerically by calculatingl2¯ as a function of C. Doing so, we typically find that
the deformability diverges at two values of compactness: a low or weak-field value and a high
or strong-field value. The values of compactness for which the tidal deformability diverges
are shown in the right panel of figure 7. Observe that the Padé resummation can only recover
the weak-field branch of the divergences, becoming highly inaccurate when C 0.35.
Observe also that l2¯ diverges only when  p l p- -1.14 0.8BL for  C0 0.43 and
when  p l p- -1.14 0.9BL for  C C0.43 max . In particular, these divergences are
absent in the BH limit, i.e.as l p -2BL and C CBH.
What is the physical meaning of these divergences? The right panel of figure 7 shows
that there exists a range of lBL for which the tidal deformability diverges for all values of
compactness, and the same would be true of the divergences in the quadrupole moment. This
suggests that we can write the location of the divergences as l = lf CQBLdiv , 2 ( )¯ ¯ , for some
Figure 7. (Left) Compactness dependence of Q¯ (top) and l2¯ (bottom) for
incompressible stars in the compactness range ÎC 0.1, 0.3( ) for various anisotropy
parameters. We also show the Padé-resummed, analytic Love-C relation for various
values of lBL (equation (B.2)) with solid curves in the bottom panel. (Right) Stellar
compactness at which the Love-C relation diverges as a function of lBL, as derived
from the Padé resummation of the tidal deformability (solid). The dashed red branch
( l -1.14BL ) is an artifact of the Padé resummation. Points are obtained from
numerical calculations in the weak-field (red crosses) and strong-field (blue pluses)
regimes. Observe that divergences occur only for  p l p- -1.14 0.8BL in the
weak-field regime and  p l p- -1.14 0.9BL in the strong-field regime, and the two
branches merge as lBL becomes more negative. Observe also that the PM
approximation becomes invalid for C 0.35. We also show the maximum
compactness in figure 3 as a black dotted–dashed curve, which crosses the blue
points at l p~ 0.9BL .
7 We do not show the Padé resummation when l p= -0.8BL , since then the PM Taylor series diverges.
functions lf CQ, 2 ( )¯ ¯ defined only for ÎC C0, max( ). In the PM regime, we can asymptotically
expand this function about zero compactness
åp= - +l l
=
¥
f C a C
4
5
. 68Q
n
n
Q n
,
wf
1
,
2
2( ) ( )¯ ¯ (
¯ ¯ )
Similarly, in the strong-field regime, we can expand about C bmax0
( ) with l p= » -b 9 10BL 0
(where the strong-field branch in blue crosses the maximum compactness curve in black in the
right panel of figure 7)
å= + -l l
=
¥
f C b b C C . 69Q
n
n
Q b n
,
sf
0
1
,
max2
2 0( ) ( ) ( )¯ ¯ (
¯ ¯ ) ( )
In both cases we have factored out the weak-field and the strong-field limits. Recall from
section 3.1 that in the weak-field limit ( C 0), the divergence in the multipole moments at
l p= -4 5BL was shown to be an artifact of the slow-rotation approximation (see
equation (37)). Therefore, the weak-field divergences captured by the asymptotic expansion
in equation (68) are also due to the slow-rotation approximation in the Q¯ case and the small-
tide approximation in the l2¯ case. Moreover, since lfQ,wf 2¯ ¯ and lfQ,sf 2¯ ¯ are two different
asymptotic representations of the same functions lf CQ, 2 ( )¯ ¯ , the strong-field divergences
captured by the asymptotic expansion in equation (69), and in fact all of the divergences
captured by lf CQ, 2 ( )¯ ¯ , are due to the use of these approximations. Obviously, these
approximations become invalid around the region where Q¯ and l2¯ diverge, since then the
metric perturbations become larger than the background and higher-order contributions can
no longer be neglected. These divergences would be absent if we had found solutions without
the slow-rotation or small-tide approximations.
Let us now study the compactness dependence of I¯ , l2¯ , Q¯ and S3¯ in the strong-field
regime. Figure 8 presents these quantities as a function of the stellar compactness in the range
 C0.46 0.5. Observe that I¯ (l2¯ ) monotonically decreases (increases) to approach the BH
limit. Observe also that the analytic I–C relation in the strongly anisotropic limit
l p= -2BL( , see equation (61)) validate the numerical results, approximating the latter very
Figure 8. Compactness dependence of I¯ (top left), Q¯ (top right), l2¯ (bottom left) and S3¯
(bottom right) for the strongly anisotropic, incompressible stars in the strongly
relativistic regime. The solid black curve in the top left panel corresponds to the
analytic relation (equation (61)) for l p= -2BL . Dashed lines on the right panel
corresponds to the BH value for Q¯ and S3¯ . Observe how each observable approaches
the BH result (black cross) as one increases the compactness.
closely even whenl p= -1.9BL . On the other hand, the behavior of Q¯ and S3¯ as the BH limit
is approached is quite different for some values oflBL. Indeed, these quantities first overshoot
the BH values, and then decrease toward the BH limit for l p p= - -1.2 , 1.4BL , and p-1.5 .
When l p= -1.9BL , Q¯ and S3¯ both increase monotonically to approach the BH limit.
In order to quantify how the properties of anisotropic compact stars approach those of
BHs, we take the difference between different observables in the strong-field limit
( C Cmax ) and the corresponding BH values. Figure 9 plots this difference as a function of
lBL, together with the analytic results of equations (53) and(57) for l C2 max¯ ( ) at l = 0BL and
p2 (green dots) to validate our numerics. Observe how rapidly each quantity approaches the
BH limit as l p -2BL , or equivalently, as C 1 2max , since recall that C 1 2max only
when l p -2BL (figure 3). Observe also that l C2 max¯ ( ), Q Cmax¯ ( ) and S C3 max¯ ( ) diverge at
l p~ -0.9 ;BL this value of lBL corresponds to that where the strong-field branch in blue
crosses the maximum compactness curve in black in the right panel of figure 7, and, as we
discussed before, they are artifacts of the slow-rotation and small-tide approximations.
Let us now investigate how fast I¯ approaches the BH limit as one takes the C 1 2
limit. In order to quantify this, we define the scaling exponent of I¯ as
tº
D
k
Id ln
d ln
, 70I
¯
( )¯
where
t º - D º -C C
C
I
I I
I
, . 71BH
BH
BH
BH
¯ ¯ ¯
¯ ( )
Figure 9. (Top) Difference between the properties (I¯ , Q¯, l2¯ and S3¯ ) of anisotropic,
incompressible stars at their maximum compactness Cmax and the properties of a non-
rotating BH, plotted as a function of l pBL . The black dashed horizontal line
corresponds to the BH result. Green dots represent analytic values of l2¯ at
l p= 0, 2BL , given in equations (53) and(57). (Bottom) Same as the top panel but
the absolute difference in a log plot. Observe how rapidly each observable approaches
the BH result as one takes the l p -2BL limit.
One can explicitly calculate kI¯ with l p= -2BL from equation (61), which we show in
figure 10. Observe that kI¯ approaches 2 as the BH limit is approached (t  0). In fact, one
can expand kI¯ as calculated from equation (61) around t = 0 to find analytically that
 t= +k 2I 1 4( )¯ . In [36], we also looked at the scaling exponents of Q¯ and S3¯ and found
that they roughly approach the BH limit linearly and quadratically respectively.
Let us finally look at the interrelations among these observables, i.e.the I-Love-Q
relations. Figure 1 already presented the I-Love and Q-Love relations for an equilibrium
sequence of anisotropic incompressible stars with increasing compactness (indicated by the
arrows) for various choices of lBL. Notice that we plot the absolute value of l2¯ and Q¯, since
these quantities can be both positive and negative. Observe that the anisotropic I-Love
relation is qualitatively very similar to the isotropic one, with both I¯ and l2∣ ¯ ∣ monotonically
decreasing as one increases the compactness, irrespective of lBL. On the other hand, the
strongly anisotropic Q-Love relation is quite different from the isotropic one, with Q∣ ¯∣
decreasing to zero at l ~ 102∣ ¯ ∣ , but then starting to increase again towards the BH limit. This
behavior is due to Q¯ changing sign at this l2∣ ¯ ∣ whenl p p= - -1.5 , 1.9BL . In the top panel of
figure 1, we also show the analytic, (3, 3)-Padé approximation to the I-Love relation of
section 3.2, which validates our numerical results, deviating from it only in the strongly
relativistic regime where the PM approximation loses accuracy.
4.2. Approaching the BH limit in dCS gravity
Let us now investigate whether the stellar quantities approach the BH values as one increases
the compactness even in theories other than GR, taking dCS gravity [42–44] as an example.
This theory is motivated from the theoretical requirement to cancel anomalies in heterotic
string theory [45], the inclusion of matter when scalarizing the Barbero-Immirzi parameter in
loop quantum gravity [46–48] and from effective theories of inflation [49].
The dCS action modifies the Einstein–Hilbert action by adding a (pseudo-) scalar field
with a canonical kinetic term, a ϑ-dependent potential and an interaction term in the form of
the product of the Pontryagin density and the scalar field. dCS gravity breaks parity at the
level of the field equations, in the sense that it introduces modifications to GR only in
Figure 10. Scaling exponent of I¯ as a function of τ with l p= -2BL , obtained
analytically using equation (61). Observe that the exponent approaches 2 (black
dashed) as one approaches the BH limit of t  0.
spacetimes that break spherical symmetry, like that of rotating compact stars. The strength of
dCS modifications are proportional to its coupling parameter α, which multiplies the inter-
action term. Since higher curvature corrections to the action are currently unknown, we treat
the theory as an effective field theory and keep terms up to leading order in the coupling
constant  a( ( ) in the scalar field and  a2( ) in the metric). This avoids problems with the
initial value formulation that arise when one insists in treating the theory as exact [50].
In this paper, we focus on how the dimensionless scalar dipole charge m¯ and the dCS
correction to I¯ approach the BH limit in the strongly relativistic regime. Following [36], we
construct a slowly-rotating anisotropic incompressible star to linear order in spin in dCS
gravity with a vanishing scalar field potential. We extract m¯ from the asymptotic behavior of
the scalar field at spatial infinity [66]:
⎛
⎝⎜
⎞
⎠⎟
*J q am c q= +R C
R
M
R
,
5
8
cos
. 723
2
3
3
( ) ¯ ( )
Imposing regularity at the horizon, one finds the BH limit of the dimensionless scalar charge
m = 8BH¯ [66]. One obtains the dCS correction to I¯ by looking at the asymptotic behavior of
the ft,( ) component of the metric. We introduce dI¯ as [36]
*d x=I
M I
I
, 73
4
CS
CS
GR
¯ ¯
¯ ( )
where I GR¯ and I CS¯ are the GR and dCS contributions to I¯ , while x paº 16CS 2. The moment
of inertia for a BH in dCS gravity is given by
= WI
S
M
, 74BH
1,BH
BH BH
3
¯ ( )
where MBH and WBH correspond to the BH mass and the BH angular velocity at the horizon
respectively. The latter is given by
⎛
⎝⎜
⎞
⎠⎟
xW = W -
M
1
709
7168
, 75BH Kerr
CS
BH
4
( )
where WKerr is the BH angular velocity for the Kerr solution. From equations (73)–(75), one
finds that the dCS corrections to the moment of inertia for a BH is
d = »I 709
1792
0.396. 76BH¯ ( )
Figure 2 shows m¯ and dI¯ as a function of C in dCS gravity for various values of lBL. The
BH limit in each panel is shown as a black cross. Observe that, unlike in the case in GR, m¯
and dI¯ do not approach the BH limit, but instead approach a different point. In order to
confirm these numerical calculations, we derived analytic, (2, 2)-Padé resummed relations
between m¯ and C within the PM approximation, as given by equation (B.12) (solid curves in
the top panel of figure 2). Observe that such analytic relations also show the feature that m¯
Table 1. Scaling exponents for dI¯ as a function of τ for incompressible anisotropic stars
in dCS gravity.
lBL p-1.2 p-1.5 p-1.9
dk I¯ −0.256 −0.311 −0.355
does not approach the BH limit in the limit C 1 2. Since dI¯ depends on m¯, one would not
expect dI¯ to approach the BH limit either, given that m¯ does not.
Following 4.1, let us finally study how the relation between dI¯ and C in dCS gravity
behaves near the =C 1 2 critical point. We estimate the scaling exponent dk I¯ of such a
relation for various values of lBL, which we show in table 1. Observe that the exponents are
negative, unlike in the relation between I¯ and C in GR. This means that dI¯ diverges in the
limit t  0 (or C 1 2), which is consistent with the bottom panel of figure 2. The BH
values for m¯ and dI¯ are obtained by imposing regularity at the horizon. Since m¯ and dI¯ for an
anisotropic compact star do not approach the BH values as one takes C 1 2, such
quantities are not regular at the surface and diverge in the limit. In this sense, such a solution
is unphysical and it is also indicative of the breakdown of the small coupling approximation.
This failure to arrive at the BH limit suggests a failure in the description of the BH limit as an
equilibrium sequence of anisotropic compact stars in dCS gravity.
5. Future directions
We have studied how the I-Love-Q relations for compact stars approach the BH limit as one
increases the stellar compactness by considering an equilibrium sequence of slowly-rotating/
tidally-deformed, incompressible stars with anisotropic pressure. We found that this sequence
approaches the BH limit in a nontrivial way, similar to the way the no-hair like relations
approach the BH limit [35]. We have also calculated the scalar dipole charge and moment of
inertia for incompressible, anisotropic stars in dCS gravity. We found that, unlike in GR,
these quantities do not approach the BH limit as one increases the stellar compactness, and in
fact, the metric diverges in this limit. These findings suggest that whether one can use a
sequence of anisotropic compact stars as a toy model to probe how the I-Love-Q relations
approach the BH limit depends on the underlying gravitational theory.
We have also carried out analytic calculations in various limits to validate and elucidate
our numerical results. In the strongly anisotropic limit, we derived the moment of inertia as a
function of the compactness and proved analytically that it reaches the BH value in the BH
limit. In the weak-field limit, we constructed incompressible, anisotropic spheroids with
arbitrary rotation that reduce to Maclaurin spheroids in the isotropic limit. We used such
spheroids to explain why strongly anisotropic, rotating stars become prolate in the weak-field
limit. We also derived analytic PM expressions for the moment of inertia and tidal deform-
ability, which accurately reproduce the numerical results for C 0.35.
Our analysis can be extended in various directions. A natural extension would be to study
how higher-order multipole moments approach the BH limit. Instead of constructing slowly-
rotating solutions, one can construct rapidly-rotating ones using e.g.the RNS code [67] and
extract higher order multipole moments. One can also calculate higher-order ( ℓ 3) tidal
deformabilities for anisotropic stars by extending the isotropic analysis of [41, 60, 68] and
study how they approach the BH limit. One should easily be able to apply some of the
analytic techniques presented here to higher-order tidal deformabilities, like the PM analysis
in section 3.2 and the calculation in the strong-field limit in section 3.3.
Another avenue for future work includes performing a stability analysis of Maclaurin-
like spheroids for anisotropic stars. Figure 5 shows that multiple values of e2 are allowed for a
given value ofW WK for a fixed value oflBL. One can study whether one of these branches is
unstable to perturbations, as we suggested in section 3.1.3 (see also [52]). One can also
construct Jacobi-like ellipsoids by relaxing axisymmetry. Then, one can again carry out a
stability analysis to see whether Maclaurin-like or Jacobi-like spheroids are energetically
favored.
Yet, another natural extension is to study universal relations for anisotropic stars in non-
GR theories other than dCS gravity. For example, Silva et al [59] constructed slowly-rotating,
anisotropic neutron star solutions to linear order in spin in scalar-tensor theories [69, 70]. One
could repeat their analysis in the strongly anisotropic regime, extract the scalar monopole
charge and the moment of inertia, and study whether these quantities approach the BH limit as
one increases the stellar compactness.
In this paper, we followed [40] and constructed Padé resummed expressions for I¯ andl2¯ .
One can extend such an analysis to Q¯ and S3¯. Our preliminary results show that the (3, 3)-
Padé approximant for Q¯ cannot capture our numerical results as well as the Padé resummation
of the PM expansion of I¯ and l2¯ even in the isotropic case. One may need to include higher-
order terms in C and construct higher-order approximants, or alternatively, one may need to
use other resummation methods or consider other expansions beyond the weak-field one.
We considered anisotropic stars as a toy model, whose compactness can reach the BH
value in the strongly anisotropic limit. An important future task would be to consider a more
realistic situation: the gravitational collapse of compact stars into BHs. One could then try to
dynamically monitor the I-Love-Q and no-hair like relations in such time-dependent situa-
tions to see how they compare to the results found in this paper and in [35]. Since the Geroch–
Hansen multipole moments [30, 31] used here and in [35] are only defined for stationary
spacetimes, one may first need to develop a generalization that is valid in dynamical situa-
tions, and yet reduces to the Geroch–Hansen moments in stationarity. One would also need to
consider how to extract such generalized moments from numerical gravitational collapse
calculations.
A final avenue for future work could be the study of a possible connection between
transitions from compact stars to BHs and second-order phase transitions in condensed matter
physics [71]. In [35], we showed that the scaling exponents for the moment of inertia (or the
current dipole moment), the mass quadrupole moment and the current octupole moment are
EoS universal to 10%( ) for isotropic stars. Such a feature may have some analog with the
universality of critical exponents in second-order phase transitions. Using the anti-de Sitter
(AdS)/conformal field theory (CFT) correspondence, [72, 73] showed that transitions from
non-rotating compact stars to BHs in AdS spacetimes correspond to phase transitions from
high density, baryonic states to thermal quark-gluon plasma states on the CFT side.
Having said this, whether this analogy can be firmly established remains to be seen. This
is because the exponents found in [35] are positive, while those in second-order phase
transitions are negative. Moreover, nobody has yet shown that scale invariance exists in the
collapse of realistic compact objects to BHs8. If this were the case, one would then have to
uncover what the analogy for the correlation length in the gravity sector is. More detailed
investigations would thus be necessary to further elucidate this intriguing line of study.
Acknowledgments
We would like to thank Steven Gubser, Jim Lattimer, Bennett Link and Anton B. Vorontsov
for useful comments, suggestions and advice. KYacknowledges support from JSPS Post-
doctoral Fellowships for Research Abroad and NSF grant PHY-1305682. NY acknowledges
8 However, see e.g. [74, 75] for scale invariance, universality and critical phenomena arising in the context of critical
gravitational collapse of scalar fields.
support from NSF CAREER Grant PHY-1250636. Some calculations used the computer
algebra-systems MAPLE, in combination with the GRTensorII package [76].
Appendix A. Tortoise coordinates for anisotropic compact stars
In this appendix, we derive the radial tortoise coordinate for non-rotating, incompressible and
anisotropic stars with l p= -2BL , which we use in section 2.2 to study their causal structure
(see figure 4). We achieve this by transforming our coordinate system to Eddington–Fin-
kelstein coordinates. The causal structure in the exterior region of an anisotropic compact star
is the same as that of a Schwarzschild BH. One introduces the null coordinate *= +v t r
ext,
where *r
ext is the radial tortoise coordinate in the exterior region. We change coordinates from
t to v in the metric and impose that *= - =g g r Rd d 1vR tt
ext( ) to find
⎛
⎝⎜
⎞
⎠⎟** * * *ò= - = + -
r
R R
R
g
R
R
C
C
R
R
1 d
2 ln
1
2
1 . A.1
tt
ext
( )
Regarding the interior region, one introduces another null coordinate *= +v t r
int, where
*r
int is the radial tortoise coordinate in the interior region. In this case, one needs to transform
not only the coordinate t to v but also the coordinate f to ψ via y f= + r¯ , where r¯ is a
function of R and θ. Such a coordinate transformation is similar to that from the Boyer–
Lindquist coordinates to Eddington–Finkelstein coordinates in the Kerr metric. Imposing
=g 1vR (and =g 0RR to further determine r¯ ), one finds
⎜ ⎟⎛⎝
⎞
⎠
⎡
⎣
⎢⎢
⎛
⎝⎜
⎞
⎠⎟
⎤
⎦⎥
*
* * *
*
= - - + - + - -
´ -- -
r
R C
C
C C C
CR
R
C
R
R
C
R
R
C
1
1
2 1 2
2 log
1
2
1
1
2 2 1 2
2
1 2
sin 2 sin 2 , A.2
int
3 2
2
2
1 1
( ) ( )
( ) ( )
where we used equation (11) for the (t, t) component of the metric and determined the
integration constant such that * * * *=r R r R
int ext( ) ( ).
Appendix B. Tables and Padé approximants for the PM analysis
In this appendix, we show some of the coefficients in the PM expressions of section 3.2, and
explain how one can construct Padé approximants following [40]. The constants cij
I( ¯), lcij 2(
¯ ) and
Table B1. Coefficients cij
I( ¯) in equation (47) for I¯ as a function of the stellar com-
pactness C within the PM approximation.
c I1,0
( ¯) c I1,1
( ¯) c I2,0
( ¯) c I2,1
( ¯) c I2,2
( ¯) c I3,0
( ¯) c I3,1
( ¯) c I3,2
( ¯) c I3,3
( ¯)
6
7
- 3
28
106
105
- 5
28
- 1
42
316
231
- 298
1155
- 155
1848
- 1
168
c I4,0
( ¯) c I4,1
( ¯) c I4,2
( ¯) c I4,3
( ¯) c I4,4
( ¯) c I5,0
( ¯) c I5,1
( ¯) c I5,2
( ¯) c I5,3
( ¯)
351872
175175
- 7225
21021
- 13747
64680
- 787
24024
- 5
3003
125632
40425
- 636578
1576575
- 2131
4550
- 421523
3603600
c I5,4
( ¯) c I5,5
( ¯) c I6,0
( ¯) c I6,1
( ¯) c I6,2
( ¯) c I6,3
( ¯) c I6,4
( ¯) c I6,5
( ¯) c I6,6
( ¯)
- 601
48048
- 1
1980
60771136
12182625
- 9147988
26801775
- 8480897
8933925
- 146991919
428828400
- 8218907
142942800
- 4207
875160
- 71
437580
lcijI 2(
¯ ¯ ) in equations (47)–(49) in GR are given in tables B1–B3 respectively. Notice that
equations (47)–(49) all have the form
åa a= + +
=
y x x x1 , B.1k
i
i
i
0
1
6
7[ ( )] ( )
where x=C for equations (47) and(48), while l= -x 2 1 5¯ for equation (49). From
equation (B.1), one can construct the (3, 3)-Padé approximant given by
⎡
⎣
⎢⎢
⎤
⎦
⎥⎥
å
åa
b
b
= +=
=
y x
x
x
x , B.2k i
i
i
i i
i
0
0
3
1
3
0
3
2
3
7( ) ( )
( )
( )
with the following identification of constants
b a a a a a a a a a a a= - - + -2 , B.3103 1 3 5 1 42 22 5 2 3 4 33 ( )( )
b a a a a a a a a a a a a a
a a a a a a a a a
= - + - + - - -
+ - - +
2
, B.4
11
3
3 5 4
2
1
2
3
3
2 4 6 3 5 2
2
4 1
2
2
6 2 3 5 2 4
2
3
2
4
( ) [ ( ) ( )]
( )
( )
b a a a a a a a a a a a a a a a
a a a a a a a a a a a a a a a
= - + + - + + -
- + + - - + + -
2
2 , B.5
12
3
3 6 4 5 1
2
2
2
6 2 4
2
3
2
4 4 6 5
2
1
2
3
5 2
2
3 4 3
3
3 6 4 5 2 3
2
5 3 4
2
( ) ( )
( ) ( )
( )
b a a a a a a a a a a a a a a a
a a a a a a a a a a a a
a a a a a a
= + - - + - + -
+ - + + + - -
+ - + B.6
2 3 2 2
2 2 2
,
13
3
6 2
3
3 5 4
2
2
2
3
2
4 1 3 6 1 5 6 4
5
2
2 3
4
1 5 6 3
2
1 4
2
4 5 3
1
2
4 6 1
2
5
2
4
3 ( )
( ) [ ( )
] ( ) ( )
( )
b b= , B.7203 103 ( )( ) ( )
b a a a a a a a a a a a a a= - + + + - - + , B.8213 2 42 1 5 32 4 1 6 2 5 3 22 6( ) ( ) ( )( )
b a a a a a a a a a a a a a a= - + + + - - , B.9223 1 52 2 4 32 5 1 4 6 2 3 6 3 42( ) ( )( )
b a a a a a a a a a a a= - + + - +2 . B.10233 2 4 6 2 52 32 6 3 4 5 43 ( )( )
In dCS gravity, we derived PM relations between the dimensionless dipole scalar charge
and the compactness, given as
⎜ ⎟
⎡
⎣⎢
⎛
⎝
⎞
⎠
⎤
⎦⎥ååm
l
p= + +
m
= =
c C C
64
5
1
7
, B.11
i j
i
ij
j
i
1
4
1
BL 5¯ ( ) ( )( ¯ )
where coefficients mcij
( ¯ ) are given in table B4. Equation (B.11) has the form of equation (B.1)
to x5( ). As we did for equation (B.2), one can then construct the (2, 2)-Padé approximant
⎡
⎣
⎢⎢
⎤
⎦
⎥⎥
å
åa
b
b
= +=
=
y x
x
x
x , B.12k i
i
i
i i
i
0
0
2
1
2
0
3
2
2
5( ) ( )
( )
( )
by identifying the constants
b a a a= - , B.13102 1 3 22 ( )( )
b a a a a a a a= + - - + , B.14112 12 3 22 4 1 2 3( ) ( )( )
Table B2. Coefficients lcij 2
( ¯ ) in equation (48) for l2¯ as a function of the stellar compactness C within the PM approximation.
lc1,02
( ¯ ) lc1,12
( ¯ ) lc1,22
( ¯ ) lc2,02
( ¯ ) lc2,12
( ¯ ) lc2,22
( ¯ )
-160
7
- 995
42
115
84
80480
441
461660
1323
82955
588
lc2,32
( ¯ ) lc2,42
( ¯ ) lc3,02
( ¯ ) lc3,12
( ¯ ) lc3,22
( ¯ ) lc3,32
( ¯ )
- 280
9
- 2825
3024
- 6598400
11319
-143668582
101871
- 2135837
2058
- 4341409
58212
lc3,42
( ¯ ) lc3,52
( ¯ ) lc3,62
( ¯ ) lc4,02
( ¯ ) lc4,12
( ¯ ) lc4,22
( ¯ )
58914551
465696
- 10915
2464
- 322075
266112
1713843200
3090087
6610428752
9270261
1547606966
27810783
lc4,32
( ¯ ) lc4,42
( ¯ ) lc4,52
( ¯ ) lc4,62
( ¯ ) lc4,72
( ¯ ) lc4,82
( ¯ )
- 351398725
7945938
234321937
1629936
2331123479
127135008
- 1030100675
18162144
- 5990225
314496
- 1097875
628992
lc5,02
( ¯ ) lc5,12
( ¯ ) lc5,22
( ¯ ) lc5,32
( ¯ ) lc5,42
( ¯ ) lc5,52
( ¯ )
8860516352
64891827
2217813273272
973377405
1550789734838
584026443
16971435941
30900870
21489165403
31783752
10043710866671
13349175840
lc5,62
( ¯ ) lc5,72
( ¯ ) lc5,82
( ¯ ) lc5,92
( ¯ ) lc5,102
( ¯ ) lc6,02
( ¯ )
- 1749233951
3302208
- 200852202065
254270016
- 130013401675
435891456
- 13797972625
290594304
- 113935625
41513472
3018408755200
9438155727
lc6,12
( ¯ ) lc6,22
( ¯ ) lc6,32
( ¯ ) lc6,42
( ¯ ) lc6,52
( ¯ ) lc6,62
( ¯ )
3742297951628752
1274151023145
55197038324995828
3822453069435
472845689999591
26003082105
24915470568822391
2184258896820
874072035313769
107864636880
- 958041289000963
499259176416
lc6,72
( ¯ ) lc6,82
( ¯ ) lc6,92
( ¯ ) lc6,102
( ¯ ) lc6,112
( ¯ ) lc6,122
( ¯ )
- 88927876896912359
6656789018880
- 2451739592866309
207484333056
- 1661453222619755
351127332864
- 698286437561675
702254665728
- 763587990875
7185604608
- 568540625
124084224
Table B3. Coefficients lcij
I 2( ¯ ¯ ) in equation (49) for I¯ as a function of the stellar compactness l2¯ within the PM approximation.
lc I1,0 2
( ¯ ¯ ) lc I1,1 2
( ¯ ¯ ) lc I1,2 2
( ¯ ¯ ) lc I2,0 2
( ¯ ¯ ) lc I2,1 2
( ¯ ¯ )
88
7
40
3
- 13
12
139616
2205
196360
1323
lc I2,2 2
( ¯ ¯ ) lc I2,3 2
( ¯ ¯ ) lc I2,4 2
( ¯ ¯ ) lc I3,0 2
( ¯ ¯ ) lc I3,1 2
( ¯ ¯ )
19127
252
- 11219
1176
3175
10584
5109760
33957
98970884
169785
lc I3,2 2
( ¯ ¯ ) lc I3,3 2
( ¯ ¯ ) lc I3,4 2
( ¯ ¯ ) lc I3,5 2
( ¯ ¯ ) lc I3,6 2
( ¯ ¯ )
1148580899
1528065
1850685733
6112260
- 30603863
2037420
- 4293193
1222452
- 10717951
9779616
lc I4,0 2
( ¯ ¯ ) lc I4,1 2
( ¯ ¯ ) lc I4,2 2
( ¯ ¯ ) lc I4,3 2
( ¯ ¯ ) lc I4,4 2
( ¯ ¯ )
2719077376
21068775
448894288
567567
73208568884
37923795
129979163705
55621566
70155024767
57047760
lc I4,5 2
( ¯ ¯ ) lc I4,6 2
( ¯ ¯ ) lc I4,7 2
( ¯ ¯ ) lc I4,8 2
( ¯ ¯ ) lc I5,0 2
( ¯ ¯ )
1602346982741
13349175840
- 1285405448827
17798901120
- 13875033761
889945056
- 15389455007
42717362688
- 507262148608
4866887025
lc I5,1 2
( ¯ ¯ ) lc I5,2 2
( ¯ ¯ ) lc I5,3 2
( ¯ ¯ ) lc I5,4 2
( ¯ ¯ ) lc I5,5 2
( ¯ ¯ )
- 7677705691856
14600661075
- 2244933136724
3369383325
1619202392194
1706570775
212537186204347
52562379870
24258802978997
5990014800
lc I5,6 2
( ¯ ¯ ) lc I5,7 2
( ¯ ¯ ) lc I5,8 2
( ¯ ¯ ) lc I5,9 2
( ¯ ¯ ) lc I5,102
( ¯ ¯ )
1882420967949367
1681996155840
- 28481247349223
120142582560
- 39174857283707
320380220160
- 1912091255831
96114066048
- 2754748757447
1922281320960
lc I6,0 2
( ¯ ¯ ) lc I6,1 2
( ¯ ¯ ) lc I6,2 2
( ¯ ¯ ) lc I6,3 2
( ¯ ¯ ) lc I6,4 2
( ¯ ¯ )
-120917752852250624
286683980207625
- 214334584423866976
57336796041525
- 2400708364625585884
172010388124575
- 43101354764736391546
1548093493121175
- 22084209085559827543
688041552498300
lc I6,5 2
( ¯ ¯ ) lc I6,6 2
( ¯ ¯ ) lc I6,7 2
( ¯ ¯ ) lc I6,8 2
( ¯ ¯ ) lc I6,9 2
( ¯ ¯ )
- 249310930642237571
19112265347175
333660111817466369461
33025994519918400
11506518464513510959
1000787712724800
8299909121934464119
2795851387929600
- 2919998240988289177
9607562042158080
lc I6,102
( ¯ ¯ ) lc I6,112
( ¯ ¯ ) lc I6,122
( ¯ ¯ )
- 8676262683332995337
35227727487912960
- 1036509722405634683
28182181990330368
- 8375090116034335069
5072792758259466240
b a a a a a a a a= - + + - -2 , B.15122 23 1 3 4 2 12 4 32( ) ( )( )
b b= , B.16202 102 ( )( ) ( )
b a a a a= - + , B.17212 1 4 2 3 ( )( )
b a a a= - . B.18222 2 4 32 ( )( )
References
[1] Steiner A W, Lattimer J M and Brown E F 2010 Astrophys. J. 722 33–54
[2] Lattimer J M and Steiner A W 2014 Astrophys. J. 784 123
[3] Miller M C 2013 arXiv:1312.0029
[4] Ozel F, Psaltis D, Guver T, Baym G, Heinke C and Guillot S 2015 arXiv:1505.05155
[5] Lattimer J M and Prakash M 2001 Astrophys. J. 550 426
[6] Lattimer J M and Prakash M 2007 Phys. Rept. 442 109–65
[7] Ozel F 2013 Rept. Prog. Phys. 76 016901
[8] Lattimer J M 2012 Ann. Rev. Nucl. Part. Sci. 62 485–515
[9] Alford M, Braby M, Paris M and Reddy S 2005 Astrophys. J. 629 969–78
[10] Alford M G, Burgio G F, Han S, Taranto G and Zappala D 2015 Phys. Rev. D 92 083002
[11] Prakash M, Cooke J and Lattimer J 1995 Phys. Rev. D 52 661–5
[12] Gendreau K C et al 2012 The Neutron star Interior Composition ExploreR (NICER): an Explorer
mission of opportunity for soft x-ray timing spectroscopy Society of Photo-Optical
Instrumentation Engineers (SPIE) Conf. Series vol 8443
[13] Psaltis D, Ozel F and Chakrabarty D 2014 Astrophys. J. 787 136
[14] Yagi K and Yunes N 2013 Science 341 365
[15] Yagi K and Yunes N 2013 Phys. Rev. D 88 023009
[16] Baubock M, Berti E, Psaltis D and Ozel F 2013 Astrophys. J. 777 68
[17] Psaltis D and Ozel F 2014 Astrophys. J. 792 87
[18] Pappas G and Apostolatos T A 2014 Phys. Rev. Lett. 112 121101
[19] Stein L C, Yagi K and Yunes N 2014 Astrophys. J. 788 15
[20] Yagi K, Kyutoku K, Pappas G, Yunes N and Apostolatos T A 2014 Phys. Rev. D 89 124013
[21] Chatziioannou K, Yagi K and Yunes N 2014 Phys. Rev. D 90 064030
[22] Majumder B, Yagi K and Yunes N 2015 Phys. Rev. D 92 024020
[23] Robinson D 1975 Phys. Rev. Lett. 34 905–6
[24] Israel W 1967 Phys. Rev. 164 1776–9
[25] Israel W 1968 Commun. Math. Phys. 8 245–60
[26] Hawking S 1971 Phys. Rev. Lett. 26 1344–6
[27] Hawking S W 1972 Commun. Math. Phys. 25 152–66
[28] Carter B 1971 Phys. Rev. Lett. 26 331–3
[29] Gurlebeck N 2015 Phys. Rev. Lett. 114 151102
[30] Geroch R P 1970 J. Math. Phys 11 2580–8
[31] Hansen R O 1974 J. Math. Phys 15 46–52
Table B4. Coefficients mcij
( ¯ ) in equation (B.11) for m¯ as a function of the stellar com-
pactness C within the PM approximation.
c I1,0
( ¯) c I1,1
( ¯) c I2,0
( ¯) c I2,1
( ¯) c I2,2
( ¯) c I3,0
( ¯) c I3,1
( ¯)
- 3
4
3
4
- 2
35
34
105
19
12
- 76
2205
1481
2205
c I3,2
( ¯) c I3,3
( ¯) c I4,0
( ¯) c I4,1
( ¯) c I4,2
( ¯) c I4,3
( ¯) c I4,4
( ¯)
186031
55440
1873
528
- 1451
121275
115371827
88288200
38537329
5045040
40345817
2882880
686755
82368
[32] Pani P 2015 Phys. Rev. D 92 124030
[33] Herrera L and Santos N O 1997 Phys. Rept. 286 53–130
[34] Bowers R L and Liang E P T 1974 Astrophys. J. 188 657
[35] Yagi K and Yunes N 2015 Phys. Rev. D 91 103003
[36] Yagi K and Yunes N 2015 Phys. Rev. D 91 123008
[37] Chandrasekhar S 1969 Ellipsoidal Figures of Equilibrium (New Haven: Yale University Press)
[38] Shapiro S L and Teukolsky S A 1983 Black Holes, White Dwarfs, and Neutron Stars: The Physics
of Compact Objects (New York: Wiley-Interscience)
[39] Poisson E and Will C M 2014 Gravity (Cambridge: Cambridge University Press)
[40] Chan T, Chan A P and Leung P 2015 Phys. Rev. D 91 044017
[41] Damour T and Nagar A 2009 Phys. Rev. D 80 084035
[42] Jackiw R 2003 and Pi S Y Phys. Rev. D 68 104012
[43] Smith T L, Erickcek A L, Caldwell R R and Kamionkowski M 2008 Phys. Rev. D 77 024015
[44] Alexander S and Yunes N 2009 Phys. Rept. 480 1–55
[45] Polchinski J 1998 String Theory Volume 2: Superstring Theory and Beyond (Cambridge:
Cambridge University Press)
[46] Alexander S H S and Gates S J Jr 2006 J. Cosmol. Astropart. Phys. JCAP06(2006)018
[47] Taveras V and Yunes N 2008 Phys. Rev. D 78 064070
[48] Calcagni G and Mercuri S 2009 Phys. Rev. D 79 084004
[49] Weinberg S 2008 Phys. Rev. D 77 123541
[50] Delsate T, Hilditch D and Witek H 2015 Phys. Rev. D 91 024027
[51] Horvat D, Ilijic S and Marunovic A 2011 Class. Quantum Grav. 28 025009
[52] Glampedakis K, Kapadia S J and Kennefick D 2014 Phys. Rev. D 89 024007
[53] Quintana H 1976 Astrosphys. J. 207 279–88
[54] Kojima Y and Hosonuma M 1999 Astrophys. J. 520 788–96
[55] Benhar O, Ferrari V, Gualtieri L and Marassi S 2005 Phys. Rev. D 72 044028
[56] Hartle J B 1967 Astrophys. J. 150 1005–29
[57] Hartle J B and Thorne K S 1968 Astrophys. J. 153 807
[58] Doneva D D and Yazadjiev S S 2012 Phys. Rev. D 85 124023
[59] Silva H O, Macedo C F B, Berti E and Crispino L C B 2015 Class. Quantum Grav. 32 145008
[60] Binnington T and Poisson E 2009 Phys. Rev. D 80 084018
[61] Hinderer T 2008 Astrophys. J. 677 1216–20
[62] Thorne K S, Price R H and MacDonald D A 1986 Black Holes: The Membrane Paradigm (New
Haven, CT: Yale University Press)
[63] Kol B and Smolkin M 2012 J. High Energy Phys. JHEP02(2012)010
[64] Chakrabarti S, Delsate T and Steinhoff J 2013 arXiv:1304.2228
[65] Quinzacara C and Salgado P 2012 Phys. Rev. D 85 124026
[66] Yagi K, Stein L C, Yunes N and Tanaka T 2013 Phys. Rev. D 87 084058
[67] Stergioulas N and Friedman J L 1995 Astrophys. J. 444 306
[68] Yagi K 2014 Phys. Rev. D 89 043011
[69] Damour T and Esposito-Farese G 1992 Class. Quantum Grav. 9 2093–176
[70] Damour T and Esposito-Farese G 1993 Phys. Rev. Lett. 70 2220–3
[71] Yeomans J M 1992 Statistical Mechanics of Phase Transitions (Oxford: Clarendon)
[72] de Boer J, Papadodimas K and Verlinde E 2010 J. High Energy Phys. JHEP10(2010)020
[73] Arsiwalla X, de Boer J, Papadodimas K and Verlinde E 2011 J. High Energy Phys. JHEP01
(2011)144
[74] Choptuik M W 1993 Phys. Rev. Lett. 70 9–12
[75] Gundlach C and Martin-Garcia J M 2007 Living Rev. Rel. 10 5
[76] GRTensorII this is a package which runs within Maple but distinct from packages distributed with
Maple. It is distributed freely on the World-Wide-Web from the address: http://
grtensor.org