NEUTRINO SIGNATURES AND THE NEUTRINO-DRIVEN WIND IN BINARY NEUTRON STAR MERGERS

The MAGMA code (Rosswog & Price 2007) is a state-of-the-art SPH code that contains physics modules that are relevant to compact binary mergers and on top implements a slew of numerical improvements over most "standard" SPH schemes. A detailed code description can be found in Rosswog & Price (2007), while results are presented in Price & Rosswog (2006) and Rosswog (2007).

We briefly summarize the most important physics modules. For the thermodynamic properties of neutron star matter we use a temperature-dependent relativistic mean-field EOS (Shen et al. 1998a, 1998b). It can handle temperatures from 0 to 100 MeV, electron fractions from Y_e = 0 (pure neutron matter) up to 0.56, and densities from about 10 to more than 10¹⁵ g cm⁻³. No attempt is made to include matter constituents that are more exotic than neutrons and protons at high densities.

The MAGMA code contains a detailed multiflavor neutrino-leakage scheme. An additional mesh is used to calculate the neutrino opacities that are needed for the neutrino emission rates at each particle position. The neutrino emission rates are used to account for the local cooling and the compositional changes due to weak interactions such as electron captures. A detailed description of the neutrino treatment can be found in Rosswog & Liebendörfer (2003).

The self-gravity of the fluid is treated in a Newtonian fashion. Both the gravitational forces and the search for the particle neighbors are performed with a binary tree that is based on that described in Benz et al. (1990). These tasks are the computationally most expensive part of the simulations and in practice they completely dominate the CPU-time usage. Forces emerging from the emission of gravitational waves are treated in a simple approximation. For more details, we refer to the literature Rosswog et al. (2000); Rosswog & Davies (2002).

In terms of numerical improvements over "standard" SPH techniques, the code contains the following:

• 1.  
  To restrict shocks to a numerically resolvable width, artificial viscosity is used. The form of the artificial viscosity tensor is oriented at Riemann solvers (Monaghan 1997). The resulting equations are similar to those constructed for Riemann solutions of compressible fluid dynamics. In order to apply the artificial viscosity terms only where they are really needed, i.e., close to a shock, the numerical parameter that controls the strength of the dissipative terms is made time dependent, as suggested by Morris & Monaghan (1997). An extra equation is solved for this parameter, which contains a source term that triggers on the shock and a term causing an exponential decay of the parameter in the absence of shocks. Tests can be found in Rosswog et al. (2000) and an illustration of the time-dependent viscosity parameter in the context of Sod's shock-tube problem can be found in (Rosswog et al. 2008, their Figure 1).
• 2.  
  A self-consistent treatment of extra terms in the hydrodynamics equations that arise from varying smoothing lengths (the so-called grad-h terms; Springel & Hernquist 2002; Monaghan 2002; Price 2004; Rosswog & Price 2007).
• 3.  
  A consistent implementation of adaptive gravitational softening lengths as described in Price & Monaghan (2007).
• 4.  
  The option to evolve magnetic fields via so-called Euler potentials (Stern 1970) that guarantee that the constraint ∇ · B = 0 is fulfilled. The details can be found in Rosswog & Price (2007) and Rosswog & Price (2008).

Our work starts from three-dimensional SPH simulations of BNS mergers in an isolated binary system, made up of two individual neutron stars, each with a baryonic mass of 1.4 M_☉ (equivalent to a gravitational mass of ∼ 1.3 M_☉), and separated by a distance of 48 km. Their initial hydrostatic structure is computed by solving the Newtonian equations. The properties of nuclear matter (pressure, energy, entropy, etc.) are obtained by interpolation in the Shen EOS (Shen et al. 1998a, 1998b), assuming zero temperature and β equilibrium. A large number of SPH particles (N ≈ 10⁶) are then mapped onto the resulting density profiles. These particles are further relaxed to find their true equilibrium state (see, e.g., Rosswog & Price 2007).

Because the tidal synchronization time exceeds the gravitational decay time, BNSs cannot be tidally locked during their inspiral phase, and are thus expected to be close to irrotational (see, e.g., Bildsten & Cutler 1992 and Kochanek 1992); however, past investigations of the coalescence of BNSs suggest that the intrinsic spin of each neutron star is an important parameter, as it determines, for example, the width of the baryon-poor funnel above the forming central object that is thought to play a decisive role in the launching of a GRB (see, e.g., Zhuge et al. 1996; Faber et al. 2001; Oechslin et al. 2007). Therefore, to encompass the range of possible outcomes, we simulate three initial configurations: (1) no spins (no intrinsic neutron star spin), (2) co-rotating spins (each neutron star has a period equal to that of the orbit, and rotates in the same direction), and (3) counter-rotating spins (same spin period as for co-rotating spins, but with rotation in a reverse direction to that of the orbit; the spins are counter-aligned). In the rest of this paper, we broadly refer to these models as the no-spin, the co-rotating spin, and the counter-rotating spin BNS models. The co-rotating and counter-rotating cases are included here to bracket the range of possibilities, but they represent unlikely extremes. We leave to future work the consideration of higher-mass BNS mergers (Rosswog & Ramirez-Ruiz 2003) or asymmetric systems (Rosswog et al. 2000).

In practice, we post-process these three-dimensional SPH simulations by constructing azimuthal-averages, selecting a time after coalescence (when the two neutron stars come into contact) of 10.0 ms (no-spins), 12.4 ms (co-rotating spins), and 20.4 ms (counter-rotating spins), corresponding in each VULCAN/2D simulations to the time origin. At these times, these systems are evolving toward, but have not yet reached, complete axisymmetric configurations. This is a compromise made to model with VULCAN/2D the epoch of peak neutrino brightness. At these initial times, the amount of mass with a density greater than 10¹³ g cm⁻³ is 2.5 M_☉ in the no-spin BNS model, 2.31 M_☉ in the co-rotating spin model, and 2.43 M_☉ in the counter-rotating spin model. In the same model order, and at these times, the mass contained inside the ∼ 13 km radius of the highest mass neutron star allowed by the Shen EOS is only 1.11, 1.26, and 1.03 M_☉, respectively. The SPH simulations employed to generate these models are Newtonian, and, thus, one would expect more mass at small radii and systematically higher densities if allowance for GR effects were made, but the use of the highly incompressible Shen EOS, together with the significant amount of mass in quasi-Keplerian motion and on wide orbits, suggests that the systems should survive at least for some time before the general-relativistic gravitational instability occurs.

Starting from fluid variables constructed with azimuthal averages of the three-dimensional SPH simulations described above, we follow the evolution of three BNS merger configurations with the Newtonian hydrodynamic code VULCAN/2D (Livne 1993), supplemented with an algorithm for neutrino transport as described in Livne et al. (2004) and Walder et al. (2005). The version of the code used here is the same as that discussed in Dessart et al. (2006a) and Burrows et al. (2006a), and uses the two-dimensional MGFLD method to handle neutrino transport (see Appendix A of Dessart et al. 2006a). Additional details on VULCAN/2D are provided in Burrows et al. (2007b). Magnetic fields are ignored in this study, as is the potential role of physical viscosity, i.e., there is no viscosity included besides one of an inherently numerical nature. The gravitational potential is assumed to be Newtonian and computed using a multipole solver (Dessart et al. 2006a). Our study presents a self-consistent treatment of neutrino emission, scattering, and absorption in a multispecies MGFLD context. Associated neutrino–matter coupling source terms are included in the momentum and energy conservation equations and their relevance to the merger evolution is, thus, computed explicitly for a typical duration of ≳ 100 ms. For these dynamical calculations, we include all neutrino processes described in Burrows et al. (2006b), but neglect the secondary processes of neutrino–electron scattering. The transport solution in the MGFLD scheme we employ is solved for the three neutrino species ν_e, ${\bar \nu }_e$, and "ν_μ," (which groups together the μ and τ neutrinos) and at eight neutrino energies: 2.50, 6.87, 12.02, 21.01, 36.74, 64.25, 112.36, and 196.48 MeV. The energy spacing is constant in the log.

