Eigenvalue spectra of QCD and the fate of $U_A(1)$ breaking towards the chiral limit

The finite temperature phase diagram of QCD with two massless quark flavors is not yet understood because of the subtle effects of anomalous $U_A(1)$ symmetry. In this work we address this issue by studying the fate of the anomalous $U_A(1)$ symmetry in $2+1$ flavor QCD just above the chiral crossover transition temperature $T_c$, lowering the light quark mass towards the chiral limit along line of constant physical strange quark mass. We use the gauge configurations generated using the Highly Improved Staggered Quark (HISQ) discretization on lattice volumes $32^3\times8$ and $56^3\times 8$ to study the renormalized eigenvalue spectrum of QCD with valence overlap Dirac operator. We have implemented new numerical techniques that have allowed us to measure about $100$-$200$ eigenvalues of the gauge ensembles with light quark masses $\gtrsim 0.6$ MeV. From a detailed analysis of the dependence of the renormalized eigenvalue spectrum and $U_A(1)$ breaking observables on the light quark mass, our study suggests $U_A(1)$ is broken at $T\gtrsim T_c$ even when the chiral limit is approached.

The finite temperature phase diagram of QCD with two massless quark flavors is not yet understood because of the subtle effects of anomalous UA(1) symmetry. In this work we address this issue by studying the fate of the anomalous UA(1) symmetry in 2 + 1 flavor QCD just above the chiral crossover transition temperature Tc, lowering the light quark mass toward the chiral limit along the line of constant physical strange quark mass. We use the gauge configurations generated using the highly improved staggered quark (HISQ) discretization on lattice volumes 32 3 × 8 and 56 3 × 8 to study the renormalized eigenvalue spectrum of QCD with valence overlap Dirac operator. We implement new numerical techniques that allow us to measure about 100-200 eigenvalues of the gauge ensembles with light quark masses 0.6 MeV. From a detailed analysis of the dependence of the renormalized eigenvalue spectrum and UA(1) breaking observables on the light quark mass, our study suggests UA(1) is broken at T Tc even when the chiral limit is approached.

