In particular, the understanding of the magnetospheres of neutron stars has recently improved through time-dependent numerical simulations, which allow one to investigate the stability of steady-state analytic solutions of the pulsar equation. Existing attempts have used both ideal relativistic magnetohydrodynamics and resistive force-free electrodynamics (see  and references therein). The investigation by  permits one to rule out some stationary solutions, which are not naturally recovered by the time-dependent simulation. Further work is needed to account for resistive relativistic MHD computations.
Another remarkable investigation based on fully-relativistic MHD simulations by  has permitted us to understand the properties of the flow produced by a magnetized pulsar wind within a plerionic nebula. The results of the axisymmetric simulations have disclosed the complex dynamics of the post-shock flow, very different from the steady quasi-radial outflow assumed in earlier (spherically-symmetric) analytical models for plerions. In particular, the jet-torus pattern discovered by Chandra in the Crab nebula can be explained within the MHD approximation when the condition of spherical symmetry is no longer enforced, as shown by the simulated synchrotron X-ray images of the MHD numerical data.
General-relativistic MHD winds from rotating neutron stars have also been studied by . Magnetically-dominated winds from stars are central to the angular momentum evolution of these objects, and GRMHD studies, such as those carried out by , may help us understand the observed spin-down of rotation-powered pulsars with strong magnetic fields.
The use of relativistic hydrodynamic codes to study the stability properties of neutron stars and to compute mode frequencies of oscillations of such objects has increased in recent years (see, e.g., the Living Reviews article by Stergioulas  and references therein). In the case of GRMHD codes, the first investigations are only just starting. An obvious advantage of the hydrodynamic approach is that it includes, by construction, nonlinear effects, which are important in situations where the linearized equations (commonly employed in calculations of the mode frequencies of pulsating stars) break down. In addition, while for nonrotating stars the frequencies of normal modes can be computed with perturbative methods and a theory of gravitational wave asteroseismology has already been formulated, there exist no accurate frequency determinations for rapidly rotating stars to date, nor has the theory of gravitational-wave asteroseismology been extended to include the effects of rotation on the oscillation frequencies. Due to the advanced status of current hydrodynamics codes, however, the computation of mode frequencies in rapidly rotating relativistic stars might be easier to achieve through nonlinear numerical evolutions than by using perturbative computations (see, e.g., the results in [131, 355, 364, 389, 99]).
Hydrodynamic evolutions of polytropic models of spherical neutron stars can be used as test-bed computations for multidimensional codes. Representative examples are the simulations by , with pseudo-spectral methods, and by  with HRSC schemes. These investigations adopted radial-gauge polar-slicing coordinates in which the general relativistic equations are expressed in a simple way that resembles Newtonian hydrodynamics. Gourgoulhon  used a numerical code to detect, dynamically, the zero value of the fundamental mode of a neutron star against radial oscillations. Romero et al.  highlighted the accuracy of HRSC schemes by finding, numerically, a change in the stability behavior of two slightly different initial neutron star models: for a given EOS, a model with mass is stable and a model of is unstable. More recently, in  a method based on the nonlinear evolution of deviations from a background stationary-equilibrium star was applied to study nonlinear radial oscillations of a neutron star. The accuracy of the approach permitted a detailed investigation of nonlinear features such as quadratic and higher-order mode coupling and nonlinear transfer of energy.
Axisymmetric pulsations of rotating neutron stars can be excited in several scenarios, such as core collapse, crustquakes and corequakes, and binary mergers, and they could become detectable either via gravitational waves or high-energy radiation. An observational detection of such pulsations would yield valuable information about the EOS of relativistic stars. In recent years, the time evolution of the nonlinear equations governing the dynamics of matter and spacetime has been introduced as a promising new approach for computing mode frequencies [138, 130, 131, 390, 389, 99]. For small amplitudes, the obtained frequencies are in excellent agreement with those expected by linear perturbation theory, while two-dimensional eigenfunctions can be obtained through a Fourier transform technique.
As a first step towards the study of pulsations of rapidly-rotating relativistic stars, Font, Stergioulas, and Kokkotas  developed an axisymmetric numerical code called ToniK (see Table 1) that integrates the GRHD equations in a fixed-background spacetime. The finite-difference code is based on a state-of-the-art approximate Riemann solver  and incorporates different second and third-order TVD and ENO numerical schemes. This code can accurately evolve rapidly rotating stars for many rotational periods, even for stars at the mass-shedding limit. The test simulations reported in  show that, for nonrotating stars, small amplitude oscillations have frequencies that agree to better than 1% with linear normal-mode frequencies (radial and nonradial) in the Cowling approximation (i.e., when the evolution of the spacetime variables is neglected). Axisymmetric modes of pulsating nonrotating stars are computed in , both in Cowling and in fully-coupled evolutions. Contrary to the 2+1 approach followed by , the code used in  evolves the relativistic stars on null spacetime foliations (see Section 2.2.2).
Until very recently (see below), the quasi-radial modes of rotating relativistic stars had been studied only under simplifying assumptions, such as in the slow-rotation approximation or in the relativistic Cowling approximation. An example of the latter is presented in , where a comprehensive study of all low-order and axisymmetric modes of uniformly and rapidly-rotating relativistic stars was presented, using the code of . This was done for a sequence of appropriately-perturbed stationary rotating stars, from the nonrotating limit to the mass-shedding limit. The frequencies of the axisymmetric modes are affected significantly by rotation only when the rotation rate exceeds about 50% of the maximum allowed. As expected, at large rotation rates, apparent mode crossings between different modes appear. A comparison of the mode-frequencies computed with ToniK was carried out by  with his pizaa code (see Table 1), which offers an improved treatment of quasistationary scenarios. The agreement found was satisfactory.
In , the first mode frequencies of uniformly-rotating stars in full general relativity and rapid rotation were obtained, using the three-dimensional code gr_astro. Such frequencies were computed both in fixed spacetime evolutions and in coupled hydrodynamic and spacetime evolutions. The Cowling runs allowed a comparison with earlier results reported by , obtaining agreement at the 0.5% level. The fundamental mode frequencies and their first overtones obtained in fully coupled evolutions show a dependence on the increased rotation, which is similar to the one observed for the corresponding frequencies in the Cowling approximation .
Small-amplitude, nonlinear pulsations of both uniformly and differentially rotating neutron stars were studied by  employing the ToniK code. The centrifugal forces and the degree of differential rotation have significant effects on the mode eigenfunction and it was found that near the mass-shedding limit, the pulsations are damped due to shocks forming at the surface of the star. This damping mechanism may set a small saturation amplitude for modes that are unstable to the emission of gravitational waves. It was also found that the fundamental quasi-radial mode is split, at least in the Cowling approximation and mainly in differentially rotating stars, into two different sequences. This study was revisited by  using the CoCoNuT code, which solves the GRHD equations for a conformally-flat three-metric (see Section 4.4). Differential rotation significantly shifts mode frequencies to smaller values, increasing the likelihood of detection by current gravitational wave interferometric detectors. An extended avoided crossing between the and first overtones was observed (previously known to exist from perturbative studies), which is important for correctly identifying mode frequencies in case of detection. For uniformly rotating stars near the mass-shedding limit, the existence of the mass-shedding-induced damping of pulsations was confirmed, even though the effect is not as strong as was previously found in the Cowling approximation .
Neutron stars following a core collapse supernova are rotating at birth and can be subject to various nonaxisymmetric instabilities (see, e.g.,  for a review). Among those, if the rotation rate is high enough that the ratio of rotational kinetic energy to gravitational potential energy , , exceeds the critical value inferred from studies with incompressible Maclaurin spheroids, the star is subject to a dynamic bar-mode ( -mode) instability driven by hydrodynamics and gravity. Its study is highly motivated these days as such an instability bears important implications in the prospects of the detection of gravitational radiation from newly-born rapidly-rotating neutron stars.
Newtonian simulations of the bar-mode instability from perturbed equilibrium models of rotating stars have shown that and is quite independent of the stiffness of the EOS provided the star is not strongly differentially rotating. The relativistic simulations of  yield a value of for the onset of the instability, while the dynamics of the process closely resembles that found in Newtonian theory, i.e., unstable models with large enough develop spiral arms following the formation of bars, ejecting mass and redistributing the angular momentum. As the degree of differential rotation becomes higher and more extreme, Newtonian simulations have also shown that can be as low as .
Get Flash to see this player.
Further relativistic hydrodynamic simulations of the dynamic bar-mode instability have been performed by  using the whisky code and extending the set of models considered by . This study focused on investigating the persistence of the bar deformation once the instability has reached its saturation and on the precise determination of the parameter signalling the threshold for the onset of the instability. Generic nonlinear mode-coupling effects were found to appear during the development of the instability, which can severely limit the persistence of the bar deformation and even suppress the instability altogether. An animation of one such simulation can be seen in the movie of Figure 16.
A word of caution is needed here, as it remains unclear whether the requirements inferred from numerical simulations are at all met by core-collapse progenitors. As shown by  magnetic torques can spin down the core of the progenitor, which leads to slowly rotating neutron stars at birth ( 10 – 15 ms). The most recent, state-of-the-art computations of the evolution of massive stars, which include angular momentum redistribution by magnetic torques and spin estimates of neutron stars at birth , lead to core-collapse progenitors, which do not seem to rotate fast enough to guarantee the unambiguous growth of the canonical bar-mode instability.
Regarding the CFS instability, relativistic hydrodynamic simulations of nonlinear -modes are presented in , using a version of the gr_astro code (see also  for Newtonian simulations). The gravitational radiation reaction-driven instability of the -modes might have important astrophysical implications, provided the instability does not saturate at low amplitudes by nonlinear effects or by dissipative mechanisms. Regarding the latter, existing studies have shown that the mode amplitude could be limited by shear and bulk viscosity, by energy loss to a magnetic field driven by differential rotation, by shock waves, or by the leak of the -mode energy into some short wavelength oscillation modes (see  and references therein). The latter mechanism would dramatically limit the -mode amplitude to a small value of , much smaller than those found in the early simulations of [390, 224] for isentropic stars (see  for a complete list of references on the subject). Energy leak of the -mode into other fluid modes was also considered by  through Newtonian hydrodynamic simulations, finding a catastrophic decay of the amplitude only once it has grown to a value larger than that reported by . Detailed analysis showed that this breakdown is due to nonlinear 3-mode coupling of the -mode to two other inertial modes. The saturation amplitude is of . It is believed that saturation amplitudes of this order are still appropriate to detect gravitational waves from the -mode instability in LMXBs.
As mentioned in Section 5.1.2 the axisymmetric collapse of differentially-rotating magnetized neutron stars has been studied by [105, 360] in full general relativity. The effects of magnetic fields on the magnetic braking of a magnetized differentially-rotating relativistic star have been studied by . (Earlier works in the literature dealt with idealized, incompressible, uniform-density stars.) In order to overcome the long evolution (Alfvén) timescale involved, of the order of for the models simulated,  adopted a perturbative metric approach and considered the limit of slow rotation and weak magnetic fields (still astrophysically realistic). This approach, however, permitted  to use a sufficiently high spatial resolution to track the development of the MRI and see its effects on the dynamics. The most important result found by  is that, independent of resolution, magnetic braking indeed drives the star toward uniform rotation within the timescales considered.
Many efforts in code development in numerical relativity and relativistic astrophysics aim to simulate the coalescence of compact binaries, composed of either two black holes, two neutron stars, or a black hole and a neutron star. The solution of the black-hole–binary problem has seen major breakthroughs in recent years, to the point that long-term stable evolutions for many orbits are now possible, and for various combinations of the parameter space of the problem, namely the spins and linear momentum of the black holes and their mass ratio (see  for a review). In this section we focus on simulations of neutron-star-binaries and neutron-star–black-hole configurations, for which progress has been equally dramatic.
Neutron-star binaries are among the most promising sources of gravitational radiation to be detected by the various ground-based interferometers worldwide. The computation of the gravitational waveform during the most dynamic phase of the coalescence and plunge depends crucially on hydrodynamic, finite-size effects. This phase begins once the stars, initially in quasi-equilibrium orbits of gradually smaller orbital radius (due to the emission of gravitational waves) reach the innermost stable circular orbit (ISCO). Approximately years after formation of the binary system, the gravitational wave frequency enters the LIGO/VIRGO high frequency band. The final plunge of the two objects takes place on a dynamic timescale of a few ms. During the last 15 minutes before the stars finally merge, the frequency is inside the LIGO/VIRGO sensitivity range. About 16,000 cycles of waveform oscillation can be monitored, while the frequency gradually shifts from 10 Hz to 1 kHz. A perturbative treatment of gravitational radiation in the quadrupole approximation is valid as long as and simultaneously, being the total mass of the binary, the neutron star radius, and the separation of the two stars. As the stars approach each other and merge, both inequalities are less valid and therefore, fully-relativistic hydrodynamic calculations become necessary.
Duez et al.  carried out the analytic modelling of the inspiral of corotational and irrotational relativistic neutron-star binaries as a sequence of quasi-equilibrium configurations, and computed the gravitational wavetrain from the last phase (a few hundred cycles) of the inspiral. These authors further showed a practical procedure to construct the entire wavetrain through coalescence by matching the late inspiral waveform to the one obtained by fully-relativistic hydrodynamic simulations. Detailed theoretical waveforms of the inspiral and plunge similar to those reported by  are crucial to enhance the chances of experimental detection in conjunction with matched-filtering techniques.
The accurate simulation of a neutron-star–binary coalescence is, however, one of the most challenging tasks in numerical relativity. These scenarios involve strong gravitational fields, matter motion with (ultra) relativistic speeds, relativistic shock waves, and strong magnetic fields. The numerical difficulties are exacerbated by the intrinsic multidimensional character of the problem and by the inherent complexities in Einstein’s theory of gravity, such as coordinate degrees of freedom and the possible formation of curvature singularities (e.g., collapse of matter configurations to black holes). It is thus not surprising that a large number of the available simulations have been attempted in the Newtonian (and post-Newtonian) framework (see  for a review). Many of these studies employ Lagrangian particle methods such as SPH, and only a few have considered (less viscous) high-order finite-difference methods such as PPM . Notwithstanding the difficulties, major progress has been achieved recently in simulating neutron-star–binary mergers within a relativistic framework.
Concerning relativistic simulations, Wilson’s formulation of the hydrodynamic equations (see Section 2.1.2) was used in [416, 418, 242]. Such investigations assumed a conformally-flat 3-metric. These early simulations revealed the unexpected appearance of a “binary-induced collapse instability” of the neutron stars, prior to the eventual collapse of the final merged object. This effect was reduced, but not eliminated fully, in revised simulations , after Flanagan  pointed out an error in the momentum-constraint equation as implemented by Wilson et al. [416, 418]. A summary of this controversy can be found in . Subsequent numerical simulations with the full set of Einstein equations (see below) did not find this effect.
The effect was neither found in the new set of CFC simulations performed by [300, 299]. In the first investigation a new 3D SPH code was presented (see Section 4.4), tested, and applied to the coalescence of neutron-star binaries modelled by relativistic polytropes with , and . Not only the system dynamics was studied but also gravitational-wave signals and luminosities were computed. The second investigation improved the microphysical treatment of the neutron-star models by using several nonzero temperature EOS, such as those of Shen et al.  and Lattimer–Swesty . The results were compared with earlier investigations employing the BSSN equations performed by  (see below). It was found that the dynamics and outcome of the models depends strongly on the EOS and on the binary parameters (neutron-star masses and spins). Larger torus masses are found for asymmetric systems (up to for a neutron-star mass ratio of 0.55). Within tens of dynamic timescales only the Lattimer–Swesty EOS leads to a collapse of the postmerger remnant to a black hole. The gravitational wave emission from these simulations was analyzed in , where the authors concluded that the peak frequency of the post-merger signal is a sensitive indicator of the EOS, provided the total mass can be determined from the inspiral chirp signal.
Nakamura et al. started a program in the late 1980’s to simulate neutron-star–binary coalescence in general relativity (see, e.g., ). The group developed a three-dimensional code that solves the full set of Einstein and hydrodynamics equations through finite differences in a uniform Cartesian grid, using van Leer’s scheme  with TVD flux limiters. Shockwaves are spread out using a tensor artificial-viscosity algorithm. The hydrodynamic equations follow Wilson’s Eulerian formulation and the ADM formalism is adopted for the Einstein equations. This code has been tested with the study of the gravitational collapse of a rotating polytrope to a black hole (comparing to the axisymmetric computation of Stark and Piran ). A few results of neutron-star–binary simulations including gravitational radiation are reported in .
Get Flash to see this player.
The headon collision of two neutron stars (a limiting case of a coalescence event) was considered by Miller et al. , who performed time-dependent relativistic simulations using the gr_astro code. These simulations analyzed whether the collapse of the final object occurs in prompt timescales (a few milliseconds) or delayed (after neutrino cooling) timescales (a few seconds). In  it was argued that in a headon collision event, sufficient thermal pressure is generated to support the remnant in quasi-static equilibrium against (prompt) collapse, prior to slow cooling via neutrino emission (delayed collapse). In , prompt collapse to a black hole was found in the headon collision of two neutron stars modeled by a polytropic EOS with . The stars, initially separated by a proper distance of , were boosted toward one another at a speed of (the Newtonian infall velocity). The simulation employed a Cartesian grid of points. The time evolution of this simulation can be followed in the Quicktime movie in Figure 17. This animation simultaneously shows the rest-mass density and the internal energy evolution during the on-axis collision. The formation of the black hole in prompt timescales is signalled by the sudden appearance of the apparent horizon at ( in code units). The violet dotted circles indicate the trapped photons. The animation also shows a moderately-relativistic shockwave (Lorentz factor ) appearing at (code units; ; yellow-white), which eventually is followed by two opposite moving shocks (along the infalling direction) that propagate along the atmosphere surrounding the black hole. We note that, using the same code, critical phenomena has recently been found by  in the headon collison of two neutron stars.
Using the code of , Marronetti et al.  performed the first dynamic determination in relativistic gravity of the ISCO of a neutron-star binary (modelled as polytropes). These authors were able to bracket the location of the ISCO by distinguishing stable circular orbits from unstable plunges. It was found that, for the simplified models considered, the ISCO frequency varies with compaction, but does not depend strongly on the stellar spin.
The most comprehensive study of neutron star coalescence in full general relativity has been performed by Shibata et al. [353, 373, 374, 357, 372, 370]. The numerical code, briefly described in Section 4.4, allows the long-term simulation of the coalescences of both irrotational and corotational binaries, from the ISCO up to the formation and ringdown of the final collapsed object (either a black hole or a stable neutron star). The code also includes an apparent horizon finder along with microphysical EOS, and can extract the gravitational waveforms emitted in the collisions. The early simulations of Shibata and Uryu  accounted for a large sample of parameters of the binary system, such as the compactness of the (equal mass) neutron stars (), the adiabatic index of the -law EOS (), and the maximum density, rest mass, gravitational mass, and total angular momentum. The initial data correspond to quasiequilibrium states, either corotational or irrotational, the latter being more realistic from considerations of viscous versus gravitational-radiation timescales. These initial data are obtained by solving the Einstein constraint equations and the equations for the gauge variables under the assumption of a conformallyflat 3-metric and the existence of a helicoidal Killing vector (see  for a detailed explanation). The binaries are chosen at the innermost orbits for which the Lagrange points appear at the inner edge of the neutron stars, and the plunge is induced by reducing the initial angular momentum by 2 – 3%.
The comprehensive parameter-space survey carried out by [353, 373, 374] shows that the final outcome of the merger depends sensitively on the initial compactness of the neutron stars before plunge. Hence, depending on the stiffness of the EOS, controlled through the value of , if the total rest mass of the system is 1.3 – 1.7 times larger than the maximum rest mass of a spherical star in isolation, the end product is a black hole. Otherwise, a marginally-stable massive neutron star forms. In the latter case, the star is supported against self-gravity by rapid differential rotation. The star could eventually collapse to a black hole once sufficient angular momentum has dissipated via neutrino emission or gravitational radiation. In turn, the different outcome of the merger is imprinted in the gravitational waveforms . Therefore, future detection of high-frequency gravitational waves could help to constrain the maximum allowed mass of neutron stars. In addition, for prompt black-hole formation, a disk orbiting around the black hole forms, with a mass less than 1% the total rest mass of the system.
The input physics of those early simulations was later improved in  where hybrid EOS were adopted to mimic realistic nuclear equations of state. In this approach the EOS is divided into two parts, , where is the cold part for which a fitting formula for realistic EOS of cold nuclear matter was assigned. The simulations used, in particular, the SLy and FPS EOS for which the maximum allowed ADM mass of cold and spherical neutron stars is and , respectively. Apart from improving the neutron-star models, the simulations of  also considered mergers of unequal-mass neutron stars. The underlying motivation was to find out whether massive accretion disks form depending on the initial mass ratio. Disk formation during neutron-star–binary coalescence, a fundamental issue for cosmological models of short duration GRBs, was indeed found to be enhanced for unequal-mass neutron stars, in which the less massive star is tidally disrupted during the evolution. The simulations also showed, in broad agreement with the simplified models, that depending on the threshold of the ADM mass of the system, which is also EOS dependent, the outcome is a black hole or a hypermassive neutron star with large ellipticity. For an ADM mass larger than the threshold, a black hole forms promptly in the merger irrespective of the mass ratio.
An example of such findings is shown in Figure 18. Both panels in this figure correspond to equal-mass neutron-star–binary mergers. The left panel is a merger, while the right panel depicts the case. Likewise, both panels display the velocity field and rest-mass isodensity contours in the equatorial plane at the final time of the evolution for the corresponding simulation. In addition, the lowest area of each panel shows the topology of the lapse function at the final time of the corresponding evolution. Despite the small difference of the initial neutron-star masses, the end product of each merger is remarkably different. In the former case, a hypermassive neutron star is formed, supported against collapse by centrifugal forces. In the latter case, the end product of this particular merger is the delayed formation of a rotating black hole. Animations of such simulations are available at .
In the formation of a hypermassive neutron star with an ellipsoidal figure, quasiperiodic gravitational waves of large amplitude are emitted, with dimensionless strains of about at a distance of 50 Mpc . While these figures are within detectability for advanced laser-interferometric gravitational-wave detectors, the frequency of those waves, between 3 and 4 kHz, may, however, be too high. The eventual detection of such signals might help to further constrain the neutron star EOS.
A large number of neutron-star–binary merger simulations with an emphasis on the the black-hole formation case and on the resulting mass of the surrounding disk was presented by . As in the simulations of , hybrid EOS were adopted to mimic realistic stiff nuclear EOS. In order to determine the mass of the disks that form, the simulations must be continued after black-hole formation, which  achieve with a simple excision technique. This technique, however, allows for long enough computations to even extract the ring-down gravitational waves associated with black-hole quasinormal modes. The resulting frequencies (too large) and amplitudes (too small at 50 Mpc) make unlikely the detection of this part of the signal by advanced laser interferometric detectors. One of the main results of the simulations of  was to find that the disk mass steeply increases with decreasing mass ratio for given ADM mass and EOS, suggesting that such small mass-ratio mergers are good candidates for producing the central engine of short-duration GRBs. Regarding the end product of the simulations, we note that the results of  do not agree with those of recent similar simulations of , at least for high-mass binaries. The reasons may be due to the different choice of nuclear EOS by both groups and also, perhaps more importantly, to the suppression of radiation-reaction effects on the gravitational waves (which rapidly dissipate angular momentum of the merged object) in the (approximate) CFC simulations of .
It is also worth noting the new neutron-star–binary merger simulations performed by [13, 12], particularly the latter, for the innovative aspect of dealing with magnetized stars for the first time (which are, however, described by simple polytropes). The single simulation of a magnetized neutron star merger reported shows the feasibility of the numerical approach, leading to a strongly differentially-rotating object for the particular initial conditions used. Differences in the gravitational waveforms are found when comparing with a purely hydrodynamic merger. Those are motivated by the extremely large value of the magnetic field considered, which has a maximum magnitude of . In this respect we point out that Newtonian SPH simulations of magnetized neutron-star mergers carried out by  have revealed the amplification of the magnetic field beyond magnetar field strength. Such ultrastrong magnetic fields are caused by Kelvin-Helmholtz instabilities in the shear layer that forms between the neutron stars and lasts for a time scale of only 1 ms.
In recent years the first numerical simulations of the mixed binary system, i.e., that formed by a neutron star and a black hole, have also been presented [228, 121, 122, 375, 395, 376, 113]. Future simulations of such systems will benefit from the advances accomplished for the black-hole–binary case, and many are already adopting the same numerical treatment to achieve long-term stability in black-hole spacetimes, namely, the moving puncture approach (see  and references therein).
The simulations of  considered only the headon collision case, but thanks to the use of six levels of fixed mesh refinement, the evolutions were carried over up to times far beyond the collision, which provided time for the extraction of the gravitational wave signal and for estimating the radiated energy.
Fully three-dimensional simulations of the mixed binary system were first performed by . However, two important simplifying assumptions were made in this study: first, the black hole was forced to remain fixed in space by making its mass sufficiently large compared to that of the neutron star, and, secondly, the treatment of the gravitational field equations was approximate due to the use of the conformal-flatness approximation. The first assumption (extreme mass-ratio case) places the onset of tidal disruption of the neutron star outside the ISCO. The initial data are quasiequilibrium models for synchronized neutron-star polytropes generated as solutions of the conformal thin-sandwich decomposition of the Einstein field equations (using the LORENE code ), from which relaxed SPH configurations were built and evolved using an SPH code. The main result of this work was to show that analytic models of mass transfer are not valid for neutron stars with low compactness, as they are highly unstable. The evolutions led to the formation of accretion disks around the black hole, which were discussed in  in the context of the mechanism that triggers short-hard GRBs.
The first fully self-consistent simulations of mixed compact binaries have been performed by Shibata and Uryu [375, 376] and Shibata and Taniguchi . These works use quasicircular states as initial conditions for the simulations, as the timescale of gravitational radiation reaction is a few times longer than the orbital period, and the black hole is modeled by a nonspinning moving puncture (see  and references therein for details). Furthermore, the first two works assume a corotating velocity field for the neutron star while the latter also considers the irrotational case. A -law equation of state with is used in all studies. High-order central schemes are used for the hydrodynamics and the Einstein equations are formulated and solved in the BSSN approach. The single simulation presented in  showed that a black-hole–plus–massive-disk () system is not formed from nonspinning black-hole–neutron-star binaries, later confirmed in  for a slightly larger sample of models. A systematic study of black-hole–neutron-star–binary mergers was later reported in , focusing on the case that the neutron star is tidally disrupted. For all initial conditions, the neutron star was found to be tidally disrupted near the ISCO. The larger the size of the neutron star the more massive the resulting torus. The largest value found was for the case of a black hole of mass and a neutron star of mass .
The most recent simulations of mixed compact binaries are those of , and are also based on moving punctures. In this work, the assumptions adopted by the authors in their early work  are relaxed and the simulations are fully relativistic and self-consistent and do not only cover the extreme mass ratio case. They differ from those of [375, 376, 371] mainly in the way the initial data are built. Here, instead of quasicircular states, the conformal thin-sandwich formalism is used (see  for details). Other than this, the same BSSN formulation is used for the evolution of the gravitational field together with a conservative system for the hydrodynamics equations (that of ), which are solved using the HLL approximate-Riemann solver. The coupling between both sets of equations is based on the Cactus parallelization framework . Despite different codes and initial data being used, the results of the simulations of  are reassuringly similar to those of [375, 376, 371]; the disks that form have masses, which may be too small to produce the required total energy output of a soft-hard GRB.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.