Bayesian analysis of capture reactions $^3\mathrm{He}(\alpha,\gamma)^7\mathrm{Be}$ and $^3\mathrm{H}(\alpha,\gamma)^7\mathrm{Li}$

TL;DR

This paper uses Bayesian methods to analyze key nuclear capture reactions, comparing effective field theory models to determine the dominant contributions and providing insights into reaction cross sections relevant for astrophysics.

Contribution

It introduces a Bayesian framework for model comparison of two EFT power countings in nuclear capture reactions, highlighting the role of two-body currents at different orders.

Findings

Two-body currents are favored at leading order for $^3$He($ o$)$^7$Be with elastic scattering data.

Both models are equally favored without elastic scattering constraints.

The model with leading order two-body currents is favored for $^3$H($ o$)$^7$Li regardless of elastic data.

Abstract

Bayesian analysis of the radiative capture reactions $^3\mathrm{He}(\alpha,\gamma)^7\mathrm{Be}$ and $^3\mathrm{H}(\alpha,\gamma)^7\mathrm{Li}$ are performed to draw inferences about the cross sections at threshold. We do a model comparison of two competing effective field theory power countings for the capture reactions. The two power countings differ in the contribution of two-body electromagnetic currents. In one power counting, two-body currents contribute at leading order, and in the other they contribute at higher orders. The former is favored for $^3\mathrm{He}(\alpha,\gamma)^7\mathrm{Be}$ if elastic scattering data in the incoming channel is considered in the analysis. Without constraints from elastic scattering data, both the power countings are equally favored. For $^3\mathrm{H}(\alpha,\gamma)^7\mathrm{Li}$, the first power counting with two-body current contributions at leading order is favored with or without constraints from elastic scattering data.

Tables

Table 1. Table 1: He 3 ​ ( α , γ ) 7 ​ Be superscript He 3 superscript 𝛼 𝛾 7 Be {}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be} : S 34 subscript 𝑆 34 S_{34} and S 34 ′ superscript subscript 𝑆 34 ′ S_{34}^{\prime} at threshold (defined as E ⋆ = 60 × 10 − 3 subscript 𝐸 ⋆ 60 superscript 10 3 E_{\star}=60\times 10^{-3} keV). The second set of errors are estimated from the EFT perturbation as detailed in the text.

Fit	$S_{34}\left(E_{\star }\right)$ (keV b)	$S_{34}^{\prime }\left(E_{\star }\right)$ (${10}^{-4}$ b)
\csvreader[head to column names, late after line=
${\text{Sp}}_{-\text{Spm}}^{+\text{Spp}}\pm \text{EFTB}$

Table 2. Table 2: He 3 ​ ( α , γ ) 7 ​ Be superscript He 3 superscript 𝛼 𝛾 7 Be {}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be} : Branching ratio R 0 subscript 𝑅 0 R_{0} at threshold (defined as E ⋆ = 60 × 10 − 3 subscript 𝐸 ⋆ 60 superscript 10 3 E_{\star}=60\times 10^{-3} keV). The second set of errors are estimated from the EFT perturbation as detailed in the text.

Fit	$R_0$
\csvreader[head to column names, late after line=
\Ratio${}_{-\text{Down}}{}^{+\text{Up}}\pm \text{EFT}$

Table 3. Table 3: H 3 ​ ( α , γ ) 7 ​ Li superscript H 3 superscript 𝛼 𝛾 7 Li {}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li} : S 34 subscript 𝑆 34 S_{34} and S 34 ′ superscript subscript 𝑆 34 ′ S_{34}^{\prime} at threshold (defined as E ⋆ = 60 × 10 − 3 subscript 𝐸 ⋆ 60 superscript 10 3 E_{\star}=60\times 10^{-3} keV). The second set of errors are estimated from the EFT perturbation as detailed in the text.

Fit	$S_{34}\left(E_{\star }\right)$ (keV b)	$S_{34}^{\prime }\left(E_{\star }\right)$ (${10}^{-4}$ b)
\csvreader[head to column names, late after line=
${\text{Sp}}_{-\text{Spm}}^{+\text{Spp}}\pm \text{EFTB}$

Table 4. Table 4: H 3 ​ ( α , γ ) 7 ​ Li superscript H 3 superscript 𝛼 𝛾 7 Li {}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li} : Branching ratio R 0 subscript 𝑅 0 R_{0} at threshold (defined as E ⋆ = 60 × 10 − 3 subscript 𝐸 ⋆ 60 superscript 10 3 E_{\star}=60\times 10^{-3} keV). The second set of errors are estimated from the EFT perturbation as detailed in the text.

