Nucleosynthesis in Outflows from Black Hole–Neutron Star Merger Disks with Full GR(ν)RMHD

Along with binary neutron star mergers, the inspiral and merger of a black hole and a neutron star is a predicted site of r-process nucleosynthesis and associated kilonovae. For the right mass ratio, very large amounts of neutron-rich material (relative to the dynamical ejecta) may become unbound from the post-merger accretion disk. We simulate a suite of four post-merger disks with three-dimensional general-relativistic magnetohydrodynamics with time-dependent Monte Carlo neutrino transport. We find that within 104 GM BH/c 3 (∼200–500 ms), the outflows from these disks are very close to the threshold conditions for robust r-process nucleosynthesis. For these conditions, the detailed properties of the outflow determine whether a full r-process can or cannot occur, implying that a wide range of observable phenomena is possible. We show that on average the disk outflow lanthanide fraction is suppressed relative to the solar isotopic pattern. In combination with the dynamical ejecta, these outflows imply a kilonova with both blue and red components.


INTRODUCTION
The merger of a black hole -neutron star (BH-NS) binary is a potential multi-messenger event.If the merger results in the tidal disruption of the NS, copious quantities of neutron-rich material undergoing rapid neutroncapture (or r-process) nucleosynthesis are ejected.The radioactive decay of freshly-synthesized heavy isotopes powers an electromagnetic transient called a kilonova.Very neutron-rich outflows are typically rich in lanthanides (140 ≤ A < 176), which are highly opaque to blue light, and thus produce a red kilonova.Less neutron-rich outflows synthesize a smaller mass-fraction of lanthanides, producing an optical or blue kilonova instead.
Whether the NS is tidally disrupted by the BH depends on the properties of the progenitor binary, such as the NS and BH masses, the BH spin, and the unknown NS equation of state (EOS).While gravitational waves from two BH-NS mergers have been confidently detected so far (GW200105 and GW200115), no kilonova counterpart has yet been found, which is not unexpected given the properties of the binary components (Abbott et al. 2021).Since many detections of BH-NS mergers are expected in the coming decades, modeling these systems with greater fidelity and predicting possible nucleosynthesis outcomes is of critical importance.
However, the composition of material ejected from BH-NS mergers remains uncertain.From theoretical modeling, the two main types of ejecta are the dynami-cal ejecta and post-merger disk outflows.A part of the disrupted NS material is ejected from the system due to tidal forces, referred to as the dynamical ejecta, while the rest settles into an accretion disk around the BH.The outflows from this post-merger accretion disk can constitute a significant fraction of the total merger ejecta and represent the focus of our work.Although a number of studies have investigated such (or similar) outflows, these ejecta are not yet fully understood (Just et al. 2015;Fernández et al. 2020;Fujibayashi et al. 2020aFujibayashi et al. ,b,c, 2022;;Wanajo et al. 2022).
The details of the r-process nucleosynthesis depend on the electron fraction (Y e ) of the material, which is set by neutrino-matter interactions.The properties of the post-merger disk driving the outflow are set by the interplay of magnetohydrodynamical effects and gravity.Thus, to obtain an accurate prediction of ejecta properties and resulting nucleosynthesis, we need general relativistic neutrino radiation magnetohydrodynamics (GRνRMHD) simulations.
There is a long history of treating post-merger outflows (Ruffert et al. 1997;Popham et al. 1999;Kohri et al. 2005) with various levels of approximation for angular momentum transport and neutrino physics, and the literature has recently exploded (Siegel & Metzger 2017;Siegel & Metzger 2018;Nouri et al. 2018;Christie et al. 2019;Fernández et al. 2019;De & Siegel 2021;Just et al. 2022;Murguia-Berthier et al. 2021;Fahlman & Fernández 2022;Nouri et al. 2022).An exciting recent development is the self-consistent evolution of a binary NS merger from inspiral through the remnant phase, seconds after merger in Hayashi et al. (2022).Recently, Foucart et al. (2020) and Foucart et al. (2022) performed merger calculations with Monte Carlo neutrino transport, however, they did not follow the post-merger disk.To-date, the only full-transport GRνRMHD simulations have been performed with our code, νbhlight (Miller et al. 2019a) in Miller et al. (2019b) and Miller et al. (2020).
Here, we present the first three-dimensional full transport GRνRMHD simulations of BH-accretion disk systems produced in BH-NS mergers.We study the conditions in four post-merger disks, spanning a range of BH and disk masses, and predict nucleosynthesis in the disk outflows.

