Axion cosmology, lattice QCD and the dilute instanton gas

Axions are one of the most attractive dark matter candidates. The evolution of their number density in the early universe can be determined by calculating the topological susceptibility $\chi(T)$ of QCD as a function of the temperature. Lattice QCD provides an ab initio technique to carry out such a calculation. A full result needs two ingredients: physical quark masses and a controlled continuum extrapolation from non-vanishing to zero lattice spacings. We determine $\chi(T)$ in the quenched framework (infinitely large quark masses) and extrapolate its values to the continuum limit. The results are compared with the prediction of the dilute instanton gas approximation (DIGA). A nice agreement is found for the temperature dependence, whereas the overall normalization of the DIGA result still differs from the non-perturbative continuum extrapolated lattice results by a factor of order ten. We discuss the consequences of our findings for the prediction of the amount of axion dark matter.


Introduction
One of the greatest puzzles in particle physics is the nature of dark matter. A prominent particle candidate for the latter is the axion A [1, 2]: a pseudo Nambu-Goldstone boson arising from the breaking of a hypothetical global chiral U(1) extension [3] of the Standard Model at an energy scale f A much larger than the electroweak scale.
A key input for the prediction of the amount of axion dark matter [4][5][6] is its potential as a function of the temperature, V(A, T ). It is related to the free energy density in QCD, F(θ, T ) ≡ − ln Z(θ, T )/V, via where V is the Euclidean space-time volume. The angle θ enters the Euclidean QCD Lagrangian via the additional term involving the topological charge density with F a µν being the gluonic field strength and α s ≡ g 2 s /(4π) the fine structure constant of strong interactions. On general grounds, the free energy density and thus the axion potential has an absolute minimum at θ = A/ f A = 0. In fact, this is the reason why in this extension of the Standard Model there is no strong CP problem [3]. The curvature around this minimum determines the axion mass m A at finite temperature, in terms of the topological susceptibility, i.e. the variance of the θ = 0 topological charge distribution, where Q ≡ d 4 xq(x). Similarly, self-interaction terms in the potential, e.g. the A 4 term occurring in the expansion of the free energy density around θ = A/ f A = 0, are determined by higher moments, These non-perturbative quantities enter in a prediction of the lower bound on the fractional contribution of axions to the observed cold dark matter as follows 1 R A 10 5 θ 2 f ( θ 2 ) GeV 3 g * S (T osc )T 3 osc χ(T osc ) where θ 2 is the variance of the spatial distribution of the initial values of the axion field A/ f A before the formation of the dark matter condensate by the misalignment mechanism, which occurs, when the Hubble expansion rate gets of the order of the axion mass, m A (T osc ) = 3H(T osc ), i.e. at a temperature T osc 50 GeV The function f ( θ 2 ) is taking into account anharmonicity effects arising from the self-interaction terms in V(A, T ) and depends on the specific form of the potential. The functions g * R (T ) and g * S (T ) denote the effective number of relativistic energy and entropy degrees of freedom, respectively. What is urgently needed for axion cosmology is thus a precise determination of the topological susceptibility and higher moments of the topological charge distribution. In this context, most predictions have been entirely based on the semi-classical expansion of the Euclidean path integral of finite temperature QCD around a dilute gas of instantons -finite action minima of the Euclidean action with unit topological charge -see e.g. Ref. [8] for an early exhausting study and Ref. [9] for a recent update concerning the quark masses. A comparative study of these predictions based on the dilute instanton gas approximation (DIGA) has been carried out in Ref. [10], where also an analysis in terms of a phenomenological instanton liquid model (IILM) [11] is presented. However, up to now, in all DIGA investigations of the topological susceptibility only the one-loop expression in the expansion around the instanton background field was used. This results in a strong renormalization scale dependence and thus large uncertainties which were neglected in the previous DIGA based predictions. In fact, in the temperature range T osc ∼ GeV of interest, one expects a large uncertainty in the overall normalization due to the neglection of higher order loop effects, since in this region α s (T osc ) is not small. We will exploit in this letter both the one-loop DIGA result as well as its two-loop renormalization group improved (RGI) version in order to study the theoretical uncertainties arising from higher loop corrections. Most importantly, we compare these predictions with the outcome of our lattice based fully non-perturbative results.
Actually, there have been a number of lattice calculations of χ(T ) and b 2 (T ) at temperatures below or slightly above the QCD phase transition, mostly in quenched QCD, see e.g. [12][13][14]. Here we go beyond those lattice calculations. We extend the available temperature range and carry out a controlled continuum extrapolation for this extended range. In addition, note that there has been no quantitative investigation whether and where the lattice results turn into the DIGA results. We will present our new high-quality lattice data for the topological susceptibility in quenched QCD (i.e. neglecting the effects of light quarks) and compare them quantitatively to the DIGA result.

