The influence of neutrinos on r-process nucleosynthesis in the ejecta of black hole–neutron star mergers

Abstract

During the merger of a black hole and a neutron star, baryonic mass can become unbound from the system. Because the ejected material is extremely neutron-rich, the r-process rapidly synthesizes heavy nuclides as the material expands and cools. In this work, we map general relativistic models of black hole–neutron star mergers into a Newtonian smoothed particle hydrodynamics (SPH) code and follow the evolution of the thermodynamics and morphology of the ejecta until the outflows become homologous. We investigate how the subsequent evolution depends on our mapping procedure and find that the results are robust. Using thermodynamic histories from the SPH particles, we then calculate the expected nucleosynthesis in these outflows while varying the level of neutrino irradiation coming from the post-merger accretion disc. We find that the ejected material robustly produces r-process nucleosynthesis even for unrealistically high neutrino luminosities, due to the rapid velocities of the outflow. None the less, we find that neutrinos can have an impact on the detailed pattern of the r-process nucleosynthesis. Electron neutrinos are captured by neutrons to produce protons while neutron capture is occurring. The produced protons rapidly form low-mass seed nuclei for the r-process. These low-mass seeds are eventually incorporated into the first r-process peak at A ∼ 78. We consider the mechanism of this process in detail and discuss if it can impact galactic chemical evolution of the first peak r-process nuclei.

hydrodynamics, neutrinos, nuclear reactions, nucleosynthesis, abundances, stars: neutron

1 INTRODUCTION

black hole–neutron star (BHNS) binary mergers are a likely candidate for Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) and Advanced Virgo detections of gravitational waves (Aasi et al. 2015; Acernese et al. 2015), they may be responsible for short gamma-ray bursts (sGRBs; e.g. Lee & Ramirez-Ruiz 2007), and they may provide a significant fraction of the r-process material found in our Galaxy (e.g. Lattimer & Schramm 1976; Korobkin et al. 2012; Bauswein, Stergioulas & Janka 2014a). Within the next few years, it is likely that Advanced LIGO will detect gravitational waves from these systems and constrain the BHNS merger rate. If electromagnetic counterparts are detected, the merger-short gamma-ray burst (sGRB) connection may be confirmed and production of the r-process nuclei may be observed in situ (Metzger & Berger 2012).

The origin of the r-process nuclei has been a longstanding question in nuclear astrophysics (Burbidge et al. 1957). Core-collapse supernovae are appealing as a possible site because of galactic chemical evolution considerations (e.g. Qian 2000; Argast et al. 2004), but there is significant difficulty finding the requisite conditions for r-process nucleosynthesis in this environment (e.g. Arcones & Thielemann 2013). Conversely, it is relatively easy to find conditions neutron-rich enough for r-process nucleosynthesis in the material ejected from binary neutron star (NS) and BHNS mergers (Freiburghaus, Rosswog & Thielemann 1999). Due to the long delay time from binary formation to merger and the large amount of material ejected per merger event, it is challenging to get simple models of galactic chemical evolution, which invoke compact object mergers for r-process production to agree with the observed distribution of r-process elements in low-metallicity halo stars (Qian 2000; Argast et al. 2004). Nevertheless, recent works taking into account more complex models of galaxy formation get reasonable agreement with the observed distribution of r-process elements (Matteucci et al. 2014; Ishimaru, Wanajo & Prantzos 2015; Shen et al. 2015; van de Voort et al. 2015), and it is possible to get r-process enrichment at very low metallicity when different channels of binary formation are considered (Ramirez-Ruiz et al. 2015). Therefore, it is plausible that compact object mergers could be a significant source of the galactic r-process nuclei.

Recently, it has been recognized that weak interactions can significantly affect the final composition of binary NS outflows (Wanajo et al. 2014; Foucart et al. 2016; Goriely et al. 2015; Palenzuela et al. 2015; Radice et al. 2016; Sekiguchi et al. 2015). Likewise, the final state and remnant product of binary NS mergers has been shown to depend on several properties of the system, e.g. important roles are played by the microphysical nuclear equation of state (EOS), electromagnetic fields and neutrino effects (Neilsen et al. 2014; Palenzuela et al. 2015). In contrast to binary NS mergers, the material ejected during BHNS mergers is unlikely to undergo significant numbers of weak interactions. Electron and positron captures are suppressed relative to the rates in the shock heated ejecta of binary NS mergers due to the low entropy present in the tidal ejecta. The high outflow speeds and low neutrino luminosities encountered in these events – compared to binary NS mergers – also make it unlikely that neutrino interactions will drastically change the number of neutrons present at the onset of r-process nucleosynthesis (Foucart et al. 2014; Foucart et al. 2015). Therefore, the dynamical ejecta of BHNS mergers have been long thought to be likely sites for production of heavy r-process nucleosynthesis (Lattimer & Schramm 1976; Lattimer et al. 1977; Korobkin et al. 2012; Bauswein et al. 2014b).

Understanding how BHNS mergers contribute to galactic chemical evolution requires knowledge of the merger rate, predictions for the amount of mass ejected per merger, the kinetic energy of the ejecta, and predictions of nuclei synthesized in these outflows. Although there are no observed BHNS binaries, theoretical predictions suggest that the rate of BHNS binary mergers could be up to a tenth of the rate of double NS binary mergers (Abadie et al. 2010; Bauswein et al. 2014b). The amount of mass ejected during BHNS mergers can depend sensitively on the binary parameters, especially the black hole (BH) spin and mass (Foucart et al. 2013; Hotokezaka et al. 2013; Bauswein et al. 2014b; Kyutoku et al. 2015). More mass is ejected as the BH spin increases in the direction of the orbital angular momentum (Foucart et al. 2014). Increasing the spin decreases the radius of the innermost stable orbit and decreases gravitational binding at the radius at which the NS is tidally disrupted. Increasing the BH mass reduces the amount of material remaining outside of the BH after merger (for fixed NS properties), since the tidal radius scales as (M_BH/M_NS)^1/3R_NS while the innermost stable orbit of the BH scales as M_BH for fixed BH spin. The fraction of the mass outside the horizon which is unbound, however, also increases with the BH mass, making the relation between BH mass and unbound mass nontrivial (Kyutoku et al. 2015).

Because the mass and spin distributions of stellar mass BHs and the expected number of BHNS system in our Galaxy are not well known (e.g. The LIGO Scientific Collaboration 2010), it is difficult to estimate the contribution of these events to the r-process material found in the Galaxy (Bauswein et al. 2014b). None the less, it is timely to investigate the detailed composition of the ejecta because the merger rate is likely to soon be constrained by Advanced LIGO (The LIGO Scientific Collaboration 2015). Additionally, there are some hints that the infrared excess associated with GRB130603B (Berger, Fong & Chornock 2013; Tanvir et al. 2013) is consistent with that event being powered by the radioactive decay of r-process products in the ejecta of a BHNS merger (Hotokezaka et al. 2013). A similar excess has recently been observed in the afterglow of GRB060614 (Jin et al. 2015; Yang et al. 2015).

In this work, we investigate the long-term hydrodynamics of the BHNS ejecta and the nucleosynthesis which occurs therein. For the first time, we focus on how neutrinos might affect the detailed nucleosynthesis patterns which are produced. Even for unrealistically large neutrino luminosities, we find that the distribution of the pre-neutron capture electron fraction is not significantly altered and the second and third r-process peaks are robustly produced in almost all of the material. This is in contrast to the dynamical ejecta of binary NS mergers, where weak processing may prevent an r-process from occurring in a significant amount of the material (Wanajo et al. 2014; Goriely et al. 2015). Of course, the BHNS result is expected because the outflows happen relatively early before the remnant disc can start to emit neutrinos, there is no hypermassive NS contributing to the neutrino flux, and the tidal ejecta possesses a very high velocity. More interestingly, we find that electron neutrino captures by neutrons can provide seed nuclei for a low-mass r-process which produces material in the first r-process peak at A ∼ 78. None the less, in our models, the ratio of the first peak to the second peak is subsolar with and without the inclusion of neutrino captures. When comparing to the yields of low-metallicity halo stars with subsolar Ge abundances (Roederer et al. 2014), we find that this first peak production can bring our models closer to agreement with the observed abundances of Ge, As and Se, although the abundances are still somewhat low.

