# Bayesian Estimation of the $S$ Factor and Thermonuclear Reaction Rate for
16O(p,$\gamma$)17F

Christian Iliadis  Vimal Palanivelrajan  Department of Physics & Astronomy,
University of North Carolina at Chapel Hill, NC 27599-3255, USA
Triangle Universities Nuclear Laboratory (TUNL), Durham, North Carolina 27708,
USA  Rafael S. de Souza  Key Laboratory for Research in Galaxies and
Cosmology, Shanghai Astronomical Observatory,
Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, China

###### Abstract

The 16O(p,$\gamma$)17F reaction is the slowest hydrogen-burning process in the
CNO mass region. Its thermonuclear rate sensitively impacts predictions of
oxygen isotopic ratios in a number of astrophysical sites, including AGB
stars. The reaction has been measured several times at low bombarding energies
using a variety of techniques. The most recent evaluated experimental rates
have a reported uncertainty of about 7.5% below $1$ GK. However, the previous
rate estimate represents a best guess only and was not based on rigorous
statistical methods. We apply a Bayesian model to fit all reliable
16O(p,$\gamma$)17F cross section data, and take into account independent
contributions of statistical and systematic uncertainties. The nuclear
reaction model employed is a single-particle potential model involving a
Woods-Saxon potential for generating the radial bound state wave function. The
model has three physical parameters, the radius and diffuseness of the Woods-
Saxon potential, and the asymptotic normalization coefficients (ANCs) of the
final bound state in 17F. We find that performing the Bayesian $S$ factor fit
using ANCs as scaling parameters has a distinct advantage over adopting
spectroscopic factors instead. Based on these results, we present the first
statistically rigorous estimation of experimental 16O(p,$\gamma$)17F reaction
rates, with uncertainties ($\pm 4.2$%) of about half the previously reported
values.

††preprint: APS/123-QED

##  I Introduction

The 16O(p,$\gamma$)17F reaction ($Q$ $=$ $600.27\pm 0.25$ keV [1]) is the
slowest process among all proton-induced reactions in the CNO mass region [2,
3]. The lowest-lying resonance is located at relatively high laboratory energy
of $\approx$2.7 MeV [4]. Below this energy, the 16O(p,$\gamma$)17F reaction is
a prime example of the nonresonant direct radiative capture process, which
assumes that the proton is captured via a single-step process into a final-
state orbit outside a closed 16O core [5, 6]. This reaction has been measured
many times at low bombarding energies using a variety of techniques, including
the activation method, in-beam detection of prompt $\gamma$ rays, and
experiments in inverse kinematics. A comprehensive analysis of the most
reliable data has been presented in Refs. [7, 8]. The reported thermonuclear
reaction rate uncertainty in Ref. [7] is about 7.5% at temperatures below $1$
GK. Knowledge of the rate at a few-percent uncertainty level is desirable
because the 16O(p,$\gamma$)17F reaction influences sensitively the 17O/16O
abundance ratio and, to a lesser degree, the 18O/16O ratio. This information
directly impacts the interpretation and paternity of oxygen isotopic ratios
measured in presolar stardust grains [9, 10, 7, 11].

Previous evaluations of the 16O(p,$\gamma$)17F rate [12, 7, 2] were performed
with methods that were conventionally employed at the time. The data from
different experiments were fitted independently because it was not clear,
within the $\chi^{2}$ method used, how to perform a common fit across several
different data sets. Also, it was neither clear how to treat independent
contributions from statistical and systematic uncertainties, nor how to
combine total cross section data in the analysis with data on individual
transitions. For these reasons, a number of ad hoc assumptions were made by
Ref. [12, 7, 2] that were not rigorous in a statistical sense.

The advent of thermonuclear rates based on hierarchical Bayesian models has
improved this situation significantly. The method was first presented in Ref.
[13] and subsequently applied to reactions of interest to Big Bang
Nucleosynthesis (BBN). In the simplest cases, the Bayesian models employed
either polynomial fitting functions or predictions from microscopic nuclear
reaction models [14, 15, 16]. The method was extended in Refs. [17, 18] to
implement a one-level R-matrix approximation into the Bayesian fitting. Here,
we report on the implementation of a single-particle potential model into the
Bayesian framework.

As will be seen below, the hierarchical Bayesian model solves a number of
problems that plagued previous work (see, e.g., Ref. [7]): (i) it allows for
the straightforward implementation of the total cross section data of Ref.
[19] and the single datum of Ref. [20] for the first-excited-state transition
(both disregarded previously); (ii) it facilitates a combined fit of all data
(results from different experiments were previously fitted independently);
(iii) it fully accounts for the independent contributions of statistical and
systematic uncertainties (which were previously combined into a single
uncertainty); (iv) it makes no ad hoc assumptions on how to combine fits from
different data sets or $\gamma$-ray transitions.

A primary goal of this work is to include the parameters of the single-
particle potential model in the random sampling, i.e., we will be exploring
the sensitivity of the $S$ factor to these quantities. It will be demonstrated
that the uncertainty of the fitted $S$ factor is relatively large when the
Bayesian analysis is performed using spectroscopic factors, but is
significantly reduced when asymptotic normalization coefficients (ANCs) are
employed instead.

The data selection is discussed in Sec. II. The nuclear reaction formalism is
given in Sec. III and the hierarchical Bayesian model is presented in Sec. IV.
Fits of the $S$ factor are found in Sec. V and thermonuclear rates calculated
in Sec. VI. Section VII provides a concluding summary.

##  II Data selection

Consistent with previous Bayesian reaction rate estimates [13, 14, 17, 15, 18,
16], we will consider only data for which statistical and systematic
uncertainties can be estimated separately. This is the case for four data
sets: Hester, Pixley and Lamb [19]; Chow, Griffiths and Hall [21]; Becker et
al. [20]; and Morlock et al. [22]. The first work [19] only measured the total
cross-section, while the third [20] reported only the cross section for the
transition into the first-excited state ($E_{x}$ $=$ $495.33\pm 0.10$ keV [4])
at a single bombarding energy. The data of Refs. [19, 20] were not taken into
account in the previous 16O(p,$\gamma$)17F rate evaluation [7], because, at
the time, it was neither clear how to fit the total cross section [19]
together with those for the individual transitions, nor how to reliably
include a data set consisting of a single data point only [20]. As will be
seen in Sec. IV, all of these data can be rigorously included in a
hierarchical Bayesian model. We discuss below the four data sets individually.

Hester, Pixley and Lamb [19] measured the total cross-section of the
16O(p,$\gamma$)17F reaction at six center-of-mass energies between $132$ and
$160$ keV. These represent the lowest-energy data points among all the data
sets taken into account in the present work. The reported statistical
uncertainties range from 14% to 40% for the highest and lowest energy,
respectively. The cross sections have been corrected using modern stopping
powers, as discussed in [7], and we adopt these corrected values in the
present work. From their quoted uncertainties in the measured beam current
(6%), counter efficiency (7%), and stopping power (10%), we estimate a total
systematic uncertainty of 14%.

Chow, Griffiths and Hall [21] measured the cross section for the transition to
the ground state at four center-of-mass energies between $1288$ and $2404$
keV, and for the transition to the first-excited state at seven energies
between $795$ and $2404$ keV. The statistical uncertainties range from 3% to
12%. The systematic effects in their measurement were discussed by Ref. [7],
including $\gamma$-ray efficiency (3%), escape peak detection (1%), angular
distribution correction (1%), and effective bombarding energy (3%).
Consequently, we adopt a value of 5% for the combined systematic uncertainty.

