{'fraction_non_alphanumeric': 0.032586584102045776, 'fraction_numerical': 0.015915865968961503, 'mean_word_length': 4.463963649073751, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accurate evolutions of unequal-mass neutron-star binaries: properties of the torus and short GRB engines

We present new results from accurate and fully general-relativistic simulations of the coalescence of unmagnetized binary neutron stars with various mass ratios. The evolution of the stars is followed through the inspiral phase, the merger and prompt collapse to a black hole, up until the appearance of a thick accretion disk, which is studied as it enters and remains in a regime of quasi-steady accretion. Although a simple ideal-fluid equation of state with \Gamma=2 is used, this work presents a systematic study within a fully general relativistic framework of the properties of the resulting black-hole--torus system produced by the merger of unequal-mass binaries. More specifically, we show that: (1) The mass of the torus increases considerably with the mass asymmetry and equal-mass binaries do not produce significant tori if they have a total baryonic mass M_tot>~ 3.7 M_sun; (2) Tori with masses M_tor ~ 0.2 M_sun are measured for binaries with M_tot ~ 3.4 M_sun and mass ratios q ~ 0.75-0.85; (3) The mass of the torus can be estimated by the simple expression M_tor(q, M_tot) = [c_1 (1-q) + c_2](M_max-M_tot), involving the maximum mass for the binaries and coefficients constrained from the simulations, and suggesting that the tori can have masses as large as M_tor ~ 0.35 M_sun for M_tot ~ 2.8 M_sun and q ~ 0.75-0.85; (4) Using a novel technique to analyze the evolution of the tori we find no evidence for the onset of non-axisymmetric instabilities and that very little, if any, of their mass is unbound; (5) Finally, for all the binaries considered we compute the complete gravitational waveforms and the recoils imparted to the black holes, discussing the prospects of detection of these sources for a number of present and future detectors.


Introduction
The numerical investigation of the coalescence and merger of binary neutron stars (NSs) within the framework of general relativity is receiving increasing attention in recent years (see e.g. [1,2,3,4,5,6,7,8] and references therein). Drastic improvements in the simulation front regarding mathematics (e.g. formulation of the equations), physics (e.g. incorporation of equations of state (EOSs) from nuclear physics), and numerical methods (e.g. use of high-resolution methods and adaptive mesh refinement), along with increased computational resources, have allowed to extend the scope of the early simulations (e.g. [9]). Larger initial separations have recently started being considered and some of the existing simulations have expanded the range spanned by the models well beyond black-hole (BH) formation [3,8,7]. This is allowing the computation of the entire gravitational waveform from the early inspiral to the decaying tail of the late ringing of the formed BH. The construction of such waveform templates is still one of the driving motivations to perform binary NS simulations, as such events are among the most promising sources of detectable gravitational radiation for laser interferometric detectors. The current estimate for the detection rate relative to the firstgeneration interferometric detectors is ∼ 1 event per ∼ 40 − 300 years, increasing to an encouraging ∼ 10 − 100 events per year for the advanced detectors [10]. The second major incentive behind this type of simulations is establishing whether the end-product of the merger can act as the underlying mechanism operating at the central engine of short-hard gamma-ray bursts (SGRBs) [11,12]. The consensus emerging from the existing simulations suggests the formation, depending on the suitability of the initial parameters of the simulated model, of a BH of stellar mass surrounded by a hot disk. Driven by neutrino processes and magnetic fields, such a compact system may be capable of launching a relativistic fireball with an energy of ∼ 10 48 erg on a timescale of 0.1 − 1 s [13].
This paper is dedicated in particular to investigating the late-time dynamics of the torus formed after the merger of unequal-mass NS binaries. As we describe below, all but one model of our initial sample have mass ratio different from one. Previous simulations have shown that the key parameter controlling the amount of mass left in the disk for a given initial mass in the system and EOS is the NS mass ratio [14,1]. Broadly speaking the general trend is simple: the larger the departure from equal-mass ratio, the more important tidal effects become in the less massive star, resulting in its tidal disruption. Because this takes place when the separation is still comparatively large, the angular momentum of the matter is still large and it results in larger-size and more massive disks. Early and low-resolution simulations with an ideal-gas EOS [14] have been shown to yield a disk mass of several percents of the total mass of the system for a mass ratio of ∼ 0.85. Improved simulations by [1] which adopted a hybrid EOS to mimic realistic, stiff nuclear EOS, indicate that the mass of the disk is ∼ 0.01 M ⊙ or slightly larger when the merger does not result in prompt collapse to a BH but in the formation of a hypermassive NS of large ellipticity (which later collapses to a BH following angularmomentum transport caused by emission of gravitational radiation). Similar disk masses, as large as ∼ 0.02 M ⊙ , are also reported in the latest simulations of [8], in which the initial orbital separation of the two stars is larger than in previous works.
Observational data seems to indicate that the total gravitational masses of the known galactic NS binary systems are in a narrow range ∼ 2.65 − 2.85 M ⊙ [15] and there is also evidence indicating that the masses of the two NSs are nearly equal, with the baryonic mass ratio q ≡ M 1 /M 2 being between 1 and ∼ 0.7 (Hereafter we will refer to q simply as the "mass ratio" but we report in table 1 also the ratio in the ADM masses q ADM ; q and q ADM do not coincide because of the nonlinear relation between baryonic and gravitational mass). However, there is no theoretical reason to assume that unequal-mass NS binaries are not produced in nature as often as the seemingly prevailing almost exactly equal-mass systems, particularly for twin giant progenitors [16,17]. There are, indeed, recent computations and observations that contradict the predominance of almost exactly equal-mass systems [18,17]. On the one hand, binary population synthesis computations performed by [18] show two peaks in the observability-weighted distribution of double NSs. One of these peaks is around q ∼ 1 and appears when both masses are close to 1.4 M ⊙ . The second peak is around considerably smaller mass ratios and depends on the assumed maximum mass of a single NS (which is in turn dependent on the EOS considered): the higher this mass (the stiffer the EOS) the more significant the second peak [18,19]. However, the crucial parameter determining the shape of the distribution is the inclusion of hypercritical accretion onto the compact object during quantitative measures of the accuracy of our numerical methods. We use a spacelike signature (−, +, +, +) and a system of units in which c = G = M ⊙ = 1 (or in other units whenever more convenient). Greek indices are taken to run from 0 to 3, Latin indices from 1 to 3 and we adopt the standard convention for the summation over repeated indices.

Mathematical and Numerical Setup
All the details on the mathematical and numerical setup used for producing the results presented here are discussed in depth in [32,3]. In what follows, we limit ourselves to a brief overview.

Einstein and Hydrodynamics equations
The evolution of the spacetime is obtained using the CCATIE code, a three-dimensional finite-differencing code providing a solution of a conformal traceless formulation of the Einstein equations with a "1 + log" slicing condition and a "Gamma-driver" shift condition (the interested reader is addressed to [32,33] for a detailed discussion of the equations and gauges used). The general-relativistic hydrodynamics equations are instead solved using the Whisky code presented in [34,35,33], which adopts a flux-conservative formulation of the equations as presented in [36] and high-resolution shock-capturing schemes. The Whisky code implements several reconstruction methods, such as Total-Variation-Diminishing (TVD) methods, Essentially-Non-Oscillatory (ENO) methods [37], and the Piecewise Parabolic Method (PPM) [38]. Also, a variety of approximate Riemann solvers can be used, starting from the Harten-Lax-van Leer-Einfeldt (HLLE) solver [39], over to the Roe solver [40], and the Marquina flux formula [41] (see [34,35] for a more detailed discussion). All the results reported hereafter have been computed using the Marquina flux formula [42] and a PPM reconstruction. We stress again (as already done in [3,7]) that the use of high-order methods and high-resolution is essential to be able to draw robust conclusions on the inspiral and merger. Lower-order methods in the reconstruction and low resolution may yield convergent and apparently reasonable results which however contain a large truncation error. Specific examples of this type of problem are presented in Appendix A of [3] and in Figure 4 of [7]. Also, a measure of our overall accuracy is presented in Appendix A below and shows that by employing such methods we are able to conserve energy and angular momentum to ∼ 1% over a timescale of ∼ 140 ms.
The system of hydrodynamics equations is closed by an EOS. As discussed in detail in [3], the choice of the EOS plays a fundamental role in the post-merger dynamics and significantly influences the survival time, against gravitational collapse, of the hyper-massive neutron star (HMNS) likely produced by the merger. It is therefore important that special attention is paid to use EOSs that are physically motivated, as done in [43] within a conformally flat description of the fields and a simplified treatment of the hydrodynamics. Because we are here mostly concerned with drawing a first qualitative picture of the properties of the torus in a space of parameters that is as vast as computationally affordable, we have employed the commonly used "ideal-fluid" EOS, in which the pressure p is expressed as p = ρ ǫ(Γ − 1), where ρ is the rest-mass density, ǫ is the specific internal energy, and Γ is the adiabatic exponent. Such an EOS, while simple, provides a reasonable approximation and we expect that the use of realistic EOSs would not change the main results of this work. Furthermore, it was shown in [44] that the inspiral of equal-mass binaries of NSs described by realistic EOSs can be reproduced quite well by studying NSs with the same mass and radii but constructed as polytropes with Γ = 2.
As in [3], the gravitational-wave signal is extracted using two methods. The first method uses the Newman-Penrose formalism so that the gravitational-wave polarization amplitudes h + and h × are then related to Ψ 4 by simple time integrals [45] where the double overdot stands for the second-order time derivative and the curvature scalar is defined as a particular component of the Weyl curvature tensor, C αβγδ , projected onto a given null frame {l, n, m,m} (see [32] for details). The second and independent method is instead based on the measurements of the non-spherical gauge-invariant perturbations of a Schwarzschild BH (see refs. [46,47,48] for some applications of this method to Cartesiancoordinate grids). In particular, the gravitational-wave polarization amplitudes are in this case expressed in terms of gauge-invariant metric perturbations [49] h where −2 Y ℓm are the s = −2 spin-weighted spherical harmonics and Q × ℓm , Q + ℓm the (gaugeinvariant) odd-parity (or axial) current multipoles and even-parity (or polar) mass multipoles of the perturbed metric, respectively. In practice, these multipoles are computed on a set of 2-spheres of fixed coordinate radius r iso = 200 (i.e. ∼ 300 km).
The two wave-extraction methods yield results whose differences are smaller than 1%, and hereafter we will concentrate only on the one using the gauge-invariant perturbations as it reduces the number of integration constants to be determined when computing the gravitational-wave strain.

Adaptive Mesh Refinements and Grid setup
The grid hierarchy is handled by the Carpet mesh refinement driver [50]. It implements vertex-centered mesh refinement, also known as the box-in-box method, and allows for regridding during the calculation as well as multiple grid centers. With mesh refinement, a small number of grids with different resolutions, called refinement levels, overlay each other, and are nested in a way that the coarsest grid has the largest extent and the finest grid the smallest extent. While the refined grids in the interior allow for an increased resolution where it is most desired, the outer boundary can at the same time be kept at a large distance.
The timestep on each grid is set by the Courant condition (expressed in terms of the speed of light) and so by the spatial grid resolution for that level; the typical Courant coefficient is set to be 0.35. The time evolution is carried out using 4th-order-accurate Runge-Kutta integration algorithm. Boundary data for finer grids are calculated with spatial prolongation operators employing 3rd-order polynomials and with prolongation in time employing 2ndorder polynomials. The latter allows a significant memory saving, requiring only three timelevels to be stored, with little loss of accuracy due to the long dynamical timescale relative to the typical grid timestep.
For the inspiral phase of the system of binary NSs, two grid centers {r c,i : i = 1, 2} are defined, with one grid center located at the grid point where the rest-mass density reaches its maximum ρ max = max(ρ) and the other grid center located at the π-symmetric point (i.e. the grid point correspondent to ρ max and rotated by 180 degrees around the z-axis). The grid hierarchy is composed of six refinement levels and a 2 : 1 refinement factor for successive levels. Once the condition ρ max = max(ρ max,i ) ≥ 1.2 ρ max,initial is satisfied, which is known from experience to occur during the merger phase and well before collapse, the grid hierarchy is reduced to a single grid center fixed at the origin of the grid. At the initial time, the finest grids cover each star completely. Later, during the merger phase, matter outflows cross the boundary to the second finest grid and subsequently to the other coarser refinement levels. The grid resolution varies from ∆ 6 = 0.15 (i.e. ∼ 221 m) for the finest level to ∆ 1 = 4.8 (i.e. ∼ 7.1 km) for the coarsest level, whose outer boundary is at 240 in our units (i.e. ∼ 360 km). Initially, the number of grid points across the linear dimension of a star is of the order of 100. The torus surrounding the BH after collapse is usually not contained within the finest grid, but its high-density region is covered by the second finest grid with resolution The whole grid is set up to be symmetric with respect to the (x, y) plane for both unequalmass binaries and equal-mass binaries. The boundary conditions are chosen to be radiative for the metric, in order to prevent gravitational waves from scattering back into the grid, and static for the hydrodynamical variables. Note that the above setup is identical to that adopted in [3]. Table 1. Properties of the binary NS initial data. From left to right the columns show: the name of the model (assembled from its rounded total baryonic mass preceded by the letter M and its mass ratio preceded by the letter q), the total baryonic mass Mtot of the system, the total ADM mass M ADM of the system, the ratio of the baryonic masses of the two stars q = M 2 /M 1 , the ratio of the ADM masses of the two stars, the total angular momentum J, the initial orbital frequency ν orb , the initial maximum rest-mass density ρmax, the mean radiusr i of each star, the axis ratioĀ i of each star, the individual ADM mass M ∞ i of each star as considered in isolation at infinity, and the compactness C ∞ i of each star as considered in isolation at infinity. The mean radius is defined asr i ≡ (r ⊢ + r ⊣ + r ⊥ + r pol )/4, where r ⊢ and r ⊣ are the radii of the star parallel to the line connecting the stars, r ⊥ is the radius in the equatorial plane perpendicular to that line, and r pol is the radius perpendicular to the equatorial plane. The axis ratio is defined as the ratio between the mean radius parallel to the line connecting the stars, and the mean radius in the plane perpendicular to that line, namelȳ A i ≡ (r ⊥ + r pol )/(r ⊢ + r ⊣ ). All values except ρmax are provided by the output of the LORENE code, and the accuracy of Mtot and J is the one at which the Whisky code is able to reproduce them for the present setup.

Initial data
We use quasi-equilibrium initial data generated with the multi-domain spectral-method code LORENE developed at the Observatoire de Paris-Meudon [51]. For more information on the code and its methods, the reader is referred to the LORENE web pages [52]. In particular, because the binaries are not expected to be corotating [53], we use irrotational configurations, defined as having vanishing vorticity, and obtained under the additional assumption of a conformally flat spacetime metric [51]. Some of the models investigated in this paper are publicly available on servers of the Meudon group [52]. Others have been created by us specifically for the unequalmass simulations presented here. The models of the lowest mass ratios have been kindly provided by Dorota Gondek-Rosińska). The EOS assumed for the initial data is in all cases the polytropic EOS p ≡ K ρ Γ with an adiabatic index Γ = 2 and a polytropic coefficient K = 0.0332 ρ nuc c 2 /n Γ nuc = 123.6 (in units where c = G = M ⊙ = 1), where ρ nuc and n nuc refer to the nuclear rest-mass and number densities, respectively. For this particular EOS, the allowed maximum baryonic mass for an individual stable NS is ∼ 2.00 M ⊙ . The initial coordinate separation of the stellar centers in all cases is d = 45 km.
The models used as initial data include both equal-mass models and most importantly unequal-mass models. As mentioned in the Introduction, we here concentrate on the dynamics of the massive tori resulting from the merger. The formation of the torus is enhanced for smaller mass ratios. Although the observational evidence is not very firm [15,54] there is also no theoretical reason ruling out their possible existence [18,17]. A full list of all considered models together with a selection of physical quantities defining them, e.g. baryon and ADM mass, orbital frequency and initial angular momentum, etc., is given in table 1. To distinguish simply the different binaries we adopt the following naming convention: any initial data binary is indicated as M%q#, with % being replaced by the rounded total baryonic mass M tot of the binary neutron-star system and # by the mass ratio q. As an example, M3.4q0.80 is the binary with total baryonic mass M tot ≃ 3.4 M ⊙ and mass ratio q = 0.80.

