Giovanni Corvino

Luciano Rezzolla

Sebastiano Bernuzzi

Roberto De Pietri

Bruno Giacomazzo

On the Shear Instability in Relativistic Neutron Stars

Keywords: numbers: 04.40.Dg, 95.30.Lz, 95.30.Sf 97.60.Jd

. We provide evidence that rapidly and differentially rotating stars that are below the expected threshold for the dynamical bar-mode instability, βc ≡ T /|W | ≃ 0.25, do nevertheless develop a shear instability on a dynamical timescale and for a wide range of values of β. This class of instability, which has so far been found only for small values of β and with very small growth rates, is therefore more generic than previously found and potentially more effective in producing strong sources of gravitational waves. Overall, our findings support the phenomenological predictions made by Watts, Andersson and Jones [2] on the nature of the low-T /|W | instability as the manifestation of a shear instability in a region where the latter is possible only for small values of the β. Furthermore, our results provide additional insight on shear instabilities and on the necessary conditions for their development.

Introduction

Non-axisymmetric deformations of rapidly rotating bodies are rather generic phenomena in Nature and can appear in a wide class of systems. Particularly interesting within an astrophysical context are those deformations taking place in fluids that are self-gravitating and the literature on this has a long history dating back to the work of [START_REF] Chandrasekhar | Ellipsoidal Figures of Equilibrium[END_REF] on incompressible Newtonian uniformly rotating bodies. Since then, the study of these instabilities has continued over the years both in Newtonian gravity and in full general relativity.

Special attention has traditionally been paid to the study of m = 2 instabilities, which are characterized by the exponential growth of m = 2 deformations, where m parametrizes the azimuthal dependence e imφ in a standard mode decomposition in spherical harmonics. Most of the interest in this type of deformation in compact stars stems from the fact that it has the shortest growth time and leads to the emission of a strong gravitational-wave signal.

The development of non-axisymmetric instabilities is commonly analyzed in terms of the quantity β ≡ T /|W | (i.e. the ratio between the kinetic rotational energy T and the gravitational potential energy W ), that provides a dimensionless measure of the amount of angular momentum that can be tapped to feed the development of the instabilities. This parameter plays an important role in what is possibly the most celebrated of the nonaxisymmetric instabilities: the so-called dynamical bar-mode instability. This is an m = 2 instability which takes place when the parameter β is larger than a critical one, β c . In the case of a Newtonian incompressible self-gravitating polytrope, for instance, the dynamical bar-mode instability develops β ≥ β d = 0.2738 [START_REF] Chandrasekhar | Ellipsoidal Figures of Equilibrium[END_REF] and is only weakly dependent on the considered polytropic index or whether the fluid is compressible. Post-Newtonian (PN) studies [START_REF] Saijo | [END_REF] or fully general-relativistic ones [5] correct this results only slightly, by reducing the threshold to somewhat lower values of the instability parameter. As an example, for a polytropic relativistic star with polytropic index Γ = 2, the accurate calculations reported in [6] reveal that the critical value is β c ∼ 0.245 and that a simple dependence on the stellar compactness allows one to track this threshold from the Newtonian limit over to the fully relativistic one [7].

The onset and development of the bar-mode instability has been traditionally studied by means of nonlinear 3D simulations of Newtonian stars that are either unmagnetized [8,9,10,11,12] and, more recently, also magnetized [13]. In addition, PN and fully relativistic simulations have been performed and highlighted, for instance, that the persistence of the bar is strongly dependent on the degree of overcriticality and is generically of the order of the dynamical timescale. Furthermore, generic nonlinear mode-coupling effects between the m = 1 and the m = 2 mode appear during the development of the instability and these can severely limit the persistence of the bar deformation and eventually suppress the bar deformation [6]. These results have been recently confirmed by the perturbative calculations in [14].

Besides dynamical instabilities, which are purely hydrodynamical, secular instabilities are also possible in rotating compact stars and these are instead triggered by dissipative processes, such as viscosity or radiation emission. If, in particular, the dissipative mechanism is the emission of gravitational radiation, then the secular instability is also known as Chandrasekhar-Friedman-Schutz or CFS instability [15,16]. Contrary to what their name may suggest, secular instabilities do not necessarily develop on secular timescales (although they normally do) and are characterized by having a much smaller threshold for the instability. Once again, in the case of a Newtonian polytrope, the critical secular instability parameter is as small as β c ∼ 0.14 and thus much more easy to attain in astrophysical circumstances.

Although widely observed in numerical simulations, the physical conditions leading to a dynamical bar-mode instability are difficult to be encountered in standard astrophysical scenarios. Such large values of the instability parameter, in fact, cannot be easily attained in old and cold neutron stars, which have been brought into uniform rotation and thus to rather small values of β. However, more recently these pessimistic prospects have been changed when a new m = 2 instability has been discovered in differentially rotating Newtonian stars [17] for values of β ≈ 0.01, therefore well below the expected values for a dynamical bar-mode instability. The most salient aspect of this new instability is that it appears in stars with a large degree of differential rotation and that it grows on a timescale which is longer but comparable with the dynamical one. This instability has been referred to as the "low-T /|W | instability" and its dependence on the polytropic index and on the degree of differential rotation has been studied in [18]. Since then, the instability has been observed or discussed in a number of related studies [19,20,21,22,[START_REF] Ott | Classical and Quantum[END_REF]24,25], all of which have highlighted the possible occurrence of this type of instability during the collapse of a massive stellar core.

Despite the abundant numerical evidence on the development of this instability, the nature of these low-T /|W | instabilities is still matter of debate and, most importantly, a sufficient criterion for its onset has not been derived yet. This instability has been studied in great detail by Watts and collaborators [26,2], who have made a number of phenomenological predictions either using a toy shell-model first introduced in [27,28], or for a stellar model in Newtonian gravity. Overall, the work of Watts and collaborators (but see also [29]) recognizes the low-T /|W | instabilities as the manifestation of a more generic class of instabilities, the shear instabilities [2] , that is unstable oscillations that do not exist in uniformly rotating systems and are associated to the existence of a corotation band [30,31]. Watts, Andersson and Jones suggest, in particular, that a necessary condition for the development of the instability is "corotation", that is the presence of a point at which the star rotates at the same pattern speed of the unstable mode [2]. An alternative suggestion on the necessary conditions has been made also by Ou and Tohline [20], who instead associate the development of the instability to the presence of a minimum in the vortensity profile of the star. This minimum can then drive unstable not only the corotating m = 2-modes but also the odd modes such as the m = 1 and m = 3-modes [32]. In this interpretation, the growth time of the instability is proportional to distance between the corotation radius, i.e. the radial position at which the unstable mode corotates with the star, and the minimum of the vortensity.

The purpose of this work is to shed some light on the development of shear instabilities and, in particular, to validate one prediction made, although not explicitly, by Watts, Andersson and Jones. More specifically, we show that, for sufficient amounts of differential rotation, shear instabilities develop for any value of the instability parameter β and also below the expected critical value for the dynamical bar-mode instability. We therefore provide evidence that the low-T /|W | instability is not a new instability but rather the manifestation of a shear instability in a region where the latter is possible only for small values of β.

Our analysis proceeds via the simulation in full general relativity of sequences of neutron star models having constant rest-mass and constant degrees of differential rotation, but with different amounts of rotation, i.e. with different values of β. The neutron-star matter is described by a realistic equation of state (EOS) defined by the unified SLy prescription [1] and we study the development of the non-axisymmetric instabilities from their linear growth up to the fully nonlinear development and suppression. Interestingly, we find that depending on the degree of differential rotation, the shear instability leads either to the growth of a single modes (for the low-β models) or to the simultaneous presence of up to three unstable modes (for the high-β models), which produce beatings in the growth of the overall m = 2 deformation. Special attention is also paid to the properties of the unstable modes and to their position within the corotation band or the vortensity profiles. In this way we are able to confirm both the necessary conditions proposed so far for the onset of the instability. In particular, we show that all the unstable modes are within the corotation band of the progenitor axisymmetric model (cf. [2]) and that all of the unstable models have vortensity profiles with a local minimum (cf. [20]).