The choice of the spatial grid is determined primarily by the very aspherical density distribution of the merger and the very steep density drop-off at the neutron star surface in the polar direction. We present in Figures 2 and 3 (top row) the initial temperature (logarithmic scale and in MeV), electron fraction, and the density distribution (white-line contours overplotted for every decade between 10⁷ and 10¹⁴ g cm⁻³) for each BNS merger configuration that we map onto the VULCAN/2D grid. Note that fast rotation in the inner region, although sub-Keplerian, displaces the density maximum by 8 km from the center and along the equator in the no-spin BNS model. These density peaks also correspond to extrema in the electron fraction of ∼ 0.1 at this time. To setup the hybrid VULCAN/2D grid, we position the transition radius between the inner Cartesian and the outer spherical polar regions at 12 km. The density at this location is at all times in excess of 10¹⁴ g cm⁻³, and thus—due to the stiffness of the Shen EOS in that regime—away from the regions of large density gradient. The resolution in this inner region is typically 200 × 200 m², corresponding to 12/0.2 = 60 zones in each (cylindrical) direction r and z from the center to the transition radius. Beyond the transition radius, we use a logarithmically increasing radial grid spacing with Δr/r = 1.85%, using 301 zones to 3000 km. The grid covers one hemisphere, from the rotation axis to the equator (both treated as reflecting boundaries), and the computation assumes axisymmetry, i.e., the azimuthal gradients are zero. We fill the grid outside the merger with material having a low density of 5 × 10³ g cm⁻³ and a low temperature of 4 × 10⁸ K. Below Y_e = 0.05, we compute thermodynamic variables, as well as opacities/emissivities, by adopting an electron fraction of 0.05. Given the smooth variation for the corresponding quantities (pressure, entropy, mean-free path, etc.) at this level of neutron richness, this approximation is expected to be quite good.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For example, the difference in any of these quantities between a neutron fraction (1 − Y_e) of 95% and 100% (the latter would be the case for Y_e = 0.0) is approximately 5% × (the neutron contribution − the proton contribution), corrected for the electron contribution difference, which is minuscule for cross sections. For the pressure, since the proton contribution at high densities is within ∼ 30% of the neutron contribution and at low densities is almost equal to the neutron contribution (ideal gas),⁷ the error is less than ∼ 1.5–2.0%. For the scattering opacities, the universality of the neutral current makes this net error even smaller. For the electron–neutrino absorptive opacities, a high neutron fraction is accompanied by a low final-state proton blocking. For the anti-electron–neutrino absorption, a low proton fraction is accompanied by a high final-state neutron blocking factor. The net result is an absorption cross section at high densities that is only a few percent different for neutron fractions of 95% and 100%.

After the MGFLD radiation-hydrodynamics simulations are completed, we post-process individual snapshots keeping the hydrodynamics frozen, and relaxing only the radiation quantities. This operation is performed for each model at 10 ms intervals over the whole evolution using the multiangle, S_n, radiation-transport module of Livne et al. (2004) and described more recently in Ott et al. (2008). We refer the reader to those papers for details on the method and the nomenclature. This S_n variant yields the explicit angular dependence of the neutrino-specific intensity, and, thus, represents a more accurate modeling of the strongly anisotropic neutrino luminosity in such highly aspherical systems (Ott et al. 2008). The S_n calculation is done with the same number of energy groups (i.e., eight) and on the same spatial grid. In the S_n algorithm, the transition at depth to MGFLD described in Ott et al. (2008) is natural here, since it occurs at 12 km, and, thus, in regions where the densities are nuclear. However, the spatial resolution is not fully satisfactory at the neutron star surface and along the polar direction, a region where the density decreases by few orders of magnitude in just a few kilometers. Thus, this represents a challenge for radiation transport on our Eulerian grid. In the first instance, we relax the MGFLD results with S_n, adopting 8 ϑ-angles. In the S_n approach, n such ϑ-angles translate into n(n + 2)/2 directions mapping the unit sphere uniformly. To obtain the steady-state angle-dependent radiation field, we apply the S_n algorithm to initial conditions resulting from the converged MGFLD model. The solution is found by iterating on the radiation field at each time step until it does not change by more than one part in a thousand at any grid location. This newly found radiation field is then evolved in time through a number of consecutive time steps, until the solution has propagated through the entire grid. In the optically thick regions, the S_n and the MGFLD solutions are inherently close, so convergence is reached after a few ms, although about 10 iterations are required for convergence at each time step. In the optically thin regions, and in the regions intermediate between these two regimes, the convergence is obtained at each time step after 3–4 iterations, but the two solvers yield quite different results, so the S_n run has to be evolved for at least 10 ms, allowing the new multiangle radiation field to "propagate out" and fill the entire grid. Hence, after having reached 10 ms, the S_n run is stopped whenever the radiation field does not change by more than one part in a thousand at any location and from one time step to the next. Once converged, we perform an accuracy check by remapping the angle-dependent neutrino radiation field from 8 to 16 ϑ-angles, and then relaxing this higher angular-resolution simulation. Due to the increased computational costs, we use the high (16 ϑ-angles) angular-resolution configuration only for snapshots at 10, 60, and 100 ms, but for all three BNS merger configurations studied here. Importantly, the explicit angle dependence of the neutrino-radiation field allows us to compute its various moments, which we employ in a novel formulation of the neutrino–antineutrino annihilation rate (Section 4.2). Surprisingly, and as documented in Table 1 (see also Section 4.1), this sophisticated, i.e., multidimensional, multigroup, multiangle, approach that we follow yields results that agree favorably with those obtained using the leakage scheme of Ruffert et al. (1996, 1997), which by contrast is gray and makes rough estimates of directional diffusion times, scale heights, etc. (the expression of the diffusion time is itself, by essence, an approximation, based on random-walk arguments). Our computation of the angle dependence of the radiation field offers, however, a number of interesting new insights, which we explore in Section 4.

Table 1. $\nu _i{\bar \nu }_i$ Annihilation Rates Using the Leakage and the S_n Schemes

	No Spins		Co-Rotating Spins		Counter-Rotating Spins
	∫dVQ+(${\nu _e}{\bar \nu }_e$)	∫dVQ+("${\nu _\mu }{\bar \nu }_\mu$")	$\int dV Q^+({\nu _e}{\bar \nu }_e$)	∫dVQ+("${\nu _\mu }{\bar \nu }_\mu$")	$\int dV Q^+({\nu _e}{\bar \nu }_e$)	∫dVQ+("${\nu _\mu }{\bar \nu }_\mu$")
Time (ms)	(B s−1)		(B s−1)		(B s−1)
	Leakage Scheme
10	1.56(−1)	3.28(−5)	2.05(−1)	1.01(−4)	2.18(−1)	1.40(−4)
20	9.28(−2)	5.99(−6)	4.79(−2)	2.38(−6)	1.05(−1)	1.01(−5)
30	3.80(−2)	1.15(−6)	3.70(−2)	1.43(−6)	5.95(−2)	3.35(−6)
40	3.11(−2)	8.54(−7)	4.67(−2)	1.08(−6)	3.65(−2)	1.36(−6)
50	2.82(−2)	6.23(−7)	5.67(−2)	1.35(−6)	2.16(−2)	5.81(−7)
60	1.78(−2)	2.57(−7)	4.43(−2)	8.74(−7)	1.50(−2)	3.00(−7)
70	1.32(−2)	1.63(−7)	2.69(−2)	3.78(−7)	1.18(−2)	1.79(−7)
80	1.22(−2)	1.67(−7)	3.55(−2)	8.16(−7)	1.01(−2)	1.21(−7)
90	1.47(−2)	2.94(−7)	2.96(−2)	6.64(−7)	7.64(−3)	7.64(−8)
100	1.60(−2)	3.73(−7)	2.54(−2)	5.41(−7)	6.48(−3)	5.47(−7)
Sn Scheme
10	1.81(−1)	1.44(−2)	5.97(−1)	9.19(−3)	3.63(−1)	4.76(−2)
20	6.82(−2)	6.41(−3)	9.22(−2)	1.15(−3)	8.92(−2)	1.70(−2)
30	3.95(−2)	3.96(−3)	4.59(−2)	4.59(−4)	4.08(−2)	1.01(−2)
40	2.71(−2)	2.58(−3)	4.13(−2)	2.82(−4)	2.77(−2)	7.83(−3)
50	2.18(−2)	1.78(−3)	4.19(−2)	1.89(−4)	1.92(−2)	5.84(−3)
60	1.82(−2)	1.39(−3)	3.59(−2)	1.34(−4)	1.28(−2)	4.27(−3)
70	1.47(−2)	1.02(−3)	2.31(−2)	8.70(−5)	1.04(−2)	3.76(−3)
80	1.11(−2)	6.94(−4)	2.30(−2)	6.22(−5)	8.49(−3)	3.28(−3)
90	1.01(−2)	5.25(−4)	1.84(−2)	4.28(−5)	7.13(−3)	2.83(−3)
100	9.67(−3)	3.93(−4)	1.59(−2)	3.02(−5)	6.21(−3)	2.46(−3)

