# Revisiting proton-proton fusion in chiral effective field theory

Bijaya Acharya<EMAIL_ADDRESS>Physics Division, Oak Ridge National Laboratory,
Oak Ridge, TN 37831, USA  Laura Elisa Marcucci<EMAIL_ADDRESS>Department of Physics “E. Fermi”, University of Pisa, Pisa, I-56127, Italy
Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, Pisa, I-56127, Italy
Lucas Platter<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Tennessee, Knoxville, TN 37996, USA  Physics Division, Oak Ridge
National Laboratory, Oak Ridge, TN 37831, USA

###### Abstract

We calculate the $S$-factor for proton-proton fusion using chiral effective
field theory interactions and currents. By performing order-by-order
calculations with a variety of chiral interactions that are regularized and
calibrated in different ways, we assess the uncertainty in the $S$-factor from
the truncation of the effective field theory expansion and from the
sensitivity of the $S$-factor to the short-distance axial current determined
from three- and four-nucleon observables. We find that $S(0)=(4.100\pm
0.019\mathrm{(syst)}\pm 0.013\mathrm{(stat)}\pm 0.008(g_{A}))\times
10^{-23}~{}\mathrm{MeV\,fm}^{2}\,,$ where the three uncertainties arise,
respectively, from the truncation of the effective field theory expansion, use
of the two-nucleon axial current fit to few-nucleon observables and variation
of the axial coupling constant within the recommended range.

##  I Introduction

Nuclear reaction rates are among the main sources of systematic uncertainty in
stellar evolution models Vinyoles et al. (2017). The proton-proton ($pp$)
fusion reaction is the rate-determining step of the $pp$ chains that power the
Sun and lighter stars. Available experimental techniques are not able to
directly measure the rate of this process with sufficient precision at
energies relevant for stellar burning, and values predicted by nuclear theory
provide critical inputs to astrophysical simulations Adelberger et al. (2011).
Smaller or more reliable nuclear-theory uncertainty estimates on this reaction
will help solar models resolve inconsistencies in observed data such as those
that exist between spectroscopic determinations of solar abundances and
helioseismology Bahcall et al. (2005); Serenelli et al. (2009). In concert
with improvements in other physics inputs such as radiative opacities,
rigorous constraints on this rate can also help us use solar data as a
powerful probe of new physics Suliga et al. (2021).

The calculations of this process were traditionally performed using potential
models Bethe and Critchfield (1938); Bahcall and May (1969); Kamionkowski and
Bahcall (1994); Schiavilla et al. (1998). Over the last few decades, methods
based on effective field theory (EFT) techniques Epelbaum et al. (2009) that
allow us to obtain theoretical predictions with reliable uncertainty
estimates, have been employed. EFTs provide a simplified yet rigorous
description of the process under study using only those degrees of freedom
that are relevant at low energies. The calculations are organized as
systematic expansions in the ratio of the typical momentum scale $Q$ of the
process to a large momentum scale $\Lambda_{b}$, beyond with the EFT expansion
breaks down. The undetermined parameters—the so-called low-energy constants
(LECs)—that appear up to a given order in this $Q/\Lambda_{b}$ expansion are
first fixed, e.g., by fitting to experimental data, and then predictive
calculations are performed for other observables. The first applications of
EFT on the $pp$ fusion process were based on hybrid approaches Park et al.
(1998, 2003) that employed wave functions obtained from phenomenological
nuclear potentials along with the nuclear electroweak current operators
derived in chiral EFT, which employs pions and nucleons as dynamical degrees
of freedom. Complete chiral EFT calculations, using potentials and currents
both derived consistently within the chiral EFT framework, have been carried
out, first in Ref. Marcucci et al. (2013) and then in Ref. Acharya et al.
(2016). However, as discussed further below, these calculations require
important updates and corrections.

Over the past few decades, a number of studies have also been performed in
pionless EFT, which uses nucleons as the only dynamical degrees of freedom
Kong and Ravndal (2001); Butler and Chen (2001); Ando et al. (2008); Chen et
al. (2013); De-Leon and Gazit (2022); Acharya et al. (2019). At energies
relevant for astrophysics, the process lies well within the domain of
convergence of pionless EFT; however, the uncertainty is currently larger than
in chiral EFT due to limitations of experiments and Lattice QCD to
sufficiently constrain the LEC $L_{1,A}$ that parametrizes the short-distance
two-body axial current Acharya and Bacca (2020); Davoudi and Kadam (2022).

The goal of this work is to present a state-of-the-art calculation of the $pp$
fusion cross section in chiral EFT. In addition to correcting Marcucci et al.
(2018) the erroneous treatment of the relationship between the three-nucleon
force parameter $c_{D}$ and the two-nucleon axial current parameter
$\hat{d}_{R}$, we also improve and expand upon the work of Ref. Acharya et al.
(2016) by using recently developed chiral EFT interactions Reinert et al.
(2018); Wesolowski et al. (2021), and by accounting for uncertainty sources
not previously considered in Ref. Acharya et al. (2016). We also review and
update the work of Ref. Marcucci et al. (2013) by taking particular care of
the convergence Acharya et al. (2017) in the expansion basis used to calculate
the nuclear wave functions, as well as the range of the integration in the
axial current matrix element, by also correcting the $c_{D}$-$\hat{d}_{R}$
relation, and by using the most recent values for the fundamental constants,
which leads to a reassessment of the Gamow-Teller matrix element of tritium
$\beta$-decay. The same nuclear interaction of Ref. Marcucci et al. (2013) is
implemented, limiting the present study to the case of cutoff $\Lambda=500$
MeV. This calculation also allows us to perform a benchmark study between the
approaches used in Ref. Acharya et al. (2016) and the one of Ref. Marcucci et
al. (2013). The two approaches will be labelled LS and VM, respectively, since
the former used the Lippmann-Schwinger equation to solve the two-body problem,
while the latter used a variational method Marcucci  _et al._ (2012) as
discussed below.

The paper is organized as follows: we outline the theoretical framework in
Section II, present our results in Section III, and conclude with a summary
and outlook in Section IV.

##  II The $S$-factor

