of free neutrons Nn, whose typical values for the s and r processes are 10 7 cm -3 < ∼ Nn < ∼ 10 11 cm -3 and Nn > ∼ 10 20 cm -3 , respectively (e.g., Thielemann et al. 2011;Käppeler et al. 2011). Cowan & Rose (1977) suggested that some n-capture nucleosynthesis might also take place in stars at intermediate values of Nn ∼ 10 13 -10 15 cm -3 . Such an i process can occur during a He-shell flash, when convection driven by He burning entrains H from its surrounding H-rich envelope, and neutrons are produced in the 12 C(p,γ) 13 N(e + ν) 13 C(α,n) 16 O reaction chain, e.g. in thermally-pulsing asymptotic giant branch (AGB) and post-AGB stars (Cowan & Rose 1977;Malaney 1986;Jorissen & Arnould 1989;Cristallo et al. 2009;Cristallo et al. 2016;Karinkuzhi et al. 2021).

The first strong observational evidence of an ongoing i process in a post-AGB star was the discovery of a shortterm (on a timescale of several months) enhancement of the surface abundances of the first peak n-capture elements Rb, Sr, Y and Zr, those with neutron numbers near magic N = 50, in Sakurai's object (V4334 Sagittarii) by Asplund et al. (1999). Herwig et al. (2011) explained that enhancement as a result of a very late thermal pulse of a He shell that had led to an H ingestion by He-shell convection and an i-process nucleosynthesis that was interrupted by a split of the He convective zone. That split was found to be a result of the hydrodynamic feedback from the violent energy release in the reaction 12 C(p,γ) 13 N operating on timescales of ≈ 15 min in the convective-reactive environment (Herwig et al. 2014b).

The first study of the impact of (n,γ) reaction rate uncertainties on the predicted abundances of the elements with the atomic number in the range 52 ≤ Z ≤ 63, in many aspects similar to our studies based on one-zone Monte Carlo simulations, results of which are presented in this paper, was done by Bertolli et al. (2013). Later, Dardelet et al. (2014) showed that the peculiar surface abundances of the elements from Sr to Ir found in the CEMP-r/s stars, the latter representing a sub-class of the carbon enhanced metal-poor stars seemingly polluted by products of both s and r process (Bisterzo et al. 2012;Abate et al. 2016), could in fact be surprisingly well fitted with i-process nucleosynthesis predictions. The results of Dardelet et al. (2014) were based on one-zone nucleosynthesis simulations in which the initial conditions, such as H abundance and temperature had been chosen to facilitate the large intermediate neutron density that is characteristic for the i process.

The i-process site creates by its nature a regime in which convective and nuclear burning time scales, specifically the time scale of the 12 C(p,γ) 13 N reaction, are equal somewhere inside the He-burning convection zone (e.g., see Cristallo et al. 2009 andappendix B of Herwig et al. 2011). Therefore, any physical model of the i-process conditions must incorporate both the mixing and the nuclear processes. Onezone models cannot accomplish this. Simple one-zone models aiming to explore the heavy-element nucleosynthesis at sufficiently high neutron density that adopt a constant temperature and create neutrons through H-burning and the 13 C(α,n) 16 O reactions then face a dilemma (Bertolli et al. 2013). At too low temperature the (α,n) reaction proceeds too slow to generate high Nn no matter how much 13 C is available, while at too high T the 13 N(p,γ) reaction is faster than the β decay and thereby limiting the formation of 13 C and preventing to reach high Nn. One simple trick is to just turn off the 13 N(p,γ) reaction and thereby obtain the high neutron densities that naturally emerge in the more realistic multi-zone simulations, as was done in the one-zone models with neutron creation via charged-particle reactions in Dardelet et al. (2014) and McKay et al. (2020) foot_0 . However, compared to multi-zone models (e.g., stellar-evolution based i-process models of Denissenkov et al. 2019) or the new 3Dbased advective-two-stream models of Stephens et al. (2020), one-zone models cannot provide any meaningful insight into the actual formation process of neutrons through convectivereactive H burning, and therefore the better solution is to just work with a constant neutron density approach as done in Hampel et al. (2016), Hampel et al. (2019) and also in McKay et al. (2020). Hampel et al. (2016) and Hampel et al. (2019) used a larger number of CEMP-r/s stars and a different fitting technique to confirm the conclusions of Dardelet et al. (2014). In all of those works the actual stellar site of the i process remained unknown, therefore they employed a one-zone nucleosynthesis model with a single set of specified density and temperature or a single value of Nn considered as free parameters.

