Axion phenomenology and $\theta$-dependence from $N_f = 2+1$ lattice QCD

We investigate the topological properties of $N_f = 2+1$ QCD with physical quark masses, both at zero and finite temperature. We adopt stout improved staggered fermions and explore a range of lattice spacings $a \sim 0.05 - 0.12$ fm. At zero temperature we estimate both finite size and finite cut-off effects, comparing our continuum extrapolated results for the topological susceptibility $\chi$ with predictions from chiral perturbation theory. At finite temperature, we explore a region going from $T_c$ up to around $4\, T_c$, where we provide continuum extrapolated results for the topological susceptibility and for the fourth moment of the topological charge distribution. While the latter converges to the dilute instanton gas prediction the former differs strongly both in the size and in the temperature dependence. This results in a shift of the axion dark matter window of almost one order of magnitude with respect to the instanton computation.


Introduction
Axions are among the most interesting candidates for physics beyond the Standard Model. Their existence has been advocated long ago [1][2][3][4] as a solution to the so-called strong-CP problem through the Peccei-Quinn (PQ) mechanism. It was soon realized that they could also explain the observed dark matter abundance of the visible Universe [5][6][7]. However, a reliable computation of the axion relic density requires a quantitative estimate of the parameters entering the effective potential of the axion field, in particular its mass and self-couplings as a function of the temperature T of the thermal bath.
The purpose of this study is to obtain predictions from the numerical simulations of Quantum Chromodynamics (QCD) on a lattice. Our results, which are summarized at the end of this section, suggest a possible shift of the axion dark matter window by almost one order of magnitude with respect to instanton computations. This shift is a consequence of the much slower decrease of the axion mass with the temperature in comparison to the dilute instanton gas prediction. Our present simulations are however limited to a range of temperatures not exceeding 600 MeV: the main obstruction is represented by the freezing of the topological modes on fine lattices, which afflicts present lattice QCD algorithms. For a more complete understanding of axion dynamics at finite T , in the future a major effort must be undertaken to reach higher temperatures.

