Constraining Inputs to Realistic Kilonova Simulations through Comparison to Observed r-process Abundances

Kilonovae, one source of electromagnetic emission associated with neutron star mergers, are powered by the decay of radioactive isotopes in the neutron-rich merger ejecta. Models for kilonova emission consistent with the electromagnetic counterpart to GW170817 predict characteristic abundance patterns, determined by the relative balance of different types of material in the outflow. Assuming that the observed source is prototypical, this inferred abundance pattern in turn must match r-process abundances deduced by other means, such as what is observed in the solar system. We report on analysis comparing the input mass-weighted elemental compositions adopted in our radiative transfer simulations to the mass fractions of elements in the Sun, as a practical prototype for the potentially universal abundance signature from neutron star mergers. We characterize the extent to which our parameter inference results depend on our assumed composition for the dynamical and wind ejecta and examine how the new results compare to previous work. We find that a dynamical ejecta composition calculated using the FRDM2012 nuclear mass and FRLDM fission models with extremely neutron-rich ejecta (Y e = 0.035) along with moderately neutron-rich (Y e = 0.27) wind ejecta composition yields a wind-to-dynamical mass ratio of M w /M d = 0.47, which best matches the observed AT2017gfo kilonova light curves while also producing the best-matching abundance of neutron capture elements in the solar system, though, allowing for systematics, the ratio may be as high as of order unity.