THE METHODS
We use the publicly available code νbhlight (Miller et al. 2019a), which uses operator splitting to couple GRνRMHD via finite volume methods with constrained transport to Monte Carlo neutrino transport.We solve the equations of general relativistic ideal magnetohy-drodynamics, closed with the SFHo EOS, described in Steiner et al. (2013) and tabulated in O' Connor & Ott (2010).
Neutrinos can interact with matter via emission, absorption, or scattering.For emission and absorption, we use the charged and neutral current interactions as tabulated in Skinner et al. (2019) and summarized in Burrows et al. (2006).Neutrino scattering is implemented as described in Miller et al. (2019a).
We use a radially logarithmic, quasi-spherical grid in horizon penetrating coordinates with N r × N θ × N φ = 192 × 168 × 66 grid points with approximately 3 × 10 7 Monte Carlo packets.Although our code is Eulerian, we track approximately 1.5 × 10 6 Lagrangian fluid packets ("tracer particles") of which approximately 5 × 10 5 become gravitationally unbound.Our tracer particles are initialized within the disk so that they uniformly sample disk material by volume.
We run our simulations for approximately 10 4 GM BH /c 3 , which allows us to observe the disk in a quasistationary turbulent state.Depending on BH mass, this translates to a different amount of physical time, ranging from ∼ 200 -500ms.
To ensure adequate resolution, we utilize the MRI quality factor defined in Sano et al. (2004) and the radiation quality factor defined in Miller et al. (2019b).We find both to be adequate for all disks at all times.For more information on our code implementation and verification, see Miller et al. (2019a).

THE MODELS
We simulate four accretion disks representing possible systems formed in the aftermath of BH-NS mergers.In the case of high component spins, a very large fraction of mass can end up outside the BH in either tidal or dynamical ejecta (Krüger & Foucart 2020).While the current population of BHs measured by LIGO-Virgo indicates a small projection of component spins onto the angular momentum axis, it is important to explore what is possible to interpret future observations.We therefore also include models with high component spins, resulting in large disk masses.The properties of our initial BH-disk configurations are listed in Table 1.Note that the disks vary not just in disk mass but with respect to several initial properties, such as BH mass and spin.
To form the accretion disk, we begin with a constant entropy and constant Y e torus in hydrostatic equilibrium around a Kerr BH, as described by Fishbone & Moncrief (1976).Our torus starts with a single poloidal magnetic field loop with a minimum ratio of gas to magnetic pressure β of 100.As the system evolves, the magnetorotational instability (MRI, Balbus & Hawley 1991) self- consistently drives the disk to a turbulent state, which provides the turbulent viscosity necessary for the disk to accrete.

