Towards a non-perturbative calculation of Weak Hamiltonian Wilson coefficients

We propose a method to compute the Wilson coefficients of the weak effective Hamiltonian to all orders in the strong coupling constant using Lattice QCD simulations. We perform our calculations adopting an unphysically light weak boson mass of around $2~\mathrm{GeV}$. We demonstrate that systematic errors for the Wilson coefficients $C_1$ and $C_2$, related to the current-current four-quark operators, can be controlled and present a path towards precise determinations in subsequent works.


I. INTRODUCTION
Weak decays of hadrons, and in particular of mesons, play an important role in our understanding of the fundamental forces and having precise theoretical predictions to compare against the experimental results can either strenghten the solidity of the Standard Model or lead to discoveries of new physics [1].
The large scale separation between the mesons, strongly bounded particles with masses of order Λ QCD , and the weak mediators with masses around 100 GeV, is used to simplify theoretical predictions of these processes in the framework of effective field theories. By integrating out the heavier degrees of freedom, specifically the weak bosons and heavy quarks, from the Standard Model, it is possible to define a new effective Hamiltonian with new operators and new coupling constants usually called the Wilson coefficients.
The coefficients capture the effect of the weak bosons and heavy quarks that are absent from the effective field theory (EFT), making them well-suited for a perturbative calculation. Instead, matrix elements of the operators involving mesonic external states require a nonperturbative calculation. In the last decade, thanks to algorithmic and computational advances, the Lattice QCD community has been able to cover a wide range of processes involving two mesons (see e.g. Ref. [2]) and also to complete the first two-body final state decays of K → ππ [3][4][5][6].
On the other hand perturbative calculations of the Wilson coefficients have been successfully carried out up to NLO (for a comprehensive review see Ref. [7]) and in some cases up to NNLO [8]. In this work we explore the possibility of a non-perturbative method to compute the Wilson coefficients to address the perturbative uncertainty of the analytic calculations. The perturbative truncation error is traded with the statistical and systematic errors usually present in lattice calculations and the purpose of this paper is to define a methodology to obtain a precise determination of the Wilson coefficients where all uncertainties have been addressed.
The problem of defining the weak effective Hamiltonian non-perturbatively has already had some initial considerations by the authors of Ref. [9]: by separating two hadronic weak currents and studying their dependence as a function of this distance it is pos-sible to define properly normalized operators, where the effect of the Wilson coefficients has been fitted away by using their perturbative expansion.
Our approach differs from the one considered in Ref. [9] as we plan to directly determine the Wilson coefficients using gauge fixed external quark states, rather than mesons, in momentum space and not in coordinate space.
The rest of the manuscript is organized as follows: in the next section we give an overview of the main features of the EFT and we describe our strategy to measure the Wilson coefficients; in section III we show our results and address the various uncertanties of the calculation; in section IV we report our determination of the Wilson coefficients in the MS scheme and discuss the comparison against the known perturbative results, and finally we conclude and present further directions for this project.

II. COMPUTATIONAL METHOD
Before introducing the lattice observables and our main strategy, we review the most important features of the EFT.
For concreteness let us restrict to transitions among hadrons. When the weak bosons and the heavy quarks are integrated out, the leading effective field theory that arises is based on operators of dimension 6 with fourquark vertices. This can be easily seen by considering the first term in the expasion of the weak propagator in the limit m W → ∞. The full expansion however contains other terms that can be related to operators with higher dimensions, in fact the most general form of the effective Hamiltonian is In eq. (1) G F represents the dimensionful Fermi constant which is related to the SU (2) L coupling of the Standard Model g 2 , according to G F / √ 2 = g 2 2 /(8m 2 W ). V CKM i denotes a generic product of the usual CKM matrix elements that depends on the flavor structure of the process and consequently of the operators. Note that among the heavy particles that are removed from the theory there arXiv:1711.05768v1 [hep-lat] 15 Nov 2017 is also the top quark, whose presence affects the Wilson coefficients [7,8].
A reliable computation of the Wilson coefficients C i requires a good suppression of the higher dimensional terms, which is obtained by performing the calculation at energy scales sufficiently small compared to m W . Once the operators in the second sum of eq. (1) can be neglected, the Wilson coefficients are obtained by equating the amplitudes computed in the EFT with the the ones computed in the full theory, which in our case differs from the Standard Model as explained below. Previous perturbative calculations [7] were based on amplitudes computed with off-shell external massless quarks and in our calculation we proceed along these lines. As ultraviolet quantities, the Wilson coefficients are expected to be independent from the infrared regulators of the theory and the external states used in the amplitudes; checking this explicitly will be an important task of our work. Similarly we will also test the gauge invariance of our results in the limit where the amplitudes go on-shell, where any depedence on the weak gauge fixing parameter disappears.
The amplitudes of the EFT require additional renormalization conditions, due to the new divergences that appear, leading to the mixing of the renormalized operators among themselves. Hence, it is important to consider a complete basis of operators that closes under renormalization, such that The basis of operators depends on the details of the process considered. For instance, a 2 → 2 transition where the four external quarks have different flavors (e.g. c → sud) requires only two current-current operators. Instead, e.g., a basis of 7 independent operators is needed in the 3-flavor EFT (i.e. H ∆S=1 eff ) to describe the K → ππ process in the zero isospin channel, involving also disconnected contributions that are typically more difficult to compute on the lattice. Therefore, in this exploratory study we will focus only on the simpler current-current operators and consequently only on the Wilson coefficients C 1 and C 2 .
A third remark concerns the running of the Wilson coefficients, as some care is required in computing them at low scales. In fact, if we naively match the two theories at scales µ m W we may encounter large logarithms in the form of log(m 2 W /µ 2 ). Therefore, in order to cancel these terms, the matching is performed exactly at µ = m W , thus defining the so-called initial conditions of the Wilson coefficients, and their values at scales lower than m W are obtained by solving their corresponding renormalization group equations. This leads to a resummation of all logarithms of the form α n α log m 2 W µ 2 k for any power k at fixed loop order n.
The step-scaling matrix involved in the running of the Wilson coefficients is given by the ratio of Z at two different scales and is a well known and studied problem on the lattice (see e.g. Refs. [10,11]). Instead, in our work we focus on the initial conditions of the Wilson coefficients. More precisely, we compute the matching between a theory with 3 light dynamical quarks in the sea, playing the role of our 3-flavor EFT, with the full theory where we also include the W boson exchange. In this study we do not consider the problem of removing a heavy quark from the theory, which is relevant to match the Standard Model, with a top quark, onto an EFT with 5 dynamical quarks.
Finally, our last remark concerns the feasibility of this study. When we introduce the lattice spacing a as a regulator for our theory we are explicitly introducing an ultraviolet cutoff of order a −1 . The Wilson coefficients, in essence, encode the information of momenta around and above m W , making them potentially very sensitive to discretization effects. We address this question by varying both m W and the lattice spacing. However, given the current limitations on the availability of fine lattice spacings, we perform our calculations in an unphysical scenario, where we take m W ≈ 2 GeV, but where we can control the other systematic uncertainties. Nevertheless we discuss, before concluding, how our results may have an impact on the determination of the Wilson coefficients with physical values of the weak boson mass.
Once the theory is discretized the path integral can be solved using numerical simulations, limited by the present computational technologies: this translates into being able to simulate only finite quark masses and finite lattice boxes, which constitute the infrared regulators of the theory. Therefore in order to compute the Wilson coefficients the following limits need to be fulfilled where p represents the typical momentum of the external states used in the evaluation of the amplitudes. If we now consider the limit where p goes to zero, in the infinitevolume theory with mass-less quarks, contributions from higher dimensional operators should vanish. However, due to dimensional transmutation, the strong interactions possess a low intrisic scale, Λ QCD , responsible for the creation of condensates that could still contribute to eq. (1) through some operators Q (dj ) j . Nevertheless in Nature, where Λ QCD m W , these contributions should be suppressed and eventually one may neglect them. In our exploratory study we achived only Λ QCD /m W ≈ 0.2 and in the next sections we describe a strategy to quantify remaining non-perturbative contaminations and to take them into account in the systematic uncertanties.
In this work we have adopted a momentum subtraction scheme as a prescription to subtract the divergences of the theory. Alternatively, position-space techniques [12,13] could be used in a similar way to directly estimate the Wilson coefficients. In principle the window problem sketched in eq. (3) can be partly circumvented with finite-size techniques [14]: in a finite and small box mass-less quarks can be simulated and the renormalization scale, in our case p, can be identified with the box size itself, thus removing the left-most inequality of eq. (3). In addition to that, very fine lattice spacings are accessible with present resources if the physical volume is small, thus imposing a large hierarchy between m W and Λ QCD . Furthermore, imposing boundary conditions on the gauge and quark fields in the temporal direction, such as Schrödinger Functional or Twisted BC [15][16][17][18][19], has the additional advantage of producing a well-behaved perturbative expansion. Calculations with this formalism (see e.g. Refs. [20,21]) have been carried out successfully for different quantities where all systematic uncertainties have been taken into account. These approaches constitute a valid alternative to further advance this project towards a physical value of the weak boson mass.

