The impact of (n,$\gamma$) reaction rate uncertainties of unstable isotopes near $N=50$ on the i process nucleosynthesis in He-shell flash white dwarfs

The first peak s-process elements Rb, Sr, Y and Zr in the post-AGB star Sakurai's object (V4334 Sagittarii) have been proposed to be the result of i-process nucleosynthesis in a post-AGB very-late thermal pulse event. We estimate the nuclear physics uncertainties in the i-process model predictions to determine whether the remaining discrepancies with observations are significant and point to potential issues with the underlying astrophysical model. We find that the dominant source in the nuclear physics uncertainties are predictions of neutron capture rates on unstable neutron rich nuclei, which can have uncertainties of more than a factor 20 in the band of the i-process. We use a Monte Carlo variation of 52 neutron capture rates and a 1D multi-zone post-processing model for the i-process in Sakurai's object to determine the cumulative effect of these uncertainties on the final elemental abundance predictions. We find that the nuclear physics uncertainties are large and comparable to observational errors. Within these uncertainties the model predictions are consistent with observations. A correlation analysis of the results of our MC simulations reveals that the strongest impact on the predicted abundances of Rb, Sr, Y and Zr is made by the uncertainties in the (n,$\gamma$) reaction rates of $^{85}$Br, $^{86}$Br, $^{87}$Kr, $^{88}$Kr, $^{89}$Kr, $^{89}$Rb, $^{89}$Sr, and $^{92}$Sr. (abridged)


