Self-interacting dark matter cusps around massive black holes

0.49 MB·

by  Stuart L. Shapiro

Self-interacting dark matter cusps around massive black holes
Stuart L. Shapiro1,∗ and Vasileios Paschalidis1
1Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA
(Dated: February 4, 2014)
We adopt the conduction fluid approximation to model the steady-state distribution of matter
around a massive black hole at the center of a weakly collisional cluster of particles. By “weakly
collisional”wemeanaclusterinwhichthemeanfreetimebetweenparticlecollisionsismuchlonger
thanthecharacteristicparticlecrossing(dynamical)timescale,butshorterthantheclusterlifetime.
Whenappliedtoastarcluster,wereproducethefamiliarBahcall-Wolfpower-lawcuspsolutionfor
thestarsboundtotheblackhole. Herethestardensityscaleswithradiusasr−7/4 andthevelocity
dispersion as r−1/2 throughout most of the gravitational well of the black hole. When applied to a
relaxed, self-interacting dark matter (SIDM) halo with a velocity-dependent cross section σ∼v−a,
4 thegasagainformsapower-lawcusp,butnowtheSIDMdensityscalesasr−β,whereβ =(a+3)/4,
1
whileitsvelocitydispersionagainvariesasr−1/2. ResultsareobtainedfirstinNewtoniantheoryand
0
then in full general relativity. Although the conduction fluid model is a simplification, it provides
2
a reasonable first approximation to the matter profiles and is much easier to implement than a

n full Fokker-Planck treatment or an N-body simulation of the Boltzmann equation with collisional
a perturbations.
J

