Abstract

We consider the influence of potential quark matter existing at high densities in neutron star (NS) interiors on gravitational waves (GWs) emitted in a binary NS merger event. Two types of equations of state (EoSs) at zero temperature are used – one describing pure nuclear matter and the other nuclear matter with a phase transition to quark matter at very high densities. Binary equilibrium sequences close to the innermost stable circular orbit (ISCO) are calculated to determine the GW frequencies just before the merger. It is found that the effects of the EoSs begin to play a role when gravitational masses are larger than M≃ 1.5 M. The difference in the GW frequency at the ISCO increases by up to ≃10 per cent for the maximum mass permitted by the EoSs. We then perform three-dimensional hydrodynamic simulations for each EoS while varying the initial mass and determine the characteristic GW frequencies of the merger remnant. The implications of the presence of quark matter show up mainly in the collapse behaviour of the merger remnant. If the collapse does not take place immediately after the merger, we find a phase difference between the two EoSs in the post-merger GW signal. We also compare the GW frequencies emitted by the remnant of the merger to values obtained from simulations using a polytropic EoS and find an imprint of the non-constant adiabatic index of our EoSs. All calculations are based on the conformally flat approximation to general relativity and the GW signal from the merger simulation is extracted up to quadrupole order.

Introduction

Binary neutron star (NS) mergers belong to the strongest gravitational wave (GW) sources for interferometer-type GW observatories, for example LIGO (Abramovici et al. 1992), VIRGO (Bradaschia et al. 1990), GEO600 (Lück 1997) and TAMA (Ando et al. 2001). After a long inspiral process which lasts millions of years, the final merger phase takes place over a time-scale of just milliseconds. The onset of the merger begins when the two companions become dynamically unstable near the innermost stable circular orbit (ISCO) and mass transfer starts. As the process has no definite symmetries, three-dimensional (3D) hydrodynamic simulations are necessary. In addition, owing to the compactness of the system, the effects of general relativity must be considered and dealt with as accurately as possible.

It has been pointed out that the GW frequency, fGW, is related to the compactness, (M/R), of a NS by fGW∼ (M/R)3/2 near the ISCO, meaning that the GW frequency just before a merger carries information to the radius of the NS. Hence, the GW signal may constrain the equation of state (EoS) of high-density matter, (see e.g. Lai & Wiseman 1996; Rasio & Shapiro 1999). It has also been found recently that a NS formed after a merger may be supported by rapid and differential rotation even if its mass exceeds 60 per cent of the maximum mass of a single non-rotating NS, and GWs are emitted due to non-axisymmetric and quasi-radial oscillations of the remnant for longer than the dynamical time-scale (Shibata & Uryū, 2000, 2002; Shibata, Taniguchi & Uryū 2003). If these oscillations persist for a long time, an integrated GW may be detectable, carrying information on the high-density matter.

In addition to general relativistic gravity, various types of physics should be taken into account when dealing with the binary NS merger problem, for example the nuclear forces summarized in the EoS, neutrino physics and magnetic fields. So far, however, owing to the complexities involved, investigations have concentrated on either the relativistic aspects of the problem using a simple EoS or the microphysical aspects of treating gravity in a Newtonian framework. Pioneering work in this field has been carried out by Oohara & Nakamura and is detailed in (Oohara & Nakamura 1999, and references therein). Relativistic aspects have been considered using the post-Newtonian approximation (Ayal et al. 2001; Faber & Rasio 2002, and references therein), the conformally flat (CF) approximation (Mathews & Wilson 2000; Oechslin, Rosswog & Thielemann 2002; Faber, Grandclément & Rasio 2003, and references therein) and a fully general relativistic treatment (Shibata & Uryū, 2000, 2002; Shibata et al. 2003, and references therein). Microphysical improvements have been applied by (Ruffert & Janka 2001, and references therein) and (Rosswog & Davies 2002, and references therein) using different EoSs (Lattimer & Swesty 1991; Shen et al. 1998) and a leakage scheme to account for neutrino emission after the merger. For a review on the topic see Rasio & Shapiro (1999).

In this paper, we focus on two EoSs for high-density NS matter – one describing pure nuclear matter and the other nuclear matter with a transition to the quark phase at very high densities. The reason for investigating these EoSs is that quark phase transition may be one of the most dramatic phenomena that changes the compactness of NSs, and therefore its influence may be clearly observed in the GW spectrum of an inspiralling binary system. Following our previous work (Oechslin et al. 2002), we consider the merger problem in the CF approximation.

The paper is organized as follows. In Section 2, we summarize the formalism to solve the Einstein and the relativistic hydrodynamic equations, and describe the numerical implementation, the choice of the EoS and the initial conditions. In Section 3, we present our results and, finally, we draw our conclusions in Section 4.

Fundamental Ingredients

In this section we describe the numerical methods and the physics on which our simulations are based. We consider a general relativistic fluid whose internal properties are described by a given EoS. Two sets of equations have to be solved simultaneously. On the one hand, we consider the relativistic hydrodynamic equations governing fluid motion and on the other hand the Einstein equation of general relativity determining the space–time metric and therefore the gravitational interaction. On top of this, an EoS, which closes the system of hydrodynamic equations, has to be given as the input.

Relativistic hydrodynamics and an approximation for general relativistic gravity

In the 3 + 1 decomposition of space–time, the metric ds2=gμν dxμ dxν can be written as 

1
formula
where α is the lapse function, βi is the shift vector and γij is the spatial metric. Space–time quantities are decomposed with respect to foliation using the hypersurface normal defined as nμ=−α∂μt with nμnμ=−1 and a projection tensor defined as γμν=gμν+nμnν, whose spatial component agrees with the spatial metric.

The Einstein field equations can be written as a set of two evolution equations, i.e. 

