Logarithmic corrections to O(a) and O(a2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^2$$\end{document}) effects in lattice QCD with Wilson or Ginsparg–Wilson quarks

We derive the asymptotic lattice spacing dependence an[2b0g¯2(1/a)]Γ^i\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^n[2b_0\bar{g}^2(1/a)]^{\hat{\Gamma }_i}$$\end{document} relevant for spectral quantities of lattice QCD, when using Wilson, O(a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{O}(a)$$\end{document} improved Wilson or Ginsparg–Wilson quarks. We give some examples for the spectra encountered for Γ^i\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\Gamma }_i$$\end{document} including the partially quenched case, mixed actions and using two different discretisations for dynamical quarks. This also includes maximally twisted mass QCD relying on automatic O(a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{O}(a)$$\end{document} improvement. At O(a2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{O}(a^2)$$\end{document}, all cases considered have miniΓ^i≳-0.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\min _i\hat{\Gamma }_i\gtrsim -0.3$$\end{document} if Nf≤4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\textrm{f}}\le 4$$\end{document}, which ensures that the leading order lattice artifacts are not severely logarithmically enhanced in contrast to the O(3) non-linear sigma model (Balog et al. in Nucl Phys B 824:563–615, 2010; Balog et al. in Phys Lett B 676:188–192, 2009). However, we find a very dense spectrum of these leading powers, which may result in major pile-ups and cancellations. We present in detail the computational strategy employed to obtain the 1-loop anomalous dimensions already used in Husung et al. (Phys Lett B 829:137069, 2022).


Introduction
Today's lattice QCD simulations are reaching a point, where statistical uncertainties of several quantities obtained from the lattice become of O(1%) and below while smaller lattice spacings a 0.04 fm become accessible. Under these conditions predictions for precision physics are within reach, but systematic errors must be kept under control. We focus here on the systematic error due to the continuum extrapolation. The method of choice to describe the approach of the lattice theory to the continuum theory as a 0 is Symanzik Effective theory (SymEFT) [1][2][3][4], see also [5, p. 39ff.]. SymEFT allows to take the quantum corrections into account, that modify the leading a n min lattice artifacts a e-mail: n.husung@soton.ac.uk (corresponding author) expected from classical field theory, where n min is a positive integer and depends on the chosen lattice discretisation. In an asymptotically free theory like QCD the leading lattice artifacts from the lattice action then have the generic form a n min [ḡ 2 (1/a)]γ i , whereḡ(1/a) is the renormalised coupling, γ i = (γ 0 ) i /(2b 0 ), (γ 0 ) i are the 1-loop anomalous dimensions of (4+n min )-dimensional operators and b 0 is the 1-loop coefficient of the β-function. These higher dimensional operators form a minimal on-shell basis describing all lattice artifacts originating from the lattice action, which can formally be written in form of the effective Lagrangian L eff (x) = L (x) + aδL (1) (x) + a 2 δL (2) (x) + · · · , (1) where L is the Lagrangian of continuum QCD = (ψ 1 , . . . , ψ N f ) T is a flavour vector, D μ (A) = ∂ μ + A μ is the continuum covariant derivative with su(N ) algebra valued gauge field A μ , F μν = D μ , D ν is the field strength tensor and M = diag(m 1 , . . . , m N f ) are the quark masses. To obtain the leading lattice artifacts we then need to expand the SymEFT around the continuum Lagrangian such that the δL (d) are treated as operator insertions in the continuum theory, which is the common strategy in Effective Field Theories. Here and in the following the superscript in δL (d) denotes the deviation of the canonical mass-dimension from the continuum field, e.g.
We focus on contributions from the lattice action. This is sufficient for spectral quantities, where corrections from local fields cancel out. For non-spectral quantities such contributions from local fields must be taken into account as well, which is beyond the scope of this paper and obviously depends on the local fields involved. 1 As shown by Balog et al. [7,8] in the 2-d O(3) sigma model quantum corrections can spoil the approach to the continuum limit with distinctly negative values forγ i . For this model they found min i (γ i ) = −3, which worsens the approach from the naive a 2 lattice artifacts (n min = 2) to something which behaves like O(a) corrections over a long range of lattice spacings due toḡ 2 (1/a) ∼ −1/ log(a ), where is the intrinsic scale of the theory. Thus computing the leading anomalous dimensions for full lattice QCD at leading order in the lattice spacing is not just a purely academic question but puts continuum extrapolations on more solid grounds by predicting the true asymptotic lattice spacing dependence. This knowledge should then be used both in the ansatz for the continuum extrapolation and to estimate uncertainties of this extrapolation through varying the leading power in the coupling in the range predicted for the competing valuesγ i .
In a previous paper [9] we discussed the special case of SU(N ) pure gauge theory with n min = 2 and Wilson's lattice QCD with n min = 1, which have a minimal basis of 2 or 1 operators respectively in the massless case. In both cases the values found forγ i are larger than zero such that discretisation errors vanish faster than the classically assumed a n min behaviour. There the general concept of SymEFT theory is discussed in more detail. It is recommended as an introduction for the reader. In this paper we will focus on the extension to lattice QCD actions with non-perturbatively O(a) improved Wilson quarks [10] and Ginsparg-Wilson quarks [11] both with n min = 2, which have a significantly larger set of operators forming the minimal basis. We also take a look at twisted mass QCD [12,13] and the connection of the minimal operator basis of untwisted QCD to the twisted case. In particular twisted mass QCD at maximal twist ensures automatic O(a) improvement due to an additional symmetry of the continuum Lagrangian [14,15]. While we discuss here the technicalities of the computation in some detail, we also recommend to take a look at the Letter published [16] alongside this paper, summarizing the overall results and consequences.

Minimal operator basis of full QCD to O(a )
To compute the variousγ i we need the complete minimal on-shell operator basis O (d) i to express 1 For commonly used local fields a similar analysis can (and should) be performed along the lines of this paper with the sole difference that total divergence operators can no longer be discarded in the additional operator bases introduced for each local field. See also [6] for a discussion of the minimal on-shell basis of local fields for n min = 1.
where ω O i (g 2 0 ) are the bare matching coefficients with bare continuum coupling g 2 0 . The expansion in g 2 0 of ω O i can be determined for any chosen lattice discretisation. Which operators must be included into the basis, depends on the symmetries realised for the lattice discretisation, i.e. the lattice action Different discretisations of the lattice gauge action S G have already been discussed in [9] and we rather focus here on the lattice fermion actions. Depending on the chosen lattice Dirac operatorD the symmetry constraints differ. We restrict considerations to Wilson quarks [17,18] and Ginsparg-Wilson quarks [11]. Starting with Wilson quarks the lattice Dirac operator readŝ where σ μν = i 2 γ μ , γ ν . Here c sw (g 2 0 ) = 1 + O(g 2 0 ) is the improvement coefficient for the Sheikholeslami-Wohlert (SW) term [19]. For the definition of ∇ μ , ∇ * μ andF μν see appendix A. The SW-term can remove O(a) lattice artifacts either perturbatively [6] or non-perturbatively [10], where the latter choice also removes double mass-dimension 5 operator insertions in the SymEFT contributing to O(a 2 ).
In contrast to Wilson quarks, Ginsparg-Wilson quarks [11] obey in the massless limit D GW , γ 5 such that they have an exact lattice chiral symmetry [20]. This symmetry ensures that terms like i/4¯ σ μν F μν are already forbidden by symmetry and the lattice artifacts of the massless theory automatically start at O(a 2 ). One particular solution to Eq. (7) are Overlap fermions [21,22] (we choose here the conventions from [23]) Another solution of Eq. (7) are Domain-Wall fermions [24,25] in the limit of infinite extent of the auxiliary 5th dimension [26]. For finite extent of the 5th dimension, chiralsymmetry violations are only exponentially suppressed as the extent of the 5th dimension increases such that flavour symmetries of Domain Wall fermions are then reduced to the ones of conventional Wilson quarks. In summary, the symmetry constraints for Wilson and Ginsparg-Wilson quarks are the following: • Local SU(N ) gauge symmetry, • invariance under charge conjugation and any Euclidean reflection, • Hypercubic H 4 symmetry as a remnant of broken O (4) symmetry, Due to being interested in the minimal on-shell basis we may further make use of the continuum equations of motion (EOM) to eliminate redundant operators [6]. Here T a denotes the generators of the su(N ) colour algebra.

Massless operator basis
For massless Wilson quarks with at most perturbative O(a) improvement, one operator is required at mass-dimension 5 describing the occurring O(a) lattice artifacts. We previously discussed this operator [9] and list it here for completeness For mass-dimension 6 we restate here the minimal basis of pure gauge theory [27,28], which is of course a subset of the operator basis of full QCD, The extension to full QCD with massless quarks at O(a 2 ) then introduces an additional fermion bilinear compatible with chiral symmetry to the minimal on-shell basis after applying the EOMs and making use of integration by parts on the basis listed in [19] Both operators O (2) 2 and O (2) 3 break O(4) symmetry and are only compatible to H 4 symmetry. Therefore their mixing under renormalisation will be fairly restricted as they cannot mix into any operator invariant under O(4) symmetryassuming use of a regulator that preserves this symmetry. At mass-dimension 6 also 4-fermion operators contribute. For Ginsparg-Wilson quarks only those compatible with chiral symmetry are allowed while Wilson quarks also require the inclusion of the chiralsymmetry violating 4-fermion operators We choose here the 4-fermion operator basis such that the flavour, colour algebra and spinor indices are contracted within¯ . . . pairs. All other ways to contract the indices compatible with the symmetry constraints are related to the basis chosen here through the following identities. Firstly through Fierz identities 2 where ψ n is a quark field with an arbitrary combination of flavour and colour index, while the spinor indices are contracted within the quark anti-quark pair. Secondly through identities from the su(N ) algebrā where ψ and η indicate different flavours, {μ} denotes a matrix of the Dirac algebra with all indices {μ} contracted with the second {μ} and A, B, C, D are contracted colour indices. This particular choice for the basis is identical to the one in [30] and, as shown there, equivalent to the choice in [19]. We choose to prepend an additional factor of g 2 0 to each 4-fermion operator motivated by the gluonic EOM (10) as well as the leading order of any kind of tree-graph leading to a 4-fermion interaction in lattice QCD. The latter happens due to the absence of terms with more than one quark-antiquark pair in the classical expansion of the lattice action in the lattice spacing a as discussed in [19], i.e. the tree-level coefficients of 4-fermion operators without the factor g 2 0 would vanish anyway.

Massive operator basis
Unless one is interested in massless quarks or quarks at very small quark masses, like up-, down-or strange-quarks, operators carrying explicit powers of quark masses should be taken into account because they will no longer be suppressed. If a hadronic renormalisation scheme is used on the lattice rather than a mass-independent scheme, such contributions can be ignored. For the general massive case one finds at mass-dimension 5 the additional linearly independent operators where the normalisation with N f is chosen to ensure that in the mass-degenerate case all operators carry the same prefactor. Similarly, the basis of mass-dimension 6 operators requires the inclusion of the larger number of additional linearly independent massive operators The way an explicit mass term added to the lattice action reduces the flavour symmetries of Ginsparg-Wilson fermions to U(1) N f V implies that only the operators O (2) 14−17 are allowed to contribute. 3 To see how this simplification arises, we first rewrite the mass term as where the subscripts indicate left-and right-handed quarks.
Promoting the mass matrix M to a spurionic field that transforms according to ensures that the action stays invariant under the flavour rotations For higher powers of the mass matrix in our SymEFT basis, we replace iteratively all occurrences of M 2 → M M † from the left. Requiring invariance under the spurionic symmetry transformation then indeed leaves only the O 14−17 as stated in the beginning. For Wilson quarks however, all massive mass-dimension 6 operators listed here are needed, because chiral symmetry was already explicitly broken by the lattice discretisation of the massless theory.

Wilson quarks with a chiral twist
When discussing Wilson quarks also Wilson quarks with a chirally twisted mass term [12,13,31] known as twisted mass QCD (tmQCD) come to mind. The lattice fermion action reads where χ are the chirally twisted flavours and the Dirac operator is defined aŝ Here τ 3 is the Pauli matrix τ 3 = diag(1, −1) acting in flavour space and μ q is the twisted mass parameter. The generalisation to mass-non-degenerate flavour doublets exists [32], but lies beyond the scope of this paper. The reasoning used here should carry over, but it relies more heavily on spurionic symmetry arguments. Also the overall renormalisation arguments for extracting the anomalous dimension matrix must be revisited. Therefore we limit considerations here to the case of multiple mass-degenerate flavour doublets. Due to the twisted mass term, invariance under parity transformation and time reversal are broken. Instead the theory of the chirally twisted flavours χ (both continuum and lattice) is invariant under [12] modified paritȳ and modified time reversal where j = 1, 2, as well as the spurionic symmetry transformations This does not have any impact on the massless operator basis, which is of course unchanged as the massless cases of Wilson QCD and tmQCD are identical. For the massive basis the spurionic symmetry transformations limit the occurrence of powers of μ q severely. In the continuum theory and thus in our SymEFT twisted (χ ) and untwisted ( ) flavours are connected via the substitution =χ exp(iωγ 5 τ 3 /2), = exp(iωγ 5 τ 3 /2)χ , as long as a regularisation is chosen, that is invariant under chiral rotations. This connection enables us to infer the anomalous dimensions of twisted operators from the untwisted operators due to, see e.g. [33], where the subscripts QCD and tmQCD denote the choice of the mass term in the action with flavours denoted as and χ respectively. Here O ext is assumed to be invariant under Eq. (26). This equivalence holds in case a mass-independent multiplicative renormalisation scheme S is used, which ensures that the mixing matrix is independent of the twist angle ω.
In contrast to the continuum action, the Wilson term in the lattice action Eq. (21) obstructs the transformation Eq. (26). This requires the inclusion of additional operators to our SymEFT. These additional operators relevant up to massdimension 6 arē dressed by an even power of the twisted mass parameter μ q to comply with the spurionic symmetry transformations in Eqs. (24) and (25). These operators are accompanied by the chirally twisted versions of¯ and i¯ σ μν , namely whereF μν = ε μνρσ F ρσ /2 is the dual field strength tensor. Both operators in Eq. (29) must be dressed by an odd power in μ q , to comply with the spurionic symmetry transformations in Eqs. (24) and (25). Any new 4-fermion operators at massdimension 6 will require an odd power in μ q , such that they eventually contribute only at mass-dimension 7, i.e. O(a 3 ), or beyond. While the reasoning in Eq. (27) is sufficient to fix the renormalisation of the massless operator basis and some mixing contributions of the operators in Eq. (29), we need additional input to take also the mixing of the additional massive operators into account. Since the computations we report here have been performed at zero twist angle we lack any information about massive mixing contributions involving μ q . The only information we have for the massive operator basis are the diagonal entries rather than the full mixing.
A special property of twisted mass QCD is automatic O(a) improvement at maximal twist, i.e. ω = π/2. This is due to the additional discrete symmetry of the continuum theory, which reads at this twist angle [14,15] While being a continuum symmetry, it is again explicitly broken on the lattice by the Wilson term, which allows for additional O(a) terms. Following the lines of [14] any operator can be split into a T 1 -even and T 1 -odd part, i.e., parts having eigenvalues ±1 under the transformation in Eq. (30) respectively. This carries over to n-point functions of operators, which can then be split into a T 1 -even and T 1 -odd part as well, where the T 1 -odd part vanishes by construction. At maximal twist all mass-dimension 5 operators parametrising the O(a) corrections are T 1 -odd. Thus any O(a) lattice artifacts vanish for quantities that are themselves T 1 -even. How-ever this does not imply that the matching coefficients ω i of the mass-dimension 5 operators are zero. Instead these operators become relevant at O(a 2 ) through T 1 -even double operator insertions but can in principle be removed through nonperturbative O(a) improvement [31] identically to untwisted Wilson quarks [10]. Without non-perturbative improvement these double operator insertions also give rise to contact terms in the SymEFT expansion to O(a 2 ).

Partially quenched QCD and mixed actions
With the results for full QCD at hand we can also consider the special case of partially quenched QCD, where the determinant of the lattice Dirac operator has been dropped for N b flavours. In perturbation theory this leads to discarding any closed fermion loops, while fermion lines contracted with local fields are kept. Conceptually this can be implemented by introducing additional bosonic fields such that closed fermion loops of the quenched flavours cancel out with closed boson loops, see e.g. [34][35][36]. Here we also introduced the short-hanď Assuming that the same lattice Dirac operatorD is used for both quenched and unquenched flavours the underlying flavour symmetries are modified to the U[ We only need to understand the constraint the graded flavour symmetry transformatioň imposes on the minimal operator basis of our effective action, wherẽ The generalisation necessary for our operators containing fermion fields can be inferred from the transformation of fermion bilinears where is flavour diagonal and, in the case of Ginsparg-Wilson fermions, preserves chiral symmetry. An analogous transformation can be given for bosonic bilinears. The special case of B =C = 0 = C =B, realises separate rotations in fermionic and bosonic flavour space but this corresponds only to a subset of the full graded flavour symmetry. We thus immediately see that fermion bilinears and 4-fermion operators must be generalised as →ˇ for the partially quenched case. This also means that there are no operators needed besides those in Eqs. (11) with the substitution to partially quenched flavours applied.
For the sake of argument mixed actions, i.e. different choices for the Dirac operator of sea (D S ) and valence (D V ) quarks, can be written as a partially quenched theory of N f unquenched and N b quenched flavours [30] The quenched flavours play the role of the valence quarks, while the unquenched flavours are the sea quarks. Due to the different discretisations of the Dirac operatorsD S and D V , the flavour symmetries are more complicated as separate flavour rotation symmetries are restricted to the quenched or unquenched flavours respectively. Thus we expect the following minimal operator basis where i;S/V are substitutes for all Dirac matrices allowed by the respective symmetry constraints with and without addi-tional colour algebra generator T a . The same holds for i;SV , but here the more restrictive symmetry constraints of either the quenched or the unquenched flavours decide, which Dirac matrices are allowed. In the massive case the minimal basis is enlarged accordingly. Instead of sea and valence quarks we could just as well introduce different discretisations for different sets of flavours. In this case the reasoning still remains the same just without the quenching.

One-loop computation of the anomalous dimension matrix
To renormalise the full QCD on-shell basis perturbatively to 1-loop order, we make use of the background field method [38][39][40][41]. There, a smooth classical background field B μ (x) is introduced as splitting the quantum field A μ (x) into the background field and the quantum fluctuations Q μ (x). Additionally the background field gauge is chosen by adding the gauge-fixing term to the continuum Lagrangian with bare gauge-fixing parameter λ 0 alongside a Faddeev-Popov term [42]. Due to the background field gauge, only the fields Q μ are gauge-fixed, while the background field itself remains unaffected. Then one-particle-irreducible (1PI) graphs with only background fields and quarks as legs, see also Fig. 1, transform manifestly gauge-covariant under gauge-transformations of the background field. This ensures absence of contributions from nongauge-invariant operators mixing into the gauge-invariant ones, which are in principle allowed due to working in the gauge-fixed theory [43,44]. Such mixing contributions are of course irrelevant when being interested in gauge-invariant observables. However, we cannot ignore contributions from gauge-invariant operators E vanishing by the equations of motion. Again such contributions vanish on-shell and the mixing of EOM vanishing operators under renormalisation is therefore restricted as where Z is the desired on-shell mixing matrix of our operator basis. For compactness we omitted here all operator indices, i.e. Z , A OE and Z E are matrices themselves. We thus need an additional minimal basis of EOM vanishing operators where the lower dimensional ones may carry additional powers of masses. With these prerequisites we can now renormalise our minimal operator basis by computing the 1PI graphs with single operator insertion as depicted in Fig. 1 using dimensional regularisation (D = 4 − 2 ) and renormalising any divergences in the modified minimal subtraction (MS) scheme [45][46][47].
Working in a mass-independent renormalisation scheme like MS ensures that the 1-loop anomalous dimensions are simply related to counterterms renormalising ultraviolet divergences. We then use the by now common trick called infrared rearrangement, see e.g. [48][49][50]. This enables us to separate the UV-divergent part of a 1-loop momentum integral from UV-finite but potentially IR divergent parts through the exact relation [48,49] 1 Here k is the loop momentum, p is an external momentum, m is a mass and we choose > 0 as a mass-scale. The second term in Eq. (42) on the r.h.s. is one power less UVdivergent. Iterating this step brings all UV-divergent terms into the standard form d D k [k 2 + ] −n k μ 1 . . . k μ l , while all UV-finite terms can eventually be dropped as they do not contribute to the 1-loop anomalous dimensions. These steps have been implemented in FORM [51] to evaluate the various 1PI graphs at 1-loop. 4 We also derived the necessary Feynman rules of the various operators in FORM and checked that we reproduce the usual Feynman rules for full QCD in the background field gauge.
Instead of inserting the flavour singlet operators, we insert variants where the flavour is kept generic and we can then build the required set of operators from these building blocks. As a consequence we are immediately able to extract the generalisations necessary for cases like partially quenched and Fig. 1 Graphical representation of all 1PI graphs of fundamental quark fields or background fields with insertion of an operator O or E considered for the renormalisation of the basis. Here the double line indicates the momentum contribution of the inserted operator, which may also be seen as an additional leg and is set to zero momentum. The wiggly lines are external background fields and the straight lines carrying arrows are quarks. The graph (e) is only needed at mass-dimension 6 to include 4-fermion operators mixed actions without too much effort. All necessary graphs are obtained using QGRAF [52,53]. The operator insertions are realised by formally introducing additional nonpropagating fields ϕ i called "anchor" and adding i ϕ i O i to the Lagrangian. These anchors correspond to the double lines in Fig. 1. For all n-point functions with a single operator insertion, except single flavour 4-fermion operators, setting options = onepi; in QGRAF is sufficient. The single flavour 4-fermion operators are implemented by splitting the 4-fermion vertex into two 3-vertices connected with an additional "mediator" (here denoted by the dotted line) according to The "mediator" ensures that relative minuses due to anticommutativity of fermions are taken care of 5  This takes into account that cutting the propagator of the "mediator" cannot imply a reducible graph and thus must be included in the set of 1PI graphs. In the same fashion one can also define (q T a q) 2 . Eventually we are interested in the flavour singlets and, for the mixed action, in combinations of flavour singlets in two different flavour subsets. Both can be easily obtained from the results for the flavoured operators.
Composite operators can mix according to the symmetry constraints, e.g. operators compatible with SU[N f ] L × SU[N f ] R flavour symmetry can mix into operators with reduced SU[N f ] V flavour symmetry but not the other way around. This mixing pattern is thus present in the mixing matrix Z as well as in the closely related anomalous dimension matrix Since we are only interested in the leading asymptotic dependence on the lattice spacing, we restrict considerations to the 1-loop coefficient γ O 0 of the anomalous dimension matrix, which has the form Notice that we hid the triangular structure of γ O 0 for the sake of compatibility to the original numbering of the operator basis. The superscripts introduced here indicate the following subsets of operators: Includes operators which are invariant under separate flavour rotations for left-and right-handed quarks.
Operators which are only invariant under joint flavour rotations of left-and right-handed quarks. These operators are in general needed for massless Wilson quarks due to less restrictive flavour symmetry constraints. • m: These operators are lower-dimensional operators multiplied by explicit powers of masses.
For the partially quenched case, where only one lattice discretisation is used for all flavours, it suffices to make the replacements where M f and M b are the diagonal mass matrices for the fermionic and bosonic flavours respectively. We checked explicitly that setting N f = 0 indeed yields the fully quenched approximation also for the 4-fermion operators. This case can be easily obtained by discarding any closed fermion loops. As a cross-check of our renormalisation procedure we also compared our mixing matrix for the non-singlet 4-fermion operators at N f = 3 with the results found in the literature [54,55] -we agree. To do so we made use of the identities Eqs. (12) - (14) to eliminate the redundant (up to Evanescent operators) single flavour 4-fermion operators via ⎛ We stress again that this equivalence does not hold in dimensional regularisation, but for the renormalisation at 1-loop order this subtlety is of no consequence as we can safely ignore any Evanescent operators to this loop order [29]. Via use of Eq. (49) also the special case of Wilson quarks with only one flavour q without quenching can be formally obtained by the change of basis introducing a set of Evanescent operators that mix only into the other operators and not the other way around at N f = 1. Those Evanescent operators take the place of the ones of the form g 2 0 (q i T a q) 2 while we keep g 2 0 (q i q) 2 in our minimal basis. The desired mixing matrix can then be obtained from the subblock without the Evanescent operators. Those renormalised Evanescent operators vanish to all orders in perturbation theory in the limit D → 4, while their mixing contributions are needed to renormalise the remaining operators. While N f = 1 is an uncommon choice, the case of having e.g. N f = 2 + 1 or N f = 3 + 1 with different lattice discretisations for the different flavour subsets is an interesting option, 6 that we want to include here as well for Wilson quarks in the single flavour subset. Notice, that the case of N f = n + 1 without quenching is special in the sense, that if we reached N f = n + 1 by having more than one quark flavour in the second set of flavours but a sufficient number of bosonic flavours, the basis is a different one and we would instead need to include the generalisation of all 4-fermion operators to the partially quenched basis.

Mass-dimension 5
In contrast to the earlier paper [9] we will also include massive operators mixing into the ones without explicit massdependence. Before discussing the operator mixing at massdimension 6 we thus give for completeness the full mixing matrix of the massive mass-dimension 5 operator basis where the corresponding minimal basis O In anticipation of the discussion in Sect. 5.5, we give here the proper basis diagonal under 1-loop renormalisation The coefficient γ B (1) 1 agrees with results found in the literature [56].

Mass-dimension 6
At mass-dimension 6 all 4-fermion operators contribute, which enlarges the mixing matrix quite a bit, and also the number of massive operators grows severely. The diagonalisation of this operator basis is then not feasible for arbitrary values of N f and N as will become clearer in Sect. 5.5. We therefore give here only the subblocks of the non-diagonal mixing matrix. For the massless case we find where we introduced 0 m×n as the m × n matrix filled with zeros. For reuse in the case of mixed actions in Sect. 5.3, we also introduced the short-hands The 13 columns and rows correspond to the non-diagonal massless operator basis O (2) 1 to O (2) 13 in the order they are numbered in Eqs. (11). The subblock A(N f = 0) has already been found for pure gauge theory [9]. The explicit dependence on the argument b is introduced to handle the occurrence of N f due to the renormalisation of the coupling differently than the N f arising from the number of flavours involved in constructing the operator. This will become important for the mixed action case.
For non-vanishing quark masses we additionally get the mixing subblocks Once a specific choice for N f > 0 and N is being made, the diagonalisation can be performed using the Mathematica notebook provided [57]. Also the spectrum of the 1-loop coefficients of the anomalous dimension matrix depends on these values and the symbolic expressions are very complicated if one does not fix N f and N to some numerical values. Therefore we refrain from trying to give here the symbolic results of a diagonal basis altogether and point again to said Mathematica notebook [57].

Generalisation to actions with two flavour subsets
Instead of considering just mixed actions, where one flavour subset is quenched, we will discuss the fully general case of two flavour subsets q and Q using different lattice fermion actions. At mass-dimension 5, this doubles the massless operator basis without introducing any new mixing. However, in the presence of massive operators the minimal basis is severely increased, but still the same three distinct eigenvalues from Eq. (53) must be considered for the spectrum. We thus skip stating the entire mixing matrix at mass-dimension 5 as it does not yield much new insight.
At mass-dimension 6 however, having two distinct flavour subsets requires the inclusion of an additional kind of 4fermion operators already for the massless case, like e.g. for the mixed action in Eq. (37). These operators give rise to new contributions in the spectrum. Allowing for the absence of quenching, we then find for the enlarged massless operator basis in Eq. (37) the nonzero subblocks of the mixing matrix 7 (58d) 7 Using γ The extension to the massive operator basis gives rise to no new values in the spectrum apart from those already found for the non-mixed action case. Of course some values in the spectrum now have a higher degeneracy or, if they are N q,Q f dependent, will split into two different values. We omit stating the full massive mixing matrix to keep the results somewhat compact. Keep in mind that this additional mixing is one-directional such that it does not affect the spectrum for γ i already found in the massless case. Also, the spectrum γ i for the full massive case is known but some matching coefficients of the massive operators are unknown, i.e. some tree-level matching coefficients may vanish thus suppressing those contributions by at least one power inḡ 2 (1/a) further.

Renormalisation of contact term interactions
When considering unimproved Wilson fermions relying on maximal chiral twist to achieve automatic O(a) improvement 8 due to the continuum T 1 -symmetry from Eq. (30), also contact terms from double insertions of the allowed massdimension 5 operators μ 2χ χ and i/4χσ μν F μν χ come into play at O(a 2 ). Since we cannot access the massive mixing in the chirally twisted theory without repeating the entire renormalisation procedure in the twisted theory, we will restrict considerations to the renormalisation of the double insertion of i/4χσ μν F μν χ in the massless case. Working in the massless theory ensures that any chiral twist does not affect the continuum QCD action and we thus can reuse the setup of the untwisted theory. For dimensional reasons, i/4χσ μν F μν χ is anyway the only double insertion that can affect the matching coefficients of the massless mass-dimension 6 operator basis. This contribution remains unchanged when going to the massive theory due to working in a mass-independent renormalisation scheme.
To determine what impact this double insertion has on the matching coefficients, we then need to renormalise the UV divergences arising from the contact interactions. Consequently we also consider the 1PI graphs in Fig. 2, which gives the additional mixing contributions δ Z The appropriate mass-dimension of the mixing operator O 1 since = (4−D)/2 and thus δ Z is dimensionless, which does not allow mixing with operators of another mass-dimension. Of course O (2) k would still contain massive operators of appropriate mass-dimension, if we considered the massive case. Fig. 2 Graphical representation of all considered 1PI graphs when renormalising the UV poles arising from contact interactions of double operator insertions. Each operator is again inserted via an "anchor" field, which is indicated via the double line. Both insertions are at zero momentum q 1 = q 2 = 0 Formally, the resulting additional mixing can be implemented by adding one row (and a trivial column, see Eq. (59)) to the existing mixing matrix. The double insertion will only get mixing contributions from the existing minimal operator basis and not introduce any new mixing beyond that. We restrict considerations to the case, where only one flavour subset contributes a double insertion of i/4Qσ μν F μν Q, e.g. when using twisted Wilson fermions as valence quarks.
For the massless case and two flavour subsets, we then find the mixing contributions where the columns correspond to the same massless operator basis used in Sect. 5.3. As we can see, 4-fermion operators both preserving and breaking chiral symmetry are required to renormalise the contact interactions, such that these operators will have non-vanishing tree-level matching coefficients in the basis in Jordan normal form as introduced in the following subsection.

Renormalisation Group invariants and the asymptotic lattice spacing dependence
For a generic Renormalisation Group invariant (RGI) spectral quantity P, the asymptotic lattice spacing dependence is, see also [9], where δP O i (1/a) is the contribution of the operator O i to the lattice artifacts renormalised at scale μ = 1/a and and thus depends only on the anomalous dimensions of our minimal operator basis. A formal solution to the RGE is where β(ḡ) = −ḡ 3 (b 0 + O(ḡ 2 )) is the β-function and Pexp denotes the path-ordered exponential with increasing x from the right to the left. The path-ordering is needed due to mixing of the operators under renormalisation, i.e. γ O is in general not diagonal. Next we switch the basis O → B such that our 1loop anomalous dimension matrix γ B 0 has Jordan normal form, which takes care of the fact that γ O 0 can be a nondiagonalisable matrix. While for the cases of N f > 0 or N q,Q f > 0 considered here this is not an issue and the Jordan normal form becomes just a diagonal matrix, both quenched and mixed actions can yield a non-diagonalisable 1-loop mixing matrix and thus require special care. To pull out γ B 0 from the path ordered exponential we need to solve another RGE, 9 see [58], with leading order solution W (μ) = 1 + O(ḡ 2 ). This allows to rewrite where we used that γ B 0 is in Jordan normal form and thus can be written as a diagonal matrix plus one matrix containing only one off-diagonal. Both matrices forming γ B 0 commute with each other, such that the path-ordering plays no role here. Eventually we can introduce the Renormalisation Group Invariants 9 We thank Agostino Patella for pointing out the issue with pathordering at subleading orders.
where we introduced The factor 2b 0 in front ofḡ 2 is the common choice for the normalisation. Finally this allows to rewrite where now all scale dependence is absorbed into the prefactor of the RGI quantity. By construction δP B k;RGI is scaleindependent.
In the cases of N f > 0 or N q,Q f > 0 considered here the 1-loop anomalous dimension matrix can be diagonalised such that the remaining exponential in Eqs. (69) and (70) becomes the identity. For both quenched and mixed actions we find terms with additional factors of ln(2b 0ḡ 2 ) modifying the simple power law [ḡ 2 (1/a)]γ . These logs will in general give the dominating contributions for operators with the given leading powerγ in the coupling. However, these particular contributions belong here to chiral symmetry violating 4-fermion operators that are typically suppressed by one power in the coupling as will be discussed in the following section.

Matching to lattice actions
With the full 1-loop mixing matrix computed in the previous section we are now able to compute the leading powers in the coupling modifying the naive a n behaviour. While in general these leading powers are greater or equal tô ii /(2b 0 ), obtained from the 1-loop coefficients of the mixing matrices listed in Sect. 5, we actually may gain additional insight from (tree-level) matching. This means we want to determine the tree-level values of the matching coef- (4) or more importantly its counterpart for the diagonal basisc B i . If a tree-level matching coefficient vanishes for an operator of our diagonalised basis, it shifts the associated leading power in the coupling by (at least) one power in the coupling further. Tree-level matching has the particularly beneficiary feature that a naive expansion of the lattice action in the lattice spacing suffices to obtain the desired tree-level matching coefficients. Beyond tree-level, the perturbative matching procedure requires l-loop compu-tations in perturbation theory for both the lattice theory and continuum SymEFT to achieve matching to l-loops. Especially the lattice side then becomes very complicated due to relying on lattice perturbation theory.
We therefore take only tree-level matching into account, when computing the spectrum of leading powers in the coupling. From tree-level matching we can only infer, whether a leading power in the coupling vanishes and if so assume the next-to-leading order (NLO) to be the truly leading power in the coupling. We then introducê i =γ i + n i (71) as the actual leading power inḡ 2 (1/a), where n i takes (n i − 1)-loop improvement into account, i.e. n i = 1 implies absence of contributions from this operator at tree-level. The various powers ofḡ 2 (1/a) for the different operators contributing to the lattice artifacts may spread over more than one power in the coupling, such that e.g. the subleading contributions of one operator compete with the leading power contributions of another operator.

Tree-level matching for O(a) improved massless action
We find as tree-level matching coefficients for the massless mass-dimension 6 operator basis O wherec O 3 is universal for the action of either Wilson, Overlap or Domain wall quarks as they eventually all employ the same Wilson Dirac operator. The reason for the suppression of most of the 4-fermion operators due to vanishing tree-level matching coefficientsc O 4,5,7−13 is explained in appendix B. Notice that this leads to an additional power inḡ 2 on top of the overall prefactor included in the definition of our 4-fermion operator basis in Eqs. (11d) and (11e). Both tree-level coefficientsc O 1 andc O 2 for the gluonic operators have been taken from [28], where e i are the coefficients for different terms in the lattice action 10 with the conventional normalisation e 0 + 8e 1 + 8e 2 + 16e 3 = 1, namely the plaquette (e 0 ), rectangle (e 1 ), chair (e 2 ) as well as twisted chair (e 3 ) Wilson loops. This general class of pure gauge actions covers a wide range of possible lattice gauge actions. A common and natural choice is to set e 2 = e 3 = 0, which reduces the computational effort as only two kinds of Wilson loops must be computed. It also setsc O 6 = 0, which is interesting as now 10 We renamed those coefficients to e i to avoid a clash with our notation.
only two coefficients remain nonzero, namelȳ A particularly useful choice for the remaining coefficient of the Wilson loops in the lattice gauge action is e 1 = −1/12, which is known as the Lüscher-Weisz action [28]. It ensures tree-level O(a 2 ) improvement of the lattice gauge action such that onlyc O 3 remains nonzero. We already discussed this for lattice pure gauge theory [9]. Having a vanishing tree-level coefficient in the basis O (2) i only ensures overall absence, if the corresponding operator does not mix into another operator that has a non-vanishing tree-level coefficient or mixes itself into another operator and so on as the (tree-level) matching coefficients of the basis in Jordan normal form can be obtained through where T is the change of basis such that the 1-loop mixing matrix is in Jordan normal form, or diagonal as in most cases discussed here. All chiral-symmetry breaking 4fermion operators have vanishing tree-level matching coefficients.
The chiral-symmetry preserving 4-fermion operators can mix into any other massless operator present in our on-shell basis, such that they will generally be present at tree-level, if any other tree-level matching coefficients are nonzero in the basis O (2) i . Knowing the TL matching coefficients has already very interesting consequences. Assuming use of only the plaquette and rectangle term in the gauge action, i.e., the matching coefficients from Eq. (73), we will find in the massless case only the two operators O (2) 2 and O (2) 3 to have nonvanishing tree-level matching coefficients, when considering non-perturbatively O(a) improved Wilson quarks as well as Ginsparg-Wilson quarks. In particular the choice of the treelevel improved Lüscher-Weisz lattice gauge action [28] will ensure that only O (2) 3 has then a non-vanishing tree-level matching coefficient. This allows either to parametrise the leading lattice artifacts in terms of this single coefficient or alternatively to perform Symanzik tree-level improvement at O(a 2 ) by adding only one additional term to the lattice fermion action like e.g.
see also [59,60], where this has already been pointed out. For a more in depth discussion of the (massive) tree-level matching coefficients and their impact on the leading lattice artifacts for O(a) improved Wilson quarks, Ginsparg-Wilson Fig. 3 Leading and subleading powers in the couplingˆ j for the minimal on-shell operator basis of GW and Wilson quarks, describing the leading order lattice artifacts for spectral quantities at N f = 0, 2, 3, 4, 8. N f = 8 is added to highlight what happens when one gets closer to the conformal window. While the solid lines correspond to the massless operator basis (possibly containing massive mixing contributions), the dash-dotted lines belong to the operators with overall mass-dependence. In case the tree-level coefficient vanishes, the first power plotted is regarded as NLO in the colour coding. Due to overlapping numerical values leading powers may be hidden by the "subleading" powers of other operators quarks and Domain-Wall fermions, we refer the reader to [16].

Tree-level matching for chirally twisted Wilson quarks
In case there are double insertions of mass-dimension 5 operators present that break chiral symmetry, as is the case for twisted mass QCD relying on automatic O(a) improvement, additional O(a 2 ) corrections occur. These double insertions then give rise to contact terms in the SymEFT, whose renormalisation affect the tree-level matching coefficients of the 4fermion operators, including chiral-symmetry violating ones, in the basis B (2) i as we can see in Eq. (60). Taking this generalisation into account is achieved in a very sloppy way by enlarging the mixing matrix and thus T in Eq. (74) by one row and column (in the massless case) formed by Eq. (60). Apart from this, the general remarks from the O(a) improved case remain the same. Regarding N f = 8 is added to highlight what happens when one gets closer to the conformal window. Again the solid lines correspond to the overall massless operator basis (possibly containing massive mixing contributions), the dash-dotted lines belong to the operators with overall massdependence. In case the tree-level coefficient vanishes the first power plotted is regarded as NLO in the colour coding. Due to overlapping numerical values leading powers may be hidden by the "subleading" powers of other operators. Notice that we do not know the tree-level matching coefficients for the massive case, which are therefore assumed to be non-vanishing. For the fully quenched case N f = 0, the lowest power may be estimated too low as it could potentially be shifted to the first massless contribution instead the matching for the massive twisted theory, we lack some information as currently only the massless mixing matrix is known. This allows to determineγ (2) i for all operators including the massive ones, but we do not know, whether any tree-level matching coefficients of the massive operators in the basis B (2) i vanish. We thus assume they do not vanish.

Leading powers in the coupling incorporating tree-level matching
Taking all of this into account, we find the leading powers and up to next-to-next-to-next-to-leading orders (N 3 LO) in the coupling, here denoted asˆ i , for the various cases as depicted in Figs. 3, 4 and 5. Going that far into the subleading powers is done only to highlight the overall spread of the leading powers from the various contributions as some of them are severely suppressed by up to three powers inḡ 2 (1/a) compared to the lowest powerˆ min , depending on N f . There also the extension to the massive case is included, where of course the same considerations for tree-level matching can and should be applied to ensure sorting out any vanishing tree-level coefficients. As we can see, the various powers in Here two distinct sets of flavours are assumed either in a mixed action with the setup described in [61] or using different lattice discretisations for light and heavy quarks, here GW and Wilson quarks respectively. While the solid lines correspond to the massless operator basis (possibly containing massive mixing contributions), the dash-dotted lines belong to the operators with overall massdependence. In case the tree-level coefficient vanishes, the first power plotted is regarded as NLO in the colour coding. Due to overlapping numerical values leading powers may be hidden by the "subleading" powers of other operators g 2 (1/a) form a dense spectrum. For use in an ansatz for the continuum extrapolation we give in Table 1 Fig. 4 with the case of Wilson in Fig. 3 at even number of flavours. Evidently, the spectrum of leading powers in the coupling is denser without having explicit improvement.
As was probably to be expected, introducing two sets of flavours, either dynamical or quenched, introduces even more operators and thus severely increases the density of the spectrum compared to having only one lattice discretisation, see e.g. Figs. 3 and 5. For the cases considered here, this will also impact the leading power in the coupling for the massless case. While e.g. at N f = 3 we only find a slight decrease fromˆ min ≈ 0.247 toˆ min ≈ 0.198, this also abandons the already weakly pronounced gap of the 1st to 2nd leading power from ˆ ≈ 0.42 to ˆ ≈ 0.05. This effect holds true for the other values for N f considered here. Notably, the case of a mixed action with two valence quarks considered here yields a non-diagonalisable mixing matrix at 1-loop order. This gives rise to an additional factor log(2b 0ḡ 2 (1/a)) modifying the pure power law in the coupling forˆ ≈ 1.586, which is however very suppressed compared to the leading powers inḡ 2 (1/a), namelyˆ min ≈ 0.230 in the massless case andˆ min ≈ −0.172 in the massive case. The spectra of mixed action and light GW + heavy Wilson quarks in Fig. 5 are very similar, but a few important differences occur for the latter case. Namely the use of Fierz identities for N f = 2 + 1 or N f = 3+1 reduces the number of 4-fermion operators and having light GW quarks also excludes some chiral-symmetry violating 4-fermion operators as explained after Eq. (37).
In general, the additional contributions to the spectrum for the massive case result again in a slightly denser spectrum, but do not affect the spectrum of the massless operator basis due to the tridiagonal form of the mixing matrix, see Eq. (47). In all cases considered, except twisted Wilson quarks purely relying on automatic O(a) improvement, the lowest power Table 1 Non-exhaustive examples for leading powers in the coupling relevant for the leading order lattice artifacts a 2 [ḡ 2 (1/a)]ˆ at N = 3 and various values for N f . Degeneracies are not listed as they cannot be distinguished numerically. While vanishing matching coefficients have been taken into account by shifting the corresponding leading power by 1, hierarchies between the various coefficients are ignored and certainly should be looked at [16]. We limit considerations to commonly used combinations and some interesting cases. For the gauge action we assume use of a variant compatible to Eq. (73). The underlined numbers belong to massive operators. If the Lüscher-Weisz action [28] is being used rather than the more general case in Eq. (73) the leading powers do not change, since the operator affected by this change of matching is already sufficiently suppressed byγ i 0.86 for 2 ≤ N f ≤ 8. The leading powers for O(a) improved tmQCD should be identical to the Wilson case, but massive tree-level matching coefficients may differ and are currently unavailable. We therefore assume those coefficients to be non-vanishing found at N f = 2, 3, 4 is always from the massive operator basis and actually slightly negative. At N f = 8 the lowest power is from the massless operator basis and negative as well. Without counting double operator insertions, there are eleven distinct massive operators at mass-dimension 6, but only three distinct powersˆ i . Taking the subleading powers into account, the spectrum of the massive operators modulo 1 has only two distinct powers. This is due to a severe degeneracy of the spectrum for these operators since additional powers of the mass only shift the overall 1-loop anomalous dimension by a constant. Introducing two sets of flavours does not affect the powers for the massive contributions as they do not depend on N Of course the degeneracy in the spectrum grows even further.
From a numerical point of view this degeneracy suggests to treat the different operators contributing with the same power as one linear combination, because one cannot distinguish them during a fit anyway (actually most of the degeneracy in the power spectrum starts at subleading orders due to vanishing tree-level matching coefficients of some of the operators). For tmQCD relying only on automatic O(a) improvement, the spectrum for the massive case has two additional values, where one of those again increases the degeneracy. Introducing the Wilson clover term with a non-perturbative improvement coefficient eliminates all truly new powers in the coupling and avoids any contact term renormalisation of our operator basis, which otherwise would affect O(a 2 ) just like in the untwisted theory. In particular, this ensures vanishing of tree-level matching coefficients for the chiral symmetry violating 4-fermion operators. Only the double insertion of one massive operator remains, that has a degenerate power in the coupling and is thus indistinguishable from the contributions of massive mass-dimension 6 operators.

Discussion
We have computed the 1-loop anomalous dimensions of all mass-dimension 5 and 6 operators relevant for the minimal basis describing lattice artifacts from either the Wilson or Ginsparg-Wilson action up to and including O(a 2 ). This also includes all relevant massive operators as well as the necessary generalisation to the quenched case, mixed actions or use of two distinct lattice discretisation for different quark flavours. These 1-loop coefficients of the anomalous dimension matrix modify the leading asymptotic lattice spacing dependence from the classical integer power a n behaviour to a n [ḡ 2 (1/a)]ˆ (n) i . Only for the quenched case and mixed actions considered here additional factors of log(ḡ 2 (1/a)) occur, see also Sect. 5.5. Hereˆ i +n i takes suppressions from matching coefficients byḡ 2n i (1/a) into account. In practice we only know whether the tree-level matching coefficients vanish for the lattice actions considered here except for the massive operators of tmQCD, where the full mixing with massive operators is currently not available. We explicitly write here the additional superscript inˆ (n) i , which has been dropped until now, to include also O(a) corrections, i.e. n = 1. We have discussed those effects already before [9] without taking the massive operators into account. We only mention those effects briefly.
For example at N f = 3 we found the lowest values of the massless spectrum to beˆ (1) min ≈ 0.074 for the O(a) lattice artifacts of unimproved Wilson quarks andˆ (2) min ≈ 0.25 for the O(a 2 ) lattice artifacts of both non-perturbatively O(a) improved Wilson quarks and Ginsparg-Wilson quarks. For typical numbers of flavours N f = 0, 2, 3, 4 all powers in the coupling from the massless operator basis improve the convergence towards the continuum limit as a 0 due to being positive. Only maximally twisted Wilson quarks without explicit O(a) improvement have e.g. at N f = 2 one slightly negative power in the couplingˆ (2) min ≈ −0.122 for the massless operator basis arising from the contribution of 4-fermion operators that violate chiral symmetry. Those contributions become worse for increasing number of flavours. In contrast to the non-perturbatively O(a) improved lattice actions those operators can arise here already at tree-level matching due to the renormalisation of contact terms from double insertions of the mass-dimension 5 operator basis in the SymEFT.
The example of twisted Wilson quarks with automatic O(a) improvement actually teaches us another lesson, namely that relying on continuum symmetries for O(a) improvement may worsen the approach to the continuum limit at O(a 2 ) and beyond compared to explicit Symanzik improvement. Performing instead explicit O(a) improvement -already treelevel improvement would suffice for that matter -ensures that double operator insertions of mass-dimension 5 operators do not occur or at least occur suppressed by at least two additional powers inḡ 2 (1/a), which then also shifts any effects from contact term renormalisation by this additional power in the coupling. Fortunately the ETMC collaboration already performs explicit O(a) improvement [31]. Under this impression the earlier attempts of obtaining a classically perfect action, see e.g. [62], may have been a worthwhile endeavour, since this argument of course occurs at any order in the lattice spacing. If there are negative powers in the coupling present at subleading powers in the lattice spacing, these will automatically be shifted up by at least one power in the coupling as will allˆ (n) i from tree-level contributions. The additional massive operators must be considered, when one is working in a mass-independent renormalisation scheme rather than a hadronic scheme. Overall they increase the density of the spectrum found forˆ (2) i even further, but introduce only three distinct powers due to a fairly degenerate spectrum. In all cases considered here apart from twisted-mass QCD relying on automatic O(a) improvement, the massive operators decrease the leading power in the coupling at N f = 2, 3, 4 slightly compared to the massless basis.
In either case all the leading powers encountered are much better behaving than the dominant power min i (ˆ (2) i ) = −3 found in the 2-d O(3) sigma model [7,8], which is good news. However, due to the large number of operators we get many different contributionsˆ (2) i to the spectrum, which in addition lie in close proximity to one another, see also Figs. 3, 4 and 5. This will make it very difficult to decide, which contributions actually dominate and therefore must be included in a proper fit-ansatz for continuum extrapolations. The density in the spectrum forˆ (2) i becomes even worse when introducing two distinct lattice discretisations, like we discussed here in form of mixed actions or different discretisations for dynamical heavy and light quark flavours. Eventually the situation is even more complicated since the sensitivity to contributions from different operators of our minimal basis will presumably depend on the (spectral) quantity to be extrapolated and also different orders of magnitude of the various LO matching coefficients will have an impact [16]. Also, there are of course always corrections subleading in the power of the lattice spacing that one must be wary about. All of this indicates that one should be very careful when doing continuum extrapolations regarding the systematic errors associated to the extracted continuum values.
Before discussing possible extensions of this work we should ask ourselves whether today's lattice simulations are at sufficiently small lattice spacings such that a perturbative description of the leading asymptotic lattice spacing dependence suffices. While there will never be absolute certainty that the leading power in the lattice spacing really dominates the picture, we may at least have a look at the running couplings α MS (1/a) = g 2 MS (1/a)/(4π) associated with the lattice spacings typically available, using perturbative 5-loop running [50]. A non-exhaustive overview of lattice spacings available in the literature is given in Table 2, where one finds that at least the smallest lattice spacings clearly start to reach the perturbative region. 11 Of course, even smaller lattice spacings would (always) be preferable, but may be too costly at the moment. Meanwhile step-scaling analysis in pure gauge theory [67] reach down to lattice spacings as small as a < 10 −3 fm. There, the absolute systematic effect on the extracted continuum value should become less severe when using the classical a n power law instead of including the logarithmic corrections discussed here.
In case one is interested in non-spectral quantities each field introduces an additional minimal basis of higher dimensional operators compatible with the transformation Table 2 Non-exhaustive overview of typical lattice spacings in today's lattice QCD simulations combined with the MS coupling at scale μ = 1/a obtained from 5-loop running [50]. The MILC HISQ ensem-bles are only added for comparison as they involve staggered quarks, which were not included in our discussion behaviour of the local field. The associated spectrum must then be taken into account as well. An analysis of such additional powers in the coupling for fermion bilinears is on its way. In general such computations aiming only at 1-loop anomalous dimensions should not be too complicated, but finding the minimal basis is somewhat tedious. As mentioned before, once the minimal basis is found, the strategy described here should be applicable also for the local fields with the sole difference that the operators must be renormalised at non-zero momentum or in other words mixing with total divergence operators must be taken into account as well. In particular for precision physics observables, performing these computations should become the norm as it offers better control over the continuum extrapolation by either reducing the error budget or providing a better estimate for this uncertainty. A more complicated issue is the generalisation to staggered quarks [68], which requires a whole new class of operators that allow for flavour changing interactions, see e.g. [69], and thus have a severely reduced flavour symmetry compared to Ginsparg-Wilson quarks, whose operators of course still contribute with the known values forˆ (2) i . For the additional set of operators we do not know whether they fall within the same range of valuesˆ (2) i or possibly undershoot them. The latter should of course give rise to concerns as there is no theoretical lower bound. Due to the prominence in the literature and low computational cost of this lattice discretisation in numerical simulations, this is probably the most pressing gap in this work. However, this is not the only theoretical concern arising when using staggered quarks. Firstly, there is no proof of perturbative renormalisability of staggered quarks to all orders in perturbation theory available, using the lattice power counting theorem generalised for staggered quarks [70]. Such a proof exists for Wilson [71] and Ginsparg-Wilson fermions [72]. Secondly, so called rooting to reduce the number of flavours, makes rooted staggered quarks a highly non-local theory, see e.g. [73]. Both aspects may invalidate the applicability of Symanzik Effective theory altogether. Nonetheless a Symanzik Effective theory analysis can still be done keeping in mind that both issues mentioned here must be resolved independently to have a solid theoretical basis.

Supplementary information
A Mathematica notebook is included in the supplementary material (matchingMixedContact.nb). It allows to obtain the leading powers for mixed actions with massless quarks for other choices of number of flavours.
Acknowledgements I am indebted to Rainer Sommer and Peter Marquard for supervising my PhD project on which this work is based on and for many helpful discussions regarding both the project and this paper. I also thank Agostino Patella for pointing out an issue with our use of the path-ordered exponential when defining RGI quantities, that would lead to wrong NLO corrections in the coupling. Finally, I thank Hubert Simma and Kay Schönwald for many discussions as well as the SHEP group at the University of Southampton to allow me some time to finish this project. The Feynman diagrams used in this paper have been generated with help of the LaTeX package TikZ-Feynman [74].

Data Availability Statement
This manuscript has data included as electronic supplementary material. The online version of this article contains supplementary material, which is available to authorized users.

Conflict of interest
The author declares that he has no financial or non-financial interests that are directly or indirectly related to the work submitted for publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.
Equivalently, we could have considered the aforementioned connected graphs of fundamental fields with some gauge-fixing and proper choices for the momenta of the external legs, i.e., analytical continuation to Minkowski space and on-shell momenta. In either case, working at tree-level ensures that we can use the naively-a-expanded lattice action to derive the relevant Feynman rules 12 for vertices and propagators of the lattice theory. This reassures us that we could have stopped the tree-level matching immediately after the naive a-expansion of the lattice action. Beyond tree-level, proper use of LPT is required.