I. INTRODUCTION
Understanding the phase diagram of strongly interacting matter described by Quantum Chromodynamics (QCD) is one of the most challenging areas of contemporary research. Apart from a theoretical understanding of how matter was formed in the early Universe, it also gives a glimpse of how strong interactions control the nature and universality class of chiral phase transitions. For two massless flavors of quarks, the QCD Lagrangian has 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 pions which are much lighter than the nucleons. The axial U A (1) subgroup is not an exact symmetry in QCD with two massless quark flavors but broken due to quantum effects [1,2]. This is essentially a non-perturbative feature of massless QCD which arises due to strong color interactions.
Though an anomalous symmetry, the U A (1) can influence the nature of the chiral phase transition of QCD with two degenerate light quark flavors. From the renormalization group studies of model quantum field theories with the same symmetries as N f = 2 QCD, it is known that the existence of a critical point at vanishingly small baryon density depends crucially on the magnitude of the U A (1) anomaly breaking near the chiral symmetry restoration temperature [3]. If the U A (1) is effectively restored as the chiral symmetry restoration occurs, then the phase transition from the hadron to the quark-gluon plasma phase is expected to be either of first order [3,4] * Electronic address: okacz@physik.uni-bielefeld.de † Electronic address: lmazur@physik.uni-bielefeld.de ‡ Electronic address: sayantans@imsc.res.in or second order of the U L (2) × U R (2)/U V (2) universality class [5,6]. On the other hand, if the magnitude of the U A (1) breaking term is comparable to its zero temperature value even at T c , then the phase transition is of second order with O(4) critical exponents [3][4][5]7].
In such model quantum field theories however, the coefficient of the U A (1) breaking term is a parameter whose magnitude can only be estimated from nonperturbative studies of QCD. Currently lattice regularization is the most practical method which can provide a reliable answer to such a question. Lattice studies in the recent years have provided some interesting initial insights about the fate of the anomalous U (1) subgroup of the chiral symmetry in 2+1 flavor QCD with physical quark masses. Using the eigenvalue spectra of the 2+1 flavor QCD Dirac operator using domainwall fermions [8][9][10], highly improved staggered quark (HISQ) discretization [11][12][13] and more recently using the twisted mass Wilson fermion discretization [14,15] and related renormalized observables, several studies have reported substantial U A (1) breaking near and above the chiral crossover region. However many recent studies for two-flavor QCD, with physical and heavier than physical light quarks (in this limit strange quark is infinitely massive), using overlap fermions [16], reweighted spectra of the domain-wall fermions [17][18][19][20], improved domain-wall fermions [21] as well as from non-perturbatively O(a) improved Wilson fermions [22] have reported effective restoration of the U A (1) near T c . Whereas the typical volume of the lattice box reported in the later set of studies which observe an effective U A (1) restoration near T c may still not be large enough to contain enough topological fluctuations responsible for the U A (1) breaking, there can be discretization effects in the QCD eigenvalue spectrum due to finite lattice spacing even when using improved versions of Wilson and staggered quarks. A more detailed understanding of the near-zero mode spec-trum of the staggered fermion (HISQ) configurations was achieved using valence overlap fermions [23] which show remarkable similarity to the pure staggered spectrum on finer lattices, closer to the continuum [24]. The near-zero modes contribute dominantly toward the U A (1) breaking [10,23]. Since neither of these studies have performed infinite volume and the continuum limit extrapolations yet, these apparent conflicting results are not surprising.
Due to a finite light quark mass m l , the singular part of the free energy which carries information about the universal critical behavior gets overwhelmed by the regular part which is an analytic function in m 2 l [25]. Thus there is no phase transition characterizing chiral symmetry restoration in QCD with physical light quark masses but is rather a smooth crossover [26][27][28][29][30][31][32][33] which should go over to an exact phase transition only in the chiral limit. However, recent lattice studies have revealed remarkable signatures of O(4) scaling in chiral observables as one lowers the light quark masses toward the chiral limit along the line of constant physical strange quark mass [34]. Within the current precision, it is possible to rule out Z(2) scaling [34], although to establish these results one ultimately needs to perform a continuum extrapolation.
Though this study indirectly hints toward a scenario where U A (1) may not be effectively restored in the chiral limit, what happens to it in the chiral limit is not yet well addressed using lattice studies. Recent results on the topological susceptibility, which quantifies the U A (1) breaking in the chiral symmetry restored phase, shows a surprising trend for two-flavor QCD as a function of the quark mass. It vanishes at some critical value of the light quark mass [35,36]. This result however needs to be confirmed in the infinite volume and in the continuum limit since it is not consistent with the other related results in the context of the N f = 2 phase transition within the so-called Columbia plot [37]. It is thus essential to look again into the fate of U A (1) breaking when the light quark masses are lowered toward the chiral limit, preferably along a different line of constant physics within the Columbia plot. This is the main aim of this work, in which we keep the strange quark mass fixed to its physical value and successively lower the light quark mass to effectively approach the massless two-flavor limit of QCD. Depending on the effective breaking or restoration of U A (1), we would either approach an O(4) or a Z(2) [may be even U L (2) × U R (2)/U V (2)] second order line respectively.
The paper is organized as follows. In the first section we describe the numerical set-up and the novel techniques we will use to measure the QCD Dirac spectrum when approaching the chiral limit. In the subsequent section, we study the eigenvalue spectrum for QCD with HISQ fermions near T c with a valence overlap Dirac operator and observe the dependence of the renormalized spectrum on light quark mass m l toward the chiral limit. In particular, we look at the quark mass dependence of the coefficient of the leading O(λ) term of the renormalized eigenvalue density in the bulk of the spectrum. In the chiral symmetry restored phase, from the chiral Ward identities it is expected [38] that this coefficient varies as m 2 l in the leading order in quark mass. We find however that the coefficient has an m l -independent contribution at T > T c and its consequences for the U A (1) breaking are discussed. We next show our results for a renormalized observable m 2 l (χ π − χ δ )/T 4 which is sensitive to U A (1), by appropriately tuning the valence and sea quark masses in order to ameliorate the discretization effects of the mixed-action set up. By performing a chiral extrapolation, we find a nonzero value of this observable, which is responsible for the U A (1) breaking. We conclude by discussing the implications and a future outlook of our present work.

II. NUMERICAL SET-UP
The gauge configurations with 2 + 1 flavors of HISQ action used in this work were generated by the HotQCD Collaboration [33,34]. We chose three different sets of gauge ensembles taken from Ref. [34] where the strange quark mass is set to its physical value and the two light quark flavors are degenerate with their mass varied such that m s /m l = 27, 40, 80. The Goldstone pion masses corresponding to these choices of the light quark mass are ∼ 135, 110, 80 MeV, respectively. The temperature range which we focus on is between T c -1.1 T c , which has been determined by representing the lattice spacing in terms of a physical scale, the kaon decay constant f K , for which we use the most recent parametrization from [13]. T c is the pseudo-critical temperature for the chiral crossover transition and is also sensitive to the pion mass. The values of T c for pion masses 135, 110, 80 MeV are T c ∼ 158, 157, 154 MeV respectively, which are extracted from the peak of the chiral susceptibility as a function of temperature [34]. We have chosen the lattice box with a temporal extent N τ = 8.The aspect ratios are chosen such that ζ = N s /N τ = 4 for the m s /m l = 27, 40 configurations and ζ = 7 for the m s /m l = 80 gauge configurations. This ensures that the corresponding lattice extent along the spatial directions is large enough m π L ∼ 2.7-3.5 to minimize the finite volume effects. The typical number of configurations analyzed at each temperature were obtained by removing first the 500-1000 configurations for thermalization and then choosing those which were separated by at least 50 hybrid Monte Carlo time steps to ensure minimal autocorrelation. The details of the configurations are given in Table I.
We use the overlap Dirac operator D ov [39] to measure the eigenvalues of the HISQ sea configurations since it has an exact index theorem on the lattice [40] and hence can resolve the small eigenvalues efficiently. Resolving the infrared eigenvalue spectrum of the HISQ configurations with a HISQ operator on relatively coarser lattices may be tricky due to the breaking of continuum flavor symmetries [12], however on finer lattices closer to the continuum, a peak of near-zero modes is observed [24]. This infrared peak can be very efficiently resolved using the overlap operator on the HISQ sea configurations even on coarser lattices [24]. We perform a proper tuning of the valence overlap quark mass to the sea HISQ quark mass in Sec. IV and then measure appropriately renormalized eigenvalue spectrum and observables sensitive to the U A (1) breaking, to ameliorate any effects of the mixed Dirac operator set-up used here. The overlap operator for a massless quark is defined as where D W is the standard massless Wilson Dirac operator but including a constant term −M which is the defect or the domain-wall height. The overlap operator was realized by calculating the sign function ε exactly in the subspace consisting of the first 50 eigenvectors of the operator D † W D W and then representing the contribution of the higher eigenvalues through a Zolotarev rational function with 25 terms. The resultant norm of the square of sign function, ε 2 deviated from unity on average by about 10 −10 . The overlap operator satisfied the Ginsparg-Wilson [41] relation to a numerical precision of magnitude 10 −9 or even lower. We have performed a detailed study of the Ginsparg-Wilson relation violation and the appearance of near-zero modes are discussed in the Appendix. The domain-wall height appearing in the overlap operator was chosen to be M = 1.8 since it gave the minimal violation of the Ginsparg-Wilson relation and approximated the sign function to the best numerical precision on the majority of the gauge configurations studied.
We then calculate the lowest eigenvalues of the overlap operator on the HISQ sea ensembles using the Kalkreuter-Simma Ritz algorithm [42]. For gauge ensembles with pion masses 135, 110 MeV we have measured the first 100 eigenvalues but increased the number of eigenvalues to 200 for configurations with pion mass of 80 MeV, because of the increasing density of the lowlying eigenvalues. The diagonalization of the Dirac operator becomes numerically quite expensive for the gauge ensembles with lighter sea-quark masses. This is due to the fact that the number of zero modes increases and they need to be calculated with very high precision. We have implemented a novel procedure to circumvent this problem which we describe in the following section.
A. Accelerating the overlap Dirac matrix diagonalization towards the chiral limit Since the overlap Dirac matrix D ov is a normal matrix, the standard procedure is to diagonalize the Hermitian operator D † ov D ov . The nonzero eigenvalues of this Hermitian operator come in degenerate pairs having opposite chiralities. The zero modes however are nondegenerate with distinct chirality and their number and the chirality depend on the topological charge of the gauge configurations. We remind here that a significant time of the diagonalization routine is spent in measuring the zero modes with a reasonable precision. We therefore projected our D † ov D ov to measure only those eigenmodes which have a chirality opposite to those of the zero modes. The corresponding eigenspace has no zero modes and, leaving them out, accelerates the diagonalization routine significantly. The zero modes do not contribute to the physical observables in the thermodynamic limit, thus measuring the eigenspectrum without zero-modes, allow for a significant speed-up of our calculations. This is especially so for the gauge ensembles with sea-quark masses in the chiral limit when the probability of occurrence of zero modes increases. We have explicitly checked on a few configurations, that for the lattice volumes we have considered, the contribution from the zero modes to the renormalized observables is negligibly small.
However, in order to project D † ov D ov onto a vector space which is devoid of its zero modes, we need to know precisely the chirality these zero modes. We estimated the chirality from the sign of the topological charge Q measured using its gluonic definition, where µνρσ is the totally anti-symmetric tensor and F µν is the non-Abelian field strength tensor. We have used an O(a 2 )-improved lattice definition of the field strength tensor [43], which greatly improves the precision of the measurement of topological charge. Before measuring the topological charge we have systematically smoothened the ultraviolet fluctuations of the gauge fields using the gradient flow [44,45]. This involves introducing a fictitious time direction denoted by t along which the fivedimensional gauge fields B µ evolve according to the following flow equations and initial conditions, The flow smoothens the fields over a region of radius √ 8t. With increasing flowtime, ultraviolet noise of the non-Abelian fields gets increasingly suppressed leaving behind the contribution of the topological modes. We have implemented the Zeuthen discretization [47] for numerically implementing the covariant (gauge) derivative, which involves an O(a 2 ) Symanzik improvement [46]. The equation of motion is integrated using a third order Runge-Kutta algorithm with an adaptive step size. The topological charge for all gauge ensembles has been measured at a flow time √ 8tT = 0.45, the results of which are shown in Fig. 1. We observe an optimal variation of the topological charge in all configurations that we have considered for different choices of the pion mass, without getting stuck in one topological sector. This provides evidence that we have chosen statistically independent configurations for our study, where Q is ergodically sampled.
When calculating the overlap eigenvalues for the m s /m l = 80 ensemble, we have observed a significant slowing down of the algorithm to converge to the desired precision. This was due to the fact that the HISQ gauge configurations tend to have significantly more small eigenvalues, some of which are localized on the scale of the lattice spacing. In order to improve the convergence, we systematically removed these ultraviolet defects by smoothening the m s /m l = 80 configurations using the Zeuthen flow technique up to a flow time of t = 0.32a 2 , before measuring their eigenvalue spectrum with overlap fermions. The smoothening of the gauge fields has been used earlier in the context of measuring the topological charge [48] and the hadron spectrum using valence overlap fermions [49,50].

