# Ratio of flavour non-singlet and singlet scalar density renormalisation
parameters in $N_{\mathrm{f}}=3$ QCD
with Wilson quarks

Jochen Heitger  Fabian Joswig  Pia L. J. Petrak  Anastassios Vladikas

###### Abstract

We determine non-perturbatively the normalisation factor $r_{\rm m}\equiv
Z_{\rm S}/Z_{\rm S}^{0}$, where $Z_{\rm S}$ and $Z_{\rm S}^{0}$ are the
renormalisation parameters of the flavour non-singlet and singlet scalar
densities, respectively. This quantity is required in the computation of quark
masses with Wilson fermions and for instance the renormalisation of nucleon
matrix elements of scalar densities. Our calculation involves simulations of
finite-volume lattice QCD with the tree-level Symanzik-improved gauge action,
$N_{\rm f}=3$ mass-degenerate ${\rm O}(a)$ improved Wilson fermions and
Schrödinger functional boundary conditions. The slope of the current quark
mass, as a function of the subtracted Wilson quark mass is extracted both in a
unitary setup (where nearly chiral valence and sea quark masses are
degenerate) and in a non-unitary setup (where all valence flavours are chiral
and the sea quark masses are small). These slopes are then combined with
$Z\equiv Z_{\rm P}/(Z_{\rm S}Z_{\rm A})$ in order to obtain $r_{\rm m}$. A
novel chiral Ward identity is employed for the calculation of the
normalisation factor $Z$. Our results cover the range of gauge couplings
corresponding to lattice spacings below $0.1\,$fm, for which $N_{\rm f}=2+1$
QCD simulations in large volumes with the same lattice action are typically
performed.

##  1 Introduction

Scalar and pseudoscalar flavour singlet and non-singlet dimension-3 bilinear
operators have the same anomalous dimension, since they belong to the same
chiral multiplet. The same is true for their renormalisation parameters,
provided that the regularisation does not break chiral symmetry. Otherwise,
the renormalisation parameters of the chiral multiplet components differ by
finite terms. This is the case for the lattice regularisation with Wilson
fermions. For example, the renormalisation parameters of the non-singlet
scalar and pseudoscalar densities (denoted as $Z_{\rm S}$ and $Z_{\rm P}$,
respectively) have a finite ratio which is a polynomial of the bare gauge
coupling $g_{0}$. This ratio can be determined by chiral Ward identities;111In
practice, distinct chiral Ward identities are used for the computation of the
ratio $Z_{\rm S}/(Z_{\rm P}Z_{\rm A})$ and $Z_{\rm A}$; the two results are
subsequently multiplied to give $Z_{\rm S}/Z_{\rm P}$. see Refs. [1, 2]. Since
$Z_{\rm P}$ and $Z_{\rm S}$ are scale dependent, imposing a renormalisation
scheme is necessary to fix one of them, and the other can be obtained using
the scheme independent ratio $Z_{\rm S}/Z_{\rm P}$.222Examples of
renormalisation schemes are ${\rm\overline{MS\kern-0.50003pt}\kern
0.50003pt}$, RI/(S)MOM [3, 4], the Schrödinger functional (SF) [5] and the
chirally rotated Schrödinger functional ($\chi$SF) [6]. In this way the
renormalised scalar and pseudoscalar densities are defined consistently in the
same scheme, with the same anomalous dimension and renormalisation group (RG)
running, and chiral symmetry is restored in the continuum limit. The ratio
$Z_{\rm S}/Z_{\rm P}$ has been computed for several gauge and Wilson fermion
actions (standard, improved etc.) in the quenched approximation [2, 7, 8, 9,
10, 11], with two dynamical quarks ($N_{\rm f}=2$ QCD) [12], and with three
dynamical quarks ($N_{\rm f}=3$ QCD) [13, 14, 15, 16].

Far less progress has been made on the computation of the ratio of the
renormalisation parameters of the non-singlet and singlet scalar densities,
$r_{\rm m}\equiv Z_{\rm S}/Z_{\rm S}^{0}$. For chirally symmetric
regularisations $r_{\rm m}=1$ holds, while for Wilson fermions $r_{\rm m}$ is
a (finite) polynomial of the gauge coupling, arising from the sea fermion
loops of the quark propagator. In the quenched approximation, $r_{\rm m}=1$.
As explained in Ref. [17], the lowest-order non-trivial perturbative
contribution to this quantity is a two-loop effect; i.e., $r_{\rm m}=1+{\rm
O}(g_{0}^{4})$. In Ref. [18] the ${\rm O}(g_{0}^{4})$ perturbative term has
been calculated for several lattice actions. Non-perturbative estimates of
this quantity have been reported in Ref. [13] at two values of the gauge
coupling for $N_{\rm f}=2+1$ QCD with the tree-level Symanzik-improved gauge
action [19] and the non-perturbatively improved Wilson-clover fermion action
[20]. This is the regularisation chosen by the CLS (Coordinated Lattice
Simulations) initiative which carries out QCD simulations with $N_{\rm f}=2+1$
flavours, on large physical volumes, for a range of bare couplings
corresponding to a hadronic regime [21, 22, 13, 23]. These CLS ensembles are
suitable for the computation of correlation functions, from which low-energy
hadronic quantities can be evaluated. In parallel, our group is performing
$N_{\rm f}=3$ simulations in the same range of bare gauge couplings, but for
small-volume lattices with Schrödinger functional boundary conditions and
nearly-chiral quark masses. These ensembles are used for the numerical
determination of the necessary renormalisation parameters and Symanzik
improvement coefficients, see Refs. [24, 15, 25, 26, 27, 28, 14] that have
various applications in lattice QCD when using this discretisation of Wilson
fermions. The present work provides high-precision estimates of $r_{\rm m}$
obtained in the same computational framework.

As seen from eq. (2.2) below, $(r_{\rm m}-1)$ contributes an ${\rm
O}(g_{0}^{4})$ term to the renormalisation of the quark masses [17]. This is
expected to be a small effect. Symanzik $\mathrm{O}(a)$ counterterms
containing $r_{\rm m}$ are often neglected in light quark mass determinations;
cf. Ref. [29]. In practical computations, however, $r_{\mathrm{m}}$ can be
relevant at $\mathrm{O}(a)$, especially when dealing with heavy flavours, and
should be taken into account in order to achieve full $\mathrm{O}(a)$
improvement; see, for example, eq. (2.13) in Ref. [30]. Another application
where $r_{\rm m}$ plays a prominent rôle is the nucleon sigma-term, which is
defined in terms of nucleon matrix elements of flavour singlet scalar
densities; see Refs. [31, 32] for example and [33, 34, 35] for more recent
works. A direct determination of $Z_{\rm S}^{0}$ is not as straightforward as
that of $Z_{\rm S}$, the former also requiring the computation of two-boundary
(“disconnected”) quark diagrams. This problem is circumvented by extracting
$Z_{\rm S}^{0}$ as the product of $Z_{\rm S}$ and $r_{\rm m}$.

Our computation of $r_{\rm m}$ is based on the relation between the current
(PCAC) mass $m$ and the subtracted quark mass $m_{\rm q}$. Close to the chiral
limit, $m(m_{\rm q})$ is a linear function with a slope that depends on the
details of the QCD model being simulated. In a unitary theory with degenerate
sea and valence quark masses, the slope of $m(m_{\rm q})$ is $Zr_{\rm m}$,
where $Z\equiv Z_{\rm P}/(Z_{\rm S}Z_{\rm A})$ and $Z_{\rm A}$ is the non-
singlet axial current normalisation. On the other hand, in a non-unitary
theory with chiral valence subtracted quark masses ($m_{\rm q}^{\rm val}=0$)
and small degenerate sea quark masses $m_{\rm q}^{\rm sea}\neq 0$, the slope
of $m(m_{\rm q}^{\rm sea})$ is $Z(r_{\rm m}-1)$. The two slopes are accessible
from two distinct sets of measurements at several common values of the bare
coupling $g_{0}$. The results are combined to give estimates of $r_{\rm
m}(g_{0}^{2})$. This approach is described in Section 2.

Alternatively, each of the two slopes $Zr_{\rm m}$ and $Z(r_{\rm m}-1)$ may be
combined with an independent estimate of $Z$, such as the results of Refs.
[14, 15]. In the present work we prefer to use a novel determination of $Z$,
relying on a chiral Ward identity which differs from the one of Ref. [15].
This identity is derived in Section 3.

In Section 4 we present our simulation setup for $N_{\rm f}=3$ QCD with
lattices of small physical volumes and Schrödinger functional boundary
conditions; these serve to numerically implement the strategies outlined in
the foregoing section. Most of our gauge field ensembles were already
generated in the context of previous works; cf. Refs. [15, 25, 26, 27, 28,
14]. Some new ensembles have also been generated, in order to cover the region
close to the origin of the function $m(m_{\rm q})$ more evenly and asses its
slope reliably.

Our results for $r_{\rm m}$, based on various combinations of $Zr_{\rm m}$,
$Z(r_{\rm m}-1)$, and $Z$ are discussed in Section 5. Different determinations
of $r_{\rm m}$ are compared, allowing us to settle for a conservative final
estimate with reliable systematic errors. Our final result is that of eq.
(5.5). In Table 5 we also list $r_{\rm m}(g_{0}^{2})$ for the
$g_{0}^{2}$-values at which CLS simulations are being performed for the
computations of hadronic quantities in $N_{\rm f}=2+1$ QCD.

In the final section we sum up our results and their uses in lattice QCD. More
detailed calculations and definitions of the correlation functions employed
can be found in Appendix A and B. Comparison of $Z$ determinations and
corresponding scaling tests can be found in Appendix C.

##  2 Wilson quark masses

In this section we recapitulate the basic quark mass definitions, namely
subtracted and current (PCAC) quark masses, and discuss how to obtain the
products $Zr_{\mathrm{m}}$ and $Z(r_{\mathrm{m}}-1)$ from relations between
the two. For any unexplained notation we refer to Ref. [14]. The starting
point is the subtracted bare quark mass of flavour $i=1,\ldots,N_{\rm f}$,

$m_{{\rm q},i}\equiv m_{0,i}-m_{\rm
crit}=\dfrac{1}{2a}\Big{(}\dfrac{1}{\kappa_{i}}-\dfrac{1}{\kappa_{\rm
crit}}\Big{)}\,,$ (2.1)

where $\kappa_{i}$ is the hopping parameter for flavour $i$, $\kappa_{\rm
crit}$ its value in the chiral limit, and $a$ is the lattice spacing. In terms
of the subtracted masses $m_{{\rm q},i}$, the corresponding renormalised quark
masses are given by

$\displaystyle m_{i,\rm R}=Z_{\rm m}\Bigg{[}m_{{\rm q},i}\,+\,(r_{\rm
m}-1)\dfrac{{\rm Tr}M_{\rm q}}{N_{\rm f}}\Bigg{]}+{\rm O}(a)\,,$ (2.2)

where $M_{\rm q}={\rm diag}(m_{{\rm q},1},\ldots,m_{{\rm q},N_{\rm f}})$ is
the $N_{\rm f}\times N_{\rm f}$ bare quark mass matrix.

We recall in passing that the renormalisation parameter $Z_{\rm
m}(g_{0}^{2},a\mu)$ depends on the renormalisation scale $\mu$ and diverges
logarithmically in the ultraviolet. It is the inverse of $Z_{\rm
S}(g_{0}^{2},a\mu)$, the renormalisation parameter of the flavour non-singlet
scalar density operator. A mass independent renormalisation scheme is implied
throughout this work. In such a scheme operator renormalisation parameters
(e.g. $Z_{\rm P},Z_{\rm m},Z_{\rm S}$), current normalisations (i.e. $Z_{\rm
A},Z_{\rm V}$) and $r_{\rm m}$ are functions of the squared bare gauge
coupling $g_{0}^{2}$. In a non-perturbative determination at non-zero quark
mass, they are affected by ${\rm O}(am_{{\rm q},i})$, ${\rm O}(a{\rm Tr}M_{\rm
q})$, and ${\rm O}(a\Lambda_{\rm QCD})$ discretisation effects, which are part
of their operational definition. As pointed out in Ref. [17], the term
$(r_{\rm m}-1)$ multiplies ${\rm Tr}M_{\rm q}$, so it arises from a mass
insertion in a quark loop. In perturbation theory it is a two-loop effect,
contributing at ${\rm O}(g_{0}^{4})$. Its non-perturbative determination is
the main purpose of this paper. An important consequence of eq. (2.2) is that
a renormalised mass $m_{i,\rm R}$ goes to the chiral limit only when all
subtracted masses $m_{{\rm q},1},\dots,m_{{\rm q},N_{\rm f}}$ vanish.

Alternatively, a bare current (PCAC) quark mass $m_{ij}$ can be defined
through the following relation:

$({\widetilde{\partial}_{\mu}})_{x}\big{\langle}(A_{{\rm
I}})^{ij}_{\mu}(x)\,\mathcal{O}^{ji}\big{\rangle}=2m_{ij}{\big{\langle}{P}^{ij}(x)\,\mathcal{O}^{ji}\big{\rangle}}\,.$
(2.3)

The quantity $m_{ij}$ is distinct from the subtracted bare quark masses, but
it is related to the mass average $(m_{{\rm q},i}+m_{{\rm q},j})/2$; see eq.
(2.9) below. The flavour non-singlet bare axial current and the pseudoscalar
density are given by

$A_{\mu}^{ij}(x)\equiv\bar{\psi}_{i}(x)\,\gamma_{\mu}\gamma_{5}\,\psi_{j}(x)\,,\qquad
P^{ij}(x)\equiv\bar{\psi}_{i}(x)\,\gamma_{5}\,\psi_{j}(x)\,,$ (2.4)

with indices $i,j$ denoting two distinct flavours ($i\neq j$). The
pseudoscalar density $P^{ij}$ and the current $(A_{{\rm I}})^{ij}_{\mu}\equiv
A^{ij}_{\mu}+ac_{\rm A}{\widetilde{\partial}_{\mu}}P^{ij}$ are Symanzik-
improved in the chiral limit, with the improvement coefficient $c_{\rm
A}(g_{0}^{2})$ being in principle only a function of the gauge coupling. In
these definitions, ${\widetilde{\partial}_{\mu}}$ denotes the average of the
usual forward and backward derivatives.333The forward derivative is defined as
$a\partial_{\mu}f(x)\equiv f(x+a\hat{\mu})-f(x)$ and the backward derivative
as $a\partial_{\mu}^{\ast}f(x)\equiv f(x)-f(x-a\hat{\mu})$. The source
operator $\mathcal{O}^{ji}$ is defined in a region of space-time that does not
include the point $x$, so as to avoid contact terms. In the ${\rm O}(a)$
improved theory, the renormalised axial current and pseudoscalar density are

$\displaystyle(A_{{\rm R}})^{ij}(x)$ $\displaystyle=Z_{\rm
A}(g_{0}^{2})(A_{\rm I})_{\mu}^{ij}(x)+{\rm O}(am_{\mathrm{q}},a^{2})\,,$
(2.5) $\displaystyle(P_{{\rm R}})^{ij}(x)$ $\displaystyle=Z_{\rm
P}(g_{0}^{2},a\mu)P^{ij}(x)+{\rm O}(am_{\mathrm{q}},a^{2})\,.$ (2.6)

The normalisation of the axial current $Z_{\rm A}(g_{0}^{2})$ is scale
independent, depending only on the squared gauge coupling $g_{0}^{2}$. The
renormalisation parameter $Z_{\rm P}(g_{0}^{2},a\mu)$ (determined, say in the
Schrödinger functional scheme of Ref. [5]) additionally depends on the
renormalisation scale $\mu$ and diverges logarithmically in the ultraviolet.
The PCAC relation, expressed by renormalised fields,

$\displaystyle({\widetilde{\partial}_{\mu}})_{x}\big{\langle}\,(A_{{\rm
R}})^{ij}_{\mu}(x)\ \mathcal{O}^{ji}\,\big{\rangle}$ $\displaystyle=(m_{{\rm
R},i}+m_{{\rm R},j})\,\big{\langle}\,(P_{\mathrm{R}})^{ij}(x)\
\mathcal{O}^{ji}\,\big{\rangle}\,,$ (2.7)

valid up to discretisation effects in the continuum, combined with eqs.
(2.3)–(2.6), implies that

$\dfrac{m_{i,\rm R}+m_{j,\rm R}}{2}=\dfrac{Z_{\rm A}}{Z_{\rm P}}m_{ij}+{\rm
O}(am_{\mathrm{q}},a^{2})\,.$ (2.8)

If we calculate the average mass $(m_{i\rm R}+m_{j\rm R})/2$ from eq. (2.2)
and equate the result to the r.h.s of eq. (2.8), we obtain an expression which
relates subtracted and PCAC bare masses:

$m_{ij}=Z\Bigg{[}\dfrac{(m_{{\rm q},i}+m_{{\rm q},j})}{2}+(r_{\rm
m}-1)\dfrac{{\rm Tr}M_{\rm q}}{N_{\rm f}}\Bigg{]}+{\rm
O}(am_{\mathrm{q}},a^{2})\,,$ (2.9)

where the product of the renormalisation parameters $Z(g_{0}^{2})\equiv Z_{\rm
P}(g_{0}^{2},\mu)/(Z_{\rm S}(g_{0}^{2},\mu)Z_{\rm A}(g_{0}^{2}))$ is scale
independent. We now exploit eq. (2.9) in two ways:

(1) In a theory with mass-degenerate quarks ($m_{{\rm q},i}=m_{{\rm q},j}={\rm
Tr}M_{\rm q}/N_{\rm f}$), it reduces to

$\displaystyle m$ $\displaystyle=Zr_{\rm m}m_{\rm q}+{\rm
O}(am_{\mathrm{q}},a^{2})$ (2.10) $\displaystyle=Zr_{\rm
m}\dfrac{1}{2a}\Big{(}\dfrac{1}{\kappa}-\dfrac{1}{\kappa_{\rm
crit}}\Big{)}+{\rm O}(am_{\mathrm{q}},a^{2})\,.$ (2.11)

In the above equation, flavour indices have been dropped from the quark masses
$m_{ij},m_{{\rm q},i}$ and the hopping parameter $\kappa_{i}$. This
simplification of notation will be adopted on most occasions below. Thus,
modelling the current quark mass $am$ as a function of $1/\kappa$ for values
of $\kappa$ close to $\kappa_{\rm crit}$, we obtain the latter as the root of
the function $am(1/\kappa)$ and the combination $Zr_{\rm m}$ as the slope of
the same curve.

(2) Once the critical hopping parameter $\kappa_{\rm crit}$ is available from
the previous step (1), we use a non-unitary setup where valence and sea quarks
of the same flavour have different bare subtracted masses $m_{{\rm q},i}^{\rm
val}\neq m_{{\rm q},i}^{\rm sea}$. In eq. (2.9), masses $m_{{\rm q},i}$ and
$m_{{\rm q},j}$ on the r.h.s. are valence quark contributions, while ${\rm
Tr}M_{\rm q}$ stands for the trace of sea quark masses; see Refs. [17, 14] for
detailed explanations. In particular, we set $\kappa^{\rm val}=\kappa_{\rm
crit}$, so as to ensure $m_{\rm q}^{\rm val}=0$ for all valence flavours.
Moreover sea quark masses are taken to be small, degenerate, and non-zero
(i.e. $\kappa^{\rm sea}\neq\kappa_{\rm crit}$, ensuring $m_{\rm q}^{\rm
sea}\neq 0$ for all sea flavours). With these conditions, the current quark
mass of eq. (2.9) reduces to

$m=Z(r_{\rm m}-1)m_{\rm q}^{\rm sea}+{\rm O}(am_{\mathrm{q}},a^{2})\,.$ (2.12)

It is remarkable that with non-zero bare subtracted sea quark masses (i.e.
$m_{\rm q}^{\rm sea}\neq 0$), all current quark masses in this setup are not
chiral (i.e. $m_{ij}\neq 0,\forall i,j$), even if all subtracted valence quark
masses vanish (i.e. $m_{{\rm q},i}^{\rm val}=0,\forall i$). From eq. (2.12) we
see that, if we compute $am$ as a function of $am_{\rm q}^{\rm sea}$ for
several sea masses, the slope of the functions gives an estimate of $Z(r_{\rm
m}-1)$.

The two slopes $Zr_{\rm m}$ and $Z(r_{\rm m}-1)$, computed in the two
different settings described above, but at the same gauge couplings
$g_{0}^{2}$, can be combined yielding estimates of $r_{\rm m}(g_{0}^{2})$; see
Subsection 5.3 for details. We stress that the above discussion concerns
relations which suffer from ${\rm O}(a)$ discretisation effects. For the quark
masses, such effects may be removed by introducing Symanzik counterterms,
leaving us with ${\rm O}(a^{2})$ discretisation errors. These counterterms
have been worked out in Refs. [36, 17]. In Ref [17] (see also eq. (2.10) of
Ref. [14]) the full ${\rm O}(am_{\rm q})$ contributions, omitted in eq. (2.9)
above, are written down explicitly. Such contributions are complicated and
taking them all into account could compromise the numerical stability of our
procedure to extract the quantities in question. We prefer a simpler and more
robust strategy, consisting of working with small quark masses so that ${\rm
O}(a)$-effects in eq. (2.9) may be safely dropped. This must of course be
checked a posteriori, by ensuring that the function $m(m_{\rm q})$ is linear
close to the origin, where our simulations are performed. The only improvement
coefficients used in this work are $c_{\mathrm{sw}}$ of the clover action and
$c_{\rm A}$, of the axial current (entering the PCAC mass).

There is an important subtlety concerning results obtained with Wilson
fermions in a Symanzik-improved setup: the bare parameters of the theory (i.e.
the gauge coupling $g_{0}^{2}$ and the $N_{\rm f}=2+1$ quark masses) are to be
varied, while staying on lines of constant physics within systematic
uncertainties of ${\rm O}(a^{2})$. In particular, if the improved bare gauge
coupling [36]

$\tilde{g}_{0}^{2}\equiv g_{0}^{2}\Big{(}1+\dfrac{1}{N_{\rm
f}}b_{g}(g_{0}^{2})a{\rm Tr}M_{\rm q}\Big{)}$ (2.13)

is kept fixed in the simulations, so is the lattice spacing, with fluctuations
being attributed to ${\rm O}(a^{2})$ effects [21]. This implies that, once
$\kappa^{\mathrm{crit}}$ has been evaluated as a function of $g_{0}^{2}$,
(re)normalisation parameters and improvement coefficients should be treated as
functions of $\tilde{g}_{0}^{2}$, rather than $g_{0}^{2}$; e.g. $Z_{\rm
A}(\tilde{g}_{0}^{2}),Z_{\rm
P}(\tilde{g}_{0}^{2},a\mu),Z(\tilde{g}_{0}^{2}),r_{\rm m}(\tilde{g}_{0}^{2})$
etc. To the extent that we are working in the chiral limit, or very close to
it (i.e. very light quark masses), this difference is immaterial. This is why
in the present work we always express our results as functions of $g_{0}^{2}$.
However, when they are to be used away from the chiral limit at low-energy
scales (see Refs. [22, 29]), this difference must be taken into account
properly. We shall elaborate further on this point when summarising our work
in Section 6.

##  3 Ward identity determination of $Z$

In the previous section we have shown how the quantities $Zr_{\rm m}$ and
$Z(r_{\rm m}-1)$ can be estimated from relations between suitably chosen
current and subtracted Wilson quark masses. They may then straightforwardly be
combined to give $r_{\rm m}$ and $Z$. The latter quantity has already been
measured in our setup ($N_{\rm f}=3$ lattice QCD with Schrödinger functional
boundary conditions) in two ways: either by using appropriate combinations of
current and subtracted quark masses with different flavours [14], or from
chiral Ward identities [15] via $Z\equiv Z_{\rm P}/(Z_{\rm A}Z_{\rm S})$. Here
we will describe yet another direct method, based on a new Ward identity, very
similar to the one of Ref. [15]. The reader is referred to that work for
details, notation etc.

We consider a product of two composite operators ${\cal O}\equiv S^{b}(y){\cal
O}^{c}$, defined as

$\displaystyle S^{b}(y)$ $\displaystyle\equiv$
$\displaystyle\mathrm{i}\bar{\psi}(y)T^{b}\psi(y)$ $\displaystyle{\cal O}^{c}$
$\displaystyle\equiv$ $\displaystyle\mathrm{i}\dfrac{a^{6}}{L^{3}}\sum_{\bf
u,v}\bar{\zeta}({\bf u})\gamma_{5}T^{c}\zeta({\bf v})\,,$ (3.1)

where $T^{b}$ and $T^{c}$ are generators of $SU(N_{\mathrm{f}})$. The former
operator is the flavour non-singlet scalar density, located in the bulk of
space-time, while the latter resides at the $x_{0}=0$ Dirichlet time boundary
of the Schrödinger functional.444For reasons of convenience, we have adopted a
slightly different notation in this section: the flavour content of operators
like $S^{b}$ or ${\cal O}^{c}$ is determined by a single flavour index $b$ or
$c$, corresponding to its flavour matrix $T^{b}$ or $T^{c}$. The fermion
fields of these operators $\psi$ and $\bar{\psi}$ are columns in flavour
space. This is to be contrasted to the notation of Section 2, where we have
introduced operators like $P^{ij}$ and $\mathcal{O}^{ji}$, which have explicit
indices, referring to the flavour of fields $\psi_{j},\bar{\psi}_{i}$ etc. The
Ward identity of interest is obtained by performing axial variations on ${\cal
O}$ in a region $R$, chosen to be the space-time volume between the hyper-
planes at $t_{1}$ and $t_{2}$ where $t_{1}<t_{2}$. With ${\cal O}^{c}$ lying
outside $R$, we have $\delta_{\rm A}{\cal O}=[\delta_{\rm A}S^{b}(y)]{\cal
O}^{c}$ and

