Dynamical evolution of black hole-neutron star binaries in general relativity: Simulations of tidal disruption

We calculate the first dynamical evolutions of merging black hole-neutron star binaries that construct the combined black hole-neutron star spacetime in a general relativistic framework. We treat the metric in the conformal flatness approximation, and assume that the black hole mass is sufficiently large compared to that of the neutron star so that the black hole remains fixed in space. Using a spheroidal spectral methods solver, we solve the resulting field equations for a neutron star orbiting a Schwarzschild black hole. The matter is evolved using a relativistic, Lagrangian, smoothed particle hydrodynamics (SPH) treatment. We take as our initial data recent quasiequilibrium models for synchronized neutron star polytropes generated as solutions of the conformal thin-sandwich (CTS) decomposition of the Einstein field equations. We are able to construct from these models relaxed SPH configurations whose profiles show good agreement with CTS solutions. Our adiabatic evolution calculations for neutron stars with low-compactness show that mass transfer, when it begins while the neutron star orbit is still outside the innermost stable circular orbit, is more unstable than is typically predicted by analytical formalisms. This dynamical mass loss is found to be the driving force in determining the subsequent evolution of the binary orbit and the neutron star, which typically disrupts completely within a few orbital periods. The majority of the mass transferred onto the black hole is accreted promptly; a significant fraction ($\sim 30\%$) of the mass is shed outward as well, some of which will become gravitationally unbound and ejected completely from the system. The remaining portion forms an accretion disk around the black hole, and could provide the energy source for short-duration gamma-ray bursts.

Figure 1

The critical mass ratio $q=M_{*}/M_{\mathrm{BH}}$ as a function of the secondary’s compaction $\mathcal{C}=M_{*}/R_{*}$, for which mass transfer begins at the ISCO, taken here as $a=6M_{\mathrm{BH}}$ (solid line). For systems above the curve, mass transfer begins while the orbit is stable; for those below, the secondary may plunge into the BH before being tidally disrupted. Dashed vertical curves show characteristic compactions of a WD or NS; dotted horizontal curves show typical mass ratios for a $1M_{\odot }$ compact object orbiting a $10M_{\odot }$ stellar-mass BH, a ${10}^3{M}_{\odot }$ IMBH, or a ${10}^6{M}_{\odot }$ SMBH. Dot-dashed curves show where $\beta$, the ratio of the light crossing time the viscous time scale of the NS, equals unity, and where ${\alpha }_{\mathrm{Vis}}$, the nondimensional turbulent viscosity assumes its largest reasonable value for the turbulent viscosity [See Eqs. (7, 9)]. Only configurations to the left of these curves can synchronize prior to merger.

Figure 2

A schematic representation of the spectral methods computational domains used to solve the NS components of the field equations, Eqs. (24, 25). Initially, a triaxial ellipsoid with surface $R_1(\theta ,\varphi )$ and origin at the NS center-of-mass is fitted around all SPH particles for which $({\rho }_{*}{)}_i>{\rho }_{\mathrm{crit}}$ (top panel). Two annular “shell-like” domains with spherical outer boundaries are also laid down, with radii $R_2$ and $R_3$, twice and 3 times the maximum value of $R_1$. During the evolution (bottom panel), we use the same procedure to fit $R_1(\theta ,\varphi )$, keep the value of $R_3$ fixed, and calculate $R_2$ as the mean of $R_3$ and the maximum of $R_1$. Some particles first leave the innermost domain, and then the entire computational volume, particularly those accreted by the black hole, shown as a circle of radius $0.5M_{\mathrm{BH}}$ centered at the origin.

Figure 3

From top to bottom, the values of ${\rho }_{*}$, ${\beta }^y$, $\psi$, and $\alpha$ along the x-axis for the SPH (solid line) and grid-based (dashed line) data representing configuration B3, a NS with a $\Gamma =2$ polytropic EOS and an initial binary separation $a_0/M_{\mathrm{BH}}=11.1$. The agreement is generally to within $1-2\%$.

Figure 4

The values of ${\rho }_{*}$, ${\beta }^y$, $\psi$, and $\alpha$ along the x-axis for the SPH (solid line) and grid-based (dashed line) data representing configuration A5, a NS with a $\Gamma =1.5$ polytropic EOS and an initial binary separation $a_0/M_{\mathrm{BH}}=11.8$. Conventions are the same as Fig. 3.

Figure 5

Value of the SPH expression for the integrated Euler equation constant $h/u^0$ as a function of the particle’s radius from the NS center-of-mass for configurations B3 (top) and A5 (bottom), featuring a NS EOS with $\Gamma =2$ and $\Gamma =1.5$, respectively. Note that the proper value differs between the two cases. The standard deviation in both cases is $<0.01\%$, with maximum variation $<0.1\%$.

Figure 6

SPH particle configurations at times $t=0$, 286, 566, and 693, projected into the orbital plane, for run A5 (for which the orbital period is $T=287M_{\mathrm{BH}}$). These configurations correspond to the initial configuration, and after $1$, $2$, and 2.5 full orbits, respectively. The BH is shown as a circle at a radius $r=0.5M_{\mathrm{BH}}$. We see that once mass transfer begins, the NS begins to expand, forming a low-density single-armed spiral pattern. Note that particles have different masses, and that those in the center of the NS are in most cases significantly more massive than those originally in the outer layers.