General dynamics
In a previous work [3], we have investigated the dynamics of the coalescence and merger of equal-mass binary NSs for models with total baryonic mass M tot = 2.912 M ⊙ and M tot = 3.250 M ⊙ . It was found that for any of the two EOSs considered, binaries with (initial) total baryonic mass below a certain limit do not collapse promptly to a BH but rather yield an oscillating HMNS, which undergoes delayed collapse to a BH. Independently of the mass ratio, all of the binaries under consideration here have masses higher than those considered in [3] and all collapse promptly never leading to a HMNS even if the EOS used here is a non-isentropic one (see discussion in [3] on the different qualitative behaviour between an isentropic and a non-isentropic EOS). This absence of a HMNS, however, is the direct consequence of the chosen initial data rather than a feature of unequal-mass mergers and it has been here exploited simply to reduce the computational costs by shortening the time of the collapse to a BH. Figure 1 shows a selection of representative isodensity contours on the equatorial plane for the equal-mass binary M3.6q1.00. At the initial time, the stars are in their quasiequilibrium configuration at a coordinate separation of 45 km. The binary progressively speeds up while inspiralling. After slightly more than two orbits have been completed (namely after about 5-6 ms), the stars merge and, about 2-3 ms later, an apparent horizon (which we search with the code of [55]) is found. The ideal-fluid EOS employed in the simulations allows for shock-heating and an increase of the specific internal energy ǫ, as shown in [3]; this, in turn, causes some matter to be ejected from the rotating central object and to propagate into the surrounding artificial atmosphere. The evolution of model M3.6q1.00 shows that matter is ejected in small amounts during the inspiral phase and in larger amounts during the merger phase, when the shocks are much stronger. Therefore, while small spiral arms can certainly be observed in the outer regions during the merger phase (see the last two snapshots of figure 1), they do not have sufficient angular momentum to reach distances as large as in the unequal-  mass models (see discussion below). Instead, the spiral arms wind around the rapidly rotating central object formed by the two NS cores. Quantitative results regarding the BH spin and the mass and angular momentum of the remaining disk will be discussed in subsequent sections.
To contrast the evolution of an unequal-mass binary, figure 2 shows the same selection shows the onset of the merger, the last two frames (at times t = 6.620 ms, t = 7.414 ms) show the behaviour of the system during the collapse to a BH. Note that the computational domain is much larger than what is shown and extends to ∼ 360 km of isodensity contours on the equatorial plane as represented in figure 1, only now for the M3.4q0.70 model, which has the smallest mass ratio considered in this work. The asymmetry of the binary system is already apparent at the initial time. The heavier star is much more compact than its extended less massive companion, which is deformed already at the initial distance by tidal forces. In addition, the center of mass does not coincide with the point halfway between the centers of the stars, but it is shifted toward the more massive star. During the inspiral phase, the heavier and more compact star is only slightly affected by its companion, whereas the latter is decompressed rapidly while being accreted onto the heavier star. This is visible in the three intermediate panels of figure 2. The tidal disruption of the lower-mass NS when it still retains a large fraction of its angular momentum results in an extended tidal tail, which, unlike what happens in the equal-mass case, transfers angular momentum outwards in a much more efficient way. This leads to the formation of large spiral arms extending well beyond the domain shown in figure 2 and ultimately to a more rapid ejection of matter. Gravitationally bound matter travelling along the spiral arms away from the central object will form a more massive accretion torus around the central BH than that formed in the case of an equal-mass, symmetric binary system. It should be noted that, although the rest-mass density of the matter in these spiral arms is much smaller than the central one, it has nevertheless densities ρ 10 10 g/cm −3 and thus well in a general-relativistic regime.