The astrophysical $S$-factor $S(E)$ at the $pp$ center-of-mass energy $E$ is
defined as

$S(E)=\sigma(E)\,E\,e^{2\pi\eta}\,,$ (1)

where $\eta=\sqrt{m_{p}/E}\,\alpha/2$ is the Sommerfeld parameter, $m_{p}$ is
the proton mass, $\alpha=1/137.036$ is the fine-structure constant, and
$\sigma(E)$ is the $pp$ fusion cross section at energy $E$, which can be
written as

$\displaystyle\sigma(E)=\int$
$\displaystyle\frac{\mathrm{d}^{3}p_{e}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}p_{\nu}}{(2\pi)^{3}}\frac{1}{2E_{e}}\frac{1}{2E_{\nu}}\,2\pi\delta\left(E+2m_{p}-m_{d}-\frac{q^{2}}{2m_{d}}-E_{e}-E_{\nu}\right)$
$\displaystyle\frac{1}{v_{rel}}\,F(Z,E_{e})\,\frac{1}{4}\sum|\langle
f|\hat{H}_{W}|i\rangle|^{2}\,.$ (2)

Here $p_{e,\nu}$ are the positron and neutrino momenta, $E_{e,\nu}$ their
energies, $m_{d}$ is the deuteron mass, $v_{rel}$ is the $pp$ relative
velocity, and $q$ is the momentum of the recoiling deuteron. The function
$F(Z,E_{e})$ accounts for the distortion of the positron wave function due to
the Coulomb field of the deuteron. Its classical expression, which can be
found in Ref. Feenberg and Trigg (1950), needs to be augmented with radiative
corrections. We enhance $F(Z,E_{e})$ by 1.62% to account for the $\gamma W$
box diagram involving one nucleon, which has been explicitly evaluated in Ref.
Kurylov et al. (2003), and ignore the diagram involving both nucleons, which
has not yet been calculated. The summation in Eq. (II) runs over the spin
projections of all the initial- and the final-state particles. The initial
state $|i\rangle$ and the final state $|f\rangle$ are products of leptonic and
nuclear states and the weak interaction Hamiltonian $\hat{H}_{W}$ can be
written in terms of the leptonic weak current $j^{\mu}$ and the nuclear weak
current $J^{\mu}$ as

$\hat{H}_{W}=\frac{G_{V}}{\sqrt{2}}\int\mathrm{d}^{3}x\left[j^{\mu}(\mathbf{x}){J_{\mu}}^{\dagger}(\mathbf{x})+\mbox{h.c.}\right],$
(3)

where $G_{V}$ is the vector coupling constant. Unless otherwise stated, the
calculations presented below adopt the value $1.149589\times
10^{-11}\,\mathrm{MeV}^{-2}$, that corresponds to the CKM matrix element
$V_{ud}=0.9737$ and the process independent radiative corrections
$\Delta_{R}^{V}=0.02454$ as obtained in a recent reanalysis of the
superallowed beta decay rates Hardy and Towner (2020).

The matrix element of the leptonic weak current operator $j^{\mu}$ between the
leptonic wave functions can be obtained from Dirac algebra and we refer the
reader to Refs. Marcucci et al. (2013); Acharya et al. (2016) for details of
the derivation of the matrix elements of the nuclear weak current operator
$J^{\mu}$ between the nuclear wave functions. The incoming $pp$ system can be
considered exclusively in an $s$-wave, as it has been shown that higher
partial-wave channels are suppressed by several orders of magnitude (see Ref.
Acharya et al. (2019) and also the erratum of Ref. Marcucci et al. (2013)).
Then, the leading contribution to $J^{\mu}$ comes from the one-body (1B)
Gamow-Teller (GT) operator and the leading two-body (2B) corrections occur at
$\mathcal{O}\left([Q/\Lambda_{b}]^{3}\right)$ relative to 1B, which correspond
to next-to-next-to-leading order (NNLO) in the chiral expansion. It should be
noted that these differ slightly from the currents used in Ref. Marcucci et
al. (2013): besides the leading-order GT 1B contribution, relativistic
corrections to the 1B GT term were also included since they are
$\mathcal{O}\left([Q/\Lambda_{b}]^{2}\right)$ in the power counting adopted in
that work. Although extremely small, these terms are retained in the updates
to the study of Ref. Marcucci et al. (2013) presented in this work.

The 2B current we use are the same as the ones used in Refs. Marcucci et al.
(2013); Acharya et al. (2016), but with the corrected Marcucci et al. (2018)
$c_{D}$-$\hat{d}_{R}$ relation (see also Refs. Gazit et al. (2009, 2019)),

$\hat{d}_{R}=-\textstyle{\frac{m_{N}}{4g_{A}\Lambda_{\chi}}}c_{D}+\frac{m_{N}}{3}c_{3}+\frac{2m_{N}}{3}c_{4}+\frac{1}{6}\,,$
(4)

where $m_{N}$ is the nucleon mass, $c_{3,4}$ are pion-nucleon ($\pi N$)
coupling constants and $\Lambda_{\chi}\approx 700$ MeV is the breakdown scale
of chiral perturbation theory. Unless otherwise stated, we use the latest
Particle Data Group (PDG) recommended value of the axial coupling constant
$g_{A}=1.2754\pm 0.0013$ Workman and Others (2022) in the calculations
presented below.

The nuclear wave functions are calculated either by solving the Lippmann-
Schwinger (LS) equations as in Ref. Acharya et al. (2016), or by applying the
variational method (VM) of Refs. Marcucci  _et al._ (2012, 2013). The chiral
EFT Hamiltonians include nuclear as well as electromagnetic potentials. In
this work, relativistic and radiative corrections to the Coulomb interaction,
discussed in Ref. Carlsson et al. (2016), are included explicitly in the
calculation of the wave function, as in Ref. Marcucci et al. (2013). This is
in contrast to Ref. Acharya et al. (2016), where we first calculated $S(E)$
using $pp$ wave functions with nuclear plus Coulomb interaction and applied a
phenomenological correction to take such higher-order electromagnetic effects
into account. Finally, in Ref. Acharya et al. (2016), we used in Eq. (II) the
deuteron mass $m_{d}$ obtained from the calculated value of the binding
energy, which gave a negligible numerical error since all of the interactions
used in Ref. Acharya et al. (2016) were fit to the experimental 2H binding
energy. As we discuss below, since some of the interactions we use in this
work do not reproduce this energy very well, we use the experimental value of
$m_{d}=1875.61294257$ MeV Tiesinga et al. instead. Note that for the
interaction models which properly reproduce the 2H binding energy, the
theoretical and experimental value of $m_{d}$ obviously coincide to sufficient
precision.

