Eigenvalues of the QCD Dirac matrix with improved staggered quarks in the continuum limit

We calculate the eigenmodes of the Highly Improved Staggered Quark (HISQ) matrix near the chiral crossover transition in QCD with $2+1$ flavors with the aim to gain more insights into its temperature dependence. On performing the continuum extrapolation, we do not observe any gap opening up in the infrared part of the eigenvalue density of the QCD Dirac operator; instead we observe a peak. The existence of the peak and oscillations of the infrared eigenmodes can be understood in terms of an interacting ensemble of instantons. From the properties of the continuum extrapolated eigenspectrum we further show that the anomalous $U_A(1)$ part of the chiral symmetry is not effectively restored simultaneously along with its non-singlet counterpart. We provide an explanation for this observation, further showing interesting connections between the anomalous $U_A(1)$ restoration and the change in the infrared part of the eigenvalue distribution.


I. INTRODUCTION
The eigenvalue spectrum of the quark Dirac operator contains valuable information about the fundamental properties of Quantum Chromodynamics (QCD).The chiral condensate which acts as an (pseudo) order parameter for the chiral (crossover) transition in QCD is related to the density of near-zero eigenvalues [1].In fact it was shown from very general considerations that the formation of the chiral condensate is related to the occurrence of small eigenvalues that scale proportional to the volume [2].The breaking of the non-singlet part of chiral symmetry, i.e., SU A (2) × SU V (2) → SU V (2) of QCD with physical quark masses at the crossover temperature T c = 156.5 ± 1.5 MeV [3] can also be explained in terms of modifications in the deep infrared part of the eigenvalue density.The flavor-singlet U A (1) part of the chiral symmetry on the other hand, is anomalous yet is believed to play an important role in determining the nature of the chiral phase transition [4][5][6].The temperature dependence of the amount of U A (1) breaking near the chiral crossover transition in QCD can be only determined using non-perturbative lattice techniques and is a topic of contemporary interest in lattice QCD; see, for e.g., Ref. [7,8] for recent reviews.Whereas there are some very compelling evidence that show U A (1) remains effectively broken in 2+1 flavor QCD with physical quark mass m [9][10][11][12][13][14][15], even when m → 0 [16], there are lattice studies which also favor an effective restoration at T c [17][18][19][20][21][22].
The eigenvalue spectrum of the QCD Dirac matrix also encodes within it some remarkable universal properties.It was shown that the route toward achieving the thermodynamic limit for the infrared modes of the Dirac op- * Electronic address: rshanker@imsc.res.inerator is universal [23], for any number of light quark flavors.The existence of a non-zero chiral condensate leads to a sum rule involving the sum of inverse squares of these small eigenvalues [2].These sum rules are universal irrespective of the details of the nature and type of gauge interactions [23,24] and could be derived from chiral random matrix theory [25].A good agreement was demonstrated for the distribution of the small eigenvalues and the spectral density of lattice QCD Dirac operator and chiral random matrix theory at zero temperature on small lattice volumes [26].In fact universal correlations between higher order spectral functions in a random matrix theory has been derived [27], and its connection to QCD was discussed.At finite temperature the universal features of infrared eigenvalues can be also accounted for within a random matrix theory [28][29][30].Additionally the infrared eigenvalue spectrum of QCD has more subtle features.A near-zero peak of localized eigenvalues has been observed for finite lattices, mixing with but very different from the delocalized bulk modes whose spectral density follows random matrix statistics [7,31].Whether or not such a feature survives in the continuum limit is yet to be ascertained.Previous studies of quark Dirac spectrum in an instanton liquid ensemble [29,32] at zero temperature have observed a similar peak-like feature.
With increasing temperature the localized modes start separating out from the random bulk modes leading to the opening up of a mobility edge [31].The corresponding temperature where a finite mobility edge separates the bulk modes from the localized one was initially estimated from lattice studies to be identical to T c in dynamical [33][34][35][36][37][38][39][40][41] as well as in quenched QCD [42,43], reminiscent of an Anderson-like transition that is observed in disordered semi-metals [44].However independent lattice studies do discuss another possible scenario where the opening of a finite mobility edge may occur at temperatures higher than T c [45], with an intermediate phase consisting of scale-invariant infinitely extended infrared modes [46,47] strongly interacting with the bulk modes leading to a singularity at the mobility edge.
Most of the previous lattice QCD studies were either performed in the quenched limit or with dynamical quarks but away from the physical point and for finite lattice spacings.On a finite lattice, the most often used lattice discretization, i.e., the staggered fermions only has a remnant of the continuum chiral symmetry group due to mixing of spin and flavor degrees of freedom.Furthermore the anomalous part of the chiral symmetry in the continuum is not realized exactly by the staggered/Wilson quarks and is expected to be recovered only in the continuum limit.We, for the first time study the properties of the eigenvalue spectrum of (highly) improved dynamical staggered Dirac operator in large volume lattices by carefully performing a continuum extrapolation.We show that the deep infrared spectrum of the QCD Dirac operator has indeed a peak of near-zero modes which survives in continuum.These are distinct from other infrared modes which have a linearly rising density and a quadratic level repulsion similar to a certain class of random matrix theories.These so-called bulk modes are delocalized in volume as compared with the near-zero modes, and they tend to distinctly disentangle from each other at a temperature ∼ 1.15 T c , which is also where U A (1) is effectively restored.In the subsequent sections we discuss our results and also provide a unified physical explanation of these phenomena we observe.

