Topological susceptibility of $N_f=2+1$ QCD from staggered fermions spectral projectors at high temperatures

We compute the topological susceptibility of $N_f=2+1$ QCD with physical quark masses in the high-temperature phase, using numerical simulations of the theory discretized on a space-time lattice. More precisely we estimate the topological susceptibility for five temperatures in the range from $\sim200$ MeV up to $\sim600$ MeV, adopting the spectral projectors definition of the topological charge based on the staggered Dirac operator. This strategy turns out to be effective in reducing the large lattice artifacts which affect the standard gluonic definition, making it possible to perform a reliable continuum extrapolation. Our results for the susceptibility in the explored temperature range are found to be partially in tension with previous determinations in the literature.


Introduction
The study of the topological properties of QCD at high temperatures is of utmost importance not only to provide better insight into the non-perturbative regime of this theory, but also because of its phenomenological implications for axion physics and cosmology. This justifies the interest in this topic, which has been the subject of several recent Lattice QCD investigations [1][2][3][4][5][6].
The axion is an hypothetical particle whose existence is predicted by the Peccei-Quinn solution of the strong CP-problem [7][8][9][10], that was also early recognized as a possible Dark Matter candidate. Since the axion is directly coupled to the topological charge Q defined by (F µν is the QCD field strength) Q = 1 32π 2 ε µνρσ Tr{F µν (x)F ρσ (x)}d 4 x , (1.1) the axion square mass is proportional to the topological susceptibility χ = Q 2 /V (V is the four-dimensional space-time volume). More precisely m 2 a = χ/f 2 a , where f a is the a priori unknown axion coupling scale. An important peculiarity of axion physics is that it is possible, modulo some general cosmological assumptions, to put an upper bound on the energy scale f a [11][12][13]. The specific value of this upper bound is fixed by the present time Dark Matter abundance and by the temperature dependence of the axion effective potential, i.e., mainly by the temperature dependence of the QCD topological susceptibility χ. Since an upper bound on f a provides a lower bound for the axion mass, the study of the QCD topological observables at finite temperature provides an essential input for current and future experimental axion searches (see [14] for a recent review of the experimental bounds).
When the temperature is asymptotically high, one can compute χ(T ) by combining semiclassical methods and perturbation theory. Assuming instantons to be well separated and thus approximately not interacting with each other (Dilute Instanton Gas Approximation, DIGA for short), and performing a one loop computation in an instanton background, it is possible to obtain the result [15,16]: where a logarithmic dependence of χ on the temperature has been implied and N c , N f are the number of colors and light flavors respectively. Being based on a combination of semiclassical and perturbative approximations, this result is expected to be trustworthy only for T ≫ Λ QCD ≈ T c . Nevertheless, due to the absence of more reliable computations, this result has been routinely used to estimate the cosmological axion relic density needed to constrain the axion coupling f a . In performing this computation the expression of χ(T ) has however to be used starting from T ≫ T c but reaching temperatures of the order of 1 ÷ 10 GeV (see, e.g., Ref. [17]), where non-perturbative deviations from the asymptotic high-T regime could be relevant. To avoid introducing systematic errors in the computation of the axion coupling bound, it was thus suggested to use first principles Lattice QCD results for χ(T ) instead of the corresponding DIGA expression [18]. When exploring the high temperature regime of QCD with lattice simulations, there are however several nontrivial numerical problems that have to be faced. Among them the most notable are: The topological susceptibility drops very rapidly with the temperature (DIGA predicts c ≈ 8 for 3 light flavors, cf. Eq. (1.2)), thus the probability of visiting configurations with Q = 0 is rapidly suppressed as the temperature is increased. This sampling problem has a physical origin and can be understood in terms of the vanishing of the variance of the topological charge probability distribution P (Q): since Q 2 = χV ≪ 1 on affordable volumes, we have P (Q = 0) ≫ P (|Q| > 0). From the numerical point of view, this implies the necessity of collecting very large statistics in order to observe a sufficient number of fluctuations above zero to reliably compute χ.
(ii ) Explicit breaking of chiral symmetry and large lattice artifacts According to the index theorem, the appearance of zero-modes in the spectrum of the continuum massless Dirac operator / D is related to the topological charge Q of the background gauge field: where n + (n − ) is the number of zero modes with positive (negative) chirality. For this reason, a finite volume configuration with non-zero topological charge enters the continuum path integral with a weight that is suppressed in the chiral limit as a power of the light quark mass. On the lattice, however, if the sea quark discretization explicitly breaks the chiral symmetry, no zero mode is present and no mode is chiral, being the smallest eigenvalues shifted by cut-off effects. Such Would-Be Zero Modes (WBZMs) make the suppression of Q = 0 configurations less efficient compared to the continuum. As a consequence the numerical determination of χ is affected by large lattice artifacts, and its continuum extrapolation requires particular care in the analysis of systematic uncertainties. In principle this issue could be solved by simulating extremely fine lattice spacings, which is in practice prevented by the topological critical slowing down problem.
(iii ) Topological critical slowing down Local updating algorithms, such as the Rational Hybrid Monte Carlo (RHMC) [19,20], become less and less effective in changing the topological charge of the configurations as the lattice spacing is decreased. While the critical slowing down generically affects all observables, it is particularly severe for the topological observables, for which the increase of the autocorrelation time is consistent with an exponential in the inverse lattice spacing [21][22][23][24][25][26]. In practice, if the lattice spacing is small enough ergodicity is lost and the whole Monte Carlo simulation remains trapped in the same topological sector; as a matter of fact, this problem is also referred to as the "freezing problem". This makes extremely difficult to reach, for the temperatures studied in this work, lattice spacings of the order or smaller than 0.01 fm.
Overcoming these issues is mandatory to provide reliable lattice computations of χ(T ) at high temperatures, and several different strategies have been proposed in literature to this end. For instance, in Ref. [3] (see also [2]), the authors compute χ(T ) bypassing issues (i ) and (iii ) simultaneously by restricting only to Q = 0 and |Q| = 1, where neglecting contributions from higher-charge sectors is justified on the basis of the DIGA itself. The computational problem (ii ), related to the adoption of staggered fermions in the lattice action, was instead addressed in Ref. [3] by reweighting configurations according to their expected continuum lowest eigenvalues of / D, thus trying to correct a posteriori for the inaccurate sampling of chiral modes. A crucial point to perform such a reweighting is the determination of WBZMs, that are not so easily distinguishable from non-chiral modes without some extra assumptions due to the explicit breaking of chiral symmetry on the lattice.
The main goal of this work is to make progress towards an independent determination of χ(T ) in full QCD from the lattice without any extra ad hoc hypothesis. In particular, we will compute this quantity in N f = 2 + 1 QCD at the physical point combining the two following strategies: (a) Topological susceptibility from staggered fermions spectral projectors In order to reduce the magnitude of lattice artifacts, we adopt a definition of the topological susceptibility based on the Spectral Projectors (SP) method [27][28][29][30][31], introduced for staggered fermions in Ref. [32]. This discretization, based on a lattice version of the index theorem, consists of defining Q as the sum of the chiralities of all the modes lying below a certain threshold M . The SP topological susceptibility is a theoretically well-defined quantity, 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. Our approach will allow us to have a better control on the systematic uncertainties related to the continuum extrapolation of χ already in the lattice spacing range employed in typical simulations, thus alleviating the necessity for extremely fine lattice spacings to reduce lattice artifacts, whose use would require a specific strategy to deal with the issue (iii ).
(b) Multicanonical algorithm The multicanonical algorithm [33] was introduced to study strong first order phase transitions, and was adapted in [5,25,34,35] to deal with the problem of the dominance of the Q = 0 sector at high temperatures (issue (i ) above). The idea is to add to the lattice action a fixed (as opposed to the case of metadynamics [36]) bias topological potential in order to enhance the probability of visiting suppressed topological sectors, thus effectively enhancing the fluctuations of Q during the Monte Carlo evolution. Expectation values with respect to the original path-integral distribution are then exactly reproduced by using standard reweighting. With respect to the standard approach, in this way a much better statistical accuracy is achieved in the determination of the relative weights of the topological sectors.
We here anticipate that our results for χ(T ) obtained from spectral projectors show some difference with respect to those reported in Ref. [3], as we observe a ∼ 2 − 3 standard deviation tension in a temperature range between 300 and 400 MeV. We also observe, in the same temperature range, a milder ∼ 2 − 2.5 standard deviation tension with respect to the determinations of Ref. [4], which however had to be mass-extrapolated to be compared with our results, as will be discussed in the following. This paper is organized as follows: in Sec. 2 we present our numerical setup, discussing in particular the spectral discretization of the topological susceptibility (a) and our implementation of the multicanonical algorithm (b). In Sec. 3 we present continuumextrapolated results for the topological susceptibility in N f = 2 + 1 QCD at the physical point. First, we apply the SP approach to the T = 0 case, which is used as a test-bed to compare this method with the standard gluonic approach. Then, we adopt the same SP method in combination with the multicanonical approach to compute the topological susceptibility as a function of the temperature in the range 200 MeV T 600 MeV. Our results are then compared with other determinations in the literature and with DIGA predictions. Finally, in Sec. 4 we draw our conclusions and discuss future outlooks of this work.

