Topological Susceptibility in $N_f=2$ QCD at Finite Temperature

We study the topological charge in $N_f=2$ QCD at finite temperature using M\"obius domain-wall fermions. The susceptibility $\chi_t$ of the topological charge defined either by the index of overlap Dirac operator or a gluonic operator is investigated at several values of temperature $T\ (>T_c)$ varying the quark mass. A strong suppression of the susceptibility is observed below a certain value of the quark mass. The relation with the restoration of $U_A(1)$ is discussed.


Introduction
Topological susceptibility in QCD at finite temperature has acquired much attention recently due to its phenomenological interest. Mass of the QCD axion, one of the candidates of dark matter, is given by the topological susceptibility, and its dependence on temperature determines the abundance of the axion in the universe. A quantitative estimate can in principle be provided by lattice QCD, and was one of the topics of the panel discussion of this year's lattice conference [1][2][3][4]. This study is not meant to provide some quantitative results at phenomenologically important temperatures 500 T 1000 MeV [1], but rather to understand the nature of the phase transition in two-flavor QCD [4].
The fate of the U A (1) symmetry at and above the phase transition for vanishing u and d quark masses is one of the long standing and fundamental questions in QCD. While at any temperature the U A (1) chiral anomaly exists, manifestation of the U A (1) breaking is only possible if the gauge field configurations with non-trivial topology actually have non-vanishing contribution. The nontrivial QCD configurations also produce the topological susceptibility. Thus there is naturally a link in between these two physical quantities.
One powerful theoretical approach for these problems is to use the properties of the spectrum of the Dirac operator [5][6][7][8]. Along this line Aoki, Fukaya and Taniguchi (AFT) revisited the problem assuming the overlap fermions for the UV regulator for quarks [9]. They claim that the U A (1) symmetry in two flavor (N f = 2) QCD is recovered in the chiral limit for temperatures at and above the critical one. Furthermore, the derivatives of the topological susceptibility with respect to the quark mass m Speaker, e-mail: yasumichi.aoki@kek.jp vanish at any order. It means that the susceptibility, which is zero at the chiral limit, stays zero in the vicinity of m = 0. As the susceptibility is non-zero for infinitely heavy quarks, there must be a critical mass which divides the regions with zero and non-zero topological susceptibility.
The relation of the spectrum of the Dirac operator with U A (1) was also studied by Kanazawa and Yamamoto (KY) more recently [10]. Assuming the U A (1) breaking they derived a relation between the U A (1) susceptibility, which is a measure of the U A (1) breaking, and the topological susceptibility through a low energy constant. According to their study, the topological susceptibility should be proportional to the squared quark mass, thus, should exhibit quite different mass dependence to that of AFT. Kanazawa-Yamamoto claims the assumption that the spectral density is analytic near the origin in AFT needs to be abandoned to have the U A (1) breaking. The analyticity, however, seems intact in the simulations with overlap fermions [11] and domain wall fermions with overlap-reweighting [12][13][14], which have exact chiral symmetry.
Studying the topological susceptibility in depth would add another dimension for the understanding of the nature of the finite temperature transition in N f = 2 QCD, especially if it is done in conjunction with the direct measurement of the U A (1) breaking. Also, understanding the fate of the U A (1) breaking should be important for the computation of the topological susceptibility to a required precision necessary for phenomenology.
Chiral symmetry plays a crucial role for the study of U A (1) [12,13]. We use the Möbius domain wall fermion and reweighting method to the overlap fermion ensemble. In this report and the one for the U A (1) breaking [14], the main lattice spacing used is finer than we have used in [12,13]. This helps to reduce the residual chiral symmetry breaking of the domain wall fermions and to make the reweighting efficient.
The AFT scenario suggests a critical mass m c > 0 which divides the regions of topological charge zero and non-zero. If this is true it is consistent with the first-order phase transition [9,15], which was suggested by Pisarski and Wilczek [16] for the case of the U A (1) restoration. This could, then, change the widely believed phase diagram, called the Columbia plot at the upper-left corner. If similar dynamics exists at the physical strange quark mass point, it would affect the nature of the transition of the physical point depending on the value of m c .
This report is organized as follows. In Sec. 2, the calculation set-up and methods are described. Starting with a discussion on the sampling of the topological charge, an elaborate estimate of the error for the topological susceptibility is explained in Sec. 3, followed by our main results. Sec. 4 is devoted to summary and outlook. We use a = 1 units throughout. All the results reported here are preliminary.