2
formula
and 
3
formula
and two constraint equations, 
4
formula
and 
5
formula
for the dynamic variables γij and Kij, the extrinsic curvature of the hypersurface (see e.g. Baumgarte & Shapiro 2003). Here, Rij is the Ricci tensor and ∇ the covariant derivative associated with γij. The stress energy tensor Tμν is decomposed into ρE=nμnνTμν, the matter energy density, jiiμnνTμν, the matter momentum density and SijiμγjνTμν, the spatial projection of the stress energy tensor. Finally, R, S and K denote the traces of Rij, Sij and Kij, respectively.

In the following discussion, we consider a perfect fluid with a matter stress energy tensor 

6
formula
Here, ρ refers to the rest-mass density, h= 1 +p/ρ+ε to the specific relativistic enthalpy, ε to the specific internal energy, uμ to the four velocity and p to the fluid pressure. Then, 
7
formula
and 
8
formula
The Lorentz factor Wu0 can be calculated using the normalization condition uμuμ=−1 and so 
9
formula

Isenberg proposed a waveless approximation to general relativity, in which he truncates some terms in the Einstein equation written in the Arnowitt, Deser and Misner (ADH) formalism to deduce an elliptic-type formalism (Isenberg 1978). Later, the same set of equations were rediscovered by (Wilson, Mathews & Marronetti 1996) and widely used to solve single or binary NSs problems (e.g. Oechslin et al. 2002) as well as binary black hole (BH) systems (Grandclement, Gourgoulhon & Bonazzola 2002). The Isenberg–Wilson–Mathews theory has its own Hamiltonian as pointed out by Isenberg himself (also discussed in Friedman, Uryū & Shibata 2002).

In this approximation, the CF conditionγij4δij for the spatial geometry and K=∂tK= 0 are imposed, where δij is the flat three-metric. It may be viewed that the second condition is a choice between a temporal gauge condition and the maximal slicing condition, while the first condition is an approximation and partly a spatial gauge choice (details can be found in Baumgarte et al. 1998). The remaining metric variables (ψ, α, βi) do not satisfy all components of the Einstein field equation consistently. We pick the constraints and the trace of the evolution equation for Kij as equations for these five variables, which leads to the same set of equations derived from Isenberg's Hamiltonian mentioned already.

The approximation leads to a considerable simplification of the Einstein equations, since all metric equations can be written in elliptic form.

The trace of the evolution equation for Kij together with the maximal slicing condition, K= 0, leads to the following equation for the lapse: 

10
formula

The trace-free part of the evolution equation for γij together with the CF condition, γij4δij, gives 

11
formula
where ∂i is a derivative with respect to the coordinate associated with δij and graphic relate to the shift βi through the flat metric as graphic. The Hamiltonian constraint provides an equation for the conformal factor, i.e. 
12
formula
Finally, with the momentum constraint we obtain an equation for the shift vector 
13
formula
where ∂iijj. Using the definition 
14
formula
equation (13) splits into two simpler parts 
15
formula
and 
16
formula
While the trace-free components of the evolution equation for γij (2) are used to relate Kij and the metric variables as in equation (11), its trace part and the trace-free components of the evolution equation for Kij (3) are dropped by the CF approximation. The trace of equation (2) can be used to check the accuracy and reliability of the CF approximation (Dimmelmeier, Font & Müller 2002; Oechslin 2003)

The relativistic hydrodynamic equations can be written as 

17
formula
 
18
formula
 
19
formula
where 
20
formula
and 
21
formula
In the case of conformal flatness, the system reduces to 
22
formula
 
23
formula
and 
24
formula
with 
25
formula
 
26
formula
and 
27
formula
Here, we have cast the system into a Lagrangian formulation which is better suited for an implementation on a computer using a Lagrangian scheme like the smoothed particle hydrodynamics (SPH) method (cf. Section 2.2). The conserved hydrodynamic variables are graphic whereas the primitive variables are (ρ, vi, ε).

Numerical implementation

To solve the system of hydrodynamic equations (22)–(24) we use the SPH method (Monaghan, & Gingold 1983; Benz 1990), which is widely used in astronomical and astrophysical simulations. The original form of SPH that solves the Newtonian hydrodynamic equations is presented in Benz (1990). We list here the generalized version which is appropriate for the already mentioned relativistic hydrodynamic equations. For details, see Oechslin et al. (2002). The continuity equation turns into a relationship between ρ* and the rest masses of the individual particles 

28
formula
where mb is the rest mass of particle b and Wab=W(|rarb|, hSPH) denotes the weight given by the standard spherical spline Kernel function W(r, hSPH). The rest masses mb and the smoothing length hSPH are initially chosen to fit an initial spatial distribution of ρ*. The pressure gradient contained in the momentum equation (24) is calculated as in standard SPH but replacing ρ by ρ*, i.e. 
29
formula
The energy equation (25) contains a velocity divergence term similar to Newtonian SPH, i.e. 
30
formula
The additional terms in the momentum and energy equation that arise from the gravitational interaction are evaluated by first solving the equations for the metric variables α, ψ and βi on a overlaid grid (see discussion below) and by mapping back those quantities on to the particles by second-order interpolation. The total time derivative of ln (αu0ψ6)a in the energy equation is evaluated using second-order finite differencing in time.

On top of these equations, an artificial viscosity (AV) scheme is implemented with time-dependent viscosity parameters (Morris & Monaghan 1997) which have only significant values in the presence of shocks. The AV produces an additional viscous pressure which is added to the physical fluid pressure (Siegler & Riffert 2000).