This paper is organized as follows. In Section 2, we present the BHNS systems we have simulated, explain how the ejected material is mapped into our smoothed particle hydrodynamics (SPH) code and describe our nuclear reaction network. Then, in Section 3.1, we discuss the effect of weak interactions on the electron fraction distribution in the ejecta. In Section 3.2, we present the integrated nucleosynthesis from our models and discuss neutrino-induced production of the first r-process peak. In Section 3.4, we discuss uncertainties in the results from our nucleosynthesis calculations and their possible implications for galactic chemical evolution and for abundance observations in low-metallicity halo stars.

2 METHODS

2.1 Relativistic merger simulations and binary systems

The BHNS merger simulations used in this work have been described in detail in our previous papers (Deaton et al. 2013; Foucart et al. 2014). Here, we review the major features and error estimates of the merger simulations, referring readers to Foucart et al. (2014) for details. The fully relativistic Einstein-hydrodynamics system is evolved with the Spectral Einstein Code (spec) (SXS Collaboration 2000). Neutrino cooling and lepton number evolution are incorporated through a neutrino leakage scheme (Deaton et al. 2013).

To model the NS, we employ the Lattimer–Swesty EOS (Lattimer & Swesty 1991) with an incompressibility K₀ = 220 MeV and a symmetry energy S_ν = 29.3 MeV (hereafter LS220), using the table available at http://www.stellarcollapse.org and described in O'Connor & Ott (2010). This EOS yields an NS radius which lies within the allowed range of radii, as determined by Hebeler et al. (2013) from nuclear theory constraints and the existence of NSs of mass ∼2 M_⊙ (Demorest et al. 2010; Antoniadis et al. 2013). For LS220, a 1.2 (1.4) M_⊙ NS has a radius R_NS of 12.8 (12.7) km and a compactness C = GM_NS/(R_NSc²) of 0.139 (0.163).

During the spec simulations, the dynamical ejecta is tracked for only about 5 ms before it exits the computational grid. However, during this time, the specific energy (u_t) of fluid elements becomes nearly constant, so it is often possible to confidently identify unbound material. Convergence of our spec simulations was observed to be faster than the second order. Assuming second-order convergence gives a conservative relative error of up to 60 per cent in the mass and kinetic energy of ejected material. Even if the true error were this large, which is unlikely, it would not affect the results of the present investigation. As we will see, variations of ejecta properties between different binary systems, which are of similar magnitude, have negligible effect on the final nuclear abundances.

In the simulations of Deaton et al. (2013) and Foucart et al. (2014), we considered BHNS binary systems with multiple masses and spins. The BH mass M_BH was taken to be 5.6, 7 or 10 M_⊙, covering most of the estimated mass distribution for stellar mass BHs (Özel et al. 2010; Farr et al. 2011). The NS gravitational mass M_NS was taken to be 1.2 or 1.4 M_⊙, which is typical for NSs (Kiziltan et al. 2013). For these masses, ejecta is produced only for at least moderately high BH spins, meaning that for most cases the Kerr spin parameter must be χ_BH > 0.7 (Foucart 2012).

For this study, we use the ejecta from three systems. The first, called ‘M12-7-S9’, with parameters M_NS = 1.2 M_⊙, M_BH = 7 M_⊙, χ_BH = 0.9, produces a very large ejecta mass of 0.16 M_⊙. The second, ‘M14-7-S8’, with M_NS = 1.4 M_⊙, M_BH = 7 M_⊙, χ_BH = 0.8, has ejecta mass 0.06 M_⊙, one of our lower ejecta mass cases. The third case, ‘M14-5-S9’, has parameters M_NS = 1.4 M_⊙, M_BH = 5.6 M_⊙, χ_BH = 0.9 and ejects a mass of 0.084 M_⊙.

2.2 SPH evolution of ejecta

After ∼5 ms, the ejecta has detached from the merger remnant and is moving for the most part ballistically. However, the outflow is not yet homologous. Also, it is possible that pressure forces will subsequently become important again because of recombination heating or collision of streams of matter (although this turns out not to be the case). Therefore, we continue the hydrodynamic evolution of the outflow using an SPH code, starsmasher (Gaburov, Lombardi & Zwart 2010; Ponce, Faber & Lombardi 2012). The SPH code is Newtonian, but since the flow is only mildly relativistic (v/c ≈ 0.2), and from the beginning somewhat far from the BH (>10M_BH), this is probably adequate for our purposes. (See check on this in the following.)

The starsmasher code is the successor to the earlier starcrash code (Lombardi et al. 2006). It represents fluids in the standard SPH way, using a finite number of fluid elements or ‘particles’. In its current implementation, the particles may have different masses (Gaburov et al. 2010), which simplifies the construction of initial data from finite volume representations. starsmasher uses variable smoothing lengths to maximize resolution, using a formalism derived consistently from a particle-based Lagrangian to ensure proper energy and entropy evolution (Monaghan 2002; Lombardi et al. 2006; Springel & Hernquist 2002).

Stable shock evolution is achieved using artificial viscosity with a Balsara switch (Balsara 1995) to suppress artificial viscosity in shear layers; fortunately, accurate shock evolution is not important for our application. Self-gravity forces are neglected, so the gravitational force is simply a function of position given by the BH potential and it is implemented in the Newtonian and Paczyński–Wiita approximations (Ponce et al. 2012; Paczynski & Wiita 1980). In order to avoid small time steps due to rapid motion, particles are removed if they come too close to the BH. These particles would eventually fall into the BH anyway, so this procedure does not affect the ejecta properties.

As initial data to the SPH simulation, hydrodynamic data from a snapshot of the spec merger simulation (taken after tidal disruption but before the tidal tail hits the outer boundary) are output on a uniform Cartesian mesh. starsmasher reads these data, reflects them to add the lower hemisphere not evolved by spec, interpolates to an hexagonal close-packed lattice and assigns a particle of appropriate density to each non-vacuum lattice point. The evolution is then continued in starsmasher using the LS220 EOS with no neutrino effects. The electron fraction Y_e of each particle is constant during the SPH evolution, and no neutrino cooling or absorption is considered. If a particle falls below the LS220 density or temperature table range, the entropy S is henceforth set to be constant, and an S ∝ ρT³, P ∝ ρ^4/3 extrapolation of the EOS is used. This only happens when pressure is negligible, and the entropy evolution in starsmasher is not used in our post-processing nucleosynthesis calculations (see the following).

Although relativistic effects are not expected to be important, the translation from relativistic to Newtonian physics must account for two subtleties. First, the late-time behaviour of an ejecta fluid element is most sensitive to its energy, especially whether it is bound or unbound, so it is important that this be appropriately translated. We therefore rescale the velocity vector so that the specific kinetic plus potential energy of each particle in the Newtonian framework is equal to its relativistic specific energy −u_t − 1 in the spec simulation. Secondly, there is no a priori guarantee that the coordinate system in which the numerical relativity simulation evolves will be close to any known coordinates. Fortunately, spec's ‘damped harmonic’ coordinates lead the space–time to settle nearly in harmonic coordinates, so we transform in starsmasher to Schwarzschild coordinates (ignoring the BH's spin, whose effects will not be important far from the hole), with a simple radial transformation r → r + M_BH. In the region of interest, the numerically evolved space–time is nearly Minkowski, and the deviation from Minkowski is mostly Schwarzschild and so can be adequately modelled by a Paczyński–Wiita potential. Lastly, one must distinguish between the rest-frame baryon density used in the EOS and the mass integrand density used to assign the mass of the SPH particle, which is the rest-frame baryon density times a Lorentz factor and a metric determinant factor. Because SPH particle mass is a constant, the mass integrand density only needs to be calculated and integrated over at the initial time.

