Study of axial U(1) anomaly at high temperature with lattice chiral fermions

We investigate the axial U(1) anomaly of two-flavor QCD at temperatures 190--330 MeV. In order to preserve precise chiral symmetry on the lattice, we employ the Mobius domain-wall fermion action as well as overlap fermion action implemented with a stochastic reweighting technique. Compared to our previous studies, we reduce the lattice spacing to 0.07 fm, simulate larger multiple volumes to estimate finite size effect, and take more than four quark mass points, including one below physical point to investigate the chiral limit. We measure the topological susceptibility, axial U(1) susceptibility, and examine the degeneracy of U(1) partners in meson and baryon correlators. All the data above the critical temperature indicate that the axial U(1) violation is consistent with zero within statistical errors. The quark mass dependence suggests disappearance of the U(1) anomaly at a rate comparable to that of the SU(2)_L x SU(2)_R symmetry breaking.

may have an impact on our understanding of the phase diagram of QCD [6]. In the standard effective theory analysis of the chiral phase transition, we expect only the pion and its chiralpartner scalar particle to be the light degrees of freedom that govern the low-energy dynamics of QCD near the phase transition. When the U (1) A symmetry is effectively restored, the η, isospin singlet pseudoscalar, and its chiral partner may play a nontrivial role [7][8][9]. Since the effective potential can have more complicated structure as the degrees of freedom increase, it was argued that the chiral phase transition likely becomes first order when U (1) A symmetry is recovered (we refer the readers to [10] for a different aspect of the first-order scenario from 't Hooft anomaly matching), though other scenarios are theoretically possible [11][12][13][14][15][16].
How much the U (1) A anomaly contributes to the dynamics has a significant importance on cosmology. The topological susceptibility is related to the mass and decay constant of the QCD axion, which is a candidate of the dark matter. Its temperature dependence influences the relic abundance of the axion [17][18][19][20][21][22][23].
For the nontrivial question of how much the U (1) A anomaly remains and affects the chiral phase transition, only lattice QCD can give a quantitative answer. Good control of chiral symmetry on the lattice [24][25][26][27] is necessary in order to precisely discriminate between the lattice artifact and the physical signal that survives in the continuum limit. Our previous studies [28][29][30] demonstrated that the signal of topological susceptibility is sensitive to the violation of chiral symmetry at high temperature, and the reweighting of the Möbius domainwall fermion determinant to that of the overlap fermion is essential when the lattice spacing is coarse a > ∼ 0.1 fm. We also found that the use of the overlap fermion only in the valence sector [31][32][33] makes the situation worse, since the lattice artifact due to the mixed action, which is unphysical, is strongly enhanced.
After removing the lattice artifact due to the violation of the Ginsparg-Wilson relation at high temperature, we observed that chiral limit of the U (1) A susceptibility is consistent with zero [30]. The disappearance of the U (1) A anomaly (at around 1.2 T c ) was also reported by other groups simulating non-chiral fermions [34,35] 1 . In [37], it was found that the U (1) A symmetry shows up at 1.3 T c but not around T c 2 .
1 In [36], they gave a preliminary results showing that the anomaly studied in [34] looks enhanced as the volume size increases. 2 A recent work [38] reported that U (1) A symmetry is still broken at 1.6 T c using highly improved staggered quarks with lattice spacings 0.06-0.12 fm.
In this work, we reinforce the conclusion of [30] by 1) reducing the lattice spacing to an extent where we observe consistency between the overlap and Möbius domain-wall fermions, 2) simulating different lattice volume sizes in a range 1.8 ≤ L ≤ 3.6 fm, and 3) simulating more quark mass points, including one below the physical point, to investigate the chiral limit.
In order to study possible artifact due to topology freezing, we also apply the reweighting from a larger quark mass, where topology tunneling is frequent, down to the mass where the topological susceptibility is consistent with zero. We measure the topological susceptibility, axial U (1) susceptibility, meson correlators, and baryon correlators. Some of the results were already reported in our contributions to the conference proceedings [39][40][41][42][43][44][45].
All the data at temperatures above 190 MeV show that the axial U (1) anomaly is consistent with zero, within statistical errors. Its quark mass dependence indicates that the disappearance of the U (1) anomaly is at a rate comparable to that of the SU (2) L × SU (2) R symmetry. At higher temperature, we have also observed a further enhancement of symmetry [46][47][48][49].
The rest of the paper is organized as follows. We describe our lattice setup and how to implement the chiral fermions in the simulations in Sec. II. In Sec. III, the numerical results for the Dirac spectrum, topological susceptibility, axial U (1) susceptibility, and meson/baryon screening masses are presented. The conclusion is given in Sec. IV.