Notes. Upper half: listing of the volume-integrated rates for the $\nu _e{\bar \nu }_e$ and "$\nu _\mu {\bar \nu }_\mu$" annihilation rates using the approach of Ruffert et al. (1996, 1997), as summarized in Section 4.1. The time for each calculation refers to the time since the start of the VULCAN/2D simulation. Numbers in parenthesis correspond to powers of 10. The strong decrease of the energy-deposition rates reflects the fading of the neutrino luminosity due to the cooling of the SMNS. A density cut above 10¹¹ g cm⁻³ is applied for both emission and deposition sites. Lower half: same as the upper half, but this time using the 8 ϑ-angle S_n calculations performed through a post-processing of the MGFLD radiation-hydrodynamics snapshots computed with VULCAN/2D. The same hydrodynamics background (ρ, T, Y_e distribution) is used for the leakage and the S_n-schemes.

Download table as:  ASCIITypeset image

Neutrino processes of emission and absorption/scattering in BNS merger simulations have been treated using leakage schemes by Ruffert et al. (1996, 1997), Rosswog & Liebendörfer (2003), Setiawan et al. (2004, 2006), and Birkl et al. (2007). In this approach, the difficult solution of the multiangle, multigroup Boltzmann transport equation is avoided by introducing a neutrino-loss timescale, for both energy and number, and applied to optically thick and optically thin conditions. In addition, an estimated rate of electron-type neutrino loss allows the electron fraction to be updated. Overall, benchmarking of these leakage schemes using more sophisticated one-dimensional core-collapse simulations that solve the Boltzmann equation suggests that the resulting neutrino luminosities are within a factor of a few (at most) of what would be obtained using a more accurate treatment.

The work presented here offers a direct test of this. Moreover, because the neutrino source terms are included in the momentum and energy equations solved by VULCAN/2D, we can model the birth and subsequent evolution of the resulting neutrino-driven wind. Our simulations being two-dimensional and less CPU intensive than three-dimensional SPH simulations can also be extended to ≳ 100 ms, and, thus, we can investigate the long-term evolution of the merger. In particular, we directly address the evolution of the neutrino luminosity and confront this with the results using the expedient of a steady state often employed in work that is also limited to a short time span of 10–20 ms after the merger.

In the leakage schemes of Ruffert et al. (1996, 1997) and Rosswog & Liebendörfer (2003), the neutrino emission processes treated are electron and positron capture on protons (yielding ν_e) and neutrons (yielding ${{\bar \nu }}_e$), respectively, and electron–positron pair annihilation and plasmon decay (each yielding electron-, μ-, and τ-type neutrinos). For electron-type neutrinos, the dominant emission processes are the charged-current β-processes. For neutrino opacity, Ruffert et al. (1996, 1997) include neutrino scattering on nucleons. Rosswog & Liebendörfer (2003) also treat electron-type neutrino absorption on nucleons.

The neutrino emission, absorption, and scattering processes used in VULCAN/2D simulations are those summarized in Burrows et al. (2006b), and include all the above, plus electron neutrino absorption on nuclei and nucleon–nucleon bremsstrahlung. The latter is the dominant emission process of ν_μ and ν_τ neutrinos in PNSs (Thompson et al. 2000). In VULCAN/2D, the ν_μ and ν_τ neutrinos are grouped together and referred to as "ν_μ" neutrinos. We present in Figure 7 the evolution of the neutrino luminosity for the BNS merger models with initially no spins (top), co-rotating spins (center), and counter-rotating spins (bottom) over a typical timespan of ≳ 100 ms, and for an adopted radius of 200 km along the equatorial (polar) directions in black (red). In each case, we plot the ν_e (solid line), the ${{\bar \nu }}_e$ (dashed line), and the "ν_μ" (dash-dotted line) neutrinos. For comparison with the predictions of the leakage scheme of Ruffert et al. (1996), we have implemented their formalism, following exactly the description given in their paper in their Appendices A and B. We plot the corresponding results at 10 ms intervals for each case (diamonds: ν_e; triangles: ${\bar \nu }_e$; squares: "ν_μ"; see also Section 4.1).

Zoom In Zoom Out Reset image size

For the BNS merger models with initially no spins, co-rotating spins, and counter-rotating spins, the neutrino signal predicted by VULCAN/2D follows a qualitatively similar evolution. The initial fast rise represents the time it takes the neutrino emission processes to ramp up to their equilibrium rates, modulated by the diffusion time out of the opaque core (and out of the not-so-opaque surface layers) and the light-travel time of ∼ 1 ms to the radius of 200 km where the luminosity is recorded. The peak luminosity occurs at about 5 ms after the start of the VULCAN/2D simulations, and is dominated by ${{\bar \nu }}_e$ and "ν_μ" neutrinos. The ν_e neutrino luminosity is initially systematically subdominant compared with the other two neutrino species. This is primarily a result of the high neutron richness of the merged object. Among the models, the stronger neutrino signal, in particular for the "ν_μ" neutrinos, is produced in the counter-rotating model, which achieves the largest central temperature (∼ 17 MeV), followed by the no-spin model (∼ 15 MeV). The co-rotating spin model has the weakest neutrino signal of all three, mainly because the tidal locking leads to a very smooth merger and correspondingly low temperatures (the central temperature is about 10 MeV). Specifically, the angle-integrated, species-integrated, peak neutrino luminosity for the BNS model with no initial spin is 1.5 × 10⁵³ (2.2 × 10⁵², 7.0 × 10⁵², and 6.7 × 10⁵²), for the co-rotating spin model 1.2 × 10⁵³ (2.2 × 10⁵², 5.7 × 10⁵², and 3.9 × 10⁵²), and for the counter-rotating spin model 2.2 × 10⁵³ erg s⁻¹(2.9 × 10⁵², 8.2 × 10⁵², and 1.2 × 10⁵³). The values given in parentheses correspond in each case to the angle-integrated peak luminosity for the ν_e, the ${\bar \nu }_e$, and the "ν_μ" neutrinos, respectively. Apart from a significant discrepancy with the "ν_μ"-neutrino luminosities, the peak neutrino luminosities using our MGFLD approach compare well with those published in the literature based on a leakage scheme (Ruffert et al. 1997; Rosswog & Liebendörfer 2003). However, in all cases, cooling of the BNS merger causes a decrease of all neutrino luminosities after the peak (i.e., rather than a long-term plateau), with a faster decline rate in the no-spin and counter-rotating spin cases. The luminosity decay rate is faster for the "ν_μ" neutrinos than for the electron-type neutrinos. This translates into a strong decrease of the neutrino–antineutrino annihilation rate due to its scaling with the square of the neutrino luminosity (see Section 4). In agreement with past work that focused on the initial 10–20 ms after the onset of coalescence, we find that the electron antineutrino luminosity dominates over the electron neutrino luminosity, but in contrast with the work of Ruffert et al. (1996) and Rosswog & Liebendörfer (2003), our "ν_μ" neutrino luminosities are at least one order of magnitude higher. This discrepancy may in part stem from their neglect of the nucleon–nucleon bremsstrahlung opacity and emission processes. Taking out the contribution from nucleon–nucleon bremsstrahlung in our VULCAN/2D BNS no-spin model leads to a general reduction of the "ν_μ" luminosity by a factor of 2–3, which then no longer dominates. However, it is still not as low as that predicted with our implementation of the leakage scheme of Ruffert et al. (1996, 1997). Moreover, in the no-spin and co-rotating spin cases, the relationships between the fluxes of the various neutrino flavors changes significantly with time, partly because they are differently affected by neutron star cooling and changes in neutron richness. While the decrease of the luminosity predicted by VULCAN/2D makes physical sense, it contrasts with the findings of Setiawan et al. (2004, 2006) in their study of neutrino emission from torus disks around ∼ 4 M_☉ black holes. In their models, they incorporate physical viscosity in the disk, together with the associated heat release, which keeps the gas temperature high. In our models, neutrino emission is a huge energy sink, not compensated by α-disk viscosity heating.

The MGFLD treatment we use permits the calculation of the energy-dependent neutrino spectrum and its variation with angle. For the BNS merger models with initially no spins (left), co-rotating spins (middle), and counter-rotating spins (right), we show such spectra along the equatorial (black) and polar (red) directions in Figure 8, using a reference radius of 200 km. We can also explicitly compute average neutrino energies, here defined as

$$\begin{eqnarray} &&\sqrt{\big\langle \varepsilon _{\nu }^2\big\rangle } \equiv \left[ \frac{\int d\varepsilon _{\nu } \varepsilon _{\nu }^2 F_{\nu }(\varepsilon _{\nu },R)}{\int d\varepsilon _{\nu } F_{\nu }(\varepsilon _{\nu },R)} \right]^{\frac{1}{2}}\,. \end{eqnarray} \tag{ 1 }$$

