Microscopic Encoding of Macroscopic Universality: Scaling Properties of Dirac Eigenspectra near QCD Chiral Phase Transition

Macroscopic properties of the strong interaction near its chiral phase transition exhibit scaling behaviors, which are the same as those observed close to the magnetic transition in a 3-dimensional classical spin system with $O(4)$ symmetry. We show that the universal scaling properties of the chiral phase transition in Quantum Chromodynamics (QCD) at the macroscale are, in fact, encoded within the microscopic energy levels of its fundamental constituents, the quarks. We establish a connection between the cumulants of the chiral order parameter, i.e., the chiral condensate, and the correlations among the energy levels of quarks, i.e., the eigenspectra of the massless QCD Dirac operator. This relation elucidates how the fluctuations of the chiral condensate arise from the correlations within the infrared part of the energy spectra of quarks, and naturally leads to a generalization of the Banks-Casher relation for the cumulants of the chiral condensate. Then, through (2+1)-flavor lattice QCD calculations with varying light quark masses near the QCD chiral transition, we demonstrate that the correlations among the infrared part of the Dirac eigenvalue spectra exhibit same universal scaling behaviors as expected of the cumulants of the chiral condensate. We find that these universal scaling behaviors extend up to the physical values of the up and down quark masses. Our study reveals how the hidden scaling features at the microscale give rise to the macroscopic universal properties of QCD.

