SYSTEMATICS OF DYNAMICAL MASS EJECTION, NUCLEOSYNTHESIS, AND RADIOACTIVELY POWERED ELECTROMAGNETIC SIGNALS FROM NEUTRON-STAR MERGERS

A. Bauswein; S. Goriely; H.-T. Janka

doi:10.1088/0004-637X/773/1/78

1. INTRODUCTION

Neutron star (NS) merger events are among the most promising candidates for the first direct measurement of a gravitational-wave signal with the upcoming Advanced LIGO and VIRGO interferometric instruments (Acernese et al. 2006; Harry & LIGO Scientific Collaboration 2010), and they are considered as the likely origin of short gamma-ray bursts and their afterglows as a consequence of ultrarelativistic, collimated outflows (see, e.g., Soderberg et al. 2006; Nakar 2007; Berger 2011; Kann et al. 2011; Fong et al. 2012). Moreover, they are possible sources of different kinds of electromagnetic signals in the precursor of the merging and in its aftermath as a consequence of magnetohydrodynamical effects, magnetospheric interactions, relativistic matter outflows, or NS crust phenomena (Lipunov & Panchenko 1996; Vietri 1996; Li & Paczyński 1998; Hansen & Lyutikov 2001; Troja et al. 2010; Shibata et al. 2011; Nakar & Piran 2011; Tsang et al. 2012; Kyutoku et al. 2012; Zhang 2013; Gao et al. 2013; Piro 2012; Lai 2012; Palenzuela et al. 2013; Metzger & Berger 2012). Thermal emission produced by hot ejecta gas, for example, may cause potentially observable optical transients (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010b; Metzger & Berger 2012), and the interaction of the ejecta cloud with the circumstellar medium is expected to create radio flares that might be detectable for periods of years (Nakar & Piran 2011; Piran et al. 2013; Rosswog et al. 2013). Observations of such signals could help to pinpoint the exact celestial locations of NS mergers (thus, e.g., supporting the analysis of data taken by gravitational-wave detectors), and repeated measurements of signals that can be unambiguously linked to NS mergers would help to constrain the still highly uncertain rate of such events in the local universe.

During the merging of two NSs a small fraction of the system mass, typically 0.1%–1%, can become gravitationally unbound and can be ejected on the dynamical timescale of milliseconds (Ruffert et al. 1997; Rosswog et al. 1999, 2000; Ruffert & Janka 2001; Oechslin et al. 2007; Rosswog et al. 2013; Piran et al. 2013; Rosswog 2013; Hotokezaka et al. 2013). Because such material is likely to possess a high neutron excess, it has been proposed as a possible site for the creation of the heaviest, neutron-rich elements, which are formed by the rapid neutron capture process (r-process; Lattimer et al. 1977; Eichler et al. 1989) (similarly, also NS–black hole mergers were suggested as sources of r-process matter; Lattimer & Schramm 1974, 1976). The radioactive decay of these freshly synthesized r-process nuclei should heat the ejecta and thus lead to an optical transient (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010b; Roberts et al. 2011; Goriely et al. 2011). The properties of such events depend on the fraction of the material that can be converted to radioactive species. Moreover, the peak luminosity, the timescale to reach the emission peak, and the effective temperature at the radiation maximum, as well as the radio brightness that accompanies the deceleration of the expelled gas during its coasting in the stellar environment, depend sensitively on the ejecta mass and expansion velocity. Detailed hydrodynamical merger models are needed to calculate these quantities and to determine the nucleosynthesis conditions in the unbound material.

Concerning their role as sources of heavy elements, binary NSs collisions have recently moved into the focus of interest because the astrophysical sources of the r-process elements have not been identified yet and core-collapse supernova simulations continue to be unable to yield the extreme conditions for forming the heaviest neutron-rich nuclei (Arcones et al. 2007; Hoffman et al. 2008; Janka et al. 2008; Roberts et al. 2010; Hüdepohl et al. 2010; Fischer et al. 2010; Wanajo et al. 2011; Arcones & Janka 2011; Arcones & Martínez-Pinedo 2011). (For reviews on r-process nucleosynthesis and an overview of potential sites, see, e.g., Arnould et al. 2007; Thielemann et al. 2011; Banerjee et al. 2011; Winteler et al. 2012.) In contrast to the situation for supernovae, investigations with growing sophistication have confirmed NS merger ejecta as viable sites for strong r-processing (Freiburghaus et al. 1999; Goriely et al. 2005; Arnould et al. 2007; Metzger et al. 2010b; Roberts et al. 2011; Goriely et al. 2011; Korobkin et al. 2012).

However, despite this promising situation a variety of aspects need to be clarified before the question can be answered whether NS mergers are a major source or even the dominant source of heavy r-process elements. On the one hand the merger rate and its evolution during the Galactic history are still subject to considerable uncertainties (see, e.g., Abadie et al. 2010 for a compilation of recent estimates), and it is unclear whether NS mergers can explain the early enrichment of the Galaxy by r-process elements as observed in metal-deficient stars (Argast et al. 2004). On the other hand it remains to be determined how much mass is ejected in merger events depending on the binary parameters and, in particular, depending on the incompletely known properties of the equation of state (EOS) of NS matter. It also needs to be understood which fraction of the ejecta is robustly converted to r-process material and whether the final abundances are always compatible with the solar element distribution, which agrees amazingly well with the r-process abundance pattern in metal-poor stars for atomic numbers Z ∼ 55–90 (see, e.g., Sneden et al. 2008).

Newtonian as well as relativistic studies showed that the mass ratio has a significant effect on the amount of matter that can become unbound (Janka et al. 1999; Rosswog et al. 1999, 2000, 2013; Ruffert & Janka 2001; Oechslin et al. 2007; Roberts et al. 2011; Goriely et al. 2011; Piran et al. 2013; Korobkin et al. 2012; Rosswog 2013; Hotokezaka et al. 2013). Such investigations, however, were performed only with a few exemplary models for high-density matter in NSs (Rosswog et al. 2000; Oechslin et al. 2007; Goriely et al. 2011; Hotokezaka et al. 2013) or even only with a single NS EOS (Roberts et al. 2011; Piran et al. 2013; Rosswog et al. 2013; Korobkin et al. 2012; Rosswog 2013), although the importance of the nuclear EOS for a quantitative assessment of the dynamical mass ejection can be concluded from published calculations (e.g., Goriely et al. 2011). These calculations, however, also suggest that the nuclear abundance pattern produced by r-processing in the ejecta may be largely insensitive to variations of the conditions in the ejecta.

It is important to note that quantitatively reliable information on the ejecta masses and their dependence on the binary and EOS properties require general relativistic (GR) simulations. Newtonian results in the literature (Rosswog et al. 1999; Janka et al. 1999; Ruffert & Janka 2001; Roberts et al. 2011; Korobkin et al. 2012; Rosswog et al. 2013; Piran et al. 2013; Rosswog 2013) exhibit significant quantitative and qualitative differences compared to relativistic models (Oechslin et al. 2007; Goriely et al. 2011; Hotokezaka et al. 2013). Newtonian calculations tend to overestimate the ejecta masses in general (Rosswog et al. 1999, 2013; Janka et al. 1999; Ruffert & Janka 2001; Roberts et al. 2011; Korobkin et al. 2012; Piran et al. 2013; Rosswog 2013). This can be understood because of several facts. First, the structure of NSs in GR is considerably more compact than that of Newtonian stars. For instance, an NS with a gravitational mass of 1.35 M_☉ described by the LS220 EOS (Lattimer & Swesty 1991) possesses a circumferential radius of 12.6 km, whereas the corresponding Newtonian star has 14.5 km. Second, GR gravity is stronger and the merging of two NSs is therefore more violent. The difference can be expressed in terms of the gravitational binding energy of a nucleon on the surface of the considered 1.35 M_☉ NSs, which is ∼200 MeV in the GR case compared to only ∼130 MeV for the Newtonian model. Third, GR forces merger remnants beyond a mass limit to collapse to black holes on a dynamical timescale. Such an effect cannot be tracked by Newtonian models. These differences are of direct relevance for the collision dynamics and the possibility to unbind matter from the inner and outer crust regions of the merging NSs.

It is the purpose of this paper to explore the influence of the high-density EOS on the ejecta properties in a systematic way, i.e., we will determine ejecta masses and the nucleosynthesis outcome for a large set of NS matter models, applying them in relativistic NS merger simulations. Most of these EOSs were already employed in our previous works (Bauswein & Janka 2012; Bauswein et al. 2012). They were chosen such that they provide as completely as possible a coverage of the possibilities for NS properties (expressed by corresponding mass–radius-relations) which are compatible with present observational constraints (e.g., the 1.97 M_☉ NS discovery of Demorest et al. 2010, and more recently 2.01 M_☉ by Antoniadis et al. 2013) and theoretical understanding (Lattimer & Prakash 2010, 2007; Steiner et al. 2010; Hebeler et al. 2010). In our study we will focus on symmetric 1.35–1.35 M_☉ systems and will compare them with asymmetric 1.2–1.5 M_☉ mergers. Because population synthesis models (Belczynski et al. 2008) and pulsar observations (Thorsett & Chakrabarty 1999; Zhang et al. 2011) suggest that the double NS population is strongly dominated by systems of nearly equal-mass stars of about 1.35 M_☉ each, the average NS merger event can be well represented by a 1.35–1.35 M_☉ configuration, and a clarification of the EOS dependence of ejecta masses, r-process yields, and properties of electromagnetic counterparts of NS mergers seems to be more important than a wide variation of binary parameters. Nevertheless, we will also present results of a more extended survey of binary mass ratios and total masses for some representative EOSs.

In our work we will exclusively concentrate on NS–NS mergers, but the discussed phenomena should play a role also for NS–black hole coalescence (Janka et al. 1999; Lee 2000; Rosswog et al. 2004, 2013; Rosswog 2005, 2013; Faber et al. 2006; Foucart et al. 2013; Piran et al. 2013) and eccentric NS mergers (East et al. 2012; East & Pretorius 2012; Rosswog et al. 2013; Rosswog 2013). However, while the existence of double NS systems is established by observations, progenitors of NS–black hole and eccentric NS mergers have not been observed yet and the rates of such types of events are even more uncertain than those of coalescing binary NSs. In investigating the latter we will only consider the phase of dynamical mass ejection between about the time when the two NSs collide until a few milliseconds later. During this phase hydrodynamical and tidal forces (shock compression, pressure forces, gravitational interaction) are responsible for the mass shedding of the merging objects. Once the remnant has formed, however, differential rotation is expected to strongly amplify the magnetic fields (e.g., Price & Rosswog 2006; Anderson et al. 2008; Liu et al. 2008; Giacomazzo et al. 2011) and viscous energy dissipation is likely to provide additional heating, enhancing the neutrino emission that accompanies the secular evolution of the post-merger configuration (Ruffert & Janka 1999; Setiawan et al. 2004; Dessart et al. 2009; Shibata & Sekiguchi 2012; Fernández & Metzger 2013). As a consequence the merger remnant will experience mass loss due to neutrino energy deposition in the near-surface regions (Ruffert & Janka 1999; Setiawan et al. 2004; Dessart et al. 2009; Wanajo & Janka 2012; Shibata & Sekiguchi 2012) (similar to the neutrino-driven wind of proto-neutron stars emerging from stellar core collapse) and due to magnetohydrodynamical outflows. Both mechanisms will add ejecta to the mass stripped during the dynamical interaction of the system components, but the details of the secular evolution and the associated mass loss will be very sensitive to the EOS-dependent stability properties of the merger remnant, i.e., to the question whether the remnant is a hypermassive NS (see Baumgarte et al. 2000 for a definition) or whether and when it collapses to a black hole–torus system. These questions lie beyond the scope of the present work.

