Topological susceptibility at high temperature on the lattice

QCD topological susceptibility at high temperature, $\chi_t(T)$, provides an important input for the estimate of the axion abundance in the present Universe. While the model independent determination of $\chi_t(T)$ should be possible from the first principles using lattice QCD, existing methods fail at high temperature, since not only the probability that non-trivial topological sectors appear in the configuration generation process but also the local topological fluctuations get strongly suppressed. We propose a novel method to calculate the temperature dependence of topological susceptibility at high temperature. A feasibility test is performed on a small lattice in the quenched approximation, and the results are compared with the prediction of the dilute instanton gas approximation. It is found that the method works well especially at very high temperature and the result is consistent with the instanton calculus down to $T\sim 2\, T_c$ within the statistical uncertainty.


I. INTRODUCTION
The standard model (SM) is not invariant under P nor CP transformation. Strangely, one of the renormalizable, CP violating terms, θ G µνGµν , exists only with an undetectably small coefficient (θ < ∼ 10 −10 ) or even is missing in the SM, where G µν is the gluon field strength tensor andG µν = ǫ µνρσ G ρσ /2. This unnatural situation is referred to as strong CP problem. Besides the solution with the vanishing up quark mass, the Peccei-Quinn (PQ) mechanism has been known to provide an elegant explanation [1][2][3][4][5][6][7][8]. In the PQ mechanism, a complex scalar field with a U(1) symmetry is introduced. After the spontaneous symmetry breaking of the U(1) symmetry at very high temperature, the radial component of the PQ field acquires the vacuum expectation value, f a , and the angular component, called axion a(x), emerges as Nambu-Goldstone (NG) boson 1 . One can rotate away the leading interaction of the axion to quarks by the U(1) chiral transformation of quark fields, and then the coefficient of the GG term in the SM, θ, is replaced by θ ′ = θ + a(x)/f a . Due to the periodicity of the moduli of the axion field, the effective potential of the axion field takes the form like χ t cos θ ′ , where χ t is the QCD topological susceptibility. This form leads to two important consequences. One is that the CP conserving vacuum (i.e. a /f a = − θ) is automatically chosen at low enough temperature, thus the strong CP problem is gone (PQ mechanism). The other is that the temperature dependent axion mass is given by the QCD topological susceptibility (and f a ) as m 2 a (T ) = χ t (T )/f 2 a . The PQ mechanism is attractive because it also provides a candidate for the dark matter of the Universe through the misalignment mechanism for the axion generation [9][10][11]. The axion abundance of the present Universe is determined by two ingredients: the axion mass as a function of T , m a (T ), and the misalignment at T = T * , θ ′ | T ≥T * , where the axion starts coherent oscillation.
In the estimate of χ t at finite temperature, the instanton picture is widely adopted [12] and predicts T * ∼ 6 T c ∼ O(1) GeV [13][14][15], where T c ∼ 150 MeV is the (pseudo-)critical temperature for chiral symmetry breaking in QCD. However, the instanton calculus is based on the perturbation theory, and hence the reliability is not very clear around T ∼ T * < ∼ 1 GeV. Furthermore, the possibility that, in two flavor QCD, χ t behaves like a step function at T = T c when the quarks are sufficiently light is argued based on reasonable assumptions [16] (see also a clarification in Ref. [17]). In such a case, a significant enhancement of the axion abundance is predicted, and even excludes the standard axion scenario if θ ′ (T * ) = O(1) [18].
Secondly, one of the full QCD calculations in Ref. [21] finds X ∼ 3, which disagrees with X ∼ 8 in the instanton calculus, while X ∼ 8 is reported in Ref. [22]. Thirdly, Q 2 | T = χ t (T ) V 4 rapidly decreases with T , where V 4 represents the four dimensional volume, and importantly existing lattice methods fail when χ t (T )V 4 ≪ 1 2 . Since χ t (T )V 4 ≪ 1 is realized above a certain temperature before reaching T * , the axion abundance becomes uncertain. Thus, methods overcoming this difficulty are desired.
To be specific, in Ref. [18], where χ t (T ) is calculated on 16 3 × 4 lattices in the quenched approximation with one of the standard methods counting the fermionic zero modes, χ t (T )V 4 is estimated to be 0.35, 0.09, 0.03 at T = 1.34 T c , 1.5 T c , 1.75 T c , respectively, and no reliable estimate is given above 2 T c .
In the present paper, we propose a novel method to determine the temperature dependence of χ t (T ) at high temperature. The method involves estimating the difference of the gauge action and the chiral condensate between two different topological sectors. In order to see how well the proposed method works, we perform a test in pure Yang-Mills theory on a small lattice. The results are found to be reasonably consistent with the instanton calculus above T ∼ 2 T c .
The paper is organized as follows. In sec. II, some ingredients of the instanton calculus relevant to subsequent sections are briefly reviewed. After the method is described in sec. III, the numerical test is presented in sec. IV. Summary and outlook are given in sec. V. To supplement the main text, three sections are put in appendix, which include the explicit form of lattice Dirac operators, the discussion of the χ t V 4 ≫ 1 case in this framework, and the analysis of the hopping parameter expansion.

