The semi-classical approximation at high temperature revisited

We revisit the semi-classical calculation of the size distribution of instantons at finite temperature in non-abelian gauge theories in four dimensions. The relevant functional determinants were first calculated in the seminal work of Gross, Pisarski and Yaffe and the results were used for a wide variety of applications including axions most recently. In this work we show that the uncertainty on the numerical evaluations and semi-analytical expressions are two orders of magnitude larger than claimed. As a result various quantities computed from the size distribution need to be reevaluated, for instance the resulting relative error on the topological susceptibility at arbitrarily high temperatures is about 5% for QCD and about 10% for SU(3) Yang-Mills theory. With higher rank gauge groups this discrepancy is even higher. We also provide a simple semi-analytical formula for the size distribution with absolute error 2 · 10. In addition we also correct the over-all constant of the instanton size distribution in the MS scheme which was widely used incorrectly in the literature if non-trivial fermion content is present.


Introduction
This paper is concerned with the semi-classical study of gauge theories. At zero temperature or large space time volume the semi-classical picture is not reliable but at high temperature and finite spatial volume or small space time volume, the femto world [1][2][3][4][5][6], it is because the renormalized coupling runs with the relevant scale, µ ∼ T or µ ∼ 1/L, and becomes small. In these regimes the path integral can accurately be computed as the saddle points associated with instantons and the perturbative fluctuations around them. At asymptotically large temperatures or asymptotically small space time volumes the higher charge sectors are suppressed relative to lower ones. If an observable only receives contributions from the non-zero sectors then the leading contribution is from the 1-instanton sector. For instance the topological susceptibility is such a quantity and the leading semi-classical result is fully determined by the size distribution of 1-instantons. This size distribution at high temperature is the object of study in the present paper.
There has been renewed interest in semi-classical results for the topological susceptibility in QCD at high temperature because of applications in axion physics; see [7] and references therein. Basically, the topological susceptibility at high temperature can be used to constrain the amount of axions as dark matter components and its mass. In order to have results from first principles lattice calculations are ideal. Unfortunately the topological susceptibility decreases fast for T > T c hence a Monte-Carlo simulation will have a hard time achieving high precision because only very few configurations fall into the non-zero sectors. Even pure Yang-Mills theory is challenging in this respect [8][9][10][11]. Beyond pure Yang-Mills, results with dynamical fermions are also available [12][13][14][15] with full QCD [14] reaching temperatures up to T /T c ∼ 13 − 14.
It is useful to compare the non-perturbative lattice results both in pure Yang-Mills and in full QCD with the semi-classical results. The relevant semi-classical expressions at finite temperature were first obtained in [16,17]. The instanton size distribution consists of two parts, the zero temperature expression and an extra factor responsible for the temperature dependence. In this work we report the correct zero temperature expression in the MS-scheme which was frequently used incorrectly once light fermions were included and secondly we also report that the factor responsible for the temperature dependence which was numerically evaluated in [17] has an uncertainty that is two orders of magnitude larger than claimed. In pure Yang-Mills theory only the latter issue is relevant and leads to an increase in the semi-classical prediction. Once fermions are included as in QCD, both issues become relevant and while the latter still leads to an increase, the former leads to a similar decrease.
The organization of the paper is as follows. In section 2 the basics of the semi-classical expansion as it relates to high temperatures is summarized, in section 3 the instanton size distribution is converted to the MS scheme and a frequent error is identified in the literature if non-trivial fermion content is present. The main result of the paper is contained in 4 where the temperature dependence of the instanton size distribution is calculated to high precision and in section 5 the results are applied to the topological susceptibility. Finally we end with some comments and conclusions in section 6.
2 Semi-classical expansion at high temperature We will consider SU (N ) gauge theory with N f flavors of light fermions in the fundamental representation. The semi-classical approximation provides a consistent and reliable expansion of the path integral at fixed spatial volume and asymptotically high temperatures similarly to the situation in small 4-volume or femto world [1][2][3][4][5][6].
It is useful to consider the theory at non-zero ϑ angle and the partition function is then given by where Z Q is the result of integrating over gauge fields with given topological charge Q. At high temperatures the Q > 1 sectors are suppressed relative to Q = 1 both exponentially in the inverse coupling as well as algebraically in the inverse temperature, similarly to the femto world. Hence we have the topological susceptibility as where V = L 3 /T is the space-time volume and L 3 is the finite spatial box. Note that here we envisage a finite and large spatial volume and the limit T → ∞. In this case all further terms in (2.2) are indeed negligible. For instance there is no need to invoke the dilute instanton gas approximation in order to obtain (2.2) at T → ∞. The subject of the present paper is Z 1 /Z 0 and our only objective is to study Z(ϑ) for T → ∞ so there is no reliance on the dilute instanton gas picture at all. In order to have a self contained presentation and clarify its relationship with the general semi-classical picture, we will nonetheless summarize the dilute instanton gas model in appendix A. Now Z 1 /Z 0 can be reliably calculated at asymptotically high temperatures as the contribution of the 1-instanton together with the perturbative fluctuations around it. The leading term is obtained by integration over the moduli space of 1-instanton solutions: (z µ , ̺) where z µ determines the position and ̺ the scale of the instanton. An integration over the arbitrary gauge orientation of the instanton is implied leading to an N -dependent over-all factor. Integration over the position gives a space-time volume factor V .
-2 -What we will be concerned with is the integration over the size ̺ or more precisely the size distribution n(̺, T ), The size distribution at finite temperature [16,17] can be expressed by the analogous quantity at T = 0 together with an explicitly temperature dependent factor S(̺, T ), The zero temperature size distribution n(̺) with N f light fermions at leading order is [18,19] with β 1 = 11/3N − 2/3N f the first coefficient of the β-function and C is a scheme-dependent constant. Its value will be detailed in the next section. The 2-loop expression for n(̺) is actually available [20] but for our purposes the leading order result is sufficient. Note that the constant C in principle depends on ̺m i given by the contribution of the non-zero mode fermion determinant [21][22][23][24] but the massless limit can be taken for the asymptotically high temperature limit. It is assumed that the theory is renormalized at T = 0 by introducing a renormalized coupling g 2 (µ) and renormalized masses m i (µ) at some renormalization scale µ in a chosen scheme. Once this is done n(̺) is finite. The presence of a finite temperature does not introduce any new divergences hence S(̺, T ) is also finite and by construction is zero at T = 0. For definiteness the renormalization scale can be chosen to be µ = cT with some O(1) constant c and the dependence on c is indicative of higher order corrections. Since S(̺, T ) is dimensionless the only dependence is in fact through the variable λ = π̺T . We have the result [17,25], where the function A(λ), encodes the temperature dependence of the 1-loop determinant in the background of the instanton. Here Π(τ, r) determines the 1-instanton solution at finite temperature, the periodic Harrington-Shepard solution [26], while Π 0 (t, r) is the corresponding function at zero temperature [27], .
Both the first and second term in (2.7) are divergent separately and the second term, formally a constant, corresponding to zero temperature is subtracted in order to make A(λ) finite and A(0) = 0. Summarizing the above, once A(λ) is found the high temperature limit of Z 1 /Z 0 and consequently that of the topological susceptibility is known.
The main object of study in the present paper is A(λ). The definition (2.7) can not be evaluated analytically. Numerically, the task is performing a 2 dimensional integral over (τ, r) and the careful subtraction of an infinite constant. This was first attempted in [17] and the numerical result was parametrized as, with α = 0.01289764 and δ = 0.15858 and a quoted absolute numerical uncertainty of at most 12 · 5 · 10 −5 = 6 · 10 −4 . Note that δ in [17] was denoted by γ but in order not to confuse it with Euler's constant later we will label it by δ. Formula (2.9) was then used in all subsequent applications [8][9][10]28] where the topological susceptibility was studied at high temperatures and compared with the semi-classical results.
In the present paper we show that the actual numerical accuracy of (2.9) is two orders of magnitude larger than claimed in [17] at around λ ∼ O(1) leading to about a 5% mismatch for the topological susceptibility in QCD or 10% mismatch for SU (3) pure Yang-Mills. For higher rank gauge groups the mismatch is even higher. Any absolute error on A(λ) needs to be scaled by N for the absolute accuracy of S(T, ̺) and this scaled absolute accuracy becomes the relative accuracy of the topological susceptibility.
Before calculating A(λ) accurately in section 4 we make a comment in the next section on the constant C appearing in (2.5). It turns out its flavor number dependence in the MS scheme was widely used incorrectly in the literature, especially recently in axion mass estimates as well as comparisons with lattice results.