##  III Results

We now present our results for $S(E)$ obtained using various chiral EFT
interactions.

###  III.1 The SMS-RS predictions

Table 1 shows the threshold values and energy-derivatives of $S(E)$ calculated
using the SMS-RS potential of Ref. Reinert et al. (2018) at regulator cutoff
$\Lambda=500$ MeV and the LS equations as in Ref. Acharya et al. (2016). In
contrast to older momentum-space chiral interactions Epelbaum et al. (2004a,
b); Entem and Machleidt (2003); Ekström et al. (2015, 2013); Entem et al.
(2017) that employed the conventional non-local Gaussian regulators along with
an additional spectral function regularization Kaiser et al. (1997) in the
two-pion exchange diagrams, the potential of Ref. Reinert et al. (2018) uses a
“semi-local” regularization scheme which consists of a local regulator for the
pion-exchange parts and non-local Gaussian regulator for the contact parts of
the $NN$ interaction in order to preserve the long-range parts that are
unambiguously determined in chiral EFT. The subleading $\pi N$ LECs in this
potential have been fixed to the precise values obtained from a Roy-Steiner
analysis Hoferichter et al. (2015) of the $\pi N$ scattering data. The
nucleon-nucleon ($NN$) LECs in the SMS-RS potential have been fit to mutually
consistent $np$ and $pp$ scattering data of the Granada 2013 database Navarro
Pérez et al. (2013). In the 2B current, we use $\hat{d}_{R}$ given by Eq. (4)
with the value of $c_{D}=-1.626$ obtained in Ref. Maris et al. (2021) by
fitting to $Nd$ scattering cross section data; however, we note that the
theoretical uncertainties in the estimation of $c_{D}$ were not explored in
that reference. It should be remarked that, thanks to the availability of the
SMS-RS potential at various orders in the EFT expansion, it is possible to
perform an order-by-order calculation consistently for both potential and
current. This is different from the study of Ref. Marcucci et al. (2013),
where the nuclear interaction order was fixed at N3LO, while the different
orders were considered only for the chiral expansion of the axial current.

We evaluate $S(E)$ in the energy interval $E=[1,30]$ keV and fit a polynomial
in $E$ to obtain the threshold $S$-factor $S(0)$ and its derivatives. The
extracted values depend on the degree of the fitted polynomial. We obtained
stable values, particularly for $S(0)$ and $S^{\prime}(0)$, for third, fourth
and fifth degree polynomials. Unless otherwise stated, we quote values from
cubic fits.

Order |  $S(0)$ [MeV fm2] |  $S^{\prime}(0)/S(0)$ [MeV-1] |  $S^{\prime\prime}(0)/S(0)$ [MeV-2] |  $S^{\prime\prime\prime}(0)/S(0)$ [MeV-3]
---|---|---|---|---
LO | 4.143 $\times 10^{-23}$ | 10.75 | 306.75 | -5150
NLO | 4.094 $\times 10^{-23}$ | 10.81 | 312.78 | -5370
NNLO | 4.100 $\times 10^{-23}$ | 10.83 | 313.72 | -5382
Table 1: The threshold $S$-factor $S(0)$, and the first and second derivatives
of the $S$-factor at threshold divided by $S(0)$ in units of MeV fm2, MeV-1,
MeV-2 and MeV-3, respectively, at different orders of the chiral EFT expansion
for both the nuclear interaction and current. The SMS-RS interaction is used.

We assess the uncertainty of the $S$-factor due to the EFT truncation by
following the method described in Refs. Furnstahl et al. (2015); Epelbaum et
al. (2015). We first express the $S$-factor as

$S(0)=S_{\rm LO}\sum_{n=0}^{3}c_{n}\left(\frac{Q}{\Lambda_{b}}\right)^{n}~{},$
(5)

where $S_{\rm LO}$ denotes the leading order (LO) result for the $S$-factor
given in Table 1, $Q$ denotes the inherent momentum scale of the problem, and
$\Lambda_{b}$ is the breakdown scale. An estimate of the truncation error is
then obtained by calculating
$(Q/\Lambda_{b})^{4}\max(|c_{0}|,|c_{2}|,|c_{3}|)$. Using the order-by-order
results of Table 1 to obtain the expansion coefficients $c_{n}$, the pion mass
for the momentum scale $Q$ and a conservative estimate of $\Lambda_{b}=500$
MeV, we obtain an uncertainty 0.019 MeV fm2 for the SMS-RS interaction,
leading to the prediction $S(0)=(4.100\pm 0.019\mathrm{(syst)})\times
10^{-23}~{}\rm{MeV\,fm}^{2}$. Here we have used the label “syst” to emphasize
that is is an estimate of the systematic uncertainty from the truncation of
the EFT expansion. We note that this estimate has been found to roughly
correspond to a 68% Bayesian credible interval for a particular choice of
Bayesian priors for the expansion coefficients $c_{n}$ Furnstahl et al.
(2015).

###  III.2 The NNLOsim updates

