Complex-energy analysis of proton-proton fusion

TL;DR

This paper provides a model-independent parametrization of the astrophysical S factor for proton-proton fusion, revealing a maximum near 0.13 MeV and linking it to the resonance pole, with implications for solar and early universe conditions.

Contribution

It introduces a new, model-independent parametrization of the matrix element Lambda for proton-proton fusion, connecting the S factor maximum to the resonance pole and enabling precise derivative calculations.

Findings

Lambda reaches a maximum near 0.13 MeV

The parametrization is valid up to about 5 MeV

Accurate derivatives of S at zero energy are obtained

Abstract

An analysis of the astrophysical $S$ factor of the proton-proton weak capture ($\mathrm{p}+\mathrm{p}\rightarrow {}^2\mathrm{H}+\mathrm{e}^++\nu_{\mathrm{e}}$) is performed on a large energy range covering solar-core and early Universe temperatures. The measurement of $S$ being physically unachievable, its value relies on the theoretical calculation of the matrix element $\Lambda$. Surprisingly, $\Lambda$ reaches a maximum near $0.13~\mathrm{MeV}$ that has been unexplained until now. A model-independent parametrization of $\Lambda$ valid up to about $5~\mathrm{MeV}$ is established on the basis of recent effective-range functions. It provides an insight into the relationship between the maximum of $\Lambda$ and the proton-proton resonance pole at $(-140-467\,\mathrm{i})~\mathrm{keV}$ from analytic continuation. In addition, this parametrization leads to an accurate evaluation of the derivatives of $\Lambda$, and hence of $S$, in the limit of zero energy.

Tables

Table 1. Table 1: Shape parameters in different potential models for the S-state two-nucleons wave functions.

	$s_d$	${\nu }_d$	$s_p$	${\nu }_p$
Av18 Wiringa et al. (1995)	$8.43$	$3.61$	$9.81$	$4.22$
Reid93 Stoks et al. (1994)	$8.45$	$3.63$	$9.55$	$3.92$
NijmII Stoks et al. (1994)	$8.78$	$4.00$	$9.59$	$3.95$

Table 2. Table 2: Numerical values of the functions f p ​ ( w ¯ ) subscript 𝑓 𝑝 ¯ 𝑤 \displaystyle f_{p}(\bar{w}) and their derivatives with respect to w ¯ ¯ 𝑤 \displaystyle\bar{w} at E = 0 𝐸 0 \displaystyle E=0 . The upper index n 𝑛 \displaystyle n is the order of the derivatives.

	$n=0$	$n=1$	$n=2$
$f_0^{\left(n\right)}$	$0.14827\left(2\right)$	$0.68187\left(8\right)$	$2.3574\left(2\right)$
$f_1^{\left(n\right)}$	$0.27417\left(4\right)$	$1.1233\left(1\right)$	$3.3682\left(2\right)$
$f_2^{\left(n\right)}$	$0.64955\left(7\right)$	$2.2505\left(2\right)$	$5.4045\left(3\right)$

Equations

$$\sigma (E)=\frac{3m_{\mathrm{e}}{c}^2(\lambda g{)}^2}{{\pi }^{2}Ek}F(E+Q){\left|{\int }_0^{\infty }{u}_{\mathrm{d}}(r)u_{\mathrm{pp}}(E,r)\mathrm{d}r\right|}^2,$$

$$\Lambda (E)=\frac{a_{\mathrm{B}}{b}^{3/2}}{\sqrt{8}}\frac{2\eta }{C_{\eta ,0}}{\int }_0^{\infty }{u}_{\mathrm{d},0}(r)u_{\mathrm{pp},0}(E,r)\mathrm{d}r,$$

$$\sigma (E)=\frac{1}{({\mathrm{e}}^{2\pi \eta }-1)E}\frac{12m_{\mathrm{e}}{c}^2(\lambda g{)}^2}{\pi a_{\mathrm{B}}{b}^3}F(E+Q){\left|\Lambda (E)\right|}^2.$$

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

$$S_{\mathrm{std}}(E)={\mathrm{e}}^{2\pi \eta }E\sigma (E),$$

$$E\le {\left(\frac{2\pi }{\ln (100)}\right)}^2\mathrm{Ry}\simeq 1.86\mathrm{Ry},$$

$$S(E)=\frac{12m_{\mathrm{e}}{c}^2(\lambda g{)}^2}{\pi a_{\mathrm{B}}{b}^3}F(E+Q){\left|\Lambda (E)\right|}^2.$$

$$u_{\mathrm{pp},0}(r)\stackrel{r\to \infty }{\to }{F}_{\eta ,0}\cos {\delta }_0+G_{\eta ,0}\sin {\delta }_0,$$

$$u_{\mathrm{d},0}(r)\stackrel{r\to \infty }{\to }A{\mathrm{e}}^{-br}.$$

$$u_{\mathrm{d},0}(r)=A(1-{\mathrm{e}}^{-s_{\mathrm{d}}br}{)}^{{\nu }_{\mathrm{d}}}{\mathrm{e}}^{-br}.$$

$$A=\sqrt{N_{\mathrm{S}}\frac{2b\Gamma (2{\nu }_{\mathrm{d}}+2s_{\mathrm{d}}^{-1}+1)}{\Gamma (2{\nu }_{\mathrm{d}}+1)\Gamma (2s_{\mathrm{d}}^{-1}+1)}},$$

$$u_{\mathrm{pp},0}(r)=(1-{\mathrm{e}}^{-s_{\mathrm{p}}br}{)}^{{\nu }_{\mathrm{p}}}(F_{\eta ,0}\cos {\delta }_0+G_{\eta ,0}\sin {\delta }_0).$$

$$u_{\mathrm{d},0}{u}_{\mathrm{pp},0}=\sum _{i,j}^{\infty }{c}_{i,j}{\mathrm{e}}^{-{\beta }_{i,j}br}(F_{\eta ,0}\cos {\delta }_0+G_{\eta ,0}\sin {\delta }_0).$$

$${\beta }_{i,j}=1+i{s}_{\mathrm{d}}+j{s}_{\mathrm{p}},$$

$$c_{i,j}=(-1{)}^{i+j}(\genfrac{}{}{0ex}{}{{\nu }_{\mathrm{d}}}{i}\left)\right(\genfrac{}{}{0ex}{}{{\nu }_{\mathrm{p}}}{j}).$$

$$F_{\eta ,\ell }(kr)=\frac{C_{\eta ,\ell }\Gamma (2\ell +2)}{(2\eta {)}^{\ell +1}}{\Phi }_{\eta ,\ell }(kr),$$

$$C_{\eta ,\ell }=\frac{(2\eta {)}^{\ell }}{\Gamma (2\ell +2)}\sqrt{\frac{2\pi \eta w_{\eta ,\ell }}{{\mathrm{e}}^{2\pi \eta }-1}}.$$

$$w_{\eta ,\ell }=\prod _{j=0}^{\ell }\left(1+\frac{j^2}{{\eta }^2}\right)\text{and}{w}_{\eta ,0}=1.$$

$$G_{\eta ,\ell }=\frac{C_{\eta ,\ell }\Gamma (2\ell +2)}{(2\eta {)}^{\ell +1}}\frac{{\mathrm{e}}^{2\pi \eta }-1}{\pi }\left(\frac{{\Psi }_{\eta ,\ell }}{w_{\eta ,\ell }}+g_{\eta ,\ell }{\Phi }_{\eta ,\ell }\right),$$

$$g_{\eta ,\ell }=\frac{\psi (\ell +1+\mathrm{i}\eta )+\psi (\ell +1-\mathrm{i}\eta )}{2}-\ln \eta ,$$

$$\begin{array}{rl}& F_{\eta ,0}\cos {\delta }_0+G_{\eta ,0}\sin {\delta }_0\\ & =\frac{C_{\eta ,0}/(2\eta )}{\left|{\Delta }_0^{+}(E)\right|}\left(\frac{a_{\mathrm{B}}}{2}{\varkappa }_0{\Phi }_{\eta ,0}+{\Psi }_{\eta ,0}\right).\end{array}$$

$${\varkappa }_0(E)=\frac{2}{a_{\mathrm{B}}}\left(\frac{\pi \cot {\delta }_0}{{\mathrm{e}}^{2\pi \eta }-1}+g_{\eta ,0}\right)\simeq \frac{-1}{{\alpha }_0}+\frac{r_0}{2}\frac{m_{\mathrm{p}}}{{\hslash }^2}E.$$

$${\Delta }_0^{+}(E)=\left(\frac{a_{\mathrm{B}}}{2}{\varkappa }_0-g_{\eta ,0}\right)-\frac{\mathrm{i}\pi }{{\mathrm{e}}^{2\pi \eta }-1}.$$