I. INTRODUCTION
For nearly half a century, neutron star binaries have been known to exist in nature, stemming from the first detection of a binary pulsar system [1].Shortly thereafter, the general relativistic prediction of gravitational radiation from a compact object binary was measured in the same system, implying the possibility of neutron star binary coalescence [2].Recently, neutron star mergers were confirmed as astrophysical sources of both gravitational wave and electromagnetic emission with the detection of the binary neutron star merger GW170817 and its transient electromagnetic counterpart AT2017gfo [3][4][5][6][7].
Around the same time as the first pulsar binary detection, compact object mergers involving neutron stars, either binary neutron star (BNS) or black-hole-neutronstar (BHNS), were predicted to be candidates for rapid neutron capture (r -process) nucleosynthesis [8][9][10][11].The nuclei synthesized in the immediate aftermath of the post-merger ejecta were thought to be heavy (A > 140), with a sizeable fraction of radioactive isotopes having short lifetimes due to their instability [12].As these nuclei decay, they release energy into the surrounding matter which would be emitted as ultraviolet, optical, and infrared thermal radiation once the ejecta becomes optically thin [13][14][15].This thermal emission is now commonly referred to as a kilonova [16] and serves as the bridge between the r -process elements produced by neutron star mergers and their resultant electromagnetic emission [17][18][19][20][21]. Aside from the transient electromagnetic kilonova emission (including a gamma-ray burst [22,23]), r -process material ejecta from neutron star binary mergers like GW170817 could produce another observable signature: relic r -process abundances such as in ancient, metal-poor stars and in our solar system.
Modeling kilonova light curves from merger events as viewed from the solar system by a distant observer requires the ejecta mass, velocity, composition, morphology, and viewing angle to be known, or otherwise supplied as model inputs.It has been conclusively demonstrated that the multiband light curves of AT2017gfo are poorly fit with single-component models, i.e., models consisting of a single type of ejecta described by fixed velocity, mass, and composition [24].Instead, the AT2017gfo light curve is better fit by two (or even three) components describing multiple types of ejecta: generally, a high-opacity "dynamical" component, and a low-opacity "wind" component [7,[24][25][26].
In two-component models of kilonovae, the low-opacity wind component typically includes elements only up to the "second r -process peak" at A ∼ 130, while the higher-opacity dynamical component includes even heavier elements [27][28][29].While these composition trends set a general opacity scale (see e.g., [30]), the full details of the composition effects on electromagnetic emission depend on the components' nuclear physics considerations as well as their physical parameterizations described, in part, by the components' masses M d , M w and velocities v d , v w , where "d" and "w" refer to the dynamical and wind components, respectively.Previous studies of kilonovae have highlighted the importance of nuclear physics inputs on r -process nucleosynthesis and the resultant effect on observed kilonova emission [31,32].In this work, we build on previous studies by considering the effect that nuclear physics uncertainties have on parameter inference from kilonova light curves.This work presents two-component nucleosynthetic yield constraints assuming r -process contribution exclusively from neutron star mergers and electromagnetic constraints assuming all neutron star mergers are phenomenologically similar to GW170817.We investigate the effects of comparing elemental abundances from kilonova simulations to solar r -process abundances under the assumption that the second (A ∼ 130) and third (A ∼ 195) r -process peaks follow universal behavior, which is justified by the robustness of the r -process pattern observed among metal-poor stars (e.g.[33][34][35][36][37][38]). We use this comparison to create a parameter estimation prior driven by explicit consideration of r -process elemental abundances in kilonova ejecta to gauge the effects on recovered ejecta properties; i.e., the masses and velocities of the ejecta components.As kilonova models improve in complexity and more observations become available for parameter estimation purposes, we can use more representative simulation abundances to hone this prior in future studies.
In this work, we will assess the extent to which our assumptions about the ejected material are simultaneously consistent with both types of aforementioned observations.Specifically, we will examine whether the abundances produced by our nucleosynthesis simulations realistically match the r -process abundances observed in the Sun while simultaneously reproducing the AT2017gfo light curve as well.In Section II A, we discuss the radiative transfer, atomic, and nuclear physics codes used to calculate the surrogate light curves, line-binned opacities, and ejecta compositions, respectively, considered in this work.In Section II B, we describe our method of comparing mass-weighted r -process abundances from our simulations with the solar abundance pattern.Section II C describes our parameter estimation framework and the effects of the r -process prior introduced in this work.In Section III we discuss whether the inclusion of the r -process prior makes a substantial difference in the parameter estimation process compared to prior work.
Our proof-of-concept analysis provides two key new approaches to multimessenger inference of BNS mergers.
On the one hand, we provide a method to quantitatively assess the hypothesis of a universal r -process origin in BNS mergers with observations of kilonovae, by requiring consistent predictions for the ejecta's electromagnetic emission and its asymptotic impact on r -process abundances.On the other hand, if the universal origin of r -process abundances is from binary neutron star mergers, our method can sharply refine our inferences about the ejected material.While in our proof-of-concept calculation we presently employ the Sun's r -process abundances as a prototype for a pristine, universal abundance signature from BNS mergers, our method should ideally be applied to ongoing and future efforts to disentangle the BNS merger's natal abundance contribution, for example from isolated metal-poor stars or from abundance principal-component correlations.

A. Simulation Setup
The aforementioned model abundances are not sufficient to create a direct link to kilonova electromagnetic emission on their own.However, they restrict which radioactive isotopes can plausibly exist and determine the radioactive heating rates powering the kilonova at different times.In this section, we describe the details of our ejecta compositions, the relevant thermalization efficiencies, and the composition-dependent ejecta opacity effects which constitute our kilonova emission model.Figure 1 schematically shows the subsequent process for using these models to perform kilonova parameter inference.
Throughout this work, we assume a two-component kilonova model composed of a lanthanide-rich "dynamical" ejecta component and a lanthanide-poor "wind" ejecta component.Each of our two ejecta components, dynamical and wind, is described by a fixed morphology and elemental composition.The morphologies are fixed to torus-shaped and peanut-shaped for the dynamical and wind ejecta, respectively (as defined in Ref. [39]).The wind component compositions, contributing to elements around and between the first (A ∼ 80) and second (A ∼ 130) r -process peaks, are fixed in this study and are described by the H5 and H1 tracers in Ref. [40] for the "wind1" and "wind2" models, respectively.We consider two different wind models with lower (wind1) and higher (wind2) neutron-richness to gauge the effects of lighter and heavier element contributions in our comparison to solar r -process abundances.The dynamical ejecta compositions, composed of the elements from the second to the third r -process peak and beyond, are varied as described in Table I; the dynamical ejecta composition used in our previous study, marked with the (*) label, is described by the model B tracer in [41].
We use the models from our previous study (see Ref. [42]) using SuperNu, a Monte Carlo code for simula- tion of time-dependent radiation transport with matter in local thermodynamic equilibrium [43].Our lightcurve simulations use radioactive power sources calculated from decaying the r -process composition from the WinNet code [44].The contributions from these power sources are weighted by thermalization efficiencies first presented in [45] (see Ref. [46] for a detailed description of the adopted nuclear heating).We utilize detailed opacity calculations from the tabulated binned opacities generated using the Los Alamos suite of atomic physics codes [47,48].Our tabulated binned opacities are not calculated for all elements; therefore, we produce opacities for representative proxy elements by combining pureelement opacities of nuclei with similar atomic properties [48].Our final SuperNu outputs are simulated kilonova spectra evaluated at 1024 equally log-spaced wavelength bins from 0.1 µm to 12.8 µm across 54 viewing angles spaced equally in cos θ for −1 ≤ cos θ ≤ 1.These spectra are then post-processed into light curves assuming a source distance of 10 parsecs.Our SuperNu simulations require discrete mass and average outflow velocity parameter inputs for the ejecta components.To sample our parameter space continuously during parameter inference, we require a continuous mapping of ejecta parameter inputs to kilonova light curve outputs in the form of a Gaussian process surrogate model.We built our surrogate model training li-brary of ∼ 450 kilonova light curve simulations using iterative simulation placement guided by Gaussian process variance minimization.In other words, we placed new simulations in regions of parameter space where our interpolation root-mean-square (RMS) uncertainty was largest.For simulation placement purposes, we only examined the uncertainty on the entire bolometric light curve rather than uncertainty at individual simulation times (see Ref. [42] for a full discussion on the creation of the simulation library).
Using Gaussian process regression interpolation in conjunction with our simulation library [42], we created a continuous mapping of our four ejecta parameter simulation inputs (M d , v d , M w , v w ) to a scalar output M AB,λ at some time t, angle θ, and wavelength λ.Because of the substantial dynamic range of our many outputs, we interpolate in AB magnitudes using the LSST grizy and 2MASS JHK bands as our reference wavelengths.Our Gaussian process uses a squared-exponential kernel and a white noise (diagonal) kernel.Unless otherwise noted, we quantify the performance of our interpolation with the RMS difference between our prediction and the true value.
Combining our surrogate light curves and parameter inference code, we generated posteriors for the ejecta parameters of GW170817 given our model assumptions [42].We perform four-dimensional Monte Carlo integration of the electromagnetic likelihood over our model's four parameters using the EM PE package1 to provide the likelihood and the RIFT adaptive Gaussian mixture model integrator to perform the integration [49].The parameter estimation in this work, discussed in Section II C, follows the same methodology as in Ref. [42], with the additional composition-based prior constraint described in Section II B.
In creating our surrogate model light-curve training library, we considered only one dynamical ejecta composition and one wind ejecta composition, indicated by the (*) label in Table I.With regard to constructing a composition-informed mass prior, the wind compositions in this work are unchanged from previous studies [39,42,46] while new considerations for dynamical ejecta compositions with different nuclear physics inputs were included in addition to the original dynamical composition used in prior studies.As a result of its ejection process, dynamical ejecta typically has a much lower electron fraction Y e ≡ (n p )/(n p + n n ), where n p is the number of protons and n n is the number of neutrons in the ejecta.A low electron fraction results in a higher availability of neutrons for capture during r -process nucleosynthesis and thus generally the creation of heavier elements such as lanthanides and actinides [50,51].Due to the dynamical ejecta's dominance on the elemental abundance pattern compared to the relatively minimal contribution of the wind ejecta, we only present newly calculated compositions for dynamical ejecta in this work.
The compositions presented in this work were generated using two nuclear network codes: WinNet and PRISM.The dynamical and wind ejecta compositions considered here and in previous work (i.e.[42,46]) were generated using WinNet.Specifically, these are the wind1 and wind2 compositions and the dynamical ejecta model using the Panov+ (2010) fission model in Table I.
The varied dynamical ejecta compositions new to this work were generated using PRISM.PRISM is a single-zone nuclear reaction network code that evolves an initial seed abundance of nuclei along a time-temperature-density thermodynamic trajectory, while allowing full flexibility with the input nuclear data [57].In this work, we use state-of-the-art nuclear reaction and decay rates that are calculated to be self-consistent with the nuclear mass model.Following from the thermodynamic trajectories of dynamical ejecta from neutron star merger simulations presented in Ref. [20], all of our PRISM runs begin in nuclear statistical equilibrium at a temperature of 10 GK in the thermodynamic trajectory.All dynamical ejecta models presented in Table I with a fission model different from Panov+ (2010) were generated using PRISM.
The mass fractions of all the compositions considered in this work are shown in Figure 2. Figure 2 highlights the main difference between our wind ejecta compositions; the wind1 composition has very low mass fractions at the second r-process peak (A ∼ 130) while a significant portion of the wind2 composition consists of elements around this peak.

B. Ejecta Prior Implied by r -process Observations
We seek to compare the combined mass fractions X sim to the seemingly universal pattern of elements between the 2nd and 3rd r -process peaks (the "main" r -process, [21]) observed among some of the oldest stars.This "rprocess universality" has been noted for iron-poor (or "metal-poor") stars that show enhancements in the main r -process elements relative to their iron content in excess of ten times the equivalent Solar ratio.However, observationally derived abundances in metal-poor stars are necessarily elemental since the abundances are derived from atomic transitions in stellar spectra that are overall insensitive to atomic mass number.Except for a handful of elements, the detailed isotopic distributions of r -process elements in metal-poor stars is observationally unknown.As a proxy for a representative example of the universal r -process, we use the well-studied solar isotopic abundance pattern X ⊙ , relying on previously-published projections of the high-A elements into different neutroncapture process contributions.Specifically, the r -process fractions presented in Ref. [58] are used in conjunction with the total abundances from Ref. [59] to isolate the contribution to the solar abundances by the r -process.
Figure 3 shows the [X/Fe] abundances of six metal-poor stars with r -process enhancements.The "[X/Fe]" notation means that each elemental ratio log ϵ(X/Fe) is compared to the same elemental ratio in the solar abundance pattern. 2 Stars with [X/Fe] > 0 are considered "enhanced" in that element relative to the solar system abundance.For many metal-poor stars, elements with Z ≥ 37 have an enhanced abundance compared to the Sun.For this work, guided by the enhancement seen in elements Z ≥ 37 in Figure 3 and with the assumption that iron was created during supernova nucleosynthesis, we assume that elements with Z ≥ 37 originate exclusively from neutron star mergers.The trends of elements with Z < 37 are less clear; they are not uniformly enhanced in stars that are enhanced with the Z ≥ 37 elements, likely pointing to multiple (non-merger) origins for these elements.
In our SuperNu simulations, we adopt a twocomponent compositional model and vary the mass ratio of the two components: the dynamical (M d ) and wind (M w ) ejecta masses.Each component has a fixed isotopic abundance, computed via nucleosynthesis network [44].Due to the fixed nature of the compositions, we weight each component's composition, represented by mass fractions X d and X w , by the mass of the respective ejecta component, dynamical M d and wind M w , to introduce composition variation as a function of component mass in the combined simulation mass fraction For every isotope, the combined mass fraction X sim is simply the mass-weighted sum of its mass fractions in the constitutive components.We varied our dynamical and wind component masses over a grid between −3 ≤ log(M d,w /M ⊙ ) ≤ −1, encompassing the most realistic ejecta masses predicted by numerical relativity simulations of neutron star mergers [29,40,[60][61][62][63].
To account for isotopes of actinides with short decay timescales, we rescale the solar mass fractions of actinides X ⊙,Ac to what they would have been at 1 day to better match the kilonova-timescale mass fractions used in our simulations.The rescaling is achieved by setting the 1day solar actinide mass fractions to values that would decay to present day values after 4.5 Gyr.The rescaled solar mass fractions are also mapped into a subset of representative elements used for SuperNu light-curve modeling as described in Section II A. Hereafter, any mention of the solar mass-fraction pattern X ⊙ refers to the 1day rescaled and mapped mass fractions using data from Refs.[58] and [59].
To get X ⊙ and X sim on the same relative scale, we introduce an offset C scale,Z that shifts X ⊙ down to comparable values for X sim by matching the two mass-fraction values at some element Z.To minimize how the choice of C scale,Z affects our results, we integrate over the range of Mass ratios were determined by calculating the mean mass ratio of the bottom 2nd percentile of all residuals so as to eliminate outliers.The residuals were calculated as in Equation 1 and the minimum residual was identified as the smallest residual across all the mass pairs considered for a given composition.The two wind1 and wind2 trajectories are described in detail in Ref. [46].The two nuclear mass models considered are FRDM2012 and HFB24 [52,53].The two nuclear fission models considered in our study are FRLDM [54] and "50/50," a simple symmetric assumption that fission yields split into two identical nuclei.5 of [56].We also consider the HFB27 mass model in place of HFB24.
possible values of C scale,Z introduced by scaling X ⊙ and X sim to matching values at different elements Z.This integration marginalizes over our uncertainty in C scale,Z , adopting a Gaussian prior on log(C scale,Z ) with mean µ = 0.08 and variance σ 2 = 0.695.After integrating over C scale,Z , we are left with a single required choice of a new C value that sets the constrained scale at some element Z choice such that log X sim,Z choice = log X ⊙,Z choice − log C. We chose Z choice = 46 as it is one of the elements present in all the dynamical, wind, and solar mass fractions.Shifting X ⊙ to be on the same relative scale as X sim and choosing a specific value of C at some Z choice are both done solely for the purpose of calculating well-behaved residuals.While the scaling by C removes the constraint that Z X ⊙,Z = 1, this has no detrimental effects on our analysis as we are interested exclusively in our composi-tions' relative abundances.
As part of our assumption that elements Z ≥ 37 were synthesized exclusively in neutron star mergers (see Figure 3), we only consider elements from Z = 37 up to and including Z = 103 when computing the residual r(M d , M w ) for all available elements in the solar abundance pattern Z ∈ Z ⊙ given a simulation with component masses M d , M w : FIG. 3. Elemental mass fraction ratios relative to iron of a sample of r -process enriched metal-poor stars.[X/Fe] > 0 implies enhanced abundance of element X compared to the solar system with respect to iron.We assume all elements that are significantly enhanced compared to iron to have been introduced post-supernova, exclusively from neutron star mergers.The region of enhanced elements (Z ≥ 37, highlighted in blue) is the focus of our comparison to solar composition.The iron-peak elements and supernova r -process are not strongly enhanced compared to solar (highlighted in yellow).Stellar elemental abundances obtained from JINAbase [64] with the respective stars reported in references [65][66][67][68][69][70].
where r(M d , M w ) is the residual for the given dynamical and wind mass pair used to calculate  4. Recreated Ye distribution for dynamical ejecta as presented in Fig. 5 of Ref. [56].The analytic fit to the distribution, presented in Equation 2, is overlaid as a red line.The analytic fit was used to generate the Ye-distribution compositions presented in Table II.
ement Z in both components (if present), σ is the uncertainty on log X sim , σ C is the uncertainty introduced by integrating out C scale,Z , X ⊙ is the average solar mass fraction across all elements, X sim is the average simulation mass fraction across all elements, and N is the total number of elements considered.
We compute the residual between each of our composition models X sim from Table I and the solar massfraction pattern X ⊙ using Equation 1.We consider component mass weights across a log-spaced grid with −3 ≤ log(M d,w /M ⊙ ) ≤ −1, with 50 masses for each component resulting in a total of 2500 residuals per composition model.The scaled residual values r − r min are shown in the top panel of Figures 5 and 6, with r representing the residual calculated for each mass pair and r min the lowest residual for all mass pairs considered for a given model.
Guided by numerical-relativity simulations which suggest a distribution of Y e values in neutron-star ejecta (e.g., [29,56,63,71]), we also analyze a selection of compositions derived from the Y e distribution for dynamical ejecta presented in Figure 5 of Ref. [56].To calculate the mass-weights for each Y e value, we recreate the Y e distribution with a piece-wise analytic fit ∆M/M = 0.05C(Y e /0.04) −1.7 e −0.5((Ye−0.17I).Bottom: Mass fractions of individual ejecta components compared to the best-fit mass fraction Xsim obtained from comparison to X⊙.The line colors represent the same quantities as in Figure 5.The minimum residual was identified as the smallest residual across all the mass pairs considered for a given composition Xsim.
to deviate from the distribution at Y e = 0.3; this is of little concern as less than 3% of the total ejecta mass is described by Y e > 0.3.
For each Y e -distribution composition, we run 10 nucle-osynthesis simulations, evenly spaced between 0.035 ≤ Y e ≤ 0.35.The mass-weight for each single-Y e composition is calculated using Equation 2. Once all the weights are calculated, they are normalized such that their net contribution describes the total ejecta mass.The final Y e -distribution composition is the weighted sum of the abundances from the single-Y e nucleosynthesis simulations.After the net Y e -distribution composition is calculated, we repeat the same methodology as for the models presented in Table I to calculate mass ratios and residuals for each composition model.

C. Parameter Inference
As in many previous applications of Bayesian inference to infer parameters of kilonovae [26,[72][73][74][75][76][77][78], we seek to compare the observed magnitudes x i at evaluation points i (denoting a combination of band and time) to a continuous model that makes predictions m(i|θ) which depend on some model parameters θ.Bayes theorem expresses the posterior probability p(θ) in terms of a prior probability p prior (θ) for the model parameters θ and a likelihood L(θ) of all observations, given the model parameters, as Unless otherwise noted, for simplicity we assume that the source sky location, distance, and merger time are known.We adopt a uniform prior on the ejecta velocity v/c ∈ [0.05, 0.3] and the two-dimensional prior discussed in Section II B on the ejecta masses m/M ⊙ ∈ [0.001, 0.1].
We assume the observations have Gaussian-distributed magnitude errors with presumed known observational (statistical) uncertainties σ i , convolved with some additional unknown systematic uncertainty σ, so that our log-likelihood is (4) where the sum is taken over every data point in every band used in the analysis.For inference using our Gaussian process surrogate models, we set σ to the estimated Gaussian process model error.For a full discussion of our parameter inference considerations, see Ref. [42].

III. RESULTS
For our two-component models, assuming a single source like GW170817 dominates the observed solar rprocess abundances, the inferred abundances from such mergers only depend on the mass ratio M w /M d .In other words, since in our study we emphasize only the relative and not absolute r -process abundances, motivated by considerable uncertainty in the binary neutron star merger rate, we therefore only use and constrain the abundance ratios.The relative abundances from a single channel only depend on the relative proportions of this channel; for our two-component model, this is simply dependent on M w /M d .Thus for each set of initial assumptions-the composition of the dynamical ejecta (represented by the electron fraction Y e ), the presumed nuclear mass and fission model, and other details-our comparison with solar abundances necessarily constrains M w /M d narrowly around a preferred value unique to that model.We note that the abundances we are considering are effectively frozen out for the processes we're interested in at times later than O(1) second.
Tables I and II provide a list of models and their preferred M w /M d , in the sense that they minimize the residual mismatch with the solar abundances as calculated in Equation 1.With a few exceptions, most models prefer M w /M d substantially lower than order unity.In other words, most of our abundance comparisons suggest that less wind than dynamical ejecta would be required for GW170817-like mergers to reproduce the solar r -process abundances.These results are at odds with those found from other contemporary modeling [74,77,79,80] as well as numerical relativity results [40,63,[81][82][83], which typically predict more post-merger (i.e.wind) mass ejection.
One particularly interesting result is that the lowestresidual single-Y e composition in Table I agrees with X ⊙ better than the Y e -distribution composition with similar nuclear physics inputs.This result suggests that a single-Y e approximation is adequate for computational simplicity in the context of nucleosynthesis calculations.
The top panel of Figure 5 shows the mass-pair residuals for the composition and morphology assumptions considered in previous work [42], denoted by the (*) label in Table I.The yellow stripe indicating the lowest residual region highlights the best-fitting ratio of windto-dynamical mass implied by the calculated residuals.This corresponds to the "Mass Ratio (M w /M d )" value recorded in Table I.
The top panel of Figure 6 shows the best-fitting windto-dynamical mass ratio for the lowest residual composition model presented in this work: dynamical ejecta with a composition characterized by the FRDM 2012 mass model, FRLDM fission model, electron fraction Y e = 0.035 and wind ejecta corresponding to the wind2 model.The bottom panels of Figures 5 and 6 show the solar mass-fraction pattern X ⊙ in gray, the dynamical and wind ejecta fixed composition mass-fraction patterns X d and X w in red and blue, respectively, and the 90% confidence interval for the mass-fraction pattern of the relevant composition model X sim as the light-blue shaded region.The 90% confidence interval was calculated for the spread of possible mass-fraction patterns which arose when scaling the fixed component compositions X d (X w ) by the component mass M d (M w ).Based on the good agreement shown in Figure 6, we adopt the associated two-dimensional likelihood versus M d and M w as a prior constraint on ejecta masses.In other words, we assume   6 with no electromagnetic information provided during sampling besides constraints on the opening angle (see Ref. [84] and references therein).We apply a scale factor in our likelihood calculation during inference to prevent underflow from the large residual values.Note the recovery of the yellow band of lowest-residual mass pairs from Figure 6 in the Mw vs. M d panel as well as the flat velocity posteriors stemming from the lack of velocity constraints introduced by our mass-focused prior.Residual small-scale substructure in the one-and two-dimensional marginal distributions reflects sample size artifacts.
GW170817 was produced from a representative member of a single population of kilonovae, which alone are responsible for the solar r-process abundance).
For each set of initial assumptions, the inferred constraint on M w /M d therefore also strongly constrains the ingredients powering the associated kilonova.For example, Figure 7 shows the results of inferring the parameters of GW170817, using only our prior constraints on M w /M d from the top panel of Figure 6 (and weak constraints on the binary orientation relative to our line of sight).Figure 8 shows how these constraints propagate into joint electromagnetic inference.The solid black contours show inferences derived without using constraints on M w /M d ; the red contours show inferences supplemented with this insight, for a specific set of initial assumptions.Figure 9 shows the light curves associated with the recovered posterior distributions presented in Figure 8.The inclusion of the composition prior results in much tighter model uncertainties compared to the lightcurve fits in Ref. [42].
Each set of our input assumptions about ejecta composition and physics makes a prediction about r -process Posterior distributions for samples generated when using the grizyJHK bands considered in [42] (black) and samples generated using the same bands along with the rprocess prior from Figure 6 (red).The values reported at the top of each posterior distribution represent the inference results from this study.The composition prior effect is most evident in the wind mass posterior shift to lower ejecta mass.6) reduces model prediction uncertainty compared to previous predictions [42].
abundances.As shown by the last column in Tables I and II, some of our input assumptions fit better than others.Given substantial systematic uncertainties associated in the many assumptions in our study, we approach these nominal residuals with considerable cautions.However, the minimum residuals presented in Tables I and II suggest that the wind2 model is a notably better fit to the solar mass-fraction pattern, consistent with similar findings in previous studies [85].The nearly distinct separation of the two wind models' lowest residuals implies that the wind1 model is less indicative of r -process nucleosynthesis from neutron star mergers; however, our work neglects consideration of lighter r-process elements (Z ≤ 37) which disfavors compositions with higher Y e like wind1.More importantly, the separation between the models also implies that new models for the wind ejecta composition need to be considered in comparison to the wind2 model.The results of Tables I and II indicate the need for further studies involving updated wind ejecta composition modeling informed by GRMHD disk simulations [29].
The results of Tables I and II also depend strongly on the assumption that neutron star mergers are the dominant r -process mechanism for the creation of elements with Z ≥ 37 (see Figure 3 and Section II B), which may not be the case; see, e.g., Ref. [86] and references therein.
Our method as stated also assumes that a narrow distribution of mergers in total mass M tot = M 1 + M 2 dominates nucleosynthesis yields.While self-evidently consistent with the binary neutron star population inferred from the merging Galactic neutron star binaries, this assumption could even still hold for a wider binary neutron star population as suggested by gravitational wave observations, if ejecta are (as expected) suppressed for the most massive mergers with large M tot .
Another caveat that presents limitations to our results is that we only incorporate very specific wind1 and wind2 compositions.There can be a broad variety of compositions permitted for electron fractions Y e > 0.20 due to varying hydrodynamic conditions.An extensitve study of these compositions, along with the tests of how much they can be considered "representative" of their respective components, is beyond the scope of this work.
Our results can further be improved by incorporating the observed higher variability of the lighter r -process abundances between the first and the second peak, compared to the universal pattern between the second and third r -process peaks.The lighter r -process as observed in metal-poor stars, exhibits variation on the order of 1 dex, while the "strong" r -process pattern varies by only about 0.3 dex [21].An investigation with more accurate numbers based on careful statistical analysis of observations will be the subject of future studies [87].
Finally, we note that we cannot ignore the bias introduced by the dynamical ejecta composition of the surrogate kilonova light curves presented in [42].While the constraints imposed by the r -process abundance prior indeed shift the recovered parameters as in Figure 8, there still remains some contribution to the parameter estimation stemming from the surrogate models having been trained on a different dynamical ejecta composition.In other words, our surrogate light curves were trained using the ejecta compositions in Table I labeled with (*).Although the primary contribution to the parameter inference in this work comes from the prior discussed in Section II B, some bias from the surrogates' original compositions is unavoidable.

IV. CONCLUSION
We have presented an approach for incorporating nuclear physics-based composition effects as a prior for our kilonova parameter inference framework.Identifying a self-consistent electromagnetic and isotopic signature from kilonova models enables us to make sharper conclusions about any specific kilonova's ejecta.Moreover, our calculations provide a Bayesian evidence, assessing how well both observations can be fit independently and together.Our approach may therefore provide a new avenue to directly test whether BNS mergers are the primary source of r -process enrichment, using any potential star as a prototype (e.g., the Sun or a metal-poor, r -process enriched example).
While our self-consistent approach will remain impactful as pertinent inputs improve, of course the quantitative numerical results shown here are merely illustrative, given substantial systematics.Our calculations rely on physical models of BNS mergers, ejecta, and the r -process abundance signature in solar and metalpoor stars, all of which have substantial and widelyinvestigated systematics.To highlight one important potential source of systematics in this approach, we considered a range of models with varying nuclear physics inputs, and, given the assumptions discussed above, the best-fitting model appears to be the one with FRDM2012 nuclear mass, FRLDM fission, Y e = 0.035 (extreme neutron richness in dynamical ejecta), moderately neutronrich wind ejecta, producing inferred dynamical and wind ejecta masses of M d ∼ 0.021 and M w ∼ 0.012 and corresponding to a relatively low mass ratio: M w /M d = 0.47.Our consideration of additional dynamical ejecta compositions, when compared to solar abundances, has indicated that the mass ratio between the two ejecta components is larger than what was implied by previous inference (M w /M d = 0.24).However, our conclusions should be taken with care, since the number of input compositions considered were quite limited.
We have also shown that the inferred mass ratio stemming from a comparison of r-process elemental abundances is highly sensitive to the input nuclear physics.For our preferred wind2 model, variations in the dynamical ejecta composition can change the recovered windto-dynamical ejecta mass ratio by a factor of ∼ 4.3.For the wind1 model, the inferred mass ratio can change by a factor of ∼ 4689, although this is largely due to the assumptions made in this work.While we focused on selected nuclear physics systematics, we have identified several areas meriting further study, including propagating different choices for nuclear physics uncertainties into our parameter inferences; examining whether more complex composition (e.g., Y e ), angular, or velocity distributions in the outflow can mimic these effects; and employing better prototypes than the Sun for a potential r -process signature.
Even allowing for extremely conservative systematic uncertainties on our inputs (e.g., assuming M w /M d 's optimal value is well-localized between 0.1 and 10), these prior abundance constraints should still provide useful insight into kilonova ejecta modeling.For instance, this framework of kilonova surrogates with abundance priors can be used as a constraint to identify merger simulations that produce consistent properties.
With the introduction of this composition-based prior, we are able to continue using our existing kilonova surrogate model and parameter inference frameworks while updating our inference priors to match contemporary results in the literature.The ability to update our mass prior using the underlying properties of kilonova models, without requiring expensive simulations (outside of nuclear network outputs), allows us to inexpensively and rapidly update our parameter inference results.
In this work, we have considered fiducial initial conditions for the outflow, including composition, without allowing for correlations induced by the fact that both the composition and outflows are initialized by binary neutron star mergers.In future work, we will explore self-consistent initialization from merger properties, in particular exploring the effects of binary mass ratio and neutron star remnant lifetime, which should have significant impact on the ejecta amount and composition [88][89][90].III.Sensitivity study results using a leave-one-out approach with specific elements, the lanthanides, the actinides, and both lanthanides and actinides.In each case, the elements in the Z column are removed from consideration during the residual calculation.Interesting cases include Z = 86 which significantly reduces the calculated residual and yields the largest ratio of Mw/M d , as well as Z = 100 which yields the lowest ratio.We find that we are particularly sensitive to the actinides in our compositions, as the mass ratio dramatically changes when they are removed from the residual calculation.

Appendix A: Sensitivity Study
Due to the variability in mass ratios as nuclear physics assumptions change (see Tables I and II), we conduct a sensitivity study to identify if certain elements have significant impact on the recovered mass ratio.We modify our best-fitting case with mass model FRDM2012, FRLDM fission, a fixed-Y e of 0.035 and a wind2 composition by choosing specific elements to be omitted during the residual calculation.The results are shown in Table III.Most notably, the inclusion (or lack thereof) of actinides strongly influences the preferred mass ratio, dropping it from an average of M w /M d ≈ 0.3 to M w /M d ≈ 0.1.We also note elements of interest Z = 86 and Z = 100 which, when removed during the residual calculation, yield the highest and lowest mass ratios, respectively.Our sensitivity to the actinides indicates that our choice of mass ratio during inference should be used with caution.Specifically, the mass prior can provide additional parameter constraints when used selfconsistently with the input model nuclear physics, but can have detrimental effects when applied incorrectly.Despite some sensitivity to the actinide composition, our central conclusion that M w /M d < 1 remains robust.

FIG. 2 .
FIG.2.Unscaled mass-fractions X as a function of atomic number Z for all single-Ye compositions considered in this work.The labels pertaining to the dynamical ejecta compositions considered in this work indicate the nuclear mass model, fission model, and electron fraction Ye used to generate the respective composition.The remaining labels indicate the solar and wind compositions.The wind1 and wind2 compositions do not extend to higher atomic numbers Z due to their higher electron fractions Ye = 0.37 and Ye = 0.27, respectively.

FIG. 5 .
FIG.5.Top: 2-D distribution of residuals calculated as in Equation 1 by comparing log X⊙ to log Xsim using the dynamical and wind compositions matching those in[42], represented by the (*) label in TableI.The residual grid is composed of 50 mass values equally log-spaced between −3 ≤ log M⊙ ≤ −1 for both dynamical and wind mass.The black lines indicate contours of equally-spaced residual values.The dashed red lines indicate a wind-to-mass ratio of 1 and serve purely as a visual aid.Bottom: Mass fractions of individual ejecta components compared to the best-fit mass fraction Xsim obtained from comparison to X⊙.The red and blue lines are the initial unweighted dynamical (X d ) and wind (Xw) ejecta mass fractions, respectively, scaled by C to match the solar mass fraction at Z = 46.The gray line is the solar mass fraction X⊙ and the blue shaded region is the 90% confidence interval for all the mass-weighted mass fractions log Xsim.The dynamical ejecta mass fraction only exceeds log X = 0 due to the scale matching at Z = 46.

FIG. 6 .
FIG.6.Top: Same 2-D distribution as described in Figure5except with the compositions which yielded the lowest residual in comparison to the solar abundance pattern X⊙ (top row of TableI).Bottom: Mass fractions of individual ejecta components compared to the best-fit mass fraction Xsim obtained from comparison to X⊙.The line colors represent the same quantities as in Figure5.The minimum residual was identified as the smallest residual across all the mass pairs considered for a given composition Xsim.
FIG.8.Posterior distributions for samples generated when using the grizyJHK bands considered in[42] (black) and samples generated using the same bands along with the rprocess prior from Figure6(red).The values reported at the top of each posterior distribution represent the inference results from this study.The composition prior effect is most evident in the wind mass posterior shift to lower ejecta mass.

FIG. 9 .
FIG.9.Broadband light-curve predictions for the ejecta parameters recovered in Figure8.The inclusion of the composition-based mass prior (top panel of Figure6) reduces model prediction uncertainty compared to previous predictions[42].

TABLE I .
Wind-to-dynamical mass ratios sorted by increasing minimum residuals for each dynamical composition considered.

TABLE II .
[55]fission rates for the simulations performed in previous work, labeledPanov+ (2010), are taken from Ref.[55].The (*) indicates the compositions used in creating the surrogate light curves used during parameter estimation (see Section III).The reported Ye values describe the neutron-richness of just the dynamical ejecta component.Same as TableIexcept considering dynamical ejecta compositions derived from the Ye distribution presented in Figure V. ACKNOWLEDGMENTS ROS, MR, and EMH acknowledge support from NSF AST 1909534.ROS acknowledges support from NSF AST 2206321.EMH acknowledges additional support for this work by NASA through the NASA Hubble Fellowship grant HST-HF2-51481.001, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.CJF, CLF, MRM, OK, RW, and TMS were supported by the US Department of Energy through the Los Alamos National Laboratory.Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).Research presented in this article was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20190021DR.This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001.