In Ref. Acharya et al. (2016), we employed the NNLOsim family of 42
interactions Carlsson et al. (2016) to assess the uncertainty in the $pp$
fusion rate due to the statistical uncertainties in the LECs, the systematic
uncertainty due to the chiral EFT cutoff dependence, and the systematic
variations in the database used to calibrate the $NN$ interaction. At each of
the 7 different cutoff values $\Lambda=[450,475,500,525,550,575,600]$ MeV, the
26 LECs in the $NN$ and $\pi N$ sectors were simultaneously fit to 6 different
pools of input data, leading to 42 interactions. Specifically, the input data
consisted of $NN$ scattering data at different truncations of the maximum
scattering energy $T_{\mathrm{lab}}^{\mathrm{max}}$, $\pi N$ cross sections,
the binding energies and charge radii of 2,3H and 3He, the one-body quadrupole
moment of 2H, as well as the comparative $\beta$-decay half-life of 3H. The
NNLOsim interactions have been refit Acharya et al. (2018) to account for the
update in the $c_{D}$-$\hat{d}_{R}$ relation given in Eq. (4) Marcucci et al.
(2018), which enters the calculation of the 3H GT matrix element fit to
tritium $\beta$-decay. We obtain $4.091\times 10^{-23}~{}\rm{MeV\,fm}^{2}$ for
the average $S(0)$ value for the refitted NNLOsim interactions. This is an
upward revision, by about $1.1\%$, partly driven by both the refitting of
$c_{D}$ and the change in the $g_{A}$ value recommended by the PDG Workman and
Others (2022), from the result $S(0)=4.047^{+0.024}_{-0.032}\times
10^{-23}~{}\rm{MeV\,fm}^{2}$ obtained in Ref. Acharya et al. (2016), where the
quoted uncertainty mainly reflected the sensitivity of the $S$-factor to input
data used to calibrate the LECs of the EFT and to the short-distance behavior
of the $NN$ interactions. While a full reanalysis of the uncertainty estimates
of Ref. Acharya et al. (2016) is beyond the scope of this work, it is
instructive to compare the EFT uncertainty obtained above using SMS-RS
interaction to a similar estimate using one of the 42 interactions from the
NNLOsim family. To this end, we choose the interaction with
$T_{\mathrm{lab}}^{\mathrm{max}}=290$ MeV and $\Lambda=500$ MeV, and obtain
$13.537\times 10^{-23}$, $4.869\times 10^{-23}$ and $4.092\times
10^{-23}~{}\rm{MeV\,fm}^{2}$ for $S(0)$ at LO, NLO and NNLO. Using the
procedure discussed above to assess the EFT truncation error, we obtain
$S(0)=(4.092\pm 0.178\mathrm{(syst)})\times 10^{-23}~{}\rm{MeV\,fm}^{2}$. We
note that this large uncertainty is a result of the LO value being rather
large because the deuteron properties are not well reproduced by this
interaction at this order. We therefore consider the estimate of the
systematic EFT truncation error obtained from the SMS-RS interaction to be
more reliable for this problem.

###  III.3 Uncertainty in calibrating the axial contact current

An additional uncertainty that has not been included in the SMS-RS result
quoted above, $S(0)=(4.100\pm 0.019\mathrm{(syst)})\times
10^{-23}~{}\rm{MeV\,fm}^{2}$, is the uncertainty arising from the axial 2B
contact current with $\hat{d}_{R}$ determined using Eq. (4). To estimate this,
we consider the Bayesian posterior probability distribution function (PDF) of
$c_{D}$-$c_{E}$ obtained by Ref. Wesolowski et al. (2021). In that study, as
in the SMS-RS interaction, the $\pi N$ LECs were fixed to central values
determined by the Roy-Steiner analysis in Ref. Hoferichter et al. (2015). The
$NN$ LECs at LO, NLO, and NNLO were then fixed by performing a new fit to the
$np$ and $pp$ scattering data with $T_{\mathrm{lab}}^{\mathrm{max}}=290$ MeV
gathered from the Granada 2013 database Navarro Pérez et al. (2013) as well as
the empirical $nn$ effective range parameters. Constraints on $c_{D}$ and
$c_{E}$ were obtained from fitting to binding energies of 3H and 4He, the
charge radius of 4He, and the GT matrix element of the 3H extracted from
tritium $\beta$-decay followed by marginalization of all other parameters with
both the uncertainty of the chiral EFT Hamiltonian and the uncertainties in
the experimental measurements taken into account.

It was found that, at the older PDG value of $g_{A}=1.2724$ Tanabashi and
Others (2018), the joint PDF of $c_{D}$-$c_{E}$ was well described by a
multivariate $t$-distribution $t_{\nu}(m,S)$ with $\nu\approx 2.8$ degrees of
freedom, a mean vector $m=[-0.0047,-0.1892]$ and a scale matrix of

$S=\begin{bmatrix}0.250&0.043\\\ 0.043&0.008\end{bmatrix}$ (6)

at one standard deviation ($1\sigma$). At least for the energy range
considered in this work, $S(E)$ has approximately linear dependence on the
values of $c_{D}$ within the $1\sigma$ range for this joint PDF. The
uncertainty on $S(E)$ can therefore be very well approximated by sampling
$c_{D}$ from its marginal $t$-distribution, which spans the approximate
$1\sigma$ range $c_{D}\in[-0.615,0.615]$, and computing $S(E)$ at the
$1\sigma$ margins. This gives $S(0)=(4.155\pm 0.013(\mathrm{stat}))\times
10^{-23}$ MeV fm2. Here we used the label “stat” to indicate that this is an
estimate of the statistical uncertainty from the $c_{D}$ probability
distribution. The rather large $S(0)$, in spite of the smaller value of 1.2724
adopted for $g_{A}$, is a result of slow convergence of the deuteron
properties in this EFT expansion up through NNLO. This indicates that the
fitting procedure adopted for this interaction results in large EFT truncation
error even at NNLO for extreme low-energy observables. Nevertheless, we
consider the quoted uncertainty from $c_{D}$ variation to be a reasonable
estimate of the uncertainty in $S(0)$ from fitting the two-body contact axial
current to few-body observables.

###  III.4 Updates to Idaho-N3LO results

We now turn our attention to the study of the $pp$ fusion performed using the
VM to calculate the deuteron and $pp$ wave functions as in Ref. Marcucci et
al. (2013) from the Idaho-N3LO potential with $\Lambda=500$ MeV Entem and
Machleidt (2003). With respect to Ref. Marcucci et al. (2013), we have put
particular attention to the convergence on the basis expansion for the
deuteron wave function, and on the integration range used to calculate the
transition operator matrix element. In particular, we have verified that not
only the deuteron binding energy, but also the asymptotic normalization
constants and the so-called $D/S$-state ratio $\eta=A_{D}/A_{S}$ are well
reproduced. Furthermore, we have verified that the integral range of
$r_{max}=50$ fm is sufficient to reach convergence of the results. In fact, by
going from $r_{max}=50$ fm to $r_{max}=60$ fm, the change in $S(0)$ is beyond
the third decimal digit. Therefore, the results we are going to present are at
convergence at least up to the the third decimal digit of $S(0)$. We note that
such an accuracy was not the primary goal of the work of Ref. Marcucci et al.
(2013), where the results were obtained with a theoretical accuracy of the
order of 1%.