$${\left|{\Delta }_0^{+}(E)\right|}^2={\left(\frac{a_{\mathrm{B}}}{2}{\varkappa }_0-g_{\eta ,0}\right)}^2+{\left(\frac{\pi }{{\mathrm{e}}^{2\pi \eta }-1}\right)}^2,$$

$${\varphi }_{\beta ,0}(ϵ)$$

$${\psi }_{\beta ,0}(ϵ)$$

$$\phi_{\beta,0}(\epsilon)=\frac{\chi}{\beta^{2}+\epsilon}\exp\!\big{(}\chi\tfrac{\arctan(\sqrt{\epsilon}/\beta)}{\sqrt{\epsilon}}\big{)}\>,$$

$${\psi }_{\beta ,0}(ϵ)\simeq \Gamma (-1,\frac{\chi }{\beta }){\varphi }_{\beta ,0}(ϵ),$$

$$\tilde{\Lambda }=\frac{A/\sqrt{2b}}{\chi \left|{\Delta }_0^{+}\right|}\sum _{i,j}^{\infty }{c}_{i,j}\left(\frac{a_{\mathrm{B}}}{2}{\varkappa }_0{\varphi }_{{\beta }_{i,j},0}+{\psi }_{{\beta }_{i,j},0}\right).$$

$$\tilde{\Lambda }=\frac{A/\sqrt{2b}}{\chi \left|{\Delta }_0^{+}\right|}\left[\left(\frac{a_{\mathrm{B}}}{2}{\varkappa }_0+\gamma \right){\varphi }_{1,0}(ϵ)-\tilde{C}(ϵ)\right],$$

$$\tilde{C}(ϵ)=-\sum _{i,j\ne 0,0}^{\infty }{c}_{i,j}\left(\frac{a_{\mathrm{B}}}{2}{\varkappa }_0{\varphi }_{{\beta }_{i,j},0}+{\psi }_{{\beta }_{i,j},0}\right).$$

Full text

I Introduction

The proton-proton fusion reaction ($\displaystyle\mathrm{p}+\mathrm{p}\rightarrow{}^{2}\mathrm{H}+\mathrm{e}^{+}+\nu_{\mathrm{e}}$), also known as the proton-proton weak capture, is a fundamental process in nuclear astrophysics. It is the starting point of the proton-proton chain for stellar nucleosynthesis in hydrogen-burning stars. Its cross section $\displaystyle\sigma(E)$ is usually expressed in terms of the astrophysical factor $\displaystyle S(E)$ at the two-proton center-of-mass energy $\displaystyle E$. Unfortunately, this cross section is so small at typical astrophysical temperature ($\displaystyle E\lesssim 0.01~{}\mathrm{MeV}$), that a reliable measurement cannot be achieved with enough statistics, even above the Coulomb barrier (about $\displaystyle 0.2~{}\mathrm{MeV}$). A theoretical prediction of $\displaystyle S$ is therefore required.

The first calculation of $\displaystyle S$ at zero energy was proposed by Bethe and Critchfield Bethe and Critchfield (1938). They also introduced the dimensionless weak capture matrix element $\displaystyle\Lambda$ at zero angular momentum from which they deduced $\displaystyle S$. Thereafter, the accuracy of $\displaystyle S(0)$ was improved by Salpeter Salpeter (1952) and Bahcall and his coworkers Bahcall and May (1969) using effective-range theory. In the 1990s, several authors calculated $\displaystyle S(0)$ from the nucleon-nucleon wave functions computed in potential models Bahcall and Pinsonneault (1992); Kamionkowski and Bahcall (1994). More recently, systematic computations of $\displaystyle S(0)$ were performed in pionless effective field theory from next-to-leading order (NLO) of the momentum expansion Kong and Ravndal (2001); Ando et al. (2008) up to N4LO Butler and Chen (2001). In parallel, efforts were made in chiral effective field theory to reduce the uncertainty on $\displaystyle S(0)$ by adding two-body corrections to the Gamow-Teller operator adjusted with data for tritium $\displaystyle\beta$-decay Carlson et al. (1991); Park et al. (2003); Schiavilla et al. (1998); Marcucci et al. (2013); *Marcucci2014a. As the uncertainty on $\displaystyle S(0)$ has diminished since early works, the small contribution of its energy derivatives $\displaystyle S^{\prime}(0)$ and $\displaystyle S^{\prime\prime}(0)$ has become important for nuclear astrophysics Adelberger et al. (2011); Marcucci et al. (2013); *Marcucci2014a; Chen et al. (2013); Acharya et al. (2016). In this regard, Adelberger et al. recommended a calculation of $\displaystyle S^{\prime\prime}(0)$ to be undertaken Adelberger et al. (2011).

The calculation of these derivatives raises the question of whether $\displaystyle S$ is analytic in the neighborhood of zero energy. Such an analysis is still missing in the literature. Yet, some peculiarities of the proton-proton scattering are known. In 1980, Kok highlighted the presence of a sub-threshold resonance pole at about $\displaystyle(-140-467\,\mathrm{i})~{}\mathrm{keV}$ in the $\displaystyle{}^{1}S_{0}$ channel Kok (1980). Up to a complex phase, this pole lies in an energy range corresponding to early Universe temperatures below the nucleosynthesis freeze-out point ($\displaystyle E\lesssim 1~{}\mathrm{MeV}$) Tanabashi et al. (2018). Therefore, the derivatives of $\displaystyle S$ are likely to be influenced by the relative closeness of this pole. In addition, it is known that the Coulomb interaction between two protons generates a sequence of poles at complex energy which accumulate to the zero-energy point Kok (1980); Gaspard and Sparenberg (2018). These poles may affect the polynomial extrapolation used in Refs. Schiavilla et al. (1998); Marcucci et al. (2013); *Marcucci2014a; Acharya et al. (2016) to determine $\displaystyle S(0)$, $\displaystyle S^{\prime}(0)$, and $\displaystyle S^{\prime\prime}(0)$ from the numerical computation of $\displaystyle S(E)$ on an energy range that does not include $\displaystyle E=0$. In this context, the issue is to understand the influence of all these structures on the astrophysical $\displaystyle S$ factor.

The purpose of this work is to obtain an efficient model-independent parametrization of $\displaystyle S$, valid on a large energy range ($\displaystyle E\lesssim 5~{}\mathrm{MeV}$), and able to impose constraints on its series expansion at $\displaystyle E=0$. Our parametrization must also describe the resonance pole at $\displaystyle(-140-467\,\mathrm{i})~{}\mathrm{keV}$ and the Coulomb poles. To do so, we resort to a recently introduced effective-range function (ERF), namely the $\displaystyle\Delta$ function Ramírez Suárez and Sparenberg (2017); Gaspard and Sparenberg (2018); Baye and Brainis (2000), which has only been considered useful for heavier systems until now Blokhintsev et al. (2017); *Blokhintsev2018a; *Blokhintsev2018b. This approach is motivated by the efficiency of the ERFs at describing the energy-dependent shape of the proton-proton wave function up to a few MeVs. It is aimed to reach a much better accuracy on the values of $\displaystyle S(0)$, $\displaystyle S^{\prime}(0)$, and $\displaystyle S^{\prime\prime}(0)$ than with the polynomial extrapolation performed in Refs. Schiavilla et al. (1998); Marcucci et al. (2013); *Marcucci2014a; Acharya et al. (2016). Finally, our results will be verified with the nucleon-nucleon wave functions computed in different local potential models, namely Av18 Wiringa et al. (1995), Reid93 Stoks et al. (1994), and NijmII Stoks et al. (1994). In particular, these potentials are intended to validate the model independence of our parametrization.

This paper is organized as follows. Section II presents the parametrization of the astrophysical $\displaystyle S$ factor, first from analytical approximations of the nucleon-nucleon wave functions in Subsec. II.2, and then in a model-independent way in Subsec. II.3. The analysis of $\displaystyle S$ at complex energies is discussed in Sec. III and supplemented by graphical illustrations. The numerical values of the logarithmic derivatives of $\displaystyle S$ are shown in Sec. IV. Section V is devoted to a conclusion. Detailed calculations of the Fermi phase-space integral and the Coulomb integrals are given in Appendices A and B, respectively.

II Proton-proton weak capture

II.1 Astrophysical $\displaystyle S$ factor

Let us start by defining the astrophysical $\displaystyle S$ factor studied in this paper. After integrating out the emitted leptons, the proton-proton weak capture cross section is known to be given to a good approximation by Salpeter (1952); Bahcall and May (1969); Baye (2013)

