Light-cone distribution amplitudes of the baryon octet

We present results of the first ab initio lattice QCD calculation of the normalization constants and first moments of the leading twist distribution amplitudes of the full baryon octet, corresponding to the small transverse distance limit of the associated S-wave light-cone wave functions. The P-wave (higher twist) normalization constants are evaluated as well. The calculation is done using $N_f=2+1$ flavors of dynamical (clover) fermions on lattices of different volumes and pion masses down to 222 MeV. Significant SU(3) flavor symmetry violation effects in the shape of the distribution amplitudes are observed.


Introduction
Understanding the structure of matter in terms of quark and gluon degrees of freedom is the ultimate goal of the physics of strong interactions. For this purpose, the intuitive quantummechanical representation of a hadron as a superposition of Fock states with different numbers of partons in the infinite momentum frame (or using light-cone quantization) is very useful in order to develop the underlying physics picture. It also provides a good basis for theoretical modeling. Although a priori there is no reason to expect that, e.g., the nucleon wave function components with, say, 100 partons (quarks and gluons) are suppressed relative to those with only three valence quarks, the phenomenological success of quark models suggests that in many cases only the first few Fock components are important. Also the analysis of hard exclusive reactions involving large momentum transfer from the initial to the final state baryon within QCD perturbation theory [1][2][3] suggests that such processes are dominated by the overlap of the valence light-cone wave functions at small transverse separations, usually referred to as light-cone distribution amplitudes (DAs).
The DAs can be viewed as light-cone wave functions integrated over the quark transverse momenta [1]. They are fundamental nonperturbative functions that are complementary to conventional parton distributions, but are more elusive because their relation to experimental observables is less direct compared to quark parton densities. DAs are scaledependent and for asymptotically large scales they are given by simple expressions, the so-called asymptotic DAs [1,3]. There are many indications, however, that these asymptotic expressions poorly approximate the real DAs for the range of momentum transfers accessible in present experiments.
The theoretical description of DAs is based on the relation of their moments, i.e., integrals over DAs weighted by powers of momentum fractions, to matrix elements of local operators. Such matrix elements can be estimated using nonperturbative techniques and the DAs can be reconstructed as an expansion in a suitable basis of polynomials in the momentum fractions.
First estimates of the first and the second moments of the baryon DAs have been obtained more than 30 years ago using QCD sum rules (QCDSR) [4][5][6][7]. These early calculations suggested large deviations from the corresponding asymptotic values and were used extensively for model building of DAs [4][5][6][7][8] that allowed for a reasonable description of the experimental data available at that time within a purely perturbative framework, see the review [3].
Despite its phenomenological success, this approach has remained controversial over the years. In particular the QCD sum rules used to calculate the moments have been criticized as unreliable, see, e.g., ref. [9]. Also, nowadays it is commonly accepted that perturbative contributions to hard exclusive reactions at accessible energy scales must be complemented by the so-called soft or end-point corrections. Estimates of the soft contributions using QCD sum rules [10], quark models [11] and light-cone sum rules [12][13][14] favor nucleon DAs that deviate only mildly from the asymptotic expressions.
It has now become possible to calculate moments of the DAs from first principles using lattice QCD. The first quantitative results for the nucleon have been obtained by the QCDSF collaboration [15,16] using two flavors of dynamical (clover) fermions, followed by [17], where a much larger set of lattices was used including ensembles at smaller pion masses, close to the physical point. The latter paper also contained an exploratory study of the DAs of negative parity nucleon resonances, see also ref. [18].
In this work we extend the analysis of ref. [17] to the full J P = 1 2 + baryon octet.
In addition to theoretical completeness, our study is motivated by applications to weak decays of heavy baryons such as Λ b and Λ c . Such baryons are produced copiously at the LHC. As more data are collected, studies of rare b-baryon decays involving flavor-changing neutral currents offer interesting insights into the quark mixing matrix and, potentially, may reveal new physics beyond the Standard Model. In particular the Λ b → Λµ + µ − decays are receiving a lot of attention, see, e.g., refs. [19,20] and references therein. B-meson decays into baryon-antibaryon pairs are also interesting. The pattern of SU(3) flavor symmetry breaking in weak decays is known in general to be nontrivial. In particular the large asymmetry observed in the decay Σ + → pγ has been fueling a lot of discussions over many years and remains poorly understood at the parton level (see, e.g., ref. [21]). Another motivation comes from the emerging possibilities to study the transition form factors for the electroexcitation of nucleon resonances at large photon virtualities, planned for the JLAB 12 GeV upgrade [22], with the hope that similar transition form factors for hyperon production, e.g., in large-angle πN scattering, will also become accessible in the future. This perspective already stimulated several theory studies, see, e.g., refs. [23,24].
In this first study we will mainly address the development of the necessary formalism and methodical issues. Studies of hyperon DAs have a long history [6], however, we found that the definitions existing in the literature are not very convenient to study effects of SU (3) breaking and that the standard notation is, in part, contradictory. Therefore, we explain our notation and provide the necessary definitions in the introductory section 2. The physical interpretation of the DAs in terms of light-cone wave functions is considered as well. The related appendix A explains the phase conventions for the flavor wave functions used in this work.
Section 3 is devoted to the lattice formulation of the problem at hand, the definition of correlation functions used in our analysis of the couplings and the first moments of baryon octet DAs, and the strategy to approach the physical limit of small pion mass. Calculations in this work are performed on a set of ensembles provided by the coordinated lattice simulations (CLS) effort [25]. These are obtained using the tree-level Symanzik improved gauge action and 2+1 dynamical Wilson (clover) quark flavors. In our calculation we start at the flavor symmetric point, where all quark masses are equal, and approach the real world in such a way that u and d quark masses decrease and simultaneously the s quark mass increases so that the average mass is kept (approximately) constant [25,26]. Section 4, complemented by appendices B and C, explains our renormalization procedure. We employ a nonperturbative method based on the well-known RI ′ SMOM scheme, combined with matching factors calculated in continuum perturbation theory to convert our results to the MS scheme. The renormalization of flavor-octet operators turns out to be more complicated than the nucleon case and is discussed in some detail.
Section 5 contains a discussion of chiral extrapolation and SU(3) symmetry breaking in the framework of three-flavor baryon chiral perturbation theory (BChPT). Our presentation is based on the recent analysis in ref. [27]. One result is a simple relation between the DAs of the Σ and Ξ hyperons which has the same theory status as the famous Gell-Mann-Okubo sum rule for baryon masses and is satisfied to high accuracy for our lattice data.
In section 6 our final results are presented and compared with the existing lattice (for the nucleon) and QCD sum rule calculations. We find that deviations of the baryon DAs from their asymptotic form at hadronic scales are small, up to an order of magnitude smaller than in old QCD sum rule calculations. The SU(3) breaking in the corresponding shape parameters is, on the contrary, much larger than anticipated. Section 7 is reserved for a summary and conclusions.
It has to be said that while our calculation provides the first qualitative insight into SU(3) breaking of octet baryon DAs from lattice QCD, it has not yet reached a quantitatively mature state. This is mainly due to the lack of a continuum extrapolation, which can have a significant impact on DAs (cf. ref. [17]). Actually, the whole CLS strategy to simulate with open boundary conditions is motivated by the fact that, at presently used lattice constants, discretization errors are significant. Moreover, on coarse lattices the lattice spacing will depend on the observable employed for scale setting. As our study is an exploratory one with significant systematic uncertainties, this fact is irrelevant in the present context. We use ensembles with the lattice spacing a = 0.0857 (15) fm, which is determined from the Wilson flow method as described in ref. [25]. The dimensionless flow time t 0 a 2 was extrapolated to the physical point and the lattice spacing was then assigned by matching to the continuum limit value √ t 0 = 0.1465(21)(13) fm determined in ref. [28]. In the future we intend to include finer lattices, which are currently being generated within the CLS effort. This will then allow us to take the continuum limit and also eliminate scale setting ambiguities related to the nonzero lattice spacing.