Properties of the black hole
As mentioned above, because of the large initial mass of the system and irrespective of the mass ratio, the merged object rapidly collapses to a BH. Its mass and angular momentum have been computed making use of the dynamical-horizon formalism [56,57], which provides a simple and accurate measure of the BH properties also when this is subject to the inflow of mass and angular momentum [35]. In the case of the equal-mass binary, because the disk resulting from the merger has comparatively small mass, the BH settles rapidly to an approximately stationary configuration, and the mass and spin of the BH measured at formation, i.e. M = 2.56 M ⊙ and a ≡ J/M 2 = 0.745, respectively, do not vary significantly throughout the subsequent evolution of the system. On the other hand, when considering unequal-mass binaries, the mass and the spin of the BH show, on the timescale of the simulations, a variation in time of ∼ 5% and ∼ 2%, respectively, because of the continued and intense accretion of both mass and angular momentum. Table 2 shows the corresponding parameters for all models at the final time of the evolution, which is not the same for the different binaries considered. We note that finding and tracking the apparent horizon in the case of binaries with small mass ratio is far from being simple since the asymmetry in the merger dynamics leads to a noticeable motion of the "center-of-mass" of the system. Hence, the location of the trial surface for the apparent horizon cannot be simply associated to a pre-existing black hole (as in the case of BH binaries [32]) or to a pre-determined coordinate location (as in the case of the collapse of a rotating star [35]). The end-result of this complication is that the apparent horizon could not be tracked successfully in all the models under consideration. This was the case of models M3.4q0.80 and M3.4q0.70, for which it was not possible to measure the mass and spin of the corresponding BH. Furthermore, in the case of the binary M3.5q0.75 the measurements were not made with the dynamical-horizon formalism but rather by using the ratio of the polar to equatorial circumference of the apparent horizon as discussed in detail in [35]. Cross-checking the two measures (i.e. apparent-horizon distortion and dynamicalhorizon formalism) in the cases where both are possible shows that they are equally reliable (see also the extended discussion in [35]). Overall, the data available suggests the existence of a local maximum of a for q ∼ 0.9, but more data is clearly necessary to confirm this.
Interestingly, when inspecting carefully the apparent horizon in the low-q model M3.5q0.75 it is possible to appreciate that its appearance precedes the time when the two stellar cores merge and is in contrast with what happens with models with high q. By comparison, we believe the same happens also for the binary M3.4q0.70, although in this case we were not able to detect an apparent horizon. Of course these considerations have little physical importance as the interior of the apparent horizon is causally disconnected with what is astrophysically observable; nevertheless this result provides another interesting example of the rich phenomenology relative to the appearance and dynamics of trapped surfaces (see, for instance, the discussion in section 4 of [58]). Table 2. Columns 2 − 5 report the properties of the final BH, i.e. mass, angular momentum, spin parameter, and kick velocity, while columns 6 − 7 report the measured torus masses Mtor and those inferred from relation (6), Mtor. Also shown in columns 8 and 9, respectively, are the numerical error ǫtor ≡ |1 − (Mtor) HR /(Mtor) MR | as computed by comparing different resolutions (medium, i.e. ∆ 6 = 0.19, and high, i.e. ∆ 6 = 0.15) for each model and the relative error ǫ fit ≡ | Mtor − Mtor|/Mtor of the phenomenological expression for the mass of the torus with respect to the numerical data. Clearly, the binaries with high mass ratio are not well described by relation (6) even though their numerical error is not very large.
Model Finally, also reported in table 2 is the recoil velocity imparted to the BH at the end of the inspiral and computed using the gravitational-wave emission as discussed in [59,32] for binary BHs. We recall, in fact, that together with energy and angular momentum, gravitational radiation also carries away linear momentum. If the binary system has a degree of asymmetry (either in the mass or in the spin) then the trajectories of the two bodies will be (slightly) different (e.g. with the smaller body moving more rapidly and, hence, being more efficient in beaming its emission) and the momentum loss in one direction will not be balanced by an equal loss in the diametrically opposite direction. This effect is well-known in binary BHs, where the recoils from quasi-circular inspirals can be as larger as ∼ 4000 km (see [60] for a recent review), but it has never been reported before for binary NSs. The recoil velocities reported in table 2 are clearly much smaller than those measured for binary BHs. However, they could still yield astrophysically interesting results being comparable or larger than the escape velocity from the core of a globular cluster that is v esc ∼ 50 km [61]. Furthermore, and possibly surprisingly, the values reported here for irrotational binaries which have very little initial spin, are not much smaller than those computed for non-spinning binary BHs (see, e.g. , [62] for a recent update) and have a local maximum for q ∼ 0.9.