A straightforward evolution of the fluid equations produces a generally realistic evolution but with some clearly unphysical artefacts. Namely, matter on the upper and lower surfaces of the ejected tidal tail blow away from the equator, something unexpected given the overall weakness of pressure forces and not indicated in the spec evolution. This vertical expansion has no influence on the energy distribution or nucleosynthesis results, but it does affect the shape of the outflow. Convergence tests show that it is not a transient caused by an insufficient number of particles, so it is likely an artefact of the transition to Newtonian physics. It can be removed by reducing pressure forces near the BH. In our simulations labelled ‘P1’, we turn off pressure forces within 10 M_BH of the BH, while within 100 M_BH, pressure is reduced by a factor varying linearly with distance between 0 (at 10 M_BH) and 1 (at 100 M_BH). Within this range, the specific entropy is held constant, because otherwise the pressure reduction would keep the fluid from adiabatically cooling. Simulations labelled ‘P2’ have full pressure forces everywhere.

We check that our evolved results are insensitive to the time at which we transition from spec/relativistic to starsmasher/Newtonian physics by starting from two different snapshots 1.6 ms apart and finding negligible variation in the evolved energy histograms, which are shown in Fig. 1. In fact, even the initial energy histograms are not very different, so after the first ≈2 ms the influence of pressure on the kinematics of the ejecta is negligible. The codes’ main function is to provide the density evolution as particles follow their ballistic trajectories.

Figure 1.

The distribution of specific kinetic plus potential energy in the unbound post-merger matter, shown for system M12-7-S9 at a time shortly after the disruption of the NS (t = 0) and 500 ms later, long after the distribution has settled. For each energy bin, we integrate the density of all particles with energy inside that bin, giving a Newtonian mass for each energy bin. We show two resolutions, ‘Res1’ and ‘Res2’, corresponding to around 79 000 and 175 000 particles, respectively. We evolve using two methods: ‘P1’ turns off pressure forces and imposes adiabatic internal energy evolution within a radius of about 100M_BH. ‘P2’ includes pressure forces everywhere but removes bound particles after 10 ms. Another SPH run begun 1.6 ms earlier in the merger has nearly stationary energy distribution if evolved with P1. A simulation using P2 with a Paczyński–Wiita potential gives results almost identical to P2 with the standard Newtonian point potential.

Our interest is only in unbound matter. The bound matter for the most part orbits the black hole in an accretion disc or is ‘eaten’ when it comes within the prescribed distance from the central point mass. Both because of the exclusion of full relativity and the lack of a transport process to drive accretion, the disc evolution cannot be regarded as believable. We find that, for P2 evolutions, if we allow the disc to evolve for long periods of time, some fraction of the mass becomes weakly unbound. This is not perhaps incorrect given the physics included, but it cannot be regarded as physical, so we remove this contamination by eliminating bound particles after 10 ms of SPH evolution. For P1 evolutions, this removal is not necessary.

Our standard evolutions use roughly 75 000 particles. The mass profile of the ejecta M12-7-S9 is shown in Fig. 2. Fig. 3 shows snapshots of the SPH particles and fluid density after ≈0.5 s of starting the SPH evolution.

Figure 2.

The mass (computed by density integral) interior to a given cylindrical radius r_cyl or vertical height |z| for binary ejecta M12-7-S9. Profiles are computed at a time 500 ms after merger, by which point the ejecta profile has settled and will thereafter spread nearly homologously. The vertical interior mass appears to asymptote to a nonzero value on the left, indicating that a significant number of particles remain near the equator. We show results for two resolutions with two ways of handling pressure forces. Simulations with pressure forces completely turned off give profiles nearly the same as P1 profiles.

Figure 3.

Outflow's profile at a representative snapshot from the SPH evolution of M12-7-S9. Upper/bottom-left panels: xy/xz-projection of SPH-particles; upper/bottom-right panels: density (in log-scale [g cm⁻³]) views in the xy/xz-planes, respectively.

2.3 Nuclear reaction network and weak interactions

To calculate the composition of a Lagrangian fluid element in the ejecta, we require the evolution of its density as a function of time as well as its initial composition and entropy. To allow evolution at very late times, we extrapolate the density histories of the particles taken from the SPH simulation assuming homology, ρ ∝ t⁻³, which accurately describes the long-term evolution of the flow. In addition to the density, we extract the entropy and electron fraction along these trajectories. The extracted electron fraction is constant due to the neglect of weak reactions during the SPH evolution. The LS220 EOS is only valid for baryon densities and temperatures above ρ = 10⁸ g cm⁻³ and T ≈ 1 GK, which does not include the entire region in the temperature density plane in which neutron capture occurs. Since corrections to the EOS due to nuclear interaction become negligible below ρ ∼ 10¹² g cm⁻³, we switch from LS220 to a multispecies non-degenerate ideal gas EOS consistent with the nuclei included in our network along with the electron EOS of Timmes & Arnett (1999) for our network evolutions. When keeping the entropy fixed and assuming an initial nuclear statistical equilibrium (NSE) composition – with the modified Helmholtz EOS described in Lippuner & Roberts (2015) – we find temperature differences less than a few per cent between the two EOS in the region where they overlap. Because the single nucleus approximation of LS220 predicts different nuclei than a full NSE calculation, there is a mismatch between the total internal energies of ∼0.1 MeV/baryon when switching between the two EOSs. This level of error is unlikely to significantly impact the nucleosynthesis calculations because the total energy released per baryon during the nucleosynthesis is of the order of 10 MeV.

Once the density evolution of the Lagrangian particles has been extracted and extrapolated, we evolve the composition of the particles using the nuclear reaction network code skynet (Lippuner & Roberts 2015) and the network described therein. The entropy generated via nuclear transmutations is self-consistently included in the evolution, similarly to Freiburghaus et al. (1999). At 3 ms after merger – the time at which the spec simulations are mapped to starsmasher – the particles are typically at temperatures over 10 GK and densities are below 10¹² g cm⁻³. At these temperatures and densities, NSE holds, but weak interactions are generally far from equilibrium over the short time-scales encountered during the merger. To follow changes in the electron fraction at high temperature, skynet includes an NSE evolution mode where strong interactions are assumed to be in equilibrium and only weak interactions are tracked. This mode is used until the temperature drops below 7 GK, at which point the full nuclear reaction network is evolved. Because inverse strong reactions are calculated via detailed balance, the transition between the two skynet evolution modes is smooth.

To track the potential importance of neutrino irradiation of the ejecta, electron neutrino capture, electron antineutrino capture, electron capture and positron capture by free nucleons are included in both evolution modes. The neutrino capture rates are given by

\begin{eqnarray} \lambda _\nu &=& \frac{G_{\rm F}^2 (1 + 3g_{\rm A}^2)}{2 \pi ^2 \hbar ^7 c^6} \nonumber \\ &&\times \int _{\tilde{Q}}^{\infty } d\epsilon _{\rm e} p_{\rm e} \epsilon _{\rm e} (\epsilon _{\rm e} - Q)^2 \skew4\bar{f}_\nu (\epsilon _{\rm e} - Q) (1 - f_{\rm e}(\epsilon _{\rm e})), \end{eqnarray}

(1)