A. The observables
As mentioned earlier we restrict ourselves to the simpler current-current operators Q 1 and Q 2 differing only in the color index routing (i, j) between the two weak currents, with ( The fields in eq. (4) are understood to reside all at the same space-time point x.
In the following we introduce a compact notation which slightly differs from the one present in the literature, e.g. Ref. [7]. Let us define the generic operators O i (x, y) as with f a generic real-valued function, γ L µ ≡ γ µ (1 − γ 5 ) and Einstein summation rule on repeated indices. In this language the choice f µν (x, y) = δ µν δ(x − y) reproduces exactly the operators Q 1 and Q 2 in eq. (4). Then let us introduce the following four-point function with a single insertion of these operators with greek and roman symbols denoting spin and color indices respectively. By inserting the operators in eq. (5) and computing the Wick contractions we end up with a Green's function that we later transform to momentum space, thus obtaining the diagram reported in Figure 1.
In our notation all the four momenta are in-coming. To simplify the notation we will omit the flavor indices in the rest of the paper since we will use degenerate quarks; following the subscript 1, 2, 3, 4 of the y coordinate or of the momenta will allow the reader to trace the flavors back.
To compute the quark propagators we invert the Dirac operator D(x, y) on plane waves with momentum p (we later refer to them as momentum sources) By saturating the x-dependence with the appropriate phase factor, we obtain the propagator in momentum space and its expectation valuẽ The periodicity of the lattice in the temporal and spatial directions imposes a constraint on the accessible momenta p µ = 2πn µ /L µ , with n µ a positive integer. Moreover the breaking of the group of continuous rotations to the hypercubic ones leads to additional discretization effects that spoil the smoothness of the propagators as a function of p. To overcome both issues we employ twisted boundary conditions (BC) in the valence sector [22][23][24], that allow to access a dense set of momenta 1 For simplicity we give below the explicit expression of Γ(O 2 ), where we omit the index contractions inside the square brackets and where we use the γ 5 -hermiticity of the Dirac operator The amputated amplitudes Λ are easily obtained by inverting the expectation value of quark propagators S(p i ) in eq. (9) At this point we define the amputated amplitudes Λ i ≡ Λ(Γ(O i )) with f µν (x, y) = δ µν δ(x−y). On the full theory side, only a single operator with color diagonal structure is needed to describe this process, namely O 2 in eq. (5) with f µν (x, y) = W µν (x, y), the tree-level propagator of the weak charged bosons in coordinate space (Euclidean metric), which we obtain by fourier transforming 2 Eqs. (6,11,12) hold in this case as well and we define Λ SM as the amputated amplitude obtained from O 2 with f replaced by the W boson propagator given above. The choice of the unitary gauge simplifies the calculation in the full theory to a single diagram. In the next section we present results also for the Feynman gauge (R ξ gauge with ξ = 1) where the contribution of the Goldstone boson needs to be included.
The δ(x − y) function in the insertion of the local operators Q i simplifies the double sum over x and y in eq. (11) to a single one. Choosing plane waves at the source of the propagators allows us to perform such a sum over the entire volume, thus sampling the operators Q i on the full lattice and reducing the statistical fluctuations of the final amplitudes. When f is replaced by the W boson propagator, the double sum needs to be performed. Also in this case the usage of momentum sources gives us the freedom (at the sink of the propagators) to evaluate both sums over x and y explicitly, thus significantly reducing the noise of Λ SM .
The only drawback of the combined use of momentum sources and twisted BCs is that a separate inversion is necessary for each momentum configuration that is considered, up to a total of 4 per configuration if we choose the four external legs with different momenta. Therefore, 2 In our calculation we use the lattice momenta apµ = 2 sin(apµ/2) in eq. (13).
in an effort to balance the cost of the inversions against the benefit of the volume average, we have explored a second strategy to sample the W boson propagator. From the simple observation that a point-source propagator can be fourier transformed to any continuous momentum, up to small finite-size errors, we have studied the possibility to stochastically sample W µν (x, y), with x and y being the sources of the inversions rather than the sinks as in the case of momentum sources. The sampling technique is essentially borrowed from the calculation of the lightby-light contribution to the anomalous magnetic moment of the muon by the RBC/UKQCD collaboration [25] and it proceeds in two steps: • for a fixed x, the sum over y in eq. (11) is achieved by randomly sampling W µν with a probability distribution falling off rapidly for large separations |x − y|; to recover the flat sum over y the appropriate reweighting factor is applied and to further decrease the cost, the hypercubic symmetries are taken into account by randomly sampling only one element per equivalence class, defined by all the points y with the same distance from x, |x − y|; this procedure defines a cloud of points stochastically sampled around the center x; • to reduce the noise of the observables a second sum over the center of the cloud, x, is performed, which is again stochastic and with a flat distribution.
Even though a continuous set of momenta is now accessible through the usage of the point sources (exceptional and non-exceptional kinematics can also be explored simultaneously), the statistical noise grows quickly: in our test we have used 40 points inside the cloud and this led to a controlled approximation of the sum over y (for the different values of the input W mass considered in this work); however the second sum over x, for which we used 16 different clouds, turns out to be the crucial one in further reducing the noise. Although the presence of a finite correlation length in our system decreases the number of useful points to reduce the noise, we have verified that the stochastic sampling leads to results that are at least 10 times noisier compared to the momentum source method, for approximately the same cost. For this reason we leave all the details of this secondary approach to the Appendix A and we concentrate in the rest of the manuscript on the analysis of the results from momentum sources. Nevertheless this technique has been very useful in an early stage of the work when we explored a vast range of momenta and kinematic configurations, on which we based our decision for the final strategy.

B. The Wilson coefficients
The appearance of new divergences in the EFT requires the introduction of additional renormalization conditions. In our study we adopt a variety of regularization independent schemes, called RI/MOM and RI/SMOM, which were introduced in Refs. [26][27][28], that are entirely defined by the choice of the external states and projectors: this translates into the momentum combination used in the calculation of the propagators and the projectors that we apply on the amputated Green's functions to obtain definite spin-color states.
In this paper we have explored two combinations of the four momenta on the external legs: the nonexceptional case called RI/SMOM where (p i + p j ) = 0 for all pairs i = j, and the exceptional one (RI/MOM) where at least one linear combination of the external momenta vanishes.
If the amplitude under study possesses the same symmetries of the light mesons, at momenta comparable with their Compton wave length, the smooth perturbative behavior of the amplitude may be spoiled by non-perturbative contaminations, referred to as Goldstone-pole contaminations [26,27]. Using non-exceptional kinematics (RI/SMOM) suppresses these unwanted effects [27,28] and we extensively discuss them in the next section.
To define specific spin-color states we use two projectors P i (i = 1, 2) that we apply to the amputated amplitudes to define the matrix such that M is invertible. To achieve this we fix the color structure of the projectors to one of the operators and we explore two options for the Dirac part, with different parities (even and odd) Eq. (16) defines the so-called γ (or γ µ ) projectors. Alternatively, the replacement γ µ → / q/|q| and γ µ γ 5 → / qγ 5 /|q| defines the / q projectors. Computing the Wilson coefficients using both parities and γ or / q Dirac structures turns out to be a crucial test of the calculation, as explained in the next section. In the rest of the paper we will refer to the two parities as VV + AA and VA + AV, with V and A labeling vector and axial Dirac structures.
The renormalization conditions that we impose on the four-quark operators read as follows with lat indicating bare lattice quantities and M tree defined by replacing the amplitudes in eq. (14) with their tree-level counterparts.
In the full theory no additional renormalization conditions are required for the projected amplitude, which we denote with W i = Tr (P i Λ SM ), beyond the usual wave function renormalization introduced above. It is important to note that since we are considering the weak theory at tree level, m W does not renormalize, since self-energy diagrams of the W boson propagator appear at higher orders in the weak coupling constant, g 2 . In the continuum theory, vector and axial current conservations (for mass-less quarks) guarantee that the same is true for g 2 as well. On the lattice however the usage of local vector and axial currents dictates the presence of the finite renormalization factors Z V and Z A , thus leading to g R 2 = Z V g 2 . Note that the V − A current can be renormalized with either Z A or Z V due to the excellent chiral properties of the Domain Wall formulation and to the employment of non-exceptional kinematics, as described later. Eventually we opt for Z V to avoid additional chiral symmetry breaking effects and a larger quark mass dependence compared to Z A . Now we have all the basic ingredients to impose the matching between the two theories We have already simplified the usual CKM factors that appear in the same form on both sides of the equation. The weak coupling constant g 2 simplifies as well with the Fermi constant G F , leaving only a factor m 2 W on the right hand side.
By expanding M RI and W RI with their corresponding bare lattice counterparts we obtain the definition of the Wilson coefficients Eq. (20) nicely separates two basic ingredients and consequently two different problems: • the bare lattice Wilson coefficients C lat i which can be used to inspect the size of the higher dimensional operators and other systematic errors, and that we compute at small momenta for a variety of external states; • the renormalization factors that we compute at high momenta and use to renormalize C lat i and eventually connect to the MS scheme, by means of the one-loop conversion matrix Z RI→MS given in the Appendix B.
Our analysis proceeds along these two steps: we first examine the dependence of the bare Wilson coefficients on the momentum scale, the quark mass and the finite box size. Then we briefly describe the results for  [34,35]. The three ensembles have been generated with the same Domain Wall discretization for the sea quarks, the Shamir formulation with fifth dimension length Ls = 16, and Iwasaki gauge action, with bare couplings g 2 0 = 6/β reported in the second column. In the third column we quote the lattice spacings measured in Ref. [36] and in the last ones we provide the lattice dimensions and the bare sea quark masses.
the renormalization matrices and discuss discretization errors, and finally present the comparison against the known perturbative results in the MS scheme.

III. RESULTS
In our study we have used the Domain Wall formalism, which retains excellent chiral properties, even at finite lattice spacing [29,30] (which become exact in the limit of infinite 5th dimension [31,32]), simplifying the renormalization pattern of the theory and suppressing the mixing among operators belonging to different representations of the chiral group [33].
We have measured the amplitudes described in the previous section on three different ensembles, reported in Table I, with 2+1 Domain Wall Fermions (DWF) in the sea. We have used a unitary setup by promoting the same discretization (Shamir DWF) also to the valence sector. Ensembles 16I and 24I have been used to study the dependence of the Wilson coefficients on the volume of the lattice. The additional ensemble 32I allows us to take the continuum limit, with two lattice spacings differing approximately by a factor of 2 in a 2 . In all our calculations we have fixed the (QCD) gauge to the Landau gauge.
Our measurements require the calculation of quark propagators for which we have used a mass (in lattice units) of 0.04 on the 24I and 16I ensemble, and 0.03 on the 32I. We utilize 27 independent configurations for both 24I and 32I, separated by 100 and 200 Molecular Dynamics Units. In our analysis we use the jackknife method with bin size of 1, after checking the stability of the error of the Wilson coefficients for larger bin sizes. On our coarser ensembles, the 16I and 24I, we measure our propagators up to momenta of O(0.8 GeV), with data points evenly spaced in p 2 . On each configuration we perform 4 different inversions at fixed p 2 for the four different legs required in non-exceptional kinematics. Moreover on the 16I and 24I we also measure the same amplitudes with momentum injected explicitly along the time direction to test finite volume effects as explained later. On the 24I we repeat those measurements for three values of the quark mass, with am = 0.02 , 0.04 and 0.08. Finally on the 32I we compute the amplitudes for four different momenta up to 0.4 GeV. Our calculations of Λ SM cover a range for am W that goes from 0.6 to 1.334 on the 24I and 0.6 < am W < 1.0 on the 32I: the masses on the 32I have been tuned to match those on the 24I according to the ratio of lattice spacings computed in Ref. [36].
When the momentum of the external quark states becomes comparable to m W , the four-quark EFT is expected to deviate from the full theory, due to the lack of higher dimensional operators that eventually become relevant. Therefore for different choices of the external states in our definition of the amplitudes, we expect to observe different behaviors with respect to the dominant scale p 2 . However, in the limit where p 2 /m 2 W 0 they should all agree and give a consistent and unique value for the Wilson coefficients, up to O(Λ QCD ) contributions. This suggests that we should be able to fit C lat i with a polynomial function of p 2 /m 2 W and we turn to this now. In the next sections we address the problems related to the usage of small momenta such as finite volume and finite quark mass effects and possible remaining non-perturbative effects of order Λ QCD /m W .

A. Fitting strategy
In order to address the size of higher dimensional operators we study the momentum dependence of the C lat i s for several choices of the external states: we adopt a nonexceptional momentum configuration given by with p 2 i = p 2 , ∀i, and transfer momentum q 2 = 2p 2 . In eq. (21) x denotes a continuous parameter obtained according to eq. (10). The exceptional momentum configuration is easily obtained by re-using the same propagator with one of the four momenta in eq. (21) on all the four legs. The momentum configuration proposed in eq. (21) leads to momentum conservation in the amplitudes and to exact equivalence between γ and / q projectors.
In Figure 2 we show the results of C lat 1 and C lat 2 as a function of p 2 /m 2 W from the 24I ensemble. We observe a good convergence of the two sets of data points (exceptional and non-exceptional) for small values of the expansion variable p 2 /m 2 W . Nevertheless higher dimensional operators are quite sizable, given the high accuracy that we are able to achieve, which forces us to explore a particularly small range of momenta and eventually The two sets of data points correspond to the exceptional and non-exceptional kinematics described in eq. (21) measured on the 24I ensemble with am = 0.04. We perform combined polynomial fits where the constant term is constrained by the two data sets, which have been obtained from / q parity odd projectors.
consider extrapolations to p 2 /m 2 W → 0. However, data obtained with non-exceptional kinematics show a milder dependence on the external momentum p 2 .
To extract the bare values of the Wilson coefficients showed in Figure 2, we adopt combined polynomial fits to the exceptional and non-exceptional points with a common constant term We base our decision for a combined fit on the fact that independent fits to the two data sets reproduce results for C lat i well compatible within 1 standard deviation (on a given fit range). To estimate the systematic uncertainties associated with these extrapolations we vary the upper limit of the fit range and the degree of the polynomial (N = 1, 2, 3) and we consider the fits with good χ 2 per d.o.f: we take as our final value the result from the highest polynomial, whose larger statistical error covers the discrepancy among the several fits 3 .
With the lattice spacings that we have studied in this work, higher dimensional operators seems to be sufficiently suppressed only for momenta around and below Λ QCD . Therefore studying the dependence on the infrared regulators of the lattice Wilson coefficients is important to control their extraction. In general we expect them to be small in C lat i due to their cancellation between W and M , especially if the quark momentum is well below m W . However to what degree this is realized in practical non-perturbative simulations is not clear a priori and we study that below.

B. Finite Volume Errors
The first infrared regulator that we investigate is the box size. For this purpose we have measured the lattice Wilson coefficients on the 16I and 24I ensembles using exceptional kinematics. In Figure 3 we present the results for C lat 2 : note that the two plots, corresponding to the two ensembles, share the same y-axis to facilitate a visual comparison. Both plots contain four sets of points obtained by combining γ and / q projectors with amplitudes where the momentum is injected along xy and time directions. The parity of the projectors corresponds to VA + AV.
For the 24I ensemble the various measurements converge to a unique point at small momenta, thus remarking the universality of the Wilson coefficients in that limit. However for the 16I ensemble a different behavior is observed: in this case only the combination of / q projectors with amplitudes with momentum along the time direction agrees with the correponding measurements in the larger volume, while the other sets of points converge to a different value. We interpret such a behavior as a finite volume error, that is largely suppressed when the momentum is injected along the time direction, which is twice as big as the spatial ones. From our numerical analysis we have noticed that the finite volume effect measured in the 16I ensemble decreases with m −1 W . In agreement with other observations presented later, this may be a consequence of a dimension-7 operator involving a non-perturbative condensate, whose strong dependence on the box volume eventually generates the spread that we observed.
A similar behavior is observed also in the other Wilson coefficient, C lat 1 , where the measurements on the 16I do not agree at small momenta, in constrast with the 24I ensemble.

C. Non-perturbative effects
In our study we have not been able to achieve a large separation between the strong and weak scales, since Λ QCD /m W ≈ 0.2. This means that our calculation, based on RI/MOM techniques, might potentially suffer from non-perturbative contaminations often re- ferred to as Goldstone-pole contaminations [26,27]: in general varying the quark mass in the measurements of the amplitudes is useful to address this issue and has been previously used to non-perturbatively subtract them away [37]. Such contaminations are present mostly in observables that share the same quantum numbers of the lighter mesons, such as the axial or pseudo-scalar bilinear operators (see for example Ref. [38]). If we extend the discussion to fermionic formulations that do not preserve chiral symmetry, such as Wilson fermions, also the mixing with wrong chiralities is allowed and at small momenta can lead to pole behavior as well. For our study we have checked how well chiral symmetry is realized by repeating some measurements of the Wilson coefficients with a larger separation of the Domain Walls along the fifth dimension: we have obtained an excellent agreement for all combinations of projectors between results obtained with the Shamir formulation with L s = 32 and L s = 16 (the latter being the same one used for the sea quarks). Hence, we expect any mixing with wrong chiralities to be predominantly of infrared origin, due to the spontaneous breaking of chiral symmetry via quark condensates, and to vanish at high momenta. Based on CPS symmetry arguments [39][40][41], we also expect the parity odd projectors to show less contaminations compared to the parity even case.
To quantify this problem we start from the difference between the Wilson coefficients computed with parity even and parity odd projectors, defined in eq. (16), for several values of the valence quark mass. In Figure 4 we show the results for the lattice Wilson coefficient C lat 2 from the exceptional momentum configuration, where we vary the quark mass and the parity of the (γ) projectors used to compute M lat ij and W lat i . For parity odd projectors we observe an excellent agreement for fixed input W mass at 1.8 GeV on the 24I ensemble. Crosses refer to parity even projectors and dots represent the opposite parity, all in the γ scheme. No difference is observed if γ projectors are replaced with the corresponding / q ones. The VV + AA projectors produce results very sensitive to the quark mass, contaminated by infrared effects proportional to 1/p 2 . A small quark mass dependence is observed in C lat 1 for parity odd projectors as well.
among the different quark masses for all points, down to the smallest momentum. On the other hand parity even projectors lead to a strong dependence of C lat 2 on the quark mass, well compatible with the expectation of a Goldstone-pole contamination: in fact, by increasing the quark mass from 0.02 to 0.08 the data is driven towards the points obtained from parity odd projectors, due to the suppression of non-perturbative effects which we expect to be proportional to (p 2 + m 2 ) −1 , with m the mass of the light state coupling to the amplitude. Changing from γ to / q projectors does not affect the general trend showed in Figure 4, as well as turning to the non-exceptional kinematics in eq. (21). For C lat 1 we draw similar conclusions on the behavior of the parity even results, with the only exception that a quark mass depedence is visible also for parity odd projectors: given the different nature and size of this Wilson coefficient it turns out to be large, on the 10% level, but only 0.015 in absolute units (with C lat 1 ≈ 0.15). We account for this effect in our systematic error.
To better understand the origin of the discrepancy between the two parities, we have projected our amputated Green's functions on the so-called wrong chiralities. For odd parity no significant mixing has been found. Instead we have measured strong contributions at small momenta from the projectors with form 1 ⊗ 1, γ 5 ⊗ γ 5 and σ µν ⊗ σ µν , which quickly vanish above 1 GeV. This provides another indication on the pollution of infrared effects in the parity even sector, as expected from CPS symmetry. In Appendix C we describe an alternative approach consisting of changing the definition of the projectors to suppress the overlap with the unwanted chiralities.
Checking the gauge independence of our calculation is a second approach that we exploit to quantify nonperturbative contributions of O(Λ QCD ). Specifically, we have computed the W lat i amplitudes in Feynman gauge. In this case the weak boson propagator simplifies to the diagonal form but the amplitude, even at tree-level, requires also the contribution of the Goldstone boson arising after the Electro-Weak Symmetry Breaking.
Goldstone boson propagator and coupling to the quark current. In the latter we have already assumed degenerate quark masses and simplified the contribution proportional to the identity.
In our case, where all the external legs have the same mass, the vertex between the Goldstone and two quarks simplifies to the form reported in Figure 5. In the limit of mass-less quarks it should vanish identically thus leaving only eq. (23) to contribute to the amplitude. This may be spoiled again by non-perturbative effects, which may be very pronounced due to the γ 5 ⊗ γ 5 structure of this diagram. Therefore studying the two contributions separately could shed some light on the size of these effects and we examine this in Figure 6, which shows the Wilson coefficients C lat i computed only from the Goldstone boson diagram and with parity even projectors. Note that the full Wilson coefficients are obtained simply by adding 4 the result of Figure 6 to the C lat i s computed from the W exchange alone. For parity odd projectors the contribution from the Goldstone boson is negligible, numerically below 10 −6 and compatible with zero within 1 sigma. For parity even projectors, in contrast, we can measure a signal from the Goldstone boson diagram. However the small prefactor (2m/m W ) 2 makes this contribution negligible also in this case, as plotted in Figure 6, although it moves both C lat 1 and C lat 2 towards the corresponding parity odd data.
Next we compare the Wilson coefficients measured in Unitary and Feynmann gauge with odd projectors exclusively, for which the Goldstone boson diagram can be neglected. The difference that we observe between the sets of data points comes from the gauge depedent part of the weak propagator in eq. (13) proportional to p µ p ν . The dependence on the weak gauge fixing condition is expected to vanish for on-shell quanti-ties, which we approach by injecting smaller and smaller momenta in the external quark legs of our amplitudes. As seen already above, in this limit non-perturbative effects produce contaminations that in this case let gauge-dependent terms survive: results between Feynman and Unitary gauge do not agree at the 2% level, even for the preferred parity odd projectors, and we may consider this difference as one source of systematic uncertainties δC gauge However, a second systematic error that we include in our calculation is taken from the difference between the even and odd projectors in Unitary gauge, δC proj In general, all the systematic uncertanties quoted in eq. (24) are estimated from the combined fits defined in eq. (22) where we exclude the points with p 2 < p 2 cut . As we increase p cut in our fit ranges, Λ QCD -type contributions are expected to vanish and this is reflected by our systematic uncertainty, which decreases also for larger values of m W . We demonstrate both behaviors in the two panels of Figure 7: the left one confirms essentially the left-most inequality in eq. (3) and the fact the natural variable to control non-perturbative effects is precisely p cut ; the right plot also confirms the general expectation that uncertainties are reduced with a larger separation between the QCD and weak scales. An additional indication of the non-perturbative origin of these contributions resides in the linear behavior in 1/m W showed in Figure 7: the presence of a condensate (a pure nonperturbative effect) could be captured by a higher dimensional operator in the form of whose functional form describes well the data.
To summarize, we have collected evidence, from the spread between different projectors and gauge fixing conditions to the volume dependence, that lead to infrared contaminations vanishing with 1/m W . In the right plot of Figure 7 we also demonstrate how a future calculation with m W ≈ 4 GeV would significantly improve the precision well below the percent level.
Finally, for the central values of the Wilson coefficients we use a p cut of 0.3 GeV for the 24I ensemble and 0.24 GeV for the 32I ensemble and we perform cubic and quadratic fits respectively according to eq. (22). We use amplitudes computed in unitary gauge projected on the odd sector with / q type projectors. No significant effect, beyond the ones already considered, is obtained from the difference with γ projectors.

IV. NON-PERTURBATIVE DETERMINATION OF THE WILSON COEFFICIENTS
The remaining systematic uncertanties that we need to address are related to discretization errors for which we need the renormalization factors. In the following we present non-perturbative results for the following values of the W boson mass: 1.4, 1.8, 2.1 and 2.4 GeV.

A. Renormalization factors
As before twisted boundary conditions turn out to be very useful, as they give us the possibility to tune the momentum to the desidered value. In our computation of the renormalization factors we have adopted the conventional RI/SMOM scheme described in Ref. [42]. Here a momentum 2p leaves the operator, as opposed to eq. (21), and only two inversions per configuration are required with momentum p 1 = (x, −x, 0, 0) and p 2 = (0, x, −x, 0) (with p 3 = p 1 and p 4 = p 2 ). Since the renormalization conditions in eq. (17) are imposed in the chiral limit we repeated the measurements for two values of the quark mass (am = 0.04 and 0.02) and we extrapolate linearly to zero quark masses. We computed the quark propagators such that |ap| ≈ am W , to obtain the Wilson coefficients renormalized at a scale coinciding with the mass of the weak boson utilized in their definition: in this way we are reproducing the so-called initial conditions, where large logarithms are avoided as explained in section II.
Finally, we deal with the renormalization of the local currents on the lattice. Therefore from the same set of propagators we also compute the following quantities After the usual amputation with the inverse propagators,  GeV   FIG. 7. Left: δC2, as a function of pcut, decreases as more points at small momenta are excluded from the fit, suppressing non-perturbative effects. Right: measured difference of C2 from odd and even projectors in the γ scheme as function of the W boson mass. The lines are linear fits forced to pass through the origin with excellent χ 2 . The non-perturbative contaminations which we quantify with δC proj i given in eq. (24) seem to vanish only with one inverse power of mW. In the plot we present the data for the exceptional kinematics, but practically no difference is observed in the non-exceptional case.
we impose the renormalization conditions In principle the appropriate renormalization condition should involve the Γ V −A µ , obtained by replacing γ µ with γ L µ in the first line of eq. (26). However we have explicitly checked from our measurements that Tr [Λ A µ γ µ ] 10 −3 (the same holds if A → V and γ µ → γ µ γ 5 ) which means that for the choice of projectors in eqs. (27,28) only Z V alone (or Z A ) renormalizes the weak coupling g 2 , thus leading to eq. (19). Replacing γ µ → ( / qq µ )/q 2 and γ µ γ 5 → ( / qγ 5 q µ )/q 2 in eqs. (27,28) defines the / q scheme as before. For simplicity, we do not combine Z q and the four-quark Z matrix from two schemes.
Let us briefly examine the properties of the renormalization factor of the Wilson coefficients in eq. (20) which we explicitly expand in terms of lattice observables for the reader's convenience (in the γ scheme) Firstly the quark mass dependence is negligible, as the four elements of the matrixZ slightly differ when the quark mass is changed from 0.04 to 0.02. Nevertheless we take this into account by linearly extrapolating to the chiral limit. Secondly we examine the difference between the renormalization matrixZ computed from two intermediate schemes, defined by the two classes of projectors, γ and / q. As we change the renormalization scale from 1.4 to 2.4 GeV the conversion factor Z RI→MS , which is only perturbative, becomes more precise. Nonetheless the ratioZ MS,γ [Z MS,/ q ] −1 significantly deviates from the identity matrix, the signature of missing terms of size α 2 s . The largest departures from 1 are observed in the diagonal terms, specifically from 23 to 18% for the (1, 1) element (that contributes mostly to C MS 1 ), and from 5 to 4% for the (2, 2) element that defines C MS 2 . These higher order effects are relatively large and can be reduced only by pushing the calculation to higher renormalization scales, where α s becomes smaller. Thirdly, at these relatively high momenta, both parity even and odd projectors give results well compatible with each other.

B. Discretization Errors
The lattice cutoff is the main limitation of the current results preventing the weak boson mass from being well separated from lower QCD scales. As we increase it, larger discretization errors are expected and values of m W of the order of the cutoff may sound already dangerous: in Figure 8 we demonstrate that the parameter space explored in this work is still in a region where discretization errors are reasonably under control.
By plotting the continuum extrapolations of C RI 1 for the four values of m W considered and normalized at the finer ensemble 32I, we show how the size of the a 2 coefficient slightly changes with m W , making the m Wdependence practically negligible compared to the overall size of cutoff effects. The scaling violations that we observe for C RI 1 range from 10 to 17%; instead they are only at 1% level for C RI 2 , where the differences among the several values of m W are also irrelevant. The different magnitudes can be easily understood from the fact that in the free theory C 1 = 0 and C 2 = 1, meaning that discretization errors start at order α s , which on our lattices is approximately 0.3. Nevertheless we perform extrapolations in a 2 rather than α s (a)a 2 , since the numerical change is irrelevant.

C. Final results and discussions
The last aspect of our work consists of changing renormalization scheme from RI to the more common MS. This is the only step where perturbation theory enters in our calculation since the conversion matrices are known to 1-loop and reported in Appendix B. In Figure 9 we plot our results together with the known 1-loop analytic formulae that we take from Ref. [7]. For the latter we use α s from the 1 and 4-loop β function with 3 flavors and Λ MS N f =3 = 341(12)MeV taken from Ref. [20]. The error from the Λ parameter is propagated to the analytic Wilson coefficients and it is represented by the two shaded regions.
The non-perturbative results include the various systematic errors described in the previous sections. The systematic uncertainty associated with the perturbative error of the RI → MS conversion, which is an O(α 2 s ) effect, can be read off from the difference between the two intermediate RI schemes studied in the paper, the γ and / q scheme. Such a difference turns out to be relatively large, following the previous discussion on the renormalization matrixZ. We remark that the RI → MS conversion matrix for the / q scheme given in Appendix B contains a large one-loop coefficient, whereas the γ scheme shows a much better convergence with a very small correction from RI to MS.
In Figure 9 we compare our results for the C MS i s against the perturbative ones as a function of the W mass. Recall that we are focusing on the initial conditions of the Wilson coefficients, which means that their renormalization scale coincides with m W . The discrepancy between our values and the 1-loop results decreases for larger values of the weak mass as expected from the running of the strong coupling constant, since the perturbative expansion of the Wilson coefficients reads with tree-level values k Given the limitations of the RI → MS conversion factor to one-loop, even with higher precision, fitting our MS results beyond O(α s ) would be incorrect. In fact, we turn now to our non-perturbative determination by sticking to the RI scheme, where we can convert the analytic results of the Wilson coefficients without loss of generality. In Table II we report various results from different polynomial fits to the γ scheme, following eq. (31), which we also augumented with a term k Λ /m W .    Table II show the potential of the method: the possibility to extract higher loop coefficients can be used to bound the error of C i (80 GeV) by considering, for example, the difference between the known 1-loop result and our 2-loop prediction. This approach has the potential to reduce the current systematic uncertainties where these Wilson coefficients are used and are relevant, such as the real part of the amplitudes of K → ππ for isospin 0 and 2 channels. One caveat that we need to remember is that any prediction of coefficients beyond one-loop depends on the number of flavors used in the simulations. For this reason we intend to continue our study by reusing the methodology developed in this paper on the finer ensembles generated by the RBC/UKQCD Collaboration [43] with 3 and 4 active flavors in the sea. Such a calculation will provide multiple benefits, from the possibility to push m W up to 4 GeV and substantially reduce the systematic uncertainties (see Figure 7), to controlling the flavor dependence of the coefficients of higher loops and reducing the systematic error from intermediate RI scheme, thus providing a solid prediction in the MS scheme. Finally, future extensions of this work will include the top quark contribution.

V. CONCLUSIONS
In this paper we have presented a method to compute the initial conditions of the Wilson coefficients to all orders in the strong coupling constant and leading order in g 2 . We have described the limitations of our exploratory study, mostly related to the presence of large infrared and non-perturbative effects. By looking at different observables we quantified those effects and took them into accout in our systematic uncertainties, that dominate the final errors.
Nevertheless, we have demonstrated how precise statistical results can be achieved with the combined use of momentum sources and twisted boundary conditions. Therefore we expect to obtain excellent results with the next iteration of this calculation repeated on finer lattices, where the systematic errors will significantly decrease.
Despite the limitations imposed by the lattice cutoff, we observed reasonable scaling violations for the values of m W that we explored. Moreover we discussed a strategy to extend the relevance of our study to place a bound on the error of the Wilson coefficients for the physical scenario with m W ≈ 80 GeV.

VI. ACKNOWLEDGEMENTS
We would like to thank our RBC and UKQCD collaborators for helpful discussions and support. M.B. is particularly indebeted to N. Christ for a critical reading of the manuscript and to T. Izubuchi and P. Boyle for many stimulating discussions. Computations for this work were carried out at the Fermilab cluster pi0 as part of the USQCD Collaboration, which is funded by the Office of Science of the U.S. Department of Energy. M.B., C.L. and A.S. were supported by the United States Department of Energy under Grant No. de-sc0012704. In addition C.L. is supported in part through a DOE Office of Science Early Career Award. In this Appendix we present further details on the alternative stochastic sampling method for the W boson propagator described in the main text. The key ingredient is the possibility to approximate the propagator in momentum space x e ip(y0−x) D −1 (x, y 0 ) with pointsource propagators transformed to any continuous momentum, up to finite-size errors It is easy to check that if the momentum is an allowed fourier mode the equation above amounts to a simple translation of the source to the origin. However if the momentum is not quantized, S app approximates well the propagator in infinte volume: the mass gap of QCD ensures that finite volume errors are exponentially small.
The second feature of the stochastic sampling relies on the approximation of the sum over the W propagator. Let us fix one end to the origin and call r the second end. Due to the fast fall-off of the weak boson propagator, only small separations contribute to the signal. Hence, we start by dividing all distances r in classes defined by their absolute value |r| with multiplicities d |r| .
Then we randomly choose one representative per class and we cover all distances up to a certain cutoff R inner . Up to lattice symmetries this sum is exact.
Finally we sample the remaining classes starting from R inner with probability p(r) ∝ |r| −3 and we evaluate the appropriate reweighting probability w(r) to obtain the flat sum again. With am W = 1.0 on the 24I ensemble, sampling 30 classes exactly (below R inner ) and 10 classes stochastically (above R inner ) produces controlled approximations with stochastic errors around 1%. This sampling strategy defines what we call a "cloud" of points around the origin.
As explained already in the text, the problem resides in the second sum over different clouds, each centered around a randomly chosen point. Although the noise of the final amplitudes and Wilson coefficients scales with the number of clouds, compared to the momentum sources, which amount to summing all of them, it is still from 5 to 10 times larger, when going from high (1 GeV) to small momenta. The access to all momenta has the additional advantage of averaging over different orientations, such as (p, 0, 0, 0), (0, p, 0, 0), (−p, 0, 0, 0), etc. . . , an effect already included in the numerical factors quote above.
To reduce the cost of this method we have also employed an All-Mode-Averaging [44] strategy by comput-ing the point-source propagators of the clouds sloppily and adding the corresponding correction term, which we tuned to be well below the statistical error throughout the entire range of momenta. In our final measurements we have used up to 16 total clouds, each containing 40 points as described above. For the same cost we have obtained the data points showed in Figure 2 and in the right panel of Figure 3. µ u i )(ū j γ L µ d j ) , Q Then we substitute u with c whenever it appears in eq. (B1) and we compute the Green's functions, similarly to eq. (6), i )] αβγδ abcd = s α a u γ c Q With the substitution u → c advocated above, we are explicitly eliminating all the disconnected diagrams of the three-flavor theory. At this point it is easy to demonstrate that Γ(Q where the operators on the r.h.s. correspond to eq. (4). In the three-flavor theory the 10 operator basis is redundant, leading to a smaller basis, with linear independent operators, often called the chiral basis (we consider again solely the linear combinations involving the three operators in eq. (B1)) 3 , Q 2 = 1 5 (2Q 1 − 2Q 3 ) , Q 3 = 1 5 (−3Q 1 + 3Q 3 ) .