General framework
Given the strong bounds on its couplings, the axion field can be safely treated as a nondynamical external field. Its potential is completely determined by the dependence of the QCD partition function on the θ-angle, which enters the pure gauge part of the QCD Euclidean Lagrangian as where is the topological charge density. The θ-dependent part of the free energy density can be parametrized as follows where χ(T ) is the topological susceptibility at θ = 0, is the global topological charge and V = V /T ), while s(θ, T ) is a dimensionless even function of θ such that s(0, T ) = 1. The quadratic term in θ, χ(T ), is proportional to the axion mass, while non-linear corrections in θ 2 , contained in s(θ, T ), provide information about axion interactions. In particular, assuming analyticity around θ = 0, s(θ, T ) can be expanded as follows [8] s(θ, T ) = 1 + b 2 (T )θ 2 + b 4 (T )θ 4 + · · · , (1.5) where the coefficients b n are proportional to the cumulants of the topological charge distribution. For instance b 2 , which is related to quartic interactions terms in the axion potential, can be expressed as (1.6) The function F(θ, T ), related to the topological properties of QCD, is of non-perturbative nature and hence not easy to predict reliably with analytic methods. This is possible only in some specific regimes: chiral perturbation theory (ChPT) represents a valid approach only in the low temperature phase; at high-T , instead, a possible analytic approach is the Dilute Instanton Gas Approximation (DIGA). DIGA predictions can in fact be classified in two groups: those that make only use of the DIGA hypothesis itself (i.e. that just relies on the existence of weakly interacting objects of topological charge one), and those that exploit also perturbation theory, the latter being expected to hold only at asymptotically high values of T . Using only the dilute gas approximation one can show that the θ-dependence of the free energy is of the form (see e.g. [9,10]) and thus b DIGA 2 = −1/12, b DIGA 4 = 1/360 and so on. Using also perturbation theory it is possible to obtain an explicit form for the dependence of the topological susceptibility on the temperature. To leading order, for N (l) f light quark flavors of mass m l , one obtains (see e.g. [9,10]) Only part of the NLO corrections to this expression are known (see [11] or [12] for a summary, [13] for the N f = 0 case).
As an alternative, a fully non-perturbative approach, which is based completely on the first principles of the theory, is represented by lattice QCD simulations. In fact, extensive studies have been carried out regarding the θ-dependence of pure gauge theories. It was shown in Ref. [14], and later confirmed in Refs. [13,15,16], that the form of the free energy in Eq. (1.7) describes with high precision the physics of the system for T 1.15 T c , while for T T c everything is basically independent of the temperature, thus strengthening the conclusion χ(T < T c ) ≈ χ(T = 0) obtained in previous studies [17][18][19][20][21]. In Refs. [13,22] it was also shown that the temperature dependence of the topological susceptibility is correctly reproduced by Eq. (1.8) for temperatures just above T c , even if the overall normalization is about a factor ten larger than the perturbative prediction.
A realistic study of θ-dependence aimed at being relevant to axion phenomenology requires the numerical simulation of lattice QCD including dynamical quarks with physical masses. Apart from the usual computational burden involved in the numerical simulation of light quarks, that represents a challenge from at least two different but interrelated points of view.
Because of the strict connection, in the presence of light fermions, between the topological content of gauge configurations and the spectrum of the fermion matrix (in particular regarding the presence of zero modes), a reliable study of topological quantities requires a discretization of the theory in which the chiral properties of fermions fields are correctly implemented. For standard discretizations, such properties are recovered only for small enough lattice spacings, so that a careful investigation of the continuum limit becomes essential. Indeed only recently it was possible to measure the dependence of the topological susceptibility on the quark masses to a sufficient accuracy to be compared with the prediction of chiral perturbation theory [23][24][25][26][27].
On the other hand, it is well known that, as the continuum limit is approached, it becomes increasingly difficult to correctly sample the topological charge distribution, because of the large energy barriers developing between configurations belonging to different homotopy classes, i.e. to different topological sectors, which can be hardly crossed by standard algorithms [28][29][30][31][32]. That causes a loss of ergodicity which, in principle, can spoil any effort to approach the continuum limit itself.
Combining these two problems together, the fact that a proper window exists, in which the continuum limit can be taken and still topological modes can be correctly sampled by current state-of-the-art algorithms, is highly non-trivial. In a finite temperature study, since the equilibrium temperature is related to the inverse temporal extent of the lattice by T = 1/(N t a), the fact that one cannot explore arbitrarily small values of the lattice spacing a, because of the above mentioned sampling problem, limits the range of explorable temperatures from above.
In the present study we show that, making use of current state-of-the-art algorithms, one can obtain continuum extrapolated results for the θ-dependence of QCD at the physical point, in a range of temperatures going up to about 4 T c , where T c ∼ 155 MeV is the pseudocritical temperature at which chiral symmetry restoration takes place. Then we discuss the consequences of such results to axion phenomenology in a cosmological context.
Our investigation is based on the numerical simulation of N f = 2 + 1 QCD, adopting a stout staggered fermions discretization with physical values of the quark masses and the tree level improved Symanzik gauge action. First, we consider simulations at zero temperature and various different values of the lattice spacing, in a range ∼ 0.05 − 0.12 fm and staying on a line of constant physics, in order to identify a proper scaling window where the continuum limit can be taken without incurring in severe problems with the freezing of topological modes. Results are then successfully compared to the predictions of chiral perturbation theory. We also show that, for lattice spacings smaller than those explored by us, the freezing problem becomes severe, making the standard Rational Hybrid Monte-Carlo (RHMC) algorithm useless.
For a restricted set of lattice spacings, belonging to the scaling window mentioned above, we perform finite temperature simulations, obtaining continuum extrapolated results for χ and b 2 in a range of T going up to around 600 MeV. These results are then taken as an input to fix the parameters of the axion potential in the same range of temperatures and perform a phenomenological analysis.

Summary of main results and paper organization
Our main results are the following. Up to T c the topological susceptibility is almost constant, and compatible with the prediction from ChPT. Above T c the value of b 2 rapidly converges to what predicted by DIGA computations; on the contrary the dependence of the topological susceptibility on T shows significant deviations from DIGA. This has a significant impact on axion phenomenology, in particular it results in a shift of the axion dark matter window by almost one order of magnitude with respect to the instanton computation. Since the T -dependence is much milder than expected from instanton calculations, it becomes crucial, for future studies, to investigate the system for higher values of T , something which also claims for the inclusion of dynamical charm quarks and for improved algorithms, capable to defeat or at least alleviate the problem of the freezing of topological modes.
The paper is organized as follows. In Section 2, we discuss the setup of our numerical simulations, in particular the lattice discretization adopted and the technique used to extract the topological content of gauge configurations. In Section 3 we present our numerical analysis and continuum extrapolated result for the θ-dependence of the free energy density, both at zero and finite temperature. Section 4 is dedicated to the analysis of the consequences of our results in the context of axion cosmology. Finally, in Section 5, we draw our conclusions and discuss future perspectives.  Table 1. Bare parameters used in this work, from [36,37] or spline interpolation of data thereof. The systematic uncertainty on the lattice spacing determination is 2 − 3% and the light quark mass is fixed by using m s /m l = 28.15.