[TABLE]

where $\displaystyle k=\sqrt{m_{\rm p}E}/\hbar$ is the proton-proton wave number, and $\displaystyle u_{\rm d}$ and $\displaystyle u_{\rm pp}$ are the radial $\displaystyle S$-wave components of the deuteron and proton-proton wave functions respectively. The contribution from higher-order partial waves can be neglected in the low-energy approximation Acharya et al. (2019). The parameter $\displaystyle\lambda=1.2724(23)$ is the weak axial/vector ratio, $\displaystyle g=G_{\rm F}\left|V_{\rm ud}\right|(m_{\rm e}c^{2})^{2}/(\hbar c)^{3}=2.96707(64)\times 10^{-12}$ is the dimensionless weak coupling constant for neutrons, $\displaystyle G_{\rm F}$ is the Fermi constant of muon decay, and $\displaystyle V_{\rm ud}$ is the first element of the CKM quark mixing matrix. All the fundamental constants are taken from Ref. Tanabashi et al. (2018). The quantity $\displaystyle F(E+Q)$ in Eq. (1) is known as the Fermi phase-space integral, which accounts for the electric repulsion of the emitted positron in a relativistic framework. It is discussed in more detail in Appendix A.

It turns out that the overlap integral between $\displaystyle u_{\rm d}$ and $\displaystyle u_{\rm pp}$ in Eq. (1) vanishes in the zero energy limit ($\displaystyle E\rightarrow 0$). This vanishing behavior originates from the cancellation of the proton-proton radial wave function $\displaystyle u_{\rm pp}$ at $\displaystyle E=0$ due to the Coulomb barrier. This behavior can be factored out of the overlap integral in defining the dimensionless orbital matrix element $\displaystyle\Lambda(E)$ as Salpeter (1952); Bahcall and May (1969)

[TABLE]

where $\displaystyle C_{\eta,0}=\sqrt{2\pi\eta/(\mathop{}\!\mathrm{e}^{2\pi\eta}-1)}$ is the Coulomb normalization coefficient, $\displaystyle\eta=1/(a_{\mathrm{B}}k)$ is the Sommerfeld parameter, and $\displaystyle a_{\mathrm{B}}=\hbar c/(\alpha m_{\rm p}c^{2}/2)=57.6398~{}\mathrm{fm}$ is the proton-proton Bohr radius. The other constants are the deuteron binding wave number $\displaystyle b=\sqrt{2m_{\rm pn}B_{\rm d}}/\hbar=0.231606~{}\mathrm{fm}^{-1}$ Bahcall and May (1969); Adelberger et al. (2011), the proton-neutron reduced mass $\displaystyle m_{\rm pn}=m_{\rm p}m_{\rm n}/(m_{\rm p}+m_{\rm n})$, and the deuteron binding energy $\displaystyle B_{\rm d}=2.22457~{}\mathrm{MeV}$. In contrast to the overlap integral, $\displaystyle\Lambda(E)$ has a finite limit at $\displaystyle E=0$. With these notations, the cross section (1) becomes after some simplifications

[TABLE]

This expression (3) suggests the most natural definition of the astrophysical $\displaystyle S$ factor, that is

[TABLE]

This definition (4) of $\displaystyle S$ is assumed in this work. Strictly speaking, the definition (4) does not reduce to the original definition due to Salpeter Salpeter (1952)

[TABLE]

and nowadays considered as standard in stellar astrophysics. We point out that Eq. (4) is considered in Sec. II of Bahcall’s and May’s paper Bahcall and May (1969), even though they actually defined $\displaystyle S$ by Eq. (5).

Anyway, at sufficiently low energy, the definitions (4) and (5) coincide. We notice indeed that the approximation $\displaystyle\mathop{}\!\mathrm{e}^{2\pi\eta}-1\simeq\mathop{}\!\mathrm{e}^{2\pi\eta}$ is valid within less than $\displaystyle 1\%$ for

[TABLE]

where $\displaystyle\mathrm{Ry}$ is the nuclear Rydberg energy which is equal to $\displaystyle\alpha^{2}m_{\rm p}c^{2}/4=12.4911~{}\mathrm{keV}$ in the proton-proton system. Since we consider energies much higher than the Rydberg energy in this work, definition (4) is preferred. Thus, we understand Eq. (5) as the low-energy approximation of Eq. (4).

Furthermore, it should be noted that definition (4) does not affect any of the values $\displaystyle S(0)$, $\displaystyle S^{\prime}(0)$, and $\displaystyle S^{\prime\prime}(0)$ presented in the main text compared to $\displaystyle S_{\rm std}$. Indeed, the relative error between the two definitions displays an essential singularity at $\displaystyle E=0$ which cancels all its derivatives in the limit $\displaystyle E\xrightarrow{>}0$. Therefore, the derivatives of $\displaystyle S$ at $\displaystyle E=0$ obtained in this work are necessarily equal to the derivatives of $\displaystyle S_{\rm std}$ at $\displaystyle E=0$.

Finally, using Eqs. (3) and (4), one finds the expression

[TABLE]

The currently recommended value of $\displaystyle S(0)$ is $\displaystyle 4.01(4)\times 10^{-23}~{}\mathrm{MeV}\,\mathrm{fm}^{2}$ Adelberger et al. (2011).

II.2 Wave function-based parametrization

Most of the energy dependence in the overlap integral (2) comes from the proton-proton wave function $\displaystyle u_{\rm pp}$. This function is normalized such that it tends to a sine wave of unit amplitude for $\displaystyle r\rightarrow\infty$. Its asymptotic behavior reads

[TABLE]

where $\displaystyle F_{\eta,\ell}(kr)$ and $\displaystyle G_{\eta,\ell}(kr)$ are the standard Coulomb wave functions Olver et al. (2010), and $\displaystyle\delta_{0}$ is the proton-proton ${}^{1}S_{0}$ phase shift. Similarly, the asymptotic behavior of the deuteron ${}^{3}S_{1}$ wave function is

[TABLE]

The normalization coefficient $\displaystyle A$ is found to be $\displaystyle 0.8850(5)~{}\mathrm{fm}^{-1/2}$ using three different potential models (Av18, Reid93, and NijmII) Wiringa et al. (1995); Stoks et al. (1994). It turns out that the asymptotic regime in Eqs. (8) and (9) is already reached for $\displaystyle r\gtrsim 2~{}\mathrm{fm}$. Since the spatial extent of the deuteron is much larger than $\displaystyle 2~{}\mathrm{fm}$, it is relevant Salpeter (1952); Bahcall and May (1969) to calculate (2) from these asymptotic behaviors.

Inserting the asymptotic behaviors (8) and (9) into Eq. (2) leads to the Laplace transforms of the Coulomb functions, which are known analytically in terms of Gauss hypergeometric functions. However, this calculation does not consider the short-range behavior of the nucleon-nucleon wave functions due to the nuclear potential.

One efficient way of including the short-range contribution in $\displaystyle\Lambda$ is to use analytical approximations of the nucleon wave functions. This method was used by many authors, especially to approximate the deuteron wave function on the basis of a series of exponential functions. Such approximations are known as Hulthén-type wave functions Moravcsik (1958); Kottler and Kowalski (1964); McGee (1966); Humberston and Wallace (1970); Oteo (1988).

We propose to use the following approximation of the ${}^{3}S_{1}$ bound state of the deuteron:

[TABLE]

The shape parameters $\displaystyle s_{\rm d}$ and $\displaystyle\nu_{\rm d}$ are fitted to the deuteron wave function. It is worth noting that $\displaystyle A$ is related to the shape parameters $\displaystyle s_{\rm d}$ and $\displaystyle\nu_{\rm d}$ in Eq. (10) by the normalization condition

[TABLE]

where $\displaystyle N_{\rm S}$ is the $\displaystyle{}^{3}S_{1}$-state probability given by $\displaystyle\left\langle u_{{\rm d},0}\middle|u_{{\rm d},0}\right\rangle=94.30(6)\%$ for the three potential models (Av18, Reid93, NijmII). The fitted values of the shape parameters subject to the constraint (11) are shown in Table 1. As one can see in Fig. 1(a), the approximation (10) of the deuteron wave function is remarkably accurate. The root-mean-square deviation from the Reid93 $\displaystyle{}^{3}S_{1}$ wave function is about $\displaystyle 0.015$. This accuracy is good enough for our needs.

The same kind of parametrization can be applied to the proton-proton wave function:

[TABLE]