$\delta_{\rm
A}S^{b}(x)=\epsilon^{a}\Big{[}d^{abe}P^{e}(x)+\dfrac{\delta^{ab}}{N_{\rm
f}}\bar{\psi}(x)\psi(x)\Big{]}\,.$ (3.2)

In what follows we simplify matters by always working with $a\neq b$, so as to
eliminate the second contribution on the r.h.s. of the above expression. In
analogy to the derivation exposed in Ref. [15], we arrive at the formal
continuum Ward identity

$\displaystyle\int\mathrm{d}^{3}{\bf y}\int\mathrm{d}^{3}{\bf
x}\Big{\langle}\Big{[}A_{0}^{a}(t_{2};{\bf x})-A_{0}^{a}(t_{1};{\bf
x})\Big{]}S^{b}(y_{0};{\bf y}){\cal O}^{c}\Big{\rangle}$
$\displaystyle-2m\int\mathrm{d}^{3}{\bf y}\int\mathrm{d}^{3}{\bf
x}\int_{t_{1}}^{t_{2}}\mathrm{d}x_{0}\langle P^{a}(x_{0};{\bf
x})S^{b}(y_{0};{\bf y}){\cal O}^{c}\rangle$ (3.3)
$\displaystyle=-d^{abe}\int\mathrm{d}^{3}{\bf y}\,\,\langle P^{e}(y){\cal
O}^{c}\rangle\,.$

Next we adapt the previous formal manipulations to the lattice regularisation
with Schrödinger functional boundary conditions. The pseudoscalar operator
${\cal O}^{c}$ is defined on the $x_{0}=0$ time boundary. Ward identity (3.3)
then becomes:

$\displaystyle Z_{\rm A}Z_{\rm S}a^{6}\Bigg{\\{}\sum_{{\bf x},{\bf
y}}\,\left\langle\Big{[}(A_{\rm I})^{a}_{0}(t_{2};{\bf x})-(A_{\rm
I})^{a}_{0}(t_{1};{\bf x})\Big{]}\,S^{b}(y_{0};{\bf y})\,{\cal
O}^{c}\right\rangle$ $\displaystyle-2am\sum_{{\bf x},{\bf
y}}\sum_{x_{0}=t_{1}}^{t_{2}}w(x_{0})\,\langle P^{a}(x_{0};{\bf
x})\,S^{b}(y_{0};{\bf y})\,{\cal O}^{c}\rangle\Bigg{\\}}$ (3.4)
$\displaystyle=-d^{abe}Z_{\rm P}\,\,a^{3}\sum_{\bf y}\langle\,P^{e}(y)\,{\cal
O}^{c}\rangle+{\rm O}(am,a^{2})\,.$

In this expression repeated flavour indices $e$ are summed, as usual. The
weight factor is $w(x_{0})=1/2$ for $x_{0}\in\\{t_{1},t_{2}\\}$ and
$w(x_{0})=1$ otherwise. It is introduced in order to implement the trapezoidal
rule for discretising integrals. Quark masses are degenerate and $m$ is the
current quark mass.

The last step is to perform the Wick contractions in Ward identity (3.4). How
this is done is explained in Appendix B; eventually, flavour factors drop out
and we are left with a Ward identity that translates into traces of products
of quark propagators and $\gamma$-matrices, graphically depicted in Fig. 1.
Solving for $Z$ we get

$\displaystyle Z\equiv\dfrac{Z_{\rm P}}{Z_{\rm A}Z_{\rm S}}={-}\dfrac{f_{\rm
AS}^{\mathrm{I}}(t_{2},y_{0})-f_{\rm
AS}^{\mathrm{I}}(t_{1},y_{0})-2am\tilde{f}_{\rm PS}(t_{2},t_{1},y_{0})}{f_{\rm
P}(y_{0})}+{\rm O}(am,a^{2})\,,$ (3.5)

where dependencies are suppressed on the l.h.s. Assuming that we work in the
chiral limit (or with nearly-vanishing quark masses, so that ${\rm O}(am)$
effects may be safely neglected), the above Ward identity is valid up to ${\rm
O}(a^{2})$ discretisation errors in lattice QCD with Wilson quarks. In this
spirit, terms proportional to Symanzik $b$-coefficients may also be safely
ignored.555This is even true for light (up/down, strange) non-chiral quark
masses, as explicitly demonstrated in Ref. [29], using the $b$-coefficients of
Ref. [14]. The renormalisation factor of the external source ${\cal O}^{c}$ is
not taken into consideration, as it cancels out in the ratio (3.5). The term
proportional to the current quark mass $m$ may also be dropped close to the
chiral limit, but since we are working with masses which are not strictly
zero, it could be advantageous to keep it in practice. In fact, it was found
in Refs. [26, 15] that this term stabilizes the chiral extrapolation leading
to smaller errors. This turns out to be true also in our case, as we will show
in Subsection 5.2 and Fig. 5.

(a) Diagram $f_{\rm P}$

(b) Diagram $f_{\rm AS;1}$

(c) Diagram $f_{\rm AS;2}$

Figure 1: The trace diagrams contributing to the expectation values of $f_{\rm
P}$, defined in eq. (B.1) (diagram (a)) and $f_{\rm AS}$, defined in eq. (B.3)
(diagrams (b) and (c)). The wall represents the time slice $x_{0}=0$ with a
$\gamma_{5}$ Dirac matrix between circles. The squares in the bulk represent
either the insertions of a pseudoscalar operator $P(y)$ (diagram (a)) or a
scalar operator $S(y)$ (diagrams (b) and (c)). The diamonds stand for an axial
operator $A_{0}(x)$. The open circles correspond to the boundary fields
$\zeta$, while the filled circles denote $\bar{\zeta}$. The diagrams
schematically represent traces, formed by starting from any point and
following the lines (quark propagators) until we close the loop. The time
ordering of points $x$ and $y$ is left unspecified in these diagrams.

It is interesting to compare Ward identity (3.3) with those of Ref. [15]:

  * •

In Ref. [15] the flavour factors gave rise to a multitude of identities, which
were combined in order to increase the signal-to-noise ratio, while here we
only have one identity. On these grounds one could expect that the numerical
results of Ref. [15] are more precise than the ones from the Ward identity
introduced here.

  * •

On the other hand, the identities of Ref. [15] involved: (i) correlation
functions with one operator insertion in the bulk of the lattice and one wall
source at each time slice; cf. Fig. 1 in that work; (ii) correlation functions
with two operator insertions in the bulk and one wall source at each time
slice; cf. Fig. 2 in that work. Here we have: (i) a correlation function with
one operator insertion in the bulk and one wall source; (ii) correlation
functions with two operator insertions in the bulk and one wall source. These
somewhat simpler correlation functions illustrated in Fig. 1 above are
expected to have less statistical fluctuations. From this point of view, the
results of the present work are expected to gain in accuracy.

Thus, one of our aims is to establish which of the two approaches leads to
more accurate results. This is discussed in Subsection 5.2 and Appendix C.

##  4 Numerical setup

$(L/a)^{3}\times T/a$ | $\beta$ | $\kappa$ | #REP | #MDU | ID |  $a$ (in fm)
---|---|---|---|---|---|---
$12^{3}\times 17$ | 3.3 | 0.13652 | 20 | 20480 | A1k1 | $0.1045(18)$
|  | 0.13648 | 5 | 6876 | A1k3 |
|  | 0.13650 | 20 | 96640 | A1k4 |
$12^{3}\times 18$ |  | $0.13612$ | 4 | 41600 | A3k1 |
|  | $0.13627$ | 4 | 41600 | A3k2 |
|  | $0.13593$ | 4 | 41600 | A3k3 |
|  | $0.136444$ | 4 | 41600 | A3k4 |
|  | $0.136575$ | 4 | 41600 | A3k5 |
|  | $0.136385$ | 4 | 41600 | A3k6 |
$14^{3}\times 21$ | 3.414 | 0.13690 | 32 | 38400 | E1k1 | $0.08381(68)$
|  | 0.13695 | 48 | 57600 | E1k2 |
$14^{3}\times 20$ |  | 0.13656 | 18 | 60480 | E2k1 |
|  | 0.13675 | 18 | 60480 | E2k2 |
$16^{3}\times 23$ | 3.512 | 0.13700 | 2 | 20480 | B1k1 | $0.06954(43)$
|  | 0.13703 | 1 | 8192 | B1k2 |
|  | 0.13710 | 2 | 16384 | B1k3 |
|  | 0.13714 | 1 | 27856 | B1k4 |
$16^{3}\times 24$ |  | 0.13677 | 1 | 25904 | B3k1 |
$20^{3}\times 29$ | 3.676 | 0.13680 | 1 | 7848 | C1k1 | $0.05170(42)$
|  | 0.13700 | 4 | 15232 | C1k2 |
|  | 0.13719 | 4 | 15472 | C1k3 |
$24^{3}\times 35$ | 3.810 | 0.13711875582 | 5 | 8416 | D1k1∗ | $0.04175(70)$
|  | 0.13701 | 2 | 6424 | D1k2 |
|  | 0.137033 | 8 | 85008 | D1k4 |
Table 1: Simulation parameters $L$, $T$, $\beta$, $\kappa$, the number of
replica #REP and the number of molecular dynamics units #MDU for the ensembles
labelled by ID. Ensembles highlighted in italics were newly generated for this
study while the remaining ones were already used in previous investigations
(see, for example Ref. [14]). The ensemble D1k1 marked by an asterisk is only
used for the determination of the PCAC masses. The lattice spacings $a$ are
obtained by interpolating the results of Ref. [22] with a polynomial fit. All
configurations are separated by 8 MDU’s except for the ensembles A1k3 (4
MDU’s) and D1k4 (16 MDU’s).

We employ the tree-level Symanzik-improved gauge action and $N_{\mathrm{f}}=3$
mass-degenerate $\mathrm{O}(a)$ improved Wilson fermions. For the
corresponding improvement coefficient $c_{\mathrm{sw}}$ we use the non-
perturbative determination of Ref. [37]. As already indicated, we impose
Schrödinger functional boundary conditions at the temporal boundaries of the
lattice. The Schrödinger functional setup is highly suitable for massless
renormalisation schemes, since nearly-vanishing quark masses are accessible in
numerical calculations due to the spectral gap of the Dirac operator. This gap
is imposed by the boundaries, so that the quark mass dependence can be mapped
out reliably in the vicinity of the chiral point. The generation of the gauge
field configurations is performed with the openQCD code [38] which employs the
RHMC algorithm [39, 40] for the third quark.

All gauge field ensembles used in this study are summarized in Table 1 and lie
on a line of constant physics (LCP), defined by a fixed spatial extent of
$L\approx 1.2\,$\mathrm{f}\mathrm{m}$$ and $T/L\approx 3/2$. The tuning was
guided by the two-loop beta-function; see Ref. [27]. Provided that this
perturbative approximation is satisfactory in the case at hand, this ensures
that our estimates of $r_{\rm m}$ and $Z$ become smooth functions of the
lattice spacing, with higher-order ambiguities vanishing monotonically. In
Ref. [25] it was explicitly shown that $L$ is constant up to ${\rm O}(a)$ cut-
off effects across the coupling range also considered in the present work. We
thus expect our final results for $r_{\rm m}$ and $Z$ to only be affected by
${\rm O}(a^{2})$ effects. These are beyond the order we are interested in and
they are treated as an ambiguity that extrapolates to zero in the continuum
limit.666More precisely, the results on scale setting for our lattice action
from Ref. [22] have been used in Ref. [25] in order to demonstrate numerically
that the deviation from a constant value of $L$ in physical units is
proportional to the lattice spacing $a$. As the latter work uses the
configuration ensembles and range of non-perturbative bare couplings used also
in the present paper, our simulation parameters define a LCP up to ${\rm
O}(a)$ lattice artefacts, so that the discretisation effects of $r_{\rm m}$
and $Z$ are ${\rm O}(a^{2})$ in the ${\rm O}(a)$ improved theory.

