Dispersion relation analysis of the radiative corrections to $g_A$ in the neutron $\beta$-decay

We present the first and complete dispersion relation analysis of the inner radiative corrections to the axial coupling constant $g_A$ in the neutron $\beta$-decay. Using experimental inputs from the elastic form factors and the spin-dependent structure function $g_1$, we determine the contribution from the $\gamma W$-box diagram to a precision better than $10^{-4}$. Our calculation indicates that the inner radiative corrections to the Fermi and the Gamow-Teller matrix element in the neutron $\beta$-decay are almost identical, i.e. the ratio $\lambda=g_A/g_V$ is almost unrenormalized. With this result, we predict the bare axial coupling constant to be {$\mathring{g}_A=-1.2754(13)_\mathrm{exp}(2)_\mathrm{RC}$} based on the PDG average $\lambda=-1.2756(13)$


I. INTRODUCTION
The recent emergence of an apparent deficit in the top-row Cabibbo-Kobayashi-Maskawa (CKM) matrix unitarity [1], |V ud | 2 + |V us | 2 + |V ub | 2 = 0.9985 (5) (1) has triggered a renewed interest in precise experimental studies of various β-decay processes giving access to |V ud |. Superallowed 0 + − 0 + nuclear decays have long been regarded as the best avenue for such purpose. Recent works, however, pointed out that current theory uncertainties in the nuclear-structure corrections may have been significantly underestimated [2][3][4][5]. Reducing these uncertainties requires novel ab-initio nuclear theory calculations that are not yet available. As a consequence, the role of alternative channels such as the β-decays of the free neutron, mirror nuclei and pion becomes increasingly important. With the future improvement in the experimental precision, these other β decay processes will offer competitive determinations of |V ud | and complementary sensitivity to possible beyond standard model (BSM) signals. Free neutron β-decay currently provides the second best determination of |V ud | through the following master formula [6,7]: where the uncertainty in the numerator arises from the Standard Model (SM) theory input. The two required experimental inputs are the neutron lifetime τ n , and the decay parameter λ ≡ g A /g V which is the ratio between the neutron axial and vector coupling constant. This parameter is renormalized by electroweak radiative corrections (RCs), and these latter are the primary focus of this article. The parameter λ can be measured either via the P-even correlation p e · p ν (the a coefficient), or the P-oddê s · p e (the A coefficient) andê s · p ν (the B coefficient) ones, withê s the unit 3-vector along the neutron polarization. The current best measurement reported by the PERKEO III collaboration λ = −1.27641(45) stat (33) sys , with a 0.04% precision [8] (in this paper we pick the sign convention λ < 0, also adopted in the Particle Data Group (PDG) review). However, the current PDG average reads λ = −1.2756 (13) [1], where the much larger uncertainty is due to a scale factor of 2.6 that accounts for the large discrepancy between the results before [9][10][11][12] and after 2002 [8,[13][14][15][16] (see Ref. [17] for more discussions). Future improvements are expected from the Nab [18] and PERC [19,20] collaborations, both aiming at an accuracy level of 10 −4 .
The exact value of λ not only serves for extracting V ud , but is also interesting in itself. The "bare" (i.e. without electroweak corrections) neutron axial couplingg A is one of the simplest hadronic matrix elements and has received much attention. Unlike its vector counterpartg V which remains non-renormalized due to the conserved vector current (CVC), the bare axial coupling is not protected and must be calculated, e.g. using lattice Quantum Chromodynamics (QCD) . The most recent FLAG average [43] reads: N f = 2 + 1 + 1 :g A = −1.251 (33) N f = 2 + 1 :g A = −1.254 (16)(30) N f = 2 :g A = −1.278 (86) , but individual calculations have achieved higher precision. For instance, Ref. [40] reported a percent-level determination ofg A = −1.271(10)(7) using an unconventional method inspired by the Feynman-Hellmann theorem, and follow-up works are aiming for sub-percent precision [42]. Such a rapid development makes λ a powerful tool for searching for redBSM physics. By comparing first-principles calculations ofg A to the experimental results for λ one thereby constraints the strength of possible BSM contributions that could modify g A , in particular the right-handed currents [44][45][46][47].
When the lattice precision reaches 10 −3 , a valid comparison betweeng A and λ will require precise O(α em /π) RCs that bringg A /g V to g A /g V . In particular, we need to deal with sizable hadronic uncertainties originating from the γW -box diagram. The latter can be written as a Q 2 -integral, and performing operator product expansion (OPE) it is well-known that the corrections to g V and g A coming from large Q 2 (which carries a large electroweak logarithm) are the same. Therefore, it was believed that the difference betweeng A and λ is numerically small [48][49][50][51]. However, for a long time there was no serious attempt to understand the RC to g A from the low-Q 2 part of the integral which is also of the order 10 −3 , comparable to the high-Q 2 contribution. It includes the elastic contributions that are fixed by the nucleon form factors, as well as the inelastic contributions that are governed by non-perturbative QCD. The first attempt for a complete analysis was performed recently in Refs. [52,53]. Anticipating the results of this work, we found that those Refs. originally contained algebraic mistakes in the computation of the elastic contribution, invalidating their numerical results. These mistakes were later corrected in the published version of Ref. [52]. Additionally, the inelastic contribution residing at low Q 2 was obtained based on a holographic QCD model, following Ref. [54] where the RC to g V was addressed. The model-dependent nature of this approach makes a rigorous estimation of the theoretical uncertainty complicated.
In this article we improve on both points. We perform a novel analysis of the RC to g A based on the dispersion relation (DR) approach. It is a powerful tool which has proved successful in the treatment of the RCs to the Fermi amplitude in the free neutron and superallowed β-decays [2][3][4]7]. In this formalism, the γW -box diagram is expressed as a dispersion integral over structure functions that are directly or indirectly related to experimental data. This ultimately allows for a fully data-driven analysis of this RC. For g A , the required input relies on the spin-dependent structure functions g 1 and g 2 , well-studied quantities in deep inelastic scattering (DIS) experiments. We utilize highprecision world data on g 1,2 to evaluate the dispersion integral, and fix the forward γW -box diagram correction to g A to an unprecedented precision better than 10 −4 . We observe that the RCs to g V and g A are numerically very close and largely cancel in the ratio, which practically removes any distinction betweeng A and λ down to 2 × 10 −4 .
The contents in this paper are arranged as follows. In Section II we define our notation and introduce the starting point for the discussion of the RC. We introduce the γW -box diagram in Section III, and derive its dispersive representation in Section IV. The elastic (Born) and inelastic contributions to the box diagram are computed in Section V and VI respectively. The final results and discussions are presented in Section VII.