The i-process conditions are found whenever H is entrained into a convective He-burning shell (Campbell et al. 2010;Herwig et al. 2011). Such H-ingestion flashes (HIFs) or proton-ingestion events (PIEs)foot_1 may occur at various stellar evolution phases in such objects as low-Z AGB stars (e.g. Fujimoto et al. 2000;Herwig 2003;Iwamoto et al. 2004;Campbell & Lattanzio 2008;Cristallo et al. 2009;Suda & Fujimoto 2010;Cristallo et al. 2016;Karinkuzhi et al. 2021), post-AGB stars (e.g. Iben Jr et al. 1983;Herwig et al. 1999a;Miller Bertolami et al. 2006;Herwig et al. 2011), super-AGB stars (Jones et al. 2016) or massive stars (Ritter et al. 2018;Banerjee et al. 2018;Clarkson et al. 2018;Clarkson & Herwig 2020). However, to reach the second neutron-magic peak around Ba and also obtain the large ratio hs/ls of the second-(heavy s-process) to the first-peak (light s-process) n-capture elements as observed in CEMPr/s stars requires to sustain the intermediate neutron density for long enough to reach high neutron exposures. Such conditions have not yet been self-consistently obtained in neither of these sites. This is mostly related to the fact that the ingested protons lead often to very energetic nuclear burning through reacting with 12 C, and this energetic feedback will disrupt the He-convection zone as the H ingestion energy forms its own separate convection zone in 1D stellar evolution models (e.g., Cristallo et al. 2009;Herwig et al. 2011). 3D stellar hydrodynamic simulations have confirmed that catastrophic instabilities can be triggered by the H in-gestion (Herwig et al. 2014b). Thus, all of these i-process candidates remain uncertain and not completely supported by simulations, especially when the goal is to reproduce the i-process abundances in CEMP-r/s stars. This includes models of low-mass, low-metallicity AGB star HIF stellar evolution simulations that have obtained peak neutron densities up to Nn ∼ 10 15 cm -3 but too small neutron exposures (Cristallo et al. 2009;Cristallo et al. 2016;Ritter et al. 2018;Karinkuzhi et al. 2021). We will briefly discuss these models in Section 4.

There are, as far as we know, only two exceptions. The H ingestion into the He-core flash of a 1M stellar evolution model with [Fe/H] = -6.5 yield a high enough neutron exposure, so that second-peak elements are produced at that low metallicity (Campbell et al. 2010). However, CEMP-r/s stars that carry the i-process signature in the second-peak region have typically higher Fe abundances ([Fe/H] > -3, e.g. Abate et al. 2016), and it has not been shown yet that the He-core flash can produce second-peak i-process abundances at such higher Fe abundances. This leaves the rapidly-accreting white dwarfs (RAWDs, Denissenkov et al. 2017Denissenkov et al. , 2019) ) the only stellar evolution phase for which self-consistent stellar evolution simulations predict the correct conditions for second-peak i-process abundances as observed in CEMP-r/s stars. Therefore, at this time, the most likely origin of the CEMP-i stars (those CEMP-r/s stars that show a very good match with the i-process predictions) are the RAWDs. They occur in close binary systems where a white dwarf is accreting H-rich material at a rate Ṁacc ∼ 10 -7 M yr -1 . The accreted H is stably burning on the white dwarf surface, which results in an accumulation of a He shell. When its mass has reached a critical value, the He shell experiences a thermal flash that leads to an i-process, like in the case of Sakurai's object. The main differences from Sakurai's object are that the H ingestion in RAWDs can last much longer (nearly a month instead of hours), without being interrupted by a split of the He convective zone, and that, following the expansion of the accreted envelope caused by the He-shell flash and its loss from the system either via a super-Eddington luminosity wind or through a common-envelope interaction, the accretion resumes, and the entire cycle can repeat many times (Denissenkov et al. 2017).

Unlike the s process, the i process involves a large number of neutron-rich isotopes ≈ 2 to 8 mass numbers outside the valley of stability. Only theoretical predictions are available for the (n,γ) reaction rates of these unstable species, e.g. those provided by Rauscher & Thielemann (2000) from Hauser-Feshbach model computations. Sensitivity studies where the impact of (n,γ) rate uncertainties on the final i-process abundances is determined are critical to guide experimental efforts to better constrain (n,γ) rates on unstable nuclei (Larsen et al. 2019;Nunes et al. 2020). A particular high priority are rates that have a strong impact on elemental abundances that are critical for the comparison with observed stellar abundances. Such studies of the impact of (n,γ) reaction rate uncertainties of unstable isotopes on the i-process nucleosynthesis have recently been presented for the first peak elements in Sakurai's object by Denissenkov et al. (2018) and for the elements with 32 ≤ Z ≤ 48 in the metal-poor star HD94028 by McKay et al. (2020) (hereafter, Paper I and Paper II).

In the present work we extend those studies to i-process nucleosynthesis of the elements with 56 ≤ Z ≤ 74, from Ba to W (hereafter referred to as the selected range of elements), using their abundances in some of the CEMP-i stars to constrain the relevant neutron densities. In Section 2, we briefly describe the analysis steps and methods developed in Papers I and II that help to identify both the i-process physical conditions and the unstable isotopes whose (n,γ) reaction rate uncertainties have the strongest impact on the predicted abundances of heavy elements at these conditions. In order to put the results of this impact study in context we discuss in Section 4 the features of low-mass, low-Z AGB stellar models as an alternative site of the i-process nucleosynthesis and their challenges to account for the abundances observed in CEMP-i stars. Sections 3 and 5 summarize our results and conclusions.

We have selected two CEMP-i stars for this study, CS31062-050 and HE2148-1247. The first star has the observed abundances of the elements from Ba to Ir (Aoki et al. 2002;Johnson & Bolte 2004) that were successfully reproduced by Dardelet et al. (2014) with a one-zone model of i-process nucleosynthesis for log 10 Nn(cm -3 ) ≈ 14.5, while the abundance distribution of the elements from Ba to Dy in the second star (Cohen et al. 2003) resembles those observed in the other two stars included in the analysis of Dardelet et al. (2014) which they fitted with the predicted i-process elemental abundances for log 10 Nn(cm -3 ) ≈ 13.5. We begin our analysis using two one-zone nucleosynthesis models that have the same fixed density ρ = 10 4 g cm -3 and temperature T9 ≡ T /10 9 K = 0.2 in which we artificially keep neutron density at the constant values Nn = 3.16 × 10 14 cm -3 and Nn = 3.16 × 10 13 cm -3 corresponding to log 10 Nn(cm -3 ) ≈ 14.5 and log 10 Nn(cm -3 ) ≈ 13.5, respectively.