II. NUMERICAL DETAILS
In this work we use the gauge configurations for 2 + 1 flavor QCD with physical quark masses generated by the HotQCD collaboration using Highly Improved Staggered quark (HISQ) discretization for the fermions and treelevel Symanzik improved gauge action.These ensembles have been previously used to measure the equation of state of QCD both at zero and finite baryon density [3,48].The Goldstone pion mass is set to 140 MeV, and the kaon mass is 435 MeV for these configurations.We focus on five different temperatures, one below T c and others above T c .For most of these temperatures we consider three different lattice spacings corresponding to N τ = 8, 12, 16, the details of which are mentioned in Table I.The number of spatial lattice sites was chosen to be N s = 4N τ such that the spatial volume in each case was about 4 fm, which ensures that the system is close to the thermodynamic limit.We calculated the first 60, 100, 200 eigenvalues of the massless HISQ Dirac matrix for N τ = 16, 12, 8 respectively on these gauge ensembles using conjugate gradient method based algorithms.We have fixed the bin size λa = 0.001 for each N τ for measuring the eigenvalue density and performed a jack-knife analysis to remove any auto-correlation effects among the data in the bins.

III. RESULTS
A. General features of the eigenvalue spectrum of QCD using HISQ Dirac operator in continuum limit In this section we study in detail the eigenvalue density ρ(λ) of the quark Dirac operator in 2 + 1 flavor QCD by performing a continuum extrapolation of the parameters characterizing the eigenspectrum calculated on the lattice with HISQ discretization.
At zero temperature it is known from chiral perturbation theory [49] that the linearly rising part of the eigenvalue density, due to the so-called bulk modes, is expressed as The intercept of bulk eigenvalue density gives the chiral condensate.The ratio of the slope and the intercept of the density as a function of λ should be proportional to the chiral condensate.We first focus on the intercept and the slope (linear in λ) of the eigenvalue density at the lowest temperature T = 145 MeV, shown in the top left panel of Fig. 1, and compare with the expectations from Eq. 1.At this temperature we could only obtain a continuum estimate of the slope and intercept as we have data for two lattice spacings.From the continuum estimate of the intercept we obtain a chiral condensate ⟨0| ψψ|0⟩/T 3 = 18.4 using Eq. 1. From the slope we could similarly extract the square ⟨0| ψψ|0⟩ 2 and by substituting N f = 3, F π = 94.14MeV the chiral condensate (normalized by T 3 ) to be 17.3 which is consistent with the one extracted from the intercept.This demonstrates the consistency of our fit procedure.The value obtained here from the eigenvalue spectrum is also consistent with the value of ⟨0| ψψ|0⟩/T 3 = 18.8 obtained from the inversion of the HISQ Dirac operator on stochastic noise vectors and performing a continuum estimate using the N τ = 12, 16 data on a much larger set of HotQCD configurations [50].Thus we conclude here that the leading features of the eigenvalue density of QCD at 145 MeV are indeed very well represented within chiral perturbation theory.The bulk eigenvalue density in the chirally symmetric phase has been studied earlier in Ref. [51].Most generally, it can be expressed as a function of λ as Here c 1 is the coefficient that characterizes the leadingorder growth of the eigenvalue spectrum in the deep infrared and c 2 is its next-to leading order coefficient which eventually has a λ 3 dependence predicted from perturbation theory.The intercept ρ 0 gives the the chiral condensate.The coefficients c 1,2,3 can in general be a function of the temperature T and the light-quark mass m.
The results of the eigenvalue density ρ(λ)/T 3 as a function of λ for T > T c are shown in Fig. 1.On the finest available N τ = 16 lattice, we observe two distinct features in the eigenvalue spectrum, a peak of near-zero eigenvalues and the linearly rising part, due to the bulk modes, as previously mentioned.For T ≲ T c , the nearzero and the bulk eigenvalues overlap strongly making it impossible to distinguish them apart.At higher temperatures, the bulk eigenvalues separate out from the deep-infrared part of the spectrum allowing for near-zero modes to be distinctly visible.Comparing the results of different lattice spacings, we observe the same trend at each temperature above T c , i.e., near-zero peak gets smeared with the bulk for coarser lattices and becomes more prominent in the continuum limit.This is thus a physical feature of the eigenspectrum and not a lattice artifact.In order to interpret its origin we recall that in the chiral random matrix theory (cRMT) at zero temperature, the scaled eigenvalue (cλ) density of the Dirac operator for N f = 2 flavors and zero topological charge sector is distributed according to [52], To compare our data with the above formula, we take c = V ⟨0| ψψ|0⟩/T , where V is the spatial volume of the system, and the value of ⟨ ψψ⟩ at finite temperature is obtained from Ref. [50] which uses the same HotQCD gauge configurations, a subset of which is used in this work.Further we also scale the eigenvalues such that the first moment of probability distribution for lowest eigenvalues for the data matches with the first moment of Eq. 4. A comparison of near zero modes for four different temperatures, T = 145, 162, 166, 171 MeV, is shown in Fig. 2. We observe a good agreement with cRMT for T = 145 MeV, in particular, the initial few oscillations of the small eigenvalue density as a function of cλ.
Incidentally an agreement between Eq. 3 and the eigenvalue density from the instanton liquid model (ILM) at T = 0 was observed in [53].However at finite temperature the oscillations in the eigenvalue density within ILM are smeared out over a length scale ∼ 1/T [29] which is qualitatively similar to what we observe for T > T c in Fig. 2. Now focusing on the bulk modes, it was shown using chiral Ward identities that in the symmetry-restored phase, the sufficient condition for U A (1) restoration evident from the degeneracy of up to six-point correlation functions in the scalar-pseudo-scalar sector are c 1 = O(m 2 ) + ... and c 3 = c 30 + O(m 2 ) + ....The perturbative λ 3 growth in Eq. 2 can have a mass-independent coefficient which however does not lead to U A (1) breaking.We verify whether indeed it is true even non-perturbatively by performing a fit to the bulk part i.e. all eigenvalues λ > λ 0 with ρ(λ) + ρ0 T 3 .This ansatz neglects higher powers in λ which is well justified since we are in the deep infrared of the eigenspectrum, represented by O(100) eigenvalues out of a total million available on such lattice sizes.The results of the fit are discussed in Table II.The extracted slope c 1 for each temperature T > T c , at three different values of N τ then allows us to perform a continuum (∼ 1/N 2 τ ) extrapolation of this coefficient.We next study the m dependence of this continuum extrapolated coefficient c 1 (m, T ).The results of the fits are shown in Fig. 3.It is evident from the fit that it is more favorable that c 1 is proportional to T 2 (χ 2 /d.o.f=0.6) to leading order rather than c 1 is proportional to m 2 (χ 2 /d.o.f=0.1).From the fit we obtain the value of c 1 (m, T )/T 2 = 16.8(4).
The finite result for the slope in the continuum limit, i.e., the m-independent term in c 1 ensures that the U A (1) part of the chiral symmetry will remain minimally broken in the chiral limit in the symmetry-restored phase as the maximum contribution to U A (1) breaking comes from the near-zero eigenvalues, which we observe in the next section.Moreover the slope of the eigendensity for T ≲ 1.12 T c is distinctly different from the perturbative λ 3 rise implying significant non-perturbative effects.Tab.II: Lattice sizes (N 3 σ × Nτ ), temperatures (T ), the estimated values of c1/T 2 and ρ0/T 3 after the fit to the bulk modes which are defined beyond the lower cutoff at λ0/T .