The field equations (10), (12), (15) and (16) are discretized on a grid which covers the matter distribution. The evaluation proceeds generally in three steps as follows.

  • (1)

    First, we calculate the source corresponding to the potential. The hydrodynamic quantities involved and defined on the SPH particles are transferred on to the grid by assigning SPH-interpolated values to the grid points, i.e. 

    31
    formula

  • (ii)

    The potential is then obtained by solving the corresponding Poisson equation. The ψ and the αψ equations are closely linked and are therefore solved in two iteration steps using the solution from the previous time-step as a guess. The remaining four equations determining the shift vector are also closely coupled and are solved simultaneously. All Poissonian equations are solved using a full multi-grid solver.

  • (iii)

    All needed derivatives are calculated on the grid using finite differencing.

  • (iv)

    Finally, the potential and its derivatives are mapped back on to the SPH particle distribution using a triangular-shaped cloud method (Hockey & Eastwood 1992) obtained using particle mesh codes. This is equivalent to second-order interpolation.

The boundary conditions and extensions of the solution beyond the grid are obtained using a multipole expansion of the source term. For details we refer to Oechslin et al. (2002).

Like the Newtonian and the first-order post-Newtonian approximations, the CF approximation does not include gravitational radiation and its reaction by construction. We therefore have to add an additional GW extraction scheme and a radiation reaction force which accounts for the angular momentum and energy carried away by GWs. The waveform in a transverse traceless gauge is extracted at the slow-motion limit (Thorne 1980; Wilson et al. 1996) and up to quadrupole order using 

32
formula
where 
33
formula
is the mass quadrupole of the system. The radiation reaction force is chosen in a similar way to the original Burke–Thorne expression, Frri=σ∂iV with V=−1/5xixjQ(5)ij, but its strength σ is chosen such that the resulting angular momentum loss by back-reaction, graphic, reproduces the expression, 
34
formula
of the slow-motion GW extraction approximation up to a difference of the order of ≃5 per cent. Since graphic, we can obtain the angular momentum change using 
35
formula
The expressions graphic, the numerical angular momentum error and graphic, the total angular momentum change per unit time, are obtained by evaluating equation (35) before adding and after adding the back-reaction force terms, respectively. As graphic has an oscillating behaviour around zero graphic depends on the binary orbit. Therefore, to practically determine the value of σ, we take an average value, graphic, for about a half orbit, and compare it with graphic computed from the quadrupole formula (34) that is averaged for the same orbit. With model B1 we obtain graphic per cent at the ISCO, where the average value is taken over one orbital period, while at a certain period of time graphic may become comparable to graphic. As the ratio graphic is negligible, the radiation reaction force is clearly driving the inspiral process, although the predicted time to the merger must be taken with care due to the approximative radiation reaction scheme and the spatial CF assumption used in our initial conditions and simulations. Fully relativistic inspiral simulations Miller, Gressman & Suen (2003) are intrinsically consistent with respect to radiation reaction dissipation but they still depend on the outer boundary location.

Equation of state

To investigate the influence of quark matter on a NS merger, we consider two EoSs – an EoS describing pure nuclear matter (‘hadronic EoS’) and an EoS describing nuclear matter with a phase transition to quark matter at very high densities (‘hybrid EoS’).

The hadronic EoS is realized using the non-linear σ−ω model in the relativistic mean field approximation with the TM1 parameter set (Sugahara & Toki 1994) which is motivated by a least-squares fit to experimental results including stable and unstable nuclei. At densities above ρ≳ 1014 g cm−3 this is a good approximation. For lower values, when inhomogeneous nuclear matter appears and the non-linear σ−ω model is no longer valid, other approximations have to be used (e.g. Shen et al. 1998). In our case, we append the polytropic EoS, p=κρΓ, where κ and Γ are adjusted to ensure smoothness of pressure and internal energy. We choose the transition density for the polytropic EoS to be 2 × 1014 g cm−3, which leads to Γ≃ 2.86; similar to other common, realistic EoSs in this density regime. The hybrid EoS is obtained by combining the hadronic EoS with a MIT bag model (Chodos et al. 1974) using a variable pressure phase transition construction (Glendenning 1992). We use massless up and down quarks and a bag constant of 90 MeV fm−3. The resulting hybrid EoS then describes three physical phases as follows.

  • (i)

    A pure hadronic phase below ≈5 × 1014 g cm−3≈ 1.8ρ0 where the EoS coincides with the hadronic EoS. Here, ρ0 := 2.8 × 1014 g cm−3 is the nuclear saturation density. The stiffness in this region varies betweenΓ≈ 3 and Γ≈ 2.5.

  • (ii)

    A quark–hadron mixed phase between ≈5 × 1014 g cm−3 and ≈1015 g cm−3≈ 3.5ρ0 where both quark and hadrons are present. In this phase transition region, the EoS substantially softens when Γ≈ 1–1.5.

  • (iii)

    A pure quark phase above ≈1015 g cm−3 where the MIT–bag model EoS with a similar adiabatic index but a lower absolute pressure compared to the hadronic EoS is applied.

In both EoSs, we neglect any temperature effects, i.e. we set T= 0. The redundant internal energy information from the EoS is dropped and we consider the pressure as a function of the density alone. This is a good approximation in the high-density regime where thermal effects caused by pressure are small. At lower densities (e.g. in the disc around the merger remnant) the thermal component is certainly not negligible. However, in this work we are concentrating on the high-density merger remnant and thus we do claim to be accurate with this approximation although thermal effects could have an effect via neutrino losses.

Initial conditions

To construct the initial conditions, we assume the binary to be in a quasi-equilibrium state. This assumption is not very well satisfied in the very vicinity of the ISCO (Miller 2003) and we therefore choose slightly wider initial orbits. In our models the binaries initially have no radial velocity and need about half an orbit to start the inspiral. We apply the method of Uryū & Eriguchi (2000) to construct the initial configuration. This method takes an EoS as the input and then solves the hydrostatic equation for an irrotational velocity field together with the Einstein field equations in the CF approximation. A similar method has been developed by Bonazzola, Gourgoulhon & Marck (1999; see also Gourgoulhon et al. 2001, Taniguchi & Gourgoulhon 2002). We map the output of the above method, the hydrodynamic quantities and the conformal factor on to a distribution of SPH particles. Finally, this SPH distribution is relaxed, to avoid any spurious inter-particle forces, using a braking term, 

