QCD phase diagram in a magnetic background for different values of the pion mass

We investigate the behavior of the pseudo-critical temperature of $N_f = 2+1$ QCD as a function of a static magnetic background field for different values of the pion mass, going up to $m_\pi \simeq 660$ MeV. The study is performed by lattice QCD simulations, adopting a stout staggered discretization of the theory on lattices with $N_t = 6$ slices in the Euclidean temporal direction; for each value of the pion mass the temperature is changed moving along a line of constant physics. We find that the decrease of $T_c$ as a function of $B$, which is observed for physical quark masses, persists in the whole explored mass range, even if the relative variation of $T_c$ appears to be a decreasing function of $m_\pi$, approaching zero in the quenched limit. The location of $T_c$ is based on the renormalized quark condensate and its susceptibility; determinations based on the Polyakov loop lead to compatible results. On the contrary, inverse magnetic catalysis, i.e. the decrease of the quark condensate as a function of $B$ in some temperature range around $T_c$, is not observed when the pion mass is high enough. That supports the idea that inverse magnetic catalysis might be a secondary phenomenon, while the modifications induced by the magnetic background on the gauge field distribution and on the confining properties of the medium could play a primary role in the whole range of pion masses

We investigate the behavior of the pseudo-critical temperature of N f = 2+1 QCD as a function of a static magnetic background field for different values of the pion mass, going up to mπ ≃ 660 MeV. The study is performed by lattice QCD simulations, adopting a stout staggered discretization of the theory on lattices with Nt = 6 slices in the Euclidean temporal direction; for each value of the pion mass the temperature is changed moving along a line of constant physics. We find that the decrease of Tc as a function of B, which is observed for physical quark masses, persists in the whole explored mass range, even if the relative variation of Tc appears to be a decreasing function of mπ, approaching zero in the quenched limit. The location of Tc is based on the renormalized quark condensate and its susceptibility; determinations based on the Polyakov loop lead to compatible results. On the contrary, inverse magnetic catalysis, i.e. the decrease of the quark condensate as a function of B in some temperature range around Tc, is not observed when the pion mass is high enough. That supports the idea that inverse magnetic catalysis might be a secondary phenomenon, while the modifications induced by the magnetic background on the gauge field distribution and on the confining properties of the medium could play a primary role in the whole range of pion masses.