Discretization adopted
The action of N f = 2 + 1 QCD is discretized by using the tree level improved Symanzik gauge action [33,34] and the stout improved staggered Dirac operator. Explicitly, the Euclidean partition function is written as 3) The symbol W n×m i; µν denotes the trace of the n × m Wilson loop built using the gauge links departing from the site in position i along the positive µ, ν directions. The gauge matrices U (2) i,µ , used in the definition of the staggered Dirac operator, are related to the gauge links U i; µ (used in S Y M ) by two levels of stout-smearing [35] with isotropic smearing parameter ρ = 0.15.
The bare parameters β, m s and m l ≡ m u = m d were chosen in such a way to have physical pion mass m π ≈ 135 MeV and physical m s /m l ratio. The values of the bare parameters used in this study to move along this line of constant physics are reported in Tab. 1. Most of them were determined in [36,37], the remaining have been extrapolated by using a cubic spline interpolation. The lattice spacings reported in Tab. 1 have a 2 − 3% of systematic uncertainty, as discussed in [36,37] (see also [38]).

Determination of the topological content
In order to expose the topological content of the gauge configurations, we adopt a gluonic definition of the topological charge density, measured after a proper smoothing procedure, which has been shown to provide results equivalent to definitions based on fermionic operators [39][40][41][42][43]. The basic underlying idea is that, close enough to the continuum limit, the topological content of gauge configurations becomes effectively stable under a local minimization of the gauge action, while ultraviolet fluctuations are smoothed off.
A number of smoothing algorithms has been devised along the time. A well known procedure is cooling [44]: an elementary cooling step consists in the substitution of each link variable U i; µ with the SU (3) matrix that minimizes the Wilson action, in order to reduce the UV noise. While for the case of the SU (2) group this minimization can be performed analytically, in the SU (3) case the minimization is usually performedà la Cabibbo-Marinari, i.e. by iteratively minimizing on the SU (2) subgroups. When this elementary cooling step is subsequently applied to all the lattice links we obtain a cooling step.
Possible alternatives may consist in choosing a different action to be minimized during the smoothing procedure, or in performing a continuous integration of the minimization equations. The latter procedure is known as the gradient flow [45,46] and has been shown to provide results equivalent to cooling regarding topology [47][48][49][50]. Because of its computational simplicity in this work we will thus use cooling, however it will be interesting to consider in future studies also the gradient flow, especially as an independent way to fix the physical scale [46,51,52].
Since the topological charge will be measured on smoothed configurations, we can use the simplest discretization of the topological charge density with definite parity [53] where U x; µν is the plaquette located in x and directed along the µ, ν directions. The tensorǫ is the completely antisymmetric tensor which coincides with the usual Levi-Civita tensor ǫ µνρσ for positive indices and is defined byǫ µνρσ = −ǫ (−µ)νρσ and antisymmetry for negative indices. The lattice topological charge Q L = x q L (x) is not in general an integer, although its distribution gets more and more peaked around integer values as the continuum limit is approached. In order to assign to a given configuration an integer value of the topological charge we will follow the prescription introduced in [29]: the topological charge is defined as Q = round(αQ L ), where 'round' denotes the rounding to the closest integer and α is fixed by the minimization of i.e. in such a way that αQ L is 'as integer' as possible. Actually, one could take the nonrounded definition Q = Q L as well, the only difference being a different convergence of results to the common continuum limit (see Ref. [54] for a more detailed discussion on this point). The topological susceptibility is then defined by s is the four-dimensional volume of the lattice, and the coefficient b 2 has been introduced in Eq. (1.6).
We measured the topological charge every 20 cooling steps up to a maximum of 120 steps and verified the stability of the topological quantities under smoothing. Results that will be presented in the following have been obtained by using 100 cooling steps, a number large enough to clearly identify the different topological sectors but for which no significant signals of tunneling to the trivial sector are distinguishable in mean values.
For all data reported in the following we verified that the corresponding time histories were long enough to correctly sample the topological charge. In particular we checked that Q is compatible with zero and that the topological charge is not frozen. Indeed it is well known that, while approaching the continuum limit, the autocorrelation time of the topological charge increases very steeply until no tunneling events between different sectors happen anymore [28][29][30][31][32]. An example of this behavior can be observed in Fig. 1, where some time-histories for zero temperature runs are showed, for three different lattice spacings. While the general features of this phenomenon are common to all the lattice actions, the critical value of the lattice spacing at which the charge gets stuck may depend on the specific discretization adopted. In our case we were able to obtain reasonable sampling for lattice spacings down to a = 0.0572 fm, while finer lattices (in particular, with a 0.04 fm) showed severe freezing over thousands of trajectories and had to be discarded.