Becker et al. [20] measured the 16O(p,$\gamma$)17F reaction in inverse
kinematics at a single center-of-mass energy of $853$ keV. Although not
mentioned explicitly, their reported cross section refers to the transition to
the first-excited state only. The statistical uncertainty amounts to 13%. The
main sources of systematic uncertainty arise from the strength of their
adopted standard resonances in 19F(p,$\alpha_{2}$)16O and the $\gamma$-ray
efficiency in their extended gas target. We estimate an overall systematic
uncertainty of 5% for their reported cross section.

Finally, Morlock et al. [22] reported 16O(p,$\gamma$)17F cross sections below
a center-of-mass energy of $3.5$ MeV. The lowest energy measured was $365$ keV
for the ground state transition and $222$ keV for the first-excited state one.
These data are presented in Fig. 3 of Ref. [22], which displays statistical
uncertainties only. Subsequently, the Morlock et al. data [22] have been
corrected by Ref. [7] for coincidence summing, and these corrected data have
been adopted for the present analysis. More detailed information regarding the
cross section uncertainties of these data is given in the caption of Fig. 2.37
in Ref. [23], which states “…the statistical errors are with few exceptions
between 1.5% and 3%. One has to add 10% systematic uncertainty (scattering
measurement, energy determinations, error propagation)…” Consequently, in the
present work, we adopt a systematic uncertainty of 10%. The lowest-lying
16O(p,$\gamma$)17F resonance is located near a center-of-mass energy of $2.5$
MeV, and, therefore, we only took their data points below an energy of $2.4$
MeV into account. The data at higher energies are irrelevant for stellar
burning.

We note that the four measurements discussed above provide independent
estimates of the 16O(p,$\gamma$)17F cross section. In particular, only the
work of Ref. [19] relied on stopping power corrections, while the other
measurements of the direct capture cross section [21, 20, 22] were performed
relative to the Rutherford scattering yield, thus obliviating the effects of
target stoichiometry or stopping powers.

The cross sections, $\sigma$, discussed above were converted to astrophysical
$S$-factors, defined by $S(E)$ $\equiv$ $\sigma(E)Ee^{2\pi\eta}$, with $\eta$
denoting the Sommerfeld parameter. The experimental $S$-factors were then
analyzed with our Bayesian model.

##  III Nuclear reaction model

The 16O(p,$\gamma$)17F reaction cross section below $2.4$ MeV center-of-mass
energy is considered as a standard case for the direct radiative capture (DC)
model since the seminal works of Christy and Duck [24] and Rolfs [5].
Subsequent analyses using the direct capture model for 16O(p,$\gamma$)17F can
be found in Refs. [25, 26, 6, 7], and references therein. The study of Ref.
[7] demonstrated that the potential model and the R matrix model provide
nearly identical data fits at low energies. We will adopt in the present work
a single-particle potential model, as discussed below.

The potential model assumes a single-step process, where the proton is
directly captured, without the formation of a compound nucleus, into a final
bound state with the emission of a photon. The dominant E1 contribution to the
theoretical (p,$\gamma$) cross section (in $\mu$b) for capture from an initial
scattering state with orbital angular momentum $\ell_{i}$ to a final bound
state with orbital angular momentum $\ell_{f}$ and the principal quantum
number $n$ (i.e., the number of wave-function nodes), is given by [5]

$\displaystyle\sigma^{DC}_{sp}(E1,n,\ell_{i},\ell_{f})=0.0716\mu^{\frac{3}{2}}\left(\frac{Z_{p}}{M_{p}}-\frac{Z_{t}}{M_{t}}\right)^{2}\frac{E_{\gamma}^{3}}{E^{\frac{3}{2}}}\times$
$\displaystyle\frac{(2J_{f}+1)(2\ell_{i}+1)}{(2j_{p}+1)(2j_{t}+1)(2\ell_{f}+1)}(\ell_{i}010|\ell_{f}0)^{2}R_{n\ell_{i}1\ell_{f}}^{2}$
(1)
$R_{n\ell_{i}1\ell_{f}}=\int_{0}^{\infty}u_{s}(r)\mathcal{O}_{E1}(r)u_{b}(r)dr$
(2)

where $\mu$ is the reduced mass, $Z_{t}$, $Z_{p}$ and $M_{t}$, $M_{p}$ are the
charges and masses (in amu), respectively, of target and projectile; $j_{p}$,
$j_{t}$, $J_{f}$ are the spins of projectile, target and final state,
respectively; $E$ and $E_{\gamma}$ are the center-of-mass energy and the
emitted $\gamma$-ray energy, respectively; $\mathcal{O}_{E1}$ is the radial
part of the E1 multipole operator; and $u_{s}$ and $u_{b}$ are the radial wave
functions of the initial scattering state and final bound state, respectively,
where $u_{b}(r=0)$ $=$ $0$ and $\int_{0}^{\infty}{u_{b}^{2}dr}$ $=$ $1$. We
disregarded any M1 and E2 contributions, which amount to less than $0.1$%
compared to the dominant E1 $S$ factor [5, 27].

Bound-state wave functions were generated by using a potential consisting of a
Woods-Saxon term and a Coulomb term, given by

$V(r)=\frac{-V_{0}}{1+e^{(r-R)/a}}+V_{C}(r)$ (3)

where $R$ $=$ $r_{0}A_{t}^{1/3}$ and $a$ are the Woods-Saxon potential radius
and diffuseness, respectively; $A_{t}$ is the mass number of the target
nucleus; $V_{C}$ corresponds to a uniformly charged sphere of radius $R$. The
well depth, $V_{0}$, was chosen to reproduce the binding energy of the final
state. Several works (see summary in Table I of Ref. [6]) have employed bound-
state square-well potentials instead of Woods-Saxon potentials. However, Refs.
[28, 29] found that the adoption of square-well potentials over-predicts the
calculated single-particle direct capture cross sections by up to a factor of
$3$.

Scattering-state wave functions were computed using a hard-sphere nuclear
potential, which gives similar results as a zero-energy nuclear potential at
the low energies explored here [6]. The insensitivity of the
16O(p,$\gamma$)17F cross section to the choice of scattering potential at low
energies has also been reported in Ref. [30]. The hard-sphere nuclear
potential radius was set equal to the Woods-Saxon radius, $R$, of the bound
state potential. The radial integration in Eq. (2) was extended to $500$ fm,
because the integrand has a maximum located far beyond the nuclear radius at
the lowest center-of-mass energies explored here. For the same reason, we used
the exact expression for the radial part of the E1 operator,
$\mathcal{O}_{E1}$, instead of its long-wavelength approximation.

For a zero-spin target nucleus, such as 16O, the direct capture to a specific
final state proceeds via a unique orbital angular momentum, $\ell_{f}$, but
may involve several values of $\ell_{i}$. In this case, the direct capture
cross section is given by an incoherent sum,

$\sigma^{DC}=C^{2}S_{\ell_{f}}~{}\sum_{\ell_{i}}\sigma^{DC}_{sp}(E1,n,\ell_{i},\ell_{f})$
(4)

with S${}_{\ell_{f}}$ and $C$ denoting the spectroscopic factor and isospin
Clebsch-Gordan coefficient, respectively. For 16O $+$ $p$, we have $C^{2}$ $=$
$1$. To avoid confusion with other symbols, we will not mention further the
isospin Clebsch-Gordan coefficient.

For the 16O(p,$\gamma$)${}^{17}F$ reaction data considered in the present
work, the population of the 17F ground state, $DC$ $\rightarrow$ $0$ keV
($J_{f}$ $=$ $5/2^{+}$), proceeds predominantly via $E1$ radiation and orbital
angular momenta of $\ell_{i}$ $=$ $1$, $3$ and $\ell_{f}$ $=$ $2$, while the
transition to the first excited state, $DC$ $\rightarrow$ $495$ keV ($J_{f}$
$=$ $1/2^{+}$), proceeds via E1 radiation and angular momenta of $\ell_{i}$
$=$ $1$ and $\ell_{f}$ $=$ $0$. We assume that the proton is transferred into
the $1d_{5/2}$ and $2s_{1/2}$ shell-model orbitals for the transition to the
ground and first-excited state, respectively. In the following, we will label
the spectroscopic factor (or later the ANC) for a given final level by the
orbital angular momentum of the bound state alone, and suppress other quantum
numbers (such as total spin). Hence, $S_{gs}$ $\equiv$ $S_{\ell_{f}=2}$ and
$S_{fes}$ $\equiv$ $S_{\ell_{f}=0}$ for the transition to the ground and
first-excited state, respectively.

