Recent progress on QCD inputs for axion phenomenology

The properties of the QCD axion are strictly related to the dependence of strong interactions on the topological parameter theta. We present a determination of the topological properties of QCD for temperatures up to around 600 MeV, obtained by lattice QCD simulations with 2+1 flavors and physical quark masses. Numerical results for the topological susceptibility, when compared to instanton gas computations, differ both in size and in the temperature dependence. We discuss the implications of such findings for axion phenomenology, also in comparison to similar studies in the literature, and the prospects for future investigations.


Introduction
Axions have been advocated long ago [1][2][3][4] as a possible solution to the so-called strong-CP problem through the Peccei-Quinn (PQ) mechanism. It was soon realized that the axion could also be a possible source of dark matter [5][6][7]. That puts axions among the best candidates for physics beyond the Standard Model. A reliable computation of the axion relic density and of the present axion mass requires a quantitative estimate of the parameters entering the axion effective potential, i.e. its mass and self-coupling, as a function of the temperature.
Large part of the information is contained in the dependence of strong interactions on the topological parameter θ, which enters the pure gauge part of the QCD Euclidean Lagrangian as where q(x) is the topological charge density. The free energy density can be parametrized as follows where χ(T ) is the topological susceptibility at θ = 0, which is proportional to the axion mass: is the global topological charge and V = V/T is the 4-dimensional volume. The coefficients b n are proportional to the higher cumulants of the topological charge distribution, for instance they provide information about axion interactions.
At high-T , instead, a possible approach is based on the Dilute Instanton Gas Approximation (DIGA). Indeed, since instantons of size ρ ≫ 1/T are suppressed by thermal fluctuations, at asymptotically high-T one can use the perturbative 1-loop estimate for the instanton action, S e f f ≃ 8π 2 /g 2 (T ), which leads to the prediction of a path integral dominated by a dilute gas of weakly interacting objects of topological charge one. The θ-dependence of the free energy is of the form [13,14]) and, using the 1-loop effective action, one obtains [13,14]) for the theory with N f light flavors of mass m l , i.e. χ ∝ T −7.66 for N f = 2. Of course, the range of reliability of the perturbative prediction is not known a priori, and in principle one would trust it only for temperatures much larger than T c . A fully non-perturbative determination of θ-dependence can be obtained by lattice QCD simulations. The task is not easy: one needs to implement a lattice discretization of the theory which is close enough to the continuum limit to reproduce the correct chiral properties of fermion fields, however the correct sampling of topological modes is affected by critical autocorrelation times as the continuum limit is approached [15][16][17][18][19], and is even more difficult in the high T region, where topological fluctuations become very rare.
In the following we present a lattice study of the topological properties of N f = 2 + 1 QCD with physical quark masses. We adopt stout improved [20] staggered fermions and the tree level improved Symanzik action [21,22] for the pure gauge sector, exploring several lattice spacings, in a range ∼ 0.05 − 0.12 fm, and staying on a line of constant physics, which is fixed following the determinations reported in Ref. [23,24]; the topological content of gauge configurations is extracted by a gauge operator after proper smoothing of gauge fields [25][26][27][28][29][30][31]. More technical details can be found in Ref. [32].  (7)).
First we consider simulations at T = 0, in order to identify a proper scaling window and compare the continuum extrapolated results for the topological susceptibility to the ChPT prediction. At finite T , we explore a region going up to around 4 T c : we compare our results with DIGA predictions and with other recent results, discussing implications for axion phenomenology and prospects for future investigations.

Numerical Results at Zero and Finite Temperature
Results obtained for χ 1/4 at T = 0, after extrapolation to the infinite volume limit, are reported in Fig. 1 as a function of a 2 , together with the expectation from ChPT. The range of explored lattice spacings is limited on the left by freezing of topological modes that we have experienced for a < 0.05 fm, see Ref. [32].
Finite cut-off effects are significant and for a ≥ 0.1 fm the reported values are closer to the quenched case rather than to the ChPT prediction. Nevertheless, we can perform a consistent continuum extrapolation of our data: just O(a 2 ) corrections are needed to describe results at the three finest lattice spacings, while O(a 4 ) corrections must be taken into account to describe the whole range. We obtain consistent continuum extrapolated results (χ 1/4 = 73(9) MeV in the former case and χ 1/4 = 69(9) MeV in the latter), which are also in agreement with ChPT.
The large cutoff effects are likely due to the fact that, for those lattice spacings, our discretization fails to reproduce the correct chiral properties of light quarks. A typical manifestation of that, in the case of staggered quarks, is the fact that the full degenerate multiplet of (pseudo)Goldstone bosons is reproduced only in the continuum limit, while at finite lattice spacings one has one light pion plus other massive states which become degenerate only as a → 0. A confirmation of that comes as one studies the combination where m ngb (a) is one of the massive pions (in particular, we adopted the one with taste structure γ i γ µ ). The values of χ 1/4 tc (a) are reported in Fig. 1 as square points: cut-off effects are strongly reduced in this way and can be parametrized as O(a 2 ) corrections over the whole range.  Finite temperature simulations have been performed at the three smallest lattice spacings, i.e. those for which O(a 2 ) corrections suffice to describe the T = 0 data, and at various values of the temporal extent N t , with a fixed spatial extent N s = 48, after checking the absence of appreciable finite size effects. Results for χ 1/4 are reported in Fig. 2 as a function of T/T c (T c = 155 MeV) and have been fitted according to the following ansatz which is inspired to the DIGA prediction and takes into account also finite a corrections. If one takes into account the range T > 1.2 T c the ansatz works quite well (χ 2 /dof ≃ 0.7), however we obtain values for the exponent A 2 which are significantly lower (by roughly a factor 2) than the instanton gas prediction.
Results emerge more clearly as one considers the ratio χ(T, a)/χ(T = 0, a), which is reported in Fig. 3, together with the DIGA and ChPT predictions. Lattice artifacts seem to be strongly suppressed in this case, indeed the dependence on the lattice spacing is quite mild. The difference with the DIGA prediction is evident, indeed a best fit of our data to the function yields a slope coefficient D 2 ∼ −3, against a DIGA prediction D 2 ∼ −8. On the other hand, the two different continuum extrapolations, Fig. 2 and Fig. 3, lead to perfectly consistent results, which are compared in Fig. 4.

