text:
# Non-perturbative Renormalisation of Domain Wall Fermions: Quark Bilinears

## I Introduction

Renormalisation of lattice operators is an essential ingredient needed to deduce physical results from numerical simulations. In contrast with the determination of hadronic masses, physical matrix elements can be determined only if the normalisation of the appropriate lattice operators can be related to that of the corresponding continuum operators, conventionally specified perturbatively at short distances. In principle, lattice perturbation theory may be used to establish this connection. However, lattice perturbation theory converges slowly and the expansion parameter, the square of the lattice coupling evaluated at the lattice scale, $`g(a)^2`$, decreases only as an inverse power of $`\mathrm{ln}(a)`$. This makes systematic improvement of perturbative results essentially impossible. This convergence may be improved when, following ideas from continuum perturbation theory , a renormalised or ”boosted”  coupling rather than the bare coupling is used as an expansion parameter. Even so, considerable arbitrariness remains, and in general it is extremely difficult to go beyond one loop order in such calculations. To overcome these difficulties, Martinelli $`\mathrm{𝑒𝑡}.\mathrm{𝑎𝑙}.`$  have proposed a promising non-perturbative renormalisation procedure. This method has been previously used to determine renormalisation coefficients for various operators using the Wilson  and staggered actions . The purpose of this work is to study the application of this technique to the renormalisation of the quark field and flavour non-singlet fermion bilinear operators for the domain wall fermion action.

Domain wall fermions  provide an action, that at the expense of introducing a fifth dimension, has a low energy theory with excellent chiral properties while at the same time preserving exact flavour symmetry. These good chiral properties lead to a suppression of the possible dimension five terms in the long-distance effective Lagrangian implying that domain wall fermions define a lattice version of QCD which is off-shell improved to $`𝒪(a^2)`$. As we will see, these domain wall off-shell Greens functions show remarkably reduced lattice artifacts. A study of operator renormalisation coefficients for this action is useful, both because these numbers are needed for use in practical calculations of physical quantities  and because it provides an excellent test of the chiral properties of the domain wall fermion action in practical simulations. In fact, we find that domain wall fermions perform quite well for non-perturbative renormalisation with negligible contributions from explicit chiral symmetry breaking. This finding is in good agreement with recent work on the chiral limit of quenched QCD with domain wall fermions .

Careful operator normalisation is especially important for the domain wall fermion method. As is reviewed in Section III A, the interpolating field conventionally used to create and destroy the physical modes is exactly localised in the fifth-dimension on the right and left walls. Since the actual physical modes extend somewhat into the fifth dimension, the overlap between the interpolating field and the physical modes will be smaller than one. This implies a wave function renormalisation factor $`(Z_q)`$ which differs from one even in the case of free fields. For the eigenvectors corresponding to the smallest 19 Dirac eigenvalues examined in the quenched, $`\beta =6.0`$ calculation of Ref. , this overlap typically varies between 75 and 85%. Fortunately, the non-perturbative methods employed here  precisely include these effects.

We begin in Section II with a brief summary of the main issues involved in applying the non-perturbative renormalisation method. In Section III, we give the domain wall fermion action and discuss the Ward-Takahashi identities it obeys. Section IV builds on this base to constrain the ways in which explicit chiral symmetry breaking terms may enter low energy matrix elements calculated using domain wall fermions. In Section V we give the details of our lattice simulations. Section VI describes the renormalisation of the quark propagator, and in Section VII we introduce the quark bilinears and compute their renormalisation on the lattice in the regularisation independent scheme. After removing expected non-perturbative pole terms, we look for effects of explicit chiral symmetry breaking and find that they are negligible. In Section VIII, we avail ourselves of the axial Ward-Takahashi identity again to compute the quark wave function renormalisation from the conserved vector and partially conserved axial-vector currents. In Section IX we calculate the renormalisation of the non-conserved, local axial current from a ratio of its hadronic matrix element to the hadronic matrix element of the partially conserved axial current and find good agreement with the results of Section VII. In Section X we convert the renormalisation coefficients to renormalisation group invariant quantities by dividing out the renormalisation group running. In Section XI we discuss the calculation of the quark wavefunction renormalisation from the propagator.

After comparing our non-perturbative results with recent perturbative calculations in Section XII, we end with our conclusions. The details of the exact conventions and equations used for the perturbative running and matching are relegated to appendices.

## II Non-Perturbative Renormalisation

In the following the method of non-perturbative renormalisation introduced in Ref.  will be studied. This method uses a renormalisation scheme that is defined by a set of conditions that mandate the renormalised values of the operators of interest between external quark states, in a fixed gauge, at large virtualities. As such these conditions may be expressed in any regularisation scheme (and so this scheme is known as the regularisation independent (RI) scheme). In particular this allows the renormalisation factors to be defined in the lattice regularisation, opening the way for renormalisation factors to be directly calculated in numerical lattice simulations.

While calculating renormalisation factors from lattice simulations neatly avoids the need to perform analytic calculations using lattice perturbation theory, which are both challenging and poorly behaved, doing so introduces several issues that must be considered:

* Calculating the matrix elements of the operators of interest between external quark states requires a fixed gauge to be used. This allows for the appearance of Gribov copies, possibly obscuring the required comparison with continuum perturbation theory where only the trivial copy appears. Earlier studies of the size of Gribov noise in the calculation of a gauge invariant normalisation factor as a ratio of two gauge-variant amplitudes suggest this may not be an important difficulty for the parameters used here. However, in future work, this difficulty can be avoided by taking two steps: i) Impose the regularisation invariant normalisation condition in a sufficiently small physical volume that non-perturbative effects are suppressed. ii) Begin the Landau gauge fixing procedure from a configuration that is in a completely-fixed axial gauge. Taking these two precautions will insure that any effects of Gribov copies will be similar to other non-perturbative effects and will vanish as the comparison with perturbation theory is done at weaker and weaker coupling.
* Numerical simulations are performed with a finite lattice spacing. This provides a natural condition,

$$|p|\frac{1}{a},$$
(1)
over the momenta range for which a direct extraction of continuum quantities is possible.
* As the renormalisation factors are determined in a non-perturbative calculation the contributions of propagating mesons, and in particular pseudo-Goldstone bosons, must be identified and removed. These effects may be reduced by working at high momenta, with a natural condition for the absence of significant deviations being

$$\mathrm{\Lambda }_{QCD}|p|.$$
(2)

Taking the last two points together suggests that this technique relies on the existence of a “window” of momenta,

$$\mathrm{\Lambda }_{QCD}|p|\frac{1}{a},$$
(3)
for which the predictions of continuum perturbation theory should correctly describe the form of the lattice data. In practical simulations however, it has been found that the effects of deviations due the violations of both these inequalities must be taken into account .

Fortunately, near either edge of this window, the form of deviations from perturbative behaviour may be predicted. In the case of momenta too low, the initial corrections may be described by an expansion in terms of momentum-suppressed condensate terms by use of the operator product expansion (OPE). In turn, the first corrections to continuum-like behaviour may be taken into account in terms of an expansion in the lattice spacing, $`a`$.

Another trivial consequence of the restricted range of momenta available in current lattice simulations is the need for many phenomenological calculations to be composed of continuum perturbation theory calculations at high scales, that are then run down to scales accessible on the lattice and combined with the lattice result. As the majority of the existing calculations for the continuum perturbative results use renormalisation schemes that may only be defined when using dimensional regularisation (such as the $`\overline{MS}`$ scheme), perturbative matching calculations between these schemes and the ones that may be defined in the lattice regularisation need to be performed.

## III Domain Wall Fermions

In this section the domain wall fermion formulation, as used in our simulations, will be reviewed.

### A Action

The domain wall fermion (DWF) method is a promising new approach to lattice QCD introduced in Ref. , which, at the expense of introducing an extra, discrete, non-gauge dimension, provides drastically improved chiral properties at finite lattice spacing while preserving exact symmetry under vectorial flavour rotations. This is achieved by using an action in the fifth dimension that is asymmetric between the left-handed and right-handed components of the fermion field. Denoting the fifth co-ordinate as $`s`$, with

$$s0,\mathrm{},L_s1,$$
(4)
the massless action may be written as

$`S_{\mathrm{fermion}}(m_f=0)`$ $`=`$ $`{\displaystyle \underset{x,s}{}}\overline{\mathrm{\Psi }}_{x,s}\{\gamma _\mu {\displaystyle \frac{1}{2}}(_\mu ^++_\mu ^{})`$ (6)
$`+[{\displaystyle \frac{1}{2}}_\mu ^{}_\mu ^++M_5]+P_L_5^+P_R_5^{}\}\mathrm{\Psi }_{x,s},`$
with

$$Z=[dU][d\overline{\mathrm{\Psi }}d\mathrm{\Psi }]\mathrm{exp}\left(S_{\mathrm{gauge}}S_{\mathrm{fermion}}\right).$$
(7)
In Eqs. 6 and 7 $`\mathrm{\Psi }_{x,s}`$ is the fermionic field, $`U_\mu (x)`$ is the gauge field and

$$S_{\mathrm{gauge}}=\beta \underset{P}{}\left(1\frac{1}{3}\mathrm{Re}\mathrm{Tr}\left[U_P\right]\right),$$
(8)
with $`\beta =6/g_0^2`$ and $`g_0`$ is the bare lattice coupling. The projectors for the left and right-handed spinors are defined as,

$`P_L`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1\gamma ^5\right)`$ (9)
$`P_R`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1+\gamma ^5\right).`$ (10)
The notation $`_\mu ^\pm `$ has been used to denote the discrete forward/backward covariant derivatives:

$`_\mu ^+\psi _x`$ $`=`$ $`\left[U_\mu (x)\psi _{x+\mu }\psi _x\right]`$ (11)
$`_\mu ^{}\psi _x`$ $`=`$ $`\left[\psi _xU_\mu ^{}(x\mu )\psi _{x\mu }\right],`$ (12)
and $`_\mu ^\pm `$ represents the corresponding derivative with no gauge term. For the case of the derivative in the fifth dimension, $`_5^\pm `$, the domain wall is implemented by giving the derivative hard boundaries. For example a one-dimensional $`_5^+`$ acting on a space with four points may be written in matrix form as