Baryon distribution amplitudes
Baryon DAs [1][2][3] are defined as matrix elements of renormalized three-quark operators at light-like separations: where B(p, λ)⟩ is the baryon state with momentum p and helicity λ, while α, β, γ are Dirac indices, n is a light-cone vector (n 2 = 0), the a i are real numbers, µ is the renormalization scale and f, g, h are quark fields of the given flavor, chosen to match the valence quark content of the baryon B. The Wilson lines, which are needed for gauge invariance, as well as the color antisymmetrization, which is needed to form a color singlet, are not written out explicitly but always implied. Renormalization of three-quark operators using dimensional regularization and minimal subtraction requires some care; we will be using the renormalization scheme proposed in ref. [29]. Restricting ourselves to the analysis of the lowest 1 2 + multiplet, neglecting electromagnetic interactions and assuming exact isospin symmetry (m l ≡ m u = m d ), it is sufficient to consider four cases: For definiteness we choose the following flavor ordering: Here Γ DA andΓ DA are the Dirac structures corresponding to the distribution amplitude DA B (x i ), see eq. (2.9) of ref. [30], and u B (p, λ) is the Dirac spinor with on-shell momentum p (p 2 = m 2 B ) and helicity λ. This decomposition can be organized in such a way that all DAs have definite collinear twist. The scale dependence will be suppressed from now on, unless it is explicitly needed. The variables x 1 , x 2 , x 3 are the momentum fractions carried by the quarks f, g, h, respectively, and the integration measure is defined as (2.5) The factor δ(1 − x 1 − x 2 − x 3 ) enforces momentum conservation.

Leading twist distribution amplitudes
In this work we will mainly be concerned with the DAs of leading twist three. To this accuracy the general decomposition in eq. (2.4) is simplified to three terms [3]: with the charge conjugation matrix C and the notatioñ Our DAs V N , A N and T N correspond to V 1 , A 1 and T 1 in ref. [30].
The equivalent definition in terms of the right-and left-handed components of the quark fields, q ↑(↓) = 1 2 (1 ± γ 5 )q, is sometimes more convenient: In the nucleon case the combination [V − A] N appearing in the first of these equations is traditionally referred to as the leading twist nucleon DA Φ N . For the full octet we define If Φ B is given, the V B and A B components can be reconstructed due to their different symmetry properties under the exchange of the first and the second quark: Using isospin symmetry one can further show for the nucleon To fully exploit the benefits of SU(3) flavor symmetry it proves convenient to define the following set of DAs: where for the nucleon Π N = Φ N + up to isospin breaking effects. In the limit of SU(3) flavor symmetry, where m u = m d = m s (and in particular at the flavor symmetric point with physical average quark mass indicated by ⋆), the following relations hold: 1 Therefore, the amplitudes Π B (or T B ) only need to be considered when flavor symmetry is broken. In the flavor symmetric limit Φ ⋆ + and Φ ⋆ − can again be combined to form a single distribution amplitude Φ ⋆ = Φ ⋆ + + Φ ⋆ − . One can show (see the detailed discussion in section 5) that to first order in the symmetry breaking parameter, m 2 K − m 2 π ∝ m s − m l , the following relation holds: To understand the physical meaning of the DAs it is instructive to work out their relation to light-front wave functions. The leading twist approximation corresponds to taking into account S-wave contributions in which case the helicities of the quarks sum up to the helicity of the baryon (cf. refs. [31,32]). Suppressing the transverse momentum dependence one finds and (2.18) DAs can be expanded in a set of orthogonal polynomials (conformal partial wave expansion) in such a way that the coefficients have autonomous scale dependence at one loop. The first few polynomials are (see, e.g., ref. [33]) Note that all P nk have definite symmetry (being symmetric or antisymmetric) under the exchange of x 1 and x 3 . Taking into account the corresponding symmetry of the DAs, defined in eq. (2.13), a generic expansion reads In this way all nonperturbative information is encoded in the set of (scale-dependent) coefficients ϕ B nk , π B nk , which can be related to matrix elements of local operators. In each DA only polynomials of one type, either symmetric or antisymmetric under exchange of x 1 and x 3 , appear.
The leading contributions 120x 1 x 2 x 3 ϕ B 00 and 120x 1 x 2 x 3 π B≠Λ 00 are usually referred to as the asymptotic DAs. The corresponding normalization coefficients ϕ B 00 and π B≠Λ 00 can be thought of as the wave functions at the origin (in position space). In what follows we will use the notation Note that for the nucleon the two couplings coincide, f N T = f N . For the Λ baryon the zeroth moment of T Λ vanishes. The higher-order coefficients are usually referred to as shape parameters. Note that, in contrast to ref. [17], we do not separate the couplings f B and f B≠Λ T as overall normalization factors, so that our ϕ N nk correspond to f N ϕ N nk of [17]. The one-loop scale evolution of the couplings and shape parameters is given by where β 0 = 11−2N f 3 is the first coefficient of the QCD β-function. In this work we restrict ourselves to the contributions of first order polynomials P 10 , P 11 and omit all higher terms. The relevant one-loop anomalous dimensions are The scale dependence of f B and f B≠Λ T is identical and is known up to three-loop order, see refs. [29,34].

