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

Abstract

We determine non-perturbatively the normalisation factor $r_{\mathrm{m}}\equiv Z_{\mathrm{S}}/Z_{\mathrm{S}}^{0}$, where $Z_{\mathrm{S}}$ and $Z_{\mathrm{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_{\mathrm{f}}= 3$ mass-degenerate ${\mathrm{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_{\mathrm{P}}/(Z_{\mathrm{S}}Z_{\mathrm{A}})$ in order to obtain $r_{\mathrm{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_{\mathrm{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_{\mathrm{S}}$ and $Z_{\mathrm{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¹; see Refs. [1, 2]. Since $Z_{\mathrm{P}}$ and $Z_{\mathrm{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_{\mathrm{S}}/Z_{\mathrm{P}}$.² 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_{\mathrm{S}}/Z_{\mathrm{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_{\mathrm{f}}= 2$ QCD) [12], and with three dynamical quarks ($N_{\mathrm{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_{\mathrm{m}}\equiv Z_{\mathrm{S}}/Z_{\mathrm{S}}^{0}$. For chirally symmetric regularisations $r_{\mathrm{m}}=1$ holds, while for Wilson fermions $r_{\mathrm{m}}$ is a (finite) polynomial of the gauge coupling, arising from the sea fermion loops of the quark propagator. In the quenched approximation, $r_{\mathrm{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_{\mathrm{m}}= 1 + {\mathrm{O}}(g_0^4)$. In Ref. [18] the ${\mathrm{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_{\mathrm{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_{\mathrm{f}}= 2 + 1$ flavours, on large physical volumes, for a range of bare couplings corresponding to a hadronic regime [13, 21,22,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_{\mathrm{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. [14, 15, 24,25,26,27,28] that have various applications in lattice QCD when using this discretisation of Wilson fermions. The present work provides high-precision estimates of $r_{\mathrm{m}}$ obtained in the same computational framework.

As seen from Eq. (2.2) below, $(r_{\mathrm{m}}-1)$ contributes an ${\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{S}}^{0}$ is not as straightforward as that of $Z_{\mathrm{S}}$, the former also requiring the computation of two-boundary (“disconnected”) quark diagrams. This problem is circumvented by extracting $Z_{\mathrm{S}}^{0}$ as the product of $Z_{\mathrm{S}}$ and $r_{\mathrm{m}}$.

Our computation of $r_{\mathrm{m}}$ is based on the relation between the current (PCAC) mass m and the subtracted quark mass $m_{\mathrm{q}}$. Close to the chiral limit, $m(m_{\mathrm{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_{\mathrm{q}})$ is $Z r_{\mathrm{m}}$, where $Z \equiv Z_{\mathrm{P}}/(Z_{\mathrm{S}}Z_{\mathrm{A}})$ and $Z_{\mathrm{A}}$ is the non-singlet axial current normalisation. On the other hand, in a non-unitary theory with chiral valence subtracted quark masses ($m_{\mathrm{q}}^{\mathrm{val}} = 0$) and small degenerate sea quark masses $m_{\mathrm{q}}^{\mathrm{sea}} \ne 0$, the slope of $m(m_{\mathrm{q}}^{\mathrm{sea}})$ is $Z(r_{\mathrm{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_{\mathrm{m}}(g_0^2)$. This approach is described in Sect. 2.

Alternatively, each of the two slopes $Z r_{\mathrm{m}}$ and $Z(r_{\mathrm{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 Sect. 3.

In Sect. 4 we present our simulation setup for $N_{\mathrm{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. [14, 15, 25,26,27,28]. Some new ensembles have also been generated, in order to cover the region close to the origin of the function $m(m_{\mathrm{q}})$ more evenly and assess its slope reliably.

Our results for $r_{\mathrm{m}}$, based on various combinations of $Z r_{\mathrm{m}}$, $Z(r_{\mathrm{m}}-1)$, and Z are discussed in Sect. 5. Different determinations of $r_{\mathrm{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_{\mathrm{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_{\mathrm{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 Appendix   B. A 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_{\mathrm{f}}$,

$$\begin{aligned} m_{{\mathrm{q}},i}\equiv m_{0,i} - m_{\mathrm{crit}} = \dfrac{1}{2a} \Big (\dfrac{1}{\kappa _i}-\dfrac{1}{\kappa _{\mathrm{crit}}} \Big ) \, , \end{aligned}$$

(2.1)

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

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

(2.2)

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

We recall in passing that the renormalisation parameter $Z_{\mathrm{m}}(g_0^2,a\mu )$ depends on the renormalisation scale $\mu$ and diverges logarithmically in the ultraviolet. It is the inverse of $Z_{\mathrm{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_{\mathrm{P}}, Z_{\mathrm{m}}, Z_{\mathrm{S}}$), current normalisations (i.e. $Z_{\mathrm{A}},Z_{\mathrm{V}}$) and $r_{\mathrm{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 ${\mathrm{O}}(am_{{\mathrm{q}},i})$, ${\mathrm{O}}(a {\mathrm{Tr}}M_{\mathrm{q}})$, and ${\mathrm{O}}(a \Lambda _{\mathrm{QCD}})$ discretisation effects, which are part of their operational definition. As pointed out in Ref. [17], the term $(r_{\mathrm{m}}-1)$ multiplies ${\mathrm{Tr}}M_{\mathrm{q}}$, so it arises from a mass insertion in a quark loop. In perturbation theory it is a two-loop effect, contributing at ${\mathrm{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,\mathrm R}$ goes to the chiral limit only when all subtracted masses $m_{\mathrm{q},1}, \dots , m_{\mathrm{q},N_{\mathrm{f}}}$ vanish.

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

$$\begin{aligned} ({\widetilde{\partial }}^{}_{\mu })_x\big \langle (A_{\mathrm{I}})^{ij}_\mu (x)\,\mathcal {O}^{ji}\big \rangle = 2m_{ij} {\big \langle {P}^{ij}(x)\,\mathcal {O}^{ji}\big \rangle } \,. \end{aligned}$$

(2.3)

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

$$\begin{aligned} A_{\mu }^{ij}(x)\equiv & {} {\bar{\psi }}_i(x)\,\gamma _{\mu }\gamma _{5}\,\psi _j(x)\,,\nonumber \\ \qquad P^{ij}(x)\equiv & {} {\bar{\psi }}_i(x)\,\gamma _{5}\,\psi _j(x) \,, \end{aligned}$$

(2.4)

with indices i, j denoting two distinct flavours ($i \ne j$). The pseudoscalar density $P^{ij}$ and the current $(A_{\mathrm{I}})^{ij}_\mu \equiv A^{ij}_\mu +ac_{\mathrm{A}}{\widetilde{\partial }}^{}_{\mu }P^{ij}$ are Symanzik-improved in the chiral limit, with the improvement coefficient $c_{\mathrm{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.³ 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 ${\mathrm{O}}(a)$ improved theory, the renormalised axial current and pseudoscalar density are

$$\begin{aligned} (A_{\mathrm{R}})^{ij}(x)&= Z_{\mathrm{A}}(g_{0}^{2}) (A_{\mathrm{I}})_{\mu }^{ij}(x) + {\mathrm{O}}(am_\mathrm {q},a^2)\,, \end{aligned}$$

(2.5)

$$\begin{aligned} (P_{\mathrm{R}})^{ij}(x)&= Z_{\mathrm{P}}(g_{0}^{2},a\mu ) P^{ij}(x) + {\mathrm{O}}(am_\mathrm {q},a^2)\,. \end{aligned}$$

(2.6)

The normalisation of the axial current $Z_{\mathrm{A}}(g_0^2)$ is scale independent, depending only on the squared gauge coupling $g_0^2$. The renormalisation parameter $Z_{\mathrm{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,

$$\begin{aligned} ({\widetilde{\partial }}^{}_{\mu })_x \big \langle \, (A_{\mathrm{R}})^{ij}_\mu (x)\ \mathcal {O}^{ji}\, \big \rangle&= (m_{\mathrm{R},i}+m_{\mathrm{R},j})\,\big \langle \, (P_\mathrm {R})^{ij}(x)\ \mathcal {O}^{ji}\,\big \rangle \,, \end{aligned}$$

(2.7)

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

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

(2.8)

If we calculate the average mass $(m_{i\mathrm R} + m_{j\mathrm 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:

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

(2.9)

where the product of the renormalisation parameters $Z(g_0^2) \equiv Z_{\mathrm{P}}(g_0^2,\mu )/(Z_{\mathrm{S}}(g_0^2,\mu )Z_{\mathrm{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_{{\mathrm{q}},i}= m_{\mathrm{q},j}= {\mathrm{Tr}}M_{\mathrm{q}}/N_\mathrm{f}$), it reduces to

  $$\begin{aligned} m&= Z r_{\mathrm{m}}m_{\mathrm{q}}+ {\mathrm{O}}(am_\mathrm {q},a^2) \end{aligned}$$

  (2.10)
  $$\begin{aligned}&= Z r_{\mathrm{m}}\dfrac{1}{2a} \Big (\dfrac{1}{\kappa }-\dfrac{1}{\kappa _{\mathrm{crit}}} \Big ) + {\mathrm{O}}(am_\mathrm {q},a^2)\,. \end{aligned}$$

  (2.11)

  In the above equation, flavour indices have been dropped from the quark masses $m_{ij}, m_{{\mathrm{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 _{\mathrm{crit}}$, we obtain the latter as the root of the function $am(1/\kappa )$ and the combination $Z r_{\mathrm{m}}$ as the slope of the same curve.

• (2) Once the critical hopping parameter $\kappa _{\mathrm{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_{{\mathrm{q}},i}^{\mathrm{val}} \ne m_{{\mathrm{q}},i}^{\mathrm{sea}}$. In Eq. (2.9), masses $m_{{\mathrm{q}},i}$ and $m_{\mathrm{q},j}$ on the r.h.s. are valence quark contributions, while ${\mathrm{Tr}}M_{\mathrm{q}}$ stands for the trace of sea quark masses; see Refs. [14, 17] for detailed explanations. In particular, we set $\kappa ^{\mathrm{val}}= \kappa _{\mathrm{crit}}$, so as to ensure $m_{\mathrm{q}}^{\mathrm{val}}=0$ for all valence flavours. Moreover sea quark masses are taken to be small, degenerate, and non-zero (i.e. $\kappa ^{\mathrm{sea}}\ne \kappa _{\mathrm{crit}}$, ensuring $m_{\mathrm{q}}^{\mathrm{sea}} \ne 0$ for all sea flavours). With these conditions, the current quark mass of Eq. (2.9) reduces to

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

  (2.12)

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

The two slopes $Z r_{\mathrm{m}}$ and $Z (r_{\mathrm{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_{\mathrm{m}}(g_0^2)$; see Sect. 5.3 for details. We stress that the above discussion concerns relations which suffer from ${\mathrm{O}}(a)$ discretisation effects. For the quark masses, such effects may be removed by introducing Symanzik counterterms, leaving us with ${\mathrm{O}}(a^2)$ discretisation errors. These counterterms have been worked out in Refs. [17, 36]. In Ref. [17] (see also Eq. (2.10) of Ref. [14]) the full ${\mathrm{O}}(am_{\mathrm{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 ${\mathrm{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_{\mathrm{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_{\mathrm{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_{\mathrm{f}}= 2 + 1$ quark masses) are to be varied, while staying on lines of constant physics within systematic uncertainties of ${\mathrm{O}}(a^2)$. In particular, if the improved bare gauge coupling [36]

$$\begin{aligned} {\tilde{g}}_0^2 \equiv g_{0}^{2}\Big ( 1 + \dfrac{1}{N_{\mathrm{f}}} b_g(g_{0}^{2}) a {\mathrm{Tr}}M_{\mathrm{q}}\Big ) \end{aligned}$$

(2.13)

is kept fixed in the simulations, so is the lattice spacing, with fluctuations being attributed to ${\mathrm{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_{\mathrm{A}}(\tilde{g}_0^2), Z_{\mathrm{P}}(\tilde{g}_0^2,a\mu ), Z(\tilde{g}_0^2), r_{\mathrm{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 Sect. 6.

3 Ward identity determination of Z

In the previous section we have shown how the quantities $Z r_{\mathrm{m}}$ and $Z(r_{\mathrm{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_{\mathrm{m}}$ and Z. The latter quantity has already been measured in our setup ($N_\mathrm{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_{\mathrm{P}}/(Z_{\mathrm{A}}Z_{\mathrm{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 ${\mathcal{O}}\equiv S^b(y) {\mathcal{O}}^c$, defined as

$$\begin{aligned} S^b(y)\equiv & {} \mathrm {i} {\bar{\psi }}(y) T^b \psi (y) \nonumber \\ {\mathcal{O}}^c\equiv & {} \mathrm {i} \dfrac{a^6}{L^3} \sum _\mathbf{u,v} {\bar{\zeta }}(\mathbf{u}) \gamma _5 T^c \zeta (\mathbf{v})\,, \end{aligned}$$

(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.⁴ The Ward identity of interest is obtained by performing axial variations on ${\mathcal{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 ${\mathcal{O}}^c$ lying outside R, we have $\delta _{\mathrm{A}} {\mathcal{O}}= [\delta _{\mathrm{A}} S^b(y)]{\mathcal{O}}^c$ and

$$\begin{aligned} \delta _{\mathrm{A}} S^b(x) = \epsilon ^a \Big [ d^{abe} P^e(x) + \dfrac{\delta ^{ab}}{N_\mathrm{f}} {\bar{\psi }}(x) \psi (x) \Big ] \, . \end{aligned}$$

(3.2)

In what follows we simplify matters by always working with $a \ne 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

$$\begin{aligned}&\int \mathrm {d}^3\mathbf{y} \int \mathrm {d}^3\mathbf{x} \Big \langle \Big [A_0^a(t_2;\mathbf{x}) - A_0^a(t_1;\mathbf{x}) \Big ] S^b(y_0;\mathbf{y}) {\mathcal{O}}^c \Big \rangle \nonumber \\&\qquad - 2m \int \mathrm {d}^3\mathbf{y} \int \mathrm {d}^3\mathbf{x} \int _{t_1}^{t_2} \mathrm {d}x_0 \langle P^a(x_0;\mathbf{x}) S^b(y_0;\mathbf{y}) {\mathcal{O}}^c \rangle \nonumber \\&\quad = - d^{abe} \int \mathrm {d}^3\mathbf{y} \,\, \langle P^e(y) {\mathcal{O}}^c \rangle \, . \end{aligned}$$

(3.3)

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

$$\begin{aligned}&\!\!\! Z_{\mathrm{A}}Z_{\mathrm{S}}a^6 \Bigg \{ \sum _{\mathbf{x},\mathbf{y}} \, \left\langle \Big [ (A_{\mathrm{I}})^a_0(t_2; \mathbf{x}) - (A_{\mathrm{I}})^a_0(t_1; \mathbf{x}) \Big ]\, S^b(y_0;\mathbf{y}) \, {\mathcal{O}}^c \right\rangle \nonumber \\&\!\!\! \quad - 2 a m \sum _{\mathbf{x},\mathbf{y}} \sum _{x_0 = t_1}^{t_2} w(x_0) \, \langle P^a(x_0;\mathbf{x}) \, S^b(y_0;\mathbf{y}) \, {\mathcal{O}}^c \rangle \Bigg \} \nonumber \\&\!\!\! \qquad = - d^{abe} Z_{\mathrm{P}}\,\, a^3 \sum _\mathbf{y} \langle \, P^e(y) \, {\mathcal{O}}^c \rangle + {\mathrm{O}}(am,a^2) \,. \end{aligned}$$

(3.4)

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

$$\begin{aligned} Z\equiv & {} \dfrac{Z_{\mathrm{P}}}{Z_{\mathrm{A}}Z_{\mathrm{S}}} \nonumber \\= & {} {-} \dfrac{f_\mathrm{A S}^\mathrm {I}(t_2,y_0) - f_\mathrm{A S}^\mathrm {I}(t_1,y_0) - 2 am {\tilde{f}}_{\mathrm{P S}}(t_2,t_1,y_0)}{f_\mathrm{P}(y_0)} \nonumber \\&+ \,\, \mathrm O(am,a^2)\,, \end{aligned}$$

(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 ${\mathrm{O}}(am)$ effects may be safely neglected), the above Ward identity is valid up to ${\mathrm{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.⁵ The renormalisation factor of the external source ${\mathcal{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. [15, 26] 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 Sect. 5.2 and Fig. 5.

Fig. 1

The trace diagrams contributing to the expectation values of $f_\mathrm{P}$, defined in Eq. (B.1) (a) and $f_\mathrm{A S}$, defined in Eq. (B.3) (b, 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) (a) or a scalar operator S(y) (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

Fig. 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

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 Sect. 5.2 and Appendix   C.

4 Numerical setup

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 \, \hbox {fm}$ 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_{\mathrm{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 ${\mathrm{O}}(a)$ cut-off effects across the coupling range also considered in the present work. We thus expect our final results for $r_{\mathrm{m}}$ and Z to only be affected by ${\mathrm{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.⁶

The gauge ensembles highlighted in italics were newly generated for this study, while the remaining ones were already used in previous investigations; see Refs. [14, 15, 24,25,26,27,28].⁷ 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 Sect. 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_{\mathrm{t}}, \widetilde{c}_{\mathrm{t}}$ values) as if the time extents were even. The fact that an odd time extent alters the tree-level value of $c_{\mathrm{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_{\mathrm{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 [13, 21, 23].

5.1 Quark mass slopes

As described in Sect. 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

$$\begin{aligned} m(x_0) = \frac{{\widetilde{\partial }}^{}_{0} f_\mathrm{A}^{ij}(x_0)+ac_\mathrm {A} \partial _0^*\partial _0^{}f_\mathrm{P}^{ij}(x_0)}{2f_\mathrm{P}^{ij}(x_0)}\,, \end{aligned}$$

(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_\mathrm{A}^{ij}(T-x_0)$ and $g_\mathrm{P}^{ij}(T-x_0)$, which are constructed from the same operators $(A_{\mathrm{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.

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} \ne \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.

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 L m \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. 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.

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} \ne \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^{\lbrace T/3 \rbrace }}$ and ${Z^{\lbrace T/4 \rbrace }}$ refer to different choices of time slices with operator insertions in the correlation functions (see text for details)

Fig. 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

Fig. 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

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 ${\mathrm{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 ${\mathrm{O}}((a m_{\mathrm{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 $a m_{\mathrm{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.

5.2 Renormalisation constant Z

As the next step in our analysis, we extract the renormalisation constant $Z \equiv (Z_{\mathrm{P}}/Z_{\mathrm{S}}Z_{\mathrm{A}})$ from the ratio (3.5), using the subset of gauge field ensembles listed in Table 1 which are not emphasised in italic font.⁸ 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^{\lbrace T/3 \rbrace }}$. Alternatively, choosing $t_1 \approx T/4$ and $t_2 \approx 3T/4$ yields a second Z estimate denoted as ${Z^{\lbrace T/4 \rbrace }}$. When T/3 and T/4 are not integers, $t_1$ and $t_2$ are rounded up/down to the nearest integer.

In the left part of Fig. 5, we depict ${Z^{\lbrace T/3 \rbrace }}$ 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^{\lbrace T/3 \rbrace }}$ and ${Z^{\lbrace T/4 \rbrace }}$ are not compatible, indicating the presence of lattice artefacts that also differ noticeably. We consider ${Z^{\lbrace T/3 \rbrace }}$ 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” $2 am {\tilde{f}}_{\mathrm{P S}}(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^{\lbrace T/3 \rbrace }}$ and ${Z^{\lbrace T/4 \rbrace }}$ 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. [13,14,15] can be found in Appendix   C.

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

Fig. 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]

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}}^{\lbrace {\mathrm {u,nu}}\rbrace }}$, defined as

$$\begin{aligned} {r_{\mathrm {m}}^{\lbrace {\mathrm {u,nu}}\rbrace }}&= \bigg (1-\left[ \frac{Z(r_\text {m}-1)}{Zr_\text {m}}\right] \bigg )^{-1} \,, \end{aligned}$$

(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}^{\lbrace \mathrm {u;}Z\rbrace }}$, defined as

$$\begin{aligned} {r_\mathrm {m}^{\lbrace \mathrm {u;}Z\rbrace }}&= \frac{Zr_\text {m}}{Z} \,. \end{aligned}$$

(5.3)

As mentioned above, this comes in two versions, ${r_\mathrm {m}^{\lbrace \mathrm {u;}Z,T/3\rbrace }}$ and ${r_\mathrm {m}^{\lbrace \mathrm {u;}Z,T/4\rbrace }}$. Moreover, from the second and third result we gain ${r_\mathrm {m}^{\lbrace \mathrm {nu;}Z\rbrace }}$ given by

$$\begin{aligned} {r_\mathrm {m}^{\lbrace \mathrm {nu;}Z\rbrace }}&= \frac{Z(r_\text {m}-1)}{Z} + 1 \,, \end{aligned}$$

(5.4)

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

In principle, the different estimates can differ by $\mathrm {O}(a^2)$ ambiguities. In Fig. 6 (left) the three determinations ${r_{\mathrm {m}}^{\lbrace {\mathrm {u,nu}}\rbrace }}$, ${r_\mathrm {m}^{\lbrace \mathrm {u;}Z,T/3\rbrace }}$, and ${r_\mathrm {m}^{\lbrace \mathrm {nu;}Z,T/3\rbrace }}$ 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}}^{\lbrace {\mathrm {u,nu}}\rbrace }}$ by an improved version of the derivative ${\widetilde{\partial }}^{}_{0}$ in Eq. (5.1).⁹ 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.

Table 4 Results for $r_\mathrm {m}$, obtained via Eqs. (5.2)–(5.4)

Fig. 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}^{\lbrace \mathrm {nu;}Z,T/3\rbrace }}$ 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. [25, 49]). 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).

Fig. 7

Non-perturbative determination of ${r_\mathrm {m}^{\lbrace \mathrm {nu;}Z,T/3\rbrace }}$ (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}^{\lbrace \mathrm {nu;}Z,T/3\rbrace }}$. Our data is best described by a Padé ansatz, constrained to the two-loop prediction of Ref. [18] for small couplings, of the form

$$\begin{aligned} 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\} \,,\nonumber \\ \end{aligned}$$

(5.5a)

where

$$\begin{aligned} \quad c_i = \left( -7.86078, 5.49175, -0.54078\right) \,, \end{aligned}$$

(5.5b)

and

$$\begin{aligned}&\mathrm {cov}(c_i,c_j) \nonumber \\&= \left( \begin{array}{l@{\quad }l@{\quad }l} 3.699760 &{}-2.198586 &{}-1.476913 \times 10^{-3}\\ -2.198586 &{} 1.306512 &{} 8.776569 \times 10^{-4}\\ -1.476913 \times 10^{-3} &{} 8.776569 \times 10^{-4} &{} 5.895710 \times 10^{-7} \end{array}\right) \,,\nonumber \\ \end{aligned}$$

(5.5c)

which is also displayed in Fig. 7. The fit function describes our data with $\chi ^2/\mathrm{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 [13, 21, 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.

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_{\mathrm{m}}\equiv Z_{\mathrm{S}}/Z_{\mathrm{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_{\mathrm{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_{\mathrm{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 ${\mathrm{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_{\mathrm{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 ${\mathrm{O}}(a)$ improved theory.

Our central results, which hold for a lattice discretisation of QCD with three flavours of non-perturbatively ${\mathrm{O}}(a)$ improved Wilson-clover sea quarks and tree-level Symanzik-improved gluons, are the continuous parameterisation of $r_{\mathrm{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 [13, 21,22,23].

Along with the numerical implementation of our strategy to extract $r_{\mathrm{m}}$, we have also developed a new method to determine the scale independent combination $Z=Z_{\mathrm{P}}/(Z_{\mathrm{S}}Z_{\mathrm{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 Sects. 2 and 3, the renormalisation factor Z is actually required to isolate $r_{\mathrm{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_{\mathrm{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 Fig. 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 Sect. 2: away from the chiral limit, the dependence of (re)normalisation parameters should be $Z({\tilde{g}}_0^2), r_{\mathrm{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

$$\begin{aligned} 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_{\mathrm{f}}}b_g(g_{0}^{2}) g_{0}^{2}a {\mathrm{Tr}}M_{\mathrm{q}}\Big ] \, , \end{aligned}$$

(6.1)

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

Appendices

Appendix A: Schrödinger functional correlation functions

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

$$\begin{aligned} f_\mathrm {P}^{ij}&= -\frac{1}{2} \frac{a^9}{L^3} \sum _\mathbf{x, u, v} \left\langle {\bar{\psi }}_i(x) \gamma _5 \psi _j(x) \cdot {\bar{\zeta }}_j(\mathbf{v}) \gamma _5 \zeta _i(\mathbf{u})\right\rangle \,, \end{aligned}$$

(A.1)

$$\begin{aligned} g_\mathrm {P}^{ij}&= -\frac{1}{2} \frac{a^9}{L^3} \sum _\mathbf{x, u, v} \left\langle {\bar{\psi }}_i(x) \gamma _5 \psi _j(x) \cdot {\bar{\zeta }}'_j(\mathbf{u}) \gamma _5 \zeta '_i(\mathbf{v})\right\rangle \,, \end{aligned}$$

(A.2)

$$\begin{aligned} f_\mathrm {A}^{ij}&= -\frac{1}{2} \frac{a^9}{L^3} \sum _\mathbf{x, u, v} \left\langle {\bar{\psi }}_i(x) \gamma _0 \gamma _5 \psi _j(x) \cdot {\bar{\zeta }}_j(\mathbf{u}) \gamma _5 \zeta _i(\mathbf{v})\right\rangle \,, \end{aligned}$$

(A.3)

$$\begin{aligned} g_\mathrm {A}^{ij}&= -\frac{1}{2} \frac{a^9}{L^3} \sum _\mathbf{x, u, v} \left\langle {\bar{\psi }}_i(x) \gamma _0 \gamma _5 \psi _j(x) \cdot {\bar{\zeta }}'_j(\mathbf{u}) \gamma _5 \zeta '_i(\mathbf{v})\right\rangle \,,\end{aligned}$$

(A.4)

$$\begin{aligned} F_1^{ij}&= -\frac{1}{2} \frac{a^{12}}{L^6} \sum _\mathbf{u', v', u, v} \left\langle {\bar{\zeta }}'_i(\mathbf{u'}) \gamma _5 \zeta '_j(\mathbf{v'}) \cdot {\bar{\zeta }}_j(\mathbf{u}) \gamma _5 \zeta _i(\mathbf{v})\right\rangle \,. \end{aligned}$$

(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(\mathbf{u})$ and $\zeta _i(\mathbf{v})$ live on time slice $x_0=0$, while ${\bar{\zeta }}'_j(\mathbf{u}')$ and $\zeta '_i(\mathbf{v}')$ 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

$$\begin{aligned}&- d^{abe} \, a^3 \sum _\mathbf{y} \langle P^e(y) {\mathcal{O}}^c \rangle \nonumber \\&\quad = - d^{abe} {\mathrm{Tr}}[T^e T^c] \dfrac{a^{9}}{L^3} \sum _\mathbf{y} \sum _\mathbf{u,v} \nonumber \\&\qquad \times \Bigg \langle \,\hbox {tr}\,\bigg \{ [\psi (y) {\bar{\zeta }}(\mathbf{u})]_{\mathrm{F}} \gamma _5 [\zeta (\mathbf{v}) {\bar{\psi }}(y)]_{\mathrm{F}} \gamma _5 \bigg \} \Bigg \rangle \nonumber \\&\quad =\,d^{abc} f_\mathrm{P}(y_0) \,, \end{aligned}$$

(B.1)

where the second equality implicitly defines $f_\mathrm{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:

$$\begin{aligned}&a^6 \sum _\mathbf{x,y} \langle A_0^a(x) S^b(y) {\mathcal{O}}^c \rangle \nonumber \\&\quad = \dfrac{\mathrm {i} a ^{12}}{L^3} {\mathrm{Tr}}[T^a T^b T^c] \sum _\mathbf{x,y} \sum _\mathbf{u,v} \nonumber \\&\qquad \times \Bigg \langle \,\hbox {tr}\,\bigg \{ \gamma _0 \gamma _5 [\psi (x) {\bar{\psi }}(y)]_{\mathrm{F}} [\psi (y) {\bar{\zeta }}(\mathbf{u})]_{\mathrm{F}} \gamma _5 [\zeta (\mathbf{v}) {\bar{\psi }}(x)]_{\mathrm{F}} \bigg \} \Bigg \rangle \nonumber \\&\qquad +\dfrac{\mathrm {i} a ^{12}}{L^3} {\mathrm{Tr}}[T^c T^b T^a] \sum _\mathbf{x,y} \sum _\mathbf{u,v} \nonumber \\&\qquad \times \Bigg \langle \,\hbox {tr}\,\bigg \{ \gamma _0 \gamma _5 [\psi (x) {\bar{\zeta }}(\mathbf{u})]_{\mathrm{F}} \gamma _5 [\zeta (\mathbf{v}) {\bar{\psi }}(y)]_{\mathrm{F}} [\psi (y) {\bar{\psi }}(x)]_{\mathrm{F}} \bigg \} \Bigg \rangle \nonumber \\&\quad = \, \mathrm {i} {\mathrm{Tr}}[T^a T^b T^c] f_{\mathrm{A S;1}} (x_0,y_0) + \mathrm {i} {\mathrm{Tr}}[T^c T^b T^a] f_{\mathrm{A S;2}} (x_0,y_0) \nonumber \\&\quad =\, \dfrac{1}{2} \bigg [ -d^{abc} \mathrm{Re}\,f_{\mathrm{A S;1}} (x_0,y_0) + f^{abc} \mathrm{Im}\,f_{\mathrm{A S;1}} (x_0,y_0) \bigg ] \,. \end{aligned}$$

(B.2)

The second in the above string of equations implicitly defines the two traces of quark propagators (devoid of flavour structure) as $f_{\mathrm{A S;1}} (x_0,y_0)$ and $f_{\mathrm{A S;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

$$\begin{aligned}&a^6 \sum _\mathbf{x,y} \langle A_0^a(x) S^b(y) {\mathcal{O}}^c \rangle = - \dfrac{1}{2} d^{abc} \mathrm{Re}\,f_{\mathrm{A S;1}} (x_0,y_0) \nonumber \\&\quad = - d^{abc} f_\mathrm{A S}(x_0,y_0) \,, \end{aligned}$$

(B.3)

which implicitly defines $f_\mathrm{A S}(x_0,y_0)$. The correlation functions $f_\mathrm{P}$ and $f_\mathrm{A S}$ are schematically drawn in Fig. 1. Analogously, from the mass dependent term of the Ward identity we also define ${\tilde{f}}_{\mathrm{P S}}(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).

Table 6 Comparison of our preferred Z determination with results from Ref. [14] (de Divitiis et al.) and the combination of results from Refs. [24, 49] (Heitger et al.)

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^{\lbrace T/3 \rbrace }}$ 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. [15, 49] (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 ${\mathrm{O}}(a^2)$ or higher between different determinations. Agreement is generally better between our results and those of Ref. [14] (de Divitiis et al.).

Fig. 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^{\lbrace T/3 \rbrace }}$ (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^{\lbrace T/3 \rbrace }}$ (excluding the coarsest lattice spacing from the fit). The vertical dotted lines correspond to the bare couplings used in CLS simulations

Fig. 9

Left: ratio of our results ${Z^{\lbrace T/3 \rbrace }}/{Z^{\lbrace T/4 \rbrace }}$, fitted as a function of the lattice spacing. Right: ratios of ${Z^{\lbrace T/3 \rbrace }}$ 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

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^{\lbrace T/3 \rbrace }}$, obtained by differentiating Eq. (C.1a) as well as $\partial \ln r_{\mathrm{m}}/\partial g_0^2$ via differentiating Eq. (5.5)

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 ${\mathrm{O}}(a^2)$, leading cut-off effects of ${\mathrm{O}}(a^3)$ can be identified in the ratio ${Z^{\lbrace T/3 \rbrace }}/{Z^{\lbrace T/4 \rbrace }}$, 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 ${\mathrm{O}}(a^2)$ ambiguities or higher which by virtue of the imposed line of constant physics decrease monotically towards the continuum limit.

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

$$\begin{aligned} 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}\,, \end{aligned}$$

(C.1a)

where

$$\begin{aligned} \quad d_i = \left( -0.34504, -0.52309\right) \,, \end{aligned}$$

(C.1b)

and

$$\begin{aligned} \mathrm {cov}(d_i,d_j) = \left( \begin{array}{lll} 2.798505 \times 10^{-7}\,\,\,\, &{}6.577477 \times 10^{-7}\\ 6.577477 \times 10^{-7}\,\,\,\,&{}1.545940 \times 10^{-6} \end{array}\right) . \nonumber \\ \end{aligned}$$

(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^{\lbrace T/3 \rbrace }}$ data point: the systematic error of ${Z^{\lbrace T/3 \rbrace }}(\beta = 3.4)$ is the statistical error of ${Z^{\lbrace T/3 \rbrace }}(\beta = 3.414)$ and that of ${Z^{\lbrace T/3 \rbrace }}(\beta = 3.85)$ is the statistical error of ${Z^{\lbrace T/3 \rbrace }}(\beta = 3.81)$. These systematic errors are added to the statistical ones in quadrature.

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

$$\begin{aligned}&\mathrm {cov}(d_i,c_j) \nonumber \\&\quad = \left( \begin{array}{lll} 3.502738 \times 10^{-4}\,\,\,&{}-2.149363 \times 10^{-4}\,\,\,&{}-1.303033 \times 10^{-7}\\ 1.464606 \times 10^{-3}\,\,\,&{}-8.876895 \times 10^{-4}\,\,\,&{}-5.079411 \times 10^{-7} \end{array}\right) ,\nonumber \\ \end{aligned}$$

(C.2)

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