Exploring a new approach to Hadronic Parity Violation from Lattice QCD

The long-range, parity-odd nucleon interaction generated by single pion exchange is captured in the parity-odd pion-nucleon coupling $h^1_\pi$. Its calculation in lattice QCD requires the evaluation of 4-quark operator nucleon 3-point functions. We investigate a new numerical approach to compute $h^1_\pi$ based on nucleon matrix elements of parity-even 4-quark operators and related to the parity-violating electro-weak theory by PCAC and chiral perturbation theory. This study is performed with 2+1+1 dynamical flavors of twisted mass fermions at pion mass $m_{\pi} \approx 260\,\text{MeV}$ in a lattice box of $L \approx 3 \,\text{fm} $ and with a lattice spacing of $a \approx 0.091 \,\text{fm}$. From a calculation excluding fermion loop diagrams we find a bare coupling of $h^1_\pi = 8.08 \,(98) \cdot 10^{-7}$.


I. INTRODUCTION
Determining the effects of hadronic parity violation (HPV) in nucleon-nucleon interaction is a challenging task, both in experiment and theory.HPV amplitudes based on parity symmetry breaking are small deviations against a large QCD background.The long-range, single pion exchange interaction is captured at the hadronic level by the parity-violating Lagrangian of proton (p), neutron (n) and the pion triplet (⃗ π) which defines the coupling h 1 π .It originates from flavor-conserving, neutral currents at the electro-weak scale and is a promising channel to study the parity-odd pion-nucleon coupling [1].The first experimental determination [2] of the associated pion-nucleon coupling related to the ∆I = 1 effective electro-weak Lagrangian has recently sparked new interest in the theoretical Standard Model (SM) prediction of h 1 π .The available theoretical estimates of h 1 π from SM physics are predominantly based on effective field theory and model calculations, apart from one exploratory lattice calculation.Starting point for model calculations was the scheme for describing parity-nonconserving nuclear forces from Desplanques, Donoghue, Holstein [3].Continuing on that basis Dubovik and Zenkin found in Ref. [4] a best value estimate of 1.3 • 10 −7 .Ref. [5] extended the quark model picture to include the weak interaction effects from the ∆ baryon and estimated 2.7 • 10 −7 .Kaiser and Meißner started investigations with a chiral soliton model [6][7][8], and Meißner and Weigelt used a three-flavor Skyrme model and calculated the coupling in the range (0.8 − 1.3) • 10 −7 in Ref. [9].A chiral quark-soliton model was used by Ref. [10] with an estimate of the coupling 0.874 • 10 −7 .Ref. [11] applied the operator product expansion to the nucleon 2-point function in an external pion field and based on QCD-sum rules found a value 3.4 • 10 −7 .Ref. [12] studied the parity-odd couplings in the nucleon-nucleon interaction with the 1/N c -expansion and esitmated for the sin 2 (θ w )/N c -suppressed h 1 π a range 0.8 (0.3) • 10 −7 .de Vries et al. used chiral effective field theory in Ref. [13,14] to compute the neutron capture on the proton process and matched to experimental data, resulting in an estimate 1.1 (1.0) • 10 −6 .The first attempt at an estimate from first principles of the strong interaction was carried out by Wasem in Ref. [15].We come back to comparing our present work to this reference and only collect its final estimate here 1.099 (0.505) • 10 −7 .The significant experimental result by the NDPgamma collaboration in Ref. [2] of 2.6 (1.2) • 10 −7 is another milestone in this timeline.The same experimental data was subsequently re-analyzed with chiral effective field theory in Ref. [16], which estimated the coupling at 2.7 (1.8) • 10 −7 .More recently, Ref. [17] used a factorization ansatz for the matrix element of the parityviolating electro-weak Hamiltonian, together with nonperturbative lattice QCD data for the nucleon quark charges to find an estimate of 3.06 (1.72) • 10 −7 .The above results are summarized in Fig. 1.We find this situation of scattered results highlights the need for systematic study and improved ab-initio theoretical determinations.The aforementioned first ab-initio lattice QCD determination and the non-perturbative estimate of the nucleon matrix elements with the parity-violating (PV) effective Lagrangian has been presented in Ref. [15].In this work the actual transition matrix elements N π LPV −→ N for a transition of a nucleon-pion state to a nucleon state mediated by the ∆I = 1, parity-violating Lagrangian was considered.Though this calculation is pioneering, it is also exploratory in many regards as also discussed in detail in Refs.[1,18].Challenges are the rigorous treatment of the pion-nucleon state in finite volume and energy nonconservation between initial and finial state on the lattice, which were circumvented in Ref. [15].Apart from FIG. 1. Collection of estimates of the nucleon-pion coupling h 1 π ; extension of Fig. 4 in [2] and Table I in [1].The labels match to references: DZ [4], FCDH [5], Quark model [10], QCD sum rule [11], Skyrme [9], 1/Nc [12], χEFT [13], χEFT + exp [16], NPDGamma [2], LQCD [15], Factor+LEC [17].The estimate for χEFT 2015 extends beyond the shown h this the calculation considered also only a certain quark flow diagram topology, arguing that the neglected diagrams are expected to have a contribution only within the statistical accuracy.Moreover, renormalization of the 4-quark operators was not included.Still, the obtained value on a coarse lattice with a heavier-thanphysical pion m π ≈ 390 MeV is consistent with the recent NPDGgamma experimental analysis.An alternative theoretical ansatz has been put forward anew in Refs.[1,19], by proposing a joint effort of chiral effective field theory (χEFT) and lattice QCD.Based on the PCAC relation the transition via the parity-violating interaction Lagrangian with a soft pion in the initial or final state is equivalent to a transition via a parityconserving Lagrangian without a soft pion.Details on the relevant PCAC relation are provided in Refs.[1,19].The PCAC relation Eq. ( 2) is, however, exact only in the limit of exact chiral symmetry.At non-zero pion mass it receives higher order corrections in χEFT.But these corrections can be argued to be numerically small [1], at the level of O (1%) at the physical pion mass, and can in principle also be calculated by studying the pion mass dependence with lattice QCD, based on known low-energy constants from meson and heavy-baryon chiral perturbation theory [19].This alternative theoretical ansatz leads to a major simplification in the lattice computation: one now considers a transition amplitude between single nucleon states N LPC −→ N with a parity-conserving (PC) Lagrangian, which from a numerical point of view is more straightforward to handle in a lattice calculation.In particular, the complication arising from the pion-nucleon state is absent since the matrix element is computed for single nucleon initial and final states.
In this work we investigate the computational concepts proposed in Ref. [1] in practice and propose a concrete numerical implementation to evaluate the nucleon 3-point functions with the 4-quark operator insertions of L PC .Of course the ensuing Wick contractions still comprise fermion loop diagrams, which were neglected in the first work [18].We will argue that these diagrams and the renormalization procedure are intricately linked: in the lattice calculation these particular fermion loop diagrams generate power-divergent mixing with lower-dimensional operators, and we add an initial discussion about such power divergent terms and our future strategy for renormalizing the 4-quark operators.
Preliminary results of this work have been reported in [20,21].