Our one-zone nucleosynthesis simulations employ the same ppn code from the NuGrid framework (Pignatari et al. 2016) that was used in Papers I and II. It has a small modification that allows to keep Nn constant. For each of the two selected stars, we call a benchmark simulation the one that provides the best fit to its observed heavy-element abundances for the NuGrid default reaction rates the full list of references for which is provided in Paper I. Because these abundances are strongly enhanced compared to the Fe abundance (black circles with error bars in Figures 1 and2), we can assume that their distributions have reached equilibrium states in which abundance ratios for neighbouring elements are mostly determined by nuclear rather than stellar physics. This assumption justifies our use of the one-zone nucleosynthesis model with a constant neutron density.

Results of reaction rate uncertainty studies are based on Monte-Carlo (MC) simulations in which selected reaction rates in the benchmark models are randomly varied within their uncertainty ranges. To minimize the computational time of the MC simulations, we choose the minimum 1, but for the one-zone benchmark model with Nn = 3.16 × 10 13 cm -3 whose predicted elemental abundances are compared with those observed in the CEMP-i star HE2148-1247. In this case, the maximum integration timestep is chosen to be 760 (at the age of 54 days).

possible integration time-steps for the one-zone benchmark models at which their predicted distributions of the heavyelement abundances have approached the assumed equilibrium states observed in the selected stars (curves in Figures 1 and2). The predicted abundance distributions in Figures 1 and2 have been scaled to the observed [Ce/Fe] ratio using the pinning method from Appendix A1 of Paper II. The first step of our analysis is the adjustment of the constant neutron density and the integration time-step for our one-zone benchmark nucleosynthesis model. The initial chemical composition for the simulations is the Asplund et al. (2009) solar abundance mixture scaled to the metallicity [Fe/H] = -2.6, which is close to the metallicities of the two CEMP-i stars, with abundances of α-elements enhanced to the mean value of [α/Fe] = +0.4. However, unlike the weak i process studied in Paper II, the results in this work do not depend strongly on the initial composition as the abundance distribution reaches equilibrium. We suppress additional production of neutrons in the reactions 12 C(p,γ) 13 N(e + ν) 13 C(α,n) 16 O by setting the initial hydrogen abundance to X0(H) = 0, with its mass fraction being distributed among C and O with X0(C) = 0.5.

To validate the one-zone benchmark models, we have included in our analysis a multi-zone i-process nucleosynthesis model based on the convective He-shell structure of the RAWD model G (Denissenkov et al. 2019) with the same metallicity [Fe/H] = -2.6. The abundances computed for this model using the NuGrid code mppnp (Pignatari et al. 2016) also fit all of the abundances in the selected element range in the star CS31062-050 well, except Ba and Pr whose abundances are under-and over-produced, respectively (Figure 12 in Denissenkov et al. 2019). The Ba discrepancy can be caused by the fact that the maximum neutron density in the He convective zone of the RAWD model G only reaches a value of log 10 Nn,max (cm -3 ) ≈ 13.5 (blue curve in Figure 3), which is close to the value of Nn in Figure 2 with a lower Ba abundance. However, because Nn,max is proportional to the H-ingestion rate, a value of log 10 Nn,max (cm -3 ) ≈ 14.5, like Nn in Figure 1 with a higher Ba abundance, could be reached if that rate in the i-process pollution source of the star CS31062-050 were ten times higher than in the RAWD model G (orange curve in Figure 3).

We cannot exclude a possibility that in some RAWDs the H-ingestion rate Ming is higher or lower than in the RAWD model G, given the uncertainties of the H-ingestion rate estimates and the lack of a detailed parameter-space study of RAWD models. Figure 5 in Denissenkov et al. (2019) compares the values of Ming estimated for six 1D stellar evolution RAWD models of nearly equal masses MWD ≈ 0.73M , central temperatures log 10 Tc ≈ 7.2 and mass accretion rates Macc ≈ 1.6×10 -7 M yr -1 but different metallicities, -2.6 ≤ [Fe/H] ≤ -0.7, with the dependence of Ming on the He-shell luminosity LHe obtained from the 3D hydrodynamic simulations of H ingestion that used the solar-metallicity 0.65M RAWD model A from Denissenkov et al. (2017) for the simulations setup. It shows that for the values of LHe at the beginning of H ingestion five of the six 1D H-ingestion rates agree with their extrapolated 3D counterparts within factors of less than 2.3 and only for the RAWD model B with [Fe/H] = -0.7 the difference in 1D and 3D values of Ming reaches a factor of 4.4. From this comparison we infer that although the maximum difference in Ming for the six RAWD models is as large as a factor of ∼ 10 it is unlikely to strongly depend on the initial metallicity of the WD progenitor, because neither Figure 5 nor Table 1 of Denissenkov et al. (2019) reveal any systematic dependence of Ming on [Fe/H]. On the other hand, like the peak H-burning temperatures and luminosities in novae (e.g., Prialnik & Kovetz 1995), the He-shell luminosity and Ming that is nearly proportional to LHe in the RAWD models may be higher for larger MWD and lower Tc. Therefore, when considering the case of log 10 Nn,max (cm -3 ) ≈ 14.5 we assume that this higher neutron density could be attained in a more massive and probably initially less luminous (colder) RAWD than the RAWD model G.

