# 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 ﬁnite temperature phase diagram of QCD with two massless quark ﬂavors is not yet under-stood because of the subtle eﬀects 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 ﬂavor 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 conﬁgurations 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 ﬂavors 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 eﬀects [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 inﬂu-ence the nature of chiral phase transition of QCD withtwo degenerate light quark ﬂavors. From the renormal-ization group studies of model quantum ﬁeld 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 ﬁrstorder [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 ﬁeld theories however, thecoeﬃcient 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 ﬂavor QCDwith physical quark masses. The eigenvalue spectraof the 2+1 ﬂavor 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 ﬂavor QCD,with physical and heavier than physical light quarks (inthis limit strange quark is inﬁnitely 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 eﬀective 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 eﬀective U A (1) restoration near T c may still notbe large enough to contain enough topological ﬂuctua-tions responsible for the U A (1) breaking, there can be dis-cretization eﬀects in the QCD eigenvalue spectrum due toﬁnite 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) conﬁgurations was achieved usingoverlap fermions [31] which show remarkable similarityto the pure staggered spectrum on ﬁner lattices, closerto the continuum [32]. The near-zero modes contributedominantly towards the U A (1) breaking [19, 31]. Sinceneither of these studies have performed inﬁnite volumeand the continuum limit extrapolations yet, these appar- a r X i v : . [ h e p - l a t ] F e b ent conﬂicting results are not surprising.Due to a ﬁnite 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 eﬀectively 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 quantiﬁes the U A (1)breaking in the chiral symmetry restored phase, showsa surprising trend for two ﬂavor QCD as a function ofthe quark mass. It vanishes at some critical value of thelight quark mass [33, 34]. This result however needs tobe conﬁrmed in the inﬁnite 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 diﬀerent line of constant physics within theColumbia plot. This is the main aim of this work, inwhich we keep the strange quark mass ﬁxed to its phys-ical value and successively lower the light quark mass toeﬀectively approach the two ﬂavor limit of QCD. Depend-ing on the eﬀective 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 ﬁrst 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 coeﬃcient 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 coeﬃcientvaries as m l in the leading order in quark mass. We ﬁndhowever that the coeﬃcient 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 eﬀects of themixed fermion action. By performing a chiral extrapola-tion, we ﬁnd 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 conﬁgurations with 2+1 ﬂavors of HISQ ac-tion used in this work were generated by the HotQCD col-laboration [15, 16]. We chose three diﬀerent sets of gaugeensembles where the strange quark mass is set to its phys-ical value and the two light quark ﬂavors 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 conﬁgurations and ζ = N s /N τ = 7 forthe m s /m l = 80 gauge conﬁgurations. 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 conﬁgurations since it hasan exact index theorem on the lattice [38] and hence canresolve the small eigenvalues eﬃciently. Resolving theinfrared eigenvalue spectrum of the HISQ conﬁgurationswith a HISQ operator on relatively coarser lattices maybe tricky due to the breaking of continuum ﬂavor sym-metries [21], however on ﬁner lattices closer to the con-tinuum, a peak of near-zero modes is observed [32]. Thisinfrared peak can be very eﬃciently resolved using theoverlap operator on the HISQ sea conﬁgurations 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 eﬀects of the mixed Dirac opera-tor set-up used here. The overlap operator for a masslessquark is deﬁned 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 ﬁrst 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 satisﬁed 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 ﬁrst 100 eigenvalues but increased the numberof eigenvalues to 200 for conﬁgurations 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 conﬁgura-tions. We remind here that a signiﬁcant 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 routinesigniﬁcantly. The zero-modes do not contribute to thephysical observables in the thermodynamic limit, thusmeasuring the eigenspectrum without zero-modes, allowfor a signiﬁcant 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 conﬁgurations, 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 deﬁnition, 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 ﬁeld strength tensor. We have usedan O ( a )-improved lattice deﬁnition of the ﬁeld strengthtensor [40], which greatly improves the precision of themeasurement of topological charge. Before measuring thetopological charge we have systematically smoothenedthe ultra-violet ﬂuctuations of the gauge ﬁelds using thegradient ﬂow [41, 42]. This involves introducing a ﬁc-titious time direction denoted by t along which the ﬁvedimensional gauge ﬁelds b µ evolve according to the fol-lowing ﬂow equations and initial conditions, ∂B µ ∂t = D Bν G Bνµ , B µ ( t = 0 , x ) = A µ ( x ) . The ﬂow smoothens the ﬁelds over a region of radius √ t .With increasing ﬂow-time, ultra-violet noise of the non-Abelian ﬁelds 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 ﬂow-time √ tT = 0 .

