Lattice QCD determination of neutron-antineutron matrix elements with physical quark masses

Matrix elements of six-quark operators are needed to extract new physics constraints from experimental searches for neutron-antineutron oscillations. This work presents in detail the first lattice quantum chromodynamics calculations of the necessary neutron-antineutron transition matrix elements including calculation methods and discussions of systematic uncertainties. Implications of isospin and chiral symmetry on the matrix elements, power counting in the isospin limit, and renormalization of a chiral basis of six-quark operators are discussed. Calculations are performed with a chiral-symmetric discretization of the quark action and physical light quark masses in order to avoid the need for chiral extrapolation. Non-perturbative renormalization is performed, including a study of lattice cutoff effects. Excited-state effects are studied using two nucleon operators and multiple values of source-sink separation. Results for the dominant matrix elements are found to be significantly larger compared to previous results from the MIT bag model. Future calculations are needed to fully account for systematic uncertainties associated with discretization and finite-volume effects but are not expected to significantly affect this conclusion.

Matrix elements of six-quark operators are needed to extract new physics constraints from experimental searches for neutron-antineutron oscillations. This work presents in detail the first lattice quantum chromodynamics calculations of the necessary neutron-antineutron transition matrix elements including calculation methods and discussions of systematic uncertainties. Implications of isospin and chiral symmetry on the matrix elements, power counting in the isospin limit, and renormalization of a chiral basis of six-quark operators are discussed. Calculations are performed with a chiral-symmetric discretization of the quark action and physical light quark masses in order to avoid the need for chiral extrapolation. Non-perturbative renormalization is performed, including a study of lattice cutoff effects. Excited-state effects are studied using two nucleon operators and multiple values of source-sink separation. Results for the dominant matrix elements are found to be significantly larger compared to previous results from the MIT bag model. Future calculations are needed to fully account for systematic uncertainties associated with discretization and finite-volume effects but are not expected to significantly affect this conclusion.

I. INTRODUCTION
In the contemporary theory of particles and fields, there is no fundamental reason for baryon number B to be conserved. Quantum effects in the Standard Model (SM) can lead to B violation, and at temperatures above the electroweak phase transition sphaleron processes can efficiently convert baryons into antileptons while preserving (B − L), where L is lepton number. Low-temperature B-violating effects have not been observed experimentally, and their existence would have significant implications for the stability of nuclear matter. However, the observed baryon-antibaryon asymmetry of the universe cannot be explained within the SM, which fulfills Sakharov's conditions for baryogenesis [1] but does not contain enough baryon number and CP violation to reproduce the observed baryon asymmetry of the universe [2][3][4][5]. Moreover, while (B − L) symmetry is preserved in the SM, it likely has to be violated in its extensions (BSM theories) aimed at explaining baryogenesis, since electroweak sphaleron transitions would otherwise "wash out" any net baryon number generated by (B −L)-conserving interactions in the early universe.
Baryon number violation might be experimentally observed in proton decays [6] or neutron-antineutron oscillations [7][8][9][10]. The implications of these two hypothetical processes are fundamentally different: proton decay changes baryon number by |∆B| = 1 unit and involves (anti)leptons, while neutron-antineutron oscillations change baryon number by |∆B| = 2 units and do not involve leptons. Proton decay, even if observed, does not necessarily violate (B − L) and may be insufficient to explain baryogenesis.
Despite decades of searches, neither process has been observed, constraining the strength of B-violating interactions. In particular models of baryogenesis, this may require higher level of CP violation, which is in turn constrained by searches for the electric dipole moments of neutrons, nuclei, and atoms. However, excluding theories of baryogenesis using results from these experiments requires knowledge of nucleon matrix elements of B-and CP -violating effective interactions expressed in terms of fundamental fields, quarks and gluons. For neutron-antineutron transitions, these calculations have previously been performed using nucleon models [11]. Modern lattice QCD methods permit modelindependent calculation of these matrix elements. This paper reports the first completely nonperturbative calculation of the neutron-antineutron transition matrix elements computed in lattice QCD with physical quark masses and chiral symmetry. In particular, we find that lattice QCD calculations result in substantially larger n-n matrix elements compared to nucleon model calculations. Our findings imply that n-n oscillation experiments should observe 1-2 orders of magnitude more oscillation events than was previously expected for the same BSM physics parameters. This paper describes in detail our methodology for computing neutron-antineutron matrix elements of operators changing baryon number by |∆B| = 2 units, which have already been reported in a short publication [12]. In particular, the operator definitions, symmetry properties of their matrix elements, and their impact on phenomenology within SU (2) L × U (1)-symmetric extensions are discussed in Sec. II. The setup for our calculation of these matrix elements on a lattice is described in Sec. III. Extraction of ground-state matrix elements from lattice correlation functions and analysis of potential excited state contaminations are performed in Sec. V. Nonperturbative renormalization and matching to the MS scheme are described in Sec. V. The final results for n-n matrix elements and their uncertainties are provided in Sec. VI. In Section VII, we discuss briefly the impact of our results in light of other potential sources of systematic uncertainties that are not controlled in our present calculation.

