Chiral phase structure of three flavor QCD in a background magnetic field

We investigate the chiral phase structure of three flavor QCD in a background $U(1)$ magnetic field using the standard staggered action and the Wilson plaquette gauge action. We perform simulations on lattices with a temporal extent of $N_\tau=4$ and four spatial extents of $N_\sigma = 8,16, 20$ and 24. We choose a smaller-than-physical quark mass in lattice spacing as $am = 0.030$ such that there exists a crossover transition at vanishing magnetic fields, and adopt two values of magnetic field strength in lattice spacing $a \sqrt{ e{B}}\simeq 1.5$ and 2. We find that the transition becomes stronger in the presence of a background magnetic field, and turns into a first order as seen from the volume scaling of the order parameter susceptibility as well as the metastable states in the time history of the chiral condensate. On the other hand, the chiral condensate and transition temperature always increase with $B$ even within the regime of a first order phase transition. This suggests that the 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 staggered fermions and those using physical pions with improved staggered fermions is mainly due to lattice cutoff effects.


I. INTRODUCTION
Strong magnetic fields have been shown to have many significant impacts on the properties of systems governed by strong interaction, and they may have observable consequences in heavy-ion collision experiments as well as magnetized neutron stars [1,2]. One of the interesting impacts is the behavior of transition temperature of strong-interaction system and chiral condensate as a function of the magnetic field strength B. Early lattice studies of N f = 2 QCD with standard staggered fermions and larger-than-physical pions on N τ = 4 lattices found the so-called magnetic catalyses, which means that the chiral condensate increases monotonically with increasing B [3]. As the chiral symmetry breaking is enhanced it 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 improved staggered fermions and physical pions [4] it turns out that the transition temperature has the opposite behavior in magnetic field as compared to earlier studies in Ref. [3]. I.e. the transition temperature decreases with increasing B. The chiral condensate, on the other hand, first increases and then decreases with B in the transition regime [4]. This nonmonotonic behavior of chiral condensate in B, which is called inverse magnetic catalysis, has also been observed in further lattice studies [5][6][7].
The discrepancy of lattice results in the behavior of chiral condensates and transition temperature in B reported in Refs [3] and [4][5][6][7] is probably due to either large quark masses or possible large cutoff effects present in the first study. To understand the role of quark masses in the intricate relation between the (inverse) magnetic catalysis and the reduction of 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 largerthan-physical values of pions (370 MeV m π 700 MeV) on N τ =6 lattices. It is found that the reduction of transition temperature always holds, however, the inverse magnetic catalysis disappears at a certain value of pion mass. It is suggested in Ref. [8] that the inverse magnetic catalysis is more like a deconfinement catalysis [10] as it is not necessarily associated with the reduction of the transition temperature as a function of B.
Despite the discrepancy mentioned above the strength of the QCD transition always becomes larger in the stronger magnetic field as presented in Refs. [3,11]. And it is speculated that the strength could be sufficiently large to turn the crossover transition into a first order phase transition which can then generate a critical point in the QCD phase diagram in the plane of temperature and magnetic field [11,12]. However, no true phase transition of QCD in 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 structure of N f = 3 QCD in a background magnetic field. At the vanishing magnetic field the true first order phase transition is not yet observed, and state-of-the-art estimates on the critical pion mass m c π based on lattice QCD simulations are m c π 50 MeV using improved staggered fermions [13,14] and m c π 110 MeV using clover-improved Wilson fermions [15,16]. Since the background magnetic field always enhances the strength of the transition one may wonder whether 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 magnetic fields with a smaller-than-physical quark mass. In our lattice simulations we use standard staggered fermions for a testbed towards 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 with a small quark mass also renders us to understand whether the discrepancy in the behavior of chiral condensate and transition temperature as a function of the magnetic field strength in [3,11] is ascribed to the lattice cutoff effects.
The paper is organized as follows. In Section II, we introduce our lattice setup and observables to be investigated. In Section III, we mainly discuss the order of the transition based on results of chiral condensates, Polyakov loops as well as their susceptibilities and Binder cumulants. In Section IV, we conclude our work. The preliminary results have been reported in [18].

II. LATTICE SETUP AND OBSERVABLES
We perform our simulations on N 3 σ ×N τ lattices with 3 mass-degenerate flavors of standard staggered quarks and the Wilson plaquette gauge action by employing the rational hybrid Monte Carlo algorithm [19]. Simulations are performed on lattices with a fixed value of temporal extent N τ = 4 and four different values of spatial size N σ = 8, 16, 20 and 24. Since the critical quark mass, where the first order chiral phase transition starts for this lattice setup is am c =0.027 [20], we choose the value of quark mass in lattice spacing am to be 0.03 in 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 will be reviewed in the following subsection II A. The relevant observables will be introduced in subsection II B.

A. Magnetic fields on the lattice
Magnetic fields couple to quarks with their electric charges and through covariant derivative in the continuum, where g is the SU(3) gauge coupling and q is the electric charge of a quark. A µ and a µ are gauge fields for SU(3) and U(1), respectively. On the lattice the background U (1) magnetic field is introduced by substituting the SU(3) link U µ by its product with the U(1) link u µ in the lattice Dirac operator, In our simulations the magnetic field only points along the z-direction. Since the x-y plane has 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-trivial conditions 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 µ − 1 (µ = x, y, z, t). The background magnetic field pointing along the z-direction B = (0, 0, B) is described by the link variable u µ (n) of the U(1) field, and u µ (n) is expressed as follows in the Landau gauge [4,21], u z (n x , n y , n z , n t ) = u t (n x , n y , n z , n t ) = 1.
The periodic boundary condition for U(1) links is applied for all directions except for the xdirection. One-valuedness of the one-particle wave function along with a plaquette requires the following quantization, where N b ∈ Z is the number of magnetic flux through the unit area for the x-y plane. The ultraviolet 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. with e the electric charge of electron. Thus to satisfy the quantization condition Eq. 4 we take the electric charge q to be that of down quark type, i.e. q = q d,s = − e 3 . To simplify the notation we useb expressed as follows to denote the magnetic field strengtĥ We choose certain values of N b such thatb are the same in physical units among various lattice sizes which are listed in Table I. Our simulation parameters and statistics are summarized in Tab. III -VI. It is worth mentioning that for our largest lattices, i.e. with N σ = 24, we have performed simulations with a hot start (configuration with random elements) and a cold start (configuration with unit elements) to check metastable states around the transition temperature to overcome small tunnelling rates.