Introduction
Most of the solar system chemical elements heavier than iron were produced in the slow (s) and rapid (r) neutron capture processes in previous generations of stars and stellar explosions (Käppeler et al., 2011;Thielemann et al., 2011). Cowan and Rose (1977) proposed that, under certain physical conditions, an intermediate (i) process with a neutron density N n ∼ 10 15 cm −3 intermediate between those typical for the s process (N n ≤ 10 11 cm −3 ) and r process (N n ≥ 10 20 cm −3 ) might be triggered in stars. Like for the main s process, neutrons for the i process are released in the reaction 13 C(α,n) 16 O. The higher neutron density in the i process is achieved at typical Heburning temperatures T ∼ 2 -3 × 10 8 K due to efficient replenishment of 13 C from the reaction 12 C(p,γ) 13 N followed by the decay 13 N(e + ν) 13 C. This combination of He and H burning reactions in one process occurs in a He-flash convective zone, where there is a plenty of He and 12 C, provided that a small amount of H is ingested into it when the upper convective boundary reaches the H-rich envelope layer. In this situation, the reactions 12 C(p,γ) 13 N and 13 C(α,n) 16 O are spatially separated, each taking place at its own favorable conditions, and 13 N with a half life of 9.96 min decays into 13 C while being carried down by convection with a comparable turnover timescale of ∼ 15 min.
It was not until the late 1990s that the first strong evidence of the i process in stars was detected. Asplund et al. (1999) investigated the evolution of the surface elemental abundances in the post-asymptotic giant branch (post-AGB) star Sakurai's object (V4334 Sagittarii) that at the time was experiencing a nova-like transient outbreak. A post-AGB star is a CO core of an AGB star that has lost most of its H-rich envelope as a result of an AGB super wind or a common-envelope event if it was in a close binary system. Such a star still has a He shell surrounded by a thin H envelope atop the CO core. When the star leaves the AGB it first moves almost horizontally to the left (to the higher effective temperatures) on the Hertzsprung-Russel diagram, as long as its luminosity is maintained by H burning in a shell, then, when the H-shell burning dies away, the star settles down on a white-dwarf cooling track (Paczyński, 1971). Depending on its thermodynamic properties and mass, the He shell can experience a very late thermal pulse (VLTP) when the star begins to cool down (Herwig, 2001b;Miller Bertolami et al., 2006). The VLTP is a He-shell flash similar to those occurring in thermally-pulsing AGB stars. During it, the temperature at the bottom of the He shell rises to ∼ 2 -3 × 10 8 K, and the entire He shell becomes convective. When the Heflash convection reaches the bottom of the H-rich envelope it begins to ingest H and this triggers the i process. VLTPs are experienced by a quarter of all post-AGB stars (Herwig, 2005). The observed outbreak of Sakurai's object was soon interpreted as a VLTP (Herwig, 2001a). The observations covered a time interval of about six months, during which the abundances of light s-process elements, such as Rb, Sr, Y and Zr, increased by approximately 2 dex. These abundance signatures can not normally arise during the preceding regular AGB evolution, and therefore provide evidence for active stellar nucleosynthesis during the He-shell flash. According to the stellar evolution models, a lag time of a few years (corresponding to the thermal time scale of the remaining thin envelope layer) is expected between the He-shell flash with its associated nucleosynthesis and the abundance observations on the stellar surface. Therefore the observed abundance changes over the six-month observational period probably signal mixing events rather than ongoing nucleosynthesis. Herwig et al. (2011) suggested that both the observed changes of the surface chemical composition of Sakurai's object and its light curve resulted from a He-shell flash that led to the physical conditions required for the i process, in particular, the vigorous convection that entrains H from a thin H-rich envelope to transport it into the He burning shell. To reproduce the observed abundances within the one-dimensional modeling approach adopted by Herwig et al. (2011), it was necessary to prolong i-process nucleosynthesis by adjusting the time when the He convective zone in Sakurai's object splits into two separate convective zones, the upper zone being driven by burning of the entrained H. It was assumed that the split occurred at a later time than predicted by the mixing-length theory (MLT) employed in the stellar evolution models. Recent three-dimensional simulations of H ingestion in Sakurai's object have revealed global horizontal oscillations of material near the top of the He convective zone caused by the violent burning of ingested H that resulted in a transient split of the convective zone , but this is not the split that quenches the i process, because not enough H has been burned at that time. Further investigations of the complex three-dimensional hydrodynamic nature of the convective-reactive He-shell flash with H ingestion in post-AGB stars are needed to provide 1D stellar evolution and postprocessing nucleosynthesis simulations with adequate properties of the convective Heshell split (Herwig et al., 2011;Stancliffe et al., 2011;Herwig et al., 2014).
Recently, a number of new abundance observations have revealed possible iprocess signatures in a broad variety of stars. This may indicate that i-process activation is not limited to post-AGB stars. Among these observations are the peculiar abundances of elements heavier than Fe in some of the carbon enhanced metal poor stars seemingly enhanced in both, r-and s-process elements (CEMP-r/s stars) (Beers and Christlieb, 2005;Masseron et al., 2010;Lugaro et al., 2012;Bisterzo et al., 2012). For some of these CEMP-r/s stars, e.g. CS 31062-050, the observed elemental abundances are consistent with predictions from i-process models, while the traditional explanation of a superposition of yields from the s and r processes shows discrepancies (Dardelet et al., 2014). The striking result of Dardelet et al. (2014) was that the observed abundances of several CEMP-r/s stars could be reproduced using simple one-zone i-process nucleosynthesis simulations, with constant temperature and density. A reasonable fit was obtained with a typical i-process neutron number density of N n ∼ 10 15 cm −3 and neutron exposure of τ ∼ 10 -50 mbarn −1 . These results have been recently confirmed by Hampel et al. (2016).
A possible i-process contribution has also been invoked in low-mass, low-metallicity AGB stars, to account for the Pb deficiency observed in post-AGB stars of low metallicity (Lugaro et al., 2015), and in young open clusters, to explain the largest [Ba/La] ‡ ratios that are neither compatible with the s process nor with a combination of s and r processes (Mishenina et al., 2015). Anomalous isotopic abundances possibly originating from an i process were also found in some presolar graphite grains (Jadhav et al., 2013), as well as in presolar SiC grains of type AB (Fujiya et al., 2013) and of type mainstream (Liu et al., 2014). Roederer et al. (2016) also proposed that the unusual [As/Ge] ratio (+0.99 ± 0.23) and the enhanced [Mo/Fe] and [Ru/Fe] ratios in the metal-poor star HD 94028 are the signatures of an i-process contribution with an even lower neutron exposure compared to the case of Sakurai's object.
With such a diverse set of observational indications it is conceivable that the i process is activated in several additional types of stars at various evolutionary stages. Stellar evolution models have pointed to various possible i-process nucleosynthesis sites including the He-core flash in low-mass stars (Campbell et al., 2010), the Heshell flash in a metal-poor AGB star (Fujimoto et al., 2000;Herwig, 2003;Iwamoto et al., 2004;Cristallo et al., 2009) or super-AGB stars . In addition, so-called rapidly accreting white dwarfs (RAWDs) have recently been proposed as possible site for an i-process (Denissenkov et al., 2017). RAWDs can occur in close binary systems where a white dwarf accretes H-rich material from a main-sequence, subgiant-or red-giant-branch star at rates of the order of 10 −7 M yr −1 . At these high accretion rates, H burns stationary resulting in a growing He shell underneath a thin H envelope. When the He shell mass reaches ∼ 0.01 − 0.03 M , a recurrent thermal He flash with H ingestion occurs. The nucleosynthesis is expected to be similar to Sakurai's object, and chemical evolution studies indicate that RAWDs could contribute significantly to the origin of Sr, Y, Zr, Nb, and Mo (Côté et al., 2017).
In order to understand the role of the i process in the origin of the elements at various stellar sites, observed abundances must be compared with astrophysical models for validation, to guide model improvements, and to constrain the underlying physical processes and conditions. This approach is limited by the observational errors and the uncertainties in the nuclear physics that link the stellar conditions to a particular abundance signature. Observational uncertainties are known and readily available as part of the published data. The goal of this paper is to also quantify the nuclear uncertainties, to enable a meaningful comparison of model predictions and observations that reveals information about the astrophysical i-process models. This is important, as current i-process models approximate complex mixing processes and make strong simplifying assumptions that need to be validated. For example, the abundance predictions from the post-AGB star model for Sakurai's object of Herwig et al. (2011) still exhibit discrepancies with observations, especially for Sr and Zr, but without quantifying the nuclear physics uncertainties it is impossible to conclude whether this indicates issues in the underlying astrophysical model or the assumed stellar parameters, such as the object's mass (Miller Bertolami and Althaus, 2007). Bertolli et al. (2013) investigated the impact of nuclear uncertainties on predicted ‡ We use the standard spectroscopic notation [A/B] = log 10 (N (A)/N (B)) − log 10 (N (A)/N (B)), where N (X) and N (X) are abundances of an element X in a star and the Sun.
Ba, La, and Eu abundances using a simple one-zone i process model and found significant uncertainties of the order of up to 1 dex. However, they did not identify the most important nuclear reactions, but instead carried out their analysis in terms of the underlying nuclear theory assumptions used in predicting reaction rates. In this paper we focus on the nuclear physics uncertainties in predicting the Rb, Sr, Y and Zr abundances in a full multi-zone post-AGB star model of relevance for Sakurai's object and RAWDs, and identify the key reaction uncertainties that need to be improved in the future. The paper is organized as follows. In section 2, we describe the methods and models used. In section 3, we discuss the nuclear physics uncertainties and the resulting reaction rate variation factors. In section 4 we present results, in particular an overall quantification of current nuclear physics uncertainties in the predicted elemental abundances, and an identification of the most significant sources of these uncertainties. The goal is to provide uncertainty information and to guide future nuclear physics work aiming at improving these uncertainties. A summary and conclusions are given in section 5.

