The topological susceptibility in finite temperature QCD and axion cosmology

We study the topological susceptibility in 2+1 flavor QCD above the chiral crossover transition temperature using Highly Improved Staggered Quark action and several lattice spacings, corresponding to temporal extent of the lattice, $N_\tau=6,8,10$ and $12$. We observe very distinct temperature dependences of the topological susceptibility in the ranges above and below $250$ MeV. While for temperatures above $250$ MeV, the dependence is found to be consistent with dilute instanton gas approximation, at lower temperatures the fall-off of topological susceptibility is milder. We discuss the consequence of our results for cosmology wherein we estimate the bounds on the axion decay constant and the oscillation temperature if indeed the QCD axion is a possible dark matter candidate.


Introduction
The visible matter consisting of protons and neutrons derive their masses due to strong interactions described by Quantum Chromodynamics (QCD). Advances in non-perturbative techniques, mainly Lattice gauge theory has lead to unraveling of the phase diagram of QCD at high temperatures and first principles determination of hadron masses and decay constants ( see e.g. Refs. [1,2,3] for recent reviews). One aspect of QCD which still lacks a comprehensive understanding is its intricate vacuum structure.
The QCD vacuum is believed to consist of energy degenerate but topologically distinct minima which are separated by potential barriers of height Λ QCD /α s , where α s is the coupling constant for strong interactions [4,5]. Each minima is characterized by a distinct integer or Chern-Simons number which represent the third homotopy class of the mapping between the Euclidean spacetime R 4 and the gauge manifold of SU (3). Tunneling between adjacent minima results in change of Chern-Simons number by unity. Such solutions are called instantons. They minimize the classical action density and are topologically protected [6]. At finite temperature, much richer topological structures are expected to arise. Due to periodicity along the Euclidean time-like direction, multi-instanton solutions can exist [7] which are called the calorons [8,9]. Calorons are characterized by a definite topological charge ±1 as well as the holonomy, which is the the Polyakov loop at spatial infinity. In presence of a non-trivial holonomy, the caloron solutions in SU (N c ) gauge theory can be decomposed into N c monopoles carrying both electric and magnetic charges, called the dyons [8,10]. On the other hand at temperatures much higher than T c , Debye screening would only allow instantons of radius ≤ 1/T to exist and the QCD partition function written in terms of non-interacting instantons, is well defined and does not suffer from infrared singularities [11]. This is the dilute instanton gas approximation (DIGA).
The outstanding issues are to understand the mechanism of confinement and chiral symmetry breaking in terms of these constituents. Some aspects of chiral symmetry breaking can be understood within the so called interacting instanton liquid model [12], but it cannot explain confinement. There are hints that dyons could be a promising missing link in our understanding of confinement mechanism [10] and preliminary signals for the existence of dyons have been reported [13] from lattice studies. The other important problem is to understand from first-principles whether DIGA can indeed describe QCD vacuum at high temperatures. Lattice studies in pure SU (2) and SU (3) gauge theories have observed the onset of a dilute instanton gas at ∼ 2 T d , T d being the deconfinement temperature [14]. However in QCD, the situation is more richer as index theorem ensures that instanton contributions to the partition function are suppressed due to the fermion determinant. Presence of light fermions thus may lead to much stronger interactions between the instantons or its dyonic constituents. Due to conceptual and computational complexity of lattice fermion formulations with exact chiral symmetry, [15,16] the QCD vacuum at finite temperature in presence of dynamical fermions have been studied in de-tail only recently [17,18,19,20]. Though some of the lattice studies on the low-lying eigenvalue spectrum of the Dirac operator with nearly perfect chiral fermions report the onset of dilute instanton gas picture in QCD already near T c [19], a more recent study on larger lattice volumes and nearly physical quark masses, observe onset of the dilute instanton gas scenario at a higher temperature ∼ 1.5 T c [20]. Here T c denotes the chiral crossover transition temperature, T c = 154(9) MeV [21]. It is thus important to study a different observable which is sensitive to the topological structures in QCD, to address when can the QCD vacuum be described as a dilute gas of instantons and whether dyons indeed exist near T c . One such observable is the topological susceptibility.
Understanding the topological susceptibility, or more generally the structure of QCD vacuum at T > 0 is not only of great conceptual importance but has applications in cosmology as well. One of the topics that has received renewed interest is whether the axions could be a possible dark matter candidate. Strong bounds exist for the magnitude of QCD θ angle that breaks CP symmetry explicitly and measurement of neutron electric-dipole moment has set the value of θ < 10 −10 . One possible way to generate such a small value of the QCD θ angle is through the Peccei-Quinn (PQ) mechanism [22,23]. This involves including a dynamical field a(x) called axion, which has a global U (1) symmetry and couples to the gauge fields in the QCD Lagrangian as, Due to dynamical symmetry breaking of this anomalous U (1) PQ symmetry, the axion field a(x) acquires a vacuum expectation value at the global minima given by θ + a /f = 0. The curvature of the axion potential or equivalently the free energy, near the global minima is quantified by the topological susceptibility of QCD, thus generating a mass m a for the axions. The mass of axions is of central importance [24], if indeed the axion is a possible cold dark matter candidate [25]. At very high temperatures or at large scales when inflation sets in, the axion field would have possibly taken random values at causally disconnected spacetime coordinates. As the universe expands and cools down the axion field feels the tilt in the potential and slowly rolls down towards the global minima. At a temperature T osc , where the Hubble scale H(T osc ) = m a (T osc )/3, the axion field would start oscillating about its minima. If the amplitude of the oscillations are large, damped only due to the Hubble expansion, these largely coherent fields would be a candidate for cold dark matter [26,27,28]. Hence the abundance of axion dark matter is sensitive to the mass of the axion and thus to the topological susceptibility in QCD at the scale T ∼ T osc . The typical value for T osc ∼ GeV and it is important to understand the temperature dependence of the topological susceptibility at high enough temperatures. The topological susceptibility, χ t has been studied on the lattice both for cold and hot QCD medium. At T = 0, its continuum extrapolated value is known in quenched QCD [29,30] as well with dynamical fermions [31]. Motivated from DIGA [11], the χ t can have a power law dependence on temperature which is characterized by an exponent γ, such that χ t ∼ T −γ . Including quantum fluctuations to the leading order in α s over classical instanton action, the exponent γ = 11/3N c − 5 + 2N f /3 for QCD with N f light quark flavors [32]. In quenched QCD, high statistics study on large volume lattices for temperatures upto 2.5 T d , employing cooling techniques, has reported a γ = 5.64(4) [33] which is significantly different from a dilute instanton gas scenario. A more extensive study [34] on a larger temperature range upto 4 T d using controlled smearing through the Wilson flow, has found exponent of χ t to be γ = 7.1(4)(2) in agreement with DIGA. However a prefactor of O(10) is needed to match the one-loop calculation of χ t within DIGA to the lattice data. Extension to QCD with realistic quark masses has been discussed [35] and very recently performed [36]. It was found in Ref. [36] that the temperature dependence of χ t is characterized by γ ∼ 2.90(65), in clear disagreement with DIGA, which predicts γ = 8 for N f = 3. This hints to the fact that the χ t in QCD with realistic quark masses could have a very different temperature dependence than in pure gauge theory.
In order to have a conclusive statement about the temperature dependence of χ t in QCD with dynamical physical quarks, we approach the problem in a two fold way. We use a very large ensemble of QCD configurations on the lattice with Highly Improved staggered Quark (HISQ) discretization and measure χ t from the winding number of gauge fields after smearing out the ultraviolet fluctuations. For the staggered quarks, the index theorem is not uniquely defined hence we carefully perform a continuum extrapolation of our measured χ t for very fine lattice spacings. We also show the the robustness of our findings and indirectly study the validity of the index theorem for staggered quarks, as we approach the continuum limit, by comparing the results for χ t with those obtained from an observable constructed out of fermionic operators, namely the disconnected part of the chiral susceptibility.
We also use our results for χ t as a tool to estimate the nature and abundance of the topological structures in QCD as a function of temperature. We show for the first time that the effective exponent γ governing the fall-off of χ t changes around 1.5 T c . The change in the exponent is very robust and persists as we go to finer lattice spacings. This may indicate that the nature of topological objects contributing to χ t is different for T < 1.5 T c and T > 1.5 T c , further supporting the conclusions of Ref. [20] about the onset of diluteness in the instanton ensemble. Moreover it points to the fact that fluctuations in local topology at temperatures close to chiral crossover transition could additionally be due to topological objects like the dyons.
The brief outline of the Letter is as follows. In the next section we discuss the details of our numerical calculations. Following this are our key results. The continuum extrapolated result for χ t is shown and compared with the disconnected chiral susceptibility. The importance of performing correct continuum extrapolation to get reliable estimates is discussed in detail. Finally we compare our results with that of χ t calculated within DIGA considering quantum fluctuations up to two-loop order. We show that the dependence of χ t in QCD on the renormalization scale is mild and discuss its consequences for axion phe-nomenology. We conclude with an outlook and possible extensions of this work.

