First-principles calculation of electroweak box diagrams from lattice QCD

We present the first realistic lattice QCD calculation of the $\gamma W$-box diagrams relevant for beta decays. The nonperturbative low-momentum integral of the $\gamma W$ loop is calculated using a lattice QCD simulation, complemented by the perturbative QCD result at high momenta. Using the pion semileptonic decay as an example, we demonstrate the feasibility of the method. By using domain wall fermions at the physical pion mass with multiple lattice spacings and volumes, we obtain the axial $\gamma W$-box correction to the semileptonic pion decay, $\Box_{\gamma W}^{VA}\big|_{\pi}=2.830(11)_{\mathrm{stat}}(26)_{\mathrm{sys}}\times10^{-3}$, with the total uncertainty controlled at the level of $\sim1$\%. This study sheds light on the first-principles computation of the $\gamma W$-box correction to the neutron decay, which plays a decisive role in the determination of $|V_{ud}|$.

Introduction -The precise determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, which are fundamental parameters of the Standard Model, is one of the central themes in modern particle physics. In the CKM matrix, V ud is the most accurately-determined element from the study of superallowed 0 + → 0 + nuclear beta decays [1] which are pure vector transitions at tree level and are theoretically clean due to the protection of the conserved vector current. Going beyond tree level, the electroweak radiative corrections involving the axial-vector current become important and ultimately dominate the theoretical uncertainties.
Among various electroweak radiative corrections, the axial γW -boson box contribution ◻ V A γW contains a significant sensitivity to low-energy hadronic effects, and is a dominant source of the total theoretical uncertainty [2]. The recent dispersive analysis [3,4] reduced this uncertainty by a factor of 2 comparing to the previous study by Marciano and Sirlin [5], and the updated result of V ud raised a 4 standard-deviation tension with the first-row CKM unitarity (barring possibly underestimated nuclear effects: see Ref. [4,6]). The main difference between those works is the use of inclusive neutrino and antineutrino scattering data that Refs. [3,4] used to estimate the contribution of the intermediate momenta inside the γW loop integral, 0.1 GeV 2 ≲ Q 2 ≲ 1 GeV 2 , prone to nonperturbative hadronic effects. To further improve the determination of V ud , it requires either better-quality experimental input or the direct, precise lattice QCD calculations of the γW -box contribution.
Ref. [36] proposed to calculate the γW -box contribution using the Feynman-Hellman theorem. In this work, we opt for a more straightforward way to perform the lattice calculation. To demonstrate the feasibility of the method, we carry out the exploratory study for the case of the pion semileptonic decays. The calculation is performed at the physical pion mass with various lattice spacings and volumes, which allows us to control the systematic effects in the lattice results. Combining the results from lattice QCD together with the perturbative QCD, we obtain the axial γW -box correction to pion decay amplitude with a relative ∼ 1% uncertainty.
The γW -box contribution -In the theoretical analysis of the superallowed nuclear beta decay rates, the dominant uncertainty arises from the nucleusindependent electroweak radiative correction, ∆ V R , which is universal for both nuclear and free neutron beta decay [1]. Among various contributions to ∆ V R , Sirlin es-tablished [2] that only the axial γW -box contribution is sensitive to hadronic scales. The relevant hadronic tensor T V A µν is defined as (1) for a semileptonic decay process H i → H f eν e . Above, H i f are given by neutron and proton for the neutron beta decay, and by π − and π 0 for the pion semileptonic decay, respectively. Furthermore, J em µ = 2 3ū γ µ u − 1 3d γ µ d is the electromagnetic quark current, and J W,A ν =ūγ ν γ 5 d is the axial part of the weak charged current.
The spin-independent part of T V A µν has only one term, where T 3 is a scalar function. For the neutron beta decay, the spin-dependent contributions, denoted by the ellipses here, are absorbed into the definition of the nucleon axial charge g A , which can be measured directly from experiments. According to current algebra [2], it is this spin-independent term that gives rise to the hadron structure-dependent contribution and dominates the uncertainty in the theoretical prediction. Using T 3 as input, the axial γW -box correction to the tree-level amplitude is given as [3] Here Q 2 = −q 2 > 0 is the spacelike four-momentum square. The normalization factor F H + arises from the local matrix element ⟨H f (p ′ ) J W,V µ H i (p)⟩ = (p + p ′ ) µ F H + , with F H + = 1 for the neutron and √ 2 for the pion decay. Methodology -In the framework of lattice QCD, the hadronic tensor T V A µν in Euclidean spacetime is given by Here the Euclidean momenta P and Q are chosen as with m H the hadron mass. By multiplying µναβ Q α P β to T V A µν , we can extract the function Here I can be written in terms of H V A µν as We can average over the spatial directions for ⃗ Q and have where j n (x) are the spherical Bessel functions. A key ingredient in this approach is that once the Lorentz scalar function µνα0 x α H V A µν is prepared, e.g. from a lattice QCD calculation, one can determine T 3 (Q 0 , Q 2 ) directly.
Putting Eqs. (8) and (6) into Eq. (2) and changing the variables as ⃗ Q = Q 2 cos θ and Q 0 = Q 2 sin θ, we obtain the master formula with For small Q 2 , lattice QCD can determine the function M H (Q 2 ) with lattice discretization errors under control. For large Q 2 , we utilize the operator product expansion There are only four possible local operators at leading twist. (For the pion decay, the hadronic matrix elements for the first three operators vanish.) Multiplying µναβ Q α P β to the relation (11) we obtain where the ellipses remind us that the higher-twist contributions are not included yet. The Wilson coefficient is calculated to four-loop accuracy [37,38] C d (Q 2 ) = n c n a n s , a s = with coefficients Here α s is the strong coupling constant and n f is the number of effective quark flavors. We introduce a momentum-squared scale Q 2 cut that separates the two regimes, and split the integral in Eq. (9) into two parts With Eq. (10) we use the lattice data to determine the integral for Q 2 ≤ Q 2 cut , while with Eq. (12) we use perturbation theory to determine the integral for Lattice setup -We use five lattice QCD gauge ensembles at the physical pion mass, generated by RBC and UKQCD Collaborations using domain wall fermion [39]. The ensemble parameters are shown in Table I. Here 48I and 64I use Iwasaki gauge action in the simulation (denoted as Iwasaki in this work) while the other three ensembles use Iwasaki+DSDR action (denoted as DSDR). We calculate the correlation function We use the wall-source pion interpolating operators φ π 0 and φ † π − , which have a good overlap with the π ground state, and find the ground-state saturation for ∆t ≳ 1 fm. In practice the values of ∆t are chosen conservatively as shown in Table I. For each ensemble, we use the gauge configurations, each separated by at least 10 trajectories. The number of configurations used is listed in Table I.  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 , the number of configurations used, N conf , the number of point-source light-quark propagator generated for each configuration, Nr, and the time separation, ∆t, used for the π ground-state saturation.
There are four types of contractions for γW -box diagrams as shown in Fig. 1. We produce wall-source quark propagators on all time slices. Using the techniques described in Ref. [35] type (A) and (B) diagrams can be calculated with high precision by performing the spacetimetranslation average over L 3 ×T measurements. Under the γ 5 hermitian of the Euclidean quark propagators, one can confirm that type (B) does not contribute to the axial γW -box diagrams. Type (C) diagram is calculated by treating one current as the source and the other as the sink. We calculate point-source propagators at N r random spacetime locations. The values of N r are shown in Table I. These point-source propagators can be placed at either electromagnetic current or weak current. We thus average the type (C) correlation functions over 2 N r measurements. This is similar with the treatment taken by Ref. [40]. We neglect the disconnected contribution (D), which vanishes in the flavor SU (3) limit. Numerical results -In practice, the integral in Eq. (10) can be performed within a range of √ t 2 + ⃗ x 2 ≤ R. Taking the ensemble 64I as an example, M π (Q 2 ) as a function of the integral range R is shown in Fig. 2. We find that for all the momenta Q 2 ∈ [0, 4] GeV 2 , the integral is saturated at large R. We choose the truncation range R 0 ≃ 4 fm, which is a conservative choice for all ensembles listed in Table I. The contributions to the integral from √ t 2 + ⃗ x 2 > R 0 is negligible, indicating that the finite-volume effects are well under control in our calculation. We can further verify this conclusion by a direct comparison using the 24D and 32D data. ) as a function of Q 2 . In the left panel, the lattice results for ensembles 64I, 48I, 32D-fine, 32D and 24D are represented by turquoise, indigo, dark green, red and black bands, respectively. Taking 64I and 24D as examples, the results for type (A) diagram are also plotted. In the right panel, Iwasaki and DSDR results at the continuum limit are shown by the gray and brown bands. The orange curve shows the results from perturbation theory by decoupling the charm quark at 1.6 GeV while the magenta one is compiled using the 4-flavor perturbation theory continuously down to 1 GeV.
The lattice results of M π (Q 2 ) as a function of Q 2 are shown in Fig. 3 together with the perturbative results. Ensemble 24D and 32D have the same pion mass and lattice spacing but different volumes. The good agreement between these two ensembles indicates that the finitevolume effects are smaller than the statistical errors. At Q 2 ≳ 1 GeV 2 , the lattice discretization effects dominate the uncertainties. In the left panel of Fig. 3 an obvious discrepancy is observed at large Q 2 for the lattice results with different lattice spacings.
For the perturbation theory, the Wilson coefficient C d (Q 2 ) is determined using the RunDec package [41], where α s is calculated to four-loop accuracy. At low Q 2 the results still contain large systematic uncertainties due to the lack of higher-loop and higher-twist contributions. In Fig. 3 we show two curves from perturbation theory. One is compiled using 4-flavor theory down to 1 GeV, while the other decouples the charm quark at 1.6 GeV and uses 3-flavor theory for (1 GeV) 2 ≤ Q 2 ≤ (1.6 GeV) 2 . The discrepancy between the two curves suggests an O(14%) systematic effect in the perturbative determination of M H (Q 2 ) at Q 2 ≈ 1 GeV 2 .
Estimate of systematic effects -For ◻ V A,≤ γW the largest uncertainties arise from the lattice discretization effects. Since Iwasaki and DSDR ensembles have different lattice discretizations, we treat them separately. After the linear extrapolation in a 2 , the Iwasaki and DSDR results at continuum limit are shown in the right panel of Fig. 3. Using Q 2 cut = 2 GeV 2 we obtain We take the Iwasaki result as the central value and estimate the residual O(a 4 ) lattice artifacts using the discrepancy between Iwasaki and DSDR. For ◻ V A,> γW the largest uncertainties arise from the higher-loop and higher-twist truncation effects. We estimate the former by comparing the 4-loop and 3-loop results from perturbation theory. For the latter, unfortunately the complete information is not available. Considering the fact that type (A) diagram, which has two currents located at different quarks lines, only contains the higher-twist contributions, we use it to estimate the size of higher twist. At Q 2 cut = 2 GeV 2 we have where the central value is compiled using the 4-flavor theory (See the magenta curve in the right panel of Fig. 3). The first error indicates the higher-loop effects. The second one stands for the higher-twist effects, which are compiled from the integral of Q 2 > Q 2 cut using the type (A) data as input.
Summary of results -After combining the results of ◻ V A,≤ γW from lattice QCD and ◻ V A,> γW from perturbation theory, we obtain the total contribution of ◻ V A γW ◻ V A γW π = 2.830(11) stat (9) PT (24) a (3) FV × 10 −3 = 2.830(11) stat (26) sys × 10 −3 , where the first uncertainty is statistical, and the rest errors account for perturbative truncation and higher-twist effects, lattice discretization effects, and lattice finite-volume effects by comparing the 24D and 32D results. We add these systematic errors in quadrature to obtain the final systematic error.
We now show how our calculation reduces the uncertainty in δ. We adopt Sirlin's parameterization [2] with slight modifications: By separating the axial γW -box part into ◻ V A γW , the remaining contributions are model independent and are given as follows.
• Sirlin's functionḡ arises from a structureindepenent, UV-finite one-loop integral. It accounts for the the infrared contributions involving the vector γW -box and the bremsstrahlung corrections. It contains a 3 ln m p term that cancels the m p -dependence in 3 ln m Z mp . Here m p is the proton mass that appears just as a matter of convention. Numerically, one has αe 2πḡ = 1.051 × 10 −2 [2,47]. •ã g represents the O(α s ) QCD correction to all oneloop diagrams except the axial γW box. The integral inã g is dominated by the high-energy scale Q 2 ≃ m 2 W , where α s is small. As a consequence αe 2πã g ≈ −9.6 × 10 −5 is a small contribution [2,48].
• δ QED HO = 0.0010(3) summarizes the leading-log higher-order QED effects which can be accounted for through the running of α e . The uncertainty assignment follows Ref. [49].
Although the detailed uncertainties forḡ andã g are not given, by power counting the intrinsic precision for the terms in the square brackets (multiplied by αe 2π ) is of the order G F m 2 p ∼ 10 −5 . Combining the ◻ V A γW in Eq. (18), we now obtain which agrees well with the result from chiral perturbation theory. The uncertainty is reduced by a factor of 3 and comes almost entirely from δ QED HO . Therefore, any theoretical improvement in the future will unavoidably require a complete electroweak two-loop analysis. Consequently, the V ud determined from the pion semileptonic decay now reads: V ud = 0.9740(28) exp (1) th .
Conclusion -In this work we perform the first realistic lattice QCD calculation of the γW -box correction to the pion semileptonic decay, ◻ V A γW π . The final result combines the lattice data at low momentum and perturbative calculation at high momentum. We use multiple lattice spacings and volumes at the physical pion mass to control the continuum and infinite-volume limits and obtain ◻ V A γW π with a total error of ∼ 1%. As a result, the uncertainty of the theoretical prediction for the pion semileptonic decay rates is reduced by a factor of 3.
The combined experimental measurement of 14 nuclear superallowed beta decays [1] is 10 times more accurate than the current pion semileptonic decay experiment. On the other hand, the free neutron decay [50,51] leads to a 4.5 times better precision. In these two cases, the nonperturbative, structure-dependent γW -box contribution plays a decisive role. The technique presented in this work can be straightforwardly generalized to a lattice calculation of the nucleon γW -box corrections, which are universial for both free and bound neutron decay. The latter is the key to a precise determination of V ud and a stringent test of CKM unitarity.
Acknowledgments -X.F. and L.C.J. gratefully acknowledge many helpful discussions with our colleagues from the RBC-UKQCD Collaborations. We thank Yan-Qing Ma and Guido Martinelli for inspiring discussions.