At the second step, we select the unstable isotopes whose (n,γ) reaction rate variations are expected to affect the predicted i-process elemental abundances. The mass fractions X of stable and unstable isotopes and reaction fluxes f = (δY /δt), where Y = X/A for an isotope with atomic mass A, are shown in Figures 4 and5 for the one-zone benchmark model from Figure 1. We use these figures to identify 164 unstable isotopes, from 131 I to 189 Hf, that are involved in the synthesis of the 18 stable elements in the range from Ba to W, selected for the analysis in this work. Variations of (n,γ) reaction and beta-decay rates of the identified unstable isotopes will change the strengths of the reaction fluxes on the path map of the i process producing the selected elements and, as a result, their predicted abundances.

At the third step, for each of the identified unstable isotopes a set of theoretical values of its (n,γ) rate ri is calculated with the Hauser-Feshbach computer code TALYS 3 (Bersillon et al. 2007). These calculations are done at a temperature of T9 = 0.3 intermediate between the constant T9 = 0.2 in the one-zone model and the maximum T9 ≈ 0.35 in the multi-zone model using 20 different combinations of the nuclear level density and γ ray strength models from the references provided in Table 1 of Denissenkov et al. (2018). Each of these 20 numbers represents a possible, theoretically predicted, value of the rate, therefore we can estimate 3 http://talys.eu its uncertainty as v max i = r max i /r min i , which we define as the maximum variation factor (Figure 6). As was shown previously (Liddick et al. 2016), rate uncertainties from Hauser-Feshbach models increase from about a factor of 2 near stability to factors of 10 a few neutrons away from stability. The default values of ri in the ppn code represent other possible values, that are taken from Rauscher & Thielemann (2000). We assume that the true value of the rate is, with equal probability, somewhere between its default value divided and multiplied by v max i . To study the impact of these uncertainties on the predicted abundances in the selected element range, we perform an MC simulation with 10000 runs of the ppn code in which the default rates of the 164 unstable isotopes from Figure 6 are multiplied by the factors fi

, where p is assigned a value 0 or 1 with equal probability, and v rand i is randomly chosen from a uniform distribution between 1 and v max i .

To investigate the impact of uncertainties of temperature-dependent beta-decay rates of the same 164 isotopes on the simulated i-process nucleosynthesis, we use Equation ( 7) from Rauscher et al. (2016) that gives us estimates of their corresponding maximum variation factors in the following form:

where the ground-state and theoretical uncertainties are set to β Ug.s. = 1.3, β U th = 10, and the partition functions of the isotopes normalized to their ground state spins G(T ) = G(T )/(2J0 + 1)) are taken for the temperature T9 = 0.3 from Table V of Rauscher & Thielemann (2000).

The ground state beta decay rates of interest in this work are experimentally known and β Ug.s. = 1.3 is therefore an overestimation of the uncertainty. However, as we show later, even these conservative uncertainties are negligible.

Finally, at the fourth step we do a correlation analysis of results of the MC simulation. The MC results are m = 10000 different sets {X k,j } m j=1 of predicted mass fractions X k of selected elements. Because in the cases of interest here, unlike for Sakurai's object studied in Paper I, the i-process elements had probably been accreted by the observed stars a long time ago, we allow all unstable isotopes to decay for 1 Gyr at the end of our nucleosynthesis simulations. Therefore, Tc and Pm, that do not have stable isotopes, are not included in the analysis. The MC abundance sets differ from each other because they are obtained using different sets of reaction rates {fi,jri} n i=1 , here for the number n = 164 of the unstable isotopes from Figure 6. To find out which of the reaction rate variations fi has the strongest impact on the predicted abundance X k , we calculate the Pearson productmoment correlation coefficients

where fi = ( m j=1 fi,j)/m and X k = ( m j=1 X k,j )/m. The benchmark simulation (the one with the NuGrid default reaction rates) is assigned the index j = 0 and it uses the multiplication factors fi,0 = 1 for all reaction rates. 115 116 117 118 119 120 121 122 123 124 125 126 127 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 115 116 117 118 119 120 121 122 123 124 125 126 127 128 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 116 117 118 119 120 121 122 123 124 125 126 127 128 129 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174

In this section we will first present the main results of our reaction rate uncertainty studies based on the one-zone models from Figures 1 and2 for radiative neutron captures and beta decays (the latter only for the second model). Then we will summarize the results of our analysis of the multi-zone MC simulation based on the RAWD model G that reproduces the heavy-element abundance distribution in the CEMP-i star CS31062-050 nearly as satisfactorily as the second of our one-zone models, i.e. without matching the observed Ba enhancement (Figure 12 in Denissenkov et al. 2019).