In order to use the most recent values for the fundamental constants entering
the calculations, we, first of all, update the value for the experimental GT
matrix element in tritium $\beta$-decay $\langle GT_{exp}\rangle$. Note that
$\langle GT_{exp}\rangle$ and the $A=3$ binding energies were the observables
of choice used to fix the LEC $c_{D}$, and consequently $\hat{d}_{R}$, and the
LEC $c_{E}$, that parametrizes the three-nucleon contact interaction entering
at NNLO (see, for instance, Ref. Marcucci  _et al._ (2012)). The GT matrix
element is defined as

$\langle GT\rangle^{2}=\bigg{[}\frac{2ft_{0^{+}\rightarrow
0^{+}}}{ft_{{}^{3}{\rm H}}}-\langle
F\rangle^{2}\bigg{]}\,\frac{1}{f\,g_{A}^{2}}\ ,$ (7)

where $g_{A}$ is the single-nucleon axial coupling constant,
$f=f_{A}/f_{V}=1.00529$ is the ratio of the axial and vector Fermi functions,
$ft_{{}^{3}{\rm H}}$ and $ft_{0^{+}\rightarrow 0^{+}}$ are the $ft$-values for
tritium $\beta$-decay and for superallowed $0^{+}\rightarrow 0^{+}$
transitions. From $\langle GT\rangle$, we define $\langle GT_{exp}\rangle$ as
$\langle GT_{exp}\rangle={\langle GT\rangle}/{\sqrt{3}}$. In the present
calculation, we have used $g_{A}=1.2754\pm 0.0013$, according to the PDG, and
$ft_{0^{+}\rightarrow 0^{+}}=(3072.24\pm 1.85)$ s, according to Ref. Hardy and
Towner (2020). This value is consistent with the value of $G_{V}$ used in Eq.
(3). Furthermore, we have used $\langle F\rangle^{2}=0.99967$, as it was
obtained in Refs. Marcucci  _et al._ (2012) with the Idaho-N3LO potentials.
This value is quite different from the one obtained using the phenomenological
AV18/UIX interaction model, and used for instance in Ref. Schiavilla et al.
(1998), $\langle F\rangle^{2}=0.9987$. However, we have verified that the
impact on $\langle GT_{exp}\rangle$ of this different $\langle F\rangle^{2}$
is negligible. Finally, we have adopted for $ft_{{}^{3}{\rm H}}$ three
different values: $(1134.6\pm 3.1)$ s as obtained in Ref. Simpson (1987),
$(1129.6\pm 3)$ s as obtained in Ref. Akulov and Mamyrin (2005), and
$(1132.1\pm 4.3)$ s as used already in Ref. Marcucci  _et al._ (2012) and
obtained averaging these two values and summing the errors in quadrature. The
$\langle GT_{exp}\rangle$ obtained with these three values for $ft_{{}^{3}{\rm
H}}$ will be labelled Fit-1, Fit-2 and Fit-3, respectively. The three values
for $\langle GT_{exp}\rangle$, together with their uncertainties arising from
the experimental errors on $g_{A}$, $ft_{0^{+}\rightarrow 0^{+}}$ and
$ft_{{}^{3}{\rm H}}$, are listed in Table 2. In the table, we report also the
corresponding ranges for $c_{D}$ obtained, as mentioned above, with the
fitting procedure of Ref. Marcucci  _et al._ (2012).

| $\langle GT_{exp}\rangle$ | $c_{D}$
---|---|---
Fit-1 | $0.9488\pm 0.0019$ | 0.1836–0.7680
Fit-2 | $0.9514\pm 0.0019$ | 0.3830–0.9690
Fit-3 | $0.9501\pm 0.0024$ | 0.2833–0.8685
Table 2: Values for $\langle GT_{exp}\rangle$ as obtained using three
different values for $ft_{{}^{3}{\rm H}}$, $1134.6\pm 3.1$ s labelled Fit-1,
$1129.6\pm 3$ s labelled Fit-2, and $1132.1\pm 4.3$ labelled Fit-3,
respectively. The corresponding ranges for the LEC $c_{D}$ obtained using the
Idaho-N3LO potential with $\Lambda=500$ MeV are also listed.

With the new values for the LEC $c_{D}$, or equivalently $\hat{d}_{R}$, we
present in Table 3 the results for the zero-energy astrophysical $S$-factor
$S(0)$, and its first, second and third derivatives, $S^{\prime}(0)/S(0)$,
$S^{\prime\prime}(0)/S(0)$ and $S^{\prime\prime\prime}(0)/S(0)$. We list in
the table also the results obtained with the same constants and the same
potential using the LS equation and the axial current of Ref. Acharya et al.
(2016). By inspecting the values in the table, we can conclude that (i) the VM
and LS calculations are in very good agreement as the small differences are
consistent with the anticipated size of the relativistic corrections to the GT
1B operator which were neglected in the LS calculations but were retained in
the VM calculations; (ii) the dependence on the adopted value for
$ft_{{}^{3}{\rm H}}$ is very weak; (iii) the theoretical uncertainty arising
from the range of $c_{D}$ values allowed for each fitting procedure, and
therefore ultimately related to the experimental error on $\langle
GT_{exp}\rangle$, is extremely small for $S(0)$ and
$S^{\prime\prime}(0)/S(0)$, and beyond the second decimal digit for
$S^{\prime}(0)/S(0)$; (iv) the calculated $S(0)$ is in excellent agreement
with the $S(0)$ value obtained with the SMS-RS and NNLOsim interactions. (v)
The value for $S(0)$ is $\sim 2.2$% larger than the value of Ref. Marcucci et
al. (2013). The difference arises from the improved convergence of the nuclear
wave functions adopted in this work, as well as from the increased values for
the single-nucleon axial coupling constant $g_{A}$.

