The merger of binary white dwarf-neutron stars: Simulations in full general relativity

1.3 MB·English

by  Vasileios Paschalidis

The merger of binary white dwarf–neutron stars:
Simulations in full general relativity
Vasileios Paschalidis1, Yuk Tung Liu1, Zachariah Etienne1, and Stuart L. Shapiro1,2
1 Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801
2 Department of Astronomy and NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801
We perform fully general relativistic (GR) simulations to address the inspiral and merger of
binary white dwarf-neutron stars. The initial binary is in a circular orbit at the Roche critical
separation. The goal is to determine the ultimate fate of such systems. We focus on binaries
whose total mass exceeds the maximum mass (M ) a cold, degenerate EOS can support against
max
gravitational collapse. The time and length scales span many orders of magnitude, making fully
2
generalrelativistichydrodynamic(GRHD)simulationscomputationallyprohibitive. Forthisreason,
1
0 we model the WD as a “pseudo-white dwarf” (pWD) as in our binary WDNS head-on collisions
2 study[1]. OurGRHDsimulations of a pWDNSsystem with a 0.98M⊙ WD and a 1.4M⊙ NS show
thatthemergerremnantisaspinningThorne-Zytkow-likeObject(TZlO)surroundedbyamassive
n
disk. ThefinaltotalrestmassexceedsM ,buttheremnantdoesnotcollapsepromptly. Toassess
max
a
whether the object will ultimately collapse after cooling, we introduce radiative thermal cooling.
J
Wefirstapplyourcoolingalgorithm toTZlOsformedinpWDNShead-oncollisions, andshowthat
3
these objects collapse and form black holes on the cooling time scale, as expected. However, when
1
we cool the spinning TZlO formed in the merger of a circular-orbit pWDNS binary, the remnant

