University of Southern Denmark Non-perturbative improvement of the axial current in N_f = 3 lattice QCD with Wilson fermions and tree-level improved gauge action

The coefficient cA required for O(a) improvement of the axial current in lattice QCD with Nf = 3 flavors of Wilson fermions and the tree-level Symanzik-improved gauge action is determined non-perturbatively. The standard improvement condition using Schrödinger functional boundary conditions is employed at constant physics for a range of couplings relevant for simulations at lattice spacings of ≈ 0.09 fm and below. We define the improvement condition projected onto the zero topological charge sector of the theory, in order to avoid the problem of possibly insufficient tunneling between topological sectors in our simulations at the smallest bare coupling. An interpolation formula for cA(g 0) is provided together with our final results. © 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3. E-mail address: heitger@uni-muenster.de (J. Heitger). http://dx.doi.org/10.1016/j.nuclphysb.2015.05.003 0550-3213/© 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3. 556 ALPHA Collaboration / Nuclear Physics B 896 (2015) 555–568


Introduction
Wilson fermions [1] are an attractive and popular fermion discretization for lattice QCD simulations. However, while the addition of the Wilson term lifts unwanted fermion 'doubler' modes, it also leads to O(a) cutoff effects. At small a, these cutoff effects are described by the (continuum) Symanzik effective theory [2,3], and the procedure of their systematic removal is called 'O(a) improvement'.
In order to eliminate these O(a) cutoff effects, a single dimension-five 'clover' term must be added to the Lagrangian [4] (with coefficient c sw ), and additional counterterms have to be added to local operators. To this end, the (unrenormalized) improved axial current is given by

1)
A a μ (x) =ψ(x)T a γ μ γ 5 ψ(x), P a =ψ(x)T a γ 5 ψ(x), (1.2) where T a is an SU(N f ) generator acting in flavor space and ∂ μ , ∂ * μ are the lattice forward and backward derivatives, respectively. Therefore, in order to improve matrix elements of the axial current, the coefficient c A must be specified as well as c sw . The non-perturbative determination of c A for a range of bare couplings is the subject of this work.
Matrix elements of the axial current are of particular importance, since they enter the computation of pseudoscalar meson decay constants and quark masses. These quantities are not only of great phenomenological interest in their own right, but it has also been demonstrated [5] that the kaon decay constant f K provides a precise determination of the lattice scale in physical units, thus affecting all dimensionful quantities. For instance, in the two-flavor theory and at lattice spacings of ≈ 0.1 fm, the c A -term discussed above typically contributes ≈ 10-15% to pseudoscalar meson decay constants and quark masses [6,7].
In this work we consider lattice QCD with N f = 3 mass-degenerate flavors of Wilson fermions and the tree-level Symanzik-improved gauge action [8]. This gauge action is demonstrably preferable in the pure gauge theory, where its cutoff effects are smaller than other standard actions [9]. In preparation for dynamical simulations in this setup (see [10] for a first report), the parameter c sw multiplying the dimension-five 'clover' term in the action has been tuned nonperturbatively [11]. Like c sw , the improvement coefficient c A has been determined in perturbation theory at 1-loop in Ref. [12]. However, previous non-perturbative determinations of c A [6,7] deviate roughly 300-400% from 1-loop perturbation theory at the largest bare couplings, so that a non-perturbative determination is required.
To determine c A non-perturbatively, we employ the (by now) standard improvement condition for dynamical fermions [6,7], which imposes the PCAC relation at constant physics to ensure a removal of O(a) effects in on-shell quantities and, at the same time, a smooth behavior of vanishing O(a 2 ) effects as the bare coupling is varied. However, for technical reasons related to the topology freezing of our simulations at the finest lattice spacing, we project the necessary correlation functions onto the trivial topological sector. The main result of this work is the interpolation formula for c A (g 2 0 ) given in Eq. (4.1), with the coefficients of Eq. (4.2), which is valid for lattice spacings of a ≈ 0.09 fm and below. A statistical error of 4% should be assigned to this formula near the largest simulated bare couplings, whereas 7% is more appropriate at the finest lattice spacing.
We detail the improvement condition as well as the projection onto the trivial topological sector in Section 2, and our simulation setup is discussed in Section 3. Numerical results together with the final interpolation formula are presented in Section 4, and conclusions are drawn in Section 5.

Improvement condition
We determine c A via a variant of the improvement condition originally introduced in quenched QCD in [13] and applied in its present form to the theory with dynamical fermions first in [6] for the two-flavor case. We briefly review this condition here, adopting the notation of Ref. [6].
Improvement conditions are typically based on imposing the PCAC relation, which in its continuum form is an operator identity, at finite lattice spacing. A consequence of this relation is that the PCAC quark mass, defined as is independent of the space-time position x as well as the external states α and β. It is this property of m PCAC , required to hold on the lattice up to O(a 2 ) cutoff effects, which we exploit to determine c A .
Employing the above relation with the improved axial current of Eq. (1.1), we can decompose As mentioned above, if the continuum form of the PCAC relation holds, m(x; α, β) = m is independent of α, β and x. If we demand this property and fix x, the quark mass m can be eliminated by using two different pairs of external states α, β and γ, δ to obtain our definition of c A : In practice we employ a finite physical system size (L ≈ 1.2 fm) with Schrödinger functional boundary conditions [14,15] in time and a periodic torus in space. Furthermore (using the standard notation), we employ a vanishing background gauge field (φ i = φ i = 0) and periodic spatial boundary conditions for the fermion fields (θ i = 0). With the tree-level Symanzik-improved gauge action, there are several possibilities for implementing boundary O(a) improvement in the Schrödinger functional. As in Ref. [11], we resort to 'choice B' of Ref. [12], which possesses the desirable property that the classical minimum of the action can be expressed analytically. Our boundary O(a) improvement is implemented at tree-level according to this choice. An additional boundary O(a) improvement possibility for this action in the Schrödinger functional has been proposed recently in Ref. [16].
Correlation functions in the Schrödinger functional setup may involve boundary quark fields, which are adopted here. We construct r and s defined above from the following correlation functions: where the pseudoscalar operator O a is composed of the boundary quark fields at x 0 = 0 (ζ and ζ ), while the 'wavefunction' ω( x) will be discussed shortly. Defining the corresponding operator O a (ω ) at x 0 = T , we also employ the boundary-to-boundary correlation function As outlined before, we aim at probing the PCAC relation with two different pairs of external states. To achieve this, we construct approximate wavefunctions of the ground and first excited state in the pseudoscalar channel, viz.
where ω π (0) (ω π (1) ) is constructed to maximize the overlap with the ground (first excited) state. These approximate wavefunctions are superpositions of trial wavefunctions with coefficients η (i) j , given by the eigenvectors of the matrix f ij = f 1 (ω i , ω j ) corresponding to the largest and next-to-largest eigenvalues for ω π (0) and ω π (1) , respectively. As our spatial trial (hydrogen-like) wavefunctions at the boundaries we takē where N i is a normalization factor chosen to ensure a 3 x ω 2 i ( x) = 1 and a 0 = L/6. Finally, our operational definition of c A now reads , (2.12) which at the same time defines r(x 0 ) and s(x 0 ). In this way, the wavefunctions ω π (0) (ω π (1) ) determine the states β (δ), cf. Eq. (2.4), while for both states α and γ the plain Schrödinger functional boundary state i 0 with vacuum quantum numbers is inserted. The choice of x 0 will be discussed later.
Unlike previous studies, we employ T = 3L/2 lattices for the generation of our dynamical gauge field ensembles, with the intention of re-using them for a determination of the axial current renormalization constant Z A according to the method of Ref. [17] (see [18] for a preliminary report). While in a quark mass independent improvement scheme as applied here one ideally would impose the improvement condition of Eq. (2.11) at zero quark mass, in practice it was found in Ref. [6] and is also confirmed by our data that c A is rather insensitive to fairly small deviations from this constraint. However, that is not the case for Z A , so in general we endeavor to tune the quark mass to values close to zero.
Like periodic temporal boundary conditions, Schrödinger functional boundary conditions are 'closed', and disconnected topological sectors emerge in the continuum limit. However, for small physical volumes, non-trivial topological sectors receive a small weight in the partition sum. Unfortunately, to keep O(a 2 ) effects under control, and with an eye toward using our ensembles for Z A , our physical volume (L ≈ 1.2 fm) is large enough that non-trivial topological sectors are not completely suppressed. Whereas these sectors are sampled sufficiently at coarser lattice spacings, ensembles at our finest lattice spacing (L/a = 24, β = 3.81) are effectively frozen in the sector with topological charge Q = 0.
The issue of topology freezing in the Schrödinger functional has been investigated recently in the pure gauge theory in Ref. [19]. There it was suggested that quantities projected to the zero topological sector have a smooth approach to the continuum limit. For the case at hand (c A as well as Z A ), quantities defined in this way differ from their conventional counterparts only by irrelevant cutoff effects. Formally, we perform this projection for all observables of interest. Adopting the notation of Ref. [19], we define where O is an arbitrary observable in the theory (such as the correlation functions in Eqs. (2.5) and (2.7)), and the topological charge Q is defined using the Wilson (gradient) flow [20,21] at flow time t given by √ 8t/L = c with c = 0.35. Since at finite lattice spacing the topological charge may take non-integer values, we define all configurations with |Q| ≤ 0.5 as the trivial topological sector. Thus (as in Ref. [19]) we make the replacement δ Q,0 → θ(Q + 0.5)θ (0.5 −Q). For ease of notation we shall henceforth take c A to mean the one projected onto the trivial sector in this manner, with the exception of Table 2, which directly compares results in all sectors (where available) with those restricted to Q = 0. Let us anticipate already here that these two kinds of analyses yield consistent results for c A as expected, because the Ward identities underlying the PCAC relation and thereby our improvement strategy hold in any topological sector.
An alternative to address topology freezing in the Schrödinger functional has recently appeared [16], namely the use of 'half-open' boundary conditions. While not considered here, these boundary conditions may help the problem in future calculations provided an improvement condition which does not require boundary-to-boundary correlation functions is devised.