We finally note that we have also used a second-degree polynomial to fit the
$S(E)$ obtained with the Idaho-N3LO interaction. With the data obtained using
the VM method this gives, $S(0)=4.105(1)\times 10^{-23}$ MeV fm2,
$S^{\prime}(0)/S(0)=11.44(1)$ MeV-1, and $S^{\prime\prime}(0)/S(0)=233.8$
MeV-2 for the three fits employed in this work.

| Method | $S(0)$ | $S^{\prime}(0)/S(0)$ | $S^{\prime\prime}(0)/S(0)$ | $S^{\prime\prime\prime}(0)/S(0)$
---|---|---|---|---|---
Fit-1 | VM | 4.115(4) | 10.60 | 347.1 | -6908
| LS | 4.101(4) | 10.83 | 313.8 | -5382
Fit-2 | VM | 4.118(4) | 10.60 | 347.1 | -6907
| LS | 4.104(4) | 10.83 | 313.8 | -5381
Fit-3 | VM | 4.117(4) | 10.60 | 347.1 | -6908
| LS | 4.104(4) | 10.83 | 313.8 | -5382
Table 3: Values for the zero-energy astrophysical $S$-factor $S(0)$ (in
$10^{-23}$ MeV fm2), its first, second and third derivatives
$S^{\prime}(0)/S(0)$ (in MeV-1), $S^{\prime\prime}(0)/S(0)$ (in MeV-2) and
$S^{\prime\prime\prime}(0)/S(0)$ (in MeV-3) obtained with the Idaho-N3LO
potential with $\Lambda=500$ MeV, using the variational method (VM) of Ref.
Marcucci et al. (2013), or the Lippmann-Schwinger (LS) equation as in Ref.
Acharya et al. (2016). The results labelled Fit-1, Fit-2, and Fit-3 are
obtained consistently with the values of $c_{D}$ listed in Table 2. The
numbers in parentheses are the theoretical error arising from the range of
$c_{D}$ values allowed for each fitting procedure. For $S^{\prime}(0)/S(0)$
these errors are beyond the quoted digits.

##  IV Summary and Outlook

We calculated the proton-proton fusion rate using various chiral interactions:
the Idaho-N3LO interaction at regulator cutoff of 500 MeV Entem and Machleidt
(2003), the NNLOsim family of 42 interactions fit at 7 different regulator
cutoffs to 6 different pools of input data, and two additional sets of
interactions—one non-local Wesolowski et al. (2021) and another “semi-local”
Reinert et al. (2018)—in which the pion-nucleon coupling constants are fixed
to precise values determined from Roy-Steiner analysis. We obtained an
estimate for the uncertainty from the truncation of the EFT expansion and also
assessed the uncertainty from fixing the LEC $\hat{d}_{R}$ in the two-body
axial current using $c_{D}$ fitted to $A\geq 3$ data. Furthermore, by using
the Idaho-N3LO potential, we have performed a benchmark calculation between
the two different approaches first applied to the $pp$ reaction in Refs.
Marcucci et al. (2013) and Acharya et al. (2016). We have used recent values
of the fundamental constants and have used the corrected relationship between
the axial current LEC $\hat{d}_{R}$ and the three-body force parameter
$c_{D}$.

The threshold $S$-factor $S(0)$ obtained from the various EFT interactions
indicate an upward revision from the recommendation made by Ref. Adelberger et
al. (2011). It is remarkable that the $S$-factor value obtained using the SMS-
RS interactions Reinert et al. (2018) is in excellent agreement with the
NNLOsim result although its $c_{D}$ value was obtained without fitting to
electroweak data. We also show that variation in $c_{D}$ within the 68%
Bayesian credible interval obtained from the joint $c_{D}-c_{E}$ distribution
calculated in Ref. Wesolowski et al. (2021) translates to an uncertainty of
$0.013~{}\mathrm{MeV}~{}\mathrm{fm}^{2}$ in $S(0)$, which is smaller than the
uncertainty of $0.019~{}\mathrm{MeV}~{}\mathrm{fm}^{2}$ stemming from the
truncation of the chiral expansion in the SMS-RS potential. We do note,
however, that the potential of Ref. Wesolowski et al. (2021) underbinds the
deuteron at all orders (by about 75% at LO, 15% at NLO and 1% at NNLO) which
causes $S(E)$ to be overpredicted, even with just one-body current. We
encounter this issue also at LO and NLO in the NNLOsim interactions, which
again impacts the extraction of truncation error. This highlights the
importance of ensuring that the low-energy properties of the $NN$ system, e.g.
the binding energy and the asymptotic normalization coefficient of the
deuteron, are reasonably well reproduced while fitting the LECs of chiral EFT.
Given the difficulty in assessing the truncation error of the EFT expansion
due to slow order-by-order convergence in the non-local chiral interactions
used in this work, we consider the central value obtained using the SMS-RS
interaction to be more reliable and the corresponding truncation error
estimate to be better calibrated. Combining the systematic EFT truncation
error with the statistical uncertainty that arises from the PDF of $c_{D}$ and
the range obtained by varying $g_{A}$ within the PDG recommendation of
$1.2754\pm 0.0013$, we recommend

$S(0)=(4.100\pm 0.019\mathrm{(syst)}\pm 0.013\mathrm{(stat)}\pm
0.008(g_{A}))\times 10^{-23}~{}\mathrm{MeV\,fm}^{2}\,.$ (8)

Note that this value is consistent with the ones obtained using the NNLOsim
and Idaho-N3LO interactions. Furthermore, we find with the SMS-RS potential
$S^{\prime}(0)/S(0)$ and $S^{\prime\prime}(0)/S(0)$ values of
$10.83~{}\mathrm{MeV}^{-1}$ and $313.72~{}\mathrm{MeV}^{-2}$, respectively.
The values of these energy derivatives of $S(E)$ depend on the interaction
used, calculation method and degree of the fitted polynomial. We refrain from
performing a detailed analysis of uncertainties for them as they are not
relevant at currently achievable precision in the energy range of interest for
solar conditions.