At 50 ms after the start of the VULCAN/2D simulations, such average neutrino energies in the no-spin BNS model and along the pole are 13.4 (ν_e), 16.7 (${\bar \nu }_e$), and 25.2 MeV ("ν_μ"), while they are 10.7, 15.6, and 21.2 MeV, respectively, along the equatorial direction. In the co-rotating model and along the polar direction, we have 14.0 (ν_e), 17.0 (${\bar \nu }_e$), and 22.0 MeV ("ν_μ"), and along the equatorial direction, we have, respectively 9.9, 15.2, and 17.5 MeV. In the counter-rotating model and along the polar direction, we have 12.9 (ν_e), 17.1 (${\bar \nu }_e$), and 27.5 MeV ("ν_μ"), and along the equatorial direction, we have 11.5, 16.7, and 25.1 MeV, respectively. The oblateness of the compact massive BNS merger, with its lower temperature and its more gradual density fall-off along the equator, systematically places the neutrinospheres in lower temperature regions, translating into softer neutrino spectra.

Zoom In Zoom Out Reset image size

The aspherical density/temperature distributions naturally lead to a strong anisotropy of the radiation field, whose latitudinal dependence is shown in Figure 9 for the three models (we use the same left to right ordering), and at a time of 60 ms after the start of the VULCAN/2D simulations. Note the stronger neutrino flux along the polar direction (see also Rosswog & Liebendörfer 2003, their Figure 12), resulting from the larger radiating surface seen from higher latitudes, as well as the systematically larger matter temperatures at the decoupling region along the local vertical (also yielding a harder neutrino spectrum). In Figure 9, we include in red the corresponding neutrino fluxes computed with the 16 ϑ-angle S_n scheme, which shows comparable fluxes along the equator, but considerably larger ones at higher latitudes. We also show in Figure 10 the energy- and species-integrated specific neutrino energy-deposition and loss rates in the no-spin BNS model at 60 ms after the start of the simulation, for both the MGFLD and the S_n calculations. Note how much more emphasized the anisotropy is with the multiangle, S_n, solver, and how enhanced is the magnitude of the deposition along the polar direction. This is a typical property seen for fast-rotating PNSs simulated with such a multiangle solver (Ott et al. 2008).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 11, we show the energy-dependent neutrinosphere radii for the ν_e (left), ${\bar \nu }_{\rm e}$ (middle), and "ν_μ" (right) neutrinos, and for the BNS merger models with initially no spins (top row), co-rotating spins (middle row), and counter-rotating spins (bottom row). These decoupling radii manifest a strong variation with latitude, but also with energy due to the approximate ε²_ν dependence of the material opacity to neutrinos. The diffusion of neutrinos out of the opaque core, which leads to global cooling of the neutron star, is also energy dependent. The location-dependent diffusion timescale $t^{\rm diff}_{\nu _i}(r,z)$ is given by (Mihalas & Mihalas 1984; Ruffert et al. 1996)

$$\begin{eqnarray} &&t^{\rm diff}_{\nu _i}(r,z) \approx 3 \tau _{\nu _i} \Delta R / c, \end{eqnarray} \tag{ 2 }$$

where c is the speed of light, ΔR is the local density scale height, $\tau _{\nu _i}$ is the optical depth (species and energy dependent), and ν_i is the neutrino species (i.e., ν_e, ${\bar \nu }_e$, or "ν_μ"). For the ν_e neutrinos, we find an optical depth in the core that varies from a few ×10² for 10 MeV neutrinos to a few ×10⁴ for 100 MeV neutrinos. These optical depths translate into diffusion times that are on the order of 10 ms to 1 s. The quasi-Keplerian disk that extends from 20–30 to 100 km is moderately optically thick to neutrinos at the peak of the energy distribution (i.e., at 10–20 MeV). Heat can thus leak out in the vertical direction over a typical diffusion time of a few tens of milliseconds (see Figure 12). This is comparable to the accretion timescale of the disk itself (Setiawan et al. 2004, 2006), which suggests that this neutrino energy is available to power relativistic ejecta, and is not advected inward with the gas, as proposed by Di Matteo et al. (2002).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Compared to PNSs formed in the "standard" core-collapse of a massive star, SMNSs formed through coalescence have a number of contrasting properties. (1) most of the material is at nuclear density, unlike in PNSs where newly accreted material deposited at the PNS surface radiates its thermal energy as it contracts and cools. (2) most of the material is neutron rich, thereby diminishing the radiation of electron-type neutrinos in favor of antielectron-type neutrinos. In PNSs, electron-type neutrinos are an important cooling agent, but have large opacities. In SMNSs, the material is strongly deleptonized and electron-type neutrinos are nondegenerate. They suffer relatively lower opacity and can therefore diffuse out more easily. The opacity of other neutrino species is even lower and therefore more prone to energy leakage through radiation. (3) energy gain through accretion is modest since there is merely a few 0.1 M_☉ of material outside of the SMNS. In our three merger models, the internal energy budget of such shear-heated SMNSs is ≳ 10⁵³ erg, half of which is thermal energy. With a total neutrino luminosity of ∼5 × 10⁵² erg s⁻¹, the cooling timescale for these SMNS should be on the order of a second. This is significantly, but not dramatically, shorter than the 10–30 s cooling timescale of hot PNSs, i.e., the time it takes these objects to radiate a total of ∼3 × 10⁵³ erg of internal energy.

A fraction of the diffusing core energy is absorbed through charge-current reactions (with associated volume-integrated energy $\dot{E}(\mathrm{cc})$) and turned into internal energy in a "gain" layer at the surface of the SMNS, primarily on the pole-facing side and at high latitudes (see Figure 10). This neutrino energy deposition is the origin of a thermally driven wind whose properties we now describe.

A few milliseconds after the start of our simulations and onset of the burst of neutrinos from the BNS merger, a neutrino-driven wind develops. It relaxes into a quasi-steady state by the end of the simulations at 100 ms. Note that in this quasi steady-state configuration, and as shown in Figure 10, the energy gain/losses due to neutrino in the polar direction is mostly net heating, with no presence of a sizable net cooling region. In Figure 13, we show this evolution for the no-spin BNS merger model for four selected times: 10, 30, 60, and 100 ms. One can see the initial transient feature, which advects out and eventually leaves the grid. This transient phase is not caused by neutrino energy deposition, but is rather due to infall at small radii and in the polar regions of the ambient medium that we placed around the merged object, as well as due to the sound waves generated by the core during the initial quakes. One can then see the strengthening neutrino-driven wind developing in the polar funnel. The side lobes of the SMNS confine the ejecta at radii ≲ 500 km to a small angle of ∼ 20° about the poles, but this opening angle grows at larger distances to ∼ 90°. The confinement at small radii also leads to the entrainment of material from the side lobes, overloading the wind and making it choke. The radial density/velocity profile of this wind flow is thus kinked, with variations in velocity that can be as large as the average asymptotic velocity value of ∼ 30,000 km s⁻¹ (Figure 4). As shown in Figure 13, the mass loss associated with the neutrino-driven wind is on the order of a few 10⁻³ M_☉ s⁻¹ str⁻¹ at a few tens of milliseconds, but decreases to 10⁻³–10⁻⁴ M_☉ s⁻¹ str⁻¹ at 100 ms. The wind is also weaker along the poles than along the 70°–80° latitudes. The associated angle-integrated mass loss summed over 100 ms approaches 10⁻⁴M_☉ and is made up in part of high Y_e material (∼ 0.5) along the pole, but mostly of low Y_e material (∼ 0.1–0.2) along midlatitudes. Thus, "r-process material" will feed the interstellar medium through this neutrino-driven wind. This latitudinal variation of the electron fraction at large distances is controlled by the relative strength of the angle- and time-dependent ν_e and ${\bar \nu }_e$ neutrino luminosities, the expansion timescale of the "wind" parcels, and the neutron richness at their launching site (Ott et al. 2007). Along the pole and at the SMNS surface, the low-density, high neutron richness, and relatively stronger ν_e neutrino luminosities at late times in the no-spin and co-rotating spin models lead to a high asymptotic Y_e value (high proton richness). Despite the relatively high resolution employed in our simulations (∼ 300 m in the radial direction at ∼ 20 km), higher resolution would be needed to resolve this region accurately. Although we expect this trend would hold at higher resolution, it would likely yield lower asymptotic values of the electron fraction along the pole.

Zoom In Zoom Out Reset image size

