In particular, the understanding of the magnetospheres of neutron stars has recently improved through timedependent numerical simulations, which allow one to investigate the stability of steadystate analytic solutions of the pulsar equation. Existing attempts have used both ideal relativistic magnetohydrodynamics and resistive forcefree electrodynamics (see [203] and references therein). The investigation by [203] permits one to rule out some stationary solutions, which are not naturally recovered by the timedependent simulation. Further work is needed to account for resistive relativistic MHD computations.
Another remarkable investigation based on fullyrelativistic MHD simulations by [205] 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 postshock flow, very different from the steady quasiradial outflow assumed in earlier (sphericallysymmetric) analytical models for plerions. In particular, the jettorus 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 Xray images of the MHD numerical data.
Generalrelativistic MHD winds from rotating neutron stars have also been studied by [66]. Magneticallydominated winds from stars are central to the angular momentum evolution of these objects, and GRMHD studies, such as those carried out by [66], may help us understand the observed spindown of rotationpowered 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 [388] 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 gravitationalwave 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 testbed computations for multidimensional codes. Representative examples are the simulations by [157], with pseudospectral methods, and by [339] with HRSC schemes. These investigations adopted radialgauge polarslicing coordinates in which the general relativistic equations are expressed in a simple way that resembles Newtonian hydrodynamics. Gourgoulhon [157] used a numerical code to detect, dynamically, the zero value of the fundamental mode of a neutron star against radial oscillations. Romero et al. [339] 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 [383] a method based on the nonlinear evolution of deviations from a background stationaryequilibrium 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 higherorder 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 highenergy 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 twodimensional eigenfunctions can be obtained through a Fourier transform technique.
As a first step towards the study of pulsations of rapidlyrotating relativistic stars, Font, Stergioulas, and Kokkotas [138] developed an axisymmetric numerical code called ToniK (see Table 1) that integrates the GRHD equations in a fixedbackground spacetime. The finitedifference code is based on a stateoftheart approximate Riemann solver [102] and incorporates different second and thirdorder TVD and ENO numerical schemes. This code can accurately evolve rapidly rotating stars for many rotational periods, even for stars at the massshedding limit. The test simulations reported in [138] show that, for nonrotating stars, small amplitude oscillations have frequencies that agree to better than 1% with linear normalmode 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 [378], both in Cowling and in fullycoupled evolutions. Contrary to the 2+1 approach followed by [138], the code used in [378] evolves the relativistic stars on null spacetime foliations (see Section 2.2.2).
Until very recently (see below), the quasiradial modes of rotating relativistic stars had been studied only under simplifying assumptions, such as in the slowrotation approximation or in the relativistic Cowling approximation. An example of the latter is presented in [130], where a comprehensive study of all loworder and axisymmetric modes of uniformly and rapidlyrotating relativistic stars was presented, using the code of [138]. This was done for a sequence of appropriatelyperturbed stationary rotating stars, from the nonrotating limit to the massshedding 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 modefrequencies computed with ToniK was carried out by [188] with his pizaa code (see Table 1), which offers an improved treatment of quasistationary scenarios. The agreement found was satisfactory.
In [131], the first mode frequencies of uniformlyrotating stars in full general relativity and rapid rotation were obtained, using the threedimensional 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 [130], 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 [130].
Smallamplitude, nonlinear pulsations of both uniformly and differentially rotating neutron stars were studied by [389] 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 massshedding 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 quasiradial mode is split, at least in the Cowling approximation and mainly in differentially rotating stars, into two different sequences. This study was revisited by [99] using the CoCoNuT code, which solves the GRHD equations for a conformallyflat threemetric (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 massshedding limit, the existence of the masssheddinginduced damping of pulsations was confirmed, even though the effect is not as strong as was previously found in the Cowling approximation [389].
Neutron stars following a core collapse supernova are rotating at birth and can be subject to various nonaxisymmetric instabilities (see, e.g., [388] 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 barmode ( 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 newlyborn rapidlyrotating neutron stars.
Newtonian simulations of the barmode 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 [358] 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 barmode instability have been performed by [29] using the whisky code and extending the set of models considered by [358]. 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 modecoupling 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 corecollapse progenitors. As shown by [384] 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, stateoftheart computations of the evolution of massive stars, which include angular momentum redistribution by magnetic torques and spin estimates of neutron stars at birth [172], lead to corecollapse progenitors, which do not seem to rotate fast enough to guarantee the unambiguous growth of the canonical barmode instability.
Regarding the CFS instability, relativistic hydrodynamic simulations of nonlinear modes are presented in [390], using a version of the gr_astro code (see also [224] for Newtonian simulations). The gravitational radiation reactiondriven 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 [26] 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 [388] for a complete list of references on the subject). Energy leak of the mode into other fluid modes was also considered by [160] through Newtonian hydrodynamic simulations, finding a catastrophic decay of the amplitude only once it has grown to a value larger than that reported by [26]. Detailed analysis showed that this breakdown is due to nonlinear 3mode 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 differentiallyrotating 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 differentiallyrotating relativistic star have been studied by [114]. (Earlier works in the literature dealt with idealized, incompressible, uniformdensity stars.) In order to overcome the long evolution (Alfvén) timescale involved, of the order of for the models simulated, [114] adopted a perturbative metric approach and considered the limit of slow rotation and weak magnetic fields (still astrophysically realistic). This approach, however, permitted [114] 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 [114] 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 blackhole–binary problem has seen major breakthroughs in recent years, to the point that longterm 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 [323] for a review). In this section we focus on simulations of neutronstarbinaries and neutronstar–blackhole configurations, for which progress has been equally dramatic.
Neutronstar binaries are among the most promising sources of gravitational radiation to be detected by the various groundbased interferometers worldwide. The computation of the gravitational waveform during the most dynamic phase of the coalescence and plunge depends crucially on hydrodynamic, finitesize effects. This phase begins once the stars, initially in quasiequilibrium 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, fullyrelativistic hydrodynamic calculations become necessary.
Duez et al. [104] carried out the analytic modelling of the inspiral of corotational and irrotational relativistic neutronstar binaries as a sequence of quasiequilibrium 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 fullyrelativistic hydrodynamic simulations. Detailed theoretical waveforms of the inspiral and plunge similar to those reported by [104] are crucial to enhance the chances of experimental detection in conjunction with matchedfiltering techniques.
The accurate simulation of a neutronstar–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 postNewtonian) framework (see [329] for a review). Many of these studies employ Lagrangian particle methods such as SPH, and only a few have considered (less viscous) highorder finitedifference methods such as PPM [342]. Notwithstanding the difficulties, major progress has been achieved recently in simulating neutronstar–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 conformallyflat 3metric. These early simulations revealed the unexpected appearance of a “binaryinduced 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 [242], after Flanagan [124] pointed out an error in the momentumconstraint equation as implemented by Wilson et al. [416, 418]. A summary of this controversy can be found in [329]. 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 neutronstar binaries modelled by relativistic polytropes with , and . Not only the system dynamics was studied but also gravitationalwave signals and luminosities were computed. The second investigation improved the microphysical treatment of the neutronstar models by using several nonzero temperature EOS, such as those of Shen et al. [351] and Lattimer–Swesty [214]. The results were compared with earlier investigations employing the BSSN equations performed by [375] (see below). It was found that the dynamics and outcome of the models depends strongly on the EOS and on the binary parameters (neutronstar masses and spins). Larger torus masses are found for asymmetric systems (up to for a neutronstar 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 [298], where the authors concluded that the peak frequency of the postmerger 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 neutronstar–binary coalescence in general relativity (see, e.g., [280]). The group developed a threedimensional code that solves the full set of Einstein and hydrodynamics equations through finite differences in a uniform Cartesian grid, using van Leer’s scheme [403] with TVD flux limiters. Shockwaves are spread out using a tensor artificialviscosity 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 [386]). A few results of neutronstar–binary simulations including gravitational radiation are reported in [304].
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. [261], who performed timedependent 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 [349] it was argued that in a headon collision event, sufficient thermal pressure is generated to support the remnant in quasistatic equilibrium against (prompt) collapse, prior to slow cooling via neutrino emission (delayed collapse). In [261], 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 restmass density and the internal energy evolution during the onaxis 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 moderatelyrelativistic shockwave (Lorentz factor ) appearing at (code units; ; yellowwhite), 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 [187] in the headon collison of two neutron stars.
Using the code of [109], Marronetti et al. [234] performed the first dynamic determination in relativistic gravity of the ISCO of a neutronstar 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 longterm 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 [374] 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 gravitationalradiation 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 3metric and the existence of a helicoidal Killing vector (see [374] 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 parameterspace 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 marginallystable massive neutron star forms. In the latter case, the star is supported against selfgravity 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 [374]. Therefore, future detection of highfrequency gravitational waves could help to constrain the maximum allowed mass of neutron stars. In addition, for prompt blackhole 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 [372] 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 neutronstar models, the simulations of [372] also considered mergers of unequalmass neutron stars. The underlying motivation was to find out whether massive accretion disks form depending on the initial mass ratio. Disk formation during neutronstar–binary coalescence, a fundamental issue for cosmological models of short duration GRBs, was indeed found to be enhanced for unequalmass 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 equalmass neutronstar–binary mergers. The left panel is a merger, while the right panel depicts the case. Likewise, both panels display the velocity field and restmass 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 neutronstar 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 [352].
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 [357]. While these figures are within detectability for advanced laserinterferometric gravitationalwave 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 neutronstar–binary merger simulations with an emphasis on the the blackhole formation case and on the resulting mass of the surrounding disk was presented by [370]. As in the simulations of [372], 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 blackhole formation, which [370] achieve with a simple excision technique. This technique, however, allows for long enough computations to even extract the ringdown gravitational waves associated with blackhole 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 [370] was to find that the disk mass steeply increases with decreasing mass ratio for given ADM mass and EOS, suggesting that such small massratio mergers are good candidates for producing the central engine of shortduration GRBs. Regarding the end product of the simulations, we note that the results of [370] do not agree with those of recent similar simulations of [297], at least for highmass 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 radiationreaction effects on the gravitational waves (which rapidly dissipate angular momentum of the merged object) in the (approximate) CFC simulations of [297].
It is also worth noting the new neutronstar–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 differentiallyrotating 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 neutronstar mergers carried out by [324] have revealed the amplification of the magnetic field beyond magnetar field strength. Such ultrastrong magnetic fields are caused by KelvinHelmholtz 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 blackhole–binary case, and many are already adopting the same numerical treatment to achieve longterm stability in blackhole spacetimes, namely, the moving puncture approach (see [323] and references therein).
The simulations of [228] 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 threedimensional simulations of the mixed binary system were first performed by [122]. 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 conformalflatness approximation. The first assumption (extreme massratio case) places the onset of tidal disruption of the neutron star outside the ISCO. The initial data are quasiequilibrium models for synchronized neutronstar polytropes generated as solutions of the conformal thinsandwich decomposition of the Einstein field equations (using the LORENE code [227]), 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 [121] in the context of the mechanism that triggers shorthard GRBs.
The first fully selfconsistent simulations of mixed compact binaries have been performed by Shibata and Uryu [375, 376] and Shibata and Taniguchi [371]. 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 [323] 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. Highorder central schemes are used for the hydrodynamics and the Einstein equations are formulated and solved in the BSSN approach. The single simulation presented in [375] showed that a blackhole–plus–massivedisk () system is not formed from nonspinning blackhole–neutronstar binaries, later confirmed in [376] for a slightly larger sample of models. A systematic study of blackhole–neutronstar–binary mergers was later reported in [371], 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 [113], and are also based on moving punctures. In this work, the assumptions adopted by the authors in their early work [122] are relaxed and the simulations are fully relativistic and selfconsistent 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 thinsandwich formalism is used (see [395] 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 [108]), which are solved using the HLL approximateRiemann solver. The coupling between both sets of equations is based on the Cactus parallelization framework [244]. Despite different codes and initial data being used, the results of the simulations of [113] 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 softhard GRB.
http://www.livingreviews.org/lrr20087 
This work is licensed under a Creative Commons AttributionNoncommercialNo Derivative Works 2.0 Germany License. Problems/comments to 