In contrast to the deuteron wave function, this wave function depends on the energy of the incoming protons. Therefore, the shape parameters $\displaystyle s_{\rm p}$ and $\displaystyle\nu_{\rm p}$ are likely to vary with the energy. However, we will neglect these variations on the considered energy range because of the great depth of the nuclear potentials. The fitted values of the shape parameters are shown in Table 1 for different potential models.

The advantage of the expressions (10) and (12) is the reduction to the known Laplace transforms of the Coulomb functions, if the shape factors $\displaystyle(1-\mathop{}\!\mathrm{e}^{-s_{\rm d}x})^{\nu_{\rm d}}$ and $\displaystyle(1-\mathop{}\!\mathrm{e}^{-s_{\rm p}x})^{\nu_{\rm p}}$ are expanded in binomial series. The product function then becomes

[TABLE]

The integer indices $\displaystyle i$ and $\displaystyle j$ in the expansion (13) run over the terms of the binomial expansion of $\displaystyle(1-\mathop{}\!\mathrm{e}^{-s_{\rm d}x})^{\nu_{\rm d}}$ and $\displaystyle(1-\mathop{}\!\mathrm{e}^{-s_{\rm p}x})^{\nu_{\rm p}}$ respectively. The variable $\displaystyle\beta_{i,j}$ takes the values

[TABLE]

and the corresponding coefficient is

[TABLE]

One important issue with the overlap integral (2) is the strongly vanishing behavior of the proton-proton wave function when the energy decreases. The wave function $\displaystyle u_{{\rm pp},0}(r)$, normalized according to Eq. (8), tends to zero as $\displaystyle\eta^{-1/2}\mathop{}\!\mathrm{e}^{-\pi\eta}$ Baye (2013). This behavior is due to the asymptotic normalization of the Coulomb wave functions $\displaystyle F_{\eta,\ell}$ and $\displaystyle G_{\eta,\ell}$ that are constrained to sine wave of unit amplitude. In Eq. (2), this cancellation is compensated by the prefactor $\displaystyle 2\eta/C_{\eta,0}$. In order to factorize the cancellation of $\displaystyle u_{{\rm pp},0}(r)$ out of the overlap integral, we rewrite the Coulomb wave functions $\displaystyle F_{\eta,\ell}$ and $\displaystyle G_{\eta,\ell}$ of Eq. (8) in terms of the modified Coulomb functions $\displaystyle\Phi_{\eta,\ell}$ and $\displaystyle\Psi_{\eta,\ell}$ introduced in Refs. Gaspard and Sparenberg (2018); Gaspard (2018). Contrary to the standard Coulomb functions, these functions have the advantage of being analytic in the complex plane of the energy. In particular, they tend toward nonzero functions of $\displaystyle r/a_{\mathrm{B}}$ at zero energy. The function $\displaystyle\Phi_{\eta,\ell}$ is related to $\displaystyle F_{\eta,\ell}$ by Gaspard (2018)

[TABLE]

where $\displaystyle C_{\eta,\ell}$ is the general Coulomb normalization coefficient that reads Olver et al. (2010); Gaspard (2018)

[TABLE]

In Eq. (17), the function $\displaystyle w_{\eta,\ell}$ is defined by

[TABLE]

The function $\displaystyle G_{\eta,\ell}$ can be expressed in terms of $\displaystyle\Phi_{\eta,\ell}$ and $\displaystyle\Psi_{\eta,\ell}$ as Gaspard (2018)

[TABLE]

where $\displaystyle g_{\eta,\ell}$ is the Bethe function given by

[TABLE]

and $\displaystyle\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$ is the digamma function Bethe (1949); Gaspard and Sparenberg (2018). Combining Eqs. (16) and (19), the asymptotic behavior of the proton-proton wave function (8) becomes

[TABLE]

In contrast to $\displaystyle F_{\eta,\ell}$ and $\displaystyle G_{\eta,\ell}$, the modified Coulomb functions $\displaystyle\Phi_{\eta,\ell}$ and $\displaystyle\Psi_{\eta,\ell}$ are free of singularity at zero energy. This property is crucial in the analysis of $\displaystyle\Lambda$ at complex energy, especially near $\displaystyle E=0$. In Eq. (21), $\displaystyle\varkappa_{0}$ is the standard ERF of the proton-proton $\displaystyle{}^{1}S_{0}$ scattering. Its first-order expansion in $\displaystyle E$ provides an accurate parametrization of the phase shift $\displaystyle\delta_{0}$ over a large energy range ($\displaystyle E\lesssim 5~{}\mathrm{MeV}$)

[TABLE]

The parameters $\displaystyle\alpha_{0}$ and $\displaystyle r_{0}$ are respectively the scattering length and the effective range. The modulus of the modified ERF $\displaystyle\Delta^{+}_{0}$ Gaspard and Sparenberg (2018); Hamilton et al. (1973) also appears in Eq. (21). This phase-shift-dependent function is related to $\displaystyle\varkappa_{0}$ by

[TABLE]

However, in contrast to $\displaystyle\varkappa_{0}$, the function $\displaystyle\Delta^{+}_{0}$ is singular at $\displaystyle E=0$ mostly because of $\displaystyle g_{\eta,0}$. The bracket in Eq. (23) is also called the $\displaystyle\Delta_{0}$ function in Ref. Gaspard and Sparenberg (2018). In practical computation, $\displaystyle\Delta^{+}_{0}$ can be evaluated from the knowledge of the effective-range function $\displaystyle\varkappa_{0}$. The square modulus of $\displaystyle\Delta^{+}_{0}$ can be expressed as

[TABLE]

as it will also enter the parametrization of $\displaystyle\Lambda$. The expression (24) is the analytic continuation of $\displaystyle\left|\Delta^{+}_{0}\right|^{2}$ to the complex plane of the energy. Finally, all these ERFs are depicted in Fig. 2.

Using the expansion (13) of the product function, the overlap integral (2) splits into a series of Laplace transforms of the modified Coulomb functions $\displaystyle\Phi_{\eta,0}$ and $\displaystyle\Psi_{\eta,0}$. These Laplace transforms are defined as follows

[TABLE]

where $\displaystyle x=br$ is the dimensionless radial coordinate, $\displaystyle\epsilon=k^{2}/b^{2}=E/B$ is the dimensionless proton energy, and $\displaystyle B=2m_{\rm pn}B_{\rm d}/m_{\rm p}$ is the deuteron binding energy corrected for the neutron-proton mass difference. As shown in Appendix B, the regular Coulomb integral is exactly given by

[TABLE]

where the dimensionless constant $\displaystyle\chi=2/(a_{\mathrm{B}}b)=0.149816$ is due to Bahcall and May Bahcall and May (1969). The Laplace transform of $\displaystyle\Psi_{\eta,0}$ cannot be obtained in a simple form. However, it is explained in Appendix B that, as far as $\displaystyle\chi$ is small with respect to $\displaystyle 1$, the following approximation holds for $\displaystyle E\lesssim 5~{}\mathrm{MeV}$

[TABLE]

where $\displaystyle\Gamma(a,z)$ denotes the upper incomplete gamma function Olver et al. (2010).

Therefore, combining Eqs. (13) and (21) into Eq. (2), we get the expansion

[TABLE]

The tilde over $\displaystyle\Lambda$ means that the result (28) assumes the analytical approximations (10) and (12).

It is useful to separate the first term of the expansion (28), for which $\displaystyle\beta_{0,0}=1$, because it corresponds to the contribution of the far-field part of the nucleon wave functions. One thus expects this contribution to be larger than the short-range part of the wave functions Bahcall and May (1969). We find convenient to introduce an energy-dependent function, denoted $\displaystyle\tilde{C}$, which gathers all the contributions of the short-range part of the wave functions. The expansion (28) is thus written as

[TABLE]

where $\displaystyle\gamma=\Gamma(-1,\chi)=4.28065$ is independent of the energy, and

[TABLE]

The minus sign in front of the series of Eq. (30) makes $\displaystyle\tilde{C}$ a positive function.

The shape parameters of the nucleon wave functions in Table 1 provide an estimation of $\displaystyle\tilde{C}$. At zero energy, we get from Eq. (30) the estimate $\displaystyle\tilde{C}(0)\simeq 0.253$. The function $\displaystyle\tilde{C}$ built on Eq. (30) is shown as curve (a) in Fig. 3.

II.3 Model-independent parametrization

Instead of relying on analytical approximations of the nucleon-nucleon wave functions, we can establish a general model-independent parametrization of $\displaystyle\Lambda$ inspired by Eq. (29)

[TABLE]

Indeed, this expression comes from the separation of the far-field part of the nucleon-nucleon wave functions from its short-range part $\displaystyle C$, and this far-field part is not supposed to change from a potential model to another. Consequently, the parametrization (31) is model independent, as far as the function $\displaystyle C$ can be fitted. The function $\displaystyle C$ is defined by inverting Eq. (31)