The structure of the paper is as follows: in section 2 we describe the numerical setting of our simulations, the EOS we used and the initial models we generated. In section 3 we describe the quantities and tools we used to monitor the evolution of the instability. In section 4 we report the results of the simulations and section 5 is dedicated to conclusions and discussion. We use a spacelike signature (-, +, +, +) and a system of units in which c = G = M ⊙ = 1 (or in cgs 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.

Numerical Setup and Initial Models

In what follows we provide a brief overview of the numerical setup used in the simulations, of the realistic EOS adopted and on the procedure followed for the construction of the initial axisymmetric models.

Numerical Setup

We solve numerically the full set of Einstein equations

G µν = 8πT µν , (1)
where G µν and T µν are the Einstein tensor and the stress-energy tensor, respectively. The equations are solved within the "3+1" decomposition of spacetime, in which the 4dimensional metric g µν is decomposed into the spatial metric γ ij , the lapse function α and the shift vector components β i . The field equations, which then also provide an evolution for the extrinsic curvature tensor K ij , are then coupled to those of general relativistic hydrodynamics

∇ µ T µν = 0 ; ∇ µ (ρu µ ) = 0 , (2)
where, in the case of a perfect-fluid, the stress-energy tensor is given by

T µν = ρ 1 + ǫ + p ρ u µ u ν + pg µν . (3)
Above u µ is the fluid 4-velocity, p is the fluid pressure, ǫ the specific internal energy and ρ the rest-mass density, so that e = ρ(1 + ǫ) is the energy density in the rest frame of the fluid. The set of hydrodynamics equation is then closed by a prescription for the properties of the matter in the form of a relation between the pressure and other quantities in the fluid, e.g. p = P (ρ, ǫ), and for which we have chosen a cold and realistic EOS which will be discussed in the following section. The evolution of the spacetime was performed using the CCATIE code, a threedimensional finite-differencing code providing a solution of a conformal traceless formulation of the Einstein equations (see [33] for the explicit expressions of the equations solved in the code and also [34] for a more recent and improved implementation). The relativistic hydrodynamics equations, on the other hand, were solved using the Whisky code, which adopts a flux-conservative formulation of the equations as presented in [35] and highresolution shock-capturing schemes or HRSC (see [START_REF] Baiotti | Computational Astrophysics in Italy: Methods and Tools[END_REF][START_REF] Baiotti | [END_REF]38] for the explicit expressions of the equations solved in the code and also [39] for a more recent extension of the code to MHD). The Whisky code implements several reconstruction methods, such as Total-Variation-Diminishing (TVD) methods, Essentially-Non-Oscillatory (ENO) methods [40] and the Piecewise Parabolic Method (PPM) [41]. Also, a variety of approximate Riemann solvers can be used, starting from the Harten-Lax-van Leer-Einfeldt (HLLE) solver [42], over to the Roe solver [43] and the Marquina flux formula [44] (see [START_REF] Baiotti | Computational Astrophysics in Italy: Methods and Tools[END_REF][START_REF] Baiotti | [END_REF] for a more detailed discussion). All the results reported hereafter have been computed using the Marquina flux formula and a PPM reconstruction.

Both the Einstein and the hydrodynamics equations are solved using the vertex-centered adaptive mesh-refinement (AMR) approach provided by the Carpet driver [45]. Our rather basic form of AMR consists of box-in-box structures centered on the origin of the coordinate system and with the finest grid covering the whole star at all times. The simulations reported here make use of 4 levels of refinement, with the finest having a resolution of 221 m and the coarsest one a resolution of 1.77 km. The outer boundary was set relatively close to the star and at a distance of ≃ 159.5 km, i.e. at about ≃ 10 times the size of the star. A reflection symmetry across the (x, y) (equatorial symmetry) plane was used to reduce the computational costs, but not a rotational one around the z-axis (π-symmetry) as it would have artificially prevented the growth of odd-m modes (see discussion in [6]).

Realistic Equation of State

As mentioned above, the system of hydrodynamics equations needs to be closed by an EOS relating the pressure with the other primitive variables, e.g. the rest-mass density. Previous studies of the bar-mode instability, both in Newtonian gravity and in general relativity, have been focused on the use of ideal fluids and analytic EOSs, either in the form a of a "polytropic" (and isentropic) EOS p = p(ρ), or of an "ideal-fluid" and (non-isentropic) EOS p = p(ρ, ǫ) (cf. discussion in [6,7]). While these two descriptions are expected to provide results that are qualitatively correct, a more accurate modelling of these instabilities in compact stars necessarily requires a more physically-motivated description of the neutron-star matter.

It is in this spirit that we have here considered a realistic EOS, namely the unified SLy EOS [1], which models high-density and cold (i.e. zero temperature) matter via a Skyrme effective potential for the nucleon-nucleon interactions. The SLy EOS, which describes via a single effective Hamiltonian the neutron star's interior, is supplemented with the HP94 EOS [46] to describe the crustal matter and with the BPS EOS [47] for lower density regions. This prescription results in a one-parameter EOS in the form p = p(e(ρ)) = p(ρ), where the SLy EOS is used for ρ > 4.979 × 10 10 , the HP94 EOS is used for 10 8 ρ 4.979 × 10 10 and the BPS EOS for ρ 10 8 (see also figure 1 of [48]). In addition, at even lower densities the EOS becomes temperature dependent (and thus no longer a simple barotropic EOS), but because these these regions are well below the threshold for the artificial atmosphere, we do not consider an additional prescription for ρ 10 8 g cm -3 .

We recall, in fact, that our HRSC methods require the use of a tenuous atmosphere which fills the regions of the computational domain not occupied by the compact star. The threshold value for the rest-mass density of the atmosphere is chosen to be several orders of magnitude smaller than the maximum value and in our simulation a fluid element is considered to be part of the atmosphere if its rest-mass density ρ satisfies ρ/ max (ρ) ≤ ×10 -8 . When this happens the fluid element is treated as a non-dynamical cold fluid described by a polytropic EOS, p(ρ) = Kρ Γ , with Γ = 2 and its velocity is set to zero (see [49] for a more detailed discussion on the use of the atmosphere in the Whisky code).

The practical implementation of the realistic EOS can take place in a number of different ways. The simplest is to use standard interpolation techniques, e.g. based on Lagrangian polynomials, on the values of the published tables. While straightforward, the interpolations in this approach do not guarantee in general that the thermodynamics relations are fulfilled (see [50] for a thermodynamical preserving interpolation). In addition, the derivatives of the fields, e.g. of the pressure to evaluate the sound speed, are typically not available in the tables. Furthermore, the use of high-order interpolation and/or finite differences can lead to undesirable spurious oscillations.

A second approach that removes all of these problems, uses analytic fits that have been proposed for the pressure. As an example, ref. [48] suggested to fit the specific internal energy of the unified SLy EOS table with the expression ‡

ǫ = p 1 ρ p2 + p 3 ρ p4 (1 + p 5 ρ) 2 f 0 {-p 6 (log(ρ) + p 7 )} + ρ 8 × 10 -6 + 2.1ρ 0.585 f 0 {p 6 (log(ρ) + p 7 )} ,( 4
)
where f 0 {x} = 1/(e x + 1), ρ and ǫ are in cgs units, and the coefficients p i are p i = {0.320, 2.17, 0.173, 3.01, 0.540, 0.847, 3.581} (see Table 2 of [48]). Equation ( 4) is obtained from Eq. ( 15) of [48] after using ρ = m B n, where n is the baryon number density and m B = 1.66 × 10 -24 g is the mass of the nucleons. As discussed in [51], it is then possible ‡ Note that a different notation is used in [48] for some primitive variables.

to compute the pressure from the value of ǫ using the first principle of thermodynamics at

T = 0 p = ρ 2 dǫ dρ , (5)
and thus to have an evaluation of the pressure which is thermodynamically consistent. The differences between the fit and the table are typically less then 2%. Unfortunately, although apparently very convenient, the evaluation of the fitting formulas containing several exponential and logarithmic functions, turns out to be computationally rather expensive even if done in a optimized way. As a third approach, which combines the efficiency of a table search with the thermodynamical consistency of an analytic fit, consists of performing a simple linear interpolation among the tabulated values constructed from the analytic fit. Besides being highly efficient, a linear interpolation also eliminates the spurious oscillations that arise, for instance, in the derivative of the pressure if high-order interpolation formulas are used. In this case, the interpolation error can be reduced simply by populating the analytically constructed tables with a large number of entries, e.g. ∼ 600 in place of the ∼ 150 which are typically available in published tables. This third approach is the one actually implemented in Whisky and provides a speed up of about 20% with respect to the evaluation of the pressure via the analytic fits and with comparable accuracy.

Initial Data

The initial data for our simulations are prepared as stationary and axisymmetric equilibrium solutions for rapidly rotating relativistic stars [52]. Adopting spherical quasi-isotropic coordinates, the line element of the corresponding spacetime is ds 2 = -e µ+ν dt 2 + e µ-ν r 2 sin 2 θ(dφ -ωdt) 2 + e 2ξ (dr 2 + r 2 dθ 2 ) , (6) where µ, ν, ω and ξ are functions of r and θ. Moreover we assume the usual relativistic j-constant law of differential rotation and that amounts to assume an angular-velocity distribution of the form

Ω c -Ω = 1 Â2 r 2 e (Ω -ω)r 2 sin 2 θe -2ν 1 -(Ω -ω) 2 r 2 sin 2 θe -2ν , (7)
where r e is the coordinate equatorial stellar radius and the coefficient  provides a measure of the degree of differential rotation. Expression (7) represents the general-relativistic equivalent of the simpler Newtonian j-constant law [27]

Ω c -Ω = Ω c r 2 sin 2 θ ( Â2 r 2 e + r 2 sin 2 θ) . (8)
Clearly,  → ∞ corresponds to a star in uniform rotation, while  → 0 corresponds to a star with increasing degree of differential rotation. As a reference,  = 1 yields a star with an angular-velocity profile which varies of a factor ∼ 3 between the center and the surface of the star (cf. left panel of figure 2).