Simulation details
Our simulations are performed using Schrödinger functional boundary conditions. We use the openQCD code 1 of Ref. [22], which was also used for the simulations to calculate c sw in Ref. [11].
The bare gauge couplings are chosen to approximately satisfy a constant physics condition, fixing L ≈ 1.2 fm. In this way it is ensured that any O(a) ambiguities in c A disappear smoothly toward the continuum limit. For a thorough an more general discussion of the idea and virtues of imposing improvement (and renormalization) conditions at constant physics, see, e.g., Refs. [23,24]. Similarly to earlier work [6,7], we fix the physical volume by beginning with a particular pair of g 2 0 and L/a (β = 6/g 2 0 = 3.3 at L/a = 12 in the present case) and choose the bare couplings for subsequent smaller lattice spacings according to the universal 2-loop β-function. 2 The range of lattice spacings covered in this way extends from a ≈ 0.09 fm to a ≈ 0.045 fm. As already mentioned, at each bare coupling we tune the bare quark mass so that the PCAC mass is Table 1 Summary of simulation parameters, number of replica and total number of molecular dynamics units of our gauge configuration ensembles labeled by 'ID'. kept approximately constant and close to zero. 3 At several lattice spacings we confirm that our determination of c A is insensitive to variations of the (small) quark mass. Information about our ensembles, consisting of several replica per parameter set in most cases, can be found in Table 1.
Due to practical reasons discussed in the openQCD documentation, 4 our lattices have temporal extents T = 3L/2 − a. Since this offset itself scales with a, one expects its influence on the determination of c A to be of O(a 2 ) and thus to be small; still, we assess it on our coarsest lattice where it is largest. There we simulate 12 3 × 17 as well as 12 3 × 19 ensembles.
Although it was found previously that small deviations from the constant physics condition have little effect on c A [6], we estimate this deviation by measuring the scale-dependent renormalized coupling ḡ 2 GF , defined in Ref. [25]. Results for this coupling are shown in Table 2. We also test the dependence of c A on L in physical units (and thereby on violations of the constant physics condition) directly by simulating an additional bare coupling at L/a = 16.
We now briefly summarize the simulation algorithm used for the generation of these gauge field ensembles. While two of the (mass-degenerate) pseudo-fermion fields can be simulated in the usual way, the RHMC algorithm [26] is employed for the third. Even-odd preconditioning is used for all fermion determinants, whereas mass preconditioning [27] with two additional pseudo-fermion fields is used for the degenerate doublet. We use a hierarchical integration scheme [28], where the gauge force is integrated on the innermost level and the remaining fermion forces on the second level. For the L/a = 16, 24 and part of the L/a = 20 lattices, the lowest poles of the RHMC are integrated on a third level. For the two inner levels, a fourthorder OMF integrator [29] is used, while the third level (when present) uses a second-order OMF integrator. A single step is used for the inner integrators, and the number of steps for the outer level is tuned to achieve an acceptance rate of ≈ 90%. For the coarsest L/a = 12 lattice ensembles A1k2 and A2k1, we adopt (type I) twisted mass reweighting [30]. The conjugate gradient solver is employed for most fermion forces, while the multi-shift variant is typically used for 3 Based on the experience from the two-flavor theory, for which the multiplicative quark mass renormalization factor only varies slowly with a, we can safely neglect it for the tuning purposes here, too. 4 For this work we employ openQCD version 1.2. This issue has been corrected in the latest version (1.4). Table 2 Summary of results for c A . The (unrenormalized) PCAC quark mass am PCAC is computed from the correlation functions projected to the approximate ground state, using the 1-loop result for c A (g 2 0 ) from [12], while ḡ 2 GF denotes the gradient (resp. Wilson) flow coupling mentioned in the text. Recall that quantities with the explicit subscript label '0' here refer to results from the analysis restricted to the sector of vanishing topological charge, whereas in the text we loosely suppress the '0'. Numbers for ensemble D1k1 (L/a = 24) are not quoted ('n.q.') for the case of covering all charge sectors in the partition sum, because our simulations are not able to sufficiently sample all the sectors and a reliable error estimation is thus not possible. Results from ensembles in italics enter into the final interpolation formula for c A (g 2 0 ), Eq. most of the RHMC poles. For the lightest mass-preconditioned field and RHMC poles on the L/a = 16, 24 lattices, we employ the SAP-preconditioned GCR algorithm [31].