where f_e is the electron distribution function, G_F is the Fermi coupling constant, g_A is the weak axial vector coupling constant, ε_e is the electron energy, p_e is the electron momentum, Q is the energy transfer from the nucleons to the final state electron, and |$\bar{f}_\nu$| is the angle-averaged neutrino distribution function. The Q-value is defined in the direction of electron or positron capture and |$\tilde{Q} = \max (Q,m_{\rm e} c^2)$|⁠. The electron and positron capture rates, |$\lambda _{{\rm e}^+}$| and |$\lambda _{{\rm e}^-}$|⁠, are calculated from similar expressions with the distribution functions interchanged. This expression assumes that there is no momentum transfer to the nucleons and neglects weak magnetism corrections. Although these corrections are potentially significant in the case of neutrino-driven winds (Horowitz 2002), they are unlikely to significantly affect the evolution of the electron fraction in BHNS merger ejecta. The α-effect locks free protons in heavy nuclei and thereby prevents significant competition from electron antineutrino capture (Fuller & Meyer 1995). We assume that the neutrino distribution has a Fermi–Dirac shape in energy space and neutrinos of all energies are emitted from a single spherical surface, which results in the distribution function

\begin{equation} f_\nu (\epsilon , \mu , r) = \frac{\theta (\mu - \mu _0(r))}{\exp (\epsilon /T_\nu ) + 1}, \end{equation}

(2)

where μ is the cosine of the angle of neutrino propagation relative to the radial direction, |$\mu _0 = \sqrt{1-(r_\nu /r)^2}$|⁠, T_ν defines the neutrino spectral temperature, θ is the Heaviside step function, ε is the neutrino energy and r_ν is the radius of neutrino emission. Inside of r_ν, μ₀ is assumed to smoothly approach negative one over a tenth of r_ν. The value of r_ν can be fixed by choosing a neutrino luminosity, L_ν, and spectral temperature. This model is crude, considering the disc-like geometry of the neutrino emitting region, but it is sufficient for this study given that we are parametrizing the properties of the neutrino field anyway. In the following sections, we consider models with fixed electron neutrino luminosities of |$L_{\nu _{\rm e}} = \lbrace 0, 0.2, 1, 5, 25\rbrace \times 10^{52} \, \mathrm{erg s}^{-1}$|⁠. The electron antineutrino luminosity is always fixed to be |$L_{\bar{\nu }_{\rm e}}=1.5 L_{\nu _{\rm e}}$|⁠, but our results are insensitive to this choice due to the α-effect. These values are in the range found in the simulations of Foucart et al. (2016) and the difference between the values accounts for re-leptonization of the disc. Since only charged current interactions are included in the nuclear network, the properties of the heavy flavoured neutrino fields do not affect our results. We employ constant luminosities to reduce the number of parameters affecting our nucleosynthesis calculations.

We perform nucleosynthesis post-processing for all of the ejected SPH trajectories. The nuclear evolution is followed until 10¹³ s after the merger, which allows for the decay of all but a handful of long-lived unstable isotopes.

3 RESULTS AND DISCUSSION

3.1 The electron fraction of the ejecta

The electron fraction of the material ejected during the BHNS merger is the most important parameter in determining the nucleosynthesis which occurs within the outflow (e.g. Lippuner & Roberts 2015). Given the short dynamical time-scales and the lack of a hypermassive central NS after the merger, it has often been assumed that the electron fraction of the dynamical ejecta from BHNS mergers is set solely by the initial beta-equilibrium electron fraction of the NS from which the material was ejected (Just et al. 2015). If there are not a substantial number of weak interactions during and after the merger, the electron fraction will be low enough that an r-process involving a significant number of fission cycles will occur: the outer layers of an NS have Y_e < 0.1 and the critical value for producing r-process material at low entropy is Y_e ≈ 0.25 (e.g. Kasen, Fernández & Metzger 2015; Lippuner & Roberts 2015). Neutrinos can impact the electron fraction of the ejecta of binary NS mergers (Wanajo et al. 2014; Foucart et al. 2016; Goriely et al. 2015; Palenzuela et al. 2015; Radice et al. 2016). In binary NS mergers, a large fraction of the prompt ejecta comes from the shock-heated material in the interaction region of the two NSs (Palenzuela et al. 2015). The increased temperatures and the large neutrino fluences near this material increases Y_e significantly and can sometimes drastically alter the character of nucleosynthesis in the outflow. In the BHNS case, there is no interaction region during the tidal disruption of the NS, and matter ejection when the tidal stream self-intersects is very subdominant (Foucart et al. 2015). The case M14-5-S9 has the most massive ejection from the tidal stream collision (Deaton et al. 2013) of these BHNS, but even for this case the imprint of this secondary ejecta source on the overall outflow composition is small. Therefore, the ejected material has a lower average entropy and electron fraction than neutron star–neutron star (NSNS) merger ejecta and there is no significant neutrino emission until a disc has formed around the BH. Here, we consider the extent to which neutrino interactions can alter the distribution of Y_e just before r-process nucleosynthesis begins in the ejecta.

We estimate the effect of neutrino captures on the BHNS outflows by considering the maximum disc neutrino luminosities found by Foucart et al. (2014). The neutrino luminosity coming from the disc in both electron neutrinos and antineutrinos is around 10⁵³ ergs⁻¹. Although the simulations of Foucart et al. (2014) used a grey leakage approximation, we can get some estimate of the average neutrino energies from the temperature of the emission region which was around 5 MeV, which suggests average neutrino energies around ε_ν ≈ 3.15T ∼ 15 MeV (e.g. Foucart et al. 2015). We can then estimate the neutrino processing time-scale as

\begin{equation} \tau _\nu (r) \approx 67.8 \, \mathrm{ms} \, \left(\frac{r}{250 \, \mathrm{km}} \right)^{2} L_{\nu _{\rm e},53}^{-1} T_{\nu _{\rm e},5}^{-1} \,, \end{equation}

(3)

where r is the radius of the fluid element, |$L_{\nu _e,53}$| is the electron neutrino luminosity in units of 10⁵³ erg s⁻¹ and |$T_{\nu _e,5}$| is the electron neutrino spectral temperature in units of 5 MeV. Electron antineutrinos are unlikely to contribute significantly to the neutrino interaction time-scale. This is because in the low entropy outflows of BHNS mergers almost all protons are locked in heavy nuclei and thus have very low neutrino capture cross-sections.

The change in Y_e due to neutrino interactions can be estimated by assuming that the tidal ejecta has a constant velocity v, the neutrino luminosity is constant, electron and positron capture are unimportant, protons are locked into heavy nuclei, and there is a finite time after merger at which neutrinos start being emitted from the disc. With these assumptions, the evolution of Y_e as a function of radius is given by

\begin{equation} \frac{{\rm d} Y_{\rm e} }{{\rm d}r} = \frac{\theta (r-vt_{\nu ,{\rm on}})}{v \tau _\nu (r) Y_{e,{\rm eq}}} \left( Y_{e,{\rm eq}} - Y_{\rm e} \right), \end{equation}

(4)

where Y_e,eq = 〈Z〉_nuclei/〈A〉_nuclei, and t_ν,on is the time after merger at which the neutrino luminosities reach their saturation value.

Assuming a constant average proton and neutron numbers of the heavy nuclei, this can easily be integrated to large radius to find the final electron fraction

\begin{eqnarray} Y_{e,f} &\approx& Y_{e,{\rm eq}} \left[1 - \exp \left(-\frac{r_0 }{v \tau _\nu (r_0) Y_{e,{\rm eq}}}\right) \right] \nonumber \\ &&+\, Y_{e,i} \exp \left(-\frac{r_0}{v \tau _\nu (r_0) Y_{e,{\rm eq}}}\right), \end{eqnarray}

(5)

where r₀ = t_ν,onv. Using the outflow velocity and neutrino luminosities calculated in the M12-7-S9 model of Foucart et al. (2014) (v ≈ 0.25 c, |$L_{\nu _{\rm e}} \approx 10^{53} \, \mathrm{erg\,s}^{-1}$| and t_ν,on ≈ 3 ms), we find that the post neutrino interaction electron fraction is Y_{e, f} ≈ 0.07 if the Y_e,eq is close to a half. Given that the r-process is robustly produced below Y_e ≈ 0.25, this suggests that neutrino interactions are much less likely to play a significant role in determining the composition of the ejecta in BHNS mergers relative to binary NS mergers, although this estimate is sensitive to t_ν,on and the velocity of the outflow.