Our paper is organized as follows. In Section 2 a brief summary of the numerical methods and microphysics ingredients of our NS merger simulations is given. In Section 3 we present our results for the relation between dynamical mass loss and NS (EOS) properties, provide a detailed description of the mass-loss dynamics in our relativistic models (drawing comparisons to Newtonian results), discuss the influence of an approximate treatment of thermal effects in the EOS, and evaluate the mass ejection for three selected, representative EOSs in merger simulations for a wider space of binary masses and mass ratios in order to determine the population-integrated mass loss. In Section 4 we describe results of nuclear network calculations performed for a subset of our merger models and draw conclusions on the Galactic merger rate and the production of long-lived radioactive species (²³²Th, ²³⁵U, ²³⁸U) used for stellar nucleocosmochronometry. Finally, we present values for the heating efficiency of the merger ejecta by radioactive decays of the nucleosynthesis products and apply them in Section 5 to estimate the properties (peak luminosity, peak timescale, effective temperature at the maximum luminosity) of the optical transients that can be expected from the expanding merger debris. We also briefly discuss the implications of our simulations for radio flares. Finally, a summary and conclusions follow in Section 6.

2. NUMERICAL MODEL AND EQUATIONS OF STATE

The simulations of our study are performed with a relativistic smooth particle hydrodynamics (SPH) code, i.e., the hydrodynamical equations are evolved in a Lagrangian manner (Oechslin et al. 2002, 2007; Bauswein et al. 2010a). The Einstein field equations are solved imposing conformal flatness of the spatial metric (Isenberg & Nester 1980; Wilson et al. 1996), and a gravitational-wave backreaction scheme is used to account for energy and angular momentum losses by the emission of gravitational radiation (Oechslin et al. 2007). The code evolves the conserved rest-mass density ρ*, the conserved specific momentum $\tilde{u}_i$ , and the conserved energy density τ, whose definitions evolve the metric potentials and the "primitive" hydrodynamical quantities, i.e., the rest-mass density ρ, the coordinate velocity v_i, and the specific internal energy epsilon . The system of relativistic hydrodynamical equations is closed by an EOS, which relates the pressure P = P(ρ, T, Y_e) and the specific internal energy epsilon = epsilon (ρ, T, Y_e) to the rest-mass density ρ, the temperature T, and the electron fraction Y_e. The temperature is obtained by inverting the specific internal energy epsilon = epsilon (ρ, T, Y_e) for given ρ and Y_e. Changes of the electron fraction are assumed to be slow compared to the dynamics (see, e.g., Ruffert et al. 1997), and the initial electron fraction, which is defined by the neutrinoless beta-equilibrium of cold NSs, is advected according to (dY_e/dt) = 0 (d/dt defines the Lagrangian, i.e., comoving, time derivative).

The EOS of NS matter is only incompletely known, and numerical studies rely on theoretical prescriptions of high-density matter. This work surveys a representative sample of 40 microphysical EOSs, which have been derived within different theoretical frameworks and make different assumptions about the composition of high-density matter and the description of nuclear interactions. Most of the employed EOSs are listed in Bauswein et al. (2012), where details can be found, while some new models are introduced below. Because of the one-to-one correspondence between the EOS and the mass–radius relation of nonrotating NSs, it is convenient to characterize EOSs by the resulting stellar properties. Stellar quantities as integral properties of an EOS are in particular useful to classify the dynamics of NS mergers and the accompanying gravitational-wave signals (Bauswein & Janka 2012; Bauswein et al. 2012). For this reason we will adopt the same approach also for this investigation. Considering for instance NSs with a gravitational mass of 1.35 M_☉, the stellar radii R_1.35 vary from 10.13 km to 15.74 km for the different EOSs of our sample. The maximum mass M_max of nonrotating NSs obtained for these EOSs ranges from 1.79 M_☉ to 3.00 M_☉. In terms of their stellar properties the employed EOSs of our study show a large variation (see the mass–radius relations in Figure 4 of Bauswein et al. 2012). Note that we do not apply any selection procedure for choosing the EOSs, except that we require a maximum mass above ≈1.8 M_☉. This limit is chosen because of the firm discovery of pulsars with gravitational masses of 1.97 ± 0.04 M_☉ (Demorest et al. 2010) and 2.01 ± 0.04 M_☉ (Antoniadis et al. 2013). EOSs which yield a maximum mass below this limit are practically excluded by this observation. Nevertheless, we accept them (at least down to M_max ≈ 1.8 M_☉) for our investigation because we expect that at densities relevant in a typical NS merger these models still provide a viable description of high-density matter (see Bauswein et al. 2012). Note that compared to our previous study in Bauswein et al. (2012) we extend our EOS survey by including also the models TM1, TMA, NL3, DD2, SFHO, and SFHX of Hempel & Schaffner-Bielich (2010), Hempel et al. (2012), and Steiner et al. (2012), relying on the interactions described in Sugahara & Toki (1994), Toki et al. (1995), Lalazissis et al. (1997), Typel et al. (2010), and Steiner et al. (2012). Moreover, we include the BSk20 and BSk21 EOSs of Goriely et al. (2010). The maximum masses M_max resulting for these EOSs are 2.21 M_☉, 2.02 M_☉, 2.79 M_☉, 2.42 M_☉, 2.06 M_☉, 2.13 M_☉, 2.16 M_☉, and 2.28 M_☉ (order as listed above), while the radii of cold 1.35 M_☉ NSs are 14.49 km, 13.86 km, 14.75 km, 13.21 km, 11.92 km, 11.98 km, 11.74 km, and 12.54 km, respectively. From our sample of EOSs in Bauswein et al. (2012) we do not consider the SKA EOS (because of its restriction to densities above 1.7 × 10⁹ g cm⁻³) and EOSs which are not compatible with the pulsar observation of Demorest et al. (2010) and Antoniadis et al. (2013) and directly form a black hole after merging. We also exclude absolutely stable strange quark matter. We refer to Bauswein et al. (2009) for the particular implications of ejecta from strange quark star mergers.

Only 12 out of the considered 40 EOSs describe thermal effects consistently and provide the dependence of thermodynamical quantities on the temperature and the electron fraction. Instead, the majority of models considers matter at zero temperature and in equilibrium with respect to weak interactions (i.e., for beta-equilibrium for neutrinoless conditions). Because temperature effects become important during the merging of the binary components and during the subsequent evolution, we employ an approximate treatment of thermal effects for those EOSs which are given as barotropic relations. This procedure supplements the pressure by an additional ideal-gas component to mimic thermal pressure support, and it requires to choose a corresponding ideal-gas index Γ_th. Appropriate values for Γ_th are in the range of 1.5–2 for high-density matter (Bauswein et al. 2010a). The uncertainties connected to the use of this approximate temperature description and the choice of the ideal-gas index were examined in Bauswein et al. (2010a), where also details about the exact implementation can be found.

From population synthesis studies (Belczynski et al. 2008) and in agreement with pulsar observations (Thorsett & Chakrabarty 1999; Zhang et al. 2011) it is expected that binaries with two NSs with gravitational masses of about M₁ ≈ M₂ ≈ 1.35 M_☉ are the most abundant systems. For this reason we focus in our EOS survey on such equal-mass binaries, albeit we also explore the influence of a system asymmetry by considering 1.2–1.5 M_☉ binaries. For a selected subset of EOS models the full range of possible binary parameters is investigated, varying the single component masses from 1.2 M_☉ to approximately the maximum mass of NSs.

Because of energy and angular momentum losses by gravitational radiation the orbits of NS binaries shrink and the binary components merge after an inspiral period, which lasts roughly 100–1000 Myr for the known systems (Lorimer 2008). The typical outcome of the coalescence of a 1.35–1.35 M_☉ binary system is the formation of a differentially rotating object, potentially a hypermassive NS (i.e., an NS that is more massive than the maximum-mass rigid-rotation configuration and that is stabilized temporarily by differential rotation (Baumgarte et al. 2000)). The merger remnant is surrounded by an extended halo structure of low-density material. Only four EOSs of our sample lead to the prompt formation of a black hole within about one millisecond after the collision because the remnant cannot be supported against the gravitational collapse. For a description of the general dynamics and a more thorough discussion of the collapse behavior we refer to Oechslin et al. (2007), Bauswein et al. (2010a), and Bauswein et al. (2012). In this paper only initially nonrotating NSs are investigated because viscosity is too low to yield tidally locked systems. The stars in NS binaries are therefore expected to rotate slowly in comparison to the orbital angular velocity, justifying the use of an irrotational velocity profile (Bildsten & Cutler 1992; Kochanek 1992).

In this study we analyze the material which becomes gravitationally unbound during or right after merging. In order to estimate whether a given fluid element, i.e., an SPH particle, can escape to infinity, we consider

$\begin{equation} \epsilon _{\mathrm{stationary}}=v^i\tilde{u}_i+\frac{\epsilon }{u^0}+\frac{1}{u^0}-1>0 \end{equation} \tag{ 1 }$

with the coordinate velocity vⁱ, the conserved momentum $\tilde{u}_i$ , and the time-component of the eigen-velocity u⁰ (in geometrical units). This expression can be derived from the hydrodynamical equations by neglecting pressure forces and assuming a stationary metric (Oechslin et al. 2002). The quantity epsilon _stationary is conserved ((d epsilon _stationary/dt) = 0) and at infinity it reduces to the Newtonian expression for the total energy of a fluid element. Hence, a particle with epsilon _stationary > 0 will be unbound. Equation (1) is evaluated in a time-dependent way, and SPH particles that fulfill this criterion 10 ms after merging are considered as ultimately gravitationally unbound. Note that our simulations neglect a possible (smaller) contribution to the ejecta by neutrino-driven winds or magnetically driven outflows from the secular evolution of the merger remnant (Dessart et al. 2009; Wanajo & Janka 2012).

3. EJECTA MASSES

In the following we employ Equation (1) to determine the unbound material in different merger simulations. After a first steep rise of the ejecta mass shortly after the merging of the two NSs, the mass fulfilling the ejecta criterion remains approximately constant (Figure 1). A few models, however, show a continuing, slow increase of the ejecta mass also at later times. The ejecta masses discussed below are computed 10 ms after merging.

**Figure 1.** Evolution of the minimum lapse function α (dashed line) and the amount of unbound matter (solid line) for the symmetric 1.35–1.35 M_☉ merger with the soft SFHO EOS (top panel), the intermediate DD2 EOS (middle panel), and the stiff NL3 EOS (bottom panel).
Download figure:
Standard image High-resolution image

In order to determine the influence of the high-density EOS on the ejecta, we employ approximately the same numerical resolution of about 350,000 SPH particles for all simulations. By using nonuniform SPH particle masses to model the stellar profile (more massive particles in the high-density core and lighter particles in the outer low-density layers) it is possible to achieve a better resolution of low-density regions. This results in an effective mass resolution of about 2 × 10⁻⁶ M_☉, which is comparable to an SPH simulation of about one million equal-mass particles. The influence of the numerical resolution is investigated by performing additional simulations with higher SPH particle numbers. For the TM1 EOS the ejecta masses are found to be (in solar masses) 1.67 × 10⁻³, 1.80 × 10⁻³, 1.71 × 10⁻³, 2.43 × 10⁻³, and 2.07 × 10⁻³ for calculations with about 339 × 10³, 550 × 10³, 782 × 10³, 1007 × 10³, and 1272 × 10³ SPH particles. In these simulations 521, 838, 1241, 2275, and 3080 particles are ejected. Determining the unbound matter 6 ms after merging for the APR EOS, we find ejecta masses (in solar masses) of 5.93 × 10⁻³, 6.08 × 10⁻³, 6.54 × 10⁻³, 6.69 × 10⁻³, and 6.14 × 10⁻³ in simulations with 339 × 10³, 592 × 10³, 782 × 10³, 1007 × 10³, and 1272 × 10³ SPH particles (with 2488, 4473, 6461, 9168, and 16958 particles unbound). Thus, the numerical resolution has an effect on the level of some 10%. This, however, is smaller than the impact of the EOS (see below), which is the focus of this paper. The nonmonotonic variations of the ejecta mass with increasing resolution indicate that statistical fluctuations have some influence on the ejected particle population as well. Note that a small fraction of weakly bound matter could become unbound by the heat generated in the nucleosynthesis processes (Metzger et al. 2010a), which is not taken into account in our hydrodynamical simulations.