In practice, we have computed a very large number of initial models using the SLy prescription for the EOS for which we have computed baryonic mass M b , the gravitational mass M , the angular momentum J, the rotational kinetic energy T and the gravitational Figure 1. Left panel: position in the (β, M/Re) plane of the initial models computed with  = 1, with the filled circles representing those we have evolved numerically (A similar behaviour is shown also by the models with  = 2. Indicated with solid thin lines are isocontours of constant baryon mass models while indicated with a thick dashed line is the threshold to the dynamical bar-mode instability as computed for a Γ = 2 polytrope [7]. Note that the threshold for the instability tends to increase for smaller rest-masses. Right panel: the same initial models as in the left panel but shown in (ρc, rp/re) planes with isocontours of constant β (upper part) or constant M/Re (lower part). See Table 1 for a summary of the properties of the initial models.

binding energy W defined as

M b ≡ d 3 x √ γW L ρ , M ≡ d 3 x -2T 0 0 + T µ µ α √ γ , (9)
E int ≡ d 3 x √ γW L ρǫ , J ≡ d 3 x T 0 φ α √ γ , (10)
T ≡ 1 2 d 3 x ΩT 0 φ α √ γ , W ≡ T + E int + M b -M , (11)
where √ γ is the square root of the determinant of three-dimensional metric γ ij and W L = αu 0 is the fluid Lorentz factor. We stress that the definitions ( 9)-( 11) of quantities such as J, T , W and β are meaningful only in the case of stationary axisymmetric configurations and should therefore be treated with care once the rotational symmetry is lost.

Out of this large set, we have then selected for the numerical evolution a number of models so as to build sequences of constant baryonic mass with M b = 2.5 M ⊙ , degree of differential rotation  = 1, 2 and values of the instability parameter β ranging between 0.140 and 0.250. Each model in this sequence is supramassive, namely it has a mass which is larger than the maximum mass allowed for a corresponding nonrotating model, i.e. M max | Ω=0 = 2.05 M ⊙ , M b,max | Ω=0 = 2.43 M ⊙ , although it is not hypermassive, namely it does not have mass which is above the maximum mass for a uniformly rotating model, i.e. M max | Ω=Ωmax = 2.41 M ⊙ , M b,max | Ω=Ωmax = 2.84 M ⊙ . There are two different reasons why such a large-mass model has been chosen. The first one is that we are interested in the development of a shear instability in the metastable star produced by the merger of a binary system of two neutron stars. As shown by a number of authors and most recently in [49], the product of this merger is either a supramassive or a hypermassive neutron star. Hence, our reference model has a mass which is sufficiently large so as to be a reasonable Table 1. Main properties of the simulated stellar models. Starting from the left: the name of the model, the differential rotation parameter Âthe instability parameter β, the compactness M/Re, the central rest-mass density ρc, the maximum of the density ρmax (Note that for models with  = 2, ρmax = ρc.), the ratio between the polar and the equatorial coordinate radii rp/re, the proper equatorial radius Re, the gravitational mass M , the total angular momentum J divided by the square of the gravitational mass. approximation to the product of a binary neutron star merger §. The second reason is that a sufficiently massive model is necessary in order to reach values of the instability parameter which are above the expected threshold for dynamical bar-mode instabilities as computed in [7], i.e. β max ∼ 0.25. Indeed, the maximum possible value for β within the computed sequence is around β max ≃ 0.2533 and thus just above the threshold (cf. long-dashed line in the left panel of figure 1). Note also that the threshold for the instability tends to increase for smaller rest-masses. Overall, when comparing with equilibrium models generated with a polytropic EOS with K = 100 and Γ = 2 (see table 1 of [6]) the realistic EOS models reach higher compactness (models in [6] typically have M/R e ∼ 0.1) but lower values of β (i.e. β max ∼ 0.2533 for the models considered here, while β max ∼ 0.28 for the polytropic models considered in [6]).

Model

The whole space of parameters is shown in figure 1, whose left panel reports the position in the (β, M/R e ) plane of the initial models computed with  = 1, and where the filled circles represent those we have evolved numerically. Indicated with solid thin lines are isocontours of constant baryon mass models while indicated with a thick dashed line is the threshold to the dynamical bar-mode instability as computed for a Γ = 2 polytrope [7]. The right panel reports the same initial models considered in the left one but shown in (ρ c , r p /r e ) planes with isocontours of constant β (upper part) or constant M/R e (lower part). The main properties of the simulated models are also reported in table 1, where we also introduce our naming convention. Any initial model is indicated as M.%.#, with % being replaced by the § We note that the use of a supramassive model has also the inconvenient consequence that it is not possible to construct a corresponding nonrotating model and this prevents us, for instance, from computing the frequency of the fundamental mode and compare it with the results of perturbative analyses. value of the differential-rotation parameter  and # by the instability parameter β. As an example, M.1.200 is the star with  = 1 and β = 0.200.

Finally, shown in figure 2 are the angular-velocity profiles (left panel) and the rest-mass density profiles (right panel) of some representative models, namely M.1.150, M.1.200 and M.1.250. Indicated with different symbols, which match the ones in the left panel of figure 5, are the normalized radial positions of the corotation radii, with the one for model M.1.150 being shown filled to help distinguish it from the others. Although we postpone to sections 4.1 and 4.2 the discussion of the implications of these corotation radii, two aspects of the initial data are worth emphasizing. The first one is that the amount of differential rotation for a given value of  effectively decreases when increasing the instability parameter β (cf. left panel of figure 2), thus resulting in a smaller corotation band for models with large β. The second one is that all the initial models evolved are axisymmetric but have a "toroidal-topology", namely have the maximum density ρ max that is not at the center of the star, and thus ρ c < ρ max . This toroidal deformation increases with the rotation and thus with β.

Figure 2.

Left panel: Initial angular-velocity profiles for three representative models with small, medium and high values of β and  = 1 (A similar behaviour is shown also by the models with  = 2). Indicated with different symbols, which match the ones in the left panel of figure 5 and6, are the normalized radial positions of the corotation radii, with the one for model M.1.150 being shown filled to help distinguish it. Right panel: the same as in the left panel but for the initial rest-mass density.

Methodology of the analysis

A number of different quantities are calculated during the evolution to monitor the dynamics of the instability. Among them is the quadrupole moment of the matter distribution

I jk = d 3 x √ γW L ρ x j x k . ( 12
)
which we compute in terms of the conserved density √ γW L ρ rather than of the rest-mass density ρ or of the T 00 component of the stress energy momentum tensor. Of course, the use of √ γW L ρ in place of ρ or of T 00 is arbitrary and all the three expressions would have the same Newtonian limit, though with different amplitudes for the gravitational waveforms produced (see for example [53]). However, we here adopt the form (12) because √ γW L ρ is a quantity whose conservation is guaranteed by the form chosen for the hydrodynamics equations.

The quadrupole moment ( 12) can be conveniently used to quantify both the growth time of the m = 2 instability τ 2 and the oscillation frequency once the instability is fully developed σ i 2 . (Hereafter we will indicate respectively with τ i and σ i (n) the growth time and frequencies of the m = n unstable modes and we note that, as will be discussed later on, during the simulation a number of different modes appear, thus justifying the use of the upper index "i"). In practice, we perform a nonlinear least-square fit of the xy component of the computed quadrupole I jk (t) and we generally use as fitting function a sum of N (usually three) exponentially modulated cosines

I jk (t) = N i=1 I jk 0(i) e t/τ (i) cos(2π σ i (n) t + φ (i) ) , (13)
where I jk 0 ≡ I jk (t = 0). Because we commonly have only about 10 cycles in the time interval considered for the fits, extreme care needs to be applied when computing the growth time, especially when the oscillation frequencies and the growth times are close to each other. In these cases, in fact, variations of the initial phase of the modes φ (i) can result in large variation of the growth times. In view of this, we will not report them.

Using three components of the quadrupole moment in the (x, y) plane we can define the distortion parameters η + (t) and η × (t), as well as the axisymmetric mode η 0 (t) as

η + (t) ≡ I xx (t) -I yy (t) I xx (0) + I yy (0) , η × (t) ≡ 2I xy (t) I xx (0) + I yy (0) , η 0 (t) ≡ I xx (t) + I yy (t) I xx (0) + I yy (0) , (14)
so that the modulus η(t) and the instantaneous orientation of the bar are given by

η(t) = η + (t) 2 + η × (t) 2 , φ bar (t) = tan -1 2I xy (t) I xx (t) -I yy (t) . (15)
Finally, as a useful tool to describe the nonlinear properties of the development and saturation of the instability, the rest-mass density is decomposed into its Fourier modes P m (t) as

P m (t) ≡ d 3 x ρ e imφ ,
where φ = tan -1 (x/y).

The phase φ m ≡ arg(P m ) essentially provides the instantaneous orientation of the m-th mode when the corresponding mode has a nonzero power. Note that despite their denomination, the Fourier modes (16) do not represent proper eigenmodes of oscillation of the star. While, in fact, the latter are well defined only within a perturbative regime, the former simply represent a tool to quantify, within the fully nonlinear regime, what are the main components of the rest-mass distribution. As a final comment we note that while all quantities ( 12)-( 16) are expressed in terms of the coordinate time t and are not invariant measurements, the lengthscale of variation of the lapse function at any given time is always larger than twice the stellar radius at that time, ensuring that events on the same timeslice are also close in proper time.

Results

In what follows we first describe the dynamics of the shear instability as deduced from the numerical simulations and then contrast our results with the phenomenological predictions on the necessary conditions for its development.

Dynamics of the instability

As mentioned when discussing the initial data, we have evolved numerically two sequences of constant baryonic mass, with M b = 2.5 M ⊙ . The first sequence has a higher degree of differential rotation (i.e. Â = 1) and instability parameter ranging from β = 0.140 to β = 0.250. The second sequence, on the other hand, has a smaller degree of differential rotation (i.e. Â = 2) and instability parameter between β = 0.125 and β = 0.200. For both sequences the highest value of β considered is also very close to the highest attainable with our initial data code (see figure 1). To keep the computational costs to an affordable level we have evolved of all these models for ∼ 15 ms and although this time window is in general insufficient to capture the suppression of the instability, is adequate to measure the frequencies of the unstable modes and provide a first estimate of the growth times.

After analyzing the results following the method outlined in the previous section and focusing on the properties of the distortion parameters η + , η × and η 0 , we find that all of the models show an m = 2 instability. The maximum distortions obtained during the simulation time are η ∼ 0.1 (i.e. distortions of 10% with respect to the axisymmetric progenitor), with the maximum values being reached for the models of the  = 1 sequence. The models of the  = 2 sequence, in fact, have in general much smaller distortions, with η ∼ 0.01 over the timescale over which they have been evolved. We believe this is simply the consequence of the fact that the models in this sequence have smaller growth rates (see also the discussion in the next section).

To show the general behaviour of the simulation we focus on the description of the result of the simulation for the  = 1 sequence as this is representative also of the one with  = 2. In particular, in the upper part of each panel in figure 3 we show the time evolution of the distortion parameter η + computed with Eq. ( 14) for six representative models: M.1.150, M.1.160, M.1.200, M.1.220, M.1.241 and M.1.250. Similarly, in the lower parts of each panel we show the corresponding evolution of the Fourier-modes computed from Eq. ( 16) for m = 1, 2, 3, 4.

Note that for models with lower values of β the bar-mode deformation is very similar to the one already discussed in [6,7], growing exponentially and with only one unstable mode appearing (cf. first row of figure 3). However, as the rotation is increased, the development of the instability is more complex and at least two unstable modes appear which develop in different parts of the star. These two modes have very similar but distinct frequencies and growth times, leading to a series of beatings in the evolutions of η + , whose irregular evolution makes the calculation of the growth times challenging (cf. second row). As the rotation is increased towards the maximal values of β, up to three distinct modes appear and the evolution of the instability is correspondingly more complex (cf. third row). A similar behaviour is shown also by the evolution of the Fourier modes. The two models with lower β, in fact, show a clear growth of the m = 2 mode and of the m = 1, with the latter becoming first comparable and then larger when the instability is suppressed. In models M.1.250, M.1.241, on the other hand, the m = 2 and m = 3 modes have comparable amplitude for a long period and then the m = 2 becomes the dominant one (this is more clear in model M.1.241). Finally, for larger values of β, the m = 1 mode never attains values comparable with either the m = 2 or the m = 3, which instead control the evolution.

Figure 4 reports the power spectral density (PSD) of η + and η 0 for the same 6 models (solid black lines) and allows to appreciate how the spectrum changes as the instability parameter is increased. More specifically, it is very apparent that at low β the spectra have only one peak, whose maximum is marked by the vertical blue dotted lines and that is present at all the values of beta considered; this is the mode we refer to as σ 2 2 (cf. table 2 and figure 5). As β increases, the amplitude of this peak decreases and it becomes the weakest for very large rotation rates. Starting from model M.1.200 (although a hint of a peak is Table 2. Results of the analysis of quadrupole evolutions. The frequencies σ f and σ i 2 are obtained from the position of the peaks in the PSD while the frequencies σi 2 are obtained from the nonlinear fitting of Eq. (13). Reported in the last two columns are the edges of the corotation band as expressed in terms of angular frequency at the surface Ωe/π ≡ Ω(r = Re)/π and on the axis Ωc/π. Note that depending on the rate of rotation and the degree of differential rotation some frequencies may not be present. present already in M.1.160), a second peak appears at higher frequency and is marked with magenta vertical dot-dashed lines; this is the mode we refer to as σ 3 2 and although it never becomes the largest one, its amplitude increases with β. Finally, a third peak appears at low frequency in the spectra starting from model M.1.220 (even though a hint is present also for model M.1.210) and is marked with black vertical long-dashed lines. This mode, which we refer to as σ 1 2 , becomes the dominant one at high β. Note also that in each panel of figure 4 we report with a dashed blue box the corotation band that will be further discussed in the following section.

Model

σ f σ 1 2 σ 2 2 σ 3 2 σ1 2 σ2 2 σ3 2 Ωe/π Ωc/π (kHz) (kHz) (kHz) (kHz) (kHz) (kHz) (kHz) (kHz) (kHz)
As a final remark we note the peaks in the spectrum of I xx + I yy (red dashed lines in figure 4) are related to axisymmetric m = 0 modes and most likely to f -modes oscillation excited by the development of the instability. Indeed, their frequencies σ f match reasonably well the phenomenological fit for the f -mode frequencies of nonrotating neutron stars with realistic EOSs computed in [54]. Because our sequences do not contain nonrotating configurations (the rest mass is larger than the maximum one), this association is just qualitative and a more detailed investigation of the frequency spectra of the equilibrium models is necessary to confirm this suggestion.

Necessary conditions for the instability

To support the interpretation of these instabilities as shear instabilities we have considered whether the necessary conditions suggested by [2] and by [20] are met. We recall that Watts and collaborators pointed out that an unstable configuration should have the unstable mode with a frequency within the corotation band. In the case of the Newtonian expression (8), the corotation band is simpler to compute and for a mode with frequency σ and azimuthal number m, this is simply given by

Ω c Â2 Â2 + 1 < σ m < Ω c . (17)
We have therefore checked whether any of the unstable m = 2 (bar-modes) with computed frequency σ i (2) has pattern speed velocity, σ i 2 /m = σ i 2 /2, within the corotation band. A careful and rather involved analysis has indeed confirmed the prediction of Watts and collaborators: namely, all the unstable modes are within the corotation band of the progenitor axisymmetric model. These results for the models with  = 1 are listed in table 2, where we report the edges of the general-relativistic corotation band in terms of angular frequency at the surface Ω e /π ≡ Ω(r = R e )/π and on the axis Ω c /π, as well as the frequencies of the unstable m = 2 modes obtained from the position of the peaks in the PSD (i.e. σ i 2 ) or from the nonlinear fitting of Eq. ( 13) (i.e. σi

2 ). Note that depending on the rate of rotation and the degree of differential rotation some frequencies may not be present in the corotation band.

The data in table 2 is also shown in the left panel of figure 5, where we plot the position of the non-axisymmetric frequencies σ i 2 within the corotation band as a function of β; indicated with a filled symbol, to distinguish it from the others, is the σ 1 frequency for model M.1.150 (cf. figure 2 and6). The right panel of the same figure shows the corresponding information and on the same scale but for the models with  = 2 (the inset shows instead a magnified view). Clearly, for both sequences all the unstable modes are within the band and, as the stellar rotation rate increases, more unstable modes appear for the same value of β. 2 for models with  = 1 and shown as a function of β; indicated with thick solid lines are the edges of the corotation band, i.e. the frequency interval between Ωc/π and Ωe/π. Shown with a filled symbol, to help distinguish it from the others, is the σ 1 2 frequency for model M.1.150 (cf. figure 2 and6). Right panel: the same as in the left one and with the same scale but for models with  = 2; the inset shows instead a magnified view. Clearly, in both panels all the unstable modes are within the band and up to three unstable modes appear for the  = 1 sequence with increasing β.

Watts, Andersson and Jones also give a qualitative description, summarized in their figure 2 of [2], of how the unstable and stable models should be distributed in the ( Â, β) plane. The considerations they make are particularly simple: for high degrees of differential rotation, that is at low values of Â, the corotation band is rather wide and there will be both an upper value and a lower value of β between which the shear instability can develop (these critical values of β correspond to the entrance and exit of the unstable mode in the corotation band). This situation corresponds therefore to the one commonly encountered in numerical simulations, such as those in [18,19,20,21,22,[START_REF] Ott | Classical and Quantum[END_REF]24], and for which the shear instability takes place only for very small values of the instability parameter and on timescales that are much longer than the dynamical one. When moving to larger degrees of differential rotation, that is when going to smaller values of Â, the corotation band becomes larger and larger and the shear instability can develop essentially for all values of β, merging with the dynamical bar-mode instability for β 0.25. This is exactly what has been found here. Conversely, when moving to smaller degrees of differential rotation, that is when going to higher values of Â, the corotation band becomes thinner and the shear instability can develop only for a smaller range of β. As the differential rotation is further decreased and the star tends to rotate uniformly, the corotation band width vanishes, all the models are stable to the shear instability and subject only to the dynamical bar-mode instability for β 0.25. This is indeed the case for the unstable models evolved in [6], which were purely (bar-mode) dynamically unstable and none of which had the unstable mode within the corresponding corotation band, but above it.

In essence, therefore, there should be an intermediate range of  for which the instability is absent at low β, appears at intermediate values and then disappears again at high β, thus defining an interval of values of β for which the models are unstable. To validate also this prediction we simulated a sequence of models with smaller differential rotation and  = 2. This seconds sequence has the same baryonic mass as that with  = 1, but with β in a smaller range, namely between β = 0.125 and β = 0.200 (Note that β = 0.200 is also very close to the largest value for which we could build an equilibrium model.) Unfortunately, also all of these models show an m = 2 shear instability, with the unstable frequencies falling within the corotation band, as it was for the  = 1 stars. Of course, lack of evidence is not evidence of absence and the fact that we have not found stable models within our range of values of β most likely means that these stable models have to be searched either for values of β < 0.125 (for this value of  = 2), or by moving towards higher values of Â, where the corotation band is less large.

An obvious consequence of the phenomenological scenario described by [2] and confirmed by the calculations reported here is that the idea of a low-T /|W | instability is indeed misleading and it is instead more meaningful to think of a more generic shear instability that, depending on the degree of rotation and of differential rotation, may manifest itself on timescales that are comparable with the dynamical ones (as in the cases reported here) or on much longer ones (as in the cases reported in [17]).

Another phenomenological prediction made by Watts, Andersson and Jones in [2], is that the growth times should be shorter in the center of the band and increase towards the edges (the growth times should in fact diverge at the edges). As shown in figure 5, all of the simulated models do have unstable modes well inside the corotation band and this probably explains why we observe them develop on a timescale which is comparable with the dynamical one and not on much longer timescales as in the original finding of [17]. We also find indications that the  = 2 models, which generically have longer growth times, have unstable modes which occupy regions of the corotation band that are overall more central (cf. right panel of figure 5) and thus in contrast with what expected from [2]. However, the difficulties mentioned above in computing an accurate estimate of the growth times and the intrinsic difficulties in determining what is the central part of the corotation band, prevent us from providing a more quantitative validation of this prediction.

We now switch to consider our results in the light of the other necessary condition for the onset of the shear instability discussed in [20]. As mentioned in the Introduction, it has in Figure 6. Left panel: radial profiles of the vortensity V for three representative models with small, medium and high values of β. The different symbols match the ones in figure 2 and in the left panel of figure 5, with the one for model M.1.150 reported as filled, and show the actual position of the corotation radii. Note that none of these coincides with the minimum of V; for compactness we have reported only the models with  = 1. Right panel: normalized corotation radii for the different frequencies σ i 2 presented in the left panel of figure 5 shown as a function of β; the same convention of the previous figures is used for the different symbols. fact been shown that the presence of a minimum in the profile of the vortensity is a necessary condition for a mode in corotation to be unstable [20]. The intuitive description is that the vortensity well can act as a resonant cavity inside the star, amplifying the modes that happen to lay near its minimum [32]. Indeed, the growth rate of the unstable mode is expected to depend on the location of its corotation radius with respect to the vortensity profile, being proportional to the depth of its corotation radius inside the vortensity well [20]. Note that the mere existence of a local minimum in the vortensity cannot be used as a sufficient condition for the occurrence of a shear instability. All of the unstable models considered in [6], in fact, do show a local minimum in the vortensity but do not have the unstable modes in the corotation band.

To validate whether this Newtonian condition holds also for the general-relativistic instabilities simulated here, we have computed for all the models the Newtonian vortensity, which is defined as the ratio, along the radial cylindrical coordinate, between the radial vorticity and the density

V = 2Ω + ̟Ω , ̟ ρ , ( 18
)
where ̟ is the radial cylindrical coordinate (a fully general-relativistic definition of the vortensity is also possible but more complicated to compute and not significantly different from the Newtonian one). When doing this we found that all the unstable models have vortensities with a local minimum in the star and this is shown in the left panel of figure 6, which reports the radial profiles of the vortensity V for three representative models with small, medium and high values of β. The different symbols match the ones in the left panel of figure 2 Note that none of these coincides systematically with the minimum of V, nor with the maximum of the rest-mass density. All of the corotation radii, however, do move towards larger radial positions as the rotation rate is increased, exactly as does the minimum of the vortensity and the rest-mass maximum. Although we cannot confirm that modes whose corotation radius is closer to the minimum of the vortensity have systematically shorter growth times (our data is not sufficiently accurate for this), we can compare the models with the higher degree of differential rotation  = 1 with those having a smaller degree of differential rotation  = 2. In the first case the vortensity well is considerably deeper, with ∆V/V = 1 -V min /V r=0 ∼ 0.6 -0.7, while for the latter the vortensity has shallower wells with ∆V/V = 1 -V min /V r=0 ∼ 0.02 -0.25. Since the models with  = 2 have smaller growth rates, our results indicate therefore that also the existence of a local minimum in the vortensity can be taken as a necessary condition for the development of the shear instability and that the depth of the vortensity well can be used to estimate the growth of the instability.

Conclusions

For many years the properties of rapidly rotating and self-gravitating fluids have been characterized by a complete analytic perturbative theory which provided, for instance, sufficient conditions for the development of instabilities (see, for instance, [START_REF] Chandrasekhar | Ellipsoidal Figures of Equilibrium[END_REF]55] for a collection of results). As numerical simulations have become increasingly accurate and stable on longer timescales, these predictions, both in Newtonian theory and in full general relativity, have been verified, corrected and in some cases extended. Our ability of modelling such configurations has now reached a maturity such that a number of new properties and instabilities have been "discovered" numerically, but which do not have behind fully perturbative description. The main reason for this is that most of this phenomenology is the result of physical scenarios which are much more complex than the ones investigated perturbatively in the past, e.g. non-isolated systems with high differential rotation and exchanging mass and angular momentum, which are much more difficult to treat analytically.