Methods and models used in our nucleosynthesis simulations
To characterize the nuclear physics uncertainties in post-AGB star white dwarf He flashes we use a multi-zone post-processing model that reproduces the observed abundances in Sakurai's object reasonably well. The impact of nuclear reaction rate variations on the predicted abundances is determined in a Monte Carlo framework, where reaction rates are varied randomly within their uncertainties. The resulting abundance variations define the nuclear uncertainties of the abundance predictions, and correlations with the nuclear reaction rate variations allow us to identify the most critical nuclear reaction uncertainties. To guide future work, we also explore the appropriateness of simpler, less computationally demanding approaches such as the use of single-zone postprocessing and individual, one-by-one reaction rate variations. The one-zone model is also used to narrow down the number of reaction rates varied in the Monte Carlo analysis.
2.1. Multi-zone simulations 2.1.1. Stellar post-processing model For the multi-zone simulations of the i process in Sakurai's object, we use the NuGrid multi-zone post-processing nucleosynthesis parallel code mppnp  customized for the H-ingestion flash problem (Herwig et al., 2011). The simulations take temperature T , density ρ, radius r and diffusion coefficient D conv profiles from a stellar evolution model for the He convective zone in the mass range 0.5756 ≤ M r /M ≤ 0.5987 at the time of its first contact with the H-rich envelope.
In the present work, we use the same stellar evolution model ET14 of an asymptotic giant branch star of Herwig et al. (2006) during its last thermal pulse that was used by Herwig et al. (2011) in their RUN48 post-processing nucleosynthesis simulations (hereafter, RUN48 model). Its He-shell flash properties are similar to those of Sakurai's object. Herwig et al. (2011) interpolated T , ρ, r and D conv for a smaller (81) than in the original ET14 model (218) number of mass zones in the He-flash convective zone to reduce the RUN48 computation time (figure 1). To reproduce the high ratio of the first to the second peak s-process elements in Sakurai's object, they assumed that the i process in the He convective zone was choked off by its split at the 950th minute after the beginning of H ingestion. This empirically adjusted parameter actually sets up a time period during which the i process at the bottom of the He convective zone can change the abundances above the mass coordinate of the split. After the split is inserted, the convective mixing above it will homogenize the abundance distributions there. Therefore, the final abundances at the top of the He convective zone, that are compared with the surface chemical composition of Sakurai's object, can simply be computed by simulating the i process in this zone for 950 minutes and then massaveraging the resulting abundance profiles above the split.
In the present study, RUN48 model is implemented with the following modifications. First, we take M r , T and ρ at the bottom of the He convective zone from the stellar evolution model ET14 and obtain the M r , T and ρ profiles as functions of r in the He convective zone by numerically integrating the equations of hydrostatic equilibrium and mass conservation under the assumption of a constant entropy, as it should be in a deep convective zone, in the presence of an ideal gas and radiation. This is done with a view of our future stellar physics uncertainty and sensitivity study of the models of Sakurai's object, because this approach allows to easily vary the properties of the He-flash convection zone. Figure 2 shows that the integrated profiles match the original ones very well. Second, we use a higher resolution in the He shell with 132 mass zones instead of 81. These two modifications result in 13% reduction of the i-process simulation time after which our final abundances at the top of the He convective zone perfectly match their counterparts from the original RUN48 model. Third, we average the final abundance distributions over the entire He convective zone. Because the i-process nucleosynthesis takes place near the bottom of the He convective zone, the produced (destroyed) heavy elements have profiles that are slightly (because of convective mixing) increasing (decreasing) with the depth and time. Therefore, the mass-averaging of the abundance profiles over the entire He zone allows us to use a shorter (by 26%) i-process simulation time before we reach the same final surface abundances as in the original RUN48 model. Taking into account that we use a timestep of 74.3 sec versus 63 sec used by Herwig et al. (2011), the total number of timesteps required to best match the chemical composition of Sakurai's object is reduced from 1008 to 601, which significantly speeds up our Monte Carlo simulations. Finally, we apply a decay time of 2 years to the final surface abundances of unstable isotopes, which is much shorter than the 1 Gyr used in Herwig et al. (2011) and matches the age of Sakurai's object (Herwig, 2001b).
For the prepared (integrated and then interpolated for the specified mass zones) distributions of T , ρ, r and D conv the mppnp code calculates the changes in the He- RUN48 Figure 1. The MLT diffusion coefficient (green) and temperature (red) profiles in RUN48 model of the He-flash convective zone in Sakurai's object from Herwig et al. (2011). The D conv profile has a split that is artificially inserted at a specified time (the 950th minute after the beginning of H ingestion) to choke off the i process by preventing 13 C formed in the reactions 12 C(p,γ) 13 N(e + ν) 13 C near the top of the He convective zone, where H is ingested, from reaching its bottom, where neutrons are released in the reaction 13 C(α,n) 16 O. The model is divided into 81 zones (green circles), and H is ingested into the outer 5 zones.
shell chemical composition due to nuclear reactions and accounts for mixing between the He and H zones. As in Herwig et al. (2011), the ingestion of H into the He-shell convection zone is taken into account in a parameterized way, by adding ∆X = 10 −4 of H mass fraction to the upper 4 × 10 −4 M of the He convective zone at every 74.3 seconds, which corresponds to the H entrainment rate of 5.3 × 10 −10 M s −1 estimated from stellar evolution computations. The initial abundances of the nuclear species affected by nuclear burning in the pp-chains, CNO cycle, and He-burning were taken from the stellar model calculations of Herwig et al. (1999). The dominant species are X( 4 He) = 0.35, X( 12 C) = 0.43, and X( 16 O) = 0.19 (Herwig et al., 1999), which reflects the primary production of 12 C and 16 O. For heavier elements we use the solar system abundance distribution of Asplund et al. (2005) with isotopic ratios from Lodders (2003) scaled to what would be expected for a metallicity (mass fraction of elements heavier than He) of Z = 0.01 at the star's formation.
We confirmed that the same initial composition can be obtained by taking the complete set of solar abundances of Asplund et al. (2005) with isotopic ratios from Lodders (2003) NuGrid one-zone code ppn with T = 55 × 10 6 K and ρ = 10 2 g cm −3 , to complete H exhaustion, followed by He burning at T = 150 × 10 6 K and ρ = 10 3 g cm −3 , down to X( 4 He) = 0.35. This provides a cross-check of our initial abundances, and may be useful for future explorations of the impact of stellar parameter variations.
As in Herwig et al. (2011), a key parameter for obtaining an abundance pattern that is similar to observations is the duration of i-process nucleosynthesis. In the stellar model, this duration is set by a split of the He-shell convection zone that quenches the i-process. However, as discussed in Herwig et al. (2011), the time of this split is highly uncertain given the limitations of 1D stellar models in accurately describing convective processes. Our goal here is to determine nuclear physics uncertainties for testing the hypothesis that the i-process is responsible for the observed abundances in Sakurai's object. We therefore determine the simulation time for our model (744 min) by requiring the best match of the He-shell mass-averaged abundances of Rb, Sr, Y, and Zr with observed abundances, after a 2 year decay time for long-lived radioactive isotopes that reflects the estimated lifetime of Sakurai's object (Herwig, 2001b). The resulting final abundances are shown in figure 3.

