The θ -dependence of the Yang–Mills spectrum from analytic continuation

: We study the θ -dependence of the string tension and of the lightest glueball mass in four-dimensional SU( N ) Yang–Mills theories. More precisely, we focus on the co-efficients parametrizing the O ( θ 2 ) dependence of these quantities, which we investigate by means of numerical simulations of the lattice-discretized theory, carried out using imaginary values of the θ parameter. Topological freezing at large N is avoided using the Parallel Tempering on Boundary Conditions algorithm. We provide controlled continuum extrap-olations of such coefficients in the N = 3 case, and we report the results obtained on two fairly fine lattice spacings for N = 6


Introduction
One of the most interesting features that emerges when studying the non-perturbative regime of quantum field theories (QFTs) is their θ-dependence: peculiar terms exist (the so called θ-terms) which, when added to the action, do not modify the classical equations of motion, and yet change the physical properties of the theory.The existence of θ-terms is related to the topological features of the gauge group and to the space-time dimensionality [1,2], so θ-dependence is not present in all QFTs.Nevertheless, several interesting QFTs display a nontrivial θ-dependence, ranging from Quantum Chromo-Dynamics (QCD) and four-dimensional SU(N ) Yang-Mills theories [1,3,4], to two-dimensional models like the CP N −1 models [5,6] and U(N ) Yang-Mills theories [7][8][9] (and even elementary quantum mechanical models [2,10,11]).
In this work we investigate an aspect of θ-dependence that has received far less attention: the θ-dependence of the spectrum of the theory.In QCD, close to the chiral limit, it is easy to derive the θ-dependence of the mass of the pseudo-Nambu-Goldstone bosons associated to the spontaneous breaking of chiral symmetry [13]; away from the chiral limit, we have once again to resort to lattice simulations.This is the case also for four-dimensional pure-gauge theories, which are however much simpler to simulate than QCD.For this reason here we focus on the case of four-dimensional SU(N ) Yang-Mills theories, whose Euclidean Lagrangian density is given by where and we investigate the θ-dependence of the string tension σ and of the lightest glueball mass m G .Analytical computations can be performed in two-dimensional U(N ) Yang-Mills models, which however do not seem to provide much insight on the physics of their fourdimensional counterparts, since there is no θ-dependence at all in the spectrum of these two-dimensional models in the continuum (see App. A).At θ = 0, it is established that the 0 ++ glueball ground state represents the lightest state [51][52][53][54].At θ ̸ = 0, spatial parity is explicitly broken, and cannot be used as a quantum number for glueball states.The latter are thus only characterized by their spin and charge conjugation quantum numbers.For this reason, we denote by m G the mass of the lightest glueball state, i.e., the one that tends to the 0 ++ glueball in the θ → 0 limit.Note that, since in our study we only investigate the small-θ regime, we can a priori exclude the possibility of a level crossing between different states.
Using the invariance under parity of the θ = 0 theory, we can parameterize the leading order θ-dependence of the string tension and of the lightest glueball mass by the constants s 2 and m 2 , defined as follows: where σ and m 0 ++ stand for the string tension and the lightest glueball mass computed at θ = 0, respectively.To the best of our knowledge, the only study in which an estimate of m 2 and s 2 was attempted is Ref. [46], where their values have been obtained from the computation, at vanishing θ, of the three-points correlation functions between the torelon or glueball interpolating operator and the square of the topological charge.As the calculation of these correlation functions is challenging, only results of limited accuracy could be obtained in Ref. [46].In this work we adopt an alternative approach, performing simulations at imaginary values of θ [31].In a few words, assuming analyticity around θ = 0, we can perform an analytic continuation from real to imaginary values of θ, which can be simulated without any sign problem.Since non analyticities are expected to arise only for θ ∼ π, and since we are dealing with the behavior of σ and m G at small θ, the analyticity assumption is well justified in our case and poses no theoretical issues.We can then obtain σ(iθ I ) and m G (iθ I ) directly for different values of θ I ∈ R, using the same standard algorithms used for computations of the spectrum at θ = 0, and estimate s 2 and m 2 from their small θ I behavior.
Although this approach avoids the complications involved in the evaluation of a threepoint correlation function, the computation of s 2 and m 2 remains non-trivial, especially for large values of N , where two main difficulties arise.The first is the rapid increase of the integrated auto-correlation time of the topological modes in simulations when approaching the continuum limit [55][56][57] (often referred to as the topological freezing problem), which becomes much stronger in the large-N limit.To address this problem, we employ the Parallel Tempering on Boundary Conditions (PTBC) algorithm [58], which has been shown to perform very well in four-dimensional SU(N ) gauge theories [59].The second is that, as can be seen by using standard large-N arguments [46], s 2 and m 2 are expected to scale as 1/N 2 .Hence they are generically expected to have small values, and to have a worsening of their signal-to-noise ratio when the value of N is increased.
As a final remark, we note that a reliable estimate of the coefficients s 2 and m 2 is not only important from a theoretical point of view, but is also directly useful in numerical simulations.As shown in [60,61], these coefficients describe the systematical error that would be introduced by estimating the spectrum from simulations at fixed topological charge Q.If M is the mass of a state and M (Q) is its estimate at fixed topological charge Q, then where M 2 is once again defined by is the topological susceptibility and V the space-time volume.The coefficient m 2 can thus be used to impose an upper bound on the finite size effects introduced by a fixed topological background in the computation of the m 0 ++ glueball mass.This paper is organized as follows: in Sec. 2 we present our numerical setup, discussing the discretization adopted, the update algorithm and the procedure used to evaluate σ and m G ; in Sec. 3 we present our numerical results for the coefficients m 2 and s 2 parametrizing the θ dependence of the string tension and of the lightest glueball state, discussing separately the cases N = 3 and N = 6; finally, in Sec. 4 we draw our conclusions and discuss some open problems.Two appendices report the analytic computations performed in two dimensional U(N ) models and the tables with the raw numerical data of the fourdimensional SU(N ) cases.

Lattice discretization and simulation details
We discretize the SU(N ) pure Yang-Mills theory at θ = 0 on an isotropic hypercubic lattice with L 4 sites using the standard Wilson action: where is the plaquette in position x oriented along the directions µν, a is the lattice spacing, and β is the inverse lattice bare coupling.For the discretization of the topological charge we adopt the standard clover discretization in which ε µνρσ coincides with the standard completely anti-symmetric tensor for positive values of the indices, and its extension to negative values of the indices is uniquely fixed by ε (−µ)νρσ = −ε µνρσ and anti-symmetry.This definition ensures that Q clov is odd under a lattice parity transformation.The action used to generate gauge configurations is thus where the lattice parameter θ L is related to the physical θ angle by −iθ = Z Q θ L , and Z Q is the (finite) renormalization constant of the lattice topological charge [62] To estimate the numerical value of the renormalization constant Z Q it is convenient to use smoothing algorithms such as cooling [63][64][65][66][67][68][69], smearing [70,71], or gradient flow [72][73][74], which dampen the short-scale fluctuations while leaving the global topology of the configurations unaltered.All smoothing algorithms have been shown to be equivalent for this purpose [69,75,76], and in this work we adopt cooling to define an integer-valued topological charge by using [26] with round{x} denoting the closest integer to x, and with α determined by the first nontrivial (i.e., 1 (2.5) In this way we can determine the renormalization constant using [31] Since the dependence of Z Q on the number n cool of cooling steps used to smooth the configurations is only very mild, reaching a plateau for n cool ∼ 10 − 15 for all the β values studied, we define Q using n cool = 20.
Simulations at imaginary values of the θ angle are by now recognized as a cost-effective technique to study θ-dependence on the lattice, as they have been shown to typically outperform simulations carried out at θ = 0 [31, 37-40, 49, 50, 59, 77-85].This is especially true whenever θ = 0 simulations would require the computation of higher-order (i.e., larger than two) correlators or susceptibilities, whose order can be effectively reduced by performing simulations with an external source, then studying the dependence of the results on the source strength.To determine the coefficients s 2 and m 2 at θ = 0 would require the computation of three point functions, see Ref. [46], while using simulations at θ ̸ = 0 we can estimate σ(θ) and m G (θ) as usual, from two-point functions.The values of s 2 and m 2 are then determined from the behavior of σ(θ) and m G (θ) for small θ values.It should be clear that, in this way, we are reducing the statistical errors.However, we have to pay attention not to introduce systematic ones, related to the determination of the O(θ 2 ) behavior of the observables for θ ≈ 0.More details about the computation of σ(θ) and m G (θ) are provided in the next subsection.
For simulations at N = 3, we rely on the standard local updating algorithms usually employed in pure-gauge simulations.More precisely, we adopt a 4:1 mixture of overrelaxation [86] and heat-bath [87,88] algorithms.For N = 6, instead, due to the severe topological freezing experienced by standard local algorithms already at coarse lattice spacing, we adopt the Parallel Tempering on Boundary Conditions (PTBC) algorithm, proposed for two-dimensional CP N −1 models in Ref. [58].This algorithm has indeed been shown to dramatically reduce the auto-correlation time of the topological change both in two-dimensional models [50,58,89] and in four-dimensional Yang-Mills theories [59,85,[90][91][92].
In a few words, in the PTBC algorithm N r replicas of the lattice theory in Eq. (2.3) are simulated simultaneously.Each replica differs from the others only by the boundary conditions imposed on a small sub-region of the lattice, called the defect; these boundary conditions depend on a single parameter, which is used to interpolate between periodic and open boundary conditions.In this way a single replica has periodic boundary conditions, another single replica has open boundary conditions, and the intermediate N r − 2 replicas have "mixed" boundary conditions, i.e., boundary conditions which interpolate between the two previous ones.The state of each replica is updated using heat-bath and over-relaxation local updates, and configuration swaps between different replicas are proposed during the MC evolution.These are accepted or rejected using a Metropolis step.This algorithm allows to exploit the fast decorrelation of Q achieved with open boundaries [93], avoiding at the same time the difficulties related to the lack of translation invariance associated with the presence of open boundaries.For more details on the implementation of this algorithm we refer the reader to Ref. [59], where exactly the same setup adopted here was used.

Extraction of glueball and torelon masses
The starting point to evaluate the torelon mass, needed to extract the string tension, and the lightest glueball mass is the selection of a variational basis of zero-momentumprojected interpolating operators O i .Each O i is a (sum of) gauge invariant single-trace operators of fat-links, built by applying blocking and smearing algorithms to the lattice link variables [51-53, 70, 94-100].For the computation of the lightest glueball mass we employ 4-, 6-and 8-link operators in the A 1 representation of the octahedral group, using a total of 160 operators.To evaluate the torelon mass we use instead 5 operators, built in terms of products of fat-links winding around the time direction once.
As noted before, the only discrete symmetry that can be used to classify the states for non-vanishing values of θ is the charge conjugation C, since parity is not conserved.For this reason, it would be natural to use a variational basis containing definite C only operators, without any projection on definite P representations.However, since we are just interested in the properties of the ground state at small θ values, we can safely use the same standard procedure adopted at θ = 0, i.e., use operators with definite parity.We have indeed verified that the extraction of the lowest glueball and torelon mass does not pose particular challenges and is always characterized by large enough overlaps (A G > 0.9, with A G the squared modulus of the matrix element between the ground state of the selected channel and the vacuum).
The optimal interpolating operator O = i v i O i for the ground state of the selected channel (i.e., the operator with the largest value of A G = |⟨0|O|G⟩| 2 , with |G⟩ denoting the ground state) is the one whose weights v i correspond to the components of the eigenvector of the Generalized Eigenvalue Problem (GEVP) associated to the largest eigenvalue λ (we typically used t 0 /a = 1, performing also some checks using t 0 /a = 2).If we denote by vi the components of this eigenvector, the optimal correlator can be written as The mass of the ground state is then obtained by fitting the functional form in a range where the t-dependent effective mass exhibits a pleateau as a function of the time separation t.Final errors on am G were estimated by means of a standard binned jack-knife analysis.
3 Numerical results

Results for the SU(3) Yang-Mills theory
For N = 3 we performed simulations for 5 different values of the inverse lattice bare coupling β, corresponding to lattice spacings ranging from ∼ 0.1 fm to ∼ 0.05 fm.The lattice size L was chosen large enough to have aL √ σ ≳ 3.5, in which case finite lattice size effects are expected to be negligible to our level of precision, see, e.g., Ref. [37].As a further check that finite size effects are indeed negligible, we compare our estimates of m 0 ++ / √ σ at θ = 0 with the results of Ref. [52], which have been obtained using larger lattices (with aL √ σ ∼ 4 − 5), finding perfect agreement.Using the method described in Sec.2.2, we computed the lightest glueball mass m G and the torelon ground state mass m tor for several values of the lattice parameter θ L .For each value of β and θ L we gathered a statistics of about O(60k) thermalized configurations, separated from each other by 10 updating steps (1 step = 1 heat-bath and 4 over-relaxation sweeps of the whole lattice).The string tension is extracted from the torelon ground state mass m tor by the usual formula [101]: and we explicitly verified that consistent results for s 2 (but not for σ(θ)) are obtained by using simply a 2 σ(θ) = am tor (θ)/L.For this reason, we used Eq.(3.1), including the next-to-leading correction, to compute a 2 σ for all values of θ L .The estimates of m G (θ L ) and σ(θ L ) thus obtained for all the probed values of β and θ L are reported for reference in App.B. In Fig. 1 we instead show a few examples of the obtained effective masses for the same coupling β = 6.40 and for two values of θ L = 0 and 8.The dependence of the lightest glueball mass and of the string tension on θ L can be parameterized, at leading order in θ L , as (see Eq. (1.3)): where the relation θ 2 = −Z 2 Q θ 2 L has been used, and Z Q is the finite renormalization constant introduced in Sec.2.1.The numerical value of Z Q depends on β, and we used the values reported in Ref. [37] in all but one case, namely β = 6.00, in which case Z Q has been estimated anew by using Eq.(2.6) on data at θ = 0 (as in Ref. [37]).
To extract the values of m 2 and s 2 we performed a best fit of our data for am G (θ L ) and a 2 σ(θ L ) using the fit function where A 1 and A 2 are fit parameters.Examples of these fits are displayed in Fig. 2, from which it can be clearly seen that our data are perfectly described by the leading O(θ 2 L ) behavior.To exclude the presence of systematical errors induced by the higher-order O(θ 4 L ) terms, we performed several fits, varying the upper limit of the fit range.When lowering    (10) Table 1.Summary of the results obtained for N = 3.The values of Z Q denoted by an asterisk are from Ref. [37], while the value for β = 6.00 has been computed anew in this work.
the upper limit of the fit range, the errors on the optimal fit parameters increase, but their central values remain well consistent with those obtained by using the full available range, as can be seen from the example shown in Fig. 2. For this reason we report the results obtained by fitting all the available θ L values.
Our estimates of m 2 and s 2 for N = 3 are summarized in Tab. 1.For comparison, in Ref. [46] the values of s 2 and m 2 were estimated by using simulations at θ = 0, where s 2 = −0.077(15) and m 2 = −0.07(4)were obtained at β = 6.00.Bearing in mind that the methods employed for these results are very different, they appear to be in reasonable agreement.Moreover, they are based on a roughly equivalent statistics, which shows that the improvement in accuracy is a benefit of the computational strategy used in the present work.
Since s 2 < 0 the string tension increases when using simulations at imaginary θ, hence we do not expect to observe significant finite-size effects at θ L ̸ = 0.As a further check of the absence of finite-size effects we compared our results for the ratio m 0 ++ / √ σ (extracted from a fit of Eq. (3.3)) with those obtained in Ref. [52] using significantly larger volumes.The comparison between these results is displayed in Fig. 3, from which it is clear that they are perfectly consistent with each other.In particular, assuming just O(a 2 ) corrections, we get the continuum limit to be compared with m 0 ++ / √ σ = 3.405 (21) reported in Ref. [52].The continuum extrapolations of m 2 and s 2 are displayed in Fig. 3.These results are obviously consistent with the presence of just the leading O(a 2 ) finite-a corrections.We thus obtain the continuum extrapolated values m 2 = −0.0083(23), s 2 = −0.0258(14), (continuum extrapolated).Remarkably, the continuum extrapolated value of s 2 is quite close to twice m 2 , which means that the dimensionless ratio m G (θ)/ σ(θ) is almost independent of θ.If we define g 2 by the equation m G (θ) we indeed have the continuum result1 : That the ratio m G (θ)/ σ(θ) is quite insensitive to the value of θ is also true at finite lattice spacing, as can be appreciated from the data reported in Tab. 1 and from the example displayed in Fig. 2.Although we are not aware of any physical argument implying the vanishing of g 2 , this result could suggest that all dimensionless quantities are independent of θ.Such a strong statement can be however shown to be false.In Ref. [83] (see also [84,85]), the θ-dependence of the SU(3) deconfinement critical temperature T c was studied, and it was concluded that where T c is the θ = 0 critical temperature and R = 0.0178 (5).(3.10)This can be recast in units of σ(θ) as: with (using Eq. (3.6)) which is definitely different from zero.

Results for the SU(6) Yang-Mills theory
The general strategy adopted at N = 6 is the same as for the case N = 3.However, obtaining precise results in this case has proven to be a much more challenging task.We thus focused on just two values of the bare inverse lattice coupling, namely β = 25.056 and β = 25.452,corresponding to quite fine lattice spacings.Obviously, using estimates at only two values of the lattice spacing prevents us from performing a reliable continuum extrapolation in this case.The use of the PTBC algorithm was instrumental in reducing the auto-correlation time of Monte Carlo simulations.In particular, for the two values of β considered in this work, the PTBC algorithm allows to reduce the integrated auto-correlation time of the topological modes (at constant CPU time) by a factor of ∼ 20 for β = 25.056[59], and by a factor of ∼ 60 for β = 25.452[90].In simulations performed at inverse coupling β = 25.056we produced and stored O(20k) thermalized configurations at θ = 0 and O(13k) configuration for each non-zero θ value, while for simulations at inverse coupling β = 25.452we produced and stored O(5k) thermalized configurations for each non-zero value of θ (for θ = 0 we used results from a previous study, see Ref. [90]).In all cases, measurements were performed every 10 parallel tempering steps, using the same setup already adopted in Refs.[59,90], to which we refer for more details.
Our results at N = 6 for m 2 and s 2 are summarized in Tab. 2, while all raw data for m G (θ L ) and σ(θ L ) can be found in App.B. As it is clear from data in Tab. 2, and from the example shown in Fig. 4, the θ dependence of m G (θ) and σ(θ) is much milder at N = 6 than it is at N = 3.In particular, we can only provide upper bounds for m 2 .This behavior is compatible with the expectation that m 2 and s 2 are suppressed in the large-N limit.
The accuracy of the data collected at N = 6 is not sufficient to directly test the expected [46] 1/N 2 behavior of m 2 and s 2 .However, our data are definitely consistent with this scaling law, as can be appreciated from Fig. 5, where N 2 m 2 and N 2 s 2 are plotted together for N = 3, 6 as a function of σa 2 .To describe the large-N behavior we can define the parameters s2 and m2 by: Assuming the leading order large-N scaling to be accurate already for N ≥ 3, as is the case for other O(θ 2 ) quantities [59,85], we can estimate these parameters using the results obtained for N = 3: s2 ≃ −0.(46) 0.6297(64) -0.0088(64) 0.04467(25) -0.0084 (26) Table 2. Summary of the results obtained for N = 6.The value of Z Q denoted by an asterisk is from Ref. [59], while the value for β = 25.452 has been computed anew in this work.

Conclusions
In this paper we presented a novel investigation of the θ-dependence of the spectrum of four-dimensional SU(N ) Yang-Mills theories.In particular, we focused on the leading order O(θ 2 ) dependence of the string tension σ(θ) and of the lightest glueball state m G (θ), and estimated the value of the parameters s 2 and m 2 defined in Eq. (1.3).The present study has been carried out by means of lattice simulations performed at imaginary values of the topological θ-angle, using the so called analytic continuation approach.This approach proved to be extremely effective in reducing statistical uncertainties with respect to the previously adopted Taylor expansion method.Our main results have been obtained for the N = 3 theory, for which we were able to perform a controlled continuum extrapolation of m 2 and s 2 , with the results: The case N = 6 was much more challenging, as expected a priori.Larger autocorrelation times make numerical simulations particularly demanding.Moreover, the values of m 2 and s 2 at N = 6 are expected (by large-N scaling) to be suppressed by a factor ≈ 4 with respect to their value at N = 3.To address the first problem we used the Parallel Tempering on Boundary Conditions algorithm, and yet the best we could do was to estimate the values of m 2 and s 2 for only two values of the lattice spacing, although two quite fine ones.
Due to the limited accuracy of the results at N = 6, we cannot perform a stringent test of the expected O(N −2 ) behavior of m 2 and s 2 .However, our data are definitely consistent with this expectation.Based on our previous experiences with O(θ 2 ) corrections (see, e.g., Refs.[59,85]), it is natural to expect the large-N scaling of m 2 and s 2 to be well approximated by the leading O(N −2 ) behavior already for N = 3; using this assumption we get for the coefficients s2 and m2 defined in Eq. (3.13) the estimates s2 ≃ −0.23(1) , m2 ≃ −0.075 (20) . (4.3) As noted in the introduction, the parameters m 2 and s 2 , which parameterize the leading order θ-dependence of the spectrum, also parameterize the way in which the lightest glueball mass and the string tension are affected by the topological freezing, i.e., the systematics induced by using for their estimation an ensemble of gauge configurations with fixed topological charge.Such a quantitative information is very useful.Indeed, given the very fast growth of the integrated auto-correlation time of the topological charge as the continuum limit is approached, it is quite common to perform simulations at fixed Q = 0 with large-N gauge groups.Using our results in Eq. (4.3) and the general formula Eq. (1.5), it is possible to estimate the bias induced by using a fixed topological background.It turns out that already modestly large volumes are sufficient to have a negligible bias: considering for example the case N = 3, in which case χ 1/4 is roughly equal to 1 fm −1 (see, e.g., [102]), we have for V ≃ (1.5 fm) 4 .Using larger values of N this estimate becomes drastically more favorable, since the topological susceptibility changes only slightly (see, e.g., Refs.[34,38,59]), while m 2 scales as 1/N 2 .These estimates constitute an independent confirmation of the results obtained in Ref. [90], with the advantage of providing a quantitative upper bound to the accuracy that can be achieved when using simulations at fixed topological sector.
The results presented in this paper can be extended quite naturally in several different ways.One possibility is to accurately investigate s 2 and m 2 for N > 3, in order to quantitatively asses the N dependence of these coefficients.From the previous discussion it should be clear that this is not an easy task, and some ideas are required to further improve the signal-to-noise ratio.A second possibility is to study the excited glueball spectrum.In particular, it would be very interesting to understand (both theoretically and numerically) how the O(θ 2 ) correction to the mass depends on the state considered.Indeed, these corrections can not be all independent from each other, since they are indirectly related to the θ-dependence of the free energy in the confined phase.This can be easily understood by using hadron resonance gas models like, e.g., those discussed in [103,104], where it was shown that a determination of the glueball masses is sufficient to obtain quantitatively accurate estimates of thermodynamical quantities.Finally, it would be very interesting to study models in which s 2 and m 2 or, more generally, the θ-dependence of the spectrum, can be investigated analytically (and is non-trivial).Two-dimensional CP

Figure 1 .
Figure1.Examples of t-dependent effective masses obtained for the lightest glueball (left) and torelon (right) states for N = 3, β = 6.40 and θ L = 0 and 8, computed using Eq.(2.10) and the method described in Sec.2.2.Shaded bands represent our final results for a G and am tor , obtained from a best fit of the optimal correlator according to Eq. (2.9) for t/a ≥ 3, i.e., in the range where effective masses exhibit a plateau.

2 Figure 2 .
Figure 2. Results for N = 3 with β = 6.40.Examples of fits performed by using the functional form Eq. (3.3) for several quantities: (top left panel) the lightest glueball mass am G , (top right panel) the string tension a 2 σ, (bottom left panel) the ratio m G / √ σ.In the bottom right panel we show the dependence of the estimates of m 2 and s 2 on the upper limit of the fit range θ (max) L .

2 Figure 3 .
Figure 3. Left panel: continuum scaling of our N = 3 results for m 0 ++ / √ σ, compared with data reported in Ref. [52].The dashed line is the result of a best fit of our data assuming O(a 2 ) scaling corrections.Right panel: continuum scaling of our N = 3 results for m 2 and s 2 and their continuum extrapolation assuming O(a 2 ) scaling corrections.

Figure 4 .
Figure 4. Results for N = 6 β = 25.452.Examples of fits performed by using the functional form Eq. (3.3) for several quantities: (top left panel) the lightest glueball mass am G , (top right panel) the string tension a 2 σ, (bottom panel) the ratio m G / √ σ.

Figure 5 .
Figure 5. Continuum scaling of N 2 m 2 and N 2 s 2 for N = 3 and N = 6.
(2.10)and the method described in Sec.2.2.Shaded bands represent our final results for a G and am tor , obtained from a best fit of the optimal correlator according to Eq. (2.9) for t/a ≥ 3, i.e., in the range where effective masses exhibit a plateau.

Table 4 .
N −1 models are natural candidates.between INFN and Cineca (under project INF22 npqcd), and on the Plymouth University cluster.Summary of all obtained results for N = 6.