Along the equatorial direction, low Y_e material (∼ 0.1–0.2) migrates outward, but its velocity is below the local escape speed and it is unclear how much will eventually escape to infinity.⁸ This is further illustrated in Figure 14 where we show the angular variation of the density and velocity at 2900 km in the no-spin BNS model, at 121 ms. The BNS merger is thus cloaked along the poles by material with a density in excess of 10⁴ g cm⁻³, while along lower latitudes even denser material from the side lobes obstructs the view from the center of the SMNS. Importantly, wind material will feed the polar regions for as long as the merger remnant remains gravitationally stable. Being so heavily baryon-loaded, the outflow can in no way be accelerated to relativistic speeds with high Lorentz factors. In this context, the powering of a short-hard GRBs is impossible before black hole formation.

Zoom In Zoom Out Reset image size

As seen in Figure 5, the rotational profile in the inner 100 km is strongly differential in the no-spin and counter-rotating spin cases, while it is quasi-uniform in the co-rotating spin case (see also Rosswog & Davies 2002, their Figure 17). The good conservation of specific angular momentum in VULCAN/2D is in part responsible for the preservation of the initial rotation profile throughout the simulation. In reality, such a differential rotation should not survive. Our two-dimensional axisymmetric setup inhibits the development of triaxial instabilities that arise at modest and large ratios of rotational and gravitational energies (Rampp et al. 1998; Centrella et al. 2001; Saijo et al. 2003; Ott et al. 2005, 2007), i.e., under the conditions that prevail here. Moreover, our good, but not excellent, spatial resolution prevents the modeling of the magneto-rotational instability, whose effect is to efficiently redistribute angular momentum (Balbus & Hawley 1991b; Pessah et al. 2006, 2007), and dissipate energy, and, in the present context, leads to mass accretion onto the SMNS. The magneto-rotational instability, operating on an rotational timescale, could lead to solid-body rotation within a few milliseconds in regions around ∼ 10 km, and within a few tens of milliseconds in regions around ∼ 100 km. Note that this is a more relevant timescale than the typical ∼ 100 s that characterizes the angular-momentum loss through magnetic dipole radiation (Rosswog & Davies 2002). Importantly, in the three BNS merger models we study, we find that the free energy of rotation (the energy difference between the given differentially rotating object and that of the equivalent solid-body rotating object with the same cumulative angular momentum) is very large. Despite differences of a factor of a few between models, it typically reaches ∼5 × 10⁵¹ erg inside the SMNS (regions with densities greater than 10¹⁴ g cm⁻³), but is on the order of 2 × 10⁵² erg in the torus disk, regions with densities between 10¹¹ and 10¹⁴ g cm⁻³. Similar conditions in the core-collapse context yield powerful, magnetically (and thermally) driven explosions (LeBlanc & Wilson 1970; Bisnovatyi-Kogan et al. 1976; Akiyama et al. 2003; Ardeljan et al. 2005; Moiseenko et al. 2006; Obergaulinger et al. 2006; Burrows et al. 2007a; Dessart et al. 2007). Rotation dramatically enhances the rate of mass ejection by increasing the density rather than the velocity of the flow, even possibly halting accretion and inhibiting the formation of a black hole (Dessart et al. 2008). In the present context, the magneto-rotational effects, which we do not include here, would considerably enhance the mass flux of the neutrino-driven wind. Importantly, the loss of differential rotational energy needed to facilitate the gravitational instability is at the same time delaying it through the enhanced mass loss it induces. Work is needed to understand the systematics of this interplay, and how much rotational energy the back hole is eventually endowed with.

Oechslin et al. (2007), using a conformally flat approximation to GR and an SPH code, find that BNS mergers of the type discussed here and modeled with the Shen EOS avoid the general-relativistic gravitational instability for many tens of milliseconds after the neutron stars first come into contact. Baumgarte et al. (2000), and more recently Morrison et al. (2004), Duez et al. (2004, 2006), and Shibata et al. (2006), using GR (and for some using a polytropic EOS), find that imposing even modest levels of differential rotation yields a significant increase by up to 50% in the maximum mass that can be supported stably, in particular pushing this value beyond that of the merger remnant mass after coalescence. Surprisingly, Baiotti et al. (2008), using a full GR treatment but with a simplified (and soft) EOS, find prompt black hole formation in such high-mass progenitors. Despite this lack of consensus, the existence of neutron stars with a gravitational mass around 2 M_☉ favors a high incompressibility of nuclear matter, such as in the Shen EOS, and suggests that SMNSs formed through BNS merger events may survive for tens of milliseconds before experiencing the general-relativistic gravitational instability. In particular, the presence of a significant amount of material (a few tenths of a solar mass), located on wide orbits in a Keplerian disk, reduces the amount of mass that resides initially in the core, i.e., prior to outward transport of angular momentum.

Our first approach to estimate the neutrino–antineutrino annihilation rate is based on a leakage scheme. We implemented the method presented by Ruffert et al. (1996, 1997), which comprises two steps. First, using the processes described in Section 3.1, the instantaneous rate of neutrino emission Q(ν_i) is computed for all grid cells. It is subsequently weighted by the direction-dependent factor that depends on the radiative-diffusion timescale $t^{\rm diff}_{\nu _i}$ and neutrino-emission timescale $t^{\rm loss}_{\nu _i}$ relevant for that cell. Since we are primarily interested in the energy deposition in the polar regions, we consider only the cylindrical-z-direction when estimating $t^{\rm diff}_{\nu _i}$. The effective emissivity Q^eff(ν_i) then has the form

$$\begin{eqnarray} &&Q^{\rm eff}(\nu _i) = \frac{Q(\nu _i)}{1 + \left(t^{\rm loss}_{\nu _i}\right)^{-1} t^{\rm diff}_{\nu _i}}, \end{eqnarray} \tag{ 4 }$$

where expressions for the various components are given explicitly in Appendices A and B of Ruffert et al. (1996).

In Figure 15, we show the effective emissivity resulting from this leakage scheme and with the neutrino processes of Ruffert et al. (1996) for the ν_e (left column), ${\bar \nu }_e$ (middle column), and "ν_μ" (right column) neutrinos, for the BNS merger models with initially no spins (top row), co-rotating spins (middle row), and counter-rotating spins (bottom row). For each process, emission is conditioned by the competing elements of optical depth, density, temperature, and electron fraction. In practice, it peaks in regions that are dense, hot, but not too optically thick, i.e., at the surface of the SMNS (see density contours in Figure 15, overplotted in black). Optical depths in excess of 100 for all neutrino energies make the inner region (the inner ∼ 15 km, where densities have nuclear values) a weak "effective" emitter, with their contribution operating on a 0.1–1 s timescale. The dominant emission associated with the "ν_μ" neutrinos originates from a considerably higher-density region than that associated with the electron-type neutrinos, i.e. 10¹²–10¹³ g cm⁻³ compared to 10⁹–10¹¹ g cm⁻³, with a corresponding "leakage" luminosity in all three merger models that is typically an order of magnitude smaller than predicted by VULCAN/2D (see Figure 7). This low "ν_μ" emissivity is likely caused by the neglect of nucleon–nucleon bremsstrahlung processes in the approach of Ruffert et al. (1996, 1997), which leads to a smaller decoupling radius for "ν_μ" neutrinos, and, therefore, an underestimate of the size of the radiating surface from which they emerge (as shown in Figure 11, right column). In practice, the decoupling radius is energy dependent, but here we stress the systematic reduction of the decoupling radius for all "ν_μ" neutrino energies because of the neglect of this extra opacity source. The importance of the bremsstrahlung process for "ν_μ" emissivity and opacity has been emphasized by Thompson et al. (2000) for "hot" PNSs. The SMNS that results from the merger is considerably heated by shocks and shear, and is thus also in a configuration where such bremsstrahlung processes are important and should be included.

Zoom In Zoom Out Reset image size

From the effective emissivity distribution computed for each neutrino flavor with the leakage scheme, Ruffert et al. (1997) then sum the contributions from all pairings between grid cells. For completeness, we briefly reproduce here the presentation of Ruffert et al. (1997). The integral to be computed for the energy-integrated $\nu _i{\bar \nu }_i$ (representing equivalently $\nu _{ e}{\bar \nu }_{\rm e}$ or "$\nu _{\mu }{\bar \nu }_{\mu }$") annihilation rate at position $\vec{r}$ is