The calculated single-particle cross section, $\sigma^{DC}_{sp}$, will depend
strongly on the adopted choice of the Woods-Saxon potential radius parameter,
$r_{0}$, and diffuseness, $a$. This means that the spectroscopic factor,
$S_{\ell_{f}}$, which is derived from a fit to experimental data, will also be
sensitive to these very same parameters. We will demonstrate below the
difficulty arising when the data are fitted in terms of the spectroscopic
factor using Eq. (4).

Previous works have demonstrated [5, 22] that, for the 16O(p,$\gamma$)17F
reaction at low energies, the integrand in Eq. (2) peaks far outside the
nuclear radius. For such a peripheral reaction, the single-particle radial
bound state wave function is asymptotically given by [30]

$u_{b,\ell_{f}}(r)\rightarrow b_{\ell_{f}}W_{-\eta,\ell_{f}+1/2}(2\kappa r)$
(5)

where $b_{\ell_{f}}$ is the single-particle asymptotic normalization
coefficient (spANC) and $W$ is the Whittaker function [31]; $\kappa$ is the
bound state wave number, with $\kappa^{2}$ = $2\mu E_{b}/\hbar^{2}$, where
$\mu$ is the reduced mass and $E_{b}$ $=$ $Q-E_{x}$ the binding energy of the
final state; $\eta$ $=$ $eZ_{p}Z_{t}\mu/(\kappa\hbar^{2})$ is the bound state
Coulomb parameter.

In a microscopic nuclear model, the capture cross section can be described in
terms of the overlap, $I_{B}$, of the 16O, 17F, and $p$ bound-state wave
functions and a many-body wave function for the relative motion. Assuming a
single-particle model, the radial dependence of the overlap function, which
represents the projection of the bound final 17F state onto the product bound
state wave functions of 16O and the proton, can be approximated by

$I_{B}(r)\approx\sqrt{S_{\ell_{f}}}u_{b,\ell_{f}}(r)$ (6)

At large distances between target and projectile, the consequences of the
complicated many-body effects will vanish and the radial dependence of the
overlap function becomes asymptotically

$I_{B}(r)\rightarrow C_{\ell_{f}}W_{-\eta,\ell_{f}+1/2}(2\kappa r)$ (7)

where $C_{\ell_{f}}$ is the asymptotic normalization coefficient (ANC).
Comparison of Eqs. (5) $-$ (7) yields the relationship between $S_{\ell_{f}}$,
$C_{\ell_{f}}$, and $b_{\ell_{f}}$ [30]

$S_{\ell_{f}}=\frac{C^{2}_{\ell_{f}}}{b^{2}_{\ell_{f}}}$ (8)

In this approach, $C_{\ell_{f}}^{2}$ appears as an observable quantity, while
both $S_{n\ell_{f}}$ and $b_{\ell_{f}}^{2}$ are derived quantities that depend
strongly on the parameters of the assumed single-particle potential model.

The nuclear model discussed above has three parameters: (i) the ANC,
$C_{\ell_{f}}^{2}$ (or the spectroscopic factor, $S_{\ell_{f}}$) specified by
the final 17F state; (ii) the radius parameter, $r_{0}$, of the Woods-Saxon,
Coulomb, and hard-sphere potentials; and (iii) the diffuseness, $a$. We will
treat these three as adjustable parameters in the fitting. It can be seen from
Eqs. (4) and (8) that $C_{\ell_{f}}^{2}$ or $S_{\ell_{f}}$ act as
multiplicative scaling factors of the computed single-particle potential cross
section (or $S$ factor) for a given transition.

Based on the above discussion, we expect that varying $r_{0}$ and $a$ will
impact the extracted value of $S_{\ell_{f}}$ significantly. When the fit is
performed instead using $C_{\ell_{f}}^{2}$, this dependence should be much
lessened, but may not entirely disappear because of the approximation in Eq.
(6). In either case, we expect that the fit will not constrain $r_{0}$ or $a$.
These effects will be explored in Sec. V.

##  IV Bayesian model

The hierarchical Bayesian model constructed for the present work is similar to
those presented in Refs. [13, 14, 17, 15, 18, 16]. The model takes all
relevant effects impacting the measured data into account. The inference
framework is based on Bayes’ theorem [32]

$p(\theta|y)=\frac{\mathcal{L}(y|\theta)\pi(\theta)}{\int\mathcal{L}(y|\theta)\pi(\theta)d\theta}$
(9)

where the data are represented by the vector $y$ and the complete set of model
parameters is denoted by the vector $\theta$. All factors in Eq. (9) represent
probability densities: $\mathcal{L}(y|\theta)$ is the likelihood, i.e., the
probability that the data, $y$, were obtained assuming given values of the
model parameters, $\theta$; $\pi(\theta)$ is the prior, which represents our
state of knowledge about each parameter before seeing the data; the product of
likelihood and prior defines the posterior, $p(\theta|y)$, i.e., the
probability of obtaining the values of a specific set of model parameters
given the data. Bayes’ theorem implies that the posterior represents an update
of our prior state of knowledge regarding the model parameters once new data
become available. The denominator, called evidence, is a normalization factor
and is not important for the discussion of the present work.

When the experimental $S$ factor is subject to statistical uncertainties only,
$\sigma_{stat}$, the likelihood is given by

$\mathcal{L}(S^{exp}|\theta)=\prod_{i=1}^{N}\frac{1}{\sigma_{stat,i}\sqrt{2\pi}}e^{-\frac{[S_{i}^{exp}-S(\theta)_{i}]^{2}}{2\sigma^{2}_{stat,i}}}$
(10)

where $S(\theta)_{i}$ $=$ $S_{i}^{DC}$ is the theoretical $S$ factor, as
predicted by direct capture theory (Eqs. (4) and (8)), while the product runs
over all the data points, labeled by the index $i$. The likelihood represents
a product of normal distributions, each with a mean of $S(\theta)_{i}$ and a
standard deviation of $\sigma_{stat,i}$, with the latter quantity given by the
experimental statistical uncertainty of datum $i$. In symbolic notation, we
can express Eq. (10) as

$S_{i}^{exp}\sim N(S(\theta)_{i},\sigma^{2}_{stat,i})$ (11)

where “$N$” denotes a normal (Gaussian) probability density and the symbol
“$\sim$” stands for “has the probability distribution of.”

According to Eqs. (4) and (8), the direct capture $S$ factor, $S^{DC}$, for a
given transition (i.e., to the ground or first-excited state) at the center-
of-mass energy, $E$, has the form

$S^{DC}(E)=S_{\ell_{f}}S_{sp}^{DC}(E,r_{0},a)=C^{2}_{\ell_{f}}\frac{S_{sp}^{DC}(E,r_{0},a)}{b^{2}_{\ell_{f}}(r_{0},a)}$
(12)

