Head-on collisions of binary white dwarf--neutron stars: Simulations in full general relativity
Vasileios Paschalidis, Zachariah Etienne, Yuk Tung Liu, Stuart L. Shapiro
aa r X i v : . [ a s t r o - ph . H E ] M a r Head-on collisions of binary white dwarf–neutron stars:Simulations in full general relativity
Vasileios Paschalidis , Zachariah Etienne , Yuk Tung Liu , and Stuart L. Shapiro , Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801 Department of Astronomy and NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801
We simulate head-on collisions from rest at large separation of binary white dwarf – neutron stars(WDNSs) in full general relativity. Our study serves as a prelude to our analysis of the circularbinary WDNS problem. We focus on compact binaries whose total mass exceeds the maximum massthat a cold degenerate star can support, and our goal is to determine the fate of such systems. A fullygeneral relativistic hydrodynamic computation of a realistic WDNS head-on collision is prohibitivedue to the large range of dynamical time scales and length scales involved. For this reason, weconstruct an equation of state (EOS) which captures the main physical features of NSs while, atthe same time, scales down the size of WDs. We call these scaled-down WD models “pseudo-WDs(pWDs)”. Using these pWDs, we can study these systems via a sequence of simulations where thesize of the pWD gradually increases toward the realistic case. We perform two sets of simulations;One set studies the effects of the NS mass on the final outcome, when the pWD is kept fixed. Theother set studies the effect of the pWD compaction on the final outcome, when the pWD mass andthe NS are kept fixed. All simulations show that after the collision, 14%-18% of the initial totalrest mass escapes to infinity. All remnant masses still exceed the maximum rest mass that our cold
EOS can support (1 . M ⊙ ), but no case leads to prompt collapse to a black hole. This outcomearises because the final configurations are hot . All cases settle into spherical, quasiequilibriumconfigurations consisting of a cold NS core surrounded by a hot mantle, resembling Thorne-Zytkowobjects. Extrapolating our results to realistic WD compactions, we predict that the likely outcome ofa head-on collision of a realistic, massive WDNS system will be the formation of a quasiequilibriumThorne-Zytkow-like object. PACS numbers: 04.25.D-,04.25.dk,04.40.Dg
I. INTRODUCTION
The inspiral and merger of compact binaries representsome of the most promising sources of gravitational waves(GWs) for detection by ground-based laser interferome-ters like LIGO [1, 2], VIRGO [3, 4], GEO [5], TAMA [6, 7]and AIGO [8], as well as by the proposed space-basedinterferometers LISA [9] and DECIGO [10]. Extractingphysical information from gravitational radiation emittedby compact binaries and determining their ultimate faterequires careful modeling of these systems in full generalrelativity (see [11] and references therein). Most effortto date has focused on modeling black hole–black hole(BHBH) binaries (see [12] and references therein), andneutron star–neutron star (NSNS) binaries (see [13] fora review), with some recent general relativistic work onblack hole–neutron star (BHNS) binaries [14–32].In this work we consider white dwarf–neutron star(WDNS) binaries in full general relativity. WDNS bi-naries are promising sources of low-frequency GWs forLISA and DECIGO and, as we argued in [33], possi-bly also high-frequency GWs for LIGO, VIRGO, GEO,TAMA and AIGO, if the remnant ultimately collapses toform a black hole.Like NSNS binaries, WDNS binaries are known to ex-ist. In [33] we compiled tables with 20 observed WDNSbinaries and their orbital properties. The NS massesin these systems range between 1 . M ⊙ and 2 . M ⊙ ,and their distribution is centered around 1 . M ⊙ . On the other hand, the WD masses in these systems rangebetween 0 . M ⊙ and 1 . M ⊙ , and their distribution iscentered around 0 . M ⊙ . Finally, 18 of these observedWDNS binaries have total mass greater than 1 . M ⊙ , ofwhich 8 have a WD component with mass greater than0 . M ⊙ , and 5 have total mass greater than 2 . M ⊙ . Thisis interesting because the expected TOV limiting massfor a cold, degenerate gas ranges between 1 . M ⊙ and2 . M ⊙ [34–42], depending on the equation of state (EOS)and degree of rotation, and one of the main goals of thiswork is to determine whether a WDNS merger can leadto prompt collapse to a black hole.Population synthesis calculations by Nelemans et al.[43] show that there are about 2 . × WDNS bina-ries in our Galaxy, and that they have a merger rate of1 . × − yr − . Furthermore, Nelemans et al. find thatafter a year of integration, LISA should be able to de-tect 128 WDNS pre-merger binaries and, after consider-ing the contribution of the double WD background GWnoise, resolve 38 of these. On the other hand, calcula-tions by Cooray [44] give much more conservative num-bers of resolvable WDNS binaries. In particular, Coorayfinds that the number of LISA-resolvable WDNS binariesranges between 1–10, using a WDNS merger rate between10 − yr − –10 − yr − . Finally, recent work by Thompsonet al. [45] suggests that the lower limit on the merger ratein the Milky Way, at 95% confidence, is 2 . × − yr − .Thompson et al. also suggest that the merger rate in thelocal universe is ∼ . − × Gpc − yr − . Therefore,leaving some uncertainties aside, all recent work on pop-ulation synthesis suggests that LISA should be able todetect a few WDNS pre-mergers per year.We note here that Newtonian work on binaries with aWD component has been performed analytically in [33,46–50] and via Newtonian hydrodynamic simulations in[51–56].In [33] we focused on WDNS binaries that have spiraledsufficiently close that they reach the termination pointfor equilibrium configurations. This is the Roche limitfor WDNSs, at which point the WD fills its Roche lobeand may experience one of at least two possible fates:i) stable mass transfer (SMT) from the WD across theinner Lagrange point onto the NS, or ii) tidal disruptionof the WD by the NS via unstable mass transfer (UMT).We also studied the key parameters that determinewhether a system will undergo SMT or UMT and foundthat, for a given NS mass, there exists a critical massratio q crit ≈ / q = M WD /M NS of a WDNSsystem is such that q > q crit , the WD quickly overfillsits Roche lobe, and the binary will ultimately undergoUMT. In the opposite case, q < q crit , the system will un-dergo SMT. We showed that a quasistationary treatmentis adequate to follow the evolution of an SMT binary dur-ing this secular phase and calculated the gravitationalwaveforms. We also pointed out that WDNS observa-tions suggest that there are candidates residing in bothregimes.In the case of tidal disruption (UMT), by contrast, thesystem will evolve on a hydrodynamical (orbital) timescale. In this scenario the NS may plunge into the WDand spiral into the center of the star, forming a quasiequi-librium configuration that resembles a Thorne-Zytkowobject (TZO) [57]; alternatively, the NS may be the re-ceptacle of massive debris from the disrupted WD.Depending on the details of the EOS, a cold degenerategas can support a maximum gravitational mass between1 . M ⊙ and 2 . M ⊙ against catastrophic collapse, if it isnot rotating (the TOV limit). It can support 20% moremass, if it is rotating uniformly at the mass-sheddinglimit (a “supramassive NS” [36]), and about 50% moremass, if it rotates differentially (a “hypermassive NS”[34–36]). If the total mass of the merged WDNS exceedsthe maximum mass supportable by a cold EOS, delayedcollapse to a black hole is inevitable after the remnantcools. However, the ultimate fate of the merged WDNSdepends on the initial mass of the cold progenitor stars,the degree of mass and angular momentum loss duringthe WD disruption and binary merger phases, the an-gular momentum profile of the WDNS remnant and theextent to which the disrupted debris is heated by shocksas it settles onto the NS and forms an extended, massivemantle. These are issues that require a hydrodynamicsimulation to resolve. Moreover, ascertaining whether ornot the neutron star ultimately undergoes a catastrophiccollapse (either prompt or delayed) to a black hole re-quires that such a simulation be performed in full gen- eral relativity. In fact, even the final fate of the NS inthe alternative scenario in which there is a long epochof SMT may also lead to catastrophic collapse, if theneutron star mass is close to the neutron star maximummass, and this scenario too will require a general rela-tivistic hydrodynamic simulation to track.In this paper we employ the Illinois adaptive mesh re-finement (AMR) relativistic hydrodynamics code [23, 58]to perform our first simulations of these alternative sce-narios. In particular, we study the head-on collision fromrest at large separation of a massive WD and a NS as aprelude to our investigation of the circular binary prob-lem, which we will report in a future work. We focus oncompact objects whose total mass exceeds the maximummass supportable by a cold EOS to determine whethersuch a collision leads to prompt collapse of the remnant,or a hot gaseous mantle composed of WD debris sur-rounding a central NS – a Thorne-Zytkow-like object(TZlO).The vast range of dynamical time scales and lengthscales involved in this problem make fully general rela-tivistic simulations extremely challenging. For example,a near-Chandrasekhar-mass WD has a radius R WD ≃ km and dynamical time scale of about 1s. On theother hand a typical NS has a radius of order R NS ≃ . M ⊙ ) butdifferent compactions, while the compaction and mass ofthe NS involved remain practically unchanged. In otherwords, while keeping the masses of the binary compo-nents and the NS radius fixed, we adjust the ratio of theradius of the pWD to that of the NS so that it varies from5:1 to 20:1 and then use our results to predict the out-come of the realistic case. The common feature among allversions of the piecewise EOS we employ is that the max-imum NS mass always is 1 . M ⊙ and the maximum WDmass always is 1 . M ⊙ , i.e., the Chandrasekhar mass.In addition to studying the effects of the pWD com-paction, we also study the effects of the NS mass. Weconsider NSs with masses 1 . M ⊙ . M ⊙ and 1 . M ⊙ .All simulations that we perform show that after thecollision, 14%-18% of the initial total rest mass escapes toinfinity. The remnant mass in all cases exceeds the maxi-mum rest mass that our cold EOS can support (1 . M ⊙ ),but we find that no case leads to prompt collapse to ablack hole. This outcome arises because the final con-figurations are hot . All our cases settle into a spher-ical quasiequilibrium configuration consisting of a coldNS core surrounded by a hot mantle. Hence, all rem-nants are TZlOs. Extrapolating our results to realisticWD compactions, we predict that the likely outcome ofa head-on collision from rest at large separation of a re-alistic massive WDNS system will be the formation of aquasiequilibrium TZlO.This paper is organized as follows. In Section II wereview the time scales and length scales involved in aWDNS merger and discuss why this problem presentssuch a computational challenge. In Section III we in-troduce the EOS adopted for our computations and de-scribe our pWD models. Section IV outlines our methodfor preparing initial data, and Section V summarizes ourmethods for evolving the gravitational and matter fields.We present the results of our fully relativistic hydrody-namic simulations in Section VI, and conclude in Sec-tion VII with a summary of our main findings. Through-out we use geometrized units, where G = c = 1. II. COMPUTATIONAL CHALLENGE
Simulating a WDNS merger in full general relativity isa difficult computational task. In this section we sketchin quantitative terms exactly why this is so.There are three fundamental time scales and lengthscales involved in the WDNS merger that must be re-solved. The relevant time scales are the dynamical timescale of each component of the binary and the orbital pe-riod; the relevant length scales are the NS and WD radiiand their orbital separation.Resolving the WD length scale and dynamical timescale is necessary in order to assess what happens to theWD at merger. Merger occurs on the orbital time scale,so this time scale must also be resolved. Resolving the NSdynamical time scale will enable us to assess whether theNS promptly collapses and forms black hole, or remainsinside the remnant WD, settling into a TZlO.The dynamical time scale of the NS, t d , NS , is given by t d , NS = s R M NS , (1)where R NS , and M NS are the characteristic NS radiusand mass respectively.Similarly, the dynamical time scale of the WD, t d , WD ,is given by t d , WD = s R M WD , (2)where R WD , and M WD are the characteristic WD radiusand mass respectively. Finally, the orbital time scale, T ,is given by T = 2 π s A M T , (3) where M T = M NS + M WD is the total mass, and A isthe orbital separation. Note that, at this separation, thehead-on collision time scale is T coll = π √ s A M T = T √ , (4)assuming the stars free-fall from rest as Newtonian pointmasses, and hence it is roughly the same order of magni-tude as T .By use of Eqs. (1) and (2) the ratio of the WD timescale to the NS time scale is t d , WD t d , NS = s R qR , (5)where q = M WD M NS (6)is the binary mass ratio.For M WD = 1 M ⊙ and using a cold degenerate electronEOS one finds R WD ≈ M NS = 1 . M ⊙ , R NS ≈ R WD R NS ∼ , (7)and from Eq. (5) the ratio of the dynamical time scalesis t d , WD t d , NS ∼ . (8)At the Roche limit, A is typically a few (two to five)times the WD radius [33]. Using, A ∼ R WD , and thevalues for the masses and radii used above, the orbitaltime scale becomes Tt d , NS = 2 π s A (1 + q ) R ∼ . (9)It is thus clear that there is a vast range of length scalesand time scales involved in this problem. The only wayto simulate the WDNS merger is by exploiting the powerof adaptive mesh refinement, so that resolution is highonly where required. However, even this does not sufficeto tackle the timescale problem as we explain below.Given that all current numerical relativity schemes forevolving both the spacetime and the fluid are explicit,there are strong limitations imposed on the size of thetimestep by the Courant-Friedrich-Lewy condition∆ t ∆ x = λ < C, (10)where λ is the Courant number and C a constant of orderunity that depends on the integration scheme employed.If one uses AMR, this implies that the size of the timestephas to be different in regions of different mesh size ∆ x .If we resolve the stars adequately, the mesh size will bemuch smaller near the NS than near the WD, becausein typical scenarios the NS is 500 times smaller than theWD. Eq. (10) then implies that the smallest timestepmust be in the domain of the NS. In particular, if theNS is covered by N NS = 2 R NS / ∆ x NS grid zones and theWD is covered by N WD = 2 R WD / ∆ x WD grid zones, thenfrom Eq. (10) we have∆ t WD ∆ t NS = ∆ x WD ∆ x NS = N NS N WD R WD R NS , (11)where ∆ t NS , ∆ x NS and ∆ t WD , ∆ x WD denote thetimestep and mesh size in the vicinity of the NS andWD, respectively.To assess the potential formation of a black hole re-quires at least N NS = 50 grid zones across the NS, inorder that a 2 M ⊙ BH (which is a probable mass a BHwould have after the merger of a WDNS system of to-tal mass of about 2 . M ⊙ ) would be covered by at least20 grid zones. Even if a BH does not form, coveringthe NS with 50 grid zones is necessary to reliably modelthe NS and maintain small hamiltonian and momentumconstraint violations. Resolving the WD requires about N WD = 30 grid zones to reliably model the star. If wecombine Eq. (7) and Eq. (11), we obtain∆ t WD ∆ t NS ≈ . (12)This means that for one timestep in the vicinity of theWD we would have to take about 833 timesteps in thevicinity of the NS. Even evolving the system for onlyone WD dynamical time scale would require millions oftimesteps in the vicinity of the NS. This shows how diffi-cult it is to resolve both the WD and the NS at the sametime.However, what renders the computation of the WDNSmerger in full GR prohibitive is that a realistic mergertakes place on an orbital time scale, which is equivalentto 10 NS dynamical time scales (see Eq. (9)).To make this quantitative, let us compare the orbitaltime scale with a typical timestep in the vicinity of theNS. Using N NS grid zones across the NS and combiningEq. (1) and Eq. (10) we find t d , NS ∆ t NS = N NS λ r R NS M NS ∼ , (13)where in the last step we used a typical value forthe Courant number λ = 0 . M NS , R NS , N NS we cited above. Combining Eq. (9) andEq. (13) we obtain T ∆ t NS ∼ . (14)Hence, a realistic WDNS simulation would require a min-imum of 10 timesteps in the vicinity of the NS. In fact FIG. 1: Mass (M) – central rest-mass density ( ρ ,c ) relation-ship of single TOV stars for various cold EOSs. In the plot ρ nuc = 2 × g/cm is the nuclear density. Plotted are theChandrasekhar electron-degenerate EOS for mean molecularweight per electron µ e = 2 (labeled as WD), the AP2 versionof the Akmal-Pandharipande-Ravenhall EOS [37, 60], a poly-tropic approximation of these realistic EOSs using EOS (15)(labeled as Fit) and a version of EOS (15) where the ratio ofthe isotropic radius of a 0 . M ⊙ pWD to the isotropic radiusof a 1 . M ⊙ NS is reduced to ≈
10 (labeled as 10:1). Theparameters of these EOSs are listed in Table I. this number of timesteps is an underestimate becauseextracting GWs would require a few orbits and the finalsystem would settle in equilibrium or collapse within afew orbital time scales after merger. As a result, a dy-namical, fully general relativistic hydrodynamics WDNScalculation would require of order 10 timesteps.We can give an estimate of how long such a simula-tion would be based on high resolution (192 grid pointsin the innermost refinement level) benchmark runs weperformed for a WDNS system with a 1 . M ⊙ -NS anda 1 . M ⊙ -WD at separation of about 2 . R WD (close tothe Roche limit), which has an orbital period of about1 . × M , where M = 2 . M ⊙ . Using 256 cores on theRanger cluster of the Texas Advanced Computing Cen-ter we found that the Illinois GR hydrodynamics codeadvances about 6 M per hour. Thus, the entire simula-tion (of about 10 orbital periods) would require about264 years of pure computational time.Realistic WDNS simulations are beyond the capabil-ities of current computational resources and numericalrelativity techniques. For this reason, we will tackle theproblem of WDNS mergers and head-on collisions adopt-ing an alternate strategy. We carefully construct an EOSwhich mimics a realistic cold NS EOS and, at the sametime, scales down the size of the WD to make such acalculation feasible. Using sequences of these systems,where the WD size gradually increases, we can extrapo-late our results to the realistic case. These scaled downWDs or pseudo-WDs are the subject of the following sec-tion. III. PSEUDO-WHITE DWARFS
In this section we introduce our EOS and describe re-sulting models for pWDs. Our EOS is the following 6-parameter piecewise polytropic EOS Pρ = κ ρ /n , ρ ρ κ ρ /n , ρ < ρ ρ ,κ ρ /n , ρ > ρ (15)where P is the pressure, ρ is the rest-mass densityand κ , κ , κ , n , n , n , ρ , ρ are the parameters of theEOS. Note that the parameters are 8 in number but con-tinuity requires that the following conditions be true κ = κ ρ /n − /n , κ = κ ρ /n − /n . (16)As a result, the adopted EOS (15) has 6 free parametersand throughout this paper we use κ , n , n , n , ρ , ρ tospecify an EOS. Note also that we use the polytropicindices n i instead of the adiabatic indices Γ i = 1 + 1 /n i .The freedom of our multi-parameter EOS enables usto capture the same characteristic curves and turningpoints on a TOV mass-central density plot as for a cold-degenerate realistic EOS (see [59]), as shown in Fig. 1.The figure shows that EOS (15) can provide a rea-sonable approximation to the mass-central density rela-tion of realistic compact objects, exhibiting both stable( dM/dρ ,c >
0) and unstable ( dM/dρ ,c <
0) branchesfor both WDs and NSs.Furthermore, EOS (15) allows us to adjust the size ofa pWD of given mass, thereby shifting the pWD branchto smaller radii (see Fig. 2), while keeping the NS massesand radii approximately unchanged for M NS > . M ⊙ .The shifted branches in Fig. 2 correspond to the starsthat we call pseudo-WDs in this work.Finally, note that all versions of EOS (15) consideredin this work have been carefully constructed so that themaximum gravitational mass of a NS is 1 . M ⊙ , i.e., thesame as that which for the AP2 version of the Akmal-Pandharipande-Ravenhall (APR) EOS [37, 60], and themaximum gravitational mass of a pWD is 1 . M ⊙ , i.e.,the maximum mass of a TOV WD obeying the Chan-drasekhar EOS for µ e = 2. In addition, EOS (15) isconstructed to preserve the shape of the M vs ρ ,c curve,yielding both stable and unstable branches and turningpoints appropriately. IV. INITIAL DATA
In this section we present the basic formalism for gener-ating valid general relativistic initial data for the head-oncollision in a given pWDNS system.
A. Gravitational Field Equations
We begin by writing the spacetime metric in the stan-dard 3+1 form [61] ds = − α dt + γ ij ( dx i + β i dt )( dx j + β j dt ) , (17)where α is the lapse function, β i the shift vector and γ ij the three-metric on spacelike hypersurfaces of constanttime t . Throughout the paper Latin indices run from 1to 3, whereas Greek indices run from 0 to 3.We conformally decompose the three-metric γ ij as γ ij = Ψ f ij , (18)where Ψ is the conformal factor and f ij the conformalmetric. We adopt the standard approximation of a con-formally flat spacetime, so that f ij = δ ij in Cartesiancoordinates.Since we are interested in head-on collisions betweencompact objects, we assume that initially the stars beginto accelerate towards each other starting from rest. As aresult the extrinsic curvature is initially zero and the mo-mentum constraints are identically satisfied [62]. Hence,we need only prepare initial data for Ψ.Under the aforementioned assumptions the only equa-tion we have to solve is the Hamiltonian constraint, whichbecomes ∇ Ψ = − π Ψ ρ, (19)where ∇ is the flat Laplacian operator. Here the sourceterm ρ is defined as ρ ≡ n α n β T αβ , (20)and where n α is the normal vector to a t = constant slice,and T αβ is the stress-energy tensor of the matter.The gauge is chosen so that the initial slice is maximal,i.e., K = 0 and ∂ t K = 0, and the shift vector is setequal to zero. Using the assumption of maximal slicingit is straightforward to derive an equation for the lapse[11, 63] ∇ ˜ α = 2 π ˜ α Ψ ( ρ + 2 S ) , (21)where ˜ α = α Ψ , (22)and the source term S is defined as S ≡ γ ij T ij . (23)Equations (19) and (21) are elliptic and hence have tobe supplemented with outer boundary conditions. Fol-lowing [63], we impose 1/r fall-off conditions on α − − TABLE I: Parameters for the piecewise polytropic EOS (15) used in generating different stellar models. The first columncorresponds to the name of the EOS. An EOS named M : N corresponds to a version of (15) for which the mass – radiusrelationship of TOV stars is such that the ratio of the isotropic radius of a 0 . M ⊙ pWD to that of a 1 . M ⊙ NS is M : N (seeFig. 2). The EOS named AP2 is the same as the AP2 EOS defined in [60]. Finally, the EOS named Fit is an approximate fitto the Chandrasekhar EOS (for µ e = 2) joined onto the AP2 EOS (see Fig. 1). Here ρ nuc = 1 . × − km − and κ isgiven in geometrized units.EOS Name κ n n n log( ρ /ρ nuc ) log( ρ /ρ nuc )20:1 5064.2599 1.56128 2.98418 0.714286 -3.17219 0.18047310:1 4993.0688 1.51515 2.96971 0.714286 -2.26862 0.2085025:1 6123.5567 1.51883 2.94291 0.699301 -1.2909 0.267623AP2 145414.043 0.60864 0.49652 0.514139 0.398915 0.698922Fit 4458.0491 2. 2.96736 0.716 -6.39356 0.208502 B. Matter fields
To model the matter, a perfect fluid stress-energy ten-sor is assumed: T αβ = ( ρ + ρ i + P ) u α u β + P g αβ , (24)where g αβ is the inverse of the four-metric and ρ , ρ i , P, u α are the rest-mass density, internal energydensity, pressure, and four-velocity of the fluid respec-tively.Since the initial configuration is assumed to be at restin the center of mass frame, the initial fluid four velocityis given by u α = u t (1 , , ,
0) (25)or u α = αu t n α . (26)A straightforward calculation shows that the sourceterm ρ in Eq. (20) can then be written as ρ = ρ + ρ i , (27)and S in Eq. (23) as S = 3 P. (28) C. Computational methods
To solve the elliptic equations (19) and (21) we de-veloped a fixed-mesh-refinement (FMR) finite differencecode based on the Portable, Extensible Toolkit for Scien-tific Computation (PETSc) algorithms [64–66]. The gridstructure used in our code is a multi-level set of properlynested, cell-centered uniform grids. We use a standardsecond-order finite difference stencil for the Laplacian op-erator and first order interpolation across the refinementlevel boundaries. The non-linearity of Eq. (19) is ad-dressed by performing Newton-Raphson iterations. Abrief description of our FMR implementation is summa-rized in Appendix A, to which we refer the interestedreader.
1. Diagnostics
To check the consistency of solutions obtained with ourFMR code we calculate the following diagnostic quanti-ties:The ADM mass is given by M ADM = − π I ∞ ∂ i Ψ dS i = − π Z ∇ Ψ d x, (29)where we have applied Gauss’ theorem to convert thesurface integral into a volume integral. The actual ex-pression we use to calculate the ADM mass volume inte-gral is (29) with ∇ Ψ replaced by the right-hand-side ofEq. (19).The Komar mass is given by M K = 14 π I ∞ ∂ i αdS i = 14 π Z ∇ αd x, (30)where again we have applied Gauss’ theorem in the laststep. By use of Eqs. (19) and (21) we find ∇ α = − ∇ i α ∇ i Ψ + 4 πα Ψ ( ρ + S ) . (31)The actual expression we use to calculate the Komarmass volume integral is (30) with ∇ α replaced by theright-hand-side of Eq. (31).Finally, the total baryon mass is given by [63] M = Z M ρ o αu t Ψ d x, (32)where M means that the integration is carried over thesupport of the matter.
2. Code Testing
Gauss’ theorem constitutes a strong consistency checkfor our FMR code. To demonstrate that the solutionsobtained with our elliptic code satisfy Gauss’ theoremand achieve second-order convergence we performed the
FIG. 2: Mass – radius relationship of TOV stars generatedby the piecewise polytropic EOS (15). Plotted are curves cor-responding to three versions of the EOS where the ratio ofthe isotropic radius of a 0 . M ⊙ pWD to that of a 1 . M ⊙ NS is 20:1, 10:1 and 5:1. The corresponding EOS parame-ters are given in Table I. The open points correspond to theNS models studied in this paper, which have masses 1 . M ⊙ ,1 . M ⊙ and 1 . M ⊙ . The solid points correspond to the pWDmodels considered in this work, which all have the same mass:0 . M ⊙ . following test. Employing the 10 : 1 piecewise poly-tropic EOS we constructed TOV NS solutions of vari-ous masses with a 1D TOV integrator. We then usedsecond-order polynomial interpolation to set up the rest-mass density profiles in our FMR elliptic code and solveEqs. (19) and (21) for the conformal factor and lapsefunction. We set up grids with five levels of refine-ment ( nl = 5) centered on the NS, and three differ-ent resolutions nx = ny = nz = 32 , ∆ x = 2 . nx = ny = nz = 64 , ∆ x = 1 . nx = ny = nz =128 , ∆ x = 0 . V. EVOLUTION OF WDNS SYSTEMSA. Basic Equations
The formulation and numerical scheme for our simula-tions are the same as those already reported in [22, 58,67], to which the reader may refer for details. Here weintroduce our notation and summarize our method.We use the 3+1 formulation of general relativity wherethe metric is decomposed in the form of Eq. (17), andwhere the fundamental dynamical variables for the met-ric evolution are the spatial three-metric γ ij and extrin-sic curvature K ij . We adopt the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism [11, 68, 69], inwhich the evolution variables are the conformal expo-nent φ ≡ ln( γ ) /
12, the conformal 3-metric ˜ γ ij = e − φ γ ij ,three auxiliary functions ˜Γ i ≡ − ˜ γ ij,j , the trace of theextrinsic curvature K , and the trace-free part of the con-formal extrinsic curvature ˜ A ij ≡ e − φ ( K ij − γ ij K/ γ = det( γ ij ). The full spacetime metric g µν is re-lated to the three-metric γ µν by γ µν = g µν + n µ n ν , wherethe future-directed, timelike unit vector n µ normal to thetime slice can be written in terms of the lapse α and shift β i as n µ = α − (1 , − β i ). The evolution equations of theseBSSN variables are given by Eqs. (9)–(13) in [22].We adopt standard puncture gauge conditions: an ad-vective “1+log” slicing condition for the lapse and a“Gamma-freezing” condition for the shift [70]. Thus, wehave ∂ α = − αK , (33) ∂ β i = (3 / B i , (34) ∂ B i = ∂ ˜Γ i − ηB i , (35)where ∂ ≡ ∂ t − β j ∂ j . We set the η parameter to0 . − for all simulations presented in this work.The fundamental matter variables are the rest-massdensity ρ , specific internal energy ǫ , pressure P , andfour-velocity u µ . We write the stress-energy tensor as T µν = ρ hu µ u ν + P g µν , (36)where h = 1+ ǫ + P/ρ is the specific enthalpy and ǫ is thetotal energy density. In our numerical implementation ofthe hydrodynamics equations, we evolve the “conserva-tive” variables ρ ∗ , ˜ S i , and ˜ τ . They are defined as ρ ∗ ≡ −√ γ ρ n µ u µ , (37)˜ S i ≡ −√ γ T µν n µ γ νi , (38)˜ τ ≡ √ γ T µν n µ n ν − ρ ∗ . (39)The evolution equations for these variables are given byEqs. (21)–(24) in [22].The EOS we adopt for the evolution has both a thermaland cold contribution, i.e., P = P th + P cold , (40) FIG. 3: Convergence test for the FMR elliptic code usingTOV stars. Five levels of refinement have been used for thistest. Upper panel: Fractional difference between the volumeintegral M vol and surface integral M surf for the ADM mass,calculated with our FMR elliptic code, versus the gravita-tional mass M TOV calculated with a 1D integrator of the TOVequations. The plot demonstrates satisfaction of Gauss’ the-orem and second-order convergence. Lower panel: Fractionaldifference between M vol and M TOV versus M TOV . The plotdemonstrates equality of M vol and M TOV and second-orderconvergence of our FMR elliptic code to the 1D (exact) re-sult. In both panels three resolutions are plotted: low ≡ ,med ≡ , high ≡ . Resolutions 32 and 64 have beenrescaled with a factor of 1 /
16 and 1 / . where P cold is given by Eq. (15) and the thermal pressureis given by P th = (Γ th − ρ ( ǫ − ǫ cold ) , (41)where ǫ cold = − Z P cold d (1 /ρ ) . (42)We set Γ th = 1 .
66 ( ≃ /
3) in all our simulations, i.e., weset it equal to the Γ exponent of the 10 : 1 EOS, appro-priate either for nonrelativistic cold degenerate electronsor (shock) heated, ideal nondegenerate baryons. Eq. (40)reduces to our piecewise polytropic law Eq. (15) for theinitial (cold) NS and pWD matter. B. Evolution of the metric and hydrodynamics
We evolve the BSSN equations using fourth-order ac-curate, centered finite-differencing stencils, except on shift advection terms, where fourth-order accurate up-wind stencils are applied. We apply Sommerfeld out-going wave boundary conditions on all BSSN fields, asin [22]. Our code is embedded in the Cactus paralleliza-tion framework [71], and our fourth-order Runge-Kuttatimestepping is managed by the
MoL (Method of Lines)thorn, with the CFL number set to 0.45 in all pWDNSsimulations. We use the Carpet [72] infrastructure toimplement the moving-box adaptive mesh refinement. Inall AMR simulations presented here, we use second-ordertemporal prolongation, coupled with fifth-order spatialprolongation, and impose equatorial symmetry to reducethe computational cost.We write the general relativistic hydrodynamics equa-tions in conservative form. They are evolved via a high-resolution shock-capturing (HRSC) technique [58, 67]that employs the piecewise parabolic (PPM) reconstruc-tion scheme [73], coupled to the Harten, Lax, and vanLeer (HLL) approximate Riemman solver [74]. Theadopted hydrodynamic scheme is second-order accurate.To stabilize our hydrodynamic scheme in regions wherethere is no matter, a tenuous atmosphere is maintainedon our grid, with a density floor ρ atm set to 10 − times the initial maximum density on our grid. The ini-tial atmospheric pressure P atm is set by using the coldEOS (15). Throughout the evolution, we impose limitson the pressure to prevent spurious heating and nega-tive values of the internal energy ǫ . Specifically, we re-quire P min ≤ P ≤ P max , where P max = 10 P cold and P min = 0 . P cold , where P cold is the pressure calculatedusing the cold EOS (15). Whenever P exceeds P max ordrops below P min , we reset P to P max or P min , respec-tively. Following [75] we impose the upper pressure limitsonly in regions where the rest-mass density remains verylow ( ρ < ρ atm ), but we impose the lower limit every-where on our grid.At each timestep, the “primitive variables” ρ , P , and v i must be recovered from the “conservative” variables ρ ∗ , ˜ τ , and ˜ S i . We perform the inversion numerically asspecified in [58]. We use the same technique as in [22]to ensure that the values of ˜ S i and ˜ τ yield physicallyvalid primitive variables, except we reset ˜ τ to 10 − ˜ τ , max (where ˜ τ , max is the maximum value of ˜ τ initially) wheneither ˜ S i or ˜ τ is unphysical [i.e., violate one of the inequal-ities (34) or (35) in [22]]. The restrictions are usuallyimposed only in the low-density atmosphere.It is instructive to discuss the mathematical structureof the system of hydrodynamic equations when a piece-wise polytropic EOS (15) is used. According to [76, 77],when the fluxes of the conservation laws are non-smooth,split waves and composite structures may be present inthe solutions. In these cases a numerical solution maynot converge to the correct solution. In our case thefluxes are not smooth everywhere because EOS (15) isnon-smooth (it is continuous but not differntiable) at theturning points ρ i i = 1 ,
2. Away from the turning pointsthe fluxes are smooth, therefore there may be some con-cern when non-linear waves cross the transition densities.However, according to [77] our adopted numerical schemeshould be able to handle such composite structures, ifthey ever arise, and our solutions should converge to thecorrect continuum solution.To study this effect we constructed an EOS which issimilar to Eq. (15), but where the pressure discontinuitiesare “smoothed out” at the turning points ρ i , ( i = 1 , ρ i (1 − ǫ ) , ρ i (1 + ǫ )], where ǫ >
0. We chose ǫ to be sufficiently small so that thesmoothed EOS mimics as closely as possible EOS (15),but large enough to avoid roundoff errors due to verylarge gradients. Setting up several generalized Riemmanproblems, we found that the numerical solutions obtainedusing EOS (15) converge to those obtained when usingits smooth counterpart and the two can hardly be dis-tinguished for the resolutions considered. Therefore, ournumerical schemes in conjunction with EOS (15) are ableto capture the correct solution, in that they are almostidentical to the solutions obtained with the smooth coun-terpart of (15). For the details of this analysis, we referthe interested reader to Appendix B. C. Evolution Diagnostics
During the evolution, we monitor the normalizedHamiltonian and momentum constraints as defined inEqs. (40)–(43) of [22].We also monitor the ADM mass and angular mo-mentum of the system, which can be calculated dur-ing the evolution by surface integrals at a large distance(Eqs. (37) and (39) of [22]). The equations used tocalculate the ADM mass and angular momentum withminimal numerical noise are as follows [11] M = Z V d x (cid:18) ψ ρ + 116 π ψ ˜ A ij ˜ A ij − π ˜Γ ijk ˜Γ jik + 1 − ψ π ˜ R − π ψ K (cid:19) , (43) J i = 18 π ǫ ijn Z V d x (cid:2) ψ ( ˜ A jn + 23 x j ∂ n K − x j ˜ A km ∂ n ˜ γ km ) + 8 πx j S n (cid:3) . (44)Here V is the volume within a distant surface, ψ = e φ , ρ = n µ n ν T µν , S i = − n µ γ iν T µν , ˜ R is the Ricci scalarassociated with ˜ γ ij , and ˜Γ ijk are Christoffel symbols as-sociated with ˜ γ ij .In this work we only focus on head-on collisions, sothere is no angular momentum involved. However, oursimulations are three-dimensional, so there is no guaran-tee that J i will remain 0. In order to quantify violations of J i = 0 we normalize the angular momentum, com-puted via (44), with the angular momentum a pWDNSsystem would have, if the binary components were New-tonian point masses in circular orbit at the initial sepa-ration A , i.e., J z , c = M / A / q (1 + q ) , (45)where the total mass M T is taken to be the sum of ADMmasses of the isolated stars.When hydrodynamic matter is evolved on a fixed uni-form grid, our hydrodynamic scheme guarantees that therest mass M is conserved to machine roundoff error.This strict conservation is no longer maintained in anAMR grid, where spatial and temporal prolongation isperformed at the refinement boundaries. Hence, we alsomonitor the rest mass M = Z ρ ∗ d x (46)during the evolution. Rest-mass conservation is also vio-lated whenever ρ is reset to the atmosphere value. Thisusually happens only in the very low-density atmosphere.The low-density regions do not affect rest-mass conser-vation significantly.Shocks occur when stars collide. We measure theentropy generated by shocks via the quantity K ≡ P/P cold ≥
1, where P cold is the pressure associated withthe cold EOS that characterizes the initial matter (seeEq. (15)). VI. CASES AND RESULTSA. Initial configurations
We perform a number of pWDNS head-on collisionsimulations varying the initial configurations, so that wecan study the effect of the pWD compaction and NS masson the final outcome separately. Table II outlines thephysical parameters for the cases considered in this work,and Table III presents the AMR grid structure used ineach case.To generate initial data for our cases we first choose theADM mass M NS of the NS and the ADM mass M WD ofthe pWD (see Table II) and solve the TOV equations foreach star in isotropic coordinates to prepare the rest-massdensity distribution for the NS and the pWD separately.We then use second-order polynomial interpolation to in-terpolate the rest-mass density profiles onto the nestedgrids of our FMR elliptic initial value code and solveEqs. (19) and (21) for Ψ and α . The two stars are placedat coordinate separation A and such that the Newtoniancenter of mass of the system is identified with the originof the coordinate system. Once a solution is achieved bythe FMR code for the initial metric, we map ρ , Ψ and α from the elliptic code grids onto the evolution grids using0 TABLE II: Summary of initial configurations. M NS ( M WD ) stands for the ADM mass of an isolated NS (pWD) ( a ) , R NS ( R WD )is the isotropic radius of an isolated NS (pWD), C NS ( C WD ) is the compaction of an isolated NS (pWD), where the compactionis defined as the ratio of the ADM mass of the isolated star to its areal radius. M ADM is the ADM mass of the system and A the initial binary separation in isotropic coordinates. All cases have exactly the same coordinate separation of 586 . M NS /M ⊙ M WD /M ⊙ C NS C WD R WD /R NS R WD /M ADM M ADM /M ⊙ A/R WD A1 1.4 0.98 0.111 0.010 8.88 41.18 2.413 4.000A2 1.5 0.98 0.130 0.010 9.96 39.36 2.524 4.000A3 1.6 0.98 0.151 0.010 11.15 37.46 2.652 4.000B 1.5 0.98 0.130 0.019 4.99 19.76 2.524 7.967C 1.5 0.98 0.130 0.005 20.01 79.08 2.524 1.991 ( a ) Here we list the ADM masses, isotropic radii and compactions of the isolated (TOV) stars whose rest-mass density profileswe used to generate initial data for Ψ and α for a given case.TABLE III: Grid configurations used in our simulations. Here M is the sum of the ADM masses of the isolated stars, Max. res.is the grid spacing in the innermost refinement box surrounding the NS, N NS denotes the number of grid points covering thediameter of the NS initially, and N WD denotes the number of grid points covering the diameter of the pWD initially. Thesmallest outer boundary distance corresponds to case A3 and is 1028 M .Case M/M ⊙ Grid Hierarchy (in units of M ) ( a ) Max. res. N NS N WD A1 2.38 (534.33, 267.16, 133.58, 66.79, 35.78[N/A], 19.08[N/A], 10.44[N/A], 7.156[N/A]) M/ .
71 63 35A2a 2.48 (510.62, 255.31, 127.65 , 63.83, 34.81[N/A], 18.86[N/A], 10.15[N/A], 6.890[N/A]) M/ .
52 44 28A2b 2.48 (534.33, 267.16, 133.58, 66.79, 35.78[N/A], 19.08[N/A], 10.44[N/A], 7.156[N/A]) M/ .
71 56 35A3 2.58 (467.27, 233.64,116.82, 58.41, 29.20 [N/A], 15.58[N/A], 8.518[N/A], 5.841[N/A]) M/ .
22 56 38B 2.48 (534.33, 267.16, 133.58, 66.79, 35.78[31.93], 19.08[N/A], 10.44[N/A], 7.156[N/A]) M/ .
71 56 35C 2.48 (534.33, 267.16, 133.58, 66.79[N/A], 35.78[N/A], 19.08[N/A], 10.44[N/A], 7.156[N/A]) M/ .
71 56 35 ( a ) There are two sets of nested refinement boxes: one centered on the NS and one on the pWD. This column specifies thehalf side length of the refinement boxes centered on both the NS and pWD. When the side length around the pWD isdifferent, we specify the pWD half side length in square brackets. When there is no corresponding pWD refinement box (asthe pWD is much larger than the NS), we write [N/A] for that box. second-order polynomial interpolation. We always makesure that the resolution of the initial data grids is higherthan the resolution of the evolution grids. The surfaces ofthe stars are a locus of rapidly decreasing density gradi-ents. As a result, small oscillations due to interpolationmay arise and lead to negative rest-mass density. Tocure this, we set the density ρ equal to the tenuous at-mosphere density that we maintain on the grid whenever | ρ /ρ i,c | < − , where ρ is the value of the density af-ter the interpolation and ρ i,c , i = WD , NS , is the centraldensity of the WD or NS. We do not find such oscillationswhen interpolating the gravitational fields, which is mostlikely due to the fact that these are sufficiently smooth.In all cases we require that the sum of the ADM massesof the isolated stars be larger than the maximum gravi-tational mass of 1 . M ⊙ that our cold EOS can support(see Fig. 2). There exist at least 18 observed WDNS sys-tems that satisfy this requirement. Since typical NS andWD masses in massive WDNS binaries lie in the range1 . M ⊙ − . M ⊙ and 0 . M ⊙ − . M ⊙ respectively (seeSec. I), we choose pWD rest-mass density profiles that correspond to an ADM mass of 0 . M ⊙ in isolation andkeep it fixed in all cases we study. This almost fixesthe pWD rest mass, because of the small compaction( < .
02) of the pWDs we consider. The pWD rest-massvariation from case B to case C, due to fixing the ADMmass, is 0 . . M ⊙ is that the ratio of the isotropicradius of a pWD of such mass to the isotropic radius ofa 1 . M ⊙ NS is ≈
10 for the 10:1 EOS.Another quantity that we fix in all our simulations isthe initial coordinate separation A of the two compo-nents. This almost fixes the kinetic energy of the starswhen they collide. In particular, we set A = 4 R WD , A .Here R WD , A denotes the isotropic radius of the spheri-cal 0 . M ⊙ pWD used in cases A1, A2, A3. We choosethe initial separation this way because we want the starsto be sufficiently far apart so that spherical TOV ini-tial models remain in near equilibrium, and at the sametime, simulate the collision within reasonable time scales,as the head-on collision time scale varies as ∼ A / (see1 FIG. 4: Snapshots of rest-mass density profiles at selected times for case A2. The contours represent the rest-mass densityin the orbital plane, plotted according to ρ = ρ , max − . j − . ( j = 0 , , . . . , ρ , max = 4 . ρ nuc . The last two snapshots are nearthe end of the simulation, and they demonstrate that the density contours within a radius of about 150km remain unchanged.Here M = 2 . M ⊙ = 3 . . × − s is the sum of the ADM masses of the isolated stars. Eq. (4)).If A = 4 R WD , A , the NS tidal field in the vicinity ofthe pWD is small, validating our assumption that anequilibrium pWD is nearly spherical. To see this, let uscalculate the ratio of the tidal force of the NS on thesurface of the pWD to the surface gravity of the pWD F tNS F WD = M NS M WD (cid:18) R WD , A A (cid:19) ≃ , (47)where we used M NS = 1 . M ⊙ , M WD = 1 M ⊙ , A =3 R WD , A . Hence, any deviations from sphericity shouldbe small. The assumption of sphericity for case B iseven better because in this case A ≈ R WD , B , but worsefor case C, where A ≈ R WD , C . In principle, we couldincrease the separation so that the sphericity approxima-tion becomes better for all cases, but if the final remnantdoes not collapse promptly to form a BH starting at closeseparations, it is unlikely that it will collapse if the initialseparation is larger. This is due to the fact that for largerseparations the kinetic energy at collision will be larger, generating more shock heating that will work to preventprompt collapse.To summarize, the set of cases A1, A2 and A3 probethe effect of the NS mass on the final outcome, whereasthe set of cases A2, B and C probe the effect of the pWDcompaction on the final outcome. In the following sec-tions we summarize the results of our simulations. B. Dynamics of collision and effects of the NS mass
Here we describe the effects of the NS mass on thedynamics of pWDNS head-on collisions. We find thatabout 18% of the initial total mass escapes to infinity forall cases A1, A2 and A3. Nevertheless, the initial totalmass in these cases is large enough to guarantee that thefinal total mass of the pWDNS remnant still exceeds themaximum mass that our cold EOS can support. How-ever, prompt collapse to a black hole does not take placein any of the cases studied because strong shock heat-2
FIG. 5: Normalized angular momentum vs time. J z and J z , c are given by Eqs. (44) and (45), respectively. The solid,dashed and dotted lines correspond to cases A2, B and C,respectively. Here M = 2 . M ⊙ = 3 . . × − s. ing gives rise to a hot remnant. Ultimate collapse to aBH is almost certain after the remnant has cooled. Theoutcome of the three cases A1, A2, A3 is a TZlO.Overall, cases A1, A2 and A3 are qualitatively similarand for this reason we mainly describe case A2 as repre-sentative of this class of our simulations. Furthermore,our study of case A2 with two different resolutions (seeTable III) shows the results to be qualitatively insensitiveto resolution indicating that the resolutions used in oursimulations are sufficiently high. In what follows all caseA2 plots correspond to the high resolution run of caseA2, i.e., case A2b.In general, the head-on collision of pWDNS systemscan be decomposed into three phases: the acceleration,the plunge, and the quasiequilibrium phase.During the acceleration phase, the two stars acceleratetoward one another starting from rest. The separationdecreases as a function of time and this phase ends whenthe two stars first make contact.As the separation decreases, the increasing NS tidalfield strongly distorts the pWD. This can be seen in theequatorial rest-mass density contours of Fig. 4. In theinsets of Fig. 4, the NS interior is almost unchanged dur-ing this phase. In reality, it oscillates but is not tidallydistorted by the pWD. Nevertheless, the NS atmospheredoes expand. The insets also show that due to numericalerrors the NS veers slightly off the x-axis, which is thecollision axis in our simulations. In general, in all oursimulations both the NS and the pWD wiggle around the x-axis. The amplitude of the NS wiggling motion isat most 1% of the pWD radius, while the amplitude ofthe pWD off-axis motion is less than 0.1% of its radius.Hence, the collision is practically head-on.It is likely that this lack of symmetry is due to smallasymmetries introduced when mapping the initial dataonto the evolution grids via 2nd order interpolation. Thisis because 2nd order interpolation requires the use of 3grid points (per direction) that surround the point towhich one interpolates. However, the effect is small andour results cannot change qualitatively due to this smallasymmetry.Along with this small off-axis motion, the pWDNS sys-tem acquires a small amount of spurious angular momen-tum. Fig. 5 shows the normalized z component of theangular momentum of the system and demonstrates thatit is always less than 3% (see Eq. (45)). In addition toconserving the angular momentum to within 3%, the nor-malized Hamiltonian constraint violations remain smallerthan 1% and the normalized momentum constraint vio-lations smaller than 3%. These results hold for all casesstudied in this work.During the plunge phase the NS penetrates the pWD,launching strong shocks that sweep through and heatthe interior of the pWD. The NS outermost layers arestripped when they encounter the dense central parts ofthe pWD, and the NS is compressed. We find that atmaximum compression in case A2, the NS central densityonly increases by about 8% of the initial central density.Eventually, strong shocks sweep through the entirepWD interior and then transfer linear momentum to thepWD outer layers, a large fraction of which receives suf-ficient momentum to escape to infinity. This can be seenin Fig. 6, where a snapshot is shown of the total energyper unit mass U = − u − U >
0) mat-ter covers most of the computational domain, as shownin Fig. 6. The rest-mass density of the ejected materialis of order 10 − ρ nuc , but the total mass that escapes toinfinity is large. This is demonstrated in Fig. 7, whichshows the fractional change in the rest mass as a functionof time. We find the amount of matter that escapes incases A1, A2 and A3 is ≃
18% of the initial rest masswhen we extrapolate our results to late times.The thermal energy deposited into the ejected materialis significant, with K = P/P cold >
40. As the ejectedmatter comes from the WD outermost layers, its den-sity is very low. This implies that its initial pre-shockedsound speed is small. As a result, the Mach number ofthe ejected material can be very large prior to shock heat-ing, and so shock heating is very strong, i.e., K increasesfrom 1 initially to greater than 40 (see also discussion inAppendix B of [23]).The time scale for shock heating to equilibrate mustbe of order few times the dynamical (free-fall) time scaleof the WD (see Eq. (2)), as this is the only relevant3
FIG. 6: Snapshot of total energy per unit mass U = − u − U <
U > U = U min . j ( j = 0 , , . . . , U < U min = 10 − . The size of the boundmatter area is insensitive to the choice of U min . Here M = 2 . M ⊙ = 3 . . × − s. timescale. For cases A, the WD dynamical timescaleis roughly 400M. Our computations show that it actu-ally takes about 800M-1000M for the star to equilibrate,which is consistent with the estimate above.Material that did not receive sufficient thrust to escapeto infinity starts to rain down onto the NS and pWDmatter and the plunge phase ends when this process isover.During the quasiequilibrium phase the remnant settlesinto a spherical quasiequilibrium object whose outer lay-ers oscillate. This can be seen in the two final snapshotsof Fig. 4, where we show that the inner equatorial rest-mass density contours do not change with time, whilethe outer layers change only a little. Fig. 8 plots xy, xzand yz rest-mass density contours. Notice that in theadopted gauge, the remnant appears to be spherical.The pWDNS final remnant consists of a cold NS corewith a hot mantle on top. This is demonstrated bythe plots in the last row of Fig. 8, where we quan-tify the results of shock heating by showing contours of K = P/P cold . Within a radius of 100 km from the centerof mass of the remnant, K ranges from unity to about15 for case A2. In all cases K ≃ . M ⊙ our cold EOS can support. In Fig. 9 we show the rest mass of the remnant as a function oftime and for various spatial domains. This figure demon-strates that the rest mass within a radius of 220km ac-counts for more than 90% of the final total rest mass andis greater than 1 . M ⊙ . However, the pWDNS remnantdoes not collapse promptly to form a black hole, becauseof extra support provided thermal pressure. Delayed col-lapse to a black hole is almost certain after the pWDNSremnant has cooled.Finally, we note that it has been suggested in [78, 79]that GWs may arise from shocks. Even though the dis-cussion in these studies focused on core collapse super-novae, the appearance of strong shocks in our case canalso generate GWs. However, here we do not calculatethe GW signature because what is really interesting froman observational and astrophysical point of view is theGW signature in the circular binary WDNS case, notthe head-on case we consider. In [33] we did calculateGWs from the inspiral phase of binary WDNS systems.General relativistic computations of the merger of circu-lar binary WDNSs will be the subject of a forthcomingpaper. C. Effects of the pWD compaction
Here we describe the effects of the pWD compactionon the dynamics of pWDNS head-on collisions. Overall,4
TABLE IV: Summary of pWD compaction study. Here C WD is the compaction of an isolated pWD (see Table II), K = P/P cold at the end of simulations ( a ) , T p is the peak temperature at the end of simulations, M (0) is the initial total rest mass,∆ M = | M , f − M (0) | , where M , f is the final total rest mass, M , r < is the mass enclosed within 220 km from the center ofmass of the remnant at the end of the simulations, ρ , c is the final value of the central rest-mass density ( b ) , and α min the finalvalue of the minimum of the lapse function.Case C WD K T p (10 o K) ( c ) ∆ M /M (0) M , r < /M ⊙ ρ , c /ρ nuc α min B 0.0190 [1,35] 3.7 14% 2.180 4.91 0.570A2 0.0100 [1,15] 3.2 18% 2.035 4.49 0.595C 0.0049 [1,10] 3.0 18% 1.900 4.10 0.609 ( a ) The K column lists the range of values which K obtains within a radius of 100 km from the centers of mass of theremnants. ( b ) ρ nuc = 2 × g/cm . ( c ) For realistic WDNS collisions we expect T p ∼ o K (see discussion following Eq. (51)).FIG. 7: Fraction of rest mass lost vs time. Here ∆ M = | M − M (0) | , where M (0) is the initial total rest mass.Small changes in the rest mass until approximately 5000 M for cases A1, A2, A3, B and 10000 M for case C are due tointerpolations when matter crosses refinement levels and in-accuracies in evolving the very low-density atmosphere. Atthe end of the simulations the amount of mass ejected is13.4% of the initial rest mass in case B and ranges from16.1%-16.7% of the initial rest mass for the other cases.Extrapolating the results to late t we find that in case B∆ M /M (0) ∼ M /M (0) ∼
18% in all other cases.Here M = 2 . M ⊙ = 3 . . × − s. our findings are qualitatively similar to those of case A2described in the previous section. An appreciable frac-tion of the initial total mass escapes to infinity, but thefinal total mass of the pWDNS remnant still exceeds themaximum mass that our cold EOS can support. Promptcollapse to a black hole does not take place either in case B or in case C, because strong shock heating gives rise toa hot remnant. The outcome of cases B and C is again aTZlO.The three phases of the head-on collision we describedin the previous section apply here, too. For this reasonwe now focus our discussion on describing the differencesbetween cases B, A2 and C, i.e., in order of decreasingpWD compaction.The tidal acceleration, which the pWD experiences dueto the NS tidal field, increases as the pWD compactiondecreases. This is because the initial coordinate separa-tion is the same for cases B, A2 and C. As a result, theacceleration phase is shorter for larger pWD compaction.This can be seen in Figs. 10 and 11, where equatorialrest-mass density contours are plotted.Shock heating far from the core of the remnants, asquantified by K , is somewhat less intense as the pWDcompaction decreases. The shorter acceleration phaseimplies that the relative speed of the two componentsat the beginning of the plunge phase is a little smaller,which in turn leads to weaker shocks.During the plunge phase, the NS interior is less affectedby decreasing pWD compaction. This can be seen (a) inthe insets in Figs. 10 and 11, where in case C the post-plunge structure of the NS core is almost the same as thatshowed in the pre-plunge snapshots, while this is not truefor case B, and (b) by the variation in the NS centraldensity ( ρ , c ); in particular, we find that at maximumcompression the NS central density increases by about42% in case B, 8% in case A2 and 5% in case C. Theseresults can be easily interpreted because in a sequenceof pWDNS head-on collisions where the NS is fixed andthe size of the pWD increases with fixed mass, the NSgradually encounters thinner and thinner material, andhence changes to the NS structure become less and lessimportant.Were the system mass loss to vary appreciably withpWD size, we might expect a corresponding variation in ρ , c . However, such a mass loss variation is not observed,as we discuss next.As in case A2, in both cases B and C, a large frac-tion of the initial mass eventually escapes to infinity (see5 FIG. 8: First three rows: Snapshots of rest-mass density profiles at the end of the simulation for cases A2, B, C. The contoursrepresent the rest-mass density in the YZ plane (first row), in the XZ plane (second row) and in the XY plane (third row)plotted according to ρ = ρ , max − . j − . ( j = 0 , , . . . , ρ , max = 4 . ρ nuc . Thesesnapshots demonstrate that in the adopted gauge, the final object is roughly spherical. Last row: Snapshots of K = P/P cold profiles at the end of the simulation for cases A2, B, C. The contours represent K in the XY plane plotted according to K = 10 . j ( j = 0 , , . . . , K ≃ K becomes larger than unityas we move outwards from the center of the objects, and shock heating is more intense in case B and less intense in case C.The color code used is the same as that defined in Fig. 4. Here M = 2 . M ⊙ = 3 .
662 km = 1 . × − s. FIG. 9: Post-merger rest mass as a function of time. Here M , tot is the total rest mass in the entire computational domain and M , r < r stands for the rest mass contained within a coordinate sphere of radius r in units of km, centered on the remnant’scenter of mass. In all cases M , r < accounts for more than 90% of the final total rest mass and it is always greater than1 . M ⊙ - the maximum rest mass that our cold EOS can support. Here M = 2 . M ⊙ = 3 . . × − s. Fig. 6 for case B), but we do not find strong variationsin the mass lost among cases B, A2 and C. We find thatthe amount of matter that escapes in case B is 14% ofthe initial rest mass, while case C loses 18%, (case A2loses 18%) (see Fig. 7) of the initial rest mass. Givenour earlier discussion that shocks in case B are strongerthan those in case A2, and in turn shocks in case A2 arestronger than those in case C, this last result may soundcontradictory, because one might expect that strongershocks would eject more matter to infinity. The appar-ent contradiction can be resolved, if one considers thatas the pWD compaction decreases, the pWD outer layersbecome less and less bound, and hence, it requires lessenergy to eject them to infinity.As in case A2, the remnants in cases B and C even-tually settle into spherical quasiequilibrium objects withoscillating outer layers (see Figs. 10, 11). The sphericityof the remnants (in the adopted gauge) in cases B and Cis demonstrated by the xy, xz and yz rest-mass densitycontours shown in Fig. 8.The pWDNS remnants in both case B and case C con-sist of a cold NS core with a hot mantle on top. Thus,all cases lead to the formation of a TZlO. This is againdemonstrated in the last row of Fig. 8, where contours of K = P/P cold are shown. Within a radius of about 100 kmfrom the center of mass of the final remnants, K ∈ [1 , K ∈ [1 ,
10] in case C ( K ∈ [1 ,
15] in caseA2). Thus, shock heating is overall strongest in case B,weaker in case A2 and even weaker in case C.In Fig. 9 the remnant rest masses in cases B and C areplotted as functions of time. The figure demonstratesthat in both cases, the rest masses within a radius of220km account for more than 90% of the final total restmasses in both case B and case C, and are greater than1 . M ⊙ , i.e., the maximum rest mass that our cold EOScan support. The pWDNS remnant does not collapsepromptly to form a black hole, because of the extra ther- mal pressure support. However, delayed collapse to ablack hole is almost certain after the pWDNS remnanthas cooled.Another important feature of Fig. 9 is that the amountof mass contained within a given radius from the centerof mass of the remnant is larger for smaller initial pWDs.For example, within a radius of 220km the remnant massis 2 . M ⊙ in case B, 2 . M ⊙ in case A2, and 1 . M ⊙ in case C. This in turn indicates that the higher the ini-tial pWD compaction the higher the core densities of thepWDNS remnant. This is supported by the rest-massdensity contours shown in Fig. 8 and by the values of thefinal central rest-mass density. In particular, we find thatthe final central rest-mass density is 4 . ρ nuc in case C,4 . ρ nuc in case A2, and 4 . ρ nuc in case B. Thus, thereis a variation in the final central density of 9 .
2% fromcase B to case A2, and 9 .
5% from case A2 to case C.Finally, it is also worth noting that the final minimumvalue of the lapse function, which is a good indicator ofcollapse, increases with increasing initial pWD size, too.Specifically, we find this value to be 0 .
57 in case B, 0 . .
609 in case C. All these facts seem toindicate that as the pWD size increases towards realisticWD sizes the less likely it is for the pWDNS remnant tocollapse. To demonstrate this trend more clearly we com-pile all aforementioned physical parameters of the finalconfigurations in cases B, A2, and C in Table IV.Hence, given the consistency in the results of casesB, A2 and C, i.e., the sequence of increasing pWD sizewith fixed pWD mass, we expect that as the parametersof our EOS vary, so as to describe realistic WDs, theresult of the head-on collision of a massive WDNS systemwill most likely lead to formation of a quasiequilibriumTZlO. If the initial total mass of the system exceeds themaximum mass a cold EOS can support, then the TZlOmay have mass that exceeds the maximum mass a coldEOS can support, so it will eventually undergo collapse7
FIG. 10: Snapshots of rest-mass density contours at selected times for case B. The contours are plotted in the orbital plane,according to ρ = ρ , max − . j − . ( j = 0 , , . . . , ρ , max = 4 . ρ nuc . The last two snapshots, which take place near the end of the simulation, demonstratethat the density contours within a radius of about 150km remain unchanged. Here M = 2 . M ⊙ = 3 . . × − sis the sum of the isolated stars’ ADM masses. to a black hole, but only after the remnant has cooled.To identify the dominant cooling mechanisms and/orrelevant nuclear reaction networks, one would need toestimate the temperatures of these TZlOs. We can dothis as follows.Using Eq. (40) and the definition of K we can calculatethe specific thermal energy as ǫ th = ( K − P cold (Γ th − ρ . (48)To estimate the temperature T of matter, we assumethat we can model the temperature dependence of ǫ th as ǫ th = 3 kT m n + f aT ρ , (49)where m n is the mass of a nucleon, k is the Boltz-mann constant and a is the radiation constant. Thefirst term represents the approximate thermal energyof the nucleons, and the second term accounts for thethermal energy due to relativistic particles. The fac-tor f reflects the number of species of ultrarelativistic particles that contribute to the thermal energy. When T << m e /k ∼ K, where m e is the electron mass,thermal radiation is dominated by photons and f = 1.When T >> m e /k , electrons and positrons becomeultrarelativistic and also contribute to radiation, and f = 1 + 2 × (7 /
8) = 11 /
4. At sufficiently high tem-peratures and densities ( T > K , ρ > g cm − ),neutrinos are generated copiously and become trapped,so, taking into account three flavors of neutrinos and an-tineutrinos, f = 11 / × (7 /
8) = 8.In Fig. 12 we show the temperature profiles of the rem-nants in cases B, A2 and C, where it is clear that typicaltemperatures of our TZlOs are of order 10 o K. This isexpected as the total energy available for shock heatingshould be of order M NS M WD /R WD , i.e., the gravitationalinteraction energy when the two stars first make contact.The total thermal energy, E th , is then E th ∼ ( M NS + M WD ) m n kT ∼ M NS M WD R WD . (50)From this last equation we can estimate the characteristic8 FIG. 11: Rest-mass density contours in the orbital plane at selected times for case C. The contours are plotted according to ρ = ρ , max − . j − . ( j = 0 , , . . . , ρ , max = 4 . ρ nuc . The last two snapshots, which take place near the end of the simulation, demonstrate that thedensity contours within a radius of about 150km remain unchanged. Here M = 2 . M ⊙ = 3 . . × − s is the sumof the isolated stars’ ADM masses. temperature as T ∼ C WD m n (1 + q ) k . (51)Thus, all things being equal (no mass loss, same kineticenergy at collision, etc.) characteristic TZlO tempera-tures should be directly proportional to the pWD com-paction. Taking case A2 as an example, C ≃ .
01 and q ≃ /
3, we find T ≃ . × o K, in rough agreementwith our simulations.Using this scaling argument we can extrapolate ourresults to realistic WDNS head-on collisions. For asolar-mass WD obeying the Chandrasekhar EOS C WD ≃ × − . Hence, we predict that typical temperaturesin realistic head-on collisions of massive WDNS systemswould be of order 10 o K. VII. SUMMARY AND DISCUSSION
In this work we studied the dynamics of the head-oncollision of WDNS binaries in full general relativity, aidedby simulations that employ the Illinois AMR relativistichydrodynamics code [23, 58]. This study serves as a pre-lude to the circular binary WDNS problem which will bethe subject of a forthcoming paper.Our primary focus is on compact objects whose totalmass exceeds the maximum mass that a cold EOS cansupport and our goal is to determine whether a WDNScollision leads either to prompt collapse to a black holeor the formation of a quasi-equilibrium configuration thatresembles a Thorne-Zytkow object (TZO) [57], which wecall Thorne-Zytkow-like object (TZlO). By a TZlO wemean a NS sitting at the center of a hot gaseous mantlecomposed of WD debris.Due to the vast range of dynamical time scales andlength scales involved in this problem, realistic WDNSsimulations (head-on or otherwise) are computationallyprohibitive, if one employs current numerical relativity9
FIG. 12: Temperature ( T ) profiles for cases B, A2 and C. Thetemperature is in units of 10 o K. The profiles correspondto the values of T at the end of the simulations and along thex-axis (for x > x c , where x c is the x position of the center ofmass of the remnant in each case). It is clear that typical tem-peratures are of order 10 o K. For realistic WDNS collisionswe expect T ∼ o K (see discussion following Eq. (51)). techniques and available computational resources. Forthis reason, we tackle the problem using a different ap-proach.In particular, we constructed a piecewise polytropicEOS which captures the main physical features of NSsand, at the same time, scales down the size of a WD.We call these scaled-down models of WDs pseudo-WDs(pWDs). Using these pWDs, we can reduce the rangeof length scales and time scales involved, rendering thecomputations feasible.A pWD is not a realistic model of a WD. However,with our proposed parametrized EOS we can construct asequence of pWD models with gradually increasing sizeand perform simulations that approach the realistic case.Then we can make predictions about the realistic caseby extrapolating the results from this sequence of simu-lations.Using pWDs, we performed two sets of simulations.One set of our simulations studied the effects of the NSmass on the final outcome, when the pWD is kept fixed ata mass of 0 . M ⊙ and its size fixed at 146km. We choosethree masses for the NS, namely 1 . M ⊙ , . M ⊙ , . M ⊙ (cases A1, A2 and A3, respectively). The other set ofsimulations studied the effect of the pWD compaction onthe final outcome, when the NS is kept fixed at a mass of1 . M ⊙ . In the latter set of calculations, we choose threevalues for the ratio of the pWD to the NS radius, namely 5 : 1, 10 : 1 and 20 : 1 (cases B, A2 and C, respectively).In general, the head-on collision of pWDNS systemscan be decomposed in three phases: i) acceleration, ii)plunge and iii) quasiequilibrium.During the acceleration phase the two stars acceleratetowards one another starting from rest. As the separationdecreases the NS tidal field becomes so strong that thepWD becomes highly distorted, while the NS interior isalmost unchanged. This phase ends when the NS andpWD first make contact.During the plunge phase the NS penetrates the pWD,launching strong shocks that sweep through and heatthe interior of the pWD. The NS outermost layers arestripped after encountering the dense central parts of thepWD (see Fig. 4), but the NS core is mostly unaffected,except when the compaction of the pWD is high (seeSec. VI C). We find that the strong shocks sweepingthe pWD transfer linear momentum to the outer pWDlayers, causing a large amount of pWD matter to escapeto infinity. In all calculations, we find the rest mass lossto be between 14% −
18% of the initial total rest mass.Material that did not escape to infinity accretes onto theunderlying NS and pWD matter.Finally, during the quasiequilibrium phase, the rem-nant settles into a spherical quasiequilibrium objectwhose outermost layers undergo damped oscillations.Although a large fraction of the initial mass escapes,the final total rest mass still exceeds the maximum restmass of 1 . M ⊙ that our cold EOS can support (seeFig. 9). However, the pWDNS remnant cannot collapsepromptly to form a black hole, because it is hot. This re-sult is the same in all our simulations. However, we pointout that delayed collapse to a black hole is almost cer-tain after the pWDNS remnant has cooled, but this willoccur on a timescale much larger than a hydrodynamicaltimescale.The final object consists of a cold NS core surroundedby a hot mantle. We quantified the results of shock heat-ing by the ratio of the total pressure to the cold pressure K = P/P cold . In all cases K ≃ K lies in the range[1 ,
15] in cases A1, A2 and A3, [1 ,
35] in case B and [1 , o K. Using a simple scaling argument (see Eq. (51))we find that TZlO temperatures should be proportionalto the compaction of the original pWD, so that in realisticWDNS head-on collisions typical remnant temperatureswould be of order 10 o K.Furthermore, we find that the smaller the initial pWDcompaction the smaller the core densities of the pWDNSremnant. This is supported by the rest-mass density con-0tours shown in Fig. 8 and by the values of the final centralrest-mass density. In particular, we find that the finalcentral rest-mass density decreases by 9 .
2% from case Bto case A2, and 9 .
5% from case A2 to case C. In addition,the final minimum value of the lapse function, which is agood indicator of collapse, increases with increasing ini-tial pWD size, too. Specifically, we find this value to be0 .
57 in case B, 0 .
595 in case A2, and 0 .
609 in case C (seeTable IV for a summary of physical parameters of thefinal configurations in cases B, A2, C). All these factsseem to indicate that as the pWD size increases towardsrealistic WD sizes the less likely it is for the pWDNSremnant to collapse.An important concern regards the invariance of theseresults with respect to larger initial separations. To fullyresolve this, one would need to extend the simulations towider separations, but this extension is outside the scopeof the current work. However, this work gives us somequalitative idea about what might happen with larger ini-tial separations. Larger separations imply larger kineticenergies during the plunge phase, which in turn implystronger shocks. Stronger shocks likely lead to largermass loss and more intense shock heating. Therefore, ourexpectation is that head-on collisions of pWDNS systemsstarting at larger separations will also result in the for-mation of TZlOs and that such collisions would not leadto prompt formation of a black hole.Given the consistency in the results of cases B, A2and C, we expect that as the parameters of our EOS areadjusted such that pWDs more closely resemble realis-tic WDs, WDNS head-on collisions are likely to form aquasiequilibrium TZlO. If the initial total mass of thesystem well exceeds the maximum mass that a cold EOScan support, then the TZlO will most likely have mass ex-ceeding the maximum mass supportable by a cold EOS,eventually collapsing to a black hole after the remnanthas cooled.We conclude by stressing that we cannot use the re-sults of this work to make definite predictions about ei-ther the pWDNS or the realistic WDNS circular binaryproblem. One might speculate that shock heating willbe minimized in such a scenario, and hence it may resultin prompt collapse to a black hole. However, sufficientangular momentum must be shed in the circular binarycase in order for the object to promptly form a blackhole. To resolve these issues, hydrodynamic simulationsin full general relativity must be performed and will bethe focus of a forthcoming paper.
Acknowledgments
We would like to thank Alexei M. Khokhlov andThomas W. Baumgarte for helpful discussions. Thispaper was supported in part by NSF Grants PHY06-50377 and PHY09-63136 as well as NASA GrantsNNX07AG96G and NNX10AI73G to the University ofIllinois at Urbana-Champaign. Z. Etienne gratefully ac- knowledges support from NSF Astronomy and Astro-physics Postodoctoral Fellowship AST-1002667.
Appendix A: Initial data Code description
In this appendix we describe some details of the fixedmesh refinement (FMR), finite difference code we devel-oped for generating general relativistic initial data.The grid structure we use is a multi-level set of prop-erly nested, cell-centered uniform grids. Each grid cor-responds to one level of refinement labeled by the levelnumber il = (0 , , , . . . , nl − nl is the totalnumber of levels. il = 0 corresponds to the coarsest leveland il = nl − nx, ny, nz ∈ Z in the x , y and z directions respectively. The coordinates on our grid aredefined as follows x il,i = x il, min + i · ∆ x il , i = 0 , , . . . , nx − ,y il,j = y il, min + j · ∆ y il , j = 0 , , . . . , ny − ,z il,k = z il, min + k · ∆ z il , k = 0 , , . . . , nz − , (A1)where x il, min , y il, min , z il, min are the minimum valuesof the coordinates in each direction on level il and∆ x il , ∆ y il , ∆ z il are the mesh sizes in each direction onlevel il .The mesh size between two consecutive levels differsby a factor of two so that∆ x il +1 = ∆ x il , (A2)and similarly in the y and z directions.In order for the grids to be properly nested we demandthat there exists an i ∈ [0 , nx −
1] such that x il,i = x il +1 , min −
32 ∆ x il +1 , il = 0 , . . . nl − y and z directions. This conditionensures that two consecutive levels share a common in-terface.We now borrow FMR terminology to name two typesof cells that exist on our grid. These are the split cellsand the leaf cells or leaves. A split cell is one withinwhich there exist higher level cells and a leaf cell is onewithin which there are no higher level cells. The totalnumber of cells N tot on our grid is N tot = nx · ny · nz · nl, (A4)and a straightforward calculation shows that the numberof leaves is N leaf = nx · ny · nz nl + 18 . (A5)When nl = 1, N leaf = N tot , i.e., all cells are leaves, asexpected.1We distinguish between these two types of cells be-cause our solutions are defined only on leaves. This maybe more cumbersome to implement, but has two mainadvantages.First, there is no ambiguity as to how one should inter-polate values of matter sources from fine levels on coarselevels in order to correctly calculate the gravitationalfields. To be more specific, let us assume that we haveone coarse cell which is split into 8 cells and that we knowthe values of the density on the fine cells. In Newtonianphysics, to ensure that the gravitational fields are com-puted correctly (at least far away), all we have to do isset the cell averaged density on the coarse level such thatthe total mass in the coarse cell is the same as the totalmass in the enclosed fine cells. In general relativity thedefinition of gravitational mass depends not only on thedensity, but also on the gravitational fields. Hence, thereis no straightforward way to set the density on coarsecells in GR. The ambiguity is immediately lifted, if onedefines all fields only on the finest cells.The second advantage of using only leaves is that thememory requirements are minimized and the calculationsare carried out faster because N leaf N tot = 7 nl + 18 nl ≤ , (A6)where the last inequality holds because nl ≥ ∇ u = f ( u ) χ, (A7)where f ( u ) is a nonlinear function of the variable u and χ a known scalar independent of u , our code employs astandard second-order finite difference scheme u il,i +1 ,j,k + u il,i − ,j,k − u il,i,j,k ∆ x il + u il,i,j +1 ,k + u il,i,j − ,k − u il,i,j,k ∆ y il + u il,i,j,k +1 + u il,i,j,k − − u il,i,j,k ∆ z il = f ( u il,i,j,k ) χ il,i,j,k . (A8)The finite difference stencil changes only across grid-level boundaries where we perform first order interpola-tion.To address the nonlinearity of Eq. (A7) we performNewton-Raphson iterations as follows. Let us assumethat u n is a guess at step n . We first calculate the residual R n from Eq. (A7) R n = ∇ u n − f ( u n ) χ, (A9)and then solve the linearized equation for the correction δu n on u n ∇ δu n = f ′ ( u n ) χδu n − R n , (A10) where f ′ ( u n ) = (cid:18) df ( u ) du (cid:19) u = u n . (A11)Once a solution to (A10) is found, we correct u n as u n +1 = u n + δu n , (A12)and iterate until this procedure converges and a solutionto Eq. (A7) is obtained.We solve the equations using the PETSc linear solverKrylov Space (KSP) methods. KSP methods are matrixmethods and hence we have to set up the matrix of thelinear system.To do this we define a global index that counts all cells(both leaves and split cells) on our grid as I = il + nl ( i + nx · j + nx · ny · k ) , (A13)In this way, every leaf cell corresponds to a unique index I . However, I takes values 0 , , . . . , N tot −
1, but thereare N leaf leaves on the grid with N leaf ≤ N tot . Hence I cannot be used to count leaves, unless nl = 1. For thisreason, we define another index ic , which counts onlythe leaves on our grid, as well as two mappings; from Ito ic, ic ( I ) and from ic to I, I ( ic ). Since for every cell onour grid we can find I from Eq. (A13) we set up thesemappings by defining two arrays ic of I, I of ic, of length N tot and N leaf , respectively. Looping over il, i, j, k , westore ic in the array ic of I assigning a value of − I in the array I of ic. Theindex I is used when we need the index ic of a neighborleaf cell in order to calculate derivatives or enter matrixelements.For example, let us assume that we are at a leaf cell ofindex ic which is not near a grid-level boundary, and wewant to enter the element of matrix A that correspondsto the right-x neighbor of this cell (where A representsthe Laplacian). From Eq. (A8) this matrix element mustbe 1 / ∆ x il . If ic represents the ic -th row of A we mustfind which column of A the neighbor corresponds to. Wefind this as follows.First, using the mapping from ic to I, we find the index I = I ( ic ) of the leaf. Next by use of Eq. (A13) we de-termine il, i, j, k that correspond to I using the followingsequence of operations k = int ( I/nl · nx · ny ) I = I − nl · nx · ny · kj = int ( I /nl · nx ) I = I − nl · nx · ji = int ( I /nl · nx ) il = I − nl · i, (A14)where int means the integer part of the division.In the next step the global index ( I p1 ) of the right x-neighbor is found, as I p1 = il + nl ( i +1+ nx · j + nx · ny · k ).Finally, using the mapping from I to ic we find the leaf (or2desired column) number ic p1 = ic ( I p1 ). Knowing the col-umn number of the neighbor, it is straightforward to as-sign A ic,ic p1 = 1 / ∆ x il . We use the same approach to setup all the matrix elements of the linear systems and cal-culate derivatives. The algorithm becomes slightly morecomplicated when the neighbor cell is a fictitious cell thatresides on a different level. In such cases we perform firstorder interpolation and use the same method as outlinedabove to find which matrix elements must be filled withnon-zero values. Appendix B: Validation of HRSC method for apiecewise polytropic EOS
In this appendix we analyze the effect our numericalschemes have on solutions obtained with our adoptednon-smooth EOS (15) and a smooth counterpart of thisEOS. We show that there is no essential difference. Thisresult is expected because an algorithm with finite resolu-tion cannot distinguish a smooth EOS from a non-smoothEOS, if the smoothing operation is performed below theresolution level of the computations.For simplicity we consider a non-smooth EOS with twobranches as follows P = κ ρ /n , ρ ρ κ ρ /n , ρ > ρ , (B1)and perform a smoothing operation over a density inter-val [ ρ i (1 − ǫ ) , ρ i (1 + ǫ )] as follows P = κ ρ /n , ρ ρ (1 − ǫ ) f ( ρ ) , ρ (1 − ǫ ) < ρ ρ (1 + ǫ ) ,κ ρ /n , ρ > ρ (1 + ǫ ) , (B2)where f ( ρ ) is a smooth spline fit such that the EOS iscontinuous and has continuous first or second derivative,depending on the order of the spline. Our particularchoice for the smoothing function is either a cubic splineor a quintic spline. In the former case the EOS becomes C , while in the latter case the EOS becomes C .We chose ǫ to be sufficiently small so that the smoothedEOS mimics as closely as possible EOS (B1), but largeenough to avoid round off errors due to very large gra-dients. For the cubic spline we set ǫ = 10 − , whilefor the quintic spline we set ǫ = 10 − . In all our nu-merical tests we choose k , k , n , n , ρ to correspond to k , k , n , n , ρ of the 10:1 EOS (see Table I), respec-tively. In Fig. 13 we show a plot of EOSs (B1) and (B2),where f ( ρ ) is a cubic spline.With the smooth EOS (B2) at our disposal we set upseveral 1D Riemman (or shock tube) problems in a spa-tial domain of length L = 4 . FIG. 13: Pressure vs rest-mass density for EOSs (B1) (black)and (B2) (red). Here P = 10 − km − and ρ nuc = 1 . × − km − . The inset zooms in the region, where EOS (B1) isnon-differentiable, and shows that the cubic spline fit smoothsout the discontinuity. th = 1 .
66 anduse Eq. (40) with P cold given by either EOS (B1) or EOS(B2).We have explored the parameter space of initial data( ρ L , P L , u xL ), ( ρ R , P R , u xR ) for the left and right states,and in all cases we found that the solutions obtainedwith EOS (B1) almost overlap those obtained with thesmooth EOS (B2). These results hold for both piece-wise parabolic (PPM) and monotonized central (MC) re-construction, regardless of resolution and the spline fitchoice. Furtermore, we verified that all our simulationswith the smooth EOS (B2) had high enough resolutionso that data points would sample the smoothing interval[ ρ i (1 − ǫ ) , ρ i (1 + ǫ )] every few timesteps.In Fig. 14 we plot a snapshot of the pressure profileand do a convergence test for one of the cases we sim-ulated with ( ρ R = 10 − , P R = P cold ( ρ R ) , u xR = 0),( ρ L = 5 × − , P L = 10 P R , u xL = 0). The figure showsthat all solutions overlap (left panel) and converge atthe expected order (right panel) to the very high reso-lution solution obtained with PPM in conjunction withthe smooth EOS.The solutions obtained with the smooth EOS (B2) withquintic spline smoothing, which results in a C and aconvex EOS, also overlap with those of the non-smoothEOS solutions even though the smoothing interval wechose was much larger than in the cubic spline case, andhence the data points would sample it more frequently.Note also that the coarsest resolution used in the shock3 FIG. 14: Left: Snapshot of pressure profile at t = 2 .
56 km, corresponding to a third of a sound speed crossing timescaleacross the domain. Right: Convergence plot (at t = 2 .
56 km) using as reference solution the very high resolution solutionobtained in conjunction with the smooth EOS (B2) with a cubic spline smoothing function. The labels in the plots denotethe resolution (LR, MR, HR, or VHR) and the reconstruction method (PPM or MC as subscripts). The resolutions used are:LR = 210 , MR = 420, HR = 840, VHR = 1680 grid points. PPM stands for the piecewise parabolic reconstruction, and MCstands for the monotonized central reconstruction. The plots demonstrate that all solutions overlap (left panel), regardless ofthe reconstruction method and the EOS used (smooth or non-smooth), and first-order convergence to the VHR run with thesmooth EOS, as expected. Here P = 10 − km − . tube tests is at least 20 times higher than the resolu-tion used in our WDNS simulations. Therefore, all theseresults demonstrate that our numerical methods capture the correct solution and indicate that no unphysical solu-tions are present in our simulations of the WDNS head-oncollisions. [1] B. Abbott and the LIGO Scientific Collaboration, ArXive-prints , 0704.3368 (2007).[2] D. A. Brown, S. Babak, P. R. Brady, N. Christensen,T. Cokelaer, J. D. E. Creighton, S. Fairhurst, G. Gon-zalez, E. Messaritaki, B. S. Sathyaprakash, et al., Class.Quant. Grav. , S1625 (2004).[3] F. Acernese and the VIRGO Collaboration, Class. Quant.Grav. , 635 (2006).[4] F. Bauville and the LIGO-VIRGO Working Group,ArXiv e-prints pp. gr–qc/0701027 (2007).[5] H. L¨uck and the GEO600 collaboration, Class. Quant.Grav. , S71 (2006).[6] M. Ando and the TAMA collaboration, Class. Quant.Grav. , 1409 (2002).[7] D. Tatsumi and the TAMA collaboration, Class. Quant.Grav. , 399 (2007).[8] .[9] G. Heinzel, C. Braxmaier, K. Danzmann, P. Gath,J. Hough, O. Jennrich, U. Johann, A. R¨udiger,M. Salusti, and H. Schulte, Class. Quant. Grav. , 119(2006).[10] S. Kawamura and the DECIGO collaboration, Class.Quant. Grav. , 125 (2006). [11] T. W. L. Baumgarte and S. L. Shapiro, Numerical Rela-tivity (Cambridge University Press, 2010).[12] I. Hinder, Classical and Quantum Gravity , 114004(2010), 1001.5161.[13] M. D. Duez, Classical and Quantum Gravity , 114002(2010), 0912.3529.[14] E. Rantsiou, S. Kobayashi, P. Laguna, and F. A. Rasio,Ap. J. , 1326 (2008).[15] F. L¨offler, L. Rezzola, and M. Ansorg, Phys. Rev. D ,104018 (2006).[16] J. A. Faber, T. W. Baumgarte, S. L. Shapiro,K. Taniguchi, and F. A. Rasio, Phys. Rev. D , 024012(2006).[17] J. A. Faber, T. W. Baumgarte, S. L. Shapiro, andK. Taniguchi, Ap. J. Lett. , L93 (2006).[18] M. Shibata and K. Uryu, Phys. Rev. D , 121503(R)(2006).[19] M. Shibata and K. Uryu, Class. Quant. Grav. , 125(2007).[20] M. Shibata and K. Taniguchi, Phys. Rev. D , 084015(2008).[21] T. Yamamoto, M. Shibata, and K. Taniguchi, Phys. Rev.D , 064054 (2008). [22] Z. B. Etienne, J. A. Faber, Y. T. Liu, and S. L. Shapiro,Phys. Rev. D , 084002 (2008).[23] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baum-garte, Phys. Rev. D , 044024 (2009).[24] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer,M. A. Scheel, and S. A. Teukolsky (2008), arXiv:gr-qc/0809.0002.[25] M. Shibata, K. Kyutoku, T. Yamamoto, andK. Taniguchi, Phys. Rev. D , 044030 (2009),0902.0416.[26] K. Kyutoku, M. Shibata, and K. Taniguchi, Phys. Rev.D , 124018 (2009), 0906.0889.[27] P. M. Motl, M. Anderson, M. Besselman, S. Chawla,E. W. Hirschmann, L. Lehner, S. L. Liebling, D. Neilsen,and J. E. Tohline, in Bulletin of the American Astronom-ical Society (2010), vol. 41 of
Bulletin of the AmericanAstronomical Society , pp. 295–+.[28] S. Chawla, M. Anderson, M. Besselman, L. Lehner, S. L.Liebling, P. M. Motl, and D. Neilsen, ArXiv e-prints(2010), 1006.2839.[29] M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, andS. A. Teukolsky, Classical and Quantum Gravity ,114106 (2010), 0912.3528.[30] F. Pannarale, A. Tonita, and L. Rezzolla, ArXiv e-prints(2010), 1007.4160.[31] F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukol-sky, ArXiv e-prints (2010), 1007.4203.[32] K. Kyutoku, M. Shibata, and K. Taniguchi, ArXiv e-prints (2010), 1008.1460.[33] V. Paschalidis, M. MacLeod, T. W. Baumgarte, and S. L.Shapiro, PRD , 024006 (2009).[34] T. W. Baumgarte, S. L. Shapiro, and M. Shibata, Ap. J.Lett. , L29 (2000), arXiv:astro-ph/9910565.[35] I. A. Morrison, T. W. Baumgarte, and S. L. Shapiro, Ap.J. , 941 (2004).[36] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Ap. J. , 227 (1994).[37] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall,Phys. Rev. C , 1804 (1998).[38] C. P. Lorenz, D. G. Ravenhall, and C. J. Pethick, Phys-ical Review Letters , 379 (1993).[39] R. B. Wiringa, V. Fiks, and A. Fabrocini, Phys. Rev. C , 1010 (1988).[40] V. R. Pandharipande and R. A. Smith, in Bulletin of theAmerican Astronomical Society (1975), vol. 7 of
Bulletinof the American Astronomical Society , pp. 240–+.[41] H. A. Bethe and M. B. Johnson, Nucl. Phys. A , 1(1974).[42] V. R. Pandharipande, Nucl. Phys. A , 641 (1971).[43] G. Nelemans, L. R. Yungelson, and S. F. P. Zwart, As-tron. and Astrop. , 890 (2001).[44] A. Cooray, MNRAS , 25 (2004).[45] T. A. Thompson, M. D. Kistler, and K. Z. Stanek (2009),arXiv:0912.0009.[46] S. Rappaport, P. C. Joss, and R. F. Webbink, Ap. J. ,616 (1982).[47] S. Rappaport, F. Verbunt, and P. C. Joss, Ap. J. ,713 (1983).[48] F. Verbunt and S. Rappaport, Ap. J. , 193 (1988).[49] P. Podsiadlowski, P. C. Joss, and J. J. L. Hsu, Ap. J. , 246 (1992).[50] T. R. Marsh, G. Nelemans, and D. Steeghs, MNRAS ,113 (2004).[51] W. Benz, R. Bowers, A. Cameron, and W. Press, Ap. J. , 647 (1990).[52] F. A. Rasio and S. L. Shapiro, Ap. J. , 887 (1995).[53] L. Segretain, G. Chabrier, and R. Mochkovitch, Ap. J. , 355 (1997).[54] J. Guerrero, E. Garcia-Berro, and J. Isern, Astron. andAstroph. , 257 (2004).[55] S.-C. Yoon, P. Podsiadlowski, and S. Rosswog, MNRAS , 933 (2007).[56] M. Dan, S. Rosswog, and M. Br¨uggen (2008),ArXiv:astro-ph/0811.1517.[57] K. Thorne and A. Zytkow, Ap. J. , 832 (1977).[58] Z. Etienne, Y. T. Liu, and S. L. Shapiro (2010),arXiv:astro-ph/1007.2848.[59] S. L. Shapiro and S. A. Teukolsky,
Black Holes, WhiteDwarfs, and Neutron Stars (John Willey and Sons, 1983).[60] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman,Phys. Rev. D , 124032 (2009).[61] R. Arnowitt, S. Deser, and C. Misner, Gravitation: AnIntroduction to Current Research (Cambridge UniversityPress, Wiley, New York, 1962).[62] T. W. Baumgarte and S. L. Shapiro, Phys. Rep. , 41(2003).[63] T. W. Baumgarte, G. B. Cook, M. A. Scheel, S. L.Shapiro, and S. A. Teukolsky, Phys. Rev. D , 7299(1998).[64] S. Balay, K. Buschelman, W. D. Gropp,D. Kaushik, M. G. Knepley, L. C. McInnes, B. F.Smith, and H. Zhang, PETSc Web page
Modern Software Tools in Scientific Comput-ing , edited by E. Arge, A. M. Bruaset, and H. P. Lang-tangen (Birkh¨auser Press, 1997), pp. 163–202.[67] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens,Phys. Rev. D , 024028 (2005).[68] M. Shibata and T. Nakamura, Phys. Rev. D , 5428(1995).[69] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D ,024007 (1998).[70] J. R. van Meter, J. G. Baker, M. Koppitz, and D.-I. Choi,Phys. Rev. D , 1465 (2004), arXiv:gr-qc/0310042.[73] P. Colella and P. R. Woodward, Journal of Computa-tional Physics , 174 (1984).[74] Harten, A. and Lax, P. D. and van Leer, B., SIAM Rev. , 35 (1983).[75] Z. Etienne, Y. T. Liu, S. Shapiro, and T. Baumgarte,APS Meeting Abstracts pp. 11007–+ (2009).[76] R. J. Leveque, Finite Volume Methods for HyperbolicProblems (Cambridge University Press, 2002).[77] V. Alexander, Ph.D. thesis, RWTH Aachen, Aachen,Germany (2005), online at http://darwin.bth.rwth-aachen.de/opus/volltexte/2005/1210/.[78] J. W. Murphy, C. D. Ott, and A. Burrows,The Astrophysical Journal , 1173 (2009), URL http://stacks.iop.org/0004-637X/707/i=2/a=1173 .[79] M. Saijo and I. Hawke, Phys. Rev. D80