Nature of the chiral phase transition of two flavour QCD from imaginary chemical potential with HISQ fermions

The nature of the thermal phase transition of two flavor QCD in the chiral limit has an important implication for the QCD phase diagram. We carry out lattice QCD simulations in an attempt to address this problem. Simulations are conducted with a Symanzik-improved gauge action and the HISQ fermion action. Within the imaginary chemical potential formulation, five different quark masses, $am=0.020,\, 0.018, \, 0.015, \, 0.013,\, 0.010$, and four different lattice volumes $N_s=8, \, 12,\, 16, \, 20$ with temporal extent $N_t=4$ are used to explore the scaling behavior. At each of the quark masses, the Binder cumulants of the chiral condensate on different lattice volumes approximately intersect at one point. We find that at the intersection point, the Binder cumulant $B_4(am,a\mu_c) $ is around $3$ which deviates from the $Z(2)$ universality class value 1.604. However, based on the expectations of $Z(2)$ criticality, the fitting result only with the data from the largest lattice volume $N_s=20$ agrees well with earlier result [ Phys. Rev., D90, 074030(2014) ]\cite{Bonati:2014kpa}. This fact implies that, although the finite cut-off effects could be reduced with HISQ fermions even on $N_t=4$ lattices, larger lattices with spatial extent $N_s>=20$ for such studies are needed to control finite volume effects.