Results
On most of our ensembles, we measure the correlation functions defined in Eqs. (2.5) and (2.7) on every fourth trajectory of length τ = 2 MDU so that the spacing between these measurements is 8 MDU. Only on A1k2 and A2k1, we use a measurement separation of 2τ = 4 MDU. The total statistics for all the ensembles considered here are tabulated in Table 1.
In addition to these correlation functions, we also measure 'smoothed' gauge field observables obtained from the Wilson (gradient) flow [20,21], which possess a well-defined continuum limit. These smoothed observables are useful in several ways. The smoothed gauge fields provide a renormalized definition of the topological charge, which we use to monitor the topology freezing discussed in Section 2. Even at lattice spacings where topology freezing is not a problem, the smoothed topological charge and action typically possess the largest observed autocorrelation times. Furthermore, the aforementioned renormalized (and L-dependent) coupling ḡ 2 GF of Ref. [25] is defined using the Wilson flow and may be used to monitor the deviation from the constant physics condition, as it is sensitive to the physical lattice size. Results for this coupling are given in Table 2, too.
In order to monitor the autocorrelation times in our simulations, we examine these smoothed observables at a flow time t given by √ 8t/L = c with c = 0.35. For all simulations we find that integrated autocorrelation times of these observables satisfy the bound τ max 200-250 MDU, except for our L/a = 24 simulations where the charge is frozen. The other smoothed observables turn out to still possess autocorrelation times of comparable order of magnitude in these simulations. The existence of the charge freezing at our finest lattice spacing is qualitatively Fig. 1. Histories of the smoothed Wilson plaquette action and the topological charge for single representative replica from each of the ensembles, which enter into the final analysis. Our inability to sufficiently sample all topological sectors at β = 3.81 is evident. illustrated in Fig. 1. There it is seen that the autocorrelation time of the smoothed action remains under control, whereas the autocorrelation time of the topological charge increases significantly from the L/a = 20 to the L/a = 24 ensembles. We are practically unable to sufficiently sample all topological sectors at this finest lattice spacing, necessitating the restriction of our observables to the trivial sector.
To provide further support for this projection, we have estimated the expectation value δ Q,0 ; it effectively corresponds to the fraction of statistics, to which the restriction to Q = 0 reduces the number of generated configurations. For the ensembles with the smallest quark mass at given L/a, {A1k1, B1k3, C1k3, D1k1}, we find δ Q,0 = 0.37(2) (L/a = 12), 0.44(2) (L/a = 16), 0.65(7) (L/a = 20), 0.90(5) (L/a = 24), and thereby to show large cutoff effects. These kind of cutoff effects have been observed before in the large-volume, two-flavor theory in Ref. [32], owing to the substantial suppression of the topological charge compared to the quenched approximation. In addition, our δ Q,0 at m PCAC ≈ 0 does not seem to smoothly (in the sense of monotonically and linearly or quadratically in a) approach 1 as L/a increases, which is the theoretical expectation for δ Q,0 in the large-volume continuum limit and at zero quark mass [33]. Even though we are not in a large-volume situation with m PCAC = 0 exactly, we interpret the encountered behavior of δ Q,0 as a sampling problem of the algorithm, on top of the cutoff effects. Therefore, we prefer to project our results to the Q = 0 sector, which in the end does not induce a noticeable difference in the final numbers for c A , cf. Table 2.
The estimation of the statistical error on our measured quantities is based on a singleelimination jackknife procedure, after first 'binning' the data (concatenated from different replica) such that the size of each bin is τ max , as well as on applying a full autocorrelation analysis according to Refs. [34,35]. Both yield similar error estimates; our quoted final results in the Q = 0 sector stem from the latter (without including any long-tail contributions to the autocorrelation functions, which are negligible for the quantities entering the c A -analysis).
After measuring the correlation functions, we solve for the largest two eigenvectors of the matrix given by f 1 (ω i , ω j ). These normalized eigenvectors have a well-defined continuum limit along our line of constant physics in parameter space, as long as the wavefunctions depend on physical scales only. In fact, as we do not observe any significant lattice spacing dependence for them, we fix these eigenvectors for once to the values calculated on the B1k2 ensemble (L/a = 16, β = 3.512, κ = 0.13703) and regard them as part of our choice of improvement condition. For our setup we find: η (0) = (0.5317(3), 0.5977(1), 0.6000(2)) and η (1) = (0.843 (5), −0.31(6), −0.44(6)), which are similar to those of Refs. [6,7].
To get an idea of the sensitivity of our method to c A , we examine the effective masses of the correlation function f P , after taking the inner product with the eigenvectors for the (approximate) ground and first excited states. These are shown in Fig. 2 for the L/a = 16 ensemble B1k2 of the previous paragraph. The distinctly seen signals display that indeed the eigenvectors effectively maximize the overlap with the ground and first excited states, since these states are clearly separated up to x 0 ≈ 11a. As already noted in Ref. [6], the energy of the first excited state is somewhat near the cutoff a −1 though, and this may influence the way in which residual cutoff effects are modified. Hence, residual O(a 2 ) effects may grow rapidly in smaller volumes [17], which justifies our choice of an intermediate volume with L ≈ 1.2 fm to impose the improvement condition at constant physics. Also shown in Fig. 2 are r(x 0 ) and s(x 0 ) (the latter being proportional to m 2 π (1) − m 2 π (0) in case of exact ground and excited state projections) from Eq. (2.11) for the same ensemble, further demonstrating our good sensitivity to c A .
Finally, the function c A (x 0 ) is shown in Fig. 3 for the ensembles used in the analysis. It indicates that c A is rather independent of x 0 for points sufficiently distant from the boundary. Moreover, only a rather little variation of c A (x 0 ) is visible for x 0 5a; this reveals that highenergy states, which could induce large O(a) ambiguities in the improvement condition, are reasonably suppressed in this region. As a compromise between cutoff effects and statistical errors, we take x 0 = L/2 as our final definition for c A , which besides is well within the regime where states with a distinct energy gap dominate the projected correlators. These data are plotted in Fig. 5, together with an interpolation.
Before discussing these final results for c A at each bare coupling, we assess our systematic errors. In order to estimate the effect of our finite (but small) quark masses and thus of small violations of the constant quark mass condition on c A , we study several different sets of ensembles, which are identical apart from |Lm PCAC | < 0.3 such as {A1k1, A1k2}, {B1k1, B1k2, B1k3}, and {C1k1, C1k2, C1k3}. Within each of these sets, the variation of the value of c A does not exceed more than about 1.5 standard deviations. To quantify the cutoff effect, which  results from T = 3L/2 − a, we compare {A1k1, A1k2} with A2k1, where the temporal extent is T = 3L/2 + a. We see that this results in a difference which is significant at less than 1σ at our coarsest lattice spacing where it is largest. Lastly, we assess the error due to the deviation   Table 2 one infers that the ensembles {A1k1, B1k3, C1k3, D1k1}, which enter into the final analysis, have ḡ 2 GF,0 (L) between ≈ 14-18, resulting in a ≈ 20% variation assuming no cutoff effects. We can explicitly test the sensitivity of c A to this variation by means of the B2k1 ensemble, which differs with respect to the B1 ensembles in a by ≈ 6% and in ḡ 2 GF,0 (L) by roughly the same amount as this overall ḡ 2 GF,0 (L)-variation. The value of c A determined on the B2k1 ensemble lies well within 1.5σ of the interpolating formula in Fig. 5, so we are confident that this variation of L does not result in a significant shift of c A . Note that, after all, any imperfection in the constant physics condition of this level and the systematic uncertainty in c A induced by it will only introduce O(a 2 ) effects in quantities involving the improved axial current, which are negligible compared to other sources of errors and vanish in the continuum limit by definition.
We now move to the final results. In order to guide our choice for an interpolation formula, we observe that (for our choice of the improvement condition) c A is almost linear in a over the range of bare couplings which we simulate. This approximate linearity is depicted in Fig. 4. Notice that the linear behavior can not extend all the way to a = 0, since this would be incompatible -owing to the non-polynomial relation between a and g 2 0 -with a polynomial dependence of c A on g 2 0 in the perturbative regime. In our case, a naive linear fit to these data within the simulated region does not even extrapolate to zero, which is the value predicted by perturbation theory at tree-level: c A (g 2 0 = 0) = 0 [36,37]. On the contrary, we interpret the behavior of the results according to c A = constant + slope × a within the region of our data such that the (non-vanishing) constant term removes the targeted O(a) effects in the non-perturbative regime, while the non-constant piece, describing the non-trivial dependence of c A on g 2 0 , only affects O(a 2 ) contributions in physical quantities. Regardless of their actual size, these intrinsic O(a) ambiguities suggest to employ in the computation of physical quantities the g 2 0 -dependence of c A induced by the constant physics condition rather than, e.g., the constant piece alone, because one then expects the remaining leading O(a 2 ) lattice artifacts to be smaller and vanish uniformly in the O(a) improved theory.
Motivated by the apparent linearity of c A as a function of the lattice spacing within our data region, we choose the following interpolation formula for c A as a function of g 2 0 = 6/β: which is constrained to reproduce the 1-loop perturbative result of Ref. [12] as g 2 0 goes to zero. The piece depending exponentially on the bare coupling is inspired by the perturbative expression of the lattice spacing in terms of g 2 0 , and it is inserted in order to capture the linear behavior reflected in Fig. 4. Our non-perturbative results for c A are then compiled in This formula can be used at all bare couplings below g 2 0 ≈ 1.8, together with the statistical errors on c A (i.e., on c A,0 as given in Table 2), which are ≈ 4-5% near the largest simulated bare couplings and increase to ≈ 7-8% near the smallest.