Nuclear reaction network and nuclear physics input
The NuGrid nucleosynthesis network kernel that is used in both the multi-zone and the single-zone (section 2.2) simulations extends the network dynamically based on the strength of the nucleosynthesis fluxes, and may include up to ∼ 5,200 isotopes and ∼ 67,000 nuclear reactions (Herwig et al., 2008). Two parameters determine the dynamic network evolution. Species with half lives below a minimum β − -decay lifetime of 10 −3 seconds are not included. This sets the maximum width of the network away from the valley of stability. Another parameter sets the maximum mass of species to be included. Since in this study we are interested in the first-peak elements, we save computing time by not including species heavier than Ag. The maximum number of species to be included in the network is then 441. Limiting the network in this way has no impact on the first-peak elements (figure 4).
Most reaction rates are taken from the JINA REACLIB v1.1 library (Cyburt et al., 2010) with the following exceptions. The NACRE recommended rates (Angulo et al., 1999) are used for all of the CNO cycle reactions, except 13 N(p,γ) 14 O and 14 N(p,γ) 15 O for which the JINA REACLIB rates are used. The rates of the protoninduced reactions for isotopes with A = 20 -40 are taken from Iliadis et al. (2001). For the (n,γ) reactions with isotopes located along the valley of stability, we mainly use the revision 0.3 of KADoNIS data §  to linearly interpolate cross sections from the KADoNIS tables in energy, more modern or different (n,γ) nuclear rates being used for a sample of species . The compilations of Fuller et al. (1985), Oda et al. (1994), and Goriely (1999) provide the NuGrid codes with weak reaction rates.