II. OPERATORS AND COUPLING
The matching between the parity-violating interaction in the electro-weak sector of the SM and the effective nucleon and pion degrees of freedom at energy scale Λ QCD ∼ m proton has been worked out in Ref. [22].Here, we largely follow the notation of the recent Ref. [1].The ∆I = 1, parity-conserving Lagrangian is given by In the operator list Eq. (4) q = (u, d) T is the up and down quark doublet and τ 3 is the third Pauli matrix, which gives the iso-vector combination ūu − dd.As a perturbative addition to pure QCD, the interaction Lagrangian Eq. ( 3) induces a proton-neutron mass splitting (δm N ) 4q due to the 4-quark operators, and with the PCAC relation the leading contribution to the coupling comes from the mass splitting where F π denotes the pion decay constant in the chiral limit.
In the following this work is mainly concerned with the lattice QCD estimate of the operator matrix elements and N the proton and neutron.

III. 4-QUARK OPERATOR MATRIX ELEMENTS FROM THE LATTICE
To determine the nucleon matrix elements of the 4-quark operators in Eq. ( 4) we follow the Feynman-Hellmann-Theorem technique advocated for in Ref. [24].The relevant correlation functions result from inserting the individual operators θ (f )′ i (x) into the proton or neutron 2-point function, summed over all lattice sites x.The nucleons are interpolated by the usual zeromomentum, positive parity proton and neutron 3-quark operators with (positive) parity projector P (+) = 1 2 (1 + γ 0 ) and C the charge conjugation matrix.From these we construct the 2-and 3-point functions The different types of Wick contractions following from Eq. ( 9) are depicted in Fig. 2 in diagrammatic form.We distinguish three types of diagrams: those containing quark loops, denoted B and D, and without a quark loop, denoted W .Note that the latter type is the only one included in the calculation of Ref. [18].For the quark loop diagrams, we further make a technical distinction between type D, where the fermion loop is individually spin-color traced and type B, where it is not.Quark-disconnected diagrams are neglected, since by virtue of the τ 3 flavor structure of the operators in Eq. ( 4) such diagrams cancel in SU(2) flavor symmetric QCD, which we work in.We connect the 2-and 3-point functions in Eqs. ( 8), ( 9) to the nucleon matrix element by spectral decomposition and the Wigner-Eckart-Theorem Here |0⟩ denotes the QCD vacuum state , |n, σ⟩ the nucleon ground state with zero 3-momentum and spin-1/2 component σ.In Eq. ( 11) we use the spin-independent matrix element for the Lorentz-scalar operator θ . With ellipsis we denote excited state contributions as well as contributions from different time-orderings, which are at most of order O t 0 .A detailed account of the application of the Feynman-Hellmann-Theorem (FHT) to the calculation of nucleon matrix elements can be found in Ref. [24].According to FHT, in the vacuum |λ⟩ including the perturbation λ L w PC in the action we can determine the effective mass of the nucleon state for sufficiently large t by up to excited state contamination, and m eff as well as the matrix element for (δm N ) 4q are taken in pure QCD.
Taking the derivative with respect to λ we then obtain the desired matrix element by studying the dependence on source-sink separation t as well as offset τ of the ratio again up to excited state contamination.We determine R per individual operator θ by fitting the ratio to a constant for various ranges in t and τ .