Disk Outflows
As the disk accretes, turbulent viscosity transfers angular momentum outward, carrying with it material.Some of this material will become unbound -the viscous outflow -due to energy injection by turbulent energy dissipation and nuclear recombination, without neutrino cooling.In-fall of material transforms gravitational potential energy into thermal, kinetic, and electromagnetic energy, which powers thermally-and magneticallydriven winds off the "surface" of the disk and through the corona.A magnetically powered jet clears the polar region and entrains some material, carrying it out the poles.
Our ejecta are comprised of material with Bernoulli parameter B e > 0 at a radius of roughly 250 gravitational radii.Most of the ejected material is moving at mildly relativistic speeds 0.1-0.2c.Due to computational cost, we did not run our simulations long enough to determine the total unbound mass from the disks.Thus, the ejecta masses obtained in the simulations are lower limits.
Figure 1 shows snapshots of the density and Y e of the four disks (left panels) and the evolution of the average Y e of the ejecta as a function of latitude and time (right panels).Typically, at early times, electron neutrinos get emitted near the equator and absorbed at relatively higher latitudes, driving the Y e in this region towards higher values.This imparts an angle-dependent structure to the Y e , with Y e values increasing with the angle off the equator, and this structure persists over time.At late times, the optical depths are lower, the neutrinos are basically free-streaming, and the Y e evolution is dominated by neutrino emission instead.
The properties of the disk outflow -primarily Y e , specific entropy s, and dynamical timescale τ dyn -set the r-process nucleosynthesis yields.Depending on s and τ dyn , ejecta are expected to be lanthanide-free for Y e 0.22 − 0.3 (Lippuner & Roberts 2015).We record the Y e and s of our ejecta, and compute τ dyn as follows: where we take the integral over the time it takes for the temperature of the material to drop from approximately 10GK to 1GK.
In Figure 2, we present the Y e , τ dyn , and s of the total ejecta, for all four disks.We can see that the Y e of virtually all of the ejected material across the four disks lies above ∼0.2 and the distribution peaks at values between ∼ 0.2 -0.4.The Y e distributions have a double peaked structure, with the peak at lower Y e corresponding to neutron-rich ejecta present in the midplane and the peak at higher Y e corresponding to mass ejection at higher latitudes, where Y e is driven up by neutrino absorption.We note a broadening of the distribution towards higher maximum Y e with increasing disk mass.
While the Y e varies meaningfully from one disk to another, s and τ dyn show little variation from disk to disk.The dynamical time is roughly 30ms for all four disks, with long tails.Most of the ejected material has entropies between 10-30 k B /baryon, with only the tails extending to different values for the different disks, with the exception of Disk 1, which has a significant high angle-high entropy outflow component.
With respect to nucleosynthesis, the Y e distribution in these outflows is both broad and marginal.For entropies and dynamical timescales typical for these simulations, Y e ∼ 0.24 is the lanthanide turn-off point i.e. above this value, the lanthanide mass-fraction drops below 10 −3 (Lippuner & Roberts 2015).These BH-NS disks have a substantial fraction of material close to this marginal value, consistent with a mildly suppressed r-process.

Nucleosynthesis
To compute r-process yields, the tracer particles representing the gravitationally unbound material are postprocessed using PRISM (Mumpower et al. 2017;Sprouse et al. 2021).The nucleosynthesis calculation for a tracer starts when its temperature falls below ∼10 GK and the initial composition is assumed to be the nuclear statisti-   cal equilbrium composition.The calculations are continued until 1 Gyr, assuming homologous expansion beyond the end time of a tracer.Nuclear data is implemented in our calculations following the prescription described in Zhu et al. (2018).
The final abundances computed with PRISM for all four disks are shown in Figure 3, showing a full range of r-process nuclei.However, the abundances of isotopes beyond the second r-process peak A ∼ 130 are suppressed to various extents relative to their solar values.Outflows from Disk 1 have the highest mass-fraction of elements beyond the second peak, including the lanthanides, while abundances of these isotopes are most suppressed for Disk 3, almost two orders of magnitude lower than solar values.Disk 4, which is the most massive disk, sits between these two abundance patterns.Abundances of Disk 2 outflows are very similar to those for Disk 3, with slightly higher values around the third peak.For Disks 3 and 4, we find enhanced production of nuclei with A ∼ 60-70 relative to Disks 1 and 2, associated mainly with the production of neutron-rich isotopes of Nickel and Zinc in the relatively high Y e ejecta component ( 0.4) present for these two disks.This suggests that it is possible to eject Zinc along with r-process elements for some fraction of BHNS mergers.Such an abundance signature is thus not unique to magnetorotational supernovae, which have been invoked to explain high [Zn/Fe] ratios seen in metal-poor stars (Nishimura et al. 2017;Tsujimoto & Nishimura 2018).
These abundance patterns are consistent with our expectations based on the Y e , s and τ distribution of the ejecta from the disks.We compute the lanthanide fractions in the ejecta, defined as the mass of lanthanides divided by the mass of all neutron-capture elements, and estimate the lanthanide fraction for the total ejecta (dynamical and disk) in Table 2.

