# Axion phenomenology and -dependence from lattice QCD

###### Abstract

We investigate the topological properties of QCD with physical quark masses, both at zero and finite temperature. We adopt stout improved staggered fermions and explore a range of lattice spacings fm. At zero temperature we estimate both finite size and finite cut-off effects, comparing our continuum extrapolated results for the topological susceptibility with predictions from chiral perturbation theory. At finite temperature, we explore a region going from up to around , 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.

###### Keywords:

Lattice QCD, QCD Phenomenology, CP violation## 1 Introduction

Axions are among the most interesting candidates for physics beyond the Standard Model. Their existence has been advocated long ago Peccei:1977hh (); Peccei:1977ur (); Wilczek:1977pj (); Weinberg:1977ma () 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 Preskill:1982cy (); Abbott:1982af (); Dine:1982ah (). 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 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 , in the future a major effort must be undertaken to reach higher temperatures.

### 1.1 General framework

Given the strong bounds on its couplings, the axion field can be safely treated as a non-dynamical 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

(1) |

where

(2) |

is the topological charge density. The -dependent part of the free energy density can be parametrized as follows

(3) |

where is the topological susceptibility at ,

(4) |

( is the global topological charge and ), while is a dimensionless even function of such that . The quadratic term in , , is proportional to the axion mass, while non-linear corrections in , contained in , provide information about axion interactions. In particular, assuming analyticity around , can be expanded as follows Vicari:2008jw ()

(5) |

where the coefficients are proportional to the cumulants of the topological charge distribution. For instance , which is related to quartic interactions terms in the axion potential, can be expressed as

(6) |

The function , 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-, 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 . Using only the dilute gas approximation one can show that the -dependence of the free energy is of the form (see e.g. Gross:1980br (); Schafer:1996wv ())

(7) |

and thus , 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 light quark flavors of mass , one obtains (see e.g. Gross:1980br (); Schafer:1996wv ())

(8) |