Lattice Setup
We use in this work, the 2 + 1 flavor gauge configurations generated by the HotQCD collaboration for determining the QCD Equation of State [37]. These were generated using the HISQ discretization for fermions which is known to have the smallest taste symmetry breaking effects among the different known staggered discretizations. The strange quark mass m s , is physical and the light quark mass is m s /20, which corresponds to a Goldstone pion mass of 160 MeV. The temporal extent of the lattices considered in this work are N τ = 6, 8, 10, 12 and the corresponding spatial extents N s = 4N τ , such that we are within the thermodynamic limit. For the coarsest N τ = 6 lattice, we consider a temperature range between 150-800 MeV and for the finest N τ = 12 lattice, we study between 150-407 MeV. We perform gradient flow [38] on the gauge ensembles to measure χ t . Typically we have analyzed about 500 configurations near T c and to a maximum of 13000 configurations at the higher temperatures and finer lattices. We usually perform measurements on gauge configurations separated by 10 RHMC trajectories to avoid autocorrelation effects. The number of configurations analyzed for each β = 10/g 2 value are summarized in table 1.

Details of our smoothening procedure
In order to calculate the the topological charge and its fluctuations it is necessary to get rid of the ultraviolet fluctuations of the gauge fields on the lattice. We employ gradient flow to do so. The gradient flow is defined by the differential equation of motion of V t , which using same notations as in Ref. [38], where g 0 is the bare gauge coupling and S[V t ] is the Yang-Mills action. The field variables V t (x = (x, τ ), µ) are defined on the four dimensional lattice and satisfy the condition (2) has the form of a diffusion equation, the gradient flow results in the smearing of the original field U µ (x) at a length scale f = √ 8t. However, unlike the common smearing procedures (see e.g. Refs. [41,42]), which can be applied only iteratively in discrete smearing steps, the amount of smearing in the gradient flow can be increased in infinitesimally small steps. This ensures that the length scale over which the gauge fields are smeared can be appropriately adjusted as we vary the lattice spacing and/or temperature. We use a step-size of 0.01 in lattice units for flow time t, which enables us to perform sufficiently fine smearing. At T = 0, the flow time has to be sufficiently large to get rid of ultraviolet noise but typically has to be smaller than 1/Λ 2 QCD . For non-zero temperature, additionally, one has to choose the flow time such that f is smaller than 1/T [43,44,45]. In this study we use a maximum of 200 flow steps in the calculation of topological charge. This corresponds to the value of f T = 0.333 for our finest N τ = 12 lattices. Even for the coarsest N τ = 6 lattices we ensure that f T < 1.
At each flow step, the topological charge Q was measured using the clover improved definition of F µ,νF µ,ν . Fractional values of Q which are within 20% of the next higher integer were approximated to that integer. We note that increasing flow time beyond f t = 0.1, the value of χ t first increases with increasing flow time and then stabilizes. This is because on the smoothened gauge configurations, χ t is rounded to the nearest non-zero integer more frequently. Since the flow time is never too large i.e. f T < 1, instantons are not washed due to smearing We determined the optimal number of flow steps at each temperature from the plateau of the χ t . For T < 300 MeV, the plateau in χ 1/4 t was already visible within 50 − 100 steps of Symanzik flow. At higher temperatures, more than 100 steps in flow time are required to obtain a stable plateau in χ 1/4 t . We compare the onset of the plateau-like behavior for χ 1/4 t as a function of the number step sizes for different lattice spacings in Fig. 1. From these data, we infer that 200 steps of smearing is sufficient to give a stable value of the χ 1/4 t . We also checked that for a finer N τ = 10 lattice, increasing the number of smearing steps to 300 did not change the value of χ 1/4 t . This is shown in the right panel of  In order to ensure that all the topological sectors have been adequately sampled, we measure the topological charge distribution for different topological sectors, labeled by Q. In Fig. 2 we show the histograms for the distribution of Q for our given configurations at different temperatures, upto ∼ 2.8 T c . It is evident that different topological sectors have been sampled effectively since the Q distribution is a Gaussian peaked at zero at both the temperatures. For β = 7.825, which corresponds to a temperature of ∼ 407 MeV for the N τ = 12 lattice, we still observe tunneling between different topological sectors even though the lattice spacing is fine enough, a = 0.041 fm and lattice volume is large, V ∼ (2.5 fm) 3 . However, we do not observe tunneling among the topological sectors for T > 600 MeV for N τ ≥ 8, even when we sample about O(1000) gauge configurations.