I. INTRODUCTION
The investigation of the properties of strong interactions in a magnetic background field has been the subject of numerous studies in the recent past, see, e.g., Refs. [1,2] for recent reviews. Part of the interest is connected with phenomenology, since strong background fields are expected in non-central heavy ion collisions [3][4][5][6][7][8], in some astrophysical objects like magnetars [9], and might have been produced during the cosmological electroweak phase transition [10,11]. However, the issue is interesting also from a purely theoretical point of view, since background fields such that eB Λ 2 QCD represent non-trivial probes of the non-perturbative properties of strong interactions. Lattice QCD simulations have been an essential tool to advance knowledge in this field, given the fact that no technical issues, such as a sign problem, appear when one introduces a magnetic background coupled to dynamical quark fields in the path-integral formulation of QCD.
An important aspect regards the influence of the magnetic field on the QCD phase diagram. Early lattice studies of N f = 2 QCD, adopting unimproved staggered fermions and larger-than-physical quark masses, showed a slightly increasing behavior of the pseudo-critical temperature as a function of the magnetic field [12]. Those preliminary indications however changed when new numerical simulations, exploiting an improved discretiza- * Electronic address: massimo.delia@unipi.it † Electronic address: floriano.manigrasso@pi.infn.it ‡ Electronic address: fnegro@pi.infn.it § Electronic address: sanfilippo@roma3.infn.it tion of N f = 2 + 1 QCD with physical quark masses, showed instead a substantial decrease of T c , of the order of 10-20% for |e|B ∼ 1 GeV 2 [13]. The reason for the discrepancy was ascribed to either the large quark masses, or possibly the large cut-off effects present in the first study.
Actually, while magnetic catalysis in the QCD vacuum (i.e. at T = 0) was confirmed by several lattice simulations [33][34][35][36][37][38], a new unexpected behavior was discovered for temperatures around T c , consisting in a decrease, rather than increase, of the quark condensate as a function of the magnetic field intensity [13,37]. This phenomenon, later confirmed by further lattice studies [38][39][40][41], was named as inverse magnetic catalysis, a name soon extended to indicate the decreasing behavior of T c (B) itself, thus assuming implicitly that the latter is caused by the former. Many efforts have been done since then to interpret the new phenomenology within various model approaches .
For this reason, one may ask whether the decrease of T c could not be ascribed to these effects and thus be a sort of deconfinement catalysis, with the decrease of the quark condensate being just a secondary effect, induced by the fact that the magnetic background fosters deconfinement, hence chiral symmetry restoration. Actually, in view of the strict entaglement between confinement and chiral symmetry breaking, which is not yet fully understood, the question of understanding which is the true driving phenomenon could be ill-posed.
However, from a practical point of view, we can ask how the situation changes as the quark mass spectrum is changed. In particular, adopting substantially larger than physical quark masses, i.e. approaching the quenched limit, one can explore a regime where chiral symmetry ceases to be a good symmetry of the theory, while symmetries associated to confinement, like center symmetry, become more and more relevant. Then, one can try to asnwer some some clear-cut questions: i) is the decrease of T c as a function of B still observed?
ii) in case it is, is the decrease always associated with inverse magnetic catalysis?
In this study we try to answer these questions adopting the same discretization of N f = 2 + 1 used in Ref. [13] and exploring three different lines of constant physics, corresponding to pseudo-Goldstone pion masses m π = 343, 440 and 664 MeV on lattices with N t = 6. Anticipating our main results, the answer to the first question is yes, while that to the second question is no. This is of course also important to understand the reason of the discrepancy observed in early lattice studies where a slight increase of T c , instead of a decrease, was observed. The unphysical quark mass spectrum adopted in Ref. [12] was considered as one possible reason [13]. Our present results suggest that the different behavior can be likely ascribed to large lattice artifacts induced by the coarse lattice spacing and unimproved discretization adopted in Ref. [12]. That is in agreement with recent studies [119], adopting the same unimproved discretization of Ref. [12], where the increase of T c as a function of B is observed even for smaller than physical quark masses.
The paper is organized as follows. In Sec. II we describe our lattice discretization of N f = 2 + 1 QCD with a magnetic background, as well as the physical observables and the numerical setup to work on lines of constant physics. In Sec. III we present results regarding the behavior of the pseudo-critical temperature as a function of the pion mass at B = 0, while Sec. IV contains our main results obtained in the presence of a magnetic background. Finally, in Sec. V we investigate the fate of inverse magnetic catalysis in the large quark mass limit and in Sec. VI we draw our conclusions.

II. NUMERICAL METHODS
In this study we investigated N f = 2 + 1 QCD at different values of the quarks masses, adopting a setup with two degenerate up and down flavors, m u = m d = m ℓ , and a strange-to-light mass ratio fixed at its physical value m s /m ℓ = 28.15. We adopted the stoutsmearing improved staggered quark action and the treelevel Symanzik improved gauge action [120,121]. With this choice the discretized partition function takes the form where DU is the functional integration over all the possible SU (3) gauge field configurations. The gauge action S YM is written in terms of the real part of the trace of the 1 × 1 and 1 × 2 Wilson loops (respectively denoted as P 1×1 i;µν and P 1×2 i;µν ): The Dirac operator M f st [m f , B] related to the flavor f which appears in Z(m ℓ , B) takes the form where the staggered phases are denoted by η i; ν and where U (2) i; µ is the two times stout-smeared [122] SU (3) link variable (with an isotropic smearing parameter ρ = 0.15) at position i pointing along the direction µ. The configuration of U (1) link variables u f i; µ is fixed at the beginning of each simulation in order to reproduce the desired value of the external uniform magnetic field B z . The only nontrivial phases are: where the charge of each flavor f is denoted as q f (q u = −2q d = −2q s = 2e/3). This expression has been derived to describe the four-potential of a uniform magnetic field over a manifold with periodic boundary conditions, such as the lattice discretized space we adopt. Such periodicity constrains B z to quantized values [123][124][125] where b is integer valued. The partition function is periodic in b with period N x N y . Numerical simulations have been performed using the Rational Hybrid Monte-Carlo algorithm (RHMC) [126] implemented in the NISSA code [127] and in the Open-StaPLE code for GPUs [128,129]. We have performed around 100 runs with different combinations of T and B for each value of the pion mass, with average statistics of approximately 3000 RHMC trajectories for each run.