Higher twist contributions
The general decomposition (2.4) contains 21 DAs of higher twist, which altogether involve only up to three new normalization constants (just two for N , Σ and Ξ), for details see refs. [27,30]. They can be defined as matrix elements of local three-quark twist four operators without derivatives. These twist four couplings are also interesting in a broader context, e.g., in studies of baryon decays in generic GUT models [35], and as input parameters for QCD sum rule calculations, see, e.g., refs. [24,36,37].
We use the following definitions: for the isospin-nonsinglet baryons (N , Σ, Ξ) and for the Λ baryon. The definitions are chosen such that at the flavor symmetric point cf. ref. [27]. For the nucleon the definitions in terms of chiral fields in eq. (2.24) are equivalent to the traditional definitions of λ N 1 and λ N 2 not involving chiral projections, as used in ref. [17]. Analogous definitions can also be given for the Λ baryon: The one-loop evolution for all twist four normalization constants is the same: (2.28) The corresponding anomalous dimensions are known up to three-loop accuracy [29,34]. The scale dependence of the couplings λ B 1 and λ Λ T is the same to all orders, whereas for λ B 2 it differs starting from the second loop.

Lattice formulation
In Euclidean spacetime a direct calculation of DAs is not possible, since this would require quark fields at light-like separations. However, lattice QCD allows us to access moments of the DAs, e.g., and similarly for the other functions. They are related to matrix elements of local threequark operators, whose general form reads Here we use a multi-index notation for the covariant derivatives, Dl ≡ D λ 1 ⋯D λ l . The Dirac structures that we consider, Γ Xr r andΓ Xr , are listed in table 1. 2 As sources for the baryon fields we have used the interpolating currents with an optimized number of smearing steps in the quark sources to suppress excited state contributions. The other baryons can then be obtained by means of isospin symmetry.

Correlation functions
Moments of baryon DAs can be extracted from the ground state contribution to the twopoint correlation functions. Neglecting the exponentially suppressed excited states the correlation functions can be written as with the energy E B = E B (p) = m 2 B + p 2 , where we assume the continuum dispersion relation. The momentum-dependent coupling Z B = Z B (p) describes the overlap between the smeared source operator and the physical baryon ground state and can be obtained from the correlator where γ + = (1 + kγ 4 ) 2 with k = m B * E B * suppresses the negative parity contribution [17,38]. 3

Leading twist -zeroth moments
In order to extract the leading twist normalization constants, the following linear combinations of operators are constructed such that their matrix elements do not contain any contributions of higher twist: where X can be V, A or T . The leading twist baryon couplings can be determined from the following correlation functions: where c V = c A = 1 and c T = −2. Again, X can be V, A or T . In practice we only consider the zero momentum correlators C B,000 X ,B and C B,000 X ,C as they are less noisy and, therefore, can be measured with higher accuracy. The couplings of interest are related to the calculated zeroth moments as follows: 8) where f N T = f N due to isospin symmetry. The remaining zeroth moments of the leading twist DAs V B , A B and T B vanish:

Leading twist -first moments
First moments of DAs can be calculated utilizing operators containing one covariant derivative. For l + m + n = 1 we define the leading twist combinations where the braces indicate symmetrization. For the calculation of the first moments of the leading twist DAs one can use the correlation functions (l + m + n = 1) One immediately notices that at least one nonzero component of spatial momentum is required to extract the first moments. We evaluate C B,lmn X ,A,1 , C B,lmn X ,B,1 and C B,lmn X ,C,1 with momentum in x direction (p = (±1, 0, 0)), 4 and C B,lmn X ,A,2 , C B,lmn X ,B,2 and C B,lmn X ,C,2 with momentum in y direction (p = (0, ±1, 0)). For momentum in z direction (p = (0, 0, ±1)) only the correlator C B,lmn X ,B,3 can be used. We do not consider the remaining two correlators as they require a higher number of nonvanishing momentum components, which would lead to larger statistical uncertainties. The shape parameters defined in eq. (2.20) can be expressed as linear combinations of V B lmn , A B lmn and T B lmn via eq. (2.13). For the N , Σ and Ξ baryons, where π N 11 = ϕ N 11 due to isospin symmetry. For the Λ baryon, In addition we define combinations corresponding to the sum of contributions with the derivative acting on each of the three quarks: where π N 00,(1) = ϕ N 00,(1) due to isospin symmetry. Thanks to the Leibniz product rule for derivatives, this sum can be written as a total derivative acting on a local three-quark operator without derivatives so that in the continuum corresponding to the momentum conservation condition x 1 + x 2 + x 3 = 1, see eq. (2.5). However, the Leibniz rule is violated by lattice discretization and this relation can only be expected to hold after continuum extrapolation in a renormalization scheme which respects the Lorentz symmetry. Note that under renormalization ϕ B 00,(1) and π B≠Λ 00,(1) mix with the other first moments, see section 4. It turns out that for the bare lattice values the equalities (3.15) are violated significantly. After renormalization and conversion to the MS scheme we find that the sum rules (3.15) are fulfilled to an accuracy between ≈ 96% and ≈ 98% for our value of the lattice spacing a ≈ 0.0857 fm, see tables 4 and 5. A violation of the momentum sum rule of this size is in perfect agreement with the results in ref. [17], where similar discretization effects have been observed.

Higher twist
Higher twist normalization constants can be calculated from the correlation functions Table 2. List of the ensembles used in this work, labeled by their CLS identifier. The pion and kaon masses have been obtained from two-point functions. β = 3.4 corresponds to the lattice spacing a ≈ 0.0857 fm. An in-depth description of these lattices can be found in ref. [25].

16)
Due to symmetry properties of the associated operators it follows that 18) and the corresponding correlators vanish.