A most notable example of these complex and yet ubiquitous instabilities is the socalled "low-T /|W | instability", which was initially found in [17] and then reproduced in a number of other different scenarios [18,19,20,21,22,[START_REF] Ott | Classical and Quantum[END_REF]24,25]. As the phenomenological description of this instability provided by the simulations has become richer and richer, a full understanding of the mechanisms that lead to its development has lagged behind and it is presently unclear. We are therefore in a situation in which numerical simulations can probe regimes and conditions which are not yet accessible to perturbative calculations, and can guide the latter by confirming or refuting those scenarios that although possible do not find a realization in practice.

This work has followed this spirit and has used fully general-relativistic calculations of rapidly and differentially rotating neutron stars modeled with a realistic EOS to shed some light on the development of the low-T /|W | instability. In particular we have concentrated our attention on validating an indirect prediction, made by Watts, Andersson and Jones [2], who recognized the low-T /|W | instability as the manifestation of a more generic class of instabilities associated to the existence of a corotation band [30,31], the shear instabilities, and should develop for any value of the instability parameter β when sufficient amounts of differential rotation are present. This is exactly what we have found in our simulations. More specifically, we have performed simulations of sequences of neutron-star models described by a realistic SLy EOS [1] and having constant rest-mass and degrees of differential rotation, but with different amounts of rotation. In all cases considered we have found the development of a bar-mode instability growing on a dynamical timescale, even when the initial axisymmetric models were well below the critical limit for the dynamical bar-mode instability. These results, which match well the phenomenological scenario portrayed in [2], suggest therefore that the idea of a low-T /|W | instability is indeed misleading and should be replaced by the more general one of shear instability. Depending then on the degree of rotation and of differential rotation, the instability will develop on timescales that are comparable to the dynamical one (as reported here) or on much longer ones (as reported in the first low-T /|W | instability studies).