Axion potential coefficients from the lattice
On the lattice, the topological susceptibility is measured on the torus as the second moment of the distribution of the global topological charge where Q is any renormalized discretization of the global topological charge, and V is the four-volume of the lattice. There are a lot of different fermionic and gluonic definitions of Q available. We choose a gluonic definition based on the Wilson flow [15,16], which has the correct continuum limit similarly to the fermionic definitions but is numerically a lot cheaper. In particular we evolve our gauge field configurations with the Wilson plaquette action to a flow time t and define the global topological charge as the integral over the clover definition of the topological charge density. This definition gives a properly renormalized observable when the flow time t is fixed in physical units.
We use a tree-level Symanzik improved gauge action. Our temperatures range from below T c up to 4T c . Here T c is the critical temperature, which is the quantity used for scale setting. The critical temperatures for different lattice spacings were determined in earlier work [17,18]. For the whole temperature range we keep the spatial lattice size approximately at L = 2/T c . We checked with dedicated high volume runs at 1.5T c and 2T c that this volume is sufficiently large to have negligible finite volume corrections on Q 2 . This is shown in Fig.1. Our spatial geometry is L × L × 2L to enable tests of subvolume methods which will be reported separately -here we only use the full volume. For all temperatures we have three lattice spacings (aT ) −1 = 5, 6, 8 to be able to perform an independent continuum extrapolation for every temperature. The local heatbath/overrelaxation algorithm is used for the update, one sweep consists of 1 heatbath and 4 overrelaxation steps. We found that the autocorrelation time of the topological charge depends weakly on (aT ), i.e. if the temperature is increased by decreasing the lattice spacing. The number of update sweeps between measurements was chosen in accordance with the autocorrelation time. Tab. 1 lists the simulation points with the number of sweeps.
We integrated the Wilson flow numerically to a maximum flow-time of about 8t ≈ 1/(2T 2 c ) for all temperatures. Figure 2 gives the dependence of the susceptibility on the flow time for T = 2T c . While in the continuum limit the result is independent of the choice of the flow time t, different choices have very different lattice artefacts. For small flow times the different lattice spac-T/T c N t N s N sweeps 0.9 5 Table 1: List of simulation points: temperatures, lattice sizes N t = 1/(aT ), N s = L/a, and number of sweeps are given.  ings give very different results. For larger flow times the expected plateau behavior can be observed for each lattice spacing and the lattice artefacts also decrease significantly. The choice of the flow time brings in some arbitrariness into the analysis, however the continuum result should not depend on this choice once t is fixed in physical units. But this is certainly a subleading source of error compared to the statistical error due to the rare topology tunneling events at high temperatures. In this analysis we choose a temperature dependent flow time for the evaluation of Q 2 as For low/high temperatures this means a temperature independent/dependent flow time. This choice is safely in the expected plateau region for all temperatures.
The resulting values for the susceptibility are plotted in Fig. 3. This plot also gives the result of a global continuum extrapolation using a set of temperatures and the 6-parameter power law ansatz where χ 0 , T 0 , and b are fit parameters giving the continuum limit. χ 0 , T 0 , and b are fit parameters describing the deviation from the continuum limit. The fit parameter T 0 is included as a consistency check and should give 1 in units of T c . This is satisfied by the fit result.
where the first error is the statistical, the second is the systematic. The point-wise continuum extrapolation is consistent with the global fit, evidently the latter has smaller errors for large temperatures. Note that though a controlled continuum extrapolation is possible using three lattice spacings, estimating the systematic uncertainty of this extrapolation would require at least one more lattice resolution.
In a recent analysis [13] the topological susceptibility was calculated using the techniques [19,20]. The calculation was carried out at two temporal extensions, corresponding to two lattice spacings at each temperature. The exponent b = −5.64(4) was found which differs from our value. Note however, that our temperature range is larger, thus we are closer to the applicability range of the DIGA. Furthermore, the two lattice spacings were not sufficient for a controlled continuum extrapolation thus uncertainties related to this final step are not included in the result of [13].
We have also determined the second important coefficient b 2 of the axion potential, characterizing its anharmonicity, by measuring the observable The result is plotted in Fig. 4.

