Topology in QCD and the axion abundance

The temperature dependence of the topological susceptibility in QCD, chi_t, essentially determines the abundance of the QCD axion in the Universe, and is commonly estimated, based on the instanton picture, to be a certain negative power of temperature. While lattice QCD should be able to check this behavior in principle, the temperature range where lattice QCD works is rather limited in practice, because the topological charge is apt to freezes at high temperatures. In this work, two exploratory studies are presented. In the first part, we try to specify the temperature range in the quenched approximation. Since our purpose here is to estimate the range expected in unquenched QCD through quenched simulations, hybrid Monte Carlo (HMC) algorithm is employed instead of heatbath algorithm. We obtain an indication that unquenched calculations of chi_t encounter the serious problem of autocorrelation already at T~2Tc or even below with the plain HMC. In the second part, we revisit the axion abundance. The absolute value and the temperature dependence of chi_t in real QCD can be significantly different from that in the quenched approximation, and is not well established above the critical temperature. Motivated by this fact and precedent arguments which disagree with the conventional instanton picture, we estimate the axion abundance in an extreme case where chi_t decreases much faster than the conventional power-like behavior. We find a significant enhancement of the axion abundance in such a case.


Introduction
It is widely believed that the instanton calculus in QCD [1] makes sense at high temperatures.
The asymptotic freedom ensures the perturbative expansion sensible and, most importantly, the infrared divergences from the large instanton contributions are cut-off by the Debye length, providing finite results after the integration over the instanton size [2].
In the semi-classical instanton picture, physics becomes θ-parameter dependent by the instanton contributions to the path integral [3]. The instanton calculus indicates that such a dependence is proportional to the product of the quark masses, m N f q and Λ b QCD , where b is the beta-function coefficient, b = 11N c /3 − 2N f /3. For example, the topological susceptibility, In the QCD axion model to solve the strong CP problem, the θ angle is promoted to the axion field a(x)/f a , where f a is the axion decay constant [4,5,6,7,8,9,10,11]. The topological susceptibility is directly related to the mass of the axion as χ t = m 2 a f 2 a . Therefore, the temperature dependence of χ t discussed above represents that of the axion mass, which is important for the calculation of the axion abundance in the Universe. In the misalignment mechanism for the axion generation in the early Universe, the axion number density is proportional to the axion mass at the temperature at which the axion field starts coherent oscillations [12,13,14]. The instanton based estimation of the temperature dependence is commonly used in the literature, and predicts that the axion can naturally be dark matter of the Universe when m a ∼ 10 −5 eV, whereas the astrophysical bound on m a is m a 10 −2 eV.
(See, e.g., [15].) The allowed region, 10 −5 eV m a 10 −2 eV, is called the "axion window." There have been arguments which imply the disappearance of the instanton effects at high temperatures. In Ref. [16], it has been argued in two-flavor QCD that the difference between the two-point functions of iso-singlet and iso-vector scalar operators vanishes in the chiral limit as fast as O(m 2 q ) when the temperature is higher than the critical temperature. By the Ward-Takahashi identities, this immediately means that χ t should be at most an O(m 4 q ) quantity, that is different from O(m 2 q ) from the instanton calculus. Moreover, in a recent paper [17], it is claimed that χ t is O(m N q ) with an arbitrary N and thus it is vanishing even with finite quark masses. Although these are somewhat surprising results, it is certainly possible that the instanton picture fails to describe the full quantum mechanical vacuum [18].
If that is the case, the estimation of the axion abundance is significantly affected.
The lattice simulations can in principle discriminate whether instantons make sense or not. However, the current situation is not conclusive. In Refs. [19,20], the temperature dependence of various susceptibilities are measured, which seem to be consistent with the semi-classical instanton picture and but at the same time not to exclude other possibilities.
Meanwhile, the analysis in Ref. [21] suggests the effective restoration of U A (1) symmetry right above the critical temperature, indicating the disappearance of the instanton effects.
The presence or absence of the U A (1) symmetry at the critical temperature T c can also be inferred through the nature of the chiral phase transition of massless two flavor QCD as attempted in, e.g., Ref. [22] (For a related work analyzing the system with the renormalization group flow, see Ref. [23]). See Ref. [24] for a new proposal to extract the instanton effects in lattice QCD, and Ref. [25] for the approach from holography.
One unavoidable problem on the lattice is that the range of the temperature in which one can reliably study the temperature dependence of χ t is rather limited, because the net susceptibility χ t V , where V is the volume, rapidly decreases with temperature and hence the topology tends to freeze at high temperature. This problem is fatal especially in dynamical lattice QCD simulations since one can not realize arbitrary large volumes while keeping both the light quark masses and lattice artifacts reasonably small. Then, it is natural to ask to what temperature χ t can be reliably calculated in dynamical QCD with a typical setup. From this viewpoint, determining the precise value of quenched χ t (T ) with large volumes and huge statistics would not provide useful information.
For the reliable calculation of χ t , it is preferable to use lattice chiral fermion for dynamical quarks. But such dynamical simulations at high temperatures are too costly to start with. As the second best, we choose to perform quenched simulations with hybrid Monte Carlo (HMC) algorithm and Iwasaki gauge action. The reasons for these choices are as follows. To answer the above question, the simulation environment needs to be as close to that in dynamical one as possible. Long auto-correlation of the topological charge is one of the serious bottlenecks in dynamical simulations, and the heatbath algorithm, which is available only in the quenched approximation, is faster but has totally different property from HMC in this regard. Iwasaki gauge action is known to result in relatively long autocorrelation time for the topological charge [26,27]. Thus, we expect that more information about dynamical simulations will be gained from quenched one by taking HMC and Iwasaki gauge action rather than taking heatbath and the standard plaquette gauge action.
Another issue to be addressed is the definition of topological charge Q. There is a subtlety in measuring Q if one adopts the field theoretical (or bosonic) definition. This method requires a suitable number of coolings of gauge configurations. If one applies it too much or too less, one would miss the right value of Q. Even if one ceases the cooling adequately, Q thus obtained is not an integer, and in the worst case it may be right in the middle of two integers, which can cause misidentification of Q. Misidentifications are potentially dangerous, especially at high temperatures (i.e. when Q 2 = χ t (T ) V is tiny), because it can significantly affect χ t (T ) even if it occurs only rarely. Importantly, it is not possible to know the right value of Q without comparing that obtained with the fermionic definition based on Atiyah-Singer index theorem. From this viewpoint, the use of the fermionic definition has an advantage, although it is much more demanding than the bosonic one. Studying the autocorrelation of Q usually requires high statistics. Thus, this choice of the definition seriously limits the lattice size to small.
In this paper, we first explore to which temperature quenched simulations using HMC and Iwasaki gauge action are able to obtain reliable results for χ t , where the topological charge Q in each configuration is determined by the number of zero modes of overlap Dirac operator.
Because of this choice, the lattice volumes are limited to relatively small. Our quenched calculations show that χ t is undetermined above T = 2 T c for N t = 4 and T = 1.5 T c for N t = 6, where N t is the number of lattice sites in the time direction.
This observation suggests two possibilities. One is that the result with N t = 6 is right and χ t suddenly decreases at T 1.5 T c . In the language of the semi-classical instanton picture, this behavior can happen if there is a sharp short-distance cut-off in the instanton size parameter ρ so that small instantons do not contribute to the path integral. (See [28,29] for the instanton based models and lattice results suggesting this picture.) If χ t exhibits such a sharp fall off, the estimation of the abundance of the QCD axion in the Universe may be significantly affected. However, this possibility appears to be unlikely because our value of χ t with N t = 4 is reasonably consistent with the previous calculations to T ∼ 1.75 T c . Another possibility is that, in the standard HMC with a fixed acceptance ratio, the suppression of χ t due to the long autocorrelation, which increases with a lattice volume, overcomes the enhancement by enlarging the volume. Through the analysis of the autocorrelation, we confirm the latter to be the case.
Outcome of the first part is that with HMC the reliable calculation of χ t is difficult for T 2 T c even in the quenched approximation. Since quenched simulations are always easier than dynamical one, the above outcome or even worse should hold for dynamical QCD.
The second part of the paper goes as follows. We first emphasize that even in the quenched theory it is hard to justify instanton calculus from the first principles because instanton description consists of the perturbative expansion in terms of the strong coupling constant α s (T ) and its prediction is only reliable above T = a few GeV. The same is true or the situation is even worse for dynamical QCD. Furthermore, there exist convincing arguments that the presence of light dynamical quarks could drastically change χ t above T c as stated above. Therefore, not much has been known about the unquenched theory, and there are many degrees of freedom in the choice of T dependence of χ t above T c .
We revisit the axion abundance with some extreme but yet allowed temperature dependence of χ t . We find that the axion abundance is significantly enhanced compared to the conventional estimates with non power-like behavior of χ t and, for some cases, the axion window is closed depending on how quickly χ t disappears.