Special attention has also been paid to the properties of the unstable modes and to their position within the corotation band or the vortensity profiles. In particular, we have shown that all the unstable modes are within the corotation band of the progenitor axisymmetric model (which is the necessary condition for the development of the instability proposed by [2]) and that all of the unstable models have vortensity profiles with a local minimum (which is the necessary condition suggested by [20]). Finally, by comparing the growth times among models with different degree of differential rotation we have shown that there is a correlation, although not a strong one, between the depth of the vortensity well and growth rate of the instability, with the latter being larger for models with deeper wells.

In summary, the results presented here shed some light on several aspects of shear instabilities, but they also reveal that more work is required, for instance, to distinguish between the predictions based on the corotation band and the ones based on the vortensity well, or to establish whether in effect they just represent two different ways of expressing the same physical conditions. Clarifying these aspects will requires additional analytical and numerical modelling, and this will be part of our future research.

Figure 3 .

 3 Figure 3. Summary of the dynamics of some representative models, i.e. M.1.150, M.1.160, M.1.200, M.1.220, M.1.241 and M.1.250, with increasing values of β.For each panel the upper part reports the evolution of the bar-distortion parameter η + , while the lower part shows the evolution of the power in the different modes of the Fourier decomposition of the rest-mass density Pm. See text for details. Indicated with dotted vertical lines are the windows within which the analytic fit with the trial function 13 is made. Note that we here report only some models and with  = 1 as they are representative also of those with  = 2.