A. Lines of constant physics
To study the chiral symmetry restoration crossover at different values of the pion mass m π , we need to perform, for each chosen value of m π , several simulations at different values of the temperature T = 1/(N t a), i.e. at different lattice spacings a, since we work at fixed N t = 6. This requires the preliminary knowledge of the lines in the β − m ℓ plane along which the pion mass stays constant at the chosen values (m π = 343, 440 and 664 MeV), which we refer to as lines of constant physics (LCP).
To achieve this, we performed preliminary T = 0 numerical simulations on a 24 3 × 32 lattice at 7 values of β (β = 3.45, 3.55, 3.62, 3.73, 3.85, 4.00, 4.10) and at 5 -7 values of m ℓ for each β (in the range 0.007 m ℓ 0.19), for a total of 42 simulation points. Then, we associated each simulation point with a value of the lattice spacing a by adopting the w 0 scale setting approach [130], based on the gradient flow technique [131], and assuming that w 0 is independent of m π . Moreover, we extracted the lightest pseudoscalar pion mass m π from the decay of Euclidean time correlators of the appropriate staggered quark operators. The obtained values of a and m π extend from ∼ 0.07 fm to ∼ 0.3 fm and from ∼ 300 MeV to ∼ 700 MeV respectively, covering completely the parameter region we are interested in.
At fixed β, we interpolated the pion mass m π as function of m ℓ with a 3rd order polynomial. This allows us to determine, for each β, the bare quark mass m ℓ that corresponds to the chosen value of the pion mass. The LCP are obtained by interpolating m ℓ as a function of β with a 4th order polynomial, as shown in Fig. 1.
For the determination of the lattice spacing along these lines we adopted a similar procedure. Again, at fixed β, we interpolated the dependence of the lattice spacing a on the pion mass m π with the function which approaches a finite value as we go towards the quenched limit m π → ∞. In this way we extract the value of a at the chosen value of the pion mass for all the explored β values. Then we give an estimate of the lattice spacing along the LCP by interpolating these data   by a 5th order polynomial 1 in β. The results of these interpolations are shown in Fig. 2.

B. Physical observables
For the purpose of this study and, in particular, for the discussion of the fate of (inverse) magnetic catalysis, we focused on the quark condensate and its susceptibility.
The quark condensate of the flavor f is expressed as where T is the temperature, V is the spatial volume and the trace is evaluated, as usual, by means of noisy estimators. At non-zero magnetic field, the two light quarks condensates differ because of their electric charge: they couple in a different way to the magnetic field. We introduce the light quark condensate Σ ℓ defined as This observable is affected by both additive and multiplicative renormalizations. As it has been pointed out in Ref. [13], the presence of the external magnetic field does not introduce new B−dependent divergencies. For this reason, the renormalization prescription introduced in Ref. [132], which exploits T = 0 quantities to perform the additive renormalization and the value of the bare quark mass to take care of multiplicative ones, can be extended to the B = 0 case as where the two condensates are computed at the same UV cutoff (i.e. the same bare parameters). The location of T c is usually defined, in terms of the renormalized light condensate, as the point of maximum slope, i.e. the point where Σ r ℓ has an inflection point as a function of T and the absolute value of ∂Σ r ℓ /∂T reaches a maximum. Alternatively, one can consider the behavior of the chiral susceptibility, i.e.
For χ ℓ the renormalization is perfomed in a similar way, and we look for the peak of the dimensionless ratio χ r ℓ (T, B)/m 4 π to locate T c . For the renormalization of both Σ ℓ and χ ℓ we have exploited the same zero temperature runs that we used to determine the LCP.
In Ref. [35] it was observed that the change in the quark condensate due to the magnetic field can be ascribed both to sea (dynamical in Ref. [35]) and to valence quarks effects. At T = 0, and for small enough magnetic fields, the two contributions add approximately to the total change of the condensate, with a sea contribution amounting to ∼ 30% of the total signal. The sea contribution is particularly interesting because it is only related to the modification of the gauge configurations which are mostly relevant in the path integral. Moreover, it was observed [39] that around T c , the sea contribution changes sign, possibly resulting in the observed overall decrease of the condensate known as inverse magnetic catalysis. It is then useful to introduce the sea quark condensate, defined as the standard condensate but with the external field switched off in the observable: The sea contribution to the change of the light condensate at finite B is then expressed as 0) . (12) Finally, as a further observable to locate the crossover, we introduce the unrenormalized Polyakov loop which is more directly connected to the confining properties of strong interactions and becomes an exact order parameter in the infinite mass limit that we are somehow approaching. As for the quark condensate, we will look for its inflection point to locate T c .