36
formula
where virr is the given initial irrotational velocity field of the initial mass distribution. Using this method, we produce a set of irrotational quasi-equilibrium configurations, both varying the EoS and the gravitational mass M, where M refers to the mass of a single NS in isolation. For M, we use values between M= 1.2 M which is at the lower limit of the observational range and M= 1.75 M which is close to the maximal gravitational NS mass of Mmax= 1.78 M of the hybrid EoS. The models are summarized in Table 1. The letter (A–E) in the model label indicates the mass while the number (1, 2) indicates the EoS. All orbits are taken slightly outside the ISCO, about 12 ∼ 25 per cent in coordinate separation. This results in smaller lag angles and oscillations in the initial phase since tidal forces are weaker with larger separation distances. In Table 2, we summarize the properties which result as a consequence our choice of initial parameters (detailed in Table 1).

1

Input parameters of the models under consideration using the hadronic and hybrid EoSs. M denotes the gravitational mass of one NS in isolation and C is the compactness. Both the masses and the compactness refer to one binary component in isolation. NParticle denotes the number of SPH particles and NGrid is the number of grid points. In model E1, 973 gridpoints are used, but at the onset of collapse, the resolution is enhanced to 1293. Models with a ‘1’ in their label are calculated with the hadronic EoS and models with a ‘2’ used the hybrid EoS.

1

Input parameters of the models under consideration using the hadronic and hybrid EoSs. M denotes the gravitational mass of one NS in isolation and C is the compactness. Both the masses and the compactness refer to one binary component in isolation. NParticle denotes the number of SPH particles and NGrid is the number of grid points. In model E1, 973 gridpoints are used, but at the onset of collapse, the resolution is enhanced to 1293. Models with a ‘1’ in their label are calculated with the hadronic EoS and models with a ‘2’ used the hybrid EoS.

2

Selected properties of our initial models resulting from the input parameters shown in Table 1. M0 is the rest mass of one single NS, P the orbital period, fGW,0 the corresponding GW frequency, d0 the initial orbital separation and R the radius of one single isolated NS measured in Schwarzschild coordinates. Note that the GW frequencies are lower than those in Fig. 2 which are taken at the ISCO.

2

Selected properties of our initial models resulting from the input parameters shown in Table 1. M0 is the rest mass of one single NS, P the orbital period, fGW,0 the corresponding GW frequency, d0 the initial orbital separation and R the radius of one single isolated NS measured in Schwarzschild coordinates. Note that the GW frequencies are lower than those in Fig. 2 which are taken at the ISCO.

1

Comparison of the nuclear EoSs under consideration. The hybrid EoS (dashed curve) has a phase transition region (1.8ρ0 <ρ <3.5ρ0) where the adiabatic indices are substantially lower, followed by a quark phase which has a similar stiffness to nuclear matter at such densities.

1

Comparison of the nuclear EoSs under consideration. The hybrid EoS (dashed curve) has a phase transition region (1.8ρ0 <ρ <3.5ρ0) where the adiabatic indices are substantially lower, followed by a quark phase which has a similar stiffness to nuclear matter at such densities.

Results

In this section, we consider the final evolution of the binary NS from a circular orbit configuration to one single object, a transient NS or a BH. This evolution takes place on a dynamical time-scale and is triggered by a radiation reaction. In Section 3.1, the EoS effects on the GW frequency from the final inspiral phase just before a merger is discussed and in the latter subsections the results of the dynamical evolution of the binary NS from a stable circular orbit to merger are presented.

The impact of the equation of state on the quasi-equilibrium binary

By considering the compactness (M/R) in Table 1, we observe that the properties of the models become EoS dependent for the most massive cases of D1, D2 and E1, E2, i.e. whenM≳ 1.5 M. For lower stellar masses, the NSs do not reach the phase transition density, ρt := 1.8ρ0, in the centre and an EoS effect cannot be observed. It has been shown that the orbital frequency and the GW frequency of a binary system depends on the compactness of its components (e.g. Lai & Wiseman 1996; Uryū, Shibata & Eriguchi 2000; Faber et al. 2002). Therefore we expect that an EoS effect may be seen in the GW frequency at the final inspiral stage in our models.

Using the method of Uryū & Eriguchi (2000); Uryū et al. (2000), we systematically investigated the locations of the ISCO, beyond which the binary become dynamically unstable, and the GW frequency at the ISCO. As shown in Fig. 2(a), an EoS effect becomes important above M= 1.5 M. The reason is that at this mass, the central density reaches the phase transition density of mixed matter in the hybrid case and a dense core of mixed matter is formed. For even larger masses, the transition density of pure quark matter at ρ≃ 3.5ρ0 is reached and a quark core is formed. Owing to this, a NS with the hybrid EoS becomes more compact than one with the hadronic EoS as shown in Fig. 2(b), which is reflected in the difference in the GW frequency. Near the maximal NS mass of the hybrid EoS, the GW frequency of the hybrid EoS binary is up to 10 per cent larger than that of the hadronic EoS binary. Such a change in the tendency for the GW frequency to increase with respect to increasing mass at M≳ 1.5 M may suggest a drastic change in the EoS such as a quark phase transition.

2

(a) GW frequencies depending on the gravitational mass in isolation at the ISCO. (b) Compactness of an isolated NS with respect to gravitational mass.

2

(a) GW frequencies depending on the gravitational mass in isolation at the ISCO. (b) Compactness of an isolated NS with respect to gravitational mass.