Lattice action
We discretize N f = 2 + 1 QCD on a N 3 s × N t lattice adopting rooted stouted staggered fermions for the quark sector and the tree-level Symanzik improved gauge action for the gluon sector. The partition function Z LQCD is thus given by where u = d ≡ l and s denote the two mass-degenerate light quarks and the strange quark respectively. The staggered Dirac operator D stag is defined by using the gauge links U (2) µ , obtained by applying to the gauge configuration n stout = 2 levels of isotropic stout smearing [37] with ρ stout = 0.15, thus: The tree-level Symanzik-improved Wilson action is instead expressed in terms of the nonstouted gauge links: where Π

Topological charge discretizations
Topological charge definitions can be divided into two broad groups: gluonic and fermionic ones. The simplest gluonic definition is the clover one, which is a straightforward discretization of Eq. (1.1) with definite parity [41,42]: where Π (1×1) µν (x) is the plaquette and the Levi-Civita symbol with negative entries is defined by ε µνρσ = −ε (−µ)νρσ and complete antisymmetry. When computing the susceptibility χ = Q 2 /V using Q clov , both multiplicative [43] and additive renormalizations appear 1 : Such renormalizations are due to ultraviolet (UV) fluctations at the scale of the lattice spacing and must be properly subtracted in order to recover the proper continuum scaling of χ gluo . To avoid dealing with these renormalization constants, smoothing algorithms are commonly employed to dampen UV fluctuations while leaving the topological content of the gauge fields unchanged. Several methods have been proposed, such as cooling [45][46][47][48][49][50][51], smearing [37,52] and gradient flow [53,54], all giving consistent results when properly matched [51,55,56].
To compute the gluonic susceptibility, in this work we adopt the cooling method, which is numerically very convenient. Since even after cooling the topological charge is noninteger (although its typical values get closer and closer to integers as the lattice spacing is reduced) we adopt the following prescription to assign an integer topological charge to configurations [57,58]: where "round" means that the quantity αQ (2.7) The parameter α is thus chosen in order to center the peaks of the distribution of αQ (cool) clov at integer values, and the constraint x ≥ 1 is required to exclude the trivial minimum x = 0. In the end, our gluonic susceptibility is: In our simulations we observe that after n cool ∼ 100 cooling steps χ gluo has reached a plateau for all explored lattice spacings. For this reason, we compute the gluonic susceptibility for that number of cooling steps in all cases. Slightly changing the value of n cool resulted in no change in the obtained results for χ gluo for each lattice spacing. It is also possible to adopt fermionic definitions of the topological charge, based on the properties of the spectrum of the lattice Dirac operator. In the continuum limit, according to the index theorem, it would be sufficient to compute the sum of the chiralities of the zero-modes u 0 , as On the lattice, since no exact zero-mode exists when using a non-chiral fermion discretization, every mode contributes and the trace becomes a sum over all modes. To properly define the topological charge, we introduce the projector on the eigenspace spanned by the eigenstates of iD stag [U (2) ] (i.e., the same operator we have included in our lattice action) with eigenvalues |λ| ≤ M : Our SP definition of the bare topological charge is: where the factor n t = 2 d/2 = 2 2 takes into account the taste degeneration of the staggered spectrum.
As discussed in Ref. [32], this definition is affected only by a multiplicative renormalization, since the fast decay of the spectral projector at infinity eliminates any additive renormalization. This multiplicative constant can be expressed in terms of traces of Γ 5 and of the spectral projector P M , thus, we are able to obtain a fully-spectral renormalized definition of the topological charge [32]: The SP expression of the topological susceptibility can then be written as (2.14) Spectral traces can be computed by several means. For example, in Refs. [28,29] noisy estimators are used to this end. In this work we follow the same strategy of Ref. [32] and we compute the first 200 smallest eigenvalues and eigenvectors of iD stag [U (2) ] using the PARPACK package [59], so that P M is obtained directly from Eq. (2.10). On our larger lattices, the computational cost to obtain the 200 lowest-lying eigenvalues of iD stag [U (2) ] turned out to be about a factor of ∼ 20 larger than the computational cost needed to perform a single RHMC step for the same lattice. Once the eigenvalues and the eigenvectors of iD stag [U (2) ] are obtained, spectral traces are then practically computed as follows: where ν(M ) is the total number of eigenmodes whose eigenvalues lie below M . The cut-off mass M is a free parameter, and from the index theorem we know that its specific value is irrelevant in the continuum limit. However, M should be kept constant in physical units as the continuum limit is approached in order to observe the usual continuum scaling of the susceptibility (i.e., O(a 2 ) corrections in the staggered case), as discussed in Ref. [28]. The cut-off mass renormalizes as a quark mass [32], i.e., M R = Z −1 S M , where Z S is the renormalization constant of the staggered flavor-singlet scalar fermionic density. For this reason, the ratio between M and any of the quark masses is a renormalizationgroup invariant quantity: . This strategy allows to completely avoid the computation of Z S and will be the one adopted in this work. When the continuum limit is taken at fixed M/m f , we thus expect the usual scaling: Previous studies of QCD at zero temperature, performed with twisted mass Wilson fermions and using twisted mass Wilson SP, have shown that the SP determination of the topological susceptibility displays much smaller lattice artifacts than the gluonic one [31]. Although we still do not have a complete quantitative understanding of this fact, an intriguing possible interpretation exists (see also the discussion at the end of Ref. [32]). When using a non-chiral fermion discretization we are sampling a lattice distribution that, for what concerns the low-laying modes relevant for topology, can have quite large lattice artifacts due to the explicit chiral symmetry breaking. If we use a gluonic definition of the topological charge, we have lost any direct connection with lattice chirality and the error introduced in the sampling propagates to the measures. If instead we use the SP definition built from the same Dirac operator used to weight the configurations, there is the possibility that the measure partially corrects the error in the sampling. Of course there is a priori no solid reason to exclude the possibility that the two different errors sum up instead of canceling each other, so this argument can not be considered in the present form conclusive. To understand to what extent such a cancellation exists, it would be very interesting to perform a study using different fermion discretizations in the generation of the configuration and in the construction of the SP used to estimate the topological charge.
As a final remark we recall that, as already discussed in Ref. [32], at finite temperature a possible ambiguity about the computation of the renormalization constant in Eq. (2.13) could arise. In particular, one could wonder if Z (stag) SP should be computed at finite T or at zero T . In principle, such quantity should be computed in the latter case, which would require to perform, along each finite temperature simulation, a twin run with the same parameters but at zero temperature to compute Z (stag) SP . However, as already discussed in Ref. [32] in the quenched theory and as we will verify also in the following in the presence of dynamical fermions, the computation of χ (stag) SP at finite temperature gives fully consistent results both when Z (stag) SP is obtained from the same finite-T ensemble employed for the computation of Q (stag) SP or from a corresponding ensemble at zero temperature, as long as the spectral traces appearing in Eq. (2.13) are computed from the staggered spectrum obtained imposing periodic boundaries along the temporal direction for D stag 2 . Thus, in our finite temperature simulations, we will compute the spectrum of the Dirac operator both for periodic and anti-periodic boundary conditions along the temporal direction, and in Eq. (2.14) we will adopt the former to compute Z (stag) SP and the latter to compute Q (stag) SP,bare .