[TABLE]

but computing $\displaystyle\Lambda$ with Eq. (2) on the basis of numerical wave functions. At zero energy, we find $\displaystyle C(0)=0.260(1)$; the error being due to the uncertainty on potential models.

The function $\displaystyle C$ in Eq. (32) is shown as curve (b) in Fig. 3. Although they are close within less than $\displaystyle 4\%$ on the considered energy range, curves (a) and (b) display different shapes. We note that curve (b) deviates from a straight line, in contrast to curve (a). This difference is due to the limitation of the assumption in Eq. (12) that the short-range part of the proton-proton wave function, i.e., the shaping factor $\displaystyle(1-\mathop{}\!\mathrm{e}^{-s_{\rm p}br})^{\nu_{\rm p}}$, does not depend on the energy. Therefore, the use of function $\displaystyle C$ is limited to relatively low energy ($\displaystyle\lesssim 100~{}\mathrm{keV}$) where it can be considered linear. In order to model the behavior over a large energy range ($\displaystyle\lesssim 5~{}\mathrm{MeV}$), we establish from Eqs. (26) and (31) a parametrization of $\displaystyle\Lambda$ that is more convenient to practical applications

[TABLE]

where the function $\displaystyle L$ can be related to $\displaystyle C$ with

[TABLE]

In contrast to $\displaystyle C$, the function $\displaystyle L$ is pretty close to a smoothly varying straight line, as shown in Fig. 4. Indeed, the three terms in the square brackets of Eq. (34) display very linear behaviors over a large energy range ($\displaystyle\lesssim 5~{}\mathrm{MeV}$). This important feature is independent from our first guess (12) about the proton-proton wave function. Therefore, the function $\displaystyle L$ is appropriate to model fitting over a large energy range.

Furthermore, the result (34) allows us to calculate the zero-energy values of $\displaystyle L$ and $\displaystyle\Lambda$. Using $\displaystyle A=0.8850(5)~{}\mathrm{fm}^{-1/2}$, $\displaystyle\alpha_{0}=-7.815(9)~{}\mathrm{fm}$, and $\displaystyle C(0)=0.260(1)$ computed in potential models, we find

[TABLE]

It should be noted that Bahcall and May originally obtained remarkably good estimates of $\displaystyle A$ and $\displaystyle C(0)$ from effective-range theory Bahcall and May (1969)

[TABLE]

In Eq. (36), we have used the numerical values of the proton-proton ${}^{1}S_{0}$ effective range $\displaystyle r_{\rm pp}=2.77(1)~{}\mathrm{fm}$ and the proton-neutron ${}^{3}S_{1}$ effective range $\displaystyle r_{\rm pn}=1.75(1)~{}\mathrm{fm}$ from Table XIV in Ref. Machleidt (2001).

If, on the other hand, we perform the linear fitting of $\displaystyle L$ on $\displaystyle[300,600]~{}\mathrm{keV}$ directly from Eqs. (2) and (33), then we get

[TABLE]

The uncertainties are due to the differences between the potential models. The function $\displaystyle L$ of Eq. (37) is compared to the results from potential models in Fig. 4. The actual curve of $\displaystyle L$ deviates from a straight line because of the short-range behavior of the wave functions. Despite this deviation, it turns out that $\displaystyle L$ is accurate by less than $\displaystyle 2\%$ error below $\displaystyle 5~{}\mathrm{MeV}$, and is especially good below $\displaystyle 1~{}\mathrm{MeV}$. This adequacy confirms the validity of the parametrization (33).

The matrix element $\displaystyle\Lambda$ computed in the Reid93 potential is depicted in Fig. 5. The curves obtained in other potentials, along with the result (33), are indistinguishable at this scale. These curves are quite rarely shown in the literature Rupak and Ravi (2015). At zero energy, we find the important value for stellar nucleosynthesis $\displaystyle\Lambda^{2}(0)=7.034(33)$, which is consistent with Ref. Adelberger et al. (2011). Using Eq. (7) and the numerical value of $\displaystyle F(Q)$ from Eq. (61) plus $\displaystyle 1.62\%$ to account for radiative corrections Adelberger et al. (2011); Chen et al. (2013); Marcucci et al. (2013); *Marcucci2014a; Kurylov et al. (2002); *Kurylov2003, the corresponding value of $\displaystyle S(0)$ is $\displaystyle 3.95(3)\times 10^{-23}~{}\mathrm{MeV}\,\mathrm{fm}^{2}$.

III Complex analysis of the weak capture

Remarkably, $\displaystyle\Lambda$ reaches a maximum near $\displaystyle 130~{}\mathrm{keV}$, that corresponds, through Eq. (33), to the minimum of $\displaystyle\left|\Delta^{+}_{0}\right|$ seen in Fig. 2. The actual origin of this maximum is revealed by the continuation of $\displaystyle\Lambda$ to complex energies, as provided by Eq. (33). The analytic continuation of $\displaystyle\Lambda^{2}$ based on Eqs. (33) and (37) is shown in Figs. 6 and 7. Note that the curve along the positive real semi-axis in Fig. 6 corresponds to Fig. 5. The deuteron bound state pole in Fig. 6 is due to $\displaystyle(1+\epsilon)^{-1}$ in Eq. (33).

Furthermore, in contrast to $\displaystyle\Lambda$, the function $\displaystyle L$ is holomorphic in the neighborhood of $\displaystyle E=0$, because of the properties of the modified Coulomb functions $\displaystyle\Phi_{\eta,0}$ and $\displaystyle\Psi_{\eta,0}$ in the limit $\displaystyle E\rightarrow 0$ Gaspard and Sparenberg (2018); Gaspard (2018). Therefore, it follows from Eq. (33) that any singularity of $\displaystyle\left|\Delta^{+}_{0}\right|^{-1}$ is reflected on $\displaystyle\Lambda$. In this regard, it can be shown that $\displaystyle\left|\Delta^{+}_{0}\right|^{-1}$ has two poles at $\displaystyle(-140\pm 467\,\mathrm{i})~{}\mathrm{keV}$, that are interpreted as the proton-proton ${}^{1}S_{0}$ resonance poles Kok (1980); Gaspard and Sparenberg (2018). In Fig. 6, only one of them is visible because the plot is restricted to the upper half-plane ($\displaystyle\operatorname{Im}E\geq 0$). The other one is shown in panel (a) of Fig. 8.

The function $\displaystyle\left|\Delta^{+}_{0}\right|^{-1}$ also possesses a branch cut along the negative real semi-axis due to the logarithm in the Bethe function $\displaystyle g_{0}$ in Eq. (23). This branch cut makes $\displaystyle\Lambda$ complex at negative energy, although it is real at positive energy. Being on the boundary of the plots, the branch cut cannot be seen either in Figs. 6 or 7, but only in the top views of Fig. 8.

In addition, $\displaystyle\left|\Delta^{+}_{0}\right|^{-1}$ is responsible for the accumulation of poles and zeros shown in Fig. 7. These singularities originate from the terms $\displaystyle\psi(\pm\mathrm{i}\eta)$ in the Bethe function $\displaystyle g_{0}$. The pole-zero pattern is repeated each $\displaystyle E=-\mathrm{Ry}/n^{2}$ for $\displaystyle n\in\{1,2,3,\ldots\}$. Such a structure can be understood as a set of virtual states generated by the Coulomb potential between the protons.

These considerations about the analytic properties of $\displaystyle\Lambda$ have two major consequences. First, the maximum of $\displaystyle\Lambda$ is directly related to the Coulomb potential. Especially, one sees in Fig. 6 that it results from a saddle point between the conjugated resonance poles and the low-energy Coulomb singularities.

Second, $\displaystyle\Lambda$ is not analytic at $\displaystyle E=0$ due to $\displaystyle\Delta^{+}_{0}$ in Eq. (33). Therefore, its series expansion is not expected to converge over a nonzero energy range around $\displaystyle E=0$. The common way of extracting the derivatives of $\displaystyle\Lambda$ would be to use polynomial extrapolation from data on a finite energy interval Schiavilla et al. (1998); Marcucci et al. (2013); *Marcucci2014a; Acharya et al. (2016). However, such a method is not accurate due to the non-negligible influence of the interval itself Acharya et al. (2016); Marcucci et al. (2013); *Marcucci2014a.

IV Zero-energy derivatives