One-zone simulations
We complement our multi-zone simulations with faster and simpler one-zone simulations and compare their results. For the one-zone simulations of the i process in Sakurai's object, we use the same model as in Dardelet et al. (2014) employing the NuGrid one-zone nucleosynthesis code ppn. The one-zone simulation setup assumes constant temperature and density, for which we adopt the typical He-shell condition values of T = 2 × 10 8 K and ρ = 10 4 g cm −3 .
Similar to the multi-zone case, the initial abundances for the one-zone simulations are taken from the solar system abundance distribution of Asplund et al. (2005) with isotopic ratios from Lodders (2003) (see text), respectively, and from the 84th minute of nuclear burning in the one-zone simulations that best match the multi-zone results. The small differences between the abundances from the two RUN48 models are entirely caused by the different decay times, 1 Gyr and 2 Yr, used in these models.
X( 1 H) = 0.1, which is supposed to mimic the H ingestion into the He convective zone, and X( 12 C) = 0.5 at the expense of 16 O, assuming that the new 16 O mass fraction is equal to the sum of the old mass fractions of 1 H, 12 C and 16 O minus the sum of the new mass fractions of 1 H and 12 C. The adopted H abundance is much higher than that expected in any zone of more realistic 1D or 3D multi-zone simulations. However, this initial condition in the one-zone model reproduces the neutron densities of the i process (Dardelet et al., 2014). All of our one-zone nucleosynthesis simulations of the i process in Sakurai's object were run with the same lists of isotopes and nuclear reactions and used the same reaction rates as their corresponding multi-zone simulations. We also used the same 2 year decay time after the end of the one-zone simulations.
The grey diamonds in figure 3 represent the abundances from the one-zone simulations at the 84th minute that best match the final abundances of the multizone simulations. Figure 5 shows the mass fractions and reaction fluxes in the region of the chart of nuclides that includes the stable (inside the heavy-line boxes) and unstable isotopes of Br, Kr, Rb, Sr, Y and Zr. log 10 X/X RUN48 multi-zone model implemented in this study same model with a small network up to Ag Sakurai's object in October of 1996 (Asplund et al. 1999) Figure 4. A comparison demonstrating that the multi-zone model used in this work reproduces reasonably well the global abundance distribution in Sakurai's object, in particular, the ratio of the abundances at the first and second s-process peaks. This comparison also shows that our using of a small network of 441 isotopes up to Ag does not affect the abundances of the elements lighter than Ag, while allowing to significantly reduce the computational time.
While the MC approach is suitable to identify the dominant nuclear uncertainties, there is the possibility that smaller, but still significant, sources of uncertainty are missed, especially when a large number of reactions contributes with comparable strength. We therefore complemented our MC study with a series of simulations D1, in which only one reaction rate is changed per run. For each reaction rate, two runs are carried out with the maximum multiplication factor f i = v max i ("rates up") and the minimum multiplication factor f i = (1/v max i ) ("rates down"). This approach determines to first order the partial derivatives ∂X k ∂r i in equation 1. In order to get an understanding of the significance of correlations between reaction rate changes, we also carried out a complete set of calculations D2 with variations of all possible reaction rate pairs. These calculations determine all ∂ 2 X k ∂r i ∂r j partial derivatives in equation 1.