To make this more concrete, we run nucleosynthesis calculations for the M12-7-S9 model including neutrino interactions induced by a constant neutrino luminosity, modelled as described above. Similar results are found for the other two models discussed in Section 2.1. In Fig. 4, the weak interaction rates and the electron fraction are shown for a single particle. Because our Lagrangian trajectories start at 3 ms after the merger, the initial density in the ejected material is below about 10¹⁰ g cm⁻³ and lepton captures are dominated by neutrino captures for neutrino luminosities above about 10⁵² erg s⁻¹. The neutrino interaction rates fall off as a power law in time, since this particular particle is moving away from the merger site at constant velocity in a nearly radial direction. Other particles can deviate from power-law behaviour at early times, but not strongly. As was expected from our estimates above, the neutrino interaction time-scale is long compared to the outflow time-scale and very little evolution of the electron fraction occurs during the first 10 ms. The evolution of Y_e after about 20 ms is driven by beta-decays occurring during the r-process.

Figure 4.

Evolution of the electron fraction and weak rates as a function of time for a characteristic fluid element. The electron neutrino luminosity is assumed to be 10⁵³ erg s⁻¹. Because of the relatively low entropy of the BHNS ejecta and because of the low initial density of our calculations, neutrino interaction rates dominate the electron and positron capture rates but neither have a large impact on the electron fraction of the outflow. The increase in Y_e seen after around 100 ms is due to beta-decay during the r-process.

To look at the effect of weak interactions globally, the distribution of Y_e in the material ejected in model M12-7-S9 is shown in Fig. 5 for a range of assumed neutrino luminosities. The GRHD simulations described in Section 2.1 include electron and positron captures, but do not include neutrino captures. The SPH simulations which follow the long-term evolution of the ejecta include no weak interactions. Therefore, we include weak interactions in our post-processing nucleosynthesis calculations to assess their impact on Y_e. As we expect, the ejected material is very neutron rich, but becomes slightly less neutron rich with increasing electron neutrino luminosity. The distribution of the electron fraction in the whole NS is also shown to emphasize that the ejecta in the absence of neutrinos has a significantly lower Y_e than the average Y_e of a cold 1.2 M_⊙ NS calculated using the LS220 EOS. The beta-equilibrium value of Y_e increases with density, so that the outer layers of the NS – which comprise most of the ejecta – have a lower electron fraction. The average electron fraction in the ejecta is 0.053, 0.053, 0.054, 0.062 and 0.127 for neutrino luminosities of {0, 0.2, 1, 5, 25} × 10⁵²erg s⁻¹.

Figure 5.

Mass-weighted histogram of the electron fraction in the ejecta from model M12-7-S9 assuming fixed electron neutrino luminosities of {0, 0.2, 1, 5, 25} × 10⁵² erg s⁻¹. For comparison, we also show the electron fraction histogram in a 1.2 M_⊙ LS neutron star (cyan line).

3.2 Nucleosynthesis and neutrino-induced production of the first r-process peak

We now consider the detailed nucleosynthesis in the ejecta of model M12-7-S9, both with and without neutrinos. We focus on the effect neutrinos can have on the isotopic abundances of the ejecta. In Fig. 6, the integrated nucleosynthesis from model M12-7-S9 is shown. Since the neutrino emission from the accretion torus formed after the BHNS merger is uncertain, we calculate the final nucleosynthetic yields of M12-7-S9 assuming the range of electron neutrino luminosities listed above. In all cases, the electron antineutrino luminosity is fixed at |$1.5 L_{\nu _{\rm e}}$| to very approximately account for re-leptonization of the neutrino emitting disc (Foucart et al. 2016). Because of the α-effect, the results are insensitive to the chosen electron antineutrino luminosity. The electron neutrino and antineutrino average energies are fixed at 12 MeV and 15 MeV, respectively. The results for the other two simulated binary systems are similar and they are discussed briefly below. We also find that in the case of zero neutrino luminosity, the nucleosynthesis results are not significantly altered when we use parametrized density histories with a fixed dynamical time-scale for all particles. This is not surprising, given the low electron fraction encountered in the ejecta, which gives rise to robust r-process nucleosynthesis across a wide range of dynamical time-scales. Nevertheless, the position history of the particles provided by our SPH code is necessary for estimating the amount of neutrino irradiation the particles undergo.

Figure 6.

Comparison of the integrated nuclear abundances in model M12-7-S9 assuming different fixed neutrino irradiation from the nascent accretion disc. We also include the classical scaled solar abundance r-process distribution from Arlandini et al. (1999) for comparison. The abundances for a selection of single particles from the |$L_{\nu _{\rm e}}=0$| and |$L_{\nu _m,52} = 5$| runs are also shown as light lines. For all runs, we assume |$L_{\bar{\nu }_{\rm e}} = 1.5 L_{\nu _{\rm e}}$|⁠.

In general, we confirm previous work which has shown BHNS mergers dynamically eject a large amount of r-process rich material (e.g. Roberts et al. 2011; Korobkin et al. 2012; Just et al. 2015). Both the second and third r-process peaks are robustly produced, independent of the neutrino luminosity. Given the low electron fractions found in the ejecta at the start of neutron capture, robust production of the r-process is not surprising (Lippuner & Roberts 2015). In all of the models, reactive flow proceeds past the third peak before neutron exhaustion occurs in the vast majority of the simulated fluid elements and they undergo fission cycling. We find that fission cycles occur in the ejecta and the number of cycles is weakly dependent on the neutrino luminosity (for the luminosities considered here). Therefore, the abundance pattern above mass number ∼90 is likely to be robust to variations in the total neutrino luminosity and the properties of the merging system. In all models, the third r-process peak is overproduced relative to the second peak. This is discussed further in Section 3.4.

We find that the abundance of the first r-process peak at mass number 78 depends on the neutrino luminosity, in contrast to the second and third peaks which are independent of the neutrino luminosity. None the less, in all cases it is underproduced relative to the solar abundance when normalizing to the second and third peaks. This first peak production is driven by low-mass r-process seed production after material falls out of NSE. This material is composed of heavy nuclei and free neutrons when strong equilibrium ceases to hold. Since the material is still relatively close to the accretion torus a few milliseconds after it is ejected, a significant number of electron neutrinos can be captured by the free neutrons. The produced protons then rapidly capture neutrons and form deuterium, which can then capture another deuteron to form an alpha particle. These alpha particles can then undergo a neutron-catalysed triple-alpha reaction, similar to what occurs in neutron-rich neutrino-driven winds (Delano & Cameron 1971; Hoffman, Woosley & Qian 1997), to produce low-mass seed nuclei for the r-process (Meyer, McLaughlin & Fuller 1998). This non-equilibrium neutrino-induced seed production creates a distinct set of seed nuclei which can undergo neutron capture, since the seeds produced by the NSE distribution tend to be between mass 78 and 100. A large number of the low-mass seeds do not get processed past the N = 50, Z = 28 point in the r-process path before neutron exhaustion occurs because of the long beta-decay half lives in that region of the chart of the nuclides. Therefore, these neutrino-produced seed nuclei are responsible for producing the first peak r-process nucleosynthesis seen in our simulations. This effect of neutrino irradiation of the outflow is distinct from the one discussed by Wanajo et al. (2014) and Goriely et al. (2015), where the neutrino luminosities are high enough to push the electron fraction over ∼0.25 and stop production of the second and third peaks.

3.3 Details of the first peak production mechanism

We now consider the details of the process by which abundances in the first peak are indirectly produced by electron neutrino captures by neutrons. The total number fraction of heavy nuclei produced by neutrino induced seed production can be estimated by using the results from Section 3.1 as follows. Low-mass seed production proceeds via the neutron catalysed triple alpha process, so it takes six protons to make a seed nucleus. The rate of proton production is just |$\dot{Y}_{\rm e}$|⁠, so the total number of low-mass seed nuclei produced by neutrino interactions is