Implications for Axion Cosmology
Whereas many high energy models exist which introduce the PQ symmetry and its breaking, the common simplest form of the low energy effective lagrangian involving the axion field can be written    as follows: where a is the axion field acquiring a non-zero vacuum expectation value a after breaking of the PQ symmetry. It has only derivative terms, apart from a coupling to the QCD topological charge density, which is regulated by the constant f a . This coupling links θ dependence and axion phenomenology: indeed one has a minimum for the axion effective potential corresponding to a zero effective θ; moreover, assuming that f a is very large, a can be considered as quasi-static and the mass and interactions of the axion can be derived from the θ-dependence of QCD. For instance, one has the relation which would give the value of the axion mass today (from the value of χ at T = 0) once f a is fixed.
In order to fix f a , one needs further inputs. A possibility is to relate the present abundance of axions to the amount of observed dark matter. One of the best candidates for cosmological axion production is misalignement: at the time the PQ symmetry breaks, the axion field is not aligned along the minimum, θ e f f = 0; the resulting excitation is converted into axion production.
A precise description needs the solution of the axion field evolution which, in the static limit and taking into account only terms quadratic in θ in the axion potential, reads as follows where H is the Hubble constant. Since H is a decreasing function of time, while m a is increasing (assuming that the Universe temperature T decreases with t), Eq. (11) predicts a damped solution as long as the friction term wins over the pull from the potential, and becomes oscillatory after the Universe cools below a temperature T osc for which m a (T osc ) ≈ 3H(T osc ). Shortly after that time, the solution develops an adiabatic invariant corresponding to the total number of axions. The value of T osc depends both on χ(T ) and on f a , so that, by requiring a given amount of axion relics, one can fix the value of f a . Actually, an unknown constant is the initial misalignement angle, θ 0 : that takes a fixed value throughout the Universe if the PQ symmetry breaking takes place before inflation, otherwise it can acquire all possible values within the visible horizon and one needs to integrate over them. In Fig. 5, based on the results obtained for χ(T ) from our numerical simulations, we show our prediction (see Ref. [32] for more details) for f a and for the present value of the axion mass, depending on the different possible assumptions, i.e. pre-or post-inflationary scenario and relative amount of dark matter density Ω DM due to axions from the misalignment mechanism.
Our results lead to values of f a which are tipically one order of magnitude larger than the ones obtained from the dilute instanton gas model, hence to smaller values of the present axion mass. This  Figure 6. b 2 as a function of T , together with the ChPT prediction at T = 0 (−0.022) and the DIGA prediction (−1/12). The light blue band is the result of a fit to the smallest lattice spacing data using a virial expansion (see Ref. [32] for more details). is mostly due to the different behavior of χ(T ), which shows, at least in the explored range, a much milder dependence on T with respect to DIGA predictions. This result of course might be affected by residual systematic effects, in particular we have results at three different lattice spacings only in a region which goes up to T ∼ 300 MeV, while the oscillation temperature deriving from our data is of the order of few GeVs, so that we are strongly relying on an extrapolation. The main difficulty in approaching higher temperatures is given by the critical slowing down observed for small lattice spacings and by the fact that topological fluctuations become very rare at high T .
Recent results from different lattice investigations have shown contrasting results. While the results of Ref. [33] are in qualitative agreement with ours, other studies have reported a much better agreement with the instanton gas prediction even in the region right above T c [34][35][36].
In particular, the authors of Ref. [34] manage to determine χ(T ) for T up to a few GeVs and observe a slope consistent with DIGA shortly after T c . The main differences consist in a reweighting procedure adopted in Ref. [34], based on the lowest lying Dirac operator eigenvalues, and in a new strategy to approach the high T regime (see also Ref. [37]). This approach avoids direct simulations at high T , and is based instead on simulations at fixed topology, which are used to determine the relative weight, in the path integral, of the topological sectors Q = ±1 with respect to the topological sector 0, since those are believed to be the only ones relevant at high T . Of course such a strategy assumes right from the beginning that the instanton gas is dilute enough and non-interacting.
Good quantities to check for the diluteness hypothesis are the b n coefficients appearing in Eq. (2), which are fixed by the fact that, for a non-interacting instanton gas, one has F(θ) ∝ (1 − cos θ). Our results for b 2 are reported in Fig. 6: even if the statistical uncertainties are still large, one sees that deviations from the dilute gas hypothesis could be still appreciable for T up to 2 − 3 T c . This is in contrast with the pure gauge case [38][39][40][41], where corrections become negligible shortly after T c . That confirms the non-triviality of fermionic contributions and claims for future studies, which should further check existing results. Various improvements are possible, especially to deal with the critical slowing down of topological modes by means of improved algorithms [42][43][44][45][46].