3.1. Origin of the Ejecta and Comparison with Other Calculations

As can be seen in Figure 2, most of the ejecta originate from the contact interface between the colliding binary components, which get deformed into drop-like shapes prior to the merging. For the 1.35–1.35 M_☉ binary the ejecta in the shear interface between the stars are separated into two components, each being fed (nearly) symmetrically by material from both colliding stars (top right panel and middle left panel). The matter in the cusps of the stars essentially keeps its direction of motion toward the companion, whereas the backward part of the contact interface mixes with some of the companion matter (top right panel). Both lumps of ejecta are squeezed out from the contact interface and expand on the retral side of the respective companion star, partially slipping over it (middle panels). The bulk matter of the binary components forms a rotating double-core structure where the two dense cores oscillate against each other (not visible because of the logarithmic density scale; see, e.g., the descriptions in Bauswein et al. 2010b and Stergioulas et al. 2011). A first bunch of the matter which was squeezed out from the contact sheet gets unbound in a first expansion phase of the rotating double core structure, which pushes the ejecta outward. This can be seen in Figure 2 (middle right panel and bottom left panel) and also in the evolution of the minimum lapse function α (Figure 1), which is a measure for the compactness of the central object. As the cores separate from each other and the lapse function grows out of its minimum, the ejecta mass increases. A second expansion of the double cores unbinds a smaller amount of matter. Finally, about two milliseconds after the first contact, the triaxial deformation grows into two spiral-arm-like extensions reaching out from the central remnant. These expand into the surrounding, low-density halo fed from the contact interface, push it away, and unbind additional matter (bottom right panel of Figure 2 and second mass-loss episode visible in Figure 1). For different EOSs the different dynamical mechanisms contribute to the ejecta production with different relative strengths. For soft EOSs (top panel of Figure 1) the first steep rise of the ejecta mass due to the expanding double core is much more pronounced, whereas for stiff EOSs (bottom panel of Figure 1) the first increase of the ejecta is very moderate and the late spiral arms unbind most of the ejecta. This can be seen in Figure 1 for the SFHO EOS representing a soft EOS, for the NL3 EOS as a stiff example, and the intermediate case of the DD2 EOS.

A significantly smaller fraction of the ejecta (typically below 25%) stems from the outer faces of the merging stars opposite to the contact layer (SPH particles at the outer left and outer right ends of the stellar body in the top right panel of Figure 2). During merging this matter at the rear of the star (SPH particles with nearly horizontal velocity vectors at the top and bottom merger tails in the middle left panel) lags behind the rotation of the star's center and is hit and ablated by the "nose" of the companion shortly after the snapshot shown in the middle left panel (see velocity arrows of particles with opposite color at the tips of the noses). This material gets mostly unbound in the first expansion phase of the oscillating double-core structure. Such type of ejecta is less abundant for stiff EOSs.

Different from relativistic simulations, Newtonian models find the ejecta originating mostly from the tips of tidal tails (see, e.g., Korobkin et al. 2012), in particular also in the case of symmetric binaries. Relativistic calculations (within the conformal flatness framework; Oechslin et al. 2007; Goriely et al. 2011) yield the dominant ejection from the contact interface as described above (see also the inset of Figure 1 in Goriely et al. 2011). Recently, the fully relativistic simulations of Hotokezaka et al. (2013) have provided further support for the ejecta origin from the contact interface, confirming the conclusions of Oechslin et al. (2007) and Goriely et al. (2011). This points to qualitative differences between the Newtonian and relativistic mass-loss dynamics with the important difference that in relativistic simulations all of the ejecta are shock-heated while in Newtonian calculations the cold, tidally stripped material dominates.

A quantitative comparison between Newtonian (Korobkin et al. 2012; Rosswog et al. 2013; Piran et al. 2013; Rosswog 2013) and relativistic simulations (Oechslin et al. 2007; Goriely et al. 2011) also reveals considerable discrepancies. For instance, simulations of a 1.4–1.4 M_☉ merger with the Shen EOS in Newtonian theory produce more than 10⁻² M_☉ ejecta (Korobkin et al. 2012; Rosswog et al. 2013; Piran et al. 2013; Rosswog 2013), whereas the relativistic calculations of this study and in Oechslin et al. (2007) and Goriely et al. (2011) yield only a few 10⁻³ M_☉ of unbound material for the 1.35–1.35 M_☉ binary with the same EOS. Comparing the results of our study with the likewise relativistic calculations in Hotokezaka et al. (2013) shows very good agreement for all four EOSs used in Hotokezaka et al. (2013). For example, for the APR EOS with Γ_th = 2 both groups find about 5 × 10⁻³ M_☉ of unbound matter. This is remarkable because the implementations differ with respect to the hydrodynamics scheme, which is an SPH algorithm here but a grid-based, high-resolution central scheme in Hotokezaka et al. (2013). (Note that we employ the conformal flatness approximation whereas the calculations in Hotokezaka et al. 2013 are conducted within full general relativity.) These findings provide confidence in the results on the quantitative level and point toward fundamental differences between Newtonian and relativistic treatments. Such differences are not unexpected because NSs are more compact in general relativity than in Newtonian gravity. The stronger gravitational attraction prevents the formation of pronounced tidal tails at the outer faces of the colliding stars and increases the strength of the collision.

3.2. Equation of State Dependence

Several NS EOSs have been employed in merger simulations by different groups, but a large, systematic investigation of the EOS dependence of the ejecta production is still missing in particular with a consistent description of thermal effects. For a given EOS the radius R_1.35 of a nonrotating NS with 1.35 M_☉ is a characteristic quantity specifying the compactness of NSs. Therefore, we use R_1.35 to describe the influence of the high-density EOS on the amount of NS merger ejecta.

The upper left panel of Figure 3 displays the amount of unbound material as a function of R_1.35 for all 40 EOSs used in our study (see also Table 1). Red crosses identify EOSs which provide the full temperature dependence. The black symbols correspond to barotropic zero-temperature EOSs, which are supplemented by a thermal ideal-gas component choosing Γ_th = 2 (see Section 2). Results based on the same zero-temperature EOS but with Γ_th = 1.5 are given in blue at the same radius R_1.35. Small symbols indicate results for EOSs which are excluded by the pulsar mass measurement of Antoniadis et al. (2013). Circles mark cases which lead to the prompt collapse to a black hole.

Table 1. Model Properties

Model M₁–M₂	Γ_th	R_1.35	M_ej	v	E_kin
Model M₁–M₂	Γ_th	(km)	(10⁻³ M_☉)	(c)	(10⁵⁰ erg)
NL3 1.35–1.35	Full	14.75	2.09	0.18	0.98
NL3 1.35–1.35	2	14.75	1.57	0.34	2.03
NL3 1.35–1.35	1.8	14.75	1.60	0.32	2.99
NL3 1.35–1.35	1.5	14.75	1.86	0.30	1.98
GS1 1.35–1.35	Full	14.72	2.19	0.32	3.43
Shen 1.35–1.35	Full	14.56	2.33	0.23	4.80
TM1 1.35–1.35	Full	14.49	1.67	0.16	0.74
TM1 1.35–1.35	2	14.49	1.37	0.36	2.02
TM1 1.35–1.35	1.8	14.49	1.33	0.34	1.77
TM1 1.35–1.35	1.5	14.49	1.53	0.32	1.86
TMA 1.35–1.35	Full	13.86	2.05	0.18	1.19
LS375 1.35–1.35	Full	13.56	2.58	0.30	3.39
GS2 1.35–1.35	Full	13.38	2.74	0.19	2.16
DD2 1.35–1.35	Full	13.21	3.07	0.22	2.18
DD2 1.35–1.35	2	13.21	2.57	0.34	3.31
DD2 1.35–1.35	1.8	13.21	2.26	0.32	2.61
DD2 1.35–1.35	1.5	13.21	2.72	0.30	2.90
LS220 1.35–1.35	Full	12.64	1.99	0.28	2.08
LS180 1.35–1.35 x	Full	12.14	2.26	0.29	3.02
SFHX 1.35–1.35	Full	11.98	6.16	0.22	4.36
SFHO 1.35–1.35	Full	11.92	4.83	0.23	3.61
SFHO 1.35–1.35	2	11.92	2.96	0.32	3.37
SFHO 1.35–1.35	1.8	11.92	3.26	0.34	4.18
SFHO 1.35–1.35	1.5	11.92	3.82	0.30	4.14
eosL 1.35–1.35	2	15.74	1.49	0.22	0.77
eosL 1.35–1.35	1.5	15.74	3.45	0.20	1.49
MS1 1.35–1.35	2	14.99	1.17	0.27	0.98
MS1 1.35–1.35	1.5	14.99	2.38	0.21	1.19
MS1b 1.35–1.35	2	14.59	1.67	0.25	1.26
MS1b 1.35–1.35	1.5	14.59	3.64	0.21	1.85
Glenh3 1.35–1.35 x	2	14.52	1.08	0.23	0.62
Glenh3 1.35–1.35 x	1.5	14.52	1.69	0.22	0.90
MS2 1.35–1.35 x	2	14.25	0.81	0.26	0.65
H3 1.35–1.35 x	2	13.95	1.43	0.27	1.15
H4 1.35–1.35	2	13.95	1.28	0.27	1.09
H4 1.35–1.35	1.5	13.95	1.93	0.27	1.64
Heb6 1.35–1.35	2	13.33	1.55	0.45	3.86
Heb6 1.35–1.35	1.5	13.33	3.43	0.24	2.58
eosO 1.35–1.35	2	12.85	3.52	0.28	3.04
eosO 1.35–1.35	1.5	12.85	4.62	0.25	3.27
ALF2 1.35–1.35	2	12.78	3.80	0.28	3.36
ALF2 1.35–1.35	1.5	12.78	4.49	0.27	3.80
BSk21 1.35–1.35	2	12.54	3.36	0.32	3.89
BSk21 1.35–1.35	1.5	12.54	4.37	0.27	3.81
Heb4 1.35–1.35	2	12.51	1.89	0.43	4.33
Heb4 1.35–1.35	1.5	12.51	2.41	0.39	5.13
MPA1 1.35–1.35	2	12.49	3.64	0.30	3.60
MPA1 1.35–1.35	1.5	12.49	4.48	0.29	4.35
Heb5 1.35–1.35	2	12.38	2.63	0.43	5.93
Heb5 1.35–1.35	1.5	12.38	2.90	0.38	5.89
eosC 1.35–1.35 x	1.5	12.06	3.09	0.27	2.49
ENG 1.35–1.35	2	12.05	5.29	0.29	5.01
ENG 1.35–1.35	1.5	12.05	6.32	0.26	5.30
APR3 1.35–1.35	2	12.04	4.65	0.30	4.69
APR3 1.35–1.35	1.5	12.04	6.15	0.27	5.50
Heb3 1.35–1.35	2	12.03	2.99	0.41	6.43
Heb3 1.35–1.35	1.5	12.03	3.70	0.37	6.93
BurgioNN 1.35–1.35	2	11.99	2.47	0.39	4.73
BurgioNN 1.35–1.35	1.5	11.99	2.44	0.37	4.67
SLy4 1.35–1.35	2	11.79	3.99	0.29	3.75
SLy4 1.35–1.35	1.5	11.79	6.40	0.27	5.53
BSk20 1.35–1.35	2	11.74	4.68	0.31	4.90
BSk20 1.35–1.35	1.5	11.74	7.83	0.26	6.80
ALF4 1.35–1.35 x	2	11.60	5.70	0.30	6.07
ALF4 1.35–1.35 x	1.5	11.60	7.40	0.29	7.65
Heb2 1.35–1.35	2	11.42	4.95	0.34	7.90
Heb2 1.35–1.35	1.5	11.42	5.01	0.37	9.07
APR 1.35–1.35	2	11.33	5.96	0.31	6.37
APR 1.35–1.35	1.5	11.33	7.38	0.30	7.90
BB2 1.35–1.35 x	2	11.30	4.95	0.28	4.17
eosUU 1.35–1.35	2	11.18	7.02	0.32	7.80
eosUU 1.35–1.35	1.5	11.18	9.42	0.31	10.0
Heb1 1.35–1.35	2	10.81	6.85	0.31	8.54
Heb1 1.35–1.35	1.5	10.81	10.88	0.32	12.54
eosAU 1.35–1.35	2	10.44	4.05	0.36	5.64
eosAU 1.35–1.35	1.5	10.44	1.51	0.29	1.69
NL3 1.2–1.5	Full	14.75	7.95	0.19	4.50
GS1 1.2–1.5	Full	14.72	6.43	0.22	6.55
Shen 1.2–1.5	Full	14.56	5.66	0.30	8.33
TM1 1.2–1.5	Full	14.49	8.66	0.17	3.94
TMA 1.2–1.5	Full	13.86	10.21	0.20	6.40
LS375 1.2–1.5	Full	13.56	6.66	0.27	7.60
GS2 1.2–1.5	Full	13.38	10.69	0.18	6.14
DD2 1.2–1.5	Full	13.21	8.79	0.20	4.97
LS220 1.2–1.5	Full	12.64	13.22	0.18	7.34
LS180 1.2–1.5 x	Full	12.14	18.58	0.20	12.13
SFHX 1.2–1.5	Full	11.98	14.67	0.19	7.91
SFHO 1.2–1.5	Full	11.92	13.39	0.22	8.94
NL3 1.2–1.2	Full	14.75	2.15	0.17	0.91
NL3 1.2–1.35	Full	14.75	4.25	0.21	2.74
NL3 1.2–1.6	Full	14.75	9.96	0.19	5.57
NL3 1.2–1.8	Full	14.75	15.68	0.15	5.75
NL3 1.35–1.5	Full	14.75	2.72	0.24	2.25
NL3 1.35–1.8	Full	14.75	18.81	0.21	11.31
NL3 1.5–1.5	Full	14.75	1.70	0.20	1.04
NL3 1.5–1.8	Full	14.75	8.10	0.21	4.94
NL3 1.6–1.6	Full	14.75	3.74	0.22	2.59
NL3 1.8–1.8	Full	14.75	9.08	0.24	7.25
NL3 1.35–2.0	Full	14.75	12.85	0.20	7.62
NL3 2.0–2.0	Full	14.75	1.91	0.29	2.18
DD2 1.2–1.2	Full	13.21	3.09	0.17	1.37
DD2 1.2–1.35	Full	13.21	3.17	0.20	2.06
DD2 1.2–1.6	Full	13.21	10.90	0.20	6.39
DD2 1.2–1.8	Full	13.21	17.08	0.17	6.72
DD2 1.35–1.5	Full	13.21	3.57	0.25	3.13
DD2 1.35–1.8	Full	13.21	14.85	0.21	9.48
DD2 1.5–1.5	Full	13.21	5.38	0.26	4.66
DD2 1.5–1.8	Full	13.21	18.84	0.25	15.52
DD2 1.6–1.6	Full	13.21	7.80	0.27	7.40
DD2 1.8–1.8	Full	13.21	1.37	0.26	1.63
DD2 1.35–2.0	Full	13.21	6.41	0.31	9.64
DD2 2.0–2.0	Full	13.21	0.25	0.25	0.25
SFHO 1.2–1.2	Full	11.92	1.88	0.21	1.26
SFHO 1.2–1.35	Full	11.92	5.44	0.22	3.86
SFHO 1.2–1.6	Full	11.92	16.91	0.21	11.10
SFHO 1.2–1.8	Full	11.92	5.78	0.34	10.08
SFHO 1.35–1.5	Full	11.92	18.73	0.23	13.34
SFHO 1.35–1.8	Full	11.92	11.76	0.31	16.22
SFHO 1.5–1.5	Full	11.92	4.10	0.27	4.13
SFHO 1.5–1.8	Full	11.92	6.34	0.42	14.40
SFHO 1.6–1.6	Full	11.92	1.13	0.21	1.00
SFHO 1.8–1.8	Full	11.92	0.17	0.29	0.24

