Effective field theory for radiative corrections to charged-current processes I: Vector coupling

We study radiative corrections to low-energy charged-current processes involving nucleons, such as neutron beta decay and (anti)neutrino-nucleon scattering within a top-down effective-field-theory approach. We first match the Standard Model to the low-energy effective theory valid below the weak scale and, using renormalization group equations with anomalous dimensions of $\mathcal{O}(\alpha, \alpha \alpha_s, \alpha^2)$, evolve the resulting effective coupling down to the hadronic scale. Here, we first match to heavy-baryon chiral perturbation theory and subsequently, below the pion-mass scale, to a pionless effective theory, evolving the effective vector coupling with anomalous dimensions of $\mathcal{O}(\alpha, \alpha^2)$ all the way down to the scale of the electron mass, relevant for beta decays. We thus provide a new evaluation of the ``inner"radiative corrections to the vector coupling constant and to the neutron decay rate, discussing differences with the previous literature. Using our new result for the radiative corrections, we update the extraction of the Cabibbo-Kobayashi-Maskawa matrix element $V_{ud}$ from the neutron decay.


Introduction
Low-energy processes mediated by the charged-current (CC) weak interaction provide promising ways to test the Standard Model (SM) and indirectly search for new physics, provided sufficiently high experimental and theoretical precision can be achieved. In recent years, there has been a resurgence of interest in beta decays and CC neutrino scattering on nuclei. On the one hand, the study of beta decays at the sub-permille level provides a unique window into possible new physics at the multi-TeV scale. Recent analyses [1][2][3][4][5][6][7][8] have uncovered a 3σ tension with the Standard Model interpretation of these processes in terms of the unitary Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix [7,9]. Moreover, global analyses of beta decay observables [10,11], including decay correlations, offer unique ways to probe nonstandard CC interactions with Lorentz structures different from the SM "V − A". On the other hand, the interest in the CC neutrino scattering process stems primarily from neutrino oscillation experiments [12][13][14][15][16][17][18][19][20][21][22][23][24], as precise theoretical predictions are needed to calibrate the neutrino fluxes and reconstruct the neutrino energy [25][26][27][28][29][30]. In what follows, we will focus on beta decays (n → peν e ) but our results, based on a low-energy effective theory, apply to neutrino scattering processes such asν e p → e + n and ν e n → ep • As a first application of the new framework, using the dispersive input from Refs. [1][2][3][4][5][6] as compiled in Ref. [8], we evaluate the combination of LECs that determine the "inner" corrections to the Fermi transition effective coupling g V . We combine this with the known O(α) radiative corrections to the matrix element in HBChPT [44,46]. We further resum the Coulomb-enhanced terms scaling as (πα/β) n (β ≡ p e /E e ) as well as subleading α/π(πα/β) n terms, in the nonrelativistic Fermi function, which is the natural quantity appearing in / πEFT. In practice, this amounts to replacing the relativistic Fermi function (which contains large logarithms ∼ α 2 ln(R p p e ), with the proton radius R p ∼ 1/Λ χ ) with its nonrelativistic counterpart. Finally, we study the impact on the extraction of V ud from neutron decay. For the total corrections to the neutron decay rate, we find a result that is one σ above the previous results, pointing to a correspondingly smaller value for V ud .
The paper is organized as follows. In Section 2, we provide a high-level summary of the results worked out in the rest of the paper, highlighting the connections to and differences from the previous literature. Following a top-down approach, we perform a multi-step matching to connect electroweak physics with neutron and nuclear decays. The first step, connecting the full Standard Model to the LEFT, is presented in Section 3. The second step, connecting the LEFT to HBChPT, is presented in Section 4. The resulting effective vector coupling g V (µ χ ∼ m N ) at the matching scale µ χ ∼ m N and its evolution to the scale of the decay, µ χ ∼ E 0 , is discussed in Section 5. In Section 6, we discuss the implications for neutron decay and the determination of V ud and comment on the relation to superallowed 0 + → 0 + transitions. Conclusions and outlook are presented in Section 7. Appendix A contains details about electric charge renormalization and running in the LEFT and Chiral Perturbation Theory. Appendix B discusses the factorization of the nonrelativistic Fermi function in nonrelativistic QED, while Appendix C contains details on the extraction of the O(α 2 ) anomalous dimension in LEFT and HBChPT// πEFT.