Scheme dependence -conversion to MS
Historically, there has been quite a bit of confusion about the constant C in (2.5). The topological susceptibility is of course a finite RG invariant quantity but its perturbative expansion in terms of a scheme dependent running coupling involve the scheme dependent constant C. In the original paper [18] a Pauli-Villars regulator was used for SU (2) which was extended to SU (N ) still with a Pauli-Villars regulator in [19]. The Pauli-Villars results in [18] were correct except for some trivial mistakes corrected later in an erratum, but the conversion to the MS scheme was incorrect. The correct conversion factor was later found in [29] and confirmed in [30]. It is possible to work directly within the MS scheme lending further support to the correct conversion factor [31]. The conversion between various schemes is via the Λ-parameter ratios, more precisely the constants C in two different schemes are related by For definiteness we record the ones relevant for us below [29,32], Note that wrong Λ-parameter ratios were published for instance in [33] as well as [34]. The MS scheme is one of the most frequently used schemes and the constant C can trivially be obtained from either -4 -the Pauli-Villars or the MS scheme. The result is with the derivative of the Riemann ζ-function ζ ′ (s). Even though there is now consensus on the pure Yang-Mills case, the flavor coefficient c 2 in MS was frequently quoted [35] as 2α(1/2) = 0.291746 in the notation of [18] and subsequently this value was used in applications [8][9][10]. 1 However 2α(1/2) equals −4ζ ′ (−1) − 5 36 − 1 3 log 2 which is the relevant result within the Pauli-Villars scheme, the one obtained in [18]. It appears that only the coefficient c 1 was converted to MS and the flavor coefficient c 2 was not although of course β 1 depends on N f . The mismatch in c 2 is precisely 1 33 = 2 3 · 1 22 which should be included using (3.1) and (3.2). Here 2 3 comes from the β-function and 1 22 from the ratio of the Λ-parameters.
Hence we believe that the fully correct MS expressions for the instanton size distribution including fermions are given by (2.5)