Notes. Basic properties of simulations and the employed EOSs. In the second column the thermal ideal-gas index is given if thermal effects are incorporated approximately; "full" refers to a simulation with a fully consistent treatment of thermal effects. R_1.35 specifies the radius of a nonrotating NS with 1.35 M_☉. M_ej is the amount of unbound matter. The fifth column provides the average outflow velocity. E_kin denotes the kinetic energy of the ejecta. A cross in the first column indicates EOSs which are incompatible with the observation of a 2.01 ± 0.04 M_☉ NS (Antoniadis et al. 2013).

Download table as: ASCIITypeset images: 1 2

One can recognize a clear EOS dependence of the ejecta mass, where EOSs with a high compactness of the NSs lead to an enhanced production of unbound material. The ejecta mass can be as big as about 0.01 M_☉ for symmetric mergers with a total binary mass M_tot = M₁ + M₂ = 2.7 M_☉. EOSs with relatively large NS radii lead to outflow masses of about 0.001–0.002 M_☉. For EOSs with approximately the same R_1.35 the ejecta masses show a scatter of up to 0.003 M_☉. However, considering only EOSs with a fully consistent description of thermal effects (red symbols), the variations are smaller. Only one simplified EOS (eosAU) leads to a prompt collapse of the merger remnant and yields significantly smaller ejecta masses (circles). Using the radius R_1.6 of a nonrotating NS with 1.6 M_☉ or the radius R_max of the maximum-mass Tolman–Oppenheimer–Volkoff solution to characterize an EOS results in diagrams similar to the upper left panel of Figure 3. However, no clear trend can be found for the ejecta mass as a function of the maximum mass of nonrotating NSs. We therefore conclude that the NS compactness is the crucial EOS parameter determining the ejecta mass. Indications of such a behavior were already observed in simulations for four simplified EOSs with an approximate temperature treatment (Hotokezaka et al. 2013).

The dynamics of the merger explain why small NS radii lead to higher ejecta masses. For smaller R_1.35 the inspiral phase lasts longer and the stars reach higher velocities before they collide. The relation between the impact velocity and the NS radius is clearly seen in the lower left panel of Figure 3, which displays the sum of the maxima of the coordinate velocities of the mass centers of the two binary components. The maximum of the coordinate velocity is reached shortly after the first contact, before the cores of the NSs are decelerated by the collision. The clash of more compact NSs is more violent, and more material is squeezed out from the collision interface, for which reason the negative correlation of the velocities with R_1.35 is reflected by a similar negative correlation of M_ejecta and R_1.35. Moreover, for smaller R_1.35 the central remnant consisting of the double cores rotates faster and the bounce and rebounce are stronger, i.e., the surface of the remnant moves faster and pushes away matter more efficiently.

3.3. Influence of the Approximate Treatment of Thermal Effects

A number of simulations of our survey (black and blue symbols) as well as calculations by other groups (e.g., Roberts et al. 2011; Hotokezaka et al. 2013) rely on an approximate description of thermal effects in the EOS, which requires the specification of an effective thermal ideal-gas index (see Section 2). As can be seen in the upper left panel of Figure 3, the choice of the value for this ideal-gas index has a considerable impact on the ejecta mass; the simulations with Γ_th = 1.5 (blue symbols) yield generally more unbound matter. The reason is the reduced thermal pressure support, which means that the two dense cores can approach each other more closely during the collision, which results in a more violent impact and shearing motion and thus in more material being squeezed out from the collision interface and in a more powerful oscillation of the central remnant. This can be clearly seen by following the centers of mass of the two cores or the evolution of the central lapse function.

The optimal choice of Γ_th is a priori unclear and may be different for different EOSs. To address this issue, we performed additional simulations (not shown in Figure 3, but listed in Table 1) for temperature-dependent EOSs (SFHO, DD2, TM1, NL3) after reducing them to the zero-temperature sector (with the constraint of neutrinoless beta-equilibrium) and supplementing them with the approximate description of thermal effects using Γ_th = 1.5, 1.8, and 2. The comparison with the fully consistent simulations reveals that generally a choice of Γ_th = 1.5 yields the best quantitative agreement with only a slight underestimation of about 10% (20% for the SFHO EOS), whereas the ejecta masses with Γ_th = 1.8 or Γ_th = 2 are significantly too low compared to the fully consistent models (for the tested EOSs between 15% and 40% for Γ_th = 2 and between 20% and 35% for Γ_th = 1.8).

The fact that a relatively low Γ_th reproduces the ejecta properties best contrasts the finding that a higher Γ_th (in the range between 1.5 and 2) has turned out to be more suitable for describing gravitational-wave features and the post-merger collapse behavior (see Bauswein et al. 2010a), i.e., the bulk mass motion of the colliding stars. The reason for this discrepancy is the density dependence of Γ_th, which drops from about 2 at supranuclear densities to about 4/3 for densities below ≈10¹¹ g cm⁻³ (see Figure 2 in Bauswein et al. 2010a). While the dynamics of the bulk mass of the merging objects, which is responsible for the gravitational-wave production, is fairly well captured with a choice of Γ_th ∼ 2, unbound fluid elements originating from the inner crust, where most of the ejecta stem from, encounter different density regimes. Consequently, the ejecta behavior cannot be well modeled with the Γ_th that is appropriate for high-density matter. We found the best compromise to be Γ_th ≈ 1.5, but we stress that the results of simulations using an approximate treatment of thermal effects should be taken with caution and do not need to be quantitatively reliable in all aspects.

3.4. Asymmetric Binaries

Even though most NS binaries are expected to be nearly symmetric systems with a total mass of about 2.7 M_☉, we investigate the mass ejection of an asymmetric setup to check whether also in this case an EOS dependence exists. The upper right panel of Figure 3 displays the ejecta masses for simulations of asymmetric binaries with a 1.2 M_☉ NS and a 1.5 M_☉ NS (see also Table 1). Again, the radius R_1.35 of a nonrotating NS is used to characterize different EOSs. Here we restrict ourselves to EOSs which provide the full temperature dependence. In comparison to the symmetric binary mergers the amount of unbound material is significantly larger. The ejecta masses are about a factor of two higher than for the symmetric binaries with the same total binary mass. Also for asymmetric binaries a decrease of M_ej with bigger R_1.35 is visible, but the scatter between models with similar R_1.35 is larger. The lower right panel of Figure 3 shows the sum of the maxima of the coordinate velocities of the mass centers of the two asymmetric binary components. As in the symmetric case the two stars collide with a higher impact velocity if the initial radii of the NSs are smaller.

Due to the asymmetry the dynamics of the merger proceeds differently from the symmetric case (see Figure 4). Prior to the merging the less massive binary component is deformed to a drop-like structure with the cusp pointing to the 1.5 M_☉ NS (top panels). After the stars begin to touch each other, the lighter companion is stretched and a massive tidal tail forms (middle left panel). The deformed 1.2 M_☉ component is wound around the more massive companion (middle panels). Also in the case of asymmetric mergers the majority of the ejecta originate from the contact interface of the collision, i.e., from the cusp of the "tear drop" and from the equatorial surface of the more massive companion, where the impact ablates matter (see top panels). Some matter at the tip of the cusp directly fulfills the ejecta criterion (top right panel), while the majority obtains an additional push by the interaction with the asymmetric, mass-shedding central remnant and the developing spiral arms (middle right and bottom panels). A smaller amount of ejecta of roughly 25% originates from the outer end of the primary tidal tail (particles in the lower part of the top right panel). A part of this matter becomes unbound by tidal forces (at the tip of the tidal tail in the middle left panel) and the other fraction by an interaction with the central remnant (middle left panel).