The gauge ensembles highlighted in italics were newly generated for this
study, while the remaining ones were already used in previous investigations;
see Refs. [24, 15, 25, 26, 27, 28, 14].777In a setup with heavy sea quarks and
very light valence quarks we approach a quenched-like situation in which
exceptional configurations are to be expected; cf. Ref. [41] where a similar
situation is discussed. In a careful analysis we identified only one gauge
field configuration in the ensemble E2k1, with an exceptionally small
eigenvalue of the massless Dirac operator. This leads to very large values of
the correlation functions $f_{\mathrm{P}}$ and $f_{\mathrm{A}}$. We have
discarded this exceptional configuration. These additional ensembles allow for
a more even and wider spread of bare quark masses around the chiral point for
each value of $\beta$, which enables a more precise extraction of the slopes
corresponding to $Zr_{\mathrm{m}}$ and $Z\left(r_{\mathrm{m}}-1\right)$ as
explained in Section 2. Since a newer version of the openQCD code was utilised
for the generation of the ensembles, the time extent $T/a$, which was odd in
the pre-existing ensembles, is even for the new ones. For all ensembles we use
tree-level boundary $\mathrm{O}(a)$ improvement for both the gauge and fermion
fields (i.e. the appropriate $c_{\rm t},\widetilde{c}_{\rm t}$ values) as if
the time extents were even. The fact that an odd time extent alters the tree-
level value of $c_{\rm t}$, depending on the definition of the line of
constants physics [42], affects the current quark masses below the precision
achieved here, as explicitly demonstrated in Ref. [27].

All Schrödinger functional correlation functions required for our numerical
investigations are $\mathrm{O}(a)$ improved. In this context we only require
the improvement coefficient $c_{\mathrm{A}}$, non-perturbatively known from
Ref. [27]. Since the Markov chain Monte Carlo sampling of the gauge field
configurations suffers from critical slowing down of the topological charge
for smaller lattice spacings (see Ref. [43]), we project our data to the
trivial topological sector as suggested in Ref. [44], in order to account for
the insufficient sampling of all topological sectors. For the analysis of the
statistical errors we employ the $\Gamma$-method [45]. We account for the
remaining critical slowing down of the Monte Carlo algorithm by attaching a
tail to the autocorrelation function, as suggested in Ref. [46]. The
corresponding slowest mode is estimated from the autocorrelation time of the
boundary-to-boundary correlation function $F_{1}^{ij}$, defined in Appendix A.
The error analysis is carried out with a python implementation of the
$\Gamma$-method, using automatic differentiation for the error propagation as
proposed in Ref. [47].

##  5 Analysis details and results

In the following we present our analysis which eventually leads to several
estimates for the ratio of the renormalisation parameters of the non-singlet
and singlet scalar densities, $r_{\rm m}$. We will first describe how we
obtain $Zr_{\mathrm{m}}$, $Z(r_{\mathrm{m}}-1)$, and $Z$ individually and then
discuss several ways of combining the three into $r_{\mathrm{m}}$. As a final
result we provide an interpolation formula for $r_{\mathrm{m}}$ and extract
its value at the bare couplings of large-volume CLS simulations [21, 13, 23].

###  5.1 Quark mass slopes

Figure 2: PCAC mass $am$ as a function of time $x_{0}/a$ for ensembles B3k1
(left) and D1k2 (right). Squares are results obtained in the unitary setup,
while diamonds are results obtained in the non-unitary setup. The final
estimate for $m$ is obtained by averaging results in the time interval
$[T/3,2T/3]$, indicated by the dashed vertical lines.

As described in Section 2, the quantities $Zr_{\mathrm{m}}$ and
$Z(r_{\mathrm{m}}-1)$ can be extracted from quark mass slopes. Our results are
based on the determination of the $\mathrm{O}(a)$ improved PCAC masses via

$\displaystyle m(x_{0})=\frac{{\widetilde{\partial}_{0}}f_{\rm
A}^{ij}(x_{0})+ac_{\mathrm{A}}\partial_{0}^{\ast}\partial_{0}f_{\rm
P}^{ij}(x_{0})}{2f_{\rm P}^{ij}(x_{0})}\,,$ (5.1)

where $f_{\mathrm{A}}^{ij}$ and $f_{\mathrm{P}}^{ij}$ are Schrödinger
functional correlation functions. In order to improve the signal, these
correlation functions are symmetrised with their $T$-symmetric counterparts
$g_{\rm A}^{ij}(T-x_{0})$ and $g_{\rm P}^{ij}(T-x_{0})$, which are constructed
from the same operators $(A_{{\rm I}})^{ij}_{0}(x)$ and ${P}^{ij}(x)$ in the
bulk but the pseudoscalar wall with operator $\mathcal{O^{\prime}}^{ji}$
positioned at the time boundary $x_{0}=T$. For exact definitions see Appendix
A.

ID | $am$ | ${Z^{\\{T/3\\}}}$ | ${Z^{\\{T/4\\}}}$
---|---|---|---
| $\kappa^{\mathrm{val}}_{i}=\kappa_{i}^{\mathrm{sea}}$ | $\kappa^{\mathrm{val}}_{i}=\kappa_{\mathrm{crit}}$ |  |
A3k3 | $\phantom{-}0.12143(82)$ | $\phantom{-}0.10759(77)$ |  |
A3k1 | $\phantom{-}0.07316(184)$ | $\phantom{-}0.06440(193)$ |  |
A3k2 | $\phantom{-}0.03070(85)$ | $\phantom{-}0.02588(89)$ |  |
A3k6 | $\phantom{-}0.01246(54)$ | $\phantom{-}0.01029(54)$ |  |
A3k4 | $\phantom{-}0.00465(56)$ | $\phantom{-}0.00370(57)$ |  |
A1k3 | $\phantom{-}0.00095(93)$ | $\phantom{-}0.00074(93)$ | 0.8195(93) | 0.7454(94)
A1k4 | $-0.00119(33)$ | $-0.00100(33)$ | 0.8101(43) | 0.7520(58)
A1k1 | $-0.00287(61)$ | $-0.00229(61)$ | 0.7892(67) | 0.7189(79)
A3k5 | $-0.00952(50)$ | $-0.00864(49)$ |  |
| $\phantom{-}0.0$ | $\phantom{-}0.0$ | 0.8184(77) | 0.7588(143)
E2k1 | $\phantom{-}0.02083(19)$ | $\phantom{-}0.01117(27)$ |  |
E2k2 | $\phantom{-}0.01072(16)$ | $\phantom{-}0.00592(17)$ |  |
E1k1 | $\phantom{-}0.00265(22)$ | $\phantom{-}0.00153(23)$ | 0.8990(47) | 0.8619(54)
E1k2 | $-0.00022(19)$ | $-0.00017(19)$ | 0.8987(47) | 0.8580(64)
| $\phantom{-}0.0$ | $\phantom{-}0.0$ | 0.8987(43) | 0.8583(59)
B3k1 | $\phantom{-}0.01502(16)$ | $\phantom{-}0.00552(22)$ |  |
B1k1 | $\phantom{-}0.00552(19)$ | $\phantom{-}0.00232(18)$ | 0.9972(45) | 0.9760(53)
B1k2 | $\phantom{-}0.00435(28)$ | $\phantom{-}0.00168(30)$ | 0.9963(73) | 0.9756(94)
B1k3 | $\phantom{-}0.00157(18)$ | $\phantom{-}0.00024(20)$ | 0.9839(48) | 0.9643(52)
B1k4 | $-0.00056(16)$ | $-0.00035(16)$ | 1.0004(50) | 0.9690(73)
| $\phantom{-}0.0$ | $\phantom{-}0.0$ | 0.9935(38) | 0.9654(50)
C1k1 | $\phantom{-}0.01322(17)$ | $\phantom{-}0.00304(21)$ | 1.0593(46) | 1.0446(42)
C1k2 | $\phantom{-}0.00601(11)$ | $\phantom{-}0.00148(11)$ | 1.0615(30) | 1.0517(35)
C1k3 | $-0.00110(11)$ | $-0.00029(11)$ | 1.0617(47) | 1.0542(42)
| $\phantom{-}0.0$ | $\phantom{-}0.0$ | 1.0621(36) | 1.0544(34)
D1k2 | $\phantom{-}0.00073(15)$ | $\phantom{-}0.00012(15)$ | 1.0896(89) | 1.0868(52)
D1k4 | $-0.00007(3)$ | $-0.00001(3)$ | 1.0908(12) | 1.0849(13)
D1k1 | $-0.00295(11)$ | $-0.00040(9)$ |  |
| $\phantom{-}0.0$ | $\phantom{-}0.0$ | 1.0907(13) | 1.0850(12)
Table 2: For each ensemble, identified in the first column by an ID label, we
list our results for the PCAC mass $am$ for simulations with
$\kappa^{\mathrm{val}}=\kappa^{\mathrm{sea}}$ (second column) and
$\kappa^{\mathrm{sea}}\neq\kappa^{\mathrm{val}}=\kappa_{\mathrm{crit}}$ (third
column). The last two columns contain $Z$ results obtained from the Ward
identity (3.5). The final results are those extrapolated to the chiral limit
at each $\beta=6/g_{0}^{2}$ (last line of each data grouping). The labels
${Z^{\\{T/3\\}}}$ and ${Z^{\\{T/4\\}}}$ refer to different choices of time
slices with operator insertions in the correlation functions (see text for
details).

We first determine the required correlation functions in a unitary setup,
$\kappa^{\mathrm{val}}=\kappa^{\mathrm{sea}}$. From these we can obtain
$\kappa_{\mathrm{crit}}$ as will be detailed below. In a second step we
compute the same correlation functions in a non-unitary setup where
$\kappa^{\mathrm{sea}}\neq\kappa^{\mathrm{val}}=\kappa_{\mathrm{crit}}$. In
Fig. 2 we show the temporal dependence of the current quark mass $m(x_{0})$
for both of these setups for the representative ensembles B3k1 and D1k2 and
demonstrate that they form well-defined plateaux as a function of time, away
from the Dirichlet boundaries. Our final estimate for the PCAC masses is
obtained by averaging $m(x_{0})$ over the central third of the temporal extent
of the lattice. This choice is motivated by the coarsest lattices; the
plateaux for the finer ones also extend closer to the boundary before lattice
artefacts become relevant as can be seen in Fig. 2. The plateau range is
adapted according to the time extent for each value of $\beta$, so as to
preserve the line of constant physics. Our PCAC mass estimates in both setups
are listed in Table 2 for all ensembles.

Figure 3: PCAC masses $am$ fitted linearly in $1/(2\kappa)$, for all simulated
$\beta$ values (i.e. for decreasing lattice spacings from top to bottom). Open
squares and filled diamonds are results in the unitary and non-unitary setups,
respectively. Note that horizontal and vertical axes are identical for all
values of $\beta$, so as to highlight the different ranges of $\kappa$ and the
change of $\kappa_{\mathrm{crit}}$ marked by the vertical dashed lines.

In order to extract $Zr_{\mathrm{m}}$ and $Z(r_{\mathrm{m}}-1)$ from the
slopes of the current quark masses with respect to the bare quark masses, we
plot $am$ against the inverse hopping parameter $1/(2\kappa)$ for both the
unitary and the non-unitary setup, as demonstrated in Fig. 3. We generally
observe that $m$ behaves linearly as a function of $1/(2\kappa)$ in the range
$-0.1\lesssim Lm\lesssim 0.3$. For the ensembles A3k1, A3k2, and A3k3 (not
displayed in Fig. 3), which correspond to $Lm\gtrsim 0.3$, linearity is lost.

Figure 4: PCAC masses $am$, at $\beta=3.3$, fitted with a quadratic polynomial
in $1/(2\kappa)$. Squares and diamonds are results in the unitary and non-
unitary setups, respectively. Note that the two rightmost points (A3k1, A3k3)
are not included in the fit, while A3k2 (at $1/(2\kappa)=3.6692$) is. The
vertical dashed line is positioned at $\kappa_{\mathrm{crit}}$ from the linear
fit; see Table 3. The $Zr_{\mathrm{m}}$ and $Z(r_{\mathrm{m}}-1)$ values shown
in the legend are obtained from the linear term of the quadratic polynomial.

Results from these ensembles have thus not been included in the linear fits.
The good linear behaviour of the data from the remaining ensembles is
justified a posteriori, by the small $\chi^{2}/\mathrm{d.o.f}$. of our fits,
as shown in Fig. 3.

We also probe the non-linear regime in both setups for $\beta=3.3$ by
performing a quadratic fit, in the presence of the ensembles A3k1, A3k2, and
A3k3, as displayed in Fig. 4. For both setups, fits confirm the presence of
${\rm O}((am_{\mathrm{q}})^{2})$ effects in this case. The two rightmost
points (A3k1, A3k3) have not been included in these fits. Including them would
result in a very large value of $\chi^{2}/\mathrm{d.o.f}$. This may also be
related to the fact that no clear-cut plateaux are seen in the current quark
mass data for these ensembles. This could be explained by the fact that
(boundary) cut-off effects for these comparatively large masses (in lattice
units) are substantial. Estimates of $Zr_{\mathrm{m}}$ and
$Z(r_{\mathrm{m}}-1)$, obtained as the linear coefficient of the quadratic
fits around $\kappa_{\mathrm{crit}}$, are compatible with those from linear
fits. The influence of the quadratic term on our final result is therefore
negligible. This ensures that our results are not affected by ${\rm
O}((am_{\rm q})^{2})$ systematic errors at $\beta=3.3$, which is our coarsest
lattice. The same conclusion holds for the finer lattices, since also for them
$am_{\rm q}$ is small and linear fits have small $\chi^{2}/\mathrm{d.o.f.}$

As implied by eq. (2.11), $\kappa_{\mathrm{crit}}$ and $Zr_{\mathrm{m}}$ are
assessed as the intercept and the slope of the linear fit to the unitary data.
Similarly, eq. (2.12) tells us that $Z(r_{\mathrm{m}}-1)$ can be estimated
from the slope of the linear fit to the non-unitary data. Our final findings
for $Zr_{\mathrm{m}}$, $Z(r_{\mathrm{m}}-1)$, and $\kappa_{\mathrm{crit}}$ are
listed in Table 3.

$\beta$ | $Zr_{\mathrm{m}}$ | $Z(r_{\mathrm{m}}-1)$ | $\kappa_{\mathrm{crit}}$
---|---|---|---
3.3 | 4.240(134) | 3.621(133) | 0.1364904(18)
3.414 | 2.015(24) | 1.092(29) | 0.1369478(26)
3.512 | 1.561(21) | 0.603(26) | 0.1371320(26)
3.676 | 1.383(19) | 0.329(20) | 0.1371611(25)
3.81 | 1.263(47) | 0.173(40) | 0.1370310(9)
Table 3: Results from the PCAC mass analyses. The second and fourth column
show results obtained in a unitary setup; the third column refers to the non-
unitary setup.

###  5.2 Renormalisation constant $Z$

As the next step in our analysis, we extract the renormalisation constant
$Z\equiv(Z_{\rm P}/Z_{\rm S}Z_{\rm A})$ from the ratio (3.5), using the subset
of gauge field ensembles listed in Table 1 which are not emphasised in italic
font.888As explained in Section 4, the ensembles in italics have been
generated for the purpose of performing reliable fits of the data in Fig. 3
and 4, in order to accurately measure their slopes. These extra ensembles have
not been used for the computation of $Z$, as they do not increase the accuracy
of the result. D1k1 (marked by an asterisk) is also not taken into account.
The correlation functions in eq. (3.5) are computed for two choices of $t_{1}$
and $t_{2}$. Our first choice is $t_{1}\approx T/3$ and $t_{2}\approx 2T/3$,
and the results obtained in this fashion are denoted as ${Z^{\\{T/3\\}}}$.
Alternatively, choosing $t_{1}\approx T/4$ and $t_{2}\approx 3T/4$ yields a
second $Z$ estimate denoted as ${Z^{\\{T/4\\}}}$. When $T/3$ and $T/4$ are not
integers, $t_{1}$ and $t_{2}$ are rounded up/down to the nearest integer.

Figure 5:  Left: Ward identity estimates of $Z$, plotted against time
$y_{0}/T$, for one representative ensemble for each lattice spacing (except
for $\beta=3.3$, corresponding to the coarsest lattice). The dashed vertical
lines bracket the two central time slices that determine the final value of
$Z$. Right: Chiral extrapolation of $Z$ at fixed $\beta$ obtained from the
Ward identity with the mass term (squares) and without it (diamonds). In the
massless case, a possible linear range in $am$ is illustrated by the dashed
line joining the two leftmost points. In the massive case, no significant
quark mass dependence is observed; the dashed line through the squares is a
linear fit where the slope vanishes within its uncertainty. Note that the
errors of the PCAC masses are also displayed and taken into account in the
fits via orthogonal distance regression [48].

In the left part of Fig. 5, we depict ${Z^{\\{T/3\\}}}$ as a function of
$y_{0}/T$ for several representative ensembles (we remind the reader that $T$
is approximately constant in physical units). Contrary to the PCAC masses in
Fig. 2, these local estimators of $Z$ do not exhibit plateau-like behaviour;
this was also observed for a similar Ward identity adopted to compute the
improvement coefficient of the vector current in Ref. [25]. Note, however,
that this is not problematic; since $Z$ is obtained from a Ward identity, its
value at any time slice qualifies as a well-defined estimate. We prefer to err
on the side of caution and quote the average of the two central time slices as
our best $Z$ estimate. Results for the two determinations of $Z$ are collected
in Table 2, where we see that ${Z^{\\{T/3\\}}}$ and ${Z^{\\{T/4\\}}}$ are not
compatible, indicating the presence of lattice artefacts that also differ
noticeably. We consider ${Z^{\\{T/3\\}}}$ the more reliable estimate because
the operator insertions in this case, being further from $x_{0}=0$ and
$x_{0}=T$, are expected to lead to less contamination through cut-off effects
induced by the boundaries.

Since the Ward identity (3.5) is only valid up to lattice artefacts of
$\mathrm{O}(am,a^{2})$, we have to interpolate our data to the chiral point,
in order to eliminate the $\mathrm{O}(am)$-effects and be left with
$\mathrm{O}(a^{2})$ only. As an additional cross-check we also compute $Z$
without the “mass term” $2am\tilde{f}_{\rm PS}(t_{2},t_{1})$ in the Ward
identity (3.5), where $am$ is the PCAC mass from the unitary setup discussed
in the previous section. This chiral interpolation is demonstrated for
$\beta=3.676$ in the right part of Fig. 5. While the data including the “mass
term” shows a very flat behaviour with respect to the current quark mass
(where the associated fit parameter even vanishes within its uncertainty
except for the coarsest lattice spacing), the truncated Ward identity results
in a considerably larger slope. If we exclude the rightmost data point for the
identity without the “mass term”, linear fits to both datasets still agree in
the chiral limit. This situation resembles closely what was observed in Ref.
[15], where $Z$ was measured employing a different Ward identity. We note that
the linear fit is based on the orthogonal distance regression method [48],
taking into account both the error of dependent and independent variables. The
final results for ${Z^{\\{T/3\\}}}$ and ${Z^{\\{T/4\\}}}$ at the chiral point
are also listed in Table 2. Compared to the indirect Ward identity
determination of Ref. [15], they have considerably smaller errors. This
confirms the expectation that the simpler structure of the correlation
functions building the Ward identity (3.4) is preferable from a numerical
perspective; see the discussion at the end of Section 3. On the other hand,
compared to the so-called ’LCP-0’ determination of Ref. [14], our results are
of similar accuracy across the bare couplings investigated. We will use our
results (Table 2) for a precise estimation of $r_{\mathrm{m}}$ in the
following. More details on the relative cut-off effects between the present
determination of $Z$ and the results obtained in Refs. [14, 15, 13] can be
found in Appendix C.

###  5.3 Results for $r_{\mathrm{m}}$

In the final step of our analysis we combine the values of $Zr_{\mathrm{m}}$
obtained in a unitary setup, $Z(r_{\mathrm{m}}-1)$ in a non-unitary setup, and
$Z$ from a chiral Ward identity, in order to arrive at different estimates for
$r_{\mathrm{m}}$. Combining the first two, we construct
${r_{\mathrm{m}}^{\\{\mathrm{u,nu}\\}}}$, defined as

$\displaystyle{r_{\mathrm{m}}^{\\{\mathrm{u,nu}\\}}}$
$\displaystyle=\bigg{(}1-\left[\frac{Z(r_{\text{m}}-1)}{Zr_{\text{m}}}\right]\bigg{)}^{-1}\,,$
(5.2) where the superscripts “u” and “nu” stand for “unitary” and “non-
unitary”, respectively. Combining $Zr_{\mathrm{m}}$ and $Z$, results in
${r_{\mathrm{m}}^{\\{\mathrm{u;}Z\\}}}$, defined as
$\displaystyle{r_{\mathrm{m}}^{\\{\mathrm{u;}Z\\}}}$
$\displaystyle=\frac{Zr_{\text{m}}}{Z}\,.$ (5.3) As mentioned above, this
comes in two versions, ${r_{\mathrm{m}}^{\\{\mathrm{u;}Z,T/3\\}}}$ and
${r_{\mathrm{m}}^{\\{\mathrm{u;}Z,T/4\\}}}$. Moreover, from the second and
third result we gain ${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z\\}}}$ given by
$\displaystyle{r_{\mathrm{m}}^{\\{\mathrm{nu;}Z\\}}}$
$\displaystyle=\frac{Z(r_{\text{m}}-1)}{Z}+1\,,$ (5.4)

which is again worked out for two cases,
${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/3\\}}}$ and
${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/4\\}}}$. All our results for
$r_{\mathrm{m}}$ from these different determinations just outlined are
gathered in Table 4.

$\beta$ | ${r_{\mathrm{m}}^{\\{\mathrm{u,nu}\\}}}$ | ${r_{\mathrm{m}}^{\\{\mathrm{u;}Z,T/3\\}}}$ | ${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/3\\}}}$ | ${r_{\mathrm{m}}^{\\{\mathrm{u;}Z,T/4\\}}}$ | ${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/4\\}}}$
---|---|---|---|---|---
3.3 | 6.848(569) | 5.181(172) | 5.424(169) | 5.588(207) | 5.772(199)
3.414 | 2.183(44) | 2.242(26) | 2.215(32) | 2.348(32) | 2.272(35)
3.512 | 1.629(32) | 1.571(21) | 1.607(26) | 1.617(22) | 1.625(27)
3.676 | 1.312(20) | 1.303(15) | 1.309(19) | 1.312(16) | 1.312(19)
3.81 | 1.158(37) | 1.158(44) | 1.158(37) | 1.164(44) | 1.159(37)
Table 4: Results for $r_{\mathrm{m}}$, obtained via eqs. (5.2) to (5.4).

In principle, the different estimates can differ by $\mathrm{O}(a^{2})$
ambiguities. In Fig. 6 (left) the three determinations
${r_{\mathrm{m}}^{\\{\mathrm{u,nu}\\}}}$,
${r_{\mathrm{m}}^{\\{\mathrm{u;}Z,T/3\\}}}$, and
${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/3\\}}}$ are plotted against the bare
coupling squared; to be able to distinguish between the different estimates,
the data points corresponding to the coarsest lattice spacing ($\beta=3.3$)
are omitted as they exhibit large cut-off effects and are thus well out of the
range displayed here. Results are compatible within their respective
$1\sigma$-errors. In Fig. 6 (right) we take a closer look at this behaviour by
plotting ratios of different $r_{\mathrm{m}}$ estimates as functions of the
lattice spacing squared; the corresponding lattice spacings can be found in
Table 1. Since the ratios have been computed on a line of constant physics,
and assuming that we are in a scaling region where Symazik’s effective theory
of cut-off effects applies, they are expected to be polynomials in the lattice
spacing, tending to 1 in the continuum limit. In this context we introduce an
additional determination,
$r_{\mathrm{m}}^{\\{\mathrm{u},\mathrm{nu};\mathrm{impr}\\}}$, which only
differs from ${r_{\mathrm{m}}^{\\{\mathrm{u,nu}\\}}}$ by an improved version
of the derivative ${\widetilde{\partial}_{0}}$ in eq. (5.1).999The improved
derivative is defined as
$a\partial_{\mu}f(x)\equiv\frac{1}{12}[-f(x+2a\hat{\mu})+8f(x+a\hat{\mu})-8f(x-a\hat{\mu})+f(x-2a\hat{\mu})]$
and its corresponding second derivative by
$a^{2}\partial_{\mu}^{*}\partial_{\mu}f(x)\equiv\frac{1}{12}[-f(x+2a\hat{\mu})+16f(x+a\hat{\mu})-30f(x)+16f(x-a\hat{\mu})-f(x-2a\hat{\mu})]$
as shown by eq. (B.4) in Ref. [14].  These ratios are very close to one except
for one of the data points at $\beta=3.3$, for which the ratio is
significantly larger. Even though it would be sufficient to demonstrate that
these ratios of $r_{\mathrm{m}}$ approach unity with a rate $\propto a^{2}$ or
higher in our particular line of constant physics framework, such ambiguities
appear to be nearly absent for $a<0.1\mathrm{fm}$.

We tried to model the data sets with and without the $\beta=3.3$ points, using
polynomials in the lattice spacing, constrained to one in the continuum limit.
When a linear term is included, we obtain unsatisfactory fits with
$\chi^{2}/\mathrm{d.o.f.}>3$. We thus conclude that our results are compatible
with the theoretical expectation of $\mathrm{O}(a^{2})$ lattice artefacts or
higher (see also Appendix C).

Figure 6: Left: Results for different $r_{\mathrm{m}}$ estimates as reported
in Table 4. The results for $\beta=3.3$ are not shown. Right: Ratio of
different $r_{\mathrm{m}}$ determinations as a function of the squared lattice
spacing. The dashed horizontal line indicates the expected continuum result.

As our preferred determination of $r_{\mathrm{m}}$ we advocate
${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/3\\}}}$ because of its small statistical
errors in our range of bare couplings and the poorer scaling behaviour of the
other estimators at the coarsest lattice spacing. In Fig. 7 we show this
result including the two-loop perturbative prediction of Ref. [18]. An
important observation is that the non-perturbative estimates strongly deviate
from the perturbative prediction in this region of strong couplings. A similar
behaviour was also observed in several studies of renormalisation factors for
which one-loop perturbative predictions are available (see, e.g. [49, 25]).
Here, we confirm this finding also for two-loop perturbation theory. We also
compare our results with those of other works. In Ref. [13], $r_{\mathrm{m}}$
was determined for two values of the bare coupling, from an alternative
renormalization condition. As inferred by Fig. 7 this result agrees with ours
at the smaller coupling, while it deviates notably at the larger coupling,
most likely due to $\mathrm{O}(a^{2})$ ambiguities (or higher).