II. LATTICE SETUP
The setup of our simulations is basically the same as our previous study [30], except for the choice of parameters (larger lattice sizes up to 3.6 fm, and smaller lattice spacing a ∼ 0.074 fm). The (pseudo) critical temperature of the chiral phase transition we estimated in [30] is T c ∼ 175 MeV, which was obtained from the Polyakov loop. Below, we summarize the essential part.
In the hybrid-Monte-Carlo (HMC) simulations 3 , we employ the tree-level improved Symanzik gauge action [52] for the link variables and the domain-wall fermion [53] with an improvement (the Möbius domain-wall fermions [54,55]) for the quark fields. Here we 3 Numerical works are done with the QCD software package IroIro++ [50] and Grid [51]. set the size of the fifth direction as L s = 16. The Möbius kernel is taken as where D W is the standard Wilson-Dirac operator with a large negative mass −1/a. In the following, we omit a, when there is no risk of confusion. The numbers without physical unit, e.g. MeV, are in the lattice unit. The fermion determinant thus obtained corresponds to a four-dimensional effective operator, Note that if we take the L s = ∞ limit, this operator converges to the overlap-Dirac operator [56] with a kernel operator H M . The stout smearing [57] is applied three times with the smoothing parameter ρ = 0.1 for the link variables in the Dirac operator.
The residual mass for our main runs at β = 4.30 is 0.14(6) MeV, which may look small enough for the measurements in this work. However, as we reported in our previous studies, the chiral symmetry breaking due to lattice artifacts is enhanced at high temperature, in particular for the axial U (1) susceptibility. Therefore, we employ an improved Dirac operator, exactly treating the near-zero modes of H M to compute the sign function. With the near-zero modes whose absolute value is less than 0.24 (∼ 630 MeV) treated exactly, we find that the residual mass becomes negligible ∼ 10 −5 or ∼ 0.01 MeV. Although it is different from the original definition in [56], we call this improved Dirac operator the "overlap" Dirac operator D ov (m). Note that our overlap operator has a finite L s , i.e. L s = 16. As we also found in our previous study that the lattice artifact from mixed action is large even For each ensemble, we simulated more than 20000 trajectories from which we carry out measurements on the configurations separated by 100 trajectories. We then bin the data in every 1000 trajectories, which is longer than autocorrelation lengths we observe. When we use the OV/MDW reweighting, we lose some amount of statistics due to its noise. An estimate for the effective number of statistics N rew in the table is defined as where R is the reweighting factor or the ratio of the determinant and R max is its maximal value in the ensemble. As discussed in [30], R has a negative correlation with the axial U (1) observables, and the statistical error estimated by the jackknife method are smaller than

III. NUMERICAL RESULTS
In this section, we show our numerical results on the axial U (1) anomaly. For meson and baryon correlators, we also present the tests of the SU (2) L × SU (2) R symmetry for comparison.