where $S_{sp}^{DC}$ denotes the single-particle direct capture $S$ factor of
Eq. (4). When the data analysis is performed using the ANC as scaling factor
to fit the data, the ratio $S_{sp}^{DC}/b^{2}_{\ell_{f}}$ is computed by the
nuclear model. Alternatively, when the fit is scaled using the spectroscopic
factor, the nuclear model calculates the quantity $S_{sp}^{DC}$. Both
$S_{sp}^{DC}$ and $b^{2}_{\ell_{f}}$ depend explicitly on the parameters
$r_{0}$ and $a$ of the Woods-Saxon potential. The total 16O(p,$\gamma$)17F $S$
factor is given by the incoherent sum over the individual transitions. The
results quoted in Sec. V will be based on adopting $C^{2}_{\ell_{f}}$ as the
scaling factor. We will briefly discuss the outcomes of tests performed when
adopting instead the spectroscopic factor.

Each model parameter requires a prior. Since one goal of the present work was
to explore the impact of Woods-Saxon potential parameter variations over wide
ranges, we adopted weakly-informative priors for the radius parameter,
$r_{0}$, and diffuseness, $a$,

$\displaystyle r_{0}\sim U(0.9~{}\mathrm{fm},1.5~{}\mathrm{fm})$ (13)
$\displaystyle a\sim U(0.5~{}\mathrm{fm},0.7~{}\mathrm{fm})$ (14)

where “$U$” denotes a uniform (i.e., constant) probability density between the
stated boundaries. Outside the ranges given above, the numerical integration
does not have a solution with the required number of nodes in the bound state
wave function, $u_{b}(r)$, when the depth of the Woods-Saxon potential is
adjusted to reproduce the binding energy of the experimental level.

Previous work [33, 26, 7, 34] has established ANC values near $C^{2}_{gs}$
$\approx$ $1$ and $C^{2}_{fes}$ $\approx$ $7000$ for the transitions to the
ground and first-excited states, respectively (see Sec. V). Our goal is to
predict ANCs based on directly measured 16O(p,$\gamma$)17F data alone. In
particular, we do not want to restrict the sampling ranges of the ANCs based
on prior information from either transfer measurements or theoretical model
calculations. For this reason, we will adopt non-informative priors

$\displaystyle C^{2}_{gs}\sim T(0,\infty)N(0,[10^{6}~{}\mathrm{fm^{-1}}]^{2})$
(15) $\displaystyle C^{2}_{fes}\sim
T(0,\infty)N(0,[10^{6}~{}\mathrm{fm^{-1}}]^{2})$ (16)

where “$N$” denotes a truncated normal density centered at zero with a large
standard deviation, multiplied by a truncation function, “$T$”, which
suppresses samples outside of the interval defined by its arguments.

Following Ref. [13], we will incorporate into our model the systematic
uncertainty in a given experiment as an informative, lognormal prior, given by

$\pi(f)=\frac{1}{\ln(f.u.)\sqrt{2\pi}f}e^{-\frac{[\ln
f]^{2}}{2[\ln(f.u.)]^{2}}}$ (17)

or in symbolic notation

$f\sim LN(0,[\ln(f.u.)]^{2})$ (18)

where “$LN$” denotes a lognormal probability density. In Eqs. (17) and (18),
the lognormal location parameter, $\mu$, is equal to zero, i.e., $\mu$ $=$
$\ln x_{med}$ $=$ $0$. In other words, we are assuming that the median value
of the lognormal distribution, $x_{med}$, equals unity. If this would not be
the case, we could have corrected the data for the systematic effect. The
lognormal spread parameter, $\sigma$, in Eqs. (17) and (18) is chosen as
$\sigma$ $=$ $\ln(f.u.)$, where $f.u.$ denotes the factor uncertainty, which
is given by the systematic uncertainty. For example, we pointed out in Sec. II
that the data of Hester, Pixley and Lamb [19] have a systematic uncertainty of
14%. In this case, we have chosen $\sigma$ $=$ $\ln(f.u.)$ $=$ $\ln(1.14)$.
The posterior of the “normalization factor,” $f$, will describe by how much a
specific data set deviates from the best-fit line.

In conventional $\chi^{2}$ fitting, normalization factors are viewed as a
systematic shift in the data. In Bayesian inference, on the other hand, the
reported data are never modified. Instead, the true (but unknown) S-factor is
multiplied by the normalization factor, $f$. This means that, during the
fitting, each data set pulls on the true S-factor curve with a strength
inversely proportional to the systematic uncertainty: a data set with a small
systematic uncertainty will pull the true S-factor curve more strongly toward
it compared to a set with a large systematic uncertainty. This “pulling” is
independent of the data set size: disregarding statistical uncertainties for a
moment, a set consisting of a single datum will have the same weight in the
fitting as one containing many data points, if both sets are described by the
same factor uncertainty, $f.u.$

In previous work (e.g., Ref. [17]), the effect of “extrinsic” scatter was
explicitly included in the Bayesian model to account for an additional and
unreported source of statistical uncertainty. For the case of
16O(p,$\gamma$)17F, we find no compelling reason to include this effect.
Neither the data points of Hester, Pixley and Lamb [19], nor those of Chow,
Griffiths and Hall [21] exhibit a scatter inconsistent with their reported
uncertainties. The data of Morlock et al. [22], both for the transition to the
ground and first-excited state, indeed reveal a scatter that clearly exceeds
their reported statistical and systematic uncertainties. But this affects only
the data points at their lowest measured energies, which also appear to be
systematically too high (see below). Therefore, it would be inappropriate to
include this effect as “extrinsic scatter” in the Bayesian model. To test
whether the observed low-energy scatter in the data of Ref. [22] influences
any of our results, we implemented a robust algorithm into the Bayesian model
that accounts for outliers [35, 13]. It was found that the observed low-energy
scatter does not impact the fit.

Our complete Bayesian model for a given transition, either to the ground or
first-excited state, can be summarized below in symbolic notation:

Nuclear reaction model: $\displaystyle\indent
S(\theta)_{i}=C^{2}_{\ell_{f}}\frac{S_{sp}^{DC}(E,r_{0},a)}{b^{2}_{\ell_{f}}(r_{0},a)}$
(19) Parameters: $\displaystyle\indent\theta\equiv(C^{2}_{\ell_{f}},r_{0},a)$
(20) Likelihoods: $\displaystyle\indent S^{\prime}_{i,j}=f_{j}\times
S(\theta)_{i}$ (21) $\displaystyle\indent S^{exp}_{i,j}\sim
N(S^{\prime}_{i,j},\sigma_{stat,i}^{2})$ (22) Priors: $\displaystyle\indent
C^{2}_{\ell_{f}}\sim T(0,\infty)N(0,[10^{6}~{}\mathrm{fm^{-1}}]^{2})$ (23)
$\displaystyle\indent r_{0}\sim U(0.9~{}\mathrm{fm},1.5~{}\mathrm{fm})$ (24)
$\displaystyle\indent a\sim U(0.5~{}\mathrm{fm},0.7~{}\mathrm{fm})$ (25)
$\displaystyle\indent f_{j}\sim LN(0,[\ln(f.u.)_{j}]^{2})$ (26)

The index $i$ labels individual data points, and $j$ denotes the experiment.

As mentioned earlier, Hester, Pixley and Lamb [19] only measured the total $S$
factor. These data points were simply described in our Bayesian model by the
sum of the ground and first-excited-state transitions. Therefore, they were
analyzed simultaneously with all other data sets and thereby constrain the
strengths of the transitions to the individual levels. At the same time, the
fitting correlations between the two partial $S$ factors are fully contained
in our results. We note that in the fitting we did not take into account the
$\pm 2$ keV beam energy uncertainties reported in Ref. [19]. Their effect on
the fit will be much smaller than the uncertainties of the reported cross
sections, as can be seen from their Fig. 4.

For the analysis of the Bayesian model, we employed the program JAGS (“Just
Another Gibbs Sampler”) using Markov chain Monte Carlo (MCMC) sampling [36].
Specifically, we will employ the rjags package that works directly with JAGS
within the R language [37]. Running a JAGS model refers to generating random
samples from the posterior distribution of model parameters. This involves the
definition of the model, likelihood, and priors, as well as the
initialization, adaptation, and monitoring of the Markov chains.