Numerical Results
The main purpose of our numerical study is to provide results for the θ-dependence of QCD at finite temperature, in particular above the pseudo-critical temperature, in order to take them as an input for axion phenomenology. However, in order to make sure that our lattice discretization is accurate enough to reproduce the chiral properties of light fermions, we have also performed numerical simulations at zero or low temperature, where results can be compared to reliable analytic predictions.
Indeed, at zero temperature, the full θ dependence of the QCD partition function, can be estimated reliably using chiral Lagrangians [4,55,56]. In particular, at leading order in the expansion, χ and b 2 can be written just in terms of m π , f π and z = m u /m d as  ChPT provides a prediction at finite temperature as well. In particular we have where the functions J n are defined in Ref. [59]. However, at temperatures around and above T c , the chiral condensate drops and chiral Lagrangians break down, so that the finite T predictions in Eq. (3.2) are expected to fail. In this regime non-perturbative computations based on first principle QCD are mandatory.

Zero Temperature
In Tab. 2 we report our determinations of the topological susceptibility on N 4 s lattices for several values of N s and different lattice spacings. The results on the three smallest lattice spacings are also plotted in Fig. 2. In order to extract the infinite volume limit of χ at fixed lattice spacing, one can either consider only results obtained on the largest available lattices and fit them to a constant, or try to model the behavior of χ with the lattice size using all available data. We have followed both procedures in order to better estimate possible systematics.
The topological susceptibility can be written as the integral of the two point function of the topological charge density; since the η ′ is the lightest intermediate state that significantly contributes to this two point function (three pions states are OZI suppressed), one expects the leading asymptotic behavior to be where C is an unknown constant. This form nicely fits data in the whole available range of aN s (see Fig. 2), both when m η ′ is fixed to its physical value and when it is left as a fit parameter. The results obtained using this last procedure (which gives the most conservative estimates) are reported in Tab. 2 and correspond to the entries denoted by N s = ∞. On the other hand, a fit to a constant value works well when using only data coming from lattices with aN s > 1.5 fm, and provides consistent results within errors. This analysis makes also us confident that results obtained on the lattices with a ≃ 0.1249 fm and a ≃ 0.0989 fm, where a single spatial volume is available with aN s > 3 fm, are not affected by significant finite size effects. In Fig. 3 we plot χ 1/4 , extrapolated to infinite volume, against the square of the lattice spacing, together with the ChPT prediction. Finite cut-off effects are significant, especially for a 0.1 fm, meaning that we are not close enough to the continuum limit to reproduce the correct chiral properties of light quarks. In the case of staggered quarks, such lattice artifacts originate mostly from the fact that the full chiral symmetry group is reproduced only in the continuum limit, so that the pion spectrum is composed of a light pseudo-Goldstone boson and other massive states which become degenerate only as a → 0. The physical signal vanishes in the chiral limit whereas this is not the case for its discretization effects. This means that it is necessary to work on very fine lattice spacings in order to keep these effects under control, a task which is particularly challenging with present algorithms, because of the freezing of the topological charge.
Despite the large cut-off effects, we can perform the continuum extrapolation of our data. If one considers only the three finest lattice spacings, the finite cut-off effects for this quantity can be parametrized as simple O(a 2 ) corrections and a best fit yields χ 1/4 = 73(9) MeV, while in order to describe the whole range of available data one must take into account O(a 4 ) corrections as well, i.e.
obtaining χ 1/4 = 69 (9) MeV. Both best fits are reported in Fig. 3. Therefore we conclude that the continuum extrapolation is already under control with the available lattice data, and is in satisfactory agreement with the predictions of chiral perturbation theory. In order to further inquire about the reliability of our continuum extrapolation and the importance of the partial breaking of chiral symmetry in the staggered discretization, we studied the combination  m ngb (a) → m phys π , this ratio converges to χ 1/4 as a → 0. The state with taste structure γ i γ µ was used, whose mass is close to the root mean square value of all the taste states (see Ref. [37] Fig. 2) and the values of χ 1/4 tc (a) are shown in Fig. 3 as square points. It is clear that the combination χ 1/4 tc (a) strongly reduces lattice artefacts with respect to χ 1/4 (a), moreover a linear fit in a 2 well describes data for all available lattice spacings, giving the result χ 1/4 = 77(3). Although a complete study of the systematics affecting χ 1/4 tc (a) was not performed (e.g. the dependence of m ngb (a) on the lattice size was not studied, just the largest size was used), this is a strong indication that the dominant source of lattice artefacts in χ 1/4 (a) is the chiral symmetry breaking present at finite lattice spacing in the staggered discretization.