Table 1. The strongest correlations between the (n,γ) reaction rate variations and the i-process elemental abundances found in the MC simulation based on the one-zone benchmark model with Nn = 3.16 × 10 14 cm -3 . Only up to two correlations with |r P (f i , X k /X k,0 )| ≥ 0.15 are shown for which the corresponding Spearman correlation coefficients r S are also provided. For this model, we have randomly varied only (n,γ) rates. Before calculating and analyzing correlation coefficients, we plot distributions of abundances predicted in the MC simulation for each of the elements in the selected range of Z (Figure 7). To make this plot we have used the dilution method from Appendix A2 of Paper II where we demonstrated that the pinning and dilution methods, the former having been used in Figures 1 and2, can be equivalently used in this case where the overabundances of the n-capture elements in question are so large (typically ≥ 2dex). Most elements respond to reaction-rate variations with an asymmetric distribution, in several cases even double-peaked. These are signatures of their abundances being strongly affected by variations of individual (n,γ) rates. To find out what reactions have the strongest impact on a given elemental abundance, and therefore are likely to be responsible for its asymmetric distribution, we display in Table 1 one or two reactions with the maximum Pearson correlation coefficients, such that |rP(fi, X k /X k,0 )| ≥ 0.15, for each of the elements in the selected range of Z, where X k,0 are their abundances predicted by the benchmark model. The lower limit 0.15 for the correlation coefficient allows us to reveal most of the (n,γ) reactions responsible for asymmetric distributions of the predicted elemental abundances. In the fourth column of Table 1 we provide Spearman correlation coefficients. The fact that for most of the selected elements the Spearman coefficients, examining the monotonicity of the found correlations, exceed by the absolute magnitude their corresponding Pearson coefficients, the latter estimating the linearity of the correlations, tells us that the correlations are not linear but significant. Finally, to check that our results are not very sensitive to the choice of the characteristic temperature at which the i-process is simulated in the one-zone model, we have changed it from T9 = 0.2 to T9 = 0.3, and the Pearson coefficients for the latter case are presented in the last column of Table 1. It is seen that they insignificantly differ from the Pearson coefficients for the case of T9 = 0.2. Table I is our main result and identifies the most important neutron capture rate uncertainties affecting each of the produced elemental abundances. The largest correlation coefficients found for each element identify the unstable isotopes whose reaction rate variations within their adopted uncertainties have the strongest impact on this element's predicted abundance.

The thick black lines in Figure 8 are the distributions of the Ba and Pr abundances obtained in the MC simulation based on the one-zone benchmark model with Nn = 3.16 × 10 14 cm -3 . They are compared with the abundances of Ba and Pr reported for the star CS31062-050 by Johnson & Bolte (2004) that are represented by the greenshaded Gaussian curves for which the mean values and standard deviations are the same that we used in Figure 7 to plot black circles with errorbars. The double-peak distribution of the Ba abundances in the MC simulation can be separated into two distinct single-peak distributions (the blue and yellow histograms) if we divide the results of the MC simulation into two groups, the one with the multiplication factors fi > 1 and the other with fi < 1 for the rate of the reaction 135 I(n,γ) 136 I whose variation has the strongest correlation with the predicted Ba abundance (Table 1, see also Bertolli et al. 2013 andCristallo et al. 2016). A similar result is obtained for Pr when we use the rate of the reaction 141 Ba(n,γ) 142 Ba identified for it in Table 1. We see that the Ba abundance predicted with the benchmark model (the vertical dashed line) agrees very well with the observed Ba abundance. However, this is not true for Pr. In this case the higher than the default rate of the reaction 141 Ba(n,γ) 142 Ba could help to significantly reduce the discrepancy between its predicted and observed abundances. The single-peak distributions, blue and yellow histograms, correspond to rate multiplication factors f i > 1 and f i < 1 for the reactions 135 I(n,γ) 136 I (top panel) and 141 Ba(n,γ) 142 Ba (bottom panel), as identified by their strongest correlation coefficients in Table 1 for Ba and Pr, respectively.

zone MC simulations as functions of the (n,γ) reaction rate multiplication factors fi with which they have the strongest correlations. The Pearson product-moment correlation coefficients in the four panels of these figures show how well compared data (anti-)correlate, but unfortunately they do not provide quantitative estimates of the slopes of the correlations that would also be useful to have. A comparison of these figures with similar ones for other elemental abundances shows that distinctly visible slopes in correlations like these appear only when |rP| > ∼ 0.4.

2 One-zone model with Nn = 3.16 × 10 13 cm -3

This benchmark model has the neutron density adjusted to reproduce the abundances of the elements in the selected range of Z in the star HE2148-1247 (Figures 2 and10). However, the Nn = 3.16 × 10 13 cm -3 model also reproduces most of the abundances in CS31062-050 (Figure 11). The exceptions are Ba and Pr, and the Ba case is the sole reason a higher neutron density was assumed for that star. However, the observed Ba abundance has large systematic errors that signficantly exceed the error bars shown. . The distributions of the Ba, Ce, Pr, and Lu elemental abundances predicted in the one-zone MC simulations with Nn = 3.16 × 10 14 cm -3 as functions of the reaction rate multiplication factors f i with which they have the strongest correlations. The small differences between the corresponding correlation coefficients in this plot and in Table 1 are caused by the replacement of f i with log 10 f i in the correlation analysis. 7, but for the one-zone benchmark model with Nn = 3.16 × 10 13 cm -3 compared with the observed abundances in the CEMP-i star HE2148-1247.

nificantly lower than the [Ba/Fe] = 2.8 from the same group published in Johnson & Bolte (2004) that we use here. This lower [Ba/Fe] value would be in agreement with the prediction of our lower neutron density model. Therefore, in addition to the nuclear uncertainties, observational uncertainties need to be addressed to constrain the neutron density more stringently. With the current uncertainties, both neutron density models can be considered potential i-process sources for CS31062-050. However, the significant decrease of Nn between the first and the second of our one-zone i-process nucleosynthesis simulations changes the i-process band path and therefore the important (n,γ) reactions. This is reflected in Table 2 that for most elements lists different (n,γ) reactions having the strongest impact on their predicted abundances than Table 1. For example, instead of 135 I, it is now 137 Cs that has the strongest impact on the Ba abundance, and because the benchmark value of [Ba/Fe] is now much lower than the one reported for the star CS31062-050 by Johnson & Bolte (2004), a lower rate of the reaction 137 Cs(n,γ) 138 Cs would reduce this discrepancy (the yellow histogram in the left panel of Figure 12). For Pr, it is now an increase of the 141 La(n,γ) 142 La reaction rate that is required to diminish the discrepancy between its predicted and observed abundances (the blue histogram in the right panel of Figure 12).