The MCMC sampling will provide the posteriors of all parameters. For each
simulation, we computed three MCMC chains, where each had a length of $3\times
10^{6}$ steps after disregarding the burn-in samples ($10^{5}$ steps for each
chain). This ensured that the chains reached equilibrium, that effective chain
lengths for all parameters exceeded $10000$, and that Monte Carlo fluctuations
were negligible compared to statistical and systematic uncertainties.

##  V Results

The resulting $S$-factor fits, adopting the model presented in Eqs. (IV) $-$
(26), are presented in Fig. 1. The top panel corresponds to the sum of the
ground and first-excited-state transitions, while the middle and bottom ones
depict the fits to the individual transitions 111The negative slope of the $S$
factor for the transition to the first-excited state has been variously
explained by a pole in the $S$ factor [38] or a supposed halo property of the
first-excited state in 17F [22]. The correct explanation was given by Ref.
[39], who demonstrated that, for weakly-bound final states, the upturn of the
$S$ factor towards lower energies arises naturally from the interplay of the
Whittaker function corresponding to the bound state and the regular Coulomb
wave function associated with the scattering state.. As pointed out in Sec.
II, only Chow, Griffiths and Hall [21] (red full circles) and Morlock et al.
[22] (black open circles) measured cross sections for both transitions
separately. The single data point of Becker et al. [20] (green full circle)
was measured for the first-excited-state transition and, thus, appears only in
the bottom panel. Hester, Pixley and Lamb [19] (blue open circles) measured
the total cross section and their data appear in the top panel only. The dark
and light shaded areas represent our $S$ factor predictions for 68% and 95%
coverage probabilities, respectively.

Numerical results are listed in Table 1. For the total (i.e., the sum of
ground and first-excited state transition) zero-energy S factor, we find a
value of $S(0)$ $=$ $0.01092$ MeVb ($\pm$4.0%), where the quoted uncertainties
are derived from the 16, 50, and 84 percentiles. Our result agrees with that
of Gagliardi et al. [26], although their value has a significantly larger
uncertainty (10%). The latter work determined the ANCs of the ground and
first-excited state by measuring the 16O(3He,d)17F transfer reaction as a test
case for the indirect estimation of the low-energy 16O(p,$\gamma$)17F S
factor. Our best estimate for the total S factor at $E$ $=$ $0.090$ MeV is
$S(0.090)$ $=$ $0.00786$ MeVb ($\pm$4.0%). This agrees with the average value
found by Iliadis et al. [7] for R-matrix and potential-model fits to the data
of Refs. [21, 22], but their uncertainty (7.5%) was not determined rigorously
in a statistical sense, as discussed in Sec. I. Our value also agrees with the
result obtained in the R-matrix analysis of the Morlock et al. [22] data by
Azuma et al. [34], but no uncertainties were reported in the latter work. For
future reference, we provide in Table 2 numerical values of our recommended
$S$ factor for a grid of energies.

Figure 1:  Results of Bayesian fits to experimental 16O(p,$\gamma$)17F $S$ factors [19, 21, 20, 22]: (Top) Total $S$ factor; (Middle) $S$ factor for the transition to the 17F ground state; (Bottom) $S$ factor for the transition to the 17F first-excited state. Becker et al. [20] measured only the $S$ factor for the first-excited state transition at a single energy (green datum in the bottom panel). Hester, Pixley and Lamb [19] measured only the total $S$ factor (blue data points in the top panel). The dark and light shaded regions represent coverage probabilities of 68% and 95%, respectively. Notice the logarithmic scale on the abscissa.  Table 1: Summary of ANCs and $S$ factors for 16O(p,$\gamma$)17F. | Present | Previous
---|---|---
| 111Results of fit to all 16O(p,$\gamma$)17F data adopting broad priors. See text. | 222Results of fit to all 16O(p,$\gamma$)17F data for fixed values of $r_{0}$ $=$ $1.25$ fm and $a$ $=$ $0.65$ fm. See text.  | Baye333Values calculated using a microscopic nuclear model.  [33] | Gagliardi444Results of fit to 16O(3He,d)17F data.  [26] | Iliadis [7] | Azuma777Results of R-matrix fit to the 16O(p,$\gamma$)17F data of Ref. [22]. No uncertainties were reported. [34]
$C^{2}_{gs}$ (fm-1) | 1.115 ($\pm$4.0%) |  $1.006$ ($\pm$4.0%) |  $1.19$, $0.947$ |  $1.08$ ($\pm$9.3%) |  $1.34$555Results of R-matrix fits to the 16O(p,$\gamma$)17F data of Refs. [21, 22]. The small uncertainties quoted in that work do not include any systematic uncertainties of the data. | $1.10$
$C^{2}_{fes}$ (fm-1) | 7063 ($\pm$4.0%) |  $6759$ ($\pm$4.0%) |  $8306$, $7468$ |  $6490$ ($\pm$10%) |  $6667$555Results of R-matrix fits to the 16O(p,$\gamma$)17F data of Refs. [21, 22]. The small uncertainties quoted in that work do not include any systematic uncertainties of the data. | $6512$
$S_{gs+efs}(0~{}\mathrm{MeV})$ (MeVb) | 0.01092 ($\pm$4.0%) |  $0.01040$ ($\pm$4.0%) |  | 0.0102 ($\pm$10%) |  |
$S_{gs+efs}(0.09~{}\mathrm{MeV})$ (MeVb) | 0.00786 ($\pm$4.0%) |  $0.00751$ ($\pm$4.0%) |  |  |  $0.0076$ ($\pm$7.2%)666Average result of R-matrix and potential model fits to the 16O(p,$\gamma$)17F data of Refs. [21, 22]. | $0.00807$
Table 2: Predicted total $S$ factors for 16O(p,$\gamma$)17F. The uncertainty of the $S$ factor is $4.0$% at all listed energies. E | S(E) | E | S(E)
---|---|---|---
(MeV) | (MeVb) | (MeV) | (MeVb)
0.002 | 0.01081 | 0.100 | 0.00767
0.004 | 0.01069 | 0.200 | 0.00639
0.006 | 0.01059 | 0.400 | 0.00521
0.008 | 0.01048 | 0.600 | 0.00465
0.010 | 0.01038 | 0.800 | 0.00435
0.020 | 0.00992 | 1.000 | 0.00413
0.040 | 0.00915 | 1.500 | 0.00376
0.060 | 0.00856 | 2.000 | 0.00356
0.080 | 0.00807 | 2.500 | 0.00342

Posteriors of the two ANCs, $C^{2}_{gs}$ and $C^{2}_{fes}$, including their
pairwise correlation, are presented in Fig. 2. As expected, the two quantities
are highly correlated: when one is increased, the other one decreases to fit
the total experimental $S$ factor. Our best-fit values are $C^{2}_{gs}$ $=$
$1.115$ fm-1 ($\pm$ 4.0%) and $C^{2}_{fes}$ $=$ $7063$ fm-1 ($\pm$ 4.0%). The
results agree with those of Gagliardi et al. [26]. It is not surprising that
their uncertainties ($\approx 10$%) are significantly larger than ours,
because their distorted wave Born approximation (DWBA) analysis of
16O(3He,d)17F data was subject to ambiguities in the choice of optical model
potentials for the incoming and outgoing channels. Our values are also in the
vicinity of previously reported theoretical [33] and experimental results [7,
34]. A more detailed comparison is not possible, as no rigorous uncertainties
were presented in these works. We do not present posteriors of the Woods-Saxon
radius parameter, $r_{0}$, or the diffuseness, $a$, because the fit does not
constrain these parameters, as discussed in Sec. III.