General Dynamics
In figures 3 and 4 we show color-coded contours of the rest-mass density for models M3.6q1.00 (left panels) and M3.4q0.70 (right panels), either in the (x, y) plane (upper rows) and in the (x, z) plane (lower rows). The snapshots in figure 3, in particular, correspond to the time t ∼ 10 ms when the systems enter the regime of quasi-stationary accretion (QSA, see below for definition), shortly after the formation of the BH, while those in figure 4 refer to the final time of the evolution, t ∼ 21 ms. These figures allow for a closer view of the morphological features of the disks, in particular, their spatial dimensions and thickness, and are a natural continuation of the dynamics already shown in figures 1 and 2, although they use a different colormap that has been tuned to yield a better contrast in the density profiles.
The large morphological differences between these two extreme models are clearly visible in these figures. The equal-mass model produces a highly symmetric, geometrically thin disk, similar to the ones already observed for other equal-mass initial data in [3]. The unequal-mass model, on the other hand, at the time of the onset of the regime of QSA is characterized by the presence of a large spiral arm, which has not yet been accumulated onto the central disk surrounding the formed BH. The asymmetry in the distribution of matter at this stage is also apparent from the color map of the rest-mass density. Only at the end of the simulation the disk of the unequal-mass binary acquires a more axisymmetric shape. The diameters of the disks and their heights perpendicular to the horizontal plane differ in a significant way between the two models. More specifically, at the end of the evolution and  using the ρ = 10 10 g/cm 3 isodensity contour as the reference value below which material is not considered part of the disk, our simulations yield disk diameters of ∼ 50 km for model M3.6q1.00 and ∼ 150 km for model M3.4q0. 70. The corresponding vertical scale is ∼ 5 km and ∼ 35 km, respectively. ‡ Taking into account all the models of our sample, we find that both scales increase as the mass ratio decreases. Even more worth noticing is the fact that, in the cases considered, while the tori differ in size by about a factor ∼ 3, they differ by a factor ∼ 200 in mass, while having comparable mean rest-mass densities (see further discussion in sections 4.2 and 4.8).

Rest-mass Evolution
In order to establish how the asymmetry in the mass of the two NSs in the binary leads to tori with different masses we show in figure 5 the evolution of the total rest mass, defined as normalized to its initial value and for the different models. In this equation W ≡ αu t is the Lorentz factor, α being the lapse function, and γ is the determinant of the spatial metric. All the curves in figure 5 have been shifted in time to coincide at t coll , which represents the (collapse) time at which a rapid decrease of the total rest-mass takes place following the formation of a BH. Note that, in practice, the collapse time is different for all models, ranging from around 6 ms for model M3.4q0.70 to around 11 ms for model M3.4q0.80. ‡ Of course it should be noted that the spatial dimensions reported here depend on the cut-off chosen for the restmass density. Using cut-offs smaller than ρ = 10 10 g/cm 3 would lead to considerably larger estimates for the sizes of the tori.

Figure 5.
Evolution of the total rest masses Mtot normalized to their initial values for all the models considered. The order of magnitude of the mass fraction in the accretion torus can be read off the logarithmic mass scale on the vertical axis. The curves referring to different models have been shifted in time to coincide at t coll , which represents the time when the very rapid decrease of the total rest mass takes place. Note that this time is not physically relevant (the apparent horizon is usually found earlier than t coll ) and simply corresponds to when the very large amount of rest-mass accumulated in a very few cells is numerically dissipated. Figure 5 shows that all models conserve the baryonic mass almost perfectly (i.e. with losses of 10 −6 ) up until the formation of the apparent horizon, after which most of the rest mass disappears in the singularity. One obvious result which can be deduced from figure 5 is that the mass of the resulting accretion disk is larger for smaller values of q. However, this trend is not entirely monotone in the figure as it is also influenced by the initial total baryonic mass of the binary. The particular values of the tori masses computed for all models are reported in table 2. While the equal-mass model produces a disk of barely 10 −3 M ⊙ , models M3.4q0.80 and M3.4q0.70 produce significantly more massive tori with masses of about 0.2 M ⊙ . A more detailed discussion of the mass in the tori will be made in section 4.8.
As the apparent horizon is formed, a substantial part of the rest-mass is still outside it, although it will accrete rapidly onto the BH. This makes the mere definition of what is the torus and its mass rather arbitrary and we decide therefore to define the torus mass M tor as the total rest mass outside the apparent horizon when the disk enters a regime of quasi-steady accretion (QSA), a regime which is found in all models investigated. More specifically, we compute the accretion rate aṡ and define the onset of the QSA as the point in time when the conditionṀ tot /M tot < 10 −6 (Gc −3 M ⊙ ) −2 is satisfied for the first time. In other words, we define the onset of the QSA as the time when the accretion has stabilized and the matter is moving on essentially circular orbits. This definition is again somewhat arbitrary, but has the advantage of allowing for a systematic comparison of the differences in the properties of the accretion tori produced by the several models considered in this work. Figure 6 shows the evolution of the total rest mass M tot (upper rows) and the mass accretion rateṀ tot (lower rows) in the regime of QSA for the two extreme models of our sample, M3.6q1.0 (left panel) and M3.4q0.70 (right panel). Also indicated with a vertical dashed line is the onset of the QSA and it should be noted that both M tot andṀ tot differ by almost two orders of magnitude when comparing equal and unequal-mass binaries. An aspect of the evolution of the accretion rate which is quite evident in figure 6 is the sharp difference between equal-and unequal-mass binaries. The equal-mass case, in fact, shows an accretion rate (and indeed the whole evolution of the torus) that is subject to quasi-periodic oscillations as the torus moves in and out at about the radial epicyclic frequency. The mass flux of the unequal-mass model, on the other hand, is rather constant in time and this reflects a very different distribution of angular momentum in the tori. Both of these aspects will be further discussed in the following sections.