Statement of the problem and results
Neutron decay is a low-energy process characterized by the energy scales of the neutron-proton mass difference, m n − m p ≈ 1.3 MeV, and the electron mass m e ≈ 511 keV. These scales, which we denote by q ext , are much smaller than the pion mass, m π ≈ 137 MeV, the nucleon mass, m N ≈ 939 MeV, and the W boson mass, M W ≈ 80 GeV. The existence of widely separated mass scales makes the process amenable to a description based on EFTs. In this work, we systematically implement EFT methods to study lowenergy charged-current processes such as neutron decay. We first integrate out the heavy particles (W , Z, h, t) and match the full Standard Model onto the so-called LEFT. Subsequently, we integrate out the scale of the nucleon mass, by matching the LEFT onto HBChPT [50]. We finally integrate out physics at the scale of the pion mass, following [44], by matching HBChPT onto / πEFT. The neutron decay rate is thus organized in an expansion in several small parameters (besides G F q 2 ext , which sets the overall scale): the electromagnetic coupling constant α, ϵ recoil = q ext /m N , which describes small kinematic corrections, ϵ / π = q ext /m π , which captures the radiative pion contributions, and the HBChPT expansion parameter ϵ χ = m π /Λ χ with the scale Λ χ = 4πF π ≈ 1 GeV.
The neutron decay rate is most conveniently computed starting from the / πEFT in which β decays are described by the Lagrangian [46,51,52] where N v = (p, n) T denotes the heavy-nucleon field doublet, v ρ is the nucleon velocity, and S ρ = (0, ⃗ σ/2) denotes the nucleon spin, with the Pauli matrices σ, while τ denotes Pauli matrices in the isospin space, satisfying [τ a , τ b ] = 2iε abc τ c , {τ a , τ b } = 2δ ab , and τ + = 1 2 τ 1 + iτ 2 . Higher-order terms in Eq. (1) include the contributions of weak magnetism, recoil corrections, and induced tensor coupling [44]. The couplings g V and g A themselves have an expansion in α, ϵ / π , and ϵ χ . At leading order, one has g V = 1. At O(α), g V does not receive any long-distance corrections from pion or photon loops and only picks up contributions from local O(e 2 p) operators in the HBChPT Lagrangian [44]: Here, C r β = 1 + O(α) is the Wilson coefficient of the Fermi operator in LEFT [see Eq. (9)], which captures electroweak corrections from energy scales above Λ χ . The LECs X 6 , g 9 , V 1,2,3,4 , and associated HBChPT operators will be defined below in Eqs. (36), (40), and (42). The couplings g V,A (µ χ ) depend on the renormalization scale of the / πEFT (in a way that cancels in the ratio λ = g A /g V ) and encode contributions from the weak scale all the way down to the pion mass scale.
In the following sections, we will detail the various steps needed to connect the low-energy coupling g V to the electroweak scale, following a top-down approach. Key new results of this work are as follows: (i) The expression for g V (µ χ ∼ Λ χ ) in terms of the Wilson coefficient C r β computed with anomalous dimensions of O(α, αα s , α 2 ) and a "subtracted" hadronic function, related to the traditional non-perturbative γW box contribution evaluated in the recent literature [1][2][3][4][5] (see Eq. (78) and discussion surrounding it); (ii) The use of two-loop anomalous dimensions in the RGE (88) needed to evolve the vector coupling down to g V (µ χ ∼ m e ), resumming large next-to-leading logarithms of order α 2 ln (m N /m e ). The resulting g V (µ χ ∼ m e ) is directly relevant to the calculation of neutron decay and can be used as input for the one-body contribution to nuclear decays.
In this work, we have focused on the application to neutron decay. With g V (µ χ ∼ m e ) at hand, we combined both virtual and real photon corrections to the decay rate [33,44,46] to obtain the effective phase-space correction ∆ f and the radiative correction ∆ R to the neutron lifetime, see Section 6, and the relation |V ud | 2 τ n 1 + 3λ 2 (1 + ∆ f ) (1 + ∆ R ) = 5283.321(5) s, with ∆ f and ∆ R given in Eqs. (110) and (111), respectively. Our definitions for ∆ f and ∆ R differ from the traditional approach both conceptually and numerically. Technically, the bulk of this difference is in shifting all short-distance contributions from ∆ f to ∆ R . ∆ f describes Coulomb-enhanced long-distance contributions and recoil corrections, while ∆ R includes all electroweak and HBChPT short-distance contributions along with the non-Coulomb radiative corrections in / πEFT, as specified in Eqs. (78), (89), and (114). Numerically, we find The uncertainty in ∆ f stems from an estimate of mixed recoil times Coulomb corrections. The dominant sources of uncertainty to ∆ R are given by the following: the non-perturbative hadronic contributions, associated with the "γW box" diagram in the standard approach [1][2][3][4][5][6]; contributions of O(αα 2 s ) not included in our renormalization group analysis in the LEFT; chiral corrections of αϵ 2 χ ; and residual dependence on the / πEFT renormalization scale, varied between m e / √ 2 and √ 2m e , which is an indicator of the O(α 2 ) corrections. A detailed discussion of uncertainties is presented in Sections 5.4 (for g V ) and 6.2 (for the remaining contributions to ∆ R ).
Our result for ∆ f in Eq. (5) differs from the one found in the literature ∆ f = 3.608 × 10 −2 [38] by −0.035%. This is because in the phase-space integration we use the nonrelativistic Fermi function, for the reasons discussed in Section 6.1, and neglect corrections induced by modeling the proton as a uniformly charged sphere of radius R p ≃ 1 fm [53] (this effect is at the level of 0.005%).
Our result for ∆ R in Eq. (6) exceeds the current value ∆ R = 3.983(27) × 10 −2 , compiled in Ref. [8] by combining the results of [1][2][3][4][5][6], by about twice the estimated uncertainties. The +0.061% shift in the central value is almost entirely due to the different treatment of the next-to-leading logarithmic terms at the hadronic level, i.e., the terms that scale as α 2 ln (m N /m e ). In both approaches, there is a contribution of this type coming from the cross term between the one-loop RGE correction to g V , scaling as α π ln (m N /m e ), and O (πα/β) terms in the Fermi function. In our approach, additional α 2 ln (m N /m e ) large logarithmic corrections arise entirely from the two-loop anomalous dimension contribution to the RGE (88) for the effective coupling g V (µ χ ) and produce a positive shift in ∆ R of 0.010%. In the EFT approach, there are no other sources of large logarithms of the ratio (m N /m e ) in the matrix element of the four-fermion operator (1) to O(α 2 ). In the literature, this class of effects is not associated with the running of g V , but arises through the negative correction α/(2π) × δ = −0.043%, introduced in Ref. [38] by adapting the results of Refs. [54,55]. 1 The mismatch of the two approaches produces a +0.053% shift in our results. The remaining difference is due to a combination of the following, individually smaller, effects: (i) we re-evaluate the "elastic" hadronic contribution, as discussed in Section 5.2, which leads to a −0.006% shift to ∆ R ; (ii) for the next-to-leading logarithmic corrections of O(α 2 ln(M W /m c )), our result differs from the one in Ref. [38], producing a negative shift of approximately −0.011%; (iii) we do not include O(αα 2 s ) terms in the running of our Wilson coefficient (corresponding to the "deep inelastic scattering" region of the γW box in the literature) that amounts to a net +0.007% in ∆ R ; and (iv) finally, different choices in the factorization between electroweak and m N /m e logarithms compared to Refs. [8,38] account for the remaining mismatch. Using ∆ f,R from Eqs. (5) and (6), respectively, in the master formula (4), we can extract V ud . This requires experimental input for the neutron lifetime τ n and the ratio λ of axial-vector to vector couplings. Using the PDG [56,57] averages for the experimental input, we obtain Both τ n and λ carry an inflated error due to scale factors. Following Ref. [8], if we instead use the most precise neutron lifetime measurement τ n = 877.75 (36) s from UCNτ @LANL [58] and the determination of λ from the most precise measurement of the beta asymmetry in polarized neutron decay by PERKEO-III [59,60], we obtain a very competitive extraction of V ud from neutron decay: with an uncertainty approaching the currently quoted error δV ud = 31 × 10 −5 from 0 + → 0 + nuclear beta decays [7]. Compared to the baseline correction of Refs. [1][2][3][4][5][6]8], the positive shift of +0.061% in ∆ R and the negative shift of −0.035% in ∆ f partially compensate, producing a smaller positive shift of +0.026% in the correction to the rate. This one, in turn, provides a negative shift in V ud , δV ud ≃ −13 × 10 −5 , compared to the results quoted in Ref. [8].
In the remainder of this paper, we provide details on the derivation of the results presented above. 3 Step I: matching the Standard Model to LEFT In this Section, we perform the matching of the Standard Model to the LEFT and present the RGE that control the effective couplings in the LEFT between the electroweak and QCD scales. We then introduce spurions and external sources in the LEFT to describe the electromagnetic and weak interactions of quarks [47,61], which is particularly useful in the matching of LEFT to chiral perturbation theory, to be described in subsequent sections. Throughout, we regulate the UV divergences in dimensional regularization, working in d = 4 − 2ϵ spacetime dimensions.