Methods and parameters
Our simulation is carried out using Möbius domain wall fermions for two dynamical quark flavors [12]. A particular focus is placed on N t = 12, β = 4.3 ensembles with five different masses in this report. The corresponding temperature is T 220 MeV. At a fixed β value, two different temperatures N t = 8 and 10 are examined and results are reported. As a check of finite lattice spacing effects, a coarser lattice at β = 4.1 and N t = 8 corresponding to T 220 MeV is examined. For all lattices reported here the spatial site number is L = 32.
The lattice cutoff as a function of β for these lattices is obtained with the Wilson flow scale t 0 using the zero temperature results and an interpolation [17].
Topological susceptibility is defined as where V is the four dimensional volume and Q t is the topological charge.
We examine two definitions of the topological charge. One is the space-time sum of the gluonic topological charge density after the Symanzik flow at t = 5. The other is the index of the overlap-Dirac operator [17].
As pointed out in [12,17], it is essential to reweight to overlap ensemble from domain wall where R is the reweighting factor defined on each gauge field configuration, to correctly take into account the effect of (near) zero modes of the overlap-Dirac operator. Partial quenching by the use of valence overlap operators on dynamical domain wall ensembles leads to an artificial enhancement of low modes. The topological charge defined through the zero mode counting suffers from such artificial effects, which can be eliminated by the reweighting. We investigate two definition of the topological charge on the original domain wall ensemble and on the overlap ensemble generated through the reweighting. Altogether, four values of topological susceptibility are obtained at each parameter point as shown in the next section.
We are aiming to acquire the data from 30,000 molecular dynamics time units with hybrid Monte-Carlo simulation for each ensemble. Some of the reported data here are still undergoing improvement of statistics.  The left panel of Figure 1 shows the Monte-Carlo time history of the topological charge for β = 4.3 with N t = 12 (T 220 MeV) and bare mass m = 0.00375 ( 10 MeV) sampled every 20th trajectory. One trajectory amounts to a unit time molecular dynamics evolution followed by an accept-reject step. The red line corresponds to the charge measured with the gluonic definition ("GL"), while cyan represents that with the overlap index ("OV"). The legends also show the ensemble on which the   calculations are based, which are domain wall ("DW") for both. The right panel plots the histogram of the charge from "GL-DW" and that after the reweighting to the overlap ensemble "GL-OV". The bin size used can be read from the combined size of a pair of neighboring red and yellow bars. It shows there is not much difference between the data before and after the reweighting. Figure 2(a) shows the histogram of the topological charge measured through the overlap index before (OV-DW) and after (OV-OV) the reweighting. Here the width of the distribution shrinks significantly after the reweighting. This is due to the fact that the spurious zero modes on the domain wall ensemble gets suppressed. On the other hand, since such spurious zero modes are also suppressed by gauge field smearing, there appeared less difference between the gluonic measurements before and after the reweighting. From these data we calculate the topological susceptibility from Eq. (1).

Topological charge sampling and error estimate
Special attention is required when there is no weight for the non-trivial topology, shown in Fig. 2(b) as an example. The OV-OV histogram shows that all samples fall in the Q t = 0 sector. There actually is a non-zero |Q t | = 1 sample, but far smaller than the minimum of the y axis shown because of the small reweighting factor. As a result, the topological susceptibility is consistent with zero, with a jackknife error χ t = 4.4(4.4) × 10 2 MeV 4 . One should not take this as the sign of exact zero of χ t . This situation is similar to null measurements of rare processes in experiment. We estimate the upper bound of Q 2 t by imposing the condition that one measurement out of the full sample had |Q t | = 1 value. If the number of samples is N, then the upper bound of the topological susceptibility is With a reweighting, the effective number of samples gets reduced. We use the following quantity for the number of samples after reweighting: where R max is the maximum value of the reweighting factor in the ensemble [17]. As ∆ χ t can also be regarded as a resolution of the topological susceptibility given the number of samples -even if countable |Q| > 0 sector exists as in the Fig. 2(a) -we estimate the corrected statistical error of χ t for all the cases as where ∆ JK χ t is the jackknife error of χ t . For the case of Fig. 2(b), N eff = 32 out of a total of 1326 samples measured every 20th trajectory. Now the error after this correction reads ∆χ t = 3.9 × 10 6 MeV 4 .  The left panel of Fig. 3 shows the quark mass dependence of topological susceptibility for N t = 12 with T 220 MeV. The color coding used here is the same as in the history and histogram shown in Figs. 1 and 2. As noted in the previous section, OV-DW can yield enhanced fictitious zero-modes. Indeed, the cyan points appear as outliers and the resulting χ t gets fictitious enhancements. Also, as mentioned for m 10 MeV, the histograms of GL-DW and GL-OV are similar. Because of this, χ t for GL-DW and GL-OV appear consistent. As the reweighting reduces the effective number of statistics, we use GL-DW in comparison with GL-OV.