$$\begin{eqnarray} Q^+(\nu _i{\bar \nu }_i) &=& {1\over 4}\, {\sigma _0\over c\,(m_e c^2)^2} \, \Biggl \lbrace {(C_1 + C_2)_{\nu _i{\bar \nu }_i}\over 3}\cdot \oint _{4\pi }d \Omega \,I_{\nu _i} \nonumber \\ & & \times \, \oint _{4\pi } d \Omega ^{\prime }\, I_{{\bar \nu }_i} [\langle \epsilon \rangle _{\nu _i} + \langle \epsilon \rangle _{{\bar \nu }_i} ](1 - \cos \Phi)^2\nonumber\\ && +\, C_{3,\nu _i{\bar \nu }_i}(m_e c^2)^2\cdot \oint _{4\pi }d \Omega \,I_{\nu _i} \nonumber\\ &&\times\,\oint _{4\pi }d \Omega ^{\prime }\,I_{{\bar \nu }_i}\, {\langle \epsilon \rangle _{\nu _i} + \langle \epsilon \rangle _{{\bar \nu }_i}\over \langle \epsilon \rangle _{\nu _i}\,\langle \epsilon \rangle _{{\bar \nu }_i}}\, \left(1 - \cos \Phi \right)\! \Biggr \rbrace. \end{eqnarray} \tag{ 5 }$$

$I_{\nu _i}$ is the neutrino-specific intensity. Ω and Ω' are the solid angles subtended by the cells producing the neutrino and antineutrino radiation incident from all directions. C₁, C₂, and C₃ are related to the weak coupling constants, C_A and C_V, and depend on the neutrino species (see Ruffert et al. 1997 as well as Section 4.2). σ₀ is the baseline weak interaction cross section, 1.705 ×10⁻⁴⁴ cm², c is the speed of light, m_e is the electron mass. Φ is the angle between the neutrino and antineutrino beams, entering the annihilation rate formulation through the term (1 − cos Φ) (squared or not), which thus gives a stronger weighting to larger-angle collisions. This is what makes the dumbell-like morphology of BNS mergers such a prime candidate for $\nu _i{\bar \nu }_i$ annihilation over spherical configurations (see also the Appendix). $\langle \epsilon \rangle _{\nu _i}$ and $\langle \epsilon \rangle _{{\bar \nu }_i}$ are the ν_i and ${\bar \nu }_i$ mean neutrino energies, respectively, whose values we adopt for consistency from the simulations of Ruffert et al. (1997), i.e., 12, 20, and 27 MeV for ν_e, ${\bar \nu }_{\rm e}$, and "ν_μ", respectively (our values are within a few times 10% of these, so this has little impact on our discussion). Note also that the average neutrino energies for all three species are much larger than m_ec², and therefore make the second term in the above equation negligible (it is about a factor of 1000 smaller than the first one, and also has a much weaker large-angle weighting). The total annihilation rate is then the sum $Q^+(\nu _{ e}{\bar \nu }_{\rm e}) + Q^+($"$\nu _{\mu }{\bar \nu }_\mu$").

In practice, we turn the angle integrals into discretized sums through the transformation

$$\begin{eqnarray} &&\oint _{4\pi } d \Omega \,I_{\nu } \longrightarrow \sum _k \Delta \Omega _k\cdot I_{\nu,k}, \end{eqnarray} \tag{ 6 }$$

where ΔΩ_k is the solid angle subtended by the cell k as seen from the location $\vec{r}$. Ruffert et al. (1997) applied their method to three-dimensional Cartesian simulations, while our simulations are two-dimensional and have both axial symmetry, and mirror symmetry about the equatorial plane. Making use of these symmetry properties, we remap our VULCAN/2D simulation from a two-dimensional (cylindrical) meridional slice onto a three-dimensional spherical volume which covers the space with a uniform grid of 40 zones for the 2π azimuthal direction and 20 zones for the π polar direction, while the reduced radial grid has a constant spacing in the log and uses 20 zones between 12 and 120 km. We compute the annihilation rate at locations in a two-dimensional meridional slice of this new volume, with 16 uniformly spaced angles from the pole to the equator, and 40 zones with a constant spacing in the log between 15 and 120 km.¹⁰ To estimate the annihilation rate at $\vec{r}$, one needs to estimate the flux received from all cells $\vec{r}_k$. Ruffert et al. (1997) assume that neutrino radiation is isotropic in the half-space around the outward direction given by the density gradient $\vec{n}_\rho$ at $\vec{k}$. Moreover, they approximate each emitting cell volume as a sphere, whose radius at (r_k, θ_k, ϕ_k) is

$$\begin{eqnarray} &&D_k = \left(\frac{\Delta \phi _k \Delta \mu _k \Delta r_k^3}{4 \pi } \right)^{1/3}, \end{eqnarray} \tag{ 7 }$$

where μ_k = cos θ_k. With our location dependent cell volume, we obtain

$$\begin{eqnarray} &&\sum _k \Delta \Omega _k \cdot I_{\nu,k} \approx \frac{1}{\pi } \sum _k \frac{\Delta \mu _k \Delta \phi _k \Delta r_k^3}{3 d_k^2}Q^{\rm eff}(\nu _i), \end{eqnarray} \tag{ 8 }$$

where $d_k = |\vec{r} - \vec{r}_k|$. When we sum over discrete cells, the cosine of the angle between $\vec{k}$ and $\vec{k^\prime }$ is given by $\cos \Phi _{kk^\prime } = (\vec{r} - \vec{r}_k)(\vec{r}-\vec{r}_{k^\prime }) / (|\vec{r} - \vec{r}_k||\vec{r}-\vec{r}_{k^\prime }|)$.

Following Ruffert et al. (1997), we apply selection criteria to determine whether to include a contribution. Specifically, emission and deposition sites must be in regions with a density lower than 10¹¹ g cm⁻³. Interestingly, the leakage scheme predicts a very small "effective" emissivity from high-density regions, due to the very long diffusion times from those optically thick regions. Relaxing the density cuts in the calculations leads only to a 50% enhancement in volume-integrated annihilation rate $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$, but an increase of a factor of ∼ 20 for $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$"). As pointed out earlier, the location in very high-density regions of the "ν_μ" emitting cells means that most of the emitting volume is truncated by the adopted 10¹¹ g cm⁻³ density cut, while results converge when this cut is increased to densities of ∼ 10¹³ g cm⁻³. Even with the latter, $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") is still three orders of magnitude smaller than $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$, likely a result of the neglect of bremsstrahlung, and the unfavorable $(1-\cos \Phi _{k{k^\prime }})^2$ for this compact emitting configuration.

By contrast with the rather large uncertainty introduced through the somewhat arbitrary density cuts, the annihilation rate calculation is weakly affected by increasing the resolution. In the no-spin BNS model at 50 ms, changing the number of radial-angle zones for the deposition sites {nr_edep, nt_edep} from (40, 16) to (60, 48) and for emitting sites {nr_emit, nt_emit, np_emit} from (20, 20, 40) to (30, 30, 60) increases $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$ from 2.82 × 10⁴⁹ to 3.14 × 10⁴⁹ erg s⁻¹ (a ∼ 10% increase) and increases $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") from 6.23 × 10⁴⁴ to 7.59 × 10⁴⁴ erg s⁻¹ (a ∼ 20% increase). The same test done for the no-spin BNS model at 60 ms after the start of the simulations yields a $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$ which is identical to within 1%, while $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") increases by ∼ 10% in the higher-resolution model. Hence, higher resolution does not change these values by a significant amount.

In the left half of each panel in Figure 16, we show the distribution in a two-dimensional meridional slice (it is in fact axisymmetric by construction) for the annihilation rate $Q^+(\nu _{ e}{\bar \nu }_{\rm e})$ of electron-type neutrinos, for the BNS merger models with initially no spins (left), co-rotating spins (middle), and counter-rotating spins (right). Volume-integral values $\dot{E}(\nu _{\rm i}{\bar \nu }_{\rm i})$ as a function of time for all three models are shown in Figure 17 (left panel; dashed line joined by star symbols) and given in Table 1. The deposition rate is maximum near the poles, where the large-angle collisions occur, and close to the peak emission sites (as shown in Figure 15), i.e., near the SMNS surface, with peak values on the order of 10²⁹ erg cm⁻³ s⁻¹ for all three models. In this formalism, $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") (the middle panel of Figure 17) is at least three orders of magnitude smaller than $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$. The peak volume-integrated energy deposition reaches a few 10⁵⁰ erg s⁻¹, typically two orders of magnitude below that achieved by charge-current reactions, whose associated energy deposition drives the neutrino-driven wind (Section 3; the right panel of Figure 17).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now present a formalism for the computation of the $\nu _i{\bar \nu }_i$ annihilation rate that directly exploits the neutrino-transport solution computed with VULCAN/2D in our BNS mergers, rather than using evaluations based on the leakage scheme (Ruffert et al. 1996, 1997).