1 PACSnumbers: 95.35.+d,98.62.Js,98.62.-g
3

] I. INTRODUCTION core collapse (i.e. reversing the “gravothermal catastro-
O
phe”) was pointed out by Shapiro [9] and studied nu-
C merically via time-dependent Monte Carlo simulations
. Determining the stellar distribution around a massive of the two-dimensional Fokker-Planck equation to deter-
h
black hole at the center of a virialized star cluster is an mine f(E,J;t) [10]. (For a review of some of this early
p
- interesting and well-studied problem. Originally formu- work and references see [11]). The same problem has
o lated by Peebles [1, 2], the problem was solved for the received considerable attention in recent years by many
tr profile in a spherical globular cluster by Bahcall and authors,focusingonsuchaspectsasnonsphericalclusters
s Wolf [3]. They assumed an isotropic velocity distribu- [e.g.,[12],masssegregation[13],resonantrelaxation[14],
a
[ tion function to solve numerically the one-dimensional relativistic corrections, tidal disruptions, and extreme-
Fokker-Planck equation for the steady-state stellar dis- mass ratio inspirals (EMRIs) (see [15] for a review and
1
tribution function f(E) of stars bound to the black hole. references)].
v
HereE istheenergyperunitmassofastar. Theyfound
5 Inthispaperwereturntothethesteady-statedistribu-
that the cumulative, distant, two-body encounters (i.e.
0 tion of matter around a massive black hole at the center
0 small-angleCoulombscattering)betweenstarsdrivesthe of a virialized, weakly collisional, spherical gas, such as
0 density throughout most of the gravitational well of the
a star cluster. Here we show that the conduction fluid
. blackhole(the“cusp”)toapower-lawprofilethatvaries
2 model provides a straightforward means of deriving the
with radius as ρ ∝ r−7/4, while the velocity dispersion
0 steady-state matter density and velocity dispersion pro-
4 varies as v ∝ r−1/2. The effects of velocity anisotropy filesinthecluster. Itisapproximateinthatittreatsthe
1 and the role of the stellar disruption loss cone in the
lowest order moments of the distribution function, and
: cusp were delineated by Frank and Rees [4] and Light-
v not the distribution function itself, and adopts several
i man and Shapiro [5]. They showed that while the distri- physically reasonable, but simplifying, relations to close
X bution function depends only logarithmically on J, the
the moment equations. However, the approach is much
r angular momentum per unit mass of a star, the inward
a easier to implement than a full Fokker-Planck treatment
flux of stars and stellar disruption rate by the black hole
oranN-bodysimulationoftheBoltzmannequationwith
are affected more significantly. They also applied the re-
collisional perturbations. We then use this approach to
sultstocuspsindensegalacticnucleiaswellasstarclus-
determine the matter distribution in a cusp formed well
ters. Allofthisworkmotivateddetailedtwo-dimensional
inside the homogeneous core of a self-interacting dark
Fokker-Planck treatments to determine f(E,J), both by
matter (SIDM) halo containing a massive, central black
Monte Carlo [6, 7] and by finite-difference [8] methods.
hole.
The role of the removal of bound stars by the black hole
SIDM with cross sections per unit mass in the
in heating the ambient star cluster and reversing secular
range σ/m ≈ 0.45 − 450 cm2 gm−1 [or 8 ×
10−(25−22) cm2 GeV−1] was proposed by Spergel and
Steinhardt [16] to rectify the problem that cosmologi-
∗AlsoDepartmentofAstronomyandNCSA,UniversityofIllinois cal simulations (e.g., [17–19]) of purely collisionless cold
atUrbana-Champaign,Urbana,Illinois61801,USA dark matter (CDM) exhibit dark matter halos that are
2
highly concentrated and have density cusps, in contrast central black holes. This possibility motivates the study
to the constant density cores observed for galaxies and in this paper, for which we adopt a velocity-dependent
galaxy clusters. SIDM cross section σ ∝ v−a and consider a relaxed,
TheSIDMremedyfelloutoffavorforvariousreasons, spherical,SIDMcorewithamassiveblackholeatitscen-
includingpossiblealternativeexplanationsforthisappar- ter and determine the cusp profile well inside the core.
ent density discrepancy (e.g. feedback and expulsion of We note that a cusp can also form in a purely collision-
baryonicgas),Bulletclusterobservationsthatseemedto less CDM system containing a central black hole. For
constrain SIDM models to have a cross section too small example, if a central black hole grows adiabatically (e.g.
to have a significant effect on the structure of dark mat- bygasaccretion)itwillperturbtheCDMparticleorbits
ter halos, gravitational lensing and x-ray data suggest- whileholdingtheadiabaticinvariantsofthemotioncon-
ing that the cores of galaxy clusters are denser and less stant, and this process will result in a cusp [1, 41, 42].
spherical than predicted for SIDM, as well as theoretical ForaSIDMsystem,however,particlecollisionswillwash
biases. However, the interpretation of some of the data out any initial and/or adiabatically altered particle dis-
isnotdefinitive(see[20]foradiscussionandreferences), tribution in a collisional relaxation time scale and the
while new observational discrepancies with dissipation- cusp will relax to the solution determined below.
less CDM models have become apparent, such as the In Sec. II we summarize our basic model and the as-
absence of dwarf spheroidal galaxies (dSphs) predicted sumptions that define it. In Sec. III we present our
by CDM simulations (the “missing satellite problem”). formulation of the problem in Newtonian theory and use
Moreover, recent SIDM simulations [20, 21] now suggest simple scaling to derive the power-law cusp solution for
that SIDM systems with a (velocity-independent) cross bound particles obeying a power-law velocity-dependent
section as large as σ/m ≈ 0.1 cm2 gm−1 can yield the interactioncrosssection(SIDMmatterorstars). Wealso
coresizesandconstantcentraldensitiesobservedindark calculate the energy flux conducted outward by the par-
matter halos at all scales, from clusters, to low surface ticlesboundtotheblackholeinthecuspanddiscusshow
brightness spirals (LSBs) and dSphs, and are consistent thisfluxmayeventuallyhaltandreversesecularcorecol-
with revised Bullet cluster observational constraints. All lapse. In Sec. IV we reformulate the problem in general
of the other large-scale triumphs of CDM appear to be relativity and integrate the resulting equations numeri-
matched by these SIDM simulations. In addition, simu- cally to obtain the full density and velocity profiles in
lationsinvolvingvelocity-dependentSIDMcrosssections the cusp. In Sec. V we discuss our findings, evaluate
also appear to be successful in forming small, less con- the contribution of unbound particles, and consider the
centrated cores and altering the inner density structure very different behavior expected at very late times when
ofsubhalosinawaythatiscompatiblewithobservations a SIDM core evolves to a collisional fluid. We also sum-
of dSphs but having a negligible effect on galaxy cluster marize some future work.
scales[22–26]. Suchvelocity-dependentinteractionsarise We adopt geometrized units and set G=1=c below.
in“hidden”sectorextensionstotheStandardModelthat
have been constructed to explain some charged-particle
cosmic ray observations in terms of dark matter annihi- II. BASIC MODEL AND ASSUMPTIONS
lations [27–30].
The renewed viability of SIDM models motivates our
Our assumptions apply both to stars in a star cluster
own interest in returning to a subject that we stud-
and to dark matter in a SIDM halo; in either case we
ied earlier. We previously explored the nature of the
refertothematteras“particles.” Weadoptassumptions
“gravothermalcatastrophe”(secularcorecollapse)iniso-
similartothosemadebyBahcallandWolf[3]whotreated
lated SIDM halos [31] using the conduction fluid formal-
stars in a globular star cluster containing a central black
ism(seealso[32]andreferencestherein)andthenshowed
hole:
thattheformationofablackholeinthecenterofagalac-
tic halo may be a natural and inevitable consequence of
1. Thedistributionofparticlesissphericallysymmet-
gravothermal evolution [33].
ricinspace,isotropicinvelocityandhasrelaxedto
Supermassiveblackholes(SMBHs)withmassesinthe
a near-equilibrium state.
range 106−1010M are believed to be the engines that
(cid:12)
power active galactic nuclei and quasars. There is also 2. Outside the central cusp about the black hole the
substantialevidencethatSMBHsresideatthecentersof cluster has a nearly homogeneous central core in
many, and perhaps most, galaxies [34, 35], including the which the particle rest-mass density ρ and (one-
0
Milky Way [36–38]. While SMBHs may arise from bary- dimensional) velocity dispersion v are nearly con-
0
onic processes, such as the collapse of Pop III stars [39] stant.
or supermassive stars [40], followed by mergers and gas
accretion, their origin is not yet known. 3. ThecentralSchwarzschildblackholehasamassM
Itisthereforenotunreasonabletoimaginethatbound that is much less than the mass of the cluster core
structures containing dark matter of all sizes (galaxy butdominatesthetotalmassofallparticlesbound
clusters, galaxies, satellitesystems, etc)containmassive, to it in the cusp.
3
4. The particles all have the same mass m, which is In a SIDM halo relaxation is driven by close, large-
very small compared to M. angle, elastic interactions between particles. The relax-
ationtimescaleisthemeantimebetweensinglecollisions
5. Themeanfreepathoftheparticleswithrespectto
and is given by
self-interactions is much longer than their charac-
teristic radius from the black hole. 1
t (SIDM) =
r ηρvσ
6. The relaxation time scale due to self-interactions (cid:20)(cid:16) η (cid:17)(cid:18) ρ (cid:19)
betweenparticlesisshorterthantheageoftheclus- (cid:39) 0.8×109yr
2.26 10−24g cm−3
ter.
(cid:16) v (cid:17)(cid:18) v (cid:19)1−a(cid:18) σ (cid:19)(cid:35)−1
7. Aparticleisremovedfromthedistributionatasuf- × 0 0 (3)
107cm sec−1 v 1 cm2 g−1
ficiently small radius deep inside the gravitational 0
potentialwelloftheblackhole,butoutsideitshori-
zon. where σ =σ0(v/v0)−a is the cross section per unit mass
and the constant η is of order unity. For example, η =
(cid:112)
The characteristic gravitational capture radius of the 16/π ≈ 2.26 for particles interacting elastically like
black hole r is taken to be billiard balls (hard spheres) with a Maxwell-Boltzmann
h
velocity distribution [see [47], Eqs. (7.10.3), (12.2.8) and
M
(12.2.12)]. WenotethatforaCoulomb-likecrosssection,
r = , (1)
h v2 wherea=4,t (SIDM)scalesthesamewaywithv andρ
0 r
as t (stars), up to a slowly varying log factor: t ∝v3/ρ.
The cusp is the region r <∼ rh inside the core. Particles Werexploit this equivalence below. r
whose orbits lie entirely in the cusp are bound to the
AstarofradiusRistidallydisruptedbytheblackhole
black hole.
whenever it passes within a radius
By assumptions 5 and 6 the particle distribution func-
tion satisfies the collisionless Boltzmann equation in r (cid:39)R(M/m)1/3. (4)
D
steadystatetoahighdegree, buttheparticularsolution
towhichitrelaxesisdeterminedbytheperturbationsin-
InaSIDMhalo,aparticleplungesdirectlyintotheblack
ducedbycollisions. Bytakingsuitablevelocitymoments,
hole once it passes within a radius
thecollisionlessBoltzmannequationreducestofluidcon-
servation equations for an ideal gas if the velocity distri-
r =4M, (5)
mb
bution of the particles is isotropic. The conduction fluid
modelprovidesanapproximatemeanshandlingthissys- in Schwarzschild coordinates. The radius r is the ra-
mb
tem, closing the moments to third order when collisions dius of the marginally bound circular orbit with energy
arerelevant. ThisformalismwasadoptedbyLynden-Bell (includingrest-massenergy)E/m=1. Itisalsothemin-
and Eggleton [43] to treat the gravothermal catastrophe imum periastron of all parabolic (E/m=1) orbits. Par-
instarclusters(seealso[44])andbyBalbergetal.[31]to ticles that approach the black hole from large distances
analyze the gravothermal catastrophe in isolated SIDM aretypicallynonrelativistic(i.e. v (cid:28)c, E/m≈1)and
∞
halos (see also [45], where the validity of the conduction hence arrive on very nearly parabolic orbits. Any such
fluid approximation is discussed, and [32], where it was particle that penetrates within r must plunge directly
mb
shown to agree with N-body simulations that incorpo- into the black hole (see, e.g., [42, 48]).
rated collisional perturbations to treat SIDM systems).
The net result is that we may set an inner boundary
Inastarclusterrelaxationisdrivenbymultiple,small-
to the cusp at r =r in a star cluster and at r =r in
D mb
angle gravitational (Coulomb) encounters. The local re-
a SIDM halo. At this boundary the density of particles
laxation time scale is given by (see e.g., [44, 46])
plummets, so it is a good approximation to set the den-
sityequaltozeroatthisradius. Formainsequencestars
r (cid:29) M, so the velocities and gravitational fields that
33/2v3 D
t (stars) = , determinethestellardistributioninthecuspareentirely
r 15.4mρln(0.4N)
Newtonian. By contrast, the SIDM particle distribution
(cid:16) v (cid:17)3
(cid:39) 0.7×109yr extends into a region in which the spacetime is highly
km sec
(cid:18)M pc−3(cid:19)(cid:18)M (cid:19)(cid:18) 1 (cid:19) relativistic.
× (cid:12) (cid:12) , (2)
ρ m ln(0.4N)
where v(r) is the local one-dimensional velocity disper- III. NEWTONIAN MODEL
sion, ρ(r) = mn(r) the local (rest) mass density, n(r)
is the local stellar number density, and N is the total In this section we formulate the problem in Newto-
number of stars in the cluster. The three-dimensional nian physics, which will help guide our general relativis-
velocity dispersion is v (r)=31/2v(r), by isotropy. tic treatment in the next section. The basic conduction
m
4
fluid equations required to determine the secular evolu- and ρ(r). We introduce the following nondimensional
tionofboundparticlesdrivenbycollisionalrelaxationin variables:
the cusp are given by [31, 43, 44]
ρ˜=ρ/ρ , v˜=v/v , r˜=r/r , (11)
0 0 h
∂(ρv2) Mρ
=− (6)
∂r r2 Dropping the tildes (˜) for convenience, the nondimen-
sional ODEs become
∂L = −4πr2ρ(cid:26) ∂ (cid:12)(cid:12)(cid:12) 3v2 +p ∂ (cid:12)(cid:12)(cid:12) 1(cid:27) dv = D , (12)
∂r ∂t(cid:12) 2 ∂t(cid:12) ρ dr v2−aρ2r4
M M
dρ ρ 2D
= −4πr2ρv2(cid:18) ∂ (cid:12)(cid:12)(cid:12) (cid:19)ln(cid:18)v3(cid:19) dr = −v2r2 − v3−aρr4, (13)
∂t(cid:12) ρ
M
= 0. (7) wheretheconstantsmultiplyingLhavebeenlumpedto-
getherwithLtodefineanewnondimensionalluminosity
Equation(6)istheequationofhydrostaticequilibrium, constant D ∝ L [see Eq. (19)]. The coupled equations
wherethekineticpressureP satisfiesP =ρv2. Equation abovemustbesolvedtogetherwithtwoboundarycondi-
(7)istheenergyequation(thefirstlawofthermodynam- tions that guarantee that the cusp profiles join smoothly
ics) for the rate of change of the entropy s given by onto the ambient core at some radius r (cid:29) r . The
0 h
nondimensional boundary conditions are then
(cid:18)v3(cid:19)
s=ln . (8)
ρ b.c.’s: ρ=1=v, r =r (cid:29)1. (14)
0
The time derivatives in Eq. (7) are Lagrangian, but in Athirdboundaryconditionrequiresthedensitytovanish
steady state the cluster is virialized and at rest on a at the inner edge of the cusp at radius r = r or r ,
in D mb
dynamical time scale and the mean fluid velocity is ev- depending on whether the particles are stars or SIDM.
erywhere negligible. Hence the time derivatives satisfy Typically r (cid:28)r , which translates to
in h
d ≈ ∂ and can be set equal to zero in seeking the
dt ∂t
steady-state solution (cf. [3, 49]). As a result, the lu- b.c.: ρ=0, r =r (cid:28)1. (15)
in
minosity L due to heat conduction is constant. Hence
Eq. (7) can be replaced by To satisfy condition (15) the constant D must be chosen
appropriately, making D an eigenvalue of the system.
L=constant. (9) As a trivial example, consider the case D = 0 = L.
HerethesolutiontoEqs.(12)and(13),imposingbound-
By assumption (5), L may be evaluated as a conductive
ary conditions (14), becomes
heat flux in the long mean free path limit,
v/v = constant=1,
L 3 H2∂v2 0 (cid:18) (cid:19)
=− bρ . (10) M/r−M/r
4πr2 2 t ∂r ρ/ρ = exp 0 . (16)
r 0 v2
0
In writing Eq. (10) we evaluated the kinetic tempera-
ture of the particles according to k T =mv2, where k This solution represents an isothermal profile with zero
B B
heat flux. However, it does not satisfy boundary con-
is Boltzmann’s constant. The parameter b is constant
dition (15) and must therefore be rejected. Physically,
of order unity and H is the local particle scale height.
anyplausiblesolutionshouldexhibitanincreasingveloc-
For a gas of hard spheres with a Maxwell-Boltzmann
ity dispersion as one moves deeper into the gravitational
distribution the coefficient b can be calculated to good
well of the black hole, i.e. dv/dr < 0. According to
precision from transport theory, and has the value of
√
Eq. (12) this requires that we search for solutions with
b ≈ 25 π/32 ≈ 1.38 [see [50], Chap. 1, Eq. (7.6), and
D < 0. We must fine-tune the search to find that value
Problem 3]. The scale height for typical bound parti-
of D <0 that enables us to satisfy condition (15) at the
cles in the cusp can be estimated as H ≈ r. It will
desired radius r .
turn out that the density and velocity profiles will not in
depend on any of the constant numerical coefficients as-
sociatedwiththefactorsappearingintheequations(e.g.,
η,b,H/r, etc.). The relaxation time scales t are given A. Power-law solution
r
by Eq. (2) for stars and Eq. (3) for SIDM. We will use
Eq.(3)forbothcases,settinga=4intheself-interaction Wepostponedisplayingfullnumericalsolutionstothe
cross section σ ∝ v−a to treat star clusters which relax above system of equations for the cusp profiles until the
via Coulomb encounters. nextsection, wherewewillformulateandsolvethesame
Using Eqs. (3) and (10), Eqs. (6) and (9) yield two problem in general relativity. Here, however, we demon-
coupled, ordinary differential equations (ODEs) for v(r) strate that Eqs. (12) and (13) admit power-law solutions
5
thatshouldapplyinthecuspinterior,wellawayfromits This kinetic heat flux is transported out into the core
inner and outer boundaries. surrounding the cusp by particle scattering. This en-
Letv ∝r−α andρ∝r−β andsubstituteintoEqs.(12) ergy can have a significant impact on the evolution of
and(13). Equatingpowersofronbothsidesofthesetwo the cluster, eventually halting and ultimately reversing
equations yield two equations relating α and β, which secular core collapse (the gravothermal catastrophe) in
when solved simultaneously yield α =1/2 and β =(3+ an isolated system, as pointed out for star clusters [9]
a)/4, or and confirmed by detailed Monte Carlo simulations [10].
The asymptotic rate of reexpansion of the core can be
v ∼r−1/2, ρ∼r−(3+a)/4. (17) estimated by equating the heating rate emerging from
the cusp into the core as given by Eq. (21) to the rate
We note that for the case a = 4 applicable to star of increase of the core energy E˙ , where E ∼−M2/R ,
clusters we recover the Bahcall-Wolf scaling laws: v ∼ 0 0 0 0
M is the mass of the core, assumed constant, R is the
r−1/2, ρ∼r−7/4. Moreoverthenumericalintegrationof cor0e radius, v ∼ (M /R )1/2 and σ ∼ v−a. The0 result
the coupled ODEs with their boundary conditions yields 0 0 0 0 0
is
a profile that is in excellent agreement with the fitting
functions provided by Bahcall and Wolf to their numer- R →t2/(7−a) (22)
0
ical solutions of the Fokker-Planck equation in steady
state (see Sec. IVC). Equation (22) agrees with [9] for the reexpansion of a
The above power-law solution can also be obtained star cluster (a = 4) containing a massive central black
by applying the simple scaling argument of Shapiro and hole: R0 →t2/3.
Lightman [51], who argued that in steady state the en-
ergy flux of bound particles must be constant indepen-
dent of radius throughout the cusp and be transported IV. GENERAL RELATIVISTIC MODEL
on a relaxation time scale (cf. Eqs. (9) and (10)):
To treat the cusp we invoke Birkhoff’s theorem
N(r)E(r) to ignore the exterior core and halo and adopt the
L∼ = constant, (18)
t Schwarzschild metric to describe the spherical spacetime
r
in the cusp:
where N ∼ ρr3 ∼ r−β+3 is the number of particles be-
tween r and 2r, E(r) ∼ v2 ∼ r−1 is the characteristic dr2
ds2 =−(1−2M/r)dt2+ +r2dΩ2. (23)
energy of a bound particle orbiting in this layer, and 1−2M/r
t ∼ 1/(σρv) ∼ 1/r(a−1)/2−β. Inserting the factors in
r
Eq. (18) and solving for β yields the same result found In using the above (vacuum) metric we neglect the small
above, Eq. (17). Evaluating the result for a velocity- contribution to the stress-energy tensor of the particles
independent cross section with a = 0 gives ρ ∼ r−3/4, a orbiting in the cusp about the black hole M, in accord
result found previously [52] using the Shapiro-Lightman with assumption (3).
scaling argument. To determine the particle profiles in the cusp requires
the relativistic generalizations of Eqs. (6) and (7). Hy-
drostatic equilibrium becomes
B. Energy flux
dP dln|ξaξ | ρ+P M
=−(ρ+P) a =− , (24)
Once the eigenvalue D is determined numerically dr dr 1−2M/r r2
the outward kinetic energy flux conducted by particles
where ξa = ∂/∂t is the time Killing vector, so that
throughout the cusp can be evaluated explicitly. Restor-
|ξaξ | = |g |, P is the kinetic pressure of the particles
ing units we have a 00
and ρ is their total mass-energy density. The evolution
L=12πDηbσ ρ2M3v−3 =constant. (19) of the entropy per particle s is now governed by the first
0 0 0 law of thermodynamics together with energy conserva-
We may recast Eq. (19) as tion along fluid worldlines, ua∇bTab = 0, where Tab is
thetotalstressenergyofthesystem(fluidplusheat)and
N(r )E(r ) ua is the fluid four-velocity. As a result, Eqs. (7)–(10)
L=18Db h h , (20)
t (r ) now become
r h
ρ+P ds
where dρ/dτ − dn/dτ =nT =−∇ qa−a qa =0,
n dτ a a
N(r )E(r ) (4πr3ρ /3m))(mv2/2) (25)
h h = h 0 0 , (21)
t (r ) (ηρ v σ )−1 where τ is proper time, n is the proper particle num-
r h 0 0 0
ber density, T is the particle kinetic temperature, aa is
thereby justifying the scaling argument leading to theirfour-acceleration,andqa istheheatfluxfour-vector
Eq. (18). (see [53], Eq. 5.103). We will relate P and T to the
6
particle rest-mass density and velocity dispersion below. Intermsofthesevariables,butagaindroppingtildes(˜),
The last equality in Eq. (25) is imposed by our seeking a thetwonondimensionalODEsthatdeterminetheprofiles
steady-state solution. become
The only nonzero component of qa for a virialized gas
dv D fa−1(v ) v2
thatisatrestinastationary,sphericalgravitationalfield N = N − 0 (33)
is qr, which can be calculated from dr vN2−aρ2Nr4(1− 2vr02)3/2 2r2(1− 2vr02)1/2
κ d(cid:0)T|g00|1/2(cid:1) dρN = −ρNγ − 2D fa−1(vN) . (34)
qr = −|g |1/2 dr dr vN2 r2 v3−aρNr4(1− 2v02)3/2
00 r
(cid:18) (cid:19)
dT T M
= −κ + (26) Here the function f(v ) inverts Eq. (31) to give v in
dr 1−2M/r r2 N
terms of v ,
N
where κ is the thermal conductivity [see [54], Eq.
1
(22.16j)]. Similarly, the only nonzero component of aa v2 =v2 f2(v ), f2(v )= .
N N N (1+9v4v4 /4)1/2+3v2v2 /2
is given by 0 N 0 N
(35)
dln|ξaξ |1/2 1 M Appearing in the function f(v ) above, the velocity v
a = a = . (27) N N
r dr 1−2M/r r2 isagainnormalizedasinEq.(32),withv0 inunitsofthe
speedoflight. Equations(33)and(34)forthecuspmust
Using Eq. (27) in Eq. (25) and integrating gives be solved subject to the two outer boundary conditions
setbytheambientclustercore,whichinnondimensional
C
r2qr = , C =constant. (28) form become [cf. Eq. (14)]
(1−2M/r)1/2
b.c.’s: ρ =1=v , r =r (cid:29)1. (36)
Equation (28) is the relativistic analog of Eq. (9). N N 0
We determine the effective conductivity κ taking the
Inaddition,athirdboundaryconditionattheinneredge
Newtonian limit of Eq. (26) and equating the resulting
of the cusp must be satisfied [cf. Eq. (15)],
expression for qr to L/4πr2 given by Eqs. (10) and (3).
In this limit qr ≈ −κdT/dr, ρ ≈ mn, and T ≈ mv2/kB, b.c.: ρN =0, r =rin (cid:28)1, (37)
where k is Boltzmann’s constant. Matching yields
B
where, nondimensionally,
3
κ=Av1−an2r2, A= ηbmk σ va =constant. (29)
2 B 0 0 r = v2(R/M)(M/m)1/3 (stars),
in 0
= 4v2 (SIDM), (38)
CombiningEqs.(26)and(28)givesanequationforthe 0
temperature profile T(r):
[cf. Eqs. (4), and (5)].
dT C T M Once again, the third boundary condition is imposed
= − . (30) by finding the appropriate eigenvalue D <0.
dr κr2(1−2M/r)3/2 1−2M/r r2
Equations (24) and (30) provide the two coupled equa-
tions needed to determine the steady-state particle pro- A. GR vs Newtonian equations
files in the clusters. To apply them we need to relate the
quantities P(r) and T(r) to the density and velocity dis- FromtheaboveODEsitisevidentthatthenondimen-
persioninthesystem. Foroursimplifiedtreatmentbased sionalGRprofilesaredeterminedbyonefreephysicalpa-
on integrating the moments of the Boltzmann equation rameter, v , thecorevelocitydispersion. Thevalueofv
0 0
in the fluid conduction approximation it is adequate to sets the degree to which the ambient core is relativistic,
modeltheparticlesasaperfect,nearlycollisionless,rela- as well as the dynamic range between the outer bound-
tivisticgaswherealltheparticleshavethesamespeedlo- ary of the cusp, r and its inner boundary, r . In the
h in
callybutmoveisotropically. Ateachradiuswethenmay Newtonian limit, Eqs. (33) and (34) reduce to Eqs. (12)
set P ≡ nkBT = ρv2, where v is the one-dimensional and(13)andareindependentofthevalueofv0,although
velocity dispersion, ρ = γmn and γ = 1/(1−3v2)1/2, thelatterparameterenterstheinnerboundarycondition
which gives kBT =γmv2. through rin.
To cast the two ODEs in a form that most closely When v (cid:28) 1 the Newtonian equations suffice to de-
0
resembles the Newtonian Eqs. (12) and (13) we define termine the cusp solution for a cluster of normal stars,
the quantities as the stars never move at relativistic velocities or in a
strong-field region before being disrupted. For a SIDM
ρ ≡mn, v2 ≡k T/m=γv2, (31)
N N B cusp, by contrast, while the bulk of an SIDM cusp re-
sides in the Newtonian regime when v (cid:28) 1, the inner-
and the nondimensional variables 0
most particles have orbits which enter the high-velocity,
ρ˜ =ρ /ρ , v˜ =v /v , r˜=r/r , etc. (32) strong-field region outside the central black hole before
N N 0 N N 0 h
7
beingcaptured. TheGRequationsarethennecessaryto
obtain an accurate solution.
When v0 <∼ 1/3 the core is relativistic and the GR
equations must be used to treat the cusp everywhere,
both for a star or a SIDM cluster. We point out that as
v →1/3 in the core (i.e. three-dimensional velocity dis-
0
persion v →1) a (nearly) collisionless core-halo cluster
m
inevitably becomes dynamically unstable and undergoes
catastrophiccollapseonadynamical(free-fall)timescale,
formingamassive, centralblackhole withintheambient
halo (see, e.g., [53, 55] for reviews of and references to
bothanalytictheoryandsimulations). Subsequentgrav-
itational encounters (stars) or collisions (SIDM) will es-
tablishthecuspdescribedhereonarelaxationtimescale
t following the collapse.
r
B. Energy flux
The kinetic heat flux can be calculated from
L
4πr2 =|qaqa|1/2. (39) FtoIGth.e1:coPrreofidleenssoitfythρec,ufsoprrseeslte-cmtacshsodiecnessitoyfρtNhe,npoorwmearl-ilzaewd
0
a in the interaction cross section σ ∝ v−a. Starting from
The relativistic analog of Eq. (19) now becomes
the bottom and moving up, the solid curves show numerical
resultsfora=0,1,2,and4. Thedashedcurvesexhibitcrude
L= 12πDηbσ0ρ20M3v0−3 = constant . (40) analytic fits to the numerical profiles that apply to all radii
1−2M/r 1−2M/r r >∼ rin = 10−3rh [Eq. (43)], where rh is the cusp radius
[Eq. (1)].
Hence, as the radius from the black hole increases, the
outward energy flux flowing through the cusp and mea-
suredbyalocal,staticobserverdecreasesfasterthanr−2. for star clusters, and
The factor (1−2M/r) now appearing on the right-hand
side of Eq. (40) accounts for the gravitational redshift v =4.7×103r1/2 km/s (42)
0 in,3
(both in energy and time) as the heat flux propagates
outwards. for SIDM halos. Here we define M ≡ M/103M and
3 (cid:12)
r ≡103(r /r ).
in,3 in h
TheresultsoftheintegrationsareplottedinFig.1for
C. Numerical results therest-massdensityprofileρ inthecusp. Thenumer-
N
ical profiles confirm that the power-law density profiles
We integrate Eqs. (33) and (34), subject to bound- given by Eq. (17) apply to the bulk of the cusp, i.e. the
ary conditions Eqs. (36) and (37), for several different region well inside the inner and outer boundaries. The
choicesofthevelocitypowerlawscharacterizingthepar- profiles can be fit reasonably well to the general analytic
ticle interaction cross section, σ ∝ v−a. We take the expression
nondimensional outer boundary radius to be r = 10r ,
0 h
which puts it well outside the cusp and into the homo- ρN/ρ0 = 1+ξ(rh/r)(3+a)/4, r >∼rin, (43)
geneous core. We assign the inner boundary radius to = 0, r <∼rin,
be r = 10−3r for illustrative purposes. This radius is
in h
small enough to give the cusp sufficient dynamic range whereξ isoforderunity. WeplotfourcasesinFig.1cor-
to confirm the power-law radial dependence anticipated responding to velocity-dependent interaction cross sec-
by Eq. (17) and to establish scaling laws applicable to tions with power-law parameter a = 0,1,2, and 4. For
other choices of rin. For typical star clusters and SIDM these four cases shown in Fig. 1 we plot four analytic
systems, the typical values of rin are much smaller. In- curves for r >∼rin, setting (a,ξ)=(0,1.5),(1,1.2),(2,1),
verting Eq. (38), the core velocity v0 is related to rin and (4,1) in Eq. (43). These analytic curves are seen to
by match the numerical solutions reasonably well.
The corresponding velocity dispersion profiles for the
(cid:18) m (cid:19)1/6(cid:18)R (cid:19)1/2 fourcasesareplottedinFig.2. Theyalsoagreewiththe
v =140 (cid:12) M1/3r1/2 km/s (41)
0 M R 3 in,3 power-law (∼ Keplerian) profile predicted by Eq. (17).
(cid:12)
8
here (see Sec. II). The analytic fitting formulas defined
in Eqs. (43) and (44), with the very coefficients for these
formulasquotedabovefora=4,wereinfactconstructed
by Bahcall and Wolf to fit their own numerical solution
for star clusters. It is reassuring to discover that these
formulas also provide good fits to our numerical solution
forthiscase,basedasitisonamuchsimplerformulation
and calculation.
Interestingly, the same a = 4 case is also rele-
vant for some new models proposed for SIDM systems.
A velocity-dependent, non-power-law, elastic scattering
cross section, motivated by a Yukawa-like potential in-
volving a new gauge boson of mass m , (and similar to
φ
“screened” Coulomb scattering in a plasma) has been
proposed recently [56, 57]. N-body simulations adopt-
ing this interaction apparently can explain the observed
coresindwarfgalaxieswithoutaffectingthedynamicsof
larger systems with larger velocity dispersions, such as
clusters of galaxies [25, 26].
Nowthepresenceofamassive,centralblackholeinthe
core of a SIDM system will develop a cusp, as we have
discussed. Foranyinteractioncross-section,powerlawor
FIG.2: Profilesofthecuspvelocitydispersionv,normalized
not,thedensityandvelocityprofilesinsuchacuspwould
to the core velocity dispersion v , for select choices of the
0 be straightforward to calculate by the method formu-
power-law parameter a in the interaction cross section σ ∝
v−a. Startingfromthetopandmovingdown,thesolidcurves lated in this paper. However, we can already anticipate
show numerical results for a = 0,1,2, and 4. The dashed theresultsforsystemsinteractingviatheaboveYukawa-
curve is a crude analytic fit to the numerical profiles that like potential, provided their core velocity dispersions v
applies to all radii r >∼ rin = 10−3rh [Eq. (44)], where rh is satisfy βc ≡ π(vmax/v)2 <∼ 0.1. In this case the cross
the cusp radius [Eq. (1)]. sections vary with velocity according to σ ∼ β2 ∼ v−4,
c
up to a slowly varying logarithmic factor. Here v is
max
the velocity at which the momentum-transfer weighted
For comparison, the analytic fit scattering rate (cid:104)σv(cid:105) peaks at a cross section value of
σmax = 22.7/m2. When v is comparable to the ve-
4 φ max
v/v0 =1+ 11(rh/r)1/2, r >∼rin, (44) lgoacliatxyyd,itshpeerNsio-bnoidnythsiemcuolraetioofnas ttyhpaitcarelsduwltarafresprheeprooritdeadl
togiveprofilesinbetteragreementwithobservationsfor
is also shown on the plot.
theseobjectsthancollisionlessCDMpredictions. Adopt-
The solution of Eqs. (33) and (34) is very sensitive
ingsuchavalueforv impliesthatthroughoutthebulk
to the eigenvalue D. For each choice of a and r , max
in
of a cusp in dwarf galaxies, as well as in any larger sys-
high integration accuracy and iteration are neces-
tem, the interaction cross section resides in the n = 4
sary to find the value of D to generate a solution
power-law regime, for which the Bahcall-Wolf solution
that matches the required inner boundary condi-
applies in a first approximation.
tion (37). The four cases examined above yield (a,D)=
(0, 0.7118010817),(1, 0.58832494372),(2, 0.50056792714)
and (4, 0.3847292732806). The eigenvalue D determines
A. Unbound particles
the outgoing heat flux through Eq. (40).
In addition to the bound particles that orbit in the
cusp there are unbound particles (E/m > 1, where E
V. DISCUSSION
again includes rest-mass energy) from the ambient core
that penetrate into the cusp. Treating the core as an in-
The case a = 4 (σ ∝ v−4) is particularly interest-
finite,collisionlessandmonoenergeticbathofsuchparti-
ing, both for star clusters and SIDM halos. This value
cles with rest-mass density ρ moving isotropically with
corresponds to a Coulomb scattering cross section and 0
(one-dimensional) velocity dispersion v yields a density
therefore yields the familiar cusp profiles for a star clus- 0
profileinthegravitationalwelloftheblackholegivenby
ter containing a massive, central black hole: ρ ∝r−7/4
N
and v ∝r−1/2. Our numerical results are in close agree- ρ /ρ =(1+2r /3r)1/2, (45)
N 0 h
ment with the steady-state profile found by Bahcall and
and a velocity profile
Wolf [3], who integrated the Fokker-Planck equation for
a star cluster satisfying the same assumptions adopted v/v =(1+2r /3r)1/2, (46)
0 h
9
assumingNewtoniangravitation[48,58](forthesolution haloes and Fig. 7 for subhaloes). Inside a cusp around a
in general relativity, see [59]). The key point is that a central black hole the ratio becomes
comparison of Eqs. (17) or (43) with (45) shows that the
λ 1 λ (cid:18)M (cid:19)(cid:16)r (cid:17)(1+a)/4
dtheenscituyspofisunmbuocuhndsmpaalrlteirclethsapnentehteradteinngsittyhethinatterbiourildosf r = ρσr ≈ R00 M0 rh , r <∼rh, (50)
up in bound particles for all plausible (a > 0) velocity-
dependent interactions. where we have used Eq. (17), and thus increases as r
decreases below r for all a > −1. Thus a typical
h
SIDM cluster with a central black hole is characterized
by λ/r (cid:29) 1 everywhere and the matter resides in the
B. Collisional fluid regime
weakly, and not the fluid, collisional regime.
NowSIDMclusterscanevolvefromaweaklycollisional
As shown in [31], gravothermal evolution in a SIDM
gastoastronglycollisionalfluid,transformingthenature
cluster characterized by a velocity-independent self-
of the cusp profiles accordingly. However, SIDM models
interaction inevitably drives the contracting core to suf-
constructed to explain observed clusters are specifically
ficiently high density that the particle mean free path
designed not to have undergone such evolution as yet.
to collisions becomes smaller than the local scale height
In particular, these models employ interaction cross sec-
in the innermost regions. A similar situation presum-
tionswhosemagnitudesaresufficientlysmallthatsecular
ablyarisesforarangeofvelocity-dependentinteractions.
core collapse – the gravothermal catastrophe – will not
When this situation occurs, the particles in the inner-
have occurred in a Hubble time. In these cases there
most regions will behave as a collisional (hydrodynam-
is insufficient time for the cores to have evolved to a
ical) fluid. A massive black hole at the center of such
fluid state, which requires a significant fraction of a core
a core would then be expected to accrete particles via
collapse time scale (≈ 290t for isolated clusters with
r
steady-state, adiabatic, spherical Bondi [60] flow. The
a=0; [31]). This conclusion is evident from the relation
cusp in such a case will fill with gas that has a density
t /t ≈ λ /R , where t ≈ 1/ρ1/2 is the charac-
profile of the form r dyn 0 0 dyn 0
teristic dynamical time scale in the core. We thus have
t ≈5×107(λ /R )(0.1M pc−3/ρ )1/2yr,wherebysub-
ρ /ρ ≈1+χ(r /r)3/2, χ∼O(1), (47) r 0 0 (cid:12) 0
N 0 h stituting Eq. 49, yields a typical core collapse time scale
significantly longer than the Hubble time. Moreover,
and an inward radial velocity approaching the free-fall
corecollapseissuppressedin(nonisolated)clusterswhen
velocity v ∼ (M/r)1/2 inside r . The temperature, or
ff h merging from hierarchical formation is included [45].
velocity dispersion of the particles, γv2 ∼ k T/m, will
B
satisfy
T/T ≈(ρ /ρ )Γ−1. (48) C. Future work
0 N 0
HeretheadiabaticindexΓis5/3foranonrelativisticgas The possibility of a massive black hole at the centers
(k T/m(cid:28)1)and4/3forarelativisticgas(k T/m(cid:29)1) ofSIDMclustersandtheformationofacusparoundthe
B B
and varies with radius from the black hole (see [48] for black hole motivates several interesting questions that
a review of spherical Bondi flow, including a relativistic we are now investigating. One question is: assuming
treatment). For a nonrelativistic core temperature T , Γ the SIDM particles experience weak (WIMP-like) inelas-
0
willbenear5/3throughoutmostofthecusp, decreasing tic interactions in addition to the elastic interactions as-
below 5/3 as the fluid approaches the horizon. In this sumed here, what are the observable signatures of such
case Eqs. (47) and (48) show that the temperature of a cusp? Will, for example, the perturbative inelastic in-
the fluid at the horizon will be marginally relativistic, teractions produce detectable radiation or high-energy
climbing to k T /m∼1, independent of T . particles above the values expected from a homogeneous
B horz 0
Theratioofthecollisionmeanfreepathtothecharac- SIDM core without a cusp [61]? Another question is,
teristic scale height in the core of a typical SIDM model given that the cusp serves as a source of kinetic energy
constructed to explain currently observed clusters satis- conducted into the ambient core, is the energy sufficient
fies to reverse secular core collapse and cause the cluster to
reexpand, as discussed in Sec. IIIB [62]? These are just
λ 1
0 = some of the issues that we will pursue in future studies
R0 ρ0σ0R0 of SIDM clusters.
≈ 50(cid:18)0.1 M(cid:12)/pc3(cid:19)(cid:18)1 cm2/g(cid:19)(cid:18)1 kpc(cid:19)(,49) Acknowledgments: ItisapleasuretothankS.Balberg,
ρ σ R B. Fields, A. Peter and P. Shapiro for useful discussions.
0 0 0
This paper was supported in part by NSF Grants No.
and this ratio increases for r >∼ R0. The above inequal- PHY-0963136 and No. PHY-1300903 and NASA Grants
ity is consistent with, e.g., the velocity-dependent SIDM No. NNX11AE11G and No. NN13AH44G at the Uni-
models constructed in [26] (see their Fig. 4 for main versity of Illinois at Urbana-Champaign. V. P. grate-
10
fully acknowledges support from a Fortner Fellowship at UIUC.
[1] P. J. E. Peebles, General Relativity and Gravitation 3, N. Weiner, Phys. Rev. D 79, 015014 (Jan. 2009)
63 (Jun. 1972) [29] P.J.FoxandE.Poppitz,Phys.Rev.D 79,083528(Apr.
[2] P. J. E. Peebles, Astrophys. J. 178, 371 (1972) 2009)
[3] J. N. Bahcall and R. A. Wolf, Astrophys. J. 209, 214 [30] J. L. Feng, M. Kaplinghat, and H.-B. Yu, Physical Re-
(Oct. 1976) view Letters 104, 151301 (Apr. 2010)
[4] J. Frank and M. J. Rees, Mon. Not. R. Astro. Soc. 176, [31] S. Balberg, S. L. Shapiro, and S. Inagaki, Astrophys. J.
633 (Sep. 1976) 568, 475 (Apr. 2002)
[5] A. P. Lightman and S. L. Shapiro, Astrophys. J. 211, [32] J. Koda and P. R. Shapiro, Mon. Not. R. Astro. Soc.
244 (Jan. 1977) 415, 1125 (Aug. 2011)
[6] S. L. Shapiro and A. B. Marchant, Astrophys. J. 225, [33] S. Balberg and S. L. Shapiro, Physical Review Letters
603 (Oct. 1978) 88, 101301 (Mar. 2002)
[7] A. B. Marchant and S. L. Shapiro, Astrophys. J. 234, [34] D. Richstone, E. A. Ajhar, R. Bender, G. Bower,
317 (Nov. 1979) A.Dressler,S.M.Faber,A.V.Filippenko,K.Gebhardt,
[8] H. Cohn and R. M. Kulsrud, Astrophys. J. 226, 1087 R.Green,L.C.Ho,J.Kormendy,T.R.Lauer,J.Magor-
(Dec. 1978) rian,andS.Tremaine,Nature(London) 395,A14(Oct.
[9] S. L. Shapiro, Astrophys. J. 217, 281 (Oct. 1977) 1998)
[10] A. B. Marchant and S. L. Shapiro, Astrophys. J. 239, [35] L. Ho, in Observational Evidence for the Black Holes in
685 (Jul. 1980) the Universe, Astrophysics and Space Science Library,
[11] S.L.Shapiro,inDynamicsofStarClusters,IAUSympo- Vol. 234, edited by S. K. Chakrabarti (1999) p. 157
sium,Vol.113,editedbyJ.GoodmanandP.Hut(1985) [36] R. Genzel, A. Eckart, T. Ott, and F. Eisenhauer, Mon.
pp. 373–412 Not. R. Astro. Soc. 291, 219 (Oct. 1997)
[12] E. Vasiliev and D. Merritt, Astrophys. J. 774, 87 (Sep. [37] A.Eckart,R.Genzel,T.Ott,andR.Scho¨del,Mon.Not.
2013) R. Astro. Soc. 331, 917 (Apr. 2002)
[13] T. Alexander and C. Hopman, Astrophys. J. 697, 1861 [38] A. M. Ghez, S. Salim, S. D. Hornstein, A. Tanner, J. R.
(Jun. 2009) Lu, M. Morris, E. E. Becklin, and G. Duchˆene, Astro-
[14] K. P. Rauch and S. Tremaine, New Astronomy 1, 149 phys. J. 620, 744 (Feb. 2005)
(Oct. 1996) [39] P. Madau and M. J. Rees, Astrophys. J. Lett. 551, L27
[15] P. Amaro-Seoane, ArXiv e-prints(May 2012), (Apr. 2001)
arXiv:1205.5240 [astro-ph.CO] [40] M. Shibata and S. L. Shapiro, Astrophys. J. Lett. 572,
[16] D.N.SpergelandP.J.Steinhardt,PhysicalReviewLet- L39 (Jun. 2002)
ters 84, 3760 (Apr. 2000) [41] P.GondoloandJ.Silk,PhysicalReviewLetters83,1719
[17] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astro- (Aug. 1999)
phys. J. 490, 493 (Dec. 1997) [42] L. Sadeghian, F. Ferrer, and C. M. Will, Phys. Rev. D
[18] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, 88, 063522 (Sep. 2013)
A. V. Kravtsov, A. A. Klypin, J. R. Primack, and [43] D.Lynden-BellandP.P.Eggleton,Mon.Not.R.Astro.
A.Dekel,Mon.Not.R.Astro.Soc.321,559(Mar.2001) Soc. 191, 483 (May 1980)
[19] R. H. Wechsler, J. S. Bullock, J. R. Primack, A. V. [44] L.Spitzer,Jr.,DynamicalEvolutionofGlobularClusters
Kravtsov, and A. Dekel, Astrophys. J. 568, 52 (Mar. (Princeton, NJ, Princeton University Press, 1987)
2002) [45] K.AhnandP.R.Shapiro,Mon.Not.R.Astro.Soc.363,
[20] M.Rocha,A.H.G.Peter,J.S.Bullock,M.Kaplinghat, 1092 (Nov. 2005)
S. Garrison-Kimmel, J. On˜orbe, and L. A. Moustakas, [46] A. P. Lightman and S. L. Shapiro, Reviews of Modern
Mon. Not. R. Astro. Soc. 430, 81 (Mar. 2013) Physics 50, 437 (Apr. 1978)
[21] A.H.G.Peter,M.Rocha,J.S.Bullock,andM.Kapling- [47] F.Reif,FundamentalsofStatisticalandThermalPhysics
hat, Mon. Not. R. Astro. Soc. 430, 105 (Mar. 2013) (New York, McGraw-Hill, 1965)
[22] C.Firmani,E.D’Onghia,V.Avila-Reese,G.Chincarini, [48] S. L. Shapiro and S. A. Teukolsky, Black Holes, White
and X. Herna´ndez, Mon. Not. R. Astro. Soc. 315, L29 Dwarfs,andNeutronStars: ThePhysicsofCompactOb-
(Jul. 2000) jects (New York, Wiley, 1983)
[23] J.L.Feng,M.Kaplinghat,H.Tu,andH.-B.Yu,Journal [49] Thesameapproachisadoptedinconstructingaahydro-
of Cosmology & Astroparticle Physics 7, 4 (Jul. 2009) static equilibrium model of a slowly evolving star.
[24] M. R. Buckley and P. J. Fox, Phys. Rev. D 81, 083522 [50] E.M.LifshitzandL.P.Pitaevskii,Physicalkinetics (Ox-
(Apr. 2010) ford, Pergamon Press, 1981)
[25] A. Loeb and N. Weiner, Physical Review Letters 106, [51] S.L.ShapiroandA.P.Lightman,Nature(London) 262,
171302 (Apr. 2011) 743 (Aug. 1976)
[26] M. Vogelsberger, J. Zavala, and A. Loeb, Mon. Not. R. [52] A. Peter, unpublished project report (2004)
Astro. Soc. 423, 3740 (Jul. 2012) [53] T. W. Baumgarte and S. L. Shapiro, Numerical Relativ-
[27] M. Pospelov, A. Ritz, and M. Voloshin, Physics Letters ity: SolvingEinstein’sEquationsontheComputer (Cam-
B 662, 53 (Apr. 2008) bridge, Cambridge University Press, 2010)
[28] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and [54] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravi-