Chiral phase structure of three flavor QCD in a background magnetic field
Heng-Tong Ding, Christian Schmidt, Akio Tomiya, Xiao-Dan Wang
AAugust 19, 2020
Chiral phase structure of three flavor QCDin a background magnetic field
Heng-Tong Ding, ∗ Christian Schmidt, † Akio Tomiya, ‡ and Xiao-Dan Wang § Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Fakult¨at f¨ur Physik, Universit¨at Bielefeld, D-33615 Bielefeld, Germany RIKEN/BNL Research center, Brookhaven National Laboratory, Upton, NY, 11973, USA
Abstract
We investigate the chiral phase structure of three flavor QCD in a background U (1) magneticfield using the standard staggered action and the Wilson plaquette gauge action. We performsimulations on lattices with a temporal extent of N τ = 4 and four spatial extents of N σ = 8 , , am = 0 .
030 with corresponding pion massestimated as m π ∼
280 MeV such that there exists a crossover transition at vanishing magneticfields, and adopt two values of magnetic field strength in lattice spacing a √ eB (cid:39) . eB/m π ∼
11 and 20, respectively. We find that the transition becomes strongerin the presence of a background magnetic field, and turns into a first order as seen from the volumescaling of the order parameter susceptibility as well as the metastable states in the time historyof the chiral condensate. On the other hand, the chiral condensate and transition temperaturealways increase with B even within the regime of a first order phase transition. This suggests thatthe discrepancy in the behavior of chiral condensates and transition temperature as a function of B between earlier lattice studies using larger-than-physical pion masses with standard staggeredfermions and those using physical pions with improved staggered fermions is mainly due to latticecutoff effects. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] a r X i v : . [ h e p - l a t ] A ug . INTRODUCTION Strong magnetic fields have been shown to have many significant impacts on the propertiesof systems governed by strong interaction, and they may have observable consequences inheavy-ion collision experiments as well as magnetized neutron stars [1, 2]. One of theinteresting impacts is the behavior of transition temperature of strong-interaction systemand chiral condensate as a function of the magnetic field strength B . Early lattice studiesof N f = 2 QCD with standard staggered fermions and larger-than-physical pions on N τ = 4lattices found the so-called magnetic catalyses, which means that the chiral condensateincreases monotonically with increasing B [3]. As the chiral symmetry breaking is enhancedit is expected and observed that the transition temperature consequently increases with B [3].However, based on continuum-extrapolated lattice results of N f = 2+1 QCD using improvedstaggered fermions and physical pions [4] it turns out that the transition temperature hasthe opposite behavior in magnetic field as compared to earlier studies in Ref. [3]. I.e. thetransition temperature decreases with increasing B . The chiral condensate, on the otherhand, first increases and then decreases with B in the transition regime [4]. This non-monotonic behavior of chiral condensate in B , which is called inverse magnetic catalysis,has also been observed in further lattice studies [5–7].The discrepancy of lattice results in the behavior of chiral condensates and transitiontemperature in B reported in Refs [3] and [4–7] is probably due to either large quark massesor possible large cutoff effects present in the first study. To understand the role of quarkmasses in the intricate relation between the (inverse) magnetic catalysis and the reductionof transition temperature, authors of Refs. [3] and [4] have performed lattice studies of N f =2+1 [8] and N f =3 QCD [9] using improved staggered fermions with various larger-than-physical values of pions (370 MeV (cid:46) m π (cid:46)
700 MeV) on N τ =6 lattices. It is foundthat the reduction of transition temperature always holds, however, the inverse magneticcatalysis disappears at a certain value of pion mass. It is suggested in Ref. [8] that theinverse magnetic catalysis is more like a deconfinement catalysis [10] as it is not necessarilyassociated with the reduction of the transition temperature as a function of B .Despite the discrepancy mentioned above the strength of the QCD transition always be-comes larger in the stronger magnetic field as presented in Refs. [3, 11]. And it is speculatedthat the strength could be sufficiently large to turn the crossover transition into a first order2hase transition which can then generate a critical point in the QCD phase diagram in theplane of temperature and magnetic field [11, 12]. However, no true phase transition of QCDin a background magnetic field has been observed in lattice QCD simulations.One of the main motivations of the current study is to explore the chiral phase structureof N f = 3 QCD in a background magnetic field. At the vanishing magnetic field the truefirst order phase transition is not yet observed, and state-of-the-art estimates on the criticalpion mass m cπ based on lattice QCD simulations are m cπ (cid:46)
50 MeV using improved staggeredfermions [13, 14] and m cπ (cid:46)
110 MeV using clover-improved Wilson fermions [15, 16]. Sincethe background magnetic field always enhances the strength of the transition one may wonderwhether it could enlarge the first order chiral phase transition region in N f =3 QCD, i.e.having a larger value of the critical pion mass.In this paper we investigate the transition of N f = 3 QCD in background magneticfields with a quark mass corresponding to pion mass estimated as ∼
280 MeV at vanishingmagnetic field. In our lattice simulations we use standard staggered fermions for a testbedtowards to probe a first-order phase transition with magnetic fields using improved fermions,e.g. Highly Improved Staggered Quarks [17]. The usage of standard staggered fermions witha small quark mass also renders us to understand whether the discrepancy in the behaviorof chiral condensate and transition temperature as a function of the magnetic field strengthin [3, 11] is ascribed to the lattice cutoff effects.The paper is organized as follows. In Section II, we introduce our lattice setup andobservables to be investigated. In Section III, we mainly discuss the order of the transitionbased on results of chiral condensates, Polyakov loops as well as their susceptibilities andBinder cumulants. In Section IV, we conclude our work. The preliminary results have beenreported in [18].
II. LATTICE SETUP AND OBSERVABLES
We perform our simulations on N σ × N τ lattices with 3 mass-degenerate flavors of standardstaggered quarks and the Wilson plaquette gauge action by employing the rational hybridMonte Carlo algorithm [19]. Simulations are performed on lattices with a fixed value oftemporal extent N τ = 4 and four different values of spatial size N σ = 8 , ,
20 and 24. Sincethe critical quark mass, where the first order chiral phase transition starts for this lattice3etup is am c =0.027 [20], we choose the value of quark mass in lattice spacing am to be 0.03in our simulations such that the QCD transition is a crossover at vanishing magnetic field.The temperature, T = 1 / ( aN τ ) is varied through the relation between the lattice spacing a and the inverse gauge coupling β , and specifically temperature increases with the value of β . The background U (1) magnetic field is implemented in a conventional way [21] and willbe reviewed in the following subsection II A. The relevant observables will be introduced insubsection II B. A. Magnetic fields on the lattice
Magnetic fields couple to quarks with their electric charges and through covariant deriva-tive in the continuum, D µ = ∂ µ − igA µ + iqa µ , (1)where g is the SU(3) gauge coupling and q is the electric charge of a quark. A µ and a µ aregauge fields for SU(3) and U(1), respectively. On the lattice the background U (1) magneticfield is introduced by substituting the SU(3) link U µ by its product with the U(1) link u µ inthe lattice Dirac operator, D [ U ] → D [ uU ] . (2)In our simulations the magnetic field only points along the z -direction. Since the x – y planehas boundaries for a finite system size, appropriate boundary conditions need to be imposed.Besides, the magnetic field is realized as a U(1) plaquette, and it introduces non-trivialconditions to the magnitude of the magnetic field as will be depicted next.Let us denote the lattice size ( N x , N y , N z , N t ) and coordinate as n µ = 0 , · · · , N µ − µ = x, y, z, t ). The background magnetic field pointing along the z -direction (cid:126)B = (0 , , B )is described by the link variable u µ ( n ) of the U(1) field, and u µ ( n ) is expressed as follows inthe Landau gauge [4, 21], u x ( n x , n y , n z , n t ) = exp[ − iqa BN x n y ] ( n x = N x − u y ( n x , n y , n z , n t ) = exp[ iqa Bn x ] , (3) u z ( n x , n y , n z , n t ) = u t ( n x , n y , n z , n t ) = 1 . x -direction. One-valuedness of the one-particle wave function along with a plaquette requiresthe following quantization, a qB = 2 πN b N x N y , (4)where N b ∈ Z is the number of magnetic flux through the unit area for the x - y plane. Theultraviolet cutoff a also introduces a periodicity of the magnetic field along with N b . Namely,a range 0 ≤ N b < N x N y /
4, represents an independent magnitude of the magnetic field B .In our simulations the 3 mass-degenerate flavors are of up, down and strange quark type,and thus two of the flavors have electrical charge of q d,s = − e and the rest one has q u = e with e the electric charge of electron. Thus to satisfy the quantization condition Eq. 4 wetake the electric charge q to be that of down quark type, i.e. q = q d,s = − e . To simplify thenotation we use ˆ b expressed as follows to denote the magnetic field strengthˆ b ≡ a √ eB = (cid:115) πN b N x N y . (5)We choose certain values of N b such that ˆ b are the same in physical units among variouslattice sizes which are listed in Table I. N b b N σ = 8 N b b N σ = 16 N b b N σ = 20 N b b N σ = 24TABLE I. Number of magnetic flux N b used in our simulations and the approximated values ofˆ b ≡ a √ eB for N σ = 8, 16, 20 and 24. Based on an earlier estimate of the pion mass in N f = 3 QCD [25] we make a estimateof the scale at the pseudo-critical temperature T pc := T pc (ˆ b = 0). We find m π ∼
280 MeVand m π /T pc (cid:39) .
77 at ˆ b = 0. The values of magnetic field strength for ˆ b = 1 . eB/ ( T pc ) (cid:39)
36 and 64, i.e. eB ∼ m π and eB ∼ m π , respectively.We further note that the temperature in the explored β -intervall [5.13,5.19] may change byapproximately 25 MeV, which would lead to corresponding changes in the magnetic field.Our simulation parameters and statistics are summarized in Tab. III - VI. It is worthmentioning that for our largest lattices, i.e. with N σ = 24, we have performed simulationswith a hot start (configuration with random elements) and a cold start (configuration with5nit elements) to check metastable states around the transition temperature to overcomesmall tunnelling rates. B. Observables
An expectation value for an operator O for given gauge coupling β , mass and magneticflux N b is defined by, (cid:104) O (cid:105) β,m,N b = (cid:90) D U P β,m,N b [ U ] O [ U ] , (6) P β,m,N b [ U ] = 1 Z β,m,N b det (cid:2) D [ U, N b ] + m (cid:3) / det (cid:2) D -1/3 [ U, N b ] + m (cid:3) / e − S gauge [ U ; β ] , (7)where D is the standard staggered Dirac operator and its subscript indicates the electriccharge of the related quark. The fractional power, e.g. 1/4 and 2/4, is due to the fourthroot of staggered fermions in our simulations.We measure the chiral condensate for a down type quark , (cid:10) ψψ (cid:11) = 14 N σ N τ (cid:28) Tr 1 D − / + m (cid:29) , (8)where Tr [ · · · ] is a trace over color and space-time, and a factor of 4 in the denominatoradjusts for taste degrees of freedom. Its disconnected susceptibility is given by χ disc = 116 N σ N τ (cid:34)(cid:42)(cid:18) Tr 1 D − / + m (cid:19) (cid:43) − (cid:28) Tr 1 D − / + m (cid:29) (cid:35) . (9)While chiral condensate and its susceptibility are related to chiral symmetry we alsocompute the Polyakov loop which is related to the deconfinement transition in a pure gluesystem P = 1 V (cid:42)(cid:88) (cid:126)x Tr (cid:89) t U ( t, (cid:126)x ) (cid:43) , (10)and its susceptibility χ P .We will also analyse the Binder cumulants B M of order parameters M , e.g. the chiralcondensate[22] and the Polyakov loop as well as the constructed order parameter from amixture of the chiral condensate and the gauge action (cf. Eq. 13). B M is defined as follows B M ( β, N b ) = (cid:104) ( δM ) (cid:105)(cid:104) ( δM ) (cid:105) , (11) In principle, we can average up, down and strange quark condensate, or up and down. However, wechoose down condensate to avoid such arbitrariness. δM = M − (cid:104) M (cid:105) gives the deviation of M from its mean value on a given gauge fieldconfiguration. From different distributions of M in different phases the value of B M can beobtained and used to distinguish phase transitions. Given that the chiral condensate is theorder parameter of the phase transition , for a first order phase transition, B ψψ = 1 [23];for a crossover B ψψ (cid:39) B ψψ (cid:39) . N eff given by N eff = τ int , (12)where τ int is the integrated autocorrelation time for the chiral condensate. Obtained resultsof τ int are listed in Appendix A. Hereafter we will show our results of various quantities ina dimensionless way, i.e. all in units of lattice spacing a . III. RESULTS
In the first subsection, we discuss the observables obtained from a fixed volume N σ = 16and the history of the chiral condensate at vanishing and nonzero values of ˆ b . In the secondsubsection, we study the volume dependences of chiral observables and the Polyakov loopas well as their susceptibilities at ˆ b = 0 and 1 . N σ =8, 16, 20and 24. In the third subsection, we show results based on a appropriate order parameterconstructed from a combination of the chiral condensate and the gluon action. A. a √ eB dependence of observables obtained on N σ = 16 lattices We show the chiral condensate and its disconnected susceptibility for ˆ b = 0 , . N σ =16 as a function of the inverse gauge coupling β in the leftand right plot of Fig. 1, respectively. We recall that the temperature is an increasing functionof β . As seen from the left plot the value of the chiral condensate is enhanced by the magneticfield in the whole temperature region. Namely only the magnetic catalysis is observed. One The same holds true for other order parameters. The autocorrelation time and its error are given in Appendix A where τ int are rounded in integer forsimplicity. β value where the chiral condensate drops most rapidly becomeslarger as the magnetic field strength increases. This indicates that the transition temperatureincreases as the magnetic field becomes stronger. Meanwhile as the strength of the magneticfield increases the dropping of chiral condensates becomes more rapidly. This means thatthe transition becomes stronger with stronger magnetic fields. These two observations arealso visible in the behaviour of disconnected chiral susceptibilities as shown in the right plotof Fig. 1. I.e. the peak location of disconnected chiral susceptibility shifts to a larger valueof β and the peak height of disconnected chiral susceptibility increases as the strength ofmagnetic field increases. . . . . . . . . . . .
28 5 .
13 5 .
14 5 .
15 5 .
16 5 .
17 5 .
18 5 .
19 5 . .
21 5 . h ψ ¯ ψ i β ˆ b =0ˆ b =1.5ˆ b =2.0 .
13 5 .
14 5 .
15 5 .
16 5 .
17 5 .
18 5 .
19 5 . .
21 5 . χ disc β ˆ b =0ˆ b =1.5ˆ b =2.0 FIG. 1. The chiral condensate (left) and disconnected chiral susceptibility (right) at ˆ b = 0 , N σ = 16 as a function of β . . . . . . . . . . .
22 5 .
13 5 .
14 5 .
15 5 .
16 5 .
17 5 .
18 5 .
19 5 . .
21 5 . h P i β ˆ b =0ˆ b =1.5ˆ b =2.0 .
13 5 .
14 5 .
15 5 .
16 5 .
17 5 .
18 5 .
19 5 . .
21 5 . χ P β ˆ b =0ˆ b =1.5ˆ b =2.0 FIG. 2. Same as Fig. 1 but for the Polyakov loop and its susceptibility.
8e show similar plots as Fig. 1 for the Polyakov loop and its susceptibility at ˆ b = 0 , . β , i.e. higher transitiontemperature with stronger magnetic field. This is similar to the observation from chiralobservables. Besides that transition temperatures signalled by the Polyakov loop and itssusceptibility are close to those obtained from chiral condensates and disconnected chiralsusceptibility. What’s more, it is interesting to see that the peak height of the Polyakovloop susceptibility also becomes higher in a stronger magnetic field, although the Polyakovloop is an order parameter for the confinement-deconfinement phase transition of a pure gluesystem while the chiral condensate is the order parameter of QCD transitions with vanishingquark massless. This could be due to the fact that neither of these two quantities are thetrue order parameter but a part of the true order parameter in N f = 3 QCD [14, 25].At vanishing magnetic field the transition is a crossover, and as the magnetic field strengthbecomes larger the QCD transition becomes stronger. In particular the jumping behaviorof the chiral condensate and the Polyakov loop at ˆ b =2 is seen from left plots of Fig. 1 andFig. 2. This might indicate a first order phase transition. We investigate the nature of thetransition by further looking into the time history of chiral condensates. We show in Fig. 3the time history of chiral condensates near the transition temperature at ˆ b = 0 (Top left),ˆ b = 1 . b = 2 (Bottom left) on lattices with N σ = 16. As expected that atvanishing magnetic field (ˆ b =0) there does not exist any metastable behavior, while in thecase of ˆ b =1.5 and 2 the metastable behavior becomes obvious. To confirm the metastablebehavior seen from the volume of N σ =16, we also study the case of ˆ b =1.5 with a largervolume, i.e. on 24 × Note that the peak height of Polyakov loop susceptibility also increases as the system approaches thefirst order phase transition region with smaller values of quark masses for the case of zero magnetic fieldstrength ˆ b = 0 [25]. . . . .
25 5000 15000 25000 35000 h ψ ¯ ψ i MC time β =5.1400 β =5.1438 β =5.1475 . . . . . . .
35 20000 50000 80000 110000 140000 h ψ ¯ ψ i MC time β =5.1660 β =5.7000 β =5.7500 . . . . . . . . h ψ ¯ ψ i MC time β =5.1780 β =5.1826 β =5.1828 β =5.1850 . . . . . . h ψ ¯ ψ i . . . . . .
35 10000 30000 50000 70000 90000 h ψ ¯ ψ i β = 5 . , hot β = 5 . , hot β = 5 . , hot β = 5 . , hotMC time coldcoldcoldcold FIG. 3. Monte Carlo time history of chiral condensates. The x axes are in units of trajectoriesof molecular dynamics. Top left: with ˆ b = 0 and N σ = 16. Top right: with ˆ b = 1 . N σ = 16.Bottom left: with ˆ b = 2 and N σ = 16. Bottom right: with ˆ b = 1 . N σ = 24. transition. The former is observed for the case of β = 5 . β = 5 . b ≥ B. Volume dependences of observables with ˆ b = 0 and . In this subsection, we discuss more on the volume dependence of observables at ˆ b = 0and 1 . b ≥ .
5. InFig. 4 we show chiral condensates and their disconnected susceptibilities as a function of β obtained from lattices with four different spatial sizes at ˆ b = 0 (Top two plots) and ˆ b =1.5(Bottom two plots). At vanishing magnetic field chiral condensates obtained from N σ = 16 , . . . . . . . . .
125 5 .
13 5 .
135 5 .
14 5 .
145 5 .
15 5 .
155 5 .
16 5 . h ψ ¯ ψ i β N σ =8N σ =16N σ =20N σ =24 051015202530355 .
125 5 .
13 5 .
135 5 .
14 5 .
145 5 .
15 5 .
155 5 .
16 5 . χ disc β N σ =8N σ =16N σ =20N σ =240 . . . . . . . . . .
155 5 .
16 5 .
165 5 .
17 5 .
175 5 .
18 5 . h ψ ¯ ψ i β N σ =8N σ =16N σ =20N σ =24 0204060801001205 .
155 5 .
16 5 .
165 5 .
17 5 .
175 5 .
18 5 . χ disc β N σ =8N σ =16N σ =20N σ =24 FIG. 4. Volume dependences of chiral condensates (left column) and its disconnected suscepti-bilities (right column) at vanishing magnetic field (top two plots) and ˆ b =1.5 (bottom two plots).
20 and 24 are consistent among each other, while large deviations are seen for the resultsobtained from N σ =8. In the case of ˆ b =1.5 the finite size effect seems to appear at β smallerand close to β c starting with a larger volume, i.e. N σ =16. This could be due to the factthat the system with the presence of a magnetic field tends to have a stronger transition andconsequently more statistics is needed to get robust results , and that the correlation lengthin the system becomes longer in the proximity of the true phase transition. Nevertheless,the point where chiral condensates drops most rapidly are consistent among various volumesfor both vanishing and nonzero magnetic fields. The is also reflected in the peak locations ofdisconnected susceptibilities shown in the right column of Fig. 4. In the case of a first orderphase transition the disconnected chiral susceptibility should grow linearly in volume. It ismore or less the case for disconnected susceptibilities at ˆ b =1.5 as seen from bottom right Although we generated two times more configurations at ˆ b =1.5, the effective configurations, however, istwo times less than that at ˆ b =0. . . . . .
125 5 .
13 5 .
135 5 .
14 5 .
145 5 .
15 5 .
155 5 .
16 5 . B ψ ¯ ψ β N σ =8N σ =16N σ =20N σ =24crossover1st order2nd order . . . . .
16 5 .
165 5 .
17 5 .
175 5 . B ψ ¯ ψ β N σ =8N σ =16N σ =20N σ =24crossover1st order2nd order FIG. 5. Binder cumulants of chiral condensate at zero magnetic field (left) and nonzero magneticfield of ˆ b =1.5 (right). We show Binder cumulants of chiral condensates at ˆ b = 0 and 1.5 in the left and rightplot of Fig. 5, respectively. One can see that at the critical β where χ disc peaks the Bindercumulant reaches to its minimum, which is almost independent of volume at ˆ b = 0 and onlystarts to saturate with N σ ≥
20 at ˆ b =1.5. In nonzero magnetic fields the minimum valuesof Binder cumulant become smaller, i.e. shifting from about 1.6 in the vanishing magneticfield to about 1. This indicates that the transition becomes stronger and tends to becomea first order phase transition region at ˆ b =1.5. This is consistent with what we found fromchiral condensates and its susceptibilities.We show similar figures as Fig. 4 and Fig. 5 for Polyakov loops in Fig. 6 and Fig. 7. Theobservation is similar to that from chiral observables. The volume dependences of Polyakovloops and corresponding susceptibilities are similar to those of chiral observables. Besides,the value for the Binder cumulant has the same tendency. The peak heights of susceptibilitiesand the minimum values of Binder cumulants for the chiral condensate and the Polyakovloop are summarized in Tab. II. 12 . . . . . . . .
18 5 .
13 5 .
135 5 .
14 5 .
145 5 .
15 5 .
155 5 . h P i β N σ =8N σ =16N σ =20N σ =24 0510152025303540 5 .
13 5 .
135 5 .
14 5 .
145 5 .
15 5 .
155 5 . χ P β N σ =8N σ =16N σ =20N σ =240 . . . . . . . . .
16 5 .
165 5 .
17 5 .
175 5 . h P i β N σ =8N σ =16N σ =20N σ =24 020406080100120140 5 .
16 5 .
165 5 .
17 5 .
175 5 . χ P β N σ =8N σ =16N σ =20N σ =24 FIG. 6. Same as Fig. 4 but for the Polyakov loop and its susceptibility. . . . . .
125 5 .
13 5 .
135 5 .
14 5 .
145 5 .
15 5 .
155 5 .
16 5 . B P β N σ =8N σ =16N σ =20N σ =24crossover1st order2nd order 00 . . . . .
16 5 .
165 5 .
17 5 .
175 5 . B P β N σ =8N σ =16N σ =20N σ =24crossover1st order2nd order FIG. 7. Same as Fig. 5 but for the Polyakov loop.
C. Binder cumulants and disconnected susceptibilities of the order parameter
In the vicinity of a critical point in 3-flavor QCD, the chiral condensate itself is nota true order parameter, but is part of a mixture of operators that defines the true orderparameter M [14, 25]. In three flavor QCD, such an order parameter can be constructed as13 σ a √ eB β c χ disc ( β c ) B ψψ ( β c ) χ P ( β c ) B P ( β c )8 0 5.1450 4.1(1) 1.65(3) 5.1(1) 1.87(3)16 0 5.1438 13.0(9) 1.90(9) 15(1) 2.03(9)20 0 5.1435 21(2) 1.7(1) 25(2) 1.7(1)24 0 5.1438 29(4) 1.88(4) 33(4) 1.96(2)8 1.5 5.1710 6.2(2) 1.57(3) 7.1(2) 1.70(3)16 1.5 5.1690 34(3) 1.62(5) 38(4) 1.57(4)20 1.5 5.1700 68(7) 1.24(4) 75(8) 1.26(4)24 1.5 5.1699 105(11) 1.15(2) 115(13) 1.173(3)TABLE II. Values of disconnected chiral susceptibility χ disc , Polyakov loop susceptibility χ P , andthe Binder cumulants for the chiral condensate and the Polyakov loop at ˆ b = 0 and 1.5. β c is thevalue of β where χ disc ( β ) peaks. a combination of the plaquette action and the chiral condensate as follows [14, 25] M ( β, s ) = ψψ ( β ) + s N σ N τ S gauge ( β ) , (13)where S gauge = 6 N σ N τ ˜ P and ˜ P is the plaquette. Correspondingly, its susceptibility is, χ mixed = 116 N σ N τ (cid:2)(cid:10) ( M ( β, s )) (cid:11) − (cid:104) M ( β, s ) (cid:105) (cid:3) . (14)In this work, we do not intend to determine the mixing parameter s and just use the valueobtained in the previous work for a √ eB = 0 in [20], i.e. s = − .
8. As will be seen next theresults obtained using s = − . s = 0 qualitatively (cf. Table II).We show in the left panel of Fig. 8 the order parameter susceptibility χ mixed divided bythe spatial volume at ˆ b = 0 and 1.5. If the phase transition is of a first order χ mixed dividedby the spatial volume should be a constant in volume. It is clearly seen that at ˆ b =1.5 datapoints obtained from lattices with different volumes all overlap at about 0.015 which holdstrue also towards the infinite volume limit of ( N τ /N σ ) →
0, while it is not the case forˆ b = 0. The Binder cumulant of the order parameter is shown in the right panel of Fig. 8.Data points for ˆ b = 0 are close the Z line, so it seems to belong to a weak second order orcrossover transition in the thermodynamic limit. On the other hand, data points for ˆ b = 1 . N τ /N σ ) →
0. This is again consistent withthe behavior in the first-order phase transition.14 . . . . . . . . .
018 0 0 .
02 0 .
04 0 .
06 0 .
08 0 . .
12 0 . χ mixed /N σ (N τ / N σ ) ˆb=0ˆb=1.5 11 . .
53 0 0 .
02 0 .
04 0 .
06 0 .
08 0 . . minmixed (N τ / N σ ) ˆb = 0ˆb=1.5crossover1st order2nd order FIG. 8. Left: χ mixed divided by spatial volume (left) and Binder cumulant of M at β = β c obtainedwith s = − . b = 0 and 1.5 as a function of the inverse spatial volume represented by ( N τ /N σ ) . In summary, we conclude that the system with am = 0 .
03 with magnetic fields a √ eB ≥ eB (cid:38) m π ) exists a first order phase transition. IV. CONCLUSION
We have investigated the chiral phase structure of three flavor QCD in the magnetic field.The simulations are performed with 3-mass degenerate flavors of standard staggered quarkand the Wilson plaquette gauge action on N τ = 4 lattices. We started with simulations ata fixed quark mass of am = 0 .
03 at vanishing magnetic field where a crossover transition isobserved. After turning on the magnetic field we studied dependences of chiral observables,i.e. chiral condensates and corresponding susceptibilities as well as Polyakov loops on themagnetic field strength. We found that chiral condensates increase with the magnetic fieldstrength, namely magnetic catalyses in the whole temperature window. And transitiontemperatures determined from both chiral condensates and Ployakov loops as well as theirsusceptibilities increase with increasing magnetic field strength. From the dropping behaviorof chiral condensates and Polyakov loops near the transition temperature and the peakheights of susceptibilities we found that the strength of transition increases with increasingmagnetic field strength. We thus checked the time history of the chiral condensate, andthe metastable behavior of chiral condensates is observed at a √ eB ≥ .
5. This indicates afirst order phase transition occurring to systems at a √ eB ≥ .
5. We further investigate a15ore appropriate order parameter which is constructed from the chiral condensate and theplaquette gauge action. We studied the volume dependence of order parameter susceptibilityand the Binder cumulant of the order parameter, and confirmed that there exists a first orderphase transition at a √ eB ≥ .
5. Based on the earlier estimate we find that a √ eB ≥ . eB/m π (cid:38)
11 with the pion mass at zero magnetic field m π ∼
280 MeV.We do not find any signs of inverse magnetic catalysis and the reduction of transitiontemperature as a function of the magnetic field strength, and this is consistent with theresults obtained from lattice studies of N f = 2 QCD with the standard staggered quarks andlarger-than-physical pions [3, 26]. Since our findings of magnetic catalysis and increasing oftransition temperature with the magnetic field strength even holds in the regime occurring afirst order phase transition this suggests that the discrepancy from results using the improvedstaggered fermions with physical pions is mainly due to the lattice cutoff effects.It is worth recalling that in the case of lattice studies using improved staggered fermionsthe strength of transition also increases with increasing magnetic field strength and a firstorder phase transition has not yet been observed in such studies. Although we find a firstorder phase transition in the current study using the standard staggered fermions, we do notintend to provide a precise determination of a critical magnetic field in which the transitionturns into a first/second order phase transition in the current discretization scheme. Wewill leave it for future studies using improved staggered fermions, such as HISQ fermionsto achieve more realistic results on the chiral phase structure of N f = 3 QCD [17]. As incurrent studies severe critical slowing down is expected for simulations in the vicinity of phasetransitions with larger volumes (see Appendix A), the autocorrelation length will probablybecome even longer in simulations with smaller lattice spacings or improved actions, in whichthe multicanonical method [27, 28] or other extended Monte Carlo method reviewed in [29]might help. ACKNOWLEDGEMENT
We thank Swagato Mukherjee for the early involvement of the work as well as enlighten-ing discussions. This work was supported by the NSFC under grant no. 11535012 and no.11775096, RIKEN Special Postdoctoral Researcher program, Deutsche Forschungsgemein-schaft project number 315477589 – TRR 211 and European Union H2020-MSCA-ITN-2018-1613942 (EuroPLEx). The numerical simulations have been performed on the GPU clusterin the Nuclear Science Computing Center at Central China Normal University (NSC ). Appendix A: Autocorrelation time
Here we explain the autocorrelation time to make this paper self-contained. A sequence ofconfigurations are affected by the autocorrelation, which is evaluated by the autocorrelationfunction. However, the autocorrelation function itself is a statistical object, so we cannotdetermine it exactly. Instead we calculate the approximated autocorrelation function [30, 31]defined by, Γ( τ ) = 1 N trj − τ N trj (cid:88) c ( O c − O )( O c + τ − O ) , (A1)where O c = O [ U ( c ) ] is the value of operator O for the c -th configuration U ( c ) and τ is thefictitious time of HMC. N trj is the number of trajectories. Conventionally, the normalizedautocorrelation function ρ ( τ ) = Γ( τ ) / Γ(0) is used.The integrated autocorrelation time τ int approximately quantifies effects of autocorrela-tion. This is given by, τ int = 12 + W (cid:88) τ =1 ρ ( τ ) . (A2)We can regard two configurations separated by 2 τ int as independent. In practice, we de-termine a window size W as a first point W = τ , where Γ( τ ) < τ . Thestatistical error of integrated autocorrelation time is estimated by the Madras–Sokal formula[31, 32], (cid:10) δτ (cid:11) (cid:39) W + 2 N trj τ , (A3)and we use the square root of (A3) for an estimate of the error on the autocorrelation time, (cid:112) (cid:104) δτ (cid:105) ≡ ∆ τ int . Appendix B: Critical slowing down
The estimated autocorrelation length for the chiral condensate and the number of in-dependent configurations are listed in Tab. IV and Tab. VI. It can be observed that the17 A u t o c o rr e l a t i o n l e n g t h N σ a √ eB =0.0 a √ eB =1.5 c N z σ , z =1.94 c N z σ , z =2.50 FIG. 9. The largest autocorrelation as a function of the spatial size N σ at ˆ b ≡ a √ eB = 0 and 1.5. autocorrelation time becomes longer in the presence of a magnetic field. This is due to thefact that the system is close to the critical region as discussed in the main text.In Fig. 9, we estimate the dynamic critical exponent of the system by using the longestautocorrelation length in each volume. A fit ansatz of log τ int ( N σ ) = z log N σ + c with fitparameters z and c is adopted and this fit ansatz is the same form as τ int ( N σ ) ∼ N zσ , whichis the definition of the dynamic critical exponent. One can see that the system is affectedby severe critical slowing down at ˆ b ≡ √ a eB = 1 .
5. It is thus practically challengingto investigate the first-order phase transition by using conventional direct simulations withHMC in the thermodynamic limit.
Appendix C: Summary of statistics
We summarize our simulation parameters and statistics in Tab. III - VI. If the statistics issufficient the binning size for the Jackknife analysis is taken as the Bin-size (cid:38) τ int , otherwisethe Bin-size is adjusted such that it is the divisor of number of trajectories by around 10.18n the latter case, the error might be underestimated. N σ β N b Bin-size τ int ∆ τ int Traj. N eff N σ β N b Bin-size τ int ∆ τ int Traj. N eff N σ = 8. Trajectories for thermalization are already discarded. Here onetrajectory denotes one time unit in the molecular dynamics.[1] D. E. Kharzeev, K. Landsteiner, A. Schmitt, and H.-U. Yee, Lect. Notes Phys. , 1 (2013),arXiv:1211.6245 [hep-ph].[2] V. A. Miransky and I. A. Shovkovy, Phys. Rept. , 1 (2015), arXiv:1503.00732 [hep-ph].[3] M. D’Elia, S. Mukherjee, and F. Sanfilippo, Phys. Rev. D82 , 051501 (2010), arXiv:1005.5365[hep-lat].[4] G. S. Bali, F. Bruckmann, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg, A. Schafer, and K. K.Szabo, JHEP , 044 (2012), arXiv:1111.4956 [hep-lat].[5] E. M. Ilgenfritz, M. Muller-Preussker, B. Petersson, and A. Schreiber, Phys. Rev. D89 ,054512 (2014), arXiv:1310.7876 [hep-lat]. σ β N b Bin-size τ int ∆ τ int Traj. N eff
16 5.1300 00 500 24 4 9500 20116 5.1350 00 500 28 5 9500 17216 5.1400 00 500 193 46 41500 10716 5.1412 00 2000 636 270 38000 3016 5.1425 00 1000 447 150 41000 4616 5.1438 00 1000 396 117 40000 5016 5.1450 00 500 190 47 40500 10716 5.1475 00 500 97 19 37500 19316 5.1500 00 500 42 6 38500 45316 5.1550 00 500 25 5 7000 14216 5.1600 00 500 22 5 7000 163 N σ β N b Bin-size τ int ∆ τ int Traj. N eff
16 5.1600 32 500 40 8 29500 37316 5.1640 32 500 36 4 47500 65916 5.1650 32 500 60 12 36000 29816 5.1660 32 500 55 8 41000 37516 5.1670 32 500 100 25 42500 21216 5.1673 32 500 94 16 97500 51816 5.1675 32 500 136 21 97500 35816 5.1680 32 2000 630 169 102000 8116 5.1685 32 4000 1600 722 96000 30 N σ β N b Bin-size τ int ∆ τ int Traj. N eff
16 5.1690 32 7000 3635 1925 140000 1916 5.1695 32 3000 1406 635 57000 2016 5.1700 32 8000 4481 3263 144000 1616 5.1710 32 500 265 56 97500 18416 5.1720 32 500 92 26 27500 15016 5.1730 32 500 63 13 27500 21916 5.1740 32 500 60 17 27500 23016 5.1750 32 500 38 6 27500 36016 5.1800 32 500 55 16 27500 249 N σ β N b Bin-size τ int ∆ τ int Traj. N eff
16 5.1700 56 100 21 8 1500 3616 5.1750 56 100 26 13 1500 2916 5.1760 56 100 32 7 9000 14016 5.1770 56 100 47 12 9000 9616 5.1780 56 150 53 5 99000 93916 5.1790 56 150 70 9 81450 58016 5.1800 56 500 70 12 45500 32516 5.1810 56 500 195 35 95500 24516 5.1820 56 2500 1166 412 95000 4116 5.1822 56 2000 142 51 42000 14816 5.1824 56 7000 4641 3642 91000 1016 5.1826 56 8000 3890 2340 88000 1116 5.1828 56 9000 4820 3223 90000 916 5.1830 56 8000 3923 2693 72000 916 5.1840 56 100 865 393 49000 2816 5.1850 56 100 59 24 8900 7516 5.1900 56 100 23 8 1900 41TABLE IV. Same as Table III but for N σ = 16. σ β N b Bin-size τ int ∆ τ int Traj. N eff
20 5.1400 00 200 106 18 50800 24020 5.1430 00 500 237 68 29500 6220 5.1435 00 1500 764 353 30000 2020 5.1440 00 1300 682 358 29900 2220 5.1450 00 1000 537 244 30000 2820 5.1460 00 300 184 51 30600 83 N σ β N b Bin-size τ int ∆ τ int Traj. N eff
20 5.1650 50 200 68 9 95200 69620 5.1690 50 5000 8327 7256 90000 520 5.1695 50 3000 1750 1157 39000 1120 5.1698 50 7000 5680 6450 35000 320 5.1700 50 5000 4919 4950 40000 420 5.1710 50 150 291 78 66750 11520 5.1750 50 100 39 10 7800 99TABLE V. Same as Table III but for N σ = 20. N σ β N b Bin-size τ int ∆ τ int Traj. N eff
24 5.1400 00 500 168 129 7500 2224 5.1412 00 500 105 41 6500 3124 5.1425 00 1000 537 280 46900 4424 5.1438 00 4000 2171 1269 78000 1824 5.1450 00 500 176 38 74500 21124 5.1475 00 500 82 34 6500 40 N σ β N b Bin-size τ int ∆ τ int Traj. N eff
24 5.1650 72 500 79 17 29500 18724 5.1680 72 500 102 17 47500 23224 5.1690 72 1000 215 28 180000 41924 5.1695 72 1500 174 26 102000 29324 5.1698 72 4000 4387 3022 76000 924 5.1699 72 2000 2535 2917 14600 324 5.1700 72 15000 8542 5236 186000 1124 5.1710 72 1000 111 25 42500 19124 5.1720 72 500 77 11 69500 45324 5.1730 72 500 111 29 51500 23224 5.1750 72 500 50 6 48000 481TABLE VI. Same as Table III but for N σ = 24.[6] V. G. Bornyakov, P. V. Buividovich, N. Cundy, O. A. Kochetkov, and A. Sch¨afer, Phys. Rev. D90 , 034501 (2014), arXiv:1312.5628 [hep-lat].[7] G. S. Bali, F. Bruckmann, G. Endr¨odi, S. D. Katz, and A. Sch¨afer, JHEP , 177 (2014),arXiv:1406.0269 [hep-lat].[8] M. D’Elia, F. Manigrasso, F. Negro, and F. Sanfilippo, Phys. Rev. D98 , 054509 (2018),arXiv:1808.07008 [hep-lat].
9] G. Endrodi, M. Giordano, S. D. Katz, T. G. Kov´acs, and F. Pittler, JHEP , 007 (2019),arXiv:1904.10296 [hep-lat].[10] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro, A. Rucci, and F. Sanfilippo, Phys.Rev. D94 , 094007 (2016), arXiv:1607.08160 [hep-lat].[11] G. Endrodi, JHEP , 173 (2015), arXiv:1504.08280 [hep-lat].[12] T. D. Cohen and N. Yamamoto, Phys. Rev. D89 , 054029 (2014), arXiv:1310.2234 [hep-ph].[13] G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo,
Proceedings, 25th International Sympo-sium on Lattice field theory (Lattice 2007): Regensburg, Germany, July 30-August 4, 2007 ,PoS
LATTICE2007 , 182 (2007), arXiv:0710.0998 [hep-lat].[14] A. Bazavov, H. T. Ding, P. Hegde, F. Karsch, E. Laermann, S. Mukherjee, P. Petreczky, andC. Schmidt, Phys. Rev.
D95 , 074505 (2017), arXiv:1701.03548 [hep-lat].[15] Y. Nakamura, Y. Kuramashi, H. Ohno, and S. Takeda, in (2019)arXiv:1909.05441 [hep-lat].[16] Y. Kuramashi, Y. Nakamura, H. Ohno, and S. Takeda, Phys. Rev.
D101 , 054509 (2020),arXiv:2001.04398 [hep-lat].[17] A. Tomiya, H.-T. Ding, X.-D. Wang, Y. Zhang, S. Mukherjee, and C. Schmidt,
Proceedings,36th International Symposium on Lattice Field Theory (Lattice 2018): East Lansing, MI,United States, July 22-28, 2018 , PoS
LATTICE2018 , 163 (2019), arXiv:1904.01276 [hep-lat].[18] A. Tomiya, H.-T. Ding, S. Mukherjee, C. Schmidt, and X.-D. Wang,
Proceedings, 35th In-ternational Symposium on Lattice Field Theory (Lattice 2017): Granada, Spain, June 18-24,2017 , EPJ Web Conf. , 07041 (2018), arXiv:1711.02884 [hep-lat].[19] M. A. Clark,
Proceedings, 24th International Symposium on Lattice Field Theory (Lattice2006): Tucson, USA, July 23-28, 2006 , PoS
LAT2006 , 004 (2006), arXiv:hep-lat/0610048[hep-lat].[20] D. Smith and C. Schmidt,
Proceedings, 29th International Symposium on Lattice field theory(Lattice 2011): Squaw Valley, Lake Tahoe, USA, July 10-16, 2011 , PoS
LATTICE2011 , 216(2011), arXiv:1109.6729 [hep-lat].[21] M. Al-Hashimi and U.-J. Wiese, Annals Phys. , 343 (2009), arXiv:0807.0630 [quant-ph].[22] K. Binder, Physical Review Letters , 693 (1981).
23] S. Ejiri, Phys. Rev.
D77 , 014508 (2008), arXiv:0706.3549 [hep-lat].[24] H. W. J. Blote, E. Luijten, and J. R. Heringa, J. Phys.
A28 , 6289 (1995), arXiv:cond-mat/9509016 [cond-mat].[25] F. Karsch, E. Laermann, and C. Schmidt, Phys. Lett.
B520 , 41 (2001), arXiv:hep-lat/0107020[hep-lat].[26] M. D’Elia and F. Negro, Physical Review D , 114028 (2011).[27] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. , 9 (1992), arXiv:hep-lat/9202004 [hep-lat].[28] C. Bonati, M. D’Elia, G. Martinelli, F. Negro, F. Sanfilippo, and A. Todaro, JHEP , 170(2018), arXiv:1807.07954 [hep-lat].[29] Y. Iba, Int. J. Mod. Phys. C12 , 623 (2001), arXiv:cond-mat/0012323 [cond-mat].[30] U. Wolff (ALPHA), Comput. Phys. Commun. , 143 (2004), [Erratum: Comput. Phys.Commun.176,383(2007)], arXiv:hep-lat/0306017 [hep-lat].[31] N. Madras and A. D. Sokal, J. Statist. Phys. , 109 (1988).[32] M. Luscher, Comput. Phys. Commun. , 199 (2005), arXiv:hep-lat/0409106 [hep-lat]., 199 (2005), arXiv:hep-lat/0409106 [hep-lat].