Calculation of A(λ)
It is seemingly a relatively straightforward task to evaluate (2.7). However analytical evaluation is apparently impossible and numerical integration is tricky primarily because of the subtraction of the infinite constant and a need for high arithmetic precision.
The 2 dimensional integral responsible for the T -dependence is over (τ, r) where τ is a periodic variable with period 1/T . First, we will perform this integral over S 1 analytically and only the second integral over r will be left numerically. Let us rescale both τ and r by 1/(2πT ) so that they become dimensionless and τ periodic with period 2π. Then we have where prime denotes differentiation with respect to r. Now the τ -integral can be done analytically for instance via the residue theorem and we obtain, where the X and Y derivatives should be taken at fix a, b, c and in the end a, b, c, X, Y are all simple functions of r. 1 In an earlier work 0.153 was reported [36] based on [37] which is also incorrect.
Now we have both terms in (2.7) in terms of a dimensionless r and arrive at, The integrals with r 2 I(r) or r 2 I 0 (r) are separately divergent and the divergence comes from the origin. For finite and fixed λ they both behave as 2/r + O(r) for r ≪ 1 hence the difference is finite. For r ≫ 1 both terms are integrable separately. More specifically we have, (4.4)

Hence (4.3) defines a convergent integral for A(λ).
Numerical evaluation is straightforward although high arithmetic precision is required because of large cancellations. First, the O(1/r) cancellation between r 2 I(r) and r 2 I 0 (r) in (4.3) needs to be resolved and also the evaluation of I(r) itself via (4.2) involve large cancellations especially for small λ. As we will see the discrepancy between our result and [17] is precisely here in the small λ region. We have found that keeping O(100) significant digits is nonetheless sufficient. In order to have as much analytic control as possible it is useful to derive the small λ and large λ asymptotics which will be used to check the numerical results.

Small and large λ asymptotics
It is relatively straightforward to obtain the λ → 0 asymptotics. The simplest is to change the integration variable r → λr in (4.3), expand the integrand in λ → 0 and perform the r-integral. The first 3 terms are, where the O(λ 7 ) terms may include fractional powers and/or logarithms and it turns out that the log on the right hand side has the same first 3 Taylor coefficients as the ones obtained from the integral. These 3 terms will be referred to as LO, NLO, NNLO in the following. The λ → ∞ regime is more difficult to obtain. Starting with the original expression (4.3) one may split the integral to two intervals, (0, r 0 ) and (r 0 , ∞) with some fixed r 0 . In the λ → ∞ limit the first term is finite whereas the second term is logarithmically divergent.
More specifically, the first term can be straightforwardly expanded in λ → ∞ leading to a constant and subleading O(1/λ) expressions. In the integral (r 0 , ∞) all exponentials e −r may be dropped, the integration over r can then be performed and the result can be expanded in λ → ∞. Finally in the sum of the two integrals, order by order in λ, r 0 → ∞ can be taken in order to justify the dropping -6 - of all exponentials. This procedure leads to, − γ + log π = 1.25338375 The divergent − log(λ 2 ) piece is originating from the second integral (r 0 , ∞) in the above split. The 4 terms above will be referred to as LO, NLO, NNLO and N 3 LO in the following, similarly to the small λ expansion. Clearly, the parametrization (2.9) satisfies the small λ expansion (4.5). However, even though the small λ expansion of the second term of (2.9) is O(λ 12 ) its coefficient is rather large 12α/δ 8 ≈ 4 · 10 5 and distorts the small λ behavior. Quite remarkably the NLO large λ expansion (4.6) matches (2.9) since the constant term in (2.9) is log(3) + 12α which agrees with C 1 to 6 significant digits. However the NNLO and N 3 LO terms of the large λ expansions do not match (2.9) but for the application to the topological susceptibility this matters much less than the mismatch at λ = O(1).