Dynamical evolution of the neutron star merger

We now turn to the merger phase that follows the quasi-stationary inspiral when the two NSs cross the ISCO and become dynamically unstable. To first illustrate the overall dynamics, we plot in Fig. 3 the density distribution for model B1 together with the velocity field in the corotating frame, whose angular velocity is determined as 

37
formula
where i runs over all particles and ωi denotes the angular velocity of the individual particles. The stars are counterrotating relative to the corotating frame and tidal lag angles are developing (snapshot 1). At merger, this counterrotation is first turned into a shear motion at the contact layer and then into a pattern of vortex rolls owing to the growing Kelvin–Helmholtz instability (snapshots 2 and 3). Finally, we end up with a differentially rotating, non-axisymmetric merger remnant (snapshot 4). In Fig. 4, we concentrate on the evolution of the merger remnant during the first few milliseconds by linearly plotting the density contour; this gives better visualization of the high-density central parts (the lowest density contour is at 5 × 1013 g cm−3.) In this figure we can see that ∼2 ms after merger, the core actually consists of two small subcores which are leftovers of the original cores of the NSs. Also visible is the fact that the two subcores still carry a small counterrotating motion – a direct consequence of the irrotating setup – which steadily dissipates until we end up in with a nearly axisymmetric, differentially rotating configuration. It is reported from finite difference simulations (e.g. Shibata & Uryū 2002; Shibata et al. 2003) that this twin-core pattern persists for a longer time than our simulations last. This may be a consequence of the numerical viscosity in our code.

3

Snapshots at selected times of the baryonic density distribution and the velocity field of model B1 during pre-merger and merger evolution. The density is plotted logarithmically in units of g cm−3 and the velocity field is plotted in the corotating frame.

3

Snapshots at selected times of the baryonic density distribution and the velocity field of model B1 during pre-merger and merger evolution. The density is plotted logarithmically in units of g cm−3 and the velocity field is plotted in the corotating frame.

4

Same as Fig. 3 but for post-merger evolution. The density is plotted linearly in units of 1013 g cm−3.

4

Same as Fig. 3 but for post-merger evolution. The density is plotted linearly in units of 1013 g cm−3.

Equation of state differences in maximal density

To investigate the different consequences of the EoSs in detail, we first consider the evolution of the maximum density during the merger phase. We split the whole merger event into pre-merger evolution where the binary still consists of two tidally stretched objects and post-merger evolution where a merger remnant, a NS or BH is forming. As the SPH particle density ρa has no direct physical meaning, the actual density ρ(x) has to be calculated as a statistical average over the various densities of the particles ρa. To find the maximum density, we first determine density values at selected gridpoints and then look for their maximum. Both statistical averaging and finite grid spacing introduce a small amount of noise. In general, the evolution ofρmax, plotted in Fig. 5 for our simulations, shows a slow decrease during the pre-merger phase when the two stars become tidally stretched and reaches its minimum during the actual merger. This is the case when the GW luminosity either approaches its maximum and the two NSs are maximally tidally stretched or just a bit later when a new merger remnant has already formed but a twin-core structure is still present (see Section 3.2). After the minimum, we can distinguish three possible evolutionary scenarios of ρmax, depending on the collapse behaviour of the remnant. If the remnant can be stabilized by pressure and centrifugal forces, ρmax slowly increases and then becomes constant. The second scenario is the delayed collapse where ρmax first slowly increases on several dynamical time-scales. This fragile equilibrium finishes with a final collapse, which is shown by the steep increase in ρmax. The final scenario is the immediate collapse just after the merger. Here, ρmax increases without delay on a dynamical time-scale.

5

Evolution of the maximal density measured in units of the nuclear saturation density ρ0 := 2.8 × 1014 g cm−3. The origin of the time axis has been shifted to the GW luminosity maximum. The solid lines correspond to the hadronic models and the dashed lines to the hybrid models.

5

Evolution of the maximal density measured in units of the nuclear saturation density ρ0 := 2.8 × 1014 g cm−3. The origin of the time axis has been shifted to the GW luminosity maximum. The solid lines correspond to the hadronic models and the dashed lines to the hybrid models.

Differences in the merger dynamics caused by the EoS, appear when ρmax reaches the phase transition density ρt. At this point the hybrid EoS and therefore the dynamical evolution of the hybrid models separates from that of the hadronic EoS. If the individual stars are more massive than M≃ 1.5 M, their central density exceeds the phase transition density ρt. In this case, EoS differences happen before the merger. This makes a very small effect which even disappears close to the merger due to tidal distortion if M is just at the threshold of ≃1.5 M as in model D1/D2. Alternatively, if M is close to the maximal gravitational mass of the hybrid EoS as in model E1/E2, the EoS difference in ρmax increases to up to ∼60 per cent. Differences can also be measured for global quantities such as the angular velocity of the binary, the compactness and especially the GW frequency (cf. Section 3.1). If the initial mass is less than M≃ 1.5 M a change in the EoS can only be seen after the merger as part of a different evolutionary scenario of the merger remnant. Typically, the remnant collapses immediately or after a couple of dynamical time-scales in the hybrid cases while it settles down to a transient NS in the hadronic cases. Among the hadronic models, only the very massive model E1 collapses immediately. In model D1, the maximal density slowly increases but does not really settle down – a sign that this object might eventually collapse. However, this might well be a numerical effect as angular momentum is transported from the core to the outer layers and differential rotation is slowly converted into uniform rotation. The amount which is transported is of the order of 40–50 per cent in the very centre during the whole evolution, despite a very good overall angular momentum conservation, by either numerical viscosity or gravitational interaction between the non-axisymmetric core and the outer layers. However, it is still possible to stabilize a remnant with a gravitational mass of nearly 3 M with an EoS having a maximal gravitational mass of 2.2 M for about 12 ms ≃ 820 Mtot, where Mtot= 2.98 M is the initial gravitational mass for this model. This stabilization effect due to differential rotation has been pointed out in Shibata & Uryū (2000); Baumgarte, Shapiro & Shibata (2000). The less massive models C1, B1 and A1 lead to transient NSs which do not collapse on a hydrodynamic time-scale. With the hybrid EoS models, the models collapse over different time-scales. While models E2 and D2 collapse immediately, C2 collapses with a delay of ∼9 ms. Model B2 does not collapse within the simulation time but it also shows a continuously growing maximum density. However, these collapse time-scales are likely to be dominated by the previously mentioned angular momentum transport.