II. GENERAL FRAMEWORK
We start by defining the hadronic currents relevant to the β-decay of the free neutron: Their single-nucleon matrix elements are given by: where N = p, n and M = (M p + M n )/2 ≈ 939 MeV. All the form factors above are functions of −(p i − p f ) 2 . We may also define the isospin combinations F S 1,2 = F p 1,2 + F n 1,2 and F V 1,2 = F p 1,2 − F n 1,2 . The values of the vector and the axial charged weak form factors at zero momentum transfer define the "bare" vector and axial coupling constants: F W 1 (0) =g V , G A (0) =g A , which represent the Fermi and Gamow-Teller matrix element in neutron β-decay respectively. In particular,g V = 1 from isospin symmetry, and the correction due to the strong isospinbreaking effects is negligible due to the Behrends-Sirlin-Ademollo-Gatto theorem [55,56]. On the other hand,g A is not protected by any exact symmetry. We do not include the isospin-breaking correction tog A separately because it is already included in the respective first-principles calculations.
The nucleon mass difference ∆ = M n − M p ≈ 1.3 MeV and the electron mass m e ≈ 0.511 MeV are much smaller than M . Therefore, the tree-level amplitude of the decay process n(p n ) → p(p p )e(p e )ν e (p ν ) is given by: where G F = 1.1663787(6) × 10 −5 GeV −2 is the Fermi constant measured from the muon decay [1], p = (p n + p p )/2 is the average nucleon momentum, and L λ =ū(p e )γ λ (1 − γ 5 )v(p ν ) is the lepton piece. The recoil corrections scale as ∆/M ∼ 10 −3 , which are small but important in precision physics. They were studied in detail with both conventional methods and effective field theory (EFT) [57][58][59][60][61][62], and will not be discussed here. RCs of the order O(α em /π) must be included for a precise extraction of the weak coupling parameters. In the usual nomenclature, they are divided into the "outer" and "inner" corrections The former is a function of {E e , m e } calculable within Quantum Electrodynamics (QED) and independent of details of strong interaction. The latter is instead a constant in E e but depends on details of the hadronic structure. The squared amplitude for the decay of a polarized neutron (to unpolarized final states) after the inclusion of the O(α em /π) RCs reads: with Here, E e = p e · p/M is the electron energy, E m = (M 2 n − M 2 p + m 2 e )/(2M n ) ≈ M n − M p is the electron end-point energy, and β = 1 − m 2 e M 2 /(p e · p) 2 is the electron speed, all in the nucleon's rest frame. With these notations, the functions δ (1,2) describe the outer corrections [48,63]: whereas F (β) ≈ 1 + α em π/β is the Fermi's function that incorporates the Coulomb interaction between the final-state proton and the electron [64]. The function δ (1) is also known as Sirlin's function g(E e , E m ). The axial to vector coupling constants' ratio, parameter λ = g A /g V , is understood as fully renormalized by the inner RCs. In near-degenerate semileptonic β-decay processes, the inner RCs are most conveniently studied in Sirlin's representation [65] (see also Refs. [66][67][68] for a detailed account). In this formalism, most of the O(α em /π) electroweak RCs are either exactly known from current algebra, or give rise to the outer corrections in Eq.(9) and the Fermi's function. As a result, the renormalized vector and axial coupling constants read: whereã g is a pQCD correction factor and δ QED HO summarizes the leading-log higher-order QED effects [69,70]. One observes that the fractional corrections tog V andg A are mostly identical and cancel in the ratio g A /g V . The only exceptions are the constants V γW and A γW that describe the inner RCs originated from the γW -box diagrams (see Fig.1), which are the focus of this paper 1 . With the above, we obtain:

III. γW -BOX DIAGRAM
The γW -box correction is of the natural size α em /π ∼ 10 −3 . Taking into account recoil corrections ∼ ∆/M, m e /M on top of the overall α em /π factor would bring us to accounting for effects in the 10 −6 range that exceed the precision goal by two or three orders of magnitude. This defines the level of the detalization that is needed in our analysis. We will consistently set ∆ = m e = 0 throughout the calculation below, as well as the proton recoil. This approximation also leads to the neglect of the pion pole due to the partially-conserved axial current (PCAC) hypothesis: the pion pole contribution, when contracted with the lepton tensor, results in lepton mass terms which, as stated are neglected. This precision level is supported by the fact that the lowest hadronic state is separated by the pion mass ∼ 140 MeV which is about hundred times larger than ∆. Notice however that this approximation may not be as safe for nuclear β-decay where the available Q-values may be as large as [15][16][17][18][19][20] MeV which are comparable to the energy level of nuclear excitations.
The part of the γW -box diagram amplitude that contributes to the inner correction must involve an antisymmetric tensor that stems from the lepton spinor structure. It reads: where with 0123 = −1 in our convention. The forward generalized Compton tensor describing the W + n → γp process, is defined as: To extract V γW and A γW , we use following identities: where the spin vector S µ is fixed by S 2 = −M 2 and S · p = 0, we obtain: In what follows, we useg A ≈ λ = −1.2756 as a normalization in the second expression. Since A γW ∼ 10 −3 , the error induced by the ambiguity ofg A is well below our precision goal.
Only those components in T µν γW that contain an antisymmetric tensor contribute to Eq. (16). These are The spin-independent, parity-violating amplitude T 3 and spin-dependent, parity-conserving amplitudes S 1,2 are functions of two invariants, ν = (p · q)/M and Q 2 = −q µ q µ . Plugging Eq. (17) into Eq. (16) gives: where we have used the following identities: that hold for any Lorentz scalar function F (ν, Q 2 ).
To evaluate the loop integrals, we need to discuss the symmetry properties of the Compton amplitudes. We start by considering the isospin structure of amplitudes T i , S i . Electromagnetic interaction does not conserve isospin and contains both isoscalar (I = 0) and isovector (I = 1) components. Therefore, each amplitude A = T i , S i can be decomposed into components contributed by the isoscalar and isovector electromagnetic current respectively: The two isospin amplitudes have a different behavior under ν → −ν: with ξ i . It is easy to show that ξ (0) i = −1 for T 3 and S 1,2 , so only the I = 0 component of these amplitudes survives in the integrals in Eq. (18).