(B4)
If we now combine eq. (B4) with eq. (B3) we can relate the chiral basis to our two operators as follows In Ref. [42] the conversion matrices Z RI→MS are given in the chiral basis, which is the reason behind the introduction of the matrix R. We stress that Q 1 transforms under the (27,1) representation of the chiral group, whereas Q 2 and Q 3 under the (8,1) representation. This prevents any mixing between the two sectors once the penguin operators are discarded, as in our case. Finally let us introduce three ad hoc matrices T such that To relate the renormalization factors of Ref. [42] to our specific case we need to recall the renormalization conditions in eq. (17) and the definition of the matrix M . Similarly to eq. (14) we can define a matrix M for the chiral basis M ij = Tr (P j Γ(Q i )) = (RM U T ) ij , with the matrix U being the rotation of the projectors P i = U ij P j . Similar to the matrices T (i) , we introduce three ad hoc matrices S (i) such that their product with U returns the 2 × 2 identity. Starting from the usual renormalization conditions for M in the RI sense, with the help of eq. (B8) and a few algebraic steps we finally obtain The universality of γ scheme is such that for any choice of U (and consequently of S (i) ) eq. (B9) is always valid for i = 1, 2 and 3 and we verified that. Instead for / q projectors the situation is different, since a naive choice of projectors can lead to mixing between the (27,1) and (8,1) operators even with the explicit suppression of the penguin diagrams. An accurate choice of projectors that forbids this is provided in the Appendix B of Ref. [42]. For simplicity we consider eq. (B9) for T (2) which corresponds to selecting Q 2 and Q 3 alone.
The final results for the conversion matrices reported below have been obtained from Table III of Ref. [42] for the γ scheme and from Table XIII for the / q scheme where in both cases we have set the penguin contributions to zero As explained in the main text, the bare lattice Wilson coefficients can be obtained from any reasonable choice of projectors in the definition of M and W . Below we present an alternative definition of the projectors defined in eq. (16) which holds both for γ and / q types. Since the aim of this digression is to reduce the discrepancy of the results obtained in the parity even and odd sectors, we focus on the former, the most problematic one. Let us introduce the "subtracted projectors" with i, j = 1, 2 referring as usual to the color diagonal and mixed case and the index α labelling the wrong chiralities that could potentially mix in the parity even sector, namely SS ± PP, VV − AA and TT. The coefficients b α ij are fixed by solving the following system Tr (P i Q α j ) = 0 , α = SS ± PP, VV − AA, TT , which minimizes the overlap of the projectors on the operators with the corresponding wrong chiralities (the index i in the operator refers again to the color structure). From a numerical experiment we found that the Dirac structures SS − PP and VV − AA contribute only to the noise without producing any effect. Instead the remaining two chiralities have a beneficial effect on the Wilson coefficients: in particular the SS + PP alone reduces the discrepancy for both C lat 1 and C lat 2 by 20-30 %; the inclusion of the tensor structure further increases the agreement with the odd projectors to less than 1 sigma for C lat 2 , but has the opposite effect on C lat 1 . In conclusion, the best result for the subtracted even projectors has been obtained with α = SS + PP alone. Since the benefit was marginal, in our estimate of δC proj i we decided to use the un-subtracted projectors, giving a more conservative and safer error.