Figure 4 .

 4 Figure 4. Power spectral density (filtered using Hanning windowing) in arbitrary units of the evolution of η + for the models described in figure3(black continuous line) and of the axisymmetric mode η 0 (red dashed line). Indicated with dotted vertical lines (black, blue and magenta) are the peak frequencies reported in the left panel of figure5within the corotation band and which are used in figure6to mark the position of the corotation radii. Finally, shown with a shaded rectangular area is the corotation band for all models.

Figure 5 .

 5 Figure 5. Left panel: Position of the non-axisymmetric frequencies σ i2 for models with  = 1 and shown as a function of β; indicated with thick solid lines are the edges of the corotation band, i.e. the frequency interval between Ωc/π and Ωe/π. Shown with a filled symbol, to help distinguish it from the others, is the σ 1 2 frequency for model M.1.150 (cf. figure2 and 6). Right panel: the same as in the left one and with the same scale but for models with  = 2; the inset shows instead a magnified view. Clearly, in both panels all the unstable modes are within the band and up to three unstable modes appear for the  = 1 sequence with increasing β.

and 5 ,

 5 with the one for model M.1.150 reported as filled, and show the actual positions of the corotation radii. The right panel of the same figure, on the other hand, shows the normalized corotation radii for the different frequencies σ i 2 presented in the left panel of figure 5 as a function of β.

  Â β M/R e ρ c /10 15 ρ max /10 15 r p /r e

				R e	M	J/M 2
		(g/cm 3 )	(g/cm 3 )	(Km) (M ⊙ )
	M.1.140 1 0.140 0.245	0.989	1.011	0.562 13.20 2.19	0.681
	M.1.150 1 0.150 0.240	0.941	0.968	0.536 13.48 2.20	0.708
	M.1.160 1 0.160 0.237	0.893	0.926	0.511 13.81 2.21	0.735
	M.1.180 1 0.180 0.226	0.802	0.847	0.461 14.49 2.22	0.789
	M.1.200 1 0.200 0.217	0.712	0.773	0.413 15.27 2.24	0.844
	M.1.210 1 0.210 0.212	0.668	0.737	0.389 15.70 2.25	0.873
	M.1.220 1 0.220 0.207	0.618	0.703	0.365 16.16 2.26	0.902
	M.1.230 1 0.230 0.200	0.576	0.669	0.340 16.67 2.26	0.933
	M.1.241 1 0.241 0.194	0.522	0.633	0.312 17.29 2.27	0.968
	M.1.244 1 0.244 0.193	0.506	0.624	0.305 17.47 2.28	0.978
	M.1.247 1 0.247 0.190	0.490	0.614	0.297 17.66 2.28	0.988
	M.1.250 1 0.250 0.188	0.474	0.604	0.289 17.85 2.28	0.999
	M.2.125 2 0.125 0.241	1.143	1.143	0.642 13.22 2.17	0.668
	M.2.150 2 0.150 0.227	1.039	1.039	0.578 14.15 2.17	0.739
	M.2.175 2 0.175 0.210	0.947	0.947	0.512 15.33 2.18	0.801
	M.2.200 2 0.200 0.184	0.865	0.865	0.435 17.34 2.15	0.878

Acknowledgments

It is a pleasure to thank Anna Watts for enlightening discussions and for sharing her notes with us, as well as Nils Andersson, Ian Jones, Joel Tohline for useful comments. We are grateful to Nikolaous Stergioulas for the code to compute the initial data used in these calculations and to E. Schnetter and all the Carpet/Cactus developers. The computations were mainly performed on the Damiana cluster at the AEI and, in a part, on the Parma's ALBERT2 cluster and on the bcx cluster at CINECA. This work was supported in part by the DFG grant SFB/Transregio 7 "Gravitational-Wave Astronomy", by "CompStar", a Research Networking Programme of the European Science Foundation, and by the "Della Riccia" Foundation.


Stephan Rosswog
email: s.rosswog@jacobs-university.de

Relativistic smooth particle hydrodynamics on a given background space-time

We review the derivation of fixed-metric, relativistic smooth particle hydrodynamics (SPH) from the Lagrangian of an ideal fluid. Combining the Euler-Lagrange equations with the first law of thermodynamics, we explicitely derive evolution equations for the canonical momentum and energy. This new set of SPH equations also accounts for corrective terms that result from derivatives of the SPH smoothing kernel and that are called "grad-h" terms in non-relativistic SPH. The new equations differ from earlier formulations with respect to these corrective terms and the symmetries in the SPH particle indices while being identical in the gravitational terms.

Introduction

Relativity is a crucial ingredient in a variety of astrophysical phenomena, both due to velocities approaching the speed of light, as for example in AGN jets, and due to strong gravity shaping space-time geometry, say, near a black hole. In the recent past, substantial progress has been made in the development of special-and general relativistic numerical tools [START_REF] Marti | Numerical Hydrodynamics in Special Relativity[END_REF][START_REF] Font | Numerical Hydrodynamics in General Relativity[END_REF][START_REF] Baumgarte | Numerical relativity and compact binaries[END_REF] that can tackle a variety of pressing astrophysical problems. While most work on numerical relativistic gas dynamics has been performed in a Eulerian framework, a couple of Lagrangian smooth particle hydrodynamics (SPH) approaches exist. The first relativistic SPH formulations were developed by Kheyfets et al. [START_REF] Kheyfets | Covariant smoothed particle hydrodynamics on a curved background[END_REF] and Mann [START_REF] Mann | A relativistic smoothed particle hydrodynamics method tested with the shock tube[END_REF][START_REF] Mann | Smoothed particle hydrodynamics applied to relativistic spherical collapse[END_REF]. Shortly after, Laguna et al. [START_REF] Laguna | Smoothed particle hydrodynamics near a black hole[END_REF] developed a 3D, general-relativistic SPH code that was subsequently applied to the tidal disruption of stars by massive black holes [START_REF] Laguna | Tidal disruptions by supermassive black holes -Hydrodynamic evolution of stars on a Schwarzschild background[END_REF]. Their SPH formulation is complicated by several issues: the continuity equation contains a gravitational source term that requires SPH kernels for curved space-times. Moreover, owing to their choice of variables, the equations contain time derivatives of Lorentz factors that are treated by finite difference approximations and restrict the ability to handle shocks to only moderate Lorentz factors. The Laguna et al. formulation has recently been extended by Rantsiou et al. [START_REF] Rantsiou | Mergers of Black Hole -Neutron Star binaries. I. Methods and First Results[END_REF] and applied to neutron star black hole binaries. In a separate approach, Chow and Monaghan [START_REF] Chow | Ultrarelativistic sph[END_REF] obtained accurate specialrelativistic test results even for flow problems with large Lorentz factors. They evolve the total energy rather than the thermal energy and by applying an artificial viscosity scheme that borrows concepts from Riemann solvers [START_REF] Monaghan | SPH and Riemann Solvers[END_REF]. More recently, Siegler and Riffert [START_REF] Siegler | Smoothed Particle Hydrodynamics Simulations of Ultrarelativistic Shocks with Artificial Viscosity[END_REF] and Siegler [START_REF] Siegler | Entwicklung und Untersuchung eines Smoothed Particle Hydrodynamics Verfahrens für relativistische Strömungen[END_REF] have presented a set of equations for both the special-and general-relativistic case. Based on the conservative form of Lagrangian hydrodynamics their choice of variables bypassed many of the complications that have plagued earlier relativistic SPH formulations. They were able to simulate accurately some test problems with very large Lorentz factors. An elegant approach that is based on the discretized Lagrangian of a perfect fluid was suggested in Monaghan and Price [START_REF] Monaghan | Variational principles for relativistic smoothed particle hydrodynamics[END_REF]. Provided that the discretized Lagrangian posseses the correct symmetries, in such an approach nature's conservation laws are hard-wired into the resulting SPH equations. Moreover, the Euler-Lagrange equations determine the evolution of the fluid and do not leave much room for arbitrariness in the derivation. We are not aware of existing numerical implementations of these equations, but a different derivation that lead to a similar equation set has been tested in [START_REF] Muir | 3D Relativistic SPH[END_REF]. All of the above approaches used a fixed background metric and it was only recently that SPH has been applied to the study of flows with strong self-gravity. In these approaches Post-Newtonian approximations were used [START_REF] Ayal | Post-Newtonian Smoothed Particle Hydrodynamics[END_REF][START_REF] Faber | Post-Newtonian SPH calculations of binary neutron star coalescence: Method and first results[END_REF][START_REF] Faber | Post-Newtonian smoothed particle hydrodynamics calculations of binary neutron star coalescence. II. Binary mass ratio, equation of state, and spin dependence[END_REF][START_REF] Faber | Post-Newtonian SPH calculations of binary neutron star coalescence. III. Irrotational systems and gravitational wave spectra[END_REF], and, more recently, the conformal flatness approximation [START_REF] Isenberg | Waveless approximation theories of gravity[END_REF][START_REF] Wilson | [END_REF] has been implemented [START_REF] Oechslin | Conformally flat smoothed particle hydrodynamics application to neutron star mergers[END_REF][START_REF] Faber | Mergers of irrotational neutron star binaries in conformally flat gravity[END_REF][START_REF] Faber | Dynamical evolution of black hole-neutron star binaries in general relativity: Simulations of tidal disruption[END_REF][START_REF] Bauswein | Discriminating Strange Star Mergers from Neutron Star Mergers by Gravitational-Wave Measurements[END_REF]. In this paper we will review the derivation of the general-relativistic SPH equations from a Lagrangian, similar to the work presented in [START_REF] Monaghan | Variational principles for relativistic smoothed particle hydrodynamics[END_REF]. Extending this work, we also account for corrective extra terms due to kernel derivatives, that are called "grad-h" terms in non-relativistic SPH [START_REF] Springel | Cosmological smoothed particle hydrodynamics simulations: the entropy equation[END_REF][START_REF] Monaghan | SPH compressible turbulence[END_REF].