Finite Temperature
Finite temperature simulations can in principle be affected by lattice artifacts comparable to those present at T = 0. For that reason, we have limited our finite T simulations to the three smallest lattice spacings explored at T = 0, i.e. those in the scaling window adopted   for the extrapolation to the continuum limit with only O(a 2 ) corrections. At fixed lattice spacing the temperature has been varied by changing the temporal extent N t of the lattice and, in all cases, we have fixed N s = 48. This gives a spatial extent equal or larger than those explored at T = 0 and an aspect ratio N s /N t ≥ 3 for all explored values of N t . The absence of significant finite volume effects has also been verified directly by comparing results with those obtained on N t × 40 3 lattices. In Tab. 3 we report the numerical values obtained for the topological susceptibility, for the ratio χ(T, a)/χ(T = 0, a) (where for χ(T = 0, a) the infinite volume extrapolation has been taken) and for the cumulant ratio b 2 .
Results for χ 1/4 as a function of T /T c , where T c = 155 MeV, are reported in Fig. 4 for the three different lattice spacings. The dependence on a is quite strong, as expected from  the T = 0 case. Inspired by the instanton gas prediction, Eq. (1.8), we have performed a fit with the following ansatz which also takes into account the dependence on a and nicely describes all data in the range T > 1.2 T c with χ 2 /dof ≃ 0.7. In Tab. 4 we report the best fit values obtained performing the fit in the region T > T cut for some T cut values; best fit curves are reported in Fig. 4 as well, together with a band corresponding to the continuum extrapolation.  It is remarkable that most of the lattice artifacts disappear when one considers, in place of χ itself, the ratio χ(T, a)/χ(T = 0, a), whose dependence on the lattice spacing is indeed quite mild. That is clearly visible in Fig. 5. Also in this case we have adopted a fit ansatz similar to Eq. (3.7) which well fits all data in the range T > 1.2 T c with a χ 2 /dof ≃ 1. The best fit curves and the continuum extrapolation are shown in Fig. 5. The best fit parameters, again for   different T cut values, are reported in Tab. 5. It is important to note, in order to assess the reliability of our continuum limit, that the two different extrapolations, Eqs. (3.7) and (3.8), provide perfectly consistent results: that can be appreciated both from Fig. 6, where we compare the two coefficients D 2 and 4A 2 , and from Fig. 7, where we directly compare the two continuum extrapolations. As a further check, we have also verified that by fitting [χ(a, T )/χ(a, T = 0)] 1/4 we obtain results in perfect agreement with the ones obtained by  fitting χ(a, T )/χ(a, T = 0). As an example, using in the fit the values corresponding to T > 1.2T c one obtains for the exponent the value −0.674 (38), to be compared with D 2 /4 from the first line of Tab. 5 and A 2 from the first line of Tab. 4. However in the following analysis we will refer to results obtained through the ratio χ(T )/χ(0), which is the quantity showing smaller finite cut-off corrections.
Let us now comment on our results for the topological susceptibility. For temperatures below or around T c , the temperature dependence is quite mild and, for temperatures up to T ≃ 1.2 T c , even compatible with the prediction from ChPT, which is reported in Fig. 5 for comparison. Then, for higher values of T , a sharp power law drop starts. However, it is remarkable that the power law exponent is smaller, by more than a factor two, with respect to the instanton gas computation, Eq. (1.8), which predicts D 2 ≃ −8 in the case of three light flavors (the dependence on the number of flavor is however quite mild). DIGA is expected to provide the correct result in a region of asymptotically large temperatures, which however seems to be quite far from the range explored in the present study, which goes up to T ≃ 4 T c . This is in sharp contrast with the quenched theory, where the DIGA power law behavior sets in at lower temperatures [13,22]. In order to allow for a direct comparison, we have reported the DIGA prediction in Fig. 5, after fixing it by imposing χ(T c ) = χ(0) as an overall normalization. As we discuss in the following section, the much milder drop of χ as a function of T has important consequences for axion cosmology.
A determination of the topological susceptibility has been presented in Ref. [68], based on the Domain Wall discretization, in the range T /T c ∼ 0.9 -1.2. Since the data have been produced with different lattice spacings at different temperatures, it is however difficult to compare their results with ours. An extended range of T has been explored in Ref. [69] in the presence of twisted mass Wilson fermions, reporting a behavior similar to that observed in this study, although with larger values of the quark masses (corresponding to a pion mass m π ∼ 370 MeV). The comparison is performed in Fig. 8, and shows a reasonable agreement if results from Ref. [69] are rescaled according to the mass dependence expected from DIGA 2 , i.e. χ(T ) ∼ m 2 q ∼ m 4 π .
Let us now turn to a discussion of our results for the coefficient b 2 (defined in Eq. (1.5)), which is related to the non-gaussianity of the topological charge distribution and gives information about the shape of the θ-dependent part of the free energy density. Data for b 2 are reported in Fig. 9. For T < T c this observable is too noisy to give sensible results but a reasonable guess, motivated by what happens in the quenched case and by the ChPT prediction, is that b 2 is almost temperature independent for T T c , like the topological susceptibility. In the high temperature region b 2 can instead be measured with reasonable accuracy and, due to the peculiar dependence of the noise on this observable on χV, data on the smaller lattice spacing turned out to be significantly more precise than the others (see the discussion in Ref. [54]). While data clearly approach the dilute instanton gas prediction b DIGA 2 = −1/12 at high temperature, deviations from this value are clearly visible on all the lattice spacings for T ≃ 1.3 T c and, for the smallest lattice spacing, also up to T ∼ 2.5 T c . This is in striking contrast to what is observed in pure gauge theory, where deviations from −1/12 are practically absent for T 1.15 T c , with a precision higher than 10% [13][14][15][16]. As discussed in the introduction deviations from b DIGA 2 cannot be simply ascribed to a failure of perturbation theory (like e.g. for the behavior of χ(T )) but are instead unambiguous indications of interaction between instantons. Another difference with respect to the quenched case is that in the pure gauge theory (with both gauge groups SU (N ) and G 2 , see [14,15]) the asymptotic value is approached from below, while in the present case it is approached from above. These peculiar features can be related to a different interaction between instantons mediated by light quarks, as it is clear from the following discussion.
To describe deviations from the instanton gas behavior it is convenient to use the parametrization of F(θ, T ) introduced in Ref. [14]: Indeed it is not difficult to show that every even function of θ of period 2π can be written in this form, and the main advantage of this form is that the value of coefficient c 2n influences only b 2j with j ≥ n: in particular (3.10) Since the coefficients c 2n parametrize deviations from the instanton gas that manifest themselves only in higher-cumulants of the topological charge, it is natural to interpret Eq. (3.9) as a virial expansion, in which the role of the "density" is played by the first coefficient (i.e. the topological susceptibility χ), and to introduce the dimensionless coefficients d 2n by where χ(T = 0) was used just as a dimensional normalization and one expects a mild dependence of d 2(n−1) (T ) on the temperature, since the strongly dependent component χ(T ) n have been factorized. To lowest order we thus have , (3.12) which nicely describes the b 2 data for the smallest lattice spacing using d 2 = 0.80 (16), see Fig. 9, where the expression χ(T )/χ(T = 0) ≃ (T c /T ) 2.7 was used.
In the spirit of a virial expansion interpretation of Eq. (3.9), the coefficient d 2 can be considered as the lowest order interaction term between instantons. In particular, a positive value of d 2 corresponds to an attractive potential, which is in sharp contrast with the pure gauge case, where a repulsive, negative value of d 2 is observed. This peculiar difference can be surely interpreted in terms of effective instanton interactions mediated by light quarks, which are likely also responsible for the much slower convergence, with respect to the quenched case, to the DIGA prediction.