Figure 2:  One- and two-dimensional projections of the posterior probability
distributions of two physical model parameters, i.e., the ANCs, $C^{2}_{gs}$
and $C^{2}_{fes}$, and four data set parameters, $f_{j}$, when adopting broad
priors in the Bayesian fit. These results correspond to those given in column
2 of Table 1 (boldface) and in Table 3). The dark and light shaded areas
correspond to 68% and 95% coverage probabilities, respectively.

Posteriors of the data-set normalization factors, $f$, are also depicted in
Fig. 2. Numerical values for both $f.u.$ and $f$ are given in Table 3. The
former determine the spreads of the lognormal priors used to implement
systematic uncertainties into the Bayesian model, while the latter are derived
from the posteriors and represent our best estimate for the deviations between
the data of a given experiment and the best-fit line. It can be seen that the
$f$ values for the data of Refs. [21, 20, 22] are consistent with unity for a
coverage probability of 64%. For the data of Ref. [19], the $f$ value is
consistent with unity for a coverage probability of 80%. Overall, our Bayesian
model predicts that the reported systematic uncertainties are consistent with
the best-fit solution, i.e., we find no evidence for significant unreported
systematic effects.

Table 3: Normalization factors, $f$, of 16O(p,$\gamma$)17F data sets. Experiment |  $f.u.$111Systematic factor uncertainties of the experimental data (see Sec. II). The values of $f.u.$ determine the spread of the lognormal probability density adopted for the priors. For example, a value of “1.14” represents a systematic uncertainty of “14%”. |  $f$222Results correspond to those given in boldface in column 2 of Table 1, i.e., by assuming broad priors. The uncertainties correspond to 68% coverage probabilities of the posterior densities.
---|---|---
Hester, Pixley & Lamb [19] | $1.14$ | $0.924\pm 0.070$
Chow, Griffiths & Hall [21] | $1.05$ | $1.022\pm 0.041$
Becker et al. [20] | $1.05$ | $0.986\pm 0.045$
Morlock et al. [22] | $1.10$ | $0.994\pm 0.041$

The small uncertainty of $4.0$% in our results (column 2 of Table 1), despite
the fact that we varied the Woods-Saxon potential parameters over broad
ranges, reflects the insensitivity of the ANCs, $C^{2}_{\ell_{f}}$, and fitted
$S$ factor to the choice of potential parameters, as discussed in Sec. III and
previously demonstrated by Gagliardi et al. [26] in the analysis of
16O(3He,d)17F data. As an additional test of this assertion, we repeated the
calculation by fixing the bound state potential parameters at values of
$r_{0}$ $=$ $1.25$ fm and $a$ $=$ $0.65$ fm. The results are listed in column
3 of Table 1. It can be seen that the values and their uncertainties obtained
for the ANCs and S factors overall agree with those obtained by sampling over
broad ranges of the potential parameters.

We will now turn to the discussion of spectroscopic factors. If we assume the
conventional Woods-Saxon potential parameter values of $r_{0}$ $=$ $1.25$ fm
and $a$ $=$ $0.65$ fm, we find $b^{2}_{gs}$ $=$ $0.8840$ fm-1 and
$b^{2}_{fes}$ $=$ $6899$ fm-1 for the single-particle ANCs. From Eq. (8) and
the results listed in column 2 of Table 1, we obtain for the spectroscopic
factors $S_{\ell_{f}=2}$ $=$ $1.26\pm 0.05$ and $S_{\ell_{f}=0}$ $=$ $1.02\pm
0.04$ for the transitions to the ground and first-excited state, respectively.
These agree with previous theoretical and experimental analyses (see, e.g.,
Refs. [40, 6, 7]).