One way to address the issue of the extrapolation to $\displaystyle E=0$ is to expand $\displaystyle\Lambda$ in power series directly from Eq. (33) taking advantage of the flatness of $\displaystyle L$ at low energy. The result (33) provides a constraint on the derivatives of $\displaystyle\Lambda$ with respect to the energy, especially the first derivative at zero energy: $\displaystyle\Lambda^{\prime}(0)$. This value plays a significant role in the proton-proton fusion at solar energies Adelberger et al. (2011). From Eq. (33), the logarithmic derivative of $\displaystyle\Lambda$ can be easily calculated

[TABLE]

where the function $\displaystyle\left|\Delta^{+}_{0}\right|$ has been replaced by $\displaystyle\Delta_{0}$ because they share the same asymptotic expansion at $\displaystyle E=0$. This is due to the fact that the term $\displaystyle(\mathop{}\!\mathrm{e}^{2\pi\eta}-1)^{-2}$ in Eq. (24) is negligible since all of its derivatives are zero. The first terms in the asymptotic expansion of $\displaystyle\Delta_{0}$ are

[TABLE]

However, it should be noted that this expansion does not converge at $\displaystyle E=0$ because of the Coulomb singularities in $\displaystyle\Delta_{0}$. It remains nevertheless valid for $\displaystyle E\ll\mathrm{Ry}$ Gaspard and Sparenberg (2018); Baye and Brainis (2000). In addition, the expansion of the regular Coulomb integral is given by

[TABLE]

Combining these results in Eq. (38), we find that the logarithmic derivative of $\displaystyle\Lambda$ at $\displaystyle E=0$ mostly depends on the effective-range parameters

[TABLE]

The prime over $\displaystyle\Lambda$ refers to the derivative with respect to $\displaystyle E$. This novel result was not obtained by Bahcall and May, although its numerical value is given in their paper Bahcall and May (1969). From $\displaystyle L_{0}$ and $\displaystyle L_{1}$ in Eq. (37), the result (41) yields

[TABLE]

This result and all the following ones are expressed in units of MeV and fm.

It should be noted that, according to Eq. (34), $\displaystyle L_{0}$ and $\displaystyle L_{1}$ also depend on effective-range parameters. In this regard, the expansion (42) is incomplete. However, it is not possible to extract the full dependence of $\displaystyle L$ in the effective-range parameters since it would be necessary to modify the potential models accordingly.

Besides, one notices in Eq. (42) that the uncertainty of the effective range $\displaystyle r_{0}$ only slightly influences the result. Therefore, it is useful to re-express Eq. (42) in the neighborhood of $\displaystyle\left|\alpha_{0}\right|=7.815~{}\mathrm{fm}$ with $\displaystyle r_{0}=2.77(1)~{}\mathrm{fm}$. We obtain the following result

[TABLE]

around the scattering length $\displaystyle\alpha_{0}=-7.815(9)~{}\mathrm{fm}$. The central value

[TABLE]

is compatible with Refs. Bahcall and May (1969); Chen et al. (2013). It turns out that the term $\displaystyle L_{1}/L_{0}$ in Eq. (41), which contains the short-range behavior of the wave functions, only contributes to about $\displaystyle 2.6\%$. Therefore, the uncertainty of $\displaystyle L_{1}/L_{0}$ marginally affects the overall uncertainty of $\displaystyle\Lambda^{\prime}(0)/\Lambda(0)$, which is primarily due to the effective-range parameters $\displaystyle\alpha_{0}$ and $\displaystyle r_{0}$.

It is worth noting that Eq. (33) now determines all the derivatives of $\displaystyle\Lambda$ at $\displaystyle E=0$. Indeed, the higher-order derivatives of $\displaystyle L$ are negligible in the expansion of $\displaystyle\Lambda$ compared to the other terms. The reason is that the derivatives of $\displaystyle\Lambda$ are dominated by $\displaystyle\Delta^{+}_{0}$. In particular, the second derivative of $\displaystyle\Lambda$ can also be calculated analytically from Eq. (33). We have

[TABLE]

where the primes refer to derivatives with respect to the main variable: either $\displaystyle E$ for $\displaystyle\Delta_{0}(E)$, or $\displaystyle\epsilon$ for $\displaystyle\phi_{1,0}(\epsilon)$. All the functions in Eq. (45) are implicitly evaluated at $\displaystyle E=0$. It turns out that the second derivative $\displaystyle L_{2}$ can be neglected, as it contributes to only $\displaystyle 0.01\%$, far below the uncertainty of the other terms. The terms of Eq. (45) containing $\displaystyle\Delta_{0}$ and its derivatives are dominating the other ones, especially the term $\displaystyle-\Delta_{0}^{\prime\prime}/\Delta_{0}$ which is about $\displaystyle 29~{}\mathrm{MeV}^{-2}$. Inserting the expansion of $\displaystyle\phi_{1,0}$ and $\displaystyle\Delta_{0}$ at $\displaystyle E=0$ into Eq. (45), but without replacing the effective-range parameters $\displaystyle\alpha_{0}$ and $\displaystyle r_{0}$ by their numerical value for now, we find

[TABLE]

The last term in Eq. (46) comes from the term $\displaystyle(\Delta_{0}^{\prime}/\Delta_{0})^{2}$ in Eq. (45). The uncertainties in this term are negligible, as it solely depends on accurately known physical quantities ($\displaystyle\hbar$, $\displaystyle\alpha$, and $\displaystyle m_{\rm p}$).

As previously, if we focus on the neighborhood of $\displaystyle\left|\alpha_{0}\right|=7.815~{}\mathrm{fm}$ assuming $\displaystyle r_{0}=2.77(1)~{}\mathrm{fm}$, we get from Eq. (46) the expression

[TABLE]

Note that the remainder term in Eq. (47) is $\displaystyle 0.0785(1)\,(\left|\alpha_{0}\right|-7.815)^{2}$ in unit $\displaystyle\mathrm{MeV}^{-2}$. The central value

[TABLE]

is in accordance with Ref. Chen et al. (2013). The uncertainty on $\displaystyle\Lambda^{\prime\prime}(0)/\Lambda(0)$ is mainly due to $\displaystyle\alpha_{0}=-7.815(9)~{}\mathrm{fm}$. This result (48) is considerably more accurate than what we get from direct fitting on $\displaystyle\Lambda$. In fact, the direct computation of the second derivative of $\displaystyle\Lambda$ depends too much on the energy interval chosen for the fitting, hence degrading its accuracy. Similar effective-range constraints can be derived for higher-order derivatives of $\displaystyle\Lambda$ from Eq. (33). The expected accuracy of this approach does not exceed about $\displaystyle 0.2\%$ as it is limited by the uncertainty on $\displaystyle\alpha_{0}$.

Finally, we deduce the zero-energy derivatives of the astrophysical $\displaystyle S$ factor from Eqs. (43) and (47). Knowing from Eq. (7) that $\displaystyle S(E)$ is proportional to $\displaystyle F(E+Q)\left|\Lambda(E)\right|^{2}$, the logarithmic derivatives read Baye (2013)

[TABLE]

and

[TABLE]

Using the values of the Fermi phase-space integral from Eq. (61), we get the results

[TABLE]

The central values $\displaystyle 11.25(1)~{}\mathrm{MeV}^{-1}$ and $\displaystyle 169.5(3)~{}\mathrm{MeV}^{-2}$, obtained by setting $\displaystyle\left|\alpha_{0}\right|=7.815(9)~{}\mathrm{fm}$, are compatible with Refs. Bahcall and May (1969); Adelberger et al. (2011); Chen et al. (2013); Marcucci et al. (2013); *Marcucci2014a; Angulo et al. (1999); Baye (2013). These results are obtained with an unprecedented high accuracy. In the literature, most of the uncertainties are due to the polynomial extrapolation of $\displaystyle S(E)$ which is highly sensitive to the chosen energy interval Marcucci et al. (2013); *Marcucci2014a; Acharya et al. (2016). Our method is based instead on the fitting of $\displaystyle L$, as suggested by the analytic structure of $\displaystyle\Lambda$ at low energy. Consequently, the results (51) are not affected by the uncertainty of $\displaystyle 1\%$ reported for $\displaystyle S(0)$ Adelberger et al. (2011).

V Conclusion