Figure 7: Non-perturbative determination of
${r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/3\\}}}$ (open circles), compared to the
results of Ref. [13] (filled diamonds) and those of two-loop perturbation
theory [18] (horizontal dotted line). The dashed line is the interpolation
(5.5) and the vertical dotted lines correspond to the bare couplings used in
CLS simulations.

Our final result consists of a continuous interpolation formula for
$r_{\mathrm{m}}={r_{\mathrm{m}}^{\\{\mathrm{nu;}Z,T/3\\}}}$. Our data is best
described by a Padé ansatz, constrained to the two-loop prediction of Ref.
[18] for small couplings, of the form

$r_{\text{m}}(g_{0}^{2})=1.0+0.004630\,g_{0}^{4}\times\left\\{\frac{1+c_{1}\,g_{0}^{2}+c_{2}\,g_{0}^{4}}{1+c_{3}\,g_{0}^{2}}\right\\}\,,$
(5.5a) where $\quad c_{i}=\left(-7.86078,5.49175,-0.54078\right)\,,$ (5.5b)
and
$\mathrm{cov}(c_{i},c_{j})=\left(\begin{array}[]{lll}\phantom{-}3.699760\phantom{\times
10^{-1}}&-2.198586\phantom{\times 10^{-1}}&-1.476913\times 10^{-3}\\\
-2.198586\phantom{\times 10^{-1}}&\phantom{-}1.306512\phantom{\times
10^{-1}}&\phantom{-}8.776569\times 10^{-4}\\\ -1.476913\times
10^{-3}&\phantom{-}8.776569\times 10^{-4}&\end{array}\right)\,,$ (5.5c)

which is also displayed in Fig. 7. The fit function describes our data with
$\chi^{2}/{\rm d.o.f.}=1.37$ and provides errors of a size comparable to the
fitted data points.

The interpolation formula can now be used in order to determine
$r_{\mathrm{m}}$ at the couplings used in CLS simulations for the computation
of hadronic quantities [21, 13, 23]. Since the CLS coupling $\beta=3.85$ lies
outside the range of our $r_{\mathrm{m}}$ computations, we perform a short
extrapolation in order to provide a value for $r_{\mathrm{m}}(\beta=3.85)$. A
systematic error, estimated as the difference between the lower error bar of
our data point at $\beta=3.81$ and the extrapolated value at $\beta=3.85$
($\sigma_{\mathrm{syst}}=0.027$) is added to the statistical error
($\sigma_{\mathrm{stat}}=0.018$) in quadrature. Our final $r_{\mathrm{m}}$
results at the CLS couplings are collected in Table 5.

$\beta$ | 3.4 | 3.46 | 3.55 | 3.7 | 3.85
---|---|---|---|---|---
$r_{\mathrm{m}}$ | 2.335(31) | 1.869(19) | 1.523(14) | 1.267(16) | 1.149(18)(27)[33]
Table 5: Values for $r_{\mathrm{m}}$ at the couplings used in CLS simulations,
obtained from the interpolation formula (5.5). As mentioned in the text, an
additional systematic error was added to the $\beta=3.85$ result. The errors
are displayed in this way:
$(\sigma_{\mathrm{stat}})(\sigma_{\mathrm{syst}})[\sigma_{\mathrm{total}}]$.

##  6 Summary

With the non-perturbative computation of the ratio of the renormalisation
constants of non-singlet and singlet scalar densities, $r_{\rm m}\equiv Z_{\rm
S}/Z_{\rm S}^{0}$, presented in this paper we have addressed a quantity, which
not only enters the renormalisation pattern of quark masses in lattice QCD
with Wilson fermions, but also constitutes an important ingredient in
calculations of renormalised nucleon (and other baryon) matrix elements of
singlet scalar densities, known as sigma terms.

Our strategy to calculate $r_{\rm m}$ merges the functional dependences of the
PCAC quark mass in terms of the subtracted quark mass, evaluated in a unitary
as well as a non-unitary setting with respect to the choice of sea and valence
quark masses. In the vicinity of the chiral limit, these dependences are found
to be linear, so that $r_{\rm m}$ can be obtained through the associated quark
mass slopes with confidence and superior control of statistical and systematic
errors. The finite-volume numerical simulations of ${\rm O}(a)$ improved QCD
with Schrödinger functional boundary conditions that enter the analysis
realise a line of constant physics by working in a volume of spatial extent
$L\approx 1.2\,$fm and thereby fixing all other relevant length scales in
physical units. This guarantees that $r_{\rm m}$ becomes a smooth function of
the bare gauge coupling as the lattice spacing is varied, where any
potentially remaining intrinsic ambiguities disappear monotonically towards
the continuum limit at a rate that stays beyond the sensitivity of the ${\rm
O}(a)$ improved theory.

Our central results, which hold for a lattice discretisation of QCD with three
flavours of non-perturbatively ${\rm O}(a)$ improved Wilson-clover sea quarks
and tree-level Symanzik-improved gluons, are the continuous parameterisation
of $r_{\rm m}$ as a function of the squared bare gauge coupling
$g_{0}^{2}=6/\beta$ in eq. (5.5), as well as its values in Table 5 at the
specific strong-coupling $\beta$ values of large-volume CLS simulations [21,
22, 13, 23].

Along with the numerical implementation of our strategy to extract $r_{\rm
m}$, we have also developed a new method to determine the scale independent
combination $Z=Z_{\rm P}/(Z_{\rm S}Z_{\rm A})$ of renormalisation parameters
of quark bilinears in the pseudoscalar, (non-singlet) scalar and axial vector
channel, respectively. It relies upon a Ward identity that, according to our
knowledge, has not yet appeared explicitly in the literature. Since, as
explained in Sections 2 and 3, the renormalisation factor $Z$ is actually
required to isolate $r_{\rm m}$ from the unitary and non-unitary quark mass
slopes, we have employed the estimates on $Z$ from this approach in our final
results of $r_{\rm m}$. However, this was primarily done for practical reasons
and served the purpose of demonstrating the feasibility of the Ward identity
method for $Z$. In fact, it is apparent from the discussion in Appendix C and
Figure 8 that these new values for $Z$ are fully compatible with the earlier
determinations available from Refs. [14, 15] and are neither superior in
statistical precision nor in systematics regarding lattice artefacts.
Nevertheless we give an interpolation formula for the present $Z$ (Table 6)
for completeness.

Finally we recall the subtlety discussed in Section 2: away from the chiral
limit, the dependence of (re)normalisation parameters should be
$Z(\tilde{g}_{0}^{2}),r_{\rm m}(\tilde{g}_{0}^{2})$, with $\tilde{g}_{0}^{2}$
defined in eq. (2.13). In order to be able to combine our results with CLS
low-energy quantities such as those of Refs. [22, 29], we should use the
expansion

$Z(\tilde{g}_{0}^{2})=Z(g_{0}^{2})\Big{[}1+\dfrac{\partial\ln
Z(g_{0}^{2})}{\partial g_{0}^{2}}\dfrac{1}{N_{\rm
f}}b_{g}(g_{0}^{2})g_{0}^{2}a{\rm Tr}M_{\rm q}\Big{]}\,,$ (6.1)

see also Ref. [50], and similarly for $r_{\rm m}(\tilde{g}_{0}^{2})$. At
present, $b_{g}$ is only known in perturbation theory [51]; $b_{g}=0.012N_{\rm
f}g_{0}^{2}$. The correction $\partial\ln Z/\partial g_{0}^{2}$ as well as
$\partial\ln r_{\rm m}/\partial g_{0}^{2}$ , computed at CLS values of the
inverse coupling $\beta$, can be found in Table 7.

Acknowledgements. This work is supported by the Deutsche
Forschungsgemeinschaft (DFG) through the Research Training Group “GRK 2149:
Strong and Weak Interactions – from Hadrons to Dark Matter” (J. H., F. J. and
P. L. J. P.). We acknowledge the computer resources provided by the WWU IT,
formerly ‘Zentrum für Informationsverarbeitung (ZIV)’, of the University of
Münster (PALMA-II HPC cluster) and thank its staff for support.

##  Appendix A Schrödinger functional correlation functions

The Schrödinger functional correlation functions employed in this work are
defined as

$\displaystyle f_{\mathrm{P}}^{ij}$
$\displaystyle=-\frac{1}{2}\frac{a^{9}}{L^{3}}\sum_{\bf
x,u,v}\left\langle\bar{\psi}_{i}(x)\gamma_{5}\psi_{j}(x)\cdot\bar{\zeta}_{j}({\bf
v})\gamma_{5}\zeta_{i}(\bf u)\right\rangle\,,$ (A.1) $\displaystyle
g_{\mathrm{P}}^{ij}$ $\displaystyle=-\frac{1}{2}\frac{a^{9}}{L^{3}}\sum_{\bf
x,u,v}\left\langle\bar{\psi}_{i}(x)\gamma_{5}\psi_{j}(x)\cdot\bar{\zeta}^{\prime}_{j}({\bf
u})\gamma_{5}\zeta^{\prime}_{i}(\bf v)\right\rangle\,,$ (A.2) $\displaystyle
f_{\mathrm{A}}^{ij}$ $\displaystyle=-\frac{1}{2}\frac{a^{9}}{L^{3}}\sum_{\bf
x,u,v}\left\langle\bar{\psi}_{i}(x)\gamma_{0}\gamma_{5}\psi_{j}(x)\cdot\bar{\zeta}_{j}({\bf
u})\gamma_{5}\zeta_{i}(\bf v)\right\rangle\,,$ (A.3) $\displaystyle
g_{\mathrm{A}}^{ij}$ $\displaystyle=-\frac{1}{2}\frac{a^{9}}{L^{3}}\sum_{\bf
x,u,v}\left\langle\bar{\psi}_{i}(x)\gamma_{0}\gamma_{5}\psi_{j}(x)\cdot\bar{\zeta}^{\prime}_{j}({\bf
u})\gamma_{5}\zeta^{\prime}_{i}(\bf v)\right\rangle\,,$ (A.4) $\displaystyle
F_{1}^{ij}$ $\displaystyle=-\frac{1}{2}\frac{a^{12}}{L^{6}}\sum_{\bf
u^{\prime},v^{\prime},u,v}\left\langle\bar{\zeta}^{\prime}_{i}({\bf
u^{\prime}})\gamma_{5}\zeta^{\prime}_{j}({\bf
v^{\prime}})\cdot\bar{\zeta}_{j}({\bf u})\gamma_{5}\zeta_{i}(\bf
v)\right\rangle\,.$ (A.5)

They refer to the general case of two distinct, i.e. not necessarily mass-
degenerate quark flavours $i,j$. Summation over the indices $i$ and $j$ is not
implied. The space-time point $x$ lies in the lattice bulk; i.e. $0<x_{0}<T$.
The Dirichlet boundary fields $\bar{\zeta}_{j}(\bf u)$ and $\zeta_{i}(\bf v)$
live on time slice $x_{0}=0$, while $\bar{\zeta}^{\prime}_{j}(\bf u^{\prime})$
and $\zeta^{\prime}_{i}(\bf v^{\prime})$ live on time slice $x_{0}=T$; the
boundary fields are introduced in Ref.[36].

##  Appendix B Wick contractions of correlation functions

In this appendix we briefly explain how to obtain eq. (3.5) from eq. (3.4).
The idea is to perform the Wick contractions of the correlation functions,
arriving at expressions which are traces of flavour matrices, multiplying
traces of products of quark propagators and $\gamma$-matrices. This procedure
has been described in full detail in Ref. [15], which deals with more
complicated Ward identities; we refer the reader to that work for unexplained
notation. Here we will only present the main features of the proof.

We start with the r.h.s. of eq. (3.4). The Wick contractions result in

$\displaystyle-d^{abe}\,a^{3}\sum_{\bf y}\langle P^{e}(y){\cal O}^{c}\rangle=$
$\displaystyle-d^{abe}{\rm Tr}[T^{e}T^{c}]\dfrac{a^{9}}{L^{3}}\sum_{\bf
y}\sum_{\bf u,v}\Bigg{\langle}\,\hbox{tr}\,\bigg{\\{}[\psi(y)\bar{\zeta}({\bf
u})]_{\rm F}\gamma_{5}[\zeta({\bf v})\bar{\psi}(y)]_{\rm
F}\gamma_{5}\bigg{\\}}\Bigg{\rangle}$ $\displaystyle=$
$\displaystyle\,d^{abc}f_{\rm P}(y_{0})\,,$ (B.1)