Late Time
Accretion, and thus outflow, has not finished during our νbhlight simulations.By the final simulation time, 2.8%, 2.2%, 1.6% and 3.1% of the original disk mass has become gravitationally unbound for Disks 1-4, respectively.The fraction of the original disk mass that has accreted during the simulation is 37.5%, 27.7%, 41.5%, and 37.9% respectively.Dividing the unbound mass by the accreted mass, we find that the percentage of the accreted disk mass that has become unbound during the course of the simulation is 7.5%, 8.3%, 3.7% and 8.3% for Disks 1-4, respectively.We thus estimate that after a disk has fully accreted, a total of roughly 8% of the disk mass will become gravitationally unbound.However, other studies suggest values as high as 40% (Siegel & Metzger 2017;Fernández et al. 2019;Christie et al. 2019).We provide the estimated disk ejecta mass range in Table 2.
As the disk drains and density falls, neutrino absorption becomes subdominant to emission, and eventually the electron fraction in the disk will be driven towards its value at emission equilibrium (where the rate of electron fraction increase from the emission of electron antineutrinos is balanced exactly by the rate of decrease from the emission of electron neutrinos).As a lower bound on the electron fraction of the outflow that might become unbound after the simulation time ends, we report the spherically averaged, density-weighted emission equilibrium Y e (Miller et al. 2020), where Y em e is computed by solving for the Y e at a given density and temperature that equilibrates the emission rates of electron neutrinos and their antiparticles.
We find that Y e em is 0.39, 0.28, 0.24 and 0.21 for Disk 1, Disk 2, Disk 3 and Disk 4 respectively.This estimate indicates that perhaps the most massive disks may produce lanthanide-rich outflow at later times.However, we emphasize that the emission equilibrium value is a lower bound, and for the most massive disks in particular, absorption will matter longer due to high neutrino opacities, and likely raise Y e in the outflow.