B. Observables
An expectation value for an operator O for given gauge coupling β, mass and magnetic where D is the standard staggered Dirac operator and its subscript indicates the electric charge of the related quark. The fractional power, e.g. 1/4 and 2/4, is due to the fourth root of staggered fermions in our simulations.
We measure the chiral condensate for a down type quark 1 , where Tr [· · · ] is a trace over color and space-time, and a factor of 4 in the denominator adjusts for taste degrees of freedom. Its disconnected susceptibility is given by While chiral condensate and its susceptibility are related to chiral symmetry we also compute the Polyakov loop which is related to the deconfinement transition in a pure glue and its susceptibility χ P .
We will also analyse the Binder cumulants B M of order parameters M , e.g. the chiral condensate [22] and the Polyakov loop as well as the constructed order parameter from a mixture of the chiral condensate and the gauge action (cf. Eq. 13). B M is defined as follows where δM = M − M gives the deviation of M from its mean value on a given gauge field configuration. From different distributions of M in different phases the value of B M can be obtained and used to distinguish phase transitions. Given that the chiral condensate is the order parameter of the phase transition 2 , for a first order phase transition, B ψψ = 1 [23]; for a crossover B ψψ 3 [23]; for a second order transition belonging to a Z(2) universality class, B ψψ 1.6 [24].
Since our lattice is rather coarse and simulated systems are in the proximity of transitions, we also estimate the effective number of independent configurations N eff given by where τ int is the integrated autocorrelation time 3 for the chiral condensate. Obtained results of τ int are listed in Appendix A. Hereafter we will show our results of various quantities in a 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 σ = 16 and the history of the chiral condensate at vanishing and nonzero values ofb. In the second subsection, we study the volume dependences of chiral observables and the Polyakov loop as well as their susceptibilities atb = 0 and 1.5 obtained from lattices with N σ =8, 16,20 and 24. In the third subsection, we show results based on a appropriate order parameter constructed 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 forb = 0, 1.5 and 2 obtained from lattices with N σ =16 as a function of the inverse gauge coupling β in the left and right plot of Fig. 1, respectively. We recall that the temperature is an increasing function of β. As seen from the left plot the value of the chiral condensate is enhanced by the magnetic field in the whole temperature region. Namely only the magnetic catalysis is observed. One 2 The same holds true for other order parameters. 3 The autocorrelation time and its error are given in Appendix A where τ int are rounded in integer for simplicity.
can also see that the critical β value where the chiral condensate drops most rapidly becomes larger as the magnetic field strength increases. This indicates that the transition temperature increases as the magnetic field becomes stronger. Meanwhile as the strength of the magnetic field increases the dropping of chiral condensates becomes more rapidly. This means that the transition becomes stronger with stronger magnetic fields. These two observations are also visible in the behaviour of disconnected chiral susceptibilities as shown in the right plot of Fig. 1. I.e. the peak location of disconnected chiral susceptibility shifts to a larger value of β and the peak height of disconnected chiral susceptibility increases as the strength of magnetic field increases.  We show similar plots as Fig. 1 for the Polyakov loop and its susceptibility atb = 0, 1.5 and 2 in Fig. 2. The peak location of the Polyakov loop susceptibility as well as the inflection point of the Polyakov loop shifts to larger values of β, i.e. higher transition temperature with stronger magnetic field. This is similar to the observation from chiral observables. Besides that transition temperatures signalled by the Polyakov loop and its susceptibility are close to those obtained from chiral condensates and disconnected chiral susceptibility. What's more, it is interesting to see that the peak height of the Polyakov loop susceptibility also becomes higher in a stronger magnetic field, although the Polyakov loop is an order parameter for the confinement-deconfinement phase transition of a pure glue system while the chiral condensate is the order parameter of QCD transitions with vanishing quark massless. This could be due to the fact that neither of these two quantities are the true order parameter but a part of the true order parameter in N f = 3 QCD 4 [14,25].
At vanishing magnetic field the transition is a crossover, and as the magnetic field strength becomes larger the QCD transition becomes stronger. In particular the jumping behavior of the chiral condensate and the Polyakov loop atb=2 is seen from left plots of Fig. 1 and vanishing magnetic field (b=0) there does not exist any metastable behavior, while in the case ofb =1.5 and 2 the metastable behavior becomes obvious. To confirm the metastable behavior seen from the volume of N σ =16, we also study the case ofb =1.5 with a larger volume, i.e. on 24 3 × 4 lattices. Since in the first order phase transition the tunnelling rate between two metastable states becomes smaller in the larger volume, here we rather investigate on the time history of the chiral condensate obtained from two different kinds of streams, i.e. one starting from a unit configuration (cold) and the other one starting from a random configuration (hot). If there is no first order phase transition chiral condensates obtained from the cold and hot starts will always overlap after thermalization, and if there exist a first order phase transition the two steams from cold and hot starts will stay apart and tunnel from one to the other as the two metastable states in the first order phase 4 Note that the peak height of Polyakov loop susceptibility also increases as the system approaches the first order phase transition region with smaller values of quark masses for the case of zero magnetic field strengthb = 0 [25].   transition. The former is observed for the case of β = 5.1689 while the latter is clearly seen for β = 5.1698.
In a short summary, above observations suggest that the transition tends to be a first order for magnetic fields ofb ≥1.5.
B. Volume dependences of observables withb = 0 and 1.5 In this subsection, we discuss more on the volume dependence of observables atb = 0 and 1.5 to further confirm the onset of the first order phase transition atb ≥ 1.5. In   20 and 24 are consistent among each other, while large deviations are seen for the results obtained from N σ =8. In the case ofb=1.5 the finite size effect seems to appear at β smaller and close to β c starting with a larger volume, i.e. N σ =16. This could be due to the fact that the system with the presence of a magnetic field tends to have a stronger transition and consequently more statistics is needed to get robust results 5 , and that the correlation length in 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 volumes for both vanishing and nonzero magnetic fields. The is also reflected in the peak locations of disconnected susceptibilities shown in the right column of Fig. 4. In the case of a first order phase transition the disconnected chiral susceptibility should grow linearly in volume. It is more or less the case for disconnected susceptibilities atb=1.5 as seen from bottom right plot of Fig. 4 and Table II. And as expected that in the case of vanishing magnetic field the disconnected susceptibility grows slower than linearly in volume as there exists a crossover transition.
FIG. 6. Same as Fig. 4 but for the Polyakov loop and its susceptibility.

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 not a true order parameter, but is part of a mixture of operators that defines the true order parameter M [14,25]. In three flavor QCD, such an order parameter can be constructed as a combination of the plaquette action and the chiral condensate as follows [14,25] M (β, s) = ψψ(β) + s 1 where S gauge = 6N 3 σ N τP andP is the plaquette. Correspondingly, its susceptibility is, In this work, we do not intend to determine the mixing parameter s and just use the value obtained in the previous work for a √ eB = 0 in [20], i.e. s = −0.8. As will be seen next the results obtained using s = −0.8 does not change from s = 0 qualitatively (cf. Table II).
We show in the left panel of Fig. 8    In summary, we conclude that the system with am = 0.03 with magnetic fields a √ eB ≥1.5 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 quark and the Wilson plaquette gauge action on N τ = 4 lattices. We started with simulations at a fixed quark mass of am = 0.03 at vanishing magnetic field where a crossover transition is observed. After turning on the magnetic field we studied dependences of chiral observables, We do not find any signs of inverse magnetic catalysis and the reduction of transition temperature as a function of the magnetic field strength, and this is consistent with the results obtained from lattice studies of N f = 2 QCD with the standard staggered quarks and larger-than-physical pions [3,26]. Since our findings of magnetic catalysis and increasing of transition temperature with the magnetic field strength even holds in the regime occurring a first order phase transition this suggests that the discrepancy from results using the improved staggered 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 fermions the strength of transition also increases with increasing magnetic field strength and a first order phase transition has not yet been observed in such studies. Although we find a first order phase transition in the current study using the standard staggered fermions, we do not intend to provide a precise determination of a critical magnetic field in which the transition turns into a first/second order phase transition in the current discretization scheme. We will leave it for future studies using improved staggered fermions, such as HISQ fermions to achieve more realistic results on the chiral phase structure of N f = 3 QCD [17]. As in current studies severe critical slowing down is expected for simulations in the vicinity of phase transitions with larger volumes (see Appendix A), the autocorrelation length will probably become even longer in simulations with smaller lattice spacings or improved actions, in which the multicanonical method [27,28] or other extended Monte Carlo method reviewed in [29] might help. In Fig. 9, we estimate the dynamic critical exponent of the system by using the longest autocorrelation length in each volume. A fit ansatz of log τ int (N σ ) = z log N σ + c with fit parameters z and c is adopted and this fit ansatz is the same form as τ int (N σ ) ∼ N z σ , which is the definition of the dynamic critical exponent. One can see that the system is affected by severe critical slowing down atb ≡ √ a 2 eB = 1.5. It is thus practically challenging to investigate the first-order phase transition by using conventional direct simulations with HMC in the thermodynamic limit.

Appendix C: Summary of statistics
We summarize our simulation parameters and statistics in Tab. III -VI. If the statistics is sufficient the binning size for the Jackknife analysis is taken as the Bin-size 2τ int , otherwise the Bin-size is adjusted such that it is the divisor of number of trajectories by around 10.
In the latter case, the error might be underestimated.