Topological susceptibility on the lattice
In the continuum, the topological susceptibility, χ t , is defined as This quantity is related to the axion mass, m a , through where f a is the axion decay constant. Through the Atiyah-Singer index theorem, one can also find where the topological charge, Q = n + − n − , is the difference between the numbers of the left and the right-handed zero modes in the eigenvectors of the Dirac operator.
On the lattice, the definition of gauge field strength tensor, F µν , is not unique, and so is χ t if one simply follows Eq. (1). With the definitions of F µν commonly used in the literature, takes non-integer values in general. Thus, a technique called "cooling" is usually applied to guess the right integer for Q.
Similarly to the continuum, an alternative is to count the zero eigenvalues of the lattice Dirac operator. This requires lattice Dirac operators to satisfy the Ginsparg-Wilson relation [30]. The overlap Dirac operator, D ov , is known as such an operator [31], with which one can express Q as where Here, As stated in the introduction, the bosonic definition has the chance to misidentify Q.
While, below and around the (pseudo-)critical temperature, previous works have revealed that the bosonic and fermionic definitions of Q give consistent results for χ t , it is not well established above 1.5 T c yet. Especially, we have to keep in mind that a tiny amount of misidentification of Q may bring significant effects to the resulting χ t at high temperatures as follows.
Whether the instanton picture are correct, it is certainly true that χ t decreases with the temperature and χ t V as well. Since χ t V represents a net width of the fluctuation of Q, the fraction of configurations with non-zero Q becomes much smaller than that with the vanishing Q when χ t V 1. Suppose that N conf configurations are generated on the lattice volume of N V = N 3 s × N t and that q of them belong to either Q = +1 or Q = −1 while the others to Q = 0. This yields a 4 χ t = Q 2 /N V = (q + δ mis )/(N conf N V ), where δ mis /N conf represents the rate of misidentification. Then, it should be noted that, even if δ mis /N conf is tiny, a 4 χ t may significantly deviate from a real value because the fraction of nonzero Q configurations, q/N conf , is also tiny at high temperatures. To avoid the misidentification, we adopt the fermionic definition in Eq. (5) for the evaluation of χ t .