For this benchmark model we have also done an MC simulation in which only beta-decay rates of the 164 unstable isotopes from Figure 6 have been varied. The results of this simulation plotted in Figure 13 show that these variations have a much weaker effect on the predicted abundances of the elements in the selected range of Z as compared to the variations of the (n,γ) rates for the same isotopes. This justifies our decision to not include them into the present reaction rate uncertainty studies.

If we used the elemental abundances predicted by the one-zone benchmark model with Nn = 3.16 × 10 13 cm -3 for the CEMP-i star HE2148-1247 from Figure 2 for comparison with the abundances in the star CS31062-050 we would find discrepancies between the predicted and observed abundances of Ba and Pr for the second star. A comparison of the surface abundances of Ba and Pr predicted for the two values of Nn in Figures 1 and2 reveals that both elements have enhanced abundances in the case with the higher Nn. Therefore, the higher neutron density could eliminate the Ba abundance discrepancy for the star CS31062-050 in the case of the lower Nn, and the assumption of a higher rate for the reaction 141 Ba(n,γ) 142 Ba would significantly reduce the Pr abundance discrepancy at the same time. The latter assumption would not affect the predicted Pr abundance for the star HE2148-1247 because at the lower neutron density, that better matches its elemental abundance distribution in the selected range of Z, [Pr/Fe] has the strongest anti-correlation with a different reaction, namely 141 La(n,γ) 142 La. 9, but for the one-zone model with Nn = 3.16 × 10 13 cm -3 .

Multi-zone benchmark models of the i-process nucleosynthesis better simulate the physical conditions at which it is believed to occur in stars. So far, only two such models have been investigated -for the H ingestion by the He-shell convection in Sakurai's object (Herwig et al. 2011) and in RAWDs (Denissenkov et al. 2017(Denissenkov et al. , 2019)). Here, we use the RAWD model G from the latter paper to simulate the i-process nucleosynthesis of the elements in the selected range of Z in the He convective zone and to study its reaction rate uncertainties. In this model, neutrons are produced near the bottom of the He shell in the reaction 13 C(α,n) 16 O for which fresh 13 C is supplied by the positron decay of 13 N that is made in the reaction 12 C(p,γ) 13 N in the middle of the convective zone using protons entrained from the H-rich envelope. The evolution of the maximum neutron density in the He zone of this model is shown in Figure 3 (the blue curve), and the results of the MC simulation based on this multi-zone benchmark model are plotted in Figure 15 and listed in Table 3. By comparing Figures 11 and15 and Tables 2 and3, we conclude that, in spite of its simplicity compared to the multi-zone model, the one-zone model with the constant neutron density Nn = 3.16 × 10 13 cm -3 close to the peak value of Nn,max in the multi-zone model can reliably predict both, distributions of the elemental abundances in the MC simulation and the identified key reactions for most of the elements (the only exception is Sm, for Tb, Tm, Lu, and W the reactions with the largest and the second largest correlation coefficient having swapped their places). This means that one-zone MC simulations can be used instead of much more computationally time-consuming multi-zone MC simulations for i-process reaction rate uncertainty studies, provided that the former use neutron densities similar to the maximum ones found in the latter.

Our paper is mainly focused on identifying the (n,γ) reaction rates whose uncertainties have the strongest impact on the predicted abundances of the elements from Ba to W for the i-process neutron densities and necessary neutron exposures that are encountered without any ad-hoc assumptions in RAWD stellar evolution models and that agree with elemental abundances observed in the CEMP-i stars CS31062-050 and HE2148-1247. A potential alternative i-process site are low-mass and low-Z AGB thermal pulse stars (Cristallo et al. 2016;Choplin et al. 2021). The fact that H ingestion can take place in low-Z AGB stars is well established, and that i-process conditions can be achieved is also suggested by several models (for another example see Ritter et al. 2018, Sections 3.2.1 and 4.6). The question is if i process in low-Z AGB stars can produce abundance patterns observed in CEMP-i stars. Karinkuzhi et al. (2021) concluded, based on comparing observations with models of Choplin et al. (2021), that AGB stars with M = 1 M and -3 ≤ [Fe/H] ≤ -2 are sufficient in explaining i-process abundances in CEMP-i stars and that RAWDs may not be needed for this. The challenge of reproducing the i-process abundance patterns in CEMP-i stars is that models need to satisfy two key observational properties of CEMP-i stars. First, they must reach a high enough neutron density to produce the ratios of abundances of nearby elements at and just beyond the second peak, most notably [Ba/Eu]. Second, they must achieve high enough neutron exposures in order to reach the high hs/ls ratios of CEMP-i stars. As already noted by Cristallo et al. (2016), unless one makes ad-hoc assumptions on continued convective mixing, when standard stellar evolution predicts a split of the convection zone due to the energy input from convective-reactive 12 C(p,γ) burning, the H ingestion and the simultaneous burning and mixing cannot be sustained for long enough to continue the convective-reactive generation of high neutron densities.