Evaluation of diagrams
We evaluate the Wick contractions by a combination of point-to-all, stochastic and sequential quark propagators.The point-to-all propagators ψ (xi,α,a) result from solving the lattice Dirac equation for a spin-color diluted source with support at a single lattice site From the point-to-all propagators the nucleon 2-point functions are evaluated in the usual way.The quark loop in diagrams B and D in Fig. 2 is constructed by a fully time, spin and color diluted stochastic timeslice propagator, where η(t, ⃗ x) ∈ {±1} are independent and identically with zero mean and unit variance To apply the FHT method we must sum the 3-point function with insertion of the 4-quark operator at each lattice site.We realize this summed simultaneous insertion by using the sequential inversion method: to that end we construct the two sequential sources for B-and D-type Here Γ is one of the relevant Dirac matrices Γ = γ µ or Γ = γ µ γ 5 .The Lorentz index µ is actually summed over at this stage.
By repeated inversion of the Dirac operators on these sources we obtain the sequential propagators T (B) and T (D) for B and D diagram, respectively, given by The ensuing contractions for C 3pt are analogous to those for C 2pt , using T (B) and T (D) .For the W diagram in Fig. 2 we again use the stochastic sequential propagator technique to split the four quark lines connecting at the insertion point into two pairs.The relevant term of propagators through the insertion point then reads with binary noise vector η(x) as in Eq. ( 17), and subscript 1, 2 denoting the quark propagator flavor.We thus generate a set of independent binary noise sources η r , r = 1, . . ., N r as in Eq. ( 17), and the corresponding sequential sources and propagators and by Eq. ( 21) the product of two such sequential propagators from independent noise sources produces in the expectation value the four quark lines connected at a single site, which is summed over the lattice.

Fierz rearrangement
The operators in Eq. ( 4) fall into two classes with respect to their spin-color structure: θ have cross-linked color and spinor indices.We refer to those as "color-crossed" operators for short.
These matrix elements of the color-crossed operators can be computed in two different ways.The first one is to compute the contractions corresponding to the colorcrossed operators.This is achieved by using color dilution of the sequential sources S (W ),r in Eq. ( 22): adding color indices a, b, c we thus have The second way is to apply a Fierz rearrangement in order to transform the color-crossed operators into products of standard quark bilinear factors, thereby avoiding the need for color dilution.Then, the analogous methods discussed above for the non color-crossed operators are applied.For instance θ is equivalently represented as where by [u ↔ d] we denote the same term as written explicitly, but up replaced by down quark flavor.The strange operators are treated accordingly.