The gravitational wave signal

A second indicator of EoS effects is the GW signal. As the waveform is sensitive to dynamical mass motions we expect that all the above maximal density differences are reflected in the GW signal mainly in the form of different frequencies. In Fig. 6 we plot the waveforms of all models sorted according to initial mass. Model E1/E2, the most massive one, is the only one which shows a significant pre-merger EoS difference in the waveform. This is because the above-mentioned maximal density difference translates via the compactness and the binary angular velocity into a GW frequency difference which amounts to ∼10 per cent at the ISCO (cf. Fig. 2) while it disappears during the actual merger phase. The next less massive model, D1/D2, does not show any EoS differences in GW frequency at the ISCO and further during the pre-merger phase, but we expect, for larger binary distances and therefore smaller tidal interaction that GW frequency differences could be seen. The more interesting part of this model is the different collapse behaviour and therefore the totally different waveform. While the hybrid model D2 only produces a short, high-frequency burst before the collapse, the hadronic model emits a long wavetrain which decreases slowly in amplitude. Both models C1/C2 and B1/B2 emit a quasi-periodic (QP) GW signal which is characteristic of the rotation and oscillation mode in the merger remnant. EoS differences in the waveform do not appear until a considerable mixed matter core has formed in the hybrid case. As a consequence, the first part of the post-merger GW signal, which is the strongest in our simulations, will not be affected by any differences in the EoS. However, when the mixed matter core becomes large enough, an accumulating phase shift in the waveform becomes clearly visible. In model C1/C2 the shift adds up to more than half a period before the hybrid remnant collapses, whereas in the B1/B2 model the shift is only very small, but visible.

6

Gravitational waveforms of all models sorted according to their initial mass. The origin of the time axis has been shifted to the GW luminosity maximum. The solid lines correspond to the hadronic models and the dashed lines to the hybrid models.

6

Gravitational waveforms of all models sorted according to their initial mass. The origin of the time axis has been shifted to the GW luminosity maximum. The solid lines correspond to the hadronic models and the dashed lines to the hybrid models.

Gravitational wave spectra

The Fourier spectrum of the GW waveform is plotted in Fig. 7. The most massive model, E1/E2, shows its EoS differences at pre-merger, i.e. in the frequency domain around 1 kHz. Here the spectrum produced by the pre-merger hybrid waveform is significantly stronger than that of the hadronic waveform, since the amplitude of the wave strain h scales roughly as (r/M) h∼ (M/R), where r is the distance from the source. This is consistent with simulations using quasi-equilibrium sequences (Faber et al. 2003). There is no high-frequency contribution from this model as both remnants collapse immediately after the merger. Note that we miss the waveform emitted by the ringing of the resulting BH. However, the expected frequency of this signal lies in the range of 5–10 kHz (Shibata & Uryū 2002; Leaver 1985) and we can clearly separate out the BH waveform.

7

Fourier spectrum of the waveform h+ for models B1/B2 to E1/E2 sorted according to their initial mass. The frequency is in kHz. The solid lines correspond to the hadronic models, the dotted lines to the spectra of the corresponding truncated waveforms (see text) and the dashed lines to the hybrid models.

7

Fourier spectrum of the waveform h+ for models B1/B2 to E1/E2 sorted according to their initial mass. The frequency is in kHz. The solid lines correspond to the hadronic models, the dotted lines to the spectra of the corresponding truncated waveforms (see text) and the dashed lines to the hybrid models.

The other models show their EoS differences via different frequency peaks which result from the QP waveform emitted post-merger. The difference is most obvious in model D1/D2. For hadronic EoS model D1 a distinct Fourier peak is clearly visible, whereas hybrid EoS model D2 only leads to a much weaker and broader Fourier peak as the remnant collapses very soon after merging producing a much shorter post-merger GW signal. In model C1/C2, both the hadronic and the hybrid models have strong QP peaks and the peak of model C2 is slightly shifted to higher frequencies owing to the more compact remnant. The same can be seen for model B1/B2. However, the effects here are smaller and measurement may be difficult. The EoS effects are only small because the spectrum is dominated by the first GW burst after merging, which is insensible to EoS differences.

We may compare the obtained QP frequencies fQP to values from a fully relativistic calculation (Shibata & Uryū 2002) where an irrotational initial velocity field with a polytropic EoS is considered. Two peaks are found in the spectrum, 

38
formula
and 
39
formula
for an adiabatic index of Γ= 2.25; and 
40
formula
and 
41
formula
for Γ= 2. Here, fQE is the GW frequency of the binary at the ISCO. Shibata & Uryū (2002) point out, that fQP2/fQE is not very dependent on the initial compactness of the stars but the higher peak, fQP2/fQE, depends on Γ so that softer EoSs lead to higher values for fQP2/fQE. As we use realistic EoSs with a non-constant Γ in our simulations, we expect the fQP2/fQE ratio to vary with increasing compactness (Zang, Centrella & McMillan 1996; Faber & Rasio 2001; Oechslin et al. 2002). Indeed our values for fQP2/fQE are strongly dependent on the initial mass M and compactness as can be seen in Table 3. From the behaviour of the fQP2/fQE ratio we can deduce that the EoS softens for more compact models, i.e. the EoS softens at higher densities. It might be interesting to calculate the fQP2/fQE− (M/R) relationship for various other realistic EoSs on the market. Evaluating this relationship might give a detailed insight into the Eos of the NS in the nuclear regime above 1014 g cm−3.