To conclude, we have derived an accurate parametrization of the energy dependence of the weak capture matrix element $\displaystyle\Lambda$ valid up to a few MeVs, that is based on recent effective-range functions Ramírez Suárez and Sparenberg (2017); Gaspard and Sparenberg (2018). This result provides the analytic continuation of $\displaystyle\Lambda$ to complex energies, and highlights the relationship between its maximum near $\displaystyle 0.13~{}\mathrm{MeV}$, the broad proton-proton resonance, and the Coulomb sub-threshold singularities. In addition, it leads to a remarkably accurate determination of the logarithmic derivatives of the astrophysical $\displaystyle S$ factor at $\displaystyle E=0$ in terms of effective-range parameters. Our method bypasses the issue Acharya et al. (2016) of the energy-range dependence in the polynomial fitting of $\displaystyle S$ by means of the function $\displaystyle L$, that is analytic at low energy, in contrast to $\displaystyle S$. In this regard, the gain in accuracy on $\displaystyle S(0)$, $\displaystyle S^{\prime}(0)$, and $\displaystyle S^{\prime\prime}(0)$ using our method is expected to be similar if corrections, such as the two-body current terms Schiavilla et al. (1998); Park et al. (2003); Marcucci et al. (2013); *Marcucci2014a, are taken into account. Finally, the new parametrization (33) is appropriate for use in stellar and Big-Bang astrophysics as it covers a large energy range up to the binding energy of the deuteron.

Acknowledgements.

This work was supported by the European Union’s Horizon 2020 (Excellent Science) research and innovation program under Grant Agreement No. 654002.

Appendix A Fermi phase-space integral

When calculating the proton-proton weak capture cross section, we are led to integrate the Dirac delta of energy-momentum conservation over the momenta of the three outgoing particles: the deuteron, the positron, and the electronic neutrino. The resulting integral is known as the Fermi phase-space integral and reads in first approximation Bahcall (1966); Wilkinson (1982); Chen et al. (2013)

[TABLE]

as long as the recoil of the deuteron is neglected. The released energy $\displaystyle Q=2m_{\rm p}c^{2}-m_{\rm d}c^{2}-m_{\rm e}c^{2}$ is found to be $\displaystyle 0.420236(17)~{}\mathrm{MeV}$ with the masses from Ref. Tanabashi et al. (2018). The variable $\displaystyle w$ in Eq. (52) is the positron energy divided by its mass. With this notation, $\displaystyle\sqrt{w^{2}-1}$ is to be understood as the positron momentum divided by its mass. From energy conservation, the upper bound denoted as $\displaystyle\bar{w}$ is equal to $\displaystyle(E+Q+m_{\rm e}c^{2})/(m_{\rm e}c^{2})$. It means that the Fermi integral $\displaystyle F$ also depends on the proton-proton energy $\displaystyle E$. The purpose of this Appendix is to calculate the low-energy dependence of $\displaystyle F$ on $\displaystyle E$.

The Coulomb factor $\displaystyle P$ in Eq. (52), accounting for the distortion of the positron wave function in the electric field of the deuteron, is given by Bahcall (1966); Wilkinson (1982); Chen et al. (2013)

[TABLE]

where $\displaystyle\nu$ is equal to $\displaystyle\sqrt{1-\alpha^{2}}$ with the fine-structure constant $\displaystyle\alpha\simeq 1/137.036$, and $\displaystyle\rho=Rm_{\rm e}c^{2}/(\hbar c)$ is the dimensionless radius of the deuteron. In the following calculations, we will assume $\displaystyle R=2.14~{}\mathrm{fm}$ Chen et al. (2013). In Eq. (53), the Sommerfeld parameter of the emitted positron $\displaystyle\eta_{\rm e}=\alpha w/\sqrt{w^{2}-1}$ must be positive, as it is repelled by the nucleus. Conversely, in a $\displaystyle\beta^{-}$ decay, the Sommerfeld parameter $\displaystyle\eta_{\rm e}$ should take a minus sign. It should be noted that in the nonrelativistic limit ($\displaystyle\nu\rightarrow 1$), the Coulomb distortion factor $\displaystyle P(w)$ becomes

[TABLE]

The Fermi integral (52) cannot be analytically calculated in a simple form. However, very efficient approximations exist. One way is to expand the Coulomb factor (53) in series of the fine-structure constant $\displaystyle\alpha$. We find

[TABLE]

where $\displaystyle\gamma=0.5772\ldots$ is the Euler-Mascheroni constant. In this work, we limit ourselves to the order $\displaystyle\alpha^{2}$, as it is enough to obtain at least five decimal places in the final results. The same approach is followed in Ref. Wilkinson (1982) up to $\displaystyle\alpha^{3}$. Now, we just have to calculate one Fermi integral for each term in the expansion (55). The advantage is that the integrals of the form

[TABLE]

which will come into play, can be expressed in terms of elementary functions for $\displaystyle p\in\mathbb{Z}$. Such expressions can be obtained by expanding the last factor $\displaystyle(\bar{w}-w)^{2}$ in Eq. (56). The results read for $\displaystyle p=0$

[TABLE]

for $\displaystyle p=1$

[TABLE]

and for $\displaystyle p=2$

[TABLE]

We notice that, according to Eqs. (57), (58), and (59), the Fermi integral is expected to behave as $\displaystyle\operatorname{\mathcal{O}}(E^{5})$ at relatively large energy ($\displaystyle E\gg m_{\rm e}c^{2}$). Therefore, $\displaystyle F(E+Q)$ dominates the $\displaystyle\operatorname{\mathcal{O}}(E^{-2})$ behavior of $\displaystyle\Lambda^{2}$ in Eq. (7).

The factor $\displaystyle\sqrt{w^{2}-1}$ in the logarithmic term of expansion (55) can be neglected because it remains of the order of $\displaystyle 1$ except at large proton-proton energies ($\displaystyle E\gg m_{\rm e}c^{2}$). Therefore, using Eq. (56), the Fermi integral (52) is approximated by

[TABLE]

This expression allows us to find at least five decimal places without requiring numerical integration. Another advantage is the computation of the derivatives of the Fermi integrals with respect to $\displaystyle E$. In this work, we need the first two derivatives of $\displaystyle F(E+Q)$ at zero proton energy ($\displaystyle E=0$). This can be easily achieved with the derivatives of $\displaystyle f_{p}(\bar{w})$ with respect to $\displaystyle\bar{w}$ that are obtained directly from Eqs. (57), (58), and (59). The derivatives of $\displaystyle F(E+Q)$ with respect to $\displaystyle E$ have thus essentially the same expressions as Eq. (60) by replacing $\displaystyle f_{p}(\bar{w})$ with the derivatives with respect to $\displaystyle\bar{w}$, denoted as $\displaystyle f_{p}^{(n)}(\bar{w})$. Note the change of variable $\displaystyle\partial_{E}^{n}F=(m_{\rm e}c^{2})^{-n}\partial_{\bar{w}}^{n}F$ in the manipulation.

The numerical values of the functions $\displaystyle f_{p}(\bar{w})$ and their derivatives at $\displaystyle E=0$, that is for $\displaystyle\bar{w}=(Q+m_{\rm e}c^{2})/(m_{\rm e}c^{2})=1.82238(3)$, are given in Table 2. Inserting the numerical values of Table 2 in the approximation (60) for the different derivative orders ($\displaystyle n=0,1,2$) leads to the results

[TABLE]

These results have also been checked by numerical integration in Wolfram Mathematica Wolfram (1999). The uncertainties in Table 2 and Eq. (61) come from the released energy $\displaystyle Q$. Finally, the low-energy behavior of the Fermi integral can be written as

[TABLE]

with the numerical values of Eq. (61).

It should be noted that the third derivative of the Fermi function (52) with respect to $\displaystyle E$ is devoid of integral and can be expressed exactly in terms of $\displaystyle P(\bar{w})$. We have

[TABLE]

from which the numerical value $\displaystyle F^{(3)}(Q)=40.498(2)~{}\mathrm{MeV}^{-3}$ at $\displaystyle E=0$ is easily found. Our approach avoids using numerical derivatives, as they are ill-conditioned in finite precision arithmetic, especially for high-order derivatives. This also ensures the accuracy of the results (61).

Appendix B Laplace transforms of the modified Coulomb functions

In this Appendix, we present the derivation of the Laplace transforms of the modified Coulomb wave functions $\displaystyle\Phi_{\eta,\ell}(\rho)$ and $\displaystyle\Psi_{\eta,\ell}(\rho)$ defined in Ref. Gaspard (2018). More explicitly, we are looking for analytical expressions of the integrals

[TABLE]

and

[TABLE]

where $\displaystyle x=br$ is the dimensionless radial coordinate, and $\displaystyle\sqrt{\epsilon}=k/b$ is the dimensionless wave number. Although we only need the result for $\displaystyle\ell=0$, we have made our derivation more general. The reason is that we use the connection formula between $\displaystyle\Phi_{\eta,\ell}(\rho)$ and $\displaystyle\Psi_{\eta,\ell}(\rho)$ developed in Ref. Gaspard (2018) to calculate $\displaystyle\psi_{\beta,\ell}(\epsilon)$ on the basis of $\displaystyle\phi_{\beta,\ell}(\epsilon)$ for any integer $\displaystyle\ell$.

