Lattice QCD Calculation of Electroweak Box Contributions to Superallowed Nuclear and Neutron Beta Decays

We present the first lattice QCD calculation of the universal axial $\gamma W$-box contribution $\square_{\gamma W}^{VA}$ to both superallowed nuclear and neutron beta decays. This contribution emerges as a significant component within the theoretical uncertainties surrounding the extraction of $|V_{ud}|$ from superallowed decays. Our calculation is conducted using two domain wall fermion ensembles at the physical pion mass. To construct the nucleon 4-point correlation functions, we employ the random sparsening field technique. Furthermore, we incorporate long-distance contributions to the hadronic function using the infinite-volume reconstruction method. Upon performing the continuum extrapolation, we arrive at $\square_{\gamma W}^{VA}=3.65(8)_{\mathrm{lat}}(1)_{\mathrm{PT}}\times10^{-3}$. Consequently, this yields a slightly higher value of $|V_{ud}|=0.97386(11)_{\mathrm{exp.}}(9)_{\mathrm{RC}}(27)_{\mathrm{NS}}$, reducing the previous $2.1\sigma$ tension with the CKM unitarity to $1.8\sigma$. Additionally, we calculate the vector $\gamma W$-box contribution to the axial charge $g_A$, denoted as $\square_{\gamma W}^{VV}$, and explore its potential implications.


Introduction:
The Cabibbo-Kobayashi-Maskawa (CKM) matrix plays a crucial role as a fundamental ingredient of the standard model (SM) where its unitarity is expected.Recent findings indicate a 2.1σ tension in the examination of the first row's unitarity [1] |V ud | 2 + |V us | 2 + |V ub | 2 = 0.9985 (6) Depending on the choice of inputs, tensions can be 3σ or larger [2].Further efforts to accurately determine the CKM matrix elements could unveil new signals and deepen our understanding of the underlying physics.
The top-left corner element |V ud | is most precisely extracted from neutron and superallowed nuclear beta decays.The latter currently claims the highest experimental accuracy, but the former is catching up thanks to improvements in experimental precision [3][4][5][6].On the theory side, superallowed nuclear decays suffer from extra nuclear-structure-dependent uncertainties that are recently under careful scrutiny [7][8][9], while neutron decay is theoretically cleaner.Nevertheless, a major source of uncertainty common to both cases is the single-nucleon γW -box diagrams (see Fig. 1) that renormalize the neutron vector and axial charges [10].They take the follow-ing forms [11] [12]: ϵ µναλ q α p λ (q 2 ) 2 T γW µν , ϵ µναλ q α S λ (q 2 ) 2 T γW µν , T γW µν = ∫ d 4 xe iq⋅x ⟨p(p, S)| T {J em µ (x)J W ν (0)} |n(p, S)⟩ , (2) where p and S are the momentum and spin vectors of the nucleon states.gA is the nucleon axial charge, where the symbol ˚indicates that it is defined in the isospin limit.J em µ is the electromagnetic current and J W µ the weak current with J W,V µ and J W,A µ its vector and axialvector part, respectively.Specifically, in 0 + → 0 + superallowed beta decay, only ◻ V A γW comes into play, while in + free neutron beta decay, both ◻ V A γW and ◻ V V γW would contribute.Various approaches to computing ◻ V A γW range from earlier works by Marciano and Sirlin [13,14] to more recent dispersive analyses by Seng et al [7,15].The latter, in particular, improved the nonperturbative contribution for loop momentum square Q 2 ≤ 2 GeV 2 and unveiled the tension with the first-row CKM unitarity, which was also observed in several follow-up works [16][17][18][19].In the meantime, while the radiative correction to the axial charge g A does not directly affect the |V ud | extraction, it is necessary for comparing the experimentally measured g A with that computed using lattice QCD.The study of ◻ V V γW has so far included estimations inspired by the holographic QCD model [19] and dispersion relations [11].
Lattice QCD offers a direct nonperturbative approach to compute the box correction ◻ V A γW , especially for Q 2 ≤ 2 GeV 2 .First lattice calculations of ◻ V A γW were successfully conducted in the pion [20] and kaon channel [21,22], and have recently been confirmed by an independent lattice calculation [23].The data reported in [20] were also used for a joint lattice QCD -dispersion relation analysis [17].This letter extends this calculation to the neutron decay channel, which entails a direct computation of the nucleon four-point function at the physical pion mass.We also briefly discuss our numerical result of ◻ V V γW , and its implication to the radiative correction to axial charge.