doesnotcollapse,demonstratingthatdifferentialrotationalsupportissufficienttopreventcollapse.
]
E Given that the final total mass exceeds Mmax for our cold EOS, magnetic fields and/or viscosity
mayredistributeangularmomentum,ultimatelyleading todelayedcollapse toaBH.Weinferthat
H
the merger of realistic massive WDNS binaries likely will lead to the formation of spinning TZlOs
h. that undergo delayed collapse.
p
- PACSnumbers: 04.25.D-,04.25.dk,04.40.Dg
o
r
t
s I. INTRODUCTION binaries and their measured orbital properties. The
a NS masses in these systems range between 1.26 M⊙
[
During inspiral and merger, compact binaries emit a and 2.08 M⊙, and their distribution is centered around
2 large flux of gravitational waves (GWs), making them 1.5 M⊙. On the other hand, the WD masses in these
7v among the most promising sources for GWs detectable systems range between 0.125 M⊙ and 1.3 M⊙, and their
7 by ground-based laser interferometers like LIGO [2, 3], distribution is centered around 0.6 M⊙. Eighteen of
these observed WDNS binaries have total mass greater
1 VIRGO [4, 5], GEO [6], TAMA [7, 8] and AIGO [9], as
5 well as by proposed space-based interferometers such as than 1.65 M⊙, 8 of which have a WD component with
9. LISA [10] and DECIGO [11]. Extracting physical infor- massgreaterthan0.8M⊙,and5havetotalmassgreater
0 mation about these binaries from their GWs may shed than 2.2 M⊙. This is interesting because the expected
Tolman-Oppenheimer-Volkoff(TOV) limiting mass for a
1 light on determining their ultimate fate, but requires
1 carefulmodelingofthesesystemsinfullgeneralrelativity cold, degenerate gas must be larger than 1.97 M⊙ [35]
v: (see [12] for review and references therein). Most effort and may reach 2.2 M⊙ [36–44], depending on the equa-
tion of state (EOS). One of the main goals of this work
i to date has focused on modeling black hole–black hole
X is to determine whether a WDNS merger remnant will
(BHBH) binaries (see [13] and references therein), and
r neutron star–neutron star (NSNS) binaries (see [14] for undergo prompt collapse to a black hole.
a
a review), with some recent work on black hole–neutron
Population synthesis calculations [45, 46] show that
star (BHNS) binaries in full GR [15–33].
there are about 2.2 106 WDNS binaries in our Galaxy,
As a follow-up to our investigation of binary WDNS andthattheyhavea×mergerratebetween10−6 yr−1 and
head-on collisions [1], in this work we perform fully gen- 1.4 10−4 yr−1. Furthermore,thesestudies find thataf-
eral relativistic simulations of circular-orbit WDNS bi- ter×ayearofintegration,LISA-likeinterferometersshould
naries through inspiral and merger. Throughout we call be able to detect 1 40 WDNS pre-mergerbinaries. Re-
this the “inspiral case” to distinguish it from the “head- centworkbyThom−psonetal. [47]suggeststhatthelower
on”collisioncase. WDNS binariesarepromisingsources limit on the merger rate of binary WDNSs in the Milky
of low-frequency GWs (for LISA and DECIGO) and, as Way,at95%confidence,is2.5 10−5 yr−1. Thompsonet
we arguedin [34], possiblyalso high-frequencyGWs (for al. alsosuggestthatthemerge×rrateinthelocaluniverse
LIGO,VIRGO,GEO,TAMAandAIGO),iftheremnant is 0.5 1 104Gpc−3 yr−1. Therefore, ignoring some
ultimately collapses to a black hole. un∼certain−ties×,allrecentpopulationsynthesiscalculations
Like NSNS binaries, WDNS binaries are known to ex- suggest that LISA-like projects should be able to detect
ist. In [34] we compiled tables of 20 observed WDNS a few WDNS pre-mergers per year.
2
A. Previous WDNS work track.
In [1] we employed the Illinois adaptive mesh refine-
In [34] we focused on possible evolutionary scenarios ment (AMR) relativistic hydrodynamics code [24, 62] to
for circular WDNS binaries that have inspiraled suffi- performthefirstsimulationsofthesesystemsinfullGR.
ciently close that they reach the termination point for In particular, we studied the head-on collision from rest
equilibrium configurations. This is the Roche limit for at large separation of a massive WD and a NS. We fo-
WDNSs, at which point the WD fills its Roche lobe and cused on compact objects whose total mass exceeds the
may experience one of at least two possible fates: i) sta- maximum mass supportable by a cold EOS in order to
ble mass transfer (SMT) from the WD across the inner determine the outcome of such collisions.
LagrangepointontotheNS,orii)tidaldisruptionofthe The vast range of time and length scales involved in
WD by the NS via unstable mass transfer (UMT). the WDNS problem make fully general relativistic sim-
Note that once an UMT binary reaches the critical ulations extremely challenging. In [1] we demonstrated
Rocheseparation,furtherinspiralandmergerisgoverned that the length scales span four orders of magnitude, as
bytidaleffects andhydrodynamicalinteractionsandnot measured in neutron star radii, and that the associated
GW emission. time scales span six orders of magnitude in M, the to-
We also studied key parameters that determine talsystemmass. Currentnumericalrelativitytechniques
whether a system will undergo SMT or UMT and found andavailablecomputationalresourcesmakesuchsimula-
that, foragivenNSmass,there existsacriticalmassra- tions prohibitive. For this reason, we tackled this prob-
tioq 0.5thatseparatestheUMTandSMTregimes. lem using a different strategy.
crit
If the m≈ass ratio q = M /M of a WDNS system is Inparticular,weconstructedasix-parameterpiecewise
WD NS
such that q > q , the WD quickly overfills its Roche polytropic EOS which mimics realistic NS EOSs while,
crit
lobe,andthebinarywillultimatelyundergoUMT.Inthe at the same time, scales down the size of the WD. We
opposite case, q < q , the system will undergo SMT. callthesescaled-downWDs“pseudo-WDs(pWDs)”. We
crit
We showedthat a quasistationarytreatmentis adequate chose all of the piecewise EOSs in such a way that the
tofollowtheevolutionofanSMTbinaryduringthissecu- maximumNSmassis1.8M⊙[63],andthemaximumWD
larphaseandcalculatedthegravitationalwaveforms. We massis1.43M⊙, i.e.,the Chandrasekharmass. Further-
also pointed out that WDNS observations suggest that more, we made sure these EOSs preserve the qualitative
there are known candidates residing in both regimes. shape of the central density–mass curves as well as the
In the case of tidal disruption (UMT), by contrast, numberofstableandunstableNSandWDbranches(see
the system will evolve on a hydrodynamical (orbital) Figs. 1 and 2 in [1]). Moreover,the scaling is performed
time scale. In this scenario the NS may plunge into so that all the length-scale and time-scale inequalities of
the WD and spiral into the center of the star, forming a the realistic problem are left unchanged. For a given set
quasiequilibrium configuration that resembles a Thorne- of EOS parameters, a realistic WD has a counterpart
Zytkow object (TZO) [48]; alternatively, the NS may be pWD whichhas the same mass butis smallerinsize. As
the receptacle of massive debris fromthe disrupted WD. a result, for every realistic WDNS system, we can con-
WDNS mergers may also give rise to gamma ray bursts struct a pWDNS counterpart which involves the same
[49, 50]. (realistic) NS and the pWD counterpart of the WD.
The ultimate fate of the merged WDNS depends on UsingpWDs weperformeda sequenceofhead-onsim-
(1) the initial mass of the cold progenitor stars, (2) the ulations in which the EOS is changed so that the pWDs
degree of mass and angular momentum loss during the have the same mass (0.98 M⊙) but decreasing com-
WD disruption and binary merger phases, (3) the angu- pactions, while the compaction and mass of the NS in-
lar momentum profile of the WDNS remnant, and (4) volved remains practically unchanged. More precisely,
the extent to which disrupted debris is heated by shocks while keeping the masses of the binary components and
and/or nuclear reactions as it settles onto the NS and the NS radius fixed, the pWD compaction was modified
forms an extended, massive mantle. These are issues so thatthe pWD:NS radiusratio variedbetween 5:1and
that requirea hydrodynamicsimulationto resolve. Note 20:1. We then scaled the results of our simulations to
that Newtonian workon binarieswith a WD component predict the outcome in the realistic case: 500:1.
has been performed analytically in [34, 51–55] and via In addition to studying the effects of the pWD com-
Newtonian hydrodynamic simulations in [56–61]. How- paction, we also studied the effects of NS mass. We con-
ever, ascertaining whether or not the neutron star ulti- sidered NSs with masses 1.4 M⊙, 1.5 M⊙, and 1.6 M⊙.
matelyundergoescatastrophiccollapse(eitherpromptor All head-on collision simulations that we performed
delayed) to a black hole requires that such a simulation showedthatafterthecollision,14%-18%oftheinitialto-
be performed in full general relativity. In fact, even the talrestmassescapestoinfinity. Inallcases,theremnant
final fate of the NS in the alternative scenario in which restmassexceededthemaximumrestmassthatourcold
there is a long epoch of SMT may also lead to catas- EOS can support (1.92 M⊙), and no case led to prompt
trophic collapse, if the neutron star mass is close to the collapse to a black hole. This outcome arises because
neutron star maximum mass. This scenario too will re- the final configurations become hot, due to shock heat-
quire a general relativistic hydrodynamic simulation to ing. Allourcasessettleintoasphericalquasiequilibrium
3
configurationconsistingof a coldNS coresurroundedby we demonstrate that it is thermal pressure alone that
a hot mantle. Hence, all remnants are Thorne-Zytkow- preventsthese objects fromundergoing promptcollapse,
likeObjects(TZlOs). ScalingourresultstorealisticWD sinceangularmomentumsupportiscompletelyabsentin
compactions,wepredictthatarealistichead-oncollision head-on collisions. Delayed collapse occurs on a cooling
will form a quasiequilibrium TZlO. time scale in all cases, providing a consistency check on
Although the head-on collision simulations appear to our cooling implementation.
lead to a consistent result (the formation of a TZlO),
these results cannot be used to predict the final fate of
WDNS systems in circular orbit. On the one hand, one Our pWDNS merger calculations show that the inspi-
might expect that the remnant in the inspiraling case ralremnantisaspinningTZlOwhichissurroundedbya
will collapse to a black hole, because shock heating is massive, extended, hot disk. In contrast to our head-on
not likely to be as intense as in the head-on case. On collisions, we do not find any outflows in the inspiral-
the other hand, the large amount of angular momen- ing case. Therefore, the final total mass is greater than
tum in the inspiraling binary case may work to pre- the maximum mass supportable by our cold EOS and
vent prompt collapse. Therefore, to predict whether the manynuclearEOSs. However,theremnantdoesnotcol-
merged WDNS remnant will collapse, promptly or fol- lapse promptly to a BH. We find that the remnant is
lowing cooling, we need to perform fully general rela- both thermally and centrifugally supported. To deter-
tivistic simulations of WDNS binaries through inspiral mine whether centrifugal forces alone can support the
and merger. remnant we incorporate cooling and find that the object
does not collapse to a black hole. Therefore, the extra
support provided by rotation is sufficient for preventing
B. Goals and objectives the collapse.
The purpose of the current work is threefold:
a) We simulate the late inspiral and merger of a WDNS Eventhoughthe TZlOdoes notcollapseafter cooling,
systemconsistingofa 1.4M⊙ NS anda0.98M⊙ WDini- we expect delayed collapse ultimately because the final
tially in circular orbit and at the Roche limit. As in our total rest mass ( 2.5 M⊙) is larger than the maximum
∼
head-on collision studies, we employ the pWD approxi- possible mass supportable by our cold EOS (and many
mationtomakethecomputationsfeasible[64]. ThepWD nuclearEOSs), evenallowingfor uniform rotation. (The
approximation is useful for predicting the ultimate fate maximumgravitationalmassofauniformlyrotatingstar
of a realistic WDNS merger using scaling. In particular, withouradoptedEOSis 2.1M⊙). Weexpectthatcol-
≃
thecollisionvelocity(v )andthepre-shockedWDsound lapse to a BH will take place after viscosity or magnetic
c
speedc bothscaleas (M/R )1/2. Thisimpliesthat fields redistribute the angularmomentum, as in the case
s WD
∼
the Mach number ( = v /c ) is invariant under scal- of a hypermassive neutron star [65–67]. his conclusion
c s
M
ing of R and so is the degree of shock heating. The will be true in the case of realistic WDNS mergers, un-
WD
thermal energy, as well as the rotational kinetic energy less the true nuclear EOS supports a uniformly rotating
(T) and the gravitational potential energy (W) all scale starwitharestmassexceedingtheremnantmass. Many
as M2/RWD, when the binary merges. Thus, T/W viableEOSsdonotsupportrestmassesaslargeas2.5M⊙
∼ | |
is also invariant under scaling of R . These consider- [37], the remnant rest mass in our simulations.
WD
ations simply mean that with respect to gravity the rel-
ative importance of thermal and rotational support in a
WDNSmergerremnantisapproximatelyinvariant,when This paper is organized as follows. In Sec. II the
the masses of the binary components are fixed and the pWD approximation and the EOS adopted in our simu-
only quantity that changes is the WD radius. As a con- lations are briefly reviewed. Section III outlines the ini-
sequence, the results obtained when adopting pWDNS tial data generation technique. Sec. IV summarizes the
systems can be scaled up to realistic WDNS systems. methods used for evolving the gravitational and matter
Notethatourcompactionstudyin[1]confirmstheabove fields. Sec. V introduces our radiative cooling formal-
scaling with the Mach number. ism, which is then applied to the TZlO remnants from
b) We introduce a radiative cooling prescription and our pWDNS head-on collision simulations in Sec. VI.
modify our adiabatic simulations by allowing for cooling We present the results of our fully relativistic hydrody-
to determine whether the merger remnant will collapse namicsimulationsofthebinarypWDNSlateinspiraland
without thermal support, if it fails to collapse promptly. merger in Sec. VII, and turn on cooling in Sec. VIIC. In
Otherwise, angular momentum provides sufficient sup- Sec. VIII we discuss possible effects of nuclear reactions
port to prevent collapse. in realistic WDNS mergers and give estimates of realis-
c) We allow cooling to occur in the TZlOs formed in tic cooling and angular momentum redistribution time
our WDNS head-on collision simulations [1] to confirm scales. Sec. IX concludes with a summary of the main
that these remnants collapse to a black hole when the findings. Throughout this work, geometrized units are
excess thermal energy is radiated away. In other words, adopted, where G=c=1, unless otherwise specified.
4
II. EQUATION OF STATE A. Gravitational Field Equations
We employ the following 6-parameter piecewise poly- The spacetimemetric in the standard3+1form[70]is
tropic cold EOS: written as
ds2 = α2dt2+γ (dxi+βidt)(dxj +βjdt), (3)
κ ρ1/n1, ρ 6ρ − ij
1 0 0 1
ρP0 =κ2ρ01/n2, ρ1 <ρ0 6ρ2 , (1) wtthimheeerteht.rαeTei-shmrthoetuergilchapoousnetfstuphnaecctpeilaoipknee,rβhLiyapttheirensusihnrfidfaticcevesescortfuorncoafnnrosdtmaγni1jt
κ ρ1/n3, ρ >ρ to 3, and Greek indices run from 0 to 3.
 3 0 0 2 The three-metric γij is then conformally decomposed