B. The fate of UA(1) breaking in the continuum limit
Since the flavor singlet part of the chiral symmetry is anomalous it has no corresponding order parameter.Hence to measure whether this singlet part of the chiral symmetry is simultaneously (and effectively) restored along with the non-singlet part, it has been suggested [54] to look at the degeneracies of the integrated correlators of mesons i.e., χ π − χ δ .In the continuum, the integrated meson correlators are related to each other through the following relations, χ δ = χ σ − 4χ disc and χ π = χ η + 4χ 5disc .These integrated meson correlators are defined as χ π = d 4 x ⟨π i (x)π i (0)⟩, χ σ = d 4 x ⟨σ(x)σ(0)⟩, χ δ = d 4 x ⟨δ i (x)δ i (0)⟩ and χ η = d 4 x ⟨η(x)η(0)⟩ where i = 1, 2, 3. We measure (χ π −χ δ )/T 2 at the four different temperatures above T c , and perform a ∼ 1/N 2 τ continuum extrapolation at each temperature, the results of which are shown in Fig. 4. For the highest temperature, we have only two data points available corresponding to N τ = 8, 12 for performing the continuum extrapolation.We hence assigned 40% and 20% error to the values for the slope and the intercept respectively, similar to that obtained from a fit to the T = 171 MeV data.This observable receives 99% contribution from the near-zero eigenvalues for N τ = 16.Performing continuum estimates with only two data points corresponding to finer lattice sizes N τ = 16, 12 at each temperature, gives a higher intercept than the corresponding extrapolation considering all three N τ values.Hence the finiteness of this observable is quite robust and we conclude that U A (1) does not get effectively restored at T c .Motivated from the fact that the major contribution to (χ π −χ δ )/T 2 comes from the near-zero modes, we expect a 1/T 2 dependence to this quantity if the corresponding eigenvalue density can be characterized as a well-defined peak.Furthermore from the chiral perturbation theory at finite temperature one expects a similar 1/T 2 dependence [55] near T ≳ T c .We thus fit the continuum extrapolated values of (χ π − χ δ )/T 2 at each temperature T > T c , i.e., the intercept of the fits shown in Fig. 4 to the ansatz a + b/T 2 .After performing the fit, shown in Fig. 5 we extract the temperature T /T c = 1.147 (25) beyond which this U A (1)-breaking observable drops to zero.
Earlier analytic works based on the properties of integrated two-point meson correlators argued that the U A (1) breaking comes from the eigenmodes of the Dirac operator at λ = 0 [56] or close to zero [57].Whereas exact zero modes do not contribute in the thermodynamic limit, the density of near-zero modes at T > T c was expected to be zero [57] in the chiral symmetry-restored phase of QCD, whereas we observe a robust presence of the near-zero modes in the continuum limit which dominantly contributes to U A (1) breaking for T ≲ 1.15 T c .
We next compare our result with the earlier observation of U A (1) (effective) restoration temperature of ∼ 200 MeV obtained from the continuum extrapolated results for the integrated screening correlators (χ π − χ δ )/T 2 in 2 + 1 flavor QCD using HISQ discretization [14] with heavier than physical light quark mass, corresponding to m π = 160 MeV.The corresponding pseudo-critical temperature is also higher by ∼ 4 MeV; hence the restoration temperature comes out to be around 1.2 T c which agrees with us within the error bar.Furthermore a recent work [15] has reported breaking of U A (1) due to a m 2 δ(λ) feature in the eigenvalue spectrum at about 1.5 T c which survives even in the chiral limit.We note here that U A (1) breaking due to this specific feature of the Dirac spectrum is expected to survive even at asymptotically high temperatures where the QCD vacuum can be explained in terms of a dilute gas of instantons [12,13].We however discuss here how a non-trivial breaking of the U A (1) part of chiral symmetry can arise due to features in the infrared part of the eigenvalue spectrum unlike the dilute instanton gas regime and show that beyond 1.15 T c such a contribution gets effectively restored which then may transition into the temperature regime studied in Ref. [15].We next verify that the chiral Ward identities are satisfied by the HISQ action.