Figure 5 displays the distribution of the ejecta in a plane perpendicular to the binary orbit for the symmetric merger (left panel) compared to the asymmetric merger (right panel) for the last timesteps shown in Figures 2 and 4, respectively. A considerable fraction of the ejected matter is expelled with large direction angles relative to the orbital plane. For a timestep about 5 ms later the ejecta geometry is visualized (azimuthally averaged) in Figure 6 excluding the bound matter. For both mergers the outflows exhibit a (torus or donut-like) anisotropy with an axis ratio of about 2:3. The velocity fields also show a slight dependence on the direction.

**Figure 6.** Ejecta geometry visualized by the rest-mass density (color-coded and logarithmically plotted in g cm⁻³) excluding matter of the bound central remnant for the 1.35–1.35 M_☉ merger (left) and the 1.2–1.5 M_☉ merger (right) with the DD2 EOS. Density contours are obtained by azimuthal averaging. Arrows represent the coordinate velocity field where an arrow length of 200 km corresponds to the speed of light. The time of the snapshots is given below the color bar of each panel. The visualization tool SPLASH was used to convert SPH data to grid data (Price 2007).
Download figure:
Standard image High-resolution image

3.5. Binary Parameter Dependence

The exploration of the full space of possible binary parameters is interesting for the determination of the highest and lowest possible ejecta mass for a given EOS and to understand the influence of the binary setup on the ejecta production. Such an investigation has been conducted only for one EOS (Shen) by Newtonian calculations (Korobkin et al. 2012; Rosswog et al. 2013; Rosswog 2013) and within a relativistic framework (Oechslin et al. 2007), which revealed quantitative differences between both approaches. Other surveys using different EOSs have been restricted to a limited variation of the binary masses (Roberts et al. 2011; Goriely et al. 2011; Hotokezaka et al. 2013). Here we present the dependence of the ejecta mass on the mass ratio q = M₁/M₂ and the total binary mass M_tot = M₁ + M₂ for a subset of EOSs employed in our study. The NL3, DD2, and SFHO EOSs are chosen because they are representative for the full set of possible EOSs: while the NL3 EOS is relatively stiff, resulting in R_1.35 = 14.75 km, the soft SFHO EOS produces rather compact NSs with R_1.35 = 11.74 km, and the DD2 represents an intermediate case with R_1.35 = 13.21 km. Figure 7 displays the amount of unbound matter as a function of the mass ratio q and total binary mass M_tot for these three EOSs (see also Table 1). The simulated binary configurations that form a differentially rotating NS are indicated in the figures by crosses, whereas systems leading to prompt black hole formation (within about one millisecond after the first contact) are marked with circles (the ejecta properties are extracted 10 ms after the merging).

**Figure 7.** Ejecta mass in M_☉ as function of the mass ratio q and the total binary mass M_tot for the soft SFHO EOS (upper panel), the intermediate DD2 EOS (center panel), and the stiff NL3 EOS (bottom panel). In all panels the simulated binary setups are marked by symbols. Crosses indicate the formation of a differentially rotating NS remnant, while circles identify configurations which lead to a direct gravitational collapse.
Download figure:
Standard image High-resolution image

All three EOSs show qualitatively the same behavior. A clear trend of increasing ejecta masses with larger binary asymmetry is visible. For symmetric binaries that do not undergo gravitational collapse there is also a slight tendency of higher total binary masses leading to more unbound material. The increase with M_tot is more pronounced for asymmetric systems. The occurrence of a prompt collapse results in a significant drop of the ejecta mass. This is an important qualitative difference to Newtonian calculations, which cannot determine and follow the relativistic gravitational collapse. The threshold for the prompt collapse depends sensitively on the EOS, and soft EOSs lead to a collapse for relatively small M_tot.

For all EOSs the maximum ejecta masses (in all cases slightly below 0.02 M_☉) are 4–10 times higher than the amount of unbound matter of the symmetric 1.35–1.35 M_☉ binaries. Here, the absolute differences between the maximum and the minimum ejected mass for non-collapsing cases are larger for stiff EOSs like the NL3 and less pronounced for soft EOSs like the SFHO. Soft EOSs yield steeper gradients of the ejecta mass in the binary parameter space, i.e., a certain variation in the binary parameters leads to a larger change of the ejecta masses than is the case for a stiff EOS. To a good approximation and ignoring cases with a prompt black hole formation, the setup with two stars of about 1.35 M_☉ is the system that produces the smallest amount of ejecta for the majority of investigated EOSs. The SFHO EOS is one of the exceptions, for which, e.g., the 1.2–1.2 M_☉ binary yields a factor two to three less ejecta than the 1.35–1.35 M_☉ setup.

Based on Newtonian calculations a fit formula for the ejecta mass as a fraction of M_tot was proposed as a function of η = 1 − 4M₁M₂/(M₁ + M₂)² in Korobkin et al. (2012) and Rosswog (2013). Reviewing our data (even without the prompt collapse cases), we find a more complicated behavior and can neither confirm the validity of the suggested fit formula nor find a generalization of it. This is not unexpected in view of the quantitative and qualitative differences between Newtonian and relativistic simulations discussed above.

3.6. Folding with Binary Populations

The dependence of the ejecta mass on the binary parameters is essential to determine the total amount of ejecta produced by the binary population within a certain time and thus to estimate the average amount of ejecta per merger event. The properties of the NS binary population are provided by theoretical binary evolution models, which still contain considerable uncertainties in many complexities of single star evolution and binary interaction. Using the standard model of Dominik et al. (2012), the folding of our results with the binary population yields an average ejecta mass per merger event of about 3.6 × 10⁻³ M_☉ for the NL3 EOS, 3.2 × 10⁻³ M_☉ for the DD2 EOS, and 4.3 × 10⁻³ M_☉ for the SFHO EOS. Therefore, the ejecta masses of the 1.35–1.35 M_☉ binary mergers give numbers for the three cases which approximate the average amount of ejecta per merger event quite well (within 70% for NL3, 3% for DD2, 11% for SFHO). This finding is simply a consequence of the fact that the binary distribution is strongly peaked around nearly symmetric systems with M_tot ≈ 2.5 M_☉ so that the average ejecta mass is not sensitive to the larger ejecta production of asymmetric systems in the suppressed wings of the binary distribution.

4. NUCLEOSYNTHESIS

4.1. r-process Abundances

The potential of NS mergers to produce heavy r-process elements in their ejecta has been manifested by several studies based on hydrodynamical simulations (Freiburghaus et al. 1999; Metzger et al. 2010b; Roberts et al. 2011; Goriely et al. 2011; Korobkin et al. 2012). These investigations have considered only a few high-density EOSs (two EOSs were used in Goriely et al. 2011). Since the NS EOS affects sensitively the dynamics of NS mergers and thus the properties of the ejecta (amount, expansion velocity, electron fraction, temperature), we explore here the influence of the NS EOS on the r-process nucleosynthesis in a systematic way.

For a selected, representative set of EOSs we extract the thermodynamical histories of fluid elements which get gravitationally unbound. For these trajectories nuclear network calculations were performed as in Goriely et al. (2011), where details on the reaction network, the temperature postprocessing, and the density extrapolation beyond the end of the hydrodynamical simulations can be found. The reaction network includes all 5000 species from protons up to Z = 110 lying between the valley of β-stability and the neutron-drip line. All fusion reactions on light elements, as well as radiative neutron captures, photodisintegrations, β-decays, and fission processes, are included. The corresponding rates are based on experimental data whenever available or on theoretical predictions otherwise, as prescribed in the BRUSLIB nuclear astrophysics library (Xu et al. 2013).

Figure 8 shows the final nuclear abundance patterns for the 1.35–1.35 M_☉ mergers described by the NL3 (blue), DD2 (red), and SFHO (green) EOSs. For every model about 200 trajectories were processed, which roughly correspond to about 1/10 of the total ejecta. Comparing the final abundance distributions of the DD2 EOS for about 200 and the full set of 1000 fluid-element histories reveals a very good quantitative agreement, which proves that a properly chosen sample of about 200 trajectories is sufficient to be representative for the total amount of unbound matter.

The scaled abundance patterns displayed in Figure 8 match closely the solar r-process composition above mass number A ≈ 140. In particular the third r-process peak around A = 195 is robustly reproduced by all models. Above mass number A ≈ 100 the results for the different NS EOSs hardly differ. For all three displayed models the peak around A ≈ 140 is produced by fission recycling, which occurs when the nuclear flow reaches fissioning nuclei around ²⁸⁰No at the end of the neutron irradiation during the β-decay cascade. The exact shape and location of this peak are therefore strongly affected by the theoretical modeling of the fission processes (including in particular the fission fragment distribution of the fissioning nuclei), which are still subject to large uncertainties (Goriely et al. 2009b). Hence, the deviations from the solar abundance pattern between A ≈ 130 and A ≈ 170 are not unexpected, while the third r-process peak around A = 195 is a consequence of the closed neutron shell at N = 126, which is robustly predicted by theoretical models. Very similar results were obtained for NS merger models performed with the LS220 and Shen EOSs in Goriely et al. (2011).

In Figure 9 the normalized abundance patterns are shown for asymmetric 1.2–1.5 M_☉ mergers employing the same representative EOSs as in Figure 8. Again a very good agreement between the solar r-process abundances and the calculated element distributions above A ≈ 130 is found for all three high-density EOSs. This confirms earlier findings that the binary mass ratio has a negligible effect on the abundance yield distribution (Goriely et al. 2011). It also confirms that the ejected abundance distribution is rather insensitive to the adopted EOS.

Besides the three temperature-dependent EOSs considered above, we conducted network calculations also for merger models computed with zero-temperature EOSs supplemented by an approximate treatment of thermal effects with a Γ_th = 2 ideal-gas component (BSk20, BSk21). In this case the temperature is estimated following a procedure described in Etienne et al. (2008), which converts the specific internal energy to temperature values. Doing so, it is assumed that the energy of the thermal ideal-gas component is composed of the thermal energy of an ideal nucleon gas and a contribution from ultrarelativistic particles (photons, possibly electrons, positrons, and neutrinos).

The network calculations for the BSk20 and BSk21 EOSs yield an abundance pattern above A ≈ 130 very similar to the other fully temperature-dependent EOSs (see Figure 10). Differences between the fully consistent models and the simulations with approximate temperature treatment are found below mass number A ≈ 50, where the calculations with the BSk EOSs yield a lower amount of elements with 5 < A < 50 but a higher mass fraction of hydrogen, deuterium, and helium. The reason is the higher temperatures found with the BSk EOSs at the beginning of the network calculations which lead to a reduced recombination of nucleons and α-particles and consequently a smaller production of heavier nuclei. In this respect the conditions in the outflows of these models resemble the situation in the neutrino-driven winds of core-collapse supernovae but for significantly higher neutron excesses.

**Figure 10.** Nuclear abundance pattern for the 1.35–1.35 M_☉ mergers with the BSk20 (red) and BSk21 (blue) EOSs compared to the solar r-process abundance distribution (black).
Download figure:
Standard image High-resolution image

Overall, it is reassuring that the BSk models yield a similar abundance pattern of r-process elements, although these calculations rely on an approximate incorporation of thermal effects and a rough estimate of temperatures. It is an important finding of our work that r-process elements are robustly produced for a representative, diverse sample of high-density EOSs and that the outcome is insensitive to the exact initial temperature conditions and the binary setup.

4.2. Merger Rates