IV. LATTICE COMPUTATION
For our numerical study we use a gauge field ensemble from the Extended Twisted Mass Collaboration [25] with dynamical up, down, charm and strange quark.
The ensemble has a lattice volume of 32 3 × 64, a lattice spacing of a ≈ 0.091 fm and, thus, a spatial lattice extend of L = 3.1 fm.The pion mass value is m π = 261(1) MeV, m π • L = 4 and the nucleon mass value is m N = 1028(4) MeV.The strange quark mass is tuned to its physical value.
The simulated action features a light mass degenerate quark doublet of twisted mass fermions at maximal twist guaranteeing O(a) improvement [26], amended by a Sheikholeslami-Wohlert "clover" term included to reduce residual a 2 lattice artifacts.The charm-strange doublet is again of twisted mass type, including a quark mass splitting term [27].The heavy doublet action is not flavor diagonal, which needlessly complicates the calculation of correlation functions involving strange quarks.For the present study we thus use a mixed-action approach, by the addition of a doublet of Osterwalder-Seiler (OS) strange quarks [28], analogous to the light quark doublet.The bare OS strange quark mass value, which is not critically important yet for this exploratory investigation, has been tuned such that the Ω baryon mass assumes its physical value.The lattice action determines symmetry properties in our computation and these are relevant for the discussion of renormalization and mixing.We thus reprint the detailed formulas for sea and valence quark action in the App. A. We summarize the statistics produced with the twisted mass gauge field ensemble for the evaluation of the individual diagrams from the 4-quark operator insertion in Tab.I.

Contributions from B and D diagram type
For later discussion it will be valuable to consider the contribution to the nucleon 3-point functions for the individual operators from the combined B and D diagram type and the W -type separately, motivated by the fermion loop present in B and D-type diagrams, but not in Wtype diagrams.
To determine the estimate for the matrix element we fit the ratios R(t, τ ) to a constant in various ranges t min ≤ t ≤ t max together with various sets {τ 1 , . . ., τ n } of joint data sets.
The fits are correlated with a block-diagonal covariance matrix, where we neglect the correlation between different τ -values, i.e. we set cov (R(t, τ ) , R(t ′ , τ ′ )) = 0 for all pairs τ ̸ = τ ′ .Including cross-τ elements renders the covariance matrix near-singular and distorts the fit with bad estimates of said elements.
To each fit we assign an Akaike Information Criterion (AIC) weight following the procedure in Ref. [29], which is given by with χ 2 based on the block-diagonal covariance matrix, N param the number of fit parameters (one for our fit to constant), and N data the number of data points ( R(t, τ ) values ) entering the fit.Moreover, the fits are bootstrapped and we obtain the fit parameter uncertainty from the variance over bootstrap samples.The ranges for t and τ applied in the fits are given by The boundaries are based on observation, where a meaningful fit is accessible, and given our current accuracy the above choice covers all such ranges.Note in addition, with the symmetric ratio in Eq. ( 14), the maximal addressed ratio value involves data at t max + τ .3. Estimates of B and D diagram contribution to the matrix element from the ratios Eq. ( 14), for various pairs (t, τ ) with the weighted fit as gray band.
Finally, we restrict the set of parameters in our fits to a single constant (for the matrix element).At present level of statistical uncertainty per R(t, τ ) data point, in most t/τ ranges we cannot model excited state contamination in our data with any statistical significance.From all available fits with best fit parameter µ i and error σ i and AIC weight w i we build the combined distribution function [29] N µ,σ is the normal distribution with mean µ and variance σ 2 .Based on CDF(x) we quote the median of the distribution function as the central value and the 16% and 84% quantiles as the uncertainty interval.
In Fig. 3 we present the ratio R (ℓ) ′ (t, τ ) as a function of t for different values of τ for the B and D-type operators.
In addition we include the result of the AIC weighting procedure as gray bands.The distribution functions and quantiles corresponding to the AIC procedure are shown in Fig. 7 in App.B.

Contributions from W diagram type
The analysis of the W diagram contribution proceeds analogously to the B + D case.We show the ratio data R(t, τ ) per operator in Fig. 4 with our estimate for the matrix element as the gray band.Fig. 4 (App.B) correspondingly justifies this estimate at the level of the cumulative distribution function.