Finally, we would like to remark that the $S$-factor is not accurately
predicted by the chiral EFT interactions of Refs. Carlsson et al. (2016);
Wesolowski et al. (2021) at low chiral orders if the deuteron bound-state
properties are not adequately reproduced. This highlights the importance of
calibrating the interactions by fitting to experimental data using strategies
that also lead to systematic convergence pattern for very low-energy
observables such as $pp$ fusion at solar conditions. Bayesian model mixing
between chiral and pionless EFT to ensure that the deuteron and the $pp$
effective-range parameters are well reproduced appears to be a promising
strategy for reliably predicting solar $pp$ fusion rate. Furthermore, we
believe that a thorough study of the $pp$ fusion reaction should be performed
using the largest possible variety of chiral EFT interactions, local or non-
local, regularized in co-ordinate or momentum space, with or without the
inclusion of $\Delta$-isobar degrees of freedom, possibly available a
different chiral orders, as, for instance, those of Refs. Piarulli et al.
(2015); Nosyk et al. (2021); Entem et al. (2017); Saha et al. (2022). Such a
study is beyond the scope of this work, but is definitely highly recommended.
Work along this line, on the footsteps of Refs. Ceccarelli et al. (2023);
Gnech et al. (2023), is currently underway.

## Acknowledgments

We would like to thank Andreas Ekström, Rocco Schiavilla and Michele Viviani
for useful suggestions and comments. This work has been supported by the
National Science Foundation under Grant Nos. PHY-1555030 and PHY-2111426, by
the Office of Nuclear Physics, U.S. Department of Energy under Contract No.
DE-AC05-00OR22725, and by the Office of High Energy Physics, U.S. Department
of Energy under Contract No. DE-AC02-07CH11359 through the Neutrino Theory
Network Fellowship awarded to BA.