3

Ratio of fQP2 to fQE for all models where a transient NS forms or a delayed collapse happens. The values for fQE are determined from quasi-equilibrium models and are taken at the ISCO.

3

Ratio of fQP2 to fQE for all models where a transient NS forms or a delayed collapse happens. The values for fQE are determined from quasi-equilibrium models and are taken at the ISCO.

8

Relationship between the ratio f QP2/ f QE and the compactness C = (M/R). The diamonds correspond to hadronic EoS models A1 − D1 and the values for Γ =2 and Γ =2.25 are taken from Shibata & Uryū (2002).

8

Relationship between the ratio f QP2/ f QE and the compactness C = (M/R). The diamonds correspond to hadronic EoS models A1 − D1 and the values for Γ =2 and Γ =2.25 are taken from Shibata & Uryū (2002).

One might expect that fQP2 is overestimated in our simulations due the already mentioned numerical viscosity effects in the code which lead to a slowly contracting merger remnant and therefore to the emission of GWs with higher frequencies. To check that this is not the case, we have also calculated the spectra of the truncated wavesignals which only contain the first few oscillations of the post-merger signals. The resulting QP peaks are reduced in strength but their positions remains unchanged within the uncertainty of the peak width (see Fig. 7).

Conclusions

In this work, we have considered the impact of the EoS on the inspiral and merger dynamics of a binary NS coalescence. We have performed equilibrium sequence studies to investigate the GW frequencies around the ISCO. We then carried out dynamical simulations of the merging event using a 3D SPH code. We chose a hadronic EoS based on the relativistic mean field approximation and a hybrid EoS obtained by combining the hadronic EoS with a MIT–bag model EoS in the high-density regime above ρt= 1.8ρ0. The GW frequency at the ISCO depends strongly on the initial mass of the binary. For a generic mass of 2 × 1.4 M, we find a GW frequency of ∼900 Hz for both EoSs. EoS differences become important for models with initial gravitational masses larger than M∼ 1.5 M. At this mass range, the central density of the individual stars reach the phase transition ρt included in the hybrid EoS model. At masses close to the maximum gravitational mass of the hybrid EoS, M= 1.78 M, the relative difference in fGW becomes as large as ∼10 per cent.

The maximum density evolution during the merger depends sensitively on the fact whether the phase transition density is reached and at what time this happens. For very massive models such as model E1/E2, which has an initial mass of M= 1.75 M, this is already fulfilled by the individual companion stars in the pre-merger phase. Less massive models with masses below M≲ 1.5 M cross this threshold later in the post-merger phase when the matter is contracting to form a transient NS or a BH. The hadronic EoS models all form a transient NS except the very massive model, E1, which collapses immediately. In model D1, the maximum density does not converge after 12 ms but rises slowly and continuously. We think that this is the effect of considerable angular momentum transportation from the remnant core to the outer layers which implies continuous conversion of the differential rotation pattern into uniform rotation. All hybrid EoS models collapse on time-scales which are highly dependent on the model mass. While the massive models E2 and D2 collapse immediately, models C2 and B2 collapse after a contraction of the remnant on several dynamical time-scales and a considerable decrease in the degree of differential rotation. This indicates that the collapse is driven here by angular momentum transport caused, for example, by numerical viscosity. The collapse dynamics of these models will be investigated in future work. The emitted gravitational waveforms are very dependent on the initial mass of the model. The most massive model E1/E2 shows a large difference prior to merger both in frequency and amplitude. This is a consequence of the different levels of compactness in models E1 and E2 owing to the different EoSs involved. The same effect is not visible in the less massive models because the difference in compactness owing to the EoSs vanishes. A second aspect of the GW signal is the waveform emitted by the remnant. If the remnant does not collapse immediately to a BH, a QP waveform with a frequency of 2-3 kHz is emitted. The ratio of this frequency to the GW frequency at the ISCO, fQP2/fQE, depends sensitively on the model and therefore on the mass or compactness, (M/R). As fQP2/fQE is fairly constant for an EoS with constant stiffness, we interpret this result as being a consequence of the varying stiffness of the EoS used. Hence, the relationship fQP2/fQE− (M/R) is characteristic for the behaviour of the stiffnesses of EoSs. A measurement of this quantity would provide important additional information on EoSs.

It has been suggested in Hughes (2002) that measurements of high-frequency GWs require a network of broadband detectors combined with narrowband detectors that have a good level of sensitivity in the high-frequency domain. Such measurements of high-frequency GWs may become feasible in the future and confirm or deny the existence of a phase transition in NS star matter at high densities.

Acknowledgments

Computations were done at the Department for Physics and Astronomy at the University of Basel. RO would like to thank T. Janka for discussions and helpful comments. KU would like to thank M. Shibata and J. L. Friedman for discussions.

RO, GP and FKT are funded by the Schweizerischer Nationalfonds under grant no 2000-061822.00.

References

Abramovici
A.
et al
,
1992
,
Sci
 ,
256
,
325
Abramovici
A.
et al
,
1996
,
Phys. Lett. A
 ,
218
,
157
Ando
M.
et al
,
2001
,
Phys. Rev. Lett.
 ,
86
,
3950
Arnowitt
R.
Deser
S.
Misner
C. W.
,
1962
, in
Witten
L.
,
Gravitation: An Introduction to Current Research
 .
Wiley
, NY , p.
226
Ayal
S.
Piran
T.
Oechslin
R.
Davies
M.
Rosswog
R.
,
2001
,
ApJ
 ,
