Dirac-mode expansion of quark number density
and its implications of the confinement-deconfinement transition
We investigate the quark number density at finite imaginary chemical potential by using the Dirac-mode expansion. In the large quark mass region, it is found that the quark number density can be expressed by the Polyakov loop and its conjugate in all order of the large quark mass expansion. Then, there are no specific Dirac-modes which dominantly contribute to the quark number density. In comparison, the small quark mass region is explored by using the quenched lattice QCD simulation. We found that the absolute value of the quark number density strongly depends on the low-lying Dirac-modes, but its sign does not. This means that the existence of the Roberge-Weiss transition which is characterized by the singular behavior of the quark number density is not sensitive to low-lying Dirac-modes. This property enables us to discuss the confinement-deconfinement transition from the behavior of the quark number density via the quark number holonomy.
Understanding the phase structure of Quantum Chromodynamics (QCD) at finite temperature () and real chemical potential is one of the interesting and important subjects in nuclear and elementary particle physics. The chiral and confinement-deconfinement transitions are key phenomena for this purpose, but the confinement-deconfinement transition is not yet fully understood comparing with the chiral transition. Although the chiral transition can be described by the spontaneous breaking of the chiral symmetry, but we can not find any classical order-parameters of the confinement-deconfinement transition in the system with dynamical quarks.
The Polyakov loop respecting the gauge-invariant holonomy becomes the exact order-parameter of the spontaneous symmetry breaking and then it has direct relation with the confinement-deconfinement transition in the infinite quark mass limit. However, it is no longer the order parameter if the quark mass is finite. Other quantity such as the dual quark condensate Bilgici:2008qy (); Fischer:2009wc (); Kashiwa:2009ki (); Benic:2013zaa (); Xu:2011pz () which has been proposed to characterize the confinement-deconfinement transition still shares the same problem since it is based on the spontaneous symmetry breaking. Therefore, we need some extension of ordinary determinations to clearly discuss and investigate the confinement-deconfinement transition in the system with dynamical quarks. One attempt is calculating the Polyakov-loop fluctuation Lo:2013hla (); Doi:2015rsa () which is originally introduce to attack the renormalization content of the Polyakov-loop, but now it is considered as the good indicator of the confinement-deconfinement transition comparing with the ordinary analysis of the Polyakov loop.
Recently, the behavior of the quark number density at finite imaginary chemical potential () has been used to determine the confinement-deconfinement transition based on topological properties of QCD. In Refs. Kashiwa:2015tna (); Kashiwa:2016vrl (); Kashiwa:2017yvy (), it has been proposed that the topological change of QCD thermodynamics at finite can be used to determine the confinement-deconfinement transition. This determination is based on the analogy of the topological order discussed in the condensed matte physics Wen:1989iv () and QCD at Sato:2007xc (). In discussions of the topological order at , the ground-state degeneracy plays a crucial role to determine the topologically ordered and dis-ordered phases. The nontrivial free-energy degeneracy in QCD induced by the Roberge-Weiss (RW) transition at finite Roberge:1986mm () can be considered as the analog of the ground-state degeneracy at zero temperature Kashiwa:2015tna ().
Based on the non-trivial free-energy degeneracy, the quark number holonomy which is defined by the contour integral of the quark-number susceptibility of has been proposed as the quantum order-parameter for the confinement-deconfinement transition. Since the quark number holonomy counts gapped points of the quark number density along and thus it becomes non-zero (zero) in the deconfined (confined) phase. Particularly, the quark number density at is quite important for the topologically determined confinement-deconfinement transition because the RW transition is expected to be happen here. Below, we use the term topological confinement-deconfinement transition when we determine the confinement-deconfinement transition by using the quark number holonomy. There are some papers which investigated the quark number density at finite DElia:2002tig (); DElia:2004ani (); Takahashi:2014rta (), but properties of the quark number density are not well understood yet. Therefore, we investigate the quark number density and its implications of the topological confinement-deconfinement transition in this article.
To investigate the quark number density and its implications of the topological confinement-deconfinement transition, the Dirac-mode expansion is a powerful and convenient tool. The chiral condensate and also the Polyakov loop are already expressed in terms of Dirac eigenvalues to investigate the relation between the chiral symmetry breaking and confinement Gongyo:2012vx (); Iritani:2013pga (); Suganuma:2014wya (). Then, we can see that dominant Dirac-modes in the chiral condensate and Polyakov loop are quite different; the dominant modes are low-lying Dirac-modes in the chiral condensate, but there are no specific modes in the Polyakov loop in the case with light quark masses. This fact may indicate that the behavior of dominant Dirac-modes can used to clarify which quantities are sensitive to the confinement-deconfinement transition. Therefore, it is interesting to analyze which behaviors is realized in the quark number density at finite .
This paper is organized as follows. In the next section, we show the the heavy quark mass (hopping parameter) expansion of the quark number density. The Dirac-mode expansion of the quark number density are discussed in Sec. III. Then, we estimate the quark number density in the large and also small quark mass regions by using the analytic calculation and the quenched lattice QCD simulation, respectively. Section IV is devoted to summary and discussions.
Ii Heavy-mass expansion of quark number density
In this study, we consider the lattice QCD on the standard square lattice. We denote each sites as and link-variables as . We impose the temporal periodic boundary condition for link-variables to generate configurations in the quenched calculation to manifest the imaginary-time formalism.
On the lattice, the quark number density is defined as
where denotes the functional trace and is taken over spinor and color indices. The operator is the Dirac operator and expresses the quark mass. In this article, we use the Wilson-Dirac operator with the chemical potential in the lattice unit as
where is the identity matrix and the link-variable operator is defined by the matrix element
with and with . The derivative is explicitly written as
It should be noted that is the dimensionless chemical potential on the lattice and it relates to as .
In the imaginary-time formalism, the temporal anti-periodicity for should be imposed to manifest the anti-periodic boundary condition of quarks. To that end, we add a minus sign to the matrix element of the temporal link-variable operator at the temporal boundary of :
In this notation, the Polyakov loop is expressed as
The minus sign stems from the additional minus on in Eq.(5).
ii.1 Leading-order contribution
In the heavy quark mass region, the quark number density (1) can be expressed by using the quark mass expansion as
where we define the effective mass and the operator . In Eq. (8), with smaller are relevant because the effective mass is supposed to be large here. Since the operator and consists of linear terms in link-variables , the -th order contribution can have many terms
where each is a product of link-variables; see Fig. 1 as an example of paths on the lattice. It should be noted that many of them become exactly zero because of Elitzur’s theorem Elitzur:1975im (); only the gauge-invariant terms corresponding to closed loops are nonzero. Moreover, spatially closed loops which do not wind the temporal length are canceled out each other and have no contribution to the quark number density in total. Noting these important facts, it is confirmed that the -th order contribution is constituted of the gauge-invariant loops with the length which winds the temporal direction. In particular, the nontrivial leading term in the expansion (8) is the -th order term which relates to the Polyakov loop () and its complex conjugate ().
Above result can be also understood from the -dependence of the quark number density. Due to at , -independent terms in the expansion (8) must vanish. As shown in Fig. 1, in the spatially closed loop, the -dependence are canceled out between and . Thus, the nonzero terms must wind the temporal direction using the temporal periodic boundary condition. The other -dependent terms up to order cannot form the gauge-invariant closed-loop and thus these should be vanished by Elitzur’s theorem. Therefore, the terms with which are proportional to and are the dominant contributions for the quark number density because is the shortest loop which leads to the -dependence.
The actual contributions of the leading term, , to the quark number density at finite take the form;
Since , equivalently , can be translated into the temporal boundary condition of quarks, and can feel -effects: at exists in the trivial center region (), but it can stay in the non-trivial center region () at with . For example, at sufficiently high , we can obtain in , in and in because of the RW transition; see Ref. Roberge:1986mm () as an example. The RW transition is also characterized by the gap of the quark number density; the quark number density has opposite sign at both side of the line.
In the heavy quark mass region, the system is almost the quenched system and then we may generate gauge configurations by only using the pure gauge action. However, if we perfectly take the quenched limit, the RW periodicity should vanish because -effects can not modify configurations. Thus, we should take into account the discrete -hopping, for the variation of as a perturbation to reproduce the RW periodicity and transition. This treatment is implicitly used in the holographic QCD calculations in the probe limit Aarts:2010ky (); Rafferty:2011hd (); Isono:2015uda (). However, in the low region, the discrete -hopping treatment may provide wrong results since is smoothly rotating with varying of . However, is exactly in the quenched calculation at low and thus there is no need to care the actual value of . In the following discussion, we consider the region and thus the system is in the trivial center region.
ii.2 Higher-order contributions
In the heavy quark-mass expansion of the quark number density (8), there are higher order terms beyond the leading terms. Noting again the important point in the previous subsection, each higher order term corresponds to a loop which winds the temporal direction and has the longer length . Specifically, the sub-leading contributions are order’s terms. For example, they include a term proportional to the quantity
The sub-leading terms correspond to the closed paths which wind the temporal length and make a detour in the spatial direction. As another example, loops winding the temporal direction twice or more can be contribute to the expansion (8). For example, a loop winding the temporal direction twice
is a possible contribution to Eq. (8) as the ()-th order term. The examples of the higher order terms are shown in Fig. 2. In the next section, we show that all the terms in the heavy quark-mass expansion (8) can be analytically represented in terms of the Dirac eigenmodes as well as the leading terms, namely the Polyakov loop and its conjugate.
Iii Dirac-mode expansion of the quark number density
In this section, we consider the quark number density in terms of the Dirac mode. In large quark mass region, we analytically show the statement in all order of the large quark mass expansion (8). In small quark mass region, we show that from numerical analysis by performing the quenched lattice QCD simulation.
The Wilson-Dirac eigenvalues are obtained from the eigenvalue equation as
where is the Wilson-Dirac eigenstate. For example, the chiral condensate is defined as
where is the four-dimensional volume. Considering the Wilson-Dirac mode expansion of the chiral condensate, the low-lying eigenmodes of the operators have dominant contribution to the chiral condensate known as Banks-Casher relation Banks:1979yr (); Giusti:2008vb ().
iii.1 large quark mass region
We start the leading term to express it in terms of the Wilson-Dirac modes. The leading contribution of the quark number density in large quark mass region (10) is expressed by the Polyakov loop and its complex conjugate. It is already known that the Polyakov loop can be expressed in terms of the eigenmodes of the naive-Dirac operator which corresponds to the case of Suganuma:2014wya (); Doi:2014zea () and the Wilson-Dirac operator Suganuma:2016lnt (); Suganuma:2016kva (). In the following, we derive a different form of the Dirac spectral representation of the the Polyakov loop using the operator on the square lattice with the normal non-twisted periodic boundary condition for link-variables, in both temporal and spatial directions. We firstly define the following key quantity,
This quantity is defined by changing a temporal link-variable of to the Wilson-Dirac operator . Substituting the definition (2) of the Wilson-Dirac operator , the quantity can be calculated as
Thus, the quantity is proportional to . Note that (other terms) vanish because of the Elitzur’s theorem and the trace over the Dirac indecies. On the other hand, since in Eq. (15) is defined through the functional trace, it can be expressed in the basis of Dirac eigenmodes as
The term arises because Wilson-Dirac operator is not normal due to the Wilson term and the completeness of the Wilson-Dirac eigenstates has the error:
However, this error is controllable and can be ignored in close to the continuum limit.
with the error. This relation is a Dirac spectral representation of the Polyakov loop. From the formula (19), it is analytically found that the low-lying Dirac modes with have negligible contribution to the Polyakov loop because the eigenvalue plays as the damping factor and the quantity has the finite value . It is also numerically shown that there is no dominant contribution in the Dirac modes to the Polyakov loop because the contributions of the Dirac modes whose eigenvalues are almost same are canceled due to the positive/negative symmetry of the Dirac-matrix element of the link-variable Doi:2014zea (). Thus, one can find the Dirac spectrum representation of the quark number density in the leading order and the low-lying Dirac modes have little contribution to the quark number density.
The above discussion on the leading term can be applicable to the Dirac spectrum representation of the higher order terms. For example, a sub-leading term can be expressed as
by considering the functional trace
instead of Eq. (15). Moreover, the loop which winds the temporal direction twice can be expressed as
by considering the different functional trace
In the same way, all the terms in the expansion (8) can be expressed in terms of the Wilson-Dirac modes. Thus, it is analytically found that the quark number density is insensitive to the density of the low-lying Wilson-Dirac modes in the all-order. However, this fact is only valid in the sufficiently large quark mass region since other contributions which can not be expressed by and can appear in the small region. Actually, every quark bilinears including the chiral condensate show the same behavior of the quark number density in terms of Dirac modes in the large quark mass region. However, we already know that the dominant contributions of the chiral condensate are low-lying Dirac modes in the small quark mass region.
iii.2 Small quark mass region
In the small quark mass region, the quark number density requires contributions which can not be expressed by and and thus we perform the lattice QCD simulation to investigate the quark number density. In this study, we perform the quenched calculation with the ordinary plaquette action and then fermionic observables are evaluated by using the Wilson-Dirac operator (2) with the imaginary chemical potential . Our calculation is performed in both the confinement phase and the deconfinement phase. In the confinement phase, we consider lattice with and which corresponds to fm and MeV. In the deconfinement phase, we consider lattice with and which corresponds to fm and MeV. Both values of correspond to . In both cases, we set the quark mass as in the lattice unit, which is equivalent to the hopping parameter , for the calculation of the eigenmodes of the Wilson-Dirac operator in the small quark mass region Aoki:1999yr ().
In the lattice QCD simulation, there is the smearing of the phase transition because of the finite size effect. Of course, the RW transition which plays a crucial role in the behavior of the quark number holonomy should be smeared in the system with dynamical quarks. However, in the quenched lattice QCD simulation, the RW transition can reproduced by considering the discrete -hopping and then the smearing of the phase transition does not matter. In addition, we can clarify the critical temperature of the confinement-deconfinement transition from previous quenched lattice QCD simulations because the RW endpoint and the ordinary critical temperature determined by the Polyakov loop are consistent with each other in the quenched lattice QCD simulation. This knowledge enable us to safely clarify which temperature is above the critical temperature when we investigate the quark number density.
The eigenstates of the Wilson-Dirac operator do not form the complete system, and then we calculate the quark number density as
This form trivially takes pure imaginary value up to the error. Each contribution, , to the quark number density of the Dirac mode with can be defined as
and then the quark number density becomes
Top (bottom) panel of Fig. 4 shows each Dirac-mode contribution to the quark number density, , with and . We here only show results with one particular configuration.
By comparing with both figures, in the distribution of at , the positive-negative symmetry Doi:2014zea () is almost realized, which was originally found in the confined phase at . On the other hand, in the case of , the symmetry is broken. This means that the violation of the positive-negative symmetry leads to at finite .
To investigate the quark number density in terms of Dirac modes, we show the infra-red (IR) cutted quark number density with the cutoff defined by
at in Fig. 5. In the evaluation of Eq. (27), the configuration averaging is basically possible but it misses a physical meaning because the averaging well works after summing over all Dirac-modes. Thus, we here show in one particular configuration.
We can see that the absolute value of the quark number density seems to depend on the low-lying Dirac modes, but its sign does not. This tendency can be found in almost all our configurations. It means that the absolute value of the quark number density shares a same property in terms of Dirac modes with the chiral condensate, while its sign shares the property with the Polyakov loop.
Finally, we discuss the topological confinement-deconfinement transition from Dirac-mode analysis. The order-parameter of the topological confinement-deconfinement transition can be expressed Kashiwa:2016vrl () as
where is the dimensionless quark number density such as . It counts gapped points of the quark number density along direction and thus it becomes zero (non-zero) in the confined (deconfined) phase. Equation (28) can be expressed as
when the RW endpoint which is the endpoint of the RW transition line becomes the second-order point at . In Eq. (29), characterizes and thus shares the same property about the Dirac modes with . The important point is that the absolute value of does not have so much meaning even if it is non-zero, but its sign is important since the sign flipping at characterizes the gapped points along -direction. From our quenched lattice QCD data, the sign of the quark number density are insensitive to the low-lying Dirac-modes and this behavior is similar to the Polyakov loop. Equation (19) holds also in the small quark mass regime because its derivation does not depend on the quark mass. It should be noted that recent lattice QCD data D'Elia:2009qz (); Bonati:2010gi () predict that the RW endpoint seems to be the triple-point where three first-order transition lines meet. We should modify the expression (29) in this case, but we can expect that the dependence of the Dirac modes is same with that of the second-order RW endpoint scenario; can be expressed by using where is the endpoint value of of the triple line at .
The flavor lattice QCD simulation predicts that the RW endpoint is MeV Bonati:2016pwz () and thus it becomes the critical temperature of the topological confinement-deconfinement transition when the RW endpoint is second order. On the other hand, the chiral pseudo critical temperature is about MeV Aoki:2006br (); Aoki:2009sc (); Borsanyi:2010bp (); Bazavov:2011nk (). This result strongly supports that the topologically determined confinement-deconfinement transition does not have the exact one-to-one correspondence with the chiral phase transition.
Iv Summary and discussion
In this paper, we have discussed properties of the quark number density at finite temperature () and imaginary chemical potential () in terms of Dirac-modes. In the heavy quark mass region, we use the heavy quark mass expansion and then the analytic form is discussed. On the other hand, we employ the quenched lattice QCD simulation in the small quark mass region.
From the heavy quark mass expansion with the Dirac mode expansion, we found that low-lying Dirac modes do not dominantly contribute to the quark number density in all order of the heavy quark mass expansion. Some other quark bilinears should also be insensitive to low-lying Dirac modes. This result is valid if the quark number density can be well expressed by the Polyakov loop () and its conjugate . If some other contributions which can not be expressed by and appear, low-lying Dirac modes can become dominant modes.
The small quark mass region has been explored by using the quenched lattice QCD simulation. We found that the absolute value of the quark number density strongly depends on low-lying Dirac modes, but its sign does not. Therefore, the sign of the quark number density shares similar properties with the Polyakov loop in terms of Dirac eigenvalues. This result means that the RW transition at is insensitive to low-lying Dirac modes where .
The RW transition plays a crucial role in the determination of the topological confinement-deconfinement transition and thus we can discuss the transition from the viewpoint of Dirac eigenvalues. The order-parameter of the topological confinement-deconfinement transition is the quark number holonomy () which is defined by the contour integral of the quark number susceptibility along . This quantity becomes non-zero if the quark number density has the gap along with fixed and then the deconfined phase is realized. From our quenched lattice QCD data, it is found that the quark number holonomy is sensitive to the confinement properties of QCD because the sign of the quark number density shares similar properties with the Polyakov loop. Therefore, our results support that the quark number holonomy is the good quantum order parameter for the confinement-deconfinement transition.
It is interesting to compare the quark number holonomy with other quantities which relates with the confinement-deconfinement transition. The dual quark condensate defined with the twisted boundary condition is one of sensitive probes for the quark-deconfinement Bilgici:2008qy (); Fischer:2009wc (); Kashiwa:2009ki (); Benic:2013zaa (); Xu:2011pz (). In the calculation of the dual quark condensate, we must break the RW periodicity because it should be zero if the RW periodicity exists. It is usually done by imposing the twisted boundary condition on the Dirac operator, while configurations are created under the anti-periodic boundary condition in the system with dynamical quarks. However, this procedure is not unique. Therefore, there is the uncertainty in the determination of the dual quark condensate. Also, it is well known that this quantity is strongly affected by other phase transitions Benic:2013zaa (); Marquez:2015bca (); Zhang:2015baa (). Another famous one is the QCD monopole. The QCD monopole which is also the order-parameter of the chiral symmetry appears by fixing the maximally Abelian gauge and plays important role in the dual-superconductor picture for the mechanism of confinement Suganuma:1993ps (). In fact, after removal of the QCD monopole from the QCD vacuum generated in the lattice QCD, both chiral restoration and quark deconfinement occur Miyamura:1995xn (); Woloshyn:1994rv (). The quark number holonomy has advantages over the QCD monopole and the dual quark condensate in the view point of the gauge invariance and the temporal boundary condition. While the QCD monopole is gauge-variant, the quark number holonomy is gauge-invariantly defined. Also, on the one hand, the dual quark condensate can be defined with only the twisted boundary condition. On the other hand, the quark number holonomy can be defined with arbitrary boundary condition, including the periodic boundary condition, which is needed for finite-temperature formalism. Also, the dual quark condensate has uncertainty in the actual calculation process, but the quark number holonomy does not. Thus, the quark number holonomy is superior order parameter of the confinement-deconfinement transition. Moreover, the quark number holonomy is better than the Polyakov loop because the Polyakov loop works as the order parameter only at the large quark-mass regime.
Acknowledgments: T.M.D. is supported by the Grantin-Aid for JSPS fellows (No.15J02108) and the RIKEN Special Postdoctoral Researchers Program.
- (1) E. Bilgici, F. Bruckmann, C. Gattringer, and C. Hagen, Phys.Rev. D77, 094007 (2008), arXiv:0801.4051 [hep-lat]
- (2) C. S. Fischer, Phys.Rev.Lett. 103, 052003 (2009), arXiv:0904.2700 [hep-ph]
- (3) K. Kashiwa, H. Kouno, and M. Yahiro, Phys.Rev. D80, 117901 (2009), arXiv:0908.1213 [hep-ph]
- (4) S. Benič, Phys.Rev. D88, 077501 (2013), arXiv:1305.6567 [hep-ph]
- (5) F. Xu, H. Mao, T. K. Mukherjee, and M. Huang, Phys.Rev. D84, 074009 (2011), arXiv:1104.0873 [hep-ph]
- (6) P. M. Lo, B. Friman, O. Kaczmarek, K. Redlich, and C. Sasaki, Phys. Rev. D88, 074502 (2013), arXiv:1307.5958 [hep-lat]
- (7) T. M. Doi, K. Redlich, C. Sasaki, and H. Suganuma, Phys. Rev. D92, 094004 (2015), arXiv:1505.05752 [hep-lat]
- (8) K. Kashiwa and A. Ohnishi, Phys. Lett. B750, 282 (2015), arXiv:1505.06799 [hep-ph]
- (9) K. Kashiwa and A. Ohnishi, Phys. Rev. D93, 116002 (2016), arXiv:1602.06037 [hep-ph]
- (10) K. Kashiwa and A. Ohnishi(2017), arXiv:1701.04953 [hep-ph]
- (11) X. G. Wen, Int.J.Mod.Phys. B4, 239 (1990)
- (12) M. Sato, Phys.Rev. D77, 045013 (2008), arXiv:0705.2476 [hep-th]
- (13) A. Roberge and N. Weiss, Nucl.Phys. B275, 734 (1986)
- (14) M. D’Elia and M.-P. Lombardo, Phys. Rev. D67, 014505 (2003), arXiv:hep-lat/0209146 [hep-lat]
- (15) M. D’Elia and M. P. Lombardo, Phys. Rev. D70, 074509 (2004), arXiv:hep-lat/0406012 [hep-lat]
- (16) J. Takahashi, H. Kouno, and M. Yahiro, Phys. Rev. D91, 014501 (2015), arXiv:1410.7518 [hep-lat]
- (17) S. Gongyo, T. Iritani, and H. Suganuma, Phys. Rev. D86, 034510 (2012), arXiv:1202.4130 [hep-lat]
- (18) T. Iritani and H. Suganuma, PTEP 2014, 033B03 (2014), arXiv:1305.4049 [hep-lat]
- (19) H. Suganuma, T. M. Doi, and T. Iritani, PTEP 2016, 013B06 (2016), arXiv:1404.6494 [hep-lat]
- (20) S. Elitzur, Phys. Rev. D12, 3978 (1975)
- (21) G. Aarts, S. P. Kumar, and J. Rafferty, JHEP 07, 056 (2010), arXiv:1005.2947 [hep-th]
- (22) J. Rafferty, JHEP 09, 087 (2011), arXiv:1103.2315 [hep-th]
- (23) H. Isono, G. Mandal, and T. Morita(2015), arXiv:1507.08949 [hep-th]
- (24) T. Banks and A. Casher, Nucl. Phys. B169, 103 (1980)
- (25) L. Giusti and M. Luscher, JHEP 03, 013 (2009), arXiv:0812.3638 [hep-lat]
- (26) T. M. Doi, H. Suganuma, and T. Iritani, Phys. Rev. D90, 094505 (2014), arXiv:1405.1289 [hep-lat]
- (27) H. Suganuma, T. M. Doi, K. Redlich, and C. Sasaki(2016), arXiv:1610.02999 [hep-lat]
- (28) H. Suganuma, T. M. Doi, K. Redlich, and C. Sasaki, Proceedings, 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII): Thessaloniki, Greece, EPJ Web Conf. 137, 04003 (2017), arXiv:1611.08742 [hep-th]
- (29) S. Aoki et al. (CP-PACS), Phys. Rev. Lett. 84, 238 (2000), arXiv:hep-lat/9904012 [hep-lat]
- (30) M. D’Elia and F. Sanfilippo, Phys. Rev. D80, 111501 (2009), arXiv:0909.0254 [hep-lat]
- (31) C. Bonati, G. Cossu, M. D’Elia, and F. Sanfilippo, Phys.Rev. D83, 054505 (2011), arXiv:1011.4515 [hep-lat]
- (32) C. Bonati, M. DâElia, M. Mariti, M. Mesiti, F. Negro, and F. Sanfilippo, Phys. Rev. D93, 074504 (2016), arXiv:1602.01426 [hep-lat]
- (33) Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, Phys. Lett. B643, 46 (2006), arXiv:hep-lat/0609068 [hep-lat]
- (34) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg, and K. K. Szabo, JHEP 06, 088 (2009), arXiv:0903.4155 [hep-lat]
- (35) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo (Wuppertal-Budapest), JHEP 09, 073 (2010), arXiv:1005.3508 [hep-lat]
- (36) A. Bazavov et al., Phys. Rev. D85, 054503 (2012), arXiv:1111.1710 [hep-lat]
- (37) F. Marquez, A. Ahmad, M. Buballa, and A. Raya, Phys. Lett. B747, 529 (2015), arXiv:1504.06730 [nucl-th]
- (38) Z. Zhang and Q. Miao(2015), arXiv:1507.07224 [hep-ph]
- (39) H. Suganuma, S. Sasaki, and H. Toki, Nucl. Phys. B435, 207 (1995), arXiv:hep-ph/9312350 [hep-ph]
- (40) O. Miyamura, Phys. Lett. B353, 91 (1995)
- (41) R. M. Woloshyn, Phys. Rev. D51, 6411 (1995), arXiv:hep-lat/9503007 [hep-lat]