III. EIGENVALUE SPECTRUM OF QCD WITH HISQ FERMIONS TOWARDS THE CHIRAL LIMIT
In this section we discuss the general features of the eigenvalue spectrum of the QCD Dirac operator, near and above the chiral crossover transition. The eigenvalue density ρ(λ) of the massless overlap Dirac operator measured on the HISQ configurations are shown as a function of the imaginary part of the eigenvalue of the overlap operator denoted as λ in Fig 2 for different temperatures near T c and physical quark masses. The exact zero modes are not shown in this plot. The general features are similar to that observed for the HISQ ensembles using the overlap operator [23] with heavier than physical quark masses m s /m l = 20. As a comparison we also plot the eigenvalue density for lower than physical quark masses m l = m s /40, for temperatures above T c in Fig. 3. Qualitatively the features of the eigenvalue spectrum that we observe for lighter quark masses are similar to physical or heavier than physical quark mass. We have earlier published some preliminary results on this comparison in Ref. [51].
The eigenvalue density ρ(λ) of the QCD Dirac operator can be characterized by a peak of near-zero modes which leads to a non-analytic dependence in the thermodynamic and continuum limits, followed by a regular dependence on λ which is called the bulk spectrum [23]. The presence of the near-zero peak is easily distinguishable from the bulk modes in the chiral crossover region. Remarkably this peak becomes more prominent as the temperature is increased gradually from T c since it is less contaminated by the bulk modes, whose density shifts further toward the larger eigenvalues. In this context we should remind that though the HISQ spectrum on coarser N τ = 8 lattices do not show any peak in the infrared [11], such a peak appears when one chooses finer N τ = 16 lattices [24]. The use of overlap Dirac matrix as the valence or probe operator on the HISQ sea ensembles corrects for this lack of an exact index theorem for the HISQ operator, extracting out the peak even on the coarser N τ = 8 lattices. This peak appearing in the eigenvalue spectrum is thus not a lattice artifact, as discussed earlier in the context of domain-wall fermions on relatively smaller lattice volumes [18,19]. In fact such a peak can appear just above T c due to an interacting ensemble of topological clusters [52] or in the high temperature phase due to a dilute gas of instantons [9,23,[53][54][55][56]. It has been recently argued that the existence of such localized nearzero modes in the chiral symmetry broken phase of QCD with massless quarks, can lead to the disappearance of Goldstone excitations at finite temperature [57]. We will study in detail how the slope of the bulk modes as well as the overlap between the near-zero and the bulk modes are sensitive to the change in temperature. Assuming correlations of up to 4-point functions in the pseudo-scalar and scalar meson channels to be analytic functions of m 2 l , it was earlier derived that the Dirac eigenvalue density for two flavor QCD is . [38]. The crucial assumptions that goes into this calculation are that the SU (2) L × SU (2) R is restored (hence a part of the chiral Ward identities are used) and that the eigenvalue density ρ(λ) is an analytic function of λ. This would imply that in the chiral limit the leading order behavior of the eigenvalue density is ρ(λ) ∼ λ 3 . With this constraint it was shown explicitly that U A (1) breaking is absent in at least 6-point correlation functions in the same scalar and pseudo-scalar sectors [38]. It is therefore important to investigate the bulk part of the eigenvalue spectrum as a function of the sea-quark mass non-perturbatively to understand the fate of U A (1) just above T c . Motivated from Ref. [23], we fit the eigenvalue density at different temperatures to the fit ansatz, where γ(m l ) characterizes the leading order dependence of the bulk eigenvalue density on λ and can be in general a function of m l . To extract the exponent γ, we choose a threshold λ 0 in the eigenvalues, beyond which the sensitivity to the infrared peak of near-zero eigenvalues is minimum. We have implemented this through a Heaviside step function Θ in the second term of the fit ansatz in Eq. 2. We have O(100) eigenvalues per configuration therefore we can measure only the leading behavior of the bulk eigenvalues. The results of the fit, including the values of the cut-off λ 0 , the exponent γ for different quark masses and temperatures and the corresponding goodness of the fits are summarized in Table II  For the temperature range we studied so far, the exponent γ ∼ 1 is consistent with the predictions from the chiral perturbation theory [58,59]. We do not observe any dependence of γ on the sea-quark mass at T c ≤ T ≤ 1.1 T c , consistent with the values obtained previously on HISQ fermions [23] with heavier than physical quark masses. In this context we would like to remind that in Ref. [38], it was argued that the coefficient c(m l ) in Eq. 2 goes as m 2 l in the chiral limit for two flavor QCD. This in-turn implies that one would not observe any linear dependence on λ in ρ(λ) of QCD when chiral symmetry is restored. In the context of the Columbia plot, the chiral limit in two flavor QCD is simply approached along the m s = ∞ line. ms/m l β T /Tc λ0/T γ χ 2 /d.o.f. 40 6.390 0.99 0.45 1.09 (22) 0.70 40 6.423 1.03 0.5 0.94 (23) 0.99 40 6.445 1.05 0.5 1.08 (15) 0.66 27 6.390 0.97 0.4 1.03 (18) 0.66 27 6.445 1.05 0.5 1.09 (11) 0.90 27 6.500 1.09 0.5 1.03 (12) 0.94 Tab. II: The temperature (T ), the exponent γ characterizing the leading order λ γ rise of the bulk eigenvalues λ and the goodness of the fits performed on eigenvalue densities for different choices of the light sea-quarks and physical value of strange sea-quarks.
In order to check this prediction we chose to instead approach the two flavor chiral limit along the line of constant physical value of m s and study the dependence of c(m l ) on the light quark mass m l . Since the line of first order transitions for three degenerate quark flavors is very tiny and survive for quark masses which are much smaller than the physical masses, the physics along these two lines of constant m s should not be much different. We neglect the lowest β values for light sea-quark masses m s /27, m s /40 respectively since we want to be in the temperature regime where chiral symmetry is restored. For the dimensionless ratio ρ(λ) T 3 = c(m l ) T 2 · λ T , if c(m l ) indeed goes as m 2 l to leading order in the sea-quark mass, then a fit to c(m l )/T 2 as a function of m 2 l /T 2 should not have a non-zero intercept. The quantity c(m l )/T 2 extracted from the eigenvalue densities for T T c and different light quark masses are shown in Fig. 4. Clearly it has a constant intercept c 0 = 0.82(17) which survives when chiral extrapolation is performed and even dominates over the usual O(m 2 l /T 2 ) term. The contribution of such a term in the eigenvalue density ρ(λ) = c 0 λT 2 to the chiral condensate, ψ ψ = dλ 2m l ρ(λ) λ 2 +m 2 l goes as 2m l T 2 c 0 ln(Λ/m l ), Λ being the ultraviolet cut-off to the eigenvalues, clearly vanishes in the chiral (m l → 0) limit. To summarize, this m l independent part of the bulk eigenvalue density does not break chiral symmetry but instead may contribute to U A (1) breaking. From this fit analysis it also is evident that the bulk eigenvalue density, to the leading order is O(λ) rather than O(λ 3 ) just above T c , even in the chiral limit. This is in addition to the contribution to the U A (1) breaking that comes due to the peak of small eigenvalues in the infrared. Moreover since this peak does not disappear and we do not observe any gap opening up in the infrared part of the eigenspectrum, we can conclude that U A (1) remains broken as we approach the chiral limit. In the next section we will provide a more quantitative estimate of the U A (1) breaking towards the chiral limit. Having observed little sensitivity of the exponent γ of the bulk eigenvalue density to the sea-quark mass, it is interesting to compare the spectra at different quark masses and also with the earlier results obtained with overlap fermions on HISQ configurations [23] for heavier than physical quark masses. Since the eigenvalue density is not a renormalization group invariant quantity, one has to renormalize the eigenvalue spectra for such a comparison [9,60,61]. The renormalized eigenvalue density is defined by scaling it with the valence strange quark mass. In order to do so one has to precisely estimate the valence quark masses. There is another advantage in precisely measuring the valence quark masses. The physics of the sea quarks can then be equivalently represented only in terms of the valence quarks or more specifically, in terms of the eigenvalues of the valence overlap Dirac matrix containing the exactly tuned valence quark masses. The tuning of the valence and sea quark masses can be numerically quite challenging. We have proposed to construct renormalized observables in terms of the valence and sea quarks respectively and match them in order to extract the valence quark masses, for given sea quark masses [23]. This is allowed because the renormalized observables describe the same physics. In this work, we measure the valence overlap strange quark mass by matching the renormalized quantity ∆,for the valence and the sea-quarks which is defined as, The chiral condensates ψ ψ (m s,l ) appearing in this observable, are first calculated for the valence overlap , where D m = D ov (1 − am/2M ) + am is the overlap Dirac operator with a (valence) quark mass m, the condensates can be calculated from the overlap eigenvalues, 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 of this expression which arises due to the zero modes, since it is a finite volume artefact to the above observable. We have checked explicitly, that in all our gauge ensembles this term due to the zero modes provides negligibly small correction to ∆, implying that the finite volume effects are under control in our tuning procedure. We then measure the ∆ for the HISQ sea by the exact calculation of the trace of the inverse of the HISQ operator using stochastic sources on the HISQ ensembles. Finally we compare the values of ∆ obtained for the valence overlap and the HISQ sea, and further tune the valence m s keeping the ratio m l /m s fixed for both the valence and sea quarks, until a perfect match of these two values of ∆ is obtained. This gives a matching valence overlap quark mass.
Once the valence m s is tuned to its sea value, we can equivalently describe the physics of the underlying sea quarks using only the valence overlap fermions with the tuned valence quark masses. The results for the tuned strange quark masses for different ensembles are tabulated in Table III. Tab. III: The valence strange quark masses obtained by matching the observable ∆ measured using valence overlap eigenvalues to that measured by inversion of the sea HISQ Dirac operator.
Subsequently, the comparison of the renormalized eigenvalue density m s ρ(λ)/T 4 as a function of λ/m s is shown in Fig. 5. While calculating the renormalized density, the total number of bins in λ/m s for different ensembles were kept fixed. We observe that the nearzero peak of the renormalized spectrum show very little sensitivity to the change in the quark mass. The bulk modes again show a linear rise. As mentioned earlier the spectrum for light quark mass m l = m s /80 has increasing density of these small eigenmodes hence the spectrum is shown only till λ/m s ∼ 2.5 with the first 200 eigenvalues we have measured.
With the tuned m s we next proceed to estimate whether the U A (1) is effectively restored above the crossover transition and its sensitivity to the light quark mass. However since U A (1) is not a symmetry there is no unique observable that is sensitive to its restoration. One can construct meson two-point correlators integrated over the spacetime volume and look for degeneracy between specific quantum number channels. For QCD with two light quark flavors, the difference of the integrated two-point correlators of isospin-triplet pion (iψτ 2 γ 5 ψ) and delta (ψτ 2 ψ) mesons, ω = χ π − χ δ is one such observable that was proposed as a measure of U A (1) breaking [62]. This is because (π, δ) transform as a doublet under U A (1) hence their correlators should be degenerate when U A (1) is restored. In fact one needs to further look at the degeneracy between higher point correlation functions for different meson quantum num- ber channels [38,63,64]. As a first test, we focus on this specific two-point correlation function. Due to chiral Ward identities χ π = ψ ψ (m l )/m l hence the expression for light quark chiral condensate can be used to measure it. Similarly χ δ is just the connected part of the scalar correlator defined as T V ∂ ∂m tr(D −1 m ∂Dm ∂m ) . The trace and inverses can be expressed in terms of the eigenvalues of the valence overlap Dirac matrix. The observable ω in terms of the eigenvalues of the overlap Dirac matrix is defined as, We have measured this quantity in terms of the first O(100) eigenvalues of the overlap Dirac operator at the tuned values of the valence quark masses. Chiral Ward identities ensure that χ π − χ δ = χ disc , where χ disc is the disconnected part of the integrated iso-singlet scalar meson correlator [66]. In order to verify the quality of our quark mass tuning we check whether this Ward identity is satisfied. We use the previously measured data for χ disc for physical quark masses from Ref. [33], obtained by the inversion of the HISQ Dirac operator using stochastic sources on N τ = 8, 12, 16 lattices and perform a continuum extrapolation of this observable. We then compare our continuum estimate of m 2 l χ disc /T 4 , to the observable m 2 l ω/T 4 that we have calculated using the eigenvalues of the valence overlap Dirac operator on the same HISQ ensembles using the tuned valence quark masses. The results for this comparative study are shown in Fig. 6. A reasonably good agreement of these renormalized quantities are observed giving us further confidence on our quark mass tuning procedure. Cont. m l 2 χ disc /T 4 Fig. 6: Comparison of χπ − χ δ measured using the overlap eigenvalues to the continuum estimates of χ disc using data from Ref. [33], shown for physical quark masses.
Next we study the quark mass dependence of the appropriately renormalized U A (1) breaking observable m 2 l ω/T 4 . It is important to note here that we have calculated only the first O(100) of the total millions of eigenvalues of the QCD Dirac operator, nonetheless these infrared eigenvalues contribute significantly to the U A (1) breaking. If ω ∼ m 2 l , i.e., as expected for a free quark gas and also in the perturbative regime then U A (1) is restored in the chiral limit. On the other hand if the leading order behavior of ω ∼ O(m 0 l ), then U A (1) is broken and its effective magnitude can be estimated. When chiral symmetry is effectively restored, we can use the following fit ansatz for our data on m 2 l ω/T 4 , corresponding to these two scenarios, where the former denotes U A (1) breaking whereas the latter is valid on its effective restoration. We have calculated the m 2 l ω/T 4 at 1.05 T c for three different tuned light valence quark masses, results of which are shown in Fig. 7. The data fits quite well to the first fit ansatz from Eq. 6, shown as a red band in the same figure, with the largest contribution to the uncertainty coming from the value corresponding to the lowest quark mass. Our data disfavors the second ansatz in Eq. 6 since the corresponding χ 2 per degree of freedom of the fit is about 5, which is almost a factor 2.5 larger than that corresponding to the first ansatz. The magnitude of (χ π − χ δ )/T 2 in the chiral limit is a 1 = 156 ± 13, which is clearly finite and non-zero within the current uncertainties. Thus we conclude that for 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 constant physical value of m s .
Noting again that in the chiral symmetry restored phase the topological susceptibility in QCD is related to m 2 l ω [9,65,66], it is evident that the former observable does not vanish when approaching the chiral limit (for finite values of m l ) along the line of constant m s . In contrast, for the two-flavor case and comparatively smaller volume lattices 2.4 fm 3 , the topological susceptibility is observed to vanish at a critical value of the light quark mass, which is lower but close to the physical quark mass [36]. It needs to be checked whether this observation in survives in the larger and finer lattices.

V. CONCLUSIONS AND OUTLOOK
In this work we report on the eigenvalue spectrum and the fate of anomalous U A (1) symmetry in the chiral symmetry restored phase of QCD on large volume N τ = 8 lattices as we approach the chiral limit along the line of constant physical strange quark mass. In order to correctly measure the number and density of the nearzero eigenmodes of the QCD ensembles using the Highly Improved Staggered quark (HISQ) discretization on the lattice, we use the overlap Dirac matrix as the valence or probe operator. This was done as the HISQ operator does not realize all the continuum flavor and anomalous symmetries on a finite lattice. In order to obtain physical and renormalized results even with different valence and sea quark actions, we have reported a fairly easy procedure to tune the valence quark masses to the sea quark mass using the eigenvalue spectrum as the input.
Comparing the appropriately renormalized eigenvalue density of the QCD Dirac operator, we observe that the bulk eigenvalue density at 1.05 T c can be represented as ρ(λ) ∼ λ in the chiral limit. This is unlike the expectations in the chiral limit [38] where the leading order behavior of the eigenvalue spectrum of QCD was derived to be ρ(λ) ∼ λ 3 , based on a part of chiral Ward identities and an assumption that the eigenvalue density is analytic in λ. This results in a non-zero value of the renormalized observable m 2 l (χ π − χ δ )/T 4 , which leads us to conclude that U A (1) is broken when one approaches the chiral limit along the line of constant physical m s .
For a final conclusive evidence, we need to re-visit this study at several choices of lattice spacing and perform a continuum extrapolation of our observables, which is computationally expensive and would require several years of dedicated efforts. Furthermore as one approaches the chiral limit, the lattice volumes have to be chosen large enough so that the spatial extent is few times larger than the corresponding pion Compton wavelength. In the present study we have chosen the lattice volumes keeping this criterion in mind and the m π L ∼ 3.5 for the lightest quark mass ensembles. However, there are already quite a few remarkable implications of our preliminary results. First of all, this is one of the first studies investigating the fate U A (1) anomalous symmetry, just above the chiral crossover for light quark masses as low as ∼ 0.6 MeV. We find that the qualitative features of the (renormalized) QCD Dirac eigenvalue spectrum is similar to the one for the physical values of the quark masses. There are no discontinuities in the infrared part of the eigenvalue spectrum. Secondly, when approaching the chiral limit along the line of constant physical strange quark mass, our study would suggest that one would eventually encounter the O(4) second order line of phase transitions.

VI. ACKNOWLEDGMENTS
We would like to dedicate this work in the memory of Prof. Edwin Laermann, who initiated this research project shortly before he passed away in August 2018. The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 'Stronginteraction matter under extreme conditions'-project number 315477589 -TRR 211. S.S. acknowledges support by the Department of Science and Technology, Govt. of India through a Ramanujan fellowship. The numerical computations in this work were performed on the GPU cluster at Bielefeld University and on Piz Daint at CSCS. We acknowledge PRACE for awarding us access to Piz Daint at CSCS, Switzerland. Our GPU code is in part based on some of the publicly available QUDA libraries [67]. We thank the HotQCD collaboration, for sharing the gauge configurations with us. S.S. is grateful to Prof. Frithjof Karsch for many helpful discussions related to this work. All data from our calculations, presented in the figures of this paper, can be found in https://doi.org/10.4119/unibi/2959004. Appendix A: Near-zero modes and Ginsparg-Wilson relation In this section we provide a detailed analysis of the violation of the Ginsparg-Wilson relation due to the numerical implementation of the overlap Dirac operator and whether it leads to the appearance of spurious near-zero modes in its eigenvalue spectrum. We will discuss here the results of our study for different sea-quark masses at T = 1.05 T c , for which we have measured the renormalized Dirac eigenvalue spectra shown in Fig. 5. It is evident from the eigenvalue spectrum, that we observe a peak comprising of small eigenvalues which seem to survive even when the quark masses are successively reduced towards the chiral limit. For the lattice volumes we have studied, the small eigenvalues did not appear sporadically for a few gauge configurations, rather all of them contributed to near-zero peak of eigenvalues. Thus for a meaningful analysis we instead chose to count the number of small eigenvalues i.e., λ < λ 0 appearing for each configuration and correlate it with the magnitude of Ginsparg-Wilson relation violation suffered by the overlap Dirac operator implemented on the corresponding gauge configurations. For m s /m l = 27, 40, we counted all overlap Dirac eigenvalues for each configuration that are lower than λ 0 = 0.4T . For m s /m l = 80 ensembles, we chose λ 0 = 0.1T instead, since these have a more dense eigenspectrum in the infrared. The magnitude of Ginsparg-Wilson relation violation due to the numerical implementation of the overlap Dirac matrix on the same configurations were simultaneously measured and compared to the probability of occurrence of small eigenvalues. The results of this analysis is shown in Fig. 8. For the m s /m l = 80 ensem-bles, we have divided the number of small eigenvalues appearing for each configuration by a factor of two to fit in the scale of the Figure. We do not find any obvious correlation between violation of the Ginsparg-Wilson relation and the proliferation of the number of small eigenvalues of the overlap Dirac operator. In fact, it is evident from Fig. 8 that the overlap Dirac matrix which had a larger count of small eigenvalues for particular HISQ configurations, satisfied the Ginsparg-Wilson relation to a relatively higher precision. We have also performed a similar comparison to check if there is any correlation between the appearance of many small eigenvalues to the numerical precision of the sign function appearing in the overlap Dirac matrix defined in Eq. 1. The result of this analysis is shown in Fig. 9. Here again we do not observe any obvious such correlation. In fact, the overlap Dirac matrix has a higher count of small eigenvalues on HISQ gauge configurations on which the sign function is implemented to a comparatively better numerical precision. Both these analyses again reconfirm the fact that the appearance of very small eigenvalues of the overlap matrix is not related to its numerical imprecision in this mixed-action set-up. This gives us a quite robust check that the small QCD Dirac eigenvalues that contribute to the U A (1) breaking above T c are not unphysical lattice artifacts.