I. INTRODUCTION
The thermodynamics of matter described by QCD is characterized by a transition from the low-temperature hadronic phase with confined quarks and gluons to the high-temperature phase with deconfined quarks and gluons. This phase transition is relevant to the early Universe, compact stars, and heavy ion collision experiments. Reviews on the study of the phase diagram can be found in Refs. [1][2][3] and references therein. Mapping out the phase diagram of QCD is one of the most challenging tasks presented for theoretical physics. Although substantial progress has been achieved in determining the phase diagram of QCD at zero density, the nature of the phase transition of QCD with two massless flavors remains still open.
At the transition point, if the U A (1) symmetry is not restored, QCD with two massless flavors has the symmetry-breaking pattern [SU (2) L ×SU (2) R ]/Z(2) V → SU (2) V /Z(2) V , on the other hand, if the U A (1) symmetry is effectively and fully restored, QCD with two massless flavors has the symmetry-breaking pattern [U (2) L × U (2) R ]/U (1) V → U (2) V /U (1) V [4][5][6]. For two-flavor QCD, Pisarsky and Wilczek pointed out that, if the U A (1) symmetry is broken at transition point T c , the sys- * Corresponding author. Email address: wuliangkai@163.com tem undergoes second-order transition with O(4) scaling, although not necessarily. On the other hand, if the U A (1) symmetry is restored at T c , the system undergoes a firstorder transition. However, further studies [5,6] show that, even if the U A (1) symmetry is restored at T c , the system also may have an infrared stable fixed point, so the transition can be of either first order or second order with different critical exponents from the O(4) universality class. Reference [7] suggests the transition is of second order, but one of critical exponents is different from the standard O(4) model.
As the interaction between the quarks and gluons is inherently strong at hadronic energy scales, lattice QCD simulation is the most reliable method up to date. The standard method to address the nature of QCD with two massless flavors is to carry out simulations by successively reducing the quark mass, and in the meantime, monitoring the transition behaviour. If the transition is of second order in the chiral limit, then this second-order transition disappears immediately when finite quark masses are turned on. On the other hand, if the transition is of first order in the chiral limit, it will get weakened gradually until at a certain Z(2) point as the quark masses increase.
Apart from the conventional method which focuses on the critical exponents, the nature of the phase transition of QCD with two massless flavors can be addressed by exploring the fate of U A (1) symmetry at high temperature [20,21,[25][26][27][28][29]. Reference [18] discusses this problem from the aspect of noninteger number of flavors. In Ref. [16], a novel approach has been developed to address the nature the phase transition of QCD with two massless flavors within the staggered fermion formulation, and this approach is employed in Ref. [17] within the Wilson fermion formulation. The approach takes advantage of the fact that when the imaginary chemical potential is switched on, the second-order line which separates the first-order region from the crossover region is governed by the tricritical scaling law, and the critical exponents around am = 0 are known [16,17,30].
So far, the investigation to address the nature of the phase transition of QCD with two massless flavors using this method are implemented through standard gauge and fermion actions [16,17]. Standard KS fermions suffer from taste symmetry breaking at nonzero lattice spacing a [31]. This taste symmetry breaking can be illustrated by the smallest pion mass taste splitting which is comparable to the pion mass even at lattice spacing a ∼ 0.1f m [32].
In this paper, we aim to investigate the nature of the phase transition of QCD with two massless flavors with a one quark-loop Symanzik-improved gauge action [33][34][35][36][37][38] and the HISQ action [39]. The one quark-loop Symanzik-improved gauge action has a discretization error of O(α 2 s a 2 , a 4 ), and the HISQ action completely eliminates the O(a 2 ) error at tree level by including smeared one-link and "Naik terms" [40,41]. Moreover, the HISQ action yields the smallest violation of taste symmetry among the currently used staggered actions [31,42,43]. These improvements are significant on the N t = 4 lattice where the lattice spacing is quite large.
The paper is organized as follows. In Sec. II, we define the lattice action with imaginary chemical potential and the physical observables we calculate. Our simulation results are presented in Sec. III, followed by discussion in Sec. IV.

II. LATTICE FORMULATION WITH IMAGINARY CHEMICAL POTENTIAL
After introducing a pseudofermion field Φ, the partition function of the system can be represented as: where S g is the Symanzik-improved gauge action, and S f is the HISQ quark action with the quark chemical potential µ. Here µ = iµ I is purely imaginary. For S g , we use with P µν , R µν , and T µνσ standing for 1/3 of the real part of the trace of 1×1, 1×2 planar Wilson loops and 1×1×1 "parallelogram" loops, respectivley, The coefficents C P , C R , and C T are tadpole improved [43], with u 0 = (< P µν >) 3/4 . The HISQ action with pseudofermion field Φ is where the form of The Dirac operator / D is constructed from smeared links [43]. The fundamental gauge links are U µ (x), the fat links after a level-1 fat7 smearing are V µ (x), the reunitarized links are W µ (x), and the fat links after level-2 asqtad smearing are X µ (x). For simplicity, we use The staggered fermion phases are absorbed into the link variables.ρ and4 are the unit vectors along ρ−direction and 4 direction, respectively.
In the study to address the chiral transition, it is natural to choose the chiral condensate as our observable. The chiral condensate is defined as: N s and N t are the spatial and temporal extent of lattice, respectively. To simplify notation, we use X to represent the chiral condensate. The susceptibility of chiral condensate is defined as We also calculate the Binder cumulant of chiral con-densate which is defined as: The Binder cumulant of the chiral condensate can be expanded around aµ c as [16,17]

III. MC SIMULATION RESULTS
Before presenting the simulation results, we describe the computation details. Simulations are carried out at quark masses am = 0.020, 0.018, 0.015, 0.013, 0.010. The Rational Monte Carlo algorithm [44][45][46] is used to generate configurations. We use different molecular dynamics step sizes for the gauge and fermion parts of the   However, from Figs. 4,5, and 6, we can find that the values of B 4 on different lattice sizes approximately intersect at B 4 (am, aµ c ) ∼ 3. We think that it is because of large finite lattice effects. To gain some understanding about the result, we fit expression to the data on the N s = 20 lattice to get the critical aµ c . The results are presented in Table II. From the results in Table II, we can see that the critical aµ c 's on lattice 20 3 × 4 are approximately in a reasonable region, which should be smaller than 0.262 on the N t = 4 lattice. If we use Eq. (15) to fit the data on the smaller lattice, it can be found that the critical aµ c 's are much larger than 0.262. Moreover, the r-square values in Table II that are close to 1 show that the fit is good. All these facts imply that the smaller lattices have significant finite volume effects.
From Fig. 8 in Ref. [16], we can see that the value of B 4 at the intersection point is 1.604, which is consistent with the Z(2) universality class value. This shows that the finite lattice volume effects in Ref. [16] are very small.

IV. DISCUSSIONS
We have made a simulation in an attempt to understand the nature of the phase transition of QCD with  (13) (14) 0.100 5.968 (40) (17)   two massless flavors with the one quark-loop Symanzikimproved gauge action and the HISQ fermion action by using the method proposed in Ref. [16] at the quark masses am = 0.020, 0.018, 0.015, 0.013, 0.010. In our simulation, we found that the Binder cumulants of the chiral condensate on different lattice volumes at one quark mass intersect at one point. The value of B 4 at the intersection point B 4 (am, aµ c ) was renormalizationinvariant. At the quark masses we used , the value of B 4 at the intersection point B 4 (am, aµ c ) was around 3.
At a nonvanishing quark mass, an additive and multiplicative renormalization ofψψ was needed to define the order parameterψψ when the scaling property was under consideration [8,31,51]. Equations. (9) and (12) in Ref. [8], Eq. (36) in Ref. [51], and Eq. (30) in Ref. [31] were used to subtract the finite quark mass influence on ψψ. However, if we start from Eqs. (12) and (13) and subtract the finite quark mass influence from the chiral condensate, then put the subtracted chiral condensate into Eqs. (12) and (13), we think that multiplicative or additive renormalization ofψψ would have no effect on the value of B 4 (am, aµ c ).
If we can detect the Z(2) transition line, the value of B 4 at the intersection point should be 1.604 [16,17]. However, in our simulation, B 4 (am, aµ c ) ∼ 3 deviated from the Z(2) universality class value. If we just used Eq. (15) to fit the data on N s = 20 lattice, we found the fit was good and the aµ c 's obtained were reasonable. So, we think that B 4 (am, aµ c ) ∼ 3 is because of large finite volume effects as described in Sec. III.
Similar behavior was observed in Ref. [52], in which Wilson-type fermions were employed to determine the critical point separating the crossover from the first-order phase transition region for three-flavour QCD. In that research, the value of kurtosis of the chiral condensate at the intersection point deviated from the universality class value on the N t = 8, 10 lattice due to finite volume correction. This observation indicates that simulation with HISQ action along this direction on the N s > 20 lattice is of great importance.