Our approach is to use the angle-dependent neutrino-specific intensity calculated for snapshots at 10 ms intervals for these BNS merger models using our multiangle, S_n, scheme (Livne et al. 2004; Ott et al. 2008; Section 2.2). In the Appendix, and for completeness, we present such annihilation-rate calculations, but in the context of the postbounce phase of a 20 M_☉ progenitor (Ott et al. 2008) and of the accretion-induced collapse (AIC) of a massive and fast-rotating white dwarf (Dessart et al. 2006b). The formalism presented here applies to all cases equivalently.

Following Burrows et al. (2006b) and Janka (1991), the expression for the local energ-deposition rate $Q^+({\nu _i{{\bar \nu }}_i})$ at a position $\vec{r}$ due to annihilation of $\nu _i{{\bar \nu }}_i$ pairs into electron–positron pairs is given to leading order by

$$\begin{eqnarray} Q^+({\nu _i{{\bar \nu }}_i}) &=& \frac{1}{6} \frac{\sigma _0 \big(C_A^2 + C_V^2\big)_ {\nu _i {{\bar \nu }}_i}}{c(m_e c^2)^2} \int _0^{\infty } d\varepsilon _\nu \int _0^{\infty } d\varepsilon _\nu ^{\prime } (\varepsilon _\nu + \varepsilon _\nu ^{\prime }) \nonumber \\ &&\times\,\oint _{4\pi } d\Omega \oint _{4\pi } d\Omega ^{\prime } I_{\nu _i} I^{\prime }_{{{\bar \nu }}_i} (1 - \cos \Phi)^2, \end{eqnarray} \tag{ 9 }$$

where $I_{\nu _i} \equiv I_{\nu _i}(\varepsilon _\nu,\vec{r},\vec{n},t)$ is the ν_i-neutrino-specific intensity at energy ε_ν, location $\vec{r}$, along the direction $\vec{n}$, and at time t. The primes denote the antiparticle. The angle Φ is the angle between the directions $\vec{n}$ and $\vec{n}^{\prime }$ of the neutrino and antineutrino, i.e., $\cos \Phi = \vec{n}\cdot \vec{n}^{\prime }$.

The formulation we present applies equally to all neutrino species, as long as the values of the weak coupling constants C_A and C_V are appropriately set. We have C_V = 1/2 + 2sin²θ_W for the electron types, C_V = −1/2 + 2sin²θ_W for the ν_μ and ν_τ types, and C_A = ±1/2, with sin²θ_W = 0.23. The other variables have their usual meanings. Note that this formulation neglects the phase-space blocking by the final-state electron–positron pair which is relevant only at high densities and negligible in the semi-optically thin regime where pair annihilation may contribute to the net heating.

Now, setting the flux factor $\vec{h} = \vec{H}/J$ and the Eddington-tensor factor ${{\mathsf k}}={{\mathsf K}}/J$, and by expanding the $\cos \Phi = \vec{n}\cdot \vec{n}^{\prime }$ term with an appropriate choice of the radiation unit vector (Hubeny & Burrows 2007), we obtain for the general case in three dimensions

$$\begin{eqnarray} Q^+({\nu _i{{\bar \nu }}_i}) &=& \frac{8\sigma _0\pi ^2 \big(C_A^2 + C_V^2\big)_{\nu _i{{\bar \nu }}_i}}{3c(m_e c^2)^2} \int _0^{\infty } d\varepsilon _\nu \int _0^{\infty } d\varepsilon _\nu ^{\prime } (\varepsilon _\nu + \varepsilon _\nu ^{\prime }) \nonumber \\ &&\times\,J_{\nu _i} J^{\prime }_{{{\bar \nu }}_i} \big(1-2\vec{h}_{\nu _i}\cdot \vec{h}^{\prime }_{{{\bar \nu }}_i} + \mathrm{Tr}[{{\mathsf k}}_{\nu _i}:{{\mathsf k}}^{\prime }_{{{\bar \nu }}_i}]\big), \end{eqnarray} \tag{ 10 }$$

where the term $\mathrm{Tr}[{{\mathsf k}}_{\nu _i}:{{\mathsf k}}^{\prime }_{{{\bar \nu }}_i}]$ is the trace of the matrix product of the two J_ν-normalized Eddington tensors.

Assuming the special choice of cylindrical coordinates in axisymmetry and neglecting the velocity dependence of the radiation field, we obtain

$$\begin{eqnarray} Q^+({\nu _i{{\bar \nu }}_i}) &=& \frac{8\sigma _0\pi ^2 \big(C_A^2 + C_V^2\big)_{\nu _i{{\bar \nu }}_i}}{3c(m_e c^2)^2} \int _0^{\infty } d\varepsilon _\nu \int _0^{\infty } d\varepsilon _\nu ^{\prime } (\varepsilon _\nu + \varepsilon _\nu ^{\prime }) \nonumber \\ &&\times\,JJ^{\prime } \times (1 - 2 h_r h^{\prime }_r - 2 h_z h^{\prime }_z \nonumber \\ &&+ {{\mathsf k}}_{rr}{{\mathsf k}}^{\prime }_{rr} + {{\mathsf k}}_{zz}{{\mathsf k}}^{\prime }_{zz} + {{\mathsf k}}_{\phi \phi }{{\mathsf k}}^{\prime }_{\phi \phi } + 2 {{\mathsf k}}_{rz}{{\mathsf k}}^{\prime }_{rz}), \end{eqnarray} \tag{ 11 }$$

where we have suppressed the ν_i and ${{\bar \nu }}_i$ subscripts. Note that k_rz is the only nonzero off-diagonal term in these coordinates. Using the trace condition on the Eddington tensor, ${{\mathsf k}}_{rr} + {{\mathsf k}}_{zz} + {{\mathsf k}}_{\phi \phi } = 1$, the above can be recast into

$$\begin{eqnarray} Q^+({\nu _i{{\bar \nu }}_i})&=& \frac{8\sigma _0\pi ^2 \big(C_A^2 + C_V^2\big)_{\nu _i{{\bar \nu }}_i}}{3c(m_e c^2)^2} \int _0^{\infty } d\varepsilon _\nu \int _0^{\infty } d\varepsilon _\nu ^{\prime } (\varepsilon _\nu + \varepsilon _\nu ^{\prime }) \nonumber \\ &&\times\,JJ^{\prime } \times (1 - 2 h_r h^{\prime }_r - 2 h_z h^{\prime }_z + {{\mathsf k}}_{rr}{{\mathsf k}}^{\prime }_{rr} + {{\mathsf k}}_{zz}{{\mathsf k}}^{\prime }_{zz}\nonumber \\ &&+ (1-{{\mathsf k}}_{rr}-{{\mathsf k}}_{zz})(1-{{\mathsf k}}^{\prime }_{rr}-{{\mathsf k}}^{\prime }_{zz}) + 2 {{\mathsf k}}_{rz}{{\mathsf k}}^{\prime }_{rz}).\nonumber\\ \end{eqnarray} \tag{ 12 }$$

For completeness and comparison with Janka (1991), we transform to spherical coordinates and reduce to spherically symmetric problems, which then yields

$$\begin{eqnarray} Q^+({\nu _i{{\bar \nu }}_i}) &=& \frac{8\sigma _0\pi ^2}{3c(m_e c^2)^2} \big(C_A^2 + C_V^2\big)_{\nu _i{{\bar \nu }}_i}\! \int _0^{\infty }\! d\varepsilon _\nu \int _0^{\infty }\! d\varepsilon _\nu ^{\prime } (\varepsilon _\nu + \varepsilon _\nu ^{\prime }) \nonumber \\ &&\times\,JJ^{\prime } \times \bigg(1 - 2 h h^{\prime } + k k^{\prime } + \frac{1}{2}(1-k)(1-k^{\prime })\bigg),\nonumber\\ \end{eqnarray} \tag{ 13 }$$

with scalar h and k. Equation (13) corresponds to the combination of Equations (1) and (6) of Janka (1991).

In practice, our approach is to start from snapshots of the (converged) MGFLD radiation-hydrodynamics simulations of the three BNS mergers presented in Section 3. Keeping the hydrodynamical variables frozen, the radiation variables are relaxed using the S_n algorithm and with 8 ϑ-angles. As an accuracy check, we also compare the S_n results based on simulations performed with 16 ϑ-angles. When we do this, the converged radiation field for the 8 ϑ-angle solution is remapped into the 16 ϑ-angle run (corresponding to a total of 144 angles), which is then relaxed. These simulations are quite costly, and take about 4 days with 48 processors. However, by contrast with the approach of Ruffert et al. (1996, 1997), the annihilation rates are then simply obtained from direct (local) integration of the moments of the neutrino-specific intensity.