Comparison between lattice and DIGA results
In this section, we confront the lattice results with the ones obtained from the DIGA framework. For the sake of completeness let us collect first the available formulas for the latter [21][22][23][24][25][26][27][28]. At very high temperatures, far above the QCD phase transition, it makes sense to infer the θ dependence of QCD from the DIGA, in which the partition function can be written as [22], where Z I = ZĪ is the contribution arising from the expansion of the path integral around a single instanton I (anti-instantonĪ). It follows directly that the potential has the form from which one infers This can be confronted right-away with our lattice results, cf. Fig. 4. Similar to Ref. [12] we find that the prediction from the DIGA for b 2 is reached already at surprisingly low values of T/T c 1.
The whole temperature dependence of the axion potential arises in the DIGA through the topological susceptibility, which in this case is explicitly given by (15) in terms of the instanton size distribution at zero temperature, D(ρ), and a factor G(πρT ) taking into account  (7) as determined from the lattice in Ref. [18].
finite temperature effects. The former is known in the framework of the semiclassical expansion around the instanton for small α s (µ r ) ln(ρ µ r ) and ρ m i (µ r ), where α s is the strong coupling, µ r is the renormalization scale and m i (µ r ) are the running quark masses. To one-loop accuracy, it is given by 2 ×(ρ µ r ) β 0 1 + O(α MS (µ r )) , with d MS = e 5/6 π 2 e −4.534122 ; β 0 = 11.
At finite temperature, electric Debye screening prohibits the existence of large-scale coherent fields in the plasma, leading to the factor [24,25], with and α = 0.01289764 and γ = 0.15858, in Eq. (15). This factor cuts off the integration over the size distribution in Eq. (15) at x = πρT ∼ 1 and ensures the validity of the DIGA at large temperatures, at which α s (πT ) 1.
T/T c 1.5 2 3 4 5 b (κ = 0.6) -6.04 -6.26 -6.43 -6.50 -6.55 b (κ = 1) -6.37 -6.46 -6.55 -6.59 -6.62 b (κ = 2) -6.55 -6.59 -6.64 -6.67 -6.69  Collecting all the factors, the topological susceptibility, in the one-loop DIGA, reads This result, however, still suffers from a sizeable dependence on the renormalization scale µ r , reflecting the importance of the neglected two-loop and higher order contributions. In fact, it is tamed by taking into account the ultraviolet part of the two-loop correction. The latter has been calculated in Ref. [27] and shown to have exactly the form that the gauge coupling becomes a parameter running according to the renormalization group (RG). Therefore, the ultimate, all order result for the topological susceptibility becomes independent of µ r , for µ r → ∞. At two loop, the corrections amount to a factor (ρµ r ) (β 1 −12 β 0 ) α MS (µ r )/(4π) ; in D(ρ). Therefore, this RG improvement can be taken into account by replacing the factor I in Eq. (20) bỹ In fact, exploiting the two-loop RG improvement, the µ r dependence is heavily reduced, as is obvious from Fig. 5, where we have used as a natural renormalization scale µ r = κ/ρ m = κπT/1.2, with ρ m being approximately the maximum of the integrand ofĨ, and varied the remaining free parameter κ between 0.6 and 2. The renormalization scale dependence appears to be highly reduced in the regimes > 3T c and < T c . However, in an intermediate region, ∼ T c -2T c , it is comparable in size to the one at one-loop.
We present in Table 2 the power-law behavior predicted by the two-loop RGI DIGA at various temperatures, which can be compared to the fit (11) to the continuum lattice result. As far as the overall normalization of the DIGA result for χ is concerned, one still expects a large uncertainty in the temperature range available from the lattice. In fact, at these temperatures, α s is not small, see Tab. 3. Apart from the ultraviolet part, there will be a finite part of the two-loop correction which will affect mainly the overall normalization of χ and will depend on the temperature only logarithmically. Unfortunately, this finite part is not known, yet. Therefore, when comparing to the continuum extrapolated lattice results, we allow a multiplicative factor K to account for this uncertainty, i.e. we absorb the remaining higher loop uncertainties by replacing in Eq. (20) and the corresponding two-loop RGI expression. Clearly, the K-factor should approach unity at very large T/T c . Figure 6 nicely illustrates the agreement between the DIGA and the lattice result, if a K-factor of order ten is included 3 . More precisely, fitting the lattice continuum data with the rescaled DIGA expression in the temperature range T/T c ≥ 1, one finds while a fit in the temperature range T/T c ≥ 2 yields Apparently, in the temperature range accessible to the lattice, the higher order corrections to the pre-factor of the DIGA are still appreciable, but there are indications of a trend that the K-factor gets smaller, as expected, towards larger values of T/T c . The K-factor strongly depends on the value of T c //Λ

