## 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, *f*_{GW}, is related to the compactness, (*M*/*R*)_{∞}, of a NS by *f*_{GW}∼ (*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 d*s*^{2}=*g*_{μν} d*x*^{μ} d*x*^{ν} can be written as

^{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.

and and two constraint equations, and for the dynamic variables γ_{ij}and

*K*, the extrinsic curvature of the hypersurface (see e.g. Baumgarte & Shapiro 2003). Here,

_{ij}*R*is the Ricci tensor and ∇ the covariant derivative associated with γ

_{ij}_{ij}. The stress energy tensor

*T*

_{μν}is decomposed into ρ

_{E}=

*n*

^{μ}

*n*

^{ν}

*T*

_{μν}, the matter energy density,

*j*=γ

^{i}^{i}

_{μ}

*n*

_{ν}

*T*

^{μν}, the matter momentum density and

*S*=γ

_{ij}_{iμ}γ

_{jν}

*T*

^{μν}, the spatial projection of the stress energy tensor. Finally,

*R*,

*S*and

*K*denote the traces of

*R*,

_{ij}*S*and

_{ij}*K*, respectively.

_{ij}In the following discussion, we consider a perfect fluid with a matter stress energy tensor

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, and The Lorentz factor

*W*=α

*u*

^{0}can be calculated using the normalization condition

*u*

_{μ}

*u*

^{μ}=−1 and so

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γ_{ij}=ψ^{4}δ_{ij} for the spatial geometry and *K*=∂_{t}*K*= 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 *K _{ij}* 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 *K _{ij}* together with the maximal slicing condition,

*K*= 0, leads to the following equation for the lapse:

The trace-free part of the evolution equation for γ_{ij} together with the CF condition, γ_{ij}=ψ^{4}δ_{ij}, gives

_{i}is a derivative with respect to the coordinate associated with δ

_{ij}and relate to the shift β

^{i}through the flat metric as . The Hamiltonian constraint provides an equation for the conformal factor, i.e. Finally, with the momentum constraint we obtain an equation for the shift vector where ∂

^{i}=δ

^{ij}∂

_{j}. Using the definition equation (13) splits into two simpler parts and While the trace-free components of the evolution equation for γ

_{ij}(2) are used to relate

*K*and the metric variables as in equation (11), its trace part and the trace-free components of the evolution equation for

_{ij}*K*(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)

_{ij}The relativistic hydrodynamic equations can be written as

where and In the case of conformal flatness, the system reduces to and with and 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 whereas the primitive variables are (ρ,*v*, ε).

^{i}### 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

where*m*is the rest mass of particle

_{b}*b*and

*W*=

_{ab}*W*(|

*r*_{a}−

*r*_{b}|,

*h*

_{SPH}) denotes the weight given by the standard spherical spline Kernel function

*W*(

**,**

*r**h*

_{SPH}). The rest masses

*m*and the smoothing length

_{b}*h*

_{SPH}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. The energy equation (25) contains a velocity divergence term similar to Newtonian SPH, i.e. 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 (α

*u*

^{0}ψ

^{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.

- (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

where is the mass quadrupole of the system. The radiation reaction force is chosen in a similar way to the original Burke–Thorne expression,*F*

^{rr}

_{i}=σ∂

_{i}

*V*with

*V*=−1/5

*x*

_{i}x_{j}Q^{(5)}

_{ij}, but its strength σ is chosen such that the resulting angular momentum loss by back-reaction, , reproduces the expression, of the slow-motion GW extraction approximation up to a difference of the order of ≃5 per cent. Since , we can obtain the angular momentum change using The expressions , the numerical angular momentum error and , 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 has an oscillating behaviour around zero depends on the binary orbit. Therefore, to practically determine the value of σ, we take an average value, , for about a half orbit, and compare it with computed from the quadrupole formula (34) that is averaged for the same orbit. With model B1 we obtain per cent at the ISCO, where the average value is taken over one orbital period, while at a certain period of time may become comparable to . As the ratio 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 ρ≳ 10^{14} 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 × 10^{14} 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 × 10

^{14}g cm^{−3}≈ 1.8ρ_{0}where the EoS coincides with the hadronic EoS. Here, ρ_{0}:= 2.8 × 10^{14}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 × 10

^{14}g cm^{−3}and ≈10^{15}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 ≈10

^{15}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,

where

*v*_{irr}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

*M*

_{max}= 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).

## 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. when*M*_{∞}≳ 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.

### 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

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 × 10

^{13}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.

### 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.

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 *M*_{tot}, where *M*_{tot}= 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.

#### 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.

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 *f*_{QP} 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,

*f*

_{QE}is the GW frequency of the binary at the ISCO. Shibata & Uryū (2002) point out, that

*f*

_{QP2}/

*f*

_{QE}is not very dependent on the initial compactness of the stars but the higher peak,

*f*

_{QP2}/

*f*

_{QE}, depends on Γ so that softer EoSs lead to higher values for

*f*

_{QP2}/

*f*

_{QE}. As we use realistic EoSs with a non-constant Γ in our simulations, we expect the

*f*

_{QP2}/

*f*

_{QE}ratio to vary with increasing compactness (Zang, Centrella & McMillan 1996; Faber & Rasio 2001; Oechslin et al. 2002). Indeed our values for

*f*

_{QP2}/

*f*

_{QE}are strongly dependent on the initial mass

*M*

_{∞}and compactness as can be seen in Table 3. From the behaviour of the

*f*

_{QP2}/

*f*

_{QE}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

*f*

_{QP2}/

*f*

_{QE}− (

*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 10

^{14}g cm

^{−3}.

One might expect that *f*_{QP2} 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 *f*_{GW} 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, *f*_{QP2}/*f*_{QE}, depends sensitively on the model and therefore on the mass or compactness, (*M*/*R*)_{∞}. As *f*_{QP2}/*f*_{QE} 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 *f*_{QP2}/*f*_{QE}− (*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.