B.1 Regular Coulomb integral

As shown in Ref. Gaspard (2018), the Coulomb function $\displaystyle\Phi_{\eta,\ell}$ in Eq. (64) is given by

[TABLE]

where $\displaystyle\kappa=\sqrt{\epsilon}$, and $\displaystyle\chi=2/(a_{\mathrm{B}}b)$ is the Bahcall and May constant Bahcall and May (1969). The regularized confluent hypergeometric function in Eq. (66) is defined by the series Olver et al. (2010)

[TABLE]

where $\displaystyle(a)_{n}=\Gamma(a+n)/\Gamma(a)$ is the Pochhammer symbol. The division by $\displaystyle\Gamma(b)$ in Eq. (67) eliminates the singularities of $\displaystyle{}_{1}F_{1}\!\left(\begin{subarray}{c}a\\ b\end{subarray};z\right)$ at $\displaystyle b\in\mathbb{Z}_{\leq 0}$ Olver et al. (2010). Using the definition (66), the regular Coulomb integral (64) expands as follows

[TABLE]

All the remaining integrals in Eq. (68) are given by

[TABLE]

One notices that the combination of Eqs. (68) and (69) leads to the Gauss hypergeometric function $\displaystyle{}_{2}F_{1}$, or more specifically to its regularized version Olver et al. (2010)

[TABLE]

Using the definition (70) in Eq. (68), we obtain the following result:

[TABLE]

Remarkably, this result is considerably simplified in the special case $\displaystyle\ell=0$. Indeed, the hypergeometric function in Eq. (71) is then of the form $\displaystyle{}_{2}F_{1}\!\left(\begin{subarray}{c}a,b\\ a\end{subarray};z\right)$, which reduces to $\displaystyle(1-z)^{-b}$ Olver et al. (2010) because of the simplification in the series (70). From Eq. (71), one finds

[TABLE]

This useful result is at the basis of the parametrization of $\displaystyle\Lambda(E)$ proposed in this paper.

B.2 Irregular Coulomb integral

Now, we calculate the Laplace transform (65) of $\displaystyle\Psi_{\eta,\ell}$. This calculation is significantly less straightforward than for $\displaystyle\Phi_{\eta,\ell}$, because it does not reduce to elementary functions for $\displaystyle\ell=0$. The Coulomb function $\displaystyle\Psi_{\eta,\ell}$ is defined in Ref. Gaspard (2018) as

[TABLE]

where the choice of the upper or lower sign is immaterial. In Eq. (73), $\displaystyle U\!\left(\begin{subarray}{c}a\\ b\end{subarray};z\right)$ is the Tricomi confluent hypergeometric function, and the Bethe functions $\displaystyle h^{\pm}_{\eta\ell}$ are defined as

[TABLE]

The subtraction by $\displaystyle w_{\eta\ell}h^{\pm}_{\eta\ell}\Phi_{\eta,\ell}$ in Eq. (73) is intended to compensate for the singularities of $\displaystyle U\!\left(\begin{subarray}{c}a\\ b\end{subarray};z\right)$ in the complex plane of the energy. This operation makes $\displaystyle\Psi_{\eta,\ell}(kr)$ regular for $\displaystyle k\in\mathbb{C}$ Gaspard (2018).

Performing the direct integration of Eq. (73) by means of the integral representation of $\displaystyle U\!\left(\begin{subarray}{c}a\\ b\end{subarray};z\right)$ Olver et al. (2010) leads to

[TABLE]

Note that, this function is not finite for partial waves higher than $\displaystyle S$ ($\displaystyle\ell>0$) due to the vertical asymptote of $\displaystyle\Psi_{\eta,\ell}(kr)$ at $\displaystyle r=0$. When $\displaystyle\ell=0$, the hypergeometric function in the above equation can be efficiently computed from its continued fraction expansion.

The expression (75) is quite difficult to analyze at low energy because the hypergeometric function shows an essential singularity at $\displaystyle\epsilon=0$. Although this singularity is compensated by $\displaystyle w_{\eta\ell}h^{\pm}_{\eta\ell}\phi_{\beta,\ell}$, it prevents the hypergeometric function from having a convergent low-energy expansion. This is why we propose to determine a suitable approximation to Eq. (75) from another approach.

It turns out that the function $\displaystyle\Psi_{\eta,\ell}$ is related to the regular Coulomb function $\displaystyle\Phi_{\eta,\ell}$. We have shown in Ref. Gaspard (2018) that $\displaystyle\Psi_{\eta,\ell}$ obeys the following connection formula

[TABLE]

where the dots refer to derivatives with respect to $\displaystyle\ell$. This useful property is preserved by the Laplace transforms (64) and (65). Therefore, the function $\displaystyle\psi_{\beta,\ell}(\epsilon)$ can be calculated from derivatives of $\displaystyle\phi_{\beta,\ell}(\epsilon)$ as follows Gaspard (2018)

[TABLE]

However, in order to calculate the derivatives in Eq. (77), we need to use the general expression (71) valid of $\displaystyle\phi_{\beta,\ell}$ for all $\displaystyle\ell\in\mathbb{C}$. In this regard, we have found convenient to approximate the hypergeometric function by its low-energy confluent limit

[TABLE]

The approximation (78) could be improved at $\displaystyle\kappa=0$ by the confluence expansion (20a) in Ref. Nagel (2004). However, this expansion converges so slowly for $\displaystyle\kappa>1$ that we will not use it here. The advantage of the approximation (78) is the consistency with Bahcall’s and May’s results in the zero-energy limit. It is useful in our calculation to rewrite the confluent hypergeometric function in Eq. (78) in terms of the regular Coulomb function

[TABLE]

The approximation of $\displaystyle\phi_{\beta,\ell}$ for $\displaystyle\ell\ll 1$ is thus given by

[TABLE]

From Eq. (71) to Eq. (80), we have taken the modulus of the factor $\displaystyle(\beta-\mathrm{i}\kappa)^{-\ell-2}$ because $\displaystyle\phi_{\beta,\ell}$ should still remain a positive real function after the approximation (78). The logarithmic derivative of $\displaystyle\phi_{\beta,\ell}$ with respect to $\displaystyle\ell$ can be easily calculated from Eq. (80):

[TABLE]

We neglect the logarithmic term in Eq. (81) because it is irrelevant in the $\displaystyle\operatorname{\mathcal{O}}(\kappa^{2})$ approximation of Eq. (78). Combining two expressions (81) evaluated at $\displaystyle\ell$ and $\displaystyle-\ell-1$ in Eq. (77), we get

[TABLE]

This relation can be simplified further by means of the definition (73) of $\displaystyle\Psi_{\eta,\ell}$ with the plus sign. Finally, after the elimination of the digamma functions with Eq. (74), we obtain

[TABLE]

In the special case of interest $\displaystyle\ell=0$, this result can be written as

[TABLE]

where $\displaystyle\Gamma(a,z)$ is the upper incomplete gamma function defined by

[TABLE]

When $\displaystyle\beta=1$ and $\displaystyle\chi=0.149816$, the ratio (84) evaluates to about $\displaystyle 4.28065$. The incomplete gamma function in Eq. (84) can also be related to the exponential integral $\displaystyle E_{1}(z)=\Gamma(0,z)$ as done in Bahcall’s and May’s work Bahcall and May (1969):

[TABLE]

Bahcall, however, limited his calculation to zero energy, in contrast to the property (84) valid up to a few MeVs.

Furthermore, the novel result (84) means that $\displaystyle\psi_{\beta,0}$ is nearly proportional to $\displaystyle\phi_{\beta,0}$ on a large energy range. The accuracy of this property is graphically tested in Fig. 9. As we can see, the relative error of the estimate at $\displaystyle\epsilon=1$, that corresponds to $\displaystyle E=B$, is only $\displaystyle 2.9\%$. The overall accuracy of the approximation (84) over a few MeVs is primarily due to the smallness of $\displaystyle\chi$ with respect to $\displaystyle 1$ ($\displaystyle\chi=0.149816$). In fact, it can be shown that both $\displaystyle\psi_{\beta,0}$ and $\displaystyle\Gamma(-1,\chi/\beta)\phi_{\beta,0}$ have the same neutral-charge limit:

[TABLE]

Therefore, the property (84) tends to be exact for $\displaystyle\chi\rightarrow 0$, but also for $\displaystyle\beta\rightarrow\infty$. These observations have important consequences in the parametrization of the weak capture matrix element $\displaystyle\Lambda$.