Density evolution
Once the BH is formed, an effective gravitational potential well builds up in which the torus undergoes radial oscillations. In the case of an equal-mass binary, the well is essentially axisymmetric and the dynamics of oscillating relativistic tori in equilibrium and in axisymmetry have been analyzed extensively in a series of papers [63,28,30,26] in the test-fluid approximation (where the self-gravity of the disk is neglected), with and without magnetic fields, and for the cases of Schwarzschild and Kerr BHs. These papers have shown that, upon the introduction of perturbations in the tori, a long-term oscillatory behavior is found, lasting for tens of orbital periods. These oscillations correspond to axisymmetric The rapid drops take place well after an apparent horizon has been formed and are caused by the numerical methods which are no longer able to resolve the very large gradients in the very central grid cells. The two insets provide a magnified view of the evolution of the density in the torus and help to contrast the periodic accretion produced in the case of equal-mass binaries and the QSA for the unequal-mass binaries.
p-mode oscillations whose lowest-order eigenfrequencies appear in the harmonic sequence 2:3. This harmonic sequence is present with a variance of ∼ 10% for tori with a constant distribution of specific angular momentum and of ∼ 20% for tori with a powerlaw distribution of specific angular momentum. More recently, those studies have been extended by [64], where systems formed by a BH (in the puncture framework) surrounded by (marginally stable) self-gravitating disks have been evolved in axisymmetry. Even in this case, the ratio of the fundamental oscillatory mode and the first overtone also shows approximately the 2:3 harmonic relation found in earlier works [28,30,26].
The dynamics of the BH-torus system produced by the merger of model M3.6q1.00 is considerably more complicated than that considered in the test-fluid studies, for which initial configurations in stable equilibrium could be found. However, despite the fact that these systems are studied ab-initio as the end-products of highly dynamical events, it is remarkable that so much of the phenomenology studied and reported in [63,28,30,26] continues to apply also here. Unfortunately, although the simulation extends to ∼ 26 ms, the timeseries is much too short to provide a firm evidence of the presence of the 2:3 harmonic relation, although the spectral analysis of the data indicates that excess power is present at such frequencies.
To provide additional evidence that the harmonic behaviour is not just in the accretion rate, figure 7 shows the evolution of the maximum of the rest-mass density ρ max , normalized to the corresponding initial value, for the two extreme models of our sample, M3.6q1.00 and M3.4q0.70. The equal-mass model M3.6q1.00 (represented in the left panel of figure  7) shows the most regular and pronounced oscillatory behavior, as was already evident in the time evolution of the total rest mass and accretion rate in figure 6. The two insets in this figure magnify these features in the QSA regime and, in the case of M3.6q1.00, they highlight the presence of both maxima and minima corresponding to configurations when the torus reaches the point of closest approach to the BH (periapsis or pericenter) and of farthest excursion (apoapsis or apocenter), respectively. A similar trend can be hinted also for the M3.4q0.70 model on the right panel of Fig 7, although the quality of the oscillations is smaller in this case, most likely because in this case the enhanced tidal disruption during the merger phase leads to a more complex dynamics. Interestingly, such oscillations seem to become more regular during the final stages of the evolution, i.e. for t 17 ms, as the torus reaches a more axisymmetric configuration.
A novel technique to analyze the evolution of the tori and to gain some insight on their dynamics is that offered by spacetime diagrams for observers comoving with the black hole. This is shown in figure 8, which reports the evolution of the color-coded rest-mass density embedded in a spacetime diagram with the (x− x AH ) coordinate on the horizontal axis, where x AH is the position of the apparent horizon, and the coordinate time t on the vertical axis. The color code is indicated to the right of each plot and isodensity contours are shown for the values of ρ = 10 10 , 10 11 , 10 12 , 10 13 g/cm 3 . Note that while these values are the same for all the panels, the spatial dimensions vary considerably. For each model, the dotted horizontal line marks the onset of the regime of QSA.
By comparing the spacetime diagrams for all models it is evident that only the equal-mass model M3.6q1.00 shows a global oscillatory movement with respect to the location of the BH horizon. The movement is indeed global as all the isodensity contours plotted oscillate simultaneously and the maximum and minimum radial extensions reached by the disk (as signalled by the location of the 10 10 g/cm 3 contour) are ∼ 25 km and ∼ 15 km, respectively. It is these oscillations that produce the periodic increase in the maximum rest-mass density reported in figure 7 and it is easy to appreciate that in this case the average density in the disk is less than about 10 12 g/cm 3 .
Scrolling through the different panels in figure 8 it is possible to appreciate that the dynamics of the torus is strongly influenced by the mass ratio. More specifically, models M3.4q0.80, M3.5q0.75, M3.4q0.70, show very rapid expansions corresponding to the ejection of the large spiral waves discussed in the previous sections. As we will comment later on, most of this matter is still bound but it nevertheless reaches distances which are several hundreds of km away from the BH, leading to tori that have spatial dimensions as large as ∼ 80 km. Furthermore, noticeably higher average rest-mass densities are reached in the three low-q models. As a result, while tidal disruption sweeps away a large fraction of the external layers of the less massive star in the binary, the corresponding tori are still able to retain the inner and denser regions; this is particularly the case for models M3.4q0.80 and M3.4q0.70, where the tori reach maximum densities as high as ∼ 10 14 g/cm 3 .
We finally note that, because the matter in the spiral arms is bound, it will eventually fall back onto the tori, where it may lead to enhanced accretion and consequently to a new and delayed gamma-ray emission as the one recently observed in [65]. Determining the dynamics of the fall-back material is therefore of great importance astrophysically and the focus of our attention in future work.

Dynamical Instabilities
As mentioned in the Introduction, current models of GRBs assume that the central engine is a system consisting of a BH and a thick disk, either formed at the late stages of the coalescence of two NSs or after the gravitational collapse of a massive star. The energy supply comes from the energy released by the accretion of disk material on to the BH and from the rotational energy of the BH itself, which can be extracted, for instance, via the Blandford-Znajek mechanism [66]. This vast amount of energy (of the order of 10 53 -10 54 erg, depending on the mass of the disk and on the BH rotation and mass) is sufficient to power a GRB if the energy released can be converted into γ-rays with an efficiency of about a few percent. This scenario requires a stable enough system to survive for a few seconds. In particular, the internal-shock model [67] implies that the duration of the energy release by the source has a duration comparable with the observed duration of the GRB. Any instability which might disrupt the system on shorter timescales, such as the so-called runaway instability [68], could pose a severe problem for the accepted GRB models. The runaway instability was first pointed out in [20] and operates as follows: If the torus is initially filling its Roche lobe, transfer of mass onto the BH is possible through the cusp located at the L 1 Lagrange point. As a result of accretion, the mass of the BH increases, thus leading to a change in the gravitational field of the system and ultimately to a change in the position of the cusp. This can move either inwards (towards the BH) or outwards (away from the BH), and when the latter happens it leads to an increase in the mass transfer and hence to the runaway accretion of the torus on a timescale of a few milliseconds. The runaway instability has been investigated under different assumptions and approximations (see [22,24] and references therein). Early simplified studies based on stationary models showed that, on one hand, the self-gravity of the disk favours the instability, and, on the other hand, there are also parameters which may help to stabilize the disk, such as the rotation of the BH and the radial distribution of specific angular momentum. The first time-dependent, general relativistic hydrodynamical axisymmetric simulations of the runaway instability of tori around BHs were performed by [22,69,25,23], who treated the dynamics of the gravitational field in an approximate way and neglected the self-gravity of the torus. Overall, [22,25,23] found that tori with constant distribution of specific angular momentum were unstable while non-constant (power-law) angular momentum disks were stable. More recently, in [24] the first simulations in full general relativity of marginallystable self-gravitating tori in axisymmetry were performed with the purpose of evaluating the influence of the torus self-gravity on the runaway instability. The results of [24] indicate that the tori are indeed stable irrespective of the angular momentum distribution. It is therefore interesting that the results presented in figure 8, which are not restricted to axisymmetry but are however constrained to much shorter timescales, reach the same conclusion: Selfgravitating tori around BHs, as those produced by the merger of binary NSs, are stable at least on the dynamical timescales investigated here. Additional considerations on the stability of the tori are presented in the following section.