lattice parameters
To calculate the topological susceptibility at finite temperatures, we perform lattice simulations of SU (3) pure Yang-Mills theory around and above T c . We employ the Hybrid Monte Carlo (HMC) algorithm to study the autocorrelation of topological charge in HMC for the reason described in introduction. To avoid the ordinary finite size effects, the aspect ratio is set to three or four, which is usually considered to be safe. In addition, we have to recall that the size of fluctuation of Q is directly affected by a volume through Q 2 = χ t V , which is a finite size effect specific to and important in the calculation of χ t . The calculation is carried out on three lattice volumes using the Iwasaki gauge action, 16 3 × 4, 18 3 × 6 and 24 3 × 6.
The topological charge Q is calculated at every ten trajectories. The statistical errors are estimated by the standard jack-knife method with the bin size of 500 trajectories.
The correspondence between β and T /T c for the Iwasaki gauge action is read from Ref. [32]. We accumulated 2000 to 40000 trajectories, depending on the simulation parameters, which are summarized in Tab. 1.
Lattice volumes chosen in the present work are somewhat smaller than other pure YM calculations. The reason is as follows. In order to reduce the ambiguity associated with the definition of Q as much as possible, we decided to use the fermionic definition for Q to calculate χ t . The purpose of this work is to find the highest temperature to which the reliable extraction of χ t is feasible in HMC, or in other words, at which temperature the autocorrelation time becomes unreasonably long. Such a study usually requires high statistics.
Furthermore the fermionic definition is much more expensive than the bosonic one. Thus this choice of definition constrains lattice volumes to be small.
Here let us mention the difference of this work from Ref. [33], which studies the T dependence of χ t around the critical temperature, using the fermionic definition for Q. Two major differences are in the range of the temperature explored and rigorous application of the overlap Dirac operator. The former is originating from the different motivation. Let us explain more about the latter. In Ref. [33], the identification of Q is made using an approximate solution of the Ginsparg-Wilson equation, while we stick to the exact one. They differ in important way. To check the approximated solutions, in Ref. [33] the chirality of the zeromodes (ψ † i γ 5 ψ i ) is examined, and found that its absolute value ranges from 0.4 to 0.9. With our exact method, it only takes +1 or -1 within a rounding error. Therefore, in Ref. [33], the possibility of the misidentification cannot be excluded thoroughly. Indeed, they observe that 2 to 9 % of configurations are misidentified by comparing their approximated solutions and the exact one below the critical temperature. Importantly, it is not clear how well the method based on the approximation works at higher temperatures. In our calculation, the size of zeromodes and near-zeromodes of Hermitian overlap operator H ov differ by a factor of O(10 6 ), and we checked for all zeromodes whether the complex conjugate pair is absent.
Thus, no ambiguity exists. While 2 to 9 % of misidentifications will not significantly affect χ t at low temperatures, it does at high temperatures as χ t itself is very small there. Note that both of the two differences are essential in the study of axion. It is naively expected that Q fluctuates more frequently at larger lattices. However, the plots in Fig. 1 shows the opposite tendency, which indicates that the suppression of χ t due to the long autocorrelation in HMC with a fixed acceptance ratio overcomes the enhancement by enlarging the volume for our choice of lattice setup.