RAWD models solve both of these challenges without any ad-hoc mixing assumptions (Denissenkov et al. 2017(Denissenkov et al. , 2019)), and it is easy to see why. When the convection zone splits in low-Z AGB HIF simulations as a result of the H in- gestion this will essentially end the mixing conditions needed for i process, because the 12 C(p,γ) reaction that creates the neutron source 13 C is separated from the hot conditions where the 13 C(α,n) reaction releases neutrons fast enough to result in high neutron densities. Thus, when the split occurs (at least the way it occurs in standard 1D models) the i-process conditions are quenched. The split happens when the energy released by the 12 C(p,γ) reaction produces locally more entropy than can be carried away through convection driven by the underlying He burning. Then, at the location of the H burning inside the He-shell convection zone the entropy added by H burning will cause an entropy step as seen in stellar evolution models and in 3D hdyrodynamic simulations, such as Figure 3 in Herwig et al. (2014b). The entropy step of course implies a radiative separation layer between the hot He-burning layer below and the cooler H-burning layer above, constituting the split that occurs roughly when LH > ∼ LHe (for example Herwig et al. 1999b, Figure 4). Of course LH depends on the entrainment rate of H which in turn depends on the He luminosity driving the convection (see Figure 5 in Denissenkov et al. 2019). RAWD thermal pulses have a key difference compared to AGB thermal pulses where the RAWD mass is the same as the AGB core mass. Due to their formation history the WD in a rapidly accreting system is cold and therefore the He-shell flashes are stronger. As a result, the entropy produced from burning of entrained H can continuously be distributed accross the convection zone due to the more efficient convection driven by higher LHe, although the entrainment rate and therefore LH is also somewhat larger. As a result, He-shell flashes in RAWDs have generally no split of the convection zone, and maybe just very weak events trailing the prolonged main H ingestion such as shown in Figure 3 of Denissenkov et al. (2017). This absence of a split during H-ingestion phases is a natural behaviour of RAWD models and does not require any forcing or ad-hoc assumptions and it leads to models with high neutron density and high neutron exposure.

The formation of the split in the models used by Karinkuzhi et al. (2021) is well documented in Figure 2 of Choplin et al. (2021). As explained in the figure caption the time span between when the peak neutron density is reached and the time of the split is 17.5 hr. As a side we note that this is about or even less than the convective turn-over time of the He-shell flash convection zone in a 1 M low-Z AGB thermal pulse according to the NuGrid model database (Figure 15 in Ritter et al. 2018). However, in the stellar evolution models the energy transport and mixing properties of a convective region are approximated in 1D with the mixing length theory (MLT, Cox & Giuli 1968) that is based on averaged quantities over space and time. Certainly, time averages to describe a convective phenomenon that lasts just one convective turn-over become a poor approximation of the 3D hydrodynamic response to a dramatic energy input in a HIF that, when integrated over a convective turn-over time, typically reaches a substantial fraction of the binding energy of the convection zone (Jones et al. 2016;Ritter et al. 2018;Clarkson & Herwig 2020). But if the stellar evolution simulation is to be taken seriously the split will generally end the i-process nucleosynthesis without additional ad-hoc assumptions (Cristallo et al. 2016). This split will therefore limit the neutron exposure and along with it the hs/ls ratio that can be reached in low-Z AGB i-process models.

This deficiency is quite evident from the model to observations comparison in Karinkuzhi et al. (2021). Our Figure 16 is made to resemble Figure 13 of Karinkuzhi et al. (2021), but we keep in it only their observed CEMP-r/s stars and abundance ratio dilution curves of the 1 M AGB model with [Fe/H] = -2.5 because this is close to both the metallicities of our RAWD F and G models, [Fe/H] = -2.3 and [Fe/H] = -2.6, respectively, and the mean [Fe/H] ratio of the CEMP-i stars. For a comparison we have plotted the dilution curves of our RAWD G model (they are the same for the RAWD F model) and we have also added the stars CS31062-050 and HE2148-1247 studied in this paper. The elemental abundance ratios plotted along the vertical axes measure the relative heights of the second-and first-peak regions, while the [La/Eu] ratio plotted along the horizontal axes is a standard diagnostic, similar to [Ba/Eu], used to classify a CEMP star either as CEMP-s, for [La/Eu] > 0.5, or CEMP-r/s, for 0 < [La/Eu] < 0.5 (Beers & Christlieb 2005). This comparison of the low-Z AGB and RAWD models shows the effect of the split limiting the neutron exposure and thus the maximum hs/ls (such as Ba or La over Y or Zr) in the low-Z AGB model vs the much higher neutron exposure and hs/ls values in agreement with observations that can be reached in RAWD models.

We also agree with Choplin et al. (2021) that their models have a too small dilution, i.e. a lot of low-Z AGB material would have to be accreted. From their equation ( 10), the ratio of the accreted mass with i-process yields to the envelope mass of the accreting star is Macc/Menv = 1/f -1, where f is the dilution coefficient. For their f ≈ 0.2 Macc/Menv ≈ 4 which is impossible if Menv is the mass of the convective envelope of a CEMP-i red giant (11 out of 14 of their considered CEMP-i stars are red giants). For RAWD values of f ≈ 0.996 Macc/Menv ≈ 0.004 which instead represents a more reasonable value.

In view of these challenges the conclusion of Karinkuzhi et al. (2021) that i process responsible for the abundance patterns observed in CEMP-i stars takes place in low-mass, low-metalliticy AGB stars, and that other sites, such as the RAWDs, may not be needed is premature. This is in agreement with Cristallo et al. (2016) who were the first to report heavy-element abundance distributions obtained as a result of n-capture nucleosynthesis at the peak value of Nn ∼ 10 15 cm -3 reached during the PIE in the He convective zone of a low-mass (1.3M ) low-Z ([Fe/H] = -2.85) AGB model experiencing its first thermal pulse that, unlike the AGB star model of Karinkuzhi et al. (2021), then followed the standard AGB evolution with further TPs, s-process nucleosynthesis, and third dredge-ups. Just as demonstrated in our Figure 16, Cristallo et al. (2016) concluded that the PIEs in AGB stars had unlikely been the source of the i-process abundances in CEMP-i stars because their models could not attain the high second-to first-peak n-capture elemental abundance ratios observed in these stars.