## References

  * Vinyoles et al. (2017) N. Vinyoles, A. M. Serenelli, F. L. Villante, S. Basu, J. Bergström, M. C. Gonzalez-Garcia, M. Maltoni, C. Peña Garay, and N. Song, Astrophys. J. 835, 202 (2017), arXiv:1611.09867 [astro-ph.SR] .
  * Adelberger et al. (2011) E. G. Adelberger et al., Rev. Mod. Phys. 83, 195 (2011), arXiv:1004.2318 [nucl-ex] .
  * Bahcall et al. (2005) J. N. Bahcall, S. Basu, M. Pinsonneault, and A. M. Serenelli, Astrophys. J. 618, 1049 (2005), arXiv:astro-ph/0407060 .
  * Serenelli et al. (2009) A. Serenelli, S. Basu, J. W. Ferguson, and M. Asplund, Astrophys. J. Lett. 705, L123 (2009), arXiv:0909.2668 [astro-ph.SR] .
  * Suliga et al. (2021) A. M. Suliga, S. Shalgar, and G. M. Fuller, Journal of Cosmology and Astroparticle Physics 2021, 042 (2021).
  * Bethe and Critchfield (1938) H. A. Bethe and C. L. Critchfield, Phys. Rev. 54, 248 (1938).
  * Bahcall and May (1969) J. N. Bahcall and R. M. May, Astrophys. J. 155, 501 (1969).
  * Kamionkowski and Bahcall (1994) M. Kamionkowski and J. N. Bahcall, Astrophys. J. 420, 884 (1994), arXiv:astro-ph/9305020 .
  * Schiavilla et al. (1998) R. Schiavilla et al., Phys. Rev. C 58, 1263 (1998), arXiv:nucl-th/9808010 .
  * Epelbaum et al. (2009) E. Epelbaum, H. W. Hammer, and U.-G. Meissner, Rev. Mod. Phys. 81, 1773 (2009), arXiv:0811.1338 [nucl-th] .
  * Park et al. (1998) T.-S. Park, K. Kubodera, D.-P. Min, and M. Rho, Astrophys. J. 507, 443 (1998), arXiv:astro-ph/9804144 .
  * Park et al. (2003) T. S. Park et al., Phys. Rev. C67, 055206 (2003), arXiv:nucl-th/0208055 .
  * Marcucci et al. (2013) L. E. Marcucci, R. Schiavilla, and M. Viviani, Phys. Rev. Lett. 110, 192503 (2013), [Erratum: Phys.Rev.Lett. 123, 019901 (2019)], arXiv:1303.3124 [nucl-th] .
  * Acharya et al. (2016) B. Acharya, B. D. Carlsson, A. Ekström, C. Forssén, and L. Platter, Phys. Lett. B 760, 584 (2016), arXiv:1603.01593 [nucl-th] .
  * Kong and Ravndal (2001) X. Kong and F. Ravndal, Phys. Rev. C 64, 044002 (2001), arXiv:nucl-th/0004038 .
  * Butler and Chen (2001) M. Butler and J.-W. Chen, Phys. Lett. B 520, 87 (2001), arXiv:nucl-th/0101017 .
  * Ando et al. (2008) S. Ando, J. W. Shin, C. H. Hyun, S. W. Hong, and K. Kubodera, Phys. Lett. B 668, 187 (2008), arXiv:0801.4330 [nucl-th] .
  * Chen et al. (2013) J.-W. Chen, C. P. Liu, and S.-H. Yu, Phys. Lett. B 720, 385 (2013), arXiv:1209.2552 [nucl-th] .
  * De-Leon and Gazit (2022) H. De-Leon and D. Gazit, (2022), arXiv:2207.10176 [nucl-th] .
  * Acharya et al. (2019) B. Acharya, L. Platter, and G. Rupak, Phys. Rev. C 100, 021001 (2019), arXiv:1903.02590 [nucl-th] .
  * Acharya and Bacca (2020) B. Acharya and S. Bacca, Phys. Rev. C 101, 015505 (2020), arXiv:1911.12659 [nucl-th] .
  * Davoudi and Kadam (2022) Z. Davoudi and S. V. Kadam, Phys. Rev. D 105, 094502 (2022), arXiv:2111.11599 [hep-lat] .
  * Marcucci et al. (2018) L. E. Marcucci, A. Kievsky, S. Rosati, R. Schiavilla, and M. Viviani, Phys. Rev. Lett. 121, 049901 (2018).
  * Reinert et al. (2018) P. Reinert, H. Krebs, and E. Epelbaum, Eur. Phys. J. A 54, 86 (2018), arXiv:1711.08821 [nucl-th] .
  * Wesolowski et al. (2021) S. Wesolowski, I. Svensson, A. Ekström, C. Forssén, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. C 104, 064001 (2021), arXiv:2104.04441 [nucl-th] .
  * Acharya et al. (2017) B. Acharya, A. Ekström, D. Odell, T. Papenbrock, and L. Platter, Phys. Rev. C 95, 031301 (2017), arXiv:1608.04699 [nucl-th] .
  * Marcucci  _et al._ (2012) L. E. Marcucci, A. Kievsky, S. Rosati, R. Schiavilla, and M. Viviani, Phys. Rev. Lett. 108, 052502 (2012).
  * Feenberg and Trigg (1950) E. Feenberg and G. Trigg, Reviews of Modern Physics 22, 399 – 406 (1950).
  * Kurylov et al. (2003) A. Kurylov, M. J. Ramsey-Musolf, and P. Vogel, Phys. Rev. C 67, 035502 (2003), arXiv:hep-ph/0211306 .
  * Hardy and Towner (2020) J. C. Hardy and I. S. Towner, Phys. Rev. C 102, 045501 (2020).
  * Gazit et al. (2009) D. Gazit, S. Quaglioni, and P. Navrátil, Phys. Rev. Lett. 103, 102502 (2009).
  * Gazit et al. (2019) D. Gazit, S. Quaglioni, and P. Navrátil, Phys. Rev. Lett. 122, 029901 (2019).
  * Workman and Others (2022) R. L. Workman et al. (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Carlsson et al. (2016) B. D. Carlsson, A. Ekström, C. Forssén, D. F. Strömberg, G. R. Jansen, O. Lilja, M. Lindby, B. A. Mattsson, and K. A. Wendt, Phys. Rev. X 6, 011019 (2016), arXiv:1506.02466 [nucl-th] .
  * (35) E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, CODATA internationally recommended 2018 values of the fundamental physical constants. NIST Standard Reference Database http://physics.nist.gov/constants (2019) .
  * Epelbaum et al. (2004a) E. Epelbaum, W. Gloeckle, and U.-G. Meissner, Eur. Phys. J. A19, 125 (2004a), arXiv:nucl-th/0304037 .
  * Epelbaum et al. (2004b) E. Epelbaum, W. Gloeckle, and U.-G. Meissner, Eur. Phys. J. A19, 401 (2004b), arXiv:nucl-th/0308010 .
  * Entem and Machleidt (2003) D. R. Entem and R. Machleidt, Phys. Rev. C68, 041001 (2003), arXiv:nucl-th/0304018 .
  * Ekström et al. (2015) A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, Phys. Rev. C 91, 051301 (2015).
  * Ekström et al. (2013) A. Ekström, G. Baardsen, C. Forssén, G. Hagen, M. Hjorth-Jensen, G. R. Jansen, R. Machleidt, W. Nazarewicz, T. Papenbrock, J. Sarich, and S. M. Wild, Phys. Rev. Lett. 110, 192502 (2013).
  * Entem et al. (2017) D. R. Entem, R. Machleidt, and Y. Nosyk, Phys. Rev. C 96, 024004 (2017), arXiv:1703.05454 [nucl-th] .
  * Kaiser et al. (1997) N. Kaiser, R. Brockmann, and W. Weise, Nucl. Phys. A 625, 758 (1997), arXiv:nucl-th/9706045 .
  * Hoferichter et al. (2015) M. Hoferichter, J. Ruiz de Elvira, B. Kubis, and U.-G. Meißner, Phys. Rev. Lett. 115, 192301 (2015), arXiv:1507.07552 [nucl-th] .
  * Navarro Pérez et al. (2013) R. Navarro Pérez, J. E. Amaro, and E. Ruiz Arriola, Phys. Rev. C 88, 064002 (2013), [Erratum: Phys.Rev.C 91, 029901 (2015)], arXiv:1310.2536 [nucl-th] .
  * Maris et al. (2021) P. Maris et al., Phys. Rev. C 103, 054001 (2021), arXiv:2012.12396 [nucl-th] .
  * Furnstahl et al. (2015) R. J. Furnstahl, N. Klco, D. R. Phillips, and S. Wesolowski, Phys. Rev. C 92, 024005 (2015).
  * Epelbaum et al. (2015) E. Epelbaum, H. Krebs, and U.-G. Meißner, Phys. Rev. Lett. 115, 122301 (2015).
  * Acharya et al. (2018) B. Acharya, A. Ekström, and L. Platter, Phys. Rev. C 98, 065506 (2018), arXiv:1806.09481 [nucl-th] .
  * Tanabashi and Others (2018) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Simpson (1987) J. J. Simpson, Phys. Rev. C 35, 752 (1987).
  * Akulov and Mamyrin (2005) Y. Akulov and B. Mamyrin, Phys. Lett. B 610, 45 (2005).
  * Piarulli et al. (2015) M. Piarulli, L. Girlanda, R. Schiavilla, R. Navarro Pérez, J. E. Amaro, and E. Ruiz Arriola, Phys. Rev. C 91, 024003 (2015), arXiv:1412.6446 [nucl-th] .
  * Nosyk et al. (2021) Y. Nosyk, D. R. Entem, and R. Machleidt, Phys. Rev. C 104, 054001 (2021), arXiv:2107.06452 [nucl-th] .
  * Saha et al. (2022) S. K. Saha, D. R. Entem, R. Machleidt, and Y. Nosyk, (2022), arXiv:2209.13170 [nucl-th] .
  * Ceccarelli et al. (2023) L. Ceccarelli, A. Gnech, L. E. Marcucci, M. Piarulli, and M. Viviani, Front. Phys. 10, 1049919 (2023), arXiv:2209.09762 [nucl-th] .
  * Gnech et al. (2023) A. Gnech, L. E. Marcucci, and M. Viviani, in preparation. (2023).