Relativistic SPH from a variational principle

Several authors [START_REF] Gingold | Kernel estimates as a basis for general particle methods in hydrodynamics[END_REF][START_REF] Speith | Untersuchung von Smoothed Particle Hydrodynamics anhand astrophysikalischer Beispiele[END_REF][START_REF] Springel | Cosmological smoothed particle hydrodynamics simulations: the entropy equation[END_REF][START_REF] Monaghan | SPH compressible turbulence[END_REF] have derived non-relativistic SPH equations from a variational principle. Recently, such a variational approach has also been applied to the general-relativistic case assuming a given background space-time metric [START_REF] Monaghan | Variational principles for relativistic smoothed particle hydrodynamics[END_REF]. In the following review we will adhere to a similar strategy, but also account for the generalrelativistic "grad-h"-terms. We assume that a prescribed metric g µν is known as a function of the coordinates and that the perturbations that the fluid induces to the space-time geometry can be safely neglected. We further use units in which the speed of light is equal to unity, c = 1, we adopt a metric with signature (-,+,+,+). We reserve greek letters for space-time indices from 0...3 with 0 being the temporal component, while i and j refer to spatial components and SPH particles are labeled by a, b and k. Contravariant spatial indices of a vector quantity w at particle a are denoted as w i a , while covariant ones will be written as (w i ) a . The line element and proper time are given by ds 2 = g µν dx µ dx ν and dτ 2 = -ds 2 and the proper time is related to a coordinate time t by

Θdτ = dt, (1)
where we have introduced a generalization of the Lorentz-factor

Θ ≡ 1 √ -g µν v µ v ν with v α = dx α dt . (2)
This relates to the four-velocity U ν by

v µ = dx µ dt = dx µ dτ dτ dt = U µ Θ = U µ U 0 , (3)
which is normalized to U µ U µ = -1.

The relativistic fluid Lagrangian

The Lagrangian of a relativistic fluid is given by [30]

L = -T µν U µ U ν √ -gdV, (4)
where g = det(g µν ) and T µν denotes the energy-momentum tensor of an ideal fluid without viscosity and conductivity

T µν = (ρ + P )U µ U ν + P g µν . ( 5
)
Here ρ is the energy density of the fluid as measured in the local rest frame and P is the fluid pressure. For clarity, we will write out explicity factors of c in the following lines. The energy density can be split up in a term associated with the rest mass and one given by the thermal energy contribution:

ρ = ρ rest + uρ rest /c 2 = nm 0 c 2 (1 + u/c 2 ) ( 6
)
Here n is the baryon number density in the local fluid rest frame, m 0 is the baryon mass ‡ and u = u(n, s) the specific energy, with s being the specific entropy. From now on, we will measure all energies in units of m 0 c 2 (and use again c = 1). With this convention the energy momentum tensor reads

T µν = {n(1 + u) + P } U µ U ν + P g µν . (7)
With the normalization of the four-velocity the Lagrangian can be written as

L = -n(1 + u) √ -gdV . (8)

The SPH discretization

In the following we apply the reasoning and principles behind modern SPH formulations to the general-relativistic, fixed-metric case. For a detailed account on modern SPH we refer the interested reader to a recent review [START_REF] Rosswog | Astrophysical smooth particle hydrodynamics[END_REF]. To perform practical simulations we give up general covariance and choose a particular frame ("computing frame") in which the computations are carried out. This requires a suitable transformation between computing frame quantities and quantities that are evaluated in the local rest frame of a fluid particle.

To find a SPH discretization in terms of a suitable density variable, we follow an approach similar to [START_REF] Siegler | Smoothed Particle Hydrodynamics Simulations of Ultrarelativistic Shocks with Artificial Viscosity[END_REF]. Local baryon number conservation, (U µ n); µ = 0, can be expressed as

∂ µ ( √ -gU µ n) = 0, (9)
or, more explicitely, as

∂ t (N ) + ∂ i (Nv i ) = 0, ( 10
)
where we have made use of Eq. ( 3) and have introduced the computing frame number density

N = √ -gΘn. (11)
The total conserved baryon number can then be expressed as a sum over fluid parcels with volume ∆V b located at r b , where each parcel carries a baryon number

ν b N = NdV b N b ∆V b = b ν b . (12)
Eq. ( 10) looks like the Newtonian continuity equation and we will use it for the SPH discretization process. Similar to standard SPH, one can approximate the continuum by fluid parcels ("particles"), so that a quantity f can be approximated as by

f ( r) b f b ν b N b W ( r -r b , h), (13)
where the subscript b indicates that a quantity is evaluated at a position r b and W is a suitable smoothing kernel whose width is determined by the so-called smoothing length h. If we keep all ν b constant in time, exact baryon number conservation is guaranteed and no continuity equation needs to be solved (this can be done, though, if desired).

For the kernel W , we assume that it has a compact support, so that sums only contain a local set of particles, and that it is radial, i.e.

W ( r -r b , h) = W (| r -r b |, h), so that ∂W (| r -r b |, h)/∂ r = -∂W (| r -r b |, h)/∂ r b .
If Eq. ( 13) is applied to the baryon number density N , one finds §

N a = N ( r a ) = b ν b W ( r a -r b , h a ), (14)
where we have chosen to evaluate the density estimate at r a with the local smoothing length h a . Note that we are evaluating locally smoothed quantities with flat-space kernels which assumes that the local space-time curvature radius is large in comparison to the local fluid resolution length. Such an approach is very convenient, but (more involved) alternatives to this approach also exist in the literature [START_REF] Kheyfets | Covariant smoothed particle hydrodynamics on a curved background[END_REF][START_REF] Laguna | Smoothed particle hydrodynamics near a black hole[END_REF]. The Newtonian mass density estimate can be recovered by the replacements

N b → ρ b and ν b → m b .
It is usually desirable to adjust the smoothing length locally to fully exploit the natural adaptivity of a particle method. One convenient way to do so by using the local density is

h a = η N 1/3 a , (15)
where the parameter η controls how many neighbor particles are taken into account for density estimates. Since Eqs. ( 15) and ( 14) mutually depend on each other, an iteration is required at each time step to obtain consistent values for both. For later use in the evolution equations we also provide the derivatives of the computing frame number density (with the notation

r bk = r b -r k , r bk = | r bk | and W bk (h) = W (| r b -r k |, h)): ∂N b ∂x i a = k ν k ∂W bk (h b ) ∂r bk ∂r bk ∂x i a + ∂W bk (h b ) ∂h b ∂h b ∂N b ∂N b ∂x i a . (16)
Collecting the ∂N b /∂x i a terms on both sides one obtains

∂N b ∂x i a = 1 Ω b k ν k ∂W bk (h b ) ∂r bk ∂r bk ∂x i a = 1 Ω b k ν k ∂W kb (h b ) ∂x i b (δ ba -δ ka ), (17)
where

Ω b = 1 - ∂h b ∂N b k ν k ∂W bk (h b ) ∂h b ( 18
)
is a corrective term of order unity, often called "grad-h"-term in a non-relativistic context. In a similar way, one finds the time derivative

dN b dt = 1 Ω b k ν k v i bk ∂W bk (h b ) ∂x i b , (19)
§ From now on we are dropping the distinction between "original" and approximated quantities.

where v bk = v b -v k . Motivated by Eqs. ( 11) and ( 12), we can re-write the fluid Lagrangian, Eq. ( 8), in terms of our computing frame number density N ,

L = - 1 + u Θ NdV, (20)
which leads to the general-relativistic, discretized SPH-Lagrangian

L SPH = - b ν b 1 + u Θ b . ( 21
)