SUMMARY AND DISCUSSION
In this paper, we present the first full transport GRνRMHD simulations of post-merger accretion disks resulting from BH-NS mergers.We compute nucleosynthesis yields in the disk outflows for four different BHdisk configurations.We find copious production of rprocess elements in these ejecta, however, the production of elements beyond the second peak is suppressed by up to a factor of 10-100 relative to solar abundances.This is consistent with the ejecta Y e being driven to values above the threshold for lanthanide production due to neutrino-matter interactions.
In Table 2, we present the total ejecta mass produced during the course of νbhlight simulations along with the lanthanide fraction X La in these ejecta.The lanthanide fraction distribution of metal-poor stars, which peaks at log X La ∼ −1.8, provides an observational test for the composition of merger ejecta.While Disk 1 produces X La ∼ 10 −2 , Disks 2, 3, and 4 produce X La ∼ a few times 10 −3 .Disk outflows from the latter three disks thus represent an interesting ejecta component distinct in their nucleosynthetic signature from dynamical ejecta, which are expected to be very neutron-rich and typically have lanthanide fractions of the order of 10 −2 .Combining the early disk ejecta with estimates of the late time disk ejecta and the dynamical ejecta gives us a sense of the overall nucleosynthesis in such mergers.While X La of ejecta from disks alone lies below the peak of the distribution observed in metal-poor stars, when tied with the dynamical ejecta, BH-NS mergers can produce lanthanide fractions close to that of the typical metal-poor star.
Given the sophisticated microphysics and neutrino transport in νbhlight, we are uniquely poised to comment on the nature of BH-NS kilonovae.Typically, for X La 3 × 10 −3 , the kilonova peaks in the near-infrared J band around 1 day (Even et al. 2020).However, the total X La is less relevant here compared to how ejecta composition varies with latitude as well as its overall morphology (Korobkin et al. 2021).Since the Y e of our outflows varies significantly with angle off of the midplane, the observed character of the kilonova will depend heavily on viewing angle.In general, the polar ejecta have higher abundances of isotopes below A ∼ 120 and lower abundances of isotopes at/between the second and third peaks, relative to the equatorial ejecta.This will shift the kilonova blue-ward due to a decrease in the lanthanide and actinide abundances.Thus, an early blue wind-produced kilonova may be visible if the remnant is viewed close to the polar axis.The observation of such a kilonova from a BH-NS merger will be a particularly exciting discovery since, unlike binary NS mergers, Note-Ejecta component masses and lanthanide fractions for all four disks, both computed from νbhlight simulations (columns 3 and 4) and estimated.Dynamical ejecta mass estimates (column 5) were provided by Francois Foucart (private communication) based on Kawaguchi et al. (2016) and Krüger & Foucart (2020).The lanthanide fraction of the dynamical ejecta is assumed to be solar, XLa∼ 0.04, using solar abundances from Arnould et al. (2007) (Ji et al. 2019).Total disk ejecta masses (column 6) are estimated to be 8-40% of the disk mass and the lanthanide fraction is assumed to be the same as that of the disk ejecta produced during the course of the simulation.X La, total (column 7) combines the total dynamical and disk ejecta estimates with their corresponding lanthanide fractions to estimate an overall lanthanide fraction for the total BH-NS ejecta.
any optical kilonova component can only be produced by the post-merger disk outflows.For BH-NS mergers, the dynamical ejecta lie on the equatorial plane (due to their tidal origin) and are very neutron rich, yielding a lanthanide-rich contribution to the kilonova (Shibata & Hotokezaka 2019).
Several long gamma ray bursts-namely GRB 060614 (Della Valle et al. 2006;Gal-Yam et al. 2006;Zhang et al. 2007), GRB 211221A (Rastinejad et al. 2022;Xiao et al. 2022;Yang et al. 2022), and GRB 211227A (Lü et al. 2022)-have been associated with a claimed kilonova afterglow (Zhu et al. 2022).We note that, even with the very large disk masses produced in our models, the reservoir of mass is insufficient to sustain accretion and non-trivial jet luminosity for the long durations required to match these observations.However, accretion might be sustained by fallback (Metzger et al. 2010), as suggested in Metzger et al. (2010), Desai et al. (2019), andZhu et al. (2022).Alternatively if the magnetic field is configured to magnetically arrest the disk, accretion may be further sustained.
In this work, we have demonstrated that a large diversity of outcomes is possible from BH-NS merger remnants.To fully connect this disk modeling to observations, more models covering a denser sampling of parameter space, simulations run to later times, and full radia-tive transfer models of the kilonova are needed.These studies will be the subject of future work.

Figure 1 .
Figure 1.Left: Azimuthally averaged density (left half) and Ye (right half) of each of the four disk models going from the least massive disk (top) to the most massive disk (bottom).The snapshots are taken at approximately 50ms.Right: Average Ye of gravitationally unbound material as a function of latitude and time.The angle off the equator is given by |90 • -θ bl |, where θ bl is the Boyer-Lindquist angle.Note the different end times for the four simulations.

Figure 2 .
Figure 2. From top to bottom: Histograms representing the electron fraction, mean dynamical timescale, and entropy of the ejecta for all four BH-disk systems.

Figure 3 .
Figure 3. Mass-weighted abundance as a function of mass number for outflows from all four accretion disks.Solar abundances from Arnould et al. (2007) are shown in gray and scaled to approximately match Disk 1 around the second peak.

Table 1 .
Foucart (2012)ionsNote-Disk label and corresponding remnant BH mass and spin, disk mass, inner radius of the disk, initial entropy and initial electron fraction.Disk masses are estimated based on the fitting formulae derived inFoucart (2012), assuming binary properties that favor NS disruption with the goal of producing a wide range of disk masses.The entropy and Ye for each disk were chosen based on the typical values obtained in dynamical BH-NS merger simulations (Francois Foucart, private communication).

Table 2 .
Ejecta masses (M ) and lanthanide fractions (-) Label Disk Mass Disk Ejecta X La,disk Estimated Dynamical Ejecta Estimated Disk Ejecta X