C. Determining observables as a function of T at fixed values of the magnetic field
When moving along a LCP, one would like to keep the magnetic field strength constant as well. Since eB is given by Eq. (4), one should tune the parameter b to balance the change of the lattice spacing along the LCP; however, being b integer-valued, this is not possible in practice. We have approached this problem similarly to Ref. [13]: simulations have been performed at various temperatures and, for each T , at various values of b. Then, to reach the desired value of eB, observables have been interpolated by a cubic spline in b, arranged in such a way that its first derivative is zero at b = 0, as expected by analyticity in eB and charge conjugation symmetry.
We report in Fig. 3 the explored grid of such simulation points for m π = 664 MeV, together with the lines of constant eB (in particular, eB = 0.1, 0.4, 0.67 GeV 2 ). Such lines are characterized by the equation A typical interpolation of the quark condensate Σ ℓ is shown in Fig. 4, where data points corresponds to m π = 343 MeV at T = 141 MeV.
III. PSEUDO-CRITICAL TEMPERATURE AT B = 0 We illustrate our results starting from the determination of the pseudo-critical temperature at B = 0, which  represents the starting input to study the dependence of T c on B. In Fig. 5, the renormalized quark condensate and its susceptibility are shown as a function of T for the different explored values of m π . The peak of the chiral susceptibility, as well as the inflection point of the quark condensate, clearly shift to higher temperatures as m π increases; at the same time, a general weakening of the crossover strength is observed, which is consistent with the fact that chiral symmetry restoration is less and less relevant to strong interaction dynamics as the quenched limit is approached.
In order to determine T c , the inflection point of the renormalized condensate has been located by fitting data according to an arctangent or a cubic polynomial ansatz, while the maximum of the susceptibility peak has been obtained by fitting data according to a quadratic or a  Lorentzian function of T . In both cases, systematic errors have been estimated by varying the range of fitted data points or the fitting function. Results obtained from the inflection point of the condensate are reported in Table I. It is interesting to consider how T c depends on m π . To that aim, in Fig. 6 we plot the results reported in Table I together with a determination at the physical point adopting the same discretization on N t = 6 lattices, T c (m π = 135 MeV) = 149(3), that we have taken from Ref. [133]. On general grounds, one expects that T c approaches finite values both in the quenched limit, lim mπ→∞ T c ≡ T quench c , and in the chiral limit, lim mπ→0 T c ≡ T χ c . Moreover, approaching the chiral limit, one expects where β and δ are the critical indexes describing the critical behavior around the chiral point, even if that could be not relevant to our range of pion masses, which is far away from the chiral limit. A previous study of the dependence of T c on m π has been reported in Ref. [134], adopting a different staggered discretization and a range of large pion masses similar to ours: in that case it was found that T c (m π ) behaves more or less linearly in m π when approaching the chiral limit, a result which is not far from the prediction of Eq. (15) for various universality classes which could be relevant in the case of a second order chiral transition (e.g., O (2) or O(4)), 2/(βδ) ∼ 1.2. Inspired by these considerations and previous findings, we have tried to fit the data reported in Fig. 6 according to (16) and, fixing T quench c = 270 MeV, we obtain T χ c = 128(4) MeV and M = 763(39) MeV withχ 2 = 1.54/2. The value obtained for the critical temperature in the chiral limit, T χ c ≡ T c (m π = 0), is not unreasonable, given the preliminary estimate T χ c = 138(5) MeV in the continuum limit recently reported in Ref. [135], and given that our estimate represents a very rough extrapolation from a region of large pion masses.
The functional dependence in Eq. (16) is by no means unique: various other functions, sharing the same rough properties exposed above, fit data equally well, like for instance 2 Taking m B π with B = 1 as an argument, which could accomplish for the possible critical behavior around the chiral point, works well also for the arctan and the exponential ansatz.

