Markov Chain Monte Carlo Predictions of Neutron-rich Lanthanide Properties as a Probe of r -process Dynamics
Nicole Vassh, Gail C. McLaughlin, Matthew R. Mumpower, Rebecca Surman
DDraft version June 9, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Markov Chain Monte Carlo Predictions of Neutron-rich Lathanide Propertiesas a Probe of r -process Dynamics Nicole Vassh, Gail C. McLaughlin, Matthew R. Mumpower, and Rebecca Surman Department of Physics, University of Notre Dame, Notre Dame, IN 46556, USA Department of Physics, North Carolina State University, Raleigh, NC 27695, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA (Dated: June 9, 2020)
ABSTRACTLanthanide element signatures are key to understanding many astrophysical observables from mergerkilonova light curves to stellar and solar abundances. To learn about the lanthanide element synthesiswhich enriched our solar system, we apply the statistical method of Markov Chain Monte Carlo toexamine the nuclear masses capable of forming the r -process rare-earth abundance peak. We describethe physical constraints we implement with this statistical approach and demonstrate the use of theparallel chains method to explore the multi-dimensional parameter space. We apply our procedure tothree moderately neutron-rich astrophysical outflows with distinct types of r -process dynamics. Weshow that the mass solutions found are dependent on outflow conditions and are related to the r -processpath. We describe in detail the mechanism behind peak formation in each case. We then compare ourmass predictions for neutron-rich neodymium and samarium isotopes to the latest experimental datafrom the CPT at CARIBU. We find our mass predictions given outflows which undergo an extended(n, γ ) (cid:29) ( γ ,n) equilibrium to be those most compatible with both observational solar abundances andneutron-rich mass measurements. INTRODUCTIONUnderstanding the nucleosynthesis of lanthanides, at57 ≤ Z ≤
71, is important in order to interpret many as-trophysical observables such as the abundances of metal-poor stars (Sneden et al. 2008; Frebel 2018) and mergerkilonova light curves (Kasen et al. 2017; Barnes et al.2016). In order to use lanthanide signatures to probe theorigins of heavy element production in our solar system,it is crucial to consider abundances derived from nu-cleosynthesis calculations. Such calculations connect toproduction sites by considering the possible outflow con-ditions present in an astrophysical environment. Theseefforts are challenged by the uncertain properties of theneutron-rich nuclei synthesized during the rapid neu-tron capture process ( r process). Such nuclear physicsuncertainties generate large ranges in the r -process nu-cleosynthetic yields of key lanthanide elements such asEu ( Z = 63) (Cˆot´e et al. 2018). Therefore gathering Corresponding author: Nicole [email protected] further observational information is alone insufficient todevelop a comprehensive picture of lanthanide produc-tion since nuclear physics advancements are also needed.In this work we approach the uncertainties surround-ing the synthesis of lanthanide elements in astrophysicalenvironments by taking advantage of statistical meth-ods which consider both nuclear and observational data.Methods to gather data which will illuminate previouslyunprobed physics are making notable advancements.Nuclear physics serves as a testament of such advance-ment with facilities such as CARIBU at Argonne Na-tional Laboratory, RIKEN, and the upcoming Facilityfor Rare Isotope Beams (FRIB) expanding the precisionand reach of studied properties.To take full advantage of current and anticipated data,Bayesian and Monte Carlo statistical methods are ex-panding their influence in nuclear physics theory. Usingthese statistical methods to treat complex processes hasgained traction for a wide variety of applications (Robert1998; Berg 2004; Brooks et al. 2011). In nuclear physicstheory, these methods have been applied to determineproperties of superhadronic matter from heavy-ion colli- a r X i v : . [ nu c l - t h ] J un Vassh et al. sions (Sangaline & Pratt 2016), nuclear emission spectra(Gulam Razul et al. 2003; Barat et al. 2007), estimatethermonuclear reaction rates (Iliadis et al. 2016) andreaction rate uncertainties (deBoer et al. 2014). Quan-tum Monte Carlo methods can be coupled with effectivefield theory to examine a variety of nuclear structureand interaction properties of light nuclei (Carlson et al.2015) such as explorations of the interactions at playduring weak transitions (Pastore et al. 2018) as well asenergy levels and level ordering (Piarulli et al. 2018).Bayesian techniques have also been applied to use massdata as well as mass model predictions to extrapolateproperties of exotic nuclei such as the masses of neutron-rich species far from stability (Utama & Piekarewicz2017; Neufcourt et al. 2018) even at the neutron dripline(Neufcourt et al. 2019). Approaches such as these canprovide valuable insights to astrophysics since they caninform nucleosynthesis calculations.In addition to wealths of current and advancing ex-perimental data, there exists opportunities to take ad-vantage of observational data. For instance telescopedata has been coupled with Markov Chain Monte Carlo(MCMC) methods to infer gas circulation processesin dwarf galaxies from stellar abundances (Cˆot´e et al.2017), to derive cosmological parameters from CMBanisotropies (Planck Collaboration et al. 2018), and toinfer the probability that Type Ia Supernovae will occurin a population of stars (Strolger et al. 2020). Bayesianmethods have also been applied in examining the neu-tron star crustal composition (Utama et al. 2016) as wellas the high density equation of state applicable to neu-tron stars (Drischler et al. 2020).In the era of multi-messenger astrophysics there ex-ist some striking testaments to how observational datacan be coupled with statistical methods. Statistical ap-proaches can been used alongside LIGO/VIRGO gravi-tational wave data to learn about the neutron star tidaldeformability and mass-radius relation (Abbott et al.2018; Capano et al. 2020; Miller et al. 2020) whichare connected to the equation of state of dense mat-ter. Bayesian and MCMC methods have also been cou-pled with X-ray observations (Steiner et al. 2010; N¨attil¨aet al. 2016; Steiner et al. 2018; Goodwin et al. 2019;Raaijmakers et al. 2019; Miller et al. 2019; Riley et al.2019) to learn about the mass, radius, thermal emission,and other properties of neutron stars. A global analysiswhich takes into account both gravitational wave andX-ray data can also be performed (Zimmerman et al.2020). Further opportunities for the applications of sta-tistical methods presented by multi-messenger scienceare exemplified by studies of the GW170817 gravita-tional wave event coupled with its electromagnetic coun- terpart (Abbott et al. 2017; Abbott et al. 2017; Cow-perthwaite et al. 2017; Villar et al. 2017). For instance,MCMC methods have been applied to model the ejec-tion dynamics of this event, such as the jet structure(Wu & MacFadyen 2018). MCMC methods were alsoused in some cases to model the kilonova light curve(Cowperthwaite et al. 2017) which indicated the syn-thesis of r -process elements through the signatures ofhigh opacity lanthanides.Since interpretations of r -process observables are af-fected by the unknown properties of neutron-rich nuclei,here we apply statistical methods to invert the problemby using solar data to find the nuclear masses which areable to reproduce lanthanide abundances. We focus onthe A ∼
164 enhancement of the solar r -process residu-als referred to as the rare-earth abundance peak. Notethat although the rare-earth elements include all lan-thanides as well as scandium ( Z = 21) and yttrium ( Z =39) here we only include the lanthanides in our analy-sis. Specifically, we employ the statistical techniques ofthe Metropolis-Hastings algorithm and MCMC to con-nect nuclear properties and rare-earth abundances, asfirst introduced in Mumpower et al. (2016b, 2017). Herewe present how we have evolved this “reverse engineer-ing” approach and describe the procedure we find tobe most suited to our particular MCMC problem. Sincethe mass predictions given by this method are tied to theastrophysical environment considered, we obtain resultsgiven several distinct outflow conditions. By consideringdistinct outflows, and cross-checking against recent Pen-ning trap mass measurement data for neutron rich nucleifound by the CPT at CARIBU (Van Schelt et al. 2012,2013; Orford et al. 2018), this method has the potentialto implicate the dominant outflow conditions respon-sible for the production of rare-earth nuclei observedin our solar system. This work therefore approachesthe problem of the uncertain outflow conditions withinastrophysical environments from a nuclear physics per-spective by examining which outflow dynamics are bothconsistent with the latest experimental measurementsand able to produce not only some lanthanide elements,but the proper elemental and isotopic ratios we observein our Sun and many other stars.This paper is organized as follows: in Section 2 we de-scribe why rare-earth abundances are of particular in-terest, their connection to nuclear physics properties,and the current state of rare-earth element abundancepredictions. In Section 3, we describe the MCMC ap-proach implemented to explore neutron-rich masses andshow the diagnostic criterion used to gain confidence inour results. Section 4 describes the three astrophysi-cal outflow conditions considered here and their distinct CMC lanthanide properties and r -process dynamics THE UNCERTAIN ORIGIN OF THE R -PROCESSRARE-EARTH ABUNDANCE PEAKThere exists a well-established connection between thesecond and third peaks of the solar r -process abundanceresiduals inferred from observational data at A ∼ A ∼
195 and the magic numbers at N = 82 and N = 126. Since nuclear states at magic numbers havean enhanced stability, during nucleosynthesis r -processproduction sees a ‘pile-up’ of such nuclei, i.e. a ten-dency of nuclei to capture into such states but then waitfor longer timescales to β -decay or capture out of thesestates. A question then naturally arises as to whetherother features in r -process abundances can be linked topile-up. For instance, the origin of the rare-earth peak,the subtle enhancement in lanthanide abundances be-tween 150 (cid:46) A (cid:46)
180 with its peak at A ∼ r process depends on the astrophysical out-flow conditions (Mumpower et al. 2012, 2016b, 2017),the rare-earth peak has the diagnostic power to shedlight on the nature of outflows from r -process sites. Forthis dynamical mechanism to operate it is only requiredthat the synthesis produces a main r process, reachingthe third peak and not necessarily beyond. Should theastrophysical outflows be either of high enough entropyor neutron-rich enough for synthesis to proceed past thethird r -process peak and into the actinide region, fissiondeposition could significantly influence rare-earth abun-dances. However since local nuclear features could alsoaffect how fission daughter products settle into their fi-nal abundances, neutron-rich rare-earth properties are of relevance in a wide range of outflow conditions. Thusto focus on the influence of solely local nuclear propertiesin the formation of the rare-earth peak, in this work weconsider astrophysical outflows which produce a main r process but do not significantly populate fissioning nu-clei. The investigation of outflow conditions which seeabundances impacted by fission daughter products willbe considered in future work.The local nuclear features which can dynamically formthe peak could be reflected in the nuclear masses whichultimately determine the nucleosynthetic outcome viatheir influence on reaction and decay rates. We there-fore focus on nuclear masses in our aim to understandthe properties of rare-earth nuclei needed to accommo-date rare-earth peak abundances. In Figure 1 we showthe calculated abundances of the rare-earth peak usingfour mass models commonly considered in r -process cal-culations: Duflo-Zuker (DZ, Duflo & Zuker 1995), theFinite Range Droplet Model (FRDM1995, M¨oller et al.1997; FRDM2012, M¨oller et al. 2016), and Hartree-Fock-Bogoliubov (HFB-21, Goriely et al. 2010). Thedistinct dynamics of the two types of outflow conditionsconsidered in Figure 1, which are characterized in termsof their entropy ( s ), expansion timescale ( τ ), and elec-tron fraction ( Y e ), are described in detail in Section 4.Although the application of some models can producea rough rare-earth peak, they often do not produce astrong enough enhancement, leave the rare-earth peak inthe wrong position, and/or miss the more subtle smoothbehavior between neighboring isotopic abundances. Wenote that the masses predicted by FRDM1995 are mostsuccessful in producing a peak in roughly the correct lo-cation for both types of outflows due to a kink in oneneutron separation energies at N = 104 which createspile-up (Surman et al. 1997). FRDM2012 mass predic-tions see a significant reduction in this feature predictedby FRDM1995 and therefore not as robust of a peak isproduced with the updated model. The DZ model showsa case where predicted deformation in the lanthanides isentirely absent, as evidenced by rare-earth abundanceswhich are on average flat. Fortunately precision nu-clear physics measurements are pushing into the lan-thanide region relevant for dynamical peak formation(Van Schelt et al. 2012, 2013; Orford et al. 2018; Vilenet al. 2018) with mass data now known for neodymiumup to neutron number N = 100 and for samarium up to N = 102. Therefore efforts such as those presented inthis work aiming to understand astrophysical outflowsvia the link between mass data and rare-earth abun-dances are timely. Vassh et al.
150 155 160 165 170 175
Mass Number (A) A b un d a n c e s/ k B =30, =70 msDZ33FRDM1995FRDM2012 HFB21Solar
150 155 160 165 170 175
Mass Number (A) A b un d a n c e s/ k B =10, =3 msDZ33FRDM1995FRDM2012 HFB21Solar Figure 1.
The predicted r -process, rare-earth peak abun-dances given two distinct astrophysical outflows ( s/k B = 30, τ = 70 ms, and Y e = 0 . s/k B = 10, τ = 3 ms,and Y e = 0 . r -process calculations. The solar data and uncer-tainties considered in this work are also shown for comparison(see Appendix A). 3. METHODHere we outline our method to explore masses ca-pable of rare-earth peak formation and demonstratethe diagnostic metrics used by showing results obtainedwith the astrophysical outflow condition considered inOrford et al. (2018). Our algorithm employs MarkovChain Monte Carlo (MCMC) to perform mass correc-tions to the Duflo-Zuker (DZ) mass model. We choosethis model as our baseline due to its lack of predictednuclear deformation in the lanthanide region which re-sults in calculated rare-earth abundances which are onaverage flat, as described in Mumpower et al. (2016b,2017) as well as the previous section. We make predic-tions for mass corrections to the DZ mass model withthe following mass parameterization (Mumpower et al. 2016b, 2017) M ( Z, N ) = M DZ ( Z, N ) + a N e − ( Z − C ) / f (1)applied exclusively to the nuclear masses in the rare-earth region. We randomly vary 28 a N parameters forneutron numbers N = 93 −
120 by generating Gaussiandistributed random numbers with a relative scaling σ ∼ . −
50% which is ideal to explore our large,multidimensional parameter space (Robert 1998). Wefix C based on preliminary runs in which we float thisquantity (in this work we use either C = 60 or C =58 depending on the outflow conditions) and set f =10 to ensure only local features in the mass surface areproduced in order to avoid modifying mass trends nearstability.After the adjustments to nuclear masses, we then cal-culate the astrophysical rates for the nuclei near therare-earth region with 45 ≤ Z ≤
69 at 93 ≤ N ≤ Q -values, β -decay rates,and neutron-capture rates for the roughly 300 nucleiin this region are performed in a self-consistent man-ner as in Mumpower et al. (2015, 2016b, 2017). Wefirst calculate β -decay Q -values and β -delayed neutronemission probabilities, P n , using the code BeoH (ver-sion 3.3.3) (Mumpower et al. 2016a). The majority ofthe nuclear data calculations are spent updating P n val-ues. We note that we use experimental data for decaysfrom NUBASE2016 (Audi et al. 2017) where availablein place of the β -decay predictions derived from our MCmasses, however for separation energies we use only thevalues from our MC procedure in our nucleosynthesiscalculations.Following the updates to separation energies and β -decay rates, we update photodissociation and neu-tron capture rates at each time step before calculat-ing the corresponding abundance prediction. We utilizethe neutron capture rates predicted by the CoH code(Kawano et al. 2016) and perform the fits introduced inMumpower et al. (2016b, 2017) which apply the func-tional form λ n,γ ( Z, N ) = exp[ a ( N, T )+ b ( N, T ) S n + c ( N, T ) S n ] (2)where a , b and c parameters are dependent on tem-perature, T , and neutron number, and S n is the oneneutron separation energy, S n ( Z, N ) = M ( Z, N − − M ( Z, N ) + M n , with M n being the mass of the neutron.This gives λ n,γ in units of cm mole − sec − . Photodisso-ciation rates are calculated from neutron capture rates CMC lanthanide properties and r -process dynamics λ γ,n ( Z, N ) ∝ T / exp (cid:20) − S n ( Z, N ) k B T (cid:21) λ n,γ ( Z, N − T is the temperature and k B is Boltzmann’s con-stant. After the relevant rates have been updated, wewrite input files for our nucleosynthesis code (PRISM)(as used in Mumpower et al. 2018; Vassh et al. 2019b;Zhu et al. 2018; Sprouse et al. 2020). Updating nuclearrates to reflect the MC mass values is critical to ensureour predicted abundance pattern reflects as realistic ofa set of nuclear data as is possible.We then evaluate how well the corresponding abun-dance pattern fits observed solar data. Before we cancompare our calculated abundances to the solar values,an overall scaling must be performed to either the solardata or the predicted abundances as is standard prac-tice. Such a scaling is needed since many factors suchas site mass ejection lie between solar data and nucle-osynthetic predictions. Since the values reported for thesolar data are themselves relative numbers often scaledaccording to the observed amount of silicon, it is therelative abundances which are meaningful when com-paring to nucleosynthesis calculations. Thus here wescale calculated abundances by determining the ratio ofsummed solar abundances between A = 150 −
180 andsummed calculated abundances for the same mass num-ber range. We make use of solar abundances, Y (cid:12) ( A ),and uncertainties, ∆ Y (cid:12) ( A ), derived from those given inGoriely (1999); Arnould et al. (2007) (see Appendix A).To consider the fit to the observational abun-dance data, we use the Metropolis-Hastings algo-rithm where the agreement between calculated abun-dances, Y ( A ), and solar data, evaluated as χ = (cid:80) A =150 ( Y ( A ) − Y (cid:12) ( A )) / ∆ Y (cid:12) ( A ), guides the evolu-tion of the Markov chain. We note that the numberof degrees of freedom used to determine the χ nor-malization depends on the number of correlations intro-duced by the parameterization (Berg 2004) which de-pends on how these parameters propagate to the finalguiding data, the abundances. In our MCMC applica-tion masses are propagated to rates and decays whichwill affect abundances in a non-linear, correlated fash-ion. Therefore since the number of correlations intro-duced by our parameterization cannot be determined ina physically meaningful way, we use an unnormalized χ when evaluating the likelihood function, L ∼ e − χ / .Since it is the likelihood ratio, R = L j L i , which deter-mines the acceptance or rejection of a new step, j , rela-tive to the previous step, i , it is important to recognizethat the common factor of the χ normalization will notaffect the MCMC evolution. The calculation begins from DZ masses such that ourMC parameters evolve away from zero to then explorethe parameter space for tens of thousands of steps. Wethen take the lowest χ solution, i.e. best step, found asthe solution from a single MCMC run. Although manysolutions with a χ which is significantly lower than thatof the initial DZ prediction are readily found, steps witha χ similar to the best step are more unique. This canbe seen in the first two panels of Figure 2 which showall steps taken during the MCMC evolution colored bytheir χ for two independent MCMC runs. Since thiscase starts with the DZ masses giving a χ > χ significantlylower than this but steps which have a χ comparable tothe best step are found in a similar region of parameterspace.Since the nuclear rates depend on several MC massvalues, and the rare-earth abundances used to calculatethe χ are determined by a convolution of nuclear re-actions, our parameters become highly correlated. Thiscauses our MCMC procedure to have a long integratedautocorrelation time and slow convergence. Since it isdifficult to ensure that an individual run explores thefull parameter space, we make use of the parallel chainsmethod of MCMC which determines the final solutionby taking the average and standard deviation of severalparallel, independent runs (Brooks et al. 2011). Thismethod also has the advantage of providing well-definederrors since uncertainty estimates are not directly de-pendent on the correlations between MC parameters.Our statistics are therefore determined from the averageand standard deviation of the lowest χ configurations(best steps) of 50 independent MCMC runs. Figure 2demonstrates that with each run taking a distinct paththrough parameter space on its way toward its lowest χ solution, the parallel chains method helps to ensurethat the full parameter space is explored.Given the use of the parallel chains method of MCMC,we determine convergence to a solution by consideringhow the average and standard deviation evolves as runsare collected. Figure 3 demonstrates how the averageand standard deviation determined from a set of 20 runscompares to the final full set of 50 runs. The averagecontinues to evolve as statistics are built when consid-ering 20 to 30 runs, but only adjusts slightly after theaddition of ten more runs when moving from a result de-termined from 30 runs to a 40 run result. However thereexists no significant difference between a 40 run resultand 50 run result implying that the addition of moreruns would provide no new information as it would onlyreinforce what has already been determined from the 50run search. Additionally the convergence of the solu- Vassh et al.
Figure 2.
Every step taken by two representative MCMCruns in the a N parameter space for N = 102 and N = 97(top and middle panels) colored by the χ of a given step.The complete space searched is demonstrated in the bottompanel which shows every fourth step for all of the 50 parallel,independent runs used to find the average (red star) andstandard deviation (red outlined box) which define the finalsolution. tion found by the parallel chains method of MCMC canbe evaluated by comparing the low χ regions found bythe parameter space search of an individual run to the
94 96 98 100 102 104 106 108 110 112Neutron Number (N) M M D Z [ M e V ]
20 runs30 runs40 runs50 runs
Figure 3.
The predicted MCMC masses for neodymium( Z = 60), relative to the DZ mass model, for the astrophysi-cal outflow condition considered in Orford et al. (2018). Theyellow lines denote the average and standard deviation deter-mined from 20 parallel, independent runs, red lines consider30 runs, blue lines consider 40 runs, and the black dottedlines show that the solution is converged upon in going from40 to 50 runs. low χ regions given the full set of all parallel runs (seeFigure 2). Such an analysis can be used as a diagnosticas to whether an individual run needs to be resumedin order to continue its search and potentially leave itspreviously identified local minimum.Since we aim to ensure that our mass predictions giveresults which are consistent with established nuclearproperties, we rein in the broad parameter space searchby imposing physical constraints. Firstly we consider acomparison with measured mass data by requiring that σ rms ( M, M
AME ) ≤ σ rms ( M DZ , M AME ) when con-sidering all nuclei with 140 ≤ A ≤ θ ( σ rms ( M, M
AME )) = , if σ rms ( M, M
AME ) > σ rms ( M DZ , M AME )1 , if σ rms ( M, M
AME ) ≤ σ rms ( M DZ , M AME ) . (4) CMC lanthanide properties and r -process dynamics σ rms ( M, M
AME )) requirement,we also consider the odd-even staggering behavior ofour MC mass values. This was implemented after initialruns located a number of solutions with an inversion inthe odd-even behavior of the one neutron separation en-ergies (see Appendix B). Since we considered this to bean unphysical mechanism of rare-earth peak formation,we introduced a check regarding the odd-even behaviorof our mass solutions using the neutron pairing metric D n ( Z, N ) = ( − N +1 ( S n ( Z, N + 1) − S n ( Z, N )) . (5)As can be seen from Figure 4, this metric reveals nu-clear structure via sharp transitions between nearby nu-clei where local nuclear properties suggest a region tobe especially stable, as is clear from predictions at theshell closures N = 82 and N = 126. The D n metriccan also hint at nuclear structure effects from pairingor collective effects such as deformation via its featuresbetween shell closures. Examining our baseline massmodel of DZ demonstrates a case where there is no en-hanced stability of rare-earth masses, as evidenced bythe purely odd-even behavior of the D n metric at andaround N = 103. After calculating D n , our algorithmvetos mass surfaces with an odd-even reversal in theirseparation energies via modifying the likelihood functionto include the step function θ ( D n ( Z, A )) = , if D n ( Z, A ) ≤ , if D n ( Z, A ) > D n metric does notexceed that of the N = 82 and N = 126 shell closures(i.e. the height of the largest peaks in Figure 4). That is,our modified likelihood function also includes the stepfunction θ ( D n ( Z, A )) = , if D n ( Z, A ) ≥ D n,AME ( Z, Z + 82)0 , if D n ( Z, A ) ≥ D n,DZ ( Z, Z + 126)1 , otherwise . (7)
80 90 100 110 120
Neutron Number (N) D n [ M e V ] DZ33FRDM1995 FRDM2012HFB21 AME2012CPT 2017
Figure 4.
The one neutron pairing metric, D n , for theneodymium chain ( Z = 60) predicted by the models consid-ered in Figure 1 as compared to AME2012 (Audi et al. 2012)and CPT at CARIBU (Orford et al. 2018) data. The impact of these D n metric checks are further dis-cussed in Appendix B. The complete modified likelihoodfunction which restricts the search to physically mean-ingful parameters is then L (cid:48) = L θ ( σ rms ( M, M
AME )) θ ( D n ( Z, A )) θ ( D n ( Z, A )) . (8)Note that since we use the σ rms check against AME2012data along with the D n metric checks to reject somecombinations of parameters outright before a step istaken, we effectively explore even more of the param-eter space than would be implied from examining thesteps taken in Figure 2. DISTINCT ASTROPHYSICAL OUTFLOWSThe nuclear physics feature which our mass adjust-ments can introduce, such as a sub-shell closure, hangsup material to form the peak. The location the al-gorithm finds such a feature to be needed depends onwhich r -process nuclei are dominantly populated whenthe neutron flux becomes exhausted (freeze-out) and de-cays to stability begin to take over. Therefore peakformation is determined by two aspects: (1) the initiallocation of the r -process path, i.e. the nuclei most pop-ulated along an isotopic chain, and (2) the dynamicswhich govern how the r process proceeds after freeze-out. We therefore considered outflow conditions withdistinct behavior: ‘hot’ scenarios in which the initialpath is the equilibrium path determined by (n, γ ) (cid:29) ( γ ,n)equilibrium (i.e. the Saha equation), and for which pho-todissociation continues to play a role after freeze-out,and ‘cold’ scenarios in which (n, γ ) (cid:29) ( γ ,n) equilibriumfails before the path populates the rare-earth region,therefore seeing nuclei closer to the dripline as its ini- Vassh et al. tial population, and showing little to no influence fromphotodissociation after freeze-out. We consider suchhot and cold scenarios for parameterized outflows whichare moderately neutron-rich and low entropy that willundergo heavy element nucleosynthesis. We emphasizethat although considering the heating introduced by nu-clear reactions can sometimes make cold dynamics differfrom their behavior when such reheating is neglected,this is not the case with all cold scenarios. In fact wefind that several scenarios can retain their cold behaviorafter including the reheating during the nucleosynthesiscalculation and thus cold dynamics remain a physicallyrealizable possibility in astrophysical environments. Theoutflow conditions considered in this work are all exam-ples which find nuclear reheating to have little to noinfluence on the expansion dynamics.Guided by merger simulations, we adopt three distincttypes of outflows (Surman et al. 2008; Metzger et al.2008; Perego et al. 2014; Fern´andez et al. 2015; Justet al. 2015; Radice et al. 2018) which could take placein both accretion disk and dynamical ejecta: (1) a hotoutflow with an entropy ( s ) of 30 k B /baryon and a dy-namical timescale ( τ ) of 70 ms, (2) a cold outflow with s = 10 k B /baryon and τ = 3 ms, and (3) a ‘hot/cold’outflow with s = 20 k B /baryon and τ = 10 ms. Herewe call this a ‘hot/cold’ outflow since it starts out char-acterized by hot r -process dynamics, and therefore theinitial r -process path is the equilibrium path, but be-haves similar to a cold outflow after freeze-out. All out-flows considered here are moderately neutron-rich withan electron fraction ( Y e ) of 0.20. These outflow param-eters are summarized in Table I. We note that in Orfordet al. (2018) we investigated whether our MCMC resultgiven outflow (1) was a viable solution in cases with sim-ilar outflow properties by considering slight adjustmentsto the entropy and expansion timescale. We found thatindeed similar expansion dynamics would require similarmass predictions in order to form a rare-earth peak com-parable to the solar data. Therefore, since similar out-flows require similar masses, the differences in requiredmasses given distinct outflow conditions such as thosein Table I can be used to discern the type of outflowscapable of accommodating both peak formation and thelatest mass measurements. Additionally we note thatalthough site ejection will actually be a mass weightedmixture of outflow conditions, the aim of exploring in-dividual trajectories here is to examine the dynamicswhich are dominant since similar dynamics will requiresimilar mass solutions.The temperature evolution, as well as a snapshot ofthe r -process path just as material begins to populatethe third peak region at N = 126, for these three cases is
80 90 100 110 120
Neutron Number (N) P r o t o n N u m b e r ( Z ) Time [s] T e m p e r a t u r e [ G K ] Hot OutflowCold OutflowHot/Cold Outflow
Figure 5. (Top) A snapshot of the r -process path just afterreaching the third peak for the ‘hot’ (red), ‘cold’ (blue), and‘hot/cold’ (green) outflows considered here. The grey regiondenotes the Duflo-Zuker dripline and black squares are stablenuclei. (Bottom) The temperature evolution as a function oftime for these astrophysical trajectories. Table 1.
Ejecta Outflow Parameters.Outflow Type Entropy (s/ k B ) Timescale (ms) Y e Hot 30 70 0.2Hot/Cold 20 10 0.2Cold 10 3 0.2 shown in Figure 5. In the hot outflow, (n, γ ) (cid:29) ( γ ,n) equi-librium persists throughout the r process correspondingto a path closer to stability than in the cold outflow.In contrast, the cold outflow sees photodissociation fallout of equilibrium early and relies almost entirely on β -decay to compete with neutron capture giving an ini-tial path closer to the dripline. In this case the mostpopulated nuclei in the rare-earth region are not wellrepresented by the equilibrium path even at early times.For the hot/cold outflow, although the initial r -process CMC lanthanide properties and r -process dynamics Figure 6.
The abundance weighted timescales for neutron capture (solid lines), photodissociation (dotted lines), and β -decay(dashed lines) as a function of temperature for all three outflow scenarios considered. path is the equilibrium path determined by (n, γ ) (cid:29) ( γ ,n)equilibrium, the nuclei initially populated lie closer tothe dripline since the hot/cold outflow has a lower ini-tial entropy than the hot outflow. This implies that thehot/cold case is more dense than the hot case at a giventemperature, such that the Saha equation sets an equi-librium path further out in neutron number (at nucleiwith lower separation energies).The distinct nature of the three outflows consideredhere can be best understood by examining Figure 6which shows the abundance weighted reaction timescales( ∼ γ ) (cid:29) ( γ ,n)equilibrium dominates early time dynamics with initialabundances well represented by the equilibrium path,however this equilibrium fails in the cold case very early,before the production of A ∼
195 nuclei. These threeoutflows also find themselves at very different tempera-tures when freeze-out begins, for instance the hot/coldoutflow has a comparatively low temperature at freeze-out ( ∼ ∼ β -decay in shaping the final rare-earth peakabundances. Therefore since both the initial populationof nuclei along the r -process path and freeze-out dy-namics vary among the three outflows described here,we expect these cases to require distinct MCMC masssolutions in order to form the rare-earth peak. RESULTSFor each of the three outflow conditions discussedin the last section, we obtain 50 independent, parallelMCMC runs and perform extensive testing and analysison the MCMC solutions found. In each case, all 50 runsare compared to evaluate how low of a χ is attainablefor the specific outflow being examined. Runs which didnot attain a χ around this value are resumed for moretimesteps to further explore the parameter space. Anadditional handful of runs, typically ones which foundthe lowest χ solutions of all runs but also had goodmovement through parameter space, are also resumedfor roughly twice as long as the standard run in orderto gain further confidence that the solution found wasin fact the global rather than local solution. The resultsfrom this procedure are presented in this section alongwith detailed descriptions of the physical mechanism bywhich the solution for each of the three distinct outflowsoperates. A summary and comparison of these solutionsis presented in the next section.5.1. Hot Outflows
Mass prediction and abundance pattern results for thehot outflow are shown in Figure 7. Note that Figure 7shows an updated result relative to that in Orford et al.(2018) since here further D n metric checks were imple-mented. However this did not significantly affect theMCMC solution for this case (see Appendix B). As de-scribed in Orford et al. (2018), this solution utilizes apile-up at N = 104 in order to form the rare-earth peak.In these hot outflow conditions, pile-up occurs becausethe updated separation energies produce a kink in theequilibrium path at N = 104, due to the dip in the red0 Vassh et al.
95 100 105 110 115 120
Neutron Number (N) M M D Z [ M e V ] AME2012AME2016CPT 2017FRIB Day 1FRIB Year 2FRIB Full StrengthANL N=126 Factory
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot Outflow Duflo-ZukerHot Outflow MCMCSolar
Figure 7. (Top) The predicted MCMC masses forneodymium ( Z = 60), relative to the Duflo-Zuker massmodel, for the hot outflow (red band). The band is producedfrom the average and standard deviation of results from 50parallel, independent runs. For comparison the AME2012data used to guide the calculation is shown, along with morerecent data from AME2016 (Wang et al. 2017) and the CPTat CARIBU (Orford et al. 2018) of which the calculationwas not informed. Potential future experimental reachesare shown as vertical lines. (Bottom) The standard devia-tion of the abundance patterns given by our 50 MCMC runs(red band) as compared to the baseline prediction using DZ(dashed black line). band mass surface of Figure 7. This N = 104 dip in themass surface may be accessible by next generation ex-periments such as the N=126 Factory and FRIB at fullbeam strength. The mass surface rise at N = 102 beforethe dip is also crucial to the solution in order place thepeak center at A = 164, as can be seen in Figure 8. Theupturn at N = 102 makes these nuclei less stable than N = 104 nuclei promoting further pile-up here at earlyand late times. With only the dip from N = 103 −
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot Outflow 50 Run AvgHot Outflow Duflo-Zuker a N N =101,102,103 a N N =103,104,105 a N N =107,108,109All 3 featuresSolar Figure 8.
The abundance pattern using the average massvalues from Figure 7 as compared to when only the N =102 or N = 104 or N = 108 key features are applied inan r -process calculation. The result when only these threefeatures are combined is also shown.
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot Outflow 50 Run AvgHot Outflow Duflo-Zuker only(n, ) only S n in ( ,n) onlySolar Figure 9.
The abundance pattern from implementing theupdated neutron capture, photodissociation, and β -decayrates determined from the average mass values in Figure 7 ascompared to when adjustments to only β -decay or neutroncapture or the separation energies in the detailed balanceequation for photodissociation are applied. aration energies are most important and the modifica-tions made by our algorithm to local neutron captureand β -decay rates play a minor role. This can be seenfrom Figure 9 which compares rare-earth abundanceswhen individual pieces of nuclear data are updated toreflect the mass adjustments shown in Figure 7. Re-call that photodissociation rates are dependent on bothseparation energies and neutron capture rates. For the‘ S n in ( γ ,n) only’ case in Figure 9 we explore the role of CMC lanthanide properties and r -process dynamics Figure 10.
The time evolution (top to bottom) for the formation of the rare-earth peak when the MCMC mass adjustmentsdetermined for the hot outflow are implemented in the neutron capture, photodissociation, and β -decay rates. Left panels showthe summed abundance as a function of mass number (red line) as compared to solar data. Right panels show the r -processpath (red) of the most abundant nuclei as compared to the equilibrium path (black boxes) determined by the Saha equation.Additionally populated nuclei are color-coded by the reaction or decay channel that a given species is most influenced by at thetime considered: pink when neutron capture dominates, yellow when photodissociation dominates, and blue when β -decay isdominant. The purple outline shows the range of experimentally established β -decay half-lives. solely the separation energy dependence in the exponen-tial of the photodissociation rate. For the ‘(n, γ ) only’case in Figure 9, our network instead updates photodis-sociation to reflect solely the new neutron capture rateand leaves separation energies to be those of DZ. Thisexercise demonstrates that only adjustments to equilib-rium path dynamics, via updating the separation en-ergies considered in the detailed balance equation, areneeded to form the peak in such hot outflow conditions.However, the changes introduced in the β decay and neutron capture rates from our mass adjustments playa role in shaping the sides of the peak.In Figure 10 we bring together the discussion of keymass surface features and their influence on each reac-tion and decay channel by showing the evolution of thepile-up which ultimately forms the rare-earth peak inthese hot outflow conditions. The figure indicates thedominant reaction or decay that a given species is un-dergoing at the time considered. Neutron capture dom-inates when the flow, F (rate times abundance), obeys2 Vassh et al. δ n,γ = F ( n,γ ) − F ( γ,n ) > | δ n,γ | > F β , photodis-sociation dominates when δ n,γ < | δ n,γ | > F β ,and when β -decay is dominant | δ n,γ | < F β . For the out-flow case considered here, the r -process path (location ofthe most abundant nuclei) is well described by the equi-librium path. At early times, there is no preferentialpile-up in the rare-earth region and the path is mostlyset by odd-even effects. When the onset of freeze-outoccurs, the path, still governed by (n, γ ) (cid:29) ( γ ,n) equi-librium, begins to move back to stability. After thistime the N = 104 kink in the equilibrium path whichcauses material to accumulate here emerges and per-sists throughout the remainder of the process. Thus thematerial which will populate the center of the peak isfound at N = 104 even at early times. At late times,after the environment falls out of (n, γ ) (cid:29) ( γ ,n) equilib-rium, the faster β -decay rates of nuclei at the bottomof the N = 104 kink in the path, as compared to theslower decay rates at the top of this feature, work asa funneling mechanism which keeps nuclei piled-up inthis location. That is, the nuclei near the bottom ofthe feature will quickly β -decay to then neutron captureback to N = 104, keeping the dominant population ofnuclei at lower Z aligned with the higher Z nuclei whichare decaying more slowly. Therefore, although Figure 9demonstrates that it was only modifications to separa-tion energies which our algorithm needed to exploit inorder to form the peak in this hot outflow, and modifica-tions to β -decay were not influential on peak formation,the late-time funneling produced by standard β -decaytrends is an important aspect of maintaining pile-up.Additionally at the latest times, for which the calcu-lation relies almost entirely on experimentally knowndecay rates, β -decay works to smooth the final pattern.5.2. Cold Outflows
In the case of cold outflow conditions, at early timesthe r -process path is in the most neutron-rich regionsnear the dripline. Therefore, unlike the hot case whichwas centered at C = 60, here calculations are centeredat C = 58 since initial tests with both C = 58 and C = 60 showed C = 58 to be able to find the lowest χ solutions given these astrophysical outflow conditions.We keep f = 10 in order to prevent significant changesin the dripline and keep the mass surface effects morelocalized. Figure 11 shows the results of 50 such parallel,independent MCMC runs. This solution also makes usein part of a pile-up at N = 104 in order to form the rare-earth peak. However this is achieved via a dip in themass surface at odd-N nuclei with N = 103, a slightlydifferent location than the dip at N = 104 observed forthe hot case. However, as is shown in Figure 12, this
95 100 105 110 115 120
Neutron Number (N) M M D Z [ M e V ] AME2012AME2016CPT 2017
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Cold Outflow Duflo-ZukerCold Outflow MCMCSolar
Figure 11.
Same as Figure 7 but with results for theneodymium mass surface (top) and final abundances (bot-tom) given the cold outflow. N = 103 feature is insufficient to form the peak alone.With only the dip from N = 102 − r -process pathlying closer to the dripline requires a feature at higherneutron number in order to properly redirect material tothe peak region, in this case at N = 108. With only thefeature at N = 107 − N = 103 and N = 108 features together which fill in the right andleft sides of the peak respectively.Here the peak formation process occurs so far out-side the region of experimentally established β -decayrates that the adjustments made by the algorithm tothe theory rates have a stronger influence than in thehot outflow, as can be seen in Figure 13. When only ad-justments to the separation energies considered in thedetailed balance equation are considered, more materialaccumulates in the rare-earth region than in the base-line case, but no real peak structure is seen. Howeverwhen only the neutron capture adjustments are applied,a clear peak structure is seen even though the height CMC lanthanide properties and r -process dynamics
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Cold Outflow 50 Run AvgCold Outflow Duflo-Zuker a N N =102,103,104 a N N =107,108,109 a N N =115,116,117All 3 featuresSolar Figure 12.
The abundance pattern using the average massvalues from Figure 11 as compared to when only the N =103 or N = 108 or N = 116 key features are applied inan r -process calculation. The result when only these threefeatures are combined is also shown.
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Cold Outflow 50 Run AvgCold Outflow Duflo-Zuker only(n, ) only S n in ( ,n) onlySolar Figure 13.
The abundance pattern from implementing theupdated neutron capture, photodissociation, and β -decayrates determined from the average mass values in Figure 11as compared to when adjustments to only β -decay or neutroncapture or the separation energies in the detailed balanceequation for photodissociation are applied. is lacking due to a need for early time accumulation ofmaterial. Therefore in such cold outflows the dynamicsare a complex interplay between neutron capture, pho-todissociation, and β -decay.The evolution of the formation of the peak in this coldoutflow case is shown in Figure 14 . Since such outflowsfall out of (n, γ ) (cid:29) ( γ ,n) equilibrium early, the equilib-rium path does not define the r -process path even atearly times. Early local competition between neutron capture and photodissociation initiates the pile-ups at N = 104 and N = 108 which we find to be responsiblefor peak formation. The influence of β -decay can alsobe seen early with N = 108 material being transferredto N = 106. Throughout the rest of the calculation, it isthe competition between β -decay and neutron capturewhich will determine the structure of local pile-up fea-tures. The material found at N = 106, ultimately orig-inating from N = 108, is later transferred to N = 104where a strong pile-up persists into late times. The peakis found off-center to the left throughout the majorityof the calculation and is eventually moved into place bylate-time neutron capture.5.3. Hot/Cold Outflows
Lastly we explore outflow conditions which fall out ofequilibrium in a manner not well represented by solelyhot or cold criterion, a hot/cold outflow, as was de-scribed in Section 4. In such outflows, the path at earlytimes is in more neutron-rich regions than the hot casebut does not push all the way to the neutron dripline asdoes the cold case. Therefore since preliminary calcu-lations centered at C = 58 and C = 60 were both ableto find low χ abundance pattern solutions, we gatherstatistics using runs center at C = 60 while again keep-ing f = 10.Figure 15 shows the results of 50 such parallel, inde-pendent MCMC runs. This solution also makes use ofa late time pile-up at N = 104 in order to form therare-earth peak where, as was the case in cold outflowconditions, this is achieved via a dip in the mass surfaceat odd-N nuclei with N = 103. However, it should benoted that in this case several of the 50 runs were able toachieve low χ solutions with the dip in the mass surfacefound at N = 104 rather than N = 103. As is shownin Figure 16, this N = 103 feature is not the dominantsource of pile-up and mostly works to fill in the rightedge of the peak. Rather, as in the cold case, with the r -process path lying in more neutron rich regions thanthe hot case, a feature at a neutron number higher than N = 104, in this case at N = 106, is required. With onlythe feature at N = 105 − β -decay rates. However as can be seen fromFigure 17, since here β -decay is most influential atlate times, the peak cannot be produced via changesto β -decay alone. Similarly neutron capture adjust-ments alone are insufficient to form the peak and infact produce abundances very near the baseline result.4 Vassh et al.
Figure 14.
Same as for Figure 10 but showing the time evolution given our MCMC mass and rate adjustments determined forthe cold outflow condition.
CMC lanthanide properties and r -process dynamics
95 100 105 110 115 120
Neutron Number (N) M M D Z [ M e V ] AME2012AME2016CPT 2017
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot/Cold Outflow Duflo-ZukerHot/Cold Outflow MCMCSolar
Figure 15.
Same as Figure 7 but with results for theneodymium mass surface (top) and final abundances (bot-tom) given the hot/cold outflow.
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot/Cold Outflow 50 Run AvgHot/Cold Outflow Duflo-Zuker a N N =102,103,104 a N N =105,106,107 a N N =115,116,117All 3 featuresSolar Figure 16.
The abundance pattern using the average massvalues from Figure 15 as compared to when only the N =103 or N = 106 or N = 116 key features are applied inan r -process calculation. The result when only these threefeatures are combined is also shown.
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot/Cold Outflow 50 Run AvgHot/Cold Outflow Duflo-Zuker only(n, ) only S n in ( ,n) onlySolar Figure 17.
The abundance pattern from implementing theupdated neutron capture, photodissociation, and β -decayrates determined from the average mass values in Figure 15as compared to when adjustments to only β -decay or neutroncapture or the separation energies in the detailed balanceequation for photodissociation are applied. When only adjustments to the separation energies con-sidered in detailed balance are applied, a peak structureis clearly produced and centered in the proper location,but is however insufficient to fill in the right side of thepeak. Therefore in such outflow conditions, it is primar-ily adjustments to (n, γ ) (cid:29) ( γ ,n) equilibrium via changesto the separation energies followed by the late-time shap-ing of the peak from β -decay which are responsible forpeak structure.The evolution of the formation of the peak for thishot/cold outflow is shown in Figure 18. At early times(n, γ ) (cid:29) ( γ ,n) equilibrium is obeyed, and the equilibriumpath is correlated with the r -process path. It is throughchanges in the separation energies of the detailed bal-ance equation that early time adjustments to the equi-librium path produce the initial pile-up of material at N = 106. However here, unlike the hot case, equilibriumis not obeyed into late-times and local competition be-tween neutron capture, photodissociation, and β -decaytransfers the early pile-up of material from N = 106 to N = 104. Later β -decay works to smooth the sharplyformed peak while late-time neutron capture shifts peakmaterial from A = 162 to A = 163. COMPARING THE DISTINCT SOLUTIONSFOUND FOR DIFFERENT ASTROPHYSICALOUTFLOWSIn the last section, we presented the mass surface so-lutions found by our MCMC method for each of thethree outflow conditions explored and outlined the peak6
Vassh et al.
Figure 18.
Same as for Figure 10 but showing the time evolution given our MCMC mass and rate adjustments determined forthe hot/cold outflow condition. formation mechanism in each case. Here we commenton the statistical diagnostics of each case and comparethese solutions directly. Table 2 provides a summary ofparameterization values and diagnostic criterion, as wellas the key features introduced by our mass adjustmentswhich we find most influential to the formation of themain peak region (from A ∼ − χ for the abundances predicted by the Duflo-Zuker mass model differ among the three cases, the χ which is able to be achieved by our MCMC procedure issimilar, with the hot/cold outflow runs able to achievethe lowest χ fits to the rare-earth peak. The numberof steps taken by our Monte Carlo runs and the accep-tance rate of runs is similar in these three cases since we work to ensure the parameter space is being properlyexplored.Given the common occurrence of pile-up at N = 104,it may appear that the three mass solutions are allequally capable of peak formation in any type of astro-physical outflow. However it is important to recognizethat each solution is tied to the distinct dynamics of thetrajectory explored. This is represented in Figure 19where we consider the formation of the rare-earth peakwith the solutions presented in Figures 7, 11, and 15in all three outflow conditions considered. The solutionfound to well form the peak in hot outflows producesa peak which is shifted greatly in mass number in coldand hot/cold outflows. The N = 108 feature crucial to CMC lanthanide properties and r -process dynamics Table 2.
MCMC Results Given Three Distinct Outflows.Outflow Type
C f
Baseline χ Avg. χ Avg. N = 102 , N = 103 , N = 103 ,
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot Outflow 50 Run AvgHot Outflow Soln Cold OutflowHot Outflow Soln Hot/Cold OutflowSolar
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Cold Outflow 50 Run AvgCold Outflow Soln Hot OutflowCold Outflow Soln Hot/Cold OutflowSolar
150 155 160 165 170 175
Mass Number (A) A b un d a n c e Hot/Cold Outflow 50 Run AvgHot/Cold Outflow Soln Hot OutflowHot/Cold Outflow Soln Cold OutflowSolar
Figure 19.
Rare-earth peak abundance predictions whenthe masses from Figure 7 are applied in nucleosynthesis cal-culations for all three outflows conditions (top), as comparedto when the masses from Figures 11 (middle) and 15 (bot-tom) are applied in all outflows considered here. peak formation in the cold case is sufficient to producea peak with outflows which are less cold, but producestoo strong of peak which is off center in both hot andhot/cold outflow cases. The N = 106 pile-up producedby the hot/cold case mass surface solution gives an offcenter peak at A ∼
166 in hot outflows while in cold out-flows a peak which is on average flat is produced withthis solution.In Section 5 we showed the MCMC mass predictionsfor all outflow cases for the neodymium isotopic chain.Here experiments have probed neutron-rich isotopes upto N = 100. Since key features for our predictions ineach case begin to be distinguishable just after this neu-tron number, neodymium comparisons would suggestall cases to presently be consistent with experimentalmasses (except for perhaps the cold case since it doesnot well reproduce the exact height of the N = 100rise seen in data). Luckily, mass values are available formore neutron-rich isotopes at higher Z . In Figure 20 wecompare our solutions directly with samarium ( Z = 62)experimental mass data which reaches up to N = 102.It is clear that this comparison suggests hot outflows tobe most consistent with the latest measurements. How-ever, recall that peak formation in the cold case is mostinfluenced by nuclei near Z = 58 so measurements atlower Z would be best to evaluate whether this remainsa viable peak formation mechanism. With the hot/coldcase solution centered at Z = 60, samarium masses bearmore weight on peak formation in such outflows. How-ever recall from Section 5.3 that N = 103 was a sup-porting feature in this case, playing a relatively minorrole in comparison with the need for pile-up at N = 106.Therefore measurements at higher N are needed.Additionally, in this work we intentionally exploredoutflows which do not host fission as fission fragmentscan greatly influence abundances in very neutron-richejecta. We note that very neutron-rich dynamical ejectacan also exhibit both ‘hot’ and ‘cold’ outflow conditions,so the path and freeze-out arguments laid out for themoderately neutron-rich cases in this work can be exam-ined in such cases as well. The complexities which canpotentially be introduced by fission deposition should beconsidered alongside local deformation of rare-earth ele-ments before more conclusive statements can be made as8 Vassh et al.
95 100 105 110 115 120
Neutron Number (N) M M D Z [ M e V ] AME2012CPT 2017
Figure 20.
The mass predictions for the samarium isotopic chain for all outflow conditions considered as compared to the mostneutron-rich mass data available from the CPT at CARIBU (Orford et al. 2018). to whether hot dynamics are indeed the most favorableoutflow conditions to form the rare-earth peak. CONCLUSIONSThe usefulness and broad applicability of statisticalmethods has lead to an increase in the popularity of suchapproaches in a wide variety of scientific applications.Here we made use of statistical techniques to model theunknown nuclear physics properties which could be re-sponsible for the pile-up forming the r -process rare-earthabundance peak. Such statistical methods which exam-ine the overlap between nuclear physics properties andastrophysical observations are able to fully exploit thecapabilities of next generation experiments and advanc-ing observational techniques.We presented how we have evolved our statistical ap-proach to soundly search parameter space and generatewell-defined statistical uncertainties by utilizing the par-allel chains method of MCMC. Additionally, by imple-menting a modified likelihood function, we are able totrain the algorithm to obey known nuclear physics. Herewe instruct the algorithm to take into account knownmasses of isotopes lighter than those being modeled aswell as some general nuclear physics properties via con-sidering the D n metric. We found that the algorithmsuccessfully meets AME2012 data when it is given thisinformation in a general way by considering a root-meansquare deviation summed over many isotopes. We notethat using a modified likelihood function approach to train the algorithm could be extended to consider otherproperties as neutron-rich rare-earth nuclei are furtherprobed by future experiment.We apply our method to the r -process rare-earthabundance peak since it is a window into the astrophys-ical outflow conditions which dominated the synthesisof lanthanide elements in our solar system. For thethree distinct types of outflows explored here, we findthat the features generating the pile-up which forms therare-earth peak correlate with the type of outflow con-ditions. Outflows which push out toward the dripline,that is outflows with cold r process dynamics which seephotodissociation drop out of competition early, showa need to pile-up abundances at a higher neutron num-ber, later transferring this bulk of material to lower neu-tron numbers. For the cold and hot/cold cases exploredhere such an early pile-up is found at N = 108 and N = 106, respectively. In the case of hot outflows withdynamics governed by (n, γ ) (cid:29) ( γ ,n) equilibrium, the r -process path finds itself comprised of isotopes which areless neutron-rich than is the case for outflows with colddynamics. The hot case considered here forms the rare-earth peak via a pile-up at N = 104 which occurs earlyand persists during the decay back to stability. This N = 104 feature is tantalizingly close to the mass mea-surements of neutron-rich rare-earth nuclei recently re-ported by the CPT at CARIBU, in particular being justtwo units in neutron number away from the recent mea-surement of Sm.
CMC lanthanide properties and r -process dynamics r process, fission frag-ments can play a role in shaping lanthanide abundances,introducing a potential connection between late-time fis-sion deposition and the rare-earth peak (Goriely 2015).However the predicted distributions of fission yieldsfor neutron-rich actinides vary widely (Goriely 2015;Goriely & Mart´ınez Pinedo 2015; Eichler et al. 2015;Mendoza-Temis et al. 2015; Roberts et al. 2011; Shiba-gaki et al. 2016; Vassh et al. 2019b,a) and the fissioningnuclei which most impact abundances are well outsidethe reach of experimental facilities (Vassh et al. 2019b).Thus whether fission deposition can influence rare-earthpeak abundances is highly uncertain. Therefore progressin understanding the synthesis of these elements is bestmade with further studies of neutron-rich lanthanideproperties, such as the masses, which are within reachof future facilities such as FRIB and the N=126 Factory.Interestingly, a past study suggested that since the ratioof odd-A and even-A r -process abundances agree best at A = 130 ,
163 and 195, such behavior supports conclu-sions that rare-earth abundances are determined by theinfluence of local nuclear properties during the r -processdecay to stability as opposed to fission deposition (Marti& Suess 1988). We will explore very neutron-rich out-flow conditions for which fission deposition occurs along-side an enhanced stability of rare-earth nuclei in futureMCMC investigations.For the moderately neutron-rich outflow conditionsconsidered in this work, we find that the most neutron-rich CPT mass data for samarium ( Z = 62) at N =100 −
102 agrees best with our MCMC results given thecase of a hot r process which does not push out to thedripline. The dip in the predicted mass data at N = 103which cold type dynamics use to keep material in placeat later times is not reflected by Sm mass measurements.If the trends suggested in the Sm isotopes also persistthrough lower Z where data cannot yet reach, the lat-est mass measurements favor rare-earth peak productionto occur in hot scenarios. Lower Z trends are of impor-tance since the r -process path can be found further fromstability in cold and hot/cold scenarios as compared tothe hot scenario. With peak formation in the cold out-flow most influenced by isotopes near the Z = 58 iso-topic chain, and features needed to form the peak in thehot case concentrated near Z = 60, higher Z elementssuch as Sm have less influence on rare-earth peak abun-dances. Other caveats should be noted. The N = 103 feature found in the cold and hot/cold cases is not theprimary feature of peak formation in these outflow con-ditions since pile-up must first occur at higher neutronnumbers. Therefore measurements are needed at higherneutron numbers where predicted features are concen-trated.Our results demonstrate that the formation of therare-earth peak is intimately tied to both the outflowconditions and nuclear properties of isotopes which arenear current experimental data. Therefore future mea-surements which push further into the neutron-rich lan-thanides can discriminate between different outflow sce-narios. Although here we use our approach to explorethe role of neutron-rich masses, this method can be re-sponsive to both new experimental information and the-ory developments for individual pieces of nuclear datasuch as neutron capture or β -decay rates. It is throughsuch collaborative efforts between theory and experi-ment that the nature of lanthanide production in as-trophysics can be understood.ACKNOWLEDGMENTSN.V. would like to thank Zolt´an Toroczkai, JorgePiekarewicz, Maria Piarulli, Daniel Phillips, Ian Ver-non and Michael Grosskopf for useful discussions. N.V.would also like to acknowledge the Information andStatistics in Nuclear Experiment and Theory (ISNET)conference series for valuable insights. The work of N.V.,G.C.M., M.R.M., and R.S. was partly supported by theFission In R-process Elements (FIRE) topical collabo-ration in nuclear theory, funded by the U.S. Depart-ment of Energy. Additional support was provided by theU.S. Department of Energy through contract numbersDE-FG02-02ER41216 (G.C.M), DE-FG02-95-ER40934(R.S.), and DE-SC0018232 (SciDAC TEAMS collabora-tion, R.S.). R.S. and G.C.M also acknowledge supportby the National Science Foundation Hub (N3AS) GrantNo. PHY-1630782. M.R.M. was supported by the USDepartment of Energy through the Los Alamos NationalLaboratory. Los Alamos National Laboratory is oper-ated by Triad National Security, LLC, for the NationalNuclear Security Administration of U.S. Department ofEnergy (Contract No. 89233218CNA000001). This workwas partially enabled by the National Science Founda-tion under Grant No. PHY-1430152 (JINA Center forthe Evolution of the Elements). This work utilized thecomputational resources of the Laboratory ComputingResource Center at Argonne National Laboratory (ANLLCRC) and the University of Notre Dame Center for Re-search Computing (ND CRC). We specifically acknowl-edge the assistance of Stanislav Sergienko (ANL LCRC)and Scott Hampton (ND CRC).0 Vassh et al.
APPENDIX A. SOLAR DATA WITH SYMMETRIZED ERRORBARSEvaluations which report uncertainties for the so-lar r -process residuals are rare, and the most widelyused, Goriely (1999), reports asymmetric uncertainties, Y + a − b . However, likelihood functions involving asymmet-ric errors are a rarely applied present area of studyin Bayesian methods (Barlow 2003, 2004; Marazzi &Yohai 2004). Here we choose to symmetrize the solaruncertainties in order to make use of a likelihood func-tion more consistent with established Bayesian methods.Symmetrization of asymmetric uncertainties is neededin many practical applications, as is the case with AMEand NUBASE mass and decay data evaluations (Wanget al. 2017; Audi et al. 2017). Thus, we adopt thesymmetrization approach outlined in Appendix A ofthe NUBASE2016 data evaluation (Audi et al. 2017).Rather than simply taking a midpoint and averaging theuncertainties a and b (as has been done with AME eval-uations in the past), this method maps an asymmetricnormal distribution into a symmetric normal distribu-tion by finding the median, m , which divides the asym-metric distribution into two equal areas: m = Y + √ a erf − (cid:0) a − b a (cid:1) , a > bY + √ b erf − (cid:0) a − b b (cid:1) , b > a. (A1)where we take erf − ( z ) ≈ √ πz/
2. The variance of thisequivalent symmetric distribution centered at m is givenby σ = (1 − π )( a − b ) + ab . Since large errors in thesolar abundance data can lead to some significant differ-ences between m and the original data point Y , we donot apply this procedure when the percent difference be-fore and after symmetrization is greater than 10%. Forthese cases we instead take the larger of the two errors torepresent the variance and leave the center unchanged.Figure 21 shows explicitly how the symmetrized solardata used in this work compares to that which was re-ported in the original evaluation. B. IMPACT OF D N METRIC CHECKWe use the neutron pairing metric, D n ( Z, N ) =( − N +1 ( S n ( Z, N + 1) − S n ( Z, N )), to provide the al-gorithm with feedback as to whether mass parametershave entered an unphysical regime. It is a useful di-agnostic since it is clearly connected to nuclear struc-ture, being influenced by odd-even effects and largest atclosed shells. Additionally, a negative value for the D n metric implies a reversal in the odd-even staggering of
150 155 160 165 170 175 180
Mass Number (A) A b un d a n c e Goriely 1999After Symmetrization
Figure 21.
Comparison between the r -process solar abun-dance residuals from Goriely (1999) and the data applied inthis work with symmetrized errors. the one-neutron separation energies and such odd-evenbehavior is not supported by any nuclear physics modelsor experimental measurements to date.In order to ensure algorithm parameters maintainphysical values, our modified likelihood function en-forces that at each time step we check that D n > D n checkthus prevents computational resources from being spentin unphysical regimes. Prior to implementing this check,roughly 40 −
50% of preliminary runs located solutionswhich violated the condition that D n remain positive.One such solution is shown in Figure 22. Solutions with D n < N = 110 in thefigure.In addition to requiring that D n remain positive, weimplement the criterion that at neutron numbers be-tween N = 82 and N = 126, the D n metric cannotbe larger than the values at these closed shells. Thatis, with D being the value of this metric given byAME2012 mass data and D being the value of thismetric predicted by the Duflo-Zuker mass model, at eachstep we check D n < D , D n < D , D n − D N − 94 96 98 100 102 104 106 108 110 112 Neutron Number (N) M M D Z [ M e V ] AME2012AME2016CPT 2017 94 96 98 100 102 104 106 108 110 112 Neutron Number (N) S n [ M e V ] Z = 60 AME2012AME2016CPT 2017 80 90 100 110 Neutron Number (N) D n [ M e V ] DZ33Ex of a final runEx without D n checkAME2012CPT 2017 Figure 22. The mass surface (top), one neutron separa-tion energy (middle), and D n metric (bottom) prediction foran MCMC run which implemented the D n > S n which we took to be unphysical. ric were implemented following several solutions locatedduring preliminary runs with the cold outflow. In suchcold outflows where the r -process path lies close to theneutron dripline, the algorithm exploited the ability toeffectively produce a new shell closure between N = 82 95 100 105 110 115 120 Neutron Number (N) M M D Z [ M e V ] AME2012AME2016CPT 2017 Figure 23. Comparison between the updated result withruns which did not violate the D n height check (maroon)and the published result from Orford et al. (2018) whichrequired only D n > and N = 126, which is not supported by experimentalmeasurements.By recursively examining our runs for the hot outflowcase previously published in Orford et al. (2018), wefound that 21 of 50 runs violated the D n height checkat some neutron number. However in this case, suchviolations were minor and were not key features formingthe rare-earth peak. To produce a full 50 run resultas is presented in the main text which obeys all D n metric conditions obeyed by the cold and hot/cold cases,we collected 21 new MCMC runs with such checks inplace. A comparison of our previously published resultto the new result presented in this work is shown inFigure 23 which confirms that the D n height check didnot significantly affect the MCMC solution we find forthis hot outflow condition. C. ADJUSTING THE PLOTTING SCALE OFEACH ISOTOPIC CHAINBecause relative differences in masses between iso-topes are most important for rare-earth peak formation,the parameterization we use is able to find the relevantmass surface trends. However, for an absolute massscale, our MCMC procedure can only predict the valuesfor the isotopic chain at which we center the calculation,typically neodymium ( Z = 60). Therefore when we dis-play our mass predictions for other isotopic chains, wepin the predicted trend to AME2012 data. To do sowe add a correction term to the parameterization whichwas applied by the Monte Carlo and use the modifiedversion of M ( Z, N ) = M DZ ( Z, N ) + a N e − ( Z − C ) / f + δ ( Z ) . (C2)2 Vassh et al. 90 95 1000.750.500.250.000.250.500.75 MM D Z [ M e V ] Z = 58 90 95 1000.750.500.250.000.250.500.75 Z = 59 90 95 1000.750.500.250.000.250.500.75 Z = 60 90 95 1000.750.500.250.000.250.500.75 Z = 61 90 95 1000.750.500.250.000.250.500.75 MM D Z [ M e V ] Z = 62 AME2012CPT 2017 90 95 1000.750.500.250.000.250.500.75 Z = 63 90 95 1000.750.500.250.000.250.500.75 Z = 64 90 95 1000.750.500.250.000.250.500.75 Z = 65 Neutron Number (N) Figure 24. The mass values for nearby isotopic chains determined by our mass parameterization formula using solely theAME2012 data for neodymium ( Z = 60) (red dashed line). We use the differences between the average of the overall trend fromparameterization values and experimental values when adjusting the scale of mass predictions (as shown by the black solid line). In order to determine δ ( Z ), we use AME2012 data toevaluate how well the parameterization fits with thisexperimental data. When centered at C = 60, weset δ (60) = 0 to find the a N values, a N,exp , at which M DZ (60 , N ) + a N,exp = M AME (60 , N ). We find a ,exp = 0 . 151 MeV, a ,exp = 0 . 128 MeV, a ,exp = − . 169 MeV, a ,exp = − . 440 MeV, a ,exp = − . a ,exp = − . 314 MeV, a ,exp = − . 354 MeV,and a ,exp = − . 162 MeV. The overall adjustmentto the mass of each isotopic chain, δ ( Z ), is thenbased on the average difference between AME2012and the parameterization prediction, that is δ ( Z ) = (cid:68) M AME ( Z, N ) − ( M DZ ( Z, N ) + a N,exp e − ( Z − C ) / f ) (cid:69) . Since we model N = 93 and above, we base our ad-justments on neodymium data for N ≥ 92 whichgives δ (58) = 0 . 498 MeV, δ (59) = 0 . 164 MeV, δ (60) = 0 . δ (61) = − . 195 MeV, δ (62) = 0 . δ (63) = − . 072 MeV, δ (64) = 0 . 163 MeV, and δ (65) = − . 042 MeV. This procedure is illustrated inFigure 24 which shows the mass parameterization pre-diction (both before and after considering δ ( Z ) adjust-ments) for multiple isotopic chains given the AME2012mass data for neodymium. We note that we havechecked that applying masses which include δ ( Z ) in our r -process calculations shows the MCMC solution foundto still produce a rare-earth peak.REFERENCES Abbott, B. P., Abbott, R., Abbott, T. D., & et al. 2017,ApJL, 848, L12—. 2018, PhRvL, 121, 161101Abbott et al., B. P. 2017, PhRvL, 119, 161101Arnould, M., Goriely, S., & Takahashi, K. 2007, PhR, 450,97Audi, G., Kondev, F., Wang, M., Huang, W., & Naimi, S.2017, Chinese Physics C, 41, 030001Audi, G., Wang, M., Wapstra, A.H., et al. 2012, ChinesePhysics C, 36, 1287Barat, E., Dautremer, T., & Montagu, T. 2007, in 2007IEEE Nuclear Science Symposium Conference Record,Vol. 1, 880 Barlow, R. 2003, in Statistical Problems in ParticlePhysics, Astrophysics, and Cosmology, ed. L. Lyons,R. Mount, & R. Reitmeyer, 250Barlow, R. 2004, arXiv e-prints, physics/0406120Barnes, J., Kasen, D., Wu, M.-R., & Mart´ınez-Pinedo, G.2016, ApJ, 829, 110Berg, B. A. 2004, Markov Chain Monte Carlo Simulationsand Their Statistical Analysis (World Scientific)Brooks, S., Gelman, A., Jones, G., & , Meng, X.-L., eds.2011, Handbook of Markov Chain Monte Carlo,Handbooks of Modern Statistical Methods (Chapman &Hall / CRC Press)Capano, C. D., Tews, I., Brown, S. M., et al. 2020, NatureAstronomy, arXiv:1908.10352 [astro-ph.HE]