C. Verifying the chiral Ward identities
In the chiral symmetry restored phase, χ σ = χ π and χ δ = χ η hence one obtaines χ π − χ δ = 4χ 5,disc .Using chiral Ward identities it is known that χ 5,disc = χ t /m 2 where χ t is the topological susceptibility of QCD.This allows one to relate the U A (1) breaking parameter to the topological susceptibility through the relation 1/4(χ π − χ δ )m 2 /T 4 = χ t /T 4 .A comparison of these two observables is shown in Fig. 6.From the figure it is evident that for T > 1.05 T c , when chiral symmetry is effectively restored, the two quantities agree with each other within errors.This is particularly interesting since for staggered quarks, even though the chiral and taste symmetries are intermixed at finite lattice spacing, the symmetries of QCD and related chiral Ward identities are recovered in the continuum limit.Fig. 6: A comparison of the integrated renormalized correlator (χπ −χ δ )m 2 /4T 4 with the topological susceptibility (measured independently using gradient flow in Ref. [58]) for temperatures > Tc.

D. Distribution of the lowest eigenvalue at finite temperature
The probability distribution of the lowest eigenvalue of the QCD Dirac operator λ min has inherent information about the microscopic degrees of freedom.For a chiral random matrix ensemble for N f = 2 (at zero temperature) the lowest eigenvalue is distributed according to [52] At the lowest temperature T = 145 MeV, we calculate the probability distribution of the scaled lowest eigenvalue cλ min at different lattice spacings and perform a continuum estimate of the distributions, for which we have extracted the lowest eigenvalue from each configuration for N τ = 12, 16 for T = 145 MeV and N τ = 8, 12, 16 for T = 171 MeV and later rescaled to the dimensionless quantity cλ min , where the value of ⟨ ψψ⟩ at finite temperature is obtained from Ref. [50].Keeping the bin size constant we obtained the probability distribution of cλ min for each N τ and then performed a spline interpolation by taking appropriate weights proportional to the errors for each data point in order to have a smoother interpolating curve.Next we performed a continuum extrapolation at each value of cλ min of the interpolating function with the ansatz c + d/N 2 τ .We assigned a 15% error for T = 145 MeV, as we only had two points while performing the continuum extrapolation.In order to compare the probability distribution of the lowest eigenvalue for both the temperatures with Eq. 4 we have to match their first moments with the cRMT distribution.We have calculated the first moment for Eq. 4 and found the result to be 4.344.Next we have scaled the cλ min and the probability distribution for the lowest eigenvalue obtained at both these temperatures from our calculations in 2 + 1 flavor QCD such that the first moment is exactly 4.344 and the area under the curve is unity.The probability distribution of the scaled lowest eigenvalue is shown in Fig. 7.The continuum extrapolation of the probability distribution at T = 145 MeV shown as the orange band is close to the probability distribution of a N f = 2 chiral Gaussian unitary random matrix ensemble.In contrast, we also plot the probability distribution of the scaled lowest eigenvalue at T = 171 MeV whose continuum extrapolation is shown as a blue band in Fig. 7.It is evident that the lowest eigenvalue which is a part of the nearzero peak follows a very different statistics rather than known from a chiral RMT.