45, the results of which are shown inFig. 1. We observe a optimal variation of the topologi-cal charge in all conﬁgurations that we have consideredfor diﬀerent choices of the pion mass. This provides anevidence that we have chosen statistically independentconﬁgurations 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 signiﬁcantslowing down of the algorithm to converge to the de-sired precision. This was due to the fact that theHISQ gauge conﬁgurations tend to have signiﬁcantlymore 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 conﬁgurations usingZeuthen ﬂow technique up to a ﬂow-time of t = 0 . a ,before measuring their eigenvalue spectrum with overlapfermions. The smoothening of the gauge ﬁelds 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 conﬁgurations used inthis work. The Q is measured using a purely gluonic operatorat a ﬂow-time √ tT = 0 .

45. The x-axis shows the conﬁg-uration number. The conﬁgurations are typically separatedby 50 hybrid Monte Carlo time steps. Diﬀerent 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 conﬁgurations are shown as a function of theimaginary part of the eigenvalue of the overlap operatordenoted as λ in Fig 2 for diﬀerent 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 ﬁner 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 ﬂavor 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 conﬁgurations with m s /m l = 27, as a function of temperature, near and justabove T c . mass at diﬀerent lattice spacings to understand the fateof U A (1) just above T c . Motivated from Ref. [31], we ﬁtthe eigenvalue density at diﬀerent temperatures to the ﬁtansatz, ρ ( λ ) = ρ 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 ﬁt ansatz in Eq. 2. Wehave O (100) eigenvalues per conﬁguration therefore wecan measure only the leading analytic behavior for theinfrared part of the eigenspectrum for λ > λ . The re-sults of the ﬁt, including the values of the cut-oﬀ λ , theexponent γ for diﬀerent quark masses and temperaturesand the corresponding goodness of the ﬁts 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 coeﬃcient c ( m l )in Eq. 2 goes as m l in the chiral limit for two ﬂavor 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 ﬂavor QCD is simply approachedalong the m s = ∞ line.In order to check this prediction we chose to insteadapproach the two ﬂavor 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 conﬁgurations 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 ﬁts performed on eigenvalue densities fordiﬀerent 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 ﬁt 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 diﬀerent 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 ﬁtanalysis 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 coeﬃcient 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 ﬁtis 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 diﬀerent quarkmasses and also with the earlier results obtained withoverlap fermions on HISQ conﬁgurations [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-ﬁned 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 ﬁrst calculated in terms of the valence overlapDirac matrix. Using the deﬁnition (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 speciﬁcally neglect the ﬁrst term in the right hand side ofthis expression which arises due to the zero modes, since itis a ﬁnite 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 ﬁnite volume eﬀects are under controlin our tuning procedure. When using the observable ∆ fortuning m s , we keep the ratio m l /m s ﬁxed for both the va-lence and the sea quarks. We then match the values of the∆ calculated from the ﬁrst 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 diﬀerent 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 diﬀerent ensembles were kept ﬁxed. 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 eﬀectively 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 diﬀerence 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 diﬀerent meson quan-tum number channels [36, 56, 57]. As a ﬁrst test, we focus onthis speciﬁc two-point correlation function. In terms of theeigenvalues of the overlap Dirac operator it is deﬁned 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 diﬀerent 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 ﬁrst 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 satisﬁed. 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 conﬁdence 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 ﬁrst O (100) of the total millions of eigenvalues of the QCD Diracoperator, nonetheless these infrared eigenvalues contribute signiﬁcantly to the U A (1) breaking. U A (1) is eﬀectively 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 eﬀective magnitude or the strength of U A (1) breaking.When chiral symmetry is eﬀectively restored, we can use thefollowing ﬁt 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 eﬀective restoration. We have calculatedthe m l ω/T at 1 . T c for three diﬀerent tuned light valencequark masses, results of which are shown in Fig. 7. The dataﬁts quite well to the ﬁrst ﬁt ansatz from Eq. 6, shown as a redband in the same ﬁgure, 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 ﬁt isabout 5, which is almost a factor 2 . χ π − χ δ ) /T in the chiral limit is 156 ±

13, which is clearly ﬁnite 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-ﬂavor 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 ﬂavor and anoma-lous symmetries on a ﬁnite lattice. In order to obtain physicaland renormalized results even with diﬀerent 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 ﬁnal 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 eﬀorts. 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 conﬁgurationswith 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 fordiﬀerent 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 conﬁgurations, 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 conﬁguration and correlateit with the magnitude of Ginsparg-Wilson relation violationsuﬀered by the overlap Dirac operator implemented on thecorresponding gauge conﬁgurations.For m s /m l = 27 ,

40, we counted all overlap Dirac eigen-values for each conﬁguration 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 conﬁgurations 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 conﬁguration by a factor of twoto ﬁt in the scale of the Figure. We do not ﬁnd 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 conﬁgurations, satisﬁed 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 eachconﬁguration. The results are shown for the gauge conﬁgura-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 conﬁguration. 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. Buchoﬀ 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. Hoﬀmann, 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. Yaﬀe, 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