autocorrelation of Q
To quantify the autocorrelation, we estimate the integrated autocorrelation time, τ int .   We observe the qualitatively similar behavior for 18 3 × 6 lattices except for the T = T c data, which shows τ int larger than other results obtained at T ∼ T c .

T dependence of χ t
The figure 3 shows the T dependence of χ t . It is seen that, restricting the data points to the one in which Q shows a recognizable fluctuation, the results obtained from such fluctuations are consistent with the existing results, for example, given in Ref. [33,34,35]. It is important to note that we have only observed |Q| ≤ 1 at T = 1.50 T c and 1.75 T c , but the resulting χ t s nevertheless reasonably agree with the previous results. It is also found that χ t with different volumes shows consistency within two standard deviations as long as the data showing a recognizable fluctuation of Q are concerned. From these observations, we realized that with HMC it is difficult to obtain the reliable χ t above T = 2 T c or even 1.5 T c depending on the lattice volume, which immediately indicates the difficulty in dynamical simulations. Even if much faster computers were used, this upper bound will not change significantly. Thus, to estimate the T dependence of χ t at O(10 T c ), we have to make a long extrapolation using those obtained in such a rather limited range of T . To push the limit upward as high as possible, it is crucial to explore the HMC parameters or improve the algorithm.
Inspired by Refs. [36,37], we tried, as an attempt, to enhance the number of configurations with nonzero Q by inserting to the path integral, where N φ is a positive integer. Then, χ t is calculated through where · · · X denotes the average over the configurations generated with the extra reweighting factor X. For µ > , the insertion of X enhances the eigenvalue density in the small eigenvalue region, whereas eigenmodes with eigenvalues λ µ are left untouched. Since, when the topology changes, the smallest eigenvalue of H W passes through zero, the above factor is  expected to increase such opportunities. However, after we performed some trial calculations, we realized that this method does not always work and the fine tuning of µ, and N φ are required. Further investigations to improve the situation is in progress.

Effects of dynamical quarks
Let us discuss what would happen when we include the dynamical quarks. The naive guess would be that χ t in the Yang-Milles theory is multiplied by a factor of m There can be more drastic possibilities. If we accept the claims of the axial U (1) restoration in two-flavor QCD [16,17], the O(m 2 q ) contributions to χ t is forbidden in two-flavor QCD. Therefore, the possibility of just multiplying by m QCD is not consistent. The results of Ref. [17] even forbid contributions with any power of m q for a small m q . An extreme possibility one can consider is with c(m q ) → ∞ as m q → 0, so that χ t cannot be expanded around m q = 0. Note that the results of Ref. [17] is contained as a special case of eq. (9). Since no unquenched result of χ t is available at high temperatures, we take c(m q ) as a free parameter in the following discussion.