E. The level spacing distribution for bulk modes
In order to understand the properties of bulk modes we look at their level spacing distribution.To study the universal properties of the eigenvalue level spacing fluctuations one has to remove the system dependent mean.This is done by a method called unfolding.Let λ represent eigenvalues in the ascending sequence for any particular gauge configuration.The average density of the eigenvalues in the sequence i.e. the reciprocal of the average spacing as a function of λ is represented as ρ(λ).The eigenvalue sequence can then be unfolded using the average level-staircase function, η(λ) = λ λ0 dλ ′ ρ(λ ′ ) which tells us how many eigenvalues in this sequence are less than λ on an average.Here λ 0 labels the eigenvalue beyond which all the higher eigenvalues are bulk modes and below which are the near-zero modes.The unfolded sequence is labeled by λ uf i = η(λ i ), where the index i labels the original eigenvalue whose unfolding is performed.When appropriately normalized, the average spacing between the unfolded eigenvalues equals unity.The nearest neighbor spacing distribution is constructed by calculating the differences between consecutive unfolded eigenvalues λ uf i+1 −λ uf i and organizing them into histogram bins.This gives us a picture of how the eigenvalue spacings fluctuate about the average which we have shown in Fig. 8 for four different temperatures T = 162, 166, 171, 176 MeV and at each temperature, for the two different lattice sizes N τ = 8, 12.We then compare the nearest neighbor spacing distributions to the Wigner surmise for the Gaussian unitary ensemble (GUE) given by P (s) = 32/π 2 s 2 e −4s 2 /π shown as dotted lines in Fig. 8.It is evident that the level repulsion between the bulk modes for small s is quadratic similar to that of random matrices belonging to the GUE.We see however that as the temperature increases, the agreement to GUE with the N τ = 12 data for s > 1 is not so good, whereas the low s < 1 part agrees very well.This occurs due to the contamination of the bulk modes, which are more closely spaced than the near-zero modes which start to build up forming a peaklike structure in the infrared part of the eigenvalue spectrum.For the N τ = 16 lattices which has a clear well-defined peak of near-zero modes at T ≥ 166 MeV, the contamination with the bulk modes is expected to be even more severe.As expected the comparison of the tail of the spacing distribution of the N τ = 16 data to the GUE prediction produces not-so-good agreement.In order to account for the long tail of the spacing distribution we fit it to a distribution P (s) ∼ s 2 exp (−αs) which shows strong quadratic level repulsion at small values of s but falls off slowly at large values of s parametrized by a fit parameter α.After performing the fit of the level separation with this ansatz, we obtain the value of α = 3.02(7), 3.17 This is a generic feature of the eigenvalue spacing distribution of a strongly disordered system [59] whose bulk eigenmodes in the center of the band follows a similar behavior as RMTs except for the tails of the distribution due to contamination with the localized states.We will explain this feature in more detail in the next section.Having shown the distinct features of near-zero and bulk modes, the next question we ask is whether the near-zero modes which arise due to instanton interactions can distinctly disentangle out from the delocalized bulk modes.The QCD medium above T c consists of quarks interacting with each other as well residing in a disorder potential very similar to an interacting electron system in a background random potential studied in detail in Ref. [59].Such systems have a mixing between the localized and delocalized many-body states which is in contrast to the traditional Anderson model, consisting of non-interacting electrons in the presence of a random disorder.In the Anderson model, one-electron states with the same energy but with different localization properties cannot co-exist in three or more spatial dimensions, d ≥ 3.There exist bands of localized and extended states, and a unique energy separating two such bands for d ≥ 3 which is termed as a mobility edge.The QCD medium above T c however cannot be understood through a conventional Anderson model; it has far more interesting properties, like the possibility of the existence of a scale-invariant infrared phase above T c discussed in the recent literature [41,47].In fact we do observe a mixing between the localized states with the bulk spectrum of the QCD Dirac operator in the level spacing distribution data as shown in Fig. 9.
In order to estimate the temperature when the bulk modes separate out from the deep-infrared peak of eigenvalues, we first estimate the typical spread of the nearzero peak visible in the eigenvalue density plots corresponding to N τ = 16 in Fig. 1.Recall that we have already estimated the slope and the intercept of the bulk eigenvalue density, using which we subtract the bulk mode contribution from the total eigenvalue spectrum for the N τ = 16 data at each temperature.The nearzero peak which we get after subtracting the bulk modes has a typical spread which we estimate to be λ/T = 0.08 for all the temperatures above T c .Next, using the fact that the bulk modes have a linear-in-λ dependence with a slope c 1 /T 2 = 16.8(4) in the continuum and the nearzero and bulk modes will separate out at a particular temperature, leading to the density to drop to zero at λ/T = 0.08, we estimate the value of the bulk intercept ρ 0 /T 3 = −1.34(corresponding to λ = 0) in the continuum limit.We then take the values of the intercept ρ 0 /T 3 of bulk mode density for all T > T c from Table II and perform a continuum extrapolation with an ansatz The continuum values of the quantity ρ 0 /T 3 for T > T c , so-obtained after the fit are shown in Fig. 10.At the highest temperature T = 176 MeV, a 10% error is assigned to the data point since we could perform a continuum estimate, with data available only for two N τ values.Next, fitting the continuum extrapolated data for ρ 0 /T 3 as a function of temperature with a fit ansatz ρ 0 /T 3 = d 1 (T /T c )+d 2 we obtain the fit parameters to be d 1 = −23.1(3)and d 2 = 25.3(3)respectively.After obtaining this parametric dependence of the continuum estimates of the intercept as a function of temperature, we can now extract the temperature where the value of the intercept ρ 0 /T 3 = −1.34,i.e., when the near-zero modes distinctly emerge out from the bulk spectrum.The extracted temperature comes out to be T = 1.15(3)T c .This is within the temperature range when the U A (1) part of the chiral symmetry is effectively restored.

