Correlated Dirac Eigenvalues and Axial Anomaly in Chiral Symmetric QCD

We investigate the Dirac eigenvalue spectrum ($\rho(\lambda,m_l)$) to study the microscopic origin of axial anomaly in high temperature phase of QCD. We propose novel relations between the derivatives ($\partial^n \rho(\lambda,m_l)/\partial m_l^n$) of the Dirac eigenvalue spectrum with respect to the quark mass ($m_l$) and the $(n+1)$-point correlations among the eigenvalues ($\lambda$) of the massless Dirac operator. Based on these relations, we present lattice QCD results for $\partial^n \rho(\lambda,m_l)/\partial m_l^n$ ($n=1, 2, 3$) with $m_l$ corresponding to pion masses $m_\pi=160-55$ MeV, and at a temperature of about 1.6 times the chiral phase transition temperature. Calculations were carried out using (2+1)-flavors of highly improved staggered quarks and the tree-level Symanzik gauge action with the physical strange quark mass, three lattice spacings $a=0.12, 0.08, 0.06$ fm, and lattices having aspect ratios $4-9$. We find that $\rho(\lambda\to0,m_l)$ develops a peaked structure. This peaked structure, which arises due to non-Poisson correlations within the infrared part of the Dirac eigenvalue spectrum, becomes sharper as $a\to0$, and its amplitude is proportional to $m_l^2$. After continuum and chiral extrapolations, we find that the axial anomaly remains manifested in two-point correlation functions of scalar and pseudo-scalar mesons in the chiral limit. We demonstrate that the behavior of $\rho(\lambda\to0,m_l)$ is responsible for it.


Introduction
The Lagrangian of the (2+1)-flavor Quantum Chromodynamics (QCD) has a global symmetry (2) L × (2) R × (1) A × (1) in the classic limit and the chiral limit of → 0. The (2) L × (2) R chiral symmetry is spontaneously broken in the vacuum and the (1) A symmetry is anomalously broken on the quantum level due to the Adler-Bell-Jackiw or chiral anomaly. For the physical lattice simulations have established quite firmly that QCD transition is a rapid cross over at a pesudocritical temperature at 156 MeV [1][2][3], while in the chiral limit → 0 chiral phase transition temperature at which the (2) L × (2) R is restored is estimated as = 132 +3 −6 MeV based on the O(4) scaling analyses [4].
Conversely, the fate of the (1) A symmetry in the high temperature phase of QCD remains unclear. Although the quantum anomaly is present at any finite temperature, at some point its effects could become negligible due to the asymptotic restoration of the (1) A symmetry with the temperature, thus the (1) A symmetry would be effectively restored. The order of the chiral transition and the associated universality class is known to depend crucially on how axial anomaly manifests itself in the two-point correlation functions of light scalar and pseudoscalar mesons for ≥ . If the isotriplet scalar and the isotriplet pseudoscalar remain non-degenerate at ≥ , then we expect a second order phase transition which belongs to the three-dimensional (4) universality class [5]. But if the and become degenerate at ≥ , then the chiral phase transition can be either first [5] or second order with the symmetry breaking pattern (2) × (2) → (2) universality class [6,7]. For the physical , the and remain nondegenerate around the chiral crossover [2,8,9]. However, what happens for as → 0 remains an open question [10][11][12][13] due to the lack of state-of-the-art lattice QCD calculations with controlled continuum and chiral extrapolations.
To gain more insight about the microscopic origin of the axial anomaly we can investigate the Dirac eigenvalue spectrum ( , ). It has been shown that if ( , ) is an analytic function of 2 and then in the chiral limit (1) A breaking effects are invisible in differences of up to 6-point correlation functions of and that can be connected via a (1) A rotation [14]. However, the dilute instanton gas approximation (DIGA) [15] predicated that ∼ 2 ( ) can lead to nondegeneracy of the two-point and correlation functions even as → 0 [16][17][18]. Some lattice QCD studies have observed infrared enhancement in [8,11,16], however, whether such enhancements scale as 2 as → 0 have not been demonstrated. In other lattice QCD calculations, no infrared enhancement in was observed [12,19,20], showing the importance of controlling lattice artifacts through continuum extrapolations. On the other hand, in Ref. [21] it was argued that if and were to remain nondegenerate at ≥ , then chiral symmetry restoration demands non-Poisson correlations among the infrared eigenvalues.
In this work we propose the novel relation between / and correlation among the eigenvalues to investigate the microscopic origin of axial anomaly at high temperature phase. The rest of paper is organized as follows. We describe the basic idea of how to obtain the relation between / and correlation among the eigenvalues in section 2. In section 3 we show the setup of our lattice simulations. We then show our numerical results in Section 4. Finally we present our conclusion in Section 5. The detailed information about this work can be found in [22].