When the fit is performed using spectroscopic factors instead of ANCs as
scaling factors (see Eq. (12), the obtained uncertainties become significantly
larger. Assuming again broad priors for $r_{0}$ and $a$ (see Eqs. (13) and
(14)), the zero-energy $S$-factor uncertainty amounts to 9% (i.e., more than
twice the value discussed above). Similarly, the uncertainty of the derived
spectroscopic factors increases substantially (to 11%). These uncertainties
can only be reduced if the sampling ranges of the priors for the potential
parameters are significantly reduced. Generally, however, these parameters are
poorly constrained by experiment and the commonly employed values of $r_{0}$
$=$ $1.25$ fm and $a$ $=$ $0.65$ fm are mainly adopted by convention. Clearly,
performing the fits in terms of ANCs has a distinct advantage over using
spectroscopic factors as scaling parameters.

##  VI Thermonuclear reaction rates

The thermonuclear reaction rate per particle pair, $N_{A}\langle\sigma
v\rangle$, at a given stellar temperature, $T$, is given by [3]

$\begin{split}N_{A}\langle\sigma
v\rangle&=\left(\frac{8}{\pi\mu}\right)^{1/2}\frac{N_{A}}{(kT)^{3/2}}\\\
&\times\int_{0}^{\infty}e^{-2\pi\eta}\,S(E)\,e^{-E/kT}\,dE,\end{split}$ (27)

where $\mu$ is the reduced mass of projectile and target, $N_{A}$ is the
Avogadro constant, $k$ is the Boltzmann constant, and $E$ is the 16O $+$ $p$
center-of-mass energy.

The 16O(p,$\gamma$)17F reaction rates were computed by integrating Equation
(27) numerically. The $S$ factor is calculated from our MCMC samples (Sec. V)
and, therefore, the values of $N_{A}\langle\sigma v\rangle$ contain the
effects of statistical and systematic uncertainties, as well as all
correlations among parameters. The results are based on 5,000 random S-factor
samples, which ensures that Markov chain Monte Carlo fluctuations are
negligible compared to the reaction rate uncertainties. Our lower and upper
integration limits were set to $0.002$ MeV and $2.5$ MeV, respectively.
Reaction rates are computed for $50$ temperature grid points between $0.003$
GK and $3.5$ GK and numerical values presented in Table 4. The recommended
rates are computed from the 50th percentile of the rate probability density
function, while the factor uncertainty, $f.u.$, is obtained from the 16th and
84th percentiles [41]. The new total rate uncertainties are 4.2% for the
entire temperature range shown. Notice that this uncertainty is slightly
larger than that of the $S$ factor given in Table 2 (4.0%). As already pointed
out, the reason is that we computed the rates directly using the MCMC samples
and, therefore, our results contain the information of parameter correlations.

Table 4: Total thermonuclear reaction rates for 16O(p,$\gamma$)17F a

T (GK) | Low | Median | High | f.u.
---|---|---|---|---
0.003 | 4.164E-41 | 4.338E-41 | 4.517E-41 | 1.042
0.004 | 1.343E-36 | 1.399E-36 | 1.457E-36 | 1.042
0.005 | 2.154E-33 | 2.244E-33 | 2.336E-33 | 1.042
0.006 | 5.980E-31 | 6.230E-31 | 6.486E-31 | 1.042
0.007 | 5.332E-29 | 5.555E-29 | 5.783E-29 | 1.042
0.008 | 2.162E-27 | 2.253E-27 | 2.345E-27 | 1.042
0.009 | 4.935E-26 | 5.142E-26 | 5.353E-26 | 1.042
0.010 | 7.290E-25 | 7.595E-25 | 7.907E-25 | 1.042
0.011 | 7.668E-24 | 7.988E-24 | 8.317E-24 | 1.042
0.012 | 6.148E-23 | 6.405E-23 | 6.669E-23 | 1.042
0.013 | 3.951E-22 | 4.117E-22 | 4.286E-22 | 1.042
0.014 | 2.114E-21 | 2.203E-21 | 2.293E-21 | 1.042
0.015 | 9.699E-21 | 1.010E-20 | 1.052E-20 | 1.042
0.016 | 3.903E-20 | 4.067E-20 | 4.234E-20 | 1.042
0.018 | 4.578E-19 | 4.770E-19 | 4.966E-19 | 1.042
0.020 | 3.809E-18 | 3.968E-18 | 4.131E-18 | 1.042
0.025 | 2.639E-16 | 2.749E-16 | 2.862E-16 | 1.042
0.030 | 6.641E-15 | 6.919E-15 | 7.204E-15 | 1.042
0.040 | 7.211E-13 | 7.513E-13 | 7.822E-13 | 1.042
0.050 | 1.998E-11 | 2.082E-11 | 2.167E-11 | 1.042
0.060 | 2.497E-10 | 2.601E-10 | 2.708E-10 | 1.042
0.070 | 1.865E-09 | 1.943E-09 | 2.023E-09 | 1.042
0.080 | 9.751E-09 | 1.016E-08 | 1.058E-08 | 1.042
0.090 | 3.933E-08 | 4.098E-08 | 4.267E-08 | 1.042
0.100 | 1.304E-07 | 1.359E-07 | 1.414E-07 | 1.042
0.110 | 3.709E-07 | 3.865E-07 | 4.024E-07 | 1.042
0.120 | 9.340E-07 | 9.731E-07 | 1.013E-06 | 1.042
0.130 | 2.129E-06 | 2.219E-06 | 2.310E-06 | 1.042
0.140 | 4.472E-06 | 4.659E-06 | 4.851E-06 | 1.042
0.150 | 8.764E-06 | 9.131E-06 | 9.507E-06 | 1.042
0.160 | 1.620E-05 | 1.688E-05 | 1.757E-05 | 1.042
0.180 | 4.791E-05 | 4.991E-05 | 5.197E-05 | 1.042
0.200 | 1.216E-04 | 1.266E-04 | 1.319E-04 | 1.042
0.250 | 7.781E-04 | 8.106E-04 | 8.440E-04 | 1.042
0.300 | 3.176E-03 | 3.309E-03 | 3.446E-03 | 1.042
0.350 | 9.703E-03 | 1.011E-02 | 1.053E-02 | 1.042
0.400 | 2.426E-02 | 2.528E-02 | 2.632E-02 | 1.042
0.450 | 5.244E-02 | 5.464E-02 | 5.689E-02 | 1.042
0.500 | 1.016E-01 | 1.058E-01 | 1.102E-01 | 1.042
0.600 | 3.000E-01 | 3.126E-01 | 3.255E-01 | 1.042
0.700 | 7.083E-01 | 7.381E-01 | 7.683E-01 | 1.042
0.800 | 1.432E+00 | 1.493E+00 | 1.554E+00 | 1.042
0.900 | 2.590E+00 | 2.699E+00 | 2.809E+00 | 1.042
1.000 | 4.301E+00 | 4.482E+00 | 4.666E+00 | 1.042
1.250 | 1.180E+01 | 1.229E+01 | 1.280E+01 | 1.042
1.500 | 2.526E+01 | 2.632E+01 | 2.740E+01 | 1.042
1.750 | 4.612E+01 | 4.806E+01 | 5.005E+01 | 1.042
2.000 | 7.545E+01 | 7.863E+01 | 8.188E+01 | 1.042
2.500 | 1.618E+02 | 1.686E+02 | 1.755E+02 | 1.042
3.000 | 2.851E+02 | 2.971E+02 | 3.093E+02 | 1.042
3.500 | 4.406E+02 | 4.591E+02 | 4.779E+02 | 1.042

  * a

In units of cm3mol-1s-1. Columns 2, 3, and 4 list the 16th, 50th, and 84th
percentiles, respectively, of the total rate probability density function at
the given temperatures; $f.u.$ is the factor uncertainty of the total reaction
rate, based on Monte Carlo sampling. The total number of samples at each
temperature was 5,000.

The new rates are displayed in Fig. 3 and compared to previous results. For
better comparison, all rates have been divided by the present median (50th
percentile) values (column 3 of Table 4). Our low (16th percentile) and high
(84th percentile) rates are depicted as magenta-shaded bands centered around
unity (dashed line). The gray bands refer to the rates of Iliadis  _et al._
[7] (top panel) and Angulo  _et al._ [12] (bottom panel). It can be seen that
our median rate is about 5% and 8% higher than those of Iliadis  _et al._ [7]
and Angulo  _et al._ [12], respectively. Furthermore, the uncertainties of our
new rate ($4.2$%) are about half of those reported by Ref. [7] ($7.5$%) and
about a factor of $8$ smaller than those of Ref. [12] ($30$%). As already
pointed out above, all previously evaluated rates lack a rigorous statistical
interpretation based on coverage probabilities.

Figure 3: Comparison of present to previously evaluated 16O(p,$\gamma$)17F
reaction rates. For better comparison, all rates have been normalized to the
present median rate (see column 3 in Table 4). The magenta-shaded regions
correspond to our 68% uncertainties (see column 5 in Table 4). The gray-shaded
regions depict the rates and uncertainties reported in Ref. [7] (top panel)
and Ref. [12] (bottom panel).

##  VII Summary

We presented the first statistically rigorous evaluation of the
16O(p,$\gamma$)17F direct capture $S$ factor and thermonuclear reaction rates.
Our analysis includes four data sets that reported statistical and systematic
uncertainties separately [21, 20, 22]. The combined fit for the transitions to
the ground and first-excited states in 17F, and their sum, was performed using
a Bayesian method. The physical model was a single-particle model employing a
Woods-Saxon potential for generating the radial bound state wave function. The
fit had three adjustable parameters: the radius parameter and diffuseness of
the Woods-Saxon potential, and the asymptotic normalization coeffients (ANCs)
for scaling the theoretical direct capture $S$ factor. We found that a poor
fit is obtained when it is performed using spectroscopic factors as scaling
parameters instead. Since the 16O(p,$\gamma$)17F reaction at low bombarding
energies is of peripheral nature, the analysis of the $S$ factor in terms of
ANCs greatly reduces the sensitivity to the single-particle potential
parameters. For the ANC’s of the ground and first-excited state transitions,
we find values of $C^{2}_{gs}$ $=$ $1.115$ fm-1 ($\pm 4.0$%) and $C^{2}_{fes}$
$=$ $7063$ fm-1 ($\pm 4.0$%), respectively. The total $S$ factor at zero
energy is $S_{gs+fes}(0)$ $=$ $0.01092$ MeVb ($\pm 4.0$%). The thermonuclear
rate is computed from the Markov Chain Monte Carlo (MCMC) samples by numerical
integration and, therefore, accounts for all correlations between parameters.
The rate uncertainties between $3$ MK and $3.5$ GK is $4.2$%, about half the
previously reported rates [7].

Our recommendation for future work is to remeasure the 16O(p,$\gamma$)17F
reaction cross section at center-of-mass energies below $200$ keV, which have
so far only been reached in the experiment of Hester, Pixley and Lamb [19],
albeit with a relatively large systematic uncertainty of $14$%. Finally, we
emphasize that the Bayesian fit results presented here are based on a specific
nuclear reaction (i.e., a single-particle potential) model. It would be
interesting to compare these results to those from a future R-matrix analysis
of the same data as analyzed here.

###### Acknowledgements.

We would like to thank Art Champagne and Richard Longland for helping with
calculating the Whittaker function. We also would like to express our
gratitude Jo Moscoso for his help during the early stages of this project. The
constructive comments of Robert Janssens are highly appreciated. This work was
supported in part by the DOE, Office of Science, Office of Nuclear Physics,
under grants DE-FG02-97ER41041 (UNC) and DE-FG02-97ER41033 (TUNL).

## References

  * Wang  _et al._ [2021] M. Wang, W. Huang, F. Kondev, G. Audi, and S. Naimi, Chinese Physics C 45, 030003 (2021).
  * Iliadis  _et al._ [2010] C. Iliadis, R. Longland, A. Champagne, A. Coc, and R. Fitzgerald, Nuclear Physics A 841, 31 (2010).
  * Iliadis [2015] C. Iliadis,  _Nuclear Physics of Stars_ , 2nd ed. (Wiley-VCH, Weinheim, Germany, 2015).
  * Tilley  _et al._ [1993] D. R. Tilley, H. R. Weller, and C. M. Cheves, Nuclear Physics A 564, 1 (1993).
  * Rolfs [1973] C. Rolfs, Nuclear Physics A 217, 29 (1973).
  * Iliadis and Wiescher [2004] C. Iliadis and M. Wiescher, Phys. Rev. C 69, 064305 (2004).
  * Iliadis  _et al._ [2008] C. Iliadis, C. Angulo, P. Descouvemont, M. Lugaro, and P. Mohr, Phys. Rev. C 77, 045802 (2008).
  * Mohr and Iliadis [2012] P. Mohr and C. Iliadis, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 688, 62 (2012).
  * Stoesz and Herwig [2003] J. A. Stoesz and F. Herwig, Mon. Not. Roy. Astr. Soc. 340, 763 (2003), arXiv:astro-ph/0212128 [astro-ph] .
  * Lugaro, M.  _et al._ [2007] Lugaro, M., Karakas, A. I., Nittler, L. R., O’D. Alexander, C. M., Hoppe, P., Iliadis, C., and Lattanzio, J. C., A&A 461, 657 (2007).
  * Palmerini  _et al._ [2021] S. Palmerini, S. Cristallo, M. Busso, M. La Cognata, M. L. Sergi, and D. Vescovi, Frontiers in Astronomy and Space Sciences 7 (2021), 10.3389/fspas.2020.607245.
  * Angulo  _et al._ [1999] C. Angulo, M. Arnould, M. Rayet, P. Descouvemont, D. Baye, C. Leclercq-Willain, A. Coc, S. Barhoumi, P. Aguer, C. Rolfs, R. Kunz, J. W. Hammer, A. Mayer, T. Paradellis, S. Kossionides, C. Chronidou, K. Spyrou, S. degl’Innocenti, G. Fiorentini, B. Ricci, S. Zavatarelli, C. Providencia, H. Wolters, J. Soares, C. Grama, J. Rahighi, A. Shotter, and M. Lamehi Rachti, Nuclear Physics A 656, 3 (1999).
  * Iliadis  _et al._ [2016] C. Iliadis, K. S. Anderson, A. Coc, F. X. Timmes, and S. Starrfield, Astrophys. J. 831, 107 (2016).
  * Gómez Iñesta  _et al._ [2017] A. Gómez Iñesta, C. Iliadis, and A. Coc, Astrophys. J. 849, 134 (2017).
  * de Souza  _et al._ [2019a] R. S. de Souza, C. Iliadis, and A. Coc, The Astrophysical Journal 872, 75 (2019a).
  * Moscoso  _et al._ [2021] J. Moscoso, R. S. de Souza, A. Coc, and C. Iliadis, The Astrophysical Journal 923, 49 (2021).
  * de Souza  _et al._ [2019b] R. S. de Souza, S. R. Boston, A. Coc, and C. Iliadis, The Physical Review C 99, 014619 (2019b).
  * de Souza  _et al._ [2020] R. S. de Souza, T. H. Kiat, A. Coc, and C. Iliadis, The Astrophysical Journal 894, 134 (2020).
  * Hester  _et al._ [1958] R. E. Hester, R. E. Pixley, and W. A. S. Lamb, Phys. Rev. 111, 1604 (1958).
  * Becker  _et al._ [1982] H.-W. Becker, W. E. Kieser, H. P. Trautvetter, and W. Wiescher, Zeitschrift für Physik A Atoms and Nuclei 305, 319 (1982).
  * Chow  _et al._ [1975] H. C. Chow, G. Griffiths, and T. Hall, Canadian Journal of Physics 53, 1672 (1975).
  * Morlock  _et al._ [1997] R. Morlock, R. Kunz, A. Mayer, M. Jaeger, A. Müller, J. W. Hammer, P. Mohr, H. Oberhummer, G. Staudt, and V. Kölle, Phys. Rev. Lett. 79, 3837 (1997).
  * Morlock [1997] R. Morlock, Ph.D. thesis, Universität Stuttgart (1997).
  * Christy and Duck [1961] R. Christy and I. Duck, Nuclear Physics 24, 89 (1961).
  * Brune [1996] C. Brune, Nuclear Physics A 596, 122 (1996).
  * Gagliardi  _et al._ [1999] C. A. Gagliardi, R. E. Tribble, A. Azhari, H. L. Clark, Y.-W. Lui, A. M. Mukhamedzhanov, A. Sattarov, L. Trache, V. Burjan, J. Cejpek, V. Kroha, i. c. v. Piskoř, and J. Vincour, Phys. Rev. C 59, 1149 (1999).
  * Descouvemont and Baye [1999] P. Descouvemont and D. Baye, Phys. Rev. C 60, 015803 (1999).
  * Wiescher  _et al._ [1980] M. Wiescher, H. Becker, J. Görres, K.-U. Kettner, H. Trautvetter, W. Kieser, C. Rolfs, R. Azuma, K. Jackson, and J. Hammer, Nuclear Physics A 349, 165 (1980).
  * Powell  _et al._ [1999] D. Powell, C. Iliadis, A. Champagne, C. Grossmann, S. Hale, V. Hansper, and L. McLean, Nuclear Physics A 660, 349 (1999).
  * Mukhamedzhanov  _et al._ [2001] A. M. Mukhamedzhanov, C. A. Gagliardi, and R. E. Tribble, Phys. Rev. C 63, 024612 (2001).
  * Hebbard and Robson [1963] D. Hebbard and B. Robson, Nuclear Physics 42, 563 (1963).
  * Hilbe  _et al._ [2017] J. M. Hilbe, R. S. de Souza, and E. E. O. Ishida,  _Bayesian Models for Astrophysical Data Using R, JAGS, Python, and Stan_ (Cambridge University Press, 2017).
  * Baye  _et al._ [1998] D. Baye, P. Descouvemont, and M. Hesse, Phys. Rev. C 58, 545 (1998).
  * Azuma  _et al._ [2010] R. E. Azuma, E. Uberseder, E. C. Simpson, C. R. Brune, H. Costantini, R. J. de Boer, J. Görres, M. Heil, P. J. LeBlanc, C. Ugalde, and M. Wiescher, Phys. Rev. C 81, 045805 (2010).
  * Andreon and Weaver [2015] S. Andreon and B. Weaver,  _Bayesian Methods for the Physical Sciences: Learning from Examples in Astronomy and Physics_, Springer Series in Astrostatistics (Springer International Publishing, 2015).
  * Plummer [2003] M. Plummer, “Jags: A program for analysis of bayesian graphical models using gibbs sampling,” (2003).
  * R Core Team [2020] R Core Team,  _R: A Language and Environment for Statistical Computing_, R Foundation for Statistical Computing, Vienna, Austria (2020).
  * Jennings  _et al._ [1998] B. K. Jennings, S. Karataglidis, and T. D. Shoppa, Phys. Rev. C 58, 579 (1998).
  * Mukhamedzhanov and Nunes [2002] A. Mukhamedzhanov and F. Nunes, Nuclear Physics A 708, 437 (2002).
  * Vernotte  _et al._ [1994] J. Vernotte, G. Berrier-Ronsin, J. Kalifa, R. Tamisier, and B. Wildenthal, Nuclear Physics A 571, 1 (1994).
  * Longland  _et al._ [2010] R. Longland, C. Iliadis, A. E. Champagne, J. R. Newton, C. Ugalde, C. AA, and R. Fitzgerald, Nuclear Physics A 841, 1 (2010).