The above-mentioned variations in the production of light elements are also reflected in the fraction of the ejecta which end up as r-process elements. In Table 2 a clear difference is observed for temperature-dependent EOSs (NL3, DD2, SFHO) and the BSk EOSs with approximate temperature treatment. While in the former cases about 96%–99% of the ejecta are converted to r-process elements, only 93%–95% of the ejecta are processed to r-process material in the latter models. This variation is a consequence of the temperature history affecting the production of light nuclei, as discussed above.

Table 2. Nucleosynthesis Calculations

Model M₁–M₂	M_ej	t_peak	M_r-process/M_ej	²³²Th/²³⁸U	²³²Th/²³⁵U	Th/Eu	E_heat	f
Model M₁–M₂	(10⁻³ M_☉)	(days)	M_r-process/M_ej	²³²Th/²³⁸U	²³²Th/²³⁵U	Th/Eu	(MeV A⁻¹)	f
NL3 1.35–1.35	2.09	0.171	0.989	1.644	1.072	0.695	3.34	1.8 × 10⁻⁶
DD2 1.35–1.35	3.07	0.189	0.980	1.671	1.080	0.627	3.13	1.6 × 10⁻⁶
SFHO 1.35–1.35	4.83	0.228	0.991	1.642	1.039	0.579	3.17	1.4 × 10⁻⁶
NL3 1.2–1.5	7.95	0.338	0.964	1.670	1.112	0.714	3.36	1.3 × 10⁻⁶
DD2 1.2–1.5	8.79	0.354	0.986	1.697	1.109	0.658	3.11	1.2 × 10⁻⁶
SFHO 1.2–1.5	13.39	0.418	0.974	1.685	1.085	0.543	3.12	1.1 × 10⁻⁶
BSk21 1.35–1.35	3.36	0.162	0.948	1.660	1.023	0.704	2.97	1.6 × 10⁻⁶
BSk20 1.35–1.35	4.68	0.195	0.931	1.689	1.021	0.698	2.96	1.5 × 10⁻⁶

Notes. Selected models for which nucleosynthesis calculations were performed. M_ej is the amount of unbound matter, whereas t_peak is the peak time of an optical transient associated with an NS merger (see Section 5). The fourth column gives the fraction of the ejecta which is processed into r-process elements. The production ratios of certain elements and isotopes are provided in the fifth to seventh columns. E_heat denotes the total amount of energy released by radioactive decays (without neutrino energy). The factor f approximates the radioactive heat generation around the time of the optical peak luminosity relative to the rest-mass energy of the ejecta (see text).

Download table as: ASCII Typeset image

Despite these differences, it is justified to assume that almost the total amount of ejecta is converted to heavy r-process elements, also considering the uncertainties in the determination of the exact ejected masses from simulations. Furthermore, in Section 3 it was argued that the total amount of ejecta, and thus the r-process material produced per event by the population of NS binaries, can be well represented by the yield of the 1.35–1.35 M_☉ merger. This allows an important consistency check by comparing the theoretically expected production to the observed amount of r-process matter with A ≳ 130 in the Galaxy, which is estimated to be about 4 × 10³ M_☉ (Qian 2000). In order to produce this amount of heavy r-process elements within the Galactic history for about 10¹⁰ yr, one requires a merger rate of 4 × 10⁻⁴ yr⁻¹ if every coalescence ejects on average 10⁻³ M_☉, which in our EOS survey corresponds to the lower bound on the ejecta mass of 1.35–1.35 M_☉ mergers (see upper left panel of Figure 3). Similarly, assuming that NS mergers are the dominant source of heavy r-process elements, the upper bound of 10⁻² M_☉ for the ejected mass from 1.35–1.35 M_☉ binaries (i.e., for EOSs with small R_1.35) would be compatible with a merger rate of 4 × 10⁻⁵ yr⁻¹. These rate estimates lie in the ballpark of theoretical predictions ranging from 10⁻⁶ to 10⁻³ yr⁻¹ (Abadie et al. 2010). This implies that all EOSs of our survey are compatible with NS mergers being the dominant or a major source of r-process elements.

This conclusion on the merger rate may be tested against future observations, in particular by multiple gravitational-wave detections or frequent observations of electromagnetic counterparts. Our work emphasizes that, in addition to a more accurate merger rate, information on the high-density EOS is needed to shed light on NS mergers as a major source of r-process elements. More specifically, a merger rate of about 4 × 10⁻⁵ yr⁻¹ may imply either that nearly all heavy r-process elements are made by NS mergers in the case of a soft high-density EOS with small NS radius R_1.35, or that only 1/10 of the observed r-process material originates from NS binaries if a stiff NS EOS with large R_1.35 is confirmed.

Inversely, considering the robustness of the r-process nucleosynthesis in NS mergers, one can infer that, assuming a minimal production of ∼10⁻³ M_☉ of r-nuclei per event and a constant merger rate during the life of the Galaxy, this rate cannot be higher than roughly 4 × 10⁻⁴ yr⁻¹; otherwise, the galactic r-process material would be more than currently observed in the Galaxy. This bound is comparable to the "optimistic" limit given in Abadie et al. (2010). Thus, the r-process element content in the Galaxy establishes further, independent evidence for an upper limit on the merger rate below ∼4 × 10⁻⁴ yr⁻¹. A restriction to soft EOSs, e.g., from other physical constraints like Hebeler et al. (2010), would lead to a smaller upper limit for the event rate of NS mergers.

4.3. Actinide Production Ratios and Stellar Chronometry

Some of the heaviest long-lived radioactive nuclei produced by the r-process can be used as nucleo-cosmochronometers. In particular, the abundance ratios of thorium to europium and thorium to uranium have been proposed for estimating the age of the oldest stars in our Galaxy. More specifically, a simple comparison of the observed abundance ratio with the production ratio can provide an age estimate of the contaminated object (Butcher 1987; Francois et al. 1993; Cowan et al. 1999; Goriely & Clerbaux 1999; Goriely & Arnould 2001; Cayrel et al. 2001; Frebel et al. 2007; Sneden et al. 2008). In addition, if we consider low-metallicity stars polluted by a small number of nucleosynthetic events that took place just before the formation of the stars, the age of the star can be estimated without calling for a complex model of the chemical evolution of the Galaxy. The major difficulty of the methodology is therefore related to the theoretical estimate of the r-production ratio and the corresponding uncertainties of astrophysics and nuclear physics origin that may affect this prediction.

In this respect, the ²³²Th to ²³⁸U chronometry has been shown to be relatively robust, in particular in comparison with the Th/Eu chronometry, i.e., to be less affected by the still large astrophysics and nuclear physics uncertainties affecting our understanding of the r-process nucleosynthesis (Goriely & Clerbaux 1999; Cayrel et al. 2001; Frebel et al. 2007). In Table 2 we provide the production ratios of the ²³²Th to ²³⁸U isotopes based on our NS merger simulations. Assuming that r-process enhanced metal-poor stars were enriched by one or a few NS merger events, we can derive ages within this scenario. From the observed ratio of log (U/Th)_obs = −0.94 ± 0.09 for the metal-poor star CS31082-001 (Cayrel et al. 2001; see also Goriely & Arnould 2001 for updated values), we compute the age as Δt = 21.8[log (U/Th)₀ − log (U/Th)_obs] = 15.7 Gyr with the production ratio log (U/Th)₀ = −0.22 ± 0.01 (see Table 2). The observational uncertainty of 0.09 dex in this case dominates the error on the age estimate since it amounts to about 2.0 Gyr, while the theoretical uncertainties associated with the different EOSs (Table 2) give a 0.2 Gyr error only.

Additional uncertainties stem from the nuclear physics aspects of the r-process nucleosynthesis. In particular, the Th and U production is known to be sensitive to the β-decay, neutron capture, and fission rates adopted in the network calculation. For the specific case of the 1.35–1.35 M_☉ merger simulation based on the DD2 EOS, about 20 different abundance calculations based on different nuclear physics ingredients were performed in order to estimate the uncertainties affecting the ²³²Th-to-²³⁸U production ratio. These inputs include (1) different mass models for the estimate of the β-decay and neutron capture rates, namely the recent microscopic Skyrme (Goriely et al. 2010) and Gogny (Goriely et al. 2009a) Hartree–Fock–Bogolyubov (HFB) mass models as well as the macroscopic–microscopic finite-range droplet model (FRDM; Möller et al. 1995); (2) different neutron-nucleus optical potentials (Jeukenne et al. 1977; Koning & Delaroche 2003); (3) different reaction mechanisms, such as the resonant compound or the direct capture models, to estimate the neutron-capture rates (Xu & Goriely 2012); (4) different predictions of the fission barrier heights to estimate the spontaneous, neutron-induced, and β-delayed fission rates, namely the HFB (Goriely et al. 2007) and Thomas–Fermi (Myers & Swiatecki 1999) fission barriers; and (5) different prescriptions for the fission fragment distributions (Kodama & Takahashi 1975; Schmidt & Jurado 2010). The final ²³²Th-to-²³⁸U production ratio is found to lie between 1.41 and 2.30, leading to an age estimate of the metal-poor star CS31082-001 of 14.9 ± 2.3 Gyr. The error associated with nuclear physics uncertainties is therefore of the same order as those affecting the observation and significantly larger than those related to the EOS. The largest age of 17.2 Gyr is obtained when use is made of the Gogny D1M mass model (Goriely et al. 2009a), while the smallest estimate of 12.6 Gyr is found with the FRDM mass predictions. The Thomas–Fermi fission barriers (Myers & Swiatecki 1999) also lead to rather larger age with respect to the HFB predictions. More details on the r-process sensitivity to the nuclear physics input within the NS merger model will be given in a forthcoming paper.

The derived age of this halo star lies within the ballpark of other age estimates (Cayrel et al. 2001). For the metal-poor star HE 1523-0901 (log (U/Th)_obs = −0.86 ± 0.13; Frebel et al. 2007) our ²³²Th-to-²³⁸U production ratio implies an age of about 14.0 ± 2.8 Gyr, which is also within the range of other calculations (Frebel et al. 2007).

Table 2 also lists the production ratios of ²³²Th to ²³⁵U as well as Th to Eu. For age estimates the Th/Eu chronometer has been widely used, although it remains highly sensitive to all types of uncertainties. In this case, the stellar age is derived from Δt = 46.7 Gyr[log (Th/Eu)₀ − log (Th/Eu)_obs], so that a 25% error on the production or observed Th/Eu ratio gives rise to an uncertainty of about 5 Gyr on the stellar age. Special care of the associated uncertainties should therefore be taken when applying this chronometer pair (Goriely & Arnould 2001). Considering our production ratio of about log (Th/Eu)_obs = −0.20 ± 0.06 (see Table 2), we find an age of Δt = 17.7 ± 2.8 Gyr for HE 1523-0901 (log (Th/Eu)_obs = −0.58, with an additional 4.8 Gyr uncertainty based on observation; Frebel et al. 2007). We stress again here that as long as the r-process site remains unidentified, corresponding uncertainties of the production ratios have to be taken into account, as well as uncertainties arising from the nuclear physics input in network calculations. We also point out the possibility that there is a measurable event-to-event variation in the production ratios, for instance in the case of NS mergers caused by the unknown binary configuration (see Table 2), but also for other sites a progenitor dependence cannot be excluded.

The production ratios summarized in Table 2 are of particular importance since NS mergers may well be a major source of r-process elements and current supernova models cannot provide suitable conditions for the formation of the heaviest r-process elements (Hoffman et al. 2008; Janka et al. 2008; Roberts et al. 2010; Hüdepohl et al. 2010; Fischer et al. 2010; Wanajo et al. 2011). The relatively reliable age estimates from the ²³²Th-to-²³⁸U ratio are compatible with the age of the universe, and thus NS mergers cannot be excluded as the source of the contamination of the considered metal-poor stars.

4.4. Nuclear Heating