\begin{eqnarray} Y_{\mathrm{s},\nu } &\approx & \frac{Y_{{\rm e},f} - Y_{{\rm e},i}}{6} \nonumber \\ &=& \frac{Y_{{\rm e},{\rm eq}} - Y_{{\rm e},i}}{6} \left[1 - \exp \left(-\frac{r_0}{v \tau _\nu (r_0) Y_{{\rm e},{\rm eq}}}\right) \right]. \end{eqnarray}

(6)

This estimate implies that around 2 × 10⁻³ seed nuclei per baryon are produced by this process, assuming that the neutrino luminosity is 10⁵³ erg s⁻¹. This number is in good agreement with the values extracted from our nucleosynthesis calculations. Comparing this to the total final abundance of the first peak for a single ejecta particle shown in Fig. 7, it is clear that only about 10 per cent of this material gets trapped in the first peak.

Figure 7.

Illustration of how the first r-process peak is produced by electron neutrino captures on neutrons for a single SPH particle. This SPH particle had initial Y_e = 0.11, initial entropy s = 9.7 k_B baryon⁻¹, and an asymptotic velocity v/c = 0.5. Top panel: the solid lines show the abundance of material in the first r-process peak, Y_1st, as a function of time (i.e. material with 72 ≤ A ≤ 79), the dashed lines show the integrated number of protons produced by weak interactions after time t divided by 6, |$Y_{{\rm s},\nu } = \int _t^{\infty } {\rm d}t Y_{\mathrm{n}}/6 (\lambda _{\nu _{\rm e}} + \lambda _{{\rm e}^+})$|⁠, and the dotted lines show the neutron abundance Y_n. Y_s,ν gives the number of low-mass seed nuclei produced by neutrino interactions. The neutrino seed nuclei produced at early times are burned past the first r-process peak, but the seed nuclei produced after the time when Y_s,ν = Y_1st,final do not get burned passed the first peak before neutrons are exhausted, and so they will end up in the first peak. Bottom panel: the solid lines show the temperature of the particle as a function of time, the dashed lines show the time-scale to process material to the first peak, τ_{(6, 25)}, and the dotted lines show the destruction time-scale of the first peak, τ_(26,28), which are defined in the text. In this particle, there is no significant variation with neutrino luminosity of the temperature or the r-process path. Therefore, the two time-scales do not change with the amount of neutrino irradiation.

Seed nuclei indirectly produced by neutrinos are not processed past the N = 50 closed shell rapidly. If such rapid processing were the case, the final amount of mass in the first peak would be set by the number of seed nuclei produced after a time just before neutron exhaustion. To illustrate when the nuclei trapped in the first peak are produced, we show the total number of seed nuclei produced by neutrino interactions after time t

\begin{equation} Y_{{\rm s},\nu }(t) = \frac{1}{6}\int _t^\infty Y_{\rm n} (\lambda _{\nu _{\rm e}} + \lambda _{{\rm e}^+}) {\rm d}t \end{equation}

(7)

in Fig. 7, along with the time dependence of the first peak abundance, Y_1st, and the neutron abundance Y_n. Y_s,ν is just the number fraction of protons produced by weak interactions after time t divided by 6, since it requires six protons to produce a seed nucleus which can capture neutrons. Material will be processed through the first peak on some time-scale τ_1st. Let t_ex be the time at which neutrons are exhausted and t_prod be the time after which neutrino-produced seed nuclei get trapped in the first peak. Seed nuclei produced at times earlier than t_prod = t_ex − τ_1st will be burned past the first peak, while seed nuclei produced within a time τ_1st of neutron exhaustion will end up in the first peak. We can estimate τ_1st by looking for solutions of Y_s,ν(t_prod) = Y_1st,final. Inspecting Fig. 7, we find t_prod is 70–100 ms and t_ex is 520–600 ms for |$L_{\nu _{\rm e},52}$| ranging from 20 to 1. Thus we estimate that τ_1st, the time it takes for seed nuclei to be processed to the N = 50 closed shell of the first peak, is between 450 and 500 ms for this particular fluid element. For reference, the final abundances of this particle are shown in Fig. 8.

Figure 8.

Final abundances as a function of neutrino luminosity for the single Lagrangian particle shown in Fig. 7.

We now attempt to explain what sets this time-scale. Assuming the waiting point approximation (cf. Kratz et al. 1993), the time-scale to go from charge Z₁ to charge Z₂ is given by

\begin{equation} \tau _{(Z_1,Z_2)}(t) = \sum _{Z=Z_1}^{Z_2} \frac{\sum _N Y_{(Z,N)}}{\sum _N Y_{(Z,N)} \tau ^{-1}_{\beta ^{-},(Z,N)}}. \end{equation}

(8)

Here, |$\tau ^{-1}_{\beta ^{-}, (Z,N)}$| is the beta-decay time-scale of a nucleus with N neutrons and Z protons. When (n, γ) reactions are in equilibrium with (γ, n) reactions – such that μ_n + μ_{(Z, N)} = μ_{(Z, N + 1)} – these time-scales are only functions of the density, temperature and neutron abundance, i.e. |$\tau _{(Z_1,Z_2)} = \tau _{(Z_1,Z_2)}(\rho , T, Y_{\rm n})$|⁠. This is often the case at the high temperatures encountered during r-process nucleosynthesis in these outflows, but the equilibrium can start to break down at lower temperatures (Kratz et al. 1993). This makes it clear that changing the temperature and electron fraction of a particular fluid element can change the path of the r-process and alter the time it takes material to be processed from one charge number to another. The quantities τ_{(6, 25)} and τ_{(26, 28)} are shown in the bottom panel of Fig. 7. Note that τ_{(6, 25)} + τ_{(26, 28)} is approximately the time it takes a seed nucleus to get to the start of the first peak and then get processed through the first peak, which we called τ_1st above. We see from the bottom panel of Fig. 7 that τ_{(6, 25)} is constant throughout the period during which the r-process is occurring and its value is in good agreement with our 450–500 ms estimate for τ_1st. Because this time-scale is determined by beta decay, the final first peak abundance goes linearly with the neutrino luminosity. The lifetimes of isotopes along the N = 50 closed shell are 40, 110 and 110 ms, for the reactions ⁷⁶Fe(β⁻, n)⁷⁵Co, ⁷⁷Co(β⁻, n)⁷⁶Ni and ⁷⁸Ni(β⁻)⁷⁸Cu. These are consistent with the τ_{(26, 28)} ≈ 100 ms we find in Fig. 7. We also note that at around 600 ms into the calculation – which is after neutron exhaustion – there is a further increase in the first peak abundance. This is driven by the reaction ⁸⁰Ni(β⁻, n)⁷⁹Cu. Significant production of ⁸⁰Ni occurs just before neutron freeze-out and it has a half-life of 175 ms.

This suggests that the neutrino flux between times t₁ = t_ex − τ_{(6, 25)} − τ_{(26, 28)} and t₂ = t_ex − τ_{(6, 25)} will determine the amount of neutrino-induced first peak production which occurs. Seed nuclei produced before time t₁ will get beyond the first peak before neutrons are exhausted, while seed nuclei produced after time t₂ will not reach the first peak before neutron exhaustion occurs. Therefore, the important quantity for understanding neutrino-induced production of the first peak will be the neutrino luminosity centred at a time around 70 ms after merger, within a window of around 100 ms.