II. EFFECTIVE n-n INTERACTIONS
A. Chiral basis of n-n operators A complete basis of color-singlet, electrically-neutral six-quark operators with uudddd flavor structure can be constructed from operators of the form [11,[13][14][15][16] where quark spinor indices are implicitly contracted in the parentheses, the P L,R = 1 2 (1 ∓ γ 5 ) are chiral projectors, and the quark color tensors T are T (symm) {ij}{kl}{mn} = ε ikm ε jln + ε jkm ε iln + ε ilm ε jkn + ε jlm ε ikn = T S1S2S3 , T (asym) [ij][kl]{mn} = ε ijm ε kln + ε ijn ε klm = T A1A2S3 , with S i , A i standing for the symmetrized and antisymmetrized pairs of color indices, respectively. These operators are identical in Euclidean and Minkowski spaces with the charge-conjugation spin matrix C, 1 that satisfies the usual condition Cγ µ C † = −γ T µ . Operators involving vector diquarks (q T CP χ γ µ q) or tensor diquarks (qCP χ σ µν q) are redundant and can be related to linear combinations of the operators in Eq. (1) by spin Fierz relations. The two choices of chirality for each (qCP χ q) diquark above in O 1,2,3 provide an overcomplete basis of 18 operators. Fierz relations O 2 χχχ − O 1 χχχ = 3O 3 χχχ reduce the number of independent operators to 14. 2 All 14 independent effective six-quarks operators are electrically neutral and change the baryon number by ∆B = −2 units. However, they are not independent under isospin symmetry transformations. The electroweak (EW) symmetry SU (2) L × U (1) Y requires that all interactions are SU (2) L -singlet, which may be achieved with with additional factors of the Higgs field (see Sec.II E). Furthermore, since the chiral symmetry SU (2) L ⊗ SU (2) R is preserved exactly in the massless perturbation theory, and preserved with good precision on a lattice with chiral fermions, it is more convenient to use a basis made of operators having definite values of chiral L, R-isospin.
The operators in Eq. (1) are built from color-symmetric (6 c ) and antisymmetric (3 c ) chiral diquarks, which can be denoted as where q 1,2 = u or d and the relative signs upon quark permutation come from their anticommutation, C T = −C, and color (anti)symmetry. Using the isospin doublet ψ = (u, d) and its conjugateψ = (ψ T C iτ 2 ), the chiral isoscalar and isovector diquarks can be written as where τ a are the Pauli matrices, [τ a , τ b ] = 2i abc τ c . The details of isospin classification were given in Ref. [17], and here we list only the chiral-basis operators and their relation to the conventional basis (1). All the SU (2) L -singlet operators can be constructed from some R-diquarks and L-isoscalar diquarks, resulting in three operators belonging to the (1 L , 3 R ) irreducible representation of the chiral isospin, and one (1 L , 7 R ) operator where τ ± = 1 2 (τ 1 ± iτ 2 ). The remaining 10 independent n-n transition operators are not SU (2) L singlets. Of these additional operators, three belong to the (5 L , 3 R ) irreducible representation, The remaining seven independent operators Q P 1 , · · · , Q P 7 are obtained from Q 1 , · · · , Q 7 by parity transformation discussed below (19) and belong to the (3 L , 1 R ), (7 L , 1 R ), and (3 L , 5 R ) irreducible representations. The operators Q 1 , · · · , Q 7 , Q P 1 , · · · , Q P 7 form a complete basis of 14 linearly independent SU (3) C × U (1) EM -invariant dimension-9 operators with baryon number ∆B = −2 and isospin ∆I 3 = −1 3 .
Isospin properties of the n-n operators are summarized in Tab. I, together with relations to notations used in other papers. In the following sections, we will discuss nucleon matrix elements only of the operators Q 1,2,3,5 , and the other matrix elements can be easily obtained using symmetries discussed below. 1 To avoid confusion, throughout the paper we use Euclidean γ-matrices ( γ, γ 4 ) Euc = ( γ, γ 4 ) † Euc = (−i γ, γ 0 ) M satisfying γ † µ = γµ. 2 These Fierz relations are valid in four spacetime dimensions but are violated in dimensional regularization at two-loop order [17]. The MS scheme defined in Ref. [17] includes evanescent operator counterterms that ensure that renormalized matrix elements obey these Fierz relations. Provided that matching between BSM theory and SM effective operators is consistently performed in this MS scheme or is performed at a high enough scale that one-loop QCD corrections are negligible, these Fierz relations can be assumed for MS renormalized matrix elements. 3 The isospin of operators ∆I Q is defined here as [Q, I] = ∆ I Q Q, leading to the selection rule I i − I f = ∆I Q for the isospins of initial and final states. QI Ref. [18] Ref. [11] Ref. [19] (I, I3)R ⊗ (I, I3)L γ

B. Operator mixing
In this work, we study lattice regularized operators that have to be nonperturbatively renormalized and then perturbatively converted to the MS scheme using the one-loop matching results of Ref. [17], as described in Sec. V, The renormalization matrix Z R IJ takes especially simple form in the "chiral basis" consisting of elements Q (P) I=1...7 , because they belong to different chiral multiplets and cannot mix with each other due to chiral SU (2) L × SU (2) R symmetry of massless QCD. Although some chiral representations appear in Tab. I more than once, they are actually also prevented from mixing. Specifically, operators Q (P) 5,6,7 consist of different components ("rows") of chiral 3-and 5-multiplets, and transform differently under SU (2) L × SU (2) R . Operators Q (P) 1,2,3 cannot mix with each other for a more subtle reason. Even though they belong to the same chiral representation, they contain different numbers of left-and right-handed diquarks. While the U (1) A symmetry is violated in QCD by the ABJ anomaly, operators Q 1,2,3 do not mix in perturbative QCD because perturbative gluon exchanges preserve the U (1) A transformation properties of external quark fields in their respective Green's functions. At the diagram level, there are only quark (and no antiquark) external fields, which cannot be contracted into closed loops, and thus penguin-like diagrams do not appear. This point is discussed and illustrated by an explicit two-loop perturbative calculation in Ref. [17].
In order to avoid mixing of renormalized operators, one has to define renormalization matrix Z R IJ in a scheme respecting chiral symmetry, such as MS, and perform perturbative matching calculations in massless QCD. Likewise, to avoid mixing of bare lattice operators Q lat I , chiral symmetry must be preserved in lattice QCD regularization, which requires [Möbius] domain wall ([M]DWF) or overlap fermion discretization. The MDWF action that we use in this work has been shown to have good chiral properties [20] (see Sec. III), and our lattice results may be safely matched to perturbative QCD in the UV regime. Finally, nonperturbative effects such as spontaneous chiral symmetry breaking and U (1) A -violating topological fluctuations (instantons) in the QCD vacuum could lead to operator mixing in nonperturbative renormalization (NPR). Mixing can be also induced by the light quark masses and residual chiral symmetry violation. However, as we study NPR numerically in Sec. V, we find that this mixing is negligible (≈ O(10 −3 )) and can be safely neglected at our level of precision.

C. Isospin relations between matrix elements
Since the chiral symmetry of QCD is spontaneously broken SU (2) L ⊗ SU (2) R → SU (2) L+R , the isospin selection rules for n-n matrix elements constrain only the total isospin I L+R of the effective operators Q (P) 1,··· ,7 . The n-n transition changes the isospin by ∆I 3 = −1, therefore the L, R isospins must add as The operator in the (1 L , 7 R ) representation (8) with the total isospin I L+R = 3 cannot couple a neutron to an antineutron (I L+R = ± 1 2 ) in our calculation that is performed with SU (2) f -symmetric QCD with m u = m d , therefore n|Q 4 |n | mu=m d = 0 .
Even if the isospin-breaking effects ∼ (m u − m d ) = 0 are included, such SU (2) f -violating matrix elements will be suppressed with powers of (m u − m d )/Λ QCD relative to those of other operators.
Similarly, while the Q 5,6,7 operators introduced in Eqs. (9) are linearly independent, isospin symmetry leads to additional relations between their n-n matrix elements that make two of them redundant. The relations between them are determined by the (I, I 3 ) L+R = (1, −1) component in the product of their chiral factors, as well as their normalization. To find the latter, one can use SU (2) ladder operators to construct the full 3 R and 5 L isospin multiplets starting from (u T Cu) S R ∼ (1, +1) R and (u T Cu) Combining these components to construct Q 5,6,7 according to Eq. (13) yields their relative normalizations. Taking into account the Clebsch-Gordan coefficients for the projection (11), one obtains the relations between matrix elements which are also fulfilled in lattice contractions up to the machine precision. Additionally, one can check that these relations hold, e.g., for the results of the Bag-model calculation [11] in the form D. C, P, and T relations The discrete symmetries C, P, and T , which are conserved in QCD, imply further relations for n-n transition matrix elements. Since the form of n-n operators is identical in Minkowski and Euclidean space, we study the relations between their matrix elements in Minkowski space but using Euclidean γ-matrix conventions. From the usual transformations for the fermion fields, we obtain C,P,T -transformation properties for quark bilinears and the 6-quark operators, which are summarized in Appendix A, where η C,P,T are arbitrary complex phases accompanying the C, P, T transformations of fermion fields. These factors and the relevance of Eq. (21) for CP -violating processes are discussed further in Refs. [21][22][23]. The conjugated operators Q † I are related to Q I by the CP transformation, The CP transformation also relates the transition matrix elementsn → n and n →n , which can be shown to be real.
Parity relates nn transition matrix elements of Q P I and Q I , where the phase factor is complementary to that in Eq. (19). For the conventional choice η P = 1, it is clear that only the pseudoscalar combination (Q I − Q P I ) has nonzero matrix elements, since n →n transition changes parity. Note that in all the cases, the arbitrary phase factors η C,P,T arising from the transformations of Q I cancel with the phase factors arising from the transformations of the states.
Finally, with the help of the T -reflection, one can also show that the matrix elements do not depend on the direction of the (anti)neutron spin. Using the transformation properties of the neutron and antineutron states, All spin-flip matrix elements of Q I are trivially zero because Q I are (pseudo)scalars.
Denoting the ground-state nn transition matrix elements for each Q I by the matrix element results derived above can be summarized as In conjunction with the results from Sec. II C, this implies that in the isospin limit where Eq. (28) is valid, nn transition rates involving the 14 operators Q (P) I are given in terms of 4 real nn transition matrix elements M 1,2,3,5 .

E. nn Effective Field Theory
The |∆B| = 2 effective interactions discussed above must be generated by some extension of the Standard model at yet unknown scale Λ BSM . It is generally assumed that such extensions have higher symmetry, which is broken at scales below Λ BSM to the electroweak symmetry SU (2) L × U (1) Y , and thus the effective interactions must be EWsymmetric. From the discussion above it follows that only Q 1,2,3 are SU (2) L × U (1) Y -singlets, while Q 4 , Q 5(67) , and all Q P I operators are not. These latter operators require additional EW-charged factors to make them EW-symmetric, which affect the power counting and result in higher suppression by the Λ BSM scale.
Such factors can be easily constructed from the Higgs field doublet φ and its conjugate iτ 2 φ * to compensate for the SU (2) L -and hypercharge of the operators Q where C I (µ) are dimensionless Wilson coefficients and I L [Q I . In addition to Eq. (29), the full |∆B| = 2 Lagrangian must also include combinations of electrically charged |∆B| = 2 operators with oppositely charged Higgs fields to assure the EW symmetry above the EW scale. Such interactions can lead to n ↔p and p ↔n,p transitions, and the emitted charged Higgs bosons (e.g., decaying into leptons) would compensate for the change in the electric charge. These transitions are suppressed by at least one factor of (v 2 /Λ 2 BSM ). 4 Using the effective Lagrangian (29) and the relations derived in the previous sections, the full n-n matrix element can be written as where M (P) I are the nucleon matrix elements of operators Q (P) I . The dimensionless low-energy constants C (P) I (µ) depend on the scale µ only logarithmically and can be computed perturbatively by using C I (Λ BSM ) ∼ O(1) given by a particular BSM scenario as an initial condition for renormalization group evolution. A nonperturbative calculation of the matrix elements M I is presented in the following sections.

III. LATTICE SETUP
In this section, we fist recount the details of the lattice QCD gauge configurations and propagators used in this study, and then describe the construction of (anti)neutron correlation functions with the n-n operators. The QCD gauge field configurations were generated with the Iwasaki gauge action on a 48 3 × 96 lattice and N f = 2 + 1 flavors of dynamical Möbius Domain Wall fermions. The fermion masses are tuned to be almost exactly at the physical point [20], such that the pion mass is approximately m π = 139.2(4)MeV and the scale (the lattice spacing) is a = 0.1141(3) fm. The residual mass m res , which encapsulates the residual violation of chiral symmetry, is smaller than 50% of the input quark mass. The physical lattice size L ≈ 5.45 fm and m π L = 3.86 should be sufficient to suppress finite volume effects of the n-n matrix elements to a level below our target precision. In particular, according to chiral perturbation theory, these finite size effects are expected to be < ∼ 1% [24]. The three-point functions needed to evaluate the matrix elements of the operators Q (P) I require six quark propagators for the u and d quarks flavors; they result from Wick contractions of the six-quark operators with the (anti)neutron interpolating fields. There are no disconnected quark-loop diagrams because the operators Q I (Q † I ) contain only quarks (antiquarks). For both two-and three-point lattice correlation functions we compute propagators on 30 independent gauge field configurations separated by 40 molecular dynamics time steps. All the quark propagators required for a single sample are computed from a point source located at the operator insertion point, which is identified in the analysis with the origin x 0 = (0, 0, 0, 0) using translational invariance. To reduce stochastic uncertainty, sampling of the neutron correlation functions is enhanced by all-mode-averaging [25], in which we compute 1 exact and 81 low-precision samples evenly distributed over the 4D volume on each gauge configuration. The low-precision quark propagators are computed with low-mode deflation and the conjugate gradient algorithm truncated at 250 iterations.
The propagators are contracted at the sink into intermediate baryon blocks [26,27] with polarized nucleon and antinucleon quantum numbers to minimize the time spent in the contraction step of the calculation. (Anti)neutron source and sink interpolating operators are constructed with either point or Gaussian-smeared (anti)quarks and are denoted with n J=P,S , respectively. Final contraction at the propagator source yields an (anti)neutron two-point correlation function sample with a point source at x 0 . Thus, the polarized neutron two-point correlation function with zero spatial momentum for positive time t > 0 is and, similarly, for the polarized antineutron where the polarization matrix Γ σ(±) = 1±γ4 2 1+σγ3γ5 2 projects on the selected parity (±) and spin σ = ± 1 2 , and the interpolating operator at the source is J = P and the one at the sink can be either J = P or J = S. Neutron/antineutron two-point functions have the spectral representation where the overlap factors Z J m are identical for neutrons and antineutrons in either spin orientations. The three-point functions involve two neutron or two antineutron fields to create and annihilate states with opposite baryon numbers. Using the (anti)neutron states defined in Eq. (A14), one can express the three-point correlation function containing, for example, the n ←n transition matrix element (see Fig. 1): where n (−) −σ is nucleon interpolating field that creates an antineutron with spin σ 5 , both the (anti)neutron operators are summed over the spatial coordinate to project on zero momentum. By calculating quark propagators with point sources located at the operator insertion point and momentum-projected P and S sinks located on all time slices, the correlation functions G JJ nQ † In can be accessed for all smearing combinations P P , P S, SP , and SS, any temporal separation between the source and the sink t sep = t 1 + t 2 , and any operator separation from the source τ = t 1 . The same propagators are used to calculate P P and P S two-point correlation functions. The spectral representation for Eq. (34) analogous to Eq. (33) is given by where (M I ) m m = m | Q I |m , the ground-state matrix element of interest is M I = (M I ) 00 , and the overlap factors Z J m , Z J m are the same as in Eq. 33. We perform contractions for all combinations of point and smeared sources and sinks in the three-point functions to enhance the analysis of the ground and excited state matrix elements in the next section. To reduce stochastic uncertainties, we also average lattice matrix elements over the spins of the neutron and antineutron states. The specific combinations of (anti)neutron 4-spinor components in the three-point functions that give matrix elements M I are where the signs correspond to the conventions listed in Appendix A. As shown in Sec. II D, these matrix elements are real, and combining them with the conjugated ones is also used to enhance statistics following Eq. (27). Covariance matrix is estimated with optimal shrinkage λ * as described in the main text. Corresponding data points show the effective masses Mn(t) = ln Gnn(t) − ln Gnn(t + 1) with their statistical uncertainties.  (15)(65) TABLE II. Results of two-point function fits from different time ranges: ground-and excited-state energies, reduced χ 2 /N dof , and optimal shrinkage parameters λ * . The uncertainties in individual fits are statistical. The last line shows "fit averages" with statistical and systematic uncertainties computed as described in Appendix B 3.

IV. ANALYSIS OF MATRIX ELEMENTS
To account for excited-state contributions, we perform two-state fits to a truncation of Eq. (34), where A JJ I and B JJ I are products of overlap factors and matrix elements involving only excited states, which are discarded in our calculation. The ground-state overlap factors Z P 0 and Z S 0 required to extract matrix elements of G JJ nQ † In can be obtained independently from fits of two-point functions G P P nn and G P S nn to an analogous two-state model The energies E 0 and E 1 appear in both Eq. (37) and Eq. (38), therefore fits of G nQ † In may be simplified by fixing the state energies E 0 , E 1 to values determined from fits of two-point functions G JJ nn . In principle, the overlaps with excited neutron states Z J 1 are also determined from two-point function fits, thus the number of parameters in Eq. (37) can be reduced by factoring A JJ I , B JJ I into excited-state matrix elements and overlap factors Z J 0,1 , of which only the latter would depend on the neutron interpolating operators. It would be possible if the two-and three-point functions were saturated by contributions only from the ground and the first excited states, or their contributions could be reliably distinguished from higher-energy states omitted from Eqs. (37,38). However, as our two-point function fits in Fig. 2 show, there are higher excited-state contributions to G nn ; in particular, there is large systematic uncertainty on E 1 in (see Tab. II).
These considerations lead us to adopt the following fit strategy: first, a combined fit of G P P nn and G P S nn to Eq. (38) is used to determine the four parameters E 0,1 and Z P,S 0 as summarized in Fig. 2  and B SS I . Also, since P P three-point functions would have even large excited-state contamination and P P threepoint/two-point ratios are not close to their plateau region for the t sep used here (not shown), we do not include G P P nQ † In in our analysis.
With all-mode-averaging described in Sec. III, we obtain one unbiased sample of the two-and three-point functions per gauge field configuration. The number of gauge field configurations used in this calculation N conf = 30 is not large enough to obtain nondegenerate determination of a covariance matrix for the required number of data points 31 ≤ K ≤ 76 included in the three-point correlator fits. Therefore, spin and parity symmetries are used to increase the effective number of unbiased samples of correlation functions. Thus, G nQ †n and G n(Q P I ) †n with two polarizations are treated as four samples per gauge-field configuration, resulting in N = 120 samples for each data point after allmode-averaging bias correction. Polarized two-point functions G nn(±1/2) , Gnn (±1/2) are similarly combined to obtain a statistical ensemble of N = 120 two-point functions. Although this yields an "ensemble" with N > K samples, it is still not sufficient for reliable determination of covariance matrix, which typically requires N > ∼ K 2 . For both two-point and three-point functions, finite sample-size fluctuations may make the sample covariance matrix ill-determined and lead to a numerically unstable inverse covariance matrix required for least-squares fitting. Shrinkage [28,29] has been proposed as a method of improving the condition number of covariance matrix estimates. Denoting the sample covariance matrix by S, the corresponding covariance matrix estimate with shrinkage is given by where 0 ≤ λ ≤ 1 is the shrinkage parameter, and diag(S) is a particular "shrinkage target". Taking any λ > 0 "shrinks" the spectrum of the covariance matrix by reducing the relative size of off-diagonal correlations compared to the diagonal covariance matrix elements. This leads to a better-conditioned covariance matrix and a more robust estimate of the inverse covariance matrix used for χ 2 -minimization. Trivial λ = 0 corresponds to no shrinkage, while λ = 1 removes off-diagonal correlations completely, which is equivalent to an uncorrelated fit. Therefore, varying the parameter 0 ≤ λ ≤ 1 interpolates continuously between correlated (albeit with potentially poorly-determined covariance matrix) and uncorrelated fits. A standard prescription for choosing the optimal shrinkage parameter is to minimize the rms difference between Σ(λ) and the true covariance matrix. A sample estimator for the optimal shrinkage parameter λ * is suggested in Ref. [29] and summarized in Appendix B 1. Bootstrap covariance matrices with optimal shrinkage 6 Σ * = Σ(λ * ) are obtained by inserting λ * from Eq. (B6) into Eq. (39) with S the bootstrap covariance matrix obtained from N boot = 10, 000 samples of two-point and three-point correlation functions. The effects of shrinkage on the central values, uncertainties, and goodness-of-fit of the matrix element fits described below are explored by varying λ, and the results for one choice of fit range are shown in Figs. 3-4. For all operators, the central values and the statistical uncertainties are relatively insensitive to the value of the shrinkage parameter once λ > 0. The χ 2 /N dof values decrease sharply in a small region around λ = 0, however they are much less sensitive for larger λ values. In all cases, the optimal values λ * for the shrinkage parameter are found outside of the region of strong dependence of χ 2 on λ.
The average two-point function and the corresponding bootstrap covariance matrix with optimal shrinkage are used for nonlinear χ 2 -minimization to determine E 0 , E 1 , Z P 0 , and Z S 0 . χ 2 -minimization is reduced to a two-parameter optimization problem by variable projection(VarPro) technique [30,31] detailed in Appendix B 2. In VarPro, the products of overlap factors in Eq. (38) are found from a linear χ 2 -fit for particular values of E 0,1 and the solution is substituted back into χ 2 in order to obtain a two-parameter function χ 2 V P (E 0 , E 1 ), which is then minimized using nonlinear numerical methods. With these four parameters held fixed, the remaining six free parameters in the three-point function fit (37) can also be found from a linear χ 2 fit. The parameter covariance matrix for all 10   Tab. III. An analogous procedure is used to obtain the uncertainties of E 0 and E 1 shown in Tab. II. Results for the matrix elements from the fits that include our smallest t sep value are compared to the ratios of three-point to two-point functions (adjusting for proper overlap factors) in Figs. 5. In addition, the two-point function fits are compared to the corresponding effective masses in Fig. 2. Systematic uncertainties of our analysis procedure are studied by varying the time ranges of data included in the two-state fits. Results of fits of two-point function data G P P nn(nn) (t min P P ≤ t ≤ t max sep ) and G P S nn(nn) (t min P S ≤ t ≤ t max sep ) for a variety of t min P P and t min P S are shown in Tab. II. Results of corresponding fits of three-point function data for a variety of τ S min and τ P min are shown in Tab. III. The data for SP and P S three-point correlation functions are averaged using relation G SP nQ † In (t sep , τ ) = G P S nQ † In (t sep , t sep − τ ) to reduce the number of data points in the fits. Results for bare ground-state matrix elements from different fits are in very good agreement with each other, as shown in Fig.6.
The five fit range choices shown in Tabs. II-III result in acceptable correlated χ 2 /N dof values in fits of two-and three-point function data. These results are combined into final estimates of M I and estimates of their statistical and systematic uncertainties. Since the various fits have different N dof , we use a weighted-averaging procedure defined in Appendix B 3. For a particular fit, the weight is a combination of the likelihood that the fit describes the data (we use its p-value as the likelihood proxy) and its statistical precision, to penalize both fits that fail to describe data and fits that do not constrain the relevant parameters. The same weights are used to determine the average statistical uncertainty, which ensures that including multiple similar fits will not lead to a spurious reduction in the τ min P τ min S t max sep N dof M lat 1 × 10 5 χ 2 /N dof λ * M lat 2 × 10 5 χ 2 /N dof λ * M lat 3 × 10 5 χ 2 /N dof λ * M lat 5 × 10 5 χ 2 /N dof λ * 6 2 13 70 -4. 13 final statistical uncertainty. The weighted mean-square difference between each fit result and the weighted average is used to define the systematic uncertainty due to arbitrarines of choice of a fit window. Applying this weighted averaging procedure to the ground-state energy E 0 of the two-point function yields the result for the nucleon mass that agrees well with the physical value, where we have used the scale-setting result a = 0.1141(3) fm from Ref. [20], which has negligible uncertainty for our purposes as it is much smaller compared to other uncertainties in our calculation. Applying the same procedure to the fit results in Tab. III provides our final estimate of the bare matrix elements including statistical and fitting systematic uncertainties, These lattice regularized matrix elements can be related to renormalized matrix elements through NPR as described in the next section.

V. RENORMALIZATION OF LATTICE OPERATORS
Since the matrix elements of the 6-quark operators are computed on a lattice, they have to be converted to some perturbative scheme, e.g., MS, before they can be used in BSM phenomenology. We calculate conversion factors between lattice-regularized operators and their perturbative definitions nonperturbatively, by computing their Green's functions on a lattice and matching them to perturbative calculations. The operators Q (P) I are the lowest-dimension operators with ∆B = −2, therefore they can only either mix with each other, or get discretization corrections from higher-dimensional operators that vanish in the continuum limit. In the chiral basis, all the 14 operators transform differently under U (2) L ⊗U (2) R flavor symmetry, so they can mix only due to the spontaneous chiral symmetry breaking (SχSB) in QCD, non-perturbative U (1) A violation, or chiral symmetry violations by quark masses and discretization of the fermion action. Mixing due to quark masses and non-perturbative effects should be small if renormalization is carried out in the UV region |p| ≈ µ {Λ QCD , m q } where perturbative matching is applicable. Furthermore, effects of the explicit chiral symmetry violation by the (M)DWF fermion action on a lattice are suppressed as the "residual mass" m res < ∼ m q [20], and thus are also negligible. Therefore, we don't expect that renormalization of our results will be affected by mixing between the chiral-basis operators 7 .

A. RI-MOM amplitudes on a lattice
The lattice renormalization constants for the 6-quark operators are defined as but, as will be shown below, in the chiral-diagonal basis |Z I =J | Z II ≡ Z I , so Q R I (µ) = Z I Q lat I both on a lattice and in continuum perturbation theory. The nonperturbative renormalization and mixing of the six-quark operators is computed using a variant of the RI-MOM scheme [32] with a specific choice of momenta of the external quark states. Since the external states are not color-singlets, the gauge is fixed to the Landau gauge using the Fourier-accelerated conjugate gradient algorithm [33]. All the operators of interest with ∆B = 2 and ∆I = 1 8 can be represented in the generic form where A i = (α i , a i ) are the spin×color indices. Then their Green's functions with external plane-wave quark states, are computed on a lattice contracting six quark propagators computed with a point source at the operator location.
The same propagators are used as for the n-n three-point correlators (see Sec. III), with the only difference that prior to the contraction the propagators are Fourier-transformed at the sink. The six-quark vertex functions are obtained by "amputating" the Green's functions (44) where contraction in {B i } is implied, and the momentum-projected quark propagators are Note that the amputated Green's functions (45) are not symmetric with respect to permutation of the spin×color indices A i , unlike the tree-level vertex function Γ in Eq. (43). This is due to the fact that G I ({p i }) and Λ I ({p i }) depend on the non-equal momenta p i of the external fields. Such dependency would break the isospin symmetry and thus may mix operators from different chiral representations.
One must choose specific momenta for external quark fields in order to preserve the chiral isospin symmetry. The simplest choice p i = p would result in a large momentum p O = 6p at the operator insertion leading to large perturbative corrections in conversion to the MS scheme. To avoid that, the external quark momenta are arranged so that i p i = 0 and, specifically, p i = ±p (see Fig. 7), where p 2 = µ 2 determines the scale for perturbative RI-MOM → MS matching. In addition, the amputated amplitudes (45) must be averaged over permutations of the ±p momenta to enforce the symmetry with respect to the external quark states [18], where the factors are determined by combinatorics. All possible permutations of momenta are implicitly included by Wick contractions, and the symmetries of the color×spin indices are restored automatically. Perturbative matching at the one-loop level for this particular scheme has been computed in Ref. [17]. The lattice renormalization factors Z lat (p 2 ) (42) can be computed by imposing the condition where Z q is the lattice quark field renormalization factor The renormalization factors Z lat IJ can be expressed in terms of the amputated and symmetrized vertex functions Λ {Ai} I (p) that are projected onto the original tree-level structures Γ where the "metric tensor" g JK is diagonal in the chiral basis Q (P) I . (Approximate) chiral symmetry on a lattice is important for ensuring that Z IJ and Λ IJ are also (predominantly) diagonal in this basis. Deviations from the diagonal form are due to the nonzero quark mass and residual chiral symmetry breaking of the DWF discretization. The effect of symmetrization (47) is evident from the magnitude of the off-diagonal components, which is shown in the log scale as the matrix in Fig. 8 comparing the momentum permutation-averaged amplitude (47) to the one with a specific choice of momentum p 1 = p 3 = p 4 = −p 2 = −p 4 = −p 6 = p as in Fig. 7. These data are shown for the momentum p = 2π a 11 48 , 11 48 , 11 48 , 22.5 96 , which is close to a 4d diagonal direction (up to (π/L) along the time axis due to the antiperiodic boundary conditions) and p 2 ≈ (5 GeV) 2 . Therefore, we conclude that in the chiral basis the renormalization matrix Z IJ is diagonal, , which is definitely within our target precision, and the operators Q (P ) I may be renormalized multiplicatively in our lattice calculation. Additionally, we observe that the mixing between 6-quark operators containing different numbers of L, R-diquarks is negligible, indicating that nonperturbative chirality-changing effects due to fluctuations of topology of the QCD vacuum do not lead to mixing in excess of the 10 −3 level.
We define lattice renormalization factors in the RI-MOM scheme for the n-n operators in the chiral basis as Finally, to get rid of the quark field renormalization, we use the renormalization constant Z A for the local axial-vector current A µ =qγ µ γ 5 q. Using the value of Z A computed in Ref. [20], we can compute Z q (p) in the RI-MOM scheme from the condition where A µ qq is the amputated Green's function for the axial current computed analogously to Eq. (45). "Scaleindependent" lattice renormalization factors Z SI Γ = Z lat Γ (p)/Z RI,pert (p) for the vector, tensor, and scalar vertices are shown in Fig. 9.
The value of the lattice renormalization constants Z I (p) may depend on the orientation of the momentum p with respect to the lattice axes due to discretization effects. We compute the lattice vertex functions (51) for various orientations of lattice momenta interpolating between 3d-diagonal and 4d-diagonal orientations to study these effects in the following sections.

B. Perturbative running
In order to convert operator normalization from the RI-MOM scheme discussed above to MS, perturbative matching calculations are required. To extract lattice renormalization factors independent from the momentum subtraction point p, the lattice factors (54) are compared to the perturbative predictions for the RI-MOM scheme in some window p min ≤ |p| ≤ p max where lattice artifacts are believed to be under control. In this section, details of relevant perturbative results are summarized.
The one-loop anomalous dimensions of the operators (7-9) were computed in Ref. [16], and the MS anomalous dimensions to the O(α 2 S ) precision together with Ø(α S ) conversion factors were computed in Ref. [17]. In the chiral basis, the perturbative renormalization of the operators is diagonal (no mixing), and their independent anomalous dimensions are with the coefficients γ (0) I given in Tab. I. These anomalous dimensions are substantially different, which would com-plicate operator renormalization if chiral symmetry was violated by a lattice fermion action and mixing was allowed. We integrate the equations (56) together with an RG equation for the coupling constant α S (µ) using the 4-loop β(α S )-function. Since our lattice QCD action has N f = 2 + 1 dynamical flavors, the lattice factors (54) are matched to Z RI (µ) factors computed in N f = 3 perturbative QCD and the coupling constant α N f =3 S is matched to its physical value at µ ≤ m c . The latter is obtained from a global fit [34] and matched at the m b,c quark mass thresholds. For the reference point µ 0 = 2 GeV, its values 9 are α where Z SI I is a "scale-independent" lattice renormalization factor with a reference point µ 0 defined in the next section. The perturbative scale dependence in both the MS and RI-MOM schemes with N f = 3 and 4 flavors is shown in Fig. 10.

C. Fits of nonperturbative and discretization effects
With known perturbative running, we can separate scale-independent renormalization from lattice artifacts and nonperturbative effects. Correlation functions computed on a lattice are subject to discretization effects that may break rotational symmetry at short distances, which are relevant for the large momenta used in the nonperturbative renormalization. In addition, they may have nonperturbative contributions that complicate matching with perturbative calculations. Below we follow closely the analysis performed in Ref. [35] and extract the scale-invariant renormalization constants Z SI I from a fit where Z SI (µ 0 , a) is the momentum-independent lattice renormalization constant, Z RI,pert I (µ) is the perturbative running of Q I in the RI-MOM scheme and ∆Z disc,NP I encapsulates discretization and nonperturbative corrections. In 9 The coupling constant in the RI-MOM scheme is conveniently defined to be equal to the MS coupling constant. our calculation with O(a 2 ) improved action, the discretization effects must scale as O (ap) 2 , [4] (ap) 4 (59) where we also include the hypercubic invariant ∝ p [4] (see Fig. 11) that breaks the rotational symmetry O(4) → H(4) for n = 0, 1. Although the vertex functions (45) are computed with "exceptional" kinematics p O = 0 (see Fig. 7), they do not have "pole" contributions ∝ 1/p 2 because, unlike the pseudoscalar density operator that can couple to pions, the 6quark operators Q I can couple only to 2-baryon (B = 2) states with masses M ≥ 2m N . However, the nonperturbative contributions are added to Eq. (58) to account for effects of the dimension-2 gluon condensate [36][37][38][39][40][41][42] that may be present in the quark propagators used to amputate the Green's functions. Contributions of condensates to correlation functions are scale-dependent and should be evaluated using OPE as in, e.g., Ref. [43]. Such analysis has not been performed yet, and the correction in Eq. (61) should be regarded as a phenomenological assumption. Another potential source of ∝ 1/p 2 effects are nonperturbative infrared contributions due to potential low-momentum subdiagrams, which may appear due to the same arguments as in Ref. [44]. We perform uncorrelated fit (58) with five parameters (Z SI I , and A, B 1,2 , C) to the lattice data Z lat I (p) for varying sets of momenta p. To keep discretization errors omitted from Eq. (59) as small as possible, we include only momenta p that interpolate between the 3d-and 4d-diagonals, The lowest rotational symmetry-breaking contribution ∝ p [4] /p 4 to Eq. (59) is shown in Fig. 11. Values Z lat I (p) at H(4)-equivalent momenta p are averaged. The fit range p 2 min ≤ p 2 ≤ p 2 max is varied with p min = 1.6, 2.0 GeV and p max = 3.5, 4.0, 4.5 GeV, resulting in 27 ≤ n mom ≤ 61 lattice momentum data points that are distinct with respect to H(4) transformations. We use uncorrelated χ 2 values to evaluate goodness-of-fit and estimate systematic uncertainties from variation of the results with the fit range and the order of the perturbation theory. Although correlated fits would be preferred, we resort to uncorrelated fits, because with a small number of independent configurations N cfg = 30, it is difficult to ensure that covariance matrices of sizes n mom ∼ N cfg are estimated with uniform reliability.
The results of the fits for all fit windows using Z RI,pert from 1-and 2-loop perturbative calculations Z RI,pert are collected in Tab. IV, together with the resulting uncorrelated χ 2 values. In order to obtain the final value, we average the central values over all the fitting methods as described in Sec. B 3. In the last row of Table IV, we show the final conversion coefficients between the lattice bare and MS-renormalized operators Q I that take into account the difference between N f = 3 and N f = 4 QCD perturbative running (see Sec. V B and Eq. (57)).

VI. RESULTS
The four indepdent non-vanishing nn matrix elements in the isospin limit are given in terms of the above bare matrix elements and renormalization factors as where the first uncertainty is the combined statistical uncertainty in M lat I and Z I and the second uncertainty is the combined systematic uncertainty associated with variation in fit window described in Sec. IV, V and Appendix B 3. Quark mass effects lead to negligible systematic uncertainties because of the nearly physical pion mass used [20]. Uncertainties in the determination of the lattice spacing in Ref. [20] are negligible compared to the fitting uncertainties in Eq. (64). Finite-volume effects have been estimated in chiral pertubation theory to be < ∼ 1% effects for the volume used for this study [24]. Discretization effects are expected to be the largest unquantified systematic uncertainty that are neglected in this work. Chiral symmetry leads to O(a) improvement of the fermion action, and discretization effects on meson observables for these configurations have been seen to be percent-level [20]. Discretization effects will be studied and removed from future calculations with multiple lattice spacing.
Final results for the n-n transition matrix elements with statistical and systematic uncertainties added in quadrature and given in Tab. V. These results can be directly compared with MIT bag model results previously used to relate experimental results to BSM couplings [11] as shown in Tab. V. The electroweak-nonsinglet matrix element M 5 is more than an order of magnitude smaller than the electroweak-singlet matrix elements in LQCD. This feature is captured by the MIT bag model, although the sign of M 5 varies between bag model parametrizations. LQCD results The effective Lagrangian for n-n oscillations given in Eq. (30) can be used to parameterize the n-n vacuum transition rate for a generic BSM theory as where η = v 2 /Λ 2 BSM is the ratio of the Higgs v.e.v. and the BSM scale squared. Both the matrix elements M and the Wilson coefficients C (P) are scheme-and scale-dependent, and these dependencies must cancel in τ n-n . Below we present results with coefficients C defined in MS scheme. The Wilson coefficients in Eq. (65) are predicted to be non-zero in various BSM theories, see Refs. [45][46][47] for reviews and further refences, and are calculable at tree-level in QCD at BSM scales µ = Λ BSM . The n-n vacuum transition rate is given in terms of the above results by To make the prefactor dimensionless, we use the "reference" normalization scale of 700 TeV. Estimates based on Eq. (66) put BSM theories with scales of Λ BSM ∼ 700 TeV and O(1) matching coefficients within reach of nextgeneration experiments that will be able to detect baryon number violation with τ −1 n-n ≥ 10 9 s [48][49][50][51]. To more precisely assess the expected signatures of theories with B-violation at Λ BSM ∼ 700 TeV, the operators can be evolved to µ = Λ BSM using the results of Refs. [16,17], The n-n transition rate can be expressed in terms of the matrix elements at this scale as This result can be combined with tree-level BSM matching results for C MS I (700 TeV) to extract constraints on BSM theory parameters from experimental constraints on n-n oscillations.

VII. CONCLUSION
We have performed the first lattice QCD calculation of the renormalized neutron-antineutron transition matrix elements needed to extract BSM physics constraints from n-n oscillation experiments. The precision of our final results including statistical and most systematic uncertainties is 15 − 30% for the electroweak-singlet matrix elements M 1 , M 2 , and M 3 , which can be straightforwardly improved in future calculations. Several important sources of systematic uncertainty are under control for the first time, most importantly non-perturbative renormalization, chiral symmetry violations, excited state contamination, and quark mass dependence. The two sources of systematic uncertainty that are not completely controlled in this pioneering calculation are finite volume and discretization effects. To summarize our control of common systematic uncertainties in lattice calculations: • The (nearly exact) physical pion mass m π = 139.2(4)MeV in our calculation eliminates the need for chiral extrapolation, which would otherwise introduce systematic uncertainties associated with low-energy effective theory. In addition, the large difference of our results from the MIT bag model may have a similar origin as the strong suppression of proton decay matrix elements found in the chiral bag model [52], therefore using the realistic light quark masses in our calculation is arguably the most important systematic effect we have under control.
• The chirally symmetric Möbius domain wall fermion action used to generate these gauge field ensembles by the RBC/UKQCD collaborations [20] and compute neutron-antineutron matrix elements in this work ensures that the 14 distinct |∆B| = 2 operators do not mix with each other and renormalization and conversion of lattice operators to MS scheme is free from associated uncertainties. In particular, the nonperturbatively computed operator mixing matrix in RI-MOM scheme is diagonal up to O(10 −3 ) corrections, which are two orders of magnitude below other uncertainties and can be safely neglected. The identical action is used for valence quarks, so this is a fully unitary calculation.
• Excited-state effects are accounted for using correlated two-state fits with 10 different values of t sep and different combinations of nucleon source and sink smearing. The energy gaps are extracted from correlated fits to nucleon two-point functions. Since we have limited statistics for such a large number of τ, t sep points included in correlated fits, we use "shrinkage" estimators to obtain well-conditioned covariance matrices. We obtain systematic errors by varying the fit ranges and averaging their results weighted by the quality-of-fit figure.
• Renormalization effects are included through NPR in an RI-MOM scheme as described in Sec. V and one-loop matching to MS using the results of Ref. [17]. Some discretization effects in NPR results such as rotational symmetry breaking and (ap) 2 dependence are studied and removed by fitting the lattice data with different quark momentum scales and orientations, varying scale ranges, and comparing to 1-and 2-loop perturbative QCD running.
• Although we do not control finite-volume effects directly in this study on a single ensemble, we expect them to be small. First, finite-volume effects are suppressed with e −mπL where m π L ∼ 3.9 for the volume used for this study, which is generally considered sufficiently large for nucleon structure calculations [53]. Second, chiral perturbation theory calculations in Ref. [24] estimate that finite-volume effects lead to corrections below 1% to M I for the volume used in this study. Future lattice calculations at additional volumes could be used to test this prediction and perform an infinite-volume extrapolation.
• Discretization effects are the least-controlled systematic uncertainty in our current work. Lattice QCD calculations with finer lattice spacing(s) in the immediate future will be used to fully quantify and remove discretization effects that are not controlled in this calculation. However, it is reasonable to assume that discretization effects are small compared to our current combined uncertainty from other sources. First, the chirally-symmetric fermion action that we use is automatically O(a 2 )-improved. Second, the meson decay constants computed on this ensemble (before finite volume and discretization corrections are applied) are within 0.6% of the physical values (f π = 131.1(4), f K = 156.4(4) GeV [20] compared to PDG values f π = 130.4(2), f K = 156.2(7) GeV [54]). Finally, the nucleon effective mass and energy dependence on the momentum is in close agreement with the continuum limit [55].
Our renormalized lattice QCD results for n-n transition matrix elements provide a significant step forward in accuracy and reliability compared to previous results from quark models and preliminary lattice studies. The matrix elements predicted by QCD are found to be 4-8 times larger than the predictions of the MIT bag model for the dominant electroweak-singlet operators. This difference between our lattice results and previously available bag model results is much larger than the statistical or systematic uncertainties present in this calculation and is also much larger than the expected size of finite-volume effects that have not yet been studied directly. There is less certainty about the size of discretization artifacts; however, the automatic O(a 2 ) improvement due to the chiral symmetry as well as minuscule discretization corrections in the meson decay constants, nucleon mass and dispersion relation make large discretization effects in the n-n matrix elements very unlikely.
The difference in M I between the bag model and our lattice results leads to increased experimental sensitivity to baryon-number violating interactions that may cause n-n oscillations. Numbers of events that can be observed both in quasi-free neutron oscillation experiments and underground nuclear decay experiments are proportional to τ −2 n-n ∝ |M I | 2 , therefore the ×(4 . . . 8) larger values of the n-n matrix elements found in our work lead to ×(16 . . . 64) increase in the event rates. Since our results are obtained from ab initio QCD calculations in a model-independent way, they must be used for more precise assessments of the potential of planned n-n oscillation searches as well as stronger constraints on theories of baryon-number violation and baryogenesis in the future.

(±)
T (x) T . (A12) which are used to construct (anti)neutron states on a lattice. This construction is more natural in the standard (Dirac-Pauli) basis in which the γ 4 matrix is diagonal. It is related to the de Grand-Rossi basis commonly used in lattice calculations by the transformation The operators (A11) create the neutron and antineutron states with definiteẑ-spin as These states are used to determine the properties of the n-n matrix elements in Sec.II D and define them in terms of three-point functions in Sec. III.

Shrinkage estimation of covariance matrices
Correlated χ 2 -fits require sample covariance matrices that are difficult to estimate when the number of data samples N is limited compared to the number of data points K, as in our case. In order to estimate covariance matrices that can be safely inverted, we use the "optimal shrinkage estimator" described in Refs. [29]. Shrinkage involves replacing the covariance matrix with a linear combination of a well-conditioned "shrinkage target" and the original covariance matrix. It has been shown that expectation values of "shrunk" covariance matrices are closer to the true covariance matrix than the sample covariance matrix [28]. The condition number of the covariance matrix is also improved by shrinkage and estimates of χ 2 relying on the inverse covariance matrix are more robust. Shrinkage targets that better approximate the true covariance matrix naturally lead to better estimates of the true covariance matrix from a finite sample, but any prescription for defining the "shrinkage parameter" introduced below that leads to zero shrinkage in the infinite statistics limit will provide a consistent estimator for the true covariance matrix.
The estimator in Ref. [29] uses a shrinkage target proportional to the K ×K identity matrix I where K is the number of data points. However, correlation functions in lattice calculations vary over orders of magnitude if a wide range of t sep are used for fitting. To transform the covariance matrix into a form where the shrinkage target of Ref. [29] more closely resembles the true covariance matrix, we normalize the data by subtracting the mean and diving by the square root of the variance. For data points x i α where i = 1, · · · , N labels decorrelated statistical samples and α = 1, · · · , K labels data points (i.e. t in two-point function fits and τ , t sep in three-point function fits), define normalized data points y i α and a normalized sample correlation matrix ρ αβ as where the sample mean and covariance are defined as The correlation matrix with optimal shrinkage is given by where µ = 1 K Tr[ρ] = 1 is the mean of the spectrum of ρ and the optimal shrinkage parameter λ * is defined to minimize the expected Frobenius norm ||X|| = Tr[X X T ] of the difference E min λ ||ρ(λ) − || 2 between the estimator ρ * and the true correlation matrix . A sample estimator for the optimal shrinkage parameter is given in Ref. [29] λ * = min{b 2 , d 2 } d 2 , (B4) The quantity d 2 estimates the dispersion of the eigenvalues of the sample correlation matrix ρ, which typically has a wider spectrum and correspondingly larger (worse) condition number compared to the true correlation matrix . The optimal estimator (B3) "shrinks" the spectrum by emphasizing the diagonal elements and makes the matrix Σ * better-conditioned, resulting in more statistically stable χ 2 values in correlated fits. Multiplying both sides of Eq. (B3) by the normalization factor in Eq. (B1) yields the corresponding estimator for the covariance matrix Σ * αβ = S αα S ββ ρ * αβ = λ * S αα S ββ δ αβ + (1 − λ * )S αβ , Σ * = λ * S + (1 − λ * )diag(S) . (B7) This shrinkage prescription is therefore equivalent to an interpolation between a fully correlated fit with λ * = 0 (no shrinkage), and an uncorrelated fit with λ * = 1 (full shrinkage). Although this prescription does not provide the strictly optimal λ * minimizing the distance between Σ * and Σ, it gives a simple practical prescription for a stable and consistent choice of the shrinkage parameter. Optimal closeness between ρ * and suggests that Σ * should provide an acceptable approximation of Σ that is better-conditioned than S.