Diagrams with strange quarks
For the strange operators θ (s)′ in Eq. ( 4), only one diagram type contributes per operator.These are D-type contributions for θ (s)′ 1, 3 and B-type contributions for θ (s)′ 2, 4 , the latter only when using Fierz rearrangement for these two operators.Still, whether B or D, only strange quark loop diagrams occur in this case.Here, we circumvent the technical complication of unitary strange and charm flavor mixing by the twisted mass heavy quark action Eq.(A2) and employ the aforementioned mixed action approach with a strange quark doublet (s + , s − ), analogously to the light quark doublet.This means the strange operators matrix elements are determined similarly to the one for θ (ℓ)′ , with the replacement of the light quark loop by the strange quark  4. Estimates of W diagram contribution to the matrix element from the ratios Eq. ( 14); the meaning of symbols as Fig. 3. loop, L u (x) → L s (x).Since the two strange quark flavors in the doublet (s + , s − ) are identical in the continuum limit, we insert the strange quark loop averaged over both strange quark flavors.Using γ 5 -hermiticity, in detail we then define We show the lattice data and the matrix element fit result for the ratios R (s)′ (t, τ ) in Fig. 5.The application of the AIC cumulative distribution functions defined in Eq. ( 28) is shown in Fig. 9 in App.B.

Discussion of matrix elements from BDW diagrams
All results for the matrix elements are compiled in Tabs.II for light quark operators and III for the strange quark operators.They are given per operator flavor f = ℓ, s, operator number k and as a third label we give the diagrams contributing.
In Tab.II we observe a difference in magnitude between B + D and W diagram contributions for each individual operator by two orders of magnitude.Our explanation here mixing with operators of lower (and equal) massdimension, in the case of B and D diagrams.The latter two types contain a quark loop and we argue in Sec.VI below, that mixing of the 4-quark operator is permitted, starting with local quark-bilinear operators.By naive visual inspection, the structure of the W diagram on the other hand, does not allow for such mixing and in Sec.
V below we use solely its contribution to arrive at an estimate of h 1 π in analogy to Ref. [18].The strange quark operators θ  5. Strange operator ratios R (s)′ and matrix element fits for operators θ  that would enforce this property at the non-perturbative level.

V. BARE COUPLING FROM W -TYPE DIAGRAM
We put our present study in perspective to the work of Ref. [15] by taking an analogous approach: we only include the contribution from the W -type diagram and ignore the multiplicative renormalization and the mixing to match the lattice result to the M S scheme, while still using the Wilson coefficients from renormalized perturbation theory.Thus, we evaluate the combination of matrix elements Based on the sets of bootstrapped fits with their associated AIC weights, we construct the cumulative distribution function for M C by first building all combinations of fits for M (ℓ) i , i = 1, 2, 3, then determining the mean and error from bootstrap mean and variance of the samplewise built linear combination in Eq. (30), and finally assigning to each such combination of fits the AIC weight as the product of weights from the three individual matrix element fits.The resulting cumulative distribution function is shown in Fig. 6, together with the median (α = 0.5 quantile) and the confidence band from the 68 % CL AIC cdf median FIG. 6. Cumulative distribution function for the combined estimate of MC in Eq. ( 30).α = 0.16, 0.84 quantiles.
From the matrix element estimate M C we determine the bare coupling based on the W -type diagrams by multiplying the numerical factors from the effective Lagrangian.
In Eq. ( 31) we restored all explicit factors of the lattice spacing a to have dimensionless quantities only. 1he pion decay constant for the gauge field ensemble considered has been determined in Ref. [30] with value in lattice units af π = 0.06674 (15).The only explicit use of the lattice spacing is made to convert the Fermi constant G F as an external scale to lattice units.We use the value a = 0.09076 (54) fm from Ref. [31].
For the matrix elements combined with the Wilson coef-ficients in lattice units we find and together with the conversion factor with Standard Model parameters from PDG [32] converted to lattice units the result for the bare coupling is then The result in Eq. ( 34) as representative of a lattice estimate of h 1 π is by construction very preliminary.Conceptually, in its underlying restrictions it is similar to the first lattice determination in Ref. [18], which found at m π ≈ 389 MeV pion mass and with a coarser lattice a = 0.123 fm.We recall, that the computational ansatz of both lattice calculations differs fundamentally in using a parity-violating versus parity-conserving interaction Lagrangian.Moreover, our preliminary result for h 1 π (W, bare) is of compatible by order of magnitude with the recent experimental value