Conclusions
This paper presents lattice and DIGA calculations of the θ and temperature dependence of the free energy density of QCD in the quenched limit, i.e. for infinite quark masses. In the lattice approach the temperature ranges from 0.9T c to 4T c and thus significantly extends former results. The precise high temperature data allows to compare the lattice results with first principle calculations within the DIGA. For the evaluation of χ(T ) we use the two-loop RGI version of the instanton size distribution, which is less sensitive with respect to the renormalization scale. We find that the two different methods show nice agreement after one includes an overall correction factor in the DIGA results. This order ten correction is supposed to take into account higher loop contributions in the normalization, similar to socalled K-factors in perturbative QCD processes. These are significant, since in the considered temperature region the strong coupling constant is still large. Nevertheless, the achieved results motivate an investigation of more realistic models involving light quarks.
The possible consequences of our findings on the predictions of the amount of axion dark matter can be read off from Fig. 7. First, we have plotted the regions where axion DM would be overproduced with respect to observations and are thus excluded. The darker yellow region is excluded just by the axions produced from the misalignment mechanism, cf. Eq. (7), assuming a flat distribution of initial misalignment angles in the ob- Figure 7: Consequence of our findings for axion dark matter from quenched (top) and full QCD (bottom). The dark yellow region is excluded because R A > 1 even without a contribution from strings. The light yellow region indicates R A > 1 when string contributions are included. In the dark green region axions even including strings give a small contribution (R A < 0.1) to dark matter. The light green region indicates R A < 0.1 just from the misalignment mechanism. In order to compare the quenched results with the axion dark matter scenario we had to transform dimensionless quantities into dimensionful ones. This is not unambiguous since we do not have pure gauge theory in nature. Using different quantities to set the scale one gets 280 -310 MeV for T c . For illustrative purposes we use T c = 294 MeV which lies in the middle of the range and is a factor two larger than the full QCD transition temperature. Our lattice results are shown by the blue points and the two-loop RGI DIGA prediction by the gray band in the quenched Figure. In the unquenched case the blue and gray bands correspond to the two-loop RGI DIGA predictions with K = 1 and K = 9.22 ± 0.6, respectively. The dashed red line shows the IILM prediction of Refs. [10,11].
servable universe 4 θ ∈ [−π, π]. For the DIGA potential, χ(T )(1 − cos(θ)), we calculate θ 2 f ( θ 2 ) 6.6 from numerically evolving the axion field with different initial conditions. The exclusion extends to the light yellow region if we consider the DM axions produced from the decay of unstable axionic strings and domain walls according to Ref. [32]. The excluded region has to be compared with our lattice results (blue dots) and the two-loop DIGA calculation with the fitted K-factor, for which we have set the scale T c = 294 MeV (twice the dynamical value) for illustration purposes, see Fig.  7 (top) . Despite the arbitrariness we can derive a number of interesting conclusions. First note that our results correspond to a region where axion DM does not offer a viable cosmology because R A > 1. However, our lattice results are too short in T/T c only by a factor of 2-3. Further developments from the lattice side could in principle lead to a direct estimation of the axion mass. Meanwhile, the DIGA result allows a controlled extrapolation of our results over the interesting region, that can be used to set a constraint (here merely illustrative): 25 µeV m A 525 µeV assuming that at least ten percent of dark matter can be attributed to axions.
For the final calculation of the QCD axion abundance as dark matter candidate one needs to include the light quarks. This complicates the lattice evaluation by far. However, the DIGA approximation can still be carried out and self-consistently takes into account the running of the quark masses as well as their influence on the strong coupling constant. We present this preliminary result as a blue band in Fig. 7 (bottom). From what we learned in the quenched case, this DIGA result still misses an O(10) constant factor but the lack of lattice data in this case prevents from fitting the expected K−factor explicitly. We can tentatively estimate it to be the required factor that makes χ(T )/χ(0) = 1 at T c . Since the transition in full QCD is a crossover [33], T c is not unambiguous. Using the transition temperature defined by the chiral susceptibility, T c = 147(2)(3) MeV [34][35][36], we obtain K = 9.22±0.6. Very interestingly, this agrees with the IILM model calculation of Ref. [10,11] (red-dashed line) at the temperatures of interest. This factor of 9.22 and the previous 8 of the quenched limits do not affect very strongly the extraction of values of the axion mass, owing to the strong T -dependence of χ(T ). Using K = 9.22 and again assuming at least ten percent axion contribution to dark matter we get the range 40 µeV m A 930 µeV 4 In the scenario where inflation homogenizes one particular value of θ across our universe, a prediction of the axion DM abundance is not possible, except perhaps in an anthropic sense [4,31] while using K = 1 we would get 50 µeV m A 1100 µeV, i.e. only a ∼ 20% correction for an O(10) K factor uncertainty.