Axion abundance
We discuss the impact on the axion abundance in the Universe for the case where the dilute instanton gas approximation fails at high temperatures. One of the source of the axion energy density today is the coherent oscillations of the axion field started in the early Universe. The equation of motion for the axion field, a, is given bÿ where H is the Hubble parameter and m a (T ) is the axion mass at a temperature T . The axion mass is related to χ t as in Eq. (2). This equation of motion leads to an oscillating solution. Although the axion mass is temperature dependent, and thus time dependent, the axion number density divided by T 3 , n a /T 3 , stays constant during the oscillation and its value is given by where T * is the temperature at which the axion field starts to oscillate, and θ is the initial By using this value and setting m a (T * ) 3H(T * ), the energy fraction today is finally given by Ω a 0.2 · θ 2 · m a 10 −5 eV for θ 1, and it is about a factor of two larger for θ ∼ O(1) [38,39]. According to this formula, the axion can naturally be the dark matter of the Universe for m a ∼ 10 −5 eV which corresponds to f a = 6 × The condition is automatically satisfied when m a (T ) ∼ H(T ) and m a (T ) ∝ T n with n ∼ O(1), but not necessarily true for exponential functions.
In the model discussed in the previous section, m a (T ) is given by where the c parameter may be much larger than O(1). The oscillation starts when both sides of Eq. (14) become comparable, and thus Since m a (T ) is a steeply varying function, we expect T * ∼ T 1 , where T 1 is defined by This causes an enhancement of the axion density by a factor of about c(T 1 /T c ) 2 compared to the conventional estimate in addition to the enhancement due to the shift of the oscillation temperature.
A numerical solution of the equation of motion in Eq. (10) shows that the axion number density is approximately given by Here Here we used 3H(T c ) 5 × 10 −20 GeV, and assumed (T 1 /T c ) 2 m a (T c )/3H(T c ). Putting it all together, the axion energy density today is given by Even for c = O(1) and m a = 10 −5 eV, we find that the axion energy density is enhanced by an order of magnitude compared to the conventional estimate. The enhancement is larger for larger c. The axion window is closed when c ∼ 400. For an extremely large value of c, i.e., c 4 × 10 5 · (m a /10 −5 eV), the axion starts to oscillate at T = T c . In that case, Ω a ∼ 2 × 10 5 · θ 2 independent of c.
We have also examined the case of χ t ∝ exp(−2c T /T c ) rather than χ t ∝ exp(−2c T 2 /T 2 c ) in Eq. (9). We find a similar enhancement factor 1.0 c instead of 2.5 c in Eq. (20).
In the scenario where the Peccei-Quinn symmetry breaking takes place after the inflation, one should take an average value of θ, θ 2 π 2 . Also, in addition to the misalignment mechanism considered above, there are contributions to the axion density from the axionic strings and domain walls [41,42,43,44]. Both contributions receive an enhancement due to the delay of time scale of the axion production which leads a milder dilution due to the cosmic expansion. The enhancement factor of Ω a compared to the estimates based on For m a = 10 −5 eV, for example, the enhancement is a factor of a few for c = O(1). The value of Ω conv.

Summary
The temperature dependence of the topological susceptibility χ t essentially determines the abundance of the QCD axion in the Universe, and can be calculated using lattice QCD in model independent way. However, at high temperatures, the lattice calculation suffers from various difficulties such as the numerically demanding definition, possible misidentifications and the long autocorrelation of the topological charge.
As the first step towards the precise calculation, we have performed the lattice calculation Meanwhile, in order to see the importance of understanding the temperature dependence of χ t for the axion abundance in the Universe, we have calculated the axion energy density, Ω a , in an extreme case where χ t drops exponentially above the critical temperature. It is found that the non-adiabatic growth of the axion potential makes the axion abundance significantly larger, possibly closing the axion window. In this particular case, the behavior of χ t around T c becomes important.