Topological susceptibility from gradient flow
Our results for the topological susceptibility χ t are summarized in Fig. 3. In most of our plots we rescale the temperature axis by the chiral crossover transition temperature T c = 154 MeV, unless stated otherwise. The cut-off effects are reflected in the strong N τ dependence of the data up to the highest temperatures considered. At low temperatures, large cutoff effects in χ t can be understood in terms of breaking of the taste symmetry of staggered fermions and can be quantified within the staggered chiral perturbation theory [3,46]. However taste breaking effects are expected to be milder at higher temperatures since the lattice spacing is finer. Our results indicate that above T c , cutoff effects depend dominantly on aT = 1/N τ rather than ∼ aΛ QCD . In fact, the cutoff effects in χ t is much stronger than for any other thermodynamic quantity calculated so far with HISQ action [37,47,48]. It is also evident from our data that χ t has an interesting temperature dependence. If indeed χ t is characterized by a power law fall off such that χ 1/4 t = AT −b , then a fit to our data restricted to intervals [165 : 220] MeV and [220 : 600] MeV resulted in b ∼ 0.9-1.2 and b ∼ 1.35-1.6 respectively. In fact, it turns out that it is impossible to obtain an acceptable fit to the data with a single exponent b in the entire temperature range. This led us to model the variation of the exponent by giving it two different values for temperatures T > 1.5T c and T < 1.5T c . This is also taken into account when performing continuum extrapolation. We note that at leading order, the exponent calculated within DIGA for N f = 3 corresponds to b = γ/4 = 2. Higher order corrections would effectively reduce the value of b.   Let us compare our results with those obtained in QCD using different fermion discretizations [19,36,49]. We do not expect the values of χ t for different choice of fermion discretizations to agree with each other at non-vanishing lattice spacings. However such comparisons provide insights on the nature and severity of the discretization effects. Another source of this difference could be due to the fact that all earlier calculations have been performed for different values of quark masses. For a fair comparison with the earlier results, we scale χ t with m 2 l ∼ m 4 π . In Fig. 3 we show our results for two lattice spacings a = 0.082 fm and a = 0.06 fm and compare them with other recent results for χ t . Our results for a = 0.06 fm agree well with the ones obtained using twisted mass Wilson fermions for pion mass 370 MeV and very similar lattice spacing once rescaling with the pion mass is performed [49]. Comparing our results with the latest calculations performed using staggered fermions with the so-called stout action [36] and physical pion mass, we find that our results are ∼ 70% of the values obtained in Ref. [36] at similar lattice spacings. This could be due to the fact that discretization errors for the HISQ action considered in our work are smaller than for the stout action. Domain wall fermions give the smallest value of χ t even on relatively coarse lattice, a = 1/(8T ), which is not surprising since the chiral symmetry in the corresponding calculation is almost exact.