The momentum equation

In the following we use the Euler-Lagrange equations

d dt ∂L ∂v i a - ∂L ∂x i a = 0, ( 22
)
to derive evolution equations for the canonical momentum and the canonical energy per baryon. The canonical momentum is

(p i ) a ≡ ∂L ∂v i a = - ∂ ∂v i a b ν b 1 + u Θ b , ( 23
)
where the velocity dependence enters directly via the generalized Lorentz-factor Θ and via

∂u b ∂v i a = ∂u b ∂n b ∂n b ∂v i a = P b n 2 b N b √ -g b ∂ ∂v i a 1 Θ b , (24)
where we have used the first law of thermodynamics, ∂u b /∂n b = P b /n 2 b , and the relation between computing frame and local rest-frame density, Eq. [START_REF] Monaghan | SPH and Riemann Solvers[END_REF]. With

∂ ∂v i a 1 Θ b = -Θ b (g iµ v µ ) a δ ab ( 25
)
one finds the canonical momentum per baryon

(S i ) a ≡ 1 ν a ∂L ∂v i a = Θ a 1 + u a + P a n a (g iµ v µ ) a = 1 + u a + P a n a (U i ) a . ( 26
)
We now need ∂L/∂x i a to obtain, via Eq. ( 22), an evolution equation for (S i ) a . The resulting terms can be split up into "hydrodynamics terms" which contain density gradients and "gravitational terms" which are proportional to gradients of the metric.

To evaluate ∂[(1 + u b )/Θ b ]/∂x i a one uses ∂ ∂x i a 1 Θ b = - U µ U ν 2Θ ∂g µν ∂x i b δ ab ( 27
)
and applies once more the first law of thermodynamics together with Eq. ( 11). The hydrodynamic contribution then becomes

∂L ∂x i a hydro = - b ν b P b √ -g b N 2 b ∂N b ∂x i a = -ν a b ν b P a √ -g a Ω a N 2 a ∂W ab (h a ) ∂x i a + P b √ -g b Ω b N 2 b ∂W ab (h b ) ∂x i a , (28)
where we have applied Eq. ( 17) and the kernel property

∂W ba (h b )/∂x i b = -∂W ab (h b )/∂x i a .
The remaining gravity terms simplify on use of Eqs. ( 11) and ( 7)

to ∂L ∂x i a gravity = ν √ -g 2N T µν ∂g µν ∂x i a . (29)
Thus, the final, general-relativistic SPH momentum equation reads

d(S i ) a dt = 1 ν a ∂L ∂x i a = - b ν b P a √ -g a Ω a N 2 a ∂W ab (h a ) ∂x i a + P b √ -g b Ω b N 2 b ∂W ab (h b ) ∂x i a + √ -g 2N T µν ∂g µν ∂x i a . (30)
The gravity contribution to this equation is identical to the one in [START_REF] Siegler | Smoothed Particle Hydrodynamics Simulations of Ultrarelativistic Shocks with Artificial Viscosity[END_REF][START_REF] Monaghan | Variational principles for relativistic smoothed particle hydrodynamics[END_REF], but the hydrodynamical terms differ both due to the presence of the grad-h terms and due to a different symmetrization in the particle indices.

The energy equation

For a suitable energy variable one can start from the canonical energy

E ≡ a ∂L ∂v i a v i a -L = a ν a v i a (S i ) a + 1 + u a Θ a (31)
to identify

e a ≡ v i a (S i ) a + 1 + u a Θ a , (32)
whose evolution equation follows from straight forward differentiation

de a dt = (S i ) a dv i a dt + v i a d(S i ) a dt + d dt 1 + u a Θ a . (33)
As we will see below, the first term cancels with a corresponding term resulting from the third term. Similar to the derivation of the momentum equation, the term d[(1 + u a )/Θ a ]/dt splits up into a term involving derivatives of the density and a term which contains derivatives of the metric tensor. Applying the first law of thermodynamics, Eq.( 11), once more and

d dt 1 √ -gΘ a = - 1 2 √ -gΘ g µν dg µν dt a - 1 √ -gΘ 2 dΘ dt a , (34)
yields the evolution equation of the thermal energy

du a dt = P a √ -g a Θ a N 2 a dN a dt - P a 2n a g µν dg µν dt - P a n a Θ a dΘ a dt . (35)
If the time derivative of the generalized Lorentz factor Θ is expressed as

dΘ a dt = Θ 3 2 v µ v ν dg µν dt + Θ 3 g µν dv µ dt v ν a (36) one finds d dt 1 + u Θ a = P a √ -g a N 2 a dN a dt -(S i ) a v i a dt - √ -g 2N T µν dg µν dt a (37)
and on using Eqs. ( 30), ( 19), (37) and

dg µν dt = ∂g µν ∂x i v i + ∂ t g µν (38)
the final general-relativistic energy equation becomes

de a dt = - b ν b P a √ -g a v i b Ω a N 2 a ∂W ab (h a ) ∂x i a + P b √ -g b v i a Ω b N 2 b ∂W ab (h b ) ∂x i a - √ -g 2N T µν ∂ t g µν a . ( 39
)
Together with an equation of state, the equations ( 14), ( 30) and (39) represent our complete and self-consistently derived set of SPH equations. The gravitational terms are identical to those of [START_REF] Siegler | Smoothed Particle Hydrodynamics Simulations of Ultrarelativistic Shocks with Artificial Viscosity[END_REF][START_REF] Monaghan | Variational principles for relativistic smoothed particle hydrodynamics[END_REF], but the hydrodynamic terms differ in both the particle symmetrization and the presence of the grad-h terms. Note that the only choices in our above derivation were the h-dependence in Eq. ( 14) and how to adapt the smoothing length. The subsequent calculation contained no arbitrariness concerning the symmetry in particle indices, everything followed stringently from the first law of thermodynamics and the Euler-Lagrange equations. Another important point to note is that the derived energy equation, Eq. (39), does not contain destabilizing [START_REF] Norman | Why ultrarelativistic numerical hydrodynamics is difficult[END_REF] time derivatives of Lorentz factors on the RHS -in contrast to the thermal energy equation (35) and to earlier SPH formulations [START_REF] Laguna | Smoothed particle hydrodynamics near a black hole[END_REF]. For a practical simulation involving shocks, Eqs. ( 30) and (39) need to be augmented by artificial viscosity terms similar to the special-relativistic case [START_REF] Chow | Ultrarelativistic sph[END_REF][START_REF] Siegler | Smoothed Particle Hydrodynamics Simulations of Ultrarelativistic Shocks with Artificial Viscosity[END_REF][START_REF] Rosswog | Conservative, special-relativistic smooth particle hydrodynamics[END_REF]. In order to couple to a generic space-time evolution code, the partial derivative ∂ t g µν in Eq.(39) would need to be re-written in terms of the extrinsic curvature, see e.g. [START_REF] Alcubierre | Introduction to 3+1 Numerical Relativity[END_REF].

The special-relativistic limit

In the special-relativistic limit we can neglect the gravitational terms in Eqs. ( 30) and (39). In addition in flat space-time with Cartesian coordinates one has √ -g → 1 and Θ → γ, and Eq. ( 11) becomes N = γn, which simply expresses the increase in the computing frame number density N with respect to the local fluid rest frame density n due the Lorentz contraction. The momentum and energy equations reduce in this limit to

  d S a dt   SR = - b ν b P a Ω a N 2 a ∂W ab (h a ) ∂ r a + P b Ω b N 2 b ∂W ab (h b ) ∂ r a (40)
and

de a dt SR = - b ν b P a v b Ω a N 2 a • ∂W ab (h a ) ∂ r a + P b v a Ω b N 2 b • ∂W ab (h b ) ∂ r a (41)
which are the equations recently derived and successfully tested in a special-relativistic context [START_REF] Rosswog | Conservative, special-relativistic smooth particle hydrodynamics[END_REF].

Summary

In this article we have reviewed in detail the derivation of the general-relativistic SPH equations form a variational principle for the case of a given background metric. We have used the computing frame baryon number density in the SPH discretization process and derived evolution equations for the canonical energy and momentum. If the latter are chosen as numerical variables, no derivatives of Lorentz factors occur in the evolution equations. As in non-relativistic SPH, one has the choice between calculating (here: computing frame baryon number) densities via a particle summation or by evolving the continuity equation. The derived equation set also contains corrective terms from smoothing kernel derivatives, so-called grad-h terms. The hydrodynamic parts of the equation set differ from earlier general-relativistic SPH formulations by the presence of these corrective terms and by the symmetry in the particle indices, the contributions from the gravitational field are identical to those found in earlier work. In the specialrelativistic limit in Cartesian coordiates, the equations reduce to a recently derived equation set that has been successfully tested in large number of benchmark problems [START_REF] Rosswog | Conservative, special-relativistic smooth particle hydrodynamics[END_REF]. Future applications of this SPH formulation will include tidal disruptions of stars by massive black holes.

‡ The appropriate baryon mass depends on the neutron to proton ratio, i.e. on the nuclear composition of the considered fluid.