Fit	$R_0$
\csvreader[head to column names, late after line=
\Ratio${}_{-\text{Down}}{}^{+\text{Up}}\pm \text{EFT}$

Equations

$$\begin{gathered} \sigma (p)=\frac{1}{16\pi M^2}\frac{1}{2}[\frac{p^2+{\gamma }_1^2}{2\mu p}|M^{{(}^2{P}_{1/2})}{|}^2 \\ +\frac{p^2+{\gamma }_0^2}{2\mu p}|M^{{(}^2{P}_{3/2})}{|}^2], \end{gathered}$$

$$\begin{gathered} |M^{(\zeta )}{|}^2=(2j+1){\left(\frac{Z_{\varphi }{M}_{\psi }}{M}-\frac{Z_{\psi }{M}_{\varphi }}{M}\right)}^2 \\ \times \left[|\mathcal{A}(p){|}^2+2|Y(p){|}^2\right]\frac{64\pi {\alpha }_e{M}^2([h^{(\zeta )}{]}^2{\mathcal{Z}}^{(\zeta )})}{\mu }, \end{gathered}$$

$$C_{1,\zeta }^2=\frac{[{\gamma }^{(\zeta )}{]}^2\Gamma (2+k_C/{\gamma }^{(\zeta )})}{\pi }[h^{(\zeta )}{]}^2{\mathcal{Z}}^{(\zeta )}.$$

$${\chi }^2$$

$$P(\mathbf{\theta }|D,H)=\frac{P(D|\mathbf{\theta },H)P(\mathbf{\theta }|H)}{P(D|H)},$$

$$P(D|H)=\int \prod _{i}d{\theta }_{i}P(D|\mathbf{\theta },H)P(\mathbf{\theta }|H),$$

$$P({\theta }_1|D,H)=\int \prod _{i\ne 1}d{\theta }_{i}P(\mathbf{\theta }|D,H),$$

$$\frac{P(M_A|D,H)}{P(M_B|D,H)}=\frac{P(D|M_A,H)}{P(D|M_B,H)}\times \frac{P(M_A|H)}{P(M_B|H)}.$$

$$P(D|M_A,H)$$

$$P(D|M_B,H)$$

$$a_0\sim U(1\mathrm{fm},70\mathrm{fm}),$$

$$r_0\sim U(-5\mathrm{fm},5\mathrm{fm}),$$

$$s_0\sim U(-30\mathrm{f}{\mathrm{m}}^3,30\mathrm{f}{\mathrm{m}}^3),$$

$${\rho }_1^{(+)}\sim U(-300\mathrm{MeV},-48\mathrm{MeV}),$$

$${\rho }_1^{(-)}\sim U(-300\mathrm{MeV},-33\mathrm{MeV}),$$

$${\sigma }_1^{(\pm )}\sim U(-5\mathrm{fm},5\mathrm{fm}),$$

$$L_1^{(\pm )}\sim U(-10,10),$$

$$K\sim U(0^{\circ },10^{\circ }).$$

$$a_0\sim U(1\mathrm{fm},70\mathrm{fm}),$$

$$r_0\sim U(-5\mathrm{fm},5\mathrm{fm}),$$

$$s_0\sim U(-45\mathrm{f}{\mathrm{m}}^3,45\mathrm{f}{\mathrm{m}}^3),$$

$${\rho }_1^{(+)}\sim U(-320\mathrm{MeV},-126\mathrm{MeV}),$$

$${\rho }_1^{(-)}\sim U(-320\mathrm{MeV},-106\mathrm{MeV}),$$

$$L_1^{(\pm )}\sim U(-10,10),$$

$$K\sim U(0^{\circ },10^{\circ }).$$

$$\frac{|\mathcal{A}(p)|}{C_0({\eta }_p)}=\left|X(p)-\frac{2\pi }{{\mu }^2}\frac{B(p)+\mu J_0(p)+{\mu }^2{k}_0{L}_{E1}^{(\zeta )}}{[C_0({\eta }_p){]}^{2}p(\cot {\delta }_0-i)}\right|,$$

$$\begin{gathered} X(p)=1+\frac{2\gamma }{3}\frac{\Gamma (2+k_C/\gamma )}{C_0({\eta }_p)} \\ \times {\int }_0^{\infty }drr{W}_{-\frac{k_C}{\gamma },\frac{3}{2}}(2\gamma r){\partial }_r\left[\frac{F_0({\eta }_p)}{pr}\right], \end{gathered}$$

$$F_l({\eta }_p,\rho )$$

$$C_l({\eta }_p)$$

$$[C_0({\eta }_p){]}^{2}p(\cot {\delta }_0-i)=-\frac{1}{a_0}+\frac{r_0}{2}{p}^2+\frac{s_0}{4}{p}^4+\cdots$$

$$-2k_{C}H({\eta }_p),$$

$$H(x)=\psi (ix)+\frac{1}{2ix}-\ln (ix).$$

$$\begin{gathered} B(p)+\mu J_0(p)=\frac{{\mu }^2}{3\pi }\frac{i{p}^3-{\gamma }^3}{p^2+{\gamma }^2}+k_{C}C(p)+\Delta B(p) \\ -\frac{k_C{\mu }^2}{2\pi }\left[2H({\eta }_p)+2{\gamma }_E-\frac{5}{3}+\ln 4\pi \right]. \end{gathered}$$

$$\begin{gathered} C(p)=\frac{{\mu }^2}{6{\pi }^2(p^2+{\gamma }^2)}{\int }_0^{1}dx{\int }_0^{1}dy\frac{1}{\sqrt{x(1-x)}\sqrt{1-y}} \\ \times (x{p}^2\ln \left[\frac{\pi }{4k_C^2}(-y{p}^2+(1-y){\gamma }^2/x-i\delta )\right] \\ +p^2\ln \left[\frac{\pi }{4k_C^2}(-y{p}^2-(1-y)p^2/x-i\delta )\right] \\ +x{\gamma }^2\ln \left[\frac{\pi }{4k_C^2}(y{\gamma }^2+(1-y){\gamma }^2/x-i\delta )\right] \\ +{\gamma }^2\ln \left[\frac{\pi }{4k_C^2}(y{\gamma }^2-(1-y)p^2/x-i\delta )\right]), \end{gathered}$$

$$B_{ab}(p)=-3\int d^{3}r\left[\frac{G_{C}^{(1)}(-B;r^{\prime},r)}{r^{\prime}}\right]\Big{|}_{r^{\prime}=0}\\ \times\frac{\partial G_{C}^{(+)}(E;\bm{r},0)}{\partial r_{a}}\frac{r_{b}}{r}\,,\\ \left[\frac{G_{C}^{(1)}(-B;r^{\prime},r)}{r^{\prime}}\right]\Big{|}_{r^{\prime}=0}=-\frac{\mu\gamma}{6\pi r}\Gamma(2+k_{C}/\gamma)W_{-\frac{k_{C}}{\gamma},\frac{3}{2}}(2\gamma r)\,,\\ G_{C}^{(+)}(E;\bm{r},0)=-\frac{\mu}{2\pi r}\Gamma(1+i\eta_{p})W_{-i\eta_{p},\frac{1}{2}}(-i2pr)\,.$$

$$\begin{gathered} Y(p)=\frac{2\gamma }{3}\Gamma (2+k_C/\gamma ){\int }_0^{\infty }drr{W}_{-k_C/\gamma ,3/2}(2\gamma r) \\ \times \left(\frac{\partial }{\partial r}+\frac{3}{r}\right)\frac{F_2({\eta }_p,rp)}{rp}. \end{gathered}$$

$$\begin{gathered} -\frac{2\pi }{h^{(\zeta )2}{\mathcal{Z}}^{(\zeta )}}=\frac{1}{p}\frac{\partial }{\partial p}9[C_1({\eta }_p){]}^2{p}^3(\cot {\delta }_1-i) \\ ={\rho }_1^{(\zeta )}-4k_{C}H\left(-i\frac{k_C}{\gamma }\right)-\frac{2k_C^2}{{\gamma }^3}(k_C^2-{\gamma }^2) \\ \left[{\psi }^{\prime }\left(\frac{k_C}{\gamma }\right)-\frac{{\gamma }^2}{2k_C^2}-\frac{\gamma }{k_C}\right], \end{gathered}$$

$$\begin{gathered} 9[C_1({\eta }_p){]}^2{p}^3(\cot {\delta }_1-i)=2k_C(k_C^2-{\gamma }^2)H(-i{\eta }_{\gamma }) \\ +\frac{1}{2}{\rho }_1(p^2+{\gamma }^2)+\frac{1}{4}{\sigma }_1(p^2+{\gamma }^2{)}^2+\cdots \\ -2k_C(k_C^2+p^2)H({\eta }_p). \end{gathered}$$

Full text

Bayesian inference, radiative capture, halo effective field theory

I Introduction

Effective field theory (EFT) calculations of low-energy capture reactions provide important uncertainty estimates in cross section evaluations. Element synthesis in the Big Bang and in stellar environments involve reactions at low energies. These reaction rates are suppressed by the Coulomb repulsion at low energy, making their measurements challenging. For example, the proton-proton fusion rate is not known experimentally at the relevant solar energy Adelberger et al. (2011). Theory is needed to extrapolate measurements at higher energies to the astrophysically relevant lower energies. EFT extrapolations are attractive as they provide model-independent error estimates. The basic idea in the EFT formulation is that one only considers certain degrees of freedom as dynamical at low momentum $p$ below a cutoff $\Lambda\gg p$. The physics associated with the high-momentum scale $\Lambda$ is not modeled but contributes through dimensionful couplings in the low-momentum theory. Observables are calculated as an expansion in $Q/\Lambda\ll 1$ called the power counting. $Q\sim p$ characterizes the low-momentum physical scale. Theoretical uncertainties are estimated from powers of $Q/\Lambda$ at a given order of the perturbative expansion.

Typical radiative capture reactions such as ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ considered in this work depend on initial state interaction, electro-weak currents and the final state wave function. Information about the initial state interaction at low momentum can be incorporated by matching the EFT couplings to elastic scattering phase shifts. The final state physics is contained in the bound state energy and the asymptotic wave function normalization which can be determined from the relevant phase shift as well. The one-body electro-weak currents are well determined, and usually they give the dominant contribution. A difficulty in these calculations is that the low momentum phase shifts are poorly known which introduces large uncertainty in the construction of the EFT interactions. Moreover, in systems that involve weakly-bound states the couplings can be fine tuned such that two-body currents contribute at leading order (LO) of the perturbation Higa et al. (2018).

The low-energy capture reactions ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ are of great interest in nuclear astrophysics. The reactions {}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$$(p,\gamma)^{8}Be  provide the high energy 8Be neutrinos detected at Super-K Fukuda et al. (2001) and SNO Aharmim et al. (2007). The neutrino flux is proportional to the ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ cross section $\sigma$ near its Gamow energy $E_{G}\sim 20$ keV as $\sigma^{0.81}$ Bahcall and Ulrich (1988); Cyburt and Davids (2008). Several measurements of the ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ capture rate have been performed in recent decades, for example  Di Leva et al. (2009); Constantini et al. (2008); Kontos et al. (2013); Brown et al. (2007); Nara Singh et al. (2004). The reaction ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ is important for calculating 7Li abundances, constraining astrophysical models of Big Bang Nucleosynthesis Kawano et al. (1991); Smith et al. (1993). There are less constraints on this reaction from recent measurements Brune et al. (1994); Bystritsky et al. (2017).

In this work we use Bayesian analysis to revisit the calculation of ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ Higa et al. (2018), and also extend the analysis to the related reaction ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$. Bayesian analysis is useful in situations as these where we have a large number of parameters that are poorly constrained by lack of accurate data especially in the elastic channels. Parameter estimation is important in not only extrapolation of reaction rates to solar energies but also to develop a systematic power counting expansion in $Q/\Lambda$. This is crucial for providing theoretical error estimates. Further, Bayesian analysis allows for a quantitative comparison of EFT formulations where, for example, the two-body currents appear at different orders in the perturbation.

The expressions for the cross sections are from Ref. Higa et al. (2018) that used halo EFT in the calculations. In halo EFT Bertulani et al. (2002); Bedaque et al. (2003), nuclear clusters are treated as point like particles at momenta that are too small to cause the clusters to break apart. Capture reaction calculations are outlined in the review article Ref. Rupak (2016). We consider two power countings depending on the size of the scattering length $a_{0}$ in the incoming $s$-wave channel. In one power counting, two-body currents contribute at LO and in the other they contribute at next-to-leading order (NLO).

II Halo EFT calculation

The calculations for ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ are similar. The discussion here follows Ref. Higa et al. (2018). The ground state of 7Be has a binding energy $B_{0}^{\mathrm{(Be)}}=1.5866$ MeV and spin-parity assignment $\frac{3}{2}^{-}$. The first excited state has a binding energy $B_{1}^{\mathrm{(Be)}}=1.1575$ MeV and spin-parity assignment $\frac{1}{2}^{-}$. These binding energies are identified with the low-energy scale Higa et al. (2018) and they are treated as a bound state of 3He and $\alpha$ clusters in the calculation. The higher excited states of 7Be with the same spin-parity as the first two states lie at least about 8 MeV above the 3He-$\alpha$ threshold. These are not included in the low-energy theory. The proton separation energy of 3He ($\sim 5.5$ MeV), excited states of $\alpha$ ($\gtrsim 20$ MeV) are identified with the high-energy scale as well. Thus, we treat the 3He and $\alpha$ clusters as point-like fundamental particles. The spin-parity assignments of $\frac{1}{2}^{+}$ for the 3He and $0^{+}$ for the $\alpha$ implies that the $\frac{3}{2}^{-}$, $\frac{1}{2}^{-}$ states of 7Be are $p$-wave bound states of 3He and $\alpha$. Then radiative capture at low energy is dominated by E1 transition from $s$-wave initial states with a small contribution from the $d$ waves Higa et al. (2018). For the ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ reaction, the ground state of 7Li has a binding energy $B_{0}^{\mathrm{(Li)}}=2.467$ MeV and spin-party $\frac{3}{2}^{-}$. The first excited state has a binding energy $B_{1}^{\mathrm{(Li)}}=1.989$ MeV and spin-parity assignment $\frac{1}{2}^{-}$. These two $p$-wave states of the $\frac{1}{2}^{+}$ 3H and $0^{+}$ $\alpha$ are included in the theory. The higher excited states with the same quantum numbers lie at least about 6 MeV above the 3H-$\alpha$ threshold. These are not included in the low-energy theory. The neutron separation energy of 3H ($\sim 6$ MeV) is identified with the high energy theory, and so 3H is treated as a point-like fundamental particle along with the $\alpha$. The main differences between the EFT expressions for the ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ cross sections are the charges of the nuclei, and the relevant masses and binding energies. The separation between the low and high energy is smaller in the 3H-$\alpha$ system which would result in a slower convergence of the $Q/\Lambda$ expansion.

The total cross section for both the capture reactions is given as Higa et al. (2018):

[TABLE]

at center-of-mass (cm) relative momentum $p$ of the incoming particles. We use the spectroscopic notation ${}^{2s+1}L_{j}$ to indicate the capture to the ${}^{2}P_{3/2}$ ground state with binding momentum $\gamma_{0}=\sqrt{2\mu B_{0}}$, and the ${}^{2}P_{1/2}$ first excited state with binding momentum $\gamma_{1}=\sqrt{2\mu B_{1}}$. $M=M_{\psi}+M_{\phi}$ is the total mass, and $\mu=M_{\psi}M_{\phi}/M$ is the reduced mass of the incoming particles. $M_{\phi}$ is the scalar $0^{+}$ $\alpha$ mass, and $M_{\psi}$ is the fermion $\frac{1}{2}^{+}$ 3He or 3H mass as appropriate. The squared amplitude is given by

[TABLE]

where $\zeta={}^{2}P_{3/2}$ for the ground state ($j=3/2$) and $\zeta={}^{2}P_{1/2}$ for the exited state ($j=1/2$). $Z_{\phi}=2$ is the charge of the $\alpha$ and $Z_{\psi}$ is the charge of either the 3He or the 3H particle. $\alpha_{e}=e^{2}/(4\pi)\approx 1/137$ is the fine structure constant. The capture from initial $s$ and $d$ wave states are given by the amplitudes $\mathcal{A}(p)$ and $Y(p)$ shown in the Appendix, Eq. (12) and Eq. (19) respectively. The factor $([h^{(\zeta)}]^{2}\mathcal{Z}^{(\zeta)})$ is associated with the normalization of the outgoing 7Be or 7Li particle wave function, shown in Appendix, Eq. (20). They are related to the asymptotic normalization constant (ANC) as:

[TABLE]

The ANCs depend on the binding momentum and the effective range $\rho_{1}$ of the corresponding $p$-wave elastic channel of the outgoing 7Be or 7Li state.

The $S$-factor is calculated as a function of the c.m. incoming kinetic energy $E=p^{2}/(2\mu)$ from the cross section as $S_{34}(E)=Ee^{2\pi\eta_{p}}\sigma(p=\sqrt{2\mu E})$ where $\eta_{p}=\alpha_{e}Z_{\psi}Z_{\phi}\mu/p$ is the Sommerfeld parameter.

III Bayesian parameter estimation and model comparison

In the EFT calculation we have several unknown couplings and parameters that we wish to constrain from available data. Suppose we represent the set of unknown couplings and parameters by the vector $\bm{\theta}$ in the multi-dimensional parameter space. In the frequentist approach, the parameters $\bm{\theta}$ are determined by minimizing the $\chi^{2}$ given by

[TABLE]

where the data set $D$ consists of the $N$ measurements $y_{i}$ with corresponding measurement errors $\sigma_{i}$. $\mu_{i}(\bm{\theta})$ are the theory predictions. Minimizing $\chi^{2}$ maximizes the likelihood function $P(D|\bm{\theta},H)\propto\exp(-\chi^{2}/2)$ which gives the conditional probability of the data set $D$ to be true given the parameters $\bm{\theta}$ and the proposition $H$. Background information, theoretical assumptions, and hypothesis are included in the proposition $H$.

In Bayesian analysis, the question one wants to answer is how likely some parameter value $\bm{\theta}$ given the data $D$. The answer is contained in the posterior probability distribution given by

[TABLE]

where $P(\bm{\theta}|H)$ is the prior distribution of the parameters which is predicated on the theoretical assumptions that is made in proposition $H$. The prior distribution is where one can, for example, include EFT assumptions about the expected sizes of the unknown couplings and parameters. The evidence $P(D|H)$ is also known as the marginal likelihood because

[TABLE]

comes from marginalizing the likelihood on the parameters. We used $\int\prod_{i}d\theta_{i}P(\bm{\theta}|D,H)=1$ above for normalized probabilities.

Parameter estimation from the posterior distribution in Eq. (5) does not require evaluation of the overall normalization factor, the evidence $P(D|H)$. Markov chain Monte Carlo (MCMC) algorithms such as Metropolis-Hastings or some variation of it can be used to sample points from the posterior distribution $P(D|\bm{\theta},H)P(\bm{\theta}|H)$. The MCMC algorithm takes a prior $P(\bm{\theta}|H)$ based on some physics principle, and combines it with the data where the likelihood $P(D|\bm{\theta},H)$ is maximized to produce the posterior distribution. Each sample in the MCMC simulation represents a point in the multi-dimensional parameter space. The calculation of the marginalized probability distribution of a parameter say $\theta_{1}$

[TABLE]

is straightforward: just look at the distribution of the coordinate $\theta_{1}$ in the sample which trivially corresponds to summing over all the other coordinates Brewer (2018).

Model comparison in Bayesian analysis requires calculating the relevant evidences. In $\chi^{2}$ minimization, one might get a better fit by including extra parameters. However, in Bayesian analysis there is a cost associated with including extra parameter as it introduces another prior distribution. The probability ratio of model A ($M_{A}$) and model B ($M_{B}$), for a data set $D$, is given by the posterior odd ratio:

[TABLE]

The first factor on the right is just the ratio of evidences for $M_{A}$ and $M_{B}$. The second factor is called the prior odd. In our analysis, when we compare models we assume no prior bias towards or against any model, and accordingly set the prior odd to 1.

The evidence calculation requires numerical integration in the multi-dimensional parameter space. If we label the parameters of $M_{A}$ by the vector $\bm{\alpha}$, and the parameters of $M_{B}$ by the vector $\bm{\beta}$, then the evidences are:

[TABLE]

The difficult task of calculating the evidences in the multi-dimensional parameter space is made easier in the Nested Sampling method introduced by Skilling Skilling (2006). The basic idea is to draw a set of points randomly from the prior distribution $P(\bm{\theta},H)$, and arrange them from best to worst in terms of the likelihood function. Then one discards the worst point, and replaces it with another point from the prior distribution with the condition that it has a larger likelihood value. This guarantees that progressively higher likelihood regions of the parameter space is probed. In the original work by Skilling, at every iteration the worst point was replaced by Monte Carlo update using the Metropolis-Hastings algorithm. There has been a lot of algorithmic development since then on how to update the list of points, see review article Brewer (2018). We tried several of these methods. For our expressions, these methods disagree by about $\sim 2$ in the calculation of $\ln[P(M_{A}|D,H)/P(M_{B}|D,H)]$, and we assume this level of systematic errors in model comparisons. We decided to use the MultiNest algorithm Feroz et al. (2009) implemented in Python Barbary . We cross checked the numbers with our own Fortran code using the original Nested Sampling algorithm Skilling (2006).

IV $\bm{S}$-factor estimation for 3He$\bm{(\alpha,\gamma)^{7}}$Be and 3H$\bm{(\alpha,\gamma)^{7}}$Li

The cross section (and the related $S$-factor) expressions for ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ are given by Eq. (1). It depends on the amplitudes $\mathcal{A}(p)$ and $Y(p)$. The latter is expressed in Eq. (19) without any unknown parameters. The former depends on the initial state strong interaction $s$-wave phase shift $\delta_{0}$ that can be written model-independently in terms of scattering parameters: scattering length $a_{0}$, effective range $r_{0}$, shape parameter $s_{0}$, etc., in Eq (A). We keep only three $s$-wave scattering parameters $a_{0}$, $r_{0}$ and $s_{0}$ in our analysis. $\mathcal{A}(p)$ also depends on 2 two-body current couplings $L^{(\zeta)}_{\mathrm{E1}}$ for $\zeta={}^{2}P_{3/2},\ \zeta={}^{2}P_{1/2}$. The cross section is multiplied by the overall factor $([h^{(\zeta)}]^{2}\mathcal{Z}^{(\zeta)})$, Eq. (20), that depends on the binding momenta $\gamma$ and the effective ranges $\rho_{1}$ in the $p$-wave channels ${}^{2}P_{3/2}$ and ${}^{2}P_{1/2}$. This amounts to 7 parameters/couplings for the capture cross section as we take the binding momenta as given. We use the notation $\pm$ to indicate the ${}^{2}P_{3/2}$ and ${}^{2}P_{1/2}$ channels, respectively.

A $\chi^{2}$ fit to capture data Higa et al. (2018) for ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ does not constraint the $p$-wave effective ranges $\rho_{1}^{(\pm)}$ accurately. The wave function normalization constants $[h^{(\pm)}]^{2}\mathcal{Z}^{(\pm)}$ have poles at approximately $\rho_{1}^{(+)}=-47.4$ MeV and $\rho_{1}^{(-)}=-32.4$ MeV, respectively. A small change in the effective range values changes the normalization of the cross section significantly. Moreover, one finds that for the parameters from the $\chi^{2}$ fit, there is a large cancellation among the terms in $\mathcal{A}(p)$ such that the two-body contribution becomes important. However, the contribution from the two-body currents can be suppressed by choosing the effective range contribution near the pole in $([h^{(\zeta)}]^{2}\mathcal{Z}^{(\zeta)})$. Therefore, we find it important to include information from $p$-wave elastic phase shifts $\delta_{1}^{(\pm)}$ to constraint $\rho_{1}^{(\pm)}$ more accurately. This entails including the shape parameters $\sigma_{1}^{(\pm)}$ in the fit. $\sigma_{1}^{(\pm)}$ contributes to phase shifts $\delta_{1}^{(\pm)}$ but not to the capture cross section. The total parameters to be fit now increases to 9: $a_{0}$, $r_{0}$, $s_{0}$, $\rho_{1}^{(\pm)}$, $\sigma_{1}^{(\pm)}$, $L_{\mathrm{E1}}^{(\pm)}$.

We start with the new analysis of ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ Higa et al. (2018). The capture and elastic scattering data is more accurately known for this reaction than for ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$.

IV.1 Capture reaction 3He$\bm{(\alpha,\gamma)^{7}}$Be

We use charge $Z_{\psi}=2$, mass $M_{\psi}=2809.41$ MeV for the incoming 3He, and charge $Z_{\phi}=2$, mass $M_{\phi}=3728.4$ MeV for the incoming $\alpha$ in Eq. (1). The ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ capture data are from ERNA Di Leva et al. (2009), LUNA Constantini et al. (2008), Notre Dame Kontos et al. (2013), Seattle Brown et al. (2007), and Weizmann Nara Singh et al. (2004). It includes prompt photon, activation and recoil (ERNA) measurements. For the prompt data, the branching ratio $R_{0}$ of capture to the excited state to the ground state is available.

In the Bayesian analysis we explore fits to capture data up to c.m. energy (momentum) $E\lesssim 1$ MeV ($p\lesssim 60$ MeV) that we call region I. Fits to capture data over $E\lesssim 2$ MeV ($p\lesssim 80$ MeV) we call region II. The phase shift data for $\delta_{0}$, $\delta_{1}^{(\pm)}$ is from an old source Spieger and Tombrello (1967) that was analyzed by Boykin et al. in Ref. Boykin et al. (1972). The phase shift data starts from around $E\sim 2$ MeV ($p\sim 80$ MeV), which is at the higher end of the range of applicability of the EFT. We fit with phase shift data, and also without (indicated by an asterisk $\ast$ as appropriate). In all, we have four combinations of data sets for the fits. Region I includes 42 $S$-factor data and 20 $R_{0}$ data points whereas region II includes 70 $S$-factor data and 32 $R_{0}$ data points. The phase shifts include 10 data points in each of the three channels.

In the EFT, the momentum $p\sim 70$ MeV and final state binding momenta $\gamma_{0}\sim\gamma_{1}\sim 60\mbox{--}70$ MeV constitutes the small scale $Q$. The pion mass, excited states of 3He, $\alpha$, etc., constitutes the cutoff scale $\Lambda\gtrsim 150\mbox{--}200$ MeV. In the $s$-wave capture amplitude $\mathcal{A}(p)$, Eq. (12), the first term $X(p)$ that only contains Coulomb interaction is $O(1)$ at small momentum. The rest of the contributions, at arbitrarily small momentum $p$, scales as $a_{0}(B+\mu J_{0})/\mu^{2}$ and $a_{0}k_{0}L_{\mathrm{E1}}^{(\zeta)}$. The linear combination $(B+\mu J_{0})/\mu^{2}$ scales as $Q^{3}/\Lambda^{2}$, and the photon energy $k_{0}=(p^{2}+\gamma^{2})/(2\mu)\sim Q^{3}/\Lambda^{2}$ Higa et al. (2018, 2008). Thus at low momentum the contributions from initial state strong interaction $a_{0}(B+\mu J_{0})/\mu^{2}\sim a_{0}Q^{3}/\Lambda^{2}$. The two-body current contribution $a_{0}k_{0}L_{\mathrm{E1}}^{(\zeta)}\sim a_{0}Q^{3}/\Lambda^{2}$ for a natural sized dimensionless coupling $L_{\mathrm{E1}}^{(\zeta)}\sim 1$. The size of the $s$-wave scattering length then determines the relative contributions in the amplitude at low momentum.

$\chi^{2}$ minimization that included two-body current contribution gave a large $a_{0}\sim 20\mbox{--}30\ \mathrm{fm}\sim\Lambda^{2}/Q^{3}$, and consequently the initial state interactions $a_{0}(B+\mu J_{0})/\mu^{2}\sim 1$ and two-body currents $a_{0}k_{0}L_{\mathrm{E1}}^{(\zeta)}\sim 1$ are LO contributions. $\chi^{2}$ minimization without two-body currents give a smaller $a_{0}\sim 5-10\ \mathrm{fm}\sim\Lambda/Q^{2}$ Higa et al. (2018). This will make initial state interaction $a_{0}(B+\mu J_{0})/\mu^{2}\sim Q/\Lambda$ a NLO contribution. In EFT, there is always going to be two-body currents unless some symmetry prevents them. Then the only consideration is at what order in perturbation they contribute. Again assuming natural sized couplings $L_{\mathrm{E1}}^{(\zeta)}\sim 1$, the two-body contribution $a_{0}k_{0}L_{\mathrm{E1}}^{(\zeta)}\sim Q/\Lambda$ is also a NLO effect for $a_{0}\sim\Lambda/Q^{2}$.

We consider two power countings. The first one assumes a large $a_{0}\sim\Lambda^{2}/Q^{3}$. This we call “Model A” though EFT is a model-independent formulation. In this power counting Higa et al. (2018), all the other parameters and couplings are natural sized determined by their naive dimensions $r_{0}\sim 1/\Lambda$, $s_{0}\sim 1/\Lambda^{3}$, $\rho_{1}^{(\pm)}\sim\Lambda$, $\sigma_{1}^{(\pm)}\sim 1/\Lambda$, $L_{\mathrm{E1}}^{(\pm)}\sim 1$. The capture cross section depends on $a_{0}$, $r_{0}$, $\rho_{1}^{(\pm)}$, $L_{\mathrm{E1}}^{(\pm)}$ at LO. The NLO contribution gets an additional contribution from $s_{0}$. The $s$-wave phase shift $\delta_{0}$ depends on $a_{0}$ and $r_{0}$ at LO, and at NLO brings in $s_{0}$. A fine-tuning in the linear combination $r_{0}p^{2}/2-2k_{C}H(\eta_{p})\sim Q^{3}/\Lambda^{2}$ promotes the effective range contribution to LO Higa et al. (2018). The $p$-wave phase shifts $\delta_{1}^{(\pm)}$ depends on $\rho_{1}^{(\pm)}$ at LO and the NLO contribution comes from $\sigma_{1}^{(\pm)}$. The $d$ wave contribution $Y(p)$ is included at NLO due to a large near cancellation at LO Higa et al. (2018).

In the second power-counting that we call “Model B”, one assumes a smaller $a_{0}\sim\Lambda/Q^{2}$. The rest of the parameters have the same scaling as Model A. However, the perturbative expansion is now different. The LO capture cross section has no initial state strong interaction. It only depends on $\rho_{1}^{(\pm)}$. $a_{0}$ and $L_{\mathrm{E1}}^{(\pm)}$ contributes at NLO, and $r_{0}$ at next-to-next-to-leading order (NNLO). In phase shift $\delta_{0}$, $a_{0}$ contributes at LO, $r_{0}$ at NLO and $s_{0}$ at NNLO. The $p$-wave phase shift contributions remain unchanged from Model A above. The $d$ wave contribution $Y(p)$ is now included at NNLO.

Model A and Model B are compared by calculating the posterior odd ratio $P(M_{A}|D,H)/P(M_{B}|D,H)$ from Eq. (8). The $\chi^{2}$ for the fit without two-body currents was larger than the one with two-body currents when phase shift data, especially $\delta_{0}$, was used. We explore the possibility that the uncertainty in the phase shifts are actually larger than estimated. We consider $\sigma^{2}\rightarrow K^{2}+\sigma^{2}$ to describe an unaccounted noise. We draw $K$ from the uniform distribution $U(0^{\circ},10^{\circ})$. We also consider the possibility $\sigma^{2}\rightarrow K^{2}\sigma^{2}$ with $K$ drawn from $U(1,10)$. Both of these give similar fits so we present the analysis for $K^{2}+\sigma^{2}$ only. Note that the errors $\sigma$ in the phase shift data are estimated to be around $\sim 1^{\circ}\mbox{--}3^{\circ}$ Boykin et al. (1972). Alternatively, one could also explore the possibility of an uncertainty in the overall normalization of the phase shift measurements. As the $s$ and $p$ wave phase shifts are from the same measurement and analysis, we did not consider the possibility of introducing separate measurement errors of the different elastic scattering channels. We use the following uniform prior distributions for the parameters and couplings in the EFT expressions:

[TABLE]

Fits without phase shifts do not depend on $\sigma_{1}^{(\pm)}$ and $K$. Model B without phase shift data also does not depend on $s_{0}$. The range for each of the uniform distributions above was guided by EFT power counting estimates. The ranges were wide enough that about 95% of the posterior distributions of the parameters and couplings are not pressed against the boundaries. For certain parameters such as the $p$ wave effective ranges $\rho_{1}^{(\pm)}$ physical constraint that the ANCs be positive determines the upper bounds. All these assumptions are considered part of the background information in the proposition $H$ in the probability distributions.

In Figs. 1 and 2, we present fits in region I with and without phase shifts. In Fig. 1, the upper panels (a) and (b) are fits with phase shift data, and the lower panels (c) and (d) are without. The panels on the left in Fig. 2 are fits to phase shift data. The ones on the right are phase shift predictions from capture data. We include the $\chi^{2}$ fit of Model A to data in region I (with phase shift) for comparison Higa et al. (2018). The parameter estimates are in Table 5. The Bayesian curves are drawn from the posterior distributions of the parameters, not directly from the median values from the table. The posterior odd for the fit with phase shift is $\ln[P(M_{A}|D,H)/P(M_{B}|D,H)]\approx 2.6\pm 0.8$ slightly favoring Model A, though the odds are similar if one accounts for systematic error. For the fit without phase shift data $\ln[P(M_{A}^{\ast}|D,H)/P(M_{B}^{\ast}|D,H)]\approx 4.9\pm 0.6$ is largely in favor of model A. We remind the reader that the asterisk $\ast$ indicates fits without phase shift data. From Fig. 1 panel (a), we see Model A (with phase shift) reproduces capture data in region II better than the other fits though it was fitted only to region I data.

Figs. 3 and 4 represent a similar analysis but fitted to capture data in region II. The parameter estimates are in Table 5. The posterior odds are $\ln[P(M_{A}|D,H)/P(M_{B}|D,H)]\approx 31.6\pm 0.9$ and $\ln[P(M_{A}^{\ast}|D,H)/P(M_{B}^{\ast}|D,H)]\approx 0.2\pm 0.8$. These fits indicate that Model B is not able to reproduce both phase shift and capture data simultaneously over a large region. The large $K$ value in Table  5 for Model B II corroborates this. Without constraint from phase shift, model A∗ II and model B∗ II are equally favored.

IV.2 Capture reaction 3H$\bm{(\alpha,\gamma)^{7}}$Li

In the ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ calculation we use $Z_{\psi}=1$ and mass $M_{\psi}=2809.43$ MeV for the incoming 3H. The charge and mass for the $\alpha$ is the same as before. The data for the capture cross section are from Caltech Brune et al. (1994) and Tomsk Bystritsky et al. (2017). As in the case of 7Be, if we associate the binding momenta for the ground and 1st excited of 7Li $\gamma_{0}\sim\gamma_{1}\sim 80\mbox{--}90$ MeV with the scale $Q$, then for this system we expect the convergence to be slower compared to the previous 3He-$\alpha$ system which had smaller $\gamma_{0}$, $\gamma_{1}$.

The role of the initial state interaction and the two-body current contributions are predicated on the size of the $s$-wave scattering length $a_{0}$ as before. However, in this system a larger $a_{0}\sim\Lambda^{2}/Q^{3}$ corresponds to only $\sim 10\mbox{--}20\ \mathrm{fm}$. A smaller $a_{0}\sim\Lambda/Q^{2}\sim 1\mbox{--}5\ \mathrm{fm}$. The larger $a_{0}$ would make the initial state strong interaction and two-body current contributions LO, and the smaller $a_{0}$ would make these contributions NLO. The power-counting arguments for ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ are similar to the ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ capture. The difference is that now $Q$ is numerically larger.

The phase shift information of the initial 3H-$\alpha$ state is less well known than the 3He-$\alpha$ state. We digitized the $s$-wave phase shift $\delta_{0}$ for momentum range $p\sim 70\mbox{--}90$ MeV from Fig. 15 in Ref. Spieger and Tombrello (1967). The $p$ waves appear very noisy and we decided not to use these in the fits. The authors estimated an error of about $5^{\circ}$ in the phase shift analysis for 3He-$\alpha$ but the errors in the 3H-$\alpha$ phase shifts are not listed. We find a mention of 3% systematic error in 3H-$\alpha$ that we use. It is likely the error in this system is similar to 3He-$\alpha$. We add an additional noise as before and use $K^{2}+\sigma^{2}$ in the fits with $K\sim U(0^{\circ},10^{\circ})$. We use the following priors:

[TABLE]

The fits do not depend on the $p$-wave shape parameters $\sigma_{1}^{(\pm)}$ since we do not use the $p$-wave phase shifts. The upper limits on the effective ranges $\rho_{1}^{(\pm)}$ are constrained by requiring positive ANCs. We use 18 $S$-factor $S_{34}$ data, 17 branching ratio $R_{0}$ data and $7$ phase shift $\delta_{0}$ data points for the fits.

Figs. 5 and  6 shows the Bayesian fits. The power counting for Model A and Model B are the same as in ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$. We include a $\chi^{2}$ fit with phase shift using Model A for comparison. The upper panels (a) and (b) of Fig. 5 are for fits with $s$-wave phase shift, and the lower panels (c) and (d) are without. Fig. 6, left panel are fits to phase shift whereas the right panel are the predictions from fits to capture data. The posterior odds overwhelmingly favor Model A with $\ln[P(M_{A}|D,H)/P(M_{B}|D,H)]\approx 19.7\pm 0.7$, $\ln[P(M_{A}^{\ast}|D,H)/P(M_{B}^{\ast}|D,H)]\approx 18.4\pm 0.6$. The parameter estimations are in Table 6. The fits to phase shift suggests errors of about $5^{\circ}$ in $\delta_{0}$.

V Discussion and Summary

We performed several analysis of the EFT cross sections for ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$  and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ to draw Bayesian inference. Two EFT power counting proposals were compared. Typically, LO cross sections in EFT do not depend on two-body currents. This allows one to perform capture calculations once a few LO couplings are determined from elastic scattering data if they are known. For the two capture reactions considered here, elastic scattering data is not well measured. Available data on elastic scattering in 3He-$\alpha$ Boykin et al. (1972); Spieger and Tombrello (1967) and 3H-$\alpha$ Spieger and Tombrello (1967) systems suggest the $s$-wave scattering lengths $a_{0}$s are fine tuned to be large. This makes initial state interactions in the capture reactions dominant, and as a consequence the two-body currents could contribute at LO Higa et al. (2018). We called this EFT power counting Model A. Alternatively, one could have a power counting where two-body currents appear at higher order in the perturbation. This would require a smaller $s$-wave scattering length. We call this power counting Model B. If we ignore the poorly known elastic scattering data, then this second power counting can be used to describe the capture data. The relatively smaller initial state interaction contributions in Model B can be compensated by a larger wave function normalization constant which requires only small variation in the respective $p$-wave effective ranges as described earlier. We discuss the results of the analysis for ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$  and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ separately below.

V.1 3He$\bm{(\alpha,\gamma)^{7}}$Be

We present the 8 different fits used to draw Bayesian inference for this reaction. First we start with the fits in the smaller capture energy region I ($E\lesssim 1000$ keV). The posterior odds favored Model A both with and without phase shift data in the fits. However, if we also look at the overall trend then we see the data for $S$-factor $S_{34}$ rises upward from around $E\sim 1500$ keV. The Model A fit with phase shift, we call Model A I, best describes the capture data over the energy range $E\lesssim 2000$ keV. We note that the Model B fits in this region are not self-consistent in that they suggest an $a_{0}$ value larger than the power counting, table 5.

In the fits to capture energy region II ($E\lesssim 2000$ keV), Model A with phase shift, we call Model A II, is favored over Model B by the posterior odd. For the fits without phase shifts, both Model A∗ II and Model B∗ II are equally favored by the posterior odd ratio. The asterisk $\ast$ indicates fits without phase shift data.

Tables 1 and  2 has the $S$-factors $S_{34}$ and branching ratios $R_{0}$ for the 4 fits described above. We also include the $\chi^{2}$ fit of Model A for comparison Higa et al. (2018). We include the derivative $S^{\prime}_{34}$ as well. All the numbers were evaluated at $E_{\star}=60\times 10^{-3}$ keV. We include the estimated EFT errors. The NLO Model A results have a 10% error, and the NNLO Model B results have a 3% error. The different EFT error estimates has to do with the distinction between “accuracy and precision”. The error estimates from higher order corrections represent precision, and different power countings have different accuracy and precision.

Fig. 7 shows the posterior distributions for the 4 $S$-factors and branching ratios from Tables 1 and  2. The symmetric distributions can be described with a Gaussian form shown by the various smooth curves unlike the skewed distributions as expected. The spread in some of the quantities is related to the uncertainty in the parameter estimates, especially the $p$-wave effective ranges $\rho_{1}^{(\pm)}$ for Model A fits, see Table 5. Large magnitude $|\rho_{1}^{(\pm)}|$ makes the wave function normalization constant smaller which can be compensated by a larger two-body current $L_{1}^{(\pm)}$ as the parameter estimates indicate. The planned TRIUMF 3He-$\alpha$ elastic scattering experiments and phase shift analysis at low energies $E\gtrsim 500$ keV would be able to shed some light on this Davids (2016), and help establish the appropriate EFT power counting.

The EFT $S$-factor at threshold can be compared to other recent calculations such as – $0.593$ keV b from FMD Neff (2011), $0.59$ keV b from NCSM Dohet-Eraly et al. (2016); and evaluations such as – $[0.580\pm 0.043(\mathrm{stat.})\pm 0.054(\mathrm{sys.})]$ keV b from Cyburt-Davids Cyburt and Davids (2008), $(0.57\pm 0.04)$ keV b from ERNA Di Leva et al. (2009), $(0.567\pm 0.018\pm 0.004)$ keV b from LUNA Constantini et al. (2008), $(0.554\pm 0.020)$ keV b from Notre Dame Kontos et al. (2013), $(0.595\pm 0.018)$ keV b from Seattle Brown et al. (2007), and $(0.53\pm 0.02\pm 0.01)$ keV b from Weizmann Nara Singh et al. (2004). We also compare to a recent EFT work using Bayesian inference Zhang et al. (2018) that finds at threshold $S_{34}(0)=0.578^{+0.015}_{-0.016}$ keV b and $R_{0}(0)=0.406^{+0.013}_{-0.011}$. Ref. Zhang et al. (2018) only used capture data to draw their inferences. Looking at tables  1 and  2, it would seem the results of Ref. Zhang et al. (2018) are more aligned with Model B∗ II, though the exact power counting used the article is not clear. The “best” recommended value from the review in Ref. Adelberger et al. (2011) is: $S_{34}(0)=[0.56\pm 0.02(\mathrm{expt.})\pm 0.02(\mathrm{theory})]$ keV b.

From the various fits, we recommend the following: Model A II if we want to include phase shift information, and either Model A∗ II or Model B∗ II if no phase shift information is used.

V.2 3H$\bm{(\alpha,\gamma)^{7}}$Li

To draw our Bayesian inferences, we performed 4 different fits in this system. Here the EFT power counting of Model A (with phase shift data) and Model A∗ (without phase shift data) are favored. Moreover, the Model B and Model B∗ fits are not consistent with the power counting estimate for the $a_{0}$ values, table 6. The $S$-factors and branching ratios are in Tables 3 and 4. The corresponding posterior distributions, and the associated Gaussian forms are shown in Fig. 8. The fit without phase shift data gives a larger $S_{34}$ at threshold. We included a 16% EFT error estimate from NNLO corrections. In this channel the binding energies are larger, resulting in a larger perturbation error compared to ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$.

The fits, especially without the phase shift data, give a central $s$-wave shape parameter value $s_{0}$ that is larger than ideally expected from the power counting. Power counting consistency would bias one towards Model A (with phase shift) over Model A∗ (without phase shift).

The EFT $S$-factor calculation can be compared to recent theoretical results – $0.12$ keV b from FMD Neff (2011), $0.13$ keV b from NCSM Dohet-Eraly et al. (2016); and experimental evaluation – $[0.1067\pm 0.0004(\mathrm{stat.})\pm 0.0060(\mathrm{sys.})]$ keV b from Caltech Brune et al. (1994).

We leave for future a more careful study of the systematic errors associated with the estimation of the parameters and evidences using the various Nested Sampling methods. In our evaluations, for example, the MultiNest Feroz et al. (2009) and Diffusive Nested Sampling Brewer et al. (2011) algorithms gave similar results. The EFT parameters (tables  5 and  6) agreed within the error bars. Similarly, the capture cross sections and branching ratios (tables  1,  2,  3,  4) also agreed within the error bars from the fits.

Acknowledgements.

The authors thank Barry Davids, Renato Higa, Daniel R. Phillips, Prakash Patil and Xilin Zhang for many valuable discussions. This work was supported in part by U.S. NSF grants PHY-1615092 and PHY-1913620. The figures for this article have been created using SciDraw Caprio (2005).

Appendix A $\bm{s}$ and $\bm{d}$ wave capture

The capture from initial $s$-wave state is given by:

[TABLE]

where $k_{0}$ is the photon energy and $L^{(\zeta)}_{E1}$ for $\zeta={}^{2}P_{3/2}$ and $\zeta={}^{2}P_{1/2}$ are the 2 two-body currents. We modified some definitions compared to Ref. Higa et al. (2018) but they are equivalent expressions.

The contribution without initial state strong interaction is given by

[TABLE]

where $\gamma$ is the binding momentum of the ground or first excited state as appropriate. The momentum $k_{C}=\alpha_{e}Z_{\psi}Z_{\phi}\mu$ is the inverse of the Bohr radius. The Coulomb expressions are:

[TABLE]

with conventionally defined Whittaker functions $M_{k,\mu}(z)$ and $W_{k,\mu}(z)$. $F_{l}(\eta_{p},\rho)$ is the regular Coulomb wave function.

The initial $s$-wave scattering in Eq. (12) is contained in the Coulomb subtracted phase shift $\delta_{0}$ parameterized by the modified effective range expansion

[TABLE]

$\psi(x)$ is the digamma function, and the $\cdots$ represents terms with higher powers in $p^{2}$. The combination $B(p)+\mu J_{0}(p)$ can be evaluated as

[TABLE]

The function $C(p)$ is given by the double integral

[TABLE]

which is reduced to a single integral before evaluating numerically. The function $\Delta B(p)$ is obtained from

[TABLE]

The integral $B_{ab}(p)$ is divergent at $r=0$. However, when combined with the contribution from $J(p)\delta_{ab}$ it is finite. Thus we make the substitution $(r_{b}/r)[\partial/\partial r_{a}]=(r_{a}r_{b}/r^{2})[\partial/\partial r]\rightarrow(\delta_{ab}/3)[\partial/\partial r]$ in the integral and accordingly $B_{ab}(p)\equiv B(p)\delta_{ab}$. The finite piece $\Delta B(p)$ is obtained numerically from $B(p)$ after subtracting the zero and single photon contributions i.e. removing terms up to order $\alpha_{e}^{2}$ Higa et al. (2018).

The capture from initial $d$-wave states is given by the amplitude Higa et al. (2018):

[TABLE]

Appendix B Wave function renormalization

The wave function renormalization constant is calculated as Higa et al. (2018):

[TABLE]

where we used the relation

[TABLE]

Appendix C Parameter Estimates

We list the parameter estimates for both ${}^{3}\mathrm{He}(\alpha,\gamma)^{7}\mathrm{Be}$ and ${}^{3}\mathrm{H}(\alpha,\gamma)^{7}\mathrm{Li}$ below. For the Bayesian fits we show the median and the bounds containing 68% of the posterior distributions. Typically, a very asymmetric bound indicates either a skewed or sometimes a bi-modal distribution.