Disconnected chiral condensate and topological susceptibility
Low-lying eigenvalues of the QCD Dirac operator are known to be sensitive to topological objects of the underlying gauge configurations. Calculating fermionic observables which are sensitive to these low-lying eigenvalues provides an alternative way to study the topological fluctuations in QCD. One such observable is the disconnected chiral susceptibility. QCD action with two light quark flavors has a SU A (2) × SU V (2) × U B (1) × U A (1) symmetry which is effectively broken to SU V (2) × U B (1) at T c . The effective restoration of these symmetries can be studied through the volume integrated meson correlation functions or susceptibilities, The above quantities can be expressed in terms of quark fields which in the one-flavor normalization [19] are given as, σ = ψ l ψ l , δ i = ψ l τ i ψ l , η = iψ l γ 5 ψ l , π i = iψ l τ i γ 5 ψ l . The susceptibilities in the different quantum number channels are related as, where the χ π,δ only contain the connected part of the correlators. When chiral symmetry is effectively restored in QCD, χ σ = χ π and χ η = χ δ . This implies that the difference in the integrated correlator χ π − χ δ = 2χ disc . Effective restoration of U A (1) would imply χ π = χ δ , which results in χ disc = 0. Therefore χ disc can be considered as the measure for U A (1) breaking. Expressing the topological charge as the trace of the Dirac fermion matrix for light quarks [19], In the chiral symmetry restored phase, χ t = (m l ) 2 χ disc since χ 5,disc = χ disc . In other words χ disc is a proxy for χ t . The above considerations are valid on the lattice only near the continuum limit. In order to check how well the corresponding equalities are satisfied at finite lattice spacing, we first consider the quantity χ π − χ δ normalized by T 2 . χ π was estimated from the Chiral Ward identity χ π = ψ ψ l /m l , and χ δ from the connected part of meson correlators χ conn from Ref. [21] (Tables X-XII) and converting the numerical results to single flavor normalization. The results for χ π − χ δ for our QCD ensembles with HISQ discretization are shown in Fig. 4. It is evident from this plot that for T 1.2T c the quantity χ π − χ δ = χ disc for all different lattice spacings considered therefore χ t can be independently estimated from the disconnected 2-point chiral susceptibility χ disc . This estimate of χ t is shown in Fig. 4. We used the combined data on χ disc from Refs. [21] and [37], since the later has high statistics results for a larger temperature range. A closer inspection of the lattice data on χ disc in Ref. [37] reveals that the magnitude of the errors on each data point increases with increasing temperatures but at some higher temperature, drops abruptly. We interpret this fact to the freezing of topological charge at high temperatures which results in unreliable estimates for χ disc . This is consistent with our measurement of χ t from Symanzik flow, where we indeed observe a signal for χ t beyond T = 500 MeV for N τ = 6, 8 lattices and beyond T > 400 MeV for N τ = 10, 12 lattices but with larger systematic and statistical uncertainties. Therefore, the corresponding data points (shown as open symbols in left panel of Fig. 3) are not included in the continuum extrapolation discussed in the following section. The significant N τ dependence of the results for χ disc evident from Fig. 4 further validates our observation in the previous section that the dominant cutoff-effects for χ t are controlled by aT = 1/N τ rather than the value of the lattice spacing in some physical units like Λ QCD . On finite lattices, the values of m 2 l χ disc are considerably smaller than the values of χ t defined through gauge observables. Furthermore, the approach to the continuum for m 2 l χ disc is opposite to χ t , determined from the fluctuations of winding number of the gauge fields. This significant difference between the values of χ t measured in terms of gluonic and fermionic observables has been observed earlier [19] using Domain wall fermions, which has nearly perfect chiral symmetry on a finite lattice. Therefore this difference is unlikely due to tastebreaking effects inherent in staggered quark discretization. In the next section we show a reliable method for the continuum extrapolation of these observables even with such large lattice artifacts.