IV. WHY IS UA(1) EFFECTIVELY RESTORED
AT TEMPERATURES ABOVE Tc?
In order to interpret these results, one could visualize quarks as a many-body state moving in the background of an interacting ensemble of instantons, where the strength of the interactions changes as a function of temperature.At the microscopic level it is conjectured that the instantons remain strongly correlated below T c , subsequently transitioning to a liquid-like phase with a finite but weaker correlation length [60] just above T c , and eventually to a gas-like phase around 2 T c [13,15].Below T c the intercept of the infrared eigenvalue density quantifies the chiral condensate which corresponds to the breaking of the non-singlet part of the chiral symmetry.Owing to very strong correlations the microscopic details of the interactions are lost and the eigenvalues repel strongly similar to random matrices of a GUE.As the temperature is increased, at ∼ 171 MeV, the near-zero eigenvalues start to become prominent.These eventually separate from the bulk at ∼ 1.15 T c .Earlier studies have observed screening of inter-instanton interactions and build-up of local pockets of Polyakov loop fluctuations [39,61,62] above such temperatures.This is also the region where the constituent dyons of the closelyspaced instantons interact semi-classically and thus start to become detectable [63][64][65][66].
Incidentally this suppression of long range instanton interactions also weakens the strength of U A (1) breaking, allowing for its effective restoration at T ≳ 1.15 T c .Lattice studies [67,68] have reported a jump in the electrical conductivity around this temperature.Similar phenomena have also been reported in many-electron systems [59] in a disordered potential where the interplay between disorder and interactions causes a separation between the localized and delocalized states leading to a jump in the electrical conductivity from near-zero to a finite value.