IV. DISPERSION REPRESENTATION OF THE FORWARD COMPTON AMPLITUDES
Forward Compton amplitudes have singularities along the real axis ν: poles due to a single nucleon intermediate state in the s− and u-channels at ν = ±ν B = ±Q 2 /(2M ), respectively, and unitarity cuts at ν ≥ ν π and ν ≤ −ν π where ν π = (2M m π + m 2 π + Q 2 )/(2M ) is the pion production threshold (see Fig.2). The discontinuity of the forward Compton tensor in the s−channel (i.e. ν ≥ ν B ) is given by the generalization of the on-shell hadronic tensor to the γW -interference: where The structure functions F 3 and g 1,2 can be decomposed similarly to I = 0, 1 components just like Eq. (20). According to the crossing behavior established earlier and noticing that they cannot diverge faster than ν when ν → ∞, the amplitudes entering Eq. (18) have the following dispersion representation: where the ν-integration is extended down to 0 to include the Born contribution. Notice that in the last line we have slightly modified the DR of S (0) 2 using the Burkhardt-Cottingham (BC) sum rule [71]: where x = Q 2 /(2M ν ) is the Bjorken variable. This sum rule is a superconvergence relation and is expected to hold at all Q 2 . The benefit of this treatment will become apparent in the later section. Substituting Eq.(24) into Eq. (18) and using the following Wick rotation formula [72], we can integrate the variable u analytically to obtain our final dispersive representation of δg γW V,A as follows: with r = 1 + 4M 2 x 2 /Q 2 . As a useful crosscheck, the two-photon exchange correction to the hyperfine splitting in ordinary and muonic atoms is expressed through analogous two-fold integrals over electromagnetic spin structure functions [73,74]. The quantity V γW , relevant for the extraction of V ud from superallowed β-decays is well-studied within the dispersive approach [2,3,7,75] and is not addressed here. The main obstacle in those studies is the absence of direct experimental data of the structure function F 3 . This forces one to either rely on data of F 3 from a different isospin channel (whose relation to F (0) 3 contains a residual model-dependence), or from indirect lattice QCD data. On the other hand, despite having received much less attention, a high-precision dispersive analysis of A γW is in fact much more robust because it depends on the parity-conserving, spin-dependent structure functions g 1,2 . The isospin symmetry unambiguously relates them to g N 1,2 measured in ordinary DIS: the latter are defined via Therefore, it is possible to perform a fully data-driven analysis of A γW without introducing further model-dependence at low Q 2 . We will perform such an analysis in the sections below. Following Sirlin's notation [48], we express our result as where d B , d 1 and d 2 represent the elastic (Born) contribution, the inelastic contributions from g 2 respectively, which we will evaluate separately in the following sections.

V. ELASTIC (BORN) CONTRIBUTION
Substituting X = p into Eq.(23) and using the elastic form factors defined in Eq. (5) give us the Born contribution to the spin structure functions needed for the evaluation of the box correction: where G E ≡ F 1 − τ F 2 and G M ≡ F 1 + F 2 are the usual electric and magnetic Sachs form factors defined for both the electromagnetic and charged weak form factors, with τ = Q 2 /(4M 2 ). All the form factors above are functions of Q 2 that drop at high Q 2 , therefore one can neglect the Q 2 dependence of the W -boxon propagator ∼ Q 2 /M 2 W . With this the Born contribution reads, where r B ≡ r| x=1 . We notice that F W 1,2 = F V 1,2 by isospin symmetry, which means Eq.(32) is fully determined by the four nucleon electromagnetic form factors: {G p E , G p M , G n E , G n M }. Different parameterizations of these form factors [76][77][78][79][80][81] all give consistent results within their respective error bars. In particular, the parametrization of Ref. [79] leads to In particular, the central values of the contribution from g We observe that the latter is much smaller, which turns out to also be the case for the inelastic contributions.
We pause here to comment on the Born contribution before moving on to the inelastic contributions. One may also try to derive it by calculating the Compton amplitudes from the first two Feynman diagrams in Fig.3 (the "pole diagrams") using the form factors in Eq.(5) as effective vertex functions, and then plugging them into Eq. (18). The pole diagrams give: We split each expression into two term, where the first term contains a singularity at ν 2 = ν 2 B and vanishes as 1/ν when ν → ∞, while the second term is regular and diverges as ν when ν → ∞. It is easy to see that, retaining only the first term leads again to Eq.(32), apart from a numerically small difference originating from our accounting for the BC sum rule in the DRs, effectively redefining the g 2 contribution to d B . The regular terms in Eq. (34) lead to an extra small deviation from Eq. (32). The origin of this deviation lies in the Gerasimov-Drell-Hearn sum rule [82,83] and its extension to finite Q 2 [84] which relate the regular low-energy term to an integral over the inelastic part of g 1 .
This discussion simply means that the definition of the "elastic" contribution is not exactly the same in the diagrammatic and the dispersive representation. Of course, if we were able to calculate the full (i.e. pole + seagull) T µν γW exactly at all values of {ν, Q 2 } with the diagrammatic approach, then the outcome must be identical to the DR analysis. But since this is impossible, the dispersive representation provides a much better starting point. We want to also point out that Refs. [52,53] attempted to calculate the Born contribution from the pole diagrams in Fig.3 (let us call it d B ). 2 Should all terms be retained, they would have obtained the following result: which only differs from the DR's definition of d B in Eq.(32) by a numerically small term. For the benefit of interested readers, we also provide in Appendix A an alternative derivation of Eq.(35) without making use of the invariant amplitudes.