Continuum Extrapolation of topological susceptibility
In this section we discuss in detail our procedure for performing continuum extrapolation of χ 1/4 t . As noted earlier, the exponent controlling the power law fall-off of χ 1/4 t on temperature depends on the chosen temperature interval. We performed simultaneous fits to the lattice data obtained for different N τ with the following ansatz, As mentioned in the previous section, data at high temperatures, where statistical errors may not have been estimated correctly are excluded from the fit.
One of the important issues that needs to be clarified is whether two different estimates of χ t obtained through the gluonic operator definition for Q and through m 2 l χ disc respectively, agree in the continuum limit. We first performed fits to our numerical results of χ 1/4 t presented in section 3.1 according to Eq. 7. We considered two independent fits for the two temperature intervals 165 MeV ≤ T < 270 MeV and 240 MeV ≤ T ≤ 504 MeV, respectively. We set b 4 = b 6 = 0 and to ensure to be within the scaling regime, we omitted the N τ = 6 data. This did not change the continuum value of χ t within errors. Furthermore, when we performed fits setting a 4 = 0, the continuum extrapolated value was 15% − 25% larger. This implies that the term proportional to a 4 is important but subdominant relative to the leading term proportional to a 2 . In the first temperature region we obtain b = 1.496(73), while the corresponding value from fit to the higher temperature interval gave b = 1.85 (15). We found that continuum extrapolations with Eq. 7 in the entire temperature interval is not possible. Thus, the temperature dependence of χ 1/4 t is clearly different for these two temperature intervals with the exponent for T > 250 MeV being consistent with the prediction within DIGA. We smoothly interpolate between the two different fits at an intermediate T = 270 MeV, the resulting continuum estimate is shown in Fig. 5. l χ disc ) 1/4 (right). In the right panel, our results are compared with the continuum extrapolated results obtained in Ref. [36]. The solid orange line corresponds to a partial two-loop DIGA calculation with µ = πT , K = 1.9 and αs(µ = 1.5 GeV) = 0.336, while the band is obtained from the variation of αs by 1σ around this central value as well as variation of the scale µ by a factor two (see text).
A similar analysis was performed for (m 2 l χ disc ) 1/4 . In this case the low and high temperature regions were defined as 165 MeV ≤ T ≤ 240 MeV and 240 MeV < T ≤ 504 MeV, respectively. From the fits we get b = 1.96 (22) in the low temperature region, while in the high temperature region, the resulting b = 2.22 (27). The two fits were matched at T = 235 MeV to obtain a continuum estimate, shown as a red band in Fig. 5. It is clearly evident that the continuum estimates for χ 1/4 t from a gluonic observable and that for (m 2 l χ disc ) 1/4 coincide. This is highly non-trivial and makes us confident that the continuum extrapolation is reliable even though the cutoff effects are important. A joint fit was also performed according to Eq. 7, allowing the parameters a 2 , a 4 , b 2 , b 4 and b 6 to be different for (m 2 l χ disc ) 1/4 and χ 1/4 t since the cut-off effects in these quantities are clearly different. We considered three different ranges 165 MeV ≤ T ≤ 210 MeV, 195 MeV ≤ T ≤ 300 MeV and 230 MeV ≤ T ≤ 504 MeV. We checked different fit ansätze, setting some of the parameters a i and b i to zero and also including or excluding the N τ = 6 data. All such trials resulted in χ 2 per degree of freedom for the fitting procedure to be about one or less. Therefore we simply averaged over all such fit results in each temperature region and used the spread in the fit to estimate the final error on our continuum extrapolation. We also checked that the error estimated this way is about the same as the statistical error of the representative fits. We matched the fit results for the three different regions at T = 194 MeV and T = 277 MeV respectively to obtain our final continuum estimate for χ t , shown in right panel of Fig. 5.
Comparing our continuum fit for χ t with the continuum extrapolated results from Ref. [36], we find that the exponent b from our fit is about factor three larger than the corresponding b reported there. The possible reason for this discrepancy could be due to the fact that χ t has been calculated only with N τ = 6, 8 lattices in Ref. [36] at high temperatures. As discussed earlier, the cutoff effects depend on aT = 1/N τ therefore, lattices with N τ ≥ 8 are needed for reliable continuum extrapolation. To verify this assertion we performed fits to χ 1/4 t data calculated for specific choices of lattice spacings, a = 0.0600, 0.0737 and 0.0825 fm which are close to the values a = 0.0572, 0.0707 and 0.0824 fm measured in the earlier work. Using identical fit ansatz (Ã 0 + a 2Ã 2 )(T c /T ) b for χ 1/4 t , we get b = 0.674(42) for T > 1.3T c , comparable to b = 0.671(40) for T > 1.2T c reported in Ref. [36]. We thus conclude that continuum extrapolation procedure for χ t is not reliable unless data from finer lattice spacings are used.