Multicanonical algorithm
The multicanonical approach consists in adding a topological bias potential V topo (Q mc ) to the action in order to enhance the probability of visiting those topological sectors that would be otherwise strongly suppressed: (2.17) The quantity Q mc is a suitable discretization of the topological charge, which in general differs from the one adopted for the measurements Q gluo . While the multicanonical algorithm is stochastically exact for any choice of Q mc , to make it more efficient than the case V topo = 0 the discretization Q mc has to satisfy a couple of requirements. First, Q mc must have a reasonable overlap with the charge used in the measures, otherwise V topo (Q mc ) would not work as a bias but just as noise. Second, it is also important that Q mc is not "too peaked" at integer values, to avoid the need for very small integration steps in the Hybrid Monte Carlo. The partition function in the presence of the topological potential Z (mc) LQCD is the following: where · mc refers to expectation values computed in the presence of the topological bias.

Topological susceptibility at zero temperature
At zero temperature, the value of the topological susceptibility in QCD can be reliably computed using Chiral Perturbation Theory (ChPT) [60][61][62][63][64][65], and the value obtained in this way constitutes a useful benchmark for lattice determinations. Using the Next-to-Leading Order (NLO) expression of Ref. [64] one gets for the case m u = m d at the physical point the estimate [1]: Moreover, several gluonic determinations from the lattice have been reported in the literature, see, e.g., Refs. [1,3]. We will use the computation of the zero temperature topological susceptibility as a test to validate our implementation, and to verify that the spectral projector approach has significantly smaller lattice artifacts than the standard gluonic one.
Our zero-temperature simulations have been performed on hypercubic N 4 s lattices and the simulation parameters adopted are reported in Tab. 1. The lattice size was chosen so that L s ≡ aN s ∼ 2.6-3.2 fm, which is sufficient to keep finite-size effects within our typical statistical errors [1]. At zero temperature, there is no need to introduce a bias potential in the action, since several topological sectors are naturally explored during the Monte Carlo evolution. As an example, we show the history of the topological charge Q gluo for our finest lattice spacing in Fig. 1   To compute χ using spectral projectors we have to choose the cut-off mass M , which in the following computations will be normalized to the strange quark mass m s . Since we know that in the continuum limit only zero modes provide a non-vanishing contribution to the topological charge, a priori the optimal possibility would be to fix a value of M/m s large enough to include in the spectral sums in Eqs. (2.15) all the Would-Be Zero Modes (WBZMs) of all the configurations, but small enough to leave out all the other modes that become irrelevant in the continuum. In general, whether this "golden choice" of M/m s exists or not depends on the amount of explicit chiral symmetry breaking of the fermion discretization. However, in the T = 0 case we can already guess that such a sharp separation can not be present in the spectrum. Indeed, due to the Banks-Casher relation, the spontaneous breaking of chiral symmetry is associated to a proliferation of near-zero modes.
Using for the eigenmodes and the eigenvalues of iD stag [U (2) ] the notation u λ and λ respectively (cf. Eq. (2.10)), we can associate to each mode its (absolute) chirality r λ ≡ |u † λ Γ 5 u λ |. In the continuum r λ = 1 if λ = 0 and r λ = 0 if λ = 0, while on the lattice we generically have 0 < r λ < 1. However, r λ can be used as a figure of merit to identify WBZMs. In Fig. 2 on the right we report a scatter plot of r λ against |λ|/m s for the first 200 low-lying modes of 540 T ≈ 0 configurations (taken every 30 RHMC steps) generated for our finest lattice spacing at this temperature, a ≃ 0.0572 fm. It is evident that a sharp separation between WBZMs and non-chiral modes does not exist. Having seen that the golden strategy is unfeasible, we identified at the finest lattice spacing a range of cut-offs M/m s chosen to reasonably include in the spectral sums all the modes |λ| ≤ M that look "chiral enough"; in the following we will simply refer to this range as the "M -range". Continuum extrapolation is then performed for several values of M/m s chosen within this range, and the residual variability will be taken as a systematic of the extrapolation procedure.
The identification of the M -range has been done for the finest lattice spacing because, being the closest point to the continuum limit, the distinction between chiral and non-chiral mode is clearer. For comparison, in Fig. 2 on the left we also report the scatter plot of r λ against |λ|/m s for our intermediate lattice spacing at this temperature, a ≃ 0.0824 fm.
While in this case the identification of a reasonable range for the cut-off appears difficult, for the finest lattice spacing it is clear that for |λ|/m s ≈ 0.1 there is a change of regime. Therefore, we chose as the M -range the interval [0.05, 0.15].
Once the M -range has been determined, it is possible to extract the continuum limit at fixed value of M/m s . In Fig. 3, we compare the continuum limit obtained using the standard gluonic definition with those obtained through SP for two different values of M/m s . In all cases continuum extrapolation is performed in two different ways: a linear fit in a 2 restricted to the three smallest lattice spacings, and a quadratic fit in a 2 in the whole range, cf. Eq. (2.16). Results obtained varying the fit range and the fit function appear to be in very good agreement in all cases, but we observe that SP continuum extrapolations are less sensitive to the fitting procedure adopted, cf. Fig. 3.
As a matter of fact, if we compare the magnitude of the O(a 2 ) lattice artifacts affecting the two discretizations, we observe that SP estimates have a much faster convergence towards the continuum limit. If we denote by c SP the coefficient of the O(a 2 ) correction to χ for the SP approach (see Eq.  In order to give a final result for χ 1/4 from SP, it is necessary to correctly assess the systematic error related to the choice of M/m s . In Fig. 4 (10) MeV is obtained by choosing a confidence interval that keeps this systematic variation into account, cf. Fig. 4. As for the gluonic result, instead, we estimated the final error from the systematic variation of the extrapolation when changing the fit function and the fit range. These final results are in good agreement among themselves and also with the gluonic determination reported in Ref. [1] and with the ChPT prediction in Eq. (3.1), as shown in Fig. 4. To test the solidity of our results, we performed a further consistency check of our continuum extrapolation procedure, consisting of a common continuum extrapolation of the gluonic results and the SP determinations at fixed M . In this case it is important to include O(a 4 ) corrections to the gluonic determination and, since SP data show milder lattice artifacts with respect to the gluonic ones, the final result turns out to be practically indistinguishable from the SP continuum extrapolation at the same M value. This test constitutes a non-trivial consistency check of the adopted procedure, and we report as our final estimate The same procedure to assess the final error on χ 1/4 will also be applied at finite temperature. Finally, we show in the right plot of Fig 4 the dependence of c SP /c gluo on the cut-off M/m s within the M -range. We observe that the SP discretization is affected by smaller lattice artifacts compared to the gluonic definition, and that corrections to the continuum limit grow as M/m s is increased, getting closer to those of the gluonic discretization. This is due to the fact that, as M/m s grows, the number of irrelevant non-chiral modes, which are more affected by UV cut-off effects, included in the spectral sums grows too.

Topological susceptibility at finite temperature
We computed the topological susceptibility for five temperature values in the high temperature phase of QCD, and specifically for T ≃ 230 MeV, T ≃ 300 MeV, T ≃ 365 MeV, T ≃ 430 MeV and T ≃ 570 MeV. In this section, we will discuss the details of the T = 430 MeV case; a similar analysis has been carried out also for the other temperature values, and the results obtained in these cases are reported in Appendix A.
Simulations at T ≃ 430 MeV have been performed on N 3 s × N t lattices following the same LCP already used in the zero temperature case; simulation parameters are reported in Tab. 2. The spatial extent of the lattice was chosen to ensure an aspect ratio not smaller than 3.   Simulations at finite temperature, unlike those at T = 0, have been performed adopting the multicanonical algorithm. Our implementation of the multicanonical algorithm closely follows the one already adopted in Ref. [5]. The quantity Q mc entering the topological potential is the clover discretization (2.4) of the topological charge computed after n mc levels of stout smearing, which allows to adopt the RHMC algorithm also in the presence of the topological potential. The values of n mc used ranged from n mc = 20 for the coarsest lattice spacing to n mc = 10 for the finest one, while the isotropic smearing parameter was always fixed to ρ mc = 0.1. The functional form of the topological potential adopted was [5]: and in our simulations, Q max = 3 turned out to be sufficient to observe a dramatic enhancement in the number of fluctuations of Q gluo . The free parameters B and C have instead been tuned through short preliminary runs by requiring the Monte Carlo histories of the measured topological charge Q gluo to uniformly explore the interval [−Q max , Q max ]. The improvement obtained with the multicanonic algorithm is exemplified in Fig. 5. Without any bias potential, Q gluo assumes a non-zero value only a handful of times, while much more fluctuations are observed as the topological potential is switched on 3 . This allows to dramatically improve the accuracy with which χ can be computed with a given machine time budget. In order to compute the topological susceptibility, we follow the same procedure already discussed for the T ≃ 0 case in Sec. 3.1. The first step consists of determining a reasonable interval for M/m s by studying the scatter plot of the chiralities for the finest lattice spacing available, which is shown in Fig. 6 on the right. With respect to the T ≃ 0 case, at finite temperature the separation between high and low-chirality modes is more evident, as two almost disconnected clusters can be observed with r λ ∼ 0.7 − 0.8 and r λ 0.5. However, also in this case, the isolation of WBZMs is an ambiguous task, as we also observe a sparsely-populated group of modes in between, with 0.5 r λ 0.7.
Such an ambiguity is not a characteristic of the whole configuration ensemble, but already arises configuration by configuration. In Fig. 7, we plot the histogram of the separation between candidate WBZMs and Non-Zero Modes (NZMs) for the same sample represented in the scatter plot of Fig. 6 on the right.  The candidate WBZMs were identified from the value of the gluonic topological charge: more precisely we assumed the first n t |Q gluo | lowest-lying eigenmodes of iD stag to be the WBZMs (the factor n t = 4 takes into account the four tastes of staggered fermions). This strategy is the same adopted in Ref. [3] to identify WBZMs for their reweighting procedure, and it is justified on the basis of the index theorem for staggered fermions in the continuum n t Q = n + − n − , supplemented by the assumption that n − = 0 (n + = 0) if Q > 0 (Q < 0), which should be approximately true if the lattice volume is not too large. However, we observe that the typical relative separation between WBZMs and NZMs is of the order of 10 −2 ÷ 10 −1 , cf. Fig. 7, meaning that a sharp separation between WBZMs and NZMs cannot be unambiguously established already at the level of the single configuration.
Thus, we cautiously choose our cut-off masses in the interval M/m s ∈ [0. 3,5] in order to include in our spectral sums all modes with |λ| ≤ M and r λ 0.5. Again, the M -range has been identified for the finest lattice spacing available at this temperature, a ≃ 0.0286 fm. However, we observe that at finite temperature also intermediate lattice spacings would lead approximately to the same choice. For comparison, in Fig. 6 on the left, we also report the scatter plot of r λ against |λ|/m s for an intermediate lattice spacing at this temperature, a ≃ 0.0381 fm.
In Fig. 8 we compare results for χ 1/4 SP at a ≃ 0.0572 fm obtained, respectively, by computing the multiplicative renormalization constant Z (stag) SP in Eq. (2.13) from our finite temperature ensemble generated on a 32 3 × 8 lattice and from a corresponding zero temperature ensemble on a 32 4 lattice. We recall that, while the computation of the spectral traces appearing in the definition of Q is computed on our finite temperature ensemble (on a 32 3 × 8 lattice) or on a corresponding zero-temperature ensemble (on a 32 4 lattice). The quantity Z (stag) SP was computed in both cases using the staggered spectrum obtained imposing periodic boundaries along the temporal direction for D stag .
. As Fig. 8 shows, both ways of computing χ 1/4 SP give perfectly consistent results, in agreement with results obtained at finite temperature in the quenched theory [32]. For this reason, for all lattice spacings we computed both Z (stag) SP and Q (stag) SP from our finite temperature ensembles, where the staggered spectrum entering the spectral definition of these quantities was obtained, respectively, imposing periodic and anti-periodic boundaries along the time direction for D stag .
In Fig. 9, we compare the continuum limits of χ 1/4 obtained via the gluonic and the SP definitions, using the same fit function in Eq. (2.16) already adopted for the T = 0 case. For the SP case, we show results for the values M/m s = 0.3 and 0.5. As for the zero-temperature case, also at finite temperature we observe that the SP definition allows to achieve a sensible reduction of the magnitude of O(a 2 ) corrections with respect to the gluonic case, c SP (0.3)/c gluo ∼ 5 × 10 −2 and c SP (0.5)/c gluo ∼ 10 −1 , allowing for a better control over systematics related to the continuum extrapolation, as we observe a very good agreement in the obtained extrapolations when varying the fit range and the fit function. Moreover, we observe a good agreement among the gluonic and the SP determinations, cf. Fig. 9. obtained for this temperature in, respectively, Refs. [3,4]. The former was temperature-interpolated according to the DIGA prediction χ 1/4 ∼ T −2 and the isospin-breaking factor was removed. The latter was mass-extrapolated according to χ 1/4 ∼ m π .
Finally, in Fig. 10, we show how the continuum limit of χ SP varies as a function of M/m s within the M -range interval determined before (left plot). Also at finite temperature we observe that the continuum limit of χ 1/4 SP is quite stable within the M -range; nevertheless the residual systematic is incorporated in our final result, also shown in Fig. 10. As for the gluonic case, the final error was instead estimated from the systematic variation of the continuum extrapolation varying the fit function and the fit range, similarly to the procedure followed at zero temperature. Finally, we also checked that performing a common fit to our gluonic results and our SP determinations for a fixed value of M/m s gave consistent results with the extrapolations performed for the SP determinations alone at the same value of M/m s . As already done at zero temperature, we consider the conservative estimate χ 1/4 SP = 15(5) MeV as our final determination, where this error takes into account both the systematic and the statistical sources of uncertainty.
In Fig. 10, in the right plot, we also report how the ratio c SP /c gluo depends on M/m s . As already observed when discussing the T = 0 results, when M/m s is chosen small enough, the SP discretization suffers for much smaller lattice artifacts compared to the gluonic one. As this cut-off increases, the ratio c SP /c gluo tends to grow, approaching 1, which can be interpreted as the effect of the inclusion of more and more non-chiral modes in the spectral sums. We now want to compare our final SP result for χ 1/4 at T ≃ 430 MeV, χ 1/4 = 15(5) MeV, with other determinations at this temperature. This value is compatible with our gluonic estimation (χ 1/4 = 14(5) MeV) and is also compatible within ∼ 1.2 σ with the value that can be obtained interpolating results of Ref. [3] in Tab. S9 according to the DIGA prediction χ 1/4 ∼ T −2 and removing the isospin-breaking factor 0.88 the authors use to restore the u−d mass difference in their N f = 2+1+1 results for the susceptibility: χ 1/4 = 9.0(4) MeV. Our determination is also in agreement within the errors with the gluonic result reported in Ref. [4] for this temperature: χ 1/4 = 11(2) MeV. Note however that in this case the comparison is less strict as the latter result was obtained rescaling gluonic determinations for χ 1/4 reported in Ref. [4] by (135 MeV)/m π , as they were obtained for m π ≈ 160 MeV. Of course this procedure assumes that χ depends on the light quark mass as predicted by the DIGA: χ 1/4 ∼ m 1/2 l ∼ m π . Finally, we mention that in Ref. [4] the authors also computed χ 1/4 adopting a fermionic definition based on the disconnected chiral susceptibility, which was however fully consistent with the gluonic results used above.

χ(T ) from spectral projectors and comparison with the DIGA
In our N f = 2 + 1 setup, semiclassical and perturbative approximations predict for the topological susceptibility when T ≫ T c the scaling (see Eq. (1.2)) The aim of the present section is to study the behavior of our results for χ 1/4 as a function of T and to compare it with the DIGA prediction. Our final results for χ 1/4 as a function of T /T c , obtained both from the gluonic and the SP definitions, are reported in Tab. 3. These results were obtained following exactly the same strategy outlined in Sec. 3.2 for T ≃ 430 MeV, and more details about them can be found in Appendix A.  (5)  570 3.68 8(6) 6(6) Table 3: Results for the fourth root of the topological susceptibility as a function of T . For the crossover temperature T c we adopted the reference value T c = 155 MeV.
In Fig. 11, we show the behavior of χ 1/4 as a function of T /T c , from which it should be clear that our data are compatible with a power-law behavior of the type 3.4 in the whole explored range. If we perform a best fit of the data using the function when including in the fit all available points, we obtain the results b SP = 1.82(43), b gluo = 1.67 (51) .
If we instead exclude from the fit the lowest temperature T = 230 MeV, we get: We thus conclude that our data are in perfect agreement with a power-law behavior in the whole explored range, with an effective exponent that turns out to be well compatible with b DIGA within our errors already for T 300 MeV, i.e., already for T /T c 2. Actually, also when T = 230 MeV is included in the fit the exponent b SP turns out to be compatible within 1 standard deviation with b DIGA . However, as it can be clearly seen from Figs. 11, the inclusion/exclusion of such point visibly changes the slope of the fit. This may be an indication, albeit at present not conclusive, that the effective exponent b changes when going from T ∼ 200 MeV to T ∼ 300 MeV.
In Fig. 11 we also compare our results with previous determinations reported in Refs. [3,4]. Concerning results of Ref. [3], to make a fair comparison, we removed the isospinbreaking factor the authors use to restore the u − d mass difference in their N f = 2 + 1 + 1 results. Our spectral determinations are systematically larger than those of Ref. [3], and we find a ∼ 2.5 − 3 standard deviation tension at T = 300 MeV and T = 365 MeV. As for the temperature behavior, we observe that the best fit of the data of Ref. [3] for T 170 MeV with the fit function (3.5) yields b = 1.945(23), i.e., in this case the DIGA-like power-law seems to set in for smaller values of the temperature.  Figure 11: Behavior of χ 1/4 as a function of T /T c in double-log scale for our SP and gluonic data. Gluonic points are slightly shifted to improve readability. Starred points represent results taken from Ref. [3], where the isospin-breaking factor was removed, while the uniform-shaded area represents the gluonic determinations reported in Ref. [4], which have been mass-rescaled according to χ 1/4 ∼ m π . The left plot reports the best fits of our data according to the fit function (3.5) and performed including all available points:χ 2 /dof = 0.54/3 for gluonic data andχ 2 /dof = 2.1/3 for SP ones. The right plot reports the result of the same fits, but performed excluding our lowest temperature T ≃ 230 MeV (excluded points are empty). In this case,χ 2 /dof = 0.06/2 for gluonic data andχ 2 /dof = 0.22/2 for SP ones. Dashed, solid and dotted lines represent, respectively, best fits of our gluonic data, of our spectral data and of gluonic data of Ref. [3].
Concerning the gluonic results of Ref. [4], since they were obtained for m π = 160 MeV, we rescaled them according to the DIGA prediction χ 1/4 ∼ m π in order to make a fair comparison. Also in this case we observe that our determinations lay systematically above. Although this comparison is less strict because of the different pion masses adopted, also in this case we observe a ∼ 2 − 2.5 standard deviation tension among determinations for T = 300 MeV and T = 365 MeV.
The observed tensions between our spectral determinations and previous results reported in the literature deserve to be further investigated, for example by refining the present errors on the spectral determinations for T 400 MeV. Also probing higher temperatures with our methods would be interesting in order to extend the present comparison towards the ∼ 1 GeV region, which is also interesting in the context of axion cosmology.

Conclusions
In this work we performed a numerical lattice study of the behavior of the topological susceptibility χ(T ) in the high temperature regime of QCD with N f = 2 + 1 quarks at the physical point.
Our computational strategy relies on the discretization of the topological charge through Spectral Projectors on the eigenmodes of the staggered Dirac operator. The reason for this choice is to reduce the large lattice artifacts affecting the gluonic definition of χ at high temperatures when non-chiral quarks, like the staggered ones, are adopted to discretize the QCD action. The problem of the dominance of the Q = 0 sector, which is related to the suppression of χ at high-T , is instead addressed by the use of the multicanonical method already applied for this purpose in the recent work [5].
The spectral definition of the topological susceptibility introduces a new free parameter, as the sum over the chiralities of the eigenmodes of the lattice staggered Dirac operator is performed by including all eigenvalues with magnitude |λ| ≤ M . In principle, any value of M (kept constant in physical units as a → 0) would provide a correct continuum limit of χ; however, since isolating WBZMs is highly ambiguous, we identified a reasonable range of M -values on which to perform the continuum limit at fixed M . Any residual systematic related to the choice of M is then assessed afterwards and included in our final determination of the uncertainty affecting χ.
To test the spectral method we first of all study the T = 0 case, where χ can be reliably computed by using chiral perturbation theory. In this case the multicanonical algorithm is not needed, since the MC evolution naturally visits Q = 0 sectors, and the final result obtained by the spectral method is perfectly compatible both with the NLO ChPT result and with previous gluonic estimates.
The advantage of the spectral projectors method with respect to the gluonic one is that lattice artifacts are in this case much smaller, and the extrapolation towards the continuum limit is thus better under control. Moreover, due to the presence of the parameter M , in the spectral projectors setting we have a natural procedure to check for the presence of systematics in the continuum extrapolation procedure: to look for residual dependence on M of the extrapolated result. As a matter of fact, this systematic gives the dominant contribution to the final error on χ for some of the temperatures studied in this work.
In the high temperature regime, we explored 5 values of T , ranging from ∼ 200 MeV to ∼ 600 MeV. Also in these cases we observe good agreement between SP and gluonic data, with the spectral results that are generically more accurate than the gluonic ones; however, more important is that, as noted above, in the SP case we have a much better control of the systematics of the continuum extrapolation.
We finally investigated the behavior of our results for χ as a function of T , comparing with expectations based on the DIGA approximation. 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. This is in agreement with the results obtained in Ref. [4] and is consistent with a growing number of observations suggesting that high-temperature QCD is dominated by strong non-perturbative effects for temperatures going approximately from the chiral crossover up to ∼ 300 MeV [66][67][68][69].
We remark that our results for χ 1/4 from spectral projectors show a ∼ 2 − 3 standard deviation tension in the range 300 MeV T 400 MeV range when compared with previous determinations in Refs. [3,4]. Moreover, also when the consistency between our results and the ones in the literature is better, our spectral determinations for χ 1/4 (T ) systematically points to larger values in the whole explored range. The same behavior, when observed in the gluonic determinations of χ, could be ascribed to a problem of the continuum extrapolation, that is unable to capture the asymptotic O(a 2 ) scaling and intro-duces 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 observed disagreement seems unlikely in this case. In conclusion, the picture emerging from the comparison carried out in Fig. 11 is that a complete quantitative understanding of the behavior of χ(T ) in the high temperature regime of QCD is still missing, and further studies will be required to clarify the sources of the observed tensions between different determinations.
Several directions can be followed to improve the present study: first of all it 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, and obtain a precise and unbiased estimate of χ(T ) in the explored temperature range. For this purpose simulations with larger statistics and smaller lattice spacings are required. Other natural extensions include a systematic study of the temperature ranges T 400 MeV, where deviations from the DIGA power law should be visible, and T ∼ 1 GeV, in order to reach the region of temperatures directly relevant for axion cosmology.
However, simulations at smaller lattice spacings (needed both to improve the present estimates and to reach higher temperatures) are practically unfeasible with standard RHMC simulations, due to the severe topological critical slowing down. A promising strategy to overcome this problem is the parallel tempering on boundary conditions algorithm proposed in Ref. [70], which has already been successfully applied both to two dimensional models [70,71] and to 4d SU(N ) Yang-Mills theories without matter fields [72,73].
• scatter plot of the chiralities of the lowest eigenmodes of the staggered operator for the finest lattice spacing explored at that temperature; • comparison of the continuum limits obtained for the gluonic and the SP definitions; • assessment of the systematics related to the choice of M/m s .
The parameters of all the performed simulations are summarized in Tab. 4. For the two coarsest lattice spacings simulated at T ≃ 230 MeV, no enhancement of topological fluctuations was needed. Therefore, simulations of these points were carried on without the multicanonic algorithm.    [3,4]. The former was temperature-interpolated according to the DIGA prediction χ 1/4 ∼ T −2 and the isospin-breaking factor was removed. The latter was massextrapolated according to χ 1/4 ∼ m π . Right: c SP /c gluo for the several values of M/m s within the M -range. A straight horizontal line is set at 0.  [3,4]. The former was temperature-interpolated according to the DIGA prediction χ 1/4 ∼ T −2 and the isospin-breaking factor was removed. The latter was massextrapolated according to χ 1/4 ∼ m π . Right: c SP /c gluo for the several values of M/m s within the M -range. A straight horizontal line is set at 0.

A.4 T = 570 MeV
The shown M -range for χ SP is narrower compared to the one shown in the scatter plot of chiralities, cf. Figs. 18 and 19, because the mean maximum eigenvalue of D stag with periodic boundaries along the time direction was smaller than the anti-periodic case. Vertically-hatched band displays the value of χ 1/4 obtained for this temperature in Ref. [3], which was temperature-interpolated according to the DIGA prediction χ 1/4 ∼ T −2 and the isospin-breaking factor was removed.   [3], which was temperature-interpolated according to the DIGA prediction χ 1/4 ∼ T −2 and the isospinbreaking factor was removed. Right: c SP /c gluo for the several values of M/m s within the M -range. A straight horizontal line is set at 0.