VI. INELASTIC CONTRIBUTIONS
Extensive measurements of the structure function g N 1 were carried out in SLAC [85][86][87], CERN [88][89][90][91], DESY [92] and JLab [93][94][95][96][97]. In particular, we utilize the results from the EG1b experiment at JLab that measured the g p 1 [97] and g n 1 [96] in a wide range of {x, Q 2 }, from which the moments Γ N i (i = 1, 3, 5) were computed in bins of Q 2 , common for p and n, from 0.05 GeV 2 to 3.5 GeV 2 . Full results are available in the supplementary material of each respective paper 3 . Data on g 2 are generally scarce [94,97,[99][100][101][102] and insufficient for a fully data-based analysis. Fortunately, its contribution is generally expected to be small due to the BC sum rule that forces the first moment of g 2 to vanish identically when accounting for elastic and inelastic contributions.
Since the data do not extend to an arbitrarily large Q 2 needed to evaluate the integrals, we make use of perturbative QCD results which are well under control theoretically above a separation scale Q 2 0 , while directly using the experimental data that contains both perturbative and nonperturbative physics below that scale. Following our earlier works on the vector RC [2,3] we take Q 2 0 = 2 GeV 2 . For the vector RC case, not only does Q 2 0 = 2 GeV 2 mark the onset of pQCD regime, it also corresponds to the scale, below which the quality of data deteriorates severely leading to some sensitivity to Q 2 0 . In the case of the axial RC, this scale lies well within the region covered by data, hence shifting it to a slightly higher value (not lower because pQCD description starts to break down) does not change the result. 2 In the earlier versions of these Refs., the author made some algebraic mistakes when dealing with the symmetric loop integral of the form d 4 qqαq β F (ν, Q 2 ) (i.e. Eq.(19)). As a consequence, an incorrect analytic formula which gave an unexpectedly large value of d B = 2.64(3) was obtained. The published version of Ref. [52] corrected these mistakes, but retained only the term proportional to (5 + 4r B )/3. 3 There is a more recent measurement of g p 1 from JLab at low Q 2 [98], but unfortunately it does not measure g n 1 simultaneously.
A. Contribution of g where x π = Q 2 / (M + m π ) 2 − M 2 + Q 2 is the pion production threshold. Notice that the definition above excludes the elastic contribution at x = 1, which is a general convention adopted by most of the experimental papers. The first moment Γ p−n 1 ≡ Γ p 1 − Γ n 1 is of particular interest because it satisfies the polarized Bjorken sum rule [103,104] at Q 2 → ∞. However, at large but finite Q 2 it receives a number of corrections [105]: here the subscript "th" denotes the theory prediction (at large Q 2 ). The first term at the right hand side is the Bjorken sum rule with a pQCD correction factor 4 , while the second term summarizes the higher-twist (HT) effects starting from twist-four. The pQCD correction factor is written as: where α s is the running strong coupling constant in the MS scheme, while the coefficients {c n } are calculated at present to n = 4 [106,107]:c with n f the number of active quark flavors, and we refer the reader to Refs. [2,3] for full detail of the pQCD contribution and relevant discussions and references. In the meantime, only the twist-four term among all the HT corrections needs to be included for our precision goal. There are several recent determinations of the coefficient µ p−n 4 [108][109][110] that are largely consistent with each other. In this work we quote the value µ p−n 4 = (−0.047±0.020)M 2 in Ref. [109]. We find that at 2 GeV 2 , the inclusion of the twist-four correction reduces the size of Γ p−n 1,th by about 13%, but its total contribution to d 1 through the integral at Q 2 > 2 GeV 2 is only about 1%. Coming back to our problem, we write When Q 2 → ∞ the functionΓ p−n 1 reduces to Γ p−n 1 , but at low Q 2 the two are not identical due to target mass corrections ∼ M 2 /Q 2 contained in the factor f (x, Q 2 ) = 4(5 + 4r)/(9(1 + r) 2 ), where r = 1 + 4x 2 M 2 /Q 2 We use the following strategy to reconstruct the fullΓ p−n 1 at low Q 2 from data. We fit f (x, Q 2 ) as function of x at fixed Q 2 as and obtain the three fitting parameters a(Q 2 ), b(Q 2 ) and c(Q 2 ) by first dividing 0 < x < x π into, say, 1000 equal intervals, evaluating the 1000 respective discrete values for f (x, Q 2 ), and performing a three-parameter fit with these reconstructed from the EG1b experiment [96,97] versus the theory prediction at large Q 2 with (red band) and without (brown curve) the twist-four correction. discrete points using, e.g. Mathematica. We find that this procedure allows for a very precise reproduction of the entire curve at 0 < x < x π , where the difference between the original and the fitted curve is negligible for all practical purposes. Contrarily, a simple Taylor expansion in powers of is only applicable at high Q 2 and low x but significantly deviates for larger x. It breaks down completely at Q 2 = 0 where it becomes divergent, whereas f (x, 0) remains finite. As an illustration, the two approximated expressions evaluated at a representative value of Q 2 = 1 GeV 2 are plotted together with the analytic form in the left panel of Fig.4, and we clearly see that Eq.(41) nicely reproduces the latter for the full range of x. We thus reconstruct the fullΓ p−n 1 in each bin of Q 2 in terms of the lowest Mellin moments, where Γ p−n i (i = 1, 3, 5) are obtained from Refs. [96,97]. We find that the difference betweenΓ p−n 1 and Γ p−n 1 (i.e., the effect of higher moments) does not exceed 3.5% for Q 2 = 2 GeV 2 , which implies a negligible difference in the integral at Q 2 > 2 GeV 2 . Therefore we will not distinguish between the two above 2 GeV 2 . In contrast, the higher-twist correction due to µ p−n 4 reaches 13% and needs to be kept along. The right panel in Fig.4 shows the reconstructedΓ p−n 1 data points versus the large-Q 2 theory prediction using Eq. (37). We find that the theory and experiment match well at Q 2 > 2 GeV 2 (observe how the twist-four correction is needed to reconcile the two), which justifies our choice of Q 2 0 = 2 GeV 2 as the separation scale between the perturbative and non-perturbative regime. We therefore evaluate d 1 separately in these two regions. At Q 2 < Q 2 0 , we fit three curves that correspond to the upper bounds, central values and lower bounds of the discrete data points respectively, and evaluate the Q 2 -integral and its uncertainty by integrating these three curves. Since the uncertainties of the data points are mainly systematics, this prescription takes into account the possible positive correlation effects. The resulting uncertainty is thus a conservative one; it is likely that it can further be reduced, but this would require a dedicated study of the systematic uncertainties of the data which lies beyond the scope of the present work. Meanwhile, at Q 2 > Q 2 0 we evaluate the integral using the theory prediction in Eq. (37). The results are as follows: where we neglected the uncertainty associated with the leading twist contribution compared to those coming from the data at low Q 2 and the HT correction (i.e. the coefficient µ p−n 4 ) at high Q 2 . The total contribution of g 1 reads, We note that the integral below Q 2 ≤ 0.05 GeV 2 not covered by the data but making part of d < 1 is controlled by the isovector GDH sum rule [82,83], dΓ p−n 1 (Q 2 = 0)/dQ 2 = (κ 2 n − κ 2 p )/8M 2 , with κ p,n denoting the proton's (neutron's) anomalous magnetic moment, respectively. Connecting the GDH-fixed value at Q 2 = 0 to the lowest data point produces a negligible d low Q 2 1 0.003 contribution which is safely accommodated within the uncertainty.
In fact, we have already implemented this sum rule in the derivation of the DR of S (0) 2 . We emphasize the importance of this procedure for a reliable estimate of the contribution of g (0) 2 to A γW : since experimental data typically only cover the inelastic region, enforcing an exact vanishing of the first moment of g (0) 2 while operating with phenomenological parametrizations of different pieces can be a delicate matter. The explicit use of the BC sum rule thus precludes any numerically significant mistake caused by an imperfection of these parametrizations. As a result, the dispersive representation of A γW only contains higher moments of g 2 , in which the non-perturbative physics at small x is suppressed. Additionally, since every extra power of x 2 is accompanied by 1/Q 2 , the contribution of g (0) 2 bears no large logarithms. Due to the smallness of the g 2 contribution, we opt for an approximate treatment and include this result in the estimate of the systematic uncertainty.
The inelastic contribution coming from g 2 reads, where the isospin relation in Eq. (28) is used. Rather than relying on data on g N 2 , we decompose g 2 into twist-two and twist-three (and higher) components [1], g N 2 = g N 2,tw2 + g N 2,tw3+ , and use the Wandzura-Wilczek relation [111] for the former, Notice that the relation above should be understood to not contain the elastic contribution at x = 1, because otherwise one could take x π < x < 1 at both sides, and then the left hand side and the first term at the right hand side would vanish but the second term at the right hand side would not, which is a contradiction. With this in mind, we multiply both sides by x i−1 and integrate them at 0 < x < x π to obtain: where Γ N n is defined in Eq. (36). Therefore, we may use the available information of g 1 to evaluate the twist-two contribution to d 2 .
Again, we discuss the integral at large and small Q 2 separately. For Q 2 > Q 2 0 , we keep the leading term, where we have used Eq.(48) and set r → 1. The Q 2 dependence of the W propagator can also be safely neglected as the integral converges. There is no simple sum rule for Γ p−n 3 at large Q 2 , but we may adopt a naïve valence quark picture that assumes each valence quark carries 1/3 of the nucleon's momentum. This gives: which automatically reproduces the free polarized Bjorken sum rule. This naïve picture predicts Γ p−n 3 ≈ 0.11Γ p−n 1 , which we may check against the experimental data: At Q 2 ≈ 3.4 GeV 2 , Refs. [96,97] give Γ p−n by some 30%, an acceptable uncertainty given our precision goal. With the above, we obtain: Next we turn to the small-Q 2 region where accounting for the leading twist may not be sufficient. As before, the leading twist contribution is reconstructed by using the WW relation and the x-integral related to measured moments of g 1 by performing a two parameter fit of the kinematical function f (x, Q 2 ) for fixed Q 2 (since the first moment is removed by the BC sum rule, a two-parameter fit is already sufficient for our precision goal).
Using Eq. (48), we obtain With Γ p−n i taken from Refs. [96,97], this gives, where we assigned a conservative 100% uncertainty to the entire contribution.
To quantify higher twist contributions we recall the definition of the "color polarizability" [112,113] of which we only consider the inelastic partd 2 (Q 2 ) coming from the interval 0 ≤ x ≤ x π [114] since the elastic part is already taken into account. In terms of this polarizability and neglecting higher moments, we obtain for the contribution of twist-three and higher, For numerical estimates, we rely on the recent analysis of generalized spin polarizabilities of the nucleon in baryon chiral effective theory [114]. We obtain, assigning a conservative 100% uncertainty, Combining the various pieces we finally arrive at as our estimate of the total inelastic contribution from g 2 . We observe that it is two orders of magnitude smaller than d 1 , following the same hierarchy as in d B .
The inelastic contribution in our DR analysis then reads d 1 + d 2 = 2.19(4) data (1) HT (3) g2 . This is to be compared with 2.31(9) from Refs. [52,53], and we see that the two do not quite agree within error bars. While being identical in the large-Q 2 treatment, their estimation of the low-Q 2 contribution is largely model-based, raising questions about the reliability of the uncertainty. In contrast, in our treatment the low-Q 2 contribution is completely fixed by experimental data without any further assumption, apart from the small g 2 correction for which only few assumption were made, which will become testable as soon as new, higher-quality low-Q 2 data for g 2 will become available.