VI. COMMENTS ON MIXING AND OUTLINE OF RENORMALIZATION
The renormalization of the set of 4-quark operators in Eq. ( 4) in continuum QCD has been discussed in Refs.[22,23], and finds its expression in the Wilson coefficients C i , S (1) i which are calculated in renormalized QCD perturbation theory in the MS scheme, together with their anomalous dimension matrix.In continuum QCD with the mass-independent MS scheme and dimensional regularization, mixing with operators of lower dimension is excluded.
Here we comment on the situation in the practical lattice QCD calculation, with Wilson-type fermions and non-zero quark mass.The explicit breaking of proper Lorentz symmetry down to discrete 3-rotations and to non-equivalence of spatial and temporal direction (due to T ̸ = L) does not play a significant role for the 4-quark operators, which are invariant under a 3-and 4-dimensional rotation.
Of practical importance in this numerical calculation with Wilson-type fermions is the breaking of chiral symmetry, of up-down SU(2) flavor and parity symmetry.We focus in these introductory remarks on the light quark propagators θ (ℓ)′ .
Based on these lattice symmetries we identify the operators which are allowed to mix, i.e. which are not excluded by quantum numbers under those symmetry transformations.These quantum numbers are for the (light) 4-quark operators Operator The mixing candidate operators of mass-dimension 3 to 5 are given by dim 3 q γ 5 ⊗ 1 q 4 m ℓ q 1 ⊗ τ 3 q, q / D ⊗ τ 3 q 5 m2 ℓ q γ 5 ⊗ 1 q , m ℓ q γ 5 / D ⊗ 1 q, q γ 5 D 2 ⊗ 1 q, q σ µν Gµν q, m ℓ G G Here, 1 denotes the unit matrix in spinor and flavor space, respectively, and Gµν = ϵ µναβ G αβ the dual (lattice) gluon field strength tensor. 2 The dimension-3 operator is allowed by breaking of parity P symmetry by the twisted mass fermion action.The dimension-4 operators are allowed by chiral symmetry breaking, which in addition to having opposite R 5 parity to θ (ℓ)′ is indicated by the explicit factor of the fermion mass m ℓ .Both operators are connected by the equation of motion for the quark field.Analogous statements hold for the dimension-5 operators.
Mixing with such operator matrix elements hampers the extraction of the renormalized 4-quark operator matrix elements in the continuum limit, due to power-divergent mixing coefficients 1/a 3 and 1/a 2 for dimension-3 and -4 operators, respectively.A suitable scheme to subtract such contributions appears to be the gradient flow together with the short flow time expansion [33].The practical application to the twisted mass case is currently under investigation.We mention two other setups, which are of interest in this study of the method using parity-even 4-quark operators.

Iso-symmetric Wilson fermions
Wilson fermions observing SU(2) isospin symmetry have individual parity P and exchange symmetry S (cf.

Eqs. (A3)-(A8)
).What remains is the explicit chiral symmetry breaking by the Wilson term.In this case the dimension-3 operator is ruled out by parity symmetry.However, the dimension-4 operator m ℓ q 1⊗τ 3 q (and thus the related q / D ⊗ τ 3 q) is still allowed to mix.

Parity-odd 4-quark operators
This case was investigated in Ref. [15], with isosymmetric clover-improved Wilson fermions.In this case, using parity P, charge conjugation C and exchange symmetry S is sufficient to show, that there is no operator of dimension 3, 4 or 5, that can mix on the lattice with the operators of L w PV .
The different mixing properties in the lattice QCD calculation with Wilson-type fermion regularization appear as a major drawback of using the PCAC relation and converting to the parity-conserving Lagrangian.However, in the parity-violating case the accurate representation of the nucleon-pion state with the Lüscher method and the ensuing signal-to-noise problem from the meson-baryon interpolator pose potentially even harder problems, especially towards physical pion mass and large lattice volume.