Specific Angular-Momentum Evolution
Besides the rest-mass density, another quantity whose evolution is useful to understand the dynamics of the tori is the the specific angular momentum. This quantity plays an important role in defining the dynamics of point particles around black holes and in defining the equilibrium of non-self gravitating tori around black holes [70]. As mentioned above, we define the specific angular momentum as ℓ ≡ −u φ /u t . We also note that a similar but distinct definition of the specific angular momentum was used in [1], namely j = hu φ . The two definitions have the same Newtonian limit of j Newt = ℓ Newt = Ωr 2 , Ω being the angular velocity. However, it is important to stress that only the definition used here yields the correct zero radial epicyclic frequency for tori with constant specific angular momentum [see eqs. (43) and (45) of [28]]. Figure 9 shows the evolution of the specific angular momentum, where the different panels show the color-coded specific angular momentum for an observer comoving with the BH. The color code is indicated to the right of each plot and in addition the same isodensity contours reported in figure 8 are shown here to aid to follow the dynamics of the matter. The most striking feature to note when scrolling through the different panels in  Figure 9. The same spacetime diagrams as in figure 8 but for the evolution of the specific angular momentum ℓ = −u φ /ut. Note that the isocontours in this case refer to the rest-mass density and are the same as in figure 8. figure 9 is that the radial distribution changes radically but systematically when going from the equal-mass binary M3.6q1.00 over to the most extreme unequal-mass binary considered M3.4q0.70. In particular, while the specific angular momentum is decreasing outwards in models M3.6q1.00, M3.7q0.94, and M3.4q0.91, it is Keplerian and increasing outwards as ∼ x 1/2 for the remaining models (see also the discussion in the following section). Furthermore, the spacetime plots show that the matter located in the outer regions of the disks acquires the largest values of the specific angular momentum. This is particularly visible in the early evolution of model M3.4q0.80, in which a large spiral arm develops extending beyond the computational boundary, and also in the late evolution of model M3.4q0.70, when the corresponding disk reaches the largest radial extension (cf. right panel of figures 3 and 4). Broadly speaking, our simulations show that, in agreement with the results of [1], the smaller the value of q, the more the angular momentum is transported outwards by a torque from the non-axisymmetric object that forms after the merger.
To better highlight the different behaviour of ℓ for different mass ratios we show in figure 10 the profiles along the x-axis for the tori produced by the binaries M3.6q1.00 (blue lines extending to 20 km) and M3.4q0.70 (red lines extending up to 70 km). The profiles are computed in a frame comoving with the BH and for densities ρ > 10 10 g/cm 3 . Different line types refer either to the onset of the QSA (i.e. t ∼ 10 ms, thin solid lines) or to the end of the simulation (i.e. t ∼ 22 ms, thick dashed lines). Quite clearly, the specific angular momentum decreases outward at all times for the equal-mass binary, while it increases outward for the unequal-mass one (although it was initially decreasing at the innermost parts). At this point it is worth remarking that Rayleigh's criterion against axisymmetric perturbations of rotating fluids requires that dℓ/dx ≥ 0 for a dynamical stability [71]. While this criterion is clearly satisfied by model M3.4q0.70, it is equally-clearly violated by M3.6q1.00, which is nevertheless stable. We believe this difference is due to the fact that Rayleigh's criterion assumes that the motion is stationary and purely azimuthal. While this is essentially the case for the unequal-mass binary, which does not show clear evidence of epicyclic oscillations, it does not hold true for the unequal-mass binary, which shows instead large radial epicyclic oscillations. The nonlinear stability of M3.6q1.00, but also of M3.7q0.94 and M3.4q0.91, seems therefore to indicate that Rayleigh's criterion can and should be extended to account for fluids which are subject to large radial excursions.  Figure 11. The same spacetime diagrams as in figure 9 but for the evolution of the angular velocity Ω. Note that the isocontours in this case refer to the rest-mass density and are the same as in figure 8.

Angular-velocity Evolution
In analogy with figures 8 and 9, figure 11 shows the spacetime diagram for the evolution of the angular velocity Ω ≡ u φ /u t for all models of our sample. It is straightforward to notice that for all models the angular velocity decreases with the radial distance from the apparent horizon. While this is qualitatively in agreement with the results of [1], it is worth noting that the radial fall-off is very different as the mass ratio is varied among the different binaries. This is shown in figure 12, which reports the profiles of Ω along the x-axis for the tori produced by the binaries M3.6q1.00 (blue lines extending to 20 km) and M3.4q0.70 (red lines extending up to 70 km). As before, the profiles are computed in a frame comoving with the BH and for densities ρ > 10 10 g/cm 3 and different line types refer either to the onset of the QSA (i.e. t ∼ 10 ms, thin solid lines) or to the end of the simulation (i.e. t ∼ 22 ms, thick dashed lines). It is then clear that while the equal-mass binary has an exponentially decaying profile (cf. long-dashed line), i.e.
which does not change significantly with time, the unequal-mass binary reaches at the end of the simulation a profile which is, especially in the outer parts, essentially Keplerian, i.e. with Ω Kep ∼ x −3/2 (cf. dotted line). This feature, which is also shared by the other low-q binaries, explains the scaling of the specific angular momentum as ℓ ∼ x 1/2 and provides firm evidence that the tori produced in this case will be dynamically stable.

Matter Ejection
As a final but nevertheless important aspect of the formation and evolution of the tori, we consider whether or not a part of the rest-mass of the system is ejected during the merger t (ms) x−x AH (km)  Figure 13. The same spacetime diagrams as in figure 9 but for the evolution of local fluid energy ut. Note that the isocontours in this case refer to the rest-mass density and are the same as in figure 8. and the subsequent evolution. To determine whether a fluid particle is bound or unbound we use the covariant time component of the 4-velocity u t and recall that, in an axisymmetric and stationary spacetime, the value of u t for a particle moving along a geodesic is conserved. If the particle is unbound, it moves outwards and −u t = W > 1 at infinity, where α ≡ 1, β i ≡ 0. The local condition u t > −1 thus provides a necessary although not sufficient, condition for a fluid element to be bound; stated differently, if a particle reaches infinity it is because it has u t < −1. Furthermore, this condition is exact only in an axisymmetric and stationary spacetime, and our spacetimes attain these properties only in the final stages of the evolution. Nevertheless, this is a useful condition for a first estimate of the amount of matter ejected and a in-depth discussion on the assumptions implicit in this criterion and on how it applies if one accounts for external forces are presented in [72]. (Note that the alternative criterion for bound flows, namely hu t > −1, would yield similar results since in the relevant regions h ∼ 1.) Figure 13 shows the evolution of u t embedded in a spacetime diagram much like the ones presented before for the rest-mass density, the specific angular momentum and the angular velocity. For all models under consideration, the criterion u t > −1 is well fulfilled, namely all the matter in the tori is bounded, except for model M3.4q0.80 which clearly shows in the early stages of its evolution, that a certain amount of unbound matter is ejected before reaching the regime of QSA. Only for the outermost, very low-density regions of the tori (which are not shown in the spacetime diagrams of figure 13) values of u t ≤ −1 are encountered in the other models and are probably the manifestation of an outflowing wind caused by the very large temperatures of those regions. As a final remark we note that although the total amount of matter ejected in this way is rather small and only of the order of ∼ 10 −4 M ⊙ , it can nevertheless act as the site for the production of the neutron-rich heavy elements that are formed by rapid neutron capture (i.e. the r-process) (see [73] and references therein). Performing such calculations and thus determining to what extent binary NS mergers contribute to the whole observed r-process material in the Galaxy requires a fully developed reaction network and is outside the scope of this study, but will be the focus of our future research.

