Eigenvalue spectra of QCD and the fate of U_A(1) breaking towards the chiral limit
EEigenvalue spectra of QCD and the fate of U A (1) breaking towards the chiral limit Olaf Kaczmarek,
1, 2
Lukas Mazur, and Sayantan Sharma Fakult¨at f¨ur Physik, Universit¨at Bielefeld, D-33615 Bielefeld, Germany Key Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China The Institute of Mathematical Sciences, HBNI, Chennai 600113, India
The finite temperature phase diagram of QCD with two massless quark flavors is not yet under-stood because of the subtle effects of anomalous U A (1) symmetry. In this work we address this issueby studying the fate of the anomalous U A (1) symmetry in 2 + 1 flavor QCD just above the chiralcrossover transition temperature T c , lowering the light quark mass towards the chiral limit alongline of constant physical strange quark mass. We use the gauge configurations generated using theHighly Improved Staggered Quark (HISQ) discretization on lattice volumes 32 × × (cid:38) . U A (1) breaking observables on the light quarkmass, our study suggests U A (1) is broken at T (cid:38) T c even when the chiral limit is approached. PACS numbers: 12.38.Gc, 11.15.Ha
I. INTRODUCTION
Understanding the phase diagram of strongly inter-acting matter described by Quantum Chromodynamics(QCD) is one of the most challenging areas of contem-porary research. Apart from theoretical understandingof how matter was formed in the early universe, it alsogives a glimpse how strong interactions control the na-ture and universality class of phase transitions. Fortwo massless flavors of quarks, the QCD Lagrangianhas a U L (2) × U R (2) chiral symmetry. The subgroup SU V (2) × SU A (2) × U V (1) is spontaneously broken to SU V (2) × U V (1) in the hadron phase giving rise to pi-ons which are much lighter than the nucleons. The axial U A (1) subgroup is not an exact symmetry in QCD butbroken due to quantum effects [1, 2]. This is essentiallya non-perturbative feature of QCD which arises due tostrong color interactions.Though an anomalous symmetry, the U A (1) can influ-ence the nature of chiral phase transition of QCD withtwo degenerate light quark flavors. From the renormal-ization group studies of model quantum field theorieswith the same symmetries as N f = 2 QCD, it is knownthat the existence of a critical point at vanishingly smallbaryon density depends crucially on the magnitude ofthe U A (1) anomaly breaking near the chiral symmetryrestoration temperature [3]. If the U A (1) is approxi-mately restored as the chiral symmetry restoration oc-curs, then the phase transition from the hadron to thequark-gluon plasma phase is expected to be either of firstorder [3, 4] or second order of U L (2) × U R (2) /U V (2) uni-versality class [5, 6]. On the other hand, if the magnitudeof the U A (1) breaking term is comparable to its zero tem-perature value even at T c , then the phase transition is ofsecond order with O (4) critical exponents [3–5, 7].In such model quantum field theories however, thecoefficient of the U A (1) breaking term is a parame- ter whose magnitude can only be estimated from non-perturbative studies of QCD. Currently lattice regular-ization is the most practical method which can pro-vide a reliable answer to such a question. Lattice stud-ies in the recent years have provided some interestinginitial insights about the fate of the anomalous U (1)subgroup of the chiral symmetry in 2+1 flavor QCDwith physical quark masses. The eigenvalue spectraof the 2+1 flavor QCD Dirac operator using domainwall fermions [17–19], Highly improved staggered quark(HISQ) discretization [20, 21] and recently using thetwisted mass Wilson fermion discretization [22, 23] andrelated renormalized observables have reported substan-tial U A (1) breaking near and above the chiral crossoverregion. However many recent studies for two flavor QCD,with physical and heavier than physical light quarks (inthis limit strange quark is infinitely massive), using over-lap fermions [24], re-weighted spectra of the domain wallfermions [25–28], improved domain wall fermions [29] aswell as from non-perturbatively O ( a ) improved Wilsonfermions [30] have reported effective restoration of the U A (1) near T c . Whereas the typical volume of the lat-tice box reported in the later set of studies which ob-serve an effective U A (1) restoration near T c may still notbe large enough to contain enough topological fluctua-tions responsible for the U A (1) breaking, there can be dis-cretization effects in the QCD eigenvalue spectrum due tofinite lattice spacing even when using improved versionsof Wilson and staggered quarks. A more detailed un-derstanding of the near-zero mode spectrum of the stag-gered fermion (HISQ) configurations was achieved usingoverlap fermions [31] which show remarkable similarityto the pure staggered spectrum on finer lattices, closerto the continuum [32]. The near-zero modes contributedominantly towards the U A (1) breaking [19, 31]. Sinceneither of these studies have performed infinite volumeand the continuum limit extrapolations yet, these appar- a r X i v : . [ h e p - l a t ] F e b ent conflicting results are not surprising.Due to a finite light quark mass m l , the singular partof the free energy which carries information about theuniversal critical behavior gets overwhelmed by the reg-ular part which is an analytic function in m l [8]. Thusthere is no phase transition characterizing chiral symme-try restoration in QCD with physical light quark massesbut is rather a smooth crossover [9–15] which should goover to an exact phase transition only in the chiral limit.However, recent lattice studies have revealed remarkablesignatures of O (4) scaling in chiral observables as onelowers the light quark masses towards the chiral limitalong the line of constant physical strange quark mass[16]. Within the current precision, it is possible to ruleout Z (2) scaling [16], although to establish these resultsone ultimately needs to perform a continuum extrapola-tion.Though this study indirectly hints towards a scenariowhere U A (1) may not be effectively restored in the chirallimit, what happens to it in the chiral limit is not yetwell addressed using lattice studies. Recent results onthe topological susceptibility, which quantifies the U A (1)breaking in the chiral symmetry restored phase, showsa surprising trend for two flavor QCD as a function ofthe quark mass. It vanishes at some critical value of thelight quark mass [33, 34]. This result however needs tobe confirmed in the infinite volume and in the continuumlimit since it is not consistent with the other related re-sults in the context of N f = 2 phase transition withinthe so-called Columbia plot [35]. It is thus essential tolook again into the fate of U A (1) breaking when the lightquark masses are lowered towards the chiral limit, prefer-ably along a different line of constant physics within theColumbia plot. This is the main aim of this work, inwhich we keep the strange quark mass fixed to its phys-ical value and successively lower the light quark mass toeffectively approach the two flavor limit of QCD. Depend-ing on the effective breaking or restoration of U A (1), wewould either approach an O (4) or a Z (2) second orderline respectively.The paper is organized as follows. In the first sectionwe describe the numerical set-up and the novel techniqueswe will use to measure the QCD Dirac spectrum whenapproaching towards the chiral limit. In the subsequentsection, we study the eigenvalue spectrum for QCD withHISQ fermions near T c with a valence overlap Dirac op-erator and observe how the appropriately renormalizedspectrum changes when the light quark mass m l is re-duced towards the chiral limit. In particular, we look atthe quark mass dependence of the coefficient of the lead-ing O ( λ ) term of the renormalized eigenvalue density.In the chiral symmetry restored phase, from the chiralWard identities it is expected [36] that this coefficientvaries as m l in the leading order in quark mass. We findhowever that the coefficient has a m l -independent con-tribution at T > T c and its consequences for the U A (1)breaking is discussed. We next show our results for arenormalized observable m l ( χ π − χ δ ) /T which is sensi- T /T c m s /m l β N N τ N configs tive to U A (1), by appropriately tuning the valence andsea quark masses in order to ameliorate the effects of themixed fermion action. By performing a chiral extrapola-tion, we find a non-zero value of this observable, whichis responsible for the U A (1) breaking. We conclude bydiscussing the implications and a future outlook of ourpresent work. II. NUMERICAL SET-UP
The gauge configurations with 2+1 flavors of HISQ ac-tion used in this work were generated by the HotQCD col-laboration [15, 16]. We chose three different sets of gaugeensembles where the strange quark mass is set to its phys-ical value and the two light quark flavors are degeneratewith their mass varied such that m s /m l = 27 , ,
80. TheGoldstone pion masses corresponding to these choices ofthe light quark mass are ∼ , ,
80 MeV, respec-tively. The temperature range which we focus upon isbetween T c -1 . T c , where T c is the pseudo-critical tem-perature which is a function of the pion mass. The valuesof T c are ∼ , ,
154 MeV in the f K scale for pionmasses 135 , ,
80 MeV respectively. We have chosenthe lattice box with a temporal extent N τ = 8.The as-pect ratios are chosen such that ζ = N s /N τ = 4 for the m s /m l = 27 ,
40 configurations and ζ = N s /N τ = 7 forthe m s /m l = 80 gauge configurations. This ensures thatthe corresponding lattice extent along the spatial direc-tions are large enough m π L ∼ . . D ov [37] to measurethe eigenvalues of the HISQ sea configurations since it hasan exact index theorem on the lattice [38] and hence canresolve the small eigenvalues efficiently. Resolving theinfrared eigenvalue spectrum of the HISQ configurationswith a HISQ operator on relatively coarser lattices maybe tricky due to the breaking of continuum flavor sym-metries [21], however on finer lattices closer to the con-tinuum, a peak of near-zero modes is observed [32]. Thisinfrared peak can be very efficiently resolved using theoverlap operator on the HISQ sea configurations even oncoarser lattices [32]. We perform a proper tuning of thevalence overlap quark mass to the sea HISQ quark massand then measure appropriately renormalized eigenvaluespectrum and observables sensitive to the U A (1) break-ing, to ameliorate any effects of the mixed Dirac opera-tor set-up used here. The overlap operator for a masslessquark is defined as D ov = 1 + γ ε ( γ D W ) = 1 + D W (cid:113) D † W D W , (1)where D W is the standard massless Wilson Dirac opera-tor but including a constant term − M which is the defector the domain-wall height. The overlap operator was re-alized by calculating the sign function ε exactly in termsof the first 50 eigenvectors of the operator D † W D W andthen representing the contribution of the higher eigen-values through a Zolotarev Rational function with 25terms. The resultant norm of the square of sign func-tion, ε deviated from unity on average by about 10 − .The overlap operator satisfied the Ginsparg-Wilson re-lation upto a numerical precision of 10 − or even lower.We have performed a detailed study of the Ginsparg-Wilson relation violation and the appearance of near-zero modes discussed in the Appendix A. The domainwall height appearing in the overlap operator was cho-sen to be M = 1 . ,
110 MeV we have mea-sured the first 100 eigenvalues but increased the numberof eigenvalues to 200 for configurations with pion massof 80 MeV, because of increasing density of the low-lyingeigenvalues. The diagonalization of the Dirac operatorbecomes numerically quite expensive for the gauge en-sembles with lighter sea quark masses. This is due tothe fact that the number of zero modes increases andthey need to be calculated with very high precision. Wehave implemented a novel procedure to circumvent thisproblem which we describe in the following section.
A. Accelerating the overlap Dirac matrixdiagonalization towards the chiral limit
Since the overlap Dirac matrix, D ov is a normal ma-trix, the standard procedure is to diagonalize the Her-mitian operator D † ov D ov . The eigenvalues of this Hermi-tian operator comes in degenerate pairs having oppositechiralities. The zero modes however are non-degenerate with distinct chirality and their number and the chiralitydepends on the topological charge of the gauge configura-tions. We remind here that a significant time of the diag-onalization routine is spent in measuring the zero modeswith a reasonable precision. We therefore projected our D † ov D ov to measure only those eigenmodes which have achirality that is of opposite to those of the zero modes.The corresponding eigenspace has no zero modes and,leaving them out, accelerates the diagonalization routinesignificantly. The zero-modes do not contribute to thephysical observables in the thermodynamic limit, thusmeasuring the eigenspectrum without zero-modes, allowfor a significant speed-up of our calculations. This is spe-cially so for the gauge ensembles with sea-quark massesin the chiral limit when the probability of occurrence ofzero modes increases. We have explicitly checked on afew configurations, that for the lattice volumes we haveconsidered, the contribution from the zero modes to therenormalized observables is negligibly small.However, in order to project D † ov D ov onto a vectorspace which is devoid of its zero modes, we need to knowprecisely the chirality these zero modes. We estimatedthe chirality from the sign of the topological charge Q measured using its gluonic definition, Q = (cid:90) d x q ( x ) , q ( x ) = g π (cid:15) µνρσ Tr { F µν ( x ) F ρσ ( x ) } , where (cid:15) µνρσ is the totally anti-symmetric tensor and F µν is the non-Abelian field strength tensor. We have usedan O ( a )-improved lattice definition of the field strengthtensor [40], which greatly improves the precision of themeasurement of topological charge. Before measuring thetopological charge we have systematically smoothenedthe ultra-violet fluctuations of the gauge fields using thegradient flow [41, 42]. This involves introducing a fic-titious time direction denoted by t along which the fivedimensional gauge fields b µ evolve according to the fol-lowing flow equations and initial conditions, ∂B µ ∂t = D Bν G Bνµ , B µ ( t = 0 , x ) = A µ ( x ) . The flow smoothens the fields over a region of radius √ t .With increasing flow-time, ultra-violet noise of the non-Abelian fields gets increasingly suppressed leaving behindthe contribution of the topological modes. We have im-plemented the Zeuthen discretization [44] for numericallyimplementing the covariant derivative, which involves an O ( a ) Symanzik improvement [43]. The equation of mo-tion is integrated using an third order Runge-Kutta algo-rithm with an adaptive step-size. The topological chargefor all gauge ensembles have been measured at a flow-time √ tT = 0 .
45, the results of which are shown inFig. 1. We observe a optimal variation of the topologi-cal charge in all configurations that we have consideredfor different choices of the pion mass. This provides anevidence that we have chosen statistically independentconfigurations for our study, where Q is ergodically sam-pled.When calculating the overlap eigenvalues for the m s /m l = 80 ensemble, we have observed a significantslowing down of the algorithm to converge to the de-sired precision. This was due to the fact that theHISQ gauge configurations tend to have significantlymore small eigenvalues, some of which are localized onthe scale of the lattice spacing. In order to improvethe convergence, we systematically removed these ultra-violet defects by smoothening the configurations usingZeuthen flow technique up to a flow-time of t = 0 . a ,before measuring their eigenvalue spectrum with overlapfermions. The smoothening of the gauge fields has beenused earlier in the context of measuring the topologicalcharge [45] and the hadron spectrum using valence over-lap fermions [46, 47]. − Q T/T c = 0 . , m s /m l = 27 − Q T/T c = 1 . , m s /m l = 27 − Q T/T c = 1 . , m s /m l = 27 − Q T/T c = 0 . , m s /m l = 40 − Q T/T c = 1 . , m s /m l = 40 − Q T/T c = 1 . , m s /m l = 400 15 30 45 60 75 90 105 120 135 150 − Monte-Carlo historyQ
T/T c = 1 . , m s /m l = 80 Fig. 1: The Monte-Carlo time-history of the topologicalcharge Q as measured on the gauge configurations used inthis work. The Q is measured using a purely gluonic operatorat a flow-time √ tT = 0 .
45. The x-axis shows the config-uration number. The configurations are typically separatedby 50 hybrid Monte Carlo time steps. Different streams havebeen concatenated.
III. EIGENVALUE SPECTRUM OF QCD WITHHISQ FERMIONS TOWARDS THE CHIRALLIMIT
In this section we discuss the general features of theeigenvalue spectrum of the QCD Dirac operator, near and above the chiral crossover transition. The eigenvaluedensity ρ ( λ ) of the overlap Dirac operator measured onthe HISQ configurations are shown as a function of theimaginary part of the eigenvalue of the overlap operatordenoted as λ in Fig 2 for different temperatures near T c and physical quark masses. The exact zero modes arenot shown in this plot. The general features are simi-lar to that observed for the HISQ ensembles using theoverlap operator [31] with heavier than physical quarkmasses m s /m l = 20. As a comparison we also plot theeigenvalue density for lower than physical quark masses m l = m s /
40, for temperatures above T c in Fig. 3. Qual-itatively the features of the eigenvalue spectrum that weobserve for lighter quark masses are similar to physicalor heavier than physical quark mass. We have earlierpublished some preliminary results on this comparisonin Ref. [48].The eigenvalue density ρ ( λ ) of the QCD Dirac oper-ator can be characterized by a non-analytic dependenceon small eigenvalues which denotes the infrared peak,followed by a regular analytic dependence on λ whichis called the bulk spectrum [31]. The presence of thenear-zero peak is easily distinguishable from the bulkmodes in the chiral crossover region. Remarkably thispeak becomes more prominent as the temperature is in-creased gradually from T c since it is less contaminatedby the bulk modes, whose density shifts further towardsthe larger eigenvalues. In this context we should remindthat though the HISQ spectrum on coarser N τ = 8 lat-tices do not show any peak in the infrared [20], such apeak appears when one chooses more finer N τ = 16 lat-tices [32]. The use of overlap Dirac matrix as the valenceor probe operator on the HISQ sea ensembles corrects forthis lack of an exact index theorem for the HISQ opera-tor, extracting out the peak even on the coarser N τ = 8lattices. This peak appearing in the eigenvalue spectrumis thus not a lattice artifact, as discussed earlier in thecontext of domain wall fermions on relatively smaller lat-tice volumes [26, 27]. In fact such a peak can appear justabove T c due to an interacting ensemble of topologicalclusters [49] or in the high temperature phase due to adilute gas of instantons [18, 31, 50–52].We will study in detail how the slope of the bulk modesas well as the overlap between the near-zero and thebulk modes are sensitive to the change in temperature.Assuming correlations of upto 4-point functions in thepseudo-scalar and scalar meson channels to be analyticfunctions of m l , it was earlier derived using chiral Wardidentities, that in the chiral symmetry restored phase oftwo flavor QCD the analytic part of the Dirac eigenvaluedensity is ρ ( λ ) ∼ O ( m l ) λ + O ( m l ) λ + O ( m l ) λ + .. [36].This would imply that in the chiral limit the leading or-der behavior of the eigenvalue density is ρ ( λ ) ∼ λ . Withthis constraint it was shown explicitly that U A (1) break-ing is absent in upto six point correlation functions in thesame scalar and pseudo-scalar sectors [36]. It is there-fore important to investigate the analytic (bulk) part ofthe eigenvalue spectrum as a function of the sea-quark ρ ( λ ) / T λ /T β =6.390 β =6.445 β =6.500 Fig. 2: The eigenvalue density of the massless valence overlapDirac operator measured on HISQ sea configurations with m s /m l = 27, as a function of temperature, near and justabove T c . mass at different lattice spacings to understand the fateof U A (1) just above T c . Motivated from Ref. [31], we fitthe eigenvalue density at different temperatures to the fitansatz, ρ ( λ ) = ρ AA + λ + c ( m l )Θ( λ − λ ) | λ | γ ( m l ) , (2)where γ ( m l ) characterizes the leading order analytic de-pendence of the eigenvalue density on λ and can be ingeneral a function of m l . To extract the exponent γ , wechoose a threshold λ in the eigenvalues, beyond whichthe sensitivity to the non-analytic peak is minimum. Wehave implemented this through a Heaviside step functionΘ in the second term of the fit ansatz in Eq. 2. Wehave O (100) eigenvalues per configuration therefore wecan measure only the leading analytic behavior for theinfrared part of the eigenspectrum for λ > λ . The re-sults of the fit, including the values of the cut-off λ , theexponent γ for different quark masses and temperaturesand the corresponding goodness of the fits are summa-rized in Table II.For the temperature range we studied so far, the ex-ponent γ ∼ γ on the sea-quark mass at T c ≤ T ≤ . T c , consistent with the values obtained pre-viously on HISQ fermions [31] with heavier than physicalquark masses. In this context we would like to remindthat in Ref. [36], it was argued that the coefficient c ( m l )in Eq. 2 goes as m l in the chiral limit for two flavor QCD.This in-turn implies that one would not observe any lin-ear dependence on λ in ρ ( λ ) of QCD when chiral sym-metry is restored. In the context of the Columbia plot,the chiral limit in two flavor QCD is simply approachedalong the m s = ∞ line.In order to check this prediction we chose to insteadapproach the two flavor chiral limit along the line of con-stant physical value of m s and study the dependence of ρ ( λ ) / T λ /T β =6.390 β =6.423 β =6.445 Fig. 3: The eigenvalue density of the valence massless overlapDirac operator measured on the HISQ sea configurations with m s /m l = 40. m s /m l β T /T c λ /T γ χ / dof40 6.390 0.99 0.45 1.09(22) 0.7040 6.423 1.03 0.5 0.94(23) 0.9940 6.445 1.05 0.5 1.08(15) 0.6627 6.390 0.97 0.4 1.03(18) 0.6627 6.445 1.03 0.5 1.09(11) 0.9027 6.500 1.09 0.5 1.03(12) 0.94Tab. II: The temperature ( T ), the exponent γ characteriz-ing the leading order λ γ rise of the bulk eigenvalues λ andthe goodness of the fits performed on eigenvalue densities fordifferent choices of the light sea-quarks and physical value ofstrange sea-quarks. c ( m l ) on the light quark mass m l . We neglect the low-est β values for light sea-quark masses m s / , m s /
40 re-spectively since we want to be in the temperature regimewhere chiral symmetry is restored. For the dimensionlessratio ρ ( λ ) T = c ( m l ) T · λT , if c ( m l ) indeed goes as m l to lead-ing order in the sea-quark mass, then a fit to c ( m l ) /T asa function of m l /T should not have a non-zero intercept.The quantity c ( m l ) /T extracted from the eigenvaluedensities for T (cid:38) T c and different light quark masses areshown in Fig. 4. Clearly it has a constant intercept whichsurvives when chiral extrapolation is performed and evendominates over the usual O ( m l /T ) term. From this fitanalysis it is evident that the eigenvalue density, to theleading order is O ( λ ) rather than O ( λ ) just above above T c , even in the chiral limit. Moreover we do not observeany gap opening up in the infrared part of the eigenspec-trum, from which we can conclude that U A (1) remainsbroken as we approach the chiral limit. This is in addi-tion to the contribution to the U A (1) breaking that comesdue to the non-analytic peak in the eigenvalue spectrum.In the next section we will provide a more quantitativeestimate of the U A (1) breaking towards the chiral limit. c ( m l ) / T m l2 /T l2 /T +0.82(17) Fig. 4: The coefficient c ( m l ) of the λ dependent term of theeigenvalue density shown as a function of light sea quark masssquared m l . Both axes are normalized by appropriate powersof the temperature. The χ per degree of freedom of the fitis 0 . IV. QUANTIFYING U A (1) BREAKING IN THECHIRAL LIMIT
Having observed little sensitivity of the exponent γ of the bulk eigenvalue density to the sea-quark mass, itis interesting to compare the spectra at different quarkmasses and also with the earlier results obtained withoverlap fermions on HISQ configurations [31] for heavierthan physical quark masses. Since the eigenvalue den-sity is not an renormalization group invariant quantity,one has to renormalize the eigenvalue spectra for sucha comparison. Here we choose a renormalized eigen-value density using the valence strange quark mass. Forobtaining the quark mass, we tune the valence overlapstrange quark mass to the HISQ sea strange quark massby matching the renormalized quantity ∆ which is de-fined as, ∆ = m s (cid:104) ¯ ψψ (cid:105) l − m l (cid:104) ¯ ψψ (cid:105) s T . (3)The chiral condensates appearing in this observable,is first calculated in terms of the valence overlapDirac matrix. Using the definition (cid:10) ¯ ψψ (cid:11) ( m ) = TV (cid:68) tr (cid:16) D − m ∂∂D m (cid:17)(cid:69) , where D m = D ov (1 − am/ M ) + am is the overlap Dirac operator with a (valence) quark mass m , the condensates can be calculated from the overlapeigenvalues, a (cid:104) ¯ ψψ (cid:105) ( m ) = 1 N σ N τ (cid:20) (cid:104)| Q |(cid:105) am (4)+ (cid:42)(cid:88) ˜ λ (cid:54) =0 am (4 M − | ˜ λ | ) | ˜ λ | (4 M − ( am ) ) + 4( am ) M (cid:43) , where Q is the topological charge, and ˜ λ is the eigenvalue of D ov which is scaled by the defect height parameter M . We specifically neglect the first term in the right hand side ofthis expression which arises due to the zero modes, since itis a finite volume artefact to the above observable. We havechecked explicitly, that in all our gauge ensembles this termdue to the zero modes provides negligibly small correction to∆, implying that the finite volume effects are under controlin our tuning procedure. When using the observable ∆ fortuning m s , we keep the ratio m l /m s fixed for both the va-lence and the sea quarks. We then match the values of the∆ calculated from the first O (100) eigenvalues of the overlapDirac operator measured on the HISQ ensembles to the exactresults for ∆ obtained by the inversion of the HISQ operatorusing stochastic sources on the same ensembles in order toextract m s . Once the valence m s is tuned to its sea value,we can describe the physics of the underlying sea quarks us-ing the valence overlap fermions. The results for the tunedstrange quark masses for different ensembles are tabulated inTable III. m s /m l β m s sea m s val (∆)80 6.423 0.0670 0.02540 6.390 0.0694 0.09040 6.423 0.0670 0.05840 6.445 0.0652 0.03827 6.390 0.0694 0.09827 6.445 0.0652 0.05127 6.500 0.0614 0.032Tab. III: The valence strange quark masses obtained bymatching the observable ∆ measured using eigenvalue den-sity of the valence overlap to that measured by inversion ofthe sea HISQ Dirac operator.Subsequently, the comparison of the renormalized eigen-value density m s ρ ( λ ) /T as a function of λ/m s is shown inFig. 5. While calculating the renormalized density, the num-ber of bins for different ensembles were kept fixed. We observethat the near-zero peak of the renormalized spectrum showvery little sensitivity to the change in the quark mass. Thebulk modes again show a linear rise. As mentioned earlierthe spectrum for light quark mass m l = m s /
80 has increas-ing density of these small eigenmodes hence the spectrum isshown only till λ/m s ∼ . m s we next proceed to estimate whetherthe U A (1) is effectively restored above the crossover transi-tion and its sensitivity to the light quark mass. Howeversince U A (1) is not a symmetry there is no unique observablethat is sensitive to its restoration. The difference of the inte-grated two-point correlators of isospin-triplet pion ( ¯ ψτ γ ψ )and delta ( ¯ ψτ ψ ) mesons, ω = χ π − χ δ is one such observ-able that was proposed as a measure of U A (1) breaking [55].In fact one needs to further look at the degeneracy betweenhigher point correlation functions for different meson quan-tum number channels [36, 56, 57]. As a first test, we focus onthis specific two-point correlation function. In terms of theeigenvalues of the overlap Dirac operator it is defined as, ρ ( λ ) m s / T λ /m s m l =m s /80m l =m s /40m l =m s /27 Fig. 5: The renormalized eigenvalue density of the QCD en-sembles at T = 1 . T c generated using Highly improvedstaggered quark (HISQ) discretization, and measured usingan appropriately mass-tuned valence overlap operator. Theseare shown for three different choices of the light quarks. a ω ( m ) = 1 N σ N τ (cid:20) (cid:104)| Q |(cid:105) ( am ) (5)+ (cid:42)(cid:88) ˜ λ (cid:54) =0 am ) (4 M − | ˜ λ | ) (cid:104) | ˜ λ | (4 M − ( am ) ) + 4( am ) M (cid:105) (cid:43) . We have measured this quantity in terms of the first O (100)eigenvalues of the overlap Dirac operator at the values of thetuned valence quark masses. Chiral Ward identities ensurethat χ π − χ δ = χ disc , where χ disc is the disconnected part ofthe integrated iso-singlet scalar meson correlator. In order toverify the quality of our quark mass tuning we check whetherthis Ward identity is satisfied. We use the previously mea-sured data for χ disc for physical quark masses from Ref. [15],obtained by the inversion of the HISQ Dirac operator usingstochastic sources on N τ = 8 , ,
16 lattices and performa continuum extrapolation of this observable. We then com-pare our continuum estimate of m l χ disc /T , to the observable m l ω/T that we have calculated using the eigenvalues of thevalence overlap Dirac operator on the same HISQ ensemblesusing the tuned valence quark masses. The results for thiscomparative study are shown in Fig. 6. A reasonably goodagreement of these renormalized quantities are observed giv-ing us further confidence on our quark mass tuning procedure.Next we study the quark mass dependence of the appropri-ately renormalized U A (1) breaking observable m l ω/T . It isimportant to note here that we have calculated only the first O (100) of the total millions of eigenvalues of the QCD Diracoperator, nonetheless these infrared eigenvalues contribute significantly to the U A (1) breaking. U A (1) is effectively re-stored in the chiral limit if ω ∼ m l , i.e., if it has a trivialdependence on the light quark mass, as also expected in theperturbative regime. On the other hand if U A (1) is broken,then the leading order behavior of ω ∼ O ( m l ), which givesus the effective magnitude or the strength of U A (1) breaking.When chiral symmetry is effectively restored, we can use thefollowing fit ansatz for our data on m l ω/T , correspondingto these two scenarios, c m l2 ( χ π - χ δ )/T Cont. m l2 χ disc /T Fig. 6: Comparison of χ π − χ δ measured using the overlapeigenvalues to the continuum estimates of χ disc using datafrom Ref. [15], shown for physical quark masses. m l ( χ π − χ δ ) T = m l ωT = a m l T + a m l T , (6)= b m l T + b m l T . (7)where the former denotes U A (1) breaking whereas the lat-ter is valid on its effective restoration. We have calculatedthe m l ω/T at 1 . T c for three different tuned light valencequark masses, results of which are shown in Fig. 7. The datafits quite well to the first fit ansatz from Eq. 6, shown as a redband in the same figure, with the largest contribution to theuncertainty coming from the value corresponding to the low-est quark mass. Our data disfavors the second ansatz in Eq. 6since the corresponding χ per degree of freedom of the fit isabout 5, which is almost a factor 2 . χ π − χ δ ) /T in the chiral limit is 156 ±
13, which is clearly finite and non-zero within the current uncertainties. Thus we conclude thatfor the large volume N τ = 8 lattices we have studied so far,the U A (1) is broken above the chiral crossover temperature,even when we approach the chiral limit along the line of con-stant physical value of m s .Noting again that in the chiral symmetry restored phase thetopological susceptibility in QCD is related to m l ω [18, 58,59], it is evident that the former observable does not vanishwhen approaching the chiral limit along the line of constant m s . In contrast, for the two-flavor case, the topological sus-ceptibility is observed to vanish at a critical value of the lightquark mass, which is lower but close to the physical quarkmass [34]. V. CONCLUSIONS AND OUTLOOK
In this work we report on the eigenvalue spectrum and thefate of anomalous U A (1) symmetry in the chiral symmetryrestored phase of QCD on large volume N τ = 8 lattices aswe approach the chiral limit along the line of constant phys-ical strange quark mass. In order to correctly measure thenumber and density of the near-zero eigenmodes of the QCD m l ( χ π - χ δ ) / T m l /T Fig. 7: The renormalized U A (1) breaking parameter shown asa function of light quark mass at 1 . T c .ensembles using the Highly Improved Staggered quark (HISQ)discretization on the lattice, we use the overlap Dirac matrixas the valence or probe operator. This was done as the HISQoperator does not realize all the continuum flavor and anoma-lous symmetries on a finite lattice. In order to obtain physicaland renormalized results even with different valence and seaquark actions, we have reported a fairly easy procedure totune the valence quark masses to the sea quark mass usingthe eigenvalue spectrum as the input.Comparing the appropriately renormalized eigenvalue den-sity of the QCD Dirac operator, we observe that the eigen-value density at 1 . T c can be represented as ρ ( λ ) ∼ λ in thechiral limit. This is unlike the expectations for N f = 2 QCDin the chiral limit where the leading order behavior of theeigenvalue spectrum of QCD was derived to be ρ ( λ ) ∼ λ [36].This results in a non-zero value of the renormalized observ-able m l ( χ π − χ δ ) /T , which leads us to conclude that U A (1)is broken when one approaches the chiral limit along the lineof constant physical m s .For a final conclusive evidence, we need to re-visit thisstudy at several choices of lattice spacing and perform a con-tinuum extrapolation of our observables, which is computa-tionally expensive and would require several years of dedi-cated efforts. Furthermore as one approaches the chiral limit,the lattice volumes have to be chosen large enough so that thespatial extent is larger than the corresponding pion Comptonwavelength. In the present study we have chosen the latticevolumes keeping this criterion in mind and the m π L ∼ . U A (1) anomalous symmetry, just above the chiralcrossover for light quark masses as low as ∼ . O (4) second order line ofphase transitions. VI. ACKNOWLEDGMENTS
We would like to dedicate this work in the memory ofProf. Edwin Laermann, who initiated this research projectshortly before he passed away in August 2018. The au-thors acknowledge support by the Deutsche Forschungsge-meinschaft (DFG, German Research Foundation) through theCRC-TR 211 ’Strong-interaction matter under extreme con-ditions’– project number 315477589 – TRR 211. S.S. acknowl-edges support by the Department of Science and Technology,Govt. of India through a Ramanujan fellowship. The numer-ical computations in this work were performed on the GPUcluster at Bielefeld University and on Piz Daint at CSCS. Weacknowledge PRACE for awarding us access to Piz Daint atCSCS, Switzerland. Our GPU code is in part based on someof the publicly available QUDA libraries [60]. We thank theHotQCD collaboration, for sharing the gauge configurationswith us. S.S. is grateful to Prof. Frithjof Karsch for manyhelpful discussions related to this work.
Appendix A: Near-zero modes and Ginsparg-Wilsonrelation
In this section we provide a detailed analysis of theGinsparg-Wilson violation due to the numerical implemen-tation of the overlap Dirac operator and whether it leads tothe appearance of spurious near-zero modes in its eigenvaluespectrum. We will discuss here the results of our study fordifferent sea-quark masses at T = 1 . T c , for which we havemeasured the renormalized Dirac eigenvalue spectra shownin Fig. 5. It is evident from the eigenvalue spectrum, thatwe observe a non-analytic peak of small eigenvalues whichsurvive even when the quark masses are successively reducedtowards the chiral limit. For the lattice volumes we havestudied, the small eigenvalues did not appear sporadically fora few gauge configurations, rather all of them contributed tonear-zero peak of eigenvalues. Thus for a meaningful analysiswe instead chose to count the number of small eigenvaluesi.e., λ < λ appearing for each configuration and correlateit with the magnitude of Ginsparg-Wilson relation violationsuffered by the overlap Dirac operator implemented on thecorresponding gauge configurations.For m s /m l = 27 ,
40, we counted all overlap Dirac eigen-values for each configuration that are lower than λ = 0 . T .For m s /m l = 80 ensembles, we chose λ = 0 . T instead,since these have a more dense eigenspectrum in the infrared.The magnitude of Ginsparg-Wilson relation violation due tothe numerical implementation of the overlap Dirac matrix onthe same configurations were simultaneously measured andcompared to the probability of occurrence of small eigenval-ues. The results of this analysis is shown in Fig. 8. For the m s /m l = 80 ensembles, we have divided the number of smalleigenvalues appearing for each configuration by a factor of twoto fit in the scale of the Figure. We do not find any obviouscorrelation between violation of the Ginsparg-Wilson relationand the proliferation of the number of small eigenvalues of theoverlap Dirac operator.In fact, it is evident from Fig. 8 that the overlap Dirac ma-trix which had a larger count of small eigenvalues for partic-ular HISQ configurations, satisfied the Ginsparg-Wilson rela-tion to a relatively higher precision. We have also performed asimilar comparison to check if there is any correlation between -9 -9 -9 -9 -9 -9
10 15 20 25 30 35 40 G W p r e c i s i on No. of small eigenvaluesm s /27m s /40m s /80 Fig. 8: The magnitude by which the Ginsparg-Wilson relationis violated by the overlap Dirac operator is compared to thethe number of its small eigenvalues λ < λ measured on eachconfiguration. The results are shown for the gauge configura-tions at T = 1 , T c for m s /m l = 27 , ,
80 respectively. Thenumber of eigenvalues for m s /
80 have been scaled by 0 . -9 -9 -9
10 15 20 25 30 35 40 | ε - | No. of small eigenvaluesm s /27m s /40m s /80 Fig. 9: The magnitude of the numerical imprecision of thesquare of the overlap sign function compared to the the num-ber of small eigenvalues λ < λ of the overlap Dirac matrixmeasured on each configuration. The results are shown forthe ensembles at T = 1 , T c for m s /m l = 27 , ,
80 respec-tively. The number of eigenvalues for m s /
80 have been scaledby 0 . U A (1) breaking above T c are not unphysical lattice artifacts.[1] S. L. Adler, Phys. Rev. 177, 2426 (1969).J. Bell and R. Jackiw, Nuovo. Cim. A 60, 47 (1969).[2] K. Fujikawa, Phys. Rev. Lett. 42, 1195 (1979).[3] R. D. Pisarski and F. Wilczek, Phys. Rev. D 29, 338(1984).[4] A. Butti, A. Pelissetto, E. Vicari, JHEP 0308, 029 (2003).[5] A. Pelissetto and E. Vicari, Phys. Rev. D 88, 105018(2013).[6] Y. Nakayama and T. Ohtsuki, Phys. Rev. D 91, 021901(2015).[7] M. Grahl and D. H. Rischke, Phys. Rev. D 88, 056014 (2013).[8] S. Ejiri et. al., Phys. Rev. D 80 094505 (2009).[9] C. Bernard et. al., Phys. Rev. D 71, 034504 (2005).[10] M. Cheng et. al., Phys. Rev. D 74, 054507 (2006).[11] Y. Aoki et. al., Nature 443, 675-678 (2006).[12] A. Bazavov et. al., Phys. Rev. D 85, 054503 (2012).[13] T. Bhattacharya et. al., Phys. Rev. Lett. 113, 082001(2014).[14] A. Bazavov et. al. [HotQCD Collaboration], Phys. Rev.D 90, 094503 (2014).[15] A. Bazavov et al. [HotQCD Collaboration], Phys. Lett. B , 15 (2019).[16] H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch,A. Lahiri, S. T. Li, S. Mukherjee, H. Ohno, P. Petreczkyand C. Schmidt, et al. Phys. Rev. Lett. , no.6, 062002(2019).[17] S. Chandrasekharan, D. Chen, N. H. Christ, W. J. Lee,R. Mawhinney and P. M. Vranas, Phys. Rev. Lett. ,2463-2466 (1999)[18] A. Bazavov et. al. Phys. Rev. D 86, 094503 (2012).[19] M. I. Buchoff et. al., Phys. Rev. D 89, 054514 (2014).[20] H. Ohno, U. M. Heller, F. Karsch and S. Mukherjee, PoSLATTICE 2011, 210 (2011); PoS LATTICE 2012, 095(2012).[21] H. Ohno, U. M. Heller, F. Karsch and S. Mukherjee, PoS LATTICE2012 , 095 (2012) [arXiv:1211.2591 [hep-lat]].[22] F. Burger, E. M. Ilgenfritz, M. P. Lombardo andA. Trunin, Phys. Rev. D , no.9, 094501 (2018).[23] L. Holicki, E. M. Ilgenfritz and L. von Smekal, PoS LAT-TICE2018 , 180 (2018) [arXiv:1810.01130 [hep-lat]].[24] G. Cossu et. al., Phys. Rev. D 87, 114514 (2013).[25] A. Tomiya, G. Cossu, H. Fukaya, S. Hashimotoand J. Noaki, PoS
LATTICE2014 , 211 (2015)[arXiv:1412.7306 [hep-lat]].[26] A. Tomiya, G. Cossu, S. Aoki, H. Fukaya, S. Hashimoto,T. Kaneko and J. Noaki, Phys. Rev. D , no.3, 034509(2017).[27] K. Suzuki et al. [JLQCD], PoS LATTICE2018 , 152(2018) [arXiv:1812.06621 [hep-lat]].[28] K. Suzuki et al. [JLQCD], PoS
LATTICE2019 , 178(2020) [arXiv:2001.07962 [hep-lat]].[29] T. W. Chiu et al. [TWQCD], PoS
LATTICE2013 , 165(2014) [arXiv:1311.6220 [hep-lat]].[30] B. B. Brandt, A. Francis, H. B. Meyer, O. Philipsen,D. Robaina and H. Wittig, JHEP , 158 (2016).[31] V. Dick, F. Karsch, E. Laermann, S. Mukherjee andS. Sharma, Phys. Rev. D , no.9, 094504 (2015).[32] S. Sharma [HotQCD], [arXiv:1801.08500 [hep-lat]].[33] S. Aoki et al. [JLQCD], EPJ Web Conf. , 07024(2018).[34] S. Aoki et al. [JLQCD], [arXiv:2011.01499 [hep-lat]].[35] S. Sharma, PoS LATTICE2018 , 009 (2019)[arXiv:1901.07190 [hep-lat]].[36] S. Aoki, H. Fukaya, Y. Taniguchi, Phys. Rev. D 86,114512 (2012).[37] R. Narayanan and H. Neuberger, Phys. Rev. Lett. 71,3251 (1993); H. Neuberger, Phys. Lett. B 417, 141 (1998).[38] P. Hasenfratz, V. Laliena and F. Niedermeyer, Phys.Lett. B 427, 125 (1998).[39] T. Kalkreuter and H. Simma, Comput. Phys. Commun.93, 33 (1996).[40] S. O. Bilson-Thompson, D. B. Leinweber andA. G. Williams, Annals Phys. , 1-21 (2003)[arXiv:hep-lat/0203008 [hep-lat]].[41] R. Narayanan and H. Neuberger, JHEP , 064 (2006)[42] M. L¨uscher, JHEP , 071 (2010) [erratum: JHEP ,092 (2014)].[43] O. Kaczmarek, L. Mazur and S. Sharma, in preparation.[44] A. Ramos and S. Sint, Eur. Phys. J. C , no.1, 15(2016).[45] A. Hasenfratz and R. Hoffmann, Phys. Rev. D , 114509(2006).[46] A. Li et al. [xQCD], Phys. Rev. D , 114501 (2010).[47] M. Lujan, A. Alexandru, Y. Chen, T. Draper, W. Free-man, M. Gong, F. X. Lee, A. Li, K. F. Liu and N. Mathur,Phys. Rev. D , 014501 (2012).[48] L. Mazur, O. Kaczmarek, E. Laermann and S. Sharma,PoS LATTICE2018 , 153 (2019) [arXiv:1811.08222[hep-lat]].[49] T. Kanazawa and N. Yamamoto, Phys. Rev. D ,105015 (2015).[50] D. J. Gross, R. D. Pisarski and L. G. Yaffe, Rev. Mod.Phys. 53, 43 (1981).[51] R. G. Edwards, U. M. Heller, J. Kiskis and R. Narayanan,Phys. Rev. D 61, 074504 (2000).[52] H. T. Ding, S. T. Li, S. Mukherjee, A. Tomiya,X. D. Wang and Y. Zhang, [arXiv:2010.14836 [hep-lat]].[53] A. V. Smilga and J. Stern, Phys. Lett. B 318, 531 (1993).[54] J. J. M. Verbaarschot, Nucl. Phys. B 427, 534 (1994).[55] E. Shuryak, Comments Nucl. Part. Phys. 21, 235 (1994).[56] S. H. Lee and T. Hatsuda, Phys. Rev. D , 1871-1873(1996).[57] M. C. Birse, T. D. Cohen and J. A. McGovern, Phys.Lett. B , 137-140 (1996).[58] L. Giusti, G. C. Rossi and M. Testa, Phys. Lett. B ,157-166 (2004).[59] P. Petreczky, H. P. Schadler and S. Sharma, Phys. Lett.B , 498-505 (2016).[60] M. A. Clark, R. Babich, K. Barros, R. C. Brower andC. Rebbi, Comput. Phys. Commun.181