550
,
846
Baumgarte
T. W.
Cook
G. B.
Scheel
M. A.
Shapiro
S. L.
Teukolsky
S. A.
,
1998
,
Phys. Rev. D
 ,
57
,
7299
Baumgarte
T. W.
Shapiro
S. L.
Shibata
M.
,
2000
,
Astrophys. J.
 ,
528
,
L29
Benz
W.
,
1990
, in
Buchler
J. R.
, ed.,
The Numerical Modelling of Non-linear Stellar Pulsations: Problems and Prospects
 .
Kluwer
, Dordrecht , p.
269
Bonazzola
S.
Gourgoulhon
E.
Marck
J.-A.
,
1999
,
Phys. Rev. Lett.
 ,
82
,
892
Bradaschia
C.
et al
,
1990
,
Nucl. Instrum. Methods A
 ,
289
,
518
Chodos
A.
Jaffe
R. L.
Johnson
K.
Thorn
C. B.
Weisskopf
V. F.
,
1974
,
Nucl. Phys Rev. D
 ,
9
,
3471
Dimmelmeier
H.
Font
J. A.
Müller
E.
,
2002
,
A&A
 ,
388
,
917
Faber
J. A.
Rasio
F. A.
,
2001
,
Phys. Rev. D.
 ,
63
,
044012
Faber
J. A.
Rasio
F. A.
,
2002
,
Phys. Rev. D.
 ,
65
,
084042
Faber
J. A.
Grandclément
P.
Rasio
F. A.
Taniguchi
K.
,
2002
,
Phys. Rev. Lett.
 ,
89
,
231102
Faber
J. A.
Grandclément
P.
Rasio
F. A.
,
2003
,
Phys. Rev. D
 , submitted
Friedman
J. L.
Uryū
K.
Shibata
M.
,
2002
,
Phys. Rev. D
 ,
65
,
064035
Glendenning
N. K.
,
1992
,
Phys. Rev. D
 ,
46
,
1274
Grandclement
P.
Gourgoulhon
E.
Bonazzola
S.
,
2002
,
Phys. Rev. D
 ,
65
,
044021
Hockney
R. W.
Eastwood
J. W.
,
1994
,
Computer Simulation Using Particles
 . Inst. Phys., London
Hughes
S. A.
,
2002
,
Phys. Rev. D
 ,
66
,
102001
Isenberg
J.
Nester
J.
,
1980
, in
Held
A.
, ed.,
General Relativity and Gravitation
 , Vol.
1
.
Plenum Press
, New York , p.
23
Lai
D.
Wiseman
A. G.
,
1996
,
Phys. Rev. D
 ,
54
,
3958
Lattimer
J. M.
Swesty
F. D.
,
1991
,
Nucl. Phys., A
 ,
535
,
331
Leaver
E. W.
,
1985
,
Proc. R. Soc. London A
 ,
402
,
285
Lück
H.
,
1997
,
Class. Quant. Grav.
 ,
14
,
1471
Mathews
G. R.
Wilson
J. R.
,
2000
,
Phys. Rev. D
 ,
61
,
127304
Miller
M.
,
2003
, (gr-qc/0305024)
Miller
M.
Gressman
P.
Suen
W. M.
,
2003
,
Phys. Rev. D
 , in press (gr-qc/0312030)
Monaghan
J.
Gingold
R.
,
1983
,
J. Comp. Phys.
 ,
52
,
374
Morris
J.
Monaghan
J.
,
1997
,
J. Comp. Phys.
 ,
136
,
41
Oechslin
R.
,
2003
,
PhD thesis
 , Univ. Basel
Oechslin
R.
Rosswog
S.
Thielemann
F. K.
,
2002
,
Phys. Rev. D
 ,
65
,
103005
Oohara
K.
Nakamura
T.
,
1999
,
Prog. Theor. Phys. Suppl.
 ,
136
,
270
Rasio
F. A.
Shapiro
S. L.
,
1999
,
Class. Quant. Grav.
 ,
16
,
R1
Rosswog
S.
Davies
M. B.
,
2002
,
MNRAS
 ,
3
,
481
Rosswog
S.
Liebendörfer
M.
,
2003
,
MNRAS
 ,
342
,
673
Ruffert
M.
Janka
H. T.
,
2001
,
A&A
 ,
380
,
544
Shen
H.
Toki
H.
Oyamatsu
K.
Sumiyoshi
K.
,
1998
,
Prog. Theor. Phys.
 ,
100
,
1013
Shibata
M.
Uryū
K.
,
2000
,
Phys. Rev. D
 ,
61
,
064001
Shibata
M.
Uryū
K.
,
2002
,
Prog. Theor. Phys.
 ,
107
,
265
Shibata
M.
Taniguchi
K.
Uryū
K.
,
2003
,
Phys. Rev. D
 ,
68
,
084020
Siegler
S.
Riffert
H.
,
2000
,
ApJ
 ,
531
,
1053
Sugahara
Y.
Toki
H.
,
1994
,
Nucl. Phys. A
 ,
579
,
557
Thorne
K. S.
,
1980
,
Rev. Mod. Phys.
 ,
52
,
299
Uryū
K.
Eriguchi
Y.
,
2000
,
Phys. Rev. D
 ,
61
,
124023
Uryū
K.
Shibata
M.
Eriguchi
Y.
,
2000
,
Phys. Rev. D
 ,
62
,
104015
Wilson
J. R.
Mathews
G. J.
Marronetti
P.
,
1996
,
Phys. Rev. D
 ,
54
,
1317
Zang
X.
Centrella
J. M..
McMillan
S. L. W.
,
1996
,
Phys. Rev. D
 ,
50
,
7261

Author notes

Present address: Max-Planck Institut für Astrophysik, Karl Schwarzschild-Strasse 1, 85741 Garching, Germany.