A. Dirac spectrum
The spectral density of the Dirac operator in volume V is used as a probe of chiral symmetry breaking. Here, λ i denotes the i-th eigenvalue of the Dirac operator on a given gauge configuration and · · · is the gauge ensemble average. In the chiral limit after taking the thermodynamical limit, we obtain the Banks-Casher relation [4], which relates the chiral condensate qq to the spectrum at λ = 0: We expect that ρ(0) = 0 above the critical temperature of the chiral phase transition.
As the vacuum expectation value (vev) qq also breaks the axial U (1) symmetry, the details of ρ(λ) in the vicinity of zero is important in this work. It has been shown that if the spectrum has a finite gap at λ = 0, the anomaly becomes invisible in mesonic two-point functions [2,3]. In Ref. [5], it is argued that the SU (2) L × SU (2) R symmetry restoration requires ρ(λ) ∼ λ α with the power α > 2, which is sufficient to show the absence of the axial U (1) A anomaly in multi-point correlation functions of scalar and pseudoscalar operators.
We measure the eigenvalues/eigenfunctions of the Hermitian four-dimensional effective Dirac operator H DW (m) = γ 5 D 4D DW (m) and those of the corresponding overlap-Dirac operator. Each eigenvalue λ m is converted to the one of massless Dirac operator by Forty lowest (in their absolute value) modes are stored for gauge configurations separated by 100 trajectories. They cover a range from zero to 300-500 MeV.
In Fig. 1, we present the Dirac eigenvalue density of the overlap-Dirac operator at T = 220 MeV at L = 32. Here the OV/MDW reweighting is applied. Compared to the solid line which represents the chiral condensate at T = 0 [59], the low-modes are suppressed by an order of magnitude. We observe a sharp peak near zero, which rapidly disappears as quark mass decreases, as expected from the SU (2) L × SU (2) R symmetry restoration. Since we find no clear gap, a region of λ where ρ(λ) = 0, we are not able to conclude if the axial U (1) symmetry is recovered or not from this observable only.
We also compare these results with the Dirac spectral density of the Möbius domain-wall fermion (dashed symbols) in Fig. 1. The agreement is remarkable. On the other hand, when we switch off the OV/MDW reweighting, which we call the non-reweighted overlap fermion setup, we observe a remarkable peak at the lowest bin, as presented in Fig. 2. Similar peaks were reported in [31][32][33], with overlap fermion only in the valence sector. Our data clearly show that these peaks are due to the lattice artifact of the mixed action.
In order to grasp possible systematic effects due to the finite lattice size, we compare the accumulated histogram at three different volumes, L = 24 (1.8 fm), 32 (2.4 fm), 40 (3.0 fm) in Fig. 3. Except for L = 24 at m = 0.01, whose aspect ratio LT is 2, which is smallest among the data sets, no clear volume dependence is seen. The point m = 0.01 is the heaviest mass in our simulation set, which is expected to be least sensitive to the volume, but the Dirac low-mode density is rather high, and some remnants of spontaneous SU (2) L × SU (2) R breaking and associated pseudo Nambu-Goldstone bosons may be responsible for this volume dependence.
We summarize the results at different temperatures in Fig. 4. The higher the temperature, the stronger the suppression of the low modes is. We find a good consistency between the Möbius domain-wall and reweighted overlap results. We also find that β = 4.30 and 4.24 results are consistent. The quark mass dependence is not very strong except for the lowest bin near λ = 0.
Finally, the quark mass dependence of the eigenvalue density of the reweighted overlap operator near λ = 0 (with the bin size ∼ 10 MeV) is presented in Fig. 5. Here both the chiral zero modes and non-chiral pair of near zero-modes are included. It is remarkable that the density of near-zero modes at different temperatures show a steep decrease towards the massless limit and becomes consistent with zero already at finite quark masses. As will be shown below, this behavior of near zero modes is strongly correlated with the signals of the axial U (1) anomaly.
Quark mass dependence of the Dirac eigenvalue density at (near) zero with the bin size ∼ 10 MeV. All the data look consistent with zero before reaching the chiral limit.

B. Topological susceptibility
In order to quantify the topological excitations of gauge fields, we measure the topological charge Q in two different ways. One is the index of the overlap-Dirac operator, or the number of zero modes with positive chirality minus that with negative chirality. For this, we perform the OV/MDW reweighting to avoid possible mixed action artifacts, which is found to be quite significant. The other is a gluonic definition, measured directly on the Möbius domain-wall ensemble, using the clover-like construction of the gauge field strength  not precise enough to determine if χ t goes to zero at finite quark mass, as predicted in [5].
Gluonic on DW T=330MeV Gluonic on DW  i.e. the axial U (1) susceptibility, defined by the difference between the pseudoscalar (π) and scalar (δ) correlators integrated over spacetime, where the ensemble average is taken at a finite quark mass m. For the overlap-Dirac operator, we can express ∆(m) using the spectral decomposition: where λ m 's are the eigenvalues of H ov (m) ≡ γ 5 D ov (m).
In our previous study [30], we found that the contribution from chiral zero modes is quite noisy. As an alternative, we definē for which the chiral zero mode contributions are subtracted. This is justified because in the thermodynamical limit, while at a fixed temperature, |Q| scales as V 1/2 (= L 3/2 ), and thus the zero-mode contributions vanish in the large volume limit as 1/V 1/2 . We numerically confirm this scaling at T = 220 MeV as presented in Fig. 9. The L 3/2 scaling of |Q| looks saturated for 1/(T L) < 0.4. Therefore,∆(m) in the thermodynamical limit coincides with ∆(m).
In this work, we further refine the observable by removing the UV divergence. From a simple dimensional analysis of the spectral expression in Eq. (8), the valence quark mass m v dependence of∆(m) can be expanded as and C has a logarithmic UV divergence. What we are interested in is the divergence free piece A m 2 v + B and if it is zero or not in the chiral limit. Note that A and B contain a sea quark mass m sea dependence, and A, in particular, should be suppressed as m 2 sea , at least, to avoid possible IR divergence in the limit of m sea = m v → 0. Measuring∆(m) at three different valence masses m 1,2,3 , we extract the UV finite quantity: We compute∆(m) through the expressions in Eqs. (8) and (9) truncating the sum at a certain upper limit λ cut (around 180 − 500 MeV). We then use Eq. (11) to obtain the UV subtracted susceptibility. Figure   This strong suppression is also seen at different temperatures, as presented in