2.
/ & +1 and (1) anomaly For (2+1)-flavor QCD, the Dirac eigenvalue spectrum is given by Here ( ) is the Dirac eigenvalue spectrum for a given gauge configuration, It is defined as are the eigenvalues of the massless Dirac matrix / [U]. Note that ( ) does not explicitly depend on and the dependence is embedded in the determinant term. Furthermore, Substituting Eq. 2 in Eq. 1 and [U] it is straightforward to obtain / [22], e.g., with The difference of the integrated two-point functions in the pion and delta channel is defined as For ≥ owing to the degeneracy of and the in the chiral limit [16] where disc is the quark-line disconnected part of the isosinglet scalar meson susceptibility, The U(1) A symmetry-breaking measures − and disc are related to through [16,22] − = In the Poisson limit, In this limit the first and second order quark mass derivatives of are expressed as follows Y. Zhang In the chiral limit, this leads to Po disc = 2( − ), in clear violation of the chiral symmetry restoration condition in Eq. 7, unless both sides of the equation trivially vanish.

Lattice setup
Lattice QCD calculations were carried out at ≈ 205 MeV ≈ 1.6 for (2+1)-flavor QCD using the highly improved staggered quarks and the tree-level Symanzik gauge action. The was tuned to its physical value, and three lattice spacings = ( ) −1 = 0.12, 0.08, 0.06 fm corresponding to = 8, 12, 16, were used [22]. were done by inverting the light fermion matrix using 50 Gaussian random sources on 2000 − 10000 configurations. Apart from the data sets as shown in above which were reported in [22], in this paper we also add new results based on simulations with = /160 on = 12, 16 lattices. For each of these two parameter sets 4200-5200 configurations each separated by 10 time units are generated, and disc and − are measured by inverting the fermion matrix on these configurations.  Fig. 1 (left) shows the dependence of −1 / and 2 / 2 at ≈ 1.6 , obtained on = 8 and the largest available for that . We observe that −1 ( / ) and 2 / 2 are almost equal to each other and independent of . Also, −1 / and 2 / 2 develops a peak at → 0 and it drops rapidly toward zero for / 1. Fig. 1 (middle) depicts the lattice spacing and volume dependence of 2 / 2 and 3 / 3 for = 80 MeV. To compare these quantities across different lattice spacings we multiply with the appropriate powers of to make them renormalization group invariant and make them dimensionless by rescaling with appropriate powers of = 132 MeV. We see that the peaked structure in 2 / 2 at → 0 Y. Zhang becomes sharper as → 0, and shows little volume dependence. Moreover, 3 / 3 are found to be consistent with zero within errors. The findings −1 / ≈ 2 / 2 and 3 / 3 ≈ 0