where . Only part of the NLO corrections to this expression are known (see Morris:1984zi () or Ringwald:1999ze () for a summary, Borsanyi:2015cka () for the 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. Bonati:2013tt (), and later confirmed in Refs. Bonati:2015uga (); Borsanyi:2015cka (); Xiong:2015dya (), that the form of the free energy in Eq. (7) describes with high precision the physics of the system for , while for everything is basically independent of the temperature, thus strengthening the conclusion obtained in previous studies Alles:1996nm (); Alles:2000cg (); Gattringer:2002mr (); Lucini:2004yh (); DelDebbio:2004rw (). In Refs. Berkowitz:2015aua (); Borsanyi:2015cka () it was also shown that the temperature dependence of the topological susceptibility is correctly reproduced by Eq. (8) for temperatures just above , 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 Bazavov:2010xr (); Bazavov:2012xda (); Cichy:2013rra (); Bruno:2014ova (); Fukaya:2014zda ().

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 Alles:1996vn (); DelDebbio:2002xa (); DelDebbio:2004xh (); Schaefer:2010hu (); Kitano:2015fla (). 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 , the fact that one cannot explore arbitrarily small values of the lattice spacing , 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 , where MeV is the pseudo-critical 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 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 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 in a range of 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.

### 1.2 Summary of main results and paper organization

Our main results are the following. Up to the topological susceptibility is almost constant, and compatible with the prediction from ChPT. Above the value of rapidly converges to what predicted by DIGA computations; on the contrary the dependence of the topological susceptibility on 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 -dependence is much milder than expected from instanton calculations, it becomes crucial, for future studies, to investigate the system for higher values of , 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.

## 2 Numerical setup

### 2.1 Discretization adopted

The action of QCD is discretized by using the tree level improved Symanzik gauge action weisz (); curci () and the stout improved staggered Dirac operator. Explicitly, the Euclidean partition function is written as

(9) | ||||

(10) | ||||

(11) |

The symbol denotes the trace of the Wilson loop built using the gauge links departing from the site in position along the positive directions. The gauge matrices , used in the definition of the staggered Dirac operator, are related to the gauge links (used in ) by two levels of stout-smearing Morningstar:2003gk () with isotropic smearing parameter .

The bare parameters , and were chosen in such a way to have physical pion mass and physical 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 physline1 (); physline2 (), the remaining have been extrapolated by using a cubic spline interpolation. The lattice spacings reported in Tab. 1 have a of systematic uncertainty, as discussed in physline1 (); physline2 () (see also physline3 ()).

[fm] | ||
---|---|---|

3.750 | 0.1249 | |

3.850 | 0.0989 | |

3.938 | 0.0824 | |

4.020 | 0.0707 | |

4.140 | 0.0572 |

### 2.2 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 Neuberger:1997fp (); Hasenfratz:1998ri (); Luscher:1998pqa (); Luscher:2004fu (); Giusti:2008vb (). 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 cooling (): an elementary cooling step consists in the substitution of each link variable with the matrix that minimizes the Wilson action, in order to reduce the UV noise. While for the case of the group this minimization can be performed analytically, in the case the minimization is usually performed à la Cabibbo-Marinari, i.e. by iteratively minimizing on the 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 Luscher:2009eq (); Luscher:2010iy () and has been shown to provide results equivalent to cooling regarding topology Bonati:2014tqa (); Cichy:2014qta (); Namekawa:2015wua (); Alexandrou:2015yba (). 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 Luscher:2010iy (); Sommer:2014mea (); Fodor:2014cpa ().

Since the topological charge will be measured on smoothed configurations, we can use the simplest discretization of the topological charge density with definite parity DiVecchia:1981qi ()

(12) |

where is the plaquette located in 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 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 DelDebbio:2002xa (): the topological charge is defined as , where ‘’ denotes the rounding to the closest integer and is fixed by the minimization of

(13) |

i.e. in such a way that is ‘as integer’ as possible. Actually, one could take the non-rounded definition as well, the only difference being a different convergence of results to the common continuum limit (see Ref. Bonati:2015sqt () for a more detailed discussion on this point). The topological susceptibility is then defined by

(14) |

where is the four-dimensional volume of the lattice, and the coefficient has been introduced in Eq. (6).

We measured the topological charge every cooling steps up to a maximum of steps and verified the stability of the topological quantities under smoothing. Results that will be presented in the following have been obtained by using 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 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 Alles:1996vn (); DelDebbio:2002xa (); DelDebbio:2004xh (); Schaefer:2010hu (); Kitano:2015fla (). 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 , while finer lattices (in particular, with fm) showed severe freezing over thousands of trajectories and had to be discarded.

## 3 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 Weinberg:1977ma (); DiVecchia:1980ve (); Leutwyler:1992yt (). In particular, at leading order in the expansion, and can be written just in terms of , and as

(15) |

NLO corrections have been computed in Mao:2009sy (); Guo:2015oxa (); diCortona:2015ldu () and are of the order of percent for physical values of the
light quark masses. Using the estimate in diCortona:2015ldu () we have
MeV and for , while
MeV and for used here. Note that
for pure SU(3) Yang-Mills these numbers become MeV Alles:1996nm (); DelDebbio:2002xa (); DelDebbio:2004ns (); Durr:2006ky (); Luscher:2010ik (); Cichy:2015jra (); Ce:2015qha () and^{1}^{1}1It is interesting to
note the striking agreement between the ChPT prediction for in the case
of degenerate light flavors and the numerical results obtained for it in the
quenched theory Bonati:2015sqt (). That seems to suggest a similar form of
-dependence in the low temperature phases of the full and of the
quenched theory. At this stage we have no particular explanation for this
agreement, which might be considered as a coincidence.
Bonati:2015sqt () (see also
Refs. DelDebbio:2002xa (); D'Elia:2003gr (); Giusti:2007tu (); Panagopoulos:2011rb (); Ce:2015qha ()
for previous determinations in the literature).

ChPT provides a prediction at finite temperature as well. In particular we have

(16) | ||||

where the functions are defined in Ref. diCortona:2015ldu (). However, at temperatures around and above , the chiral condensate drops and chiral Lagrangians break down, so that the finite predictions in Eq. (16) are expected to fail. In this regime non-perturbative computations based on first principle QCD are mandatory.

### 3.1 Zero Temperature

In Tab. 2 we report our determinations of the topological susceptibility on lattices for several values of 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

(17) |

where is an unknown constant. This form nicely fits data in the whole available range of (see Fig. 2), both when 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 . On the other hand, a fit to a constant value works well when using only data coming from lattices with , and provides consistent results within errors. This analysis makes also us confident that results obtained on the lattices with fm and fm, where a single spatial volume is available with fm, are not affected by significant finite size effects.

In Fig. 3 we plot , extrapolated to infinite volume, against the square of the lattice spacing, together with the ChPT prediction. Finite cut-off effects are significant, especially for , 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 . 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 corrections

(18) |

and a best fit yields MeV, while in order to describe the whole range of available data one must take into account corrections as well, i.e.

(19) |

obtaining 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

(20) |

Here is the mass of one of the non-Goldstone pions, i.e. of a state that becomes massless in the chiral limit only if the continuum limit is performed; since in the continuum , this ratio converges to as . The state with taste structure was used, whose mass is close to the root mean square value of all the taste states (see Ref. physline2 () Fig. 2) and the values of are shown in Fig. 3 as square points. It is clear that the combination strongly reduces lattice artefacts with respect to , moreover a linear fit in well describes data for all available lattice spacings, giving the result . Although a complete study of the systematics affecting was not performed (e.g. the dependence of 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 is the chiral symmetry breaking present at finite lattice spacing in the staggered discretization.

[fm] | ||

32 | 0.1249 | |

32 | 0.0989 | |

12 | 0.0824 | |

16 | 0.0824 | |

20 | 0.0824 | |

24 | 0.0824 | |

32 | 0.0824 | |

0.0824 | ||

16 | 0.0707 | |

20 | 0.0707 | |

24 | 0.0707 | |

32 | 0.0707 | |

40 | 0.0707 | |

0.0707 | ||

20 | 0.0572 | |

24 | 0.0572 | |

32 | 0.0572 | |

40 | 0.0572 | |

48 | 0.0572 | |

0.0572 |

### 3.2 Finite Temperature

Finite temperature simulations can in principle be affected by lattice artifacts comparable to those present at . For that reason, we have limited our finite simulations to the three smallest lattice spacings explored at , i.e. those in the scaling window adopted for the extrapolation to the continuum limit with only corrections. At fixed lattice spacing the temperature has been varied by changing the temporal extent of the lattice and, in all cases, we have fixed . This gives a spatial extent equal or larger than those explored at and an aspect ratio for all explored values of . The absence of significant finite volume effects has also been verified directly by comparing results with those obtained on lattices. In Tab. 3 we report the numerical values obtained for the topological susceptibility, for the ratio (where for the infinite volume extrapolation has been taken) and for the cumulant ratio .

[fm] | |||||

16 | 0.0824 | 0.964 | -0.23(33) | ||

14 | 0.0824 | 1.102 | -0.191(95) | ||

12 | 0.0824 | 1.285 | -0.003(25) | ||

10 | 0.0824 | 1.542 | -0.083(29) | ||

8 | 0.0824 | 1.928 | -0.065(25) | ||

24 | 0.0707 | 0.749 | -0.11(12) | ||

16 | 0.0707 | 1.124 | -0.057(43) | ||

14 | 0.0707 | 1.284 | -0.019(30) | ||

12 | 0.0707 | 1.498 | -0.101(21) | ||

10 | 0.0707 | 1.798 | -0.096(26) | ||

8 | 0.0707 | 2.247 | -0.077(18) | ||

6 | 0.0707 | 2.996 | -0.076(10) | ||

16 | 0.0572 | 1.389 | -0.049(17) | ||

14 | 0.0572 | 1.587 | -0.058(16) | ||

12 | 0.0572 | 1.852 | -0.0626(90) | ||

10 | 0.0572 | 2.222 | -0.0705(74) | ||

8 | 0.0572 | 2.777 | -0.091(11) | ||

6 | 0.0572 | 3.703 | -0.0792(12) |

(MeV) | (fm) | ||
---|---|---|---|

1.2 | 76.0(5.1) | 103(17) | -0.671(40) |

1.4 | 80.2(6.9) | 100(20) | -0.728(57) |

1.6 | 86(10) | 85(24) | -0.750(80) |

2.0 | 83(18) | 97(53) | -0.718(14) |

(fm) | |||
---|---|---|---|

1.2 | 1.17(24) | -2(28) | -2.71(15) |

1.4 | 1.56(39) | -12(30) | -3.02(23) |

1.6 | 1.81(76) | -33(32) | -2.99(39) |

2.0 | 1.8(1.5) | 2(22) | -2.90(65) |

Results for as a function of , where MeV, are reported in Fig. 4 for the three different lattice spacings. The dependence on is quite strong, as expected from the case. Inspired by the instanton gas prediction, Eq. (8), we have performed a fit with the following ansatz

(21) |

which also takes into account the dependence on and nicely describes all data in the range with . In Tab. 4 we report the best fit values obtained performing the fit in the region for some 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 , 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. (21)

(22) |

which well fits all data in the range with a . The best fit curves and the continuum extrapolation are shown in Fig. 5. The best fit parameters, again for different 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. (21) and (22), provide perfectly consistent results: that can be appreciated both from Fig. 6, where we compare the two coefficients and , and from Fig. 7, where we directly compare the two continuum extrapolations. As a further check, we have also verified that by fitting we obtain results in perfect agreement with the ones obtained by fitting . As an example, using in the fit the values corresponding to one obtains for the exponent the value , to be compared with from the first line of Tab. 5 and from the first line of Tab. 4. However in the following analysis we will refer to results obtained through the ratio , 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 , the temperature dependence is quite mild and, for temperatures up to , even compatible with the prediction from ChPT, which is reported in Fig. 5 for comparison. Then, for higher values of , 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. (8), which predicts 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 .

This is in sharp contrast with the quenched theory, where the DIGA power law behavior sets in at lower temperatures Borsanyi:2015cka (); Berkowitz:2015aua (). In order to allow for a direct comparison, we have reported the DIGA prediction in Fig. 5, after fixing it by imposing as an overall normalization. As we discuss in the following section, the much milder drop of as a function of has important consequences for axion cosmology.

A determination of the topological susceptibility has been presented
in Ref. Buchoff:2013nra (), based on the Domain Wall
discretization, in the range 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 has been explored in
Ref. Trunin:2015yda () 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 MeV). The comparison is performed in
Fig. 8, and shows a reasonable agreement if results
from Ref. Trunin:2015yda () are rescaled according to the mass
dependence expected from DIGA^{2}^{2}2One might wonder why DIGA
should work for the mass dependence of and not for its
dependence on . A possible explanation is that, while the
temperature dependence stems from a perturbative computation, the
mass dependence comes from the existence of isolated zero modes in
the Dirac operator, i.e. from the very hypothesis of the existence
of a dilute gas of instantons, which seems to be verified already at
moderately low values of , see the following discussion on ,
but not at ., i.e. .

Let us now turn to a discussion of our results for the coefficient (defined in Eq. (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 are reported in Fig. 9. For 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 is almost temperature independent for , like the topological susceptibility. In the high temperature region can instead be measured with reasonable accuracy and, due to the peculiar dependence of the noise on this observable on , data on the smaller lattice spacing turned out to be significantly more precise than the others (see the discussion in Ref. Bonati:2015sqt ()). While data clearly approach the dilute instanton gas prediction at high temperature, deviations from this value are clearly visible on all the lattice spacings for and, for the smallest lattice spacing, also up to .

This is in striking contrast to what is observed in pure gauge theory, where deviations from are practically absent for , with a precision higher than Bonati:2013tt (); Borsanyi:2015cka (); Bonati:2015uga (); Xiong:2015dya (). As discussed in the introduction deviations from cannot be simply ascribed to a failure of perturbation theory (like e.g. for the behavior of ) 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 and , see Bonati:2013tt (); Bonati:2015uga ()) 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 introduced in Ref. Bonati:2013tt ():

(23) |

Indeed it is not difficult to show that every even function of of period can be written in this form, and the main advantage of this form is that the value of coefficient influences only with : in particular

(24) | ||||

Since the coefficients parametrize deviations from the instanton gas that manifest themselves only in higher-cumulants of the topological charge, it is natural to interpret Eq. (23) 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 by

(25) |

where was used just as a dimensional normalization and one expects a mild dependence of on the temperature, since the strongly dependent component have been factorized. To lowest order we thus have

(26) |

which nicely describes the data for the smallest lattice spacing using , see Fig. 9, where the expression was used.

In the spirit of a virial expansion interpretation of Eq. (23), the coefficient can be considered as the lowest order interaction term between instantons. In particular, a positive value of corresponds to an attractive potential, which is in sharp contrast with the pure gauge case, where a repulsive, negative value of 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.

## 4 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 diCortona:2015ldu ()

(27) |

where and are the present abundance and temperature of photons while is the ratio between the comoving number density and the entropy density computed at a late time such that the ratio became constant. The number density can be extracted by solving the axion equation of motion

(28) |

The temperature (and time) dependence of the Hubble parameter 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 . At high temperatures the Hubble friction dominates over the vanishing potential and the field is frozen to its initial value . As the Universe cools the pull from the potential starts winning over the friction (this happens when , defined as ) and the axion starts oscillating around the minimum. Shortly after becomes negligible and the mass term is the leading scale in Eq. (28). In this regime the approximate WKB solution has the form

(29) | ||||

where is the cosmic scale factor. Since the energy density is given by , the solution (29) implies that what is conserved in the comoving volume is not the energy density but , which can be interpreted as the number of axions Preskill:1982cy (); Abbott:1982af (); Dine:1982ah (). Through the conservation of the comoving entropy , it follows that becomes an adiabatic invariant. Hence, it is enough to integrate the equation of motion (28) in the small window around the time when . We integrated numerically Eq. (28) in the interval between the time when to that corresponding to and extract the ratio when , namely a factor a hundred since the oscillation regime begins. The value for varies from to several GeV depending on the axion decay parameter and the temperature dependence of the axion potential. More details about this standard computation can be found for example in Wantz:2009it (); diCortona:2015ldu (). In order to estimate the uncertainty in the results given below we varied the fitting parameters of the topological susceptibility , and those relative to the QCD equation of state physline2 () within the quoted statistical and systematic errors.

Given that converges relatively fast to the value predicted by a single cosine potential, we can assume for . Using the most conservative results for the fit of , i.e. , in Fig. 10 we plot the prediction for the parameter as a function of the initial value of the axion field assuming that the misalignment axion contribution make up for the whole observed dark matter abundance, Ade:2015xua (). We also plot the case where the axion misalignment contribution accounts only for part (10% for definiteness) of the dark matter abundance.

In some cases the axion field acquires all possible values within the visible horizon, therefore the initial condition to the Eq. (28) 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 ). Numerically it corresponds to choosing and the value of the PQ scale in this case is fixed to:

(30) |

where the errors correspond respectively to the uncertainties on the fit parameters , , 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. (30) is almost one order of magnitude bigger than the one computed with instantons techniques ( 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 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 in Eq. (30) the axion field starts oscillating around GeV. An even longer extrapolation is required for GeV, corresponding to , where the axion starts oscillating around 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 corresponds instead to smaller , for example for GeV, GeV while for GeV, GeV. Above these values no extrapolation is required and the corresponding results are free from the extrapolation uncertainties.

## 5 Conclusions

We studied the topological properties and the -dependence for QCD along a line of constant physics, corresponding to physical quarks masses, and for temperatures up to , where MeV is the pseudo-critical temperature at which chiral symmetry restoration takes place. We explored several lattice spacings, in a range 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 MeV, in reasonable agreement with the ChPT prediction in the case of degenerate light flavors, MeV.

At finite temperature we observe that cut-off effects are drastically reduced when one considers the ratio , which turns out to be the most convenient quantity to perform a continuum extrapolation. The agreement with ChPT persists up to around . At higher temperatures the topological susceptibility presents the power-law dependence , with . The exponent of the power-law behavior is however substantially smaller than the DIGA prediction, .

Regarding the shape of the -dependent part of the free energy density, and in particular the lowest order coefficient , a visible convergence to the instanton gas prediction is observed in the explored range. With respect to the pure gauge case, the convergence of to DIGA is slower and the deviation is opposite in sign Bonati:2013tt (); Bonati:2015uga (). This suggests a different interaction between instantons right above , namely repulsive in the quenched case and attractive in full QCD Callan:1977gz ().

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 ( 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 Borsanyi:2015cka (); Berkowitz:2015aua () 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 ( 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 deForcrand:1997fm (); openbc () and some new strategies have been put forward quite recently surfing (); endres (). In view of the exploration of higher temperatures, one should also consider the inclusion of dynamical charm quarks.

Numerical simulations have been performed on the Blue Gene/Q Fermi machine at CINECA, based on the agreement between INFN and CINECA (under INFN project NPQCD). Work partially supported by the INFN SUMA Project, by the ERC-2010 DaMESyFla Grant Agreement Number: 267985, by the ERC-2011 NEWPHYSICSHPC Grant Agreement Number: 27975 and by the MIUR (Italy) under a contract PRIN10-11 protocollo 2010YJ2NYW. FN acknowledges financial support from the INFN SUMA project.

## References

- (1) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
- (2) R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977).
- (3) F. Wilczek, Phys. Rev. Lett. 40 (1978) 279.
- (4) S. Weinberg, Phys. Rev. Lett. 40 (1978) 223.
- (5) J. Preskill, M. B. Wise and F. Wilczek, Phys. Lett. B 120 (1983) 127.
- (6) L. F. Abbott and P. Sikivie, Phys. Lett. B 120 (1983) 133.
- (7) M. Dine and W. Fischler, Phys. Lett. B 120 (1983) 137.
- (8) E. Vicari and H. Panagopoulos, Phys. Rept. 470, 93 (2009) [arXiv:0803.1593 [hep-th]].
- (9) D. J. Gross, R. D. Pisarski and L. G. Yaffe, Rev. Mod. Phys. 53, 43 (1981).
- (10) T. Schäfer and E. V. Shuryak, Rev. Mod. Phys. 70, 323 (1998) [hep-ph/9610451].
- (11) T. R. Morris, D. A. Ross and C. T. Sachrajda, Nucl. Phys. B 255, 115 (1985).
- (12) A. Ringwald and F. Schrempp, Phys. Lett. B 459, 249 (1999) [hep-lat/9903039].
- (13) S. Borsanyi et al., Phys. Lett. B 752, 175 (2016) [arXiv:1508.06917 [hep-lat]].
- (14) C. Bonati, M. D’Elia, H. Panagopoulos and E. Vicari, Phys. Rev. Lett. 110, 252003 (2013) [arXiv:1301.7640 [hep-lat]].
- (15) C. Bonati, JHEP 1503, 006 (2015) [arXiv:1501.01172 [hep-lat]].
- (16) G. Y. Xiong, J. B. Zhang, Y. Chen, C. Liu, Y. B. Liu and J. P. Ma, Phys. Lett. B 752, 34 (2016) [arXiv:1508.07704 [hep-lat]].
- (17) B. Alles, M. D’Elia and A. Di Giacomo, Nucl. Phys. B 494, 281 (1997) [Erratum Nucl. Phys. B 679, 397 (2004)] [hep-lat/9605013].
- (18) B. Alles, M. D’Elia and A. Di Giacomo, Phys. Lett. B 483, 139 (2000) [hep-lat/0004020].
- (19) C. Gattringer, R. Hoffmann and S. Schaefer, Phys. Lett. B 535, 358 (2002) [hep-lat/0203013].
- (20) B. Lucini, M. Teper and U. Wenger, Nucl. Phys. B 715, 461 (2005) [hep-lat/0401028].
- (21) L. Del Debbio, H. Panagopoulos and E. Vicari, JHEP 0409, 028 (2004) [hep-th/0407068].
- (22) E. Berkowitz, M. I. Buchoff and E. Rinaldi, Phys. Rev. D 92, 034507 (2015) [arXiv:1505.07455 [hep-ph]].
- (23) A. Bazavov et al. [MILC Collaboration], Phys. Rev. D 81, 114501 (2010) [arXiv:1003.5695 [hep-lat]].
- (24) A. Bazavov et al. [MILC Collaboration], Phys. Rev. D 87, 054505 (2013) [arXiv:1212.4768 [hep-lat]].
- (25) K. Cichy et al. [ETM Collaboration], JHEP 1402, 119 (2014) [arXiv:1312.5161 [hep-lat]].
- (26) M. Bruno, S. Schaefer and R. Sommer [ALPHA Collaboration], JHEP 1408, 150 (2014) [arXiv:1406.5363 [hep-lat]].
- (27) H. Fukaya et al. [JLQCD Collaboration], PoS LATTICE 2014, 323 (2014) [arXiv:1411.1473 [hep-lat]].
- (28) B. Alles, G. Boyd, M. D’Elia, A. Di Giacomo and E. Vicari, Phys. Lett. B 389, 107 (1996) [hep-lat/9607049].
- (29) L. Del Debbio, H. Panagopoulos and E. Vicari, JHEP 0208, 044 (2002) [hep-th/0204125].
- (30) L. Del Debbio, G. M. Manca and E. Vicari, Phys. Lett. B 594, 315 (2004) [hep-lat/0403001].
- (31) S. Schaefer et al. [ALPHA Collaboration], Nucl. Phys. B 845, 93 (2011) [arXiv:1009.5228 [hep-lat]].
- (32) R. Kitano and N. Yamada, JHEP 1510, 136 (2015) [arXiv:1506.00370 [hep-ph]].
- (33) P. Weisz, Nucl. Phys. B 212, 1 (1983).
- (34) G. Curci, P. Menotti and G. Paffuti, Phys. Lett. B 130, 205 (1983) [Erratum-ibid. B 135, 516 (1984)].
- (35) C. Morningstar and M. J. Peardon, Phys. Rev. D 69, 054501 (2004) [hep-lat/0311018].
- (36) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 0906, 088 (2009) [arXiv:0903.4155 [hep-lat]].
- (37) S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, JHEP 1011, 077 (2010) [arXiv:1007.2580 [hep-lat]];
- (38) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabo, Phys. Lett. B 730, 99 (2014) [arXiv:1309.5258 [hep-lat]].
- (39) H. Neuberger, Phys. Lett. B 417, 141 (1998) [hep-lat/9707022].
- (40) P. Hasenfratz, V. Laliena and F. Niedermayer, Phys. Lett. B 427, 125 (1998) [hep-lat/9801021].
- (41) M. Luscher, Phys. Lett. B 428, 342 (1998) [hep-lat/9802011].
- (42) M. Luscher, Phys. Lett. B 593, 296 (2004) [hep-th/0404034].
- (43) L. Giusti and M. Luscher, JHEP 0903, 013 (2009) [arXiv:0812.3638 [hep-lat]].
- (44) B. Berg, Phys. Lett. B 104, 475 (1981); Y. Iwasaki and T. Yoshie, Phys. Lett. B 131, 159 (1983); S. Itoh, Y. Iwasaki and T. Yoshie, Phys. Lett. B 147, 141 (1984); M. Teper, Phys. Lett. B 162, 357 (1985); E. -M. Ilgenfritz et al., Nucl. Phys. B 268, 693 (1986).
- (45) M. Luscher, Commun. Math. Phys. 293, 899 (2010) [arXiv:0907.5491 [hep-lat]].
- (46) M. Luscher, JHEP 1008, 071 (2010) [JHEP 1403, 092 (2014)] [arXiv:1006.4518 [hep-lat]].
- (47) C. Bonati and M. D’Elia, Phys. Rev. D 89, 105005 (2014) [arXiv:1401.2441 [hep-lat]].
- (48) K. Cichy, A. Dromard, E. Garcia-Ramos, K. Ottnad, C. Urbach, M. Wagner, U. Wenger and F. Zimmermann, PoS LATTICE 2014, 075 (2014) [arXiv:1411.1205 [hep-lat]].
- (49) Y. Namekawa, PoS LATTICE 2014, 344 (2015) [arXiv:1501.06295 [hep-lat]].
- (50) C. Alexandrou, A. Athenodorou and K. Jansen, Phys. Rev. D 92, 125014 (2015) [arXiv:1509.04259 [hep-lat]].
- (51) R. Sommer, PoS LATTICE 2013, 015 (2014) [arXiv:1401.3270 [hep-lat]].
- (52) Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi and C. H. Wong, JHEP 1409, 018 (2014) [arXiv:1406.0827 [hep-lat]].
- (53) P. Di Vecchia, K. Fabricius, G. C. Rossi and G. Veneziano, Nucl. Phys. B 192, 392 (1981).
- (54) C. Bonati, M. D’Elia and A. Scapellato, Phys. Rev. D 93, 025028 (2016) [arXiv:1512.01544 [hep-lat]].
- (55) P. Di Vecchia and G. Veneziano, Nucl. Phys. B 171 (1980) 253.
- (56) H. Leutwyler and A. V. Smilga, Phys. Rev. D 46, 5607 (1992).
- (57) Y. Y. Mao et al. [TWQCD Collaboration], Phys. Rev. D 80 (2009) 034502 [arXiv:0903.2146 [hep-lat]].
- (58) F. K. Guo and U. G. MeiÃner, Phys. Lett. B 749 (2015) 278 [arXiv:1506.05487 [hep-ph]].
- (59) G. G. di Cortona, E. Hardy, J. P. Vega and G. Villadoro, JHEP 1601, 034 (2016) [arXiv:1511.02867 [hep-ph]].
- (60) L. Del Debbio, L. Giusti and C. Pica, Phys. Rev. Lett. 94, 032003 (2005) [hep-th/0407052].
- (61) S. Durr, Z. Fodor, C. Hoelbling and T. Kurth, JHEP 0704, 055 (2007) [hep-lat/0612021].
- (62) M. Luscher and F. Palombi, JHEP 1009, 110 (2010) [arXiv:1008.0732 [hep-lat]].
- (63) K. Cichy et al. [ETM Collaboration], JHEP 1509, 020 (2015) [arXiv:1504.07954 [hep-lat]].
- (64) M. Cé, C. Consonni, G. P. Engel and L. Giusti, Phys. Rev. D 92, 074502 (2015) [arXiv:1506.06052 [hep-lat]].
- (65) M. D’Elia, Nucl. Phys. B 661, 139 (2003) [hep-lat/0302007].
- (66) L. Giusti, S. Petrarca and B. Taglienti, Phys. Rev. D 76, 094510 (2007) [arXiv:0705.2352 [hep-th]].
- (67) H. Panagopoulos and E. Vicari, JHEP 1111, 119 (2011) [arXiv:1109.6815 [hep-lat]].
- (68) M. I. Buchoff et al., Phys. Rev. D 89, 054514 (2014) [arXiv:1309.4149 [hep-lat]].
- (69) A. Trunin, F. Burger, E.-M. Ilgenfritz, M. P. Lombardo and M. Müller-Preussker, arXiv:1510.02265 [hep-lat].
- (70) O. Wantz and E. P. S. Shellard, Phys. Rev. D 82 (2010) 123508 [arXiv:0910.1066 [astro-ph.CO]].
- (71) P. A. R. Ade et al. [Planck Collaboration], arXiv:1502.01589 [astro-ph.CO].
- (72) C. G. Callan, Jr., R. F. Dashen and D. J. Gross, Phys. Rev. D 17, 2717 (1978).
- (73) P. de Forcrand, J. E. Hetrick, T. Takaishi and A. J. van der Sijs, Nucl. Phys. Proc. Suppl. 63, 679 (1998) [hep-lat/9709104].
- (74) M. Luscher and S. Schaefer, JHEP 1107, 036 (2011) [arXiv:1105.4749 [hep-lat]].
- (75) A. Laio, G. Martinelli and F. Sanfilippo, arXiv:1508.07270 [hep-lat].
- (76) M. G. Endres, R. C. Brower, W. Detmold, K. Orginos and A. V. Pochinsky, Phys. Rev. D 92, 114516 (2015) [arXiv:1510.04675 [hep-lat]].