In the right half of each panel in Figure 16, we show the distribution in a two-dimensional meridional slice for the annihilation rate $Q^+(\nu _{ e}{\bar \nu }_{\rm e})$ of electron-type neutrinos, for the BNS merger models with initially no spins (left), co-rotating spins (middle), and counter-rotating spins (right), but now computed from the results of the S_n calculation and the formalism presented above. We provide in Figure 17 volume-integral values as a function of time for all three models (the top and middle panels; solid lines) and given in Table 1. By contrast with the leakage scheme results, the deposition is more evenly spread around the SMNS, is maximum at depth (in the region where the density cut applies), and reaches peak values a factor of a few larger. We find that it is not only the large-angle collisions that favor the deposition, but also the dependence on J² (and its 1/R⁴ radial dependence in optically thin regions) which considerably weights the regions close to the radiating SMNS, making the deposition large not only in the polar regions, but also all around the SMNS. After 40 ms, and although the associated neutrino luminosities are quite comparable between the three models, the (somewhat unrealistic) co-rotating spin BNS model boasts the largest annihilation rate. This stems from the very extended high-density SMNS configuration, leading to larger neutrinosphere radii and extended regions favoring large-angle collisions. This contrasts with the fact that it is the least hot of all three configurations, and is a somewhat weaker emitter, suggesting that it is not just the neutrino luminosity, but also the spatial configuration of the SMNS that sets the magnitude of the annihilation rate. In the new formalism, which employs results from the S_n simulations, $\dot{E}$("$\nu _{\mu }{\bar \nu }_{\mu }$") (the middle panel of Figure 17) is now only one order of magnitude smaller than $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$, even for the same adopted density cut. By contrast with the leakage scheme, and to some extent because of the treatment of bremsstrahlung processes, the emission from "ν_μ" neutrinos is associated with more exterior regions of the SMNS, the ∼ 25 MeV "ν_μ" neutrinos decoupling at ∼ 18 and ∼ 100 km along the polar and equatorial directions, respectively, and thus in regions with densities on the order of, or lower than, 10¹¹ g cm⁻³ (Figure 11). The emission is, thus, not significantly truncated by the adopted density cut. Moreover, the subtended angle of the representative "ν_μ" neutrinosphere is correspondingly much bigger, favoring larger-angle collisions.

This disagreement should not overshadow the good match between the leakage- and the S_n-scheme predictions for the $\nu _{ e}{\bar \nu }_{\rm e}$ annihilation rate, which are typically within a few tens of percent, only sometimes in disagreement by a factor of ∼ 3 (see, e.g., the no-spin BNS model at 10 ms). The S_n-scheme results, thus, support the long-term decline of the annihilation rate predicted with the leakage scheme, but with a $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") value that is typically a factor of 10 smaller than $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$ at all times. The results given here using 8 ϑ-angles are already fairly converged. For the co-rotating spin BNS model at 60 ms after the start of the simulation, we obtain the following differences between the 16 and 8 ϑ-angle S_n calculations: $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$ changes to 3.54 × 10⁴⁹ from 3.59 × 10⁴⁹ erg s⁻¹ (down by 2%), $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") changes to 1.73 × 10⁴⁷ from 1.34 × 10⁴⁷ erg s⁻¹ (up by 29%), and $\dot{E}(\mathrm{cc})$ changes to 1.40 × 10⁵¹ from 1.64 × 10⁵¹ erg s⁻¹ (down by 15%).

The long time-coverage of our simulations, up to ≳ 100 ms, shows that the peak value, which is also coincident with the peak neutrino luminosity at 5–10 ms, is not sustained. The conditions at the peak of the annihilation rate do not correspond to a steady state, but instead herald the steady decrease that follows as the BNS merger radiates, contracts, and cools. In all simulations, $\dot{E}(\nu _{ e}{\bar \nu }_{\rm e})$ is down by a factor 10–100 at 100 ms compared to its peak value, while for $\dot{E}($"$\nu _{\mu }{\bar \nu }_{\mu }$") the decrease is by 2–3 orders of magnitude (this component is in any case subdominant). The J² dependence, combined with the strong fading of all neutrino emissivities, is at the origin of the steady decrease of the annihilation rate over the time span considered here. In their black hole/torus-disk configuration with α-disk viscosity (which causes significant shear heating not accounted for here), Setiawan et al. (2004, 2006) observe that the annihilation rate decreases as t^−3/2, and, thus, only slightly more gradually than our prediction. In our work, neutrino energy deposition by charge-current reactions is the dominant means to counteract the global cooling effect of neutrino emission.

We now use the results for the S_n scheme to discuss the assumption of Ruffert et al. (1996, 1997) that neutrino emission is quasi-isotropic. More precisely, they propose that the effective emission for every cell is isotropic in the half space around the outward direction given by the density gradient $\vec{n}_\rho$. By contrast, in the S_n scheme, the angular distribution of the neutrino radiation field is solved for, and we can determine how isotropic this emission is. For that purpose, we study the spatial variation of the angular term in Equation (11),

$$\begin{eqnarray} W &=& (1 - 2 h_r h^\prime _r - 2 h_z h^\prime _z \nonumber \\ &&+\, k_{rr}k^\prime _{rr} +k_{zz} k^\prime _{zz} + k_{\varphi \varphi } k^\prime _{\varphi \varphi } +2 k_{rz} k^\prime _{rz}), \end{eqnarray} \tag{ 14 }$$

and plot it in Figure 18 for the ν_e–${\bar \nu }_e$ annihilation rate at a representative neutrino energy of 12.02 MeV. Overplotting the (representative) 12.02 MeV ν_e flux vectors, together with isodensity contours, one sees that the flux is everywhere oriented predominantly in the radial direction, and peaks along the polar direction. The outward direction given by the local density gradient (given by the perpendicular-to isodensity contours, shown here as short white tick marks) is collinear to the flux vector along near-polar latitudes, while at midlatitudes they are perpendicular to each other. Hence, the assumption by Ruffert et al. (1996, 1997) of isotropic neutrino emission in the half space around the outward direction given by the local density gradient is not accurate in this instance. Note that it will likely hold more suitably once the black hole forms and only the torus disk radiates.

Zoom In Zoom Out Reset image size

In Figure 18, W has a maximum of ∼ 1.33 in the optically thick regions of the extended (equatorial) disk-like structure. There, the flux-factor terms and k'_φφ + 2k_rzk'_rz are nearly zero, while k_rrk'_rr + k_zzk'_zz + k_φφk'_φφ = 1/3. In polar regions, the transition to free streaming happens at small radii and W decreases from ∼ 1.33 at ≲ 20 km to ∼ 0.3 at 40 km. Outside ∼ 60–100 km, W ≈ 0, since h_zzh'_zz ≈ k_zzk'_zz ≈ 1 and all other terms are small. Considering now the fact that the distribution of the mean intensity, J, irrespective of energy group and species, is oblate inside the rapidly spinning SMNS and transitions to a prolate shape outside (see Figure 10), we can explain the spatial distribution of the pair annihilation rate (as shown in Figure 16). The rate peaks inside the oblate SMNS, since JJ' and W are largest there. Along the polar axis, W decreases rapidly with radius, while JJ' becomes prolate and remains large out to large radii, compensating in part for the small W. In equatorial regions, on the other hand, JJ' drops off more rapidly, but W remains near 1.33 out to a radius of ∼ 150 km, and only slowly transitions to zero as the equatorial radiation field gradually becomes forward-peaked. This systematic behavior sustains a significant $\nu _i{\bar \nu }_i$ annihilation rate for equatorial radii of up to ∼ 150 km. Note, however, that net energy deposition by pair annihilation does not occur in the equatorial charged-current loss regions.

The findings described in the above paragraph indicate that the dominant contribution to $Q^+(\nu _i{\bar \nu }_i)$ is not coming from large-angle collisions of neutrinos emitted from the disk, but rather is due to collisions at all angles and close to the high-density, high-temperature, SMNS surface. This is in distinction to the annihilation rate computed using the leakage scheme and the assumptions concerning the morphology of the neutrino radiation field made in Ruffert et al. (1997). However, once the SMNS transitions to a black hole, the inner contribution will vanish and the radiation will come exclusively from the cooling hot torus. In a future paper we will address the annihilation rates that we obtain with this BNS merger configuration.

In the future, while retaining the merits of the S_n method, we will need to improve the computation of the annihilation rates by accounting for GR effects. The compact and massive configuration of BNS mergers, with GM/Rc² on the order of 20% at 20 km, suggests that the numbers we present could be modified with this improvement, although Birkl et al. (2007) suggest the magnitude of the effect is at most a few tens of percent.