Implications for Axion Phenomenology
The big departure of the results for the topological susceptibility at finite temperature from the DIGA prediction has a strong impact on the computation of the axion relic abundance. In particular the model independent contribution from the misalignment mechanism is determined to be [59] Ω mis a = 86 33 where Ω γ and T γ are the present abundance and temperature of photons while n ⋆ a /s ⋆ is the ratio between the comoving number density n a = m a a 2 and the entropy density s computed at a late time t ⋆ such that the ratio n a /s became constant. The number density n a can be extracted by solving the axion equation of motion The temperature (and time) dependence of the Hubble parameter H is determined by the Friedmann equations and the QCD equation of state. The biggest uncertainties come therefore from the temperature dependence of the axion potential V (a). At high temperatures the Hubble friction dominates over the vanishing potential and the field is frozen to its initial value a 0 . As the Universe cools the pull from the potential starts winning over the friction (this happens when T ≈ T osc , defined as m a (T osc ) ≈ 3H(T osc )) and the axion starts oscillating around the minimum. Shortly after H becomes negligible and the mass term is the leading scale in Eq. (4.2). In this regime the approximate WKB solution has the form where R(t) is the cosmic scale factor. Since the energy density is given by ρ a ∼ m 2 a A 2 /2, the solution (4.3) implies that what is conserved in the comoving volume is not the energy density but N a = ρ a R 3 /m a , which can be interpreted as the number of axions [5][6][7]. Through the conservation of the comoving entropy S, it follows that n ⋆ a /s ⋆ becomes an adiabatic invariant. Hence, it is enough to integrate the equation of motion (4.2) in the small window around the time when T ≈ T osc . We integrated numerically Eq. (4.2) in the interval between the time when m a = H/10 to that corresponding to m a = 2400H and extract the ratio n ⋆ a /s ⋆ when m a ∼ 300H, namely a factor a hundred since the oscillation regime begins. The value for T osc varies from T c to several GeV depending on the axion decay parameter f a and the temperature dependence of the axion potential. More details about this standard computation can be found for example in [59,70]. In order to estimate the uncertainty in the results given below we varied the fitting parameters of the topological susceptibility D 2 , D 0 and those relative to the QCD equation of state [37] within the quoted statistical and systematic errors.
Given that b 2 (T ) converges relatively fast to the value predicted by a single cosine potential, we can assume V (a) = −χ(T ) cos(a/f a ) for T T c . Using the most conservative results for the fit of χ(T ), i.e. χ(T )/χ(0) = (1.8 ± 1.5)(T c /T ) 2.90±0.65 , in Fig. 10 we plot the prediction for the parameter f a as a function of the initial value of the axion field θ 0 = a 0 /f a assuming that the misalignment axion contribution make up for the whole observed dark matter abundance, Ω DM = 0.259(4) [71]. We also plot the case where the axion misalignment contribution accounts only for part (10% for definiteness) of the dark matter abundance. Figure 10. Values of the axion decay constant f a as a function of the initial field value θ 0 = a 0 /f a such that the axion misalignment contribution matches the full or a tenth of the observed dark matter abundance (red band or dotted green line respectively). When the PQ symmetry is broken only after inflation the axion abundance is reproduced by choosing θ 0 ≈ 2.2, i.e. the vertical blue dashed line.
In some cases the axion field acquires all possible values within the visible horizon, therefore the initial condition to the Eq. (4.2) needs to be integrated over. This happens if the PQ symmetry is broken only after inflation or if the PQ symmetric phase is temporarily restored after inflation (e.g. if the Hubble scale during inflation or the maximum temperature at reheating are at or above the PQ scale f a ). Numerically it corresponds to choosing θ 0 ≈ 2.2 and the value of the PQ scale in this case is fixed to: where the errors correspond respectively to the uncertainties on the fit parameters D 2 , D 0 , and to the error in the QCD equation of state. Note that, in this case, also other contributions from topological defects are expected to contribute. In case their effects are not negligible, or if axions are not responsible for the whole dark matter abundance, the value above should just be read as an upper bound to the PQ scale. The value in Eq. (4.4) is almost one order of magnitude bigger than the one computed with instantons techniques (f a 0.2 · 10 12 GeV), which in fact corresponds to the value we instead get for one tenth of the dark matter abundance.
Some few remarks are in order. While the uncertainties on the axion mass fit we used are not small, the final axion abundance is rather insensitive to them and the prediction for f a has therefore a good precision. The results above, however, rely on the extrapolation of the axion mass fit formula up to few GeV, where no lattice data is available. In particular for the value of f a in Eq. (4.4) the axion field starts oscillating around T osc = 4.3 GeV. An even longer extrapolation is required for f a = 1.67 · 10 11 GeV, corresponding to Ω mis a = 0.1 Ω DM , where the axion starts oscillating around T osc = 7.2 GeV.
Given the stability of the fit in the accessible window of temperatures (see Fig. 6) big changes in the axion mass behavior are not expected. Still, extending the analysis to even higher temperatures would be extremely useful to control better the extrapolation systematics.
Larger values of f a corresponds instead to smaller T osc , for example for f a = 10 13 GeV, T osc = 2.2 GeV while for f a = 10 14 GeV, T osc = 1.1 GeV. Above these values no extrapolation is required and the corresponding results are free from the extrapolation uncertainties.

