Near Threshold Proton-Proton Fusion in Effective Field Theory

The astrophysical S-factor for proton-proton fusion, S_11(E), is obtained with the nuclear matrix element analytically calculated in pionless effective field theory. To the third order, the zero-energy result S_11(0) and the first energy derivative S'_11(0) are found to be (3.99 \pm 0.14)* 10^-25 MeV b and S_11(0)*(11.3 \pm 0.1) MeV^-1, respectively; both consistent with the current adopted values. The second energy derivative is also calculated for the first time, and the result S"_11(0) = S_11(0)*(170 \pm 2) MeV^-2 contributes at the level of 0.5% to the fusion rate at the solar center, which is smaller than 1% as previously estimated.


I. INTRODUCTION
As the trigger of stellar hydrogen burning, the proton-proton (pp) fusion reaction, p+p → d+e + +ν e , plays a fundamental role in astrophysics. Its reaction rate, conventionally encoded by the astrophysical S-factor S 11 (E), with E the energy in the center of mass frame, is thus an important input in studies such as stellar evolution and solar neutrinos.
For such a pivotal reaction, the value of S 11 (E) at the center of the Sun (with temperature ∼ 1.55×10 7 K which yields E ∼ a few keV), unfortunately can not be determined or reliably extrapolated by terrestrial experiments. This is because they have to be performed at much higher energies to overcome the Coulomb barrier and get sensible statistics. Therefore, one has to rely on theory for predictions.
Ever since the first proposal of the pp chains by Bethe and Critchfield [1], the nuclear transition amplitude of the pp fusion, together with other solar fusion cross sections have been extensively studied. Those results were extensively reviewed first in Ref. [2] (SFI) and then most recently in Ref. [3] (SFII). For solar fusion, the temperature is much lower than typical energy scales in nuclear physics, thus only the first few terms in the expansion of S 11 (E) around E = 0 is needed for solar models. Currently, the recommended value of S 11 (0) in SFII has an error of 1%, which results in a 1% error in the pp fusion rate. The recommended value of S ′ 11 (0) also has an error of 1%, which only results in a less than 0.1% error in the fusion rate. S ′′ 11 (0) was not given in SFII. However, Bahcall and May [4] estimated its contribution to the rate to be ∼ 1%, comparable to the overall error in S 11 .
Thus, SFII recommended a modern calculation of S ′′ 11 (0) be undertaken. In this work, we calculate S 11 (E) using the nuclear effective field theory (EFT) with pions integrated out. One of the major ingredients, the square of the orbital matrix element Λ 2 (E) which S 11 (E) is proportional to, is extracted from the cross section of an analogue process: ν e + d → p + p + e − , analytically computed in Ref. [5].
This pionless EFT is applicable for low energy processes with the characteristic momentum p much smaller than the pion mass m π [6][7][8], which is the case for solar pp fusion. In this theory, pions are integrated out. All the nucleon-nucleon interactions and two-body currents are described by point-like contact interactions with a systematic expansion in powers of p/m π . A close analogy of this theory is the Fermi theory of four fermion contact interactions.
The one-and two-body contributions both depend on the momentum cut-off but the sum does not. In pionless EFT, there is only one two-body current (with coupling L 1,A ) in every weak interaction deuteron breakup process to next-to-next-to-leading (NNLO) in the p/m π expansion [5]. This two-body current is a Gamow-Teller operator. The other two-body currents are either missing due to vector current conservation or the matrix elements are suppressed because of the orthogonality of the initial and final state wave functions in the zero recoil limit. This means the universal number L 1,A encodes the two-body contributions for all low energy weak deuteron breakup processes, and it takes just one measurement to calibrate all the processes. This feature is also seen in the other complimentary approaches to the solar fusion such as potential models [9], hybrid version of EFT [10], and pionless EFT with dibaryon (see SFII for more details).