A phenomenological expression for the mass in the torus
As mentioned above, determining the amount of rest-mass in the torus may be one of the most important aspects of this research for the impact it has on the modelling of the emission in SGRBs. Table 2 also reports the mass of the torus and since the latter slightly decreases in time, we have arbitrarily chosen the time of t ∼ 17 ms as a reference (it is about the latest time for which we have data from all the simulations). Because of the importance of the information and because of the scarcity of the numerical data available, it would be valuable to derive a phenomenological expression for the mass in the torus which can be constructed on basic expectations and that can be constrained by using the numerical data.
Following this spirit, we first search for a phenomenological expression for the torus mass which will depend only on the mass ratio and on the total mass of the binary, i.e. M tor = M tor (q, M tot ). Next, we exclude the trivial case in which the total mass is larger than the maximum mass of the binary system M max (based on the maximum allowed mass for isolated stars with the given EOS); in practice we impose that M tor (q, M tot ≥ M max ) = 0 for any value of q. Finally we impose the expectation that the mass of the torus should depend, at least at lowest order, on the mass ratio (this was already noted by [1]) and yield the torus with the smallest possible mass for an equal-mass binary. Collecting all these constraints, our ansatz is where in the second expression we have written the maximum mass of the binary in terms of the maximum mass M * of an isolated nonrotating star, i.e. M max = c 3 (1 + q)M * .
Note that as introduced in expression (6), the coefficients c 2 and c 3 have a direct physical interpretation: c 2 is proportional to the mass of the torus for equal-mass binaries, while c 3 parametrises the excess of maximum mass that can be supported in the binary because of the stabilizing effect produced by the nonzero spin of the stars and of the tidal potential (i.e. c 3 is expected to be slightly larger than 1). The three coefficients, c 1 , c 2 , and c 3 can then be computed by comparing expression (6) with the numerical data reported in table 2 as well as the one computed in [3] for equal-mass binaries. The fitting procedure then yields c 1 = 1.115 ± 1.090, c 2 = 0.039 ± 0.023, c 3 = 1.139 ± 0.149, with a reduced χ 2 ≃ 2 × 10 −3 .  [3] (green squares). Also shown in the parameter space (q, Mtot) considered here is the phenomenological modelling Mtor suggested by expression (6). Note that to highlight the functional behaviour of the phenomenological expression, the x− and y−axes are shown as decreasing when moving to the right and to the left, respectively. Figure 14 shows the torus mass M tor either as measured in the simulations reported here (red crosses) or in those presented in [3] (green squares) and against the phenomenological modelling M tor suggested by expression (6) in the region where M tot ≤ M max . Note that to highlight the functional behaviour of the phenomenological expression, the x− and y−axes are shown as decreasing when moving to the right and to the left, respectively. Overall, the figure shows rather generically that: 1) The mass of the torus increases with the asymmetry in the mass ratio; 2) Such an increase is not monotonic and for sufficiently small mass ratios the tidal disruption leads to tori that have a smaller mass for binaries with the same total mass; 3) Tori with masses 0.21 M ⊙ have been measured and even more massive ones, i.e. with masses up to ∼ 0.35 M ⊙ , are possible for mass ratios q ∼ 0.75 − 0.85.
We note that somewhat similar considerations about the mass of the torus were made also in [1], where a different phenomenological expression for the mass of the torus was proposed. When applied to the data computed here, the expression suggested in [1] does not reproduce well the data and yields rather large errors. There are a number of reasons that could justify these differences and that are related to the different initial data chosen (ref. [1] has only two initial total masses which are smaller than those considered here), to the different EOSs employed (ref. [1] uses cold but realistic EOSs in contrast to the ideal-fluid chosen here), and to the different numerical techniques adopted (ref. [1] uses a uniform grid with rather coarse resolution in place of the mesh-refined grid employed here). All these differences make the comparison between the two calculations rather difficult, although they also motivate a closer comparison using at least the same initial data and the same EOSs, and which will be the subject of our future work. However, common conclusions of both calculations are that: The mass of the torus can be as large as ∼ 0.1 M ⊙ and larger; It increases with the mass asymmetry in the binary; It is larger for systems with smaller total mass. We believe these features are robust and will be also present when different initial data and EOSs are considered.
A final note of caution must be mentioned: Although figure 14 indicates a very good match between the data and expression (6), it also shows that the latter is inaccurate for q ≃ 1, where the tori masses are much smaller and the prediction leads to small but negative values; luckily, the regime where (6) is less accurate is also the least interesting one from an astrophysical point of view. Most importantly, however, it is clear that the attempt to produce a phenomenological description for the mass of the torus after having investigated only a small portion of the space of the parameters (especially with respect to the total mass of the binary) and after using as support only 8 simulations is a very demanding task and potentially a flawed one. However, because we believe that expression (6) is a reasonable description of the expected results, we foresee that it will reveal its robustness as additional simulations are performed and the coefficients will be further improved. This will indeed be the subject of our future work. Figure 15 shows the waveforms in the two polarizations of the gravitational-wave amplitude (h + ) 22 (upper panels) and (h × ) 22 (lower panels) for all the models considered and as computed from the gauge-invariant perturbations of a Schwarzschild spacetime. As predicted by the post-Newtonian approximation [74], the inspiral phase is characterized by harmonic oscillations at roughly twice the orbital frequency but that show an increase both in amplitude and frequency as the merger approaches. We note that the initial part of the inspiral of the binary M3.4q0.80 shows a comparatively larger contamination from the initial spurious burst of radiation. This is simply due to the fact that such a binary has been constructed with a comparatively larger initial violation of the constraints (i.e. the violation of the L 2 norm of the Hamiltonian constraint is ∼ 3 × 10 −6 and about 50% larger than the violation measured in other binaries). We believe that this larger initial error is also the one responsible for a longer time spent by this binary before the merger.

Gravitational-Wave Emission
As already discussed before, because of the very high total mass of the systems no transient HMNS forms, whose dynamics would have been dramatically imprinted on the waveforms (cf. the detailed comparison of the HMNS dynamics for different EOSs presented in [3]). As a result, the post-merger waveform is essentially the one corresponding to the collapse of the HMNS to a BH. Indeed, as noted above, in the low-q cases M3.5q0.75 and M3.4q0.70, for which a common apparent horizon is found almost simultaneously with the merger, the part of the waveform produced by the newly formed BH starts essentially together with the end of the one coming from the inspiral.
The ringdown part of the waveform starts increasingly early for binaries with smaller mass ratios and its signature in the waveform is also less evident. More specifically, while the ringdown of the BH created after the merger can be clearly identified in the waveform of the equal-mass model M3.6q1.00, it becomes much less clear as one scrolls down in the different panels of figure 15 and it seems almost absent in model M3.4q0.70. Indeed Figure 15. Gravitational waveforms in the two polarizations h + (upper panels) and h × (lower panels) as computed from the lowest ℓ = m = 2 multipole for all the binaries considered. For those models where it was found, the vertical dashed lines mark the time of the first detection of the apparent horizon. Note that as the mass ratio q decreases, the ringdown part of the signal starts earlier but it is also less evident because of the increasingly large accretion after the formation of the apparent horizon. Finally, shown as an aid to comparison, the panel of the binary M3.4q070 also reports with dotted lines the waveforms for the equal-mass binary M3.6q1.00.
it is necessary to examine M3.4q0.70 on a logarithmic scale in order to appreciate the presence of an exponential ringdown. We believe that this behaviour is mostly likely due to the copious mass accretion after the formation of the apparent horizon that becomes increasingly large as the mass ratio decreases. We recall, in fact, that the mass accretion rate following the BH formation is highly sensitive on the mass ratio and inversely proportional to it (see figure 5 where this is very apparent). Under these conditions of very intense mass accretion, the BH is continuously "hit" by generically nonspherical flows of matter which prevent its natural ringdown, essentially "chocking" it. A detailed analysis on the role played by mass accretion on the properties of the ringdown has already been investigated in [75], where however the BH ringdown was always observed because of the intrinsically perturbative nature of the approach. The rather different accretion regime reached in these simulations suggests therefore that the dynamics observed in figure 15 reflects a nonlinear response of the BH that was not accessible in the work of [75]. Additional work is needed to clarify the relation between hypercritical accretion and BH ringdown and will be the subject of future investigations.
A more systematic analysis of the waveforms as a function of the mass ratio is beyond the scope of this paper and will be considered elsewhere using more realistic or parametrised EOSs. Here, however, as an aid to comparison, the panel relative to the binary M3.4q070 in figure 15 also reports with dotted lines the waveforms for the equal-mass binary M3.6q1.00 and highlights that besides the different amplitude evolution, the mass asymmetry also results into a different phase evolution which is likely to provide important information on the EOS.
In addition, we show with black continuous lines in figure 16 the scaled power spectral densities (PSD) of h + , i.e.h + (f )f 1/2 , for all the binaries considered when placed at a distance of 100Mpc (see [33] for a definition ofh + (f )). Shown also here are the noise curves of the Virgo detector (dotted magenta), of the advanced LIGO detector [76] (short-dashed blue), and of the planned Einstein Telescope (ET) [77] (long-dashed red). Since the number of cycles computed is very small, the peak emission is the one corresponding to the last stages of the inspiral, around 0.6 − 0.7 kHz for all models considered. The amplitude, however, depends sensitively on the mass ratio, being maximal for the high-q binaries and above the noise curve for Virgo in these cases. As the mass ratio is decreased, in fact, the peak values of the PSD decrease and the binaries at the distances considered become then undetectable by an interferometer like Virgo (we recall that the binary M3.4q0.80 has an extended inspiral induced by the larger initial violation of the constraints; hence its PSD amplitude is spuriously increased in figure 16). New-generation detectors such as advanced LIGO will instead be able to reveal the inspiral signal in the frequency interval ∼ 0.3 − 2.0 kHz, while essentially all of the late-inspiral and merger signal would be measured by the Einstein Telescope (see ref. [78] for an introduction to the science reach of third-generation detectors such as ET, and ref [79] for a detailed discussion of the impact that gravitational waves from NSs may have on such detectors). Table 3 summarizes this by reporting the signal-to-noise ratios (SNRs) for the different detectors and clearly highlights that present detectors are unlikely to detect any of the binaries considered here if at a distance of 100 Mpc and observed only during the final part of the inspiral. On the other hand, advanced detectors will be able to reveal these sources even at such large distances (the only ones that can provide an interesting event rate) and, in the case of third-generation detectors such as ET, even measure them with significant SNRs.  Table 3. SNR for the different binaries considered as computed when placed at a distance of 100Mpc for a presently operating detector such as Virgo, as well as for detectors of second and third generation, such as advanced LIGO and ET.