Figure 7

Binary separation (top panel) and total mass lost from the NS (bottom panel) as a function of time for all runs calculated using the soft NS EOS, runs A1–A7, which have $\Gamma =1.5$. The properties of the initial configurations, which differ only in their binary separation, are shown in Table . Here, particles are defined as “lost” when they lie outside the innermost computational domain around the NS. In all cases, mass loss never quenches once it begins, eventually leading to the complete disruption of the NS.

Figure 8

Total mass lost inward toward the BH (top panel) and outward away from the BH (bottom panel) as a function of time, for the runs shown in Fig. 7. In general, the smaller the initial separation, the more mass that falls inward toward the BH. For the runs started from a larger separation, much of the initial mass loss results from spurious numerical of low-mass particles out of the Roche lobe, leading to nearly equal flows directed inward and outward. In contrast, the mass loss for the closer cases is in exactly the form expected for Roche-lobe overflow through the inner Lagrange point.

Figure 9

Binary separation (top panel) and total mass lost from the NS (bottom panel) as a function of time for the runs calculated using the stiff NS EOS without GW radiation-reaction, runs B1-B4, which have $\Gamma =2$. Conventions are as in Fig. 7. For run B1, we see periodic mass loss, as the NS gets kicked into an elliptical orbit, losing mass during every periastron passage. Runs B2 and B3 are similar to those with the softer EOS, as the NS gets disrupted during the first mass loss phase. Run B4 is essentially stable, only showing the long-term effects of numerical diffusion of particles.

Figure 10

Total mass lost inward toward the BH (top panel) and outward away from the BH (bottom panel) as a function of time, for the runs shown in Fig. 9. Conventions are as in Fig. 8. Approximately 50% more mass is lost inward for runs B2 and B3, started from near the stability limit, whereas for run B1, marked by an elliptical orbit and periodic mass transfer, nearly $4$ times much mass is lost toward the BH.

Figure 11

SPH particle configurations at times $t/M_{\mathrm{BH}}=0$, 261, 514, and 767, projected into the orbital plane, for run B3a, which includes the dissipative effects of radiation-reaction (for this configuration, $T=265M_{\mathrm{BH}}$). In terms of the binary orbit, these correspond to the initial configuration and $1$, $2$, and $3$ full orbits, respectively. Conventions are as in Fig. 6. Here, we see that radiation-reaction initially drives only an inwardly directed mass-transfer stream onto the BH, until somewhere after $t/M_{\mathrm{BH}}=700$, when the expansion of the NS becomes unstable and tidal disruption occurs.

Figure 12

Binary separation (top panel) and total mass lost from the NS (bottom panel) as a function of time for the runs calculated using the stiff NS EOS and dissipative GW radiation-reaction effects, runs B3a and B3b, which have $\Gamma =2$. Conventions are as in Fig. 7. The only difference between the two runs is the magnitude of the radiation-reaction effects; we use the quadrupole-order form in run B3a, and double its strength for run B3b. We find that doubling the radiative drag force forces the binary to a smaller separation, leading to a larger initial burst of mass loss and a more rapid increase in the binary separation. As a result, the system tidally disrupts slower than the case shown in run B3a.

Figure 13

Total mass lost inward toward the BH (top panel) and outward away from the BH (bottom panel) as a function of time, for the runs shown in Fig. 9. Conventions are as in Fig. 8. Initially, in run B3b with twice the physical GW radiation-reaction force applied, all mass lost in run B3b is directed inward. This leads to a rapid increase in the binary separation and slows the growth of the mass-transfer rate. In contrast, in run B3a, mass loss is more evenly balanced between inwardly and outwardly directed flows, and the NS disrupts more quickly.

Figure 14

Gravity wave emission angular frequency (top panel) and waveform in both polarizations, $h_{+}$ (solid curve) and $h_{\times }$ (dashed curve), defined by Eqs. (76, 77), for the merger calculated in run B3b. Initially, we see conclusion of the standard “chirp” signal, in which the binary separation decreases while the GW amplitude and frequency increases, all of which happens extremely gradually on the secular GW time scale. This lasts until the onset of mass transfer, at which point we encounter a much more rapid “reverse chirp”, as the GW amplitude and frequency rapidly decrease while the NS is tidally disrupted.

Figure 15

Critical compactness for determining the stability of mass transfer from a relativistic polytropic NS companion as a function of the binary mass ratio, defined in terms of the polytropic index $n$ (bottom), or the adiabatic index $\Gamma \equiv 1+1/n$. Mass transfer is unstable to the right of the curve, stable to the left. The square and triangle represent the positions of the two models we are evolving dynamically.

Figure 16

Critical mass ratio for determining the stability of mass transfer from a relativistic polytropic NS companion as a function of the NS compactness, defined in terms of the polytropic index $n$ (bottom), or the adiabatic index $\Gamma \equiv 1+1/n$. Mass transfer is unstable to the right of the curve, stable to the left. The heavy solid lines depict the maximum possible compaction for a given polytropic EOS. The square and triangle represent the positions of the two models we are evolving dynamically.