D. Meson correlators
In the previous subsection we studied the difference between the pseudoscalar and scalar two-point correlation function integrated over the whole lattice. In this subsection, we investigate the mesonic two point correlation function itself, which must contain more microscopic information of QCD. We measure the spatial correlator in the z direction with O Γ =q τ Γq. Here τ are the generators in the flavor space. For Γ we choose γ 5 (P S), 1(S), γ 1,2 (V ), γ 5 γ 1,2 (A), γ 4 γ 3 (T t ) and γ 5 γ 4 γ 3 (X t ). In this work, we focus on T t , and X t channels, which are related by the axial U (1) transformation, as well as the V and A channels to check the recovery of the SU (2) L × SU (2) R symmetry. We find that the S correlator is too noisy to extract the "mass" and compare it with that of P S. For other channels, we will report elsewhere.
Since this quantity represents shorter distance physics than the axial U (1) susceptibility obtained from integration over whole lattice, the violation of the Ginsparg-Wilson relation enhanced by near-zero modes found in [30] is less severe. Therefore, we employ the Möbius domain-wall fermion formalism without reweighting. In order to improve the statistics, For the channels other than S, the signal is good enough to extract the asymptotic behavior of C Γ (z) at large z, which should contain the information of the screening mass.
However, the correlator at high temperatures may not behave like a single exponential even at large z. To circumvent this, recent studies often apply multi-state fits and introduce various source types to extract ground-state values, e.g. in [37].
The asymptotic behavior depends on the structure of the spectral functions (in the spatial direction) for each channel Γ: When the spectral function ρ Γ (ω) starts from a series of delta functions, which represents isolate poles, C Γ (z) at large z is dominated by a single exponential. On the other hand, if the correlator is described by deconfined two quarks, ρ Γ (ω) is a continuous function of ω, provided that the volume is sufficiently large. Let us assume that ρ Γ (ω) becomes nonzero at a threshold 2m, wherem is a constituent screening quark mass. For large z, we can expand ρ Γ (ω) as θ(ω − 2m)(c 0 + c 1 (ω − 2m) + · · · ) (with a step function θ), which results 1/z 2 )). In the Appendix A, we show that this form of the spectral function is indeed realized in the free two quark propagators. The essential difference from the single exponential is, thus, the factor 1/z. In this study we therefore apply two types of fitting functions: the standard cosh function, and the two-quark-inspired function (2q), In the free quark limit, we obtained more complete form of the two-quark propagations for each channel [46]. Note that in this limit,m = πT , which is the lowest Matsubara mass. It is, therefore, interesting to see how much m Γ is close to 2πT at our simulated temperature.
In Fig. 13, we present some of typical effective mass plots for m Γ and m Γ at T = 220 MeV to make a comparison between the single cosh function (square symbols) Eq. It is interesting to see that the plateaux are located at a mass lower than 2πT which is indicated by dashed lines. We also note that near the middle point ∼ L/2, the effective mass substantially drops. This drop may indicate a constant contribution to the correlator, which is caused by one quark propagating normally but the other going backward to wrap the torus, which is not taken into account in our formula Eq. (15). While we plan to give a detailed analysis in another publication [62], we show the screening mass m obtained by a fit to the 2q formula Eq. (15). The results are listed in Tab. II.
In Fig. 14, we plot the difference of the fitted m between T t and X t correlators at T = 220 MeV at different volumes. They are connected by the U (1) A rotation, and therefore, the difference, denoted by ∆m screen , is a probe of the axial U (1) symmetry. For the reference, we also plot in Fig. 15 the results for the difference between A and V channels, which is an indicator for the SU (2) L × SU (2) R symmetry.
Although the U (1) data at heavier quark masses are noisier than those for SU (2) L × SU (2) R , their chiral limit looks consistent with zero, and the central values are only a few MeV, at the lightest quark mass. We note that their individual mass is ∼ 1 GeV. Therefore, the axial U (1) symmetry relation is satisfied at a sub-% level. This behavior is also seen at different temperatures, as shown in Fig. 16, except for T = 190 MeV (but they are still consistent with zero with large error bars). The disappearance of the axial U (1) anomaly is consistent with other observables obtained using the reweighted overlap fermions.
Along the lines of the two-quark-inspired function for mesons, we use a three-quarkinspired function to extract a screening mass m j for each channel by fitting the forward propagating states.
This procedure gives qualitatively the same picture as in the previous section: a more stable plateau located at a lower energy value than that from the single cosh function. We therefore   use m to compare masses of different channels. Contrary to the case of mesons, we measure the baryon correlation functions exclusively in z-direction without low-mode averaging.
The numerical results presented in Fig. 17 show the difference of m between the axial U (1) partners. While the uncertainty grows with increasing quark mass, the signal is consistent with zero for all volumes at T =220 MeV. A similar restoration pattern is seen for SU (2) L × SU (2) R symmetry, as shown in Fig. 18, albeit with less fluctuations.