where the second equality implicitly defines $f_{\rm P}$ (see also eq. (A.1)
and Appendix B of Ref. [14]). The left-hand-side consists of correlation
functions with one boundary operator and two insertions in the bulk. So the
Wick contractions of such a correlation function give:

$\displaystyle a^{6}\sum_{\bf x,y}\langle A_{0}^{a}(x)S^{b}(y){\cal
O}^{c}\rangle=$ $\displaystyle\dfrac{\mathrm{i}a^{12}}{L^{3}}{\rm
Tr}[T^{a}T^{b}T^{c}]\sum_{\bf x,y}\sum_{\bf
u,v}\Bigg{\langle}\,\hbox{tr}\,\bigg{\\{}\gamma_{0}\gamma_{5}[\psi(x)\bar{\psi}(y)]_{\rm
F}[\psi(y)\bar{\zeta}({\bf u})]_{\rm F}\gamma_{5}[\zeta({\bf
v})\bar{\psi}(x)]_{\rm F}\bigg{\\}}\Bigg{\rangle}$
$\displaystyle+\dfrac{\mathrm{i}a^{12}}{L^{3}}{\rm
Tr}[T^{c}T^{b}T^{a}]\sum_{\bf x,y}\sum_{\bf
u,v}\Bigg{\langle}\,\hbox{tr}\,\bigg{\\{}\gamma_{0}\gamma_{5}[\psi(x)\bar{\zeta}({\bf
u})]_{\rm F}\gamma_{5}[\zeta({\bf v})\bar{\psi}(y)]_{\rm
F}[\psi(y)\bar{\psi}(x)]_{\rm F}\bigg{\\}}\Bigg{\rangle}$ $\displaystyle=$
$\displaystyle\,\mathrm{i}{\rm Tr}[T^{a}T^{b}T^{c}]f_{\rm
AS;1}(x_{0},y_{0})+\mathrm{i}{\rm Tr}[T^{c}T^{b}T^{a}]f_{\rm
AS;2}(x_{0},y_{0})$ $\displaystyle=$
$\displaystyle\,\dfrac{1}{2}\bigg{[}-d^{abc}{\rm Re}\,f_{\rm
AS;1}(x_{0},y_{0})+f^{abc}{\rm Im}\,f_{\rm AS;1}(x_{0},y_{0})\bigg{]}\,.$
(B.2)

The second in the above string of equations implicitly defines the two traces
of quark propagators (devoid of flavour structure) as $f_{\rm
AS;1}(x_{0},y_{0})$ and $f_{\rm AS;2}(x_{0},y_{0})$. In the last equation we
have made use of the fact that the two traces of propagators are complex
conjugates of each other which is a consequence of the
$\gamma_{5}$-Hermiticity property of Wilson fermion propagators. Finally, the
fact that the above correlation function is invariant under charge conjugation
leads to the vanishing of the term proportional to $f^{abc}$ in the last
expression. Hence, we obtain

$\displaystyle a^{6}\sum_{\bf x,y}\langle A_{0}^{a}(x)S^{b}(y){\cal
O}^{c}\rangle=-\dfrac{1}{2}d^{abc}{\rm Re}\,f_{\rm
AS;1}(x_{0},y_{0})=-d^{abc}f_{\rm AS}(x_{0},y_{0})\,,$ (B.3)

which implicitly defines $f_{\rm AS}(x_{0},y_{0})$. The correlation functions
$f_{\rm P}$ and $f_{\rm AS}$ are schematically drawn in Fig. 1. Analogously,
from the mass dependent term of the Ward identity we also define
$\tilde{f}_{\rm PS}(x_{0},y_{0})$; the summation over all times from $t_{1}$
up to $t_{2}$ (see eq. (3.4)) is included in its definition. It is important
to note that $d^{abc}$ appears in both eqs. (B.1) and (B.3). Therefore, it
cancels out in the Ward identity, which becomes an expression between traces
of propagators, without any flavour indices. Putting everything together, we
eventually obtain eq. (3.5).

##  Appendix C Comparison of $Z$ determinations and scaling tests

In this appendix we present more details on our $Z$ results, listed in Table
2. In Fig. 8 and Table 6 our preferred determination for $Z$, namely
${Z^{\\{T/3\\}}}$ is compared to Ref. [14] (de Divitiis et al.), $Z$
determined at two values of the gauge coupling in Ref. [32] (Bali et al.) and
to a $Z$ estimate that we work out from the results of Refs. [49] and [15]
(Heitger et al.). In particular, we extract the axial current normalisation
$Z_{\mathrm{A}}$ at our couplings from the interpolation formula of Ref. [49]
and combine it with the ratio $Z_{\mathrm{S}}/Z_{\mathrm{P}}$ of the
pseudoscalar and scalar renormalisation constants from Ref. [15]. In addition,
we give an interpolation formula for our preferred determination for $Z$ (also
displayed in Fig. 8).

Our result agrees with the other determinations at weaker bare couplings,
while disagreements are seen at stronger couplings. These are attributed to
lattice artefacts associated with intrinsic ambiguities of ${\rm O}(a^{2})$ or
higher between different determinations. Agreement is generally better between
our results and those of Ref. [14] (de Divitiis et al.).

Figure 8: $Z$ results, obtained with different methods, as a function of the squared bare coupling $g_{0}^{2}$. The preferred determination of this work is $Z={Z^{\\{T/3\\}}}$ (pentagons). The squares are obtained by combining results from Refs. [49] and [24] (Heitger et al.). The two $Z$ estimates determined in Ref. [32] (Bali et al.) are depicted by triangles. The circles correspond to the $Z$ results from Ref. [14] (de Divitiis et al.). One-loop perturbation theory is illustrated by the dotted line, Ref. [14]. The dashed line shows the interpolation (C.1a) of ${Z^{\\{T/3\\}}}$ (excluding the coarsest lattice spacing from the fit). The vertical dotted lines correspond to the bare couplings used in CLS simulations.  $\beta$ |  $Z={Z^{\\{T/3\\}}}$ this work |  $Z$, LCP-0 de Divitiis et al. |  $Z$, LCP-1 de Divitiis et al. |  $1/Z_{\mathrm{A}}\cdot Z_{\mathrm{P}}/Z_{\mathrm{S}}$ Heitger et al.
---|---|---|---|---
3.3 | 0.8184(77) | 0.7462(56) | 0.7896(36) | 0.884(26)0
3.414 | 0.8987(43) | 0.8762(40) | 0.8992(26) | 0.990(12)0
3.512 | 0.9935(38) | 0.9764(33) | 0.9861(23) | 1.0396(80)
3.676 | 1.0621(36) | 1.0588(31) | 1.0611(23) | 1.0901(89)
3.81 | 1.0907(13) | 1.0882(11) | 1.0884(8)0 | 1.1029(61)
Table 6: Comparison of our preferred $Z$ determination with results from Ref.
[14] (de Divitiis et al.) and the combination of results from Refs. [49] and
[24] (Heitger et al.).

In order to confirm this claim of consistency (leaving aside higher cut-off
effects) we construct ratios of different determinations and investigate their
behaviour as a function of the lattice spacing. Interestingly, rather than
${\rm O}(a^{2})$, leading cut-off effects of ${\rm O}(a^{3})$ can be
identified in the ratio ${Z^{\\{T/3\\}}}/{Z^{\\{T/4\\}}}$, as seen in Fig. 9
(left). The scaling behaviour of our results compared to those of previous
works is shown in Fig. 9 (right). All ratios are fitted with an ansatz
$1+ca^{3}$, excluding the coarsest lattice spacing. When adding a term linear
in the lattice spacing, its fit parameter vanishes within its uncertainty in
all cases. In conclusion, these scaling tests indicate that our results for
$Z$ are in accordance with the theoretical expectation of ${\rm O}(a^{2})$
ambiguities or higher which by virtue of the imposed line of constant physics
decrease monotically towards the continuum limit.

Figure 9: Left: Ratio of our results ${Z^{\\{T/3\\}}}/{Z^{\\{T/4\\}}}$, fitted
as a function of the lattice spacing. Right: Ratios of ${Z^{\\{T/3\\}}}$ and
$Z$ computed in previous works, fitted as a function of the lattice spacing.
The coarsest lattice spacing is excluded from the fits. The dashed lines are
the fits while the horizontal dotted lines indicate the expected continuum
results.

In addition, we interpolate our $Z$ data using a Padé ansatz, constrained to
the one-loop prediction of Ref. [52] for small couplings; see eq. (C.1) and
Fig. 8. Owing to relevant higher order cut-off effects, shown in the right
panel of Fig. 9, we do not include the coarsest lattice spacing in the fit.
Other than that, it must be noted that the coarsest lattice spacing is well
outside the range of CLS couplings. We obtain

$Z(g_{0}^{2})=1+0.0703169\cdot
g_{0}^{2}\times\frac{1+d_{1}g_{0}^{4}}{1+d_{2}g_{0}^{2}}\,,$ (C.1a) where
$\quad d_{i}=\left(-0.34504,-0.52309\right)\,,$ (C.1b) and
$\mathrm{cov}(d_{i},d_{j})=\left(\begin{array}[]{lll}\phantom{-}2.798505\times
10^{-7}&\phantom{-}6.577477\times 10^{-7}\\\ \phantom{-}6.577477\times
10^{-7}&\end{array}\right).$ (C.1c)

Our $Z$ results at the CLS couplings are gathered in Table 7 and compared to
those of Ref. [14] for two different LCP conditions. The two outmost CLS
$\beta$ values ($3.4$ and $3.85$) lie outside the range of our fitted $Z$
estimates, so they are obtained by extrapolation. Their systematic errors are
estimated from the statistical uncertainty of the nearest ${Z^{\\{T/3\\}}}$
data point: the systematic error of ${Z^{\\{T/3\\}}}(\beta=3.4)$ is the
statistical error of ${Z^{\\{T/3\\}}}(\beta=3.414)$ and that of
${Z^{\\{T/3\\}}}(\beta=3.85)$ is the statistical error of
${Z^{\\{T/3\\}}}(\beta=3.81)$. These systematic errors are added to the
statistical ones in quadrature.

$\beta$ |  $Z={Z^{\\{T/3\\}}}$ interpolated, this work |  $Z$, LCP-0 de Divitiis et al. |  $Z$, LCP-1 de Divitiis et al. | $\partial\ln Z/\partial g_{0}^{2}$ | $\partial\ln r_{\mathrm{m}}/\partial g_{0}^{2}$
---|---|---|---|---|---
3.4 | $\phantom{-}0.8798(47)(43)[64]$ | $\phantom{-}0.8758(52)$ | $\phantom{-}0.8981(35)$ | $-3.241(144)$ | $\phantom{-}8.975(195)$
3.46 | $\phantom{-}0.9507(25)$ | $\phantom{-}0.9320(50)$ | $\phantom{-}0.9468(35)$ | $-1.974(62)$ | $\phantom{-}5.915(179)$
3.55 | $\phantom{-}1.0147(15)$ | $\phantom{-}0.9937(42)$ | $\phantom{-}1.0015(30)$ | $-1.104(23)$ | $\phantom{-}3.647(149)$
3.7 | $\phantom{-}1.0696(13)$ | $\phantom{-}1.0591(23)$ | $\phantom{-}1.0612(17)$ | $-0.522(7)$ | $\phantom{-}1.962(90)$
3.85 | $\phantom{-}1.0961(12)(13)[18]$ | $\phantom{-}1.0975(25)$ | $\phantom{-}1.0971(18)$ | $-0.278(3)$ | $\phantom{-}1.190(49)$
Table 7:  $Z$ values at the couplings used in CLS simulations, obtained from
the interpolation formula (C.1a), excluding the coarsest lattice spacing from
the fit. As explained in the text, an additional systematic error is added to
the $\beta=3.85$ and $\beta=3.4$ results and the errors are displayed as:
$(\sigma_{\mathrm{stat}})(\sigma_{\mathrm{syst}})[\sigma_{\mathrm{total}}]$.
In the last two columns we list $\partial\ln Z/\partial g_{0}^{2}$ for
${Z^{\\{T/3\\}}}$, obtained by differentiating eq. (C.1a) as well as
$\partial\ln r_{\rm m}/\partial g_{0}^{2}$ via differentiating eq. (5.5).

The interested reader may also use our interpolation formula for $Z$ (from eq.
(C.1)) and $r_{\rm m}$ (from eq. (5.5)) and the covariance between the fit
parameters of the two different interpolations,

$\mathrm{cov}(d_{i},c_{j})=\left(\begin{array}[]{lll}\phantom{-}3.502738\times
10^{-4}&-2.149363\times 10^{-4}&-1.303033\times 10^{-7}\\\
\phantom{-}1.464606\times 10^{-3}&-8.876895\times 10^{-4}&\end{array}\right),$
(C.2)

to construct combinations of the two such as $Zr_{\mathrm{m}}$.