Instanton Gas and Axion Phenomenology
At high temperatures, the maximum size of instantons is limited due to the Debye screening therefore, the fluctuations of the topological charge can be described within DIGA [11]. In this section, we calculate χ t in QCD within DIGA following the general outline in Refs. [32,34]. The χ t for QCD with N f quark flavors can be written in terms of the density of instantons D(ρ) of size ρ as [32,34], where we use the 2-loop renormalization group (RG) invariant form [50] of the instanton density in M S scheme Here G(x) is the cut-off function that includes the effects of Debye screening in a thermal medium [11] and β 0 and β 1 are the one and two-loop coefficients of the QCD beta function, respectively. The constant d MS has been calculated in Refs. [51,52,53]. At temperatures T > T c , we expect the strange quark to contribute to χ t so we henceforth consider N f = 3 in our calculations. At oneloop level, χ t calculated from Eq. (8) varies as ∼ T −8 in QCD with three quark flavors. With the two-loop corrections included, the temperature dependence becomes more complicated. To evaluate χ t (T ) according to Eqs. (8,9) one needs to specify α s (µ) and the strange quark mass m s (µ) in M S scheme as function of the renormalization scale µ. For the running of coupling constant, we set the starting value α s (µ = 1.5GeV) = 0.336(+12)(−0.008) obtained from lattice calculations of the static quark anti-quark potential [54]. This value is compatible with the earlier determination of α s [55] but has smaller errors. We fix the strange quark mass m s (µ = 2GeV) = 93.6 MeV, as determined in Ref. [56]. Since we are interested in comparing to our lattice results, the value of light quark mass was set to m l = m s /20. For the scale dependence of the α s (µ) and quark masses we use four-loop RG equations as implemented in the RunDeC Mathematica package [57].
The central value of the renormalization scale is chosen to be µ = πT , motivated from the fact that the instanton size is always re-scaled with πT in the screening functions G(ρπT ). To match the partial two-loop result for χ t that we calculated within DIGA for µ = πT to our lattice results, we have to scale it by a factor K = 1.90 (35). This factor is obtained by requiring that the calculations within DIGA and the continuum extrapolated lattice results agree at T = 450 MeV. The value of the K-factor in QCD is similar to that obtained in Ref. [34] for SU (3) gauge theory. The comparison between the scaled DIGA and lattice results can be seen from Fig. 5, where the uncertainty band in the DIGA prediction was estimated by varying the renormalization scale µ by a factor of two and the input value of α s (1.5 GeV) by one sigma. The strong temperature dependence of the upper bound is due to the fact that the lower value of scale µ becomes quite small when temperatures are decreased, resulting in a larger value of α s . We note that while the two-loop corrections to the instanton density reduces the uncertainties due to the choice of renormalization scale, the difference between the one and two-loop results is not large, e.g. for T = 450 MeV it is only about 1%. One may thus wonder about the origin of the large K-factor needed to match the DIGA estimates with lattice results since the uncertainties due to higher loop effects in the instanton density D(ρ) are small. We believe the reason for a large value of K-factor is inherent in the physics of the high temperature QCD plasma. It is well known that the naive perturbative series for thermodynamic quantities is not convergent [58,59,60]. The situation is particularly worse for the Debye mass, where corrections to the leading order result could be as large as the leading order result itself [61,62,63] even up to very high temperatures. Therefore, the calculation of the instanton density in Ref. [11] should be extended to higher orders within DIGA for a better comparison with lattice results. This is clearly a very difficult task. For the Debye mass, higher order corrections do not modify significantly its temperature dependence, i.e. they appear as a constant multiplicative factor to the leading order Debye mass. Therefore, including higher order effects within DIGA by a K-factor appears to be reasonable and appears to work well in the temperature range probed by the lattice calculations.
Let us finally mention that at temperatures of about 1 GeV, the effects of the charm quarks need to be taken into consideration when calculating D(ρ). It is not clear at present, how to consistently treat the charm quarks within DIGA since they are neither light nor heavy compared to this temperature. Treating charm quarks as light degrees of freedom, the value of χ t obtained within DIGA increases by ∼ 11% at T = 1 GeV compared to the N f = 3 result. This is clearly a modest effect compared to the large loop effects encoded within the K-factor.
We can now predict a bound on the axion decay constant f a , if indeed the cold dark matter abundance Ω DM h 2 in the present day universe is due to QCD axions i.e., Ω a h 2 ≤ Ω DM h 2 . Assuming that the PQ symmetry breaking and the slow roll of the axion field happened after inflation and choosing the misalignment angle θ = a /f a , at PQ scale to be averaged over different regions of the universe such that θ 2 = π 2 /3, we estimate lines of constant T osc /f a if axions constitute a fraction r = Ω a /Ω DM of the present dark matter abundance. The lines of constant T osc /f a for r = 0.1, 1 are denoted by the solid and the dashed lines in Fig. 6 respectively, considering the latest PLANCK data [64] for Ω DM h 2 = 0.1199 (27). For estimating these lines, we assumed that the temperatures are high enough for the charm quarks and τ leptons to be still relativistic, whereas the bottom quarks are not considered. Clearly for the values of T osc beyond the bottom quark mass threshold, one has to consider them as relativistic degrees of freedom as well, which however will shift the current estimates of the lines of constant T osc /f a by only a few percent. From the condition of slow roll m a (T osc ) = 3H(T osc ) and using the relation m 2 a = χ t /f 2 a , where χ t is an input from QCD calculations, we can again constrain the allowed regions in the T osc −f a plane. The dark-orange band in Fig. 6 shows the allowed values for f a and T osc within DIGA with a K-factor of 1.9 (35) which is favored by our data for χ t obtained from lattice. The background orange band shows the corresponding bound within DIGA when the K-factor is varied between 1-2.5 whereas the red band is estimated using the continuum result of χ t taken from Ref. [36] considering a 1σ spread of the exponent b. As evident from Fig.  6, the T osc is extremely sensitive to the characteristic temperature exponent b of the χ t data. The values of T osc obtained from our analysis for any fixed value of f a within the range 10 11 − 10 13 GeV, is atleast a factor five lower than the corresponding T osc extracted using data from Ref. [36]. The best estimate of the bound for f a from our analysis comes out to be f a ≤ 1.2 × 10 12 GeV with a ten percent error. The current uncertainties within the DIGA which is reflected in the magnitude of the K-factor would only mildly effect the estimates of the bound for f a , changing it by ∼ 20% when the K-factor varies between 1-2.5 as seen in Fig. 6. This range of f a is to be probed from the measurements of axion mass in the current and next generation ADMX experiments [65].