Details and strategy of the lattice simulation
In this analysis we use lattice ensembles generated within the coordinated lattice simulations (CLS) effort. These N f = 2 + 1 simulations employ the nonperturbatively order a improved Wilson (clover) quark action and the tree-level Symanzik improved gauge action. We used a modified version of the Chroma software system [39,40], the LibHadronAnalysis library and efficient inverters [41][42][43]. To enhance the ground state overlap the source interpolators are Wuppertal-smeared [44], employing spatially APE-smeared [45] transporters. A special feature of CLS configurations is the use of open boundary conditions in time direction [43,46]. This will eventually allow for simulations at very fine lattices without topological freezing. We achieve an efficient and stable hybrid Monte Carlo (HMC) sampling by applying twisted-mass determinant reweighting [43], which avoids near-zero modes of the Wilson Dirac operator. A list of the CLS ensembles used in this work is given in table 2. As schematically represented in figure 1, these ensembles are tuned such that the average quark mass reproduces (approximately) the physical value. They have rather large spatial volumes (L > 2.7 fm, with m π L ≳ 4) and high statistics. Consecutive gauge field configurations are separated by four molecular dynamics units.
Lattice calculations with the average quark mass fixed at the physical value have already been carried out for hadron masses and some form factors [26,47,48].  Figure 1. Plot showing the meson masses of the lattice ensembles used in this study. All quantities are made dimensionless using the average octet baryon mass X b , cf. section 5. Along the flavor symmetric line (blue) all pseudoscalar mesons have equal mass (m 2 K = m 2 π ), which is equivalent to equal quark masses (m l = m s ). The (green) line of physical normalized average quadratic meson mass ((2m 2 K + m 2 π ) X 2 b = phys.) corresponds to an approximately physical mean quark mass (2m l + m s ≈ phys.). The red line is defined by (2m 2 K − m 2 π ) X 2 b = phys. and indicates an approximately physical strange quark mass (m s ≈ phys.). The red dot marks the physical point. is then approached in such a way that u and d quark masses decrease and simultaneously the s quark mass increases so that their average is kept constant.
For each configuration we have carried out all measurements with 3 different source positions t src = 30a, 47a and 65a. Taking the average of correlators from all these different sources is not advisable as the open boundary conditions break translational invariance in time. Instead, we average suitable forward and backward propagating states, i.e., the forward direction from t src = 30a and the backward direction from t src = 65a as well as the forward and the backward running state from t src = 47a. The two remaining ones (backward from t src = 30a and forward from t src = 65a) are not considered in this analysis, as sink positions closer than ∼ 20 time slices to the boundary can show significant boundary effects due to the open boundary conditions. The second step of the data analysis is conducted by averaging over appropriate twopoint functions and momenta as outlined in section 3.1. For the statistical analysis we then N (47,47) N (30,65) Λ (47,47) Λ (30,65) Σ (47,47) Σ (30,65) Ξ (47,47) Ξ ( Figure 2. The data points in this plot show the effective baryon masses obtained from the two forward-backward averaged smearedsmeared correlation functions (as described in the main text) calculated on the C101 ensemble with zero three-momentum. The plateaus start at t − t src = 7a, where excited states are sufficiently suppressed. For each baryon the horizontal line represents the result of a combined fit to both correlators in the range 7a ≤ t − t src ≤ 20a. . Typical plot (from the C101 ensemble) used for the determination of the fit ranges by varying the value of the minimal source-sink distance t start . It shows the leading twist normalization constants obtained from the correlators given in eq. (3.7). A conservative choice is t start = 9a, where the results have fully saturated for all leading twist couplings. A variation of the maximal source-sink separation within reasonable bounds did not have any significant impact on the result. Here it is always set to t end = 20a. generate 1000 bootstrap samples per ensemble using a binsize of 8 to eliminate autocorrelations. For each sample we use a χ 2 -measure to simultaneously fit the two correlation functions resulting from the forward-backward averaging procedure described above.
In order to exclude contributions from excited states the choice for the lower bound of the fit range is crucial. Figure 2 demonstrates, that, with increasing source-sink distance, the excited states decay and clear plateaus in the effective masses emerge. To determine the optimum minimal source-sink distance t start we perform multiple fits with varying fit ranges for all observables. t start is chosen in such a way that fits with even larger starting times no longer show any systematic trend in the fit results. As an example, figure 3 shows the fitted leading twist coupling constants as a function of t start .

Renormalization
Bare lattice results have to be renormalized. The preferred scheme in phenomenological applications is based on dimensional regularization where, for baryons, there are subtleties due to contributions of evanescent operators that have to be taken into account, see refs. [29,34]. For simplicity, we refer to the prescription suggested in [29] as the MS scheme. In principle, lattice perturbation theory could be used to compute the required renormalization coefficients, however, such calculations are nontrivial and often poorly convergent. Therefore, we employ a nonperturbative method combined with matching factors Table 3. List of three-quark operator multiplets transforming irreducibly under H(4), sorted by operator dimension and representation. For zero derivatives all operators are listed. For the operators with derivatives only the leading twist multiplets are shown. The dots indicate the position of the remaining higher twist operators. The nomenclature follows ref. [49]. no derivatives 1 derivative 2 derivatives dimension 9/2 dimension 11/2 dimension 13/2 calculated in continuum perturbation theory to convert our results to the MS scheme in the end. The details of our renormalization procedure are described in appendix C.
Studying the renormalization of our three-quark operators, we face the problem of the reduced symmetry of the four-dimensional lattice relative to the continuum. The lattice symmetry group for fermions is known as the spinorial hypercubic group H(4), which has five irreducible spinorial representations: τ Operators of higher dimension (e.g., with derivatives) can in general mix with operators of lower dimension transforming according to the same representation. This mixing is highly undesirable as the admixture of the lower dimensional operators is proportional to negative powers of the lattice spacing a and leads to severe problems when taking the continuum limit. It can be avoided if one chooses operators from a H(4) representation where no lower dimensional multiplets exist. For single-derivative operators this happens in the case of the representation τ  (see table 3). Therefore, we use O D5−D7 (and do not use O D2−D4 ) to calculate the leading twist first moments. There exists a fourth multiplet, O D8 , in this representation which can in principle mix with these operators. However, these operators have different chirality. Hence, an admixture is a power-suppressed O(a) effect. Furthermore, in the continuum limit all octet-baryon-to-vacuum matrix elements of operators within O D8 vanish identically, even though the operators themselves are nonzero. We have verified this property numerically on the lattice. In summary, the admixture of O D8 to O D5−D7 is completely negligible and can safely be ignored in our analysis.
The classification of three-quark operators according to irreducible representations of H(4) and their behavior under renormalization is independent of the flavor structure. Indeed, in ref. [49] the irreducible multiplets have been given assuming three generic flavors. However, the mixing is further restricted by the behavior of the operators under permutations of the three quarks. Consider operator multiplets that transform irreducibly also with respect to the permutation group S 3 . Such operators are given in appendix B.2. In a flavor symmetric world the three inequivalent irreducible representations of S 3 would correspond to a decuplet of SU(3) (trivial representation of S 3 , labeled D), an SU(3) singlet (totally antisymmetric representation, labeled S ) and a doublet of octets (two-dimensional representation of S 3 , labeled O). Of course, in the real world flavor symmetry is broken, which can lead to the appearance of renormalization constants from different representations of S 3 in a single matrix element.
To be more specific, we obtain the renormalized leading twist couplings from operators belonging to the representation τ 12 1 of H(4). These can be arranged into a doublet of octet multiplets (with a common renormalization factor Z Of ) and a decuplet multiplet (with renormalization factor Z Df ), so that we end up with If f B≠Λ T = f B (as is the case for the nucleon and for the SU(3) symmetric limit), the first equation reduces to a multiplicative renormalization with one and the same factor Z Of .
As detailed above, for the higher twist normalization constants we use the H(4) representation τ 4 1 . In this case one obtains a singlet multiplet (renormalization factor Z S λ ) and two doublets of octet multiplets, which can mix under renormalization (with a 2 × 2 renormalization matrix Z Oλ mm ′ ). Thus, we have At the flavor symmetric point, where λ Λ T = λ Λ 1 , the second equation reduces to the first one, which is then valid for all octet baryons.
In the case of the first moments of the leading twist DAs we work with the H(4) representation τ 12 2 . Here all three representations of S 3 appear: one singlet multiplet (renormalization factor Z S ϕ ), four doublets of octet multiplets (renormalization matrix Z Oϕ mm ′ ) and three decuplet multiplets (renormalization matrix Z Dϕ mm ′ ). The resulting renormalization pattern is the following: In the SU(3) symmetric limit π B≠Λ 00,(1) = ϕ B 00,(1) , π B≠Λ

Chiral extrapolation and SU(3) flavor breaking
For the chiral extrapolation we use the three-flavor baryon chiral perturbation theory expressions derived in ref. [27]. All data points used in the present study have approximately physical average quark mass. Therefore, we use the simplified version of the extrapolation formulas, where the mean quark mass is kept fixed and all quantities are expanded around the flavor symmetric point. This scenario corresponds to the green line of figure 1. Using the average octet baryon mass X b ≡ (2m N +3m Σ +2m Ξ +m Λ ) 8, we define the dimensionless quantity to parametrize this path in a natural way starting from the flavor symmetric point at δm = 0 and hitting the physical point at δm phys ≈ 0.228. For the leading twist DAs, defined in eq. (2.13), the extrapolation formulas read The formulas for the higher twist normalization constants are similar: 5 (5.3c) 5 In eq. (5.3c) the subscript Ξ refers to the chiral even higher twist DAs Ξ B ±,4 5 , see ref. [27].
The functions g B DA (δm) contain the nonanalytic contributions from the leading one-loop diagrams and of the wave function renormalization. These are normalized such that g B DA (0) = 1, which means that Φ ⋆ ± are the distribution amplitudes at the flavor symmetric point given in eq. (2.14). The functional form of g B DA (δm) is known and can be extracted from eq. (5.18) of ref. [27]. 6 It is important that these nonanalytic terms entering as multiplicative factors do not depend on the quark momentum fractions. The remaining quantities Φ ⋆ ± , ∆Φ B ± , ∆Π B , λ ⋆ 1,2 , ∆λ B 1,2 and ∆λ Λ T play the role of low energy constants, meaning that they are independent of δm. However, note that Φ ⋆ ± , ∆Φ B ± , ∆Π B still depend on x 1 , x 2 , x 3 and their functional forms cannot be predicted by an effective low energy theory. Chiral perturbation theory [27] imposes, however, certain relations between the DAs ∆Φ B ± and ∆Π B which parametrize the SU(3) breaking: Analogously, the SU(3) breaking parameters of the higher twist couplings should satisfy the constraints In the following we will call fits to the lattice data constrained, if the relations (5.4) are imposed, and unconstrained otherwise. Combining eqs. (5.4) with the explicit form of g B DA (δm) one can find specific linear combinations of DAs for which all terms linear in δm cancel so that the SU(3) breaking is minimized. Similar combinations exist for the baryon masses: The first relation is the famous Gell-Mann-Okubo (GMO) sum rule for baryon masses [57], whose almost exact realization in nature is widely known. The second one cannot be checked for the physical masses since it depends on m ⋆ , the baryon mass at the flavor symmetric point, which is inherently inaccessible in experiment. The analogous expressions for the leading twist DAs read: For appropriately defined higher twist DAs (see ref. [27]) similar relations hold. For the normalization constants one has where the first and the second equation follow directly from eq. (5.6), and the remaining ones result from the corresponding relations for higher twist DAs. The similarity between the relations (5.6c) and (5.7c) is not accidental. A chiral perturbation theory analysis [27] reveals that the DAs of arbitrary twist can be categorized into classes with definite chiral behavior, and λ B 1 is the normalization of several higher twist DAs within the same class as Φ B − . Similar relations hold for the higher moments of the DAs. In particular, the relations for the first moments of the leading twist DAs are obtained from eq. (5.6) by replacing Φ + ↦ ϕ 11 , Π B≠Λ ↦ π B 11 , Φ − ↦ ϕ 10 and Π Λ ↦ π Λ 10 . To visualize the size of higher order SU(3) breaking terms it is convenient to form dimensionless expressions that vanish in the flavor symmetric limit (δm → 0): In figure 4 we show linear and quadratic fits to the data. Even though for all these combinations the expected δm dependence is quadratic, we find that a linear dependence cannot be excluded. The largest deviation at the physical point is found for δ 2 f (up to ≈ 15%). Most remarkably, the deviation from the GMO-like relation for the leading twist DAs, δ 1 f , is very small ( δ 1 f ≈ 1% at the physical point). For comparison, the violation of the GMO sum rule (5.5a) using the experimental values of baryon masses is In figures 5-10 we show constrained (left) and unconstrained (right) combined fits to the lattice data. For most of the measured quantities we find that the constraints in  and λ Λ T , as well as for the first moments of Φ B − and Π Λ (ϕ B 10 and π Λ 10 ), which have the same chiral behavior as λ B 1 and λ Λ T , one finds an extraordinarily good agreement with the lattice data (cf. figures 6 and 10). Also for the first moments ϕ B 11 and π B≠Λ 11 , which appear in Φ B + and Π B≠Λ , and are predicted to have the same chiral logarithms as the couplings f B and f B≠Λ T , the constraints are fulfilled within errors (cf. figure 9). In contrast, for the leading twist normalization constants f B and f B≠Λ T , as well as for ϕ B 00,(1) and π B≠Λ 00,(1) (which have to coincide with f B and f B≠Λ T in the continuum), these relations seem to be broken rather badly (cf. figures 5 and 8). Also for λ B 2 the agreement is not really flawless (cf. figure 7).
We can summarize that leading one-loop BChPT can qualitatively describe our data, even though in some cases the observed SU(3) breaking cannot be reproduced by the constrained fit. This might indicate that for these quantities higher order BChPT effects are particularly large. However, the observed discrepancies could also be caused by systematic errors in the lattice data, for which finite volume and discretization effects are prominent candidates. In particular lattice spacing effects have already been identified as a major source of systematic uncertainty in the two-flavor calculation of ref. [17], where it was also argued that for the leading twist normalization constants discretization effects are expected to be larger than for the higher twist couplings.
A heuristic parametrization of the leading discretization effects can be constructed by introducing a multiplicative factor into the extrapolation formulas. The leading corrections are linear in the lattice spacing, since the operators we use are not O(a) improved. At fixed mean quark mass this would yield, for instance, for the leading twist couplings: The constant C has to be equal for all baryons in the octet while the D B (T ) can be dif-   Figure 6. Constrained fit (left, 4 parameters) and unconstrained fit (right, 6 parameters) for the chiral odd higher twist normalization constants λ N 1 , λ Σ 1 , λ Ξ 1 , λ Λ 1 and λ Λ T .  Figure 9. Constrained fit (left, 4 parameters) and unconstrained fit (right, 7 parameters) for the first moments ϕ N 11 , ϕ Σ 11 , π Σ 11 , ϕ Ξ 11 , π Ξ 11 and ϕ Λ 11 of the leading twist DAs Φ B + and Π B≠Λ .   (1)(5) 105 (1)(5) ferent and are not necessarily subject to the same constraints as ∆f B (T ) . One can easily convince oneself that, at nonzero lattice spacing, terms O(aδm) can override the effect of the constraints given in eq. (5.4).
In this work we only use data at a single lattice spacing and cannot study discretization effects. Therefore, for the time being, the difference between chiral extrapolations using constrained and unconstrained fits has to be interpreted as evidence for systematic uncertainties.

Results
The results of the chiral extrapolations as shown in figures 5-10 are summarized in table 4 (constrained fit) and table 5 (unconstrained fit). For all quantities the first error refers to a combined statistical and extrapolation error, and the second error is an estimate of the uncertainty due to the renormalization procedure as described in appendix C. We do not expect significant finite volume effects [17,58,59] since all our ensembles have values of m π L ≳ 4 and at the same time L > 2.7 fm, cf. table 2. As discussed above, for some quantities the difference between constrained and unconstrained chiral extrapolations is sizable and can be viewed as part of the systematic uncertainty. Since the overall quality of the unconstrained fit is better (χ 2 per degree of freedom is smaller than 1.5 for all unconstrained fits), we present the corresponding numbers as our final results for this lattice spacing (see table 5). All further tables and figures in this section are generated using these values.
Our results for the nucleon normalization constants (at a ≈ 0.0857 fm with N f = 2 + 1) are approximately 30% larger for f N and about 20% larger in the case of λ N 1 and λ N 2 ,  (1) (1)(5) 86 (1)(4) 96 (1)(5) 101 (1)(5) in comparison to the N f = 2 lattice study [17], where a continuum extrapolation was performed. As one can see from figure 7 therein, 7 the continuum extrapolation from lattices with a ≈ 0.06 − 0.08 fm resulted in a decrease of f N by ≈ 30% and a somewhat smaller decrease for λ N 1,2 , so that our results are in fact very compatible. Given that we use a similar lattice action, we have to expect discretization effects of the same magnitude as in [17], and therefore, a thorough continuum extrapolation will be of utmost importance and is a primary goal for future studies. Note, however, that our results for the momentum sums ϕ B 00, (1) and π B 00,(1) defined in eq. (3.14) are within 5% of the corresponding couplings, cf. eq. (3.15), indicating that discretization errors in the derivatives are under control, see also figure 8 in ref. [17].
Our results for the first order shape parameters of the leading twist DA of the nucleon, ϕ N 11 = π N 11 and ϕ N 10 , agree with the results of ref. [17] within errors, 8 and also with the parameters extracted from the study of the nucleon electromagnetic form factors in lightcone sum rules [14]. Note that our ϕ N nk correspond to f N ϕ N nk in refs. [14,17]. We also confirm the approximate equality ϕ N 10 ≈ ϕ N 11 found in ref. [17]. Our results for the shape parameters of hyperons are, however, up to an order of magnitude smaller than the values obtained using QCD sum rules [6], see table 6. In ref. [17] it has already been reported that, in general, modern lattice simulations and light-cone sum rule calculations yield estimates of the first moments of the nucleon DA that are one order of magnitude smaller than in "old" phenomenology, cf. refs. [3,6]. Our measurements confirm this observation also for Table 6. Comparison of the central values of our N f = 2 + 1 results (unconstrained fit, see table 5) with the N f = 2 lattice study for the nucleon [17] and the Chernyak-Ogloblin-Zhitnitsky (COZ) model [6]. All values are given in units of GeV 2 . All quantities have been converted to the conventions established in this work and rescaled to µ 2 = 4 GeV 2 , using the three-loop evolution equation for the couplings with the anomalous dimensions calculated in ref. [34], and the one-loop equation (2.22) for the shape parameters. Note that f T Λ in ref. [6] is proportional to the first moment π Λ 10 in our nomenclature.
B work method Interestingly, the SU(3) breaking in the shape parameters of the octet baryons turns out to be very large, e.g., π Ξ 11 ≳ 3ϕ N 11 and ϕ Λ 10 ≳ 7ϕ N 10 . This effect is much stronger than in QCD sum rule calculations [6], even though the absolute values are much smaller. This large SU(3) breaking is somewhat surprising as shape parameters have autonomous scale dependence and should be viewed as independent nonperturbative parameters, and is in stark contrast to the situation for the normalization constants where the differences between octet baryons are at most 50%. As a consequence, SU(3) breaking in hard exclusive reactions that are sensitive to the deviations of the DAs from their asymptotic form can be enhanced.
The SU(3) breaking in the shape of the leading twist DAs can be represented in many ways. Consider, e.g., normalized combinations of symmetric and antisymmetric DAs all of which are equal both in the asymptotic limit, φ as ≡ φ B,as = B,as = 120x 1 x 2 x 3 , and in the limit of exact flavor symmetry, φ ⋆ ≡ φ B⋆ = B⋆ . Due to isospin symmetry N = φ N . Hence, there are seven independent functions that can be used to visualize the deviations from the DA φ ⋆ in the SU(3) flavor symmetry limit. These seven functions, φ N − φ * , etc., are shown in figure 11 together with φ ⋆ itself, which is almost (but not exactly) symmetric in x 1 , x 2 , x 3 due to small (but nonvanishing) values of ϕ ⋆ 11 and ϕ ⋆ 10 (cf. figures 9 and 10). In phenomenological applications it is more convenient to consider the standard representation of DAs in terms of [V − A] B and T B . In this way also the physical interpretation  is more straightforward as every momentum fraction can be attributed to a quark of definite helicity and flavor. [V −A] B and T B do not coincide, however, at the flavor symmetric point, so that for these DAs it is more natural to show the deviations from the asymptotic shape φ as rather than from φ * , see figure 12. The plots in the left and in the right column show normalized DAs [V − A] B and T B after the subtraction of the asymptotic DA. Note that the amplitudes T B≠Λ are symmetric under the interchange of x 1 and x 2 by construction. The approximate symmetry of [V − A] N under the exchange of x 2 and x 3 is, in contrast, nontrivial. It is due to the approximate equality of the two nucleon shape parameters ϕ N 10 ≈ ϕ N 11 mentioned above. In the nucleon Fock state u ↑ u ↓ d ↑ this is equivalent to a symmetric distribution of momentum between the second and third quark. In agreement with earlier studies [4,6,17], we observe that the "leading" u ↑ quark, which has the same helicity as the nucleon, carries a larger momentum fraction. In the u ↑ u ↑ d ↓ nucleon state, which is described by T N , the peak of the distribution is shifted towards the two u quarks in a symmetric manner. T N , however, is not an independent DA. Taking into account the isospin relation (2.12), the spin-flavor structure of the nucleon light-cone wave function (2.16) can be presented, schematically, as In this picture our result for [V − A] N corresponds to a shift of the momentum distribution towards the u ↑ quark, which carries the nucleon helicity, and the symmetry under x 2 ↔ x 3 may be interpreted as an indication for the remaining valence quarks forming a dynamical scalar "diquark", which is assumed in many models.
For the Σ baryon state d ↑ d ↓ s ↑ one sees that the maximum of the distribution is shifted from d ↓ towards s ↑ , whereas in the d ↑ d ↑ s ↓ state the s quark gathers additional momentum from both d quarks equally. The overall size of the deviations from the asymptotic distribution is, however, quite small, smilar to the nucleon case. For the Ξ baryon the deviations are slightly larger. In the s ↑ s ↓ u ↑ state, the distribution is tilted towards the s ↑ quark and leaves less momentum for the u ↑ quark. T Ξ is clearly dominated by the two s quarks. In summary, for the isospin-nonsinglet baryons one can identify two competing patterns: First, the strange quarks carry, in general, a larger fraction of the momentum. Second, in the ↑↓↑⟩ state the first quark is favored over the second, while in the ↑↑↓⟩ state the first two quarks behave identically. These rules do not apply to the Λ baryon due to its reversed symmetry properties, see eq. (2.11): In the u ↑ d ↓ s ↑ state the maximum of the distribution is shifted towards the s quark. T Λ is a special case, since it does not contain the leading asymptotic part due to the antisymmetry under exchange of x 1 and x 2 . Hence, for the Λ baryon, the Fock state u ↑ d ↑ s ↓ is expected to be highly suppressed.
In order to quantify this picture, we consider normalized first moments of which are sometimes referred to as momentum fractions in the literature. Note that this name is imprecise since the averaging is done with the DA and not a wave function squared, and, in particular, for T Λ , which has no asymptotic part, the interpretation as momentum fractions breaks down completely. The ⟨x i ⟩ can be calculated in terms of the shape 296 parameters as follows: (1) .
The results are summarized in table 7. They support the qualitative picture suggested by the discussion of figure 12. Finally, we consider the higher twist matrix elements that are related to the normalization of the P -wave light-cone wave functions and also appear as low energy constants in effective theories for generic GUT models [35]. We obtain, for the nucleon, λ N 2 ≈ −2λ N 1 , which is well known, see, e.g., refs. [16,38,60]. The same relation also holds for the Σ and Ξ hyperons but not for the Λ baryon. Instead, we find λ Λ 2 ≈ −2λ Λ T , i.e., the matrix element in eq. (2.27b) is zero within the error bars. The likely interpretation (similar to the familiar relations for isospin-nonsinglet baryons) is that the corresponding matrix elements vanish in the nonrelativistic quark model limit.

Conclusions and outlook
In this work we have performed the first N f = 2+1 lattice QCD analysis of the normalization constants and (leading twist) first moments of the octet baryon distribution amplitudes with pion masses down to 222 MeV. The results are scheme-and scale-dependent and, thus, have to be renormalized. To this end we first carried out a nonperturbative renormalization in a RI ′ SMOM scheme, followed by a conversion to the MS scheme applying continuum perturbation theory at one-loop accuracy. We extrapolated our results to the physical point using three-flavor BChPT formulas derived in ref. [27].
We find significant SU(3) flavor breaking effects for the leading twist normalization constants (7.1) and somewhat smaller symmetry breaking for the higher twist couplings where the number in parentheses gives a combined statistical and chiral extrapolation error, while the uncertainty from the renormalization procedure is negligible for these ratios. It is likely that these ratios are less sensitive to discretization effects than the couplings themselves.
Deviations from the asymptotic DAs are quantified by the values of shape parameters. They are small for all baryons in the octet, in agreement with the findings of ref. [17] for the nucleon, and much smaller than results of old QCD sum rule calculations [6]. The SU(3) breaking in the shape parameters is, however, very large, see table 5. For the isospinnonsinglet baryons one can identify two competing patterns: First, the strange quarks carry, in general, a larger fraction of the momentum. Second, in the f ↑ g ↓ h ↑ state (using our flavor conventions (2.3)) the first quark is favored over the second, while in the f ↑ g ↑ h ↓ state the first two quarks behave identically. These rules do not apply to the Λ baryon due to its reversed symmetry properties, see eq. (2.11). The interplay of these two patterns leads to the rather elaborate structure shown in figure 12.
To first order in the SU(3) symmetry breaking parameter we have derived the following relation between the DAs of the Σ and Ξ hyperons: This relation has the same theory status as the renowned Gell-Mann-Okubo relation for the masses, and is satisfied with similarly high accuracy ∼ 1% in our data. The analysis presented here, using a trajectory with fixed mean quark mass, constitutes the first half of the twofold strategy pursued by the CLS effort. It will be complemented by a second set of lattices at fixed physical strange quark mass as indicated by the red line in figure 1. The extrapolation to the physical point along this second path can be described using chiral perturbation theory with only two flavors, while any other path requires a full SU(3) treatment. The combination of these two methods will provide one with an additional tool to estimate systematic errors. Its full implementation lies beyond the scope of this work, where we have focused on the development of the necessary formalism to describe patterns of SU(3) breaking at the wave function level. Future studies will have access to a rich landscape of CLS ensembles along both trajectories, including ensembles at (nearly) physical quark masses and various lattice spacings down to a ≈ 0.04 fm, thus allowing for a reliable continuum extrapolation.

A The baryon octet
Starting with the standard representation for the quark triplet we define lowering operators T − , U − and V − for the isospin, U -spin, and V -spin, respectively, in this (fundamental) representation as so that The baryon octet is usually presented as [62] ⎛ where K B are matrices in flavor space, e.g., We further define the action of the lowering operators T − , U − and V − on the octet by the usual expressions for the adjoint representation aŝ without any additional phase factors. The above choices specify our phase conventions. Starting from the proton state, the complete octet can be constructed by applying the following transformations as illustrated in figure 13: Figure 13. Illustration of our phase conventions. The Λ baryon is not shown since one needs a linear combination for its construction, cf. eq. (A.7g). Blue arrows indicate the cases where one has to apply a Fierz transformation (see ref. [30]) to relate the distribution amplitudes at the symmetric point. An explicit calculation shows that this always yields an additional minus sign that has to be taken into account in order to reproduce eq. (A.10) and eq. (2.14).
Starting from the mixed-symmetric and mixed-antisymmetric flavor wave functions for the proton defined as our conventions also fix the relative signs of the baryon DAs. As shown in ref. [27] this choice corresponds to in the limit of exact isospin symmetry. This also fixes the relative phases at the flavor symmetric point in eq. (2.14). All phases are now unambiguously determined up to a single unphysical global phase, which is commonly fixed by the condition that f N has to be positive.  Table 9. Mixed-symmetric (MS) flavor wave functions.  Table 11. Totally symmetric (S) flavor wave functions.

B.1 Relation to H(4) operators
In the following we will relate the operators defined in (3.6) and (3.10) to those of ref. [49]. It is implied that within the generic operators appearing on the right hand side of the equations, the quark flavors f , g and h are chosen such that they agree with the convention for the baryon B, see (2.3). For the operators without derivatives we have right hand side the position of the derivative is set as mandated by the superscripts lmn: where the operators for the structure V + A (or V − A) can be obtained by replacing Similarly, the operators which are relevant for higher twist normalization constants (see eq. (3.16)) can be expressed in terms of O 1−5 . In the chiral odd sector we have In the chiral even sector (λ B 2 ) we obtain:

B.2 Operator bases for renormalization
For the purpose of renormalization it is convenient to employ operator multiplets that transform irreducibly not only with respect to the spinorial hypercubic group H(4) but also with respect to the group S 3 of permutations of the three quark flavors. The latter group has three inequivalent irreducible representations, which we label by the names of the corresponding ground state particle multiplets in a flavor symmetric world. Therefore, the one-dimensional trivial representation is labeled by D in the main text, the one-dimensional totally antisymmetric representation by S and the two-dimensional representation by O. We construct multiplets with the desired transformation properties from the multiplets defined in ref. [49]. For operators without derivatives in the representation τ (with the first multiplet being mixed-symmetric and the second one being mixed-antisymmetric) and one operator multiplet transforming trivially under S 3 : For operators without derivatives in the H(4) representation τ 4 1 we have one multiplet that is totally antisymmetric under flavor permutations, and two doublets of operator multiplets transforming according to the two-dimensional representation of S 3 : For operators with one derivative in the H(4) representation τ 12 2 we have one multiplet that is totally antisymmetric under S 3 , four doublets of operator multiplets corresponding to the two-dimensional representation of S 3 ,

C Renormalization procedure
Every local three-quark operator can be represented as a linear combination of the operators with the same multi-index notation as above. Aiming at a mass-independent renormalization scheme we assign the same mass to all flavors and eventually consider the chiral limit where this mass is sent to zero. In order to conveniently display the behavior of the operators under permutations of the three quarks we write the above operators in the form or, in an abbreviated notation, as Ψ f α (l; x). Then we have for all permutations π in the symmetric group S 3 of three elements, where From these "elementary" operators we construct the operators of interest with the help of flavor coefficients F and spinor coefficients S according to where a sum over all (multi-)indices which appear twice is implied. Under SU(3) the quark fields transform according to the fundamental representation 3 and for our three-quark operators we have the decomposition 3 ⊗ 3 ⊗ 3 = 1 ⊕ 8 ⊕ 8 ⊕ 10. The flavor-singlet (flavor-decuplet) representation corresponds to the totally antisymmetric (totally symmetric or trivial) representation of S 3 . The two flavor octets, called mixedsymmetric (MS) and mixed-antisymmetric (MA), form a basis for the two-dimensional representation of S 3 . More explicitly, we have the singlet flavor structure The spinor structures should be chosen to yield a flavor-spinor structure that is totally symmetric under simultaneous permutations of the flavor, spinor and derivative indices f a , α a andl a (a = 1, 2, 3). Starting from the operator multiplets given in ref. [49], which Here m labels the different H(4) multiplets and i labels the different members of the multiplets. Then for t ∈ {1, 2}. Therefore, we can always work with the MA flavor structure (t = 2) and assume that the flavor-spinor structure factorizes into a flavor structure and a spinor structure as in (C.5). For the singlet and decuplet baryons this factorization is trivially satisfied. Let us now describe our renormalization procedure, which is similar to the well-known RI ′ SMOM scheme. In particular, we compute the quark field renormalization factor Z q from the quark propagator as in ref. [63].
Under renormalization, multiplets of operators transforming according to the same representation of H(4) and having the same or lower dimension will in general mix. Since mixing with operators of lower dimension is particularly problematic we have chosen the operators such that this type of mixing is absent. Moreover, there is no mixing between operators transforming according to inequivalent representations of S 3 .
We write the mixing operator multiplets for a fixed flavor structure F B,f o2 in the form The renormalization and mixing coefficients Z mm ′ are fixed by the renormalization condition (C. 19) Here and in the following the superscript "Born" indicates the corresponding tree level expression (Born term). This is taken with all lattice artefacts included. More explicitly, our renormalization condition can be written as This procedure yields (matrices of) renormalization factors leading from the bare operators on the lattice to operators renormalized according to the MOM scheme introduced above. However, in the applications we need operators renormalized in the MS scheme. These are constructed with the help of (matrices of) conversion factors calculated in continuum perturbation theory, where we use the particular version of the MS scheme introduced in ref. [29]. Due to the complexity of higher-loop calculations we had to limit ourselves to one-loop accuracy. Also the anomalous dimensions of our operators are in general only known to one loop, with the exception of the operators without derivatives, for which the anomalous dimensions have been calculated to three loops [34].
Let us now give a few more technical details of the computation of the renormalization matrices. The required propagators have been evaluated on gauge field ensembles which were generated on lattices of size 32 4 for four different quark masses using periodic (antiperiodic) boundary conditions for the quark fields in spatial (temporal) direction and periodic boundary conditions for the gauge fields in all four directions. For the external momenta we have taken employing twisted boundary conditions. The extrapolation to the chiral limit is performed linearly in the square of the pseudoscalar mass. The values to be used in the analysis of the physical matrix elements are determined by interpolating the (chirally extrapolated) numerical data linearly in µ 2 to our target scale µ 2 0 = 4 GeV 2 . While the statistical errors are quite small, systematic uncertainties are far more important. In order to estimate their impact we proceed as follows. Let the value of the renormalization coefficient under study be z 0 . We consider two additional procedures to determine this coefficient. In the first procedure we take the (interpolated) data at 10 GeV 2 and evolve them down to 4 GeV 2 with the help of the perturbative renormalization group using as many loops in the anomalous dimensions as are available. Let us call the resulting number z 1 . In the second procedure we read off our result directly at 4 GeV 2 , but use for the perturbative conversion to the MS scheme one loop order less than before. (This means in our case that the conversion matrix is set to unity.) This gives the value z 2 . Then we perform the complete analysis with the three choices z 0 , z 1 , z 2 for the needed renormalization and mixing coefficients. Defining δ i as the difference between the outcome of the analysis employing z i and the outcome of the analysis employing z 0 , we estimate the systematic uncertainties due to the renormalization factors as δ 2 1 + (δ 2 2) 2 . Here we have multiplied δ 2 by 1 2, because going from one to two or more loops in the conversion factors is expected to yield a smaller change than going from zero loops to one loop.