Go to previous page Go up Go to next page

5.3 Hydrodynamic and magnetohydrodynamic evolution of neutron stars

The numerical investigation of many interesting astrophysical processes involving neutron stars, such as the rotational evolution of proto-neutron stars (which can be affected by a dynamic bar-mode instability and by the Chandrasekhar–Friedman–Schutz instability), the appearance of MHD winds and outflows from pulsar magnetospheres, or the gravitational radiation from unstable pulsation modes or, more importantly, from the catastrophic coalescence and merger of neutron-star compact binaries, requires the ability of accurate, long-term hydrodynamic and MHD evolutions employing relativistic gravity. These scenarios are receiving increasing attention in recent years, as we discuss next.

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 [203Jump To The Next Citation Point] and references therein). The investigation by [203] 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 [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 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 [66Jump To The Next Citation Point]. Magnetically-dominated 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 spin-down of rotation-powered pulsars with strong magnetic fields.

5.3.1 Pulsations and instabilities of relativistic stars

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 [388Jump To The Next Citation Point] 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 [131Jump To The Next Citation Point355364389Jump To The Next Citation Point99Jump To The Next Citation Point]).

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 [157Jump To The Next Citation Point], with pseudo-spectral methods, and by [339Jump To The Next Citation Point] 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 [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 1.94532 M ⊙ is stable and a model of 1.94518 M ⊙ is unstable. More recently, in [383] 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 [138Jump To The Next Citation Point130Jump To The Next Citation Point131Jump To The Next Citation Point390Jump To The Next Citation Point389Jump To The Next Citation Point99Jump To The Next Citation Point]. 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 [138Jump To The Next Citation Point] 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 [102] 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 [138Jump To The Next Citation Point] 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 [378Jump To The Next Citation Point], both in Cowling and in fully-coupled evolutions. Contrary to the 2+1 approach followed by [138Jump To The Next Citation Point], the code used in [378] 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 [130Jump To The Next Citation Point], where a comprehensive study of all low-order ℓ = 0,1,2 and 3 axisymmetric modes of uniformly and rapidly-rotating relativistic stars was presented, using the code of [138]. 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 [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 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 [130Jump To The Next Citation Point], 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].

Small-amplitude, nonlinear pulsations of both uniformly and differentially rotating neutron stars were studied by [389Jump To The Next Citation Point] 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 [99] 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 l = 0 and l = 4 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 [389].

Neutron stars following a core collapse supernova are rotating at birth and can be subject to various nonaxisymmetric instabilities (see, e.g., [388Jump To The Next Citation Point] for a review). Among those, if the rotation rate is high enough that the ratio of rotational kinetic energy T to gravitational potential energy W, β ≡ T∕|W |, exceeds the critical value β ∼ 0.27 d inferred from studies with incompressible Maclaurin spheroids, the star is subject to a dynamic bar-mode (l = m = 2 f-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 βd ∼ 0.27 and is quite independent of the stiffness of the EOS provided the star is not strongly differentially rotating. The relativistic simulations of [358Jump To The Next Citation Point] yield a value of β ∼ 0.24 − 0.25 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 βd can be as low as 𝒪 (0.01).

Get Flash to see this player.

Figure 16: mpeg-Movie (8995 KB) The time evolution of the bar-mode instability of one of the differentially rotating neutron-star models of [29Jump To The Next Citation Point]. The animation depicts the distribution of the rest-mass density. The exponential growth of the instability, which is followed by its saturation, the development of spiral arms, and the attenuation of the bar deformation, is all clearly visible in the movie. The last part of the dynamics is dominated by an almost axisymmetric configuration. Visualization developed at SISSA. Used with permission.

Further relativistic hydrodynamic simulations of the dynamic bar-mode 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 βd 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 16Watch/download Movie.

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 [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, 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 [172], 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 r-modes are presented in [390Jump To The Next Citation Point], using a version of the gr_astro code (see also [224Jump To The Next Citation Point] for Newtonian simulations). The gravitational radiation reaction-driven instability of the r-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 r-mode energy into some short wavelength oscillation modes (see [26Jump To The Next Citation Point] and references therein). The latter mechanism would dramatically limit the r-mode amplitude to a small value of 𝒪 (10−3), much smaller than those found in the early simulations of [390224] for isentropic stars (see [388] for a complete list of references on the subject). Energy leak of the r-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 3-mode coupling of the r-mode to two other inertial modes. The saturation amplitude is of −2 𝒪 (10 ). It is believed that saturation amplitudes of this order are still appropriate to detect gravitational waves from the r-mode instability in LMXBs.

As mentioned in Section 5.1.2 the axisymmetric collapse of differentially-rotating magnetized neutron stars has been studied by [105360] in full general relativity. The effects of magnetic fields on the magnetic braking of a magnetized differentially-rotating relativistic star have been studied by [114Jump To The Next Citation Point]. (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 104 M for the models simulated, [114Jump To The Next Citation Point] adopted a perturbative metric approach and considered the limit of slow rotation and weak magnetic fields (still astrophysically realistic). This approach, however, permitted [114Jump To The Next Citation Point] 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.

5.3.2 Neutron-star binary and black-hole–neutron-star coalescence

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 [323Jump To The Next Citation Point] 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 8 10 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 M ∕R ≪ 1 and M ∕r ≪ 1 simultaneously, M being the total mass of the binary, R the neutron star radius, and r 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. [104Jump To The Next Citation Point] 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 [104] 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 [329Jump To The Next Citation Point] 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 [342]. 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 [416Jump To The Next Citation Point418Jump To The Next Citation Point242Jump To The Next Citation Point]. 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 [242], after Flanagan [124] pointed out an error in the momentum-constraint equation as implemented by Wilson et al. [416418]. 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 [300299]. 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 Γ = 2.0,2.6, and 3.0. 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. [351] and Lattimer–Swesty [214]. The results were compared with earlier investigations employing the BSSN equations performed by [375Jump To The Next Citation Point] (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 0.3M ⊙ 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 [298], 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., [280]). 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 [403] 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 [386]). A few results of neutron-star–binary simulations including gravitational radiation are reported in [304].

Get Flash to see this player.

Figure 17: mpeg-Movie (558 KB) An animation of a headon collision simulation of two 1.4 M ⊙ neutron stars obtained with a relativistic code [137261Jump To The Next Citation Point]. The movie shows the evolution of the density and internal energy. The formation of the black hole in prompt timescales is demonstrated by the sudden appearance of the apparent horizon at t = 0.16 ms (t = 63.194 in code units), which is indicated by violet dotted circles representing the trapped photons. See [245] for download options of higher-quality versions of this movie.

The headon collision of two neutron stars (a limiting case of a coalescence event) was considered by Miller et al. [261Jump To The Next Citation Point], 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 [349] 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 [261], prompt collapse to a black hole was found in the headon collision of two 1.4M ⊙ neutron stars modeled by a polytropic EOS with Γ = 2. The stars, initially separated by a proper distance of d = 44 km, were boosted toward one another at a speed of ∘ ------- GM ∕d (the Newtonian infall velocity). The simulation employed a Cartesian grid of 1923 points. The time evolution of this simulation can be followed in the Quicktime movie in Figure 17Watch/download Movie. 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 t = 0.16 ms (t = 63.194 in code units). The violet dotted circles indicate the trapped photons. The animation also shows a moderately-relativistic shockwave (Lorentz factor ∼ 1.2) appearing at t ∼ 36 (code units; ∼ 0.09 ms; yellow-white), which eventually is followed by two opposite moving shocks (along the infalling z 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 neutron-star binary (modelled as Γ = 2 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. [353Jump To The Next Citation Point373Jump To The Next Citation Point374Jump To The Next Citation Point357Jump To The Next Citation Point372Jump To The Next Citation Point370Jump To The Next Citation Point]. 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 [374Jump To The Next Citation Point] accounted for a large sample of parameters of the binary system, such as the compactness of the (equal mass) neutron stars (0.12 ≤ M ∕R ≤ 0.16), the adiabatic index of the Γ-law EOS (1.8 ≤ Γ ≤ 2.5), 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 [374Jump To The Next Citation Point] 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%.

View Image

Figure 18: Isodensity contours in the equatorial plane for the merger of a 1.3M ⊙ − 1.3M ⊙ neutron-star binary (left panel) and a 1.35M ⊙ − 1.35 M ⊙ case. Vectors indicate the local velocity field (vx,vy). The image corresponds to the final snapshot of the evolution for the two models. The outcome of the merger for the more massive case is the delayed formation of a black hole, as signalled by the collapse of the lapse function shown in the lower panel of the image. In the less massive case, a long-lived hypermassive neutron star supported by rotation is formed. Simulations performed by [372Jump To The Next Citation Point]. (Used with permission.)

The comprehensive parameter-space survey carried out by [353373374Jump To The Next Citation Point] 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 [374]. 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 [372Jump To The Next Citation Point] where hybrid EOS were adopted to mimic realistic nuclear equations of state. In this approach the EOS is divided into two parts, P = Pcold + Pth, where Pcold 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 ∼ 2.04M ⊙ and 1.80M ⊙, respectively. Apart from improving the neutron-star models, the simulations of [372Jump To The Next Citation Point] 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 18View Image. Both panels in this figure correspond to equal-mass neutron-star–binary mergers. The left panel is a 1.3M ⊙ − 1.3 M ⊙ merger, while the right panel depicts the 1.35M − 1.35M ⊙ ⊙ 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 [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 −21 h ∼ 6 − 7 × 10 at a distance of 50 Mpc [357]. 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 [370Jump To The Next Citation Point]. 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 black-hole formation, which [370Jump To The Next Citation Point] 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 [370Jump To The Next Citation Point] 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 [370] do not agree with those of recent similar simulations of [297Jump To The Next Citation Point], 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 [297].

It is also worth noting the new neutron-star–binary merger simulations performed by [1312], 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 ∼ 1016 G. In this respect we point out that Newtonian SPH simulations of magnetized neutron-star mergers carried out by [324] 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 [228Jump To The Next Citation Point121Jump To The Next Citation Point122Jump To The Next Citation Point375Jump To The Next Citation Point395Jump To The Next Citation Point376Jump To The Next Citation Point113Jump To The Next Citation Point]. 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 [323Jump To The Next Citation Point] 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 three-dimensional simulations of the mixed binary system were first performed by [122Jump To The Next Citation Point]. 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 [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 short-hard GRBs.

The first fully self-consistent simulations of mixed compact binaries have been performed by Shibata and Uryu [375Jump To The Next Citation Point376Jump To The Next Citation Point] and Shibata and Taniguchi [371Jump To The Next Citation Point]. 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 Γ = 2 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 [375Jump To The Next Citation Point] showed that a black-hole–plus–massive-disk (∼ 1 M ⊙) system is not formed from nonspinning black-hole–neutron-star binaries, later confirmed in [376Jump To The Next Citation Point] for a slightly larger sample of models. A systematic study of black-hole–neutron-star–binary mergers was later reported in [371Jump To The Next Citation Point], 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 ∼ 0.16M ⊙ for the case of a black hole of mass ∼ 4M ⊙ and a neutron star of mass ∼ 14.7 km.

The most recent simulations of mixed compact binaries are those of [113Jump To The Next Citation Point], 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 self-consistent and do not only cover the extreme mass ratio case. They differ from those of [375Jump To The Next Citation Point376Jump To The Next Citation Point371Jump To The Next Citation Point] mainly in the way the initial data are built. Here, instead of quasicircular states, the conformal thin-sandwich 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 approximate-Riemann 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 [375376371]; the disks that form have masses, which may be too small to produce the required total energy output of a soft-hard GRB.

  Go to previous page Go up Go to next page