Another outcome of our nucleosynthesis calculations is the determination of the heating due to radioactive decays in the ejecta. This is in particular important because the radioactive heating provides the energy source for an optical counterpart associated with NS mergers (see Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010b, Section 5). Our calculations allow us to check the robustness and general behavior of the heating rate. In all models the heating rate due to radioactive decays (beta-decays, fission, and alpha-decays) looks similar (see Figure 3 in Goriely et al. 2011). For instance at the time t_peak, when the luminosity of the optical transient reaches its maximum (typically several hours; see Section 5), the heating rate varies from 3 × 10¹⁰ erg g⁻¹ s⁻¹ to 1 × 10¹¹ erg g⁻¹ s⁻¹ for the cases where detailed nucleosynthesis calculations were made. This implies that the heating efficiency $f\equiv \dot{Q}(t_{\mathrm{peak}})t_{\mathrm{peak}}/{M_{\mathrm{ej}}c^2}$ that enters the estimates of optical emission properties (see Section 5 and Li & Paczyński 1998; Metzger et al. 2010b) varies between 1.1 × 10⁻⁶ and 1.8 × 10⁻⁶ (see Table 2). One observes a moderate dependence of f on the EOS and the mass ratio. This is caused by the longer-duration t_peak of the emission peak for larger ejecta masses in models with soft EOSs or asymmetric binaries. Overall, f ≈ (1.5 ± 0.3) × 10⁻⁶ seems to be a fair approximation, which is half of the value suggested in Metzger et al. (2010b) for a longer peak time of about one day. As detailed in Section 5, we expect from our relativistic merger simulations shorter peak times in the range of 2–7 hr for symmetric binaries. Note that a fraction of the energy of radioactive decays is released in gamma-rays, of which a fraction may escape from the ejecta without efficient thermalization (see Kulkarni 2005; Metzger et al. 2010b). The assumption of a full thermalization of gamma-rays implies that the factor f is somewhat overestimated. The total amount of energy released by radioactive decays (without neutrinos) is in the range of about 3.2 ± 0.2 MeV per nucleon (see Table 2).

5. OPTICAL COUNTERPARTS AND RADIO REMNANTS

The radioactive decay of the synthesized r-process elements generates heat, which is deposited in the ejecta and powers an optical display (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010b; Roberts et al. 2011; Goriely et al. 2011). This electromagnetic counterpart of an NS merger is potentially observable with existing and upcoming optical surveys such as the Palomar Transient Factory, the Synoptic All Sky InfraRed Survey, the Panoramic Survey Telescope and Rapid Response System, and the Large Synoptic Survey Telescope (see, e.g., Strubbe & Quataert 2009 and Metzger & Berger 2012 for a compilation of certain characteristics of these facilities). We also refer to Kann et al. (2011) for a report on attempts to detect signatures of a radioactively powered transient in the light curve following a gamma-ray burst. The detection of such optical signals would provide valuable information on ejecta properties and the sky position of an event. As detailed in the next section, the peak luminosity, peak time, peak width, and effective temperature depend on the amount of ejecta and the expansion velocity. From an observation one might therefore derive these ejecta characteristics, which are otherwise only accessible by numerical simulations. Detecting radioactively powered emission and determining ejecta masses will consolidate the role of NS mergers for the enrichment of the Galaxy with heavy r-process elements. A precise localization of an unambiguously identified optical transient of an NS merger event will help improve the sensitivity of gravitational-wave detections and will provide information about the host galaxy or environment. Moreover, if the distance scale of the events can be constrained, better observational limits on the merger rate will become available.

5.1. Model

The bolometric peak luminosity of an optical transient associated with an NS coalescence can be estimated by

$\begin{equation} L_{\mathrm{peak}}\approx 5\times 10^{41} {{\rm erg\ s}^{-1}} \left(\frac{f}{10^{-6}} \right)\left(\frac{v}{0.1 c} \right)^{1/2}\left(\frac{M_{\mathrm{ej}}}{10^{-2}\ M_{\odot }} \right)^{1/2} \end{equation} \tag{ 2 }$

with the average outflow velocity v, the ejecta mass M_ej, and the heating efficiency f already introduced in Section 4.4 (Arnett 1982; Li & Paczyński 1998; Metzger et al. 2010b).

The time of the peak luminosity and the effective temperature at the time of the maximum luminosity can also be expressed as functions of the outflow velocity, the ejecta mass, and f:

$\begin{equation} t_{\mathrm{peak}}\approx 0.5\, \mathrm{days} \left(\frac{v}{0.1 c} \right)^{-1/2}\left(\frac{M_{\mathrm{ej}}}{10^{-2}\ M_{\odot }} \right)^{1/2}, \end{equation} \tag{ 3 }$

$\begin{equation} T_{\mathrm{peak}}\approx 1.4\times 10^4\, \mathrm{K} \left(\frac{f}{10^{-6}} \right)^{1/4} \left(\frac{v}{0.1 c} \right)^{-1/8}\left(\frac{M_{\mathrm{ej}}}{10^{-2}\ M_{\odot }} \right)^{-1/8} \end{equation} \tag{ 4 }$

(see Li & Paczyński 1998; Metzger et al. 2010b). Within this model the width Δt_peak of the luminosity peak is proportional to t_peak (Arnett 1996; Kasen & Woosley 2009). On the basis of the one-zone model in Goriely et al. (2011) we find that the FWHM can be very well approximated as

$\begin{equation} \Delta t_{\mathrm{peak}} \approx 2.5 t_{\mathrm{peak}}. \end{equation} \tag{ 5 }$

The above formulae can be understood by some general considerations. For larger M_ej or smaller v the ejecta need a longer time to become transparent. If the outflow gets optically thin at a later time, expansion cooling reduces the effective temperature at L_peak. Since, however, the increase of the emission radius R_peak ≈ v × t_peak dominates the decrease of the temperature in the Stefan–Boltzmann law, $L_{\mathrm{peak}}\propto T_{\mathrm{peak}}^4 R_{\mathrm{peak}}^2$ , the peak luminosity increases with larger M_ej and higher v. The simple scaling laws of Equations (2)–(4) for the properties of optical transients are confirmed by calculations in Metzger et al. (2010b), Roberts et al. (2011), and Goriely et al. (2011).

5.2. EOS Dependence

The panels on the left side of Figure 11 display the peak luminosity, the peak time, and the effective temperature of electromagnetic counterparts of 1.35–1.35 M_☉ mergers for different EOSs, which are characterized by the radii of the corresponding 1.35 M_☉ NSs. A clear dependence of the optical display on the compactness of the NSs can be seen, with soft EOSs yielding brighter transients, which peak on longer timescales with a lower effective temperature. The relatively clear relations are mainly a consequence of the strong EOS impact on M_ej (Figure 3), while the average outflow velocity varies only by a factor of three (see Figure 12 and Table 1). The average expansion velocities show the tendency of being higher for EOSs which yield smaller R_1.35 (see Figure 12). This is consistent with the reasoning used before that more compact NSs lead to more violent collisions. However, the relation between the average outflow velocity and R_1.35 is not very tight. For symmetric binaries the outflow velocities (measured 10 ms after the merging at which time the asymptotic values are fairly well determined) vary (with larger scatter) between 0.16 and 0.45 times the speed of light.

**Figure 12.** Average ejecta expansion velocity for 1.35–1.35 M_☉ mergers (with symbols analogous to Figure 3) and for 1.2–1.5 M_☉ (red squares) for different EOSs characterized by the corresponding radius R_1.35 of a nonrotating NS with a mass of 1.35 M_☉.
Download figure:
Standard image High-resolution image

Also for asymmetric 1.2–1.5 M_☉ binaries an EOS dependence of optical counterpart properties is observed (right side of Figure 11). As in the case of equal-mass mergers, EOSs with smaller R_1.35 lead to more mass ejection and therefore to more luminous events with longer t_peak and lower T_eff. Compared to symmetric binaries, asymmetric setups generally produce brighter transients, which reach their peak luminosities on longer timescales and therefore lower effective temperatures, because their ejecta masses are higher whereas the expansion velocities are comparable with those of symmetric systems described by the same EOS (see Figure 12). For all considered quantities the asymmetric models exhibit a milder EOS dependence. Plotting the peak luminosity as a function of the binary mass ratio and the total binary mass for the NL3, DD2, and SFHO EOSs reveals a qualitatively similar behavior as the ejecta masses shown in Figure 7.

The relations shown in Figure 11 suggest the possibility to constrain NS radii and thus the high-density EOSs from observations of optical transients associated with NS mergers. Optimally such a detection could be supplemented by a gravitational-wave measurement, which provides the involved binary masses, the distance, and the merger time. But even without an associated gravitational-wave signal the observable features of a transient may have the potential to yield constraints for the NS EOS. The combinations of L_peak, t_peak, and T_eff vary systematically with the NS properties. For instance a low peak luminosity, a small peak width, and a high effective temperature imply a large NS radius.

5.3. Implications for Observations

Symmetric 1.35–1.35 M_☉ binaries are predicted to be the most common configurations and thus are likely to be the ones first and most frequently observed. Unfortunately, the 1.35–1.35 M_☉ systems yield the smallest ejecta masses and thus the lowest luminosities and peak timescales. The peak widths are important to estimate the prospects of blind searches, whereas the peak time sets the scale for the response time after a gravitational-wave trigger.

The nearly complete coverage of EOS possibilities by our survey allows us to determine the possible range of signal properties of optical transients associated with NS merger events. From our survey we find that the optical peak luminosity of a 1.35–1.35 M_☉ NS merger should be expected to be between about 3 × 10⁴¹ erg s⁻¹ and 14 × 10⁴¹ erg s⁻¹, corresponding to absolute bolometric magnitudes of M = −15.0 and M = −16.7. The peak times range from only 2 hr to 7 hr, and the duration of the emission is expected to be between 4.8 hr and 18 hr depending on the high-density EOS. Note that the ballpark of our models yields fainter and shorter transients than typical estimates based on Newtonian models, which in general obtain higher ejecta masses and lower average expansion velocities (Roberts et al. 2011; Korobkin et al. 2012; Rosswog et al. 2013; Piran et al. 2013; Rosswog 2013; see also the discussion in Section 3). While for the peak luminosity these differences lead to partially compensating effects, the timescales are more strongly affected. For symmetric binaries even the maximum peak time of about 7 hr found in our sample is well below most predictions based on Newtonian models (Metzger et al. 2010b; Roberts et al. 2011; Piran et al. 2013; Rosswog et al. 2013).

As shown in Figure 6, during the early stages of the expansion the ejecta exhibit a fair asymmetry between polar and equatorial directions. The formulae used in this section for estimating the properties of optical counterparts assume spherically symmetric outflows. It remains to be explored whether the donut-like shape visible in Figure 6 persists at late times or whether the outflow becomes more symmetric when the peak of the optical display occurs. Multidimensional radiation transport calculations coupled to long-term hydrodynamical simulations are required to determine to which extent the simplified emission model provides reliable estimates of the observable properties of the electromagnetic transients in dependence on the observer direction.

The latest detailed atomic models (Kasen et al. 2013; Barnes & Kasen 2013) predict that the opacity of r-process elements may be significantly enhanced compared to iron group elements, whose opacities have been adopted in Equations (2)–(4). The effects of an increased opacity can be readily estimated because the opacity κ enters the prefactors as κ^−1/2 in Equation (2), as κ^1/2 in Equation (3), and as κ^−3/8 in Equation (4) (see Metzger et al. 2010b). While the peak luminosity of the transient is reduced and its duration stretched, the differences between different EOSs remain. By observing two signal features like, for instance, the peak luminosity and the peak width, one can in principle remove the degeneracy due to the uncertain opacity because every observable individually exhibits a specific approximate correlation with R_1.35, where the relation is known except for a constant factor.

5.4. Radio Flares