IV. PSEUDO-CRITICAL TEMPERATURE AT NON-ZERO MAGNETIC FIELD
The main result of our study can be already appreciated by looking at Figs. 7 and 8, where we plot respectively the renormalized condensate and the renormalized chiral susceptibility as a function of the temperature for different values of the magnetic background and of the pion masses. The inflection point of the condensate and the peak of the susceptibility always move to lower temperatures as B is increased. The inflection point of the unrenormalized Polyakov loop, which is plotted in Fig. 9, shows a similar behavior. It is interesting to notice, especially looking at the behavior of the chiral susceptibility peak, that at the same time the strength of the crossover transition seems to increase.
For each value of eB, we have determined the location of T c using the same procedure (also regarding the determination of the systematic error) already described above, i.e. the inflection point for the renormalized con- densate and for the Polyakov loop, the maximum of the peak for the chiral susceptibility. Results obtained for m π = 440 and for all different observables are reported in Fig. 10: the temperature decreases in a similar way in all cases. In Fig. 11 we report instead the pseudocritical temperature obtained from the inflection point of the renormalized condensate as a function of eB for the different values of the pion mass. It can be deduced that the decrease of T c with B is a general phenomenon which, from a qualitative point of view, is independent of the quark mass spectrum. However, some dependence on m π must be present at a quantitative level, because the influence of the magnetic field on the gluon field distributions takes place only through dynamical quarks, so when their masses go to infinity the magnetic field decouples and T c must become independent of B. A tendency for a less pronounced influence as m π increases is qualitatively visible in Fig. 12, where we report T c (eB)/T c (0) obtained from the renormalized quark condensate for the different pion masses. In order to check this expectation on a quantitative basis, we have tried to describe the small-B dependence of T c in terms of a curvature coefficient obtaining the results reported in Table II and in Fig. 13. The reported errors on v B take into account the systematic uncertainty related to the choice of the range over which a best fit according to Eq. (17) is performed; moreover, the value at the physical pion mass has been obtained by fitting data for T c (eB) reported for N t = 6 in Ref. [13]. The dotted line in Fig. 13 represents the result of a best fit to a behavior v B = a/m α π which yields α = 0.62 (12) and serves just to prove that data are compatible with a vanishing v B in the quenched limit, how- ever we stress that different behaviors, obtained by fixing α = 1 or choosing v B = a exp(−m π /b) work equally well and we are not able to fix the actual dependence of v B , due to the present precision of our data. Moreover, since we are considering data at fixed values of the cut-off (N t = 6), the actual observed power law might be affected by the fact that, with our lattice discretization, just one pion becomes massless as the chiral limit is approached. As a final remark, we stress that the results presented in this section show that the discrepancy between the results reported in Refs. [12] and [13] was just due to discretization effects, and not the result of a different quark mass spectrum. This is also supported by Ref. [119] where, adopting the same unimproved discretization of Ref. [12], T c continues to be an increasing function of B, even for lighter-than-physical quark masses.