Conclusions
Numerical-relativity simulations of non-vacuum spacetimes have now reached a sufficient stability and accuracy to be able to describe in a complete manner all of the stages of the inspiral, merger and post-merger of binary NSs. Determining the properties of the black-holetorus system produced by the merger represents a key aspect in the modelling of the central engine of SGRBs. Of the many different properties characterizing the torus, the total rest-mass clearly represents the most important one, since it is the torus' binding energy which can be tapped to extract the large energies necessary to power the SGRB emission. However, the rest-mass density and angular momentum distributions in the torus also represent important elements which determine its secular evolution and need to be computed equally accurately for any satisfactory modelling of SGRB engine. As a first step towards modelling ab-initio the central engine of SGRBs, we have here presented new results from accurate and fully general-relativistic simulations of the coalescence of unmagnetized binary NSs with unequal masses as these are the ones expected to yield the largest tori. The evolution of the stars has been followed through the inspiral phase, the merger and prompt collapse to a BH, up until the appearance of a thick accretion disk, which was studied as it enters and remains in a regime of quasi-steady accretion. Although we have employed a simple ideal-fluid equation of state, we have performed a systematic study of the properties of the black-hole-torus obtaining a number of results that can be summarized as follows: • The mass of the torus increases considerably with the mass asymmetry and equal-mass binaries do not produce significant tori if they have a total baryonic mass M tot 3.7 M ⊙ . Those produced have masses M tor ∼ 10 −3 M ⊙ and a radial extension of ∼ 30 km.
• Tori with masses as large as ∼ 0.2 M ⊙ have been measured with binaries having M tot ∼ 3.4 M ⊙ and mass ratios q ∼ 0.75 − 0.85. The tori in these cases are much more extended with typical sizes 120 km. • The mass of the torus can be described accurately by the simple expression M tor (q, M tot ) = [c 3 (1 + q)M * − M tot ] [c 1 (1 − q) + c 2 ], involving the maximum mass for the binaries and some coefficients, both of which can be constrained from the simulations. • Using the phenomenological expression we conclude that tori with masses as large as M tor ∼ 0.35 M ⊙ can be produced for binaries with total masses M tot ∼ 2.8 M ⊙ and q ∼ 0.75 − 0.85. • Tori from equal-mass binaries exhibit a quasi-periodic form of accretion associated with the radial epicyclic oscillations of the tori, while those from equal-mass binaries exhibit a quasi-steady form of accretion. • When analyzing the evolution of the angular-momentum distribution in the tori, we find no evidence for the onset of non-axisymmetric instabilities, that angular momentum is transported outwards more efficiently for smaller values of q thus yielding Keplerian angular-velocity distributions, and that very little of the mass of the tori is unbound. • Present gravitational-wave detectors are unlikely to detect any of the binaries considered here if at a distance of 100 Mpc and observed only during the final part of the inspiral. • Advanced detectors will be able to reveal these sources even at large distances and measure them with significant SNRs in the case of third-generation detectors such as ET.
Overall, these results indicate that large-scale tori with large masses and quasi-stationary evolutions can be produced as the result of the inspiral and merger of binary NSs with unequal masses. Hence, they may provide the energy reservoir needed to power short GRBs. Although complete and accurate, our results are also far from being realistic. Much remains to be done to improve them either by considering physically-motivated EOSs, or by including the effect of magnetic fields, or by taking into account the modifications introduced by a self-consistent treatment of the radiation transfer. All and each of these improvements will be the subject of our future research. , and by the Spanish Ministerio de Educación y Ciencia (AYA 2007-67626-C03-01).

Appendix A. On our accuracy: conservation of mass and angular momentum
In a recent work [5] we have discussed in detail the convergence properties of our numerical simulations and, in particular, the deterioration of the convergence rate at the merger and during the survival of the merged object, when strong shocks are formed and turbulence develops. In particular, in figure 3 of that work we have shown a very stringent measure of the overall conservation properties of our simulations by reporting the time evolution of the energy and angular momentum which are partially radiated during the simulation. In order to reduce the computational costs associated with the measurements made in [5], we had limited ourselves to a single configuration and in particular one that, because of the rather high mass, formed a BH soon after the merger. In particular, we had considered an equal-mass binary with a total baryonic mass of M b = 3.56 M ⊙ and a total ADM mass of M ADM = 3.23 M ⊙ as evolved from an initial (coordinate) separation of ∼ 45 km. As a result, we were able to show an overall conservation of both mass and angular momentum to a precision of ∼ 1% over a timescale of ∼ 10ms. Figure A1. Left panel: Evolution of the maximum rest-mass density normalized to its initial value. Shown with different colours are the different parts of the evolution which are then also magnified in the three insets (cf. the timescale to associate the insets to the different parts of the evolution). Right panel: The same as in the left panel but in terms of the ℓ = m = 2 mode of the h + polarization amplitude.
We here reconsider the same assessment of the conservation properties but for a far more challenging case of a binary with a small total mass and for which the HMNS survives a considerably larger time before collapsing to a BH. In particular, we examine the evolution of an equal-mass binary with a total baryonic mass of M b = 2.912 M ⊙ and a total ADM mass of M ADM = 2.694 M ⊙ as evolved from an initial (coordinate) separation of ∼ 45 km; this same binary was indicated as 1.46-45-IF in [3] and evolved there only up to 25 ms. As representative information about the binary inspiral and merger we report in the left panel of figure A1 the evolution of the maximum rest-mass density normalized to its initial value (cf. figure 15 in [3]). Shown with different colours are the different parts of the evolution and which are also magnified in the different insets (cf. the timescale to associate the insets to the three parts of the evolution). They refer to the immediate formation of the HMNS (red line), to the secular evolution of the HMNS as a contracting bar-deformed object (green line) and to the exponential growth when the threshold to black-hole formation has been crossed (blue line). The right panel of the same figure shows instead the same stages of the evolution but in terms of the ℓ = m = 2 amplitude of the + polarization. Note that the timescale over which the evolution is reported is ∼ 140 ms and thus a factor ∼ 5 larger than the one discussed in [3]. (As for the evolutions in [5], here too we have used a rotational symmetry around the z-axis to reduce computational costs.) Figure A2. Left panel: conservation of energy. The continuous black line is the ADM mass computed as an integral over the whole grid, while the long-dashed blue line is the energy carried from gravitational waves outside the grid and magnified by a factor of 10; the shortdashed red line is the sum of the two and it should be conserved. The numerical violation is at most 1.5% (cf. dot-dashed line). Right panel: the same as in the left panel but for the conservation of the angular momentum. Also in this case the violation is at most 1%.
In analogy with figure 3 of [5], we show in the left panel of figure A2 the evolution of the total mass as normalized to the initial value (cf. left panel of figure 3 of [5]). Indicated with different lines are the volume-integrated values of the ADM mass (solid black line), of the energy lost to gravitational waves magnified of a factor 10 (long-dashed blue line), and of their sum (short-dashed red line). The last quantity should be strictly constant and this is the case to a precision of ∼ 0.5% during the inspiral, but with a secular decrease that brings the total error to be ∼ 1.5% at the end of the simulation (as an aid to comparison the value at 1.015 is shown with a dot-dashed line). Similar considerations apply also to the conservation of the angular momentum as shown in the right panel of figure A2 (cf. right panel of figure 3 of [5]) which uses the same conventions as the left panel (here J vol is computed with the integral (15) in [3]). In this case the radiative losses are much larger (almost 15% of the available angular momentum is lost to gravitational waves) but the overall conservation is still accurate to ∼ 1%. Once again it is worth nothing that the timescale over which we can show accurate conservation of mass and angular momentum is a factor ∼ 14 larger than the one discussed in [5] and provides us with great confidence over the numerical accuracy of our results. Of course this does not provide us with any measure of whether such results are indeed realistic.