Nuclear physics input uncertainties
For a given neutron density evolution (which is affected by its own nuclear uncertainties from the neutron producing reactions, and nuclear uncertainties in the prior stellar evolution) i-process nucleosynthesis calculations rely on neutron capture rates and β − decay rates within the i-process band, which extends from two to eight neutrons away from the valley of stability (figure 5). In this mass range the terrestrial β − -decay rates used in this work are all experimentally well known, though corrections for the stellar environment may be significant in some cases. However, none of the neutron capture rates on unstable nuclei are experimentally constrained. We therefore focus in this study on the role of uncertainties in theoretical neutron capture rate predictions. Neutron capture rates are predicted using the Hauser-Feshbach model of statistical decay of a compound nucleus. This theoretical model describes the decay of "highly" excited nuclei with large numbers of levels per MeV, provided that a reasonable description of nuclear statistical properties (for example level density and gamma ray strength function) and nuclear potentials are available. Two sources of uncertainty need to be distinguished (1) the uncertainty due to the inherent limitations of the statistical model that describes an expected average behavior of nuclei (intrinsic uncertainty) and (2) the uncertainty due to the inability of models for nuclear potentials, level densities, and gamma-ray strength functions to describe the average behavior of nuclei away from stability where experimental data on these quantities are not available (extrapolation uncertainty).
The intrinsic model error can be constrained by comparison of rate predictions with experimental reaction rate data in mass regions where experimental data constrain global descriptions of level density, gamma ray strength functions and nuclear potentials. Currently reliable measurements of astrophysical neutron capture rates are only possible for stable or very long-lived isotopes. A comprehensive comparison of neutron capture rate predictions from statistical models with experimental data was carried out by Beard et al. (2014) and demonstrated that the intrinsic uncertainty of the statistical model approach is of the order of a factor of 2.
The extrapolation uncertainty for predictions of neutron capture rates away from stability in the i-process region is more difficult to estimate. We follow here the approach of Liddick et al. (2016) who took advantage of the fact that for each of the statistical model ingredients (nuclear potentials, level densities, and gamma-ray strength functions) there is a range of theoretical models that describe the data near stability reasonably well. The divergence of the predictions of these various models away from stability can then be taken as a lower limit of the model uncertainty, which can then be propagated to the nuclear reaction rates. Liddick et al. (2016) found that for r-process temperatures of 1.5 GK within a few mass units away from stability rate uncertainties in the Mn-Ga range quickly exceed factors of 10 or even 100. We apply here the same approach to the nuclide region and lower temperatures (2 × 10 8 K) of relevance for the i-process in Sakurai's object.
We use the statistical model code TALYS ¶ (Bersillon et al., 2007) and take advantage of various options for potentials, nuclear level densities (NLD) and γ ray strength functions (γSF). The choice of input models was restricted to the models listed in table 1. The selection criteria are discussed in more detail in Liddick et al. (2016). Comparison calculations with the Koning-Delaroche (Koning and Delaroche, 2003) and Jekeune-Lejenne-Mahaux (Jeukenne et al., 1977) optical potentials yielded no significant differences. We therefore used the Koning-Delaroche potential for all calculations.
Calculations were carried out with all combinations of the 5 NLD models and 4 γSF models in table 1. The ratio of the largest and smallest neutron capture rate is indicated in figure 6 and increases from less than a factor of 5 to more than a factor of 20 for increasingly neutron rich nuclei. We bin the ratio into ranges and use the values indicated in figure 6 as an estimate of the extrapolation uncertainty in the reaction rate predictions. As it is generally much larger than the intrinsic uncertainty of about a factor of 2, we use this ratio as the maximum variation factor v max i in our sensitivity study, v max i = 30 being used for v max i > 20. The subset of neutron capture reactions varied in this study (see figure 6) was identified using the single-zone model. We include 52 neutron capture rates that exhibit strong reaction fluxes log 10 f > −14.4 at the peak of i-process processing just prior to its quenching by the split of the convection zone (figure 5). 90 Y was included for completeness.

Results
In this section we present the simulation results for varying the radiative n-capture rates of the 52 unstable isotopes shown in figure 6 both randomly and systematically. We analyze the impact of these reaction rate variations on the predicted surface abundances of Rb, Sr, Y and Zr in models of Sakurai's object. Table 1. List of models used to describe the nuclear level density and the γ strength in the Hauser-Feshbach calculations. Calculations were performed with all 20 possible combinations of these models.
For comparison, the grey-shaded area shows the observed Zr abundance represented by a normal distribution with the mean and standard deviation from Asplund et al. (1999).