Conclusions
In summary, we have presented the continuum extrapolated results for topological susceptibility in 2+1 flavor QCD for T 3 T c from first principles lattice calculations using Symanzik flow technique. The exponent b that controls the fall-off of χ 1/4 t with temperature is found to be distinctly different in temperature intervals below and above 250 MeV. For T > 250 MeV, we show that b ∼ 2, in agreement with the calculations when considering that the topological properties of QCD can be well described as a dilute gas of instantons. The main conclusion from this work is that the χ 1/4 t in QCD determined non-perturbatively for T > 250 MeV can be very well described with the corresponding (partial) twoloop result within DIGA when the later is scaled up with a factor K = 1.90 (35). This scaling factor is similar to the recent estimate in pure gauge theory [34].
Our results are significantly lower in magnitude than the other continuum estimate of χ t in QCD with dynamical fermions [36] for T > 250 MeV. This could possibly due to the fact that the previous study uses relatively coarser lattices. We show that the cut-off effects in χ t depends on aT rather than aΛ QCD Figure 6: The allowed values of fa and Tosc, if indeed the present day dark matter abundance is due to the QCD axions. The input from QCD for determining the allowed regions is the magnitude of χt at the time of slow roll. The dark-orange band correspond to estimates within DIGA with K = 1.90(35) that is favored by our lattice data for χt whereas the red band is obtained using the data from Ref. [36]. The light-orange band corresponds to estimates of χt using DIGA with a K-factor between 1-2.5. The blue (dashed) and green lines denote the upper bounds for fa for different values of Tosc, if axions constitute the total or 10% respectively, of the present day dark matter abundance (for details see text). and hence extremely sensitive to the lattice spacing. Considering this fact, our continuum extrapolation has been carefully performed. We also calculated the disconnected part of chiral susceptibility which gives an independent estimate for χ t in the high temperature phase of QCD when the chiral symmetry is effectively restored. The continuum estimate of χ t from this observable matches with our result using the gauge definition of topological charge which gives us further confidence in our procedure.
It would be desirable to check how sensitive is the scale factor K on the matching procedure, particularly when the matching of the lattice result for χ t with the one within DIGA is performed at a temperature ∼ 1 GeV. With the present day algorithms, we do not observe any signal for χ t at finer lattice spacings for T > 450 MeV hence it is necessary to explore new algorithms that sample different topological sectors more efficiently at higher temperatures (see e.g. [66]). However for axion phenomenology, the effects due to large uncertainties in K-factor is fairly mild. A change in the magnitude from K = 2.5 to K = 1, would reduce the estimates of axion decay constant by 20% which is within the current experimental uncertainties for axion detection. With the temperature dependence of χ t for T > 1.5 T c now fairly well understood in terms of the microscopic topological constituents, it would be important to study in detail the topological properties in the vicinity of the chiral crossover transition in QCD.