Conclusions
We studied the topological properties and the θ-dependence for N f = 2 + 1 QCD along a line of constant physics, corresponding to physical quarks masses, and for temperatures up to 4 T c , where T c = 155 MeV is the pseudo-critical temperature at which chiral symmetry restoration takes place. We explored several lattice spacings, in a range 0.05 − 0.12 fm, in order to perform a continuum extrapolation of our results. Our investigation at even smaller lattice spacings has been hindered by a severe slowing down in the decorrelation of the topological charge.
At zero temperature we observe large cut-off effects for the topological susceptibility. Nevertheless we are able to perform a continuum extrapolation, obtaining from the three finest lattice spacings χ 1/4 = 73(9) MeV, in reasonable agreement with the ChPT prediction in the case of degenerate light flavors, χ 1/4 ChPT = 77.8(4) MeV. At finite temperature we observe that cut-off effects are drastically reduced when one considers the ratio χ(T )/χ(T = 0), which turns out to be the most convenient quantity to perform a continuum extrapolation. The agreement with ChPT persists up to around T c .
At higher temperatures the topological susceptibility presents the power-law dependence χ(T )/χ(0) = D 0 (T /T c ) D 2 , with D 2 ∼ −3. The exponent of the power-law behavior is however substantially smaller than the DIGA prediction, D 2 ∼ −8.
Regarding the shape of the θ-dependent part of the free energy density, and in particular the lowest order coefficient b 2 , a visible convergence to the instanton gas prediction is observed in the explored range. With respect to the pure gauge case, the convergence of b 2 to DIGA is slower and the deviation is opposite in sign [14,15]. This suggests a different interaction between instantons right above T c , namely repulsive in the quenched case and attractive in full QCD [72].
The deviations from the dilute instanton gas predictions that we found in the present study have a significant impact on axions, resulting in particular in a shift of the axion dark matter window by almost one order of magnitude with respect to estimates based on DIGA. The softer temperature dependence of the topological susceptibility also changes the onset of the axion oscillations, which would now start at a higher temperatures (T ∼ 4 GeV).
An important point is that this seems an effect directly related to the presence of light fermionic degrees of freedom: indeed, pure gauge computations [13,22] observe a power-law behavior in good agreement with DIGA in a range of temperatures comparable to those explored in the present study.
One might wonder whether a different power law behavior might set in at temperatures higher than 1 GeV. That claims for future studies extending the range explored by us. The main obstruction to this extension is represented by the freezing of the topological modes at smaller lattice spacings which would be necessary to investigate such temperatures (a < 0.05 fm). Such an obstruction could be overcome by the development of new Monte-Carlo algorithms. Proposals in this respect have been made in the past [73,74] and some new strategies have been put forward quite recently [75,76]. In view of the exploration of higher temperatures, one should also consider the inclusion of dynamical charm quarks.