Topological susceptibility in high temperature QCD: a new investigation with spectral projectors

. We compute the topological susceptibility of high temperature QCD with 2 + 1 physical mass quarks using the multicanonical approach and the spectral projector estimate of the topological charge. This approach presents reduced lattice artifacts with respect to the standard gluonic one, and makes it possible to perform a reliable continuum extrapolation.


Introduction
In recent years the study of the topological properties of QCD at high temperatures has been the subject of several Lattice QCD investigations [1][2][3][4][5][6].This is motivated by the general interest in the non-perturbative physics of high-temperature QCD, but also by its phenomenological implications for axion physics and cosmology.
For asymptotically high temperatures the topological susceptibility χ(T ) = Q 2 /V (we denote by Q the topological charge and by V the space-time volume) can be computed by semiclassical methods and perturbation theory: in the Dilute Instanton Gas Approximation (DIGA) from a one loop computation in the instanton background it is possible to obtain the result [7,8]: where N c , N f are the number of colors and light flavors respectively.This result can however be trusted only for T Λ QCD ≈ T c , and it was thus suggested in [9] to use first principles Lattice QCD results for χ(T ) to estimate the upper bound of the axion coupling constant f a [10][11][12].
The computation of χ(T ) in the high temperature regime by means of Lattice QCD simulations is however a very challenging task.This is due to the simultaneous presence of three complementary numerical obstructions that have to be faced.
The first problem is related to the fact that the topological susceptibility approaches zero quite rapidly by increasing the temperature (Eq.( 1) predicts χ(T ) ∝ T −8 ), thus in the high tempertature regime the probability of visiting configurations with Q 0 is strongly suppressed on finite volume simulations.Very large statistics are needed to observe a sufficient number of fluctuations of Q to reliably estimate χ.
The second problem is related to the presence of very large lattice artifacts.In continuum QCD the index theorem connects the topological charge of a configuration with the number of zero modes of the massless Dirac operator / D, and the presence of zero modes suppresses in the chiral limit the probability of finite volume configurations with non zero topological charge.On the lattice chiral symmetry is explicitly broken when using Wilson or staggered fermion discretizations.As a consequence no exact zero modes are present, but only Would-Be Zero Modes (WBZMs), which make the suppression of Q 0 configurations less efficient with respect to the continuum.The numerical determination of χ is thus affected by huge lattice artifacts, and its continuum extrapolation requires particular attention.
The third problem is the so called topological freezing problem: the effectiveness of the commonly used local update algorithms in varying the topological charge of the configurations is subject to a severe form of critical slowing down as the lattice spacing is decreased.Indeed the increase of the autocorrelation time for topological observables is consistent with an exponential in the inverse lattice spacing [13,14].This makes extremely difficult to use, for temperatures of a few hundred MeVs, lattice spacings of about 0.01 fm or smaller.
The use of specific strategies to overcome these issues is mandatory to estimate χ(T ) at high temperatures, and different approaches have been proposed in the literature to this end.For instance, in Ref. [3] (see also [2]) the Authors compute χ(T ) by restricting simulations to the Q = 0 and |Q| = 1 sectors, and to reduce lattice artifacts related to the use of staggered fermions they a posteriori reweight configurations by using the lowest eigenvalues of / D, assumed to be WBZMs.
In this proceeding we discuss some of the results recently obtained in [15], where some progress has been made towards a lattice determination of χ(T ) in QCD at the physical point (with N f = 2 + 1) without any extra ad hoc hypothesis.To make this possible we use two different strategies to cope with the first and the second of the problems discussed above.An important point to be stressed is that these strategies help resolving or alleviating the computational problems but they are known not to introduce any biases in the final results.
To deal with the problem of the dominance of the Q = 0 sector at high temperature we use an adaptation of the multicanonical algorithm [16], analogously to what was done in [5,[17][18][19].In short, the idea is to add to the lattice action a bias topological potential to enhance the probability of visiting the suppressed topological sectors.Expectation values with respect to the original distribution are then computed by using standard reweighting to remove the XV effect of the bias potential.We explicitly note that this reweighting has no overlap problem and can thus be carried out without introducing any systematical errors.
To reduce the size of the lattice artifacts, for the topological susceptibility we adopt a definition based on the Spectral Projectors (SP) method [20][21][22][23][24], adapted in Ref. [25] to the case of staggered fermions.In this discretization Q is defined as the sum of the chiralities of all the eigenmodes of the Dirac operator lying below a certain threshold M.This definition is theoretically well-defined, as it can be shown to converge to the correct continuum limit, and the choice of M can be used to reduce the lattice artefacts with respect to the standard gluonic definition of the topological charge.This approach increases our control of the systematic uncertainties related to the continuum extrapolation of χ, thus alleviating the necessity for extremely fine lattice spacings to reduce lattice artifacts.

Numerical setup
In our simulations we adopted a discretization of N f = 2 + 1 QCD using rooted stouted staggered fermions for the quark sector and the tree-level Symanzik improved gauge action for the gluon sector.The staggered Dirac operator is defined by using the gauge links U (2)  μ , obtained by applying to the gauge configuration 2 levels of isotropic stout smearing [26] with ρ stout = 0.15.The bare coupling β and masses ms and mu = md ≡ ml have been tuned in order to move on the Line of Constant Physics (LCP) corresponding to the physical values of the pion mass m π 135 MeV and of the ratio ms / ml = m s /m l 28.15 [27][28][29].
For the topological charge we adopt two different definitions.The first one is a standard gluonic definition, obtained by using the simplest clover discretization with definite parity of the topological charge density [30,31] on cooled configurations, in order not to explicitly deal with lattice renormalizations [32].In particular we use the definition introduced in [33] (see [15] for more details) and 100 cooling steps, since this number was checked to be sufficient to reach a plateau for all explored lattice spacings.
To introduce the SP definition of the topological susceptibility we first of all need the projector P M on the vector space spanned by the eigenstates of iD stag [U (2) ] (the same operator entering the lattice action) with eigenvalues |λ| ≤ M: The SP definition of the bare topological charge is then where the factor n t = 2 d/2 = 2 2 is needed to take into account the taste degeneracy of the staggered spectrum.As discussed in Ref. [25], this definition is affected only by a multiplicative renormalization, and the SP expression of the renormalized topological susceptibility can be written as Spectral traces have been estimated by directly computing the first 200 smallest eigenvalues and eigenvectors of iD stag [U (2) ].From the eigenvalues and the eigenvectors of iD stag [U (2) ] spectral traces are then practically computed by using ,  The cut-off mass M is a free parameter whose specific value is irrelevant in the continuum limit.It is however important to keep M constant in physical units as the continuum limit is approached.Since the cut-off M renormalizes as a quark mass [25], to keep M constant in physical units it is sufficient to keep M/ m f constant as we move m f along the LCP, where M stands for M in lattice units.Adopting this procedure we thus expect the usual continuum scaling χ where c SP (M/m f ) parameterizes the finite lattice spacing corrections.
The last ingredient we need is the multicanonical algorithm: in this algorithm a topological bias potential V topo (Q mc ) is added to the gauge action, in order to enhance the probability of visiting those topological sectors that would be otherwise strongly suppressed.The quantity Q mc is a gluonic discretization of the topological charge, and the multicanonical algorithm is stochastically exact for any choice of Q mc .However, to make the algorithm efficient, Q mc must have a reasonable overlap with the topological charge used in the measures.Moreover, to avoid the need for very small integration steps in the Hybrid Monte Carlo, Q mc can not be "too peaked" on integer values.Also the precise form of the function V topo is largely irrelevant, its only role being that of increasing the probability of the Q 0 sectors, without completely depleting the Q = 0 sector to avoid overlap problems in the reweighting procedure (see [15] for the specific form adopted).For a generic observable O, expectation values with respect to the original path-integral distribution O are finally recovered by reweighting: where • mc denotes the average computed with the topological bias in the action.

Results
To apply the SP definition of the topological susceptibility we have first of all to fix the value of the parameter M. While from a theoretical point of view any value of M can be used, from the practical point of view a wise selection of M can help reducing the lattice artifacts.Comparison of the continuum extrapolations of χ 1/4 at T 430 MeV obtained with the gluonic and the SP discretizations.Data from Refs.[3,4] are also reported for comparison.Data from [3] has been modified to remove the isospin-breaking term, while data from [4] has been massextrapolated to physical point using χ 1/4 ∼ m π .
To understand which can be a reasonable interval for M we looked at the scatter plot of the (absolute value of the) mode chirality r λ as a function of |λ| for the different lattices studied.An example of these scatter plots is shown in Fig. 1 for the case of T 430 MeV and two different values of the lattice spacing, from which the emergence of two clusters is evident when reducing the lattice spacing.The cluster of points with high chirality is likely mainly composed of WBZMs (alhough a precise identification of WBZMs is far from trivial, see [15] for a discussion of this point), which in the continuum limit are the only modes which contribute to the SP definition of Q.For this reason we selected M values in the region of |λ|/m s which separates the two cluster, varying M in the region delimited by the dashed vertical lines in Fig. 1 to check for systematics of the continuum extrapolation.
The continuum extrapolation of the results obtained at T 430 MeV is shown in Fig. 2, where we compare results obtained by the SP method with the result obtained on the same ensamble of configurations using the gluonic definition of the topological charge.It is clear that the SP method has much smaller lattice artifacts than the gluonic approach, moreover by investigating the dependence of the result on M we also have a better control on the systematics of the continuum limit, something which is not possible when using the gluonic definition.In Fig. 2 we also report, for comparison, data from [3] (modified to remove the a posteriori DIGA-like isospin-breaking effect introduced in [3]) and from [4] (extrapolated in mass by the DIGA-like scaling χ 1/4 ∼ m π , since an unphysical pion mass was used in [4]).The same procedure was followed for five temperatures between 230 and 570 MeV, and in Fig. 3 we show our final results for the topological susceptibility χ(T ) as a function of the temperature, in which error bars take into account both statistical and systematic uncertainties.Results obtained by using the SP method have a significantly smaller error then those obtained by using the gluonic definition of the topological charge, moreover, as already stressed, the error is not only smaller but also more reliable when using the SP definition.Data clearly follow, with the possible exception of T 230 MeV, a simple power law behav- Behavior of χ 1/4 as a function of T/T c in log-log scale.Gluonic points are slightly shifted to improve readability.Starred points represent results taken from Ref. [3] (with isospin-breaking factor removed), while the shaded area represents the gluonic determinations reported in Ref. [4], massrescaled according to χ 1/4 ∼ m π .
ior, and a fit performed for T > 230 MeV by using the ansatz provides for b (by using the SP determinations) the optimal value value b SP = 2.63(81), which is consistent with the DIGA prediction b DIGA = 2.
In Fig. 3 we also report for comparison data from [3,4], with the same modification noted above to make the results refer to the same physical setting.While the global behavior of the data is the same, and in particular all data scale with a power-law consistent with DIGA, some discrepancies are clearly evident: all our data appear to be systematically above the previous determinations.

Conclusions
In this proceeding we presented our results concerning the behavior of the topological susceptibility χ(T ) in the high temperature regime of QCD with N f = 2+1 quarks at the physical point [15].
To perform this computation we relied on the discretization of the topological charge through Spectral Projectors on the eigenmodes of the staggered Dirac operator.With this choice we find lattice artifacts that are much smaller then those obtained by using the standard gluonic definition of χ.The problem of the dominance of the Q = 0 sector is instead addressed by using the multicanonical method, already adopted for this purpose in [5].
The SP definition of the topological susceptibility introduces a new free parameter, which is the cut-off of the sum on the chiralites of the eigenmodes of the lattice staggered Dirac operator.Any value of M (kept constant in physical units as a → 0) would in principle provide the correct continuum limit of χ.To fix a specific value of M we first of all preliminary investigated the chirality of the staggered Dirac eigenmodes, identifying two clusters of modes, which can be roughly associated to WBZMs and to modes with nonvanishing eigenvalue , 06001 (2022) https://doi.org/10.1051/epjconf/202227406001t h Quark Confinement and the Hadron Spectrum EPJ Web of Conferences 274 XV in the continuum limit.While an unambiguous identification of a would-be-zero-mode is a nontrivial task, the purpose of this preliminary investigation was only the identification of a region of "reasonable" values for M. The systematic related to the choice of M is then investigated and included in our final uncertainty for χ.As a matter of fact, this systematic is typically the dominant contribution to the final error on χ for the temperatures investigated in this work.
We explored five different values of T in the high temperature regime, ranging from ∼ 200 MeV to ∼ 600 MeV, investigating the behavior of χ as a function of the temperature T .Comparing our numerical results with semiclassical expectations, we find that a decaying power law well describes our data in the whole explored range; in particular the effective exponent of this power law is well in agreement with the DIGA prediction for T 300 MeV, i.e., for T/T c 2.
Our results for χ 1/4 obtained by using the SP method show a ∼ 2 − 3 standard deviation tension in the 300 MeV T 400 MeV range with respect to the previous determinations reported in [3,4], with our estimates of χ 1/4 (T ) systematically pointing to larger values.The same behavior, when observed in the gluonic determinations of χ, could be ascribed to a problem of the continuum extrapolation, which due to large lattice artefacts is unable to capture the true asymptotic O(a 2 ) scaling and introduces a bias in the extrapolation.The lattice spacing dependence of the SP determinations is however much milder than that of the gluonic estimates, and such an interpretation of the disagreement does not seems likely in this case.
The final picture emerging from the comparison carried out in Fig. 3 is that we still do not have a quantitatively complete understanding of the behavior of χ(T ) in the high temperature regime of QCD, and further studies are required to clarify the sources of the observed tensions between the different determinations.
It thus seems crucial to refine our present estimates of the topological susceptibility, in order to make the comparison with the results of Refs.[3,4] more stringent.For this purpose simulations with larger statistics and smaller lattice spacings are required.Despite being essential to improve the present results and to explore larger temperatures, simulations at smaller lattice spacings are practically unfeasible with standard simulations algorithms, due to the exponential nature of the topological critical slowing down.A promising approach that can help to overcome this problem is the one introduced in [34], and already applied to two dimensional CP N−1 models [34,35] and four dimensional Yang-Mills theories [36,37].

Figure 1 .
Figure 1.Scatter plot of the (absolute value of the) chirality r λ vs |λ|/m s for the first 200 eigenvalues of configurations generated at T 430 MeV.(Left) for a 0.0381 fm, (Right) for a 0.0286 fm.The two vertical dashed lines at |λ|/m s = 0.03 and 5 delimit the M-range used to check for systematics.

, 3 SP, M/m s = 0. 5 Figure 2 .
Figure 2. Comparison of the continuum extrapolations of χ 1/4 at T430 MeV obtained with the gluonic and the SP discretizations.Data from Refs.[3,4] are also reported for comparison.Data from[3] has been modified to remove the isospin-breaking term, while data from[4] has been massextrapolated to physical point using χ 1/4 ∼ m π .

,Figure 3 .
Figure3.Behavior of χ 1/4 as a function of T/T c in log-log scale.Gluonic points are slightly shifted to improve readability.Starred points represent results taken from Ref.[3] (with isospin-breaking factor removed), while the shaded area represents the gluonic determinations reported in Ref.[4], massrescaled according to χ 1/4 ∼ m π .