IV. CONCLUSION
In this work, we simulated two-flavor lattice QCD and tried to quantify how much of the axial U (1) anomaly survives at high temperatures 190-330 MeV. We employed the Möbius domain-wall fermion action and the overlap fermion action whose determinant is obtained by a stochastic reweighting technique. We fixed the lattice spacing to 0.074 fm, and chose more than four quark masses, including one below the physical point.
We confirmed that our data are consistent with those in the previous works [30] extending statistics of ensembles at β = 4.24. We also observed a good consistency between the Möbius domain-wall and overlap fermions, except for the axial U (1) susceptibility, which is very sensitive to the violation of the chiral symmetry at T = 220 MeV. The discretization effect is therefore well under control. We also confirmed that the systematics due to finite size of the lattice is under control. Our data with various lattice sizes agree, except for those with L = 24, which has a small aspect ratio T L = 2 at T = 220 MeV.
In the Dirac spectrum we found a strong suppression of low but non-zero eigenmodes.
The higher temperature, the more suppression of the low lying modes observed. On the other hand, for the chiral zero mode, a peak is found at all four simulated temperatures but its quark mass dependence is steep and the chiral limit is consistent with zero.
As expected from the behavior of the chiral zero mode, a sharp disappearance of the topological susceptibility is found, which suggests a mass dependence starting with a power ∼ m 4 near the chiral limit. Our numerical data for the axial U(1) susceptibility, meson and baryon correlators also indicate that the axial U (1) anomaly is consistent with zero in the chiral limit. From these observations we conclude that the remaining anomaly of the axial U (1) symmetry at the physical point for T ≥ 1.1 T c is at most a few MeV level, which is ∼ 1% of the simulated temperatures.
To examine the disappearance of the U (1) A anomaly near T c , we need a simulation around the critical temperature, which is beyond the scope of this paper. In this appendix, we compute propagators for two and three non-interacting quarks in ddimensions (to show that d = 4 is special). To take the finite temperature into account, the spacetime is assumed to be an Euclidean flat continuum space with one direction compactified. Namely, we consider S 1 × R d−1 , and anti-periodic boundary conditions are imposed on the fermions. We denote the compact direction by x 0 and consider spatial propagators in the x 1 direction.
Noting that the 0-th component of p denoted by p 0 is discrete, and neglecting higher p 0 contribution except for the lowest Matsubara frequency p 0 = M = ±πT , we can use the following approximation where the (d − 2)-dimensional vector q is given by q = (p 2 , p 3 , · · · p d−1 ) and the factor two comes from the two possible signs of M . Changing the variables as P 1 = p 1 − p 1 , R 1 = (p 1 + p 1 )/2 and explicitly integrating over R 1 and P 1 , we obtain where the constant C comes from the solid angle integral. In the last line we have changed the integral variable to ω = 2 M 2 + q 2 . Note that a fractional power is absent for d = 4.
From the above integral, we can read off the spectral function for d = 4 as which supports our assumption for the fitting form Eq. (15) of the meson correlators. Here we have chosen the pseudoscalar correlators but it was confirmed in [46] by a full computation including higher Matsubara frequencies that this asymptotic form is universal in all other channels.
From the above integral, we can read off the spectral function for d = 4 as with a constant (matrix) D, which supports the asymptotic three-quark form exp(−3M x 1 )/x 2 1 corresponding to Eq. (17).