II. S 11 (E) FROM PIONLESS EFT
The astrophysical S-factor, S(E), for a nuclear reaction at kinetic energy E is related to the reaction cross section by The rapid-varying energy dependence of σ(E) due to the Coulomb barrier is mostly accounted for in the exponential term which is controlled by the Sommerfeld parameter where α = 1/137.036 is the fine structure constant and m p = 938.272 MeV is the proton mass. 1 Using the same convention and inputs as SFII, we write the S-factor for the pp fusion as The energy dependence of S 11 (E) is determined by two terms: (i) Λ(E), the orbital matrix element, which is proportional to the nuclear transition matrix element, and (ii) f R pp (E), the phase space factor in this nuclear β + process. We note the conventional way of separating the radiative correction in S 11 (E): the long-distance (so-called "outer") part, which is processdependent, is included in the phase factor f R pp (E) (annotated by a superscript "R"), while the short distance (so-called "inner") part, which is process-independent, is taken into account by the product of (g A /g V ) 2 and 1/(f t) 0 + →0 + . For a detailed account, see Ref. [11].
At very low energy, the pp fusion predominantly goes from the 1 S 0 partial wave state to the deuteron state ( 3 S 1 with some 3 D 1 mixture) through the spatial axial current operator.
In pionless EFT, it takes the form [12] A where N denotes the nucleon field; all σ's and τ 's are the Pauli matrices for spin and isospin, respectively.
The coupling constant L 1,A of the leading axial two-body current appears at next-to-leading order (NLO) and there is no new two-body current contributing until next-to-next-to-nextto-leading order (N 3 LO). The orbital matrix element Λ(p), with p = m p E, is then related to the nuclear matrix elements of A − k by where j denotes the deuteron |d polarization and is the square of the Sommerfeld factor. The effect of the deuteron recoil, with momentum q 0.4 MeV, is suppressed by a factor q 2 /γ 2 < 10 −4 and hence can be neglected. In the zero recoil limit the vector matrix element between pp and d vanishes because the wave functions are orthogonal.
In Ref. [5], various (anti)neutrino deuteron breakup processes were computed in pionless EFT up to N 2 LO with analytic expressions. The needed hadronic matrix element for the pp fusion can be extracted from the result of the similar process v e + d → p + p + e − . The square of the orbital matrix element can be simplified and expressed as where ρ d = 1.764 fm is the effective range parameter for deuteron, and the corresponding energy-dependent structure functions are The complex integral The p i components of the pp scattering amplitude in the 1 S 0 channel, A i (p), are found to be where a = −7.82 fm and r 0 = 2.79 fm are the scattering length and effective range parameter for 1 S 0 , respectively; and the function depends on momentum via the Sommerfeld parameter η(p). The renormalization scale µindependent coupling constants L ′ 1,A is defined as and when being evaluated at the pion-mass scale, µ = m π ≈ 140 MeV, the low energy constants C 0,−1 = −3.77 fm 2 and C 2,−2 = 7.50 fm 4 . 2 Note that we keep the small expansion parameter ǫ explicit in Eqs. (7a-7d) to make the order of each term transparent. For evaluation, one has to make a series expansion in ǫ up to the maximum order at which the result is valid, and then set ǫ = 1 in the end.
At zero energy, our result of Λ 2 (0) can be written in a compact form where the function yields a value of 1.4655 with χ ≡ α m p /γ = 0.1498. Up to the order of ǫ 2 , this result is in agreement with Ref. [12] in which Λ(0) is worked out to the order of ǫ 4 ; and the constant L 1,A defined therein is exactly the same as our L ′ 1,A . This result is also consistent with Ref. [13] by ignoring the L ′ 1,A term, and Ref. [14] by a redefinition of L ′ 1,A .
The phase space factor f (E) in a nuclear β process is conventionally written as (see, e.g., Ref. [15]) where m e = 0.511 MeV is the electron (and positron) mass; W and p e = W 2 − m 2 e are the relativistic energy and momentum of the emitted electron or positron; Z is the charge of the final nuclear state with the sign "±" being assigned to β ∓ emission, respectively; and 2 As we take the q 2 /γ 2 → 0 limit in F 1 , F 4 , and B 0 , and re-define G 2 and L 1,A by some normalization factors, we add a prime to remind readers about these changes from Ref. [5].
Q is the difference of the total rest mass of the initial and final nuclear states. The Fermi where γ 0 ≡ (1 − Z 2 α 2 ) 1/2 and v e = ±Z α W/p e take both the relativity of electron and the finite size of nucleus (through a spherical radius R) into account.
For pp fusion, Z = 1, R = 2.1402 fm, Q = 0.420 MeV, and the outer radiative correction evaluated to be δ out pp = 1.62% [16,17], so the "corrected" phase space factor becomes The linear term in E is in excellent agreement with [4]. A 0.1% error is assign to f R pp (0) in SFII. An error of the same size is assigned to f R pp (E) in this work.