VII. CONCLUSION
In this work we investigated the numerical implementation of a new method to calculate nucleon matrix elements of parity-even flavor-conserving 4-quark operators from lattice QCD, that pertain to the fully theoretical prediction of the long-range nucleon-pion coupling h 1 π .This constitutes the first step towards a full-fledged calculation of the coupling from a combination of chiral perturbation theory and non-perturbative lattice matrix elements.For one ensemble at pion mass 260 MeV and lattice spacing a ≈ 0.091 fm we demonstrated for the first time the calculation of all relevant (bare) matrix elements at the level of 10% combined statistical and systematic uncertainty.The specific implementation shown here is readily and feasibly scalable towards physical pion mass, continuum and infinite volume.We achieve this result due to the simplified representation of the coupling by application of the soft-pion-theorem and the implied sufficiency to calculate (single-hadron) nucleon matrix elements of parity-even operators.Of course, corrections to this leading order defining relation in χPT, though expected to be small, can in principle be computed and, therefore, the approximation is systematically improvable.However, renormalization of lattice matrix elements poses a challenge due to mixing with lower-dimensional operators.To appreciate the expected impact of this mixing we summarize our results for the bare matrix elements normalized by twice the nucleon mass in lattice units and weighted by the Wilson coefficients in Tab.IV.
In the table we keep the three contributions (ℓ) ′ i,w /(2am N ) , (ℓ) ′ i,bd /(2am N ) , (s) ′ i,bd /(2am N ) , separate, which are light operators with W -type diagrams, light operators with B + D-type diagrams and the light-strange operators with B + D-type diagrams.
From our numerical investigation we find that only W light has the expected order of magnitude, while the light and strange B + D-type diagram contributions are much larger.The latter two types contain fermionic loop subdiagrams, which we deem as the origin of the mixing.This mixing comes about due to reduced symmetries at non-zero lattice spacing and explicit breaking of chiral symmetry due to finite quark mass values.The Gradient Flow method for subtracting power-divergent mixing and renormalization appears as a promising direction to study, due to the possibility to study mixing and matching to a continuum renormalization scheme only after extrapolating lattice data to the continuum, and thus with restored symmetries.The investigation of its practical implementation for the pertinent 4-quark operators is our on-going work.
mixing strange and charm sea quark flavor by lattice artifacts.
To simplify the calculation of nucleon correlators with strange operator insertion we follow the mixed action approach in Ref. [28] and introduce another doublet of twisted mass valence strange quarks (s + , s − ).It is formally identical to the light quark doublet, except for the value of the bare twisted quark mass m ℓ → m s .
In particular it shares the critical hopping parameter (as tuned to maximal twist), and the Sheikholeslami-Wohlert parameter c SW with the light quark sector.The bare quark mass m s at maximal twist is given by twisted quark mass parameter m s = µ s and the latter is tuned, such that the mass of the Ω baryon takes the physical value.

Symmetry transformations for light quarks
We list the complete set of discrete transformations, which pertain to our identification of operator mixing for the 4-quark operators.Apart from those listed here, there are the 3-rotations, and the residual (continuous) U(1) 3 flavor symmetry, under which the lattice action is invariant.
The discrete transformations of charge conjugation C, parity P and time reversal T are given by ψ(t, ⃗ x)

FIG. 2 .
FIG. 2. Quark flow diagrams from inserting 4-quark operators into the nucleon 2-point function.The red bar denotes the single-point vertex.
from B-and D-type diagrams, albeit with the strange quark flavor running inside the fermion loop.These operators are therefore equally prone to mixing as the light quark operators.This mixing of θ(ℓ)′ k , θ (s)′ kin lattice QCD with lower dimensional operators does not come entirely unexpected, given the dimension 6 of the operators and the reduced symmetry of the lattice model.We add several comments on potential subtractions and renormalization in Sec.VI below.A second feature is the antisymmetry ofM (ℓ)′ 1 and M (ℓ)′ 3 .The opposite-equal values can be shown for the bare matrix element at tree-level of perturbation theory.Beyond that, we are currently unaware of a symmetry argument, ′ k, X /(2amN ) 1.307 +0.082 −0.094 3.93 +0.25 −0.28 −1.320 +0.097 −0.085 0.218 +0.019 −0.017 • 10 −2 1.66 +0.24 −0.21 • 10 −2 −0.218 +0.022 −0.025 • 10 −2

TABLE I .
Summary of statistics produced for the evaluation of the individual diagram types from the 4-quark operators inserted into the nucleon 2-point function.N conf gives the number of independent gauge configurations, N sample the number of stochastic samples and Nsource the number of point sources employed in the estimates.

TABLE II
. Summary of best estimates for dimensionless bare operator matrix elements per diagram type B + D ("bd") and W ("w"), divided by 2mN , from our fits and AIC analysis.

TABLE III .
Summary of best estimates for dimensionless bare strange operator matrix elements, divided by 2mN , from our fits and AIC analysis.

TABLE IV .
Summary of matrix elements per diagram type, weighted by the Wilson coefficients.