Numerical results
After making sure the arithmetic precision is high enough it is straightforward to calculate A(λ) numerically using (4.3). First, numerical integration is done on the interval (0, 8) with either the trapezoid or Simpson's rule with step size 10 −4 . The difference between the two schemes is at most 10 −6 which is the estimated accuracy. The remaining integral on (8, ∞) can be approximated by dropping all exponentials e −r in the integrand and performing the integral analytically. This second piece is then added to the numerical integral on the (0, 8) finite interval and amounts to at most O(10 −4 ). We have checked that the exact integrand at r = 8 agrees with the one where the exponentials are dropped to either O(10 −6 ) or the added term is at most O(10 −9 ). The procedure is then reliable in absolute precision to at least O(10 −6 ).
The final result is shown in figure 1. The numerical result is compared with the increasing orders of both the small and large λ expansions in figure 2 showing increasing agreement order by order in both regimes as expected.  Let us now compare with [17]. We show the difference between the parametrization (2.9) and the numerical result in the left panel of figure 3. Clearly, the difference is well outside the claimed absolute precision of 12 · 5 · 10 −5 = 6 · 10 −4 . Note that any absolute error on A(λ) will turn into a relative error on the size distribution. Unfortunately, the differences are largest for λ = O(1) which is precisely the region which contributes the most to the topological susceptibility through the integration over the instanton sizes.
Once a reliable result for A(λ) is obtained the topological susceptibility at asymptotically high temperatures can be calculated using (2.3) -(2.6). Here we are interested in studying the discrepancy between using (2.9) and our result (4.7) for the topological susceptibility for various gauge groups and fermion content. It is a simple exercise to compute the integrals over the instanton size and it turns out that there is an approximately 10% discrepancy between them for SU (3) pure Yang-Mills and about 7%, 6% and 4% for 2, 3 and 4 light flavors, respectively, with the correct results being larger by these amounts. For fixed N the largest relative difference is always at N f = 0 and for increasing N the discrepancy grows. For instance with SU (10) we have about 22% while with SU (20) it is about 40% which is quite sizable. The reason for this exploding discrepancy is that in the large-N limit the size distribution is more and more peaked around λ 0 ≈ 1.6935 which unfortunately is approximately where the discrepancy between our result and [17] is largest in absolute terms, see left panel of figure  3, and any absolute discrepancy is scaled up by N . The large absolute discrepancy is translated into a large relative discrepancy for the topological susceptibility. Since only instantons of the critical size λ 0 contribute in the large-N limit, a straightforward explicit expression can in fact be obtained for the topological susceptibility.
Recently there has been renewed interest in the non-perturbative determination of the topological susceptibility in lattice calculations [8-15, 28, 38] motivated mainly by axion physics. Even the pure Yang-Mills case is challenging and obtaining high precision results at high enough temperatures is difficult. It is nevertheless clear [9,11] that at T /T c = 2.5 the temperature is not high enough for the semi-classical topological susceptibility to agree with the non-perturbative lattice result, however at T /T c = 4.1 the two already only deviate [11] within 3σ. The lattice result [11] extrapolated to the continuum is log χ/T 4 c = −12.47 (21), whereas the semi-classical one is −13.80(10)(40) where the 2-loop result for the susceptibility itself [20] and 5-loop running [39,40] for the coupling was used (however 3-loop or 4-loop running for the coupling gives identical results). The first error estimate includes the residual dependence on the scale µ and the second (dominant) one originates from the uncertainty of T c /Λ MS = 1.26 (7) which was also obtained in lattice calculations [41]. Note that without lattice input it is only possible to obtain the semi-classical susceptibility in Λ MS units with the temperature also measured in Λ MS , however the lattice results for the susceptibility are in T c units so the uncertainty related to T c /Λ MS is unavoidable. The resulting uncertainty is rather large because the β 1 − 4 = 7 power of T c /Λ MS enters the semi-classical susceptibility.
In QCD the situation [14] is similar in that the semi-classical susceptibility is in Λ MS units and with 4 flavors Λ MS = 292 (16)MeV; the uncertainty is about 5% [42]. This uncertainty is amplified by the even higher factor β 1 − 4 + N f = 8.33 in the logarithm of the susceptibility, leading to the semiclassical result log χ/MeV 4 = 1.15(3)(46) at T = 2000 MeV where again the inherent uncertainty of the semi-classical result is negligible compared to the one associated with the uncertainty of the scale. The corresponding 3 + 1 flavor lattice result extrapolated to the continuum is 3.99(68) hence the deviation between the two is about 3.5σ; for higher temperature the deviation is decreasing [14]. The semi-classical result is consistently below the lattice result, similarly to pure Yang-Mills, and it is most likely that further perturbative O(g 2 ) corrections to the 2-loop semi-classical topological susceptibility are responsible for the current 3.5σ deviation at the relatively high temperature reachable today by lattice calculations. Note that in QCD the two issues reported in this paper (more accurate A(λ) and correct C MS ) nearly cancel each other in the final result.
-9 -In this paper we have studied the instanton size distribution in gauge theories with light fermions semiclassically. At large but finite spatial volume and asymptotically high temperatures the contribution of the 1-instanton sector is the exact result, there is no need to invoke the dilute instanton gas model. We were motivated by the recent interest in the semi-classical calculation of the topological susceptibility at high temperature which in turn was motivated by axion physics.
It turned out that there were two issues in previous semi-classical calculations which led to incorrect results. One, even at zero temperature the over-all constant of the instanton size distribution was not correctly converted to MS when N f > 0 and second, the numerically evaluated 1-loop determinant in the instanton background was two orders of magnitude less precise than originally thought. We have managed to evaluate these determinants to high precision with an absolute uncertainty of at most 2 ·10 −4 by doing the periodic integral over the temperature analytically leaving only a one dimensional integral numerically which could be handled by standard methods. Although the discrepancy with previous work is not large the integration over the instanton size is peaked at around the same point where the discrepancy is largest. Any discrepancy in absolute terms is then converted to a relative discrepancy for the topological susceptibility which can be 10% for SU (3) or larger for higher gauge groups.
Once correct semi-classical expressions are obtained a meaningful comparison with the lattice results can be carried out however the uncertainty of the scale Λ MS needs to be taken into account. Taking everything into account in SU (3) pure Yang-Mills the semi-classical and lattice results already at T /T c = 4.1 are within 3σ [11] and with 3 + 1 flavor QCD at T = 2000 MeV the deviation is below 3.5σ [14]. The remaining deviation can surely be reduced by simply considering higher temperatures which is however quite demanding on the lattice or by calculating further perturbative corrections to the semi-classical result. where I ν is the modified Bessel-function of the first kind. The probability to find a gauge field with charge Q is then, Hence the single parameter χ determines all ϑ-dependence and the full topological charge distribution. Unsurprisingly we have, i.e. the single parameter χ is indeed the topological susceptibility. We emphasize again that a priori there does not need to be any reference to weak coupling or semi-classical objects at all in order for (A.1) -(A.3) to hold with some parameter χ. Perhaps it would be more accurate to call these set of assumptions the ideal gas model rather than the dilute instanton gas model in order to make this distinction clear. This is because the following two questions are independent at a given temperature: (A) whether the semi-classical expansion at fixed order is accurate or not and (B) whether the topological charge distribution and the full ϑ-dependence follows (A.1) -(A.3) with a single parameter χ or not. As to question (A) we know that the semi-classical expansion is reliable at finite spatial volume and asymptotically high temperatures without any further assumptions. At asymptotically high temperatures (B) is also true with χ given by the semi-classical result (2.2) however at asymptotically high temperatures all the additional assumptions (A.1) -(A.3) are not necessary, the semi-classical picture alone is sufficient to calculate anything. Question (B) becomes a non-trivial proposition at finite temperature T > T c where perhaps it does hold, although there is no a priori reason for it, with some parameter χ which is however different from the semiclassical result (2.2). This is precisely what appears to be the case based on detailed lattice calculations in the pure Yang-Mills case [43,44]. Immediately above T c the topological charge distribution seems to follow (A.3) but with a χ which is non-perturbative and does not agree with the semi-classical formula.
This comparison is what we would call testing the ideal gas model. It would be very interesting to extend these results to full QCD.