III. RESULTS AND DISCUSSION
To obtain a numerical result of Λ 2 (p), the last piece of information we need is the value of L 1,A . Using the reactorν e d breakup data, L 1,A is determined to be (3.6 ±4.6) fm 3 [18,19], and using the data of solar neutrino deuteron breakups through charged and neutral currents, and v e e − elastic scattering, L 1,A is determined to be (4.0±6.3) fm 3 [19]. Both fits have taken into account the radiative corrections computed in Refs. [16,17]. Treating these two fits in two nucleon systems as independent and using the same average scheme as in Ref. [18], we obtain where a −0.3 fm 3 shift of the central value is introduced due to updating g A from 1.267 to 1.2695. This range of L 1,A is consistent with the nave dimensional analysis value |L 1,A | ∼ 6 fm 3 [20]. One can further constrain this two-body current by the tritium beta decay [3], which is a three nucleon system. This is carried out in potential model [9] and hybrid EFT [10], and both approaches yielded the determination of S 11 (0) with 1% accuracy.
The two sources of errors are correlated in the sense that when two near-threshold weak deuteron breakup processes have the same kinematics, the two processes will be governed by the same Gamow-Teller matrix element such that higher order effects can be largely absorbed by shifting the value of L 1,A . Indeed, the two nucleon processes used to constrain L 1,A and the pp fusion have similar kinematics -even though not exactly the same, so the combined error on Λ 2 (0) should be between 3% (100% correlated) to 4.2% (not correlated, added in quadrature). We will assign the combined error to be 3.6% which is slightly larger than the 3% error adopted in SFII with pionless EFT determination. Note that in SFII, a 0.9% error is assign to Λ 2 (0) by using tritium beta decay rate to constrain the two body current. Here we just use the constraint of Eq.(17) available in the two nucleon sector and study its implication in S 11 (E).
It is remarkable that c 1 and c 2 are quite insensitive to L 1,A . This is because for near threshold pp fusion, the energy dependence is dominated by the initial state pp scattering which alone gives c 1 = 3.09 and c 2 = 36.2 without L 1,A dependence. c 1 and c 2 then receive −28% and −6% L 1,A -dependent corrections, respectively, to reach the values of Eq. (20).
This explains why the errors due to the uncertainty of L 1,A in c 1 and c 2 are 23% and 4% of that in Λ 2 (0). This result along with the combined error assignment discussed above leads to error reduction in S ′ 11 (0)/S 11 (0) and S ′′ 11 (0)/S 11 (0). Combining our results of Λ 2 (E) and f pp (E), we found the S-factor S 11 (E) for the low- Our central value of S 11 (0) is the same as the pionless EFT value in SFII while the error is slightly larger as discussed above. As for S ′ 11 (0), our result is consistent with the adopted value in SFII: S ′ (SFII) 11 (0) = (11.2 ± 0.1) S 11 (0)MeV −1 , which was obtained by Bahcall and May [4]. Finally, we also report a value of S ′′ 11 (0), which was not computed before, with ∼ 1% accuracy. At the center of the sun, this higher-order term has a contribution at the level of 0.5% to the fusion rate, which is smaller than ∼ 1% that was previously estimated [3,4].