$$_5^+\left(\begin{array}{c}\mathrm{\Psi }_{x,0}\\ \mathrm{\Psi }_{x,1}\\ \mathrm{\Psi }_{x,2}\\ \mathrm{\Psi }_{x,3}\end{array}\right)=\left(\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 1& 0\\ 0& 0& 1& 1\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}\mathrm{\Psi }_{x,0}\\ \mathrm{\Psi }_{x,1}\\ \mathrm{\Psi }_{x,2}\\ \mathrm{\Psi }_{x,3}\end{array}\right).$$
(13)
It should be noted that the action in Eq. 6 is actually the hermitian conjugate of the action proposed in Ref. . This change was made for practical reasons related to compatibility with the existing Wilson operator implementation for the QCDSP machine.

In the free theory, for $`0<M_5<2`$, the effect of this is to produce a spectrum with one light fermionic mode, with exact chiral symmetry in the $`L_s\mathrm{}`$ limit, and $`16L_s1`$ heavy modes. The wavefunction of this light mode has its right-handed component concentrated on the wall at $`s=L_s1`$ and its left-handed component on the wall at $`s=0`$. This light fermion mode may be studied by introducing an interpolating operator of the form

$`q_x`$ $`=`$ $`P_L\mathrm{\Psi }_{x,0}+P_R\mathrm{\Psi }_{x,L_s1}`$ (14)
$`\overline{q}_x`$ $`=`$ $`\overline{\mathrm{\Psi }}_{x,0}P_R+\overline{\mathrm{\Psi }}_{x,L_s1}P_L.`$ (15)
The above considerations also naturally lead to the introduction of an explicit mass term to the action of the form

$$S_{\mathrm{fermion}}(m_f)=S_{\mathrm{fermion}}(m_f=0)+\underset{x}{}m_f\overline{q}q,$$
(16)
where $`m_f`$ is the bare quark mass. In the free case, this leads to a spectrum with one light fermion of mass

$$M_5\left(2M_5\right)\left(m_f+\left(1M_5\right)^{L_s}\right).$$
(17)
Note that in the $`L_s\mathrm{}`$ limit this is proportional to $`m_f`$, while for finite $`L_s`$ there remains a residual mass, $`m_{\mathrm{res}}`$, that acts as an additive renormalisation to $`m_f`$.

However, the properties of domain wall fermions in the presence of gauge fields is a much more difficult question. In particular while the form of the mass of the light mode is expected to be proportional to $`m_f+m_{\mathrm{res}}`$, the dependence of $`m_{\mathrm{res}}`$ on $`L_s`$ must be determined. Perturbative calculations  have shown that the existence of the light mode is stable to small perturbations and that this mode has all chiral symmetry breaking proportional to $`m_f`$ as $`L_s\mathrm{}`$. These studies also highlight several issues that must be considered when undertaking numerical simulations:

1. The dependence of $`m_{\mathrm{res}}`$ on $`L_s`$ may no longer be of the simple exponential form shown in Eq. 17.
2. $`M_5`$ undergoes a strong additive renormalisation. This is understandable, as the five dimensional problem has no approximate chiral symmetry to protect it.

Indeed, extensive numerical studies in the quenched approximation  have shown that the $`L_s`$ dependence of $`m_{\mathrm{res}}`$ does not fit a single exponential in the range $`L_s=1248`$ for lattices with the same lattice spacing ($`a=0.520\mathrm{GeV}^1`$) as the results in this paper. For $`L_s=16`$, the value used in this work, $`m_{\mathrm{res}}`$ was found to be $`4\mathrm{M}\mathrm{e}\mathrm{V}`$ in the $`\overline{MS}`$ scheme at $`2\mathrm{G}\mathrm{e}\mathrm{V}`$ . The strong additive renormalisation of $`M_5`$ requires that an input value be chosen numerically so that a single light mode forms and that its decay in $`L_s`$ is as rapid as possible. It has been found that for even coarser lattices than used here such a choice can be made .

### B Lattice Ward-Takahashi Identities

For the purpose of analysing the consequences of the symmetries of the action, it is convenient to introduce an extended mass term, $`M`$, with flavour structure such that the mass term reads

$$\overline{q}_LM^{}q_R+\overline{q}_RMq_L,$$
(18)
and so the mass term is invariant under a transformation of the quark fields and the mass matrix $`M`$ of the form

$`q_L`$ $``$ $`U_Lq_L`$ (19)
$`q_R`$ $``$ $`U_Rq_R`$ (20)
$`M`$ $``$ $`U_RMU_L^{}.`$ (21)
Following Ref. , on a finite lattice, an exact vector Ward-Takahashi identity may be derived by considering transformations of the 5-dimensional fermion field, $`\mathrm{\Psi }`$, such that

$`\delta _V\mathrm{\Psi }_{x,s}`$ $`=`$ $`iϵ_x^aT^a\mathrm{\Psi }_{x,s}`$ (22)
$`\delta _V\overline{\mathrm{\Psi }}_{x,s}`$ $`=`$ $`iϵ_x^a\overline{\mathrm{\Psi }}_{x,s}T^a,`$ (23)
where $`\left\{T^a\right\}`$ is the set of hermitian traceless matrices acting on $`SU(N_f)`$ flavour-space. This leads to an exact Ward-Takahashi identity that reads:

$`_\mu ^{}𝒱_\mu ^a(z)O(x_1,\mathrm{},x_n)`$ (24)
$`+\overline{q}[M,T^a]q(z)O(x_1,\mathrm{},x_n)`$ $`=`$ $`i\delta ^aO(x_1,\mathrm{},x_n),`$ (25)
where

$`𝒱_\mu ^a(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{s}{}}(\overline{\mathrm{\Psi }}_{x+\mu ,s}(1+\gamma _\mu )U_{x,\mu }^{}T^a\mathrm{\Psi }_{x,s}`$ (27)
$`\overline{\mathrm{\Psi }}_{x,s}(1\gamma _\mu )U_{x,\mu }T^a\mathrm{\Psi }_{x+\mu ,s}).`$
For the case of axial transformations the analogous choice is a transformation of the form:

$`\delta _A\mathrm{\Psi }_{x,s}`$ $`=`$ $`iϵ_{x,s}^aT^a\mathrm{\Psi }_{x,s}`$ (28)
$`\delta _A\overline{\mathrm{\Psi }}_{x,s}`$ $`=`$ $`iϵ_{x,s}^a\overline{\mathrm{\Psi }}_{x,s}T^a,`$ (29)
with