as
where P is the pressure, ρ is the rest-mass density
0
and κ ,κ ,κ ,n ,n ,n ,ρ ,ρ are the parameters of the
1 2 3 1 2 3 1 2 γ Ψ4f , (4)
EOS. The parameters in Eq. (1) are 8 in number, but ij ij
≡
continuityrequiresthatthe followingconditionsbe true
where Ψ is the conformal factor and f the conformal
ij
metric. We adopt the standard approximation of a con-
κ1 =κ2ρ11/n2−1/n1, κ2 =κ3ρ12/n3−1/n2. (2) formally flat spacetime, so that fij = δij in Cartesian
coordinates.
As a result, the adopted EOS has 6 free parameters We split the extrinsic curvature (Kij) into trace (K)
κ ,n ,n ,n ,ρ , and ρ . and tracefree parts (Aij)
3 1 2 3 1 2
Because of its multiple parameters, this EOS gives us
1
thefreedomtocapturethesamecharacteristiccurvesand Kij =Aij + γijK, (5)
turningpointsonaTOVmass-centraldensityplotasfor 3
a realistic cold-degenerate EOS (see [68]), as shown in
take the initial slice to be maximal
Fig. 1in[1]. The EOSexhibitsbothstable(dM/dρ >
0,c
0) and unstable (dM/dρ < 0) branches for both WDs
0,c K =0, (6)
and NSs, as in the realistic case.
Furthermore, this EOS allows us to adjust the size of and introduce a “conformal”, traceless extrinsic curva-
a pWD of any given mass, thereby shifting the pWD ture as
branch to smaller radii (see Fig. 2 in [1]), while keeping
the NS branch approximately unchanged. For more de- A¯ij Ψ10Aij. (7)
tails about our EOS and pWDs we refer the interested ≡
reader to [1]. Using Eqs. (4)–(7) and assuming the existence of
In this work the EOS parameters correspond to the an approximate helical Killing vector, the Hamiltonian
10:1 EOS we considered in [1]: κ = 4993, Γ = and momentum constraintequations assume the form of
3 1
1.515, Γ = 2.969, Γ = 0.714, log(ρ /ρ ) = 2.268, the conformal-thin-sandwich (CTS) equations [12]. The
2 3 1 nuc
−
log(ρ /ρ )=0.208,whereallvaluesareingeometrized Hamiltonian constraint becomes
2 nuc
unitsandΓ =1+1/n ,ρ =1.485 10−4km−2. These
i i nuc
× 1
parametersarechosensuchthattheratiooftheisotropic ¯2Ψ= Ψ−7A¯ A¯ij 2πΨ5ρ. (8)
ij
radiusofaTOV0.98M⊙ pWDtothatofaTOV1.5M⊙ ∇ −8 −
NSis10:1. Inaddition,theEOShasbeenconstructedso
where ¯2 is the flat Laplacian operator associated with
thatthemaximumgravitationalmassofaNSis1.8M⊙, ∇
f . Here the source term ρ is defined as
i.e., the same as that for the AP2 version of the Akmal- ij
Pandharipande-Ravenhall (APR) EOS [39, 69], and the
ρ nαnβT , (9)
maximum gravitational mass of a pWD is 1.43 M⊙, i.e., ≡ αβ
the maximum mass of a TOV WD obeying the Chan-
where nα is the normal vector to a t = constant slice,
drasekhar EOS for mean molecular weight µ =2.
e
and T is the stress-energy tensor of the matter.
αβ
The momentum constraint yields
1
III. INITIAL DATA ¯2βi+ ¯i(¯ βj)=2A¯ij¯ (αΨ−6)+16παΨ4ji, (10)
j j
∇ 3∇ ∇ ∇
Thissectionintroducestheformalismadoptedforgen- where the source term ji is given by
erating valid general relativistic initial data for binary
pWDNS systems in circular orbit. jα γαn Tβγ. (11)
≡− β γ
5
B. Matter fields
TABLE I. Outer boundary conditions imposed on the CTS
variables when generating WDNSinitial data.
As we argued in [34] the WD in a WDNS binary with
Variable Fall-off condition
close separation likely will be tidally locked. For this
Ψ−1 ∼1/r reason we focus on corotating WDNS systems only.
α−1 ∼1/r We assume that the matter is described by a perfect
Gx ∼y/r3 fluid stress-energy tensor:
Gy ∼x/r3
Tαβ =(ρ +ρ +P)uαuβ +Pgαβ, (18)
Gz ∼xyz/r7 0 i
B ∼xy/r3 where gαβ is the inverse of the four-metric and
ρ ,ρ ,P,uα are the rest-mass density, internal energy
0 i
density, pressure, and four-velocity of the fluid respec-
Takingthe traceofthe evolutionequationforKij (see tively. For allinitial configurations,the pressure is given
Eq. (2.106)in [12]), imposing the maximal slicing condi- by the cold EOS as specified in Eq. (1). The internal
tion Eq. (6), and combining the result with Eq. (8), we energy density can be derived by integrating
obtain an equation for the lapse [12]
ρ 1
i
d = Pd , (19)
¯2(αΨ)=αΨ 7Ψ−8A¯ijA¯ij +2πΨ4(ρ+2S) . (12) (cid:18)ρ0(cid:19) − (cid:18)ρ0(cid:19)
∇ 8
(cid:18) (cid:19)
and for Eq. (1) the integration yields
Here the source term S is defined as
n κ ρ1/n1, ρ 6ρ
1 1 0 0 1
S γijT . (13)
ij
≡ ρ 
In all equations above ρ0i = n2κ2ρ01/n2 +c2, ρ1 <ρ0 6ρ2 (20)
Ψ6 2 nκρ1/n3 +c , ρ >ρ
A¯ij = 2α (cid:18)∇¯iβj +∇¯jβi− 3fij∇¯kβk(cid:19), (14) where  0 3 0 2
and A¯ =f f A¯mn.
ij im jn c =(n n )κ ρ1/n1, c =c +(n n )κ ρ1/n2,
InsteadofsolvingEq. (10)fortheshiftvectordirectly, 2 1− 2 1 1 3 2 2− 3 2 2
itisconvenienttodecomposeβi asasumofavectorand (21)
In Cartesian coordinates we choose the orbital plane
a gradient (cf. [71])
of the binary to be the z = 0 plane, so that the fluid
1 four-velocity takes the form [12]
βi Gi ¯iB. (15)
≡ − 4∇
uα =ut(1, Ωy,Ω(x x ),0), (22)
rot
− −
Eq.(10)canthenbereplacedbythetwoequivalentequa-
tions where Ω is the constantorbitalangularvelocityand xrot
isthexcoordinateoftheaxisofrotation. Following[72],
¯2Gi =2A¯ij¯ (αΨ−6)+16παΨ4ji (16) we introduce a vector
j
∇ ∇
and ξα =(0, y,(x xrot)), (23)
− −
¯2B = ¯ Gi. (17) and rewrite the four-velocity as
i
∇ ∇
uα =ut(αnα+Ωξα+βα). (24)
Equations (8), (12), (16) and (17) form a system of 6
coupled, nonlinear elliptic equations for the 6 unknowns
Ψ,αΨ,GiandB,whichmustbesolvediteratively. These The source term ρ in Eq. (9) can then be written
equations are elliptic and hence require outer boundary
ρ +ρ +P
conditions to be specified. We impose the same fall- ρ= 0 i P, (25)
1 v2 −
off boundary conditions as in [72], except that here we −
choose the binary components to be initially lined up on
wherevisthemagnitudeofthethree-velocityofthefluid.
the x-axis and the binary rotation axis parallel to the Using uαu 1, it can be shown that v2 is given by
α
z-axis. Table I lists the full set of outer boundary condi- ≡−
tions imposed in our initial data. v2 = Ψ4 (Ωy βx)2+ Ω(x x )+βy 2+(βz)2 .
α2 − − rot
h (i26)
(cid:0) (cid:1)
6
The momentum source ji in Eq. (11) becomes corotational NS at such large separations. The spin of
a corotating NS is very small. To understand this, con-
ji = (ρ0+ρi+P)(Ωξi+βi), (27) sider the ratio of the angular velocity of the corotating
α 1 v2 NS (Ω ) to that at the mass-shedding (Ω ) limit:
cor ms
−
and S in Eq. (13) is given by 3/2 3/2
Ω M R R
cor total NS NS
= 1.3 , (29)
S =(ρ +ρ +P) v2 +3P. (28) Ωms r MNS (cid:18) AR (cid:19) ≈ (cid:18) AR (cid:19)
0 i 1 v2
− where A is the Roche limit separation. For the typi-
R
cal system we consider R /A R /3R . For re-
NS R NS WD
≈
C. Computational methods alistic massive WDs RWD/RNS 500, and for pWDs
R /R 10. Thus, Ω /Ω≈ 10−5 for realistic
pWD NS cor ms
WDNSs and≈Ω /Ω 10−2 for pW≈DNSs. Therefore,
We solve the nonlinear elliptic equations (8), (12), cor ms
≈
thecorotationspintheNSacquiresisverysmallandhas
(16) and (17) using a fixed-mesh-refinement (FMR) fi-
no physical significance.
nite difference code we developed, which is based on the
HavingprescribedtheNSandpWDrest-massdensity,
Portable, Extensible Toolkit for Scientific Computation
usingsecond-orderpolynomialinterpolation,weinterpo-
(PETSc) library [73–75]. A full description of our code
late the NS and pWD matter distribution on the grid of
may be found in [1]. Here we summarize the basic fea-
tures. our FMR elliptic code and solve the CTS equations.
The grid structure used in our FMR elliptic code is
a multi-level set of properly nested, uniform grids. We
use standard cell-centered, second-order accurate finite IV. EVOLUTION OF WDNS SYSTEMS
difference stencils for the Laplacian operator and the
derivatives of the variables, using first-order interpola- A. Basic Equations
tion across the refinement level boundaries when neces-
sary. We calculate the solution across the entire grid,
Theformulationandnumericalschemeforoursimula-
and only on leaf cells (i.e. cells within which there exist
tions are the same as those reported in [1, 23, 62, 76], to
no higher resolution cells). In [1] we performed a series
whichthereadermayreferfordetails. Hereweintroduce
of tests involving single NSs, and we demonstrated that
our notation and summarize our method.
the code converges to the expected solutions at second
We use the 3+1 formulation of general relativity, in
order.
which the metric is decomposed as in Eq. (3). In this
Giventhe matter distribution,Ω andx we solvethe
rot formalism, the fundamental dynamical variables for the
CTSequationsiteratively,addressingthenon-linearityof
metric evolution are the spatial three-metric γ and ex-
ij
Eq. (8) by performing Newton-Raphson iterations, until
trinsic curvature K . The Baumgarte-Shapiro-Shibata-
ij
the residuals of all six equations become smaller than
Nakamura(BSSN)formalism[12,77,78]isadopted. The
some set tolerance (usually set to 10−15).
BSSN evolution variables are the conformal exponent
We obtain the WD rest-mass density distribution, Ω φ ln(γ)/12, the conformal 3-metric γ˜ = e−4φγ ,
ij ij
and x at the Roche limit for equilibrium, corotating ≡
rot three auxiliary functions Γ˜i γ˜ij , the trace of the
binaryWDNSsincircularorbitobeyingourcoldEOSus- ≡ − ,j
extrinsiccurvatureK,andthetrace-freepartofthecon-
ing the unigridNewtoniancode we developedandtested
formal extrinsic curvature A˜ e−4φ(K γ K/3).
in[34]. AttheRochelimit,thebinaryseparationislarge ij ≡ ij − ij
Here γ = det(γ ). The full spacetime metric g is re-
enough so that the tidal effects on the NS are negligible, ij µν
latedtothethree-metricγ byγ =g +n n ,where
and hence the NS will be spherical to a high degree and µν µν µν µ ν
thefuture-directed,timelikeunitvectornµ normaltothe
point-likefromthepointofviewoftheWD.Thus,inthe
timeslicecanbewrittenintermsofthelapseαandshift
Newtonian code we model the NS as a point mass and
βi asnµ =α−1(1, βi). Theevolutionequationsofthese
we self-consistently solve for the WD rest-mass density −
BSSN variables are given by Eqs. (9)–(13) in [23].
distribution via the integrated Euler equation. We use
We adoptstandardpuncture gauge conditions: an ad-
the Newtonianequationsforthis step,because itiscom-
vective “1+log” slicing condition for the lapse and a
putationally simple and fast. Also, the large separation
“Gamma-freezing” condition for the shift [79]. Thus, we
at the Roche limit ensures that the WD and NS inter-
have
action lies in the Newtonian regime, so that our initial
configuration is nearly in equilibrium.
∂ α= 2αK , (30)
0
After the WD rest-mass density distribution has been −
∂ βi =(3/4)Bi , (31)
calculated, the point-mass NS is replaced by a TOV NS 0
with gravitational mass equal to that of the point-mass ∂ Bi =∂ Γ˜i ηBi , (32)
0 0
−
NS, centered at the position of the point mass. For sim-
plicity, we model the NS as corotational because there where ∂0 ∂t βj∂j. We set the η parameter to
is no essential difference between an irrotational and a 0.01 km−1≡for all−simulations presented in this work.
7
The fundamental matter variables are the rest-mass that employs the piecewise parabolic (PPM) reconstruc-
density ρ , specific internal energy ǫ, pressure P, and tion scheme [82], coupled to the Harten, Lax, and van
0
four-velocity uµ. We write the stress-energy tensor as Leer (HLL) approximate Riemman solver [83]. The
adopted hydrodynamic scheme is second-order accurate.
Tµν =ρ0huµuν +Pgµν , (33) To stabilize our hydrodynamic scheme in regions where
there is no matter, a tenuous atmosphere is maintained
whereh=1+ǫ+P/ρ0isthespecificenthalpyandǫisthe on our grid, with a density floor ρ set to 10−10 times
atm
specificinternalenergy. Inournumericalimplementation
the initial maximum density on our grid. The average
of the hydrodynamics equations, we evolve the following density of the pWD is 1011gr/cm3, and at least six or-
“conservative” variables:
dersofmagnitudelargerthanthatofthe artificialatmo-
ρ∗ √γρ0nµuµ , (34) sphere. Thus,theatmosphereposesnoprobleminevolv-
≡− ingthepWD.TheinitialatmosphericpressureP isset
S˜i ≡−√γTµνnµγνi , (35) byusingthecoldEOS(1). Throughouttheevoluattmion,we
τ˜ √γTµνnµnν ρ∗ . (36) imposelimitsonthepressuretopreventspuriousheating
≡ −
andnegativevalues ofthe internalenergyǫ. Specifically,
The evolution equations for these variables are given by
we require P P P , where P =10P and
min max max cold
Eqs. (27)–(29) in [62]. ≤ ≤
P = 0.8P , where P is the pressure calculated
min cold cold
TheEOSweadoptfortheevolutionhasbothathermal
using the cold EOS (1). Whenever P exceeds P or
max
and cold contribution, and can therefore be written
drops below P , we reset P to P or P , respec-
min max min
tively. Following[24]weimposetheupperpressurelimits
P =P +P , (37)
th cold only inregionswhere the rest-massdensity remainsvery
low(ρ <100ρ ),butweimposethelowerlimitevery-
where P is given by Eq. (1) and the thermal pressure 0 atm
cold
where on our grid. We impose the pressure floor every-
is given by
where, because numericalerrorsometimes leads ǫ ǫ
cold
−
P =(Γ 1)ρ (ǫ ǫ ), (38) slightly below zero, resulting in negative thermal pres-
th th 0 cold
− − sure. Wehavefoundexperimentallythatifthissituation
where arises, it can be avoided in the subsequent timesteps by
imposing the pressure floor.
ǫcold = Pcoldd(1/ρ0). (39) At eachtimestep, the “primitive variables”ρ0, P, and
− vi must be recovered from the “conservative” variables
Z
We set Γ = 1.66 ( 5/3) in all our simulations. That ρ∗, τ˜, and S˜i. We perform the inversion numerically as
th
≃ specifiedin[62]. Weusethesametechniqueasin [84,85]
is,wesetΓ totheΓ exponentofthe10:1EOS,appro-
th 1 toensurethatthevaluesofS˜ andτ˜yieldphysicallyvalid
priateeitherfornonrelativisticcold,degenerateelectrons i
primitive variables.
or (shock) heated, ideal nondegenerate baryons. Equa-
tion (37) reduces to our piecewise polytropic law Eq.(1)
for the initial (cold) NS and pWD matter.
C. Diagnostics
B. Evolution of the metric and hydrodynamics
During the evolution, we monitor the normalized
Hamiltonian and momentum constraints as defined in
We evolve the BSSN equations using fourth-order ac- Eqs. (40)–(43) of [23]. We also monitor the ADM mass
curate,cell-centeredfinite-differencingstencils,excepton and angular momentum of the system. The equations
shift advection terms, where fourth-order accurate up- used to calculate the ADM mass and angular momen-
wind stencils are applied. We apply Sommerfeld out- tum with minimal numerical noise are as follows [12]
going wave boundary conditions on all BSSN fields, as
1 1
in [23]. Our code is embedded in the Cactus paral- M = d3x ψ5ρ+ ψ5A˜ A˜ij Γ˜ijkΓ˜
ij jik
16π − 16π
lelization framework [80], and our fourth-order Runge- Z (cid:18)
Kutta timestepping is managed by the MoL (Method of 1 ψ 1
+ − R˜ ψ5K2 , (40)
Lines) thorn, with the Courant-Friedrichs-Lewy num- 16π − 24π
(cid:19)
ber set to 0.45 in all pWDNS simulations. We use the
Carpet [81] infrastructure to implement the moving-box 1 2
J = ǫ n d3x ψ6(A˜j + xj∂ K
adaptive mesh refinement. In all AMR simulations pre- i 8π ij n 3 n
sented here, we use second-order temporal prolongation, Z
1 (cid:2)
coupledwithfifth-orderspatialprolongation,andimpose xjA˜ km∂nγ˜ km)+8πxjSn . (41)
−2
equatorial symmetry to reduce the computational cost.
We write the generalrelativistic hydrodynamics equa- Here ψ = eφ, ρ = nµnνTµν, Si = nµγiνT(cid:3)µν, R˜ is the
−
tions in conservative form. They are evolved via a high- Ricci scalar associated with γ˜ , and Γ˜ are Christoffel
ij ijk
resolution shock-capturing (HRSC) technique [62, 76] symbols associated with γ˜ .
ij
8
TABLE II. Summary of initial configurations. M (M ) is the ADM mass of an isolated NS (pWD)(a), R (R ) the
NS WD NS WD
isotropic radius of an isolated NS (pWD), C the compaction of an isolated NS, where the compaction is the ratio of the
NS
ADMmassoftheisolatedstartoitsarealradius. AllpWDsconsideredherehaveC =0.01. M istheADMmassofthe
WD ADM
systemandAtheinitialbinaryseparationinisotropiccoordinates. CasesA1andA3arethesameasthehead-oncollisioncases
we studied in [1]. The initial coordinate separation for these cases was set to 587 km. Case A corresponds to our simulation
of a binary pWDNS in circular orbit starting at the Roche limit. For this case ΩM = 6.95×10−4. All cases have been
ADM
produced with the10:1 EOS of [1].
Case MNS/M⊙ MWD/M⊙ CNS J/MA2DM RWD/RNS RWD/MADM MADM/M⊙ A/RWD
A1 1.4 0.98 0.11 0. 8.88 41.18 2.41 4.00
A3 1.6 0.98 0.15 0. 11.15 37.46 2.65 4.00
A 1.4 0.98 0.11 2.88 8.88 40.05 2.48 3.14
(a) Here we list the ADMmasses, isotropic radii and compactions of the isolated NSstars, whose rest-mass density profiles
were used to generate initial data. Thesame holds for the pWDsin cases A1 and A3. In case A the pWDrest-mass density
profile,theRochelimitseparationandΩweregeneratedbyaNewtonianbinaryWDNScodeandthenusedinourCTSsolver.
When hydrodynamic matter is evolved on a fixed uni- The dynamics of radiation is governedby [86–88]
formgrid,ourhydrodynamicschemeguaranteesthatthe
Rαβ = Gβ, (43)
rest mass M is conserved to machine roundoff error. α
0 ∇ −
This strict conservation is no longer maintained in an whereRαβ isthe radiationstress-energytensorgivenby
AMR grid, where spatial and temporal prolongation is
performed at the refinement boundaries. Hence, we also
Rαβ = dνdΩI NαNβ, (44)
monitor the rest mass ν
Z
and Gα is the radiation four-force density given by
M0 = ρ∗d3x (42)
Z Gα = dνdΩ(χ I j )Nα. (45)
ν ν ν
−
during the evolution. Rest-massconservationis alsovio- Z
lated wheneverρ0 is resetto the atmospherevalue. This In the equations above dΩ is the solid angle, ν and
usuallyhappensonlyintheverylow-densityatmosphere. I = I (xα,Ni,ν) are the radiation frequency and spe-
ν ν
The low-density regions do not affect rest-mass conser- cific intensity of radiation at xα moving in direction
vation significantly. Nα = pα/hν, respectively. All quantities are measured
In all simulations we present in this work the nor- inthelocalLorentzframeofafiducialobserverwithfour-
malizedHamiltonianconstraintviolationsremainsmaller velocity uα , i.e.,
fid
than 0.9% and the normalized momentum constraint vi-
olations smaller than 2.1%. Rest mass is conserved to hν =−pαuαfid, (46)
within 4% and angular momentum to within 10%.
where pα is the photon four-momentum and h denotes
Shocks occur when the stars collide. We measure
Planck’s constant. The energy-momentum conservation
the entropy generated by shocks via the quantity K
≡ equation then becomes
P/P 1, where P is the pressure associated with
cold cold
≥
the cold EOS that characterizes the initial matter (see (Tαβ +Rαβ)=0 (47)
α
Eq. (1)). ∇
or after using Eq. (43)
Tαβ =Gβ. (48)
α
V. RADIATIVE COOLING ∇
After projecting this equation using the fluid four-
velocity uα, we obtain the modified energy equation:
Our binary WDNS head-on collision studies in [1]
demonstrate that the hot, quasiequilibrium TZlO rem- uα ε= (ε+P) uα uαG , (49)
α α α
nants do not collapse promptly to a black hole, even ∇ − ∇ −
though the final total mass is larger than the maximum wheretheperfectfluidstress-energytensorhasbeenwrit-
mass supportable by the cold EOS. This outcome might ten as
also arise in the case of WDNS mergers in circular or-
Tµν =(ε+P)uµuν +Pgµν, (50)
bit, and can be due to additional support provided by
thermal pressure and/or rapid rotation. In order to de-
and ε is the total energy density. Using the continuity
termine whether thermal support is dominant, we add
equation
radiativecoolingtotheGRhydrodynamicequations. We
now describe our formalism for implementing this. (ρ uα)=0, (51)
α 0
∇
9
TABLE III. Grid configurations used in our simulations. Here M is the sum of the ADM masses of the isolated stars, ∆x
is the grid spacing in the innermost refinement box surrounding the NS, N denotes the number of grid points covering the
NS
diameteroftheNSinitially,andN denotesthenumberofgridpointscoveringthe(smallest)diameterofthepWDinitially.
WD
The outer boundarydistance to the centerof mass is approximately 1020M in cases A1 and A3, and 540M in case A.
Case M/M⊙ Grid Hierarchy (in units of M)(a) ∆x NNS NWD
A1 2.38 (534.33, 267.16, 133.58, 66.79, 35.78, 19.08, 10.44, 7.156) M/6.71 63 35
A3 2.58 (467.27, 233.64,116.82, 58.41, 29.20 , 15.58, 8.518, 5.841) M/8.22 56 38
A 2.38 (270.87, 135.44, 67.72, 36.28[N/A], 19.35[N/A], 9.674[N/A], 5.744[N/A]) M/13.2 124 73
(a) There are two sets of nested refinement boxes: one centered on theNS and one on thepWD.This column specifies the
half-length of thesides of therefinement boxescentered on both theNS and pWD.When there is no corresponding pWD
refinement box (as thepWD is muchlarger than the NS),we write [N/A] for that box.
Eq. (49) becomes quasiequilibriumsphericalobjectwhichisthermallysup-
ported, and we cool it quasistatically,i.e., choose a cool-
uα ε= ε+Puα ρ uαG . (52) ing time scale much longer than the free-fall time scale
α α 0 α
∇ ρ0 ∇ − of the star, it will radiate thermal energy away and con-
tract. While the contractiongeneratesextraheat, radia-
Thetotalenergydensityisrelatedtothespecificthermal
tion tends to removeany residual thermal energy. Thus,
energy via the following equation
after cooling, TZlOs are expected either to collapse to
a BH or be supported by the residual cold pressure and
ε=ρ (1+ǫ +ǫ ). (53)
0 th cold centrifugal force.
In the optically thin regime, assuming the fluid radi-
UsingEqs.(52)and(53)wefindthatthespecificthermal
ates isotropically in its rest frame, and using uα = uα,
energy evolves as fid
the source term Gα is generally expressed as [87, 88]
P 1
uα ǫ = thuα ρ uαG , (54)
∇α th ρ20 ∇α 0− ρ0 α Gα =uα(Γ Λ)+ dν(χa +χs)Fα, (58)
− ν ν ν
Z
where we have used dǫ /dρ = P /ρ2, and P =
cold 0 cold 0 th
P P . whereΓandΛaretheheatingandcoolingterms,respec-
cold
−In the comoving reference frame uα =d/dτ, where tively, given by
α
∇
τ is the proper time. Thus, in the comoving frame Eq.
(54) becomes Γ= dνdΩχaI ,
ν ν
Z (59)
d P d 1
dτǫth = ρt2hdτρ0− ρ uαGα. (55) Λ= dνdΩjν,
0 0 Z
In order to achieve cooling, the radiation term in (55) and where χa,χs are the absorption and scattering co-
ν ν
mustbespecifiedsothatthermalenergycanberemoved. efficients, respectively, and j is the emissivity. Finally,
ν
For this reason we choose the following cooling law that Fα is the total radiation flux four-vector
ν
gives rise to exponential cooling
uαG =ǫ ρ /τ , (56) Fνα =Pαβ dΩIνNβ, (60)
α th 0 c
Z
where τc > 0 is the cooling time scale. Substituting where the projection operator Pαβ is defined as
Eq. (56) in Eq. (55) we obtain
Pα δα +uαu . (61)
β β β
d (Γ 1)dρ 1 ≡
th 0
ǫ = − ǫ , (57)
th th
dτ ρ dτ − τ If we assume that there is no absorption and no scat-
(cid:20) 0 c(cid:21)
tering, Eq. (58) becomes
where we used Eq. (38) to substitute for the thermal
pressure. Gα = uαΛ. (62)
−
The first term in brackets on the RHS of Eq. (57) re-
sults from adiabatic compression or expansion. The sec- As a result,
ond term results from cooling and radiates away ther-
mal energy exponentially. Thus, if initially we have a u Gα =Λ. (63)
α
10
FIG.1. Snapshotsofrest-mass densityprofilesat selected timesfor head-oncase A1aftercooling isturnedon. Thecontours
represent the rest-mass density in the orbital plane, plotted according to ρ0 =ρ0,max10−0.68j−0.16 (j =0,1,...,7). The color
sequence dark red, red, orange, yellow, green, light green, blue and light blue implies a sequence from higher to lower values.
This roughly corresponds to darker grey-scaling for higher values. Here ρ0,max = 4.645ρnuc, where ρnuc = 2×1014 g/cm3.
These snapshots clearly demonstrate that theentireTZlO remnant collapses once cooling is turned on. HereM =2.38 M⊙ =
3.52 km=1.17×10−5 s is thesum of the ADMmasses of the isolated stars.
AstraightforwardcomparisonofEqs.(56)and(63)shows tions, which is precisely how cooling is implemented in
that the integrated emissivity in our cooling model is our HRSC code.
given by
Note that the optically-thin approximation employed
ρ here is valid for neutrino cooling in the WDNS merger
0
Λ= ǫ . (64)
τ th scenario(head-onorotherwise). Accordingtoouranaly-
c
sisin[1],thetemperaturesanddensitiesofthehotmantle
If we project Eq. (48) using the timelike unit vector ofa TZlO are suchthat thermalneutrino emissionlikely
nα normal to spacelike hypersurfaces and the projection willbethedominantsourceofcooling. ThediffuseTZlO
operatorhαβ =δαβ+nαnβ,wefindthatthe3+1GRHD mantle composed of the WD debris is optically thin to
equations become neutrinos, justifying the above approximation.
1 Finally,notethattheself-gravityoftheradiationisne-
∂ S˜ +∂ (α√γTj )= α√γTαβg α√γu Λ, (65)
t i j i 2 αβ,i− i glected,i.e.,weassumethatradiationdoesnotaffectthe
spacetime structure, so only the perfect fluid stress en-
and
ergy tensor contributes to the BSSN source terms. This
∂tτ˜+∂i(α2√γT0i ρ∗vi)=s α2√γu0Λ, (66) is a good approximation as long as the radiation energy
− − density is subdominant (n n Rαβ n n Tαβ). This
α β α β
≪
where we haveused Eq.(62) andΛ is givenby Eq. (64). is indeed the case in a WDNS scenario because the rest
Thus,coolingentersasasourcetermintheGRHDequa- mass dominates the mass-energy as can be inferred by