Topological susceptibility at T 220 MeV
The right panel of Fig. 3 shows χ t at m 6.6 MeV and T 220 MeV as a function of squared lattice spacing a 2 , where the finer lattice results are on the measured point and the coarser lattice results are obtained by linear-interpolation from the nearest two points 1 . The GL-DW result develops a large discretization error, and it gets close to OV-OV towards the continuum limit. The OV-OV result is more stable against lattice spacing. All results suggest χ t is vanishing in the continuum limit.
Focusing on the OV-OV result in the left panel the mass dependence of the topological susceptibility indicates two regions for mass: one is 0 < m 10 MeV where the observation of continuum scaling above strongly suggests χ t = 0. Actually, χ t with OV-OV is consistent with zero in this region. The other is m 10 MeV where χ t is significantly non-zero. We note that the existence of the boundary at non-zero m is also suggested from GL-DW. While χ t > 0 for 0 < m 10 MeV, it is almost constant. For χ t 10 MeV sudden development of χ t is observed. Due to its better precision over OV-OV, GL-DW results may be useful to identify the location of the boundary.
We note that a preliminary computation of the pion mass on the zero temperature configuration leads to an estimate of the physical ud quark mass as m = 4 MeV for the bare mass, which is well inside the region where χ t = 0 is suggested.   [9] (orange) and Kanazawa-Yamamoto [10] (brown). A zero-temperature result [18] for N f = 2 + 1 is plotted as a reference (green). Figure 4 shows a magnified view of the left panel of Fig. 3 without GL-OV and OV-DW. The newly added green line shows a zero temperature reference represented as a two-flavor ChPT fit with N f = 2 + 1 results [18]. In this figure two scenarios are compared: one is Aoki-Fukaya-Taniguchi [9] (AFT), where they claim that the derivatives of χ t with respect to quark mass vanish. With χ t = 0 at m = 0 a natural solution would be χ t = 0 for m < m c . The OV-OV result is consistent with this picture with 10 m c 12 MeV. The AFT result is based on the analyticity of Dirac eigenvalue spectral density ρ(λ). On the other hand, Kanazawa-Yamamoto [10] (KY) claims that U A (1) should be violated for T > T c due to its violation in the high enough temperature claimed in [19,20]. They reported that the analyticity of ρ(λ) needs to be abandoned for the U A (1) violation. There is a KY scenario of χ t (m) given in [10]. To evaluate χ t (m), one needs to know the value of a low energy constant, which may be extracted from the U A (1) order parameter measured with fixed topology. At the lightest mass where the topological charge is practically fixed at |Q t | = 0 after the reweighting (see Fig. 2(b)), we take their proposal with our U A (1) breaking parameter ∆ π−δ [14] and obtain the brown curve (∝ m 2 ) in the figure. This curve shows how χ t (m) behaves if the U A (1) symmetry were violated in the thermodynamic limit. Comparing with our OV-OV result, it has a tension (> 2σ) at m 13 MeV.

Topological susceptibility for T 220 MeV
To check whether the jump of the topological susceptibility observed at T 220 MeV persists at other temperatures, ensembles with different N t with fixed β = 4.3 have been generated and analyzed. The additional lattices are N t = 10 and 8 with fixed L = 32 as for N t = 12. The corresponding temperatures are T 264 and 330 MeV respectively. Figure 5 shows the topological susceptibility as a function of quark mass for three different temperatures, where only GL-DW data are shown. A similar jump of χ t at finite quark mass is observed also for T 264 and 330 MeV. The position of the jump shifts toward larger mass as T is increased.

Summary and outlook
Topological susceptibility χ t in N f = 2 QCD was examined at temperatures above the critical one with Möbius domain wall fermion ensembles reweighted to the overlap fermion ensembles. A special focus is put on the T 220 MeV ensembles with N t = 12. The preliminary results suggest that for the range of bare mass 0 ≤ m 10 MeV (which includes physical ud mass m 4 MeV) χ t = 0 and for m 10 MeV a sudden development of χ t starts. It is consistent with the prediction of Aoki-Fukaya-Taniguchi [9] with U A (1) symmetry restoration in the chiral limit, thus consistent with the direct measurement of the order parameter of U A (1) [14]. If that were due to finite volume effects and eventually we were to see the breaking in the thermodynamic limit, then Kanazawa-Yamamoto [9] explains how the U A (1) order parameter at finite volume is related to χ t . The result has a > 2σ tension.
We have examined the stability of the observation of the χ t = 0 region with a comparison to the coarse lattice result at approximately the same temperature, which indeed suggests the result is robust in the continuum limit. However this comparison is done with the lattice site number in the spatial direction fixed, therefore the physical box sizes are different. We are now examining the volume dependence on the finer lattice used in this report to check this.
Higher temperatures T 264 and 330 MeV are also studied with fixed lattice spacing a −1 2.64 GeV. A sudden change of χ t as a function of the quark mass is also observed for these temperatures. The point where the change occurs shifts towards larger mass for higher temperature. To get more insight for this observation, a systematic study for these high temperatures in conjunction with the U A (1) order parameter is planned.