$$ϵ_{x,s}^a=\{\begin{array}{cc}\hfill ϵ_x^a& 0s<L_s/2\\ \hfill ϵ_x^a& L_s/2s<L_s.\end{array}$$
(30)
This leads to a Ward-Takahashi identity of the form

$`_\mu ^{}𝒜_\mu ^a(z)O(x_1,\mathrm{},x_n)`$ (31)
$`+\overline{q}\{M,T^a\}\gamma _5q(z)O(x_1,\mathrm{},x_n)`$ (32)
$`+2J_{5q}^aO(x_1,\mathrm{},x_n)`$ $`=`$ $`i\delta _A^aO(x_1,\mathrm{},x_n),`$ (33)
where

$`𝒜_\mu ^a`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{s}{}}\mathrm{sign}(s{\displaystyle \frac{L_s1}{2}})(\overline{\mathrm{\Psi }}_{x+\mu ,s}(1+\gamma _\mu )U_{x,\mu }^{}T^a\mathrm{\Psi }_{x,s}`$ (35)
$`\overline{\mathrm{\Psi }}_{x,s}(1\gamma _\mu )U_{x,\mu }T^a\mathrm{\Psi }_{x+\mu ,s})`$
$`J_{5q}^a`$ $`=`$ $`\overline{\mathrm{\Psi }}_{x,L_s/21}P_LT^a\mathrm{\Psi }_{x,L_s/2}+\overline{\mathrm{\Psi }}_{x,L_s/2}P_RT^a\mathrm{\Psi }_{x,L_s/21}.`$ (36)
Therefore, in contrast to the previous case, the axial current is not exactly conserved. This is necessary both to provide a mechanism for physical terms due to the $`U(1)_A`$ axial anomaly to enter the calculated amplitudes and also to allow for explicit chiral symmetry breaking contributions at finite $`L_s`$. The situation is analogous to that for Wilson fermions , where the role of $`J_{5q}^a`$ is played by the chiral variation of the Wilson term, except that the contributions from $`J_{5q}^a`$ are expected to tend to zero as $`L_s\mathrm{}`$ in the present case . The form of the contributions from $`J_{5q}^a`$ will be further discussed in the next section.

## IV Operator Mixing and Chiral Symmetry

The major attraction of the domain wall fermion formalism is its ability to decrease the size of chiral symmetry breaking by increasing the parameter $`L_s`$, the distance between the two four-dimensional lattice boundaries to which the left and right chiral modes are bound. However, it is often impractical or inefficient to choose such a large value of $`L_s`$ that all chiral symmetry breaking effects from mixing between these walls can be neglected. Thus, it is important to characterise the effects of this chiral symmetry breaking and in this section we will determine how it can effect the low energy physics of lattice QCD. As we will see, this can be done as either an expansion in the size of the wall-mixing effects, which for simplicity we will denote by $`𝒪(e^{\alpha L_s})`$ although the exact $`L_s`$ dependence may be different, and/or as an expansion in the lattice spacing $`a`$.

This analysis is easily made by starting with the interpretation of chiral symmetry proposed by Furman and Shamir . Here one identifies the full $`SU(N_f)_LSU(N_f)_R`$ chiral symmetry of the continuum theory as the independent $`SU(N_f)`$ rotation of the fermion fields defined on the left- and right-hand halves of the five-dimensional lattice:

$$\begin{array}{cccc}\hfill \mathrm{\Psi }(x,s)& & U_L\mathrm{\Psi }(x,s)\hfill & 0sL_s/21\hfill \\ \hfill \mathrm{\Psi }(x,s)& & U_R\mathrm{\Psi }(x,s)\hfill & L_s/2sL_s1,\hfill \end{array}$$
(37)
where $`U_L`$ and $`U_R`$ are $`N_f\times N_f`$ special unitary matrices belonging to the left and right factors of $`SU(N_f)_LSU(N_f)_R`$. From Eq. 14 it is clear that this transformation will act on the four-dimensional quark fields as a standard element of the full chiral symmetry.

Of course, the transformation in Eq. 37, whose generators are given in Eq. 29 and Eq. 23, cannot be an exact symmetry of the five-dimensional theory as the derivative terms in the fifth dimension, taken collectively, couple the left and right hand walls and prevent such independent rotations of this single, five-dimensional field. However, in the low energy sector of the theory this symmetry can be quite good. The physical, chiral modes which survive at low energy are expected to be exponentially bound to the walls with an overlap that is suppressed as $`L_s`$ increases. The higher energy modes which can propagate freely between the walls are all far off-shell with propagators which are necessarily also exponentially suppressed at long distances, especially for the large distance $`L_s`$.

In order to characterise the effects of this controlled symmetry breaking that comes from communication between the walls, we will generalise somewhat the Dirac domain wall fermion operator of Shamir given in Eq. 16. We will introduce a special-unitary, flavour matrix $`\mathrm{\Omega }`$ in the derivative term joining the four-dimensional planes $`s=L_s/21`$ and $`s=L_s/2`$. Thus, we will modify Eq. 6 by adding the term:

$$S_\mathrm{\Omega }=\underset{x}{}\left\{\overline{\mathrm{\Psi }}_{x,L_s/21}P_L\left(\mathrm{\Omega }^{}1\right)\mathrm{\Psi }_{x,L_s/2}+\overline{\mathrm{\Psi }}_{x,L_s/2}P_R\left(\mathrm{\Omega }1\right)\mathrm{\Psi }_{x,L_s/21}\right\}.$$
(38)
If we include the transformation of the matrix $`\mathrm{\Omega }`$,

$$\mathrm{\Omega }U_R\mathrm{\Omega }U_L^{},$$
(39)
this generalised domain wall Dirac operator will now possess exact chiral symmetry. Note, a comparison with Eq. 21 shows that $`\mathrm{\Omega }`$ transforms “like a mass term”.

Thus, if we examine this generalised theory that includes the chiral matrix $`\mathrm{\Omega }`$, all amplitudes will become functions of $`\mathrm{\Omega }`$ but will exactly obey the chiral symmetry described by Eq. 37 and Eq. 39. Therefore, we need only understand how the matrix $`\mathrm{\Omega }`$ will enter the low energy Green’s functions of interest to determine in a precise way the transformation properties of the chiral symmetry breaking induced by mixing between the walls.

To zeroth order in $`e^{\alpha L_s}`$, the fermion degrees of freedom will remain bound to the walls and propagation from one wall to the other can be neglected. In such circumstances, the matrix $`\mathrm{\Omega }`$ which is introduced at a point mid-way between the walls cannot enter, the amplitude will be independent of $`\mathrm{\Omega }`$ and hence naively invariant under the full $`SU(N_f)_LSU(N_f)_R`$ chiral symmetry. To the next order, $`e^{\alpha L_s}`$, we expect phenomena which involve a single propagation between $`s=0`$ and $`s=L_s1`$. Thus, the matrix $`\mathrm{\Omega }`$ should enter linearly in such amplitudes.

An important application of this analysis is to constrain the form of the effective continuum action which gives amplitudes that agree with those of the domain wall theory through a given order in the lattice spacing. To leading order in the lattice spacing, this effective Lagrangian has the standard continuum form. The above analysis requires that the mass term in this leading order effective Lagrangian must have the form:

$$Z_mm_f\overline{\psi }\psi +c\left\{\overline{\psi }\mathrm{\Omega }^{}P_R\psi +\overline{\psi }\mathrm{\Omega }P_L\psi \right\},$$
(40)
where $`c`$ is a constant with the dimensions of mass. Here the field $`\psi `$ represents a conventional continuum multiplet of quark fields and all quantities carry their physical dimensions. The first piece is the normal chiral symmetry breaking introduced by the input mass $`m_f`$. The second comes from mixing between the walls and is required by the extended symmetry of Eq. 37 and Eq. 39 to be linear in $`\mathrm{\Omega }`$. Thus, this induced mass term can occur to first order in the mixing between the walls, permitting $`ce^{\alpha L_s}/a`$. With the conventional choice of Shamir, $`\mathrm{\Omega }_{a,b}=\delta _{a,b}`$, the second term in Eq. 40 reduces to our usual residual mass term with $`am_{\mathrm{res}}10^3`$ .

In a similar fashion the $`𝒪(a)`$ effective Lagrangian will contain a clover term induced by mixing between the walls, again to first order in $`e^{\alpha L_S}`$, since it also has the permitted $`SU(N_f)_LSU(N_f)_R`$ chiral structure:

$$ac_1\left\{\overline{\psi }\sigma _{\mu \nu }F_{\mu \nu }\mathrm{\Omega }^{}P_R\psi +\overline{\psi }\sigma _{\mu \nu }F_{\mu \nu }\mathrm{\Omega }P_L\psi \right\}.$$
(41)
where $`c_1e^{\alpha L_s}`$ is $`𝒪`$($`am_{\mathrm{res}}`$). Thus, such a term is suppressed both by the lattice spacing and by the smallness of $`m_{\mathrm{res}}`$.

If we extend these considerations to $`𝒪(a^2)`$ terms in the effective Lagrangian, we can conclude that a four-fermi operator of the form

$$c_2a^2(\overline{\psi }\psi )(\overline{\psi }\psi ),$$
(42)
where $`c_2`$ is a constant, cannot occur to order $`e^{\alpha L_s}`$. Since this operator will become a chiral singlet only when contracted with two powers of the matrix $`\mathrm{\Omega }`$ or one power of $`\mathrm{\Omega }`$ and one powers of the mass matrix $`M`$ of Eq. 18, the coefficient of such an operator must contain a double suppression $`c_2e^{2\alpha L_s}`$ or a further factor of $`m_f`$.

## V Simulation Details

In the following discussions much use will be made of the momentum space quark propagator in Landau gauge. The first step in calculating this quantity is to fix the gauge. We implement Landau gauge fixing by iteratively sweeping over all lattice sites, maximising the functional

$$\underset{x,\mu }{}\left[\mathrm{Tr}\left(U_\mu (x)+U_\mu ^{}(x)\right)\right].$$
(43)
At each lattice site we determine a gauge transformation matrix, $`g(x)`$, which increases the value of Eq. 43. The maximum is achieved when

$$\underset{x}{}\mathrm{Tr}\left[B(x)B^{}(x)\right]$$
(44)
is zero, where

$$B(x)=A(x)A^{}(x)\frac{1}{3}\mathrm{Tr}(AA^{})$$
(45)
$$A(x)=\underset{\mu }{}\left[\stackrel{~}{U}_\mu ^{}(x)+\stackrel{~}{U}_\mu (x\mu )\right]$$
(46)
and

$$\stackrel{~}{U}_\mu (x)=g(x)U_\mu (x)g^{}(x+\mu ).$$
(47)
In practice we stop when the quantity in Eq. 44 is smaller than $`10^8`$.

On this gauge-fixed configuration the quark propagator, $`S(x,0)`$, from one source, denoted as $`0`$, to all possible sinks is then calculated. A discrete Fourier transform is then performed over the sink positions giving,

$$S(p,0)=\underset{x}{}\mathrm{exp}\left(ip^{\mathrm{latt}}x\right)S(x,0),$$
(48)
with

$$p_\mu ^{\mathrm{latt}}=\frac{2\pi }{L_\mu }n_\mu ,$$
(49)
where $`\mu `$ is one of $`x`$,$`y`$,$`z`$ or $`t`$ and $`n_\mu `$ may in principle lie in the range $`0L_\mu 1`$. In practice, however, only a subset of this range is used.

Unless otherwise stated all the data that will be presented is from calculations on a $`16^3\times 32\times 16`$ lattice (where the last number refers to the extent of the lattice in the fifth dimension). The simulation was performed at $`\beta =6.0`$ with $`2000`$ heatbath sweeps between every configuration and with $`2000`$ thermalisation sweeps performed at the outset. In total $`142`$ configurations were generated. For this lattice size the momentum range was restricted to those momenta for which $`n_\mu =0,1,2`$ for $`\mu =x,y,z`$ and $`n_t=0,1,2,3`$. Quark propagators for $`5`$ bare masses, $`m_f=0.01`$, 0.02, 0.03, 0.04 and 0.05 were calculated all using $`M_5=1.8`$.

The results will often be quoted against the square of the absolute momentum, where this refers to the Euclidean inner product of the momentum defined in Eq. 49. To be more specific

$$\left(ap\right)^2=\underset{\mu }{}p_\mu ^{\mathrm{latt}}p_\mu ^{\mathrm{latt}},$$
(50)
where $`p`$ is dimensionful.

## VI The Renormalised Propagator

Before we treat the more complicated situation of the fermion bilinears it is necessary to first consider the renormalisation of the quark propagator. Neglecting, for the moment, potential contributions from lattice artifacts the renormalised quark field may be defined as

$$q_{\mathrm{ren}}(x)=Z_q^{\frac{1}{2}}q_0(x).$$
(51)
If we similarly introduce a renormalised mass, defined by

$$m_{\mathrm{ren}}=Z_mm_0,$$
(52)
where $`m_0`$ represents a generic multiplicatively renormalised bare mass, then the renormalised propagator may be written

$$S_{\mathrm{ren}}(p,m_{\mathrm{ren}})=Z_qS_0(p,m_0)|_{m_0=m_{\mathrm{ren}}/Z_m}.$$
(53)
Both $`Z_q`$ and $`Z_m`$ are fixed in the RI scheme by requiring that the renormalised propagator obey the Euclidean space relations

$$\underset{m_{\mathrm{ren}}0}{lim}\frac{i}{12}\mathrm{Tr}\left(\frac{S_{\mathrm{ren}}^1}{\begin{array}{c}/\hfill \\ p\hfill \end{array}}(p)\right)_{p^2=\mu ^2}=1$$
(54)
$$\underset{m_{\mathrm{ren}}0}{lim}\frac{1}{12m_{\mathrm{ren}}}\mathrm{Tr}\left(S_{\mathrm{ren}}^1(p)\right)_{p^2=\mu ^2}=1.$$
(55)
While Eq. 54 and Eq. 55 seem to give a simple and appealing way to calculate $`Z_q`$ and $`Z_m`$ by directly applying them to the lattice propagators, the effect of both lattice artifacts and spontaneous chiral symmetry breaking must be considered.

Lattice actions with explicit chiral symmetry breaking require an additive renormalisation of the input mass, which may be taken into account for domain wall fermions by making the replacement

$$m_0m_f+m_{\mathrm{res}},$$
(56)
in the equations above. However, the effects of lattice artifacts on the correct definition of the renormalised and improved quark field are more complicated. They have been studied in Ref.  for Wilson fermions, where it is noted that there are three terms that may mix with the definition at $`𝒪`$(a) in the lattice spacing, giving rise to an expression for the improved and renormalised quark field of

$$q_{\mathrm{ren}}=Z_q^{\frac{1}{2}}\left(1+b_qma\right)\left\{1+ac_{}^{}{}_{q}{}^{}\left(\begin{array}{c}/\hfill \\ D\hfill \end{array}+m_{\mathrm{ren}}\right)+ac_{\mathrm{NGI}}\begin{array}{c}/\hfill \\ \hfill \end{array}\right\}q_0,$$
(57)
where   $`/`$  $``$  may appear because the gauge is fixed. If such extra terms appear then conditions must be found that allow them to be subtracted from the bare quark field before Eq. 54 and Eq. 55 may be applied. In the context of simulations using $`𝒪`$(a) improved Wilson action at $`\beta =6.0`$ these terms have been found to give significant contributions to the form of the propagator . In particular $`c_q^{}`$ was found to be large. However, they all break chiral symmetry and so, following the arguments of Section IV, should be suppressed by a factor of $`O(am_{\mathrm{res}})`$ for simulations using domain wall fermions. As such, studying the form of the propagator provides an excellent test of the chiral properties of domain wall fermions.

The effects of spontaneous chiral symmetry breaking on the form of the propagator are well known . The most noticeable effect is that the trace of the inverse propagator picks up an extra contribution, which at lowest order in a power expansion in $`1/p^2`$ may be described as

$$\frac{1}{12}\mathrm{Tr}\left(S_{\mathrm{ren}}^1(p)\right)=m_{\mathrm{ren}}+C_1\frac{\overline{q}q}{p^2}+\mathrm{}$$
(58)
where, at first order in perturbation theory, $`C_1=4\pi \alpha _s/3`$. Putting Eq. 58 and Eq. 57 together, the predicted form for the trace of the lattice quark propagator is

$`{\displaystyle \frac{1}{12}}\mathrm{Tr}\left(S_{latt}^1(ap)\right)`$ $`=`$ $`\mathrm{}+{\displaystyle \frac{a^3\overline{q}q}{(ap)^2}}C_1Z_q`$ (61)
$`+Z_mZ_q\left\{am_f+am_{\mathrm{res}}\right\}`$
$`+\mathrm{\hspace{0.17em}2}\left(c_{\mathrm{NGI}}Z_qc_{}^{}{}_{q}{}^{}\right)(ap)^2+\mathrm{},`$
where terms of $`𝒪`$($`mc_{\mathrm{NGI}}`$) have been neglected.

In Figure 1 we plot the left-hand side of Eq. 61 versus $`(ap)^2`$ for a variety of values for $`m_f`$. As can be seen, for domain wall fermions this quantity approaches a constant value for moderately large values of $`(ap)^2`$. Also, it is encouraging that while at low momenta the effects of spontaneous chiral symmetry are visible, there is no evidence to suggest appreciable effects from explicit chiral symmetry breaking. In particular, there is no evidence of a large additive mass renormalisation. This is visible in Figure 2, which shows the result of a linear extrapolation of the data to the point $`m_f=0`$. Figure 3 shows the slope of this extrapolation, which from Eq. 61 is $`Z_mZ_q`$ at large $`(ap)^2`$.

### A Extracting $`Z_q`$ from the Propagator

The extraction of $`Z_q`$ from the propagator via Eq. 54 is numerically challenging due to the need for a discrete derivative to be calculated. A much simpler method , is to calculate

$$Z_{}^{}{}_{q}{}^{}=i\frac{1}{12}\underset{\mu }{}\frac{ap_\mu }{(ap)^2}\mathrm{Tr}\left(\gamma _\mu S_{\mathrm{latt}}^1(ap)\right).$$
(62)
This quantity may then be related to $`Z_q`$ by a perturbative matching calculation performed in the continuum .

On the other hand, the use of $`Z_{}^{}{}_{q}{}^{}`$ to determine $`Z_q`$ introduces significant $`O(a^2)`$ errors through the choice of how the discretised momenta are defined. If we were to replace in Eq. 62 the definition of the lattice momentum in Eq. 49 by, for example,

$$\overline{p}_\mu =\frac{1}{a}\mathrm{sin}p_\mu a,$$
(63)
then the resulting $`Z_{}^{}{}_{q}{}^{}`$ would differ from that given in Eq. 62 by $`O(a^2)`$ because the trace includes an explicit factor of $`p_\mu \gamma _\mu `$. One can estimate the size of this error by using various definitions for the lattice momentum in the analysis. As will be shown later in Sec. XI, this uncertainty is roughly 10–20%. In Sections VIII and IX we introduce two methods for computing $`Z_q`$ which avoid this uncertainty.

## VII Flavour Non-Singlet Fermion Bilinears

### A Introduction

In the following the renormalisation of the flavour non-singlet fermion bilinear operators will be considered. To simplify notation explicit quark flavours ($`u`$ and $`d`$) will be used in the following equations. The most general fermion bilinear may be written as $`\overline{u}\mathrm{\Gamma }_id`$ with

$$\mathrm{\Gamma }_i\{1,\gamma _\mu ,\gamma _5,\gamma _\mu \gamma _5,\sigma _{\mu \nu }\},$$
(64)
where $`i`$ represents whatever indices the gamma matrices have. The renormalised operator is defined as

$$\left[\overline{u}\mathrm{\Gamma }d\right]_{\mathrm{ren}}=Z_\mathrm{\Gamma }\left[\overline{u}\mathrm{\Gamma }d\right]_0.$$
(65)
The factor $`Z_\mathrm{\Gamma }`$ is fixed in the RI scheme by defining the unrenormalised, amputated vertex function

$$\mathrm{\Pi }_{\mathrm{\Gamma },0}(p,q)=\frac{1}{V}d^4zd^4x_1d^4x_2e^{ip.x_1+iq.x_2}\left[\overline{u}\mathrm{\Gamma }d\right]_0(z)u_0(x_1)\overline{d_0}(x_2)_{\mathrm{AMP}},$$
(66)
the corresponding renormalised, amputated vertex function

$$\mathrm{\Pi }_{\mathrm{\Gamma },\mathrm{ren}}(p,q)=\frac{Z_\mathrm{\Gamma }}{Z_q}\mathrm{\Pi }_{\mathrm{\Gamma },0}(p,q)$$
(67)
and requiring that

$$\mathrm{\Lambda }_{\mathrm{\Gamma }_i,\mathrm{ren}}(p,p)_{p^2=\mu ^2}=\frac{1}{\mathrm{Tr}\left(_i\mathrm{\Gamma }_i\mathrm{\Gamma }_i\right)}\mathrm{Tr}\left(\underset{i}{}\mathrm{\Gamma }_i\mathrm{\Pi }_{\mathrm{\Gamma }_i,\mathrm{ren}}(p,p)\right)_{p^2=\mu ^2}=1.$$
(68)
Note, the corresponding, unrenormalised vertex amplitude is defined by

$$\mathrm{\Lambda }_{\mathrm{\Gamma },0}=\frac{1}{\mathrm{Tr}\left(_i\mathrm{\Gamma }_i\mathrm{\Gamma }_i\right)}\mathrm{Tr}\left(\underset{i}{}\mathrm{\Gamma }_i\mathrm{\Pi }_{\mathrm{\Gamma }_i,0}(p,p)\right),$$
(69)
so that

$$\mathrm{\Lambda }_{\mathrm{\Gamma }_i,\mathrm{ren}}(p,p)=\frac{Z_\mathrm{\Gamma }}{Z_q}\mathrm{\Lambda }_{\mathrm{\Gamma }_i,0}(p,p).$$
(70)
While Eq. 70 completely defines a procedure for calculating a renormalisation factor for the bilinear operator of interest, practically we need to be able to use a renormalisation condition that allows us to match to perturbative calculations. In general the value of $`\mathrm{\Lambda }_{\mathrm{\Gamma }_i,0}(p,p)`$ has contributions from intrinsically non-perturbative effects (such as those due to propagating pions) that perturbative calculations do not include. As we are interested in the value of the renormalisation factors in the perturbative regime we either apply the renormalisation condition at a high enough momenta such that the non-perturbative effects are suppressed, or remove such effects from the data and in the following that is what we will do. We will reserve the notation $`Z_i`$, $`i\{q,S,P,T,A,V\}`$ for the renormalisation factors in the perturbative regime.

### B $`Z_A`$ and $`Z_V`$

A (partially) conserved current that is normalised in a fashion which is consistent with the usual Ward-Takahashi identities will undergo no renormalisation and the corresponding $`Z_\mathrm{\Gamma }`$ will be unity. In particular, for domain wall fermions and the RI renormalisation scheme specified by Eqs. 54 and 68 and imposed at high-momentum, $`\mu \mathrm{\Lambda }_{QCD}`$, we expect

$$Z_𝒜=Z_𝒱=1.$$
(71)
However, on the lattice the (partially) conserved currents are not local and it is frequently more convenient to work with their local counterparts. Provided that these are related by a chiral transformation one still has

$$Z_A=Z_V.$$
(72)
This does not, however, mean that $`\mathrm{\Lambda }_A`$ must equal $`\mathrm{\Lambda }_V`$, and several mechanisms exist for splitting them away from each other at low energies.

Even if there are no significant effects from explicit chiral symmetry breaking, the effects of spontaneous chiral symmetry breaking must be taken into account. Consideration of the operator product expansion, to lowest order in powers of $`1/p^2`$, shows that $`\mathrm{\Lambda }_A`$ and $`\mathrm{\Lambda }_V`$ may get contributions from terms of the form

$$\frac{m_{\mathrm{ren}}^2}{p^2}$$
(73)
and

$$\frac{m_{\mathrm{ren}}\overline{q}q}{p^4}.$$
(74)
Since such terms, by their very nature, stem from chiral symmetry breaking they are not constrained to enter $`\mathrm{\Lambda }_A`$ and $`\mathrm{\Lambda }_V`$ with the same weight. At large momenta these terms are suppressed and do not effect the extraction of $`Z_A`$ and $`Z_V`$.

If the action being used explicitly breaks chiral symmetry $`Z_A`$ and $`Z_V`$ need not be equal, but their ratio will still be scale independent. This means that while $`\mathrm{\Lambda }_A`$ and $`\mathrm{\Lambda }_V`$ need not approach one another at high momenta, their ratio should become scale independent for large enough momenta.

While we want to work in the chiral limit for the extraction of the renormalisation factors it is also worthwhile to consider what the mass dependence of $`\mathrm{\Lambda }_A`$ and $`\mathrm{\Lambda }_V`$ should be (especially as we wish to extract the chiral limit from data measured in finite mass simulations). Requiring that the “generalised” symmetry introduced in Eq. 21 is satisfied, constrains the mass dependence to either be of the same form as Eq. 73 (a single power of mass multiplying something that breaks chiral symmetry - and therefore by the argument of the previous paragraph damped with momentum) or proportional to a second or higher power of mass. In the latter case any effect on our data should be negligible.

The above considerations suggest looking at the quantity

$$\mathrm{\Lambda }_A\mathrm{\Lambda }_V.$$
(75)
This is shown in Figure 4 and, as with the case for the quark propagator, while the effects of non-perturbative breaking terms are visible at low momenta, they are damped at higher momenta. There is no significant signal for effects from explicit chiral symmetry breaking, since $`\mathrm{\Lambda }_A\mathrm{\Lambda }_V`$ is tending to zero. At momenta of interest, there also seems to be no significant splitting due to non-perturbative effects with the difference between $`\mathrm{\Lambda }_A`$ and $`\mathrm{\Lambda }_V`$ being less than $`1\%`$ in the chiral limit at $`(ap)^2=0.8`$ and smaller for momenta above this. This being so it is sensible to use the quantity $`\frac{1}{2}\left(\mathrm{\Lambda }_A+\mathrm{\Lambda }_V\right)`$ for the extraction of both $`Z_A/Z_q`$ and $`Z_V/Z_q`$ to increase the statistical accuracy. This is shown in Figure 5.

### C  $`Z_S`$ and $`Z_P`$

For a theory with chiral symmetry the RI scheme preserves the well known $`\overline{MS}`$ relations

$`Z_S`$ $`=`$ $`Z_P`$ (76)
$`Z_m`$ $`=`$ $`{\displaystyle \frac{1}{Z_S}}.`$ (77)
If the potential for explicit chiral symmetry breaking is taken into account, these equalities cease to be valid, but the quantities $`Z_S/Z_P`$, $`Z_SZ_m`$ and $`Z_mZ_P`$ are expected to be scale independent. However, lattice studies using both the Wilson and staggered actions have shown that the ratio $`\mathrm{\Lambda }_P/\mathrm{\Lambda }_S`$, which in perturbation theory would be equal to $`Z_S/Z_P`$ and therefore might be expected to be momentum independent up to small corrections, is strongly momentum and mass dependent, with the bulk of this dependence arising from $`\mathrm{\Lambda }_P`$. It is instructive to consider the source of this discrepancy . We start from the continuum axial Ward-Takahashi identity which is derived by taking the axial variation of the quark propagator. This reads

$$(m_u+m_d)𝑑zu(x_1)\left[\overline{u}\gamma _5d\right](z)\overline{d}(x_2)=\gamma _5d(x_1)\overline{d}(x_2)+u(x_1)\overline{u}(x_2)\gamma _5.$$
(78)
Moving to momentum space gives

$$(m_u+m_d)u\left[\overline{u}\gamma _5d\right]\overline{d}(p,p)=\gamma _5d\overline{d}(p)+u\overline{u}(p)\gamma _5,$$
(79)
which in the $`m_um_dm`$ limit gives

$$2mu\left[\overline{u}\gamma _5d\right]\overline{d}(p,p)=\{\gamma _5,S(p)\}.$$
(80)
Amputating and tracing with $`\gamma _5`$ yields

$$m\mathrm{\Lambda }_P(p,p)=\frac{1}{12}\mathrm{Tr}\left(S^1(p)\right).$$
(81)
Neglecting all lattice artifact terms except the additive mass renormalisation (which is justified by the discussions in Section VI), this leads to an approximate expression for $`\mathrm{\Lambda }_{P,\mathrm{latt}}`$, including the first order contribution of spontaneous chiral symmetry breaking:

$$\left(m_f+m_{\mathrm{res}}\right)\mathrm{\Lambda }_{P,\mathrm{latt}}(ap,ap)=\mathrm{}+\frac{a^2\overline{q}q}{(ap)^2}C_1Z_q+Z_mZ_q\left(m_f+m_{res}\right),$$
(82)
neglecting terms of $`𝒪`$($`(ap)^2`$). While in the absence of the condensate term this equation reduces to $`Z_P=1/Z_m`$, the condensate term, which is clearly visible in Figure 2, gives rise to a pole term of the form

$$\frac{\overline{q}q}{\left(m_f+m_{res}\right)(ap)^2}C_1Z_q$$
(83)
in $`\mathrm{\Lambda }_{P,\mathrm{latt}}`$. Figure 6 shows the data for $`\mathrm{\Lambda }_{P,\mathrm{latt}}`$ with this effect clearly visible, with the rise for small $`(ap)^2`$ becoming much more pronounced as the mass decreases.

A similar argument may be put forward for $`\mathrm{\Lambda }_S`$. In this case it is the Ward-Takahashi identity arising from a vector rotation of the fields.

$$(m_um_d)𝑑zu(x_1)\left[\overline{u}d\right](z)\overline{d}(x_2)=d(x_1)\overline{d}(x_2)u(x_1)\overline{u}(x_2)$$
(84)
Moving to momentum space, this gives

$$(m_um_d)u\left[\overline{u}d\right]\overline{d}(p,p)=d\overline{d}(p)u\overline{u}(p),$$
(85)
which in the $`m_um_d`$ limit tends to

$$u\left[\overline{u}d\right]\overline{d}(p,p)=\frac{}{m}S(p).$$
(86)
Finally, amputating and tracing gives

$$\mathrm{\Lambda }_S=\frac{1}{12}\frac{\mathrm{Tr}\left[S(p)^1\right]}{m}.$$
(87)
Note that this relation should be exact for domain wall fermions (with $`m=m_f`$) for any value of $`L_s`$. Using Eq 61, and noting that both the residual mass and the renormalisation factors should be independent of $`m_f`$, gives the approximate expression

$$\mathrm{\Lambda }_{S,\mathrm{latt}}=Z_mZ_q+\frac{C_1Z_q}{(ap)^2}\frac{a^3\overline{q}q}{(am_f)}.$$
(88)
If the mass dependence of $`\overline{q}q`$, for small masses, is proportional to only positive powers of the mass, then the second term in Eq. 88 is almost certainly unimportant as it is suppressed in exactly the region of parameter space in which we are working: large momenta and small masses. (The effect might be larger than naively expected, as $`\overline{q}q`$ is quadratically divergent in the lattice spacing.) It is necessary, however, to consider the effects of fermionic zero-modes on $`\overline{q}q`$. Assuming a theory with chiral symmetry, the spectral decomposition of the quark propagator leads to an expected form for $`\overline{q}q`$, on a single configuration C, of

$`\overline{q}q_C`$ $`=`$ $`Tr\left[\left(\begin{array}{c}/\hfill \\ D\hfill \end{array}+m\right)^1\right]`$ (97)
$`=`$ $`{\displaystyle \frac{n_0}{mV}}+{\displaystyle \frac{1}{V}}{\displaystyle \underset{\lambda _n>0}{}}{\displaystyle \frac{1}{m+i\lambda _n}};n_00.`$ (98)
where $`n_0`$ is the number of fermionic zero-modes, $`V`$ is the four dimensional space-time volume and the $`\lambda _n`$ are such that $`\begin{array}{c}/\hfill \\ D\hfill \end{array}\psi _n=i\lambda _n\psi _n`$. The number of such zero modes should grow more slowly than the volume, and so the first term in Eq. 98 will vanish in the infinite volume limit. However, for the lattice parameters used for this simulation the effects of zero-modes have been found to be noticeable in both $`\overline{q}q`$ and hadronic spectrum calculations  and so must be considered for the present case. Comparing Eq. 88 and Eq. 98 shows that, as $`m0`$ for fixed momentum, $`\mathrm{\Lambda }_S`$ gets a large contribution from zero-modes of the form

$$\frac{1}{m^2}n_0\frac{C_1Z_q}{p^2},$$
(99)
that must be subtracted before $`Z_S/Z_q`$ may be calculated from Eq. 68. Figure 7 shows our data for $`\mathrm{\Lambda }_{S,\mathrm{latt}}`$. While the effect of the condensate term is smaller than that for $`Z_P`$, it is noticeable for the lighter masses.

### D Fitting the Pole Terms

Considering Eq. 99 and moving to lattice notation, the method for extracting $`Z_S`$ from $`\mathrm{\Lambda }_S`$ becomes clear. Working at a fixed momentum, $`\mathrm{\Lambda }_S`$ may be fitted to the form

$$\mathrm{\Lambda }_{S,\mathrm{latt}}=\frac{c_{1,S}}{\left(am_f+𝒪(am_{\mathrm{res}})\right)^2}+c_{2,S}+c_{3,S}(am_f)^2,$$
(100)
with $`Z_q/Z_S`$ being given by $`c_{2,S}`$. While one might naively expect the denominator in the above equation to be $`am_f+am_{\mathrm{res}}`$, as shown in Ref.  the residual chiral symmetry breaking effects that appear in the pole term in $`\overline{q}q`$ are not parameterised precisely by $`am_{\mathrm{res}}`$ since the singular behavior of the pole enhances what are expected to order $`a^2`$ variations in the quantity $`am_{\mathrm{res}}`$. One only knows that the residual chiral symmetry breaking effects are of $`𝒪`$($`am_{\mathrm{res}}`$). However, for the pole subtractions in this paper we have used precisely $`am_f+am_{\mathrm{res}}`$. This is justified since the statistical errors on our data are such that the fits are insensitive to the exact value of $`m_{\mathrm{res}}`$.

The situation for the $`Z_P`$ extraction is slightly more complicated. Examining Eq. 82 and Eq. 98, shows that $`\mathrm{\Lambda }_P`$ should have a double pole of the form

$$\mathrm{\Lambda }_{P,\mathrm{latt}}=\frac{c_{1,P}}{(am_f+am_{\mathrm{res}})^2}+\frac{c_{2,P}}{(am_f+am_{\mathrm{res}})}+c_{3,P}+c_{4,P}(am_f)^2,$$
(101)
with $`c_{3,P}`$ being equal to $`Z_q/Z_P`$ and the quadratic mass pole due to zero-mode effects in $`\overline{q}q`$. For simplicity, we have again used $`am_f+am_{\mathrm{res}}`$ in the first term of the right-hand side in Eq. 101. For practical purposes, however, the need to fit to the quadratic term may be avoided by working with $`m_f0.02`$. Good evidence that the above fitting forms are correct is shown in Figure 8. This shows the average, over all the momenta in the range $`0.5<\left(ap\right)^2<2.0`$, of the $`\chi ^2`$ per degree of freedom for a correlated fit to the above forms for $`m_f`$ dependence with the power of the pole treated as a free parameter. One sees that a single pole is favored for the $`Z_P`$ case while a double is identified in the fit for $`Z_S`$. Further evidence is provided by considering the resulting values for $`Z_S/Z_q`$ and $`Z_P/Z_q`$. Figure 9 shows a comparison between the extracted values of these two quantities. As chiral symmetry would predict for $`Z_S/Z_q`$ and $`Z_P/Z_q`$, the two quantities coincide at large momenta. This provides an excellent test of both the extraction method and the chiral properties of domain wall fermions. This can be further seen by comparing $`Z_P`$ and $`Z_m`$ (as calculated from the trace of the inverse propagator), which is shown in Figure 10. This product is clearly very close to unity.

### E  $`Z_T`$

The tensor density renormalisation, while sometimes neglected in bilinear renormalisation coefficient calculations, is a quantity of use to current lattice simulations . An extraction of its value will be postponed to Section X, but for completeness a plot of $`\mathrm{\Lambda }_T`$ is shown in Figure 11.

## VIII Extracting $`Z_q`$ from the Exact Ward Identities.

The vector Ward-Takahashi identity of Eq. 24 is exact at finite lattice spacing. As such, the renormalisation coefficient for the conserved vector current defined in Eq. 27, is equal to unity. Additionally, the considerations of Section IV show that through first order in the residual chiral symmetry breaking the extra $`J_{5q}^a`$ in Eq. 31 can be completely absorbed into the additive renormalisation of the mass so that the axial Ward-Takahashi identity, for low-energy physics, takes on the normal, continuum form. Therefore, the renormalisation factor for the axial current defined by Eq. 35 should also be unity to a good approximation.

The above facts can be used to compute the quark field renormalisation $`Z_q`$ from Eq. 68, as applied the conserved vector and axial currents. For the case of the conserved vector current Eq. 68 reads:

$$Z_q^1Z_𝒱\frac{1}{48}\mathrm{Tr}(\gamma ^\mu \mathrm{\Lambda }_𝒱^\mu (p))=1$$
(102)
which therefore implies

$$Z_q=\frac{1}{48}\mathrm{Tr}(\gamma ^\mu \mathrm{\Lambda }_𝒱^\mu (p)).$$
(103)
A similar equation also holds for the conserved axial current.

Because the formulae for the conserved currents contain fermion fields at separate points on the lattice as well as summation over all $`s`$, the calculation of the matrix elements for these operators can be very expensive. We used a random source estimator to compute the part of the sum between $`s=1`$ to $`s=L_s2`$. Also, instead of calculating all four components of $`\mathrm{\Lambda }^\mu (p)`$ for a given momentum $`p`$, we calculate $`\mathrm{\Lambda }^0(p)`$ for momenta related to $`p`$ by interchange of its $`0`$th component with each of the other three. We then use the equality

$`\mathrm{Tr}\left(\gamma ^\mu \mathrm{\Lambda }_𝒱^\mu (p^0,p^1,p^2,p^3)\right)`$ $`=`$ $`\mathrm{Tr}\left(\gamma ^0\mathrm{\Lambda }_𝒱^0(p^0,p^1,p^2,p^3)\right)`$ (107)
$`+\mathrm{Tr}\left(\gamma ^0\mathrm{\Lambda }_𝒱^0(p^1,p^0,p^2,p^3)\right)`$
$`+\mathrm{Tr}\left(\gamma ^0\mathrm{\Lambda }_𝒱^0(p^2,p^1,p^0,p^3)\right)`$
$`+\mathrm{Tr}\left(\gamma ^0\mathrm{\Lambda }_𝒱^0(p^3,p^1,p^2,p^0)\right),`$
to obtain the needed result. As the time to Fourier transform a matrix element is negligible in comparison with the calculation of the matrix element itself, this allows us to obtain the result with only a quarter of running time of direct calculation of $`\mathrm{\Lambda }_𝒱^\mu `$.

Because the volume $`16^3\times 32`$ used in the simulations is not symmetric, strictly speaking the above equation is not exact, as the third component of momentum is not related by symmetry to the first and second ones. However, this difference is suppressed by two powers of the lattice spacing and, in practice, the results obtained for the last three terms in Eq. 107 all agree within statistical errors for the momenta used. We used momenta with the first three integer components no larger than $`2`$ and the last component equal to $`0`$ or $`2`$. The two exceptions were the momenta with integer components $`(2,2,2,0)`$ and $`(2,2,2,2)`$ that were excluded since they would require usage of momenta $`(0,2,2,4)`$ and $`(1,2,2,4)`$.

Figure 12 shows the difference between $`Z_q`$ calculated using the axial and vector Ward-Takahashi identities, while Figure 13 shows the average. As can be seen from Figure 12, while for low momenta the two methods give different results, this difference is damped at large momenta as would be expected if this effect stems from spontaneous (rather than explicit) chiral symmetry breaking. Again, this provides strong evidence that the effects of explicit chiral symmetry breaking are negligible in these calculations.

## IX $`Z_A`$ from Hadronic Matrix Elements

As mentioned in Section VIII, to a good approximation the axial current defined by Eq. 35 is conserved, and therefore has a renormalisation coefficient equal to the identity. This provides a simple way to calculate the renormalisation coefficient of the local axial current operator, $`Z_A`$, directly from hadronic matrix element calculations. One method, that has been used in Ref. , is to note that the matrix element of any operator with the renormalised axial current is a well defined quantity, which will be independent of the interpolating operator for the axial current at distances above the scale of the lattice spacing. Therefore,

$$𝒜_0(t_1)\overline{q}\gamma _5q(t_2)=Z_AA_0(t_1)\overline{q}\gamma _5q(t_2).$$
(108)
for $`|t_1t_2|a1`$. A full discussion of this method and the results are given in Ref. , but it is useful to summarise the results here. Table I collects together values for $`Z_A`$ for several $`L_S`$ values on a $`\beta =6.0`$, $`16^3\times 32`$ lattice with $`M_5=1.8`$. The quoted errors are statistical only.

## X Renormalisation Group Behaviour

The previous sections have provided an extraction of the renormalisation coefficients of interest taking into account the possible effects stemming from chiral symmetry breaking (either explicit or spontaneous). In general, perturbation theory predicts that these coefficients may be logarithmically dependent on the momentum scale. Lattice artifacts may also cause the result to depend on the definition of momentum on the lattice. When small these lattice artifacts will be manifest in the data as added terms proportional to $`\left(ap\right)^2`$.

The simplest approach to using the renormalisation coefficients calculated in this paper is to take the (mass-pole subtracted) value for $`\mathrm{\Lambda }_i\left(ap\right)`$ as $`\left(Z_q/Z_i\right)\left(ap\right)`$. However, we need to address two significant uncertainties:

First, this choice assumes, without verification, that the $`𝒪`$($`(ap)^2`$) lattice artifacts are small. We should attempt to understand the momentum dependence of the amplitude $`\mathrm{\Lambda }_i\left(ap\right)`$ to determine the size of the $`𝒪`$($`(ap)^2`$) contamination and, if possible, remove it.

Second, most operators of interest in lattice calculations are ultimately defined using continuum perturbation theory. Making this connection requires the use of perturbation theory to connect the RI scheme and momentum scale used here with the renormalisation procedure and momentum scale used in the original continuum definition, typically the $`\overline{MS}`$ scheme. This final matching step between the normalised lattice and continuum operators is done at a specific momentum scale for the renormalised lattice operator. In general, both the normalisation of lattice operators and the matching coefficients will depend on this momentum scale. It is not known a priori how many loops in perturbation theory must be calculated to correctly describe the momentum range probed in current lattice calculations, or even if perturbation theory can describe the region we are studying. Comparing the momentum behaviour predicted from perturbation theory to that of the data therefore provides an important consistency check for the general framework of the method.

Our approach to comparing the known perturbative running of the quantities of interest to our numerical data will be to divide the data by the predicted renormalisation group running, with the overall normalisation set by requiring that at the point $`(ap)^2=1`$ this divisor is one. If the perturbative result correctly describes the data, and the effect of $`(ap)^2`$ terms may be neglected, the result will be completely scale independent.

There are three components that are needed to calculate these quantities:

1. The anomalous dimensions for the operators from perturbation theory.
2. The ratio of $`Z^{\mathrm{RI}}/Z^{\overline{\mathrm{MS}}}`$ must be known to perform the matching. Since the renormalisation condition that determines $`Z^{\mathrm{RI}}`$ is well defined both on the lattice and when in continuum dimensional regularisation this ratio may be calculated perturbatively using the latter regularisation. In general it can be expanded as

$$\frac{Z^{\mathrm{RI}}}{Z^{\overline{\mathrm{MS}}}}=1+\frac{\alpha _s}{4\pi }Z_0^{(1)\mathrm{RI}}+\frac{\alpha _s^2}{\left(4\pi \right)^2}Z_0^{(2)(\mathrm{RI})}+\mathrm{}.$$
(109)
For a consistent treatment this ratio need only be known to one less power of $`\alpha _s`$ than the running is known.
3. A lattice value for $`\alpha _s`$. The value of $`\alpha _s`$ affects the scale dependence of both the matching and running for this calculation. For this work the value of $`\alpha _s`$ was calculated at three loops using a lattice value of $`\mathrm{\Lambda }_{QCD}`$ taken from Ref.  as

$$\mathrm{\Lambda }_{QCD}=238\pm 19\mathrm{M}\mathrm{e}\mathrm{V}.$$
(110)
To do this consistently with the way the lattice treatment in Ref.  was performed, their value of $`r_0=0.5\mathrm{fm}`$ was taken and converted into a lattice spacing using the results of Ref. . For the dimensionful scales that we will quote, we set the physical scale through the rho mass computed with domain wall fermions , which for $`\beta =6.0`$ gives

$$a=0.520(11)\mathrm{GeV}^1.$$
(111)

Both $`Z_A`$ and $`Z_V`$ should be scale independent, but this is not the case for $`Z_q`$. Figure 14 shows both $`1/\mathrm{\Lambda }_A`$ and the scale invariant (SI) quantity calculated as described above:

$$\mathrm{\Lambda }_A^{SI}\left((ap)^2\right)=\mathrm{\Lambda }_A\left((ap)^2\right)/C_A\left((ap)^2\right).$$
(112)
The quantity $`C_A`$ is determined through three loops using the anomalous dimension coefficients calculated in Refs.  as described in Appendix B. It is normalised so that $`C_A(1)=1`$. As can be seen, in this case the renormalisation group running actually goes in the opposite direction from the data. The scale dependence of this data, either predicted or actual, is, however, very small and a plausible explanation for this is an $`(ap)^2`$ error. Indeed, when a linear fit of the SI data versus $`(ap)^2`$ is performed, for $`0.8<(ap)^2<2.0`$, the gradient is $`0.02`$. A more compelling test of the renormalisation group behaviour is provided by studying the data for $`Z_S/Z_q`$. In this case the predicted scaling behaviour over the range of momenta studied is much larger and, as Figure 15 shows, the agreement between the predicted behaviour and the data is impressive (with a gradient, in this case, of $`0.003`$). The values for $`Z_S/Z_q`$ versus momentum used here are taken after the mass-pole has been subtracted and, again, the three loop results for the running are taken from Refs. . Unfortunately, a matching calculation for $`Z_T`$ could not be found in the literature, so the data could only be compared to the one loop running (which is taken from Ref. ). The SI quantity so calculated is shown in Figure 16 and has a gradient of $`0.02`$.

Taking the interpretation that the remaining scale dependence is due to $`𝒪`$($`(ap)^2`$) effects, the correct way to extract the renormalisation coefficients is to first construct the SI quantity as described above, and then fit any remaining scale dependence  to the form

$$y=c_1+c_2(ap)^2,$$
(113)
for a range of momenta that is chosen to be “above” the region for which condensate effects are deemed to be important. Table II shows the fitted values for the RI and $`\overline{MS}`$ scheme renormalisation coefficients using a fitting rage of $`0.8<(ap)^2<2.0`$. Now that the renormalisation group running has been taken into account, it is possible to make a comparison of the various methods of calculating $`Z_q`$ and thus give final results for the renormalisation factors. Table II already gives $`Z_q`$ as calculated from the conserved currents (Figure 17 shows the momentum dependence of both the SI and bare form). Another simple way to derive this quantity is by taking $`Z_A/Z_q`$ from Table II and combining it with the value of $`Z_A`$ obtained from hadronic matrix elements. This gives $`Z_q=0.805(3)(15)`$. This, approximately $`5\%`$, difference may be taken as an indication of the size of the systematic errors.

## XI $`Z_q`$ from the Propagator – Results

Here we show results for the wave function renormalisation computed through Eq. 62 and demonstrate that this methods contains a comparatively large systematic uncertainty due to the ambiguity in defining discrete momentum. We use the perturbative matching between $`Z_{}^{}{}_{q}{}^{}`$ and $`Z_q`$, as given in Ref. .Their convention is that a given $`Z`$-factor in Ref.  is the reciprocal of ours. Then the SI $`Z_q`$ can be constructed as described above.

Figure 18 shows the SI $`Z_q`$ using $`Z_{}^{}{}_{q}{}^{}`$ defined in Eq. 62, and Figure 19 shows the SI $`Z_q`$ where the replacement

$$ap_\mu \mathrm{sin}(ap_\mu )a\overline{p}_\mu $$
(114)
is made in Eq. 62. Note that the former is plotted vs. $`(ap)^2`$ and the latter vs. $`(a\overline{p})^2`$. We use the data at $`m_f=0.02`$ since no mass dependence can be observed.

As in the previous section, we extrapolate to $`(ap)^2=0`$. We find for the data in Figure 18, $`Z_q=0.715\pm 0.007\pm 0.040`$, where the first error is statistical and the second comes from different choices for the range of momenta over which to fit. The data in Figure 19 give $`Z_q=0.733\pm 0.007\pm 0.050`$. We can further probe these discretisation uncertainties by extrapolating $`Z_q^1`$ to zero $`(ap)^2`$ or $`(a\overline{p})^2`$. This results in $`Z_q=0.732\pm 0.006\pm 0.020`$ and $`Z_q=0.772\pm 0.006\pm 0.020`$, respectively.

The spread in values of $`Z_q`$ obtained depending on momentum ranges and on the definition of discrete momenta mean that extracting a SI $`Z_q`$ in the same manner as the $`Z`$’s for bilinear operators is less precise. However, the results are in rough agreement with the more precise methods described above.

## XII Comparison with Perturbation Theory

All the renormalisation factors considered above have also been calculated in lattice perturbation theory, at the one loop level, for the domain wall fermion action in the $`L_s\mathrm{}`$ limit . As we see little evidence of explicit chiral symmetry breaking effects in our study, the fact that the perturbative calculations have been performed in the $`L_s\mathrm{}`$ limit will probably not affect this comparison. However, a more serious issue is which $`M_5`$ value to use in the perturbative formulae.

The reason for this is easy to understand. Away from the walls, the massless domain wall fermion Lagrangian, Eq. 6, may be viewed as a simple extension of the standard Wilson fermion action to five dimensions with a negative mass term, $`M_5`$. Like the four dimensional Wilson mass term, $`M_5`$ undergoes a strong additive renormalisation, the size of which perturbation theory is not good at predicting. While a more accurate prediction may be made using tadpole improved perturbation theory , a good deal of ambiguity remains in the perturbative prediction of any quantity that is strongly dependent on $`M_5`$.

Further consideration of the similarity of the domain wall and standard Wilson actions, leads to a non-perturbative estimate of the magnitude of this additive renormalisation. As argued in Ref. , in the $`L_s\mathrm{}`$ limit, the effect of this additive renormalisation may be taken into account by using

$$M_5^{\mathrm{pert}}=M_5^{\mathrm{sim}}\left(4\frac{1}{2\kappa _c}\right),$$
(115)
in the perturbative equations, where $`\kappa _c`$ is the four-dimensional critical Wilson hopping parameter and $`M_5^{\mathrm{sim}}`$ refers to the value used in the non-perturbative simulation (in our case $`M_5^{\mathrm{sim}}`$=1.8). For $`\beta =6.0`$, this ansatz leads to shift of magnitude $`0.8`$, which has been found to describe the dependence of the pion mass squared as a function of $`M_5`$ to a good degree of accuracy in a numerical simulation with $`L_s=14`$.

Eq. 4.10 in Ref.  gives the complete one loop bilinear renormalisation constants in the $`\overline{MS}`$ scheme:

$`Z_\mathrm{\Gamma }^{\mathrm{total}}`$ $`=`$ $`((1w_0^2)Z_w)^1Z_\mathrm{\Gamma }.`$ (116)
Here $`Z_w`$, $`Z_2`$, and $`Z_\mathrm{\Gamma }`$ are to be computed from Eqs. 3.30, 3.42, and 4.11 and Tables 2 and 3 in Ref. , while $`w_0=1M_5`$. In the mean-field improved case, the above relations hold with the replacements $`w_0w_0^{MF}=w_0+4(1u)`$, $`Z_wZ_w^{MF}`$, $`Z_2uZ_2^{MF}`$, and $`Z_\mathrm{\Gamma }uZ_\mathrm{\Gamma }^{MF}`$ , whose values can also be computed from Tables 2 and 3 in Ref. . The factor $`u`$ in these formulae is the mean link variable in Feynman gauge. As it is not possible to use the value of the mean link in Feyman gauge we have instead used the fourth root of the plaquette and the perturbative results of Ref.  to convert the results of Ref. .

In Figure 20 and Figure 21 we plot $`Z_A^{\mathrm{total}}`$ and $`Z_S^{\mathrm{total}}`$, respectively, as functions of the variable $`M_5`$ in naive perturbation theory, in naive perturbation theory with the variable $`M_5`$ shifted according to Eq. 115 and in the mean-field improved case. To compute $`\alpha _s`$, we used the same input values for $`\mathrm{\Lambda }_{QCD}`$ and $`a`$ as in the perturbative running calculations in Section X. We obtain $`\alpha _s\left((ap)^2=1\right)=0.20`$. These figures show appreciable $`M_5`$ dependence. Our non-perturbative result is shown as a point corresponding to the single value of $`M_5=1.8`$ that we have studied.

The naive perturbation theory curve has a significant dependence on the precise value of $`\alpha _s`$. In the mean-field improved case this problem is not as serious as the coefficient of $`\alpha _s`$ is a factor of 2-3 times smaller. Examining Figures 20 and 21, one recognises that naive perturbation theory does a poor job of determining $`Z_A`$ or $`Z_S`$ giving values nearly 2 times too small for $`M_5=1.8`$. Introducing the shift of Eq. 115 improves the situation noticeably giving values 15% too small and to within a few percent, although the perturbative result is rapidly varying with $`M_5`$ in this case. The mean-field results differ from the non-perturbative result by around 5% in both cases.

## XIII Conclusions

In this paper we have described a first study of non-perturbative renormalisation of the quark field and flavour non-singlet fermion bilinear operators in the context of domain wall fermions. We presented a theoretical argument constraining the form that explicit chiral symmetry breaking effects may take, and found that numerically these are insignificant, as might be expected from the measured size of the additive mass renormalisation, $`m_{\mathrm{res}}`$, . However, systematic effects due to spontaneous chiral symmetry breaking and zero-modes are significant, but accurately follow the expected form and can be effectively subtracted away.

Renormalisation group invariant quantities were obtained in Section X by dividing the regularisation independent scheme coefficients by the three loop renormalisation group running (where available). The residual scale dependence of these quantities is small and was treated as an $`𝒪`$($`a^2`$) error. Three different quantities were used to determine the quark renormalisation factor: the off-shell vertex functions of the conserved vector and axial currents; the trace of the product of $`p^\mu \gamma ^\mu `$ and the off-shell quark propagator; and the combination of $`Z_A`$ as determined from hadronic matrix elements with the value of $`Z_A/Z_q`$ obtained in this study from the off-shell, axial vector vertex function. The technique of obtaining this $`Z_q`$ directly from the propagator suffers from large discretisation errors, but is roughly consistent with the other two methods which gave results differing by $`5\%`$.

In the final section we compared our results against the predictions of both standard and mean-field improved one loop perturbation theory.

## Acknowledgments

The authors would like to acknowledge useful discussions with Sinya Aoki. We thank RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy for providing the facilities essential for the completion of this work.

This research was supported in part by the DOE under grant # DE-FG02-92ER40699 (Columbia), by the NSF under grant # NSF-PHY96-05199 (Vranas), by the DOE under grant # DE-AC02-98CH10886 (Dawson-Soni) and by the RIKEN-BNL Research Center (Blum-Wingate).

## A The Running of $`\alpha _s`$

In the following the definitions

$$\alpha =\frac{g^2}{4\pi },$$
(A1)
and

$$C_F=\frac{N_C^21}{2N_C},$$
(A2)
will be used. The renormalised coupling may be defined in terms of the bare coupling by,

$$\alpha _b=Z_g^2\alpha _s\mu ^{2ϵ}.$$
(A3)
As $`\alpha _b`$ is completely independent of $`\mu `$,

$$\mu ^2\frac{d\alpha _s}{d\mu ^2}=ϵ\alpha _s+\beta $$
(A4)
with,

$$\beta =\alpha _s\frac{2}{Z_g}\mu ^2\frac{dZ_g}{d\mu ^2}$$
(A5)
The results for the beta function are most easily given in terms of the $`\beta _i`$ variable:

$$\frac{\beta (\alpha _s)}{4\pi }=\beta _0\left[\frac{\alpha _s}{4\pi }\right]^2\beta _1\left[\frac{\alpha _s}{4\pi }\right]^3\mathrm{}.$$
(A6)
The values used in the current work are summarised in Table IV. They are taken from Ref.  with the number of flavours set to zero (as we are working in the quenched case) and the number of colours set to three.

Once the $`\beta _i`$’s are known the running equation,

$$\mu \frac{d}{d\mu }\alpha _s=2\beta _0\frac{\alpha _s^2}{4\pi }2\beta _1\frac{\alpha _s^3}{(4\pi )^2}\mathrm{},$$
(A7)
may be solved:

* One loop solution:

$$\frac{\alpha _s}{4\pi }=\frac{1}{\beta _0\mathrm{ln}\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}$$
(A8)
* Two loop solution :

$$\frac{\alpha _s}{4\pi }=\frac{1}{\beta _0\mathrm{ln}\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}\frac{\beta _1\mathrm{ln}\mathrm{ln}\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}{\beta _o^3\mathrm{ln}^2\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}$$
(A9)
* Three loop solution :

$`{\displaystyle \frac{\alpha _s}{4\pi }}`$ $`=`$ $`{\displaystyle \frac{1}{\beta _0\mathrm{ln}\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}}{\displaystyle \frac{\beta _1\mathrm{ln}\mathrm{ln}\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}{\beta _o^3\mathrm{ln}^2\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}}`$ (A12)
$`+{\displaystyle \frac{1}{\beta _0^5\mathrm{ln}^3\left(\mu ^2/\mathrm{\Lambda }_{QCD}^2\right)}}\{\beta _1^2\mathrm{ln}^2\mathrm{ln}(\mu ^2/\mathrm{\Lambda }_{QCD}^2)\beta _1^2\mathrm{ln}\mathrm{ln}(\mu ^2/\mathrm{\Lambda }_{QCD}^2)`$
$`+\beta _2\beta _0\beta _1^2\}`$

## B The Running of the $`Z`$-factors

As mention previously, the renormalised operators we are working with are defined as

$$Z_OO_{bare}=O_{ren}.$$
(B1)
Requiring that the bare operator is independent of the renormalisation scale gives the RG equation,

$`\mu {\displaystyle \frac{d}{d\mu }}O_{ren}`$ $`=`$ $`{\displaystyle \frac{1}{Z_O}}\mu {\displaystyle \frac{dZ_O}{d\mu }}O_{ren}`$ (B2)
$`=`$ $`{\displaystyle \frac{\gamma _O}{2}}O_{ren}.`$ (B3)
Writing the solution to this equation as

$$Z_O\left(\mu ^2\right)=\frac{C_O\left(\mu _{}^{}{}_{}{}^{2}\right)}{C_O\left(\mu ^2\right)}Z_O\left(\mu _{}^{}{}_{}{}^{2}\right).$$
(B4)
and using the notation

$`\gamma _O`$ $`=`$ $`{\displaystyle \underset{i}{}}\gamma _O^{(i)}\left({\displaystyle \frac{\alpha _s}{4\pi }}\right)^{i+1}`$ (B5)
$`\overline{\gamma }_{Oi}`$ $`=`$ $`{\displaystyle \frac{\gamma _O^{(i)}}{2\beta _0}}`$ (B6)
$`\overline{\beta }_i`$ $`=`$ $`{\displaystyle \frac{\beta _i}{\beta _0}},`$ (B7)
gives rise to solutions to the running equation of the form (where we have suppressed the subscripts identifying the particular operator O):

* One loop solution :

$$C(\mu ^2)=\alpha _s(\mu )^{\overline{\gamma }_0}.$$
(B8)
* Two loop solution:

$$C(\mu ^2)=\alpha _s(\mu )^{\overline{\gamma }_0}\left\{1+\frac{\alpha _s(\mu )}{4\pi }\left(\overline{\gamma }_1\overline{\beta }_1\overline{\gamma }_0\right)\right\}.$$
(B9)
* Three loop solution:

$`C(\mu ^2)`$ $`=`$ $`\alpha _s(\mu )^{\overline{\gamma }_0}\{1+{\displaystyle \frac{\alpha _s(\mu )}{4\pi }}(\overline{\gamma }_1\overline{\beta }_1\overline{\gamma }_0)`$ (B11)
$`+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\alpha _s(\mu )}{4\pi }}\right)^2[(\overline{\gamma }_1\overline{\beta }_1\overline{\gamma }_0)^2+\overline{\gamma }_2+\overline{\beta }_1^2\overline{\gamma }_0\overline{\beta }_1\overline{\gamma }_1\overline{\beta }_2\overline{\gamma }_0]\}.`$

Tables V to VII show the anomalous dimensions used in this work. These values were taken from Refs.  with the number of flavours set to zero and the number of colours to three.

## C Matching Coefficients

The numerical values of the matching coefficients, $`Z_0^{(1)}`$ and $`Z_0^{(2)}`$ in Eq. 109, used for $`Z_q`$ and $`Z_S`$ are collected together in Table VIII and Table IX.