We will close this section with an interesting twist to this argument which just demonstrates that the actual site of the i process for CEMP-i stars remains still not entirely settled. A similar situation to the question of whether low-Z AGB stars can or cannot account for CEMP-i star abundances exists for Sakurai's object. Herwig et al. (2011) showed that post-processing simulations of the plain stellar evolution model of the very-late thermal pulse born-again star would not provide high enough neutron exposures to reach in this case the very low hs/ls. In Sakurai's object only the first-peak elements are enhanced. Sakurai's object is a nearby Pop I object. Since the initial Fe abundance is much higher than in CEMP-i stars a similar neutron exposure is expected to reach only to the first peak. However, due to the problem of the split, just as in low-Z AGB models, the presently favoured post-AGB very-late thermal pulse model, if post-processed as is, cannot reproduce the observed abundances. Herwig et al. (2011) showed that if mixing accross the split was assumed to continue for 1000 min then the observed abundances could be reproduced. 3D hydrodynamic simulations (Herwig et al. 2014a) confirmed the early split predicted by 1D stellar evolution. It is interesting to note that RAWD models of high Z similar to that of Sakurai's object do an excellent job of explaining Sakurai's abundance pattern, just as the low-Z RAWD models do an excellent job in explaining CEMP-i stars. From a nucleosynthesis point of view Sakurai's object is favoured to be a Pop I RAWD. If Sakurai's object is instead a born-again AGB star then our 1D stellar evolution models are fundamentally unable to predict the occurrance and properties of the split correctly, and this would then also apply to the low-Z AGB stellar evolution models.

For the neutron densities Nn = 3.16 × 10 14 cm -3 and Nn = 3.16 × 10 13 cm -3 that were presumably attained in the i-process stellar sites that polluted the CEMP-i stars CS31062-050 and HE2148-1247, respectively, we have identified the (n,γ) reaction rates of the 164 unstable isotopes selected for this study (Figure 6) whose variations fi within the ranges constrained by the Hauser-Feshbach model computations have the strongest impact on the predicted mass fractions X k of 18 elements from Ba to W. This has been done by calculating the Pearson product-moment correlation coefficients rP(fi, X k /X k,0 ), where X k,0 are the abundances predicted by the benchmark i-process nucleosynthesis models for which all fi = 1. Up to two maximum coefficients with |rP| > 0.15 are presented in Tables 1, 2, and 3. For each of the 18 selected elements, they identify one or two key reactions whose rate variations within their theoretical uncertainty limits have the strongest impact on its predicted abundance and whose rates need therefore to be measured experimentally first to improve the predictive power of stellar i-process nucleosynthesis models. Among those, experimental measurements of the rates of the reactions 135 I(n,γ) 136 I, 137 Cs(n,γ) 138 Cs, 141 Ba(n,γ) 142 Ba, and 141 La(n,γ) 142 La have the high priority because these rates have the strongest impact on the predicted abundances of Ba and Pr that have been observed in the considered CEMP-i stars and these elements are particularly sensitive to neutron density. Our MC simulations provide in addition statistically meaningful nuclear uncertainties for the i-process abundance predictions. This enables a quantitative comparison with observations, e.g. compare the predicted and observed abundances of Ba and Pr in Figure 11 where their error bars are shown. As Figure 7 shows, with the currently large nuclear uncertainties our i-process model predictions agree with observations, except for Pr, and, in the cases of low Nn (Figures 11 and15) for Ba.

We have shown that the uncertainties of the temperature-dependent β-decay rates of the same 164 unstable isotopes have a much weaker effect on the predicted abundances and that, at least, uncertainty studies of the elements from Ba to W can rely on one-zone MC simulations, instead of computationally time-consuming multi-zone MC simulations, provided that the former use values of Nn consistent with those attained in the latter.

In order to underpin the astrophysical context of this impact study we have discussed the question whether recently proposed low-mass low-Z AGB star models (Karinkuzhi et al. 2021;Choplin et al. 2021) could be an alternative site of the i-process nucleosynthsis for the heavyelement abundance enhancements in CEMP-i stars. We have shown that the ratios of the second-to the first-peak ncapture elemental abundances predicted by these AGB models are not as high as those observed in CEMP-i stars that, on the other hand, are well reproduced by our low-Z RAWD models. Therefore, in agreement with the results of Cristallo et al. (2016), we conclude that the low-mass low-Z AGB stars are not as likely candidates for the source of the i-process abundances for CEMP-i stars as the RAWDs. However, the treatment of convective-reactive H-ingestion flash events is not well approximated by 1D stellar evolution, and future 3D hydrodynamic simulations will have to reveal the true nature of the formation and evolution of the He-shell convection split and the subsequent convective mixing and nucleosynthesis.

The computational and analysis tools developed during our work on Papers I, II and this paper can be used to extend our analysis to other neutron densities and to other reactions and elements.

Figures with correlations like those shown in Figure 9 for all of the elements listed in Table 1 can be found at https: //doi.org/10.5281/zenodo.4148667. Other data underlying this article can be shared on reasonable request to the corresponding author.