## References

  * [1] M. Bochicchio, L. Maiani, G. Martinelli, G. C. Rossi and M. Testa, _Chiral Symmetry on the Lattice with Wilson Fermions_ , _Nucl. Phys. B_ 262 (1985) 331.
  * [2] L. Maiani, G. Martinelli, M. L. Paciello and B. Taglienti, _Scalar Densities and Baryon Mass Differences in Lattice QCD With Wilson Fermions_ , _Nucl. Phys. B_ 293 (1987) 420.
  * [3] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, _A General method for nonperturbative renormalization of lattice operators_ , _Nucl. Phys. B_ 445 (1995) 81, [hep-lat/9411010].
  * [4] C. Sturm, Y. Aoki, N. H. Christ, T. Izubuchi, C. T. C. Sachrajda and A. Soni, _Renormalization of quark bilinear operators in a momentum-subtraction scheme with a nonexceptional subtraction point_ , _Phys. Rev. D_ 80 (2009) 014501, [0901.2599].
  * [5] S. Capitani, M. Lüscher, R. Sommer and H. Wittig, _Non-perturbative quark mass renormalization in quenched lattice QCD_ , _Nucl. Phys. B_ 544 (1999) 669, [hep-lat/9810063]. [Erratum: Nucl. Phys. B 582, 762 (2000)].
  * [6] M. Dalla Brida, S. Sint and P. Vilaseca, _The chirally rotated Schrödinger functional: theoretical expectations and perturbative tests_ , _JHEP_ 08 (2016) 102, [1603.00046].
  * [7] G. M. de Divitiis and R. Petronzio, _Nonperturbative renormalization constants on the lattice from flavor nonsinglet Ward identities_ , _Phys. Lett. B_ 419 (1998) 311, [hep-lat/9710071].
  * [8] T. Bhattacharya, S. Chandrasekharan, R. Gupta, W.-J. Lee and S. R. Sharpe, _Nonperturbative renormalization constants using Ward identities_ , _Phys. Lett. B_ 461 (1999) 79, [hep-lat/9904011].
  * [9] T. Bhattacharya, R. Gupta, W.-J. Lee and S. R. Sharpe, _Order a improved renormalization constants_ , _Phys. Rev. D_ 63 (2001) 074505, [hep-lat/0009038].
  * [10] M. Guagnelli, R. Petronzio, J. Rolf, S. Sint, R. Sommer and U. Wolff, _Non-perturbative results for the coefficients $b_{\mathrm{m}}$ and $b_{\mathrm{a}}-b_{\mathrm{P}}$ in ${\rm O}(a)$ improved lattice QCD_, _Nucl. Phys. B_ 595 (2001) 44, [hep-lat/0009021].
  * [11] T. Bhattacharya, R. Gupta, W. Lee and S. R. Sharpe, _Scaling behavior of discretization errors in renormalization and improvement constants_ , _Phys. Rev. D_ 73 (2006) 114507, [hep-lat/0509160].
  * [12] P. Fritzsch, J. Heitger and N. Tantalo, _Non-perturbative improvement of quark mass renormalization in two-flavour lattice QCD_ , _JHEP_ 08 (2010) 074, [1004.3978].
  * [13] G. S. Bali, E. E. Scholz, J. Simeth and W. Söldner, _Lattice simulations with $N_{\mathrm{f}}=2+1$ improved Wilson fermions at a fixed strange quark mass_, _Phys. Rev. D_ 94 (2016) 074501, [1606.09039].
  * [14] G. M. de Divitiis, P. Fritzsch, J. Heitger, C. C. Köster, S. Kuberski and A. Vladikas, _Non-perturbative determination of improvement coefficients $b_{\mathrm{m}}$ and $b_{\mathrm{A}}-b_{\mathrm{P}}$ and normalisation factor $Z_{\mathrm{m}}Z_{\mathrm{P}}/Z_{\mathrm{A}}$ with $N_{\mathrm{f}}=3$ Wilson fermions_, _Eur. Phys. J. C_ 79 (2019) 797, [1906.03445].
  * [15] J. Heitger, F. Joswig and A. Vladikas, _Ward identity determination of $Z_{\mathrm{S}}/Z_{\mathrm{P}}$ for $N_{\mathrm{f}}=3$ lattice QCD in a Schrödinger functional setup_, _Eur. Phys. J. C_ 80 (2020) 765, [2005.01352].
  * [16] G. Bali, S. Bürger, S. Collins, M. Göckeler, M. Gruber, S. Piemonte et al., _Nonperturbative Renormalization in Lattice QCD with three Flavors of Clover Fermions: Using Periodic and Open Boundary Conditions_ , 2012.06284.
  * [17] T. Bhattacharya, R. Gupta, W. Lee, S. R. Sharpe and J. M. Wu, _Improved bilinears in lattice QCD with non-degenerate quarks_ , _Phys. Rev. D_ 73 (2006) 034504, [hep-lat/0511014].
  * [18] M. Constantinou, M. Hadjiantonis, H. Panagopoulos and G. Spanoudes, _Singlet versus nonsinglet perturbative renormalization of fermion bilinears_ , _Phys. Rev. D_ 94 (2016) 114513, [1610.06744].
  * [19] M. Lüscher and P. Weisz, _On-Shell Improved Lattice Gauge Theories_ , _Commun. Math. Phys._ 97 (1985) 59. [Erratum: Commun. Math. Phys. 98, 433 (1985)].
  * [20] B. Sheikholeslami and R. Wohlert, _Improved Continuum Limit Lattice Action for QCD with Wilson Fermions_ , _Nucl. Phys. B_ 259 (1985) 572.
  * [21] M. Bruno et al., _Simulation of QCD with $N_{\mathrm{f}}=2+1$ flavors of non-perturbatively improved Wilson fermions_, _JHEP_ 02 (2015) 043, [1411.3982].
  * [22] M. Bruno, T. Korzec and S. Schaefer, _Setting the scale for the CLS $2+1$ flavor ensembles_, _Phys. Rev. D_ 95 (2017) 074504, [1608.08900].
  * [23] D. Mohler, S. Schaefer and J. Simeth, _CLS $2+1$ flavor simulations at physical light- and strange-quark masses_, _EPJ Web Conf._ 175 (2018) 02010, [1712.04884].
  * [24] J. Heitger, F. Joswig, A. Vladikas and C. Wittemeier, _Non-perturbative determination of $c_{V},Z_{V}$ and $Z_{S}/Z_{P}$ in $N_{f}=3$ lattice QCD_, _EPJ Web Conf._ 175 (2018) 10004, [1711.03924].
  * [25] J. Heitger and F. Joswig, _The renormalised $\mathrm{O}(a)$ improved vector current in three-flavour lattice QCD with Wilson quarks_, _Eur. Phys. J. C_ 81 (2021) 254, [2010.09539].
  * [26] J. Bulava, M. Della Morte, J. Heitger and C. Wittemeier, _Nonperturbative renormalization of the axial current in $N_{\mathrm{f}}=3$ lattice QCD with Wilson fermions and a tree-level improved gauge action_, _Phys. Rev. D_ 93 (2016) 114513, [1604.05827].
  * [27] J. Bulava, M. Della Morte, J. Heitger and C. Wittemeier, _Non-perturbative improvement of the axial current in $N_{\mathrm{f}}=3$ lattice QCD with Wilson fermions and tree-level improved gauge action_, _Nucl. Phys. B_ 896 (2015) 555, [1502.04999].
  * [28] L. Chimirri, P. Fritzsch, J. Heitger, F. Joswig, M. Panero, C. Pena et al., _Non-perturbative renormalization of $O(a)$ improved tensor currents_, _PoS_ LATTICE2019 (2020) 212, [1910.06759].
  * [29] M. Bruno, I. Campos, P. Fritzsch, J. Koponen, C. Pena, D. Preti et al., _Light quark masses in $N_{\mathrm{f}}=2+1$ lattice QCD with Wilson fermions_, _Eur. Phys. J. C_ 80 (2020) 169, [1911.08025].
  * [30] J. Heitger, F. Joswig and S. Kuberski, _Determination of the charm quark mass in lattice QCD with $2+1$ flavours on fine lattices_, 2101.02694.
  * [31] G. S. Bali et al., _The strange and light quark contributions to the nucleon mass from Lattice QCD_ , _Phys. Rev. D_ 85 (2012) 054502, [1111.1600].
  * [32] G. S. Bali, S. Collins, D. Richtmann, A. Schäfer, W. Söldner and A. Sternbeck, _Direct determinations of the nucleon and pion $\sigma$ terms at nearly physical quark masses_, _Phys. Rev. D_ 93 (2016) 094504, [1603.00827].
  * [33] S. Aoki et al., _FLAG Review 2019: Flavour Lattice Averaging Group (FLAG)_ , _Eur. Phys. J. C_ 80 (2020) 113, [1902.08191].
  * [34] K. Ottnad, _Excited states in nucleon structure calculations_ , in _38th International Symposium on Lattice Field Theory_ , 11, 2020.  2011.12471.
  * [35] J. Green, _Systematics in nucleon matrix element calculations_ , _PoS_ LATTICE2018 (2018) 016, [1812.10574].
  * [36] M. Lüscher, S. Sint, R. Sommer and P. Weisz, _Chiral symmetry and $\mathrm{O}(a)$ improvement in lattice QCD_, _Nucl. Phys. B_ 478 (1996) 365, [hep-lat/9605038].
  * [37] J. Bulava and S. Schaefer, _Improvement of $N_{\mathrm{f}}=3$ lattice QCD with Wilson fermions and tree-level improved gauge action_, _Nucl. Phys. B_ 874 (2013) 188, [1304.7093].
  * [38] M. Lüscher and S. Schaefer. http://luscher.web.cern.ch/luscher/openQCD.
  * [39] A. D. Kennedy, I. Horvath and S. Sint, _A New exact method for dynamical fermion computations with nonlocal actions_ , _Nucl. Phys. Proc. Suppl._ 73 (1999) 834, [hep-lat/9809092].
  * [40] M. A. Clark and A. D. Kennedy, _Accelerating dynamical fermion computations using the rational hybrid Monte Carlo (RHMC) algorithm with multiple pseudofermion fields_ , _Phys. Rev. Lett._ 98 (2007) 051601, [hep-lat/0608015].
  * [41] M. Lüscher, S. Sint, R. Sommer, P. Weisz and U. Wolff, _Non-perturbative O(a) improvement of lattice QCD_ , _Nucl. Phys. B_ 491 (1997) 323, [hep-lat/9609035].
  * [42] P. Perez-Rubio, S. Sint and S. Takeda, _An O(a) modified lattice set-up of the Schrödinger functional in SU(3) gauge theory_ , _JHEP_ 07 (2011) 116, [1105.0110].
  * [43] L. Del Debbio, H. Panagopoulos and E. Vicari, _$\theta$ dependence of SU(N) gauge theories_, _JHEP_ 08 (2002) 044, [hep-th/0204125].
  * [44] P. Fritzsch, A. Ramos and F. Stollenwerk, _Critical slowing down and the gradient flow coupling in the Schrödinger functional_ , _PoS_ Lattice2013 (2014) 461, [1311.7304].
  * [45] U. Wolff, _Monte Carlo errors with less errors_ , _Comput. Phys. Commun._ 156 (2004) 143, [hep-lat/0306017]. [Erratum: Comput. Phys. Commun. 176, 383 (2007)].
  * [46] S. Schaefer, R. Sommer and F. Virotta, _Critical slowing down and error analysis in lattice QCD simulations_ , _Nucl. Phys. B_ 845 (2011) 93, [1009.5228].
  * [47] A. Ramos, _Automatic differentiation for error analysis of Monte Carlo data_ , _Comput. Phys. Commun._ 238 (2019) 19, [1809.01289].
  * [48] P. T. Boggs and J. E. Rogers, _Orthogonal distance regression_ , tech. rep., National Institute of Standards and Technology, Gaithersburg, MD, 1989.  10.6028/NIST.IR.89-4197.
  * [49] M. Dalla Brida, T. Korzec, S. Sint and P. Vilaseca, _High precision renormalization of the flavour non-singlet Noether currents in lattice QCD with Wilson quarks_ , _Eur. Phys. J. C_ 79 (2019) 23, [1808.09236].
  * [50] A. Gerardin, T. Harris and H. B. Meyer, _Nonperturbative renormalization and $O(a)$-improvement of the nonsinglet vector current with $N_{f}=2+1$ Wilson fermions and tree-level Symanzik improved gauge action_, _Phys. Rev. D_ 99 (2019) 014519, [1811.08209].
  * [51] S. Sint and P. Weisz, _Further results on O(a) improved lattice QCD to one loop order of perturbation theory_ , _Nucl. Phys. B_ 502 (1997) 251, [hep-lat/9704001].
  * [52] S. Aoki, K.-i. Nagai, Y. Taniguchi and A. Ukawa, _Perturbative renormalization factors of bilinear quark operators for improved gluon and quark actions in lattice QCD_ , _Phys. Rev. D_ 58 (1998) 074505, [hep-lat/9802034].