Direct CP violation and the ΔI=1/2 rule in K→ππ decay from the Standard Model
Ryan Abbott, Thomas Blum, Peter A. Boyle, Mattia Bruno, Norman H. Christ, Daniel Hoying, Chulwoo Jung, Christopher Kelly, Christoph Lehner, Robert D. Mawhinney, David J. Murphy, Christopher T. Sachrajda, Amarjit Soni, Masaaki Tomii, Tianle Wang
CCERN-TH-2020-058, MIT-CTP/5197
Direct CP violation and the ∆ I = / ∆ I = / ∆ I = / rule in K → π π K → π π K → π π decay from theStandard Model R. Abbott, T. Blum,
2, 3
P.A. Boyle,
4, 5
M. Bruno, N.H. Christ, D. Hoying,
3, 2
C. Jung, C. Kelly, C. Lehner,
7, 4
R.D. Mawhinney, D.J. Murphy, C.T. Sachrajda, A. Soni, M. Tomii, and T. Wang (RBC and UKQCD Collaborations) Physics Department, Columbia University, New York, NY 10027, USA Physics Department, University of Connecticut, Storrs, CT 06269-3046, USA RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA Brookhaven National Laboratory, Upton, NY 11973, USA SUPA, School of Physics, The University of Edinburgh, Edinburgh EH9 3JZ, UK Theoretical Physics Department, CERN, 1211 Geneve 23, Switzerland Universit¨at Regensburg, Fakult¨at f¨ur Physik, 93040, Regensburg, Germany Center for Theoretical Physics, Massachusetts Institute of Technology, Boston, MA 02139, USA School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK (Dated: April 21, 2020) a r X i v : . [ h e p - l a t ] A p r ABSTRACTWe present a lattice QCD calculation of the ∆ I = / K → ππ decay amplitude A and ε (cid:48) , the measure of direct CP-violation in K → ππ decay, improving our 2015 calculation [1]of these quantities. Both calculations were performed with physical kinematics on a 32 × a − = . ( ) GeV. However, the current calcu-lation includes nearly four times the statistics and numerous technical improvements allowingus to more reliably isolate the ππ ground-state and more accurately relate the lattice operatorsto those defined in the Standard Model. We find Re ( A ) = . ( . )( . ) × − GeV andIm ( A ) = − . ( . )( . ) × − GeV, where the errors are statistical and systematic, re-spectively. The former agrees well with the experimental result Re ( A ) = . ( ) × − GeV. These results for A can be combined with our earlier lattice calculation of A [2] toobtain Re ( ε (cid:48) / ε ) = . ( . )( . )( . ) × − , where the third error represents omitted isospinbreaking effects, and Re ( A ) /Re ( A ) = . ( . )( . ) . The first agrees well with the experimen-tal result of Re ( ε (cid:48) / ε ) = . ( . ) × − . A comparison of the second with the observed ratioRe ( A ) / Re ( A ) = . ( ) , demonstrates the Standard Model origin of this “ ∆ I = / I. INTRODUCTION
A key ingredient to explaining the dominance of matter over antimatter in the observable uni-verse is the breaking of the combination of charge-conjugation and parity (CP) symmetries. Theamount of CP violation (CPV) in the Standard Model is widely believed to be too small to explainthe dominance of matter over antimatter, suggesting the existence of new physics not present inthe Standard Model. CPV in the Standard Model is highly constrained, requiring the presence ofall three quark-flavor doublets and described by a single phase [3]. These properties imply thatthe “direct” CPV in K → ππ decays is a highly suppressed O ( − ) effect in the Standard Model,making it a quantity which is especially sensitive to the effects of new physics in general, and newsources of CPV in particular.Direct CPV was first observed in K → ππ decays by the NA48 (CERN) and KTeV (FermiLab)experiments [4, 5] in the late 1990s, and the most recent world average of its measure is Re ( ε (cid:48) / ε ) = . ( . ) × − [6], where ε is the measure of indirect CPV ( | ε | = . ( ) × − ). However,despite the impressive success of these experiments, it was only recently that a reliable, first-principles Standard Model determination of ε (cid:48) that could be compared to the experimental valuebecame available. This is due to the presence of low-energy QCD effects that are difficult to modelreliably.Lattice QCD is the only known technique for determining the properties of low-energy QCDfrom first principles with systematically improvable errors. In this regime the high-energy physicsis precisely captured by the ∆ S = H W = G F √ V ∗ us V ud ∑ i = [ z i ( µ ) + τ y i ( µ )] Q i ( µ ) , (1)where the Fermi constant G F = . × − GeV − , V q (cid:48) q is the Cabibbo-Kobayashi-Maskawamatrix element connecting the quarks q (cid:48) and q and τ = − V ∗ ts V td / V ∗ us V ud . The quantities z i and y i arethe Wilson coefficients that encapsulate the high-energy behavior, and which have been computedto next-to-leading-order (NLO) in QCD perturbation theory and to leading order (with some im-portant NLO terms) in electroweak perturbation theory [7], in the MS scheme as a function of thescale µ . The task of the lattice calculation is to determine the matrix elements (cid:104) ππ | Q i ( µ ) | K (cid:105) ofthe weak effective operators Q i renormalized in a scheme which can be defined non-perturbatively.A further perturbative calculation is subsequently necessary to match such matrix elements tothose in the MS scheme. Conventionally, as shown in Eq. (1), the weak Hamiltonian is expressedin terms of 10 operators { Q i } ≤ i ≤ (as defined, for example, in Eqs. (4.1) – (4.5) of Ref. [7] )that are linearly dependent due to the Fierz symmetry. More convenient for our purposes is a sec-ond, 7-operator “chiral” basis [8] { Q (cid:48) j } j = , , , , , , in which the operators are linearly independentand transform as irreducible representations of SU ( ) L × SU ( ) R . The relationship between thesebases is discussed in more detail in Sec. VI B.For an isospin-symmetric lattice calculation it is convenient to formulate the K → ππ matrixelements in terms of two amplitudes, A and A , where A I = (cid:104) ( ππ ) I | H W | K (cid:105) and the subscriptindicates the isospin representation of the two-pion state. These correspond to ∆ I = / ∆ I = / ε (cid:48) can be obtained directly as ε (cid:48) ε = i ω e i ( δ − δ ) √ ε (cid:20) Im ( A ) Re ( A ) − Im ( A ) Re ( A ) (cid:21) (2)where δ I are the ππ scattering phase shifts and ω = Re ( A ) / Re ( A ) . Note that the effects ofisospin breaking and electromagnetism are not included in our simulation and are instead treatedas systematic errors as discussed in Sec. VIII D.In 2015 the RBC & UKQCD collaborations published [1] the first lattice calculation of A using 216 lattice configurations with a 32 ×
64 volume, an inverse lattice spacing of a − = . ( ) GeV, and with physical kinematics. We found Re ( A ) = . ( . )( . ) × − GeV and Im ( A ) = − . ( . )( . ) × − GeV, where the parentheses contain the statis-tical and systematic errors, respectively. Within the uncertainties, the former agrees with theexperimental result of Re ( A ) = . ( ) × − GeV, and the latter, combined with the ex-perimental value of Re ( A ) and the result of our previous calculation of A [2], gives Re ( ε (cid:48) / ε ) = . ( . )( . ) × − , which is 2 . σ below the experimental value.In order to obtain on-shell kinematics, i.e. to ensure that E ππ , the energy of the two-pionfinal state, satisfies E ππ = m K , we exploit the possibility of choosing appropriate spatial boundaryconditions. With periodic boundary conditions for all the quarks, the ground state of the two-pionfinal state corresponds to E ππ = m π , with each of the pions at rest, and the state with E ππ = m K appears as an excited state. We would therefore need to resort to multi-state fits to rather noisy datain order to obtain the physical amplitudes. The change in the finite-volume corrections inducedby modifying the boundary conditions is exponentially small [9, 10] or else accounted for by theL¨uscher and Lellouch-L¨uscher [11, 12] prescriptions with minor alterations [13].In our calculation of A [2, 14, 15] we employ antiperiodic spatial boundary conditions (APBC)for the down quark in some or all directions, which results in the charged pions also obeying cor-responding antiperiodic boundary conditions. The momenta of the charged pions are thereforediscretized in odd-integer multiples of π / L in these directions, where the spatial volume of the lat-tice is V = L . Since only the spectrum of the charged pions is changed by the APBC, we compute K + → π + π + matrix elements of operators which change I z , the third component of isospin, by3/2 and then use the Wigner-Eckart theorem to obtain the physical K + → π + π amplitude whichis proportional to A . Note that in order to ensure that E ππ = m K , L must be appropriately tuned.The technique of using APBC on the down quark naturally breaks the isospin symmetry. Forthe ∆ I = / K + → π + π + matrix element is the only doubly-charged two-pion state and there-fore cannot mix with other ππ states due to charge conservation. However, the final state inthe ∆ I = / | π + π − (cid:105) and | π π (cid:105) states. Thus the breaking of isospin symmetry at the boundaries results in different energies forthe | π + π − (cid:105) and | π π (cid:105) states and the APBC technique cannot be used. As a result, for the calcula-tion of the ∆ I = / ◦ isospin rotation about the y-axis, ˆ G = ˆ Ce π ˆ I y .The charged and neutral pions are both odd eigenstates of this operation, hence when applied asa boundary condition all pion states again become antiperiodic in the spatial boundary. Whilemore general than the APBC approach, the use of GPBC introduces a number of technical andcomputational difficulties that we discuss in Ref. [10] and below.Note that due to the ππ interaction being repulsive in the I = I = ππ energies in these two representations differ at fixed lattice size andit is therefore not possible to use the same ensemble to measure both the ∆ I = / ∆ I = / A and will combine it with the results for A given in Ref. [2].Among the necessary ingredients in the lattice calculation of the K → ππ matrix elements arethe two-pion energies E ππ and the amplitudes (cid:104) ππ | O ππ | (cid:105) , where O ππ is an interpolating opera-tor that can create the required two-pion state from the vacuum. These quantities are determinedby correlation functions describing the propagation of the two-pion state. The matrix elements (cid:104) ππ | Q i | K (cid:105) are obtained from K → ππ correlation functions in which the Euclidean time depen-dence is exponential in E ππ , and the amplitudes corresponding to the annihilation of the two-pionstate (and the creation of the kaon state) have to be removed. From the measurement of E ππ andusing the L¨uscher formula [11], we also determine the s-wave isospin 0 ππ -phase shift, δ ( m K ) ,which enters the expression relating the matrix elements to ( ε (cid:48) / ε ) , Eq. (2). The derivative of thephase-shift with respect to the energy is additionally required to determine the power-like (i.e. non-exponential) finite-volume corrections through the Lellouch-L¨uscher formula [12] (cf. Sec. VI A).In the 2015 calculation we obtained δ ( E lat ππ ≈ m K ) = . ( . )( . ) ◦ , substantially smaller thanthe dispersive result [16].The observation of a discrepancy from the predicted phase shift increased our motivation toextend the earlier calculation by increasing the statistics and using more sophisticated methodsto better analyze the I = E ππ ; the result was consistent between one-and two-state fits to our data (i.e. whether we assumed that just the ground-state was propagatingor allowed for a contribution from an excited state) and was also independent, within the uncer-tainties, of the time separation between the insertion of the creation and annihilation operators (the O ππ introduced above). Nevertheless, we considered the best explanation for the discrepancy to becontamination from one or more excited states whose contribution with increasing time is maskedby the rather rapid reduction in the signal-to-noise of our data. Therefore, in addition to increas-ing our statistics by more than a factor of 3, we have introduced two additional ππ interpolatingoperators. For our original calculation we used a ππ operator comprising two quark bilinear op-erators that create back-to-back moving pions of a particular momentum. Alongside this operator,which we label ππ ( ) , we have now added a scalar operator σ = √ ( ¯ uu + ¯ dd ) , and an operatorcreating pions with larger relative momenta that we label ππ ( ) . Here the number appearingin the parentheses of the ππ ( . . . ) operators is related to the components of the pion momentum inlattice units: ( xyz ) → ( ± x , ± y , ± z ) π / L (the total ππ momentum is zero in all cases). Here andfor the remainder of this document we will assume the lattice size L to be in lattice units unlessotherwise stated. All three operators, once suitably projected onto a state that is symmetric undercubic rotations, have the same quantum numbers as the s -wave I = × ππ two-point correla-tion functions in which the two-pion states are created or annihilated by one of these three opera-tors, results in a substantial reduction in the statistical and systematic errors. We find that, once theexcited states are taken into account, the resulting I = ππ -scattering phase shift at E lat ππ = . δ ( E lat ππ ) = . ( . )( . ) ◦ , where the errors are statistical and systematic, respectively.This significant increase in our result for δ ( E lat ππ ) brings us into much closer agreement with thedispersive prediction, which at our present value of E lat ππ is δ ( E lat ππ ) disp = . ◦ , obtained usingEqs. 17.1-17.3 of Ref. [16] with m π = . ∆ I = / K → ππ matrix elements obtained from our expanded data set of 741 measurements, using allthree ππ interpolating operators.In this analysis we also include an improved non-perturbative determination of the renormal-ization factors relating the bare matrix elements in the lattice discretization to those of operatorsrenormalized in the RI-SMOM scheme (see Sec. V). Perturbation theory is then required to matchthe operators renormalized in the RI-SMOM scheme to those in the MS scheme in which the Wil-son coefficients have been computed. This calculation utilizes step-scaling to raise the matchingscale from 1.53 GeV to 4.01 GeV, significantly reducing the systematic error associated with theperturbative matching.Throughout this document results are presented in lattice units unless otherwise stated.While the current paper is intended to be self-contained it should be viewed as the third in aseries of three closely related papers. The first of these is Ref [10] which gives a detailed discus-sion of the implementation and properties of lattice calculations which impose G-parity boundaryconditions. The second paper is Ref. [18] in which the same ensemble of gauge configurationsand many of Green’s functions used in the current paper are analyzed to study ππ scattering. Thissecond paper contains the two-pion, finite-volume energy eigenvalues from which the I = I = ππ scattering phase shifts are derived as well as the matrix elements of the ππ interpolatingoperators between the corresponding energy eigenstates and the vacuum which are used in thecurrent paper.For the convenience of the reader we summarize the primary results of this work in Tab. I. Forfurther discussion we refer the reader to Sec. VIII. It is important to stress that the results anduncertainties in Tab. I have been obtained by combining a number of elements. The major directcontribution from this work is the evaluation of the matrix elements (cid:104) ( ππ ) I = | Q i | K (cid:105) in isosym-metric QCD, with the operators renormalized in the RI-SMOM ( / q , / q ) scheme (see Tab. XXVII),with the lattice systematic uncertainties carefully estimated (see Sec. VII). These matrix elementsare combined with the perturbatively-calculated Wilson coefficients in the MS scheme and the per-turbative matching of the matrix elements from the RI-SMOM to MS schemes with estimates ofthe corresponding systematic uncertainties. If and when these perturbative uncertainties, as wellas those in the CKM matrix elements and isospin breaking, are reduced then the matrix elements Quantity ValueRe ( A ) × − GeVIm ( A ) -6.98(0.62)(1.44) × − GeVRe ( A ) / Re ( A ) ( ε (cid:48) / ε ) TABLE I: A summary of the primary results of this work. The values in parentheses give thestatistical and systematic errors, respectively. For the last entry the systematic error associatedwith electromagnetism and isospin breaking is listed separately as a third error contribution.in Tab. XXVII can be used to improve the precision in the determination of A .The layout of the remainder of this paper is as follows: In Sec. II we introduce our latticeensemble and give a general overview of our measurement techniques. In Sec. III we discuss andpresent results from fits to the single-pion, two-pion and kaon two-point correlation functions, thevalues of which are required as inputs to the fits of the K → ππ three-point correlation functionsfrom which the matrix elements of the bare lattice operators are determined. In Sec. IV we discussthe measurement of these three-point functions and provide the results from the fits. In Sec. Vwe discuss our procedure for the non-perturbative renormalization of the operators Q i , the resultsof which are combined with the matrix elements of the bare lattice operators and other inputs todetermine A and ε (cid:48) / ε in Sec. VI. We follow this by a detailed discussion of the systematic errorsin Sec. VII and present our final results for the matrix elements, decay amplitudes, and ε (cid:48) / ε ,together with a discussion of the ∆ I = / II. OVERVIEW OF MEASUREMENTS
In this section we provide an overview of the calculation, including information on the ensem-ble and the measurement techniques.
A. Gauge ensemble
For this calculation we employ a single lattice of size 32 ×
64. We utilize 2 + L s =
12 and M¨obius parameters b + c = /
12 and b − c = × − and 0.045, respectively. We use the Iwasaki+DSDRgauge action with β = .
75, corresponding to an inverse lattice spacing of a − = . ( ) GeV.The dislocation suppressing determinant ratio (DSDR) term [19] is a modification of the gaugeaction that suppresses the dislocations, or tears in the gauge field that enhance chiral symmetrybreaking at coarse lattice spacings. This enables the calculation to be performed with larger latticespacings, and hence larger physical volumes, at fixed computational cost, ensuring good controlover finite-volume systematic errors. We use G-parity boundary conditions (GPBC) in three spatialdirections in order to obtain nearly physical kinematics for our K → ππ decays.The lattice parameters are equal to those of the 32ID ensemble documented in Refs. [20, 21]except for the boundary conditions and that we now simulate with a lighter, physical pion mass of142 MeV versus the 170 MeV pion mass of the 32ID ensemble. This allows the use of existingmeasurements such as the lattice spacing, and also enables the computation of the non-perturbativerenormalization factors in a regime free of the complexities associated with GPBC.The ensemble used for our 2015 calculation comprised 864 molecular dynamics (MD) timeunits (after thermalization), upon which 216 measurements were performed separated by 4 MDtime units. Subsequent to the calculation, it was discovered [22] that an error existed in the gen-eration of the random numbers used to set the conjugate momentum at the start of each trajectory,which gave rise to small correlations between widely separated lattice sites. While the resultingeffects were determined to be two-to-three orders of magnitude smaller than our statistical errors,we nevertheless do not include these configurations in the present calculation.In the period following our previous publication, we have dramatically increased the number ofmeasurements. Configurations were generated on seven independent Markov chains originatingfrom widely seperated configurations of our original ensemble. Subsequent algorithmic improve-ments, particularly the introduction of the exact one-flavor algorithm (EOFA) [23–25] further en-hanced our rate of generation such that we have completed over 5000 additional MD time units todate.Continuing with a measurement separation of 4 MD time units, we can potentially performalmost 1300 measurements in total. For this analysis, we include measurements on ∼
60% of the0available configurations, totaling 741. We aim to provide updated results containing measure-ments on the remaining portion in a future publication. For further information on the ensembleproperties, generation algorithms and details of the configurations used for this analysis we referthe reader to Ref. [18].
B. Goodness-of-fit and error estimation
Aside from the central values of our fit parameters we must also estimate the standard error andthe goodness-of-fit. These are obtained via bootstrap resampling, specifically the non-overlappingblock bootstrap variant [26] which allows us to account for mild autocorrelation effects observedin our data. A block size of 8 is used.The bootstrap measurement of the goodness-of-fit is a technique developed specifically for thisand our companion work [18], and is detailed in Ref. [27]. To summarize, the goodness-of-fit istypically parameterized by a p-value that represents the likelihood that the data agrees with themodel, allowing only for statistical fluctuations. The p-value is computed by first measuring q = ∑ i , j ( ¯ x i − f ( (cid:126) p , i )) ( cov ) − i j (cid:0) ¯ x j − f ( (cid:126) p , j ) (cid:1) . (3)where ¯ x i are the ensemble-means of the data at coordinate i , (cid:126) p the fitting parameters, f the modelfunction, and ( cov ) ab is the covariance matrix. The value obtained for q is then compared to thenull distribution that describes how this quantity varies between independent experiments if onlystatistical fluctuations are allowed around the model. The null distribution is typically assumedto be the χ distribution, but this is inappropriate when the fluctuations in the covariance matrixbetween experiments become significant, as is the case for our ππ measurements [27]. In that workwe demonstrate that the null distribution can be estimated directly from the data through a simplebootstrap procedure, allowing for a more reliable p-value that is free from assumptions. Thisprocedure also has the benefit of allowing us to neglect the autocorrelations in the determinationof the covariance matrix on each bootstrap ensemble, which dramatically improves the statisticalerror but changes the definition of q in a subtle way that cannot be accounted for by traditionalmethods.1 C. Measurement technique
Measurements are performed using the all-to-all (A2A) propagator technique of Ref. [28],whereby the quark propagator is decomposed into an exact low-mode contribution obtained froma set of, in our case 900, predetermined eigenvectors, and a stochastic approximation to the high-mode contribution. This allows for the maximal translation of correlation functions in order to takefull advantage of each configuration, as well as easy implemention of arbitrarily smeared sourceand sink operators. We perform full spin, color, flavor and time dilution such that the stochasticsource is required only to produce a delta-function in the spatial location.For all quantities we use smeared meson sources with an exponential (1 s hydrogen wavefunction-like) structure, Θ ( | (cid:126) x − (cid:126) y | ) = exp ( −| (cid:126) x − (cid:126) y | / α ) (4)where α = (cid:126) x and (cid:126) y are the spatial coordinates of the two quarkoperators. Several technicalities must be considered when using G-parity boundary conditions,including limitations on the allowed quark momenta which has implications for the cubic rotationalsymmetry, the preservation of which is essential for producing an operator that projects onto therotationally-symmetric (s-wave) ππ state. These are detailed in Ref. [18].More specific details of the various measurements are provided in the following sections. III. RESULTS FROM TWO-POINT CORRELATION FUNCTIONS
In order to compute the K → ππ matrix elements it is necessary to measure the energies andamplitudes of the pion, kaon and ππ two-point Green’s functions. In this section we present resultsfor the kaon two-point function and summarize the results of Ref. [18] for the pion and ππ two-point functions. We also detail the determination of the energy dependence of the phase shift atthe kaon mass scale, which is used to obtain the Lellouch-L¨uscher [12] finite volume correction tothe matrix elements. A. Notation
G-parity boundary conditions mix quark flavor at the boundary, introducing additional Wickcontractions in which a quark propagates through the boundary and is annihilated by an operator2
State Fit Range
A E p-valueKaon 10-29 4 . ( ) × . ( ) × TABLE II: Fit results in lattice units, fit ranges and p-values for the pion and kaon states. Here E is the energy of the state in question, which for the kaon is equal to the kaon mass, m K .of the opposite quark flavor. In Ref. [10] we introduced a notation whereby the quark field and itsG-parity partner are placed in a two-component vector, (5) ψ l = dC ¯ u T , where C is the charge conjugation matrix. We will refer to the index of these vectors as a “flavorindex”. In this notation the propagator becomes a 2 × s (cid:48) , into which the strange quark transforms at theboundary. The corresponding field operator is (6) ψ h = sC ¯ s (cid:48) T . With the introduction of this extra quark flavor a square-root of the s / s (cid:48) determinant is required inorder to generate a 2+1 flavor ensemble [10]. B. Kaon two-point function
Following Ref. [10], a stationary (G-parity even) kaon-like state can be constructed as (7) | ˜ K (cid:105) = √ (cid:0) | K (cid:105) + | K (cid:48) (cid:105) (cid:1) , where K is the physical kaon and K (cid:48) a degenerate partner with quark content ¯ s (cid:48) u . This | ˜ K (cid:105) statecan be created using the following operator (8) O ˜ K ( t ) = i √ ∑ (cid:126) x ,(cid:126) y e i (cid:126) p · ( (cid:126) x − (cid:126) y ) ψ l ( (cid:126) x , t ) γ Θ ( | (cid:126) x − (cid:126) y | ) ( + σ ) ψ h ( (cid:126) y , t ) (cid:126) p = ( , , ) π L is the quark momentum and Θ is defined in Eq. (4). Note that in the aboveequation and for the other operators presented in this document, the projection operators ( ± σ ) appear; these are necessary to define quark field operators that are eigenstates of translation andhence have definite momentum [10].The two-point function (9) C K ( t , t ) = (cid:104) | O †˜ K ( t ) O ˜ K ( t ) | (cid:105) is measured for all t and t , and subsequently averaged over t at fixed t = t − t . The data arefolded in t , i.e. data with t = t − t are averaged with those with t = L T − ( t − t ) , where L T is thelattice temporal extent, to improve statistics. We perform correlated fits to the following function,(10) C K ( t ) = A K (cid:16) e − m K t + e − m K ( L T − t ) (cid:17) , where the second term accounts for the state propagating backwards in time through the latticetemporal boundary. The chosen fit range, p-value and the results of the fit are given in Tab. II. Inphysical units our kaon mass is 490.5(2.4) MeV, which is within 2% of the physical neutral kaonmass. C. Pion two-point function
The isospin triplet of pion states can be constructed from the operators listed in Sec. V.A. ofRef. [10]. Due to the isospin symmetry the resulting two-point functions all have the same Wickcontractions, and are most conveniently generated with the neutral pion operator, (11) O π ( (cid:126) p π , t ) = − i √ ∑ (cid:126) x ,(cid:126) y e i ( (cid:126) p · (cid:126) x + (cid:126) p · (cid:126) y ) ψ l ( (cid:126) x , t ) γ σ Θ ( | (cid:126) x − (cid:126) y | ) P (cid:126) p π ψ l ( (cid:126) y , t ) , where (cid:126) p π = (cid:126) p + (cid:126) p is the total pion momentum and P (cid:126) p π is a flavor projection operator of theform ( ± σ ) whose sign depends on the particular choice of the quark momentum, per thediscussion in Sec. IV.G. of Ref. [10]. We measure the two-point function with four differentmomentum orientations related by cubic transformations in order to improve the statistical error: (cid:126) p π ∈ { ( , , ) , ( − , , ) , ( , − , ) , ( , , − ) } in units of π / L . The corresponding choices ofquark momentum are given in Ref. [18]. The two-point function (12) C π ( (cid:126) p π ; t , t ) = (cid:104) | O † π ( (cid:126) p π , t ) O π ( (cid:126) p π , t ) | (cid:105) is again averaged over all source timeslices and also over all four momentum orientations, and thedata folded to improve statistics. Correlated fits are performed to the function, (13) C π ( t ) = A π (cid:16) e − E π t + e − E π ( L T − t ) (cid:17) , Parameter Value2-state fit 3-state fitFit range 6-15 4-15 A ππ ( ) . ( ) . ( ) A ππ ( ) . ( ) . ( ) A σ − . ( ) − . ( ) E . ( ) . ( ) A ππ ( ) . ( ) . ( ) A ππ ( ) − . ( ) − . ( ) A σ . ( ) . ( ) E . ( ) . ( ) A ππ ( ) — 0 . ( ) A ππ ( ) — 0 . ( ) A σ — 0 . ( ) E — 0 . ( ) p-value 0.314 0.092 TABLE III: Fit parameters in lattice units and the p-values for multi-operator fits to the I = ππ two-point functions. Here E i are the energies of the states and A i α represents the matrix elementof the operator α between the state i and the vacuum, given in units of √ × . The secondcolumn gives the parameters for our primary fit which uses two-states and three operators. Thethird column shows a fit with the same three operators and one additional state that is used toprobe the systematic effects of this third state on the K → ππ matrix element fits.where t = t − t as before. The chosen fit range, p-value and the results of the fit are also givenin Tab. II. In physical units, and assuming the continuum dispersion relation, our pion mass is142.3(8) MeV, approximately 5% larger than the physical value of 135 MeV. The small effect ofthis difference on our final results is expected to be negligible in comparison to our other errors.5 D. I = I = I = ππππππ two-point function Details of the strategy for measuring the I = ππ two-point function can be found in Ref. [18].In summary, we construct three operators with the quantum numbers of the I = ππ state: Thefirst and second operators, labeled ππ ( ) and ππ ( ) , comprise two single-pion operators car-rying equal and opposite momenta separated by ∆ = ( ± , ± , ± ) π / L , andthose of the latter in the set ( ± , ± , ± ) π / L and permutations thereof. We average over all non-equivalent directions of the pion momentum in order to project onto the rotationally symmetricstate. The final, σ operator corresponds to the scalar two-quark operator √ ( ¯ uu + ¯ dd ) . As men-tioned previously, the pion and σ bilinear operators are smeared with a hydrogen wavefunction(exponential) smearing function of radius 2 lattice sites in order to improve their overlap with thelowest-energy states.Two-point correlation functions are constructed from pairs of source and sink operators thus,(14) C ππαβ ( t , t ) = (cid:104) | O † α ( t ) O β ( t ) | (cid:105) − (cid:104) | O † α ( t ) | (cid:105)(cid:104) | O β ( t ) | (cid:105) , where we include an explicit vacuum subtraction. Here t specifies the earliest time in which anyfermion operator appearing in the annihilation operator O † α is evaluated, and likewise t is thelatest time appearing in the creation operator O β , such that t = t − t is the time of propagation ofthe shortest-lived pion state. We average over many t at fixed t = t − t and the data are foldedto improve statistics as follows: (15) C ππαβ ( t ) → (cid:104) C ππαβ ( t ) + C ππαβ ( L T − ∆ i − ∆ j − t ) (cid:105) where ∆ i = ππ ( ) and ππ ( ) operators and zero for the σ operator. To the matrixof correlation functions we perform simultaneous correlated fits to the functions, (16) C ππαβ ( t ) = i max ∑ i = A i α A i β (cid:16) e − E i t + e − E i ( L T − ∆ i − ∆ j − t ) (cid:17) . We will use the result obtained by uniformly fitting to the temporal range 6 −
15 with all three ππ source/sink operators and allowing for two intermediate states ( i max = ππ and kaon energies differ by 2.2(3)%,where the error is statistical only, and as such our K → ππ calculation is not precisely energy-conserving. The effect of this difference is incorporated as a systematic error on our final result,as discussed in Sec. VII B.6It is interesting to compare the statistical errors of our ππ ground-state fit parameters to those ofour 2015 analysis, which was performed using a single operator ( ππ ( ) ) and the same t min = A ππ ( ) = . ( ) E = . ( ) . Comparing these to the results of this work in Tab. III we find that the error on the ground-stateamplitude has reduced by a factor of 1 . (cid:112) / = . σ operator, which vastly enhance theresolution on the ground-state energy.The I = ππ scattering phase shift is obtained via L¨uscher’s method [11, 29] and has the value,(18) δ = . ( . )( . ) ◦ , where the errors are statistical and systematic, respectively. The procedures by which we estimateour errors are detailed in Ref. [18].Our decision to fit the ππ two-point function with two states limits the number of states that wecan include in our K → ππ matrix element analysis. In order to study the possibility of residualcontamination from a third state we repeat the analysis of the ππ two-point function with 3 states,the results of which are given in the third column of Tab. III. For a stable fit to the ππ data wefound it necessary to use t min =
4, which is lower than the t min = K → ππ matrix elements aresmall. E. Phase-shift derivative at the kaon mass
As detailed in Sec. VI A, the Lellouch-L¨uscher finite volume correction to the K → ππ matrixelements requires the evaluation of the derivative of the phase-shift with respect to the ππ energyevaluated at the kaon mass scale, or more specifically with respect to the variable q = kL / π where k = ( E ππ / ) − m π is the square of the interacting pion momentum. This derivative cannot7presently be obtained experimentally at this energy scale, and therefore an interpolating ansatz ordirect lattice measurement is required.In Ref. [18], alongside the stationary state examined above, we also compute the ππ energy atseveral non-zero center-of-mass momenta, allowing us to obtain the phase-shift at two values ofthe rest-frame energy that are lower than the kaon mass as well as a threshold determination of thescattering length. These results are also close to their corresponding dispersive predictions, albeitwith somewhat larger excited-state systematic errors. Using these results we can directly measurethe derivative of the phase-shift with respect to the energy using a finite-difference approximation,for which we obtain (19)d δ d q = . ( ) radfrom the difference with the nearest energy to the kaon mass, and (20)d δ d q = . ( ) radfrom the next-to-nearest.We can also obtain the derivative from the dispersive prediction of Colangelo et al [16]. Thederivative with respect to s = E ππ , computed at our lattice ππ energy using Eqs. 17.1-17.3 ofRef. [16] with m π =
135 MeV, is found to be (21)d δ d s = . ( ) × − rad MeV − , where the error is the statistical error arising from the uncertainty in the lattice spacing and mea-sured lattice ππ energy. Note that this result is obtained at the physical pion mass, which is 5%smaller than our lattice value. In order to estimate the impact of the difference in pion masses onthis derivative we use NLO chiral perturbation theory [16, 30] (ChPT) to estimate the derivativewith respect to energy at k = . √ s is roughly constant (which is well motivated by the dispersion theoryresult, cf. Fig. 7 of Ref. [16]) we estimate the change in d δ d s evaluated at our lattice ππ energyas 1.2%. This value is small relative to the final systematic error we assign to the derivative inSec. VII D and can therefore be neglected here. Finally, applying d s / d q = . ( ) × MeV ,where again the errors are statistical, we obtain (22)d δ d q = . ( ) rad . δ d E ππ ≈ δ E ππ − m π may also be appropriate. With this ansatz we find (24)d δ d q = . ( ) rad . Given that the derivative of the phase shift is a subleading contribution and that the above valuesare all in reasonable agreement, we expect that the Lellouch-L¨uscher factor can be obtained reli-ably. The variation in these results will be taken into account in our systematic error in Sec. VII D.In our 2015 work [1] we also considered a linear ansatz in q , (25)d δ d q ≈ δ q for which we obtain (26)d δ d q = . ( ) rad . This value is not as well motivated as the ansatz in Eq. (24) and is in disagreement with all fourof the above results. Given the good agreement between our measured phase-shifts and the aboveestimates of the derivative with the dispersive predictions, we will not include this result in oursystematic error estimate.
F. Optimal ππππππ operator
For use later in this document we define here an optimal operator that maximally projects ontothe ππ ground state relative to the first-excited state.Under the excellent assumption that the backwards-propagating component of the time depen-dence is small in the fit window, the two-point functions can be described as a sum of exponentials:(27) C ππαβ ( t ) = ∑ i A i α A i β e − E i t , where again Greek indices denote operators and Roman indices states. We wish to define anoptimized operator that projects onto the ground state: (28) O opt = ∑ α O α r α , for which (29) C ππ opt ( t ) = (cid:104) | O †opt ( t ) O opt ( ) | (cid:105)≈ [ A ] e − E t , ππ operators. Expanding the Green’s function, (30) (cid:104) | O †opt ( t ) O opt ( ) | (cid:105) = ∑ αβ r α (cid:104) | O † α ( t ) O β ( ) | (cid:105) r β = ∑ i ∑ αβ r α A i α A i β r β e − E i t = ∑ i (cid:20) ∑ α A i α r α (cid:21) e − E i t . Without loss of generality we can fix A =
1, which alongside Eq. (29) is sufficient to define r i :(31) ∑ α A i α r α = δ i , . If the number of states is equal to the number of operators this can be interpreted as a matrixequation, (32) A (cid:126) r = ˆ0 , where the row index of A is the state index i and the column index the operator index α . Here ˆ0 isa unit vector in the 0-direction, and as such (33) (cid:126) r = A − ˆ0 . which gives (34) r α = [ A − ] α , i.e. (cid:126) r is the first column of the inverse matrix.As our ππ fits include only two states, we drop the noisier ππ ( ) operator in order to form asquare matrix of correlation functions. We then obtain (cid:126) r T = ( . ( ) × − , − . ( ) × − ) (35)where the elements are the coefficients of the ππ ( ) and σ operators, respectively. In Fig. 1we compare the effective energy obtained with the optimal operator to that of the ππ ( ) and σ operators alone. We observe a marked reduction in the ground-state energy and a noticeableimprovement in the length of the plateau region resulting from the removal of excited-state con-tamination, as well as a significant improvement in the statistical error. This optimal operator willalso be used in our matrix element fits in the following section.0 t E e ff (111),(111) FIG. 1: A comparison of the effective ground-state energy obtained from the optimal operator( i.e. the optimal combination of the σ and ππ ( ) operators, and labeled “ ππ ( ) , σ ” here)with the energies obtained from the σ and ππ ( ) operators separately. (a) type1 (b) type2 (c) type3 (d) type4 FIG. 2: The four classes of K → ππ Wick contractions.
IV. RESULTS FROM THREE-POINT CORRELATION FUNCTIONS FOR ∆ I = / , K → ππ ∆ I = / , K → ππ ∆ I = / , K → ππ DECAYS
In this section we detail the measurement and fitting of the K → ππ three-point Green’s func-tions, from which the unrenormalized matrix elements (cid:104) ( ππ ) I = | Q i | K (cid:105) are obtained.1 A. Overview of measurements
On the lattice we measure the following three-point functions, (36) C i ( t , t K → snksep ) = (cid:104) | O †snk ( t K → snksep ) Q i ( t ) O ˜ K ( ) | (cid:105) , where t denotes the time separation between the kaon and four-quark operators, and t K → snksep thetime separation between the kaon and the ππ “sink” operator, O snk . As described in Ref. [31], theWick contractions of these functions fall into four categories based on their topology, as illustratedin Fig. 2.Note that here and below we take care to differentiate between the G-parity kaon state ˜ K ,which is a G-parity even eigenstate of the finite-volume Hamiltonian, and the physical kaon K that is not an eigenstate of the system. The matrix elements of the physical kaon are related tothose of the G-parity kaon by a constant multiplicative factor of √ type3 and type4 diagrams these are measured with kaon sources on every timeslice, 0 ≤ t K < L T . The farmore precise type1 and type2 contributions are measured every eighth timeslice in order to reducethe computational cost. For the remainder of this section we will assume all correlation functionsto have been averaged over the kaon timeslice where appropriate.We compute each diagram with 5 different time separations between the kaon and the ππ sink operators, t K → snksep ∈ { , , , , } , with the ∆ S = ππ operator for those cases when the ππ operator is aproduct of single-pion operators evaluated on different time slices. (This convention of specifyingthe minimum time separation from those ππ operators which are non-local in the time is followedthroughout this paper.) As these ππ operators comprise back-to-back moving pions with zero totalmomentum, we must measure each diagram for all possible orientations of the pion momenta inorder to project onto the rotationally symmetric state.The type3 and type4 diagrams both contain a light or strange quark loop beginning and endingat the operator insertion point that results in a quadratic divergence regulated by the lattice cutoff.This divergence is removed by defining the subtracted operators [31, 32], (37) Q i → Q i − α i ¯ s γ d . Q i . The coefficients α i inEq. (37) are defined by imposing the condition, (38) (cid:104) | (cid:110) ˆ Q i ( t ) − α i ( t )[ ¯ s γ d ]( t ) (cid:111) O ˜ K ( ) | (cid:105) = , where we have allowed α i to vary with time as this was found to offer a minor statistical improve-ment. Although the matrix element of this pseudoscalar operator vanishes by the equations ofmotion for energy-conserving kinematics and is therefore not absolutely necessary for our calcu-lation, the subtraction reduces the systematic error resulting from the small difference between our ππ and kaon energies while simultaneously reducing the statistical error and suppressing excited-state contamination.Due to having vacuum quantum numbers, the I = ππ operators project also onto the vacuumstate and this off-shell matrix element dominates the signal unless an explicit vacuum subtractionis performed, (39) C i ( t , t K → snksep ) → C i ( t , t K → snksep ) − (cid:104) | O †snk ( t K → snksep ) | (cid:105)(cid:104) | Q i ( t ) O ˜ K ( ) | (cid:105) . However, due to our definition of the subtraction coefficient α i in Eq. (38), the vacuum matrixelements appearing in the right-hand side vanish making this subtraction unnecessary. In practicethis cancellation is not exact in our numerical analysis for the following reason: While the ππ “bubble” (cid:104) | O †snk | (cid:105) is formally time-translationally invariant we observed a minor statistical ad-vantage in evaluating this quantity with the ππ operator on the same timeslice as it appears in thefull disconnected Green’s function that is being subtracted, such that it is maximally correlated.Therefore, for the right-most term in Eq. (39) we compute (40)1 n t K ∑ t K ∈{ t K } (cid:104) | O †snk ( t K + t K → snksep ) | (cid:105)(cid:104) | (cid:110) ˆ Q i ( t + t K ) − α i ( t )[ ¯ s γ d ]( t + t K ) (cid:111) O ˜ K ( t K ) | (cid:105) , where t K is the kaon timeslice and { t K } the set of timeslices upon which measurements wereperformed, i.e. with the product of the K → vacuum matrix element and the ππ bubble performedunder the average over the kaon source timeslice rather than after. As suggested by the above,the coefficients α i ( t ) are computed separately from the t K -averaged matrix elements and thereforethe cancellation between the two terms in brackets is exact only up to the degree to which thetime translation symmetry is realized at finite statistics. Due to our large statistics we found thedifference in the fitted Q matrix element obtained with and without the vacuum subtraction to beat the 0.1% level.3 t -3-2-101234 C ( t , t K s n k s e p = )( × ) type type type type t -8-6-4-2024 C ( t , t K s n k s e p = )( × ) type type type type FIG. 3: The contributions of the four Wick contraction topologies type1 - type4 to the C (left) and C (right) three-point functions with the ππ ( ) sink operator, plotted as a function of the timeseparation between the kaon and the four-quark operator, t , at fixed t K → snksep =
16. For clarity weplot with an inverted x-axis such that the ππ sink operator is on the left-hand side. Thesecorrelation functions include the subtraction of the pseudoscalar operator.We perform measurements with all three two-pion operators described in Sec. III D. For the K → ππ matrix elements of the four-quark operators, the full set of Wick contractions for the ππ ( ) and ππ ( ) sink operators can be found in Appendix B.1 and B.2 of Ref. [33], and thoseof the σ operator in Appendix A of this document. The Wick contractions for the K → ππ matrixelements of the pseudoscalar operator (with all three sink operators) as well as the K → vacuummatrix elements of this and the four-quark operators are provided in Appendix B of this document.In Fig. 3 we plot the contributions of the four classes of Wick contraction illustrated in Fig. 2 tothe three-point functions of the (subtracted) Q and Q operators with the ππ ( ) sink operator.As the individual topologies are not separately interpretable as Green’s functions of the QCDpath integral, their time dependence is not necessarily described by the propagation of physicaleigenstates of the QCD Hamiltonian. As such we cannot combine our data sets with different t K → snksep when generating such plots, and instead plot with a single, fixed t K → snksep =
16. Despite theinability to interpret the time dependence physically, we can look at the relative contributions ofeach topology within the central region of the plot in which the behavior of the combined data isdominated by the kaon and ππ ground-states, i.e. the region in which we perform our fits below.Our final choices of cut incorporate data from this set in the range 6 ≤ t ≤
11 (cf. Sec. IV E 4). Inthis window we observe that for both the C and C correlation functions, the contribution of thenoisy, type4 disconnected diagrams is largely consistent with zero, albeit with much larger errors4for the former. C appears dominated by the type1 and type3 diagrams, which both contributewith the same sign, with a negligible contribution from the type2 diagrams. The contribution ofthe type1 and type3 diagrams appears to behave similarly for the C three-point function, howeverhere we observe a strong cancellation between those and the type2 diagrams. B. Determination of α i α i α i The subtraction coefficients α i are computed via Eq. (38) as the following ratio of two-pointfunctions, (41) α i ( t ) = (cid:104) | ˆ Q i ( t ) O ˜ K ( ) | (cid:105)(cid:104) | [ ¯ s γ d ]( t ) O ˜ K ( ) | (cid:105) , where the average of the correlation functions over the kaon source timeslice is implicit as above.The Wick contractions for the (cid:104) | ˆ Q i ( t ) O ˜ K ( ) | (cid:105) two-point functions are identical to the com-ponents of the type4 K → ππ diagrams that are connected to the kaon. While these connectedcomponents are formally independent of the sink two-pion operator, in practice these quantitieswere computed using code that was organized differently for the ππ and σ operators. As describedin Appendix B of this paper and Appendix B.2 of Ref. [33], the factors entering the type4 diagramsthat determine the α i were constructed from two separate bases of functions of the quark propa-gators, one for the σ and the other for the ππ ( . . . ) operators, where for each basis γ hermiticitywas used in a different way. While γ hermiticity is an exact relation, the fact that we are using astochastic approximation for the high modes of the all-to-all propagator allows small differencesto arise between the values of the α i computed in these two bases. We therefore have separateresults for the α i from the ππ and σ three-point functions calculations.In Fig. 4 we plot the time dependence of the α i for all ten operators. We observe excellentagreement between the results obtained from the two different bases of contractions as expected.For a number of operators we find statistically significant but relatively small excited-state con-tamination for small t that in all cases appears to die away by t =
6. While the effects of thiscontamination are unlikely to significantly affect our final results, the cuts that we later apply toour fits nevertheless exclude data with t < t ( t ) (111) t ( t ) (111) t ( t ) (111)2 4 6 8 10 12 14 16 t ( t ) (111) t ( t ) (111) 2 4 6 8 10 12 14 16 t ( t ) (111) t ( t ) (111) t ( t ) (111) t ( t ) (111) t ( t ) (111) FIG. 4: The pseudoscalar subtraction coefficient α i as a function of time for each of the tenoperators in the following order: ˆ Q - ˆ Q on the first line, ˆ Q - ˆ Q on the second, ˆ Q - ˆ Q on the thirdand ˆ Q on the fourth. Red circles denote data obtained in the basis of correlation functions usedfor the ππ ( ) operator, and blue squares for the σ sink operator. C. (cid:104) ππππππ | ¯ s γ d ¯ s γ d ¯ s γ d | ˜ K ˜ K ˜ K (cid:105) matrix elements The K → ππ matrix elements of the pseudoscalar operator ¯ s γ d are required to perform thesubtraction of the divergent loop contribution. In this section we independently analyze thesematrix elements in order to understand their time dependence and the corresponding effect of the6subtraction on the amount of excited state contamination in the final K → ππ result.In the limit of large time separation between the source/sink operators and the four-quark opera-tor, only the lowest-energy ππ and kaon states are present. Since the pseudoscalar matrix elementsvanish by the equations of motion when the decay conserves energy and the kaon and ππ ground-state energies in our calculation differ by only 2%, we expect the subtraction to result in only anegligible shift in the central value but a marked improvement in the statistical errors in this limit.However at finite time separations, the contributions of the excited states may take a long time todie away due to the increasing magnitude of the corresponding matrix elements between initialand final states of different energies. It is this concern that prompts us to study this system morecarefully.The lattice three-point function (42) C P ( t , t K → snksep ) = (cid:104) | O †snk ( t K → snksep ) [ ¯ s γ d ]( t ) O ˜ K ( ) | (cid:105) for a generic sink ππ operator, O snk , has the following time dependence: (43) C P ( t , t K → snksep ) = ∑ i j A i in A j out M i jP exp (cid:0) − E i in t (cid:1) exp (cid:16) − E j out ( t K → snksep − t ) (cid:17) , where the subscript ‘in’ refers to the incoming kaonic state, ‘out’ to the outgoing two-pion state,and M i jP is the matrix element for the term involving in and out states i and j , respectively. It is con-venient to define an “effective matrix element” by dividing out the ground-state time dependenceand operator amplitudes, (44) M eff , snk P ( t (cid:48) , t K → snksep ) = M P + ∑ i , j (cid:54) = A (cid:48) i in A (cid:48) j out M i jP exp (cid:16) − ∆ E i in ( t K → snksep − t (cid:48) ) (cid:17) exp (cid:16) − ∆ E j out t (cid:48) (cid:17) , where (45) t (cid:48) = t K → snksep − t is the separation between the four-quark operator and the sink and (46a) A (cid:48) i in / out = A i in / out / A / out , (46b) ∆ E i in / out = E i in / out − E / out . Note that M eff , snk P is dependent on the sink operator through the terms involving the excited states,in which a ratio of ground and excited state amplitudes appears.We measure the correlation function Eq. (42) for each of our three two-pion operators. Note thata vacuum subtraction is also required here and is performed in the same way as for the four-quark7 t M e ff , s n k P ( t , t K s n k s e p ) t M e ff , s n k P ( t , t K s n k s e p ) FIG. 5: The effective pseudoscalar matrix element M eff , snk P as a function of the time separationbetween the four-quark operator and the sink, t (cid:48) . In the left pane we show the data for the ππ ( ) operator (circles) and the σ operator (squares) separately, and in the right pane we showthe same for the optimal operator. Colored data correspond to the different t K → snksep as follows: red(10), green (12), blue (14), orange (16) and mauve (18). The data for each of these differentseparations are staggered in order such that t K → snksep =
10 is the left-most point of each cluster and t K → snksep =
18 the right-most.operators. In Fig. 5 we plot M eff , snk P for the ππ ( ) and σ operators for each of the five valuesof t K → snksep . The corresponding data for the ππ ( ) operator is much noisier and has thereforebeen excluded. The form of this plot can be explained as follows: As E ≈ E we expect M P to be small. If we then assume that the dominant excited state contributions come from the terminvolving the excited kaon state and ground ππ state ( i = , j =
0) and the term with the groundkaon state and the first excited ππ state ( i = , j = M eff , snk P ( t (cid:48) , t K → snksep ) ≈ A (cid:48) M P exp (cid:16) − ∆ E t K → snksep (cid:17) exp (cid:0) + ∆ E t (cid:48) (cid:1) + A (cid:48) M P exp (cid:0) − ∆ E t (cid:48) (cid:1) . This ansatz then implies an exponentially falling contribution from the excited pion state and anexponentially growing piece from the excited kaon state, giving rise to a bowl-like shape assum-ing that A (cid:48) and A (cid:48) have the same sign, which appears to the the case here. Furthermore, theexponentially-growing piece in t (cid:48) is expected to be larger for smaller t K → snksep , and indeed we ob-serve that the turnover point at which the exponentially-growing term begins to dominate occurssooner for smaller t K → snksep .While the effective matrix elements of both sink operators initially trend towards zero, for themore precise ππ ( ) data it seems that none of the five data sets are statistically consistent with8zero at their maxima, suggesting we do not reach the limit of ground-state dominance. This isnot necessarily an issue for our calculation given that the subtraction will heavily suppress thesecontributions in our final result, and furthermore the inclusion of multiple sink operators willimprove our ability to extract the ππ ground-state matrix element. In order to disentangle thesetwo effects it is convenient to examine the three-point function for the optimized sink operatordiscussed in Sec. III F. The time dependence of M eff , snk P for this operator is also shown in Fig. 5.By definition this operator heavily suppresses A (cid:48) j out for j >
0, and indeed we find the data to bemuch flatter in the low- t (cid:48) region and also considerably closer to zero. The exponential growth and t K → snksep dependence that enters due to the excited kaon term is expected to be largely unaffected bythis transformation, however it seems that in several cases the plateaus extend much further intothe large- t (cid:48) region than previously. It is likely that is due to an accidental cancellation owing to thefact that A (cid:48) is positive for the ππ ( ) operator and negative for the σ operator (cf. Tab. III) andhence the exponentially-growing terms for these operators have opposite signs.We conclude by discussing the expected size of the excited-state contamination in the matrixelements of the subtracted four-quark operators arising from the pseudoscalar operator. In the K → ππ calculation, this dimension-3 operator is introduced to remove what in the continuumlimit would be a quadratic divergence resulting from the self-contraction between two of the fourquark operators appearing in those operators ˆ Q i with a component transforming in the ( , ) or ( , ) representations of SU ( ) L × SU ( ) R . In our lattice calculation these terms behave as 1 / a when expressed in physical units. To leading order in a this 1 / a coefficient does not depend onthe external states and is therefore removed from our (cid:104) | ( ππ ) ˆ Q i ˜ K | (cid:105) amplitude by the subtractiondefined above, even though the coefficients α i are determined from the (cid:104) | ˆ Q i ˜ K | (cid:105) matrix elementin Eq. (41). Because of the chiral structure of the ( , ) and ( , ) operators, these coefficientshave the structure: α i ∼ m s − m d a + . . . [8], where the ellipsis represents terms which are not power-divergent.Thus, the ¯ s γ d subtraction removes the leading 1 / a term in the matrix element of ˆ Q i , leavingbehind a finite piece of size ∼ ( m s − m d ) Λ ¯ s γ d . This remainder is not physical and dependson the condition chosen to define the α i . However, it will contribute to our final result if E ππ (cid:54) = m K . For the ground-state component ( i = j =
0) this term is thus heavily suppressed by thefactor ( E ππ − m K ) . However for the excited states we expect this piece to be on the order of thephysical contribution from the dimension-6 four-quark operator. As such it may result in a modestenhancement of the excited state matrix elements. Providing we are able to demonstrate that we9have the excited ππ and kaon states under control through appropriate cuts on our fitting ranges,this should pose no obstacle to our calculation. D. Description of fitting strategy
For a lattice of sufficiently large time extent that around-the-world terms in which states prop-agate through the lattice temporal boundary can be neglected, and assuming that the four-quarkoperator is sufficiently separated from the kaon source that the kaon ground state is dominant,the three-point Green’s functions C i of the weak effective operators defined in Eq. (36) have thegeneral form, C i ( t , t K → snksep ) = ∑ j √ A K A j snk e − m K t M ji e − E j ( t K → snksep − t ) , (48)where M ji = (cid:104) ( ππ ) j | Q i | K (cid:105) is the matrix element of the four-quark operator Q i with the ππ state j , with M i corresponding to the physical K → ππ matrix elements required to compute A . Thefactor of 1 / √ A K is the amplitude of the G-parity kaon operator, A j snk are the amplitudesof the sink operator with the state j , and E j is the energy of that state. These parameters are fixedto those obtained from the two-point function fits in Sec. III: A K and m K to the results given inTab. II, and A j snk and E j to the results obtained from the three-operator, two-state ππ fits given inthe second column of Tab. III.We perform simultaneous correlated fits over multiple sink operators to the form Eq. (48) in or-der to determine the matrix elements M ji , allowing for one or more states j . Independent one-statefits are also performed to the optimized sink operator defined in Sec. III F. The fits are performedto each weak effective operator separately, in the 10-operator basis (the relationship between these10 linearly-dependent operators serves as a useful cross-check of the fit results) using the strategyoutlined in Sec. II B. We apply a cut t min on the separation t between the kaon and the four-quark operator in order to isolate the ground-state kaon, and also a cut t (cid:48) min on the separation t (cid:48) = t K → snksep − t between the four-quark and sink operators. These cuts, the number of sink opera-tors, and the number of excited ππ states included in the fit are varied in order to study systematiceffects.For use below we again define an “effective matrix element” in which the ground-state ππ and0kaon amplitudes and time dependence are multiplied out, (49) M eff , snk i ( t (cid:48) ) = C i ( t , t K → snksep ) (cid:18) √ A K A e − m K t e − E ( t K → snksep − t ) (cid:19) − = M i + ∑ j A j snk A M ji e − ( E j − E ) t (cid:48) . These effective matrix elements converge exponentially to the ground-state matrix element at large t (cid:48) . Note that, unlike in Sec. IV C, we are assuming that a cut, t (cid:48) min , on the separation betwen thekaon and four-quark operators has been applied that is sufficient to isolate the contribution of thekaon ground state. As a result, these effective matrix elements can be assumed to be independentof t K → snksep and a weighted average of our five datasets of different t K → snksep can be applied to improvethe statistical resolution of the data presented in our plots. E. Fit results
In this section we examine the results of fitting various subsets of our data, with the goal offinding an optimal fit window in which systematic errors arising from both excited ππ and kaonstates are minimized.In Figs. 6 and 7 we plot the fitted ground-state matrix elements M i as a function of t (cid:48) min forvarious choices of t min , the number of sink operators and the number of states. The three-operatorfits are performed using the ππ ( ) , ππ ( ) and σ sink operators; for the two-operator fitswe drop the noisier ππ ( ) data; and for the one-operator fits we further drop the σ data. Theone-operator, one-state fits are equivalent to those performed in our 2015 work, albeit with morestatistics and more reliable ππ energies and amplitudes.The discussion below will be focused on these figures. We will first discuss general featuresaddressing the quality of the data and the reliability of the fits, and will then concentrate on search-ing for evidence of systematic effects (or lack thereof) arising from kaon and ππ excited states.Based on those conclusions we will then present our final fit results.
1. Discussion of data and fit reliability
We will first comment on the fits to the optimal operator, labeled “opt.” in the figures. Thisapproach is outwardly advantageous in that the fits are performed to a single state and the covari-1 t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. FIG. 6: Fit results in lattice units for the K → ππ ground-state matrix elements M − M as afunction of t (cid:48) min , the minimum time separation between the four-quark and sink operators that isincluded in the fit. The results have been shifted along the x-axis for clarity in order of theirappearance in the legend. The legends are given in the format × t min on the time separation between the kaon and the four-quark operators. Here “opt.” is the fit tothe optimal operator and “sys.” is used to estimate the systematic error resulting from a thirdstate.2 t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. t min M t min =63×2 t min =73×2 t min =82×2 t min =61×1 t min =6 opt.1×1 t min =63×3 t min =6 sys. FIG. 7: The extension of Fig. 6 to the ground-state matrix elements M − M .ance matrix is considerably smaller. In Fig. 8 we compare the t (cid:48) dependence of the M eff , snk2 and M eff , snk6 effective matrix elements of this optimal operator to that of the ππ ( ) and σ operatorsalone, where we note a marked improvement in the quality of the plateau. This behavior, which isalso accounted for implicitly in the multi-state fits, demonstrates the power of the multi-operatortechnique for isolating the ground state. In Figs. 6 and 7 we observe that the fit results for theoptimal operator agree very well with the multi-state fit results in all cases. While this approachdoes not appear to offer any statistical advantage, the strong agreement suggests that our complexmulti-state correlated fits are under good control.In Figs. 6 and 7 we observe for several ground-state matrix elements a trend in the fit results upto an extremum at t (cid:48) min =
7, followed by a statistically significant correction at the level of 1-2 σ for the fits with t (cid:48) min =
8. In this and Sec. VII A we present substantial evidence that the systematicerrors resulting from excited kaon and ππ states are minimal, which makes it unlikely that thisrise is associated with excited state contamination. Certainly if it were due to excited ππ stateswe would expect an improvement as more sink operators are added, but there is little evidence of3 t M e ff , s n k opt(111) t M e ff , s n k opt(111) FIG. 8: The effective matrix elements M eff , snk2 (left) and M eff , snk6 (right) for the ππ ( ) and σ sink operators and the two-operator two-state optimal sink operator (labeled “opt.” here), plottedas a function of t (cid:48) . The error-weighted average has been applied to the five different K → sinkseparations subject to a cut of t min = i P-value i P-value1 0.314 6 0.4462 0.737 7 0.8433 0.02 8 0.884 0.123 9 0.5815 0.421 10 0.545
TABLE IV: The p-values assessing how well the data with t (cid:48) ≥ C i correlation functions obtained by fitting to 3 operators and 2 states with t (cid:48) min = t min = t min cut, whereas no significant change is observed. The most likely explanation isa statistical fluctuation in our correlated data set, and indeed in Fig. 8 we see evidence of such afluctation peaking at t (cid:48) = t (cid:48) min =
5, which in all cases lie within the plateau region before this rise, are a gooddescription of the subset of data with t (cid:48) ≥
7, or in other words how likely it is that these data areconsistent with this model allowing only for statistical fluctuations. In Tab. IV we list the p-valuesfor these data using the model obtained by fitting to 3 sink operators and 2 states with t (cid:48) min = t M e ff , s n k (111) FIG. 9: The M eff , snk8 effective matrix element for the ππ ( ) (red circles) and σ (blue squares)sink operators overlaid by curves showing the function M eff , snk8 ( t (cid:48) ) predicted using the parametersobtained by fitting the data with t min = t (cid:48) min =
5. The lighter part of the band is the portionof the curve outside of the fit region. An error-weighted average over t K → snksep has been performedto the data. Recall that the effective matrix elements are defined (Eq. 49) such that the resultconverges to the ground-state matrix element at large t (cid:48) .and t min =
6, computed using the technique discussed in Sec. II B (with no free parameters). Weobserve excellent p-values in all cases bar M , and to a lesser extent M . The lower p-valuesfor these operators are common for all of the multi-operator fits and are likely associated withthe statistical fluctuations described above which are more apparent for these matrix elements (cf.Fig. 6). We expect that such unusual statistical fluctuations will be found when so many differentoperators and fitting ranges are examined. Of most importance in a calculation of Im ( A ) is M ,for which we find that the model obtained with t (cid:48) min = t (cid:48) ≥
7. The p-value is in fact little different from the value p = .
525 obtained by fitting to thesedata directly, suggesting that the models are equally good descriptions despite the tension in theground-state matrix elements.For M and M (and to a lesser extent, M ) we observe a discrepancy between the one-operatorand multi-operator results at the 1-2 σ level that persists even to large t (cid:48) min . Given the very clearplateaus in the multi-state fit results, this disagreement is likely due again to statistical effects inthese correlated data. This is evidenced for example in Fig. 9 in which we overlay the M eff , snk8 effective matrix element for the ππ ( ) and σ sink operators by the multi-operator fit curve.We observe that the fit curve for the ππ ( ) operator is completely compatible with the data5 t M e ff , s n k t M e ff , s n k FIG. 10: The t (cid:48) dependence of the M eff , snk4 (left) and M eff , snk6 (right) effective matrix elementswith the optimal sink operator. Each cluster of points, separated for clarity, shows the data for thefive different K → snk separations: 10,12,14, 16, and 18, from left to right in that order. Thedarker, filled points are those that lie within our cut of t min =
6, and the lighter, hollow points arethose excluded.but favors a value that is consistently within the upper half of the error bar, suggesting that theapparent flatness of the ππ ( ) effective matrix element represents a false plateau, and the factthat the multi-operator method is capable of resolving the behavior is a testament to its power.
2. Excited kaon state effects
We now address excited kaon state effects. Because the data rapidly becomes noisier as wemove the four-quark operator closer to the kaon operator and thus further away from the ππ operator, such effects are not expected to be significant. The simplest test is to vary the cut on thetime separation between the kaon operator and the four-quark operator, t min . The first three pointsfrom the left of each cluster in Figs. 6 and 7 show the result of varying t min between 6 and 8 atfixed t (cid:48) min . As expected we observe no statistically significant dependence on this cut.We can also test for excited kaon effects by examining the data near the kaon operator in moredetail, alongside looking for trends in the five different K → sink separations at fixed t (cid:48) . Theoptimal operator proves convenient for examining this behavior as it neatly combines the twodominant sink operators and should be flat within the fit window. In Fig. 10 we plot the data forthe M eff , snk4 and M eff , snk6 effective matrix elements with a distinction drawn between data includedand excluded by a cut on the kaon to four-quark operator time separation of t min =
6. We find no6apparent evidence of excited kaon state contamination even for data excluded by the cut, nor dowe observe any trends of the data in the K → sink separation.We therefore conclude that excited kaon effects in our results are negligible.
3. Excited ππ state effects The dominant fit systematic error is expected to be due to excited ππ states. Fortunately, giventhat we can change both the number of operators and the number of states alongside varying the fitwindow within a region where our data is most precise, there are a number of tests we can performto probe this source of error.We begin by comparing the multi-operator fits to the one-operator ( ππ ( ) ) fit, the latter be-ing equivalent to the procedure used for our 2015 work. In the majority of cases we see littleevidence of excited state contamination in the one-operator data, as evidenced by its agreementwith the multi-operator fits as well as the strong consistency between the fits as we vary the fitwindow. However for the M and M matrix elements we observe strong evidence of excited-statecontamination in these fits at smaller t (cid:48) min . Fig. 6 clearly demonstrates how these effects are sup-pressed as we add more operators: Initially the one-operator results converge with the 3-operatorresults at t (cid:48) min = t (cid:48) min = t (cid:48) min =
3. This suggests that the 5% excited-statesystematic error on our 2015 result which used t (cid:48) min = ππ ( ) data are considerably noisier than those of theother operators, and the associated ππ energy and amplitudes are less-well known, and as suchthese data contribute relatively little to the fit. Nevertheless we do observe that for the M and M matrix elements, the introduction of the third operator results in values that for low t (cid:48) min (3 or 4)are in considerably better agreement with the results for larger t (cid:48) min , suggesting that in the regimein which these data are less noisy (i.e. closer to the ππ operator) the third operator is acting toremove some residual excited-state contamination. We conclude that it is beneficial to include thethird operator.7 t M e ff , s n k (111) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 t M e ff , s n k (111)0 1 2 3 4 5 6 7 8 9 10 11 12 13 t M e ff , s n k (111) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 t M e ff , s n k (111) FIG. 11: The effective matrix element M eff , snk6 for the ππ ( ) (red circles) and σ (blue squares)sink operators, overlaid by the fit curves. The lighter part of the band is the portion of the curveoutside of the fit region. The upper panels are for the 2-state fits and the lower panels are for the3-state fits. In each case the left panel is for t (cid:48) min = t (cid:48) min =
5. All fits areperformed with 3 operators and use t min = ππ two-point function fit parametersgiven in the third column of Tab. III and the same fit ranges for t and t (cid:48) used in the three-operator,two-state fits. The results for the ground-state matrix elements are also included in Figs. 6 and 7with the label “sys.”. We find that including this third state has very little impact and the resultsagree very well with the three-operator, two-state fits. This again suggests that we have the ππ excited-state systematic error under control.A further test for excited-state contamination is to study the agreement of the fit curves withthe data outside of the fit region. To this end in Fig. 11 we plot the ππ ( ) and σ operatordata for the M eff , snk6 effective matrix element overlaid by the fit curves for the 3-operator, 2-state8fits, and for the 3-operator, 3-state fits described above, using t (cid:48) min = t (cid:48) min = ππ ( ) data at t (cid:48) = σ data at this timeslice. Fitting with t (cid:48) min = σ operator data at t (cid:48) =
4. This is consistent with the pattern of couplings of the operators to the statesin Tab. III which show a significant reduction in the couplings to higher states for the ππ ( ) operator but almost equal-sized couplings of the σ operator to all three states. The 3-operator,3-state fit with t (cid:48) min = ππ two-point datafrom this timeslice. However with t (cid:48) min = E = . ( ) (in lattice units) obtained by our fits issomewhat larger than the value of E = .
692 predicted by dispersion theory suggesting that theeffects of even higher excited states may be playing a role here. Nevertheless the strong agreementbetween the ground-state matrix elements for all of these fits suggest that the residual effects ofthe higher excited states on the 3-operator, 2-state fits are negligible.For our final result we choose to focus upon the three-operator, two-state fits. While the ma-jority of the corresponding curves in Figs. 6 and 7 are essentially flat from t (cid:48) min =
3, we opt fora conservative and uniform cut of t (cid:48) min =
4. Final fit results
As discussed above we choose the 3-operator, 2-state fit with t (cid:48) min = t min =
6. In Tab. V we present the full set of p-values and parameters forthese fits. We obtain acceptable p-values in the majority of cases, with the notable exception of the Q four-quark operator for which p = t (cid:48) min , and also that the p-value of the one-operator, one-state fit with the same fit range – withwhich our chosen value is in excellent agreement – has a p-value of 15%. The low probability9 Param Value Param Value M − . ( ) M . ( ) M . ( ) M − . ( ) p-value 0.488 p-value 0.743 M . ( ) M . ( ) M . ( ) M − . ( ) p-value 0.036 p-value 0.139 M − . ( ) M − . ( ) M . ( ) M . ( ) p-value 0.458 p-value 0.159 M . ( ) M . ( ) M − . ( ) M − . ( ) p-value 0.913 p-value 0.676 M − . ( ) M . ( ) M . ( ) M − . ( ) p-value 0.327 p-value 0.56 TABLE V: Final K → ππ matrix element results in lattice units obtained from a three-operator,two-state fit with t min = t (cid:48) min =
5. Here M ji refers to the matrix element of the Q i operatorwith ππ state j .is therefore unlikely to be associated with any systematic effect and can be attributed to low-probability statistical effects.We conclude this section with a comparison of the statistical errors of the matrix elements M and M to those determined in our 2015 analysis. Previously we obtained (50) M = . ( ) M = − . ( ) . Comparing these values to those in Tab. V we find that the errors have reduced by factors of 2.8 and2.4 for M and M , respectively. Comparing the 3-operator, 2-state fits to the 1-operator, 1-statefits in Fig. 6 we observe that the larger improvement for M can be explained by the additionaloperators, however for M these two approaches have similar errors. The fact that the error on0 M has improved considerably more than the factor of 1.9 expected by the increase in statisticscan therefore be attributed to the improved precision of the ππ two-point function fits observed inSec. III D. V. NON-PERTURBATIVE RENORMALIZATION OF LATTICE MATRIX ELEMENTS
The Wilson coefficients are conventionally computed in the MS (NDR) renormalizationscheme, and therefore we are required to renormalize our lattice matrix elements also in thisscheme. This is achieved by performing an intermediate conversion to a non-perturbative reg-ularization invariant momentum scheme with symmetric kinematics (RI-SMOM). As the namesuggests, these schemes can be treated both non-perturbatively on the lattice (provided the renor-malization scale is sufficiently small compared to the Nyquist frequency π / a ) and in continuumperturbation theory (providing the renormalization scale is sufficiently high that perturbation the-ory is approximately valid at the order to which we are working). Thus, we can use continuumperturbation theory to match our RI-SMOM matrix elements to MS, avoiding the need for latticeperturbation theory. The matching factors have been computed to one-loop in Ref. [34].In our 2015 calculation we computed the renormalization matrix at a somewhat low renormal-ization scale of µ = .
529 GeV in order to avoid large cutoff effects on our coarse, a − = .
38 GeVensemble. Due to this low scale, the systematic error associated with the perturbative RI to MSmatching was our dominant error, with an estimated size of 15%. In this paper we utilize the step-scaling procedure [35, 36] (summarized below) in order to circumvent the limit imposed by thelattice cutoff and increase the renormalization scale to 4.0 GeV at which the error arising from theuse of one-loop perturbation theory is expected to be significantly smaller. A separate step-scalingcalculation to 2.29 GeV was performed in Ref. [37] and we will utilize those results to study thescale dependence of the perturbative and discretization errors in our operator normalization.
A. Summary of approach
Due to operator mixing, the renormalization factors take the form of a matrix. This is mostconveniently expressed in the seven-operator chiral basis in which the operators are linearly inde-pendent and transform in specific representations of the SU ( ) L ⊗ SU ( ) R chiral symmetry group,an accurate symmetry of our DWF formulation even at short distances. In this basis the renormal-1ization matrix is block diagonal, with a 1 × Q (cid:48) operator that transformsin the ( , ) representation, a 4 × ( , ) operators Q (cid:48) , Q (cid:48) , Q (cid:48) and Q (cid:48) , and a 2 × ( , ) operators Q (cid:48) and Q (cid:48) .In the RI-SMOM scheme the renormalized operators are generally defined thus, (51) O RI i = Z RI ← lat i j O lat j where Einstein’s summation conventions are implied and the label “RI” is used as short-hand forthe RI-SMOM scheme. The renormalization factors are defined via (52) Z − q [ P m ] β αδ γ [ Γ RI im ] αβ γδ ( p , p ) = F im , where the index m is not summed over. Here α − δ are combined spin and color indices, Z q is thequark field renormalization, q is a four-momentum that defines the renormalization scale and P m are “projection matrices” described below. The quantities F im on the right-hand side are found byevaluating the left-hand side of the equation at tree level. Γ RI im are computed as (53) [ Γ RI im ] αβ γδ ( p , p ) = (cid:28) E m ∑ x e iqx O RI i ( x ) (cid:29) αβ γδ amp . where the sum is performed over the full four-dimensional lattice volume and q = p − p . Here E m are a set of seven four-quark operators that each create the four quark lines that connect to theweak effective operator, E = E = E = E = s ( − p ) ¯ d ( p ) u ( − p ) ¯ u ( p ) E = E = E = s ( − p ) ¯ d ( p ) ∑ q = u , d , s q ( − p ) ¯ q ( p ) , (54)where the momentum arguments indicate the incoming momenta and the quark momenta satisfy symmetric kinematics : p = p = ( p − p ) = q ≡ µ . The subscript “amp.” in Eq. (53) impliesthat the external propagators are amputated by applying the ensemble-averaged inverse propagator,such that the resulting Green’s function has a rank-4 tensor structure in the spin-color indices.These Green’s functions are not gauge-invariant, hence the procedure must be performed usinggauge-fixed configurations, for which we employ Landau gauge-fixing. The use of momentum-space Green’s functions introduces contact terms that prevent the use of the equations of motionso that additional operators, beyond those needed to determine on-shell matrix elements, mustbe introduced if all possible operator mixings are to be included, as is required if the RI-SMOMscheme is to have a continuum limit. These are discussed below.2Note that the Wick contractions of Eq. (53) result in disconnected penguin-like diagrams thatinteract only by gluon exchange; these diagrams are evaluated using stochastic all-to-all propaga-tors and are typically noisy, requiring multiple random hits and hundreds of configurations. Thepresence of disconnected diagrams also precludes the use of partially-twisted boundary conditionsand therefore limits our choices of the renormalization momentum scale to the allowed latticemomenta.The quark field renormalization Z q is also computed in the RI-SMOM scheme via the ampu-tated vertex function of the local vector current operator, ¯ q γ µ q , from which we compute Z V / Z q where Z V is the corresponding renormalization factor for the local vector current. The factor Z V isnot unity as the local vector current is not conserved, however it can be computed independentlyfrom the ratio of hadronic matrix elements containing the local and conserved (five-dimensional)vector current allowing Z q to be obtained from the above. Alternatively, Z q can also be computedfrom the local axial-vector current operator ¯ q γ µ γ q . Again the ratio Z A / Z q is determined from athree-point function evaluated in momentum space and, providing non-exceptional kinematics areused, is equivalent up to negligible systematic effects at large momentum [38]. The quantity Z A isthen determined by comparing the pion-to-vacuum matrix elements of the local and approximatelyconserved (five-dimensional) axial current.The independent projection matrices P m contract the external spin and color indices, and arechosen with a tensor structure that reflects that of the operator with the same index. For the weakeffective operators, we can choose both parity-even and parity-odd projectors, which project ontothe parity-even and parity-odd components of the amputated Green’s function, respectively, andwhich should both provide the same result due to chiral symmetry. In practice however we havefound that the parity-odd choices are better protected against residual chiral symmetry breakingeffects that induce non-zero mixings between the different SU ( ) L ⊗ SU ( ) R representations (cf.Sec. 4.5 of Ref. [39]), and so we will use the parity-odd projectors exclusively. We consider twodifferent projection schemes: the “ γ µ scheme”, for which the parity-odd projectors have the spinstructure, (55) P γ µ m = ± γ µ ⊗ ( γ γ µ ) − ( γ γ µ ) ⊗ γ µ , and the “ / q scheme” with spin structure (56) P / qm = ± / q ⊗ ( γ / q ) − ( γ / q ) ⊗ / q . For the full set of parity-odd and parity-even projectors we refer the reader to Sec. 3.3.2 ofRef. [33].3Similar choices of γ µ and / q projector exist also for the quark field renormalization. We willfollow the convention of describing our RI-SMOM schemes with a label of the form SMOM ( A , B ) where the quantities A and B in parentheses describe the choices of projector for the four-quarkoperator and Z q , respectively. In this work we consider only the SMOM ( γ µ , γ µ ) and SMOM ( / q , / q ) schemes as previous studies of the renormalization of the neutral kaon mixing parameter B K in-dicate that the non-perturbative running is better described by perturbation theory for these twochoices than for the two mixed schemes [40]. We will compare our final results obtained usingboth intermediate schemes in order to estimate the systematic perturbative and discretization errorsin computing the RI to MS matching. B. Operator mixing
The seven weak effective operators mix with several dimension-3 and dimension-4 bilinearoperators. For the parity-odd components these are S = ¯ s γ d , S = ¯ s → / D γ d and S = ¯ s ← / D γ d ,where the arrow indicates the direction of the discrete covariant derivative. These are accountedfor by performing the renormalization with subtracted operators, (57) Q (cid:48) sub , lat i = Q (cid:48) i + ∑ j = b j S lat j . The subtraction coefficients b j are obtained by applying the following conditions, (58) P β α j (cid:68) s ( − p ) ¯ d ( p ) O sub , lat i ( q ) (cid:69) αβ amp . = q . The projection operators can be found in Sec. 7.2.6 ofRef. [37]. In practice we find that the subtraction coefficients are small due to the suppression ofthe mixing by a factor of the quark mass as a result of chiral symmetry, and also the observation thatthe amputated vertex function Eq. (53) with a four-quark external state and a two-quark operatornecessarily involves only disconnected diagrams that are small at large momentum scales due tothe running of the QCD coupling.Mixing also occurs with the dimension-5 chromomagnetic penguin operator and a similar elec-tric dipole operator, conventionally labeled Q and Q , respectively [41]. These operators donot vanish by the equations of motion and therefore contribute also to the on-shell matrix ele-ments, but break chiral symmetry and as such are expected to be heavily suppressed [41, 42]. Itis therefore conventional to neglect their effects in, for example, the determination of the Wilson4coefficients [7]. In our DWF calculation the dimension-1 mixing coefficients of these dimension-5operators will be of order the input quark masses used in our RI-SMOM calculations or the DWFresidual mass — effects, when combined with the required gluon exchange, should be at or belowthe percent level. Thus, in this work we neglect these operators.In addition to the lower-dimension operators there is also mixing with both gauge-invariant andgauge-noninvariant dimension-6 two-quark operators. These operators enter at next-to-leadingorder and above, and are therefore naturally small provided we perform our renormalization atlarge energy scales.The gauge-noninvariant dimension-6 operators vanish due to gauge symmetry and in manycases also by the equations of motion, and therefore do not contribute to on-shell matrix ele-ments [43]. These operators enter the renormalization only at the two-loop level [8] and above,and given that the RI → MS matching factors are at present only available to one loop, the sys-tematic effect of disregarding these operators is likely to be much smaller than our dominantsystematic errors. Nevertheless we are presently investigating position-space renormalization [44]which does not require gauge fixing and therefore does not suffer from such mixing, and as suchwe may be able to remove this systematic error in future work.Of the gauge-invariant dimension-6 operators, (59) G = ¯ s (cid:2) D µ [ D µ , D ν ] (cid:3) γ ν ( − γ ) d is the only operator that mixes at one loop [45], with all others entering at two-loops and above. InRef. [37] we have investigated the impact of including the G operator in our RI-SMOM renormal-ization and have computed the subsequent effect on the K → ππ amplitudes. This can be achievedwithout the need for measuring matrix elements of G between kaon and ππ states by taking ad-vantage of the equations of motion to rewrite those matrix elements for on-shell kinematics interms of the matrix elements of the conventional four-quark operators, such that the entire effectof this operator is captured by changes in the values of the ( , ) elements of the renormalizationmatrix. Note that at present the results including the G operator have been computed only at the2.29 GeV renormalization scale and not the 4.0 GeV scale used for our final result. However, asdemonstrated in Ref. [37] and also in Sec. VII F, the effects of including G are at the few per-cent level as expected, implying that the resulting systematic error is small compared to our othererrors.5 C. Step-scaling
Step-scaling [36] allows for the circumvention of the upper limit on the renormalization scaleimposed by the lattice spacing through independently computing the non-perturbative running ofthe renormalization matrix to a higher scale using a finer lattice. The multiplicative factor relatingthe RI-SMOM operators renormalized at two different scales can be obtained from the ratio (60) Λ RI ( µ , µ ) = Z RI ← lat ( µ )( Z RI ← lat ( µ )) − , where µ is a renormalization scale that lies below the cutoff on the original coarser lattice while µ is a higher scale, likely inaccessible on the coarser lattice. The quantity Λ RI ( µ , µ ) is computedon finer lattices for which µ also lies below the cutoff and can be applied thus, (61) Z RI ← lat ( µ ) = Λ RI ( µ , µ ) Z RI ← lat ( µ ) in order to raise the renormalization scale to µ , giving the renormalization matrix Z RI ← lat ( µ ) which non-perturbatively converts our course-lattice operators into an RI scheme defined at a scale µ potentially much larger that the inverse of our coarse lattice spacing. We will take advantageof this technique to avoid having to match perturbatively to MS directly at the lower energy scalesallowed by our coarse, a − = .
38 GeV lattice.
D. Details and results of lattice calculation
We use the step-scaling procedure to obtain the renormalization matrix at a scale of µ = . β = . a − = . ( ) GeV (32ID) ensemble and a second,finer ensemble with β = .
37 and a − = . ( ) whose properties are described in Ref. [21]under the label “32Ifine”. These ensembles have periodic spatial boundary conditions rather thanG-parity boundary conditions, but as previously mentioned, boundary effects can be neglected forthese high-energy Green’s functions. Such quantities are also constructed to be insensitive to thequark mass scale, and therefore we can disregard the unphysically heavy 170 MeV and 370 MeVpion masses on the 32ID and 32Ifine ensembles, respectively. Note also that, although we do nottake the continuum limit of the step-scaling matrix computed on the 32Ifine ensemble, the finelattice spacing and the typically small size of discretization effects on such quantities [46] suggestthe induced error is also negligible compared to our other errors. We remind the reader that thesecalculations do not include the G operator, and its absence in our calculation is treated as a sourceof systematic error in Sec. VII.6Due to the presence of disconnected diagrams in our calculation, the choices of quark momentaare restricted to the discrete values allowed by the finite-volume. The closest match betweenallowed momenta on the 32ID and 32Ifine ensembles that can be chosen as an intermediate scaleis µ = .
531 GeV and µ = .
514 GeV, respectively. The fact that these scales differ by1.1% introduces a systematic error that, given the slow evolution of the QCD β -function, can betreated as negligible.We obtain the quark field renormalization for the 32Ifine ensemble via the vector current op-erator as described in Sec. V A. For the 32ID ensemble we use the axial-vector operator as thecorresponding renormalization factor, Z A has been measured to much higher precision than Z V (0.05% versus 1.2%, respectively) [47]. The measurements of Z A and Z V are treated as statis-tically independent from those of the amputated vertex functions and are incorporated into thecalculation using the superjackknife technique.On the 32ID ensemble we extend the calculation at µ = .
531 GeV performed in our previ-ous work and documented in Ref. [33] from 100 to 234 configurations, where for each configura-tion we have increased the number of stochastic sources used in the evaluation of the disconnecteddiagrams from 1 to 20, improving the statistical errors substantially. We measure the amputatedGreen’s function Eq. (53) with quark momentum choices (62) p = ( , , , ) π L , p = ( , , , ) π L , that satisfy symmetric kinematics p = p = ( p − p ) = ( µ ) . Combined with the followingmeasurements of the quark field renormalization coefficient in the γ µ and / q schemes at µ ,(63) Z γ µ q ( µ ) = . ( ) , Z / qq ( µ ) = . ( ) , we obtain the renormalization matrices Z RI ← lat i j for the SMOM ( γ µ , γ µ ) and SMOM ( / q , / q ) schemesgiven in Tab. VI.For the measurement of the step-scaling matrix on the 32Ifine ensemble we likewise use (64) p = ( , , , ) π L , p = ( , , , ) π L , µ = .
514 GeV and (65) p = ( , , , ) π L , p = ( , , , ) π L , at the high scale µ = .
006 GeV. The corresponding values of Z q are (66) Z γ µ q ( µ ) = . ( ) , Z / qq ( µ ) = . ( ) , at µ = .
514 GeV and (67) Z γ µ q ( µ ) = . ( ) , Z / qq ( µ ) = . ( ) , at µ = .
006 GeV.The results for the step-scaling matrix Λ ( .
006 GeV , .
514 GeV ) i j in both schemes are givenin Tab. VII. In Tab. VIII we combine these step-scaling results with the 32ID Z RI ← lat results toproduce the final renormalization matrices at 4.0 GeV, where the errors on the two independentensembles have been propagated using the super-jackknife procedure.As mentioned previously, we will also utilize step-scaled renormalization matrices computedat µ = .
29 GeV both with and without the G operator included. This calculation used anintermediate scale of µ = .
33 GeV to match between the coarse and fine ensemble. Details ofthis calculation can be found in Ref. [37]. In that work the statistical errors on Z V and Z A werenot included in the results, and Z V was used rather than Z A in the determination of Z q on the 32IDensemble. In order to match the procedure outlined above we have reanalyzed the data from thatwork, the results of which are presented in Tab. IX for µ = .
33 GeV and Tab. X for µ = . ( γ µ , γ µ ) scheme are available with G included. VI. RESULTS FOR A A A AND ε (cid:48) ε (cid:48) ε (cid:48) In this section we combine our lattice measurements with experimental inputs to obtainRe ( ε (cid:48) / ε ) . The set of Standard Model parameters and other experimental values used for these8 . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) TABLE VI: The elements of the 7 × ( γ µ , γ µ ) (upper) and SMOM ( / q , / q ) (lower)renormalization matrices Z ( . ) RI ← lat i j with renormalization scale µ = .
531 GeVcomputed on the 32ID ensemble.calculations are listed in Tab. XI and their uncertainties are accounted for as a systematic errorin the following section. In this table the value of Re ( A ) was obtained from the experimentalmeasurement of K + → π + π decays, and the value of Re ( A ) from K S → π + π − and K S → π π decays. The relationship between the isospin amplitudes and the experimental branching fractionsand decay widths is described in detail in Secs. III.A and III.B of Ref. [14].As previous mentioned, the Wilson Coefficients that incorporate the short distance physics“integrated out” from the Standard Model are known in perturbation theory in the 10-operatorbasis to NLO in the MS scheme. Partial calculations at NNLO are available in the literature [48–52], together with a preliminary study on a direct lattice determination [53]; in this manuscript weutilize the complete NLO results of Ref. [7] in the MS-NDR scheme for our central values, andthe LO predictions to assign a systematic error due to the truncation of the perturbative series.9 . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) TABLE VII: The elements of the 7 × ( γ µ , γ µ ) (upper) and SMOM ( / q , / q ) (lower)step-scaling matrices Λ ( . , . ) i j between renormalization scales µ = .
514 and µ = .
006 GeV computed on the 32Ifine ensemble.For consistency with the NLO determination of the Wilson coefficients we follow Ref. [7] inutilizing the 2-loop determination of α s given in Ref. [7] (and the 1-loop determination for the LOWilson coefficients used to estimate the systematic error) despite the fact that a 4-loop calculationis available [54]. In order to fix the parameters of the 2-loop (1-loop) calculation, a value of α s at a reference scale is required, and to minimize the perturbative truncation error it is desirablethat this scale be close to the typical scale of the physical problem, in our case O ( ) . Wetherefore utilize the 4-loop calculation of α s to run the value of α N f = s ( M Z ) given in Tab. XI downto 1.7 GeV in the 4-flavor theory, and use the result, (68) α N f = s ( . ) = . . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) TABLE VIII: The elements of the 7 × ( γ µ , γ µ ) (upper) and SMOM ( / q , / q ) (lower)renormalization matrices Z ( . ) RI ← lat i j with renormalization scale µ = .
006 GeVcomputed by applying the step-scaling matrices in Tab. VII with the renormalization matrices inTab. VI. This matrix converts the lattice matrix elements computed in this paper to theappropriate RI scheme at µ = .
006 GeV
A. Lellouch-L ¨uscher factor
The Lellouch-L¨uscher factor F [12] removes the leading power-law finite-volume correctionsto the lattice matrix element. It is defined as (69) F = π m K E ππ k (cid:18) k d δ dk + q d φ dq (cid:19) , where δ is the I = ππ scattering phase shift and φ is a known function [11] of q = Lk π , appro-priately modified for our antiperiodic pion boundary conditions [13], with k the interacting pionmomentum defined via (70) k = (cid:18) E ππ (cid:19) − m π . . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) TABLE IX: The elements of the 7 × ( γ µ , γ µ ) renormalization matrix Z ( . ) RI ← lat i j with (upper) and without (lower) the effects of the G operator included. This matrix converts thelattice matrix elements computed in this paper to the SMOM ( γ µ , γ µ ) scheme at µ = .
33 GeVNote that Eq. (69) differs by a factor of two from the corresponding equation in Ref. [12] due toour different conventions on the decay amplitude (cf. Ref. [31]).The calculation of the Lellouch-L¨uscher factor requires the derivative of the phase shift withrespect to interacting pion momentum, or correspondingly the ππ energy, evaluated at the kaonmass. The determination of this derivative is detailed in Sec. III E where we present values ob-tained both directly from the lattice and also from the dispersive prediction. Given the goodagreement between our phase shifts and the dispersive predictions [18] we will use the dispersiveresult given in Eq. (22). The variation in the results will be incorporated as a systematic error inSec. VII D.We find (71) F = . ( ) , . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) − . ( ) . ( ) − . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) TABLE X: The elements of the 7 × ( γ µ , γ µ ) renormalization matrix Z ( . ) RI ← lat i j with (upper) and without (lower) the effects of the G operator included. This matrix converts thelattice matrix elements computed in this paper to the SMOM ( γ µ , γ µ ) scheme at µ = .
29 GeVwhere the error arises primarily from the uncertainty in measured ππ energy and its small sizeresults from the small contribution of the ππ scattering phase shift relative to that of the knownfunction φ in Eq. (69). B. Renormalized physical matrix elements
The infinite-volume matrix elements of the seven chiral-basis operators Q (cid:48) Rj in a scheme R at thescale µ can be expressed without ambiguity in terms of the matrix elements M (cid:48) lat j = (cid:104) ππ | Q (cid:48) lat j | K (cid:105) of the corresponding lattice operators: (72) M (cid:48) R ( µ ) = Z R ← lat ( µ ) (cid:16) a − F M (cid:48) lat (cid:17) , where a is the lattice spacing, Z R ← lat ( µ ) a 7 × F the Lellouch-3 Quantity Value G F × − GeV − V ud V us φ ε τ | ε | ω ( A ) expt . ( ) × − GeV (†)Re ( A ) expt . ( ) × − GeV (†) m c ( m c ) m b ( m b ) m W ( m W ) m Z ( m Z ) m t ( m t ) α N f = s ( m Z ) α / . ( ) (*)sin ( θ W ) TABLE XI: Standard Model and other experimental inputs required to determine A andRe ( ε (cid:48) / ε ) from the lattice matrix elements. The parameters given in this table were obtained fromthe PDG Review of Particle Physics [6], apart from those of Re ( A ) , Re ( A ) and their ratio, ω ,which were taken from Ref. [1]. Here φ ε is the phase of the indirect CP-violation parameter ε .The CKM ratio τ = − V ∗ ts V td / V ∗ us V ud is obtained using the Wolfenstein parameterization expandedto eighth order, with parameters taken from the aforementioned review. The impact upon ourresult of the errors on those quantities marked with a ( ∗ ) is incorporated as a systematic error inSec. VII H. The errors on those quantities marked with ( † ) are included within the quotedstatistical errors on our results. The errors on the remaining quantities are neglected as theircontributions to our final error are small in comparison to our statistical error.4L¨uscher factor obtained in Eq. (71).The ten conventional, linearly-dependent operators Q i are defined in terms of the seven inde-pendent operators Q (cid:48) j as follows: (73) Q i = ∑ i T i j Q (cid:48) j , where 1 ≤ i ≤ j runs over the set { , , , , , , } and the matrix T is given by (74) T = / / /
10 0 − / − which can be found as Eqs. (58) and (59) of Ref. [34]. This relationship applies both to RI schemeand bare lattice operators.In our lattice calculation we have evaluated the matrix elements of all ten linearly-dependentoperators Q i as given in Tab. V. This gives us a consistency test of the three Fierz identities: theseidentities are obeyed to within statistical errors and with an absolute size at the 1% level, validatingour code. We do not expect the Fierz relations to be obeyed to floating point accuracy since our useof all-to-all propagators introduces a stochastic element into the inversion of the Dirac operator andour use of γ hermiticity differs between the ten operators introducing statistical noise in differentways into each evaluation.Since the Fierz identities are not obeyed exactly by the data in Tab. V, we have a choice as tohow the ten linearly-dependent matrix elements M lat i in that table are to be combined to give theseven independent matrix elements M (cid:48) lat i needed on the right-hand side of Eq. (72). To this endwe choose to treat M (cid:48) lat i as fit parameters whose best fit values are obtained by minimizing thecorrelated χ : (75) χ = ∑ i j = (cid:32) M lat i − ∑ k = T ik M (cid:48) lat k (cid:33) ( C − ) i j (cid:32) M lat j − ∑ (cid:96) = T j (cid:96) M (cid:48) lat (cid:96) (cid:33) . × − × − × − × − -0.0002789 -0.0003660 -0.0010780.001206 0.0004747 0.003463 0.002873 0.008692 0.004380 -0.0006516 -0.001387 -0.0008054 -0.00032950.0001964 0.0008078 0.003764 0.006152 0.004380 0.02195 -0.001279 -0.006099 -0.0003987 -0.0013770.0004749 0.0004188 -0.0001617 6.055 × − -0.0006516 -0.001279 0.002804 0.003961 0.001241 0.00060637.289 × − × − -0.0009426 -0.0003660 -0.0008054 -0.0003987 0.001241 0.0004238 0.002475 0.0003710-2.695 × − TABLE XII: The 10 ×
10 covariance matrix C i j between the unrenormalized, infinite-volumelattice operators in the conventional basis and physical units of GeV .The result is an optimal combination that provably minimizes the statistical error on the resulting M (cid:48) lat i . The 10 ×
10 covariance matrix C i j is estimated by studying the variation of the bootstrapmeans of the matrix elements, and is given in Tab. XII. Note that we use the same covariancematrix for the fit to each bootstrap sample (a frozen fit) and therefore do not take into account inour errors the fluctuations in the covariance matrix over bootstrap samples. However such effectsare expected to be minimal due to our large number of configurations. The results for the barematrix elements obtained by this procedure, along with those obtained by applying Eq. (73) toconvert those results back into the 10-basis, are given in Tab. XIII. These results are quoted inphysical units and incorporate the Lellouch-L¨uscher finite-volume correction.The results for the seven operators converted to the SMOM ( γ µ , γ µ ) and SMOM ( / q , / q ) schemesare given in the left two columns of Tab. XIV. The right two columns of that table show the matrixelements of the ten conventional operators in the MS scheme obtained from the left two columnsby an application of Eqs. (73) and (74). For the convenience of the reader in utilizing these resultswe also provide the covariance matrices for the SMOM ( / q , / q ) scheme matrix elements, which wewill use as our central values in Sec. VIII, and also the MS matrix elements derived from them, inTabs. XV and XVI, respectively.6 i Q (cid:48) i ( GeV ) Q i ( GeV ) . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) − . ( )
10 - 0 . ( ) TABLE XIII: The bare lattice matrix elements in the 7-operator chiral basis (second column) thatminimize the correlated χ Eq. (75), and those results converted back into the 10-operator basisby applying Eq. (73) (third column). These results are quoted in physical units and incorporatethe Lellouch-L¨uscher finite-volume correction. The errors are statistical, only.
C. Results for A A A We can now obtain A from our lattice calculation as follows: (76) A = G F √ V ∗ us V ud ∑ i = (cid:16) z MS i ( µ ) + τ y MS i ( µ ) (cid:17) M MS i ( µ ) . The Wilson coefficients have been computed to next-to-leading order in QCD and electroweak per-turbation theory in the MS scheme [7], and at µ = .
006 GeV take the values given in Tab. XVII.For the CKM matrix element ratio τ we use the value given in Tab. XI. Combining these withthe MS-renormalized matrix elements obtained in Tab. XIV we obtain the following for theSMOM ( / q , / q ) intermediate scheme, (77a)Re ( A ) = . ( ) × − GeV , (77b)Im ( A ) = − . ( ) × − GeV . and for the SMOM ( γ µ , γ µ ) intermediate scheme, (78a)Re ( A ) = . ( ) × − GeV , (78b)Im ( A ) = − . ( ) × − GeV . i SMOM ( / q , / q ) (GeV ) SMOM ( γ µ , γ µ ) (GeV ) MS via SMOM ( / q , / q ) (GeV ) MS via SMOM ( γ µ , γ µ ) (GeV )1 0 . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( )
10 - - 0 . ( ) . ( ) TABLE XIV: Physical, infinite-volume matrix elements in the SMOM ( / q , / q ) and SMOM ( γ µ , γ µ ) schemes at µ = .
006 GeV given in the 7-operator chiral basis, as well as those convertedperturbatively into the MS scheme at the same scale in the 10-operator basis. The errors arestatistical only. . . × − − . × − . − . . . . × − . − . × − . . . × − − . − . × − − . × − . . . − . × − . . . . . . − . × − − . − . . . . . − . − . . . × − − . × − − . × − − . . . . − . . − . − . . . TABLE XV: The 7 × ( / q , / q ) scheme in the chiral basis.The values of Re ( A ) agree to 4 . ( . ) % between the two schemes, and those of Im ( A ) to3 . ( . ) %. This excellent agreement suggests that the systematic errors resulting from discretiza-tion effects and the truncation of the perturbative series in the non-perturbative renormalizationare minimal at our high 4 GeV scale. In the following section a more detailed discussion of these8 . . × − . . . . . − . × − . − . . × − . . . . . . × − . − . . . . . . . . − . × − − . − . − . . . . . . . − . × − . × − − . − . . . . . . . − . × − − . − . − . . . . . . . − . − . − . − . . . × − − . × − − . × − − . × − − . . . . . − . × − . − . . × − − . − . . . . × − . . − . − . − . − . − . . . × − . . × − − . . − . − . − . − . . . . × − . TABLE XVI: The 10 ×
10 covariance matrix between the renormalized, infinite-volume matrixelements in the MS scheme in the chiral basis obtained using the SMOM ( / q , / q ) intermediatescheme. i y i z i TABLE XVII: The MS Wilson coefficients (cid:126) y and (cid:126) z at µ = .
006 GeV computed via NLOQCD+EW perturbation theory.systematic errors is presented.The contributions of each of the ten operators to the real and imaginary parts of A are given inTab. XVIII. The result for Im ( A ) is dominated by the Q matrix element with a 14(4)% cancela-tion from Q , where the errors are statistical only and the value is obtained using the SMOM ( / q , / q ) intermediate scheme to match the scheme used for the previous work. This is in contrast to the9 Re( A ) Im( A )i ( / q , / q ) ( × − GeV) ( γ µ , γ µ ) ( × − GeV) ( / q , / q ) ( × − GeV) ( γ µ , γ µ ) ( × − GeV)1 0 . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) . ( ) . ( ) . ( ) . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( . ) × − − . ( . ) × − − . ( ) − . ( )
10 2 . ( ) × − . ( ) × − − . ( ) − . ( ) Total 2 . ( ) . ( ) − . ( ) − . ( ) TABLE XVIII: The contributions of each of the ten four-quark operators to Re ( A ) and Im ( A ) for the two different RI-SMOM intermediate schemes. The scheme and units are listed in thecolumn headers. The errors are statistical, only.51(29)%-level cancelation observed in Ref. [1] and is largely due to a 5 . σ increase in the Q con-tribution from − . ( ) × − GeV to − . ( ) × − GeV (again using the SMOM ( / q , / q ) intermediate scheme). This change appears to largely result from excited-state contamination inour previous result, as we can see in Fig. 6 comparing the (larger-statistics) single-operator resultat the value of t (cid:48) min = t (cid:48) min =
5. This suggests that the 5% systematic error we formerly associated with excited-statecontamination was significantly underestimated.
D. Incorporating experimental results to improve the determination of Im( A A A ) The real and imaginary parts of A comprise different linear combinations of the same basis ofreal lattice matrix elements. As the real part of the amplitude is precisely known from experimentand is not expected to receive significant contributions from new physics, we can use this quantityto replace part of the lattice input and thereby improve the precision of the imaginary part. The0appropriate procedure is discussed in Refs. [55, 56] in the context of the conventional basis of 10non-independent operators, where the latter authors use it to eliminate the Q matrix element. Forour purpose it is more convenient to express the method in terms of the unrenormalized matrixelements in the 7-operator basis. We writeRe ( A ) = G F √ V ∗ us V ud ∑ k = Re ( w MS ← lat k ) M (cid:48) lat k (79)Im ( A ) = G F √ V ∗ us V ud ∑ k = Im ( w MS ← lat k ) M (cid:48) lat k (80)where the M (cid:48) lat j = (cid:104) ππ | Q (cid:48) k | K (cid:105) are the matrix elements of the unrenormalized lattice operators inthe 7-basis in infinite-volume and physical units, andRe ( w MS ← lat k ) = ∑ i = ∑ j = (cid:16) z MS i + Re ( τ ) y MS i (cid:17) T i j Z MS ← lat jk (81)Im ( w MS ← lat k ) = ∑ i = ∑ j = (cid:16) Im ( τ ) y MS i (cid:17) T i j Z MS ← lat jk (82)are the “lattice Wilson coefficients”. Here T i j is the 10 × Z MS ← lat is the product of the 7 × × Q (cid:48) (cid:96) from Im ( A ) if wewrite Im ( A ) = G F √ V ∗ us V ud ∑ k = Im ( w MS ← lat k ) M (cid:48) lat k + λ (cid:34) Re ( A ) − G F √ V ∗ us V ud ∑ k = Re ( w MS ← lat k ) M (cid:48) lat k (cid:35) (83)and choose (84) λ = Im ( w MS ← lat (cid:96) ) Re ( w MS ← lat (cid:96) ) In Tab. XIX we present values for Im ( A ) obtained through using this procedure to replacesuccessive lattice matrix elements. The most significant gain in statistical error is achieved byreplacing the matrix element M (cid:48) lat3 , for which we obtain the following for the SMOM ( / q , / q ) inter-mediate scheme, (85)Im ( A ) = − . ( ) × − GeV1 i SMOM ( / q , / q ) ( × − GeV) SMOM ( γ µ , γ µ ) ( × − GeV)1 − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( . ) − . ( . ) − . ( . ) . ( . ) − . ( . ) − . ( . ) . ( . ) . ( . ) TABLE XIX: Values of Im ( A ) obtained for each of the two intermediate schemes by eliminatinglattice data for the matrix element of operator Q (cid:48) (cid:96) in favor of experimental value for Re ( A ) .and for the SMOM ( γ µ , γ µ ) intermediate scheme, (86)Im ( A ) = − . ( ) × − GeVwhich have 6% smaller statistical errors.We could instead choose the parameter λ to give that result for Im ( A ) with the smallest sta-tistical error. Since the value obtained for λ from this procedure is extremely close to that neededto remove the matrix element M (cid:48) lat3 , we adopt the simpler procedure of eliminating M (cid:48) lat3 and theresults given in Eqs. (85) and (86). E. Determination of ε (cid:48) ε (cid:48) ε (cid:48) Re ( ε (cid:48) / ε ) can now be obtained via Eq. (2). We use the lattice values for the I = I = ππ scattering phase-shifts : δ is given in Eq. (18) and for δ we use (87) δ = − . ( . )( . ) ◦ , obtained from our continuum result [2]. Here the parentheses list the statistical error and an esti-mate of the excited-state systematic error, respectively.Writing ε = | ε | e i φ ε , where both | ε | and its phase φ ε can be found in Tab. XI, the overall complexphase of ε (cid:48) / ε is (88) ie i ( δ − δ ) e − i φ ε = e i ( δ − δ + π / − φ ε ) . The resulting real part of the complex phase, (89)cos ( δ − δ + π / − φ ε ) = . ( ) , ( A ) and Re ( A ) ,and use the results for Im ( A ) given in Eqs. (85) and (86) that incorporate the experimental valueof Re ( A ) . The continuum, lattice value for Im ( A ) is given in Eq. 64 of Ref. [2] and must becorrected for the 20% change of Im ( τ ) = − . ( A ) = − . ( . ) × − GeVFor the SMOM ( / q , / q ) intermediate scheme we find (91)Re ( ε (cid:48) / ε ) = . ( ) and for the SMOM ( γ µ , γ µ ) intermediate scheme, (92)Re ( ε (cid:48) / ε ) = . ( ) , where the error is statistical only.It is illustrative to break the value of Re ( ε (cid:48) / ε ) into the so-called “QCD penguin” (93)Re (cid:18) ε (cid:48) ε (cid:19) QCDP = − ω cos ( δ − δ + π / − φ ε ) √ | ε | Im A ReA and “electroweak penguin” (94)Re (cid:18) ε (cid:48) ε (cid:19) EWP = ω cos ( δ − δ + π / − φ ε ) √ | ε | Im A ReA contributions, the sum of which is equal to Re ( ε (cid:48) / ε ) . These terms have opposite sign such that thesum involves an important cancellation. For the electroweak penguin contribution we find (95)Re (cid:18) ε (cid:48) ε (cid:19) EWP = − . ( ) × − . Using the results for Im ( A ) obtained using the SMOM ( / q , / q ) intermediate scheme we find (96)Re (cid:18) ε (cid:48) ε (cid:19) QCDP = . ( ) , and likewise for the SMOM ( γ µ , γ µ ) intermediate scheme, (97)Re (cid:18) ε (cid:48) ε (cid:19) QCDP = . ( ) . ( / q , / q ) and SMOM ( γ µ , γ µ ) results, respectively. This degree of cancellationis considerably less than the 71(36)% observed in our 2015 analysis. Here the errors are statisticalonly.We can also compute a purely lattice value of Re ( ε (cid:48) / ε ) using Re ( A ) from Eqs. (77a) and (78a),Im ( A ) from Eqs. (77b) and (78b), and both Re ( A ) and Im ( A ) from Eq. 64 of Ref. [2]. Note wedo not correct Re ( A ) for the change in Re ( τ ) as its contribution is much smaller than that of theWilson coefficients z i . For the SMOM ( / q , / q ) intermediate scheme we obtain (98)Re ( ε (cid:48) / ε ) = . ( ) and for the SMOM ( γ µ , γ µ ) intermediate scheme, (99)Re ( ε (cid:48) / ε ) = . ( ) , where the errors are again statistical. Unfortunately these pure-lattice results have considerablylarger statistical errors, which suggests that there is little statistical correlation between the resultsfor Im ( A ) and Re ( A ) which would be needed to reduce the error in their ratio. Thus, we will usethe results given in Eqs. (91) and (92) for our final results. F. Origin of the change in ε (cid:48) ε (cid:48) ε (cid:48) compared to our 2015 calculation In this section we provide further insight into the origin of the significant change between our2015 result of Re ( ε (cid:48) / ε ) = . ( . )( . ) × − and our results above. Several factors maycontribute to this effect:1. The increase in the minimum time separation between the four-quark operator and the sink ππ operator from 4 to 5 in the K → ππ matrix element fitting.2. The change in the procedure for determining the derivative with respect to energy of ππ scattering phase-shift that enters the Lellouch-L¨uscher factor.3. The increase in statistics from 216 to 741 configurations.4. The addition of the ππ ( ) and σ sink operators.5. The use of step-scaling to raise the renormalization scale from 1.53 GeV to 4.01 GeV.46. The change in the value of the experimental inputs, notably that of the CKM ratio τ from0 . − . i to 0 . − . i .We first note that repeating the ππ two-point function analysis for our larger data set but with aone-state fit to a single operator ( ππ ( ) ), and a fit range 6-25 to match that of the 2015 analysis,yields a result (in lattice units), (100) A ππ ( ) = . ( ) × √ × E = . ( ) that is consistent with the results of our 2015 analysis, (101) A ππ ( ) = . ( ) × √ × E = . ( ) to 1 . σ and 1 . σ for the amplitude and energy, respectively. Furthermore, the p-value of this fitis 0.451 indicating an excellent fit to the one-state model. The ground-state energy is, however,significantly larger than the value of E = . ( ) found using three operators and two statesin Sec. III D.We next repeat the analysis of the K → ππ matrix elements but with only the ππ ( ) operatorand a one-state fit with t (cid:48) min = ππ fit parameters fromEq. (100) above. Recall t (cid:48) min is the minimum time separation between the four-quark operatorand the ππ sink for data included in the fit. We use the same input experimental parameters andother analysis strategies as in the original work, including the approach to obtaining the Lellouch-L¨uscher parameter and the same SMOM ( / q , / q ) non-perturbative renormalization factors with µ = .
529 GeV. We find, (102)Re ( ε (cid:48) / ε ) = . ( . ) × − , where the errors are statistical only. This result is completely consistent with our 2015 result,(103)Re ( ε (cid:48) / ε ) = . ( . ) × − , indicating that a 3.4 × increase in statistics is not sufficient to account for the difference.Repeating the above but with the K → ππ analysis and input parameters updated to match thatof the present work gives, (104)Re ( ε (cid:48) / ε ) = . ( . ) × − , which is slightly larger but still considerably smaller than the results in the previous section. Withthe step-scaled renormalization factors with µ = .
01 GeV we find, (105)Re ( ε (cid:48) / ε ) = . ( . ) × − . ππ and K → ππ fittingstrategies. Adopting the final fit ranges determined for the ππ and K → ππ fits in Secs. III and IV,such that the analysis now differs only in the number of ππ operators, gives (106)Re ( ε (cid:48) / ε ) = . ( . ) × − . This result is now much closer to our final result. The behavior we observe here is consistent withthat displayed in Fig. 6 where we plot the dependence of the fitted matrix elements on the cut t (cid:48) min and the number of ππ operators included in the fits to the matrix elements (the ππ two-pointfunction fits remain unchanged between the results displayed in this figure). This figure showsa significant discrepancy between the Q matrix element obtained from a one-operator, one-statefit with t (cid:48) min = t (cid:48) min = t (cid:48) min = Q matrix element is still several standard deviations larger than the strong plateau observed in themulti-operator fits.We therefore conclude that the difference in Re ( ε (cid:48) / ε ) between our present and 2015 analysisresults can be attributed primarily to unexpectedly large excited-state contamination in our previ-ous analysis masked by the rapid reduction in the signal to noise ratio, and that multiple operatorsare essential to isolate the ground-state matrix element even with large statistics. VII. SYSTEMATIC ERRORS
In this section we describe the procedure used to estimate the systematic errors on our results.We will quote the values as representative percentage errors on either the matrix elements or on A as appropriate. A discussion of the systematic errors in the ∆ I = / A. Excited state contamination
In Sec. IV E we devoted considerable effort to finding an optimal fit window in which excitedstate effects are minimal. We were unable to find evidence of such effects arising from excitedkaon states, which is to be expected given both the large relative energy of such states and also6 i Rel. diff1 − . ( ) . ( ) − . ( . ) − . ( ) . ( ) . ( ) . ( ) − . ( ) . ( ) − . ( ) TABLE XX: Relative differences between the ground-state elements obtained by fitting to 3operators and 3 states with t (cid:48) min = t (cid:48) min = ππ operator implies that the data furthest from the kaon operator dominates the fit results. Assuch we do not assign a systematic error to possible contamination from excited kaon states.As for the contribution of excited ππ states, we found little evidence for such effects evenwithin the single operator fits to the ππ ( ) data, except for the Q and Q matrix elementswhere the single-operator fits showed statistically significant deviations from the common plateauregion that did not die away until t (cid:48) =
6. We observed that by adding further sink operators andallowing for more ππ states substantially reduced the excited-state contamination such that the fitswere highly consistent even if we include data at times as low as t (cid:48) =
3. Despite this we chose aconservative uniform cut of t (cid:48) min = ππ states, we consider thecomparison of the 3-operator, 3-state fit with t (cid:48) min = t (cid:48) min = t (cid:48) min = . σ , is for the Q matrix element, which has only a negligible contribution toIm ( A ) . For the dominant Q and Q matrix elements the differences cannot be resolved withinour errors. We therefore conclude that the excited state systematic error is likely to be muchsmaller than our dominant systematic errors and can be neglected. B. Unphysical kinematics
As our values of E ππ and m K differ by 2.2(3)%, the K → ππ matrix elements are not preciselyon shell. As discussed in Sec. IV, the primary result of these unphysical kinematics is the riseof a divergent contribution from the pseudoscalar operator ¯ s γ d that vanishes when on shell bythe equations of motion. In order to suppress this error we perform an explicit subtraction of thepseudoscalar operator that leaves behind a finite, regulator-independent term that represents thedominant remaining systematic error from the unequal kaon and ππ energies. As we are close tobeing on shell we can reasonably assume a linear ansatz for the dependence of our result on theenergy difference E ππ − m K . We estimate the associated systematic error by observing the changein the Q matrix element as the kaon mass is increased by 4.5%. The measurement was performedusing 69 configurations of our original ensemble [1], with 3 different K → π time separations (10,12, and 14), and we observed a 6.9% increase in the matrix element. We scale this increase by therelative difference between our kaon and ππ energies, giving 3%.Another means of estimating this systematic error is to vary the subtraction coefficients α i by anamount consistent with the expected size of the residual contribution of the pseudoscalar operator.Given that the operator is dimension-3, its coefficient is originally O ( m s / a ) where the strangequark mass is in physical units. After the subtraction is performed, the residual term is expectedto be of size O ( m s Λ ) , which has a relative size of ∼ a Λ , or ∼ Λ QCD =
300 MeV. Increasing the subtraction coefficients α i by this amount gives riseto the differences in the unrenormalized lattice matrix elements given in Tab. XXI. The observedvariations are generally consistent with the above, but to be conservative we assign a relative sys-tematic error of 5% on the matrix elements resulting from the off-shell difference E ππ (cid:54) = m K .8 i Rel. diff1 − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) − . ( ) TABLE XXI: Relative differences in the unrenormalized lattice matrix elements of Q i as thepseudoscalar subtraction coefficients α i are uniformly increased by 5% C. Finite lattice spacing
We use the value provided in Ref. [1] as an estimate of the finite lattice spacing systematicerror. This was obtained by comparing the values of the ∆ I = / × β = .
75 (32ID) lattice.The parameters of the latter ensemble are identical to those used in this work to compute A , albeitwithout G-parity boundary conditions and with a larger-than-physical light quark mass givinga unitary pion mass of 170 MeV. The MS values for the three continuum matrix elements thatcontribute to A are obtained by combining the continuum values of those matrix elements in theSMOM ( / q , / q ) scheme (Tab. XIV of Ref. [2]) with the RI → MS renormalization matrix computedon the 32ID lattice (Eq. 66 of Ref. [2]). As such this estimate addresses only the discretizationerrors on the matrix elements and not to those on the renormalization factors (which are expectedto be small). We find the values given in Tab. XXII. Averaging the three relative errors we arriveat an estimate of 12% discretization errors on the matrix elements.9
Operator 32ID continuum rel. diff ( , ) ( , ) ( , ) mix TABLE XXII: The three ∆ I = / µ = . that contribute to A , calculated on the 32ID ensemble (Ref. [14], Eq. (31)) and inthe continuum limit (Ref [2], Tab. XIV) along with their relative difference. Only statisticalerrors are shown. D. Lellouch-L ¨uscher factor
As described in Sec. VI A, the calculation of the Lellouch-L ¨uscher factor, F , that accounts forthe power-law finite-volume corrections to the matrix element, requires an ansatz for the derivativeof the ππ phase shift with respect to energy. In Sec. III E we present values for this derivativeobtained from three methods: • The Schenk parameterization [57] of the dispersive energy dependence obtained in Ref. [16] • A linear approximation in the ππ energy above threshold, d δ dE ππ = δ E ππ − m π , which is inspiredby the dispersive low-energy dependence found in Ref. [16] and can be related to d δ / dq via Eq. (70). • A direct lattice calculation of the phase shift at energies close to and including the kaonmass.Ignoring the noisier of the two lattice determinations, the results varied between d δ d q = .
26 and1.41, a 12% spread. The resulting values of F differ by 1.5% since the dominant contributionarises from the derivative of the analytic function φ . We therefore assign a 1.5% systematic errorto the matrix elements from this source. E. Exponentially-suppressed finite volume corrections
We expect the remaining finite volume corrections to our matrix elements to be dominated bythe (exponentially-suppressed) interactions between the final state pions that are not accounted for0 i Relative difference1 − . ( ) − . ( ) . ( ) − . ( ) . ( ) . ( ) ( ) ( ) − . ( ) − . ( ) TABLE XXIII: The relative difference in MS matrix elements at µ = .
33 GeV obtained throughthe SMOM ( γ µ , γ µ ) intermediate scheme due to including the G operator.by the L ¨uscher and Lellouch-L ¨uscher prescriptions. In Refs. [2, 14] we performed an in-depthanalysis of the finite-volume errors on the matrix elements that comprise A using SU(3) chiralperturbation theory, in which the mesonic loop integrals are replaced by discrete sums over theallowed momenta. We do not expect these effects to depend strongly on the form of the four-quark operator, and indeed comparable O ( − . ) corrections were estimated for both classesof operator that enter the calculation of A . We therefore assign a representative 7% systematicerror to the matrix elements. F. Neglecting the contribution of the G G G operator In the calculation of our step-scaled non-perturbative renormalization factors with scale µ = .
01 GeV we have not incorporated the effects of the G operator. A previous lattice study [37],performed in the SMOM ( γ µ , γ µ ) scheme and utilizing step-scaling from a low-scale of µ = . .
29 GeV on a finer lattice, revealed the effects on A of including this operator to be on the order of a few percent when combined with the matrixelements measured in our 2015 work [1]. Unfortunately the statistical errors on the differencesin the renormalized matrix elements at µ = .
29 GeV with and without G included were found1to be too large to resolve the effect with any precision, and we find that this also applies to thematrix elements obtained in the present work. (The renormalization matrices with and without G at µ = .
29 GeV can be found in Tab. X.)As discussed in Ref. [37], the increase in the relative error on the bootstrap differences isassociated largely with the step-scaling matrix Λ RI that describes the running between the low andhigh energy scales. However it is reasonable to expect that the largest effects of neglecting G appear at the low energy scale in the step-scaling where the QCD coupling is larger. We thereforecompare the matrix elements renormalized at the low scale in the MS scheme in order to estimatethe size of this systematic error with greater precision. We perform this comparison using theSMOM ( γ µ , γ µ ) intermediate scheme with µ = .
33 GeV, the renormalization matrices of whichare given in Tab. IX. The relative differences of the resulting MS matrix elements are given inTab. XXIII. While the observed differences are still poorly resolved, the typical size of the effectappears to be O ( ) , and we therefore assign a 3% systematic error to the effect of neglecting G . (This estimate is quite conservative given the tiny difference in the dominant, Q operatorobserved in the table.) G. Sytematic errors in
MSMSMS operator renormalization
The most important systematic errors in determining the renormalization matrix Z MS ← lat arisefrom three sources: i) The omission of dimension-6, quark bilinear operators which vanish onshell such as G discussed above. ii) Finite lattice spacing errors that result from our large choiceof RI renormalization scale µ . iii) The perturbative truncation error introduced when one-loopQCD perturbation theory is used to relate the RI-SMOM and MS schemes. In order to estimatethese systematic errors, we examine the difference between the results in the MS scheme obtainedfrom our two different intermediate RI-SMOM schemes. Rather than examining the matrix ele-ments themselves, which can be statistically noisy and vary significantly in size and importance,it is convenient to study instead the differences between the elements of the 7 × → MSrenormalization matrix R MS ← RI ← lat1 - loop j (cid:96) ( µ ) = H MS ← RI1 - loop jk ( µ ) R RI ← lat k (cid:96) ( µ ) , (107)where H is the perturbative matching matrix. In the absence of systematic errors the matrix R MS ← RI ← lat is independent of the intermediate RI scheme. We can then study this systematic2error by examining the matrix Ξ ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) I − (cid:20) R MS ← SMOM ( / q ,/ q ) ← lat1 - loop (cid:21) − R MS ← SMOM ( γ µ , γ µ ) ← lat1 - loop (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (108)where I is the 7 × | . | implies that the absolute value of each element is taken.The ratio of R -matrices in this equation converts from the lattice scheme to MS through one inter-mediate scheme, and converts back to the lattice scheme via the other scheme, and hence becomesthe unit matrix if no systematic errors exist. The difference from the unit matrix is therefore ameasure of the size of the systematic error: Under the reasonable assumption that the systematicerrors in the two schemes are comparable in size, we expect the elements of Ξ to vary betweenzero and approximately twice the size of the systematic error present in each. We therefore assigna percentage systematic error that is one half of the largest observed element of Ξ at a scale µ .In Tab. XXIV we tabulate the non-zero elements of Ξ for various MS scales and step-scalingprocedures. Once again we observe that the effects of including or discounting the G operator,while harder to statistically resolve after passing through the step-scaling procedure, are at thepercent scale.As expected there is a general trend towards smaller values as we increase the scale that ap-pears consistent with the factor of three decrease in α s between 1.33 GeV and 4.01 GeV that isexpected to describe the scaling of the missing NNLO terms. Unfortunately the statistical errors onthe results at 4.01 GeV are too large to resolve the residual systematic effects. Nevertheless, con-sidering the results of this table and also the 3-4% differences observed in Re A and Im A betweenthe schemes in Sec. VI C, we assign a 4% systematic error to the non-perturbative renormalizationfactors. H. Parametric errors
We propagate the parametric uncertainties shown in Tab. XI to Re A and Im A . For Re A thelargest such uncertainty is the charm-mass dependence, which, however, is only a 0 .
3% effect.For Im A , the largest uncertainty is 5% from the τ parameter, 3% from α s , and less than 1% fromthe charm and top quark masses. The other uncertainties have been estimated but are negligiblecompared to those quoted. We therefore estimate a total parametric uncertainty of 6% for Im A and 0 .
3% for Re A .3 Element ( i , j ) . ( ) . ( ) . ( ) . ( ) (2,2) 0 . ( ) . ( ) . ( ) . ( ) (2,3) 0 . ( ) . ( ) . ( ) . ( ) (2,5) 0 . ( ) . ( ) . ( ) . ( ) (2,6) 0 . ( ) . ( ) . ( ) . ( ) (3,2) 0 . ( ) . ( ) . ( ) . ( ) (3,3) 0 . ( ) . ( ) . ( ) . ( ) (3,5) 0 . ( ) . ( ) . ( ) . ( ) (3,6) 0 . ( ) . ( ) . ( ) . ( ) (5,2) 0 . ( ) . ( ) . ( ) . ( ) (5,3) 0 . ( ) . ( ) . ( ) . ( ) (5,5) 0 . ( ) . ( ) . ( ) . ( ) (5,6) 0 . ( ) . ( ) . ( ) . ( ) (6,2) 0 . ( ) . ( ) . ( ) . ( ) (6,3) 0 . ( ) . ( ) . ( ) . ( ) (6,5) 0 . ( ) . ( ) . ( ) . ( ) (6,6) 0 . ( ) . ( ) . ( ) . ( ) (7,7) 0 . ( ) . ( ) . ( ) . ( ) (7,8) 0 . ( ) . ( ) . ( ) . ( ) (8,7) 0 . ( ) . ( ) . ( ) . ( ) (8,8) 0 . ( ) . ( ) . ( ) . ( ) TABLE XXIV: The non-zero elements of the matrix Ξ computed using the renormalizationmatrices obtained at µ = .
33 GeV and 1.53 GeV on the 32ID ensemble, as well as thestep-scaled renormalization matrices with µ = .
29 GeV and 4.01 GeV. We do not include the G operator here, and its absence is treated as a separate systematic error in Sec. VII F. I. Wilson coefficients
As mentioned previously we compare the NLO and LO determinations of the Wilson coeffi-cients in order to estimate the systematic error arising due to missing higher-order terms. More4specifically, we compare Im ( A ) obtained from LO and NLO Wilson coefficients, computed usingthe 1-loop and 2-loop determinations of α s , respectively, while keeping fixed the renormalized ma-trix elements in the MS scheme at 4.01 GeV obtained using the SMOM ( / q , / q ) intermediate scheme,given in Tab. XIV, together with the various input parameters, such as the quark masses and theQCD coupling constant. For the latter we use the solution of the 4-loop β function [54] to compute α N f = s ( ˆ µ ) in the 4-flavor theory, starting from the value of α s ( m Z ) in Tab. XI, and we study thedependence of the LO prediction of Im ( A ) as a function of ˆ µ , relative to the NLO result. (Asexpected, the NLO shows a mild dependence simply due to the mismatch between the running of α s from the Z pole (4 loops) and the running used in the calculation of the Wilson Coefficients (2loops).) Starting at 8% at ˆ µ ≈ m c , it increases up to 16% at ˆ µ ≈ m b ; hence for our systematic errorestimate on the Wilson coefficients, we choose the intermediate point ˆ µ = . Λ N f = leads to similarconclusions.Additionally we consider the same difference of LO vs NLO predictions for Im ( A ) , as a func-tion of the RI intermediate schemes and the scale of the RI to MS conversion, while keeping fixedall parameters, α N f = s ( ˆ µ ) included. We find that, despite varying the renormalization scale by al-most a factor of two and the use of different intermediate RI schemes, the differences in the valuesof Im ( A ) are quite consistent, in the range 11-15%. This suggests that the bulk of the observeddifference arises from the perturbative 3-to-4 flavor matching and running above the charm thresh-old, which is common to all of these determinations, and that improved theory input for the 3-to-4flavor matching could significantly reduce it. (Note that in our calculation we take the matchingscale across a flavor threshold equal to the corresponding quark mass in order to avoid large loga-rithms. Additional insights could be gained by studying the dependence on this matching scale asin Ref. [52].)In conclusion we assign a 12% systematic error on both Re A and Im A associated with theNLO determination of the Wilson coefficients. J. Error budget
We divide the systematic errors into those that affect the calculation of the matrix elementsof the MS weak operators Q (cid:48) j and those that enter when these matrix elements are combined toproduce the complex, physical decay amplitude A . The former are collected in Tab. XXV. In5 Error source ValueExcited state -Unphysical kinematics 5%Finite lattice spacing 12%Lellouch-L¨uscher factor 1.5%Finite-volume corrections 7%Missing G operator 3%Renormalization 4%Total 15.7% TABLE XXV: Relative systematic errors on the infinite-volume matrix elements of theMS-renormalized four-quark operators Q (cid:48) j . Error source ValueRe ( A ) Im ( A ) Matrix elements 15.7% 15.7%Parametric errors 0.3% 6%Wilson coefficients 12% 12%Total 19.8% 20.7%
TABLE XXVI: Relative systematic errors on Re ( A ) and Im ( A ) .order to obtain the final systematic error on Im ( A ) arising from these matrix elements we notethat the result is dominated by the Q operator with only a 20% cancellation from Q . In thiscircumstance it is reasonable simply to apply the same flat percentage error to Im ( A ) as to Q .Since Re ( A ) is similarly dominated by Q , we apply the same strategy. For A we then arrive atthe error budget given in Tab. XXVI which includes this error arising from the uncertainties in thematrix elements as well as those arising from the use of perturbation theory when computing theMS Wilson coefficients and the values of the needed Standard Model input parameters.6 i SMOM ( / q , / q ) (GeV )1 0 . ( )( ) − . ( )( ) . ( )( ) − . ( )( ) − . ( )( ) . ( )( ) . ( )( ) TABLE XXVII: Physical, infinite-volume matrix elements in the SMOM ( / q , / q ) scheme at µ = .
006 GeV given in the 7-operator chiral basis. The errors are statistical and systematicrespectively. Note that our 4% estimate of the renormalization systematic error includes bothlattice systematic errors and those associated with the truncation of the perturbative series in theRI → MS matching. While the latter are inappropriate to apply to matrix elements in thenon-perturbative schemes, due to our estimation procedure we are at present unable to isolatethese two effects and as such apply the full 4% systematic error also to these RI matrix elements.
VIII. FINAL RESULTS AND DISCUSSION
In this section we collect our final results including systematic errors and discuss the implica-tions of our results. For consistency with our previous work we will use the SMOM ( / q , / q ) interme-diate scheme for our central value. A. Matrix elements
The renormalized, infinite-volume matrix elements in the RI and MS schemes are given inTab. XIV, where the errors are statistical only. The corresponding relative systematic errors canbe found in Tab. XXV. For the convenience of the reader we have reproduced the matrix elementsin the SMOM ( / q , / q ) scheme including their systematic errors in Tab. XXVII. In order to allow thereader to compute derivative quantities from these matrix elements, the covariance matrices forthe renormalized matrix elements in the SMOM ( / q , / q ) and MS schemes at 4.01 GeV can be foundin Tabs. XV and XVI, respectively.7 B. Decay amplitude
For the real part of the decay amplitude we take the value from Eq. (77a) and apply the system-atic errors given in Tab. XXVI to obtain (109)Re ( A ) = . ( . )( . ) × − GeV , where the errors are statistical and systematic, respectively. The imaginary part is obtained like-wise from Eq. (85), giving (110)Im ( A ) = − . ( . )( . ) × − GeV . The breakdown of the contributions of each of the 10 operators to these amplitudes can be foundin Tab. XVIII. We observe that, at the scale at which we are working, the dominant contributionto Re ( A ) (97%) originates from the tree operator Q , while Q has a contribution of about 13%that is largely cancelled by that of the penguin operator [58, 59] Q . Likewise, the dominantcontribution to Im ( A ) is from the QCD penguin [58, 59] operator, Q , with a 14% cancellationfrom Q . C. A comment on the ∆ I = / ∆ I = / ∆ I = / rule The “ ∆ I = / I = K → ππ decay rate relative to that of the I = ( A ) / Re ( A ) = . ( ) . A factor of two contribution to this ratio arises from the perturbativeWilson coefficients [60–62]. While the remaining factor of ten has been viewed for some timeas a consequence of the strong dynamics of QCD, the origin of this large factor has remainedsomething of a mystery with no widely-accepted dynamical explanation.In the past [14, 15, 63], and most recently in Ref. [2], when simulating with physical pionmasses we have observed a sizeable cancellation between the two Wick contractions of the dom-inant ( , ) operator contributing to the ∆ I = / ( A ) . In these calculations we reproduced the experimental value of Re ( A ) andconcluded that this cancellation was likely to be a very significant element in the ∆ I = / ( A ) depends sensitivelyon the light quark mass and becomes much less significant as the light quark mass is increasedabove its physical value. Note also that such a cancellation is not consistent with na¨ıve factoriza-tion, which predicts that both contributions have the same sign and differ in size by a factor of8three due to color suppression, although calculations using chiral perturbation theory and the 1 / N c expansion [64, 65] have previously suggested a reduction in A .In order to obtain a quantitative, first-principles result for Re ( A ) / Re ( A ) , we also requireknowledge of Re ( A ) which we provide in Eq. (109) of the present paper. Combining this withour earlier result for A [2], we obtainRe ( A ) Re ( A ) = . ( . ) ( . ) , (111)where the errors are statistical and systematic respectively. The value in Eq. (111) agrees verywell with the experimental result, demonstrating quantitatively that, within the uncertainties, the ∆ I = / D. Result for Re( ε (cid:48) / εε (cid:48) / εε (cid:48) / ε ) For ε (cid:48) / ε we use Eq. (2), combining the lattice values for the imaginary parts of the decay am-plitudes with the experimental measurements of the real parts. The systematic error for Im ( A ) is taken from Tab. XXVI and that of Im ( A ) from Eq. 64 of Ref. [2]. The statistical and system-atic errors on Im ( A ) and Im ( A ) are combined in quadrature and are therefore enhanced by thecancellation between the two terms in Eq. (2). However, one further important systematic errorshould be addressed: that arising from the effects of electromagnetism and the isospin-breakingdifference, m d − m u , between the down and up quark masses.While for most quantities these corrections enter at the 1% level or below, for ε (cid:48) this familiarsituation does not hold. As can be seen from the formula used to compute ε (cid:48) in the Standard Modelgiven in Eq. (2), the I = I = A and A enter with equal weight. However, asis summarized by the ∆ I = / A is 22.5 times smaller than A . Thus, a 1%correction to A can introduce an O ( ) correction to A and a potential correction to ε (cid:48) of 20%or more.The effects on ε (cid:48) of electromagnetism and m d − m u have been the subject of active researchfor some time [66–68]. The most recent results are those of Cirigliano et al. [68]. They providea correction that is appropriate for our calculation in which the contribution of the electro-weakpenguin operators Q and Q has been included. Their result is parametrized by ˆ Ω eff which is9introduced into a version of Eq. (2) which incorporates these effects: ε (cid:48) ε = i ω + e i ( δ − δ ) √ ε (cid:34) Im ( A emp2 ) Re ( A ( ) ) − Im ( A ( ) ) Re ( A ( ) ) (cid:0) − ˆ Ω eff (cid:1)(cid:35) . (112)and find ˆ Ω eff = ( . + . − . ) × − . Here we are reproducing Eqs. (54) and (60) from Ref. [68],where Re ( A ( ) , ) refer to the real amplitudes in the absence of isospin breaking, Im ( A emp2 ) repre-sents the dominant contribution to Im ( A ) and arises from the electroweak penguin operators Q , ,and Im ( A ( ) ) additionally includes the effects of QCD penguin operators. At the present level ofaccuracy, our use of the experimental rates for the real amplitudes, together with small differ-ences from the definition of the isosymmetric limit in Ref. [68], do not affect the applicabilityof Eq. (112) to our calculation. (For a review of earlier work on this topic see Ref. [69].) Notealso that ω + = Re ( A + ) / Re ( A ) , where the plus (+) indicates the amplitude obtained from chargedkaon decay, is equal to the value of ω used to represent the isospin-symmetric ratio in this workand given in Tab. XI.Since a careful discussion of these corrections is beyond the scope of this paper we choose totreat these effects of isospin breaking as a systematic error whose size is given by the effect ofincluding ˆ Ω eff in Eq. (112). We find (113)Re ( ε (cid:48) / ε ) = . ( )( )( ) , where the errors are statistical and systematic, with the systematic error separated as isospin-conserving and isospin-breaking, respectively. We note that if we were to apply this negativecorrection directly to our result for Re ( ε (cid:48) / ε ) , the central value obtained, 0.00167, would nearlycoincide with the experimental value, albeit with appreciable errors.Our first-principles calculation of ε (cid:48) / ε also allows us to place a new, horizontal-band constrainton the CKM matrix unitarity triangle in the ¯ ρ − ¯ η plane. In Fig. 12 we overlay this band with con-straints arising from other sources. We find that our result is consistent with the other constraintsand does not at present suggest any violation of the CKM paradigm. For more information on howthis band was obtained, as well as the corresponding plot obtained using our 2015 results, we referthe reader to Ref. [72]. IX. CONCLUSIONS
We have described in detail a calculation which substantially enhances our 2015 lattice calcu-lation of ε (cid:48) [1]. Both the 2015 and the current calculation were performed on a single, 32 × η _ ρ _ ∆ M s / ∆ M d ε K + |V cb |sin 2 β |V ub /V cb | ε ’ FIG. 12: The horizontal-band constraint on the CKM matrix unitarity triangle in the ¯ ρ − ¯ η planeobtained from our calculation of ε (cid:48) , along with constraints obtained from other inputs [6, 70, 71].The error bands represent the statistical and systematic errors combined in quadrature. Note thatthe band labeled ε (cid:48) is historically (e.g. in Ref. [72]) labeled as ε (cid:48) / ε , where ε is taken fromexperiment.M¨obius domain wall ensemble with the Iwasaki+DSDR gauge action, with an inverse lattice spac-ing of 1.378(7) GeV and physical pion masses. G-parity boundary conditions are used in the threespatial directions which induces non-zero momentum for the ground-state pions so that the energyof the lightest two-pion state matches the kaon mass to around 2%, thereby ensuring a physical,energy-conserving decay.The new calculation reported here is based on an increase by a factor of 3.4 in the numberof Monte Carlo samples and includes two additional ππ interpolating operators, which have dra-matically improved our control over contamination arising from excited ππ states. The greaterresolution among the excited finite-volume ππ states provided by our now three interpolating op-erators has allowed us to resolve the approximately 2 σ discrepancy between our earlier latticeresult for the I = ππ scattering phase shift and the dispersive prediction, as will be detailedin Ref. [18]. These improved techniques result in a significant, 70% (2 . σ if σ is determinedfrom only the statistical error) relative increase in the size of the unrenormalized lattice value of Q , suggesting that our excited-state systematic error was previously underestimated. A detailedcomparison of our old and new result can be found in Sec. VI F.We have also included in this new calculation, an improved renormalization technique. As dis-cussed in Sec. V, the lattice matrix operators must be renormalized in the MS scheme in which theWilson coefficients that parameterize the high-energy weak interactions have been evaluated. This1is accomplished by performing an intermediate non-perturbative conversion into two RI-SMOMschemes, each of which can be matched perturbatively to MS at some high energy scale. Aswe use a somewhat coarse, a − = .
38 GeV ensemble, our renormalization scale was formerlylimited by this cutoff and µ = .
53 GeV was chosen as the momemtum scale at which our RI-SMOM schemes were converted to MS. In the new calculation reported here we have applied thestep-scaling procedure to bypass the limitation imposed by the lattice cutoff and raise our renor-malization scale to 4 .
006 GeV, thereby improving our control over the systematic error resultingfrom the perturbative matching to MS. This improved method results in a reduced discrepancy be-tween the results obtained from the two different RI-SMOM intermediate schemes and a reductionin the renormalization systematic error. In the future we expect to improve this systematic errorby further raising the renormalization scale.Unfortunately raising the renormalization scale does not result in a similar improvement forthe Wilson coefficients as their error is dominated by the use of perturbation theory at the scaleof the charm quark mass to match the effective weak interaction theory between 3 and 4 flavors.We are presently working [73] to circumvent this issue by computing the 3- to 4-flavors matchingnon-perturbatively using a position-space NPR technique [44].Finally in the current calculation we have adopted a new bootstrap method [27] to determinethe χ distributions appropriate for our calculation in which the data is both correlated and non-Gaussian. The resulting improved p -values provide better guidance in our choice of fitting rangesand multi-state fitting functions.Finite lattice spacing effects remain a significant source of systematic error as at present wehave computed ε (cid:48) at a single, somewhat coarse lattice spacing. In the future we intend to follow theprocedure used in our A calculation [2] to compute A at two different lattice spacings, allowingus to perform a full continuum limit. This is hampered by the need to generate new ensembleswith GPBC, which alongside the high computational cost of the measurements and the need forlarge statistics requires significantly more computing power than is presently available.A second important systematic error, which we plan to reduce in future work, comes from theeffects of electromagnetic and light quark mass isospin breaking. As discussed in Sec. VIII D,the small size of the amplitude A relative to A gives a potential twenty times enhancement ofsuch effects which are normally at the 1% level. The effects of electromagnetism and the quarkmass difference m d − m u have been studied in considerable detail using chiral perturbation theoryand large N c arguments, most recently in Ref. [68]. We take the size of their correction as an2important systematic error for our present result and are exploring possible methods to also uselattice techniques to determine these effects [74, 75].For our final result we obtain (114)Re ( ε (cid:48) / ε ) lattice = . ( )( )( ) . The third error here is the systematic error associated with isospin breaking and electromagneticeffects, and the first and second errors are the statistical error and the remaining systematic error.This result can be compared to the experimental value (115)Re ( ε (cid:48) / ε ) expt . = . ( ) . These values are consistent within the quoted errors.We believe that ε (cid:48) continues to offer a very important test of the Standard Model with excitingopportunities for the discovery of new physics. For this promise to be realized substantially moreaccurate Standard Model predictions are needed. Important improvements can be expected froma simple extension of the work presented here, studying a sequence of ensembles with decreasinglattice spacing so that a continuum limit can be evaluated. In addition, we are developing a second,complementary approach to the study of K → ππ decay which is based on periodic boundary con-ditions. This avoids the complexity of the G-parity boundary conditions used in the present workbut requires that higher excited ππ states be used as the decay final state [76]. More challenging isthe problem posed by the inclusion of electromagnetism where new methods [74, 75] are neededto combine the finite-volume methods of L¨uscher [11] and Lellouch and L¨uscher [12] with thelong-range character of electromagnetism. ACKNOWLEDGMENTS
We would like to thank our RBC and UKQCD collaborators for their ideas and support. Weare also pleased to acknowledge Vincenzo Cirigliano for helpful discussion and explanation ofisospin breaking effects. The generation of the gauge configurations used in this work was primar-ily performed using the IBM BlueGene/Q (BG/Q) installation at BNL (supported by the RIKENBNL Research Center and BNL), the Mira computer at the ALCF (as part of the Incite program),Japan’s KEKSC 1540 computer and the STFC DiRAC machine at the University of Edinburgh(STFC grants ST/R00238X/1, ST/S002537/1, ST/R001006/1), with additional generation per-formed using the NCSA Blue Waters machine at the University of Illinois. The majority of the3measurements and analysis, including the NPR calculations, were performed using the Cori su-percomputer at NERSC, with contributions also from the Hokusai machine at ACCC RIKEN andthe BG/Q machines at BNL.T.B. and M.T. were supported by U.S. DOE grant
Appendix A: Wick contractions for the K → ππ K → ππ K → ππ three-point function with the σσσ operator In this appendix we provide the expressions for the Wick contraction required to compute the K → ππ three-point function with the σ operator. The corresponding diagrams for the ππ ( . . . ) operators can be found in in Appendix B.1 and B.2 of Ref. [33].For this appendix we will utilize the notation described in Sec. III A whereby the quark fieldoperators are placed in two-component “flavor” vectors ψ l and ψ h for the light and heavy quarks,respectively, and the corresponding propagators are matrices also in this flavor index. In thisnotation the creation operator for the G-parity even neutral kaon analog has the form, (A1) O ˜ K = i √ ( ¯ d γ s + ¯ s (cid:48) γ u )= i √ ψ l γ ψ h , where the physical component corresponds to the usual neutral kaon operator (cf. Sec. VI.A ofRef. [10]). The σ creation operator has the form, (A2) O σ = √ ( ¯ uu + ¯ dd )= √ ψ l ψ l . For convenience we will treat the meson bilinears as point operators in which both quarks resideon the same lattice site. (In our actual lattice calculation we use more elaborate source and sink4operators but those details are not needed to specify how we evaluate the Wick contractions.) Theten effective four-quark operators Q i for i ∈ { . . . } written in the above notation are given inSec. 3.2.2 of Ref. [33]. While the exact forms are not important for this discussion, we highlightthe fact that the operators are written in terms of a common set of matrices, (A3) M µ , V ± A = F γ µ ( ± γ ) , M µ , V ± A = − F γ µ ( ± γ ) , where F i are diagonal flavor matrices that pick out either the upper (0) or lower (1) element of thevector upon which they act: (A4) F = , F = . The matrices M µ i , V ± A appear inside products of two bilinear operators and the space-time index µ is summed over implicitly. Following the notation of Ref. [33] we will suppress this index.The Wick contractions of the K → ππ three-point function with the σ operator, (A5) A i = (cid:104) | O † σ ( z ) ˆ Q i ( y ) O ˜ K ( x ) | (cid:105) , where ˆ Q i are the unsubtracted four-quark operators, are divided into three classes by their topol-ogy that we label with indices 1, 3 and 4 in homage to the conventional labeling of the ππ ( . . . ) contractions. The type3 and type4 diagrams are those that contain a quark loop at the location ofthe four-quark operator, with type4 corresponding to that subset of those diagrams that are dis-connected (i.e. for which the σ operator self-contracts). For the ππ ( . . . ) operators the remaining,connected, contractions can be subdivided based on whether the two pion bilinear operators aredirectly connected by a quark line ( type2 ) or not ( type1 ); no such distinction exists of course forthe σ sink operator. Hence we classify all remaining diagrams as type1 .As in Ref. [33] it is convenient to write the ten expressions A i in terms of a common basis of, inthis case 23, functions D α ( Γ , Γ ) where the subscript indexes the function and Γ , are spin-flavormatrices.We will first write down the expressions for the correlation functions A i in terms of thesefunctions and will conclude the section with their definition. We list the contributions for each of5the three types separately. The type1 contributions are as follows: (A6a) A type11 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) (A6b) A type12 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) (A6c) A type13 = D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A ) (A6d) A type14 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A ) (A6e) A type15 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A ) (A6f) A type16 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A ) (A6g) A type17 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A ) (A6h) A type18 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A ) (A6i) A type19 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) (A6j) A type110 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) , the type3 contributions are: (A7a) A type31 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) (A7b) A type32 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) (A7c) A type33 = D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) (A7d) A type34 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A )+ D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) (A7e) A type35 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A )+ D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) (A7f) A type36 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A )+ D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) A type37 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) (A7h) A type38 = D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) (A7i) A type39 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) (A7j) A type310 = D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) , and the type4 : (A8a) A type41 = − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) (A8b) A type42 = − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) (A8c) A type43 = − D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A )+ D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) (A8d) A type44 = − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) (A8e) A type45 = − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) (A8f) A type46 = − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) + D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) (A8g) A type47 = − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A )+ D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) A type48 = − D ( M , V − A , M , V − A ) + D ( M , V − A , M , V − A ) − D ( M , V − A , M , V + A )+ D ( M , V − A , M , V + A ) − D ( M , V − A , M , V + A ) (A8i) A type49 = − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A )+ D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) (A8j) A type410 = − D ( M , V − A , M , V + A ) + D ( M , V − A , M , V + A ) − D ( M , V − A , M , V − A )+ D ( M , V − A , M , V − A ) − D ( M , V − A , M , V − A ) . The type1 contractions are: (A9a) D ( Γ , Γ ) = tr (cid:16) Γ G ly , x γ G hx , y Γ G ly , z G lz , y (cid:17) (A9b) D ( Γ , Γ ) = tr (cid:16) G hx , y Γ G ly , x γ (cid:17) tr (cid:16) G lz , y Γ G ly , z (cid:17) (A9c) D ( Γ , Γ ) = tr s f (cid:18)(cid:104) Γ G ly , z G lz , y (cid:105) αβ (cid:104) Γ G ly , x γ G hx , y (cid:105) αβ (cid:19) (A9d) D ( Γ , Γ ) = tr s f (cid:16) G ly , x γ G hx , y Γ (cid:17) αβ tr s f (cid:16) Γ G ly , z G lz , y (cid:17) αβ (A9e) D ( Γ , Γ ) = tr s f (cid:16) tr c (cid:104) G ly , x γ G hx , y Γ (cid:105) tr c (cid:104) G ly , z G lz , y Γ (cid:105)(cid:17) , and the type3 are: (A10a) D ( Γ , Γ ) = tr (cid:16) γ G hx , y Γ G ly , z G lz , x (cid:17) tr (cid:16) G ly , y Γ (cid:17) (A10b) D ( Γ , Γ ) = tr (cid:16) G ly , z G lz , x γ G hx , y Γ G ly , y Γ (cid:17) (A10c) D ( Γ , Γ ) = tr s f (cid:18)(cid:104) Γ G ly , z G lz , x γ G hx , y (cid:105) αβ (cid:104) Γ G ly , y (cid:105) αβ (cid:19) (A10d) D ( Γ , Γ ) = tr s f (cid:16) Γ G ly , y (cid:17) αβ tr s f (cid:16) Γ G ly , z G lz , x γ G hx , y (cid:17) αβ (A10e) D ( Γ , Γ ) = tr (cid:16) G ly , z G lz , x γ G hx , y Γ (cid:17) tr (cid:16) G hy , y Γ (cid:17) (A10f) D ( Γ , Γ ) = tr (cid:16) G hy , y Γ G ly , z G lz , x γ G hx , y Γ (cid:17) (A10g) D ( Γ , Γ ) = tr s f (cid:16) tr c (cid:104) G ly , y (cid:105) tr c (cid:104) Γ G ly , z G lz , x γ G hx , y Γ (cid:105)(cid:17) (A10h) D ( Γ , Γ ) = tr c (cid:16) tr s f (cid:104) G hy , y Γ (cid:105) tr s f (cid:104) Γ G ly , z G lz , x γ G hx , y (cid:105)(cid:17) (A10i) D ( Γ , Γ ) = tr s f (cid:16) tr c (cid:104) G hy , y (cid:105) tr c (cid:104) Γ G ly , z G lz , x γ G hx , y Γ (cid:105)(cid:17) , G l and G h are light and strange quark propagators, respectively, and α , β are color indices.We indicate spin and flavor traces as tr s f and color traces as tr c ; traces over all three indices (spin,color and flavor) are denoted as tr without a subscript.For simplicity, in Eqs. (A13) given below for the type4 diagrams we do not include the discon-nected σ “bubble”, (A11) B σ = tr (cid:16) G lz , z (cid:17) . In computing the expectation values of these diagrams it is also necessary to perform a vacuumsubtraction. Thus, the expressions D ∗ i given in Eqs. (A13) can be used to obtain the completecontributions of the corresponding diagrams to the type4 amplitudes as follows: (A12) (cid:104) D i ( Γ , Γ ) (cid:105) = (cid:104) D ∗ i ( Γ , Γ ) B σ (cid:105) − (cid:104) D ∗ i ( Γ , Γ ) (cid:105) × (cid:104) B σ (cid:105) , where D ∗ are defined as: (A13a) D ∗ ( Γ , Γ ) = tr (cid:16) G ly , x γ G hx , y Γ G ly , y Γ (cid:17) (A13b) D ∗ ( Γ , Γ ) = tr (cid:16) G hx , y Γ G ly , x γ (cid:17) tr (cid:16) G ly , y Γ (cid:17) (A13c) D ∗ ( Γ , Γ ) = tr s f (cid:18)(cid:104) Γ G ly , y (cid:105) αβ (cid:104) Γ G ly , x γ G hx , y (cid:105) αβ (cid:19) (A13d) D ∗ ( Γ , Γ ) = tr s f (cid:16) G ly , x γ G hx , y Γ (cid:17) αβ tr s f (cid:16) G ly , y Γ (cid:17) αβ (A13e) D ∗ ( Γ , Γ ) = tr (cid:16) γ G hx , y Γ G ly , x (cid:17) tr (cid:16) Γ G hy , y (cid:17) (A13f) D ∗ ( Γ , Γ ) = tr (cid:16) Γ G ly , x γ G hx , y Γ G hy , y (cid:17) (A13g) D ∗ ( Γ , Γ ) = tr s f (cid:16) tr c (cid:104) G ly , y (cid:105) tr c (cid:104) Γ G ly , x γ G hx , y Γ (cid:105)(cid:17) (A13h) D ∗ ( Γ , Γ ) = tr c (cid:16) tr s f (cid:104) G ly , x γ G hx , y Γ (cid:105) tr s f (cid:104) G hy , y Γ (cid:105)(cid:17) (A13i) D ∗ ( Γ , Γ ) = tr s f (cid:16) tr c (cid:104) G hy , y (cid:105) tr c (cid:104) Γ G ly , x γ G hx , y Γ (cid:105)(cid:17) . Appendix B: Wick contractions for matrix elements required for subtraction of the vacuum andpseudoscalar operator contributions
As described in Sec. IV it is necessary to subtract a pseudoscalar operator P = ¯ s γ d from theunsubtracted weak effective four-quark operators ˆ Q i in order to remove a divergent contribution foroff-shell terms. The subtraction and the evaluation of the corresponding coeffients, α i , require themeasurement of (cid:104) O † ππ P ˜ O ˜ K (cid:105) , (cid:104) P O ˜ K (cid:105) and (cid:104) ˆ Q i O ˜ K (cid:105) correlation functions. The vacuum subtractionof the type4 diagrams also requires evaluating the (cid:104) ˆ Q i O ˜ K (cid:105) correlation functions. Here and below9we use the shorthand (cid:104) ABC . . . (cid:105) to denote the n-point Green’s functions of the operators A , B , C ,and so on, in descending time order.It is easy to see that the A vac i = (cid:104) ˆ Q i O ˜ K (cid:105) are directly proportional to the type4 , disconnectedcontributions to (cid:104) O † ππ ˆ Q i O ˜ K (cid:105) with the ππ “bubble” removed. The results are (B1a) A vac1 = √ (cid:0) C ( M , V − A , M , V + A ) − C ( M , V + A , M , V − A ) (cid:1) (B1b) A vac2 = √ (cid:0) C ( M , V − A , M , V + A ) − C ( M , V + A , M , V − A ) (cid:1) (B1c) A vac3 = √ (cid:0) C ( M , V − A , M , V + A ) + C ( M , V − A , M , V − A ) − C ( M , V + A , M , V − A ) − C ( M , V − A , M , V − A ) + C ( M , V − A , M , V − A ) − C ( M , V − A , M , V − A ) (cid:1) (B1d) A vac4 = √ (cid:0) C ( M , V − A , M , V + A ) + C ( M , V − A , M , V − A ) − C ( M , V + A , M , V − A ) − C ( M , V − A , M , V − A ) + C ( M , V − A , M , V − A ) − C ( M , V − A , M , V − A ) (cid:1) (B1e) A vac5 = √ (cid:0) C ( M , V − A , M , V − A ) + C ( M , V − A , M , V + A ) − C ( M , V − A , M , V − A ) − C ( M , V + A , M , V − A ) + C ( M , V − A , M , V + A ) − C ( M , V − A , M , V + A ) (cid:1) (B1f) A vac6 = √ (cid:0) C ( M , V − A , M , V − A ) + C ( M , V − A , M , V + A ) − C ( M , V − A , M , V − A ) − C ( M , V + A , M , V − A ) + C ( M , V − A , M , V + A ) − C ( M , V − A , M , V + A ) (cid:1) (B1g) A vac7 = √ (cid:18) C ( M , V − A , M , V − A ) − C ( M , V − A , M , V + A ) − C ( M , V − A , M , V − A )+ C ( M , V + A , M , V − A ) − C ( M , V − A , M , V + A ) + C ( M , V − A , M , V + A ) (cid:19) (B1h) A vac8 = √ (cid:18) C ( M , V − A , M , V − A ) − C ( M , V − A , M , V + A ) − C ( M , V − A , M , V − A )+ C ( M , V + A , M , V − A ) − C ( M , V − A , M , V + A ) + C ( M , V − A , M , V + A ) (cid:19) (B1i) A vac9 = √ (cid:18) C ( M , V − A , M , V + A ) − C ( M , V − A , M , V − A ) − C ( M , V + A , M , V − A )+ C ( M , V − A , M , V − A ) − C ( M , V − A , M , V − A ) + C ( M , V − A , M , V − A ) (cid:19) (B1j) A vac9 = √ (cid:18) C ( M , V − A , M , V + A ) − C ( M , V − A , M , V − A ) − C ( M , V + A , M , V − A )+ C ( M , V − A , M , V − A ) − C ( M , V − A , M , V − A ) + C ( M , V − A , M , V − A ) (cid:19) . C − C type4 contributions from the expressionsin Sec. 3.2.2 of Ref. [33] and multiplying the result by a factor of 1 / √
3. Equivalent results canalso be obtained from the type4 contributions given in Eq. (A13) by multiplying the result by afactor of √
2. When measured with A2A propagators the results computed in these two bases arenot exactly equal due to differing choices of where to employ γ -hermiticity, a symmetry that isbroken by the stochastic “high-mode” approximation and restored only in the large ensemble-sizelimit (or the large-hit limit on a single configuration). This gives rise to the small differencesobserved in Sec. IV B.In our notation the pseudoscalar operator becomes (B2) P = ¯ s γ d = ¯ ψ h γ F ψ l , where F is defined in Eq. (A4).The (cid:104) P O ˜ K (cid:105) and (cid:104) O † ππ P O ˜ K (cid:105) correlation functions with the ππ ( . . . ) and σ operators can bewritten in terms of three diagrams: (B3a)mix3 = tr (cid:16) G lz , x γ G hx , y γ F G ly , z γ σ G lz , z γ σ (cid:17) (B3b)mix3 σ = tr (cid:16) G lz , x γ G hx , y γ F G ly , z (cid:17) (B3c)mix4 = tr (cid:16) G hx , y γ F G ly , x γ (cid:17) . where x and y are the locations of the kaon source and the operator insertion, respectively. The σ sink operator is located at z , and the coordinates of the two pion bilinear operators in the ππ ( . . . ) operators are z and z .The result for A vac , P = (cid:104) P O ˜ K (cid:105) is (B4) A vac , P = − √ . The amplitudes A ππ ( ... ) , P = (cid:104) O † ππ P O ˜ K (cid:105) for the ππ ( . . . ) operators are computed as (B5) A ππ ( ... ) , P = − √ ( B mix4 + mix3 ) where (B6) B = −
12 tr (cid:16) G lz , z γ σ G lz , z γ σ (cid:17) is the ππ self-contraction “bubble” introduced in Sec. B.2 of Ref. [33]. The corresponding resultfor the σ sink operator is (B7) A σ , P = ( B σ mix4 − mix3 σ ) B σ is defined in Eq. (A11). [1] Z. Bai et al. (RBC, UKQCD), Phys. Rev. Lett. (21), 212001 (2015), 1505.07863.[2] T. Blum et al. , Phys. Rev. D91 (7), 074502 (2015), 1502.00263.[3] M. Kobayashi and T. Maskawa, Prog. Theor. Phys. , 652 (1973).[4] J. R. Batley et al. (NA48), Phys. Lett. B544 , 97 (2002), hep-ex/0208009.[5] E. Abouzaid, M. Arenton, A. R. Barker, M. Barrio, L. Bellantoni, E. Blucher, G. J. Bock, C. Bown,E. Cheu, R. Coleman, M. D. Corcoran, B. Cox, et al. , Phys. Rev. D , 092001 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.83.092001 .[6] M. Tanabashi, K. Hagiwara, K. Hikasa, K. Nakamura, Y. Sumino, F. Takahashi, J. Tanaka, K. Agashe,G. Aielli, C. Amsler, M. Antonelli, D. M. Asner, et al. (Particle Data Group), Phys. Rev. D , 030001(2018), URL https://link.aps.org/doi/10.1103/PhysRevD.98.030001 .[7] G. Buchalla, A. J. Buras, and M. E. Lautenbacher, Rev. Mod. Phys. , 1125 (1996), hep-ph/9512380.[8] T. Blum et al. (RBC), Phys. Rev. D68 , 114506 (2003), hep-lat/0110075.[9] C. Sachrajda and G. Villadoro, Phys. Lett. B , 73 (2005), hep-lat/0411033.[10] N. Christ, C. Kelly, and D. Zhang (2019), 1908.08640.[11] M. Luscher, Nucl. Phys.
B354 , 531 (1991).[12] L. Lellouch and M. Luscher, Commun.Math.Phys. , 31 (2001), hep-lat/0003023.[13] T. Yamazaki et al. (CP-PACS), Phys. Rev. D , 074513 (2004), hep-lat/0402025.[14] T. Blum, P. Boyle, N. Christ, N. Garron, E. Goode, et al. , Phys.Rev. D86 , 074513 (2012), 1206.5142.[15] T. Blum, P. Boyle, N. Christ, N. Garron, E. Goode, et al. , Phys.Rev.Lett. , 141601 (2012),1111.1699.[16] G. Colangelo, J. Gasser, and H. Leutwyler, Nucl. Phys.
B603 , 125 (2001), hep-ph/0103088.[17] G. Colangelo, E. Passemar, and P. Stoffer, Eur. Phys. J.
C75 , 172 (2015), 1501.05627.[18] N. H. Christ et al. , Lattice determination of I = and I = ππ scattering phase shifts with physicalpion masses , forthcoming (2020).[19] D. Renfrew, T. Blum, N. Christ, R. Mawhinney, and P. Vranas, PoS LATTICE2008 , 048 (2008),0902.2587.[20] R. Arthur et al. (RBC, UKQCD), Phys. Rev.
D87 , 094514 (2013), 1208.4412.[21] T. Blum et al. (RBC, UKQCD), Phys. Rev.
D93 (7), 074505 (2016), 1411.7017. [22] Z. Bai et al. B738 , 55 (2014), 1403.1683.[24] Y.-C. Chen and T.-W. Chiu (TWQCD), PoS
IWCSE2013 , 059 (2014), 1412.0819.[25] C. Jung, C. Kelly, R. D. Mawhinney, and D. J. Murphy, Phys. Rev.
D97 (5), 054503 (2018),1706.05843.[26] E. Carlstein, Ann. Statist. (3), 1171 (1986), URL https://doi.org/10.1214/aos/1176350057 .[27] C. Kelly and T. Wang, in (2019), 1911.04582.[28] J. Foley, K. Jimmy Juge, A. O’Cais, M. Peardon, S. M. Ryan, and J.-I. Skullerud, Comput. Phys.Commun. , 145 (2005), hep-lat/0505023.[29] M. Luscher, Commun. Math. Phys. , 153 (1986).[30] J. Bijnens, G. Colangelo, G. Ecker, J. Gasser, and M. E. Sainio, Nucl. Phys. B508 , 263 (1997), [Erra-tum: Nucl. Phys.B517,639(1998)], hep-ph/9707291.[31] T. Blum et al. , Phys. Rev.
D84 , 114503 (2011), 1106.2714.[32] C. W. Bernard, T. Draper, A. Soni, H. D. Politzer, and M. B. Wise, Phys. Rev.
D32 , 2343 (1985).[33] D. Zhang,
Kaon to Two Pion decay from Lattice QCD and CP violation , Ph.D. thesis, ColumbiaUniversity (2015).[34] C. Lehner and C. Sturm, Phys. Rev.
D84 , 014001 (2011), 1104.4948.[35] R. Arthur and P. A. Boyle (RBC, UKQCD), Phys. Rev.
D83 , 114511 (2011), 1006.0422.[36] R. Arthur, P. A. Boyle, N. Garron, C. Kelly, and A. T. Lytle (RBC, UKQCD), Phys. Rev.
D85 , 014501(2012), 1109.1223.[37] G. McGlynn,
Advances in Lattice Chromodynamics , Ph.D. thesis, Columbia University (2016).[38] Y. Aoki et al. , Phys. Rev.
D78 , 054510 (2008), 0712.1061.[39] Q. Liu,
Kaon to two pions decays from lattice QCD: ∆ I=1/2 rule and CP violation , Ph.D. thesis,Columbia University (2012).[40] Y. Aoki et al. , Phys. Rev.
D84 , 014503 (2011), 1012.4178.[41] S. Bertolini, J. O. Eeg, and M. Fabbrichesi, Nucl. Phys.
B449 , 197 (1995), hep-ph/9409437.[42] A. J. Buras and J.-M. Grard, JHEP , 126 (2018), 1803.08052.[43] G. Martinelli, G. C. Rossi, C. T. Sachrajda, S. R. Sharpe, M. Talevi, and M. Testa, Nucl. Phys. B611 ,311 (2001), hep-lat/0106003. [44] M. Tomii and N. H. Christ, Phys. Rev. D99 (1), 014515 (2019), 1811.11238.[45] A. J. Buras, M. Jamin, M. E. Lautenbacher, and P. H. Weisz, Nucl. Phys.
B400 , 37 (1993), hep-ph/9211304.[46] P. A. Boyle, N. Garron, R. J. Hudspith, C. Lehner, and A. T. Lytle (RBC, UKQCD), JHEP , 054(2017), 1708.03552.[47] R. Arthur et al. (RBC, UKQCD), Phys.Rev. D87 , 094514 (2013), 1208.4412.[48] A. J. Buras, P. Gambino, and U. A. Haisch, Nucl. Phys.
B570 , 117 (2000), hep-ph/9911250.[49] M. Gorbahn and U. Haisch, Nucl. Phys.
B713 , 291 (2005), hep-ph/0411071.[50] J. Brod and M. Gorbahn, Phys. Rev.
D82 , 094026 (2010), 1007.0684.[51] M. Cerd`a-Sevilla, M. Gorbahn, S. J¨ager, and A. Kokulu, J. Phys. Conf. Ser. (1), 012008 (2017),1611.08276.[52] M. Cerd`a-Sevilla, Acta Phys. Polon.
B49 , 1087 (2018).[53] M. Bruno, C. Lehner, and A. Soni, Phys. Rev.
D97 (7), 074509 (2018), 1711.05768.[54] T. van Ritbergen, J. A. M. Vermaseren, and S. A. Larin, Phys. Lett.
B400 , 379 (1997), hep-ph/9701390.[55] A. J. Buras, M. Gorbahn, S. J¨ager, and M. Jamin, JHEP , 202 (2015), 1507.06345.[56] T. Kitahara, U. Nierste, and P. Tremper, JHEP , 078 (2016), 1607.06727.[57] A. Schenk, Nucl. Phys. B363 , 97 (1991).[58] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys.
B120 , 316 (1977).[59] F. J. Gilman and M. B. Wise, Phys. Rev.
D20 , 2392 (1979).[60] M. K. Gaillard and B. W. Lee, Phys. Rev. Lett. , 108 (1974).[61] G. Altarelli and L. Maiani, Phys. Lett. , 351 (1974).[62] F. J. Gilman and M. B. Wise, Phys. Lett. , 83 (1979).[63] P. A. Boyle et al. (RBC, UKQCD), Phys. Rev. Lett. (15), 152001 (2013), 1212.1474.[64] W. A. Bardeen, A. J. Buras, and J. M. Gerard, Phys. Lett. B192 , 138 (1987).[65] A. Pich and E. de Rafael, Phys. Lett.
B374 , 186 (1996), hep-ph/9511465.[66] A. Buras and J. Gerard, Phys. Lett. B , 156 (1987).[67] V. Cirigliano, A. Pich, G. Ecker, and H. Neufeld, Phys. Rev. Lett. , 162001 (2003), hep-ph/0307030.[68] V. Cirigliano, H. Gisbert, A. Pich, and A. Rodrguez-S´anchez, JHEP , 032 (2020), 1911.01359.[69] V. Cirigliano, G. Ecker, H. Neufeld, A. Pich, and J. Portoles, Rev. Mod. Phys. , 399 (2012),1107.6001. [70] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D98 (3), 030001 (2018 and 2019 update).[71] P. A. Boyle, L. Del Debbio, N. Garron, A. Juttner, A. Soni, J. T. Tsang, and O. Witzel (RBC/UKQCD)(2018), 1812.08791.[72] C. Lehner, E. Lunghi, and A. Soni, Phys. Lett. B , 82 (2016), 1508.01801.[73] M. Tomii, PoS
LATTICE2018 , 216 (2019), 1901.04107.[74] N. Christ and X. Feng, EPJ Web Conf. , 13016 (2018), 1711.09339.[75] Y. Cai and Z. Davoudi, PoS
LATTICE2018 , 280 (2018), 1812.11015.[76] D. Hoying, PoS