VII. FINAL DISCUSSIONS
Collecting all the results from Sec.V and VI gives: where the uncertainties come from the elastic form factors, the low-Q 2 g 1 data, the HT-correction to g 1 and g 2 , respectively. We compare this to our recent update of V γW using indirect lattice inputs: V γW = 3.83(11) × 10 −3 [7]. These two numbers are very close to each other, and in fact their difference is consistent with zero: Using Eqs. (11), (61) and the PDG average λ = −1.2756 (13) [1], we obtain: which is consistent with the result from the current best lattice QCD determination. Our result indicates that there is no practical distinction between λ andg A , unless the experimental precision of the former and the lattice precision of the latter have reached 3 × 10 −4 or better. We wrap up with some discussions of the future prospects. Within the same DR framework, a much better precision is achieved for A γW than for V γW thanks to the existence of high-quality data of the structure function g 1 at Q 2 < 2 GeV 2 . On the other hand, the precision of V γW is limited by the low-quality data of the structure function F 3 from neutrino (antineutrino)-nucleus scattering experiments in the 80s [115,116]. Better-quality data may come from the Deep Underground Neutrino Experiment (DUNE) in the next decade [117,118].
It was pointed out that a direct lattice QCD calculation of V γW is a promising way to proceed at the present stage [119]. Several exploratory calculations of mesonic γW -box diagrams have shown great success [67,120] and the same technology is directly applicable to nucleon. At present, no such direct calculation on the nucleon is available yet. A more involved comparison of the DR result for V γW with the lattice computation of the respective quantity on the pion, amended with further phenomenological ingredients shows a nearly perfect agreement [7]. Even with this reassuring agreement, it is not unthinkable of that a direct lattice calculation could still disagree with the phenomenological, DR-based evaluation. Examples of such an unexpected disagreement are the pion-nucleon sigma term σ πN (see Ref. [43] and references therein) and, more recently, the hadronic vacuum polarization contribution to g µ −2 [121,122]. They show that even carefully-performed first-principles calculations or fully data-driven analysis may still contain unknown, previously unanticipated systematic effects which may seriously affect the implications of the corresponding precision experiments. Given these precedents, it is always useful to cross-check the lattice calculations with alternative methods. Our new result of A γW is perfectly up to this task as it is a solid phenomenological determination with the uncertainty very well under control.
Further effort from the DR side should be dedicated to RCs to the GT strength in nuclear mirror decays where a recent study [123] revealed inconsistencies in the previous analyses. Removing these inconsistencies led to a better agreement for the V ud extracted across mirror and superallowed nuclear decays, as well as neutron decay. However, Ref. [123] only partially accounted for the γW -box contribution. The dispersion formulation of the A γW correction developed in this work can be directly applied to mirror systems. Following Refs. [3,4], nuclear modifications of the universal free-neutron γW -box correction can be computed, and we defer this task to future work.
We start by computing T γW µν from the direct (D) and crossed (C) pole diagrams in Fig.3. Using the elastic form factors as effective vertices, we obtain: are matrices in the Dirac space. We can get rid of the nucleon spinorsū s (p), u s (p) in the expression of T γW,B µν using the following trick. First, we recall that any Dirac structure Γ can be decomposed in terms of standard Dirac basis 1, γ 5 , γ α , γ α γ 5 , σ αβ using the following identity: Next, we have the following identities when a Dirac basis is sandwiched betweenū s (p) and u s (p): u s (p)u s (p) = 2M u s (p)γ 5 u s (p) = 0 u s (p)γ α u s (p) = 2p ᾱ u s (p)γ α γ 5 u s (p) = 2S ᾱ u s (p)σ αβ u s (p) = 2 M αβρσ p ρ S σ .
where i = D, C. The right hand side is free from the nucleon spinors. The trace of the Dirac matrices can be performed using various Mathematica packages, so we do not display the explicit results here. We then plug our expression of T γW,B µν , now free from nucleon spinors, into Eq.(16) and use Eq. (19) to simplify the integrals. Finally, we perform the Wick rotation using Eq.(26) and integrate out the variable u analytically. This brings us exactly to Eq. (35).