We have shown that the production of first peak nuclei goes linearly with the electron neutrino luminosity for the specific Lagrangian particle shown in Fig. 7, but Fig. 6 shows that production of the first peak appears to saturate at luminosities above ∼5 × 10⁵² ergs⁻¹. Below this luminosity, the dependence of first peak production on luminosity is approximately linear as expected. In Fig. 9, we show a different Lagrangian particle which exhibits non-monotonic behaviour of the first peak abundance with the neutrino luminosity. The first peak abundance increases at low luminosity, decreases with luminosity around |$L_{\nu _{\rm e}} = 10^{53} \, \mathrm{erg s}^{-1}$|⁠, and then increases with luminosity again. This particle has a lower asymptotic velocity than the particle shown in Fig. 7 and therefore experiences more neutrino irradiation. Additionally, it has lower initial entropy and Y_e, which means neutrino interactions can have a larger effect on its thermodynamic state.

Figure 9.

The same as Fig. 7, except for a different thermodynamic trajectory. This SPH particle had initial Y_e = 0.05, initial entropy s = 4.33 k_B baryon⁻¹, and an asymptotic velocity v/c = 0.29. Because of the lower velocity, lower initial entropy, and lower Y_e present in this particle relative to the particle shown in Fig. 7, neutrino interactions significantly alter the thermodynamic state of the material and τ_{(6, 25)}. This causes the first peak abundance to vary non-monotonically with the neutrino luminosity.

The lower panel of Fig. 9 clearly shows that neutrinos significantly alter the thermodynamic state of the considered particle and that the low-mass r-process path is shifted by the inclusion of neutrino interactions. In particular, τ_{(26, 28)} and τ_{(6, 25)} increase with the neutrino luminosity. Fig. 10 shows how the r-process path varies with the neutrino luminosity 100 ms into the calculation, giving rise to the processing time-scales dependence on neutrino luminosity. The total number of seeds increases with initial Y_e and temperature, corresponding to larger neutrino luminosities in this fluid element. Additionally, increasing the entropy of the outflow reduces the rate at which material can bypass the A = 8 stability gap. The large difference in the first peak processing time-scale, τ_{(26, 28)}, seen in Fig. 9, is due to the r-process path shifting from being far beyond the N = 50, Z = 28 closed shells at lower temperatures (and lower neutrino luminosities) to proceeding through closed shells at higher temperatures (and higher neutrino luminosities). This significantly alters how first peak nuclei are produced throughout the calculation and breaks the linear dependence on the neutrino luminosity.

Figure 10.

The r-process path for the SPH particle shown in Fig. 9 for different neutrino luminosities at 100 ms into the calculation. The inset shows the mass summed abundances at the same time. Notice how the path differs for different neutrino luminosities.

Even in the absence of neutrinos, there is some production of first peak nuclei. As we have mentioned, this material is produced by fission of heavy nuclei. Since we are employing symmetric fission fragment distributions, it is likely that more realistic fission fragment distributions will result in a broader distribution of fission daughters and more material being left behind in this region. None the less, it seems likely that there will be at least some production of the first peak even in the very neutron-rich outflows of BHNS mergers, as long as neutrino luminosities from the post-merger remnant are above about 10⁵² ergs⁻¹ within a hundred milliseconds of the merger. We also emphasize that neutrino-induced production of the first peak does not produce enough material in our models to agree with the solar r-process abundances when they are normalized to the second peak. Instead, the abundance is around an order of magnitude too low.

3.4 Isotopic and elemental abundances, Galactic chemical evolution, and low-metallicity halo stars

In Fig. 11, the integrated abundances from the models M12-7-S9, M14-7-S8 and M14-5-S9 are shown for a fixed neutrino luminosity of 10⁵³ ergs⁻¹ (and |$L_{\bar{\nu }_{\rm e}} = 1.5 L_{\nu _{\rm e}}$|⁠). Clearly, there is little discernible difference between the predicted nucleosynthesis from these models. The electron fraction in almost all of the ejecta in all three models is below the threshold for fission cycling to occur (Lippuner & Roberts 2015) and the entropy of the ejecta is quite low. These are conditions which result in a second and third peak nuclear abundance pattern which is quite insensitive to the detailed properties of the ejecta (Lippuner & Roberts 2015, e.g.). Interestingly, the neutrino-produced first peak is also very insensitive to the binary parameters if the same neutrino luminosities are assumed. A priori, it would seem that different dynamics during the merger could give rise to different dynamics of the ejecta and alter the radius at which neutron exhaustion occurs. Of course, the important parameter will be the velocity of the ejecta. Larger velocities will result in lower local neutrino densities around the time at which neutrons are exhausted in the ejected material.

Figure 11.

The integrated nuclear abundances of the dynamical ejecta in the models M12-7-S9, M14-7-S8 and M14-5-S9 assuming an electron neutrino luminosity of 10⁵³ ergs⁻¹. The low electron fraction of all the ejecta results in production of a robust r-process, independent of the dynamics of the merger and parameters of the binary system.

Although there is reasonable qualitative agreement with the solar r-process abundance pattern above A ∼ 100 in all of our models, there are significant quantitative differences. Most clearly, the third peak is significantly overproduced and has an offset with respect to the observed solar pattern. Uncertainties in the ejecta abundance pattern can come from two sources, uncertainties in the properties of the ejecta and uncertainties in the nuclear data which serve as input for our nucleosynthesis calculations. Given how robust the r-process pattern in our models is above A ∼ 100 to variations in the binary parameters and neutrino irradiation, it seems unlikely – for our chosen nuclear data – that BHNS mergers can make a pattern which agrees exactly with the r-process pattern found in the Sun and low-metallicity halo stars. Strictly following this argument to its conclusion, BHNS mergers would be ruled out as the dominant contributor to the galactic chemical evolution of r-process elements. This would put a significant constraint on the combined merger rate, the BH spin distribution and the BH mass distribution in these binaries (Bauswein et al. 2014a). Of course, it is easy to imagine scenarios where the galactic r-process nucleosynthesis is produced by multiple types of events, so strongly ruling out a single r-process production channel on the basis of inexact agreement with the solar pattern seems premature at best.

The second possible source of uncertainty in our results is the input nuclear physics data. The nuclear masses, beta-decay rates, neutron capture rates, fission barrier heights and fission fragment distributions in the r-process path, which lies far from nuclear stability, have, on the whole, not been experimentally determined but are instead determined from models which are in part constrained by data from nuclei closer to stability (e.g. Möller, Nix & Kratz 1997; Goriely, Chamel & Pearson 2009). For instance, different nuclear mass models can give significantly different abundance patterns for the same thermodynamic histories (e.g. Arcones & Martínez-Pinedo 2011; Martin et al. 2015; Mendoza-Temis et al. 2015). By varying masses within a particular mass model within the expected uncertainty, Mumpower et al. (2015) have shown that the uncertainties in the final r-process abundance pattern solely due to nuclear physics uncertainties can be as a large as a factor of 10. We have also not included neutrino-induced fission in our nuclear network (Qian 2002; Kolbe, Langanke & Fuller 2004), which could change the nuclei which undergo fission and potentially alter the low-mass r-process pattern. Nevertheless, the neutrino irradiation in these outflows is rather weak so it would be surprising if neutrino-induced fission drastically changed our results. Therefore, given the level of agreement we find with the solar r-process isotopic abundance pattern, our results seem wholly consistent with BHNS mergers contributing to the galactic budget of heavy r-process nuclei.

For confrontation with observations of abundances in low-metallicity halo stars, it is more instructive to examine the elemental abundance pattern of the ejecta (i.e. |$Y_{\rm Z} = \sum _i \delta _{Z_i, {\rm Z}} Y_i$|⁠), since only elemental abundances are easily determined in these stars. In Fig. 12, we show the final elemental abundances in the ejecta of model M12-7-S9 for a variety of assumed neutrino luminosities. In the region 51 < Z < 81, most of the abundances agree with the solar pattern to within a factor of 3. The notable exceptions are gold and platinum in the third peak and cesium in the second peak, all of which are overproduced by a factor of around 10. Although the agreement is not perfect, our patterns above Mo are within the errors due to uncertainties in the nuclear physics input (Mumpower et al. 2015).

Figure 12.