II. INSTANTON CALCULATION
For later use, the calculation of the topological susceptibility in the dilute instanton-gas approximation (DIGA) [14] is briefly reviewed. In general, the gauge action is bounded by where the equality holds when G µν = ±G µν . This self-duality relation is realized in the BPST instanton solution [24]. In instanton calculus, the BPST instanton is taken as the classical background, and the effects of quantum or thermal fluctuations around it is incorporated perturbatively.
Using the partition functions with the topological charge Q = 0 and 1 3 , the instanton density n(ρ, T ) for the one-instanton sector is defined by where z and ρ are the position and the size of the instanton, respectively, and n(ρ, T ) is factorized into the gauge (n G ), the fermionic part (n F ) and the finite temperature effect (n T ).
After the explicit calculation of Z Q=1 [12,25], the gauge contribution to the instanton density is found to be where N is the number of colors and the renormalization scale µ is introduced. Other undefined constants are shortly given. The constant N/6 in eq. (5) does not appear in the Pauli-Villars regularization and appears in the MS scheme [26].
In the presence of N f flavors of quarks with the mass m f [12, [27][28][29][30], the fermionic contribution is given, using the Padé approximation [29,30], by The constants involved in the above equations are a 1 = −13.4138, a 2 = 2.64587, b 1 = 25 592955 21609 a 2 + 255 49 In the finite temperature QCD, the explicit form of the instanton on S 1 × R 3 is known as the HS caloron [31]. While the ρ integral diverges at zero temperature, it becomes finite at finite temperature [14] since the Debye screening exponentially suppresses the large size instanton. This effect is embedded in n T , which is known to be where λ = πρT , c 1 = 0.01289764, and c 2 = 0.15858. Collecting the above expressions, the DIGA predicts the topological susceptibility χ t (T ) at finite temperature to be Later, d ln χ t (T )V 4 /d ln T in the DIGA is numerically estimated to compare with the lattice result, where the running coupling is calculated with the four loop β function. Focusing on the temperature dependence in the high temperature limit where the DIGA is reliable, it follows from eq. (13) that where

III. METHOD
We begin with clarifying our conventions and notations. The QCD partition function on the lattice with a finite θ and in the specific topological sector can be written as respectively, where ∈ Q denotes that the integral is restricted to the configurations with a topological charge Q, andE(θ) (−π < θ ≤ π) is the internal energy density. E(θ) is periodic in θ with periodicity of 2π and symmetric under the change of the sign of θ, i.e E(θ) = E(−θ) = E(2π + θ). The lattice gauge action S g is given by where β = 6/g 2 represents the lattice gauge coupling. The action density P is given by where W P and W R denote the 1 × 1 plaquette and 1 × 2 rectangle averaged over fourdimensional lattice sites, respectively. c 0 and c 1 satisfying c 0 = 1 − 8c 1 are the improvement coefficients for the lattice gauge action. The total number of lattice sites is N site = N 3 S ×N T = V 4 (β)/a 4 (β), where N S and N T represent the number of lattice sites in the spatial and time directions, respectively. For fixed N S and N T , the physical, four dimensional volume V 4 (β) and the lattice spacing a(β) depend only on β. The temperature of the system in the physical unit is given by and hence it can be changed by adjusting either the temporal size N T or the lattice bare coupling g 2 = β/6 4 . S q is the lattice quark action, which we do not specify here. In the following, we consider N f flavors of quarks with a degenerate mass m q (m q = a(β)m q in the lattice unit) for simplicity. The extension to non-degenerate cases is straightforward. Assuming θ = 0, the expectation value of an operator O is expressed as where we have defined Thus, the topological susceptibility times four dimensional volume is written as The simplest method to calculate χ t is to generate an ensemble on the lattice and look at the distribution of Q. As is seen, for example, in Fig. 1 of Ref. [18], an update algorithm employed there only generate configurations with Q = 0 or ±1 at some high temperature 5 . Since those with Q = 0 dominates the other, Z 0 ≫ Z ±1 should hold, and it follows from eq. (22) 4 Thus, we take the mass independent scale setting prescription, where the lattice spacing a does not depend on the quark mass. 5 Note that, even in such a case, the resulting value of χ t turns out to be consistent with more extensive lattice simulations such as Refs. [19,20].
So far, the partition function, Z Q , has been written as a function of β andm q , but an arbitrary pair of arguments can be chosen as long as they fix the QCD coupling and the quark masses. In the following, we consider (T, m q ) and (w = χ t V 4 , m q ) as a pair of arguments, and fix m q to the physical quark mass as function of T or w. In this case, Z Q can be viewed as the function of T or w. Furthermore, the numbers of lattice sites in the spatial and the temporal directions, i.e. N T and N S , are also fixed.
Consider the derivative of the ratio of the partition functions of different topological sectors with respect to temperature with m q and N site fixed, d ln( . Using the chain rule, we rewrite it as d ln Then, the T dependence of w is expressed as In the following, the symbol N site is omitted for simplification. How to estimate each of two factors in the r.h.s. is described below. The first factor, d ln(Z Q 2 /Z Q 1 )/d ln T , can be calculated using lattice numerical simulations as follows, where the temperature, defined in eq. (19), is varied by changing β while N T is fixed. The differentiation of Z Q with respect to T is then given by The β derivative term in eq. (26) is found to be where we have used eq. (17), eq. (19) and the β function for the QCD coupling In perturbation theory, the first two coefficients of β g are given by For our purpose, β g has to be numerically determined as the temperature considered here is of O(T c ).
The term including the mass derivative in eq. (26) are estimated as follows. The first factor is found to be which is related to the anomalous dimension of the quark mass. The second factor is calculated to be where the explicit form of the scalar density operator, sq q , requires specifying the quark action, S q . For example, it is given by for the Wilson fermion, and for the overlap fermion. For details, see appendix A. Gathering eqs. (27), (31) and (32), eq. (26) becomes Taking the difference of eq. (35) for Q 2 and Q 1 , we obtain d ln where we have defined for later use. From eq. (36), it turns out that the differences of the gauge action and the chiral condensate between two topological sectors are required to determines the temperature dependence of Z Q 2 /Z Q 1 . Next, we turn to the second factor in eq. (25), d ln(Z Q 2 /Z Q 1 )/d ln w. In the following, the arguments of the partition functions are omitted for the sake of simplicity. When w ≫ 1, existing lattice methods should work well, and our method is not more efficient than those. However, since it is still instructive to analyze the w ≫ 1 case within the new framework, several remarks are described in the appendix B. Hereafter, we focus on the w ≪ 1 case.
Note that, although w ≪ 1, we assume that the spatial volume is still larger than the typical length scale of the system (∼ 1/T ). Now assume that Z Q can be expanded in terms of w as with an unknown coefficient a Q . While n Q for arbitrary Q is not known, previous numerical simulations tell that there is a temperature region where Z ±1 /Z 0 ≈ w/2 holds [eq. (23)], indicating n ±1 − n 0 = 1. With the assumption eq. (38), it follows that Using eqs. (36) and (39) and recalling d ln If the boundary condition for this differential equation is provided, we can determine the absolute value of χ t (T ). It should be noted that the l.h.s. of eq. (40) is independent of the choice of Q 1 and Q 2 up to O(w). By equating the r.h.s. of eq. (25) for different pairs of Q, we can numerically determine the ratio (d ln( independently of the size of w. Then, the assumption eq. (38) gives Especially, when Q 2 = Q 4 = 0 and Q 3 = 1, R (Q 1 ,0,1,0) (β) = n Q 1 − n 0 . Thus, measuring R (Q,0,1,0) (β) with various Q enables us to investigate the leading power of Z Q /Z 0 , i.e. n Q −n 0 . On the other hand, when w ≫ 1, the behavior of In this case, calculating R (Q 1 ,Q 2 ,Q 3 ,Q 4 ) (β) may serve to check whether w ≫ 1 indeed holds. 6 See eq. (B1) in the appendix B for more details.
Here let us comment on our method. If one could calculate the right hand side of eq. (36) over a wide range of T , Z Q 2 /Z Q 1 can be obtained by the numerical integration with an suitable input. By repeating this procedure for arbitrary pairs of (Q 1 , Q 2 ) and substituting Z Q 2 /Z Q 1 thus obtained into eq. (22), one can determine χ t (T ) over a wide range of T without any assumptions, in principle. If that is possible, the most of above arguments are unnecessary. However, as we will show soon, it turns out that the numerical accuracy is rather limited and the above naive procedure does not work well.
In this work, we instead focus on d ln χ t (T )/d ln T in the temperature region, where χ t (T )V 4 ≈ 2 Z ±1 /Z 0 is valid, because this quantity still provides useful information. For example, the leading powers of w in Z Q 2 /Z 0 , i.e. n Q 2 − n 0 , extracted through eq. (42) for various Q 2 (Q 1 is fixed to zero for simplicity) can be used to identify the θ dependence of the energy density 7 . Furthermore, once an integer value of n Q 2 − n 0 was determined, d ln(Z Q 2 /Z 0 )/d ln T provides an independent determination of d ln(Z ±1 /Z 0 )/d ln T through eq. (40) with n Q 1 = n 0 as we will explicitly show in the next section.

A. high temperature limit
It is instructive to see the high temperature limit of eq. (40). In this limit, the gauge action in each topological sector is expected to realize the BPST instanton solution, at least in the continuum theory, i.e. S g β,mq /β has a finite value in the high T limit, Using the perturbative expression for β g and keeping only the leading order contribution, Collecting the above yields where the O(w) contribution is omitted.
calculus for N f = 0 should also be reproduced in the heavy quark limit, in which the heavy quarks will be decoupled from the theory and hence the β-function is reduced to the one for N f = 0. By imposing that the heavy quark limit of Eq. (46) yields χ t ∼ T −7 , is obtained. The vicinity of the heavy quark limit can be analyzed by applying the hopping parameter expansion, which is described in the appendix C. When N f = 3, the instanton calculus predicts χ t ∼ T −8 , which indicates This coincides with the contributions from the fermion zero modes, −(|Q 2 | − |Q 1 |), when n Q = |Q|.

IV. TEST IN THE QUENCHED APPROXIMATION
A. lattice setup In order to see how well the method described in the previous section works, we perform a test in the quenched approximation. The configurations are generated using the renormalization group improved Iwasaki gauge action, i.e. c 1 = −0.331. The lattice volume is fixed to 16 3 × 4 in this feasibility test except one simulation, in which the calculation is repeated on 24 3 × 4 lattice to see the volume dependence. However, the number of configurations required for a fixed statistical error grows as N site , and our limited computational resources did not allow us to investigate the size dependence in detail.
We use the index theorem in defining the topological charge, Q = Index[D ov ], where D ov is the overlap Dirac operator shown in (A5). Since the configurations in a fixed topological sector is needed, we insert the topology fixing (TF) term, into the path integral [32] to fix Q during the update process. The explicit form of the Hermitian Wilson Dirac operator, H W , is found in eq. (A3). Due to this term, the appearance of the eigenvalues of H W smaller than µ, |λ H W | < ∼ µ, is suppressed, and so is the topology change. In this work, µ = 0.2 and M 0 = 1.6 are used. The standard hybrid Monte Carlo (HMC) method is applied in the configuration generation. The step size in the molecular dynamics procedure is tuned to realize the acceptance ratio of 75% to 90 %. In the preparation step, we first generate configurations at around T c without the TF term to sample the configurations with various Q values. Then, the TF term is turned on, and β is changed to a desired value. The topological charge of configurations thus generated is monitored by calculating the index of the overlap Dirac operator [see eq. (A5)] with the same value of M 0 as that in the TF term, and we checked that no transition to a different Q sector occurs within the configurations used in the analysis except in the Q = −2 sector on 24 3 × 4 lattice, where Q = −2 is changed to −1 after 1,310 trajectories.
In the following plots, we present the statistical error only, which is estimated by the standard single elimination jackknife method with the bin size of 50 trajectories. Increasing the bin size by a factor two only changes the size of uncertainty by a few %.
The theory with the TF term (49) is not rigorously equivalent to the quenched QCD, because the TF term (49) would break Z 3 symmetry. Thus, strictly speaking, the action with the TF term may not allow us to study the phase transition of the quenched QCD. Thus, our study focuses on the temperature region like T ≥ 2 T c .
It is also important to note that the presence of the TF term, in general, changes the correspondence between the simulation parameter β and temperature T . By using the fact that the spectrum of the Dirac operator is sensitive to the temperature, we see how much the correspondence between the simulation parameter β and temperature T is shifted in the presence of the TF term. The distribution of the smallest eigenvalues of the Hermitian Wilson (H W ) and overlap (H ov ) Dirac operators are shown in Fig. 1 as examples, where β = 2.450, 2.802 and 10 correspond to T ∼ 1.3 T c , 2.25 T c and 8 × 10 3 T c , respectively.
In Fig. 1 (left), the suppression of the appearance of small eigenvalues is clear at T ∼ 1.3 T c (left) while no significant difference is observed at T > 2 T c . As for the Hermitian overlap Dirac operator [ Fig. 1 (right)], while the effect of the TF term is again clear at low temperature (top) especially in the near-zero mode region, the distributions reasonably agree at high temperatures (middle and bottom). The temperature region we are interested in is T > ∼ 2 T c and in such a region the Dirac spectra with and without the TF term turn out to agree at the same β values. This observation allows us to employ the relationship between the simulation parameter β and temperature T obtained in simulations with the same gauge action but without the TF term. Although it would be possible to numerically take the µ=0 limit, we do not pursue the limit in this exploratory study.
The configurations are generated at 12 values of β ranging from T c to 10 4 T c and in four different topological sectors, Q = 0, 1, −1, −2. The configurations are stored every 10 and 5 trajectories for 16 3 × 4 and 24 3 × 4 lattices, respectively. The simulation parameters and statistics are tabulated in Tab. I. The values of T /T c in the table are obtained by using the formula provided in Ref. [34], where the lattice spacings are determined in a wide range of β using the same gauge action as ours but without the TF term.

B. numerical results
In quenched QCD, the T dependence of χ t is determined by The results of ∆S (Q,0) g (β) with Q = ±1 and −2 are shown in Fig. 2, where it is seen that the data for Q = 1 and −1 agree well within the statistical error as expected. Thus, the averaged value over Q = 1 and −1 is used in the following analysis.
The horizontal dotted lines represent the action difference in the BPST instanton solutions, or in the high temperature limit, for |Q| = 1 and 2 from top to bottom. The lattice data for Q = ±1 are on top of the corresponding BPST line down to β ∼ 2.5 (or T ∼ 1.45 T c ) and suddenly decrease at β ∼ 2.4. The similar behavior is observed for Q = −2 but the deviation from the corresponding BPST line starts at slightly larger β, β ∼ 3. The jump observed at β ∼ 2.4 may be associated with the phase transition. Studying the phase transition itself within this framework is interesting, but we focus on the high temperature region in this paper.
The large volume results are also shown in Fig. 2 (filled symbols). It is confirmed that the Q = ±1 result are consistent with that from the smaller lattice. We omit the Q = −2 result on the larger lattice from the figure because of a large uncertainty.
In order to estimate d ln χ t /d ln T using ∆ (Q,0) S g (β) with Q = ±1 or −2, we need to know n Q − n 0 . n ±1 − n 0 = 1 is empirically known 8 . We can estimate n 2 − n 0 by looking at R (2,0,1,0) (β) [see eq. (38)]. Figure 3 shows that R (2,0,1,0) (β) is consistent with two over the whole range of β we have studied, but the large statistical errors do not allow the precise determination except for the region of β ≥ 10. It is seen that, when the mean value is relatively large, the error is also large. Thus, we assume n −2 − n 0 = 2 in the following analysis.
The QCD beta function, β g , down to a low energy scale (∼ T c ) is necessary in estimating eq. (50). We use the result of Ref. [34], in which the lattice spacing is expressed as a function of the lattice gauge coupling, β, as where σ denotes the string tension and Using this expression, β g is numerically determined through At the same time, the relationship between T /T c and β is found to be where T c = T (β c ) and β c = 2.288 [34]. β g and T (β)/T c are shown as a function of the lattice gauge coupling β in Fig. 4. In the plot, we also show β 2 β g , which approaches to β 2 β g → 792/(4π) 2 ∼ 5 in the large β limit. Substituting the above results into eq. (50), d ln χ t /d ln T is calculated as shown in Fig. 5, where the two solid curves represent the prediction of the DIGA (13) with µ = πT /2 and 2πT , respectively although they can not be distinguished at this axis scale.
The results with |Q| = 1 and 2 are consistent with each other, which is expected from the observation in Fig. 3. These results are also consistent with the high temperature limit and the DIGA down to T /T c ∼ 1.5. Note that the results using the Q = −2 sector has the uncertainty smaller than those using Q = ±1 by a factor n −2 − n 0 = 2, which indicates that once the n Q − n 0 has been fixed one can obtain very accurate result by performing a simulation at large Q.
One of the concerns in this approach is the finite volume effect since the physical volume becomes extremely small at large β. Figure 5 shows that the lattice results well reproduce the high temperature limit at high temperature. From this observation, it is unlikely that the finite size effect significantly affects the lattice results, and it is natural to think that N T ≪ N S is the necessary condition for the finite volume effects to be under control. Indeed, the aspect ratio of our lattices is N S /N T = 4, and hence the above condition seems to be satisfied. Nevertheless, calculations with different lattice sizes are clearly useful to explicitly check the finite size effect and whether w ≪ 1 holds or not. However, since the uncertainty of the action value grows as √ N site , we need the statistics proportional to N site to keep the size of the uncertainty constant.
From the phenomenological point of view, χ t (T ) for T c < ∼ T < ∼ 10 T c is important. In this range of T , the statistical uncertainty is relatively large (typically ±4 for O(10, 000) trajectories), which makes the axion abundance ambiguous. It is thus important to accumulate a large number of statistics. On the other hand, if χ t (T ) behaves like a step function, our method should be able to detect such a behavior.

V. SUMMARY AND FUTURE PROSPECTS
The QCD topological susceptibility, χ t , at high temperature provides an important input for the estimate of the axion abundance in the present universe. Existing methods to calculate χ t on the lattice in the literature fail when χ t (T )V 4 ≪ 1. We proposed a novel lattice method to calculate the temperature dependence of the susceptibility, which is expected to work well especially in high temperature region where χ t (T )V 4 ≪ 1. To see how it works, we performed quenched simulations on the 16 3 × 4 lattice, and found that the results of d ln χ t /d ln T well agree with the DIGA prediction above 1.5 T c . The simulation on a slightly larger lattice confirms that there is no unexpected large finite volume effect, although keeping the statistical error constant requires statistics proportional to N site . Thus, that error may be the main source of uncertainty in future serious works.
To predict the axion abundance, we still have to include dynamical quarks with the physical masses. In order for the method to work, the difference of the chiral condensate between two topological sectors has to be precisely determined, for which the dynamical overlap fermion seems to be preferred. Then, accumulating a large number of configurations requires large amount of resources. But, if χ t (T ) behaves like a step function, a large number of statistics may not be necessary to detect such a behavior.
A possible way out is to generate configurations in large Q sectors, with which one can achieve an uncertainty smaller than that with |Q| = 1 by a factor of n Q − n 0 . Note that this requires the signal on R (Q,0,1,0) [eq. (43)] at w ≪ 1 to be sufficiently precise to unambiguously identify the integer n Q − n 0 . Knowing n Q − n 0 for various Q is also useful to put the constraints on the θ dependence of E(θ) in eq. (15), whose general form would be given by with χ t = n c n /n 2 .
The relationship between expectation values in the θ vacuum and a fixed topological sector is derived in Refs. [36][37][38], where the results are written in the form of 1/w expansion (w = χ t V 4 ). From eq. (2.38) of Ref. [38], it is read that d ln where c 2n is defined as the coefficient of θ 2n in the expansion, especially c 2 (T ) = χ t (T ). Note that c i (T ) is dynamical quantity and depends on T . No assumption is made in deriving eq. (B1) except for w ≫ 1 and Q 2 ≪ w, but in order for the expansion to be sensible, c 2n /χ t ∼ O(1) is required. It is interesting to see that eq. (B1) is proportional to Q 2 2 − Q 2 1 while the same quantity is to |Q 2 | − |Q 1 | when w ≪ 1 [see, for example, eq. (39)].
Using eq. (B1) and estimating d ln(Z Q 2 /Z Q 1 )/d ln T , one can estimate the T dependence of χ t through eq. (25). Although d ln(Z Q 2 /Z Q 1 )/d ln T can be directly estimated on the lattice as explained in the main text, we can further use the relationship between expectation values in the θ vacuum and a fixed topological sector [36][37][38] to proceed. From eq. (3.20) in Ref. [38], it follows that where x k are defined as follows, and q 4 = Q 4 θ=0 . The 1/w expansion in eq. (B3) is sensible if x i ∼ O(1). Applying the above expansions to the quenched case yields where Appendix C: hopping parameter expansion We consider the temperature dependence of χ t in the presence of N f flavors of heavy quarks (and no light quarks). In the following, we assume that the degenerate heavy mass is larger than the temperature, m q ≫ T . By introducing the hopping parameter the Wilson Dirac operator, eq. (A2), can be rewritten as Thus, when κ q ≪ 1, we can expand D W in terms of κ q . If we take the Wilson fermion as the heavy quark action, the partition function and the expectation value of operator O in a fixed topology sector can be written as respectively. Applying the hopping parameter expansion (HPE) to the heavy quarks, the determinant can be expanded as where L denotes the real part of the Polyakov loop and we have used N T = 4. Then, the expectation value, eq. (C4), is given by If O consists of quark fields like O =q x q x , the HPE is further applied. Using β,mq can be written as In summary, It turns out that, in the heavy quark region, the Polyakov loop plays an important role. In Fig. 6, the difference of L in different Q sectors is shown. We define the following quantities, where we omit the term including the anomalous dimension of quark mass. Figure 7 shows the κ q dependence of the above two quantities with N f = 2. It turns out that the contribution of ∆ (Q) Y K is much smaller than that of ∆ (Q) S K except for β = 100. Thus, omitting the term of the anomalous dimension does affect the final result by much. ∂ ln χ t /∂ ln T is plotted in Fig. 8, which shows no clear κ q dependence up to κ q = 0.1 except at β = 100. At β = 100, the critical kappa is ∼ 0.125, thus the result at β = 100 explores the relatively light quark mass region and may indicate a tendency that d ln χ t /d ln T decreases towards the chiral limit, although the convergence of the HPE has to be checked.  [15] O. Wantz and E. P. S. Shellard, "Axion Cosmology Revisited," Phys. Rev. D 82, 123508