Multi-zone Monte-Carlo simulations
Most of the abundance distributions from the 10,000 Monte Carlo variations of the 52 (n,γ) reaction rates follow closely a normal distribution. As an example, figure 7 shows the distribution of the predicted surface abundance of Zr (logarithmic and with respect to its solar value). We therefore fitted the distributions with a normal distribution to obtain standard deviation and mean. Figure 8 shows the Monte Carlo predicted abundance distributions, together with mean and standard deviations from a fit with a normal distribution, in comparison with observations from Sakurai's object. It is evident from this figure that nuclear uncertainties are significant and overall comparable to observational errors.  (Asplund et al., 1999).  explains the deviation of its resulting abundance distribution from a normal distribution (figure 8). Sr is affected by mainly 2 uncertainties -87 Kr(n,γ) and 88 Kr(n,γ). Y and Zr are affected by a larger number of reactions with comparable |r P |. Table 2 also shows that some reactions have a strong impact on elemental abundance ratios. For example, 87 Kr(n,γ) affects strongly Rb/Sr, 88 Kr(n,γ) affects strongly Sr/Y, and 89 Rb(n,γ) affects strongly Y/Zr. Some of these results can be anticipated if we compare the mass fractions of isobars from the one-zone simulation presented in upper panel of figure 5. From such a comparison, one may come to the following conclusions. The final abundance of Zr could have been most affected by the production of 90 Kr if its decay to 90 Zr had not been shielded by the long-lived 90 Sr whose half life time t 1/2 ≈ 28.8 yr considerably exceeds our used 2 yr decay time. Therefore, the Zr abundance in Sakurai's object is actually determined by the synthesis of 92 Sr. The results of our statistical analysis in table 2 confirm this conclusion. The presence of 85 Br,86 Br,87 Kr and 88 Kr with the positive correlation coefficients in the Zr, Y and Sr columns of table 2 rather indicates that the n-capture cross sections of these isotopes determine the total strength of the i process in our model. 89 Y is the only stable isotope of Y. It is mainly produced from the decay of the short-lived 89 Kr. The only unshielded stable isotope of Sr is 88 Sr. It is predominantly produced from 88 Kr that has the maximum mass fraction among the isobars with A = 88. Rb has two unshielded stable isotopes, 85 Rb and 87 Rb, but the production of 85 Rb via the β − -decay of the abundant 85 Br is hindered because it passes through the relatively long-lived (t 1/2 ≈ 10.8 yr) isotope 85 Kr. Therefore, the final Rb abundance is determined by the synthesis of 87 Rb from 87 Kr, which agrees with the results presented in table 2.
4.1.1. Multi-zone variations of individual rates (D1 and D2) We also carried out 104 multi-zone model calculations where each of the 52 reaction rates was varied up and down by the highest (f i = v max i ) and lowest (f i = 1/v max i ) variation factor, respectively (calculation set D1). Based on this simplified approach we find the sensitivities listed in tables 3 and 4. Almost all the important reaction rate uncertainties identified with the Monte Carlo study are also identified as significant sensitivities in the single rate variations when we select the reactions for which the condition log 10 |X f i =1 /X f i =1 | ≥ 0.075 is true for at least one of the final abundances of Rb, Sr, Y and Zr. There is only one additional reaction of comparable sensitivity in the single rate variations -85 Kr(n,γ) -that does not show a correlation with |r P | > 0.15 in the Monte Carlo study. This indicates that single rate variations are as sensitive in identifying sources of uncertainty as the correlation analysis of the Monte Carlo results.
We also varied all possible pair combinations of the 52 neutron capture reactions to explore correlations between uncertainties (calculation set D2). This required 5304 additional multi-zone calculations. The results of the D2 simulations are presented in the second column of table 4. As expected, the maximum changes in final abundance are larger when two reactions are changed simultaneously. However, the list of important Table 3. The list of reactions (the first column) for which the absolute magnitude of at least one of the logarithmic ratios of the abundances of Rb, Sr, Y and Zr predicted with the rate multiplication factor f i varied up and down relative to the case of f i = 1 (shown in the other columns) has exceeded 0.075 in our multi-zone D1 runs.
Reaction Rb (up/down) Sr (up/down) Y (up/down) Zr (up/down) 85 Br ( Table 4. Individual isotopes (D1) and their pairs (D2) whose maximum (n,γ) reaction rate variations have the largest impact on the predicted elemental abundances and their ratios (the absolute magnitudes are given in parenthesis) from multi-zone simulations of the i process in Sakurai's object.

D1 D2
nuclear physics uncertainties is the same, again indicating that correlations among rate changes are negligible, and that for the case of the i-process a single rate variation approach is appropriate for identifying the critical nuclear reactions.

One-zone simulations
The D1 and D2 series of simulations have been repeated using the one-zone model of the i process in Sakurai's object. The results are presented in table 5. The onezone simulations identify neutron capture rates on 85 Br, 87 Br,87 Kr, 88 Kr, and 89 Kr as the dominant uncertainties affecting the final abundances. Although this list partially overlaps with the one obtained from multi-zone simulations, it over-emphasizes the importance of 87 Br(n,γ) and overlooks 86 Br(n,γ), 89 Rb(n,γ), 89 Sr(n,γ), and 92 Sr(n,γ). We conclude that single zone models are not adequate for i-process sensitivity studies of convective-reactive processes, such as the i process in Sakurai's object.

Summary and Conclusions
We studied the impact of nuclear uncertainties on the nucleosynthesis predictions of a multi-zone i-process model for the post-AGB star Sakurai's object to determine whether this model can explain the large observed enhancements in Rb, Sr, Y, and Zr. We analyzed the impact of the uncertainties in the theoretical predictions of 52 (n,γ) reaction rates on neutron rich unstable isotopes of Br-Zr with neutron numbers around N = 46 − 61. These reactions cannot be measured directly in experiments as both reactants, the neutron rich isotope and the neutron, are unstable. We estimated the theoretical uncertainties by varying input parameters in the TALYS nuclear statistical model code and find, in agreement with previous work (Liddick et al., 2016) that uncertainties increase dramatically with distance from the valley of stability due to the lack of experimental constraints on key ingredients, such as level structure and γstrength distribution. These reaction rate uncertainties can reach more than a factor of 10 across the band of the i-process. We propagated these errors through the stellar model using a Monte Carlo approach and the NuGrid 1D multi-zone post-processing framework mppnp that enables sufficiently fast modeling with large nuclear reaction networks, while taking fully into account convective mixing processes. We find that the nuclear uncertainties in the predicted abundances can reach 0.3 dex and are thus comparable to current observational uncertainties. With both, observational and nuclear uncertainties quantified, we are now in a position to compare abundance predictions with observations and test the appropriateness of the underlying astrophysical model (see figure 3). Model predictions agree well with observations for Rb and Y. Abundance predictions are low for Sr, and high for Zr. However, they do agree with observations at the 2σ level. To determine whether these differences are significant and point to issues with the astrophysical model assumptions would require a significant reduction of uncertainties in either observations, nuclear physics, or both.
There are significant opportunities to reduce the nuclear physics uncertainties. Experimental techniques have been developed to indirectly constrain neutron capture rates on unstable nuclei. Neutron transfer reactions (Cizewski et al., 2007;Bardayan, 2016) and other surrogate techniques (Escher et al., 2012) populate compound nuclei and study their properties. The β-Oslo technique uses β − -decay to constrain level densities and γ-strength functions and has specifically been developed to reduce uncertainties in statistical model predictions of neutron capture rates on unstable nuclei (Liddick et al., 2016). Furthermore, significant increases in intensities of radioactive beams needed for these measurements are expected with new rare isotope beam facilities such as FRIB in the US or FAIR in Germany coming online in the near future.
To guide future experimental and theoretical efforts in reducing nuclear physics uncertainties in i-process models we identified the most critical reactions that are the dominant sources of uncertainties. The reactions identified using the Pearson productmoment correlation coefficient between rate and abundance changes in the Monte Carlo calculations, and reactions identified by single variation of individual reaction rates mostly agree. This demonstrates that single rate variations are an appropriate approach to identify critical reaction rates. The important reaction rate uncertainties for i-process production of Rb, Sr, Y, and Zr are neutron captures on 85 Br,86 Br,85 Kr, 87 Kr, 88 Kr, 89 Kr, 89 Rb, 89 Sr, and 92 Sr. Neutron captures on 85 Kr were only identified in the single rate variations, indicating that their uncertainties are less dominant than the others. Single rate variations using a simplified one-zone i-process led to significantly different results. This shows that one-zone approximations are not reliable for identifying critical reaction rates in convective-reactive regimes such as the post-AGB star i-process. The 1D stellar model used here has significant astrophysical uncertainties. For example, the split of the convection zone that quenches the i-process, and therefore actually determines its duration, had to be delayed compared to what has been predicted by the mixing length theory in the 1D stellar model in order to obtain a strong iprocess production of Rb, Sr, Y, and Zr (Herwig et al., 2011). The goal here was to investigate whether comparison with observations reveals any deficiencies of the underlying astrophysical model that would guide future model improvements to address astrophysical model uncertainties. We found that with current nuclear uncertainties model predictions and observations are consistent, though there are some possible hints of discrepancies. Future improvements in nuclear uncertainties are needed to provide strong guidance for model improvements. In addition, it will be important to investigate the sensitivity to astrophysical parameters such as H-ingestion rate and duration, r, T and ρ at the bottom of the convective zone, Z, D conv , or mass resolution. Should future work reveal strong significant discrepancies with observations, such sensitivity studies can then point to specific issues in the model.
Our results were obtained using an i-process model for the post-AGB star Sakurai's object. The model does reproduce the observed abundances within nuclear and observational errors. We therefore expect our nuclear sensitivity results to be applicable to any post-AGB star i-process model that reproduces the observed abundances of Rb, Sr, Y, and Zr. Our results should also be directly applicable to near-solar metallicity rapidly accreting white dwarfs, where models have indicated that an i-process occurs under similar physical conditions (Denissenkov et al., 2017). The computational method and analysis tools that we have developed and used in conjunction with the NuGrid postprocessing nucleosynthesis codes can be used for uncertainty and sensitivity studies of many other convective-reactive nucleosynthesis sites.