Wilson coefficient and RGE
The part of the LEFT Lagrangian relevant for muon and β decays just below the weak scale reads Here µ denotes the MS renormalization scale and is the scale-independent Fermi constant that is extracted from precise measurements of the muon lifetime [62][63][64][65], expressed in terms of MS Standard Model parameters (with s 2 . The function g (µ) can be found in Ref. [66] and reduces to g (µ) = 1 at tree level. The effective coupling multiplying the semileptonic operator that mediates β decays involves the same G F as the pure-leptonic term in Eq. (9), the CKM matrix element V ud , and the MS-subtracted Wilson coefficient C r β (a, µ), which reads [36,66,67] C r β (a, µ) = 1 + The finite O(α) matching coefficient depends on the scheme through B(a). We have used the Naive Dimensional Regularization (NDR) scheme for γ 5 and kept track of the additional evanescent operator scheme dependence via the parameter a, defined by [68][69][70] with an evanescent operator E(a) that has a vanishing matrix element in d = 4. Current conservation protects C β from O(α s ) corrections. Concerning the terms of O(αα s ), we only keep logarithmic contributions, as the finite matching coefficients and the corresponding three-loop anomalous dimensions are not known. The renormalized Wilson coefficient C r β (a, µ) obeys the following RGE: γ se = +1 [36,66,71], whereñ is the scale-dependent effective number of fermions, α (µ) and α s (µ) are the electromagnetic and strong running coupling constants. We have obtained γ N DR 1 (a) by adapting the QCD calculation in [68]. As far as we know, this is the first time the full two-loop anomalous dimension is worked out. 2 With appropriate rescalings of the QCD diagrams of Ref. [68], we also reproduce γ se = 1. γ 0 and γ se are scheme-independent. The scheme independence of γ se follows from the general argument given in Ref. [72], combined with the fact that there is no finite matching term nor anomalous dimension to O(α s ) for the operator under study here. On the other hand, γ 1 depends on both the treatment of γ 5 in d spacetime dimensions and on the scheme used for evanescent operators.
In our final result, we will use the numerical solution for C r β (a, µ). However, it is quite instructive to provide an approximate analytic solution, based on the perturbative treatment of the next-to-leading logarithm (NLL) terms associated with the scheme-dependent two-loop anomalous dimension γ 1 (a) = O(α 2 ) and the finite one-loop matching condition B(a). First, setting γ 1 (a) → 0 and consistently taking as an initial condition C r β (a, µ SM ) = 1, generalizing the result of Ref. [71] given at the τ -mass scale, we obtain the solutionC where we have subsequently integrated out the b quark, τ lepton, and c quark, and the strong and electromagnetic running couplings are obtained by solving the one-loop RGEs. This solution resums all the terms of O(α n ln n (µ SM /µ)) and O(αα n s ln n (µ SM /µ)). We can then perturbatively include the effects of O(α 2 ln(µ SM /µ)) due to γ 1 (a) and B(a), arriving at where and the scheme-independent combination κ is given by 3 In the equation above, β 0 = −(4/3)ñ controls the one-loop β function for α via µdα/dµ = −(β 0 /(2π))α 2 .
The scale-dependent effective number of fermions takes the valuesñ(µ < m c ) = 4,ñ(m c ) = 16/3,ñ(m τ ) = 19/3, andñ(m b ) = 20/3. Note that the scheme dependence of C r β (a, µ) in the solution (16) appears only through the initial factor involving B(a). As we will show below, this term explicitly cancels when one includes the O(α) corrections to the matrix element of the semileptonic operatorū L γ α d LēL γ α ν eL .

External sources and spurions
The matching of LEFT to HBChPT is conveniently performed by introducing classical source fieldsl µ (x) andr µ (x) for the left-and right-handed chiral currents of quarks as well as electromagnetic left q L and right q R spurions, and the weak spurion q W [47,61,76,77]. These allow one to handle the explicit chiral symmetry breaking introduced by the electromagnetic and weak interactions at the quark level in a compact way. With this motivation in mind, we write the source term for currents plus the QED and weak CC interactions of the light quarks q T = (u, d) as where A µ denotes the photon field. The Lagrangian in (20) is with L, R ∈ SU (2) L,R , provided q L,R and q W transform as "spurions" under the chiral group, namely q L,W → Lq L,W L † and q R → Rq R R † , and thatl µ andr µ transform as gauge fields under G. At the physical point, with τ + = (τ 1 +iτ 2 )/2 in terms of the Pauli matrices τ a . Note that we include the LEFT Wilson coefficient C r β in the definition of the spurion q W . With this identification, Eq. (20) reproduces the semileptonic piece of Eq. (9).
The O(e 2 ) counterterms in the LEFT Lagrangian can be written in terms of spurions as [47] L CT where g 00 is the counterterm related to the electron wavefunction renormalization, g 02 and g 03 come from the counterterm of C β , while g 23 includes contributions from both the counterterm of C β as well as divergences related to the quark wavefunction renormalization. Furthermore, are chiral covariant derivatives, expressed in terms of the fields l µ (x) and r µ (x) that combine the classical sources, the photon, the leptons, and the spurions: In the MS scheme, the g ij couplings appearing in Eq. (23) are determined by the 1/ε divergences and can be written as with h 00 = 1/2, h 23 = (1/2)(1 − α s /π), h 02 = −1 − α s /π, and h 03 = 4 − 2α s /π. 4 Step II: matching LEFT to HBChPT The goal of this section is to find a representation for the LECs appearing inĈ V , see (3), in terms of the LEFT counterterms g ij and quark correlation functions, which can then be modeled or computed via non-perturbative techniques such as lattice QCD.

The Chiral Lagrangian
The chiral representation of Eq. (23) can be built using standard spurion techniques. As in Eq. (23), we will need purely leptonic operators, purely electromagnetic operators, and operators with charged leptons and neutrinos. The corresponding chiral Lagrangians were built in Refs. [44,61,[78][79][80]. Here we extend the bases of [44,80] in order to avoid assumptions regarding q L and q R , allowing us to keep the spurions q L,R completely general. Moreover, we do not use the equations of motions to reduce the operator set in order to avoid hadronic contributions to purely leptonic LECs [47]. As we will see below, to perform the matching between LEFT and HBChPT it is convenient to introduce vector and axial-vector charge spurions and sources, which we define as It is also convenient to decompose the electromagnetic charge spurions in isovector and isoscalar components with J ∈ {L, R, V, A}. The physical values are q 0 L,R = q 3 L,R = 1 2 for the left and right spurions, q 0 V = q 3 V = 1 for the vector spurion, and q 0 A = q 3 A = 0 for the axial-vector case. 4 The chiral Lagrangians are built using the chiral covariant functions of the charges and of the corresponding covariant derivatives in Eqs. (24) and (25): At lowest order, the HBChPT Lagrangian is given by where F π and g (0) A denote the pion decay constant and the nucleon axial-vector coupling in the chiral limit. u µ and χ + are given by with the light quark masses m q and a LEC of dimension of mass B 0 , which is related to the quark condensate. We further introduced the nucleon chiral covariant derivative where by the superscript baryon we indicate that the photon couples to the nucleon via the charge q baryon V in Eq. (30). In addition to the weak and electromagnetic interactions arising from chiral covariant derivatives, Eq. (32) contains electromagnetic effects mediated by high-momentum photons via the coupling Z π , which is related to the pion-mass splitting.
The chiral Lagrangian needed at O(G F α) is given by L CT lept is a purely leptonic counterterm Lagrangian The coupling X 6 is determined by computing the electron propagator in LEFT and chiral perturbation theory, obtaining in arbitrary R ξ gauge, where µ and µ χ are the LEFT and HBChPT renormalization scales, respectively. X r 6 (µ χ , µ) denotes the renormalized coupling, after subtraction of the 1/ε pole in the MS χ scheme. Note that, following standard practice [49], in the MS χ scheme, we subtract instead of the conventional MS subtraction used in LEFT: For the electromagnetic Lagrangian L e 2 p πN , we use the construction of Ref. [80]. Only one operator is required to describe Fermi transitions, For the electroweak sector with charged leptons and neutrinos, we provide the most general weakinteraction Lagrangian in the heavy-baryon sector with one charge and one weak spurion, where we only assumed the constraint ⟨q W ⟩ = 0 [81]: where g (0) A denotes the nucleon axial-vector coupling in the chiral limit and The dimensionless low-energy coupling constants in Ref. [44] are related to the couplings in Eq. (41) by A when the spurions take physical values. In Ref. [44], the authors have used the equations of motion to eliminate some S ρ -dependent operators. In addition, they reduced operators that are bi-linear in spurions to linear expressions by exploiting the relations q L q W = (2/3)q W and q W q L = −(1/3)q W , valid for physical values of the spurions.
As realized in the mesonic sector in Refs. [47,77], we can interpret amplitudes in LEFT and in HBChPT as functionals of the same charges q 0,a (x), promoted to be spacetime-dependent external fields. The matching between the LEFT and HBChPT can then be obtained by equating functional derivatives of the effective action with respect to q 0,a (x) in both theories. As we will see, this allows us to derive an explicit representation for the LECs and to keep track of unphysical scale and scheme dependence appearing at intermediate steps of the calculation.

Electromagnetic coupling constant
We start from the electromagnetic coupling g 9 . Expanding the charge covariant derivative in Eq. (40), we obtain g 9 can then be evaluated by taking functional derivatives with respect to two isovector charges, provided that the charge carries non-zero momentum, or by taking three derivatives, two with respect to the charges and one with respect to a vector or axial-vector source. The first representation is more convenient, since, as we will see, it allows one to automatically obtain cancellations between electromagnetic and weak couplings. More precisely, we will consider the following object: Figure 1: Diagrams that contribute to Γ V V in HBChPT are shown. Double, wiggly, and dashed lines denote nucleons, photons, and pions, respectively. Dashed circles denote insertions of the sources q a,b V . The arrows denote the flow of the momentum r inserted by the sources. The first three diagrams originate from the leading-order π and πN Lagrangians, L p 2 π , L e 2 π , and L p πN [44,80,82], which are presented in Eq. (32). The last diagram denotes contributions from L e 2 p πN and is proportional to g 9 .
where k and k ′ are the nucleon momenta, σ and σ ′ denote the nucleon spins, and i, j are the nucleon isospins. We take the nucleon to be at rest, k = k ′ = m N v and use the nonrelativistic normalization for heavy-particle states πN Lagrangian. g 9 provides the only contribution to Γ V V . The loops are determined by couplings in the leading-order (LO) pion and pion-nucleon Lagrangians. In particular, the diagram with pion-mass splitting Z π is symmetric in isospin, and vanishes once contracted with the Levi-Civita tensor, so that the loop corrections are purely determined by the minimal coupling of the photon to the nucleon. In arbitrary R ξ gauge, we introduce the photon mass λ γ as an infrared regulator and obtain g r 9 (µ χ , µ) in the second line denotes the renormalized coupling, after subtraction of the 1/ε pole in the MS χ scheme. For ξ = 1, the anomalous dimension of g r 9 (µ χ , µ) agrees with the result of Ref. [80], so that Eq. (45) is independent of the scale µ χ .
In the LEFT, the same matrix element is given by Eq. (46) contains a tree-level term, proportional to the counterterm g 23 that cancels the divergences generated by loop diagrams. The loop contribution contains the hadronic tensor T µν V V (q, v), which can be expressed in terms of the two-point correlation function of quark currents. Here, we use the definition [34] The gauge-dependent term in Eq. (46) is obtained using which follows from the conservation of the vector current.
To highlight the UV structure of Eq. (47), we add and subtract the high-energy limit of the hadronic tensor provided by the operator product expansion (OPE) where for the OPE of the relevant currents we use results from Refs. [83,84], adapted to include the appropriate color factors [35]. Since our calculation is only accurate at the leading logarithm in O(αα s ), the O(α s ) correction to the OPE is computed in d = 4. Note that in Eq. (49) we have introduced an arbitrary scale µ 0 to regulate infrared divergences that appear when evaluating the convolution integrals with T OPE . Performing the relevant integrations, we obtain where T denotes the subtracted hadronic tensor, T = T − T OPE . T depends on µ 0 in such a way that the final results are µ 0 -independent. Finally, note that we are dropping terms of O(αα s ) that appear without logarithmic enhancements, because they are beyond the accuracy of our calculation. Equating Eqs. (45) and (46), we obtain a representation for g 9 : Alternatively, to control the infrared region and see a cancellation of the infrared divergences, we can introduce the combinationT = T − T IR , where T IR is the leading infrared contribution g µν T µν IR = i/ (v · q), and obtain

Electroweak coupling constants
We follow the same strategy for the determination of the electroweak coupling constants. In this case, the operators V 1 and V 2 receive contributions from the isovector component of the electromagnetic charges, while V 3 and V 4 receive from the isoscalar component. We thus define two matrix elements At the order we are working, the electron and neutrinos can be taken to be massless and to carry zero momentum.  Fig. 1. In this case, the sources inject zero momentum. The first two diagrams originate from the LO πN Lagrangian L p πN , and the last diagram denotes contributions from L e 2 p πN ℓ . Diagrams with the sources coupling to pions do not contribute at this order.
The HBChPT diagrams contributing to Γ (0,1) V W are shown in Fig. 2. The loop diagrams cancel for isoscalar electromagnetic couplings, so that we obtain In the LEFT, the isovector and isoscalar components are given bȳ The hadronic tensors with two isovector currents are defined in Eq. (47), while we define the hadronic tensor with one isoscalar vector current as As in Sec. 4.2, the UV divergences in the LEFT are determined by the operator product expansion. In NDR, the leading-order OPEs of T µν V V − T µν V A and T µν V V,0 − T µν V A,0 are proportional to the symmetric and antisymmetric combinations of Dirac matrices (γ µ / qγ ν ± γ ν / qγ µ )P L , respectively. The symmetric combination does not depend on the scheme, while the antisymmetric piece depends on the definition of the evanescent operators, in such a way as to compensate the dependence of the couplings in the LEFT. Using the OPE, we obtain The integrals of the subtracted hadronic tensorsT are convergent, so that we can perform the Dirac algebra on the leptonic leg in the d = 4 dimension. Putting everything together, we arrive at the matching equations To obtain the second line of Eqs. (62) and (63), we used the Ward identities on the subtracted tensors, the symmetry (antisymmetry) of unpolarized hadronic tensors T µν V V (0) (T µν V A(0) ) under µ ↔ ν, and, in the contractions with the Levi-Civita tensor, we replaced The non-perturbative QCD input in the LECs is encoded in the subtracted hadronic tensors T V V , T V A , T V V, 0 , and T V A, 0 . Using time reversal and crossing symmetry [2,33], we can show that the scalar functions in the matching equations (62) and (63) are odd or even under q → −q. Explicitly we have where we indicated that the functions depend only on the invariants As a consequence of Eqs. (66) and (67), T V A and T V V, 0 do not contribute to the matching, and the final expressions for the combinations of LECs V 1 + V 2 and V 3 + V 4 are Note that in this framework, the LECs depend not only on the chiral renormalization scale (µ χ ) but also on the LEFT renormalization scale (µ) and the schemes adopted for γ 5 and the evanescent operators.

Corrections to g V
In this Section, we combine the coupling constants of the heavy-baryon chiral perturbation theory into the counterterm of g V in / πEFT. We subsequently evaluate the non-perturbative inputs to the vector coupling constant, resum logarithms between the chiral and electron-mass scales, and provide numerical results for g V .

Matching at the baryon-mass scale
Having determined the electroweak coupling constants V 1 -V 4 and the electromagnetic coupling constant g 9 , we can evaluate the O(α) contribution to g V in the low-energy effective theory, cf. Eqs. (2) and (3). These corrections are known in the literature as "inner" radiative corrections.
Before getting to the final result, we can combine the LECs that depend on the V V hadronic tensor, g 9 and V 1 + V 2 , and the lepton wavefunction renormalization X 6 , obtaining which is independent of the gauge parameter ξ. T V V enters this combination of LECs multiplied by the IR regulator λ 2 γ . The only contribution to the integral can thus come from the infrared limit of T V V , where the hadronic tensor is well approximated by the elastic piece. The integral over the hadronic tensor then only leaves behind a finite piece, yielding Thus, the only contributions to − X 6 2 + 2 (V 1 + V 2 ) − g 9 are due to the different renormalization scales, µ vs µ χ , and the different subtraction schemes commonly used in HBChPT, MS χ vs MS.
The other combination of LECs V 3 + V 4 is conveniently expressed in terms of the scalar amplitude T 3 ν, Q 2 as where we defined the amplitude T 3 from the tensor decomposition of the hadronic tensor as [85][86][87][88][89][90] with the OPE-subtracted expression In the OPE, we have retained the O(α s ) correction, which is needed to cancel the µ-dependent term proportional to αα s ln(M W /µ) in C r β . To the order we are working, we can use α s (µ) at any µ where QCD is perturbative. We will use α s (µ 0 ) in what follows.
Combining the HBChPT coupling constants into the / πEFT countertermĈ V according to Eqs. (2), (3), (72), and (73), we achieve the matching condition where we resummed logarithms in the Wilson coefficient C r β (a, µ), as it is described in Section 3.1. This expression does not contain electroweak-scale parameters or artificial hadronic scales, besides the dependence contained in the coupling constant C r β (a, µ). The vector coupling g V (µ χ ) does not depend on the scale and scheme used in the LEFT at the one-loop level.
We can further simplify the expression for g V (µ χ ) and connect it to the previous literature. First, we eliminate the evanescent scheme dependence by defining the scheme-independent NLO Wilson coeffi- which can be immediately read off from Eq. (16). We then have where the non-perturbative input is in the "subtracted" hadronic contribution □ We will evaluate the non-perturbative input in Eq. (79) in Sec. 5.2. Eq. (78) encodes the so-called "inner" radiative corrections to the Fermi transitions in the EFT language in the form of a µ χ -dependent coupling g V (µ χ ), which appears in the effective Lagrangian of Eq. (1). Once all large electroweak logarithms are resummed via the RGE in C β (µ), Eq. (78) does not contain additional large logarithms when the scales µ χ , µ, and µ 0 are similar and of order Λ χ ∼ 1 GeV. As shown below, the µ χ -scale dependence in g V (µ χ ) is canceled in physical amplitudes by the µ χ dependence of the virtual photon corrections computed in the pionless theory. Since the only scale of these loops is O(m e ), we will evolve g V (µ χ ) down to the scale µ χ ∼ m e in order to avoid large logarithms, see Sec. 5.3.

Evaluation of the non-perturbative input
As shown in Refs. [1,2], the box function can be represented as a one-dimensional integral over the Q 2 > 0 variable where F (Q 2 ) = (12/Q 2 )M 3 (1, Q 2 ) and M 3 (1, Q 2 ) is the first Nachtmann moment of the structure function defined in terms of the imaginary part of T 3 (ν, Q 2 ). Following Refs. [1,2], it is useful to isolate the well-defined elastic contribution to F (Q 2 ), which we denote by F el (Q 2 ), known in terms of the nucleon isoscalar magnetic vector and axial-vector form factors, and define where F (Q 2 ) includes inelastic contributions. For Q 2 ≤ Q 2 0 = 2 GeV 2 , 6 F (Q 2 ) contains contributions from the resonance region and the so-called Regge region. Current knowledge is based on modeling [1][2][3][4][5] and lattice QCD input [6]. For Q 2 ≥ Q 2 0 = 2 GeV 2 , one enters the deep inelastic scattering region (DIS), controlled by the OPE with Wilson coefficients computed in perturbative QCD (pQCD). The OPE representation of F (Q 2 ) is known to leading order in 1/Q 2 , with coefficients known to O(α 4 s ) [3,[91][92][93]: In practice, we will use only the n = 1 term (with coefficientc 1 = 1) in ∆(Q 2 ), as higher-order terms are beyond the accuracy of our NLL LEFT analysis. Moreover, for consistency with the OPE terms that we subtract in the matching procedure, we will use ∆(Q 2 ) → ∆(µ 2 0 ) in F DIS (Q 2 ). In terms of the quantities defined above, the subtracted hadronic contribution reads Isolating the elastic contribution and separating the integration in the regions below and above Q 2 0 = ( √ 2 GeV) 2 , we find Numerically, for the non-perturbative contributions we find = (0.49(11) + 0.04(1)) × 10 −3 .
Up to negligible contributions of O(Q 2 0 /M 2 W ), the integral of F (Q 2 ) between 0 and Q 2 0 in Eq. (87b) coincides with the Q 2 ≤ Q 2 0 inelastic piece of the "box diagram", recently considered in the literature [1][2][3][4][5][6]. The result is usually written as the sum of the "Regge" plus "Resonance" contributions. The various evaluations in the literature have recently been combined by Ref. [8], leading to the numbers used in Eq. (87b). This part of our result is fully correlated with previous work and carries the dominant contribution to the error budget for the radiative corrections.

RG evolution of g V below the baryon scale
To account for higher-order perturbative logarithms, which are needed for precise predictions of β-decay rates and (anti)neutrino-nucleon scattering, we evolve the low-energy coupling constant g V (µ χ ) from the matching scale µ χ ∼ Λ χ to the physical scale µ χ ∼ m e using the one-and two-loop anomalous dimensions. The vector coupling constant evolves according to with the effective number of particlesñ, as it is described in Appendix A.2. The appropriate one-loop anomalous dimensionγ 0 has been identified in Refs. [38,44,46,51,96,97]. 7 It can also be extracted from calculations of the "heavy-light" current QCD anomalous dimension in the context of heavy quark physics, as for example in Refs. [98,99]. As discussed in Appendix C, we can exploit this analogy to extract the QED two-loop anomalous dimensionγ 1 , by adapting the results from Ref. [100] (see also Refs. [101][102][103]). The above expression forγ 1 only includes two-loop diagrams involving two virtual photons in the pionless theory. Possible contributions arising from diagrams involving pions and photons are not included. Note that the term inγ 1 proportional to π 2 can lead to contributions to the decay rate that scale as α 2 ln (m N /m e ), larger than a typical two-loop contribution. Using the expression for the evolution operator in Eq. (19), we solve the RGE in (88) and resum the leading and subleading logarithms between particle thresholds according to Below the baryon scale, we determine α from its value in the Thomson limit by evolving it up in scale with the electron, muon, and charged pion as active degrees of freedom, which leads tõ See Appendix A for details on the definition of the fine-structure constant in both LEFT and χPT. In Eq. (89), g V (Λ χ ) is obtained by evaluating Eq. (78) at µ χ = Λ χ ∼ m N . Note that bothγ 0 andγ 1 are negative, implying g V (m e )/g V (Λ χ ) > 1.

Numerical results and uncertainty estimates
We next present numerical results for the vector coupling g V (m e ) and discuss the various sources of uncertainty. We start by providing some intermediate results that illustrate the impact of corrections at various orders in our RGE analysis.
For the semileptonic Wilson coefficient C r β , we include α, αα s , and α 2 contributions to the running, as described in Section 3.1. 8 To illustrate the effect of running from the electroweak to GeV scales, we provide results for the fixed order (LO) C β (m c ) = 1 + (α(m c )/π) ln(M Z /m c ), leading logarithms (LL), next-to-leading logarithms NLL1, which includes the anomalous dimensions up to order αα s , and nextto-leading logarithms NLL2, including the anomalous dimensions up to orders αα s and α 2 . For the initial conditions, we specify After numerically solving the RGEs, we obtain the following values for the effective couplings at µ = m c : The effects of NLL1 and NLL2 resummations combine to essentially "undo" the effect of LL resummation. The final result is very close to the perturbative one. The numerical solution of the RGEs agrees with the analytic solutions provided in Section 3.1. Our result for the NLL1 correction is consistent with the finding of Ref. [71]. The impact of NLL2 corrections in our result is more than a factor of 2 larger than in Ref. [38], reflecting the difference discussed in Section 3.1.
For the running of the vector coupling constant g V , we include the O(α) and O(α 2 ) anomalous dimensions, as described in Section 5.3. 9 We provide the relative running contributions for the one-loop logarithm (LO), namely g V (m e )/g V (m p )| LO = 1 + (3/4)(α/π) ln(m p /m e ), the LL resummation, where we include onlyγ 0 in the RGE, and NLL resummation, where we also includeγ 1 in the RGE: At the level of decay rate, our NLL correction implies an increase of 1.0 × 10 −4 (roughly half of the final uncertainty of the radiative corrections). Putting together all the results obtained so far, we evaluate the vector coupling constant g V (µ χ ) in the MS renormalization scheme of χPT at the scale µ χ = m e , where O(α n ) loop corrections to the matrix elements of the Lagrangian (1) do not contain large logarithms: In contrast, the vector coupling at fixed order g 1−loop V (i.e., without resummation, without αα s corrections, and with taking the value for the electromagnetic coupling constant in the Thomson limit) takes the value In the RGE evolution, the electromagnetic and αα s effects contribute with opposite signs, resulting in a net increase of g V at the level of 0.07%. For the uncertainty estimate, we add the following dominant sources in quadrature: • 0.012%: the hadronic error for Regge, resonance, and πN contributions from Ref. [45] is added in quadrature to the uncertainty propagated from the lepton-nucleon experimental data for the elastic contribution.
• 0.004%: the higher-order αα 2 s uncertainty is estimated by including the known terms of O(α 2 s ) [1][2][3]38] in the pQCD correction ∆(µ 2 0 ) that controls the DIS region of □ V γW in Eq. (85). In our approach, this DIS contribution maps onto the αα 2 s anomalous dimension for the Wilson coefficient C β (µ) in LEFT.
All other perturbative and parametric sources of uncertainties are at the level 0.001% or even below. We conclude this section by noting that the effective coupling g V (µ χ ≈ m e ) captures the "inner" corrections to one-body weak transitions through NLL, i.e., up to and including terms of order α 2 L 2 and α 2 L (where L indicates large logarithms of M Z /m N and m N /m e ), with residual uncertainty at O(α 2 ) due to finite terms in two-loop diagrams. Importantly, g V controls both neutron decay and the one-body contribution to nuclear β decays, in combination with appropriate n → peν e and (N, Z) → (N −1, Z+1)eν e matrix elements computed to the same accuracy. For applications in neutrino and nuclear physics, in Table 1 we provide the coupling constant g V for a few values of the renormalization scale up to 50 MeV.

Corrections to neutron decay and impact on V ud
We can now use the / πEFT Lagrangian in Eq. (1) with g V (µ χ = m e ) from Eq. (89) to compute the neutron decay rate including radiative corrections. The final ingredient is the square modulus of the n → peν e and n → peν e γ matrix elements in HBChPT, evaluated at µ χ ∼ m e . To match the accuracy achieved in g V (m e ), since ln(µ χ /m e ) ∼ O(1), we will need the matrix elements to O(α) and will ignore terms of O(α 2 ) and higher. The only exceptions are "Coulomb"-enhanced terms scaling as (πα/β) n and α/π(πα/β) n , where β ≡ p e /E e , which are parametrically large, diverge for β → 0, and can be resummed in the nonrelativistic Fermi function.

"Long distance" electromagnetic corrections and differential decay rate
After including the contributions from both virtual and real photons [44,46] as well as recoil corrections [46,53], the differential decay rate dΓ n for unpolarized neutrons takes the form [33,53] where E 0 = (m 2 n − m 2 p + m 2 e )/(2m n ) is the electron endpoint energy and λ ≡ g A /g V is the ratio of effective axial-vector and vector couplings in the low-energy Lagrangian (1). The ratio λ = λ QCD (1 + δ (λ) RC ) is affected by a µ χ -independent electromagnetic correction δ (λ) RC parameterized in terms of calculable pion loops and certain chiral LECs (see Ref. [44]). λ itself can be extracted from beta decay correlation experiments, so that we do not need to know δ (λ) RC for the purpose of studying total decay rates and the extraction of V ud . δ recoil (E e ) collects recoil corrections that can be found in Ref. [46]. They are usually factorized since the impact of the product of radiative times recoil corrections is estimated to be well below 10 −4 . Finally,δ RC (E e ) represents the electromagnetic corrections arising from the matrix element squared. To O(α), one finds whereĝ(E e , E 0 ) is a "subtracted" Sirlin function defined in terms of the Sirlin function g (E e , E 0 ) of Ref. [33].ĝ (E e , E 0 ) arises naturally in the EFT calculation and does not contain any large logarithm of m N /m e . The corrections proportional to πα/β in Eq. (98) are enhanced by a factor of π 2 compared to the naive scaling of loop corrections, and are numerically dominant even for β ∼ O(1). The leading terms in the series in πα/β arise from the momentum regions of loop integrals in which the photon momentum has potential scaling, k 0 ∼ m e β 2 ≪ | ⃗ k| ∼ m e β, and they can be identified with nonrelativistic EFT methods [104][105][106][107]. Their resummation leads to the nonrelativistic Fermi function F N R (β) [108][109][110][111][112][113][114][115][116][117][118], which we include in the matrix element squared as where As we discuss in Appendix B, the factorization ansatz in Eq. (101) captures all numerically-enhanced leading and subleading terms in 1/β, and reproduces similar results for the production of two heavy quarks at threshold, derived with nonrelativistic QCD and potential nonrelativistic QCD [104][105][106][107][119][120][121][122]. At O(α 2 ), Eq. (101) gives Indeed, the first cross term −(11/4)α 2 /β corresponds to the matching coefficient of heavylight to heavy-heavy current [123] in the MS χ renormalization scheme. The second cross term (α 2 /β) (E 0 − m e ) 2 / 12m 2 e comes from the product of the Fermi function with real radiation. These Only the first two diagrams give rise to terms in theγ 1 enhanced by π 2 [100]. These diagrams also give rise to the leading π 2 α 2 /β 2 behavior captured by the nonrelativistic Fermi function.
terms are beyond the accuracy of our calculation and can be booked as O(α 2 β 3 ) in the nonrelativistic limit. In the case of neutron decay, this term provides a negligible shift of 1.6 × 10 −5 to the decay rate.
We thus arrive to our final form for the differential decay rate: Compared to state-of-the-art analyses of neutron decay in the literature (see e.g. Ref. [38]), our result (104) amounts to replacing the relativistic Fermi function [53,[109][110][111][124][125][126][127] with the nonrelativistic one, F 0 → F N R . While we arrived at this result by constructing the relevant terms of the amplitude in the EFT framework, one could also argue for this replacement along the following lines. First, recall that the leading corrections to the phase space coming from the distortion of the electron wavefunction in the Coulomb field of the proton is usually captured by the function [53] F 0 (β) = 2 1 + γ F (β) = 4(2E e βR) 2(γ−1) e πy |Γ(γ + iy)| 2 (Γ(1 + 2γ)) 2 , y = This form is obtained by solving the Dirac equation for an electron moving in the charge distribution of a uniformly charged sphere of radius R [53], but corresponds to a rescaling of the solution of the Dirac equation for a point-like proton, F (β), evaluated not at the origin, where the wavefunction diverges logarithmically, but at the "nucleon radius" R. R corresponds to a mass scale much larger than m e and effectively acts as a UV regulator. So we see that while F 0 (β) coincides with F N R (β) at the one-loop level, F 0 includes a dependence on the UV regulator via the logarithms of R that first appear at O(α 2 ). Expanding F 0 in series of α, one obtains The dependence on the UV regulator R ∼ 1/µ χ does not match the µ χ -dependence of g V (µ χ ) in the MS χ scheme presented so far. In dimensional regularization, indeed, the ln R term in Eq. (106) corresponds to a UV singularity that appears in the first two diagrams in Fig. 3, when we consider only the contribution arising from picking the two nucleon poles. This is only one piece of the full anomalous dimensionγ 1 . In order not to double-count large logarithms, one should set the logarithmic term in F 0 to zero when using the RGEs to evaluate the large logarithms as we do here. The remaining O(α 2 ) terms in Eq. (106) are incomplete and beyond the accuracy of our calculation, which allows us to drop them and replace the relativistic Fermi function F 0 by its nonrelativistic counterpart F N R .

Total decay rate and extraction of V ud
Upon performing the integration over E e in Eq. (104), the decay rate can be written as where the phase-space integral is given by with x 0 = E 0 /m e and E 0 = 1.292581 MeV, and takes the value f 0 (x 0 ) = 1.62989. Following standard practice [38,53], in Eq. (107) we have lumped the Coulomb (F N R ) and recoil terms into an effective phase-space correction ∆ f , separating the remaining radiative corrections into ∆ R . In this factorization scheme, the various corrections to the decay rate are defined by where β(x) = 1 − 1/x 2 . A few remarks are in order: • The decay rate in Eq. (107) corresponds to the usual definition adopted in the literature [38], upon identifying f ≡ f 0 (1 + ∆ f ). Therefore, the total shift in the decay rate which impacts the extraction of V ud , requires specifying both ∆ f and ∆ R . The expressions and numerical values of ∆ f and ∆ R in our EFT approach differ from the results found in the literature (see Ref. [38] and most recent calculations of ∆ R [1][2][3][4][5][6]8]). In what follows, when necessary we will discuss the origin of the differences.
• For ∆ f , which encodes Coulomb and recoil corrections, we find where we estimated the uncertainty to be of the size of Coulomb corrections times the recoil cross term. The difference from the standard result ∆ f = 3.608 × 10 −2 [38] is mainly due to the fact that we use the nonrelativistic Fermi function, for the reasons discussed above, while Ref. [38] uses the relativistic Fermi function. We also do not include the corrections induced by modeling the proton as a uniformly charged sphere of radius R p ≃ 1 fm [53]: this is a small effect shifting ∆ f by 0.005%.
• Up to the accuracy of our calculation, the remaining radiative correction ∆ R in our framework is given by where µ χ ∼ m e andĝ (E 0 ) = −9.58766 is obtained by averaging the subtracted Sirlin function g(E e , E 0 ) over the phase space, according to Eq. (111). At leading order in α, the µ χ -scale dependence in Eq (114) cancels between the coupling constant g V (µ χ ) and virtual one-loop contributions, while higher-order perturbative logarithms from virtual diagrams at scales µ χ ∼ m e are small.
• To separate hadronic and electroweak contributions to g V (µ χ ), and to make contact with some of the previous literature, we provide the fixed-order result In the above relations, the explicit dependence on µ 0 is canceled by the implicit dependence in □ V Had (µ 0 ). The hadronic physics is included in □ V Had , while the two logarithms in Eq. (115), which are proportional to the anomalous dimensions, correspond to the ratios between electroweak vs hadronic and hadronic vs beta-decay scales.
which, apart from the uncertainty coming from g V discussed in Sect. 5.4, includes a perturbative uncertainty of 0.005% obtained by varying the scale of the calculation µ χ in the range m 2 e /2 ≤ µ 2 χ ≤ 2m 2 e . Our result for ∆ R is 0.061% above the most recent evaluation [8] based on Refs. [1][2][3][4][5][6]. The sources of this difference are discussed in Section 2. Combining ∆ f and ∆ R in the factorization scheme of Eq. (107) we obtain ∆ TOT = 7.761(27)%.
• As a consistency check on the accuracy of the calculation and the size of cross terms (such as recoil × electromagnetic corrections), we have performed the phase-space integration in a different scheme that does not assume factorization of F N R and δ recoil , defined by with For the numerical values in this scheme, we find ∆ g V = 5.060(27)%, ∆ C = 3.375%, ∆ RC = −0.969%, and ∆ recoil = 0.173%, leading to ∆ TOT = 7.770%. The latter differs from the factorized result by 0.009%, consistent with its expected size of O(α 2 ) and the uncertainties quoted above.
Finally, we extract the CKM matrix element V ud from precise measurements of the neutron lifetime with our updated calculation of radiative corrections and present the results in Section 2.

Comments on radiative corrections to nuclear decays
Finally, we comment on the connection to the standard framework for the analysis of superallowed 0 + → 0 + transitions, described for example in Ref. [7]. The corrections to nuclear beta decays are combined into the quantity Ft, related to the experimental f t values as where K is a constant and δ C is the isospin-symmetry breaking contribution. The correction δ N S corresponds to the transition-dependent nuclear structure correction. ∆ V R is the transition-independent part of the radiative correction, which is related to the correction to neutron decay via δ ′ R contains the so-called "outer corrections", which depend on the transition but not on the nuclear structure, and corresponds to soft photon emissions from point-like nuclei. δ ′ R reduces to the Sirlin function at O(α) and includes a set of O(Zα 2 ) and O(Z 2 α 3 ) corrections [54,55,[128][129][130]. In addition, it contains the leading-logarithm renormalization group evolution from m N to m e , using the RGE kernel derived in Ref. [38] and discussed in Ref. [131]. Therefore, the standard breakdown of radiative corrections corresponds to evaluating the coupling g V at a scale µ χ ∼ Λ χ ∼ m N in Eq. (124), and then lumping the leading-logarithm RG evolution and the matrix element in δ ′ R , namely which agrees with the result compiled in Ref. [8].
From an EFT point of view aiming at describing nuclei starting from nucleon degrees of freedom, it is more natural to evolve the single-nucleon (and possibly two-nucleon, three-nucleon . . .) coupling g V all the way down the scale µ χ = m e , and only leave the evaluation of the fixed-order matrix element in δ ′ R . This can be achieved by defining the universal correction ∆ V R | EFT : and appropriately redefining the Fermi function, the outer correction δ ′ R , and the nuclear correction δ N S in such a way that they do not contain large logarithms. This requires a new EFT analysis of both δ ′ R and δ N S .

Conclusions and Outlook
In this paper, we developed a systematically-improvable top-down effective field theory framework for radiative corrections to the neutron β decay and low-energy (anti)neutrino-nucleon scattering. As a first step, we perform the matching at the electroweak scale of the Standard Model to the LEFT with effective four-fermion operators. We resum leading and next-to-leading large electroweak logarithms to all orders by evolving the semileptonic coupling constant in the LEFT from electroweak to GeV scale according to the LEFT RGEs. Next, we perform the matching to the heavy-baryon chiral perturbation theory at the hadronic scale and express the HBChPT low-energy constants in terms of non-perturbative correlation functions of quark currents. To avoid large logarithms in the evaluation of the matrix elements, we resum leading and subleading logarithms between the hadronic scale and electron-mass scale according to the RGEs in HBChPT. In this framework, all contributions from physics above the scale of the electron mass play the role of short-distance effects, which are captured by the Wilson coefficient g V . We compare our framework to the traditional current-algebra approach and find an agreement at the one-loop level. Contrary to the traditional approach, we employ dimensional regularization with minimal subtraction (MS) and specify the scale and scheme dependence in all steps of the calculation allowing us to consistently include the next-to-leading logarithms and their resummation. In our approach, the so-called DIS region of the γW box is mapped onto a contribution to the Wilson coefficient in LEFT.
In our new EFT framework, we determined the low-energy vector coupling constant g V (m e ) − 1 = (2.499 ± 0.012)%, which controls the neutron decay rate as well as low-energy (anti)neutrino-nucleon scattering, and provides the basis for one-body contributions to nuclear decays. We also extracted the CKM matrix element V ud from neutron decay measurements with our new values for the radiative corrections. An updated value |V ud | = 0.97402(42) based on the most precise determinations of the neutron lifetime and axial-vector to vector coupling constants ratio is smaller than previous results. The difference with respect to the previous analyses originates from the consistent inclusion of the next-to-leading logarithms and Coulomb corrections within our framework.
The effective field theory approach to radiative corrections in weak processes advocated in this paper can be extended to the analysis of the axial-vector coupling constant g A , which is a natural next step that will be presented in future work. The developed EFT approach can straightforwardly be applied to precise first-principles cross-section calculations in low-energy (anti)neutrino-nucleon scattering and can be extended to describe neutral-current processes with nucleons at low energies. The EFT framework can also be generalized to address radiative corrections to nuclear decays. In fact, one of the advantages of EFT is that the effective couplings g V and g A already determine the one-body inner corrections to nuclear decays at O(G F α). Consequently, matrix elements of the weak Lagrangian of Eq. (1) should be computed to O (α) in the low-energy nuclear many body theory. In this approach, the so-called "nuclear γW box" arises from contributions at scales smaller than or equal to the Fermi momentum k F , which are calculable in the nuclear EFT. Short-range physics is captured by g V , and potentially, by two-nucleon and/or few-nucleon weak operators in HBChPT, whose contributions at a given order in ϵ χ can be be estimated by the power counting in chiral EFT. The full analysis of radiative corrections to nuclear beta decay will require the development of the EFT framework for few-nucleon systems to O(G F αϵ n χ ).
In any theory, including LEFT and χPT, charge renormalization is studied in connection with the photon self-energy tensor Π µν (q 2 ) and the vacuum polarization function Π(q 2 ) defined by [136] Π µν (q) = g µν q 2 − q µ q ν Π(q 2 ). (127) Including resummed self-energy corrections, the amplitude for scattering of two charged particles is proportional to the physical (renormalization scale and scheme-independent) combination where the subscript "R" labels the renormalization scheme, α R denotes the renormalized fine-structure constant and Π R is the corresponding subtracted, UV-finite vacuum polarization function. For example, in the "on-shell" (OS) renormalization scheme the renormalized vacuum polarization function is defined by Π OS (q 2 ) = Π(q 2 ) − Π(0) ≡ ∆α(q 2 ). In this scheme, using observables at q 2 → 0, one extracts α OS = 1/137.036. This scheme can be implemented in any field theory, including χPT, LEFT, and the full Standard Model. Importantly, Eq. (128) allows one to relate α R (in any scheme and at any renormalization scale) to α OS , in terms of Π R (q 2 = 0). Eq. (128) can also be used to relate the electromagnetic couplings defined in any two renormalization schemes and even in two different EFTs.

A.1 Charge renormalization in LEFT
Throughout this paper, we use the notation α(µ) to indicate the electromagnetic coupling in LEFT, defined in the modified minimal subtraction scheme (MS).
At any value of µ < M W,Z , where the LEFT is applicable, α(µ) can be defined by its relation to α OS via Eq. (128), leading to where in the second equality we have expressed Π MS (0) in terms of Π MS (μ 2 ) and Π OS (μ 2 ), at the arbitrary scaleμ ≫ Λ QCD . In LEFT, the vacuum polarization receives contributions from charged fermions only. The contribution of charged leptons to both Π MS (μ 2 ) and Π OS (μ 2 ) = Π(μ 2 ) − Π(0) can be computed in perturbation theory. For quarks, the calculation of Π MS (μ 2 ) can be carried out in perturbation theory becauseμ 2 ≫ Λ 2 QCD , with each quark flavor of charge Q q contributing (to zeroth order in the QCD coupling α s ) The non-perturbative contributions are encoded in Π OS (μ 2 ) and can be evaluated via a dispersion relation where R(s) = σ e + e − →hadrons (s)/σ e + e − →µ + µ − (s), see for example Ref. [137]. The scheme matching conditions between the on-shell and MS couplings given in Eqs. (129) are conceptually clean and could be implemented at any value of µ < M W,Z . They are, however, not extremely practical, because the dispersive integral ∆α(μ 2 ) is usually given in the literature only atμ = M Z [56,137]. To circumvent this issue, we define the LEFT MS coupling α(µ) by relating it to the MS fine-structure constant in the full Standard Model with five quark flavors, denoted byα (5) (µ) in the PDG review on the electroweak theory [56].α (5) (M Z ) is related to α OS by an expression which is analogous to (129), with µ =μ = M Z , up to an additional contribution to charge renormalization due to the W boson [138]. Taking this into account, we find at µ = µ SM ∼ M W , The numerical valueα (5) (M Z ) −1 = 127.951(9) [56] impliesα (5) (M W ) −1 = 127.989 (9) and α(M W ) −1 = 127.936 (9). We use the latter value as the initial condition for the RGE where Q f denotes the fermion charge and n f is the multiplicity (n f = 1 for charged leptons, n f = N C for quarks). To the accuracy we work, α(µ) can be treated as continuous across heavy-fermion thresholds in LEFT.

A.2 Charge renormalization in χPT
Below the QCD scale Λ χ ∼ m N , we work with (baryon) χPT extended to include dynamical photons and leptons. For illustrative purposes, we consider SU (2) χPT with both muon and electron as dynamical fields. This is the field content relevant for matching to the LEFT at the scale Λ χ . We add the gaugekinetic terms for photons and leptons and couple them to mesons and baryons by appropriate shifts to the external field sources that appear in chiral covariant derivatives, see Eq. (27). Upon redefining the photon field A µ → (1/e)A µ , the electromagnetic coupling e only appears in the photon kinetic term and the relevant renormalized Lagrangian, written in terms of renormalized couplingê χ and counterterms Z A,χ − 1, reads The low-energy constant (LEC) h 2 was introduced in Ref. [49] in the context of SU (2) meson χPT. X 8 was introduced in Ref. [61], which extended χPT to include dynamical leptons, for the study of semileptonic processes. In χPT, the LECs contain a pure counterterm that subtracts the UV divergences of meson and lepton loops, and a finite, renormalized coupling that encodes contributions from heavy states not included in the EFT. Adopting dimensional regularization in d = 4 − 2ε and following Refs. [49,139], the generic LEC C i is written as where γ h 2 = 1/12, γ X 8 = −4/3 [61], and 1/ε = 1/ε − γ E + ln (4π). The renormalized couplings C r i (µ χ ) depend on the scale in such a way that, after including loops, the physical amplitudes are µ χ -independent. h 2 cancels divergences induced by pseudoscalar meson loops, while X 8 cancels divergences produced by loops with the electron and muon.
In the standard χPT scheme defined by Eqs. (133), we can study charge renormalization and vacuum polarization. The renormalized electromagnetic couplingα χ ≡ê 2 χ /(4π) and Πχ(q 2 ) are separately scaleindependent. In fact, the subtracted vacuum polarization is given by in terms of the pion loop, the lepton loop, and the counterterm contributions The divergent parts of the LECs, independently presented in Refs. [49,61], cancel the loop contributions from light charged particles, as expected. The µ χ dependence cancels in the sum of loop and LECs contributions, so that Πχ(q 2 ) does not depend on µ χ . Moreover, the relation between the known α OS and α χ depends on unknown LECs through Πχ(q 2 = 0) in the scheme-matching relation: A more convenient renormalization scheme, that more closely resembles the standard minimal subtraction, is achieved by rewriting Z A,χ in Eq. (133b) as where δz A,χ contains the divergent parts of the LECs, proportional to (1/ε + 1), and the fine-structure constant is redefined as α χ (µ χ ) ≡α Such a redefinition corresponds to a different choice (scheme) in separating the counterterm Lagrangian. The first of Eqs. (133) now reads (up to higher-order terms in α χ ) The corresponding finite vacuum polarization is given by whereΠ π,ℓ are obtained from Π π,ℓ by replacingα χ → α χ and 1/ε → −1 in (136a) and (136b). In this new scheme, the beta function for the renormalized coupling α χ (µ χ ) is similar to the standard minimal subtraction scheme. The running of α χ (µ χ ) is controlled by where the sum over charged leptons (ℓ) includes the electron and the muon. Another benefit of this renormalization scheme is that the relation to the on-shell coupling does not involve any unknown LECs: with Π χ (0) obtained, at one loop, from Eq. (138e). This implies The above formula accounts for the discontinuity of α χ (µ χ ) at the particle mass threshold, due to the fact that conventionally in χPT one subtracts 1/ε + 1 rather than 1/ε. In the main text of this manuscript, we imply α χ as soon as α has an argument µ χ . Numerically, using the boundary conditions described above and running α(µ) down from µ = M Z and α χ (µ χ ) up from µ χ = m e , we find 1/α χ (m p ) = 135.112 while 1/α(m p ) = 134.302.
B Factorization of the decay rate in the nonrelativistic limit In this appendix, we provide a justification of the factorized form for the electromagnetic corrections provided in Eq. (101). This form can be rigorously derived for the nonrelativistic electron, but, even for the relativistic electron, it captures the leading series in (πα/β) n and α/π(πα/β) n . As these terms are enhanced by a factor of π 2 with respect to naive two-loop corrections, they are relevant at the level of ∼ 10 −4 , and their estimate requires the evaluation of the diagrams in Fig. 3 and of next-to-next-to-leading order real-virtual and real emission diagrams. As we argue below, the first corrections to Eq. (101) are of order O(α 2 ). For E 0 − E e ≪ m e , which does not apply to neutron decay but would apply, for example, to triton decay, we could further integrate out the scale of the electron mass, and match / πEFT onto a theory with nonrelativistic electrons (NRQED). In this theory, the charged-current vector operator, with coupling constant g V in front, matches onto where ψ e denotes a nonrelativistic electron field. The matching coefficient at one loop can be extracted from the matching of heavy-light onto heavy-heavy currents performed in Ref. [123], and, in the MS χ scheme, it is given by The NRQED Lagrangian still contains various photon modes: soft (k 0 ∼ | ⃗ k| ∼ m e β), ultrasoft (k 0 ∼ | ⃗ k| ∼ m e β 2 ), and potential (k 0 ∼ m e β 2 ≪ | ⃗ k| ∼ m e β). We can integrate out soft and potential modes by matching onto potential NRQED (pNRQED) [140,141]. In this EFT, the proton and electron interact via non-local potentials, and via the exchange of ultrasoft photons. In momentum space, the potential is given at leading order by the Coulomb potential, In principle, the above potential should be evaluated using α in the Thomson limit, as the electron no longer contributes to the running of α within pNRQED. However, as we discuss below, the difference with the coupling in HBChPT, α(µ χ ), leads to very small effects numerically. Corrections to the potential are organized in powers of α and β. NLO (O(α) and O(β)) and N 2 LO (O(α 2 ), O(αβ), and O(β 2 )) corrections to the potential have been computed, both for systems of equal mass (tt, quarkonium, and positronium), and for systems with different masses (such as the hydrogen atom, or the B c meson). The results can be found in Refs. [107,[140][141][142]. No corrections to the potential appear at NLO. At N 2 LO, we can adapt the result of Refs. [141,142] to the case of QED, C F = 1, C A = 0, and n f = 0. We find where c D is the Darwin term, and c D = 1 at the order we are working. Finally, the interactions of ultrasoft photons with heavy quarks do not contribute up to N 3 LO [107,140]. This can for example be seen by the fact that one-loop real emission diagrams only contribute at O(αβ 2 ), three orders smaller than the leading order.
In the hypothetical case in which the electron emitted in neutron decay would be nonrelativistic, the decay rate could thus be expressed as with the pNRQED matrix element organized as |M pNRQED | 2 = n απ β n 1, α, β, α 2 , αβ, β 2 , . . . .
From the previous discussion, the LO is provided by the iteration of the Coulomb potential, which leads to the nonrelativistic Fermi function Since there are no NLO potentials and the ultrasoft modes only contribute at N 3 LO, M LO pNRQED does not receive any O(α) corrections. O(α) corrections are entirely contained in the matching coefficient C NRQED . Therefore, the leading (πα/β) n and the subleading α/π(πα/β) n terms are captured by the factorized expression This result was proved in the case of tt production at threshold in Refs. [104][105][106]119].
To make a connection with the relativistic expressions, we now notice that where we neglected real emission diagrams which give a contribution that goes as (E 0 −E e ) 2 /E 2 e ∼ O(β 4 ). By expressing the decay rate as dΓ n dE e = G 2 F |V ud | 2 (2π) 5 1 + 3λ 2 p e E e (E 0 − E e ) 2 [g V (µ χ )] 2 F N R (β) 1 + δ RC (E e , µ χ ) , our expression correctly reproduces the relativistic one-loop result, and, in addition, captures all subleading terms of α/π(πα/β) n . Numerically, the various contributions to ∆ f (see Eq. (110)) average to be as follows: α 2 β , α 2 , α 2 ln β = (8, 5, 2) × 10 −5 . To obtain a diagrammatic expansion of the Fermi function F N R in the electromagnetic coupling constant at the level of decay rate, we evaluate it with χPT value α(µ χ ). We also test the O(α 2 ) difference induced by computing F N R using α(µ χ ) versus α in the Thomson limit and obtain a result for ∆ f that is only 0.002% higher. Therefore, theoretical improvements (i.e., knowledge of the terms O(α 2 ln β)) would not induce much change in our results and uncertainty estimates.

C Details on the two-loop anomalous dimensions
As discussed in Sections 3.1 and 5.3, we obtain the two-loop O(α 2 ) anomalous dimensions, γ 1 andγ 1 , in the LEFT and χPT by adapting calculations in the literature. For γ 1 in the LEFT, we use Refs. [68,72] which consider the two-loop QCD anomalous dimension of a four-quark operator. As each diagram is given separately, their results can be modified to obtain the two-loop QED diagrams for the operator in Eq. (9). This involves replacing the QCD couplings, multiplicities, and color factors with their QED counterparts. The same procedure also allows us to reproduce the O(αα s ) anomalous dimension γ se . Although Refs. [68,72] provide results using one particular scheme for the evanescent operators, a = −1, the full a dependence can be recovered due to the fact that both the 1/ϵ 2 and 1/ϵ coefficients of the diagrams are given, leading to the result in Eq. (14).
To obtainγ 1 , relevant for the RGE in χPT, we use calculations for the heavy-light currents in heavyquark effective theory, [100], see also Refs. [101][102][103]. These calculations compute the two-loop anomalous dimension in QCD for a heavy-light current,QΓq, where Q denotes a nonrelativistic field, q is a relativistic particle, and Γ is an arbitrary Lorentz structure. As we argue below, the graphs in Ref. [100] can be adapted to obtain the relevant diagrams for the anomalous dimension of g V , see Fig. 3.
One way to see that the two cases are related to each other is by rewriting the weak operator in Eq. (1) using Fierz identities, where i, j run over the Dirac structures, c ij are coefficients determined by the Fierz relation, and E is an evanescent operator. The last bilinear on the right-hand side of Eq. (153) now takes the same form as the heavy-light current, where the proton and the charge-conjugated electron, e c , play the role of the Q and q fields. This means the diagrams will take the same form as for the heavy-quark calculation, with the heavy-light vertex Γ replaced byΓ j , while the neutral neutrino and neutron fields are irrelevant to the calculation. The final ingredient to show a correspondence is the fact that the loop diagrams do not depend on the Dirac structure of the vertex [103]. Since QED vertices, ∼ v µ , and propagators, ∼ i/(v · k), on the heavy-quark/proton line do not involve any Dirac structure, they cannot modify the original vertex. This is less obvious in the part of the diagrams involving the light-quark/electron line as it consists of a string of QED vertices, each of which comes with a propagator. However, one can show that, after performing the loop integrals, all gamma matrices will be contracted either with each other or with factors of v µ . This implies that the string of gamma matrices on the electron side also becomes proportional to the identity and leaves the original vertex unchanged. The independence of the gamma structure then allows us to adapt the results of Ref. [100] to obtain the diagrams with insertions of the right-hand side of Eq. (153). 10 The relevant replacements again involve replacing the QCD couplings, multiplicities, and color factors with their QED counterparts, where the appearing QED charges are now Q p and Q e + = −Q e − . The result of this procedure is given in Eq. (88). Consequently, the anomalous dimensions can be obtained by a simple substitution C F = 1, C A = 0, T F = 1 in the QCD calculation of Ref. [143]. This also allows us to obtainγ 2 as Even though this expression contains terms enhanced by π 4 , the numerical value ofγ 2 is of natural size, γ 2 ≲ 2 in / πEFT, and its contribution is beyond the required accuracy for neutron β decay.