Another potentially observable phenomenon connected with NS mergers is radio emission that is produced by the interaction of the outflowing material with the ambient medium (Nakar & Piran 2011; Piran et al. 2013; Rosswog et al. 2013). The ejecta properties which determine the appearance of a radio remnant are the kinetic energy and the outflow velocity. The peak flux density was computed to be proportional to the total kinetic energy E_kin of the outflow and to the (initial) outflow velocity v to a power of about 2.5 (Piran et al. 2013). For symmetric mergers we find values for the kinetic energy between 6 × 10⁴⁹ erg and 10⁵¹ erg (Figure 13 and Table 1), whereas the average outflow velocities vary from 0.16 to 0.45 times the speed of light (Figure 12). The kinetic energy scales well with the expansion velocity, i.e., the models with the highest outflow velocities yield also the highest kinetic energies, and configurations with smaller v result in lower kinetic energies. This implies that the peak flux of radio remnants is uncertain by a factor of 200 because of the variations in the kinetic energy and expansion velocity associated with the incomplete knowledge of the NS EOS. Similar values are found for asymmetric mergers. Moreover, the circumburst densities affecting the signal brightness and length are uncertain, and low densities are likely to reduce the radio detectability (Berger et al. 2005; Soderberg et al. 2006; Metzger & Berger 2012; Fong et al. 2012, 2013).

**Figure 13.** Kinetic energy of the ejecta for 1.35–1.35 M_☉ mergers (with symbols analogous to Figure 3) and for 1.2–1.5 M_☉ (red squares) for different EOSs characterized by the corresponding radius R_1.35 of a nonrotating NS with a mass of 1.35 M_☉.
Download figure:
Standard image High-resolution image

6. SUMMARY AND CONCLUSIONS

We have performed relativistic hydrodynamical simulations of NS mergers to investigate the mass ejection, the nucleosynthesis outcome, and the properties of associated optical transients. The main goal of this study was to explore the systematics of the EOS dependence of these aspects by employing a large set of candidate EOSs, while focusing mostly on binaries with two 1.35 M_☉ NSs and on asymmetric 1.2–1.5 M_☉ systems.

The unbound ejecta mass is strongly affected by the adopted EOS. We find that the NS compactness is the crucial EOS parameter determining the ejecta properties. Using the radius R_1.35 of a nonrotating 1.35 M_☉ NS to characterize different EOSs, we find that "softer" EOSs which yield more compact NSs tend to produce more ejecta. The ejecta masses are between 10⁻³ M_☉ and 1.5 × 10⁻² M_☉ depending on the EOS and binary mass ratio (Figures 3 and 7). Most of the unbound material originates from the contact interface between the colliding stars for symmetric as well as asymmetric binaries. In the latter case about 25% of the mass ejection is shed off the outer end of the spiral, arm, like tail into which the lower-mass component is stretched during its final approach to collision with the more massive companion.

A qualitative and quantitative agreement of our SPH simulations (employing the conformal flatness approximation) with fully relativistic grid-based simulations is found, whereas considerable differences compared to Newtonian models concerning the origin and the amount of ejecta are observed. The pronounced spiral arms, for example, which form in Newtonian simulations during the merging of symmetric binaries and whose mass stripping dominates the ejecta, are absent in relativistic mergers of equal-mass NSs. Newtonian models therefore tend to produce considerably higher ejecta masses.

When temperature effects are mimicked by adding a thermal ideal-gas component with a constant ideal-gas index Γ_th to EOS models which are provided as zero-temperature barotropes, the best match of the ejecta masses with fully consistent calculations is achieved for a relatively low value of Γ_th = 1.5. This is in conflict with values of 1.8 or 2 which have been widely used and work well for gravitational-wave determinations (Bauswein et al. 2012, 2010a). This can be understood from the fact that the gravitational-wave signal is produced by the bulk of the merger mass in the high-density regime, whereas the ejecta depend on the thermodynamics of lower-density matter expanding away from the colliding stars and being accelerated by pressure forces.

The binary parameters have qualitatively the same influence on the ejecta masses for all investigated EOSs. A larger binary mass asymmetry leads to a strong increase of the mass of unbound matter, whereas a higher total binary mass results in larger ejecta masses for asymmetric mergers but only a weak increase of M_ej for symmetric systems (Figure 7). The occurrence of a prompt collapse of the merger remnant is associated with a significant drop in the ejecta mass. For a given EOS (and no prompt collapse) the smallest amount of ejecta is to a good approximation produced by 1.35–1.35 M_☉ binaries, while the ejecta mass can be up to a factor 10 higher for other binary configurations. For soft EOSs the ejecta masses show steeper gradients in the binary parameter space (q, M_tot).

Nuclear network calculations show that r-process elements with mass numbers above 130 are robustly produced in the ejecta of NS mergers for all investigated EOSs. The vast majority of ejected material is fission recycled, producing a final mass-integrated abundance pattern that resembles closely the solar composition of r-process elements. The robustness with respect to variations of the high-density EOS confirms that NS mergers are a very promising source of r-process elements.

For some EOS models which do not provide the temperature dependence consistently, the hydrodynamical simulations are based on an approximate treatment of thermal effects and the nucleosynthesis calculations rely on an estimation of the temperature in the ejecta. Also in these cases r-process elements are produced with solar distribution, showing the insensitivity to the exact temperature value at the beginning of the network calculations. Moreover, the abundance patterns are not affected by asymmetries in the binary setup. The amazing insensitivity of the outcome of the nucleosynthesis processes can be understood as a consequence of the large neutron excess in the ejected inner-crust material of the merging NSs, which allows for fission recycling in essentially all considered conditions.

By folding with the binary population, we identify the results of 1.35–1.35 M_☉ mergers as a good approximation for the average ejecta mass per merger event. The main uncertainty in the average ejecta mass per merger is therefore associated with the incomplete knowledge of the NS EOS, which implies variations in the average ejecta mass of about a factor of 10. The observed abundance of r-process matter in the inventory of our Galaxy can be accounted for with an NS merger rate that is compatible with current predictions based on population synthesis and pulsar observations (Abadie et al. 2010). Final conclusions, however, require a more precise determination of the merger rate, e.g., by gravitational-wave detections or observations of electromagnetic counterparts. Our work implies that in addition to more accurate information on the merger rate also better constraints on the high-density EOS are needed to decide whether NS mergers are a major (or the dominant) source of r-process elements.

Moreover, our simulations in combination with estimates of the Galactic r-process material provide independent evidence that the Galactic merger rate cannot be higher than approximately 4 × 10⁻⁴ events per year if NS mergers should not overproduce heavy r-nuclei compared to observations.

The nucleosynthesis calculations of our survey also provide important information on the production ratios of certain isotopes which are used for nucleocosmochronometry. For instance the ratio of 232-thorium to 238-uranium is found to be about 1.65 with only small variations depending on the high-density EOS and the binary configuration (Table 2). Using this result, we derived ages of metal-poor stars which are consistent with other age estimates. This implies that NS mergers are not excluded as r-process element sources for metal-poor stars and the production ratios provided here for the first time in the NS merger context should be taken into account in stellar age estimates as long as mergers cannot be excluded as r-process sites in the early Galactic history.

Just as the ejecta masses of binary NS mergers exhibit a strong sensitivity to the properties of the nuclear EOS and thus to the radius R_1.35 of the merging stars, we also predict the optical transients powered by the radioactive energy release in the ejecta to depend on the compactness of the binary components. EOSs which lead to smaller NS radii produce more ejecta and therefore cause brighter optical counterparts, which peak on longer timescales with longer durations and with lower effective temperatures. On the basis of our extensive survey of EOSs, which suggests clear correlations between observable features (luminosity, peak timescale, effective temperature) and NS radii, we propose that optical observations of transients associated with NS mergers could yield valuable constraints on the NS EOS.

The very broad range of possibilities included in our EOS sample allows us to bracket the expected range of signal features of optical counterparts associated with NS mergers. Optical transients of 1.35–1.35 M_☉ mergers should (at least) reach an absolute bolometric peak magnitude between −15.0 and −16.7 (3 × 10⁴¹ erg s⁻¹ and 14 × 10⁴¹ erg s⁻¹). Depending on the high-density EOS, the peak times vary from 2 to 7 hr, implying durations of about 4 to 18 hr, whereas effective temperatures between 1.3 × 10⁴ K and 1.9 × 10⁴ K can be expected. We emphasize that the peak luminosities, peak times, and peak widths of the optical counterparts are found to be considerably lower in our analysis compared to earlier investigations based on Newtonian models (Metzger et al. 2010b; Roberts et al. 2011; Piran et al. 2013; Rosswog et al. 2013). The reduction is a consequence of the smaller ejecta masses especially for symmetric binaries. Because of the shorter peak time and duration of the optical transient suggested by relativistic merger results, we also find a smaller fraction f of radioactive decay energy relative to the rest-mass energy of the ejecta. While Newtonian merger models yield f ∼ 3 × 10⁻⁶ at the time of the luminosity peak (Metzger et al. 2010b), we obtain f ∼ 1.5 × 10⁻⁶ with little sensitivity to the EOS and the binary parameters.

For different EOSs considerable differences are also found in the properties of radio remnants. Based on our sample of models, we estimate an uncertainty of up to a factor of 200 for the theoretical predictions of the brightness of these events.

Future work should address a variety of issues. The hydrodynamical models of NS mergers should include magnetic fields, and the effects of neutrino interactions should be explored. Our results also need to be confirmed by fully relativistic merger simulations. The properties of emitted electromagnetic radiation should be computed for detailed multi-dimensional outflow models including the corresponding nuclear network calculations to determine the composition and heating. Radiative transfer calculations will have to be performed to study the observational appearance of the potentially anisotropic ejecta dependent on the viewing direction. This will require employing appropriate opacities of the r-processed material as recently derived by Kasen et al. (2013). The corresponding opacity increase by orders of magnitude leads to considerably lower peak luminosities and longer peak timescales than estimated in our work. However, in view of the faster expansion timescales and the lower ejecta masses of our relativistic models, the very long durations of radioactive transients calculated by Barnes & Kasen (2013) (on the basis of assumed outflow properties) appear to be on the extreme side. In this context it will also be important to determine the contribution of mass ejection from the secular evolution of the merger remnant (a black hole–torus system or hypermassive NS), which will lose mass through neutrino-driven and magnetohydrodynamical outflows. The corresponding matter will increase the ejecta mass (see, e.g., Fernández & Metzger 2013 for long-term evolution models) and will ultimately have to be taken into account for reliable predictions of the properties of electromagnetic counterparts of NS mergers. However, the differences in the composition and velocities of this additional mass outflow are predicted to lead to a second light curve peak that can be observationally discriminated from the signature of the dynamical ejecta (Barnes & Kasen 2013). Therefore, the potential of inferring EOS information from the counterpart of the dynamical ejecta is likely to remain unaffected. The robustness of the nucleosynthesis outcome has to be explored concerning variations connected to uncertainties of the nuclear reaction rates. Further work is also needed to address how NS mergers as r-process sources fit into chemical evolution scenarios of the Milky Way, which should explain the observations of r-element-enhanced metal-poor stars. Finally, the capabilities of various observational facilities have to be evaluated in view of the bounds on the observable features set by our survey.

We thank Matthias Hempel for helpful discussions and for providing his EOS tables. This work was supported by Sonderforschungsbereich Transregio 7 "Gravitational Wave Astronomy," and the Cluster of Excellence EXC 153 "Origin and Structure of the Universe" of the Deutsche Forschungsgemeinschaft, and by CompStar, a research networking program of the European Science Foundation. S.G. is F.R.S.-FNRS Research Associate. Computing resources provided by the Rechenzentrum Garching of the Max-Planck-Gesellschaft and Leibniz-Rechenzentrum Garching are acknowledged.

SYSTEMATICS OF DYNAMICAL MASS EJECTION, NUCLEOSYNTHESIS, AND RADIOACTIVELY POWERED ELECTROMAGNETIC SIGNALS FROM NEUTRON-STAR MERGERS

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL MODEL AND EQUATIONS OF STATE

3. EJECTA MASSES