V. CONCLUSIONS
In this work we have addressed a long-standing question of whether the flavor singlet U A (1) subgroup of the chiral symmetry gets effectively restored simultaneously with the non-singlet part for QCD with two light quark flavors at T c .The effective restoration of the anomalous U A (1) symmetry is a non-perturbative phenomenon driven by the deep infrared part of the QCD Dirac eigenvalue spectrum.By carefully performing the continuum extrapolation of the staggered Dirac spectrum on the lattice and studying in detail its properties, we explicitly demonstrate that U A (1) remains effectively broken in the chirally symmetric phase (T > T c ) for T ≲ 1.15 T c .We also provide arguments for why this conclusion should remain unchanged even in the chiral limit.
With the increase in temperature the strength of interactions between the instantons starts to weaken due to which the deep infrared peak of the spectrum is separated out from the bulk modes, which happens at around T ∼ 1.15 T c .The tunneling probability due to instantons also decreases with increasing temperature which results in lowering of the height of the near-zero peak of eigenvalue density.We show for the first time that both these phenomena are possibly the reason behind the U A (1) restoration, which also surprisingly happens to be around the same temperature.Observations of such rich interplay of phenomena in QCD matter above T c should be quite robust, since these are made after performing a continuum extrapolation.It will be interesting to observe further finer details of chiral transition in the massless limit with QCD Dirac operators which have exact chiral symmetry on the lattice.
All data from our calculations, presented in the figures of this paper, can be found in Ref. [69].

Fig. 5 :
Fig. 5: The continuum estimates for (χπ −χ δ )/T 2 for temperatures greater than Tc shown as points fitted to a functional form a + b/T 2 shown as a band.

Fig. 7 :
Fig. 7: The continuum extrapolated probability distribution of scaled lowest eigenvalue for T = 145, 171 MeV shown as orange and blue bands respectively and these are compared with the cRMT prediction for N f = 2.

Fig. 8 :
Fig. 8: Unfolded level spacing distributions of bulk eigenvalue modes for different temperatures shown as a function of different lattice spacings or equivalently, Nτ .The dotted curves in each plot correspond to the Wigner surmise for Gaussian unitary random matrix ensembles.

Fig. 9 :
Fig.9: A fit to the eigenvalue level spacing distribution using a mixed ansatz for Nτ = 16 at T = 162, 166, 171 MeV and compared with the prediction from a GUE of random matrices.

F
. Separating the near-zero from the bulk eigenvalues of the QCD Dirac spectrum

Fig. 10 :
Fig.10: Continuum extrapolation of the bulk intercept for eigenvalue densities at different temperatures above Tc.The horizontal line corresponds to ρ0/T 3 = −1.34 for the bulk spectrum when it completely separates from near-zero modes. ρ Fig.3: Continuum estimates for c1(m, T )/T 2 for T > Tc obtained after fitting the points with an m-independent constant (orange band) and a sum of quadratic (m 2 /T 2 ) and quartic (m 4 /T 4 ) dependent fit function (gray band).