NSPT estimate of the improvement coefficient $c_A$ to two loops

By using Numerical Stochastic Perturbation Theory (NSPT), we carry out a quenched two-loop computation of the improvement coefficient $c_A$ associated to the isovector axial current. Within the Schr\"odinger Functional formalism, we compute the bare quark mass $m$ and fix $c_A$ by requiring discretization corrections on $m$ to be of order $O(a^2)$ in the lattice spacing $a$.

Introduction -After performing a set of numerical simulations of a generic field theory, corresponding results become meaningful only if extrapolated to the physical point (Φ). This means that, given simulation parameters like -among others -lattice spacing a, volume V and, in case, quark masses m q , the limits a → 0, V → +∞ and m q → m (Φ) q have to be accurately and reliably computed. While these limits were thought to be hardly accessible for lattice QCD some time ago [1], a few decades of algorithmic developments have paved the way to controlling all the above-mentioned sources of systematic uncertainty. Nowadays, not only N f = 1 + 1 + 1 + 1 QCD simulations with small lattice spacings, large boxes and quark masses close to the physical point are being performed, but also challenging QED effects have started to be taken into account [2].
At present, the standard setup for QCD simulations features the Hybrid Monte Carlo (HMC) algorithm [3], usually supplemented with other techniques -like preconditioning [4], the Hasenbusch trick [5], multiple timescale integration [6] and smearing [7][8][9][10]. In the framework of the HMC algorithm, such techniques are useful in that they allow for a larger value of the time-step in integrating the equations of motion (thus decreasing the autocorrelation among subsequent configurations), they reduce the condition number κ(D) of the Dirac operator D (speeding up the inversion of D needed both within the HMC algorithm and at measure time) and they average out ultraviolet fluctuations, improving the signal-to-noise ratio (SNR).
In this scenario, a technique often employed to reduce discretization effects is the so-called improvement, usually in a form following the Symanzik programme [11]. To understand how the latter works, it is useful to recall that the mean value O of a generic continuum observable O computed on the lattice reads where O Lat is the lattice counterpart of O, S Lat is the lattice QCD action (given by the sum of gauge part S G and a fermionic part S F ) and where the dependence of O with respect to simulation parameters has been made explicit. According to Symanzik, close to the continuum limit, O Lat can be Taylor-expanded in the lattice spacing as where O 1 , O 2 , . . . have to be interpreted as contributions stemming from operator insertions in the continuum and must have symmetry properties consistent with O. A similar expression holds also for the gauge and the fermionic action S G and S F entering S Lat . Plugging said Taylor expansions with respect to a into Eq.(1) results in a similar expansion for O as well. In the Symanzik improvement programme, the leading correction to O -usually linear in a as in Eq.(2) -can be cancelled by adding irrelevant terms to S G , S F and to O Lat 1 . In this way, the dependence of O with respect to a is flattened and, consequently, larger values of the lattice spacing can be used to recover the continuum limit, thereby further reducing k(D) and increasing the SNR at the same time.
In general, each irrelevant term is multiplied by its own coefficient that has to be appropriately tuned: its value can be determined either non-perturbatively or within perturbation theory (PT). Obviously, among these socalled improvement coefficients, the most important ones are those improving on S G and S F since they enter in the improvement procedure of any observable. When it comes to perturbative computations, the expansions of these coefficients in the bare coupling g 0 are usually known up to very low orders only, usually at one loop.
In this paper, we investigate whether Numerical Stochastic Perturbation Theory (NSPT) [13,14] can be applied to compute the improvements coefficients in PT to orders higher than g 2 0 . We want to clearly state that the present study essentially aims at being a proof of concept, i.e., at assessing the feasibility -or not -of a similar NSPT computation in principle. Given such exploratory character, we tackle the simplest case possible, namely the two-loop computation of the improvement coefficient c A associated to the isovector axial current in the quenched approximation. In spite of its seeming minor importance, a perturbative computation of c A to a given order j allows for a determination to the same order of the much more important coefficient c SW , i.e., the improvement coefficient multiplying the irrelevant term improving on the fermionic action S F [19] which will be introduced later.
Lattice setup -In this and the next section we outline our approach which is based on the Schrödinger Func-tional formalism [15] and closely follows the strategy described in [16] and used in [19] to compute c SW and c A to one loop.
We simulate a four-dimensional lattice made up of N T × N 3 S sites, each one labelled with integer coordinates n = (n 0 , n 1 , n 2 , n 3 ) varying in the intervals [0, N T −1] and [0, N S − 1] along the time and spatial directions respectively. The lattice volume V will therefore be equal to V = L T × L 3 S , with L T = aN T and L S = aN S . A generic gauge variable U µ (n) -with µ ∈ {0, 1, 2, 3} -belongs to the SU (3) group and is associated to the link connecting site n to site n +μ,μ being a unit vector along direction µ. Lattice group variables are related to their continuum counterparts A µ (n) in the Lie algebra of SU (3) through the equation U µ (n) = exp[iA µ (n)]. We stick to the usual convention according to which U −1 µ (n) = U † µ (n). Quark and antiquark degrees of freedom are Grassmann variables -denoted as ψ(n) andψ(n) respectively -associated to the lattice sites 2 . For simplicity, we assume the presence of N f mass-degenerate flavours though, in what follows, the fermionic action S F and its irrelevant term will always be written by taking into account only one flavour to ease the notation: it is understood that there are actually N f replicae of such operators. Anyway, it is worth stressing that, while the quenched approximation implies -obviously -that fermionic degrees of freedom play no role in updating the lattice configuration and that, consequently, the results of this paper have to be considered valid for N f = 0, the setup described in the next sections is such that N f never enters into play at measurement time as well, as long as the mass degeneracy holds.
While boundary conditions are periodic along the three spatial directions, they are of Dirichlet type in the time direction. In other words, by labelling a generic spatial direction with k from now on, for the gauge fields the following equalities hold with n = (n 1 , n 2 , n 3 ) and where W k ( n) can be expressed in terms of a smooth, fixed field C k ( n) as being P the path-ordering symbol. W k ( n) is parametrized by another field C k ( n) in an analogous way.
With this setup, the lattice gauge action S G is given by the modified Wilson action where U (p) is the product of the link variables around a lattice plaquette, the sum runs on all oriented plaquettes and the weights ω(p) are equal to 1 for each plaquette, except for the spatial ones at n 0 = 0 and n 0 = N T − 1 where ω(p) = 1 2 . Note that S G is already O(a)−improved with respect to its continuum formulation: this can be shown by Taylor-expanding Eq.(5) in a. The irrelevant terms cancelling leading corrections proportional to a 2 in S G are described in [17] but they are not to be needed in this study, which aims at reaching O(a)−improvement only.
The unimproved fermionic action S F is given by m 0 being the bare quark mass and the Wilson-Dirac operator D reading where repeated indices are summed and γ µ 's are Euclidean Dirac matrices. Covariant derivatives in Eq. (9) are defined as λ 0 = 1 and λ k = exp(iaθ k /L S ) being phase factors (with −π < θ k ≤ π). For simplicity, all three angles θ k will be set to the same unique value θ = 0. Strictly speaking, Eq.(8) holds in an infinite volume. However, it remains valid also in the present setup -i.e., a box of finite size with Dirichlet boundary conditionsprovided some technical conventions are assumed: the interested reader can find more details in subsection 4.2 of [16]. We tacitly take such conventions for granted and carry on with Eq.(8) in combination with said lattice topology.
The leading discretization correction to S F is linear in a and, as mentioned in the introduction, an O(a)−improved fermionic action S imp F can be obtained by adding to Eq.(8) an irrelevant term δS V , i.e., δS V is usually referred to as clover term and its expression is given by [12] δS V = a 5 c SW where with The perturbative expansion of the c SW coefficient appearing in Eq. (12) is known up to one loop and it can be written as where c SW has been computed in several papers [18][19][20][21] yielding results slightly different but in agreement within errorbars.
Before concluding this section, it is important to stress that the bare quark mass m 0 in Eq.(8) will be set to 0 from now on and that quarks will be kept massless by subtracting the appropriate mass counterterms order by order in PT (see [22] for their calculation in infinite volume to two loops).
Methodology -The improvement coefficient c A targeted by this study is associated to the isovector axial current τ b being a Pauli matrix acting on flavour indices and γ 5 = γ 0 γ 1 γ 2 γ 3 as usual. An O(a)−improved expression is obtained by adding an irrelevant term δA b µ (n) reading where δ * µ and δ µ stand for the standard left and right derivative on the lattice while P b (n) is the isovector axial density As in Eq. (15), the improvement coefficient c A in Eq. (17) can be expanded as where c A has been determined in [19] while estimating c (2) A is the goal of this work. Following [16], we begin by relating A b µ (n) and P b (n) to the unrenormalized quark mass m by means of the PCAC relation with O the product of fields located at non-zero distance from site n and from each other. Then, we set O to be with the constraint n 0 = n 0 = 0 and with being ρ( n) andρ( n) the fields introduced in Eqs. (6)(7).
After defining the correlators f A (n) and f P (n) as again with the constraint n 0 = n 0 = 0, the unimproved quark mass m in Eq.(20) is given by Noting that, just like S G , both P b (n) and the operator O in Eq. (21) are already O(a)−improved [16] and recalling Eq. (17), the O(a)−improved quark mass m imp is given by provided that the irrelevant term in Eq. (12) is also added to S F and that both c A and c SW are correctly set. The last observation gives us a prescription to determine the improvement coefficients. Before explaining why, it is worth recalling that any quantity computed on the lattice is dimensionless until the lattice spacing is determined through scale setting. However, the NSPT practice detailed in the next section is such that any reference to a is purely formal in this study: consequently, any observable being measured in the NSPT framework eventually remains dimensionless.
Bearing this in mind and going back to Eq.(25), in PT m imp can be expanded as where coefficients m (i) imp will depend on the coefficients c A and c SW and on the kinematic parameters a, N S , N T , θ as well as on the time coordinate n 0 of site n in Eq.(25) -there is no dependence with respect to the spatial coordinates of n because of the translational invariance along the corresponding directions. m imp should also depend on the boundary fields C k , C k , ρ,ρ, ρ ,ρ : however, fermionic boundary fields will be set to zero after derivatives in Eq. (21) are computed while fields C k ( n) and C k ( n) will be fixed to 0 for every n (so that W k ( n) = W k ( n) = 1). The last choice will be motivated later on.
Since infinite-volume mass counterterms are subtracted up to two loops and since every quantity does not carry any dimension, both m (i) imp -with i = 1, 2can be expanded in a out of dimensional analysis as where dots denote terms of higher order in a and all coefficients d (i) 1,... depend on c A and c SW -such dependence will be left implicit to ease the notation. By setting N T = 2N S + 1, n 0 = N S /2 (as in [19]) and by keeping θ fixed, the resulting mathematical setup is such that the terms on the r.h.s. of the previous formula can be collected into one as By comparing Eqs. (27) and (28), it should be evident that an expansion in powers of a is equivalent to an expansion in powers of N S . Bearing this observation in mind and recalling that we aim at O(a)-improvement, coefficients c A in Eqs. (15) and (19) can be determined up to two loops as follows: with the setup outlined above, the PCAC quark mass is first measured for several values of N S , then it is fitted vs. 1/N S and c SW and c A are finally determined by requiring the coefficient d In this approach, there is actually one last issue to be solved: in fact, to a given loop i in PT, the coefficient d SW with j ≤ i, so that their effects have to be disentangled. This can be done by choosing the boundary fields C k and C k appropriately. In particular, if such fields are both set to 0 everywhere along the time boundaries as in the present setup, it can be proven [16] that the dynamical gauge degrees of freedom U are 1 throughout the whole lattice and, consequently, Eqs. (13) and (14) imply 3 that the lowest order of F µν in Eq.(12) will be proportional to g 0 . In turn, this means that, truncating any expansion in g 0 at a given loop i, only c (i) A -as well as all coefficients at loops lower than i in Eqs. (15) and (19) -will be left into play. This yields to a well-defined procedure to evaluate c have already been determined for i ≤ 1, by setting these tree-level and one-loop coefficients to their known values as well as the fields C k and C k to 0 and by truncating any perturbative expansion at the second loop in g 0 , d 1 will only depend on c (2) A and the latter coefficient can thus be fixed by fitting m (2) imp with respect to 1/N S . Though it is not the goal of this study, let us recall how c (2) SW could be evaluated. The very same setup needed to compute c (2) A is maintained but the fields C k and C k have now to be set as explained in Sect. 6.2 of [16]: fixing c (2) A to the value found as outlined in the previous paragraphs, d 1 will now depend solely on c (2) SW , so that the correct value of the latter coefficient could be determined, again by fitting m (2) imp vs. 1/N S . This overall procedure can obviously be iterated to the third loop (and higher), provided that the corresponding mass counterterm is subtracted.
Before concluding this section, it is worth recalling that, within the Schrödinger Functional formalism, also boundary irrelevant terms in a have in principle to be introduced in order to achieve O(a)−improvement, each one with its own coefficient. However, as stated in [19], these terms can be eventually dropped and remaining improvement coefficients can be determined by solely requiring the unrenormalized quark mass to be independent of the kinematic parameters, which corresponds to the strategy outlined above.
NSPT practice -In this section we describe how configurations are generated by means of Numerical Stochastic Perturbation Theory. NSPT stems from Stochastic Quantization (SQ) [24], an algorithm that allows for the computation of expectation values in quantum field theories. Being convenient from the point of view of computer simulations in that it does not require an accept/reject step (thus allowing for a fast updating process), SQ has been used in several domains of research, one of the latest being the search for solutions to the sign problemsee [25,26] and references therein. Its main drawback is the need to perform simulations for different time steps in order to extrapolate to continuous stochastic time, as explained in what follows.
To introduce the basics of SQ in a simple way, we start with a lattice scalar field theory with action S[φ]. In SQ its degrees of freedom φ(n) are updated by numerically integrating a Langevin equation reading where t is the so-called stochastic time and η(n, t) is a Gaussian noise satisfying η(n, t) η = 0 , η(n, t)η(n , t ) η = 2δ(n − n )δ(t − t ) .
The subscript "η" stands for an average over the noise. Given a generic observable O(φ), it can be shown [27] that the time averagē is equal to the path-integral mean value, i.e., Z being the partition function. After discretizing the stochastic time t as t = m (with integer m), Eq.(29) can be numerically integrated through the prescription 4 φ(n, m + 1) = φ(n, m) − f (n, m) , where the force term f (n, m) in the Euler scheme is given by with η(n, m) = √ η(n, t = m ). Since the equivalence in Eq.(31) holds only for continuous t, computer simulations with different values of have to be carried out to extrapolate to → 0, as already mentioned. The SQ setup for the scalar theory has to be modified in order to be applied to SU (3) link variables. In this respect, Eq.(29) is modified as where matrices T a are the generators of the SU (3) algebra (with normalization tr(T a T b ) = 1 2 δ ab ) and ∇ a n,µ is the Lie derivative -with respect to the algebra fields associated to variable U µ (n) -defined as [28] f being a scalar function of the group variable U and ω a 's small parameters. As for the group counterpart of Eq.(32), it reads U µ (n, m + 1) = e −i a T a f a µ (n,m) U µ (n, m) , with f a µ (n, m) = ∇ a n,µ S G [U ] − √ η a µ (n, m) .
In this framework, PT up to N L loops can be introduced by a formal expansion of each gauge field as with β = 6/g 2 0 . A few remarks are in order concerning Eqs. (38). First, the leading order of U µ (n) is 1 in the light of the previous choice of the fields C k and C k (see the comments at the end of the previous section). Second, while the fields A (i) µ (n) are elements of the Lie Algebra of SU (3), the fields U (i) µ (n)taken one by one -do not belong to the SU (3) group but U µ (n) on the l.h.s. of the second of Eqs. (38) does, up to terms of order O(β −N L +1 ) excluded. Finally, a Taylor expansion of the exponential of A µ (n) and of the logarithm of U µ (n) allows for switching between algebra and group.
Plugging the second of Eqs. (38) into a discretized version of Eq.(34) results in a hierarchical system of differential equations where the evolution of a given order U (i) µ (n) only depends on lower orders, thus allowing for a consistent truncation at the needed loop N L . In this setup, the noise η a µ (n) enters at order g 0 after a further rescaling of , as explained in [37]. This is the core of the NSPT algorithm which has been applied to lattice QCD in order to study -among others -the free energy density at finite temperature [29][30][31], renormalization constants [32][33][34][35] and renormalons [13,[37][38][39][40][41][42].
In order to prevent fluctuations associated with Gribov copies [43], the updating step in Eq. (36) has to be alternated with the so-called stochastic gauge fixing [44]. In other words, before moving from configuration U µ (n, m) to U µ (n, m + 1) in stochastic time, each field U µ (n, m) has to undergo a gauge transformation like U µ (n, m) → G(n, m)U µ (n, m)G † (n +μ, m) , where matrices G belong to the SU (3) group as well. With periodic boundary conditions along all directions, a common choice for G(n) is G(n) = exp[Ω(n)] where Ω(n) reads [45] which results in fluctuations around the Landau gauge.
With Dirichlet boundaries, this definition of the entries of matrix Ω(n) is modified for sites with n 0 = 0 or n 0 = N T − 1 as follows: where the constraint n 0 = 0 holds in the sum over n .
To carry out the computation of correlators f A (n) and f P (n) entering in the PCAC quark mass, it is necessary to contract the fermionic fields in Eqs. (23) and, therefore, to invert the Wilson-Dirac operator D in Eq. (9). This is done following [14]: given a generic operator M and its perturbative expansion M = k=0 g k 0 M (k) , its inverse M −1 can be expanded perturbatively as where the tree-level term is the inverse of the zero-order term of M while, recursively, The expression for D (0) −1 can be found in Sect. 3.1 of [19].
Finally, it is worth mentioning that, even if the bare quark mass m 0 in Eq. (8) has been set to 0, the condition θ = 0 implies a non-vanishing quark mass at tree level [19]: thus, in this setup there are no zero modes to take care of.
Data analysis and results -Simulations were performed with the following values of the parameters: θ = 1.2, ∈ {0.005, 0.007, 0.010, 0.015, 0.020} and N S ∈ {11, 12, 13, 14, 15, 16, 17, 18, 20, 24, 32}. Note that, since the Euler scheme is employed in the integration of the Langevin equation, only three values of would be actually needed to extrapolate to → 0. However, we used 5 time steps in order to increase the precision of the extrapolation.
We did not carry out simulations with N S lower than 11 since corresponding results turned out to be rather noisy. Our understanding is that, with such values of N S , low modes -which are not zero but should still be relatively small -occupy a non-negligible fraction of the configuration space and eventually spoil the SNR. Tab.(I) details the number of measurements of observables f A and f P in Eqs.(23) -and, hence, of the unrenormalized quark mass -for the different combinations of simulation parameters (N s , ). For each setup two subsequent measurements were separated by 100 Langevin updates of the lattice configuration in order to reduce the autocorrelation. It is worth stressing that these 100 updating steps usually took approximately half of the time needed to perform a single measurement of the quark mass. Fig.(1) shows a plot of the unimproved quark mass m in Eq.(24) at the first (m (1) ) and second loop (m (2) ) vs.
together with the extrapolated result. Statistical errors on the data in blue in Fig.(1) have been obtained through a jackknife procedure.
Before computing c A , we check to which extent our approach is reliable: it is understood that all quark masses referred to in the rest of this section are obtained from an extrapolation to → 0 as that shown in Fig.(1 As a first test, we compared the analytical result for the tree-level quark mass [19] with our estimate. The two values agree apart from round-off errors that are usually lower than 0.03% -the smallest tree-level quark mass we measure is equal to 2.97 · 10 −6 for N S = 32, which is also the extent where round-off errors turn out to be the largest. As a second check, we now work out the one-loop 5 The only exception to this assumption is given by tree-level quantities since they do not depend on the Langevin time . This is due to the r.h.s. of Eq.(37), whose leading order is g 0 : in fact, in the present setup, the action derivative is zero at tree level while, as already mentioned before, the Gaussian noise enters at order g 0 in PT. Consequently, the tree level of any variable Uµ(n) does not evolve with respect to the stochastic time and bears no statistical error. coefficient c (1) A and compare it to the value determined in [19]. To better explain the fit strategy, we rewrite the improved quark mass m imp in Eq. (25) as where m is defined in Eq.(24) and a has been set to 1 -let us recall once more that any dependence in the lattice spacing is purely formal. Plugging Eq. (19) into the previous formula, at one-loop level we have δm (0) corresponds to the tree-level contribution to the term (including the denominator) multiplied times c A on the r.h.s. of Eq.(25): being a quantity at tree level, it bears no error. In order to determine c A , we proceed as follows: tentative valuesc (1) A 's are assigned to c (1) imp is fitted vs. f (1/N S ) = d 1 N S and anyc (1) A is eventually retained as valid if the coefficient d is compatible with 0 within errorbars. In this way, we obtain a range [c (1) A,min ,c As for the systematic error, its main source is given by the truncation of function f (1/N S ) at the first order in 1/N S . Therefore, the best way to assess the systematics would be to repeat the procedure above with a function of higher order in 1/N S and to compare the outcome with that obtained with a linear function. In fact, as pointed out in [19], for relatively large values of θ as that used in this study, corrections going like 1/N 2 S can be comparable to the leading contribution 1/N S . Unfortunately, our data are not precise enough to support higher-order terms in f (1/N S ): any attempt of fit in this sense, either with logarithmic or power corrections, winds up with an extremely poor determination of c (1) A . In order to assess the systematic error in a somehow coarser way, we carried out the fit of m (1) imp vs. f (1/N S ) = d 1 N S as before but within a range of N S limited to the three smallest extents (i.e., N S ∈ [11,13]). Since the latter is the regime where higher-order corrections in N S should have more impact, the mean value of c (1) A obtained in this way should feature the largest deviation with respect to the mean value of c (1) A computed with the fit employing all available sizes. Such a deviation could then be considered as a rough estimate of the systematic error.
Following this strategy and converting our expansion in β −1/2 into a series in g 0 , the result we obtain for c where the first and second error are the statistical and systematic uncertainty respectively. Within errorbars, this values is in reasonable agreement with c A = −0.00756(1) quoted in [19], though the precision of the latter estimate is much higher than that obtained with NSPT.
If we now move to the second loop and to the evaluation of c where c (1) A has been set to its known value c (1) A = −0.00756. Starting from the last expression, the same procedure adopted at one-loop level can be applied to the fit of m (2) imp with respect to 1/N S , only m (1) in Eq.(43) has to be replaced withm (2) (1) . Similarly to what remarked before for δm (0) , δm (1) corresponds to the one-loop contribution to the term (including the denominator) multiplied times c A on the r.h.s. of Eq. (25).
The result we get for c (2) A is While at one loop the systematic error is slightly larger than the statistical one, at two-loop level systematic effects seem to be apparently less important. This is most likely a consequence of a worse signal-to-noise ratio at the second loop, as shown in Fig.(2) where m (1) andm (2) are plotted -together with their O(a)−improved counterparts -on the left and right panel respectively. At one loop, statistical errors are smaller and, consequently, the trend of the data is altogether better defined, thus allowing for the rough assessment of subleading corrections in 1/N S outlined above. On the contrary, at two-loop level the trend of the data is harder to be determined due to the larger errorbars: this eventually leads to a more difficult estimate of the impact of subleading contributions in 1/N S and, in practice, to an apparently small systematics.
Conclusions -By studying the cutoff dependence of the unrenormalized PCAC quark mass, we carried out an exploratory, perturbative determination of the improvement coefficient c A to two loops by using NSPT in the quenched approximation. The crosschecks at zero-and one-loop in PT reproduce -with fair precision -results previously obtained in other independent studies, thus indicating that NSPT can in principle be used in this kind of computations, even at orders where non-numerical approaches are usually too cumbersome from an algebraic point of view. Let us recall that the order N L at which the Taylor series in Eqs. (38) are truncated is in fact just a parameter and increasing it does not entail any extra work algebra-wise.
However, though viable in theory, at a practical level NSPT seems to be very demanding in terms of computer resources when it comes to producing an accurate estimate of the improvement coefficients, even in a relatively simple setup as that examined in this study, featuring the quenched approximation and the Wilson action at two loops. Indeed, in spite of having employed some millions of CPU hours, data are still rather noisy, as shown in Fig.(2): though not optimal at the first loop already (and, in fact, our estimate for c (1) A is much less precise than that contained in [19]), the deterioration of the signal at higher order in PT is manifest and gets particularly reflected in a difficult assessment of the systematic uncertainty. It might be that observables other than that examined in this study -i.e., the unrenormalized quark mass -display an intrinsically better signal-to-noise ratio but, at present, we would have no suggestion to make with respect to this issue.
Our understanding is that reducing errorbars to a level at which NSPT data can be accurately fitted (and improvement coefficients precisely determined) would most likely require considerable computer resources irrespective of the quantity being monitored, especially if higher order in PT are envisaged and/or more interesting -but also more demanding -setups are considered, as that of unquenched simulations.
More or less recently, some works have actually been published [46][47][48] proposing some technical changes to the standard NSPT approach followed in this study. Such changes yield new setups -called Instantaneous Stochastic Perturbation Theory (ISPT), Hybrid Stochas- A on the right have been set to the mean values reported in Eqs. (45) and (47) respectively.
tic Perturbation Theory (HSPT) and Kramers Stochastic Perturbation Theory (KSPT) -that seem to be highly beneficial in remeding the shortcomings of the standard version of NSPT we employed: this looks particularly true for the formulation described in [48]. It would be extremely interesting to carry out a computation of c (2) A with the latter setup in order to (hopefully) get a more precise estimate of this improvement coefficient at much lower computer costs 6 .
Nevertheless, we believe that the original result of this project, i.e., the evaluation of the coefficient c (2) A for N f = 0, can serve as a benchmark to some extent: indeed, we are confident that at least its negative sign and its order of magnitude (that is, 10 −3 ) are reliable.