Conclusions
In this work we have determined non-perturbatively c A (g 2 0 ), the coefficient required for O(a) improvement of axial current matrix elements in lattice QCD with N f = 3 flavors of Wilson quarks, non-perturbative c sw [11] and the tree-level Symanzik-improved gauge action.
The main result is the interpolation formula in Eqs. (4.1) and (4.2), obtained using the standard improvement condition, which is imposed together with a variation of boundary wavefunctions in the PCAC relation with external states and along a line of constant physics in parameter space. This implies that potentially large O(a) ambiguities in c A are avoided and that its remaining intrinsic O(a) ambiguities disappear smoothly toward the continuum limit. Eq. (4.1) for c A (g 2 0 ) (with coefficients (4.2)) is valid for bare couplings below g 2 0 ≈ 1.8 (or, equivalently, for lattice spacings a 0.09 fm). We have treated the topology freezing, encountered in our simulations at the finest lattice spacing (of around 0.045 fm), by restricting the improvement condition to the trivial topological sector for all ensembles considered.
The gauge field ensembles entering this work were generated with T = 3L/2, in order to be re-used for the determination of the axial current renormalization constant Z A , see [18] for a preliminary report. After this is completed (the results of which will appear in a future publication), axial current matrix elements such as pseudoscalar meson decay constants can be calculated precisely from large-volume ensembles of gauge configurations at lattice spacings typically employed in the context of phenomenological applications of lattice QCD.