Methodology:
The notations used in this work align with those used in [20].We define the hadronic function H V A µν within Euclidean space where H i/f represents the zero-momentum projected neutron/proton state, created by a smeared-source nucleon operator.The computation of box contribution ◻ V A γW involves a momentum integral M n (Q 2 ) is a weighted integral of the hadronic function with m W and m N the masses of the W -boson and the nucleon.The weighting function is where θ and j n (x) are the spherical Bessel functions.
To evaluate M n (Q 2 ) as prescribed in Eq. ( 5), it is necessary to extend the temporal integration range sufficiently to reduce truncation effects.However, as the time separation between the two currents increases, the lattice data tend to exhibit greater noise-to-signal ratio.
Here we employ the infinite volume reconstruction (IVR) method [24] to incorporate the long-distance (LD) contribution arising from the region where |t| > t s .Here, t s is the time slice at which the short-distance (SD) and LD contributions are separated.The IVR method, in addition to eliminating the power-law suppressed finite volume error, can also reduce the lattice statistical error in the long distance region.To elaborate, we divide the integral into SD contribution, weighted by ω(t, ⃗ x), and LD contribution, weighted by ω(t, ⃗ x) with and Here, H(t, ⃗ Q| 2 and t g is chosen to be large enough to ensure the groundintermediate-state dominance.Once t g is fixed, t s can be varied to further verify the ground-state dominance.In the final results, it is natural to choose t s = t g . Due to the factor 1/Q 2 in Eq. ( 4), we observe that ◻ V A γW encounters a notably increased noise originating from M n (Q 2 ) at small Q 2 region.To mitigate this noise, we can use the model-independent relation to substitute M LD n (Q 2 , t s , t g ) with Above, as far as ground-state dominance is satisfied, H(t g , ⃗ x) is independent of t g .μp,n denote the proton/neutron magnetic moments defined in the isospin limit.During the substitution process, we incorporate experimentally measured values for g A and µ p,n as depicted in Eq. (11).The difference is of a higher order and numerically negligible.Furthermore, ω0 is defined as Importantly, the convergence of the integral with ω−ω 0 at Q 2 → 0 is considerably faster than that with ω, enabling a more efficient control over statistical uncertainties.We refer to the calculation of M LD n using Eq. ( 8) and Eq.(11) as the "direct" and "substitution" methods, respectively.
We introduce a four-momentum squared scale Q 2 cut which separates the Q 2 integral into two regimes, = 3α e 2π ( ∫ we use lattice results as inputs.Con- , we utilize the perturbative QCD and employ the leading twist contribution from the operator product expansion [25][26][27].Further details can be found in Ref. [20].A common representative value for the scale of It is also feasible to vary this value to investigate potential systematic effects.

Numerical analysis:
We use two lattice QCD gauge ensembles at the physical pion mass, generated by RBC and UKQCD Collaborations using 2 + 1-flavor domain wall fermion [28].The ensemble parameters are outlined in Table I.Both ensemble utilize Iwasaki + DSDR action.For each configuration we produce 1024 pointsource and 1024 smeared-source propagators at random spatial-temporal locations and calculate the correlation function ⟨ψ p (t f )J em µ (x)J W,A ν (y)ψ † n (t i )⟩ with t f = max{t x , t y } + ∆t f and t i = min{t x , t y } − ∆t i using the random sparsening-field technique [29,30].Local vector and axial vector current operators are contracted with the renormalization factors quoted from Ref. [31].We calculate all the connected insertions, discarding disconnected insertions which vanish under the flavor SU(3) limit.
142.6(3) 24 64 1.023(2) 207 32D-fine 143.6(9) 32 64 1.378 (5) 69 Table I.Ensembles used in this work.For each ensemble we list the pion mass mπ, the spatial and temporal extents, L and T , the inverse of lattice spacing a −1 [32], the number of configurations used, N conf .The lattice spacing is determined using the mass of Ω baryon as input.14) as a function of ts.
Here, we employ 24D as an illustrative example.∆ti + ∆t f is fixed at 0.77 fm.
To incorporate the LD contribution, the appropriate values for t g and t s must be determined.We calculate the LD contribution to ◻ V A,≤2GeV 2 γW using M LD n (Q 2 , t s , t g ) as inputs, labelling the relevant part of the box contribution as ◻ V A,≤2GeV 2 γW (t s , t g ).For small t g values, a visible contamination from excited states is anticipated.To extend this analysis, we calculate ◻ V A,≤2GeV 2 γW (t s , t s ) with t g = t s for various t s values and construct a ratio as depicted in Fig. 2. Evidently, for t g below 0.39 fm, the ratio is not fully consistent with 1 for t s ≥ t g .However, at t g ≥ 0.58 fm, the LD contribution reconstructed using H(t s = t g , ⃗ x) and H(t s > t g , ⃗ x) agree well.This implies that at t g = 0.58 fm, the ground state begins to dominate the hadronic function, and statistical deviations between the reconstruction at t s > t g and t s = t g are negligible.
To further verify the ground-state dominance, we present in Fig. 3 the SD, LD and total contributions to ◻ V A,≤2GeV 2 γW for both the 24D and 32Dfine ensembles .The LD contribution is reconstructed utilizing the lattice data of H(t g , ⃗ x) with t g ≈ 0.6 fm.When the SD and LD contributions are combined, a discernible plateau emerges for t s ≳ 0.6 fm, indicating that the influence of excited-state contributions beyond this specific time slice is notably mild.
In Fig. 4 we present a plot of ◻ V A,≤2GeV  tion of ∆t i + ∆t f utilizing t s = t g ≈ 0.6 fm for both 24D and 32Dfine.The maximal source-sink separation extends to to ∼ 1.8 fm.In our analysis, we note that when √ Q 2 < 0.2 GeV, the substitution method's error is less than that of the direct method.Thus, we adopt the substitution method for calculating the momentum integral for ◻ V A,≤2GeV 2 γW when √ Q 2 < 0.2 GeV, and resort to the direct method for integrals exceeding this threshold.For ∆t i + ∆t f ≳ 0.7 fm lattice results manifest a discernible plateau for both 24D and 32Dfine ensembles.By fitting to a constant, we obtain ◻ V A,≤2GeV .The lattice result at the continuum limit is slightly lower than prediction from dispersive analysis.uum extrapolation, as carried out in Fig. 5, where we obtain ◻ V A,≤2GeV 2 γW = 1.490(51) × 10 −3 in the continuum limit.In the study of the pion decay [20], we observed that the continuum extrapolation derived from the two Iwasaki ensembles with finer lattice spacings, 48I/64I, deviates from the outcomes of the 24D/32Dfine ensembles by 2.4 × 10 −5 .Considering that lattice discretization effects predominantly arise from short-distance effects in the loop momentum integral, we anticipate consistency between the pion and the nucleon.In practice, we adopt a more conservative approach by assigning a twice-larger uncertainty, i.e. 4.8 × 10 −5 , for the nucleon than for the pion.A future calculation performed at finer lattice spacings would be beneficial for a more robust error estimate.The other systematic effects including finite-volume effects (FV), excited-state contaminations (ES) and contributions from disconnected digrams (disc) [33] are also estimated, with the details incorporated in the supplementary materials.Finally, we obtain This value is in agreement within 1σ with the value ◻ V A,≤2GeV 2 γW = 1.62(10) × 10 −3 obtained from dispersive analysis [7,15].
At the leading-twist level, the perturbative contribution to ◻ V A,>2GeV 2 γW aligns with that for the pion decay.Therefore, we utilize the 4-loop analysis detailed in Ref. [20] to deduce where the first error accounts for the higher-loop trunca-tion effects.This is determined by contrasting the results from 4-loop and 3-loop perturbative calculations.The second error pertains to higher-twist truncation effect and is gauged through diagrams involving two currents positioned on distinct quark lines, thus containing solely higher-twist contributions.Upon combining ◻ V A,>2GeV 2 γW and ◻ V A,≤2GeV 2 γW , the final result becomes To assess the Q 2 cut dependence, we also consider Q 2 cut = 1 GeV 2 and 3 GeV 2 , with remarkably consistent results, as is seen from Table II.Results and conclusion: According to the PDG [1], precise values of |V ud | 2 are obtained from superallowed nuclear and neutron beta decays [14,[34][35][36] as , superallowed, ( 18) For superallowed decays, the first uncertainty stems from analyzing the half-lives of 15 precisely measured decays [35].Additionally, an uncertainty associated with the nuclear structure, (54) NS , is quoted from Ref. [8].For neutron decays, the primary contributors to the uncertainty of |V ud | are of the experimental origin, namely the neutron lifetime τ n = 878.4(5)s and the ratio of axialvector to vector couplings, g A = 1.2754( 13) [1].It is worth noting that due to the SM radiative corrections g A deviates from the axial charge gA calculated in lattice QCD.∆ V R denotes a universal, nuclear-structureindependent electroweak radiative correction.In the framework proposed by Sirlin [10], ∆ V R is given by with ãg = −0.083representing the O(α s ) QCD correction to all one-loop diagrams except the axial γW box.Additionally, δ QED HO = 0.00109(10) [36] summarizes the leading-log higher-order QED effects, which can be accounted for via the running α e (see an updated discussion in Ref. [37]).Using the lattice results for ◻ V A γW , we obtain |V ud | = 0.97386 (11) exp.( 9) RC ( 27) NS (22) for superallowed beta decays.Upon entering |V ud |, the uncertainties from both lattice results and δ QED HO contribute collectively as (9) RC , with the lattice uncertainty accounting for 83%.Likewise, for free neutron decay, We note that the primary source of uncertainty is attributed to the experimental measurements of g A .In this study, we also calculate the γW box correction to g A using the formula provided in [19] g Interestingly, the vector γW -box term largely cancels the axial one, resulting in a value of ◻ V V γW − ◻ V A γW = 0.07(11) × 10 −3 that aligns with zero.This outcome also agree with the value of 0.13(13) × 10 −3 derived from dispersive analysis [11].It is essential to emphasize that the γW box correction to g A doesn't constitute the dominant contribution.Ref. [38] suggests that the most substantial radiative correction originates from vertex corrections, and may approach a percent level.To address this possibility on the lattice, one would need to include, e.g., the radiative corrections to the axial form factors ("f.f." in the equation above).This, in turn, requires computing 5-point correlation functions on the lattice.
To conclude, in this work we perform the first lattice QCD calculation of the γW -box contribution to the superallowed nuclear and free neutron beta decays.As depicted in Fig. 5, lattice results for both ensembles are consistent with the outcome obtained from dispersive analysis.The only distinction arises in the slightly smaller value exhibited by the continuum-extrapolated result in comparison to the dispersive analysis.It is always advisable to have future independent check using different discretizations of QCD.Using the updated |V ud | as provided in Eq. ( 22) and combining it with |V us | = 0.2243 (8) from PDG [1], we obtain This result exhibits 1.8σ tension with CKM unitarity.

Quark contractions for nucleon 4-point correlation functions
The connected insertions encompass 10 distinct contraction types, as illustrated in Fig. S 1. Notably, types (b) and (d) do not contribute to the axial-vector and vector γW -box diagrams.This is because the quark current associated with the W boson, altering the isospin, cannot be inserted between isospin-0 diquark blocks.Demonstration of Eq. ( 8) We follow the derivation presented in Ref. [20] to express M n (Q 2 ) as Here, the hadronic tensor T V A µν (Q) is defined as The tensor T V A µν (Q) is split into long-distance and shortdistance parts ) Utilizing Eq. (S 4), T LD µν (t s , t g , Q) is expressed as Taking the real part, we obtain Inserting T LD µν (t s , t g , Q) into Eq.(S 1), we derive The third line of Eq. (S 7) is expressed in terms of H(t g , ⃗ x) through Inserting Eq. (S 8) into Eq.(S 7), we arrive at Eq. ( 8).
Demonstration of Eq. ( 10) Here we demonstrate that once the ground-state dominance is satisfied at |t| ≥ t g , the spatial summation of the hadronic function H(t, ⃗ x) can be written in terms of gA , μp and μn .We start with the expression where In the small | ⃗ Q| limit, the above expression can be simplified as lim where is the magnetic form factor.For t ≥ t g , we have G M (0) = μp .For t ≤ −t g , the expression is similar as Eq.(S 11), albeit with G M (0) = μn .We thus have lim On the other hand, we have Combining Eqs.(S 11) and (S 14), we obtain In our calculation, we find that the lattice results for the 24D, 32Dfine ensembles, and the continuum extrapolation are all consistent with the PDG value, as depicted in Table S  x H(tg, ⃗ x) for the 24D, 32Dfine ensembles, and the continuum extrapolation contrasted with the PDG value of −3g A (µp + µn).Here, tg takes a value of ∼ 0.6 fm.

Demonstration of the necessity of using the IVR method
To demonstrate the necessity of using the IVR method in our calculation, we use the ensemble 24D as an example and present in Fig. S juxtaposed with findings from dispersive analysis [7,15] and perturbation theory [25][26][27].Two curves from perturbation theory are presented.One is constructed using 4-flavor theory down to 1 GeV, while the other decouples the charm quark at 1.6 GeV and employs 3-flavor theory for the range 1 GeV ≤ √ Q 2 ≤ 1.6 GeV.The observed discrepancy between these two curves suggests the unreliability of perturbation theory at small √ Q 2 .Fortuitously, the contribution to ◻ V A γW from this specific momentum range remains relatively small.A substantially more significant contribution emerges from the √ Q 2 < 1 GeV region, where the lattice results generally exhibit reduced uncertainties in comparison to those derived from dispersive analysis.At √ Q 2 = 0.2 GeV, we compare the lattice results of M n (Q 2 )/ √ Q 2 obtained from the substitution and direct methods.For 24D and 32Dfine ensembles, the discrepancy is 0.006(49) GeV −1 and 0.035(53) GeV −1 , indicating that the two approaches yield consistent outcomes.It's worth noting that these results are highly concordant with the findings presented in Table S I.This convergence is unsurprising, considering that the disparity between the two methods stems from the contrast between the lattice calculation of ∫ d 3 ⃗ x H(t g , ⃗ x) and the PDG value of −3g A (µ p + µ n ).

Analysis of systematics
In our calculation, the primary source of systematic effects arises from lattice discretization effects.We employ the lattice results of pion γW box contributions from our previous calculation [20] as inputs to estimate the residual lattice discretization effects after continuum extrapolation for the nucleon.
To assess the finite-volume effect, we perform the spatial integral in M SD n (Q 2 , t s ) and M LD n (Q 2 , t s , t g ) within a range of |⃗ x| ≤ R. Taking the ensemble 24D, which has better statistics, as an example, we depict ◻ V A,≤2GeV 2 γW as a function of R in Fig. S 5. The final result is calculated at maximal R ≃ 4 fm, while the plateau is already evident at R ≃ 1.8 fm, as indicated by the dashed line in

Figure 1 .
Figure 1.The γW -box diagrams for the semileptonic decay process Hi → H f eνe.

Figure 2 .
Figure2.The ratio defined in Eq. (14) as a function of ts.Here, we employ 24D as an illustrative example.∆ti + ∆t f is fixed at 0.77 fm.

Figure 4 .
Figure 4. ◻ V A,≤2GeV 2 γW for 24D and 32Dfine as a function of ∆ti + ∆t f .Here ts and tg are maintained as ts = tg ≈ 0.6 fm.
2 1.32(7) ×10 −3 2.31(2) × 10 −3 3.63(7) × 10 −3 2 GeV 2 1.49(7) × 10 −3 2.16(1) × 10 −3 3.65(7) × 10 −3 3 GeV 2 1.59(7) × 10 −3 2.06(1) × 10 −3 3.65(7) × 10 −3 acknowledges support by DOE Office of Science Early Career Award No. DE-SC0021147 and DOE Award No. DE-SC0010339.The work of M.G. is supported in part by EU Horizon 2020 research and innovation programme, STRONG-2020 project under grant agreement No 824093, and by the Deutsche Forschungsgemeinschaft (DFG) under the grant agreement GO 2604/3-1.K.F.L. and B.G.W. are partially supported by the DOE Grant No. DE-SC0013065 and No. DE-AC05-06OR23177.The work of C.-Y.S. is supported in part by the U.S. Department of Energy (DOE), Office of Science, Office of Nuclear Physics, under the FRIB Theory Alliance award DE-SC0013617, by the DOE grant DE-FG02-97ER41014, and by the DOE Topical Collaboration "Nuclear Theory for New Physics", award No. DE-SC0023663.The work of B.G.W. is supported in part by the DOE Office of Science, Office of Nuclear Physics under the umbrella of the Quark-Gluon Tomography (QGT) Topical Collaboration with Award DE-SC0023646.The research reported in this work was carried out using the computing facilities at Chinese National Supercomputer Center in Tianjin.It also made use of computing and long-term storage facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy.

Figure S 1 . 8
Figure S 1. 8 out of the 10 quark contraction types are pertinent to the neutron γW -box diagrams.The symbols of ⊗ indicate the insertions of the vector or axial-vector current.

Fig. S 4
Fig. S 4 displays the lattice results forM n (Q 2 )/ √ Q 2 as a function of √ Q 2 ,juxtaposed with findings from dispersive analysis[7,15] and perturbation theory[25- 27].Two curves from perturbation theory are presented.One is constructed using 4-flavor theory down to 1 GeV, while the other decouples the charm quark at 1.6 GeV and employs 3-flavor theory for the range 1 GeV ≤ √ Q 2 ≤ 1.6 GeV.The observed discrepancy between these two curves suggests the unreliability of perturbation theory at small √ Q 2 .Fortuitously, the contribution to ◻ V A γW from this specific momentum range remains relatively small.A substantially more significant contribution emerges from the √ Q 2 < 1 GeV region,

Figure S 4 . 2 .
Figure S 4. Mn(Q 2 )/ √ Q 2 as a function of √ Q 2 .The parameters for ∆ti +∆t f , ts and tg are maintained at the same values as those employed in Fig. S 3. The lattice results for the 24D and 32Dfine ensembles are juxtaposed with results obtained from dispersive analysis and perturbation theory.The lattice data below and above √ Q 2 = 0.2 GeV are compiled using the substitution and direct methods, respectively.
The IVR is performed at tg = 0.58 fm for 24D and tg = 0.57 fm for 32Dfine.∆ti + ∆t f is set at ∼ 0.75 fm for both ensembles.

Table II .
Utilizing the scale Q 2 cut partition the integral range, we present the contributions of ◻ I.Table S I.The lattice results of ∫ d 3 ⃗ 2 the results of M SD n (Q 2 , t s ) as a function of Q 2 for different values of t s .Notably, even when increasing t s to 1.16 fm while maintaining ∆t i +∆t f 2with various choices of ts for 24D.∆ti+∆t f is set at 0.77 fm.The error band denoted as "SD+LD" is the lattice result which incorporates the reconstruction of the LD contributions using IVR and tg = ts = 0.58 fm.