Top panel: the integrated elemental abundances in the ejecta of model M12-7-S9 for a range of imposed neutrino luminosities compared to arbitrarily scaled solar r-process abundances. We include the name of every other element to guide the eye. Bottom panel: ratio of the calculated abundances to the solar abundances (from Arlandini et al. 1999). In the lower panel, factor of 3 error bars are included on our calculated abundance patterns to approximately account for errors in the input nuclear physics. The solar r-process abundances have been transformed by an overall scaling factor to minimize the residuals in log space in the region 40 < Z < 81. We also include abundances from the low-metallicity halo stars HD 115444 and 122563 (Westin et al. 2000) and HD 108317 and 128279 (Roederer et al. 2012, 2014), with [Fe/H] of −3.0, −2.7, −2.5 and −2.5, respectively. The observational error bars are smaller than the plotted symbols and the abundances patterns have been scaled by a factor to minimize the deviation from the solar r-process in the range 40 < Z < 81.

For elements below Mo, there are a number of primary nucleosynthesis processes which can contribute to these abundances, even in low-metallicity stars (Travaglio et al. 2004; Montes et al. 2007; Qian & Wasserburg 2008; Arcones & Montes 2011; Hansen, Montes & Arcones 2014). Therefore, it is hard to rule out our models because they fail to produce certain abundances below Mo. As can be seen in the bottom panel of Fig. 12, the abundances below Mo are significantly underproduced relative to solar in all of our models. Nevertheless, neutrino-induced r-process seed production brings these abundances up much closer to the solar value.

Since many processes contribute to the solar abundances, it is more reasonable to compare our abundances with those measured in the atmospheres of low-metallicity halo stars. In Fig. 12, we show the abundances of four low-metallicity halo stars with subsolar Ge abundances. The abundance data are taken from Westin et al. (2000), Roederer et al. (2012) and Roederer et al. (2014) for the stars HD 115444, 122563, 108317 and 128279, which have metallicities of [Fe/H] = −3.0, −2.7, −2.5 and −2.5, respectively. These metallicities are low enough such that s-process contamination from low-mass stars is highly unlikely (e.g. Simmerer et al. 2004). None the less, it is possible that the s-process in massive stars may have contributed to these abundances (e.g. Woosley, Heger & Weaver 2002).

We can compare our predicted abundances with observed abundances of Ge, As and Se in low-metallicity halo stars. For these elements, neutrino-induced r-process seed production of the first peak brings the abundances of Ge, As and Se closer to agreement with the observed values for larger values of the electron neutrino luminosity, even though the abundances are still small by close to a factor of 10. Given the significant uncertainties in the level of neutrino irradiation and the properties of nuclei far from stability near and below the first peak, the actual amount of first peak material created by neutrino-induced seed production is highly uncertain. If first peak production could be increased by around an order of magnitude, certain low-metallicity halo star abundance patterns above Zn can possibly be explained by only invoking nucleosynthesis in the dynamical ejecta of BHNS mergers, although there are clearly a number of other possible ways to produce the first peak in the material ejected from the disc (e.g. Just et al. 2015). This contrasts with the results of Wanajo et al. (2014) for binary NS mergers, where they find that the first peak nuclei are produced with a ratio to the second and third peak which is close to the ratio observed in the sun and therefore are overproducing the first peak by about a factor of 10 relative to the low-metallicity halo star abundances. In contrast to Ge, As and Se, all the BHNS merger outflows underproduce Sr by orders of magnitude and the final abundances of Sr and Zr are insensitive to the level of neutrino irradiation.

4 CONCLUSIONS

In this work, we have considered nucleosynthesis in the ejecta of BHNS mergers including the effect of neutrino interactions. Starting from general relativistic hydrodynamic simulations of BHNS mergers with reasonably high BH spins, we have extracted the unbound material and mapped it into a Newtonian SPH code, starsmasher. We then evolved the ejecta over a long enough time to reach homologous expansion. Using the Lagrangian histories of the SPH particles, we then performed post-processing calculations of nucleosynthesis in the neutron-rich ejecta with skynet. In particular, we focused on the influence of neutrinos on the nucleosynthesis in these outflows by parametrizing the neutrino luminosity coming from the disc.

As is expected from previous work (e.g. Lattimer et al. 1977; Freiburghaus et al. 1999; Roberts et al. 2011; Korobkin et al. 2012; Foucart et al. 2014; Bauswein et al. 2014a), we find that the second and third r-process peaks are robustly produced in the outflow. In contrast to the case of NSNS mergers (Wanajo et al. 2014; Goriely et al. 2015; Sekiguchi et al. 2015), we find that – for reasonable luminosities of |$L_{\nu _{\rm e}} = \lbrace 0, 0.2, 1, 5, 25\rbrace \times 10^{52} \, \mathrm{erg\,s}^{-1}$| – neutrinos are unable to significantly shift the distribution of the pre-r-process electron fraction in the ejecta and almost all of the ejected mass elements produce the full r-process. This is due to the rapid outflow time-scales of the tidal ejecta and the relatively low expected neutrino luminosities of the disc formed around the remnant BH. Additionally, we find that there are negligible differences between the nucleosynthesis calculated for different BHNS binary parameters.

Although weak interactions have no gross effect on the electron fraction distribution, we find that they do alter nucleosynthesis in more subtle ways. Once r-process neutron captures have begun, electron neutrino captures can turn a small fraction of the neutrons into protons due to the α-effect (Fuller & Meyer 1995). These protons form alpha particles which rapidly combine to form ¹²C. This provides a source of low-mass seed nuclei for the r-process (Meyer et al. 1998) and results in nuclear flow below mass 78 which is not present in low entropy outflows when neutrinos are neglected. In the few hundred milliseconds before neutron exhaustion, this flow can build up a significant amount of material in the first r-process peak at mass 78. This process is likely to operate both in the dynamical ejecta of BHNS mergers and binary NS mergers, as well as in material which is ejected from the discs left over in these events, which we did not consider in this work.

Although this neutrino induced seed production produces the first peak, we find that it does not produce enough material to bring the final abundance distributions into agreement with the solar r-process distribution below mass 80. None the less, this first peak production will be interesting when comparing with the abundance patterns of low-metallicity halo stars. In particular, it can significantly increase the abundance of the elements Ge, As and Se present in these outflows.

Further work is required to determine how robust neutrino induced r-process seed production is to variations in beta-decay lifetimes and neutron capture rates around the N = 50 closed shell. Additionally, employing more realistic fission fragment distributions than the ones employed here may also effect the abundances found just above the first r-process peak. Of course, this process is sensitive to the electron neutrino flux in the outflow. Here, we have chosen to parametrize the neutrino luminosity. Therefore, better models of the electron neutrino irradiation of the outflow are required to determine if neutrino-induced r-process seed production is robust in nature. The results are likely to be somewhat sensitive to the beta-decay rates away from stability below Z = 28, so better measurements of decays along the r-process path could shed significant light on whether or not this process is important in the ejecta of BHNS mergers. It will also be interesting to investigate if this process operates in other environments with low Y_e outflows, namely the ejecta of binary NS mergers and material ejected from discs formed after compact object mergers.

Acknowledgments

LFR acknowledges Yongzhong Qian and Projjwal Banerjee for useful discussions relating to this work.

Support for this work was provided by National Aeronautics and Space Administration (NASA) through Einstein Postdoctoral Fellowship grants numbered PF3-140114 (LFR) and PF4-150122 (FF) awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. JL and CDO are partially supported by the National Science Foundation (NSF) under award nos. TCAN AST-1333520, CAREER PHY-1151197 and AST-1205732, and by the Sherman Fairchild Foundation. JCL is supported by NSF grant number AST-1313091. This work also benefitted from NSF support through award no. PHY-1430152 (Joint Institute for Nuclear Astrophysics Center for the Evolution of the Elements). MDD acknowledges support through NSF Grant PHY-1402916.