V. INVERSE OR DIRECT MAGNETIC CATALYSIS AROUND Tc
In QCD with physical quark masses, the decreasing behavior of T c as a function of B is associated with another interesting and unexpected phenomenon taking place around T c : the enhancement of chiral symmetry breaking, which is expected on general grounds and actually observed at T = 0, is reversed around T c , where instead ψ ψ becomes a decreasing function of B. This phenomenon has been given the name inverse magnetic catalysis, which has a clear correspondence with what happens. However, with a less obvious correspondence, the same name has been usually assigned also to the decreasing behavior of T c itself.
Apart from merely lexical issues, the question is whether the observed decrease of the condensate is actually the driving phenomenon leading to the decrease of T c , or if, on the contrary, some other phenomena force T c to decrease, like those related to the influence of B on the confining properties. In this scenario, the observed behavior of ψ ψ would be a secondary phenomenon: around the transition, when B is switched on while keeping the temperature T fixed, the condensate decreases just because that temperature is moving into the deconfined and chirally restored phase. However let us say that, given our present ignorance about the relation existing between confinement and chiral symmetry breaking, the question itself might be ill-posed.
Nevertheless, from a practical point of view, we can investigate if the relation between inverse magnetic catalysis and the decrease of T c is maintained also for largerthan-physical values of the pion mass. To that purpose, in Fig. 14 we report the variation of the quark condensate, measured with respect to its value at B = 0, as a function of the temperature and for two different values of the magnetic field. Inverse magnetic catalysis is well visible, in a region around T c , for m π = 343 MeV; however the effect becomes barely visible for m π = 440 MeV and completely disappears for m π = 664 MeV, where the condensate is always an increasing function of B at all temperatures. That gives evidence that one may have a decreasing behavior of T c even in absence of inverse magnetic catalysis, challenging the strict connection, or even identfication, which has been usually assumed between the two phenomena.
On the other hand, when one considers the separation of the modification of the quark condensate into a valence and sea contribution, one still continues to see a behavior consistent with inverse magnetic catalysis in the sea contribution only and for all values of the explored pion masses. This is visible in Fig. 15, where we report the sea contribution ∆Σ sea ℓ (T, B) defined in Eq. (12) for the largest pion mass and for a set of temperatures around the transition.
The fact that the sea contribution still continues to be a decreasing function of B around T c is of course related to how the gauge field distribution gets modified by B, so as to move towards (or deeper into) the deconfined/chirally restored phase. As a consequence, the distribution of eigenvalues of the B = 0 Dirac operator changes, and small eigenvalues are suppressed.

VI. DISCUSSION AND CONCLUSIONS
We have investigated the modification of the pseudocritical temperature of QCD with N f = 2 + 1 flavors as a magnetic background field is switched on. The study has been done on lattices with a fixed temporal extent, N t = 6, changing the temperature along lines of constant physics corresponding to three different values of the pseudo-Goldstone pion mass, m π = 343, 440 and 664 MeV.
We have found that the pseudo-critical temperature has always a decreasing behavior as a function of eB, even if the relative variation, measured in terms of the curvature coefficient v B defined in Eq. (17), appears to be a decreasing function of m π approaching zero in the quenched limit, as expected. At the same time, we have observed that inverse magnetic catalysis disappears at all temperatures for the largest value of the pion mass, even if it is maintained in the sea contribution only. Finally, we have observed that the magnetic field induces a strengthening of the crossover, which is agreement with previous lattice studies [12,13,119] as well as with the prediction for a first order transition at large enough magnetic field strengths [41,136]. Present results are limited to a single value of the temporal extension, N t = 6, therefore the exploration of finer lattice spacings would be welcome in the future.
The persistence of the decreasing behavior of T c (B) observed even for large values of the pion mass, where the chiral properties of the theory are not relevant, and possibly up to the quenched limit, clarifies definitely the origin of the discrepancy between Refs. [12] and [13] and sheds some light on the origin of this phenomenon. The fact that it is qualitatively independent of the quark mass spectrum is in agreement with model computations which have found evidence of it in a large-N c framework (see, e.g., Refs. [42,77]). The fact that it is not necessarily associated with a decrease of the quark condensate as a function of B would suggest to name the phenomenon as deconfinement catalysis [116] rather than inverse magnetic catalysis.
The paramagnetic behavior of the deconfined phase of strongly interacting matter has been sometimes suggested as a possible origin of this deconfinement catalysis. Actually, such paramagnetic behavior has been observed even when adopting the same unimproved lattice action leading to an increase of T c [12], however the results of Refs. [101,103] show that coarse and unimproved discretizations tend to make the paramagnetic behavior weaker, especially around the pseudo-critical temperature.
The direct effects on the confining properties of the theory, which are observed both a zero and at finite temperature, are other natural candidates to explain the deconfinement catalysis. At finite temperatures below T c the magnetic field suppresses the string tension, while strong enough magnetic fields could even lead to an anisotropic deconfinement at T = 0 [116]. It would be interesting, in the future, to investigate if, analogously to what happens for the behavior of T c , such effects on the confining properties persists also for larger-than-physical quark masses and possibly towards the quenched limit. At the same time, in view of the observed strengthening of the transition induced by the magnetic field, it would be interesting to investigate if the first order region which is found around the quenched limit is enlarged by the presence of the magnetic background.