Results
show that the peaked structure ( → 0, → 0) ∝ 2 . In Fig. 1 (right) we show the difference 2), with the Poisson approximations for / as defined in Eq. 10 and Eq. 11. The fact Δ Po < 0 shows that the repulsive non-Poisson correlation within the small gives rise to the ( → 0) peak.   (left) and disc (right) with those reconstructed (filled symbols, slightly shifted horizontally for visibility) from and / using Eq. 9.
In Fig. 2 we show that and / reproduce directly measured − and disc using Eq. 9. We checked that only the infrared / 1 parts of and / are needed for the reproductions of − and disc . Additionally, we checked that once the bin-size of in the numerical integration of left equation of Eq. 9 is chosen to reproduce directly measured − , the same bin size automatically reproduces disc and ¯ without any further tuning. In the left panel of Fig. 3 we show the quark mass dependence of chiral condensate in detail. We performed linear fits (solid lines) and a quadratic fits (dotted lines) in quark mass to the chiral Y. Zhang condensate. It can be clearly seen that the linear fits give a good description of the data and the fit result of chiral condensates at each lattice spacing vanish in the chiral limit. This is in accord with the expectation [U] is an even function of for ≥ due to the restoration of the (2) subgroup of (2) L × (2) R . This leads to the expectation that the disc should be quadratic in quark mass as → 0. As can be seen from the right panel of Fig. 3 which shows the 2 disc / 4 as a function of quark mass for = 8, 12, 16, the data indeed favors the quadratic dependence of 2 disc / 4 in quark mass as → 0. In Fig. 4 we show the continuum and chiral extrapolated results for disc and − . With the additional 2 data points at / = 160 (or = 55 MeV) on = 12 and 16, we follow the same analysis methods as in our previous studies [22]. I.e. using data for = 8, 12, 16 and ≤ 140 MeV, we performed a joint , → 0 extrapolation of the form disc ( , ) = disc (0, 0) + 1 / 2 + 2 / 4 + ( / ) 2 0 + 1 / 2 + 2 / 4 . Fits were performed on each bootstrap sample of the data set. The bootstrap samples were created by randomly choosing data from Gaussian distributions with means equal to the average values and variances equal to the errors of disc . We chose the median value as the final result (depicted by the upward triangles) and the 68% percentiles confidence interval of the resulting bootstrap distribution as the errors (the band labeled by 8,12,16 − −−−− → ∞). Since we used the so-called rooted-staggered formulation [27], we also checked that the same disc (0, 0) is obtained within errors by first carrying out the → 0 extrapolations for each and then performing the → 0 extrapolation. For this purpose, we used the = 12, 16 data for each of = /27, /40, /80, /160 to obtain disc (0, ) by fitting to the ansatz disc ( , ) = disc (0, ) + 1 / 2 . Then the chiral extrapolation was carried out using disc (0, ) = disc (0, 0) + 2 ( / ) 2 based on the continuum estimates of disc (0, ). These extrapolations were done by using the same bootstrap procedure described before and the final results are indicated with the label 12,16 − −−− → ∞. The same procedures were followed also for − to obtain its continuum and chiral extrapolated values. After carrying out continuum and chiral extrapolations we obtained that disc (0, 0) is 3.0 ± 1.1 for the sequential fit and 5.7 ± 2.3 for the joint fit, which is 2-3 away from 0, while [ − ] (0, 0) is 6.7 ± 1.1 for the sequential fit In the staggered discretization formalism the remnant chiral symmetry at nonzero lattice spacing is O(2). Y. Zhang and 7.8 ± 2.2 for the joint fit, which is 4-6 away from 0. We find that Eq. 7 is satisfied within errors, and disc and − are nonvanishing at a confidence level above 95%. These results are consistent with those obtained without the two additional data points [22].

Conclusions
In this work we establish relations between / and +1 . Based on these relations, we present direct computations of / employing state-of-the-art lattice QCD techniques. Based on these results we conclude that, in chiral symmetric (2+1)-flavor QCD at ≈ 1.6 , (i) ( → 0, ) develops a peaked structure due to repulsive non-Poisson correlations within small ; the peak becomes sharper as → 0, and its amplitude is ∝ 2 . (ii) The underlying presence of this ( → 0, ) leads to manifestations of (1) A anomaly in − and disc . (iii) Axial anomaly remains manifested in − and disc even in the chiral limit. These suggest that for ∼ 1.6 the microscopic origin of axial anomaly is driven by the weakly interacting (quasi)instanton gas motivated ( → 0, → 0) ∼ 2 ( ), and the chiral phase transition in (2+1)-flavor QCD is of the three-dimensional (4) universality class.
The above conclusions are based on the continuum extrapolated lattice QCD calculations using the (2+1) flavors of staggered fermions. Confirmations of these continuum extrapolated results using other fermion actions, especially using chiral fermions, are needed in future. Even in those future calculations it will be very difficult to directly identify a structure like 2 ( ) in itself as → 0. The formalism developed and techniques presented in this work for directly accessing / will be essential for those future studies too.