Introduction.-Criticalphenomena exhibited in the vicinity of continuous second order phase transition are ubiquitous in nature [1].Phase transitions at high temperatures involving the electroweak and strong forces of nature have given rise to the Universe as we experience today.Lattice-regularized field theory calculations have shown that the high temperature transitions in matters governed by the electroweak (see, e.g., discussions in [2]) and strong interactions [3][4][5][6] are rapid crossovers.But for small enough Higgs mass [7][8][9][10][11][12] and in the massless (chiral) limit of up and down quarks [13] the electroweak and strong forces, respectively, are expected to undergo true phase transitions.
In the vicinity of a second order phase transition macroscopic quantities related to the order parameter exhibit telltale scaling behaviors that are uniquely characterized by the dimensionality and global symmetries of the system, irrespective of the details of its microscopic degrees of freedom and interactions.Based on the global symmetries of quantum chromodynamics (QCD), the theory of strong interaction, the second order QCD chiral transition can be in the 3-dimensional O(4) universality class [13] [14].For lattice-regularized QCD using costly chiral fermions, e.g., overlap fermions, the chiral symmetry of QCD is strictly preserved for any finite value of the regulator, i.e., for nonvanishing lattice spacing, the universality class is 3-dimensional O(4); while for the case using staggered fermions the chiral symmetry is only partially preserved and the universality class falls into 3-dimensional O(2) [15][16][17][18].With present computational resources it is only feasible to carry out large scale lattice QCD simulations toward a chiral limit of up and down quarks using the staggered fermions.Thus, lattice QCD studies of chiral phase transition using staggered fermion discretizations are expected to observe the same macroscopic scaling behaviors as that in the vicinity of the liquid to superfluid λ transition in 4 He [19]; notwithstanding, the microscopic degrees of freedom for QCD are quarks and gluons governed by the strong force while for 4 He are electrons and photons interacting via the electromagnetism.Because of this universal feature, to understand and predict macroscopic properties of a system close a second order phase transition one, most often, resorts to a simplified effective theory possessing the same dimensionality and global symmetries of the original theory, ignoring its microscopic complexities [19].
However, it is unclear if and how the macroscopic universal scaling properties of the strong interaction in the vicinity of the chiral phase transition are concealed within the microscopic energy spectrum of its elementary degrees of freedom, quarks.The goal of this work is lattice QCD-based understanding of possible connections between the universal features at the macroscopic and microscopic scales of QCD.An analogous goal for quantum electrodynamics will be to comprehend how the macroscopic scaling properties near the λ transition of 4 He arise from the energy levels of electrons without resorting to an effective theory.
In this Letter, first, we will establish theoretical relations between the cumulants of the chiral condensate and correlations among the eigenvalue spectrum of the massless Dirac operator.We will then demonstrate how the O(2) scaling properties of the cumulants of the chiral condensate are reflected within the correlations of the Dirac The nth order cumulants, K n , of the order parameter ψψ(m l ) can be obtained from the generating functional as Hereafter, T is the temperature and V is the spatial volume of the system, and ⟨•⟩ 0 denotes expectation value with respect to the QCD partition function in the chiral limit, Z(0 and Z(m l )/Z(0) = exp{−m l ψψ(m l )} 0 , it is easy to see that K n are the standard cumulants of ψψ(m l ); e.g., Energy levels of a massless quark in the background of U are given by the eigenvalues, In terms of λ j [U] the probe operator can be expressed as ψψ(ϵ) ≡ 2Tr( / D[U] + ϵ) −1 = 2 j (iλ j + ϵ) −1 .Thus, Eq. 1 becomes where From Eq. 2 it is straightforward to obtain where P 1 (λ) = K 1 [P U (λ; m l )] for n = 1, and for n ≥ 2 Here, K 1 is the first order joint cumulant of n variables (X i ) defined as Eq. 5 is our main theoretical result connecting the cumulants of the order parameter to the n-point correlations of the quark energy levels ρ U (λ).The explicit expressions of K 1 [ ψψ], K 2 [ ψψ], and K 3 [ ψψ] in terms of ρ U (λ) are provided in Eq. 14 of Supplemental Material.
The chiral phase transition in the staggered lattice QCD at nonvanishing lattice spacings is expected to be in the 3-dimensional O(2) universality class.Following the expectations for a 3-dimensional O(2) spin model near criticality [20], in the vicinity of the chiral transition The scaling variable z ∝ z 0 m −1/βδ l (T − T c )/T c , where T c is the chiral phase transition temperature and z 0 is a scale parameter; both are system specific.β and δ are the universal critical exponents, and f n+1 (z) = (1/δ−n+ 1)f n (z)−zf ′ n (z)/βδ are the universal scaling functions of 3-dimensional O(2) universality class.Here, n ≥ 1 and the superscript prime denotes derivative with respect to z [21].In our work we adopted β = 0.349, δ = 4.78 and for consistency the scaling functions f n (z) of the O(2) universality class determined from Refs.[16,20].
Eq. 8 indicates that universal scaling properties of the macroscopic observables K n [ ψψ] arise from the correlations among the microscopic energy levels P n (λ).To elucidate this point, consider m l → 0. Then from Eq. 4 one finds P U (λ; Noting that K 1 of n identical variables is equivalent to K n , in the chiral limit Eq. 5 thus becomes a generalization of the Banks-Casher relation [22] expressed as follows: To the best of our knowledge this generalized relation between the higher order cumulants of chiral condensate in the chiral limit and density of the deep infrared energies of quarks is new in literature.
In the chiral limit and close to T c , K n [ ψψ] should manifest universal scaling, e.g., . According to Eq. 9 this must arise from the universal behaviors of the λindependent K n [ρ U (0)].Thus, it is natural to expect that for small m l within the scaling window the critical scaling of K n [ ψψ] in Eq. 8 arises from the universal behaviors of the amplitudes of P n (λ) at the infrared, and not from its system-specific λ dependence; i.e., Here, g n (λ) are nonuniversal functions encoding the properties of the specific system under consideration.Next, we numerically establish Eq. 10 through lattice QCD calculations.
Lattice QCD calculations.-LatticeQCD calculations were carried out between T = 135 − 176 MeV for (2+1)flavor QCD using the highly improved staggered quarks and the tree-level Symanzik gauge action, a setup extensively used by the HotQCD Collaboration [23][24][25][26][27][28].m s was fixed to its physical value with a varying m l = m s /27, m s /40, m s /80, m s /160, which correspond to the Goldstone pion masses m π ≈ 140, 110, 80, 55 MeV, respectively.The temporal extents of the lattices were N τ = 8, and spatial extents were chosen to be N σ = (4 − 7)N τ .The gauge field configurations were generated using a software suite SIMULATeQCD [29], and the same gauge ensembles were used for the determination of chiral phase transition temperature in the continuum limit [30].
Observables were calculated on gauge configurations from every tenth molecular dynamics trajectory of unit length, after skipping at least first 800 trajectories for thermalization.ρ U (λ) and P n (λ) for n = 1, 2, 3 over the entire range of λ were computed using the Chebyshev filtering technique combined with the stochastic estimate method [31][32][33][34][35][36] on about 3000 configurations.Orders of the Chebyshev polynomials were chosen to be 2 × 10 5 and 24−96 Gaussian stochastic sources were used.Exact details are provided in Tab.I of Supplemental Material.
Results.-Owing to Eq. 9 we expect the relevant infrared energy scale is λ ∼ m l for small values of m l .It is natural to express all quantities as functions of the dimensionless and renormalization group invariant λ/m l : where the dimensionless and renormalization group invariant Kn [ ψψ] = m n s K n [ ψψ]/T 4 c .In Fig. 1 we show Pn ( λ) for n = 1, 2, 3 as a function of λ in the proximity of T c = 144.2(6)MeV [17].Pn ( λ) rapidly vanishes for λ ≳ 1, and the regions where Pn ( λ) ̸ = 0 get smaller with increasing n.This reinforces that the relevant infrared energy scale turns out to be λ ∼ 1.In this infrared region Pn ( λ) at a fixed T shows clear dependences on m l , which becomes stronger for increasing n.The form of m l dependence of Pn ( λ) also changes with varying T .
As shown in Fig. 3 of Supplemental Material we have checked that integrals over the relevant nonvanishing infrared regions of Pn ( λ) (cf.Eq. 11) reproduces K n [ ψψ], independently calculated through inversions of the fermion matrices, for n = 1, 2, 3.As seen from Fig. 1, expectedly, our results become increasingly noisy with increasing n and decreasing m l .With our present statistics we cannot access correlation functions n > 3, particularly for smaller m l .
The m l and T dependence of Pn ( λ) shown in Fig. 1 can be understood in terms of the 3-dimensional O(2) scaling properties.Once the Pn ( λ) are rescaled with respective m1/δ+1−n l f n (z) the data in Fig. 1 magically collapse onto each other, see Fig. 2. The system-specific parameters T c = 144.2(6)MeV and z 0 = 1.83(9) needed to obtain f n (z) were taken from Ref. [17], where 3-dimensional O(2) scaling fits were carried out for the same lattice ensembles but using an entirely different macroscopic ob- servable, namely the m l dependence of the static quark free energy.Thus, our expectations from Eq. 10 are clearly borne out in Fig. 2, namely where ĝn ( λ) characterize the system specific of the nth order energy-level correlations.To satisfy our generalized Banks-Casher relations of Eq. 9 the ĝn ( λ) must also satisfy lim V →∞ lim a→0 lim m l →0 ĝn ( λ) → δ( λ), such that K n [ ψψ] has the correct scaling behavior in (T − T c )/T c .The values of z 0 and T c used to demonstrate the universal scaling in Fig. 2 were obtained fitting lattice results for m l dependence of the static quark free energy only for 55 MeV ≤ m π ≤ 110 MeV [17].It is noteworthy that the physical QCD with m π ≈ 140 MeV also shows the same universal scaling for 135 MeV ≤ T ≤ 145 MeV.Outside this temperature window we do not observe scaling (see Fig. 5

of Supplemental Material).
As mentioned in Ref. [17], presently {T c , z 0 } are not very well determined.By using other values for {T c , z 0 } quoted in Ref. [17] we checked that the scaling of Fig. 2 is fairly insensitive to the exact values of {T c , z 0 } (see Fig. 6 of Supplemental Material).Presumably, this is because Pn ( λ) are sensitive only to the deep infrared physics λ ∼ m l .This is in contrast to many other macroscopic operators used for detailed scaling studies that contain large contributions from the ultraviolet energies [16,30,[37][38][39][40].This suggests that it even might be advantageous to use Pn ( λ) for detailed scaling studies to determine the system-specific parameters.Our focus here is to reveal the underlying connection between the universal features and the quark spectra, and such detailed scaling studies are beyond the scope of the present work.
Conclusions.-In this Letter, we investigate how the universal critical scaling of macroscopic observables near the QCD chiral transition arises from the microscopic degrees of freedom.We have presented a theoretical connection between the nth order cumulant of the chiral order parameter and the n-point correlations of the quark energy spectra.This connection led us to a generalized Banks-Casher relation, equating the nth order cumulant of chiral condensate to the nth order cumulant of the zero mode of the quark energy in the chiral limit.These new theoretical developments establish a direct connection between the universal scaling observed at the macroscale and the microscopic energy levels of the system.Through staggered lattice QCD calculations in the vicinity of the chiral phase transition with a series of light quark masses we have discovered the hidden universality within the correlations among the quark energy spectra.We have found that these universal behaviors are also imprinted within the microscopic energy levels of QCD with physical light quark masses.
The new theoretical developments presented in this work can be applied to spin systems near criticality.Away from the phase transition and at temperature much lower compared to T c with nonvanishing chiral condensate, the newly proposed quantity P n (λ) as well as the generalized Banks-Casher relation could be interesting to be investigated in the random matrix theory [41][42][43][44][45].The numerical techniques used here can be straightforwardly carried over to lattice QCD calculations with controlled thermodynamic and continuum limits when sufficient computing power is available.
We thank Yu Zhang for the early involvement of the work, Sheng-Tai Li for technical support, Jacobus Verbaarschot for discussions on the random matrix theory, and the members of the HotQCD Collaboration for numerous interesting discussions.This material is based The gauge configurations used in this work were generated using computing resources obtained through PRACE grants at CSCS, Switzerland, and at CINECA, Italy as well as grants at the Gauss Centre for Supercomputing and NICJülich, Germany.These grants provided access to resources on Piz Daint at CSCS and Marconi at CINECA, as well as on JUQUEEN and JUWELS at NIC.Additional calculations have been performed on GPU clusters of USQCD, at Bielefeld University, the PC2 Paderborn University, and NSC 3 .I. Summary of lattice parameters; i.e., values of lattice gauge coupling β, temperature T , strange quark mass ams in unit of lattice spacing, number of gauge configurations #conf used in both the direct measurements of chiral observables and in the computation of ρU via the Chebyshev polynomial method, as well as the number of random vectors Nvec used in the Chebyshev method shown in the square brackets.The order of Chebyshev polynomials are all set to 2 × 10 5 .The pion mass mπ and lattice size N 3 σ × Nτ corresponding to different gauge ensembles are also listed.

III.A. Reproduction of Kn ψψ via Pn(λ)
We show the temperature dependence of K n ψψ with n ≤ 3 in Fig. 3.The bin-size of λ needs to be fixed firstly when computing ρ U , P n and the numerical integration of Eq. 5.In this work the bin-size of λ was chosen such that the direct measurements of K 2 ψψ via fermion inversions (cf.Eq. 13) can be reproduced via Eq.14. Similarly as in Ref. [32,33], once the bin-size of λ is chosen to reproduce directly measured K 2 ψψ , the same value can also be used to reproduce other directly measured K n ψψ without any further tuning.The consistency of K n ψψ for n = 1, 2, 3 given by the two different methods as shown in Fig. 3 verifies the reliability of the results of P n (n = 1, 2, 3).Supplementary to the results shown in Fig. 1 we show results of Pn ( λ) with n = 1, 2, 3 at higher temperatures ranging from 147 MeV to 176 MeV in Fig. 4. Supplementary to the results shown in Fig. 2 we show results of Pn ( λ) rescaled by m1/δ+1−n l f n (z) for n = 1, 2, 3 at higher temperatures ranging from 147 MeV to 176 MeV in Fig. 5. P1 / ml  Supplementary to the results shown in Fig. 2 we show Pn ( λ) rescaled by m1/δ+1−n l f n (z) for n = 1, 2, 3 but using different values of {T c , z 0 } to obtain f n (z) in Fig. 6.The values of {T c , z 0 } adopted in Fig. 6 are T c = 145.6(3)MeV and z 0 = 2.24(5) as quoted in [17], which are obtained from scaling fits to lattice data of static quark energy in a broader pion mass region of 55 MeV≤ m π ≤140 MeV.λ FIG. 6.Similar to Fig. 2 but using Tc = MeV and z0 = 2.24(5) as also quoted in [17] to obtain fn(z).
arXiv:2305.10916v2 [hep-lat] 20 Oct 2023 eigenvalue spectrum through the state-of-the-art lattice QCD calculations in the staggered discretization scheme.Theoretical developments.-Forthe lattice QCD calculations we will use (2+1)-flavor QCD with degenerate light up (u) and down (d) quarks having masses m l = m u = m d and a heavier strange quark with physical mass m s .Since the physical strange quark plays no significant role in the discussion of the O(2) critical behavior, for simplicity in this subsection we develop the main theoretical idea by considering QCD with two degenerate light quarks.Consider the Euclidean-time QCD action S[U, m l ] = S g [U]+ ψ / D[U]ψ+m l ψψ, where ψψ = ψu ψ u + ψd ψ d , S g [U] is the pure gauge action, and / D[U] is the massless QCD Dirac operator for given background SU (3) gauge field U. To probe the system in chiral limit, m l → 0, we introduce a probe operator, ψψ(ϵ) ≡ 2Tr( / D[U] + ϵ) −1 , with the background U distributed according to exp{−S[U, 0]}.The valance quark mass, ϵ > 0, is introduced to facilitate the evaluation of the probe operator.Traces over the color, spin and space-time indices are denoted by Tr.

FIG. 2 .
FIG. 2. Pn( λ) in Fig. 1 by m1/δ+1−n l upon work supported partly by the National Key Research and Development Program of China under Contract No. 2022YFA1604900; the National Natural Science Foundation of China under Grants No. 12293064 , No. 12293060 and No. 12325508, and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics through Contract No. DE-SC0012704 and within the framework of Scientific Discovery through Advance Computing (SciDAC) award Fundamental Nuclear Physics at the Exascale and Beyond.Numerical simulations have been performed on the GPU cluster in the Nuclear Science Computing Center at Central China Normal University (NSC 3 ), Wuhan supercomputing center, the facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy.

FIG. 3 .
FIG. 3. Comparisons of direct measurements (open symbols) of K1 ψψ (left), K2 ψψ (middle) and K3 ψψ (right) (cf.Eq. 13) with those reproduced (filled symbols, slightly shifted horizontally for visibility) from Pn(λ) (cf.Eq. 5).The results are shown for all available values of light quark masses and temperatures.Lines which are just connections of data points and not fits are used to guide the eye.