Magnetized baryons and the QCD phase diagram: NJL model meets the lattice

We determine the baryon spectrum of 1 + 1 + 1-flavor QCD in the presence of strong background magnetic fields using lattice simulations at physical quark masses for the first time. Our results show a splitting within multiplets according to the electric charge of the baryons and reveal, in particular, a reduction of the nucleon masses for strong magnetic fields. This first-principles input is used to define constituent quark masses and is employed to set the free parameters of the Polyakov loop-extended Nambu-Jona-Lasinio (PNJL) model in a magnetic field-dependent manner. The so constructed model is shown to exhibit inverse magnetic catalysis at high temperatures and a reduction of the transition temperature as the magnetic field grows -- in line with non-perturbative lattice results. This is contrary to the naive variant of this model, which gives incorrect results for this fundamental phase diagram. Our findings demonstrate that the PNJL model can be reconciled with the lattice findings in a systematic way, employing solely zero-temperature first-principles input.


Introduction
The impact of background electromagnetic fields on strongly interacting matter is relevant for a range of physical situations including off-central heavy-ion collisions, magnetized neutron stars and the evolution of the early universe [1]. In particular, the elementary properties of magnetized hadronic degrees of freedom are important for cold astrophysical environments. The masses of baryons and mesons enter the nuclear equation of state and influence the mass-radius relations of magnetars. Together with hadronic decay rates, these also affect stability of such compact objects and cooling mechanisms that characterize the emitted neutrino spectra [2]. For heavy-ion collisions, the magnetic field is produced in the very early stages and is expected to be short-lived [3,4], primarily affecting heavy baryons. A special role might be played by charged vector mesons that were conjectured to condense for sufficiently strong magnetic fields [5].
Besides their phenomenological importance, magnetic fields also represent external probes of strongly interacting matter i.e. of the underlying theory, quantum chromodynamics (QCD). One particular feature of the magneto-response of QCD matter that received great attention in the last decade is the phase diagram for nonzero temperatures and static, spatially uniform background magnetic fields, see, e.g., the review [6]. This phase diagram features a chiral symmetry restoration/deconfinement crossover [7,8], where the chiral condensateψψ drops towards zero and, almost simultaneously, the Polyakov loop P rises. According to lattice simulations, the pseudo-critical temperature T c , where the transition occurs, is reduced as the magnetic field strength B grows [9][10][11]. For physical quark masses (i.e. such that the pion mass is M π = 135 MeV), this behavior emerges due to the non-trivial dependence ofψψ on the temperature and on the magnetic field. On the one hand, for temperatures well below T c the magnetic field enhancesψψ (a phenomenon referred to as magnetic catalysis [12]). On the other hand, for T ≈ T c the opposite is observed andψψ is reduced by B (inverse magnetic catalysis [13]). While magnetic catalysis originates from the high degeneracy of the lowest Landau-level [14,15], inverse magnetic catalysis arises as a result of the rearrangement of gluonic configurations induced by the magnetic field -it is thus a secondary effect that can be associated to the indirect interaction between B and electrically neutral gluons via sea quark loops [13]. This mechanism is suppressed if quarks are heavy. Indeed, contrary to the situation at the physical point, for sufficiently heavy quarks (M π 500 MeV), inverse magnetic catalysis does not occur anymore [16] -nevertheless, the transition temperature is still reduced by B [17].
The above summarized results are based on first-principles lattice QCD simulations. In addition to the lattice approach, a multitude of low-energy models and effective theories of QCD were also employed to investigate the phase diagram for B > 0. Surprisingly, the initial studies observed the exact opposite of the lattice results: magnetic catalysis at all temperatures and the enhancement of the transition temperature with growing B. A prime example for this behavior was obtained in the Polyakov loop-extended Nambu-Jona-Lasinio (PNJL) model [18], but various other models resulted in the same picture, see the reviews [6,19]. The failure of these approaches was associated to the fact that gluons merely enter as a static background in these models so that the indirect mechanism behind inverse magnetic catalysis cannot be truly captured. Later it was recognized that including a B-dependence in model parameters might improve the situation and bring model calculations closer to the lattice results, see e.g., Refs. [20][21][22][23][24][25][26][27]. In a PNJL model study [21], this was performed by tuning the coupling G(B) to reproduce the transition temperature T c (B) obtained on the lattice. While this shows that the model can be made compatible with full QCD, in this example the predictive power of the effective approach is clearly lost.
Let us emphasize that effective models, albeit approximations to full QCD, are helpful for identifying the relevant degrees of freedom and interaction mechanisms, and thus guide our understanding of the physics of strongly interacting matter in extreme environments. For large baryon chemical potentials, where lattice simulations are hindered by the sign problem, low-energy models represent one of the few possibilities for the investigation of the phase diagram. Therefore it is highly desirable to test the limitations of such models in cases where importance sampling-based lattice investigations can be performed -like the phase diagram at nonzero magnetic field or at nonzero isospin density [28].
In the present paper our aim is to develop a systematic approach to fix the parameters of the PNJL model utilizing first-principles input at zero temperature. In particular we determine the baryon spectrum in three-flavor QCD using continuum extrapolated lattice simulations with physical quark masses. From this analysis, T = 0 constituent quark masses are inferred and used to set the model parameters in a magnetic field-dependent manner 1 . According to our results, the phase diagram of the so constructed lattice-improved PNJL model agrees with all features of the available lattice findings. Our method may also be extended to further low-energy models of QCD.
Besides fixing free parameters of effective descriptions, our results constitute the first lattice determination of magnetized baryon masses at the physical point. This complements earlier lattice calculations of baryon masses with heavier-than-physical quarks [30][31][32], meson masses [9,33,34] and decay rates [35], and properties of heavy quarkonia [36]. Our results might provide useful information for magnetized compact stars and the early stages of heavy-ion collisions, as pointed out above. This paper is structured as follows. In Sec. 2 we describe our numerical setup and measurement strategy and present the results for the baryon spectrum. This is followed by Sec. 3, where the definition of the constituent quark masses and the details of our PNJL model are given. The results for the magnetic field-dependent model parameters and the thermodynamics of the model is presented in Sec. 4. Finally in Sec. 5 we summarize our findings and give an outlook for potential future research.
Our numerical simulations are performed on N 3 s × N t lattices with spacing a, using the tree-level Symanzik improved gauge action and three flavors (u, d and s) of stout-improved rooted staggered quarks [37]. The quark masses m u = m d and m s are set to their physical values along the line of constant physics [38]. The magnetic field is chosen to point in the z direction and is implemented via U(1) phases satisfying periodic boundary conditions [9]. This setup gives rise to a quantized magnetic where e > 0 is the elementary charge and the quark electric charges are set as q u = −2q d = −2q s = 2e/3. The details of our lattice ensembles are listed in Refs. [9,10]. The baryon masses can be extracted from the exponential decay of baryon correlators C b (t) at large Euclidean times 2 . We employ localized corner sources. To enhance statistics we average over sources living on different time-slices as well as at different spatial locations. In addition, a sum over spatial coordinates is performed at the sink (at B = 0 summing over the x and y coordinate components achieves zero momentum projection p x = p y = 0, while for B > 0 it merely helps to reduce fluctuations [34]). For staggered quarks, single-time-slice baryon operators mix parity partners so that the correlator takes the form [39], requiring a four-parameter fit to extract the mass of the baryon, M b , and of its parity partner, M b . We consider members of the baryon octet, including baryons with strangeness S = 0, −1 and −2. In the effective model study we will only use four of the baryons b = p, n, Σ 0 and Σ + , for reasons which will become clear later.
The determination of a baryon mass M b at a certain lattice spacing a and at a certain B is as follows: an effective mass (M eff b ) as a function of the fitting region (labeled by t min ) is obtained by fitting the function (2.2) to the measured correlator data in the region [t min , N t − t min ]. A plateau is then extrapolated in t min from the acquired M eff b (t min ) as the t min → ∞ limit of a simple exponential decay. The statistical error is then estimated with the jackknife method, while a systematic error is estimated from the exponential fit to find the plateau. The example of the proton effective mass is shown in Fig. 1 along with the exponential fits and the mass estimates obtained, including statistical and systematic errors.
The continuum limit is carried out in two steps. First at vanishing magnetic field the masses M b (B = 0) are extrapolated to a = 0, then separately only the magnetic field dependence M b (B)/M b (B = 0) is extrapolated to the continuum. The latter step requires interpolation for the magnetic field dependence, since at different lattice spacings we have measurements at different physical magnetic field values. We carry out the continuum limit of the B-dependence by fitting a Taylor-expansion 2 We note that in the present study we do not aim for precision results for the magnetic moments (related to the weak magnetic field-region), but concentrate on strong magnetic fields, which will be relevant for the phase diagram, see below. Thus we do not consider spin-projected operators but look for the state that minimizes the baryon energy. effective mass average fit, type 1 average fit, type 2 mass estimate 1 mass estimate 2 final mass estimate Figure 1. Effective mass diagram of the proton at vanishing magnetic field. The parameter t min characterizes the fitting region for (2.2), the larger it is the less early time correlator data is used. We estimate the t min → ∞ limit by fitting exponential decays, and the deviation of the estimates is used as the systematic error of our method. The grey band around the final mass estimate is the combined statistical and systematic error.
with lattice spacing dependent coefficients To estimate the systematic error of our approach we redo the fits excluding the B 6 term to see how much the result changes. The statistical errors are estimated both for the B = 0 and B = 0 cases by the bootstrap method. The continuum extrapolation of the nucleon and Σ masses at B = 0 is shown in Fig. 2, comparing to their respective experimental values. Notice that at zero magnetic field isospin symmetry is present, which is reflected in our results as well. In the S = −2 channel, large lattice artefacts, together with the closeness of excited states prevent us from reaching an acceptable continuum limit for the Ξ baryons. (For precision results at B = 0 including further baryons we refer the reader to Ref. [40].) The continuum limit of the magnetic field dependence of the remaining baryon masses is shown in Fig. 3. At low magnetic fields, a few outlier points are visible, related to the fact that here the Zeeman-splitting cannot be fully resolved. For strong magnetic fields this issue is absent. Notice furthermore that the behaviour of the Σ − is completely different compared to the others. This might be explained within a simplified quark model: for all other baryons, quarks can orient their magnetic moments in an energetically favorable way with respect to the magnetic field in the lowest energy configuration (i.e. in the lowest Landau-level), however in the case of the Σ − one of the quarks is forced to be in an excited state (first Landau-level).

Construction of the PNJL model
As an application for the B-dependent baryon masses, we use them as input for the magnetic field dependent reparameterization of the two-flavor PNJL model, which in turn will be used to explore     the B − T phase diagram of strongly interacting matter. First of all, since the PNJL model can only deal with constituent quark masses and not baryons, we use a simple non-relativistic quark model (NRQM) based on Ref. [41] to define B-dependent u, d and s constituent quark masses. For this reason it is also advantageous to discuss baryons instead of mesons -the latter receive their masses substantially from explicit chiral symmetry breaking and not from constituent quarks. Our working assumption is that the baryon masses can be obtained by merely summing the masses of their constituents: with M f being the constituent quark mass for flavor f . We determine M f as a function of B by a least squares fit of the three quark masses to the results shown in Fig. 3. According to Ref. [41], in the Σ − baryon at least one quark is forced into a spin state for which the Zeeman energy is added instead of subtracted. In the other four baryons, however, all quarks can be in the energetically most favorable spin state. To avoid having to describe excited states of the constituent quarks, we therefore disregard Σ − from the least squares fit. The goodness of the fits are found to be satisfactory, χ 2 < 1 for all magnetic fields. The obtained constituent quark masses are shown in Fig. 4. The errors are propagated by bootstrap resampling, while systematic errors of the NRQM model are estimated by redoing the fits leaving out one baryon at a time.
We now briefly summarize the basic properties and equations of the PNJL model following Ref. [42], except that we use Schwinger's proper time method as the ultraviolet regularization scheme, see e.g. Ref. [43]. Errors are propagated over from the constituent quark masses to all PNJL results by bootstrap resampling. The Lagrangian of the PNJL model is where ψ is the constituent quark field coupled to the Polyakov loop P through the covariant derivative and m 0 is the bare current quark mass. The Polyakov loop potential U(P, T ) is a classical one constructed to reproduce pure gluonic lattice results for the temperature dependence of the Polyakov loop expectation value [44], We adopt the usual choice of parameters a 0 = 3.51, a 1 = −2.47, a 2 = 15.2, b 3 = −1.75, except for T 0 , which sets the transition temperature in the pure gauge theory. It is usually set to 270 MeV, however -following Ref.
[45] -we set it to T 0 = 208 MeV to include corrections induced by the two quark flavors. We use the mean-field approximation for the quarks, in which the thermodynamic potential at finite B reads , (3.5) where M = M u + M d is the dynamically generated average constituent quark mass for the two flavors, In the mean-field approximation, both M and P minimizes the thermodynamic potential. We solve the model at every B and T by numerically searching for the two dimensional minimum of Ω. Once the minimum is found, the quark condensate can be obtained by evaluating Note that we follow the convention, where ψ ψ is positive and therefore our (3.8) contains an extra minus sign compared to most NJL studies. The potential Ω depends on three model parameters: the bare current quark mass m 0 , the fourfermion coupling G and the cutoff scale Λ of the theory (for different regularizations of the NJL model, see Ref. [46]). In mapping out the B − T phase diagram we first fix m 0 and Λ at B = T = 0 by setting the predictions of the NJL model for the pion mass m π , lattice result lattice-improved PNJL standard PNJL Figure 5.
Left: The magnetic field dependent coupling inferred from the baryon masses of Fig. 3. The significant deviation from a constant already signals a strong effect on this level. Right: The average quark condensate at T = 0 as a function of B compared to lattice results of Ref. [10] and to a standard PNJL calculation with B independent coupling. and for the pion decay constant f π , to their physical value, that is 138 MeV and 93 MeV, respectively, following Ref. [43]. To fully fix the parameters of the NJL model we prescribe M (B, T = 0) to take the value which is consistent with the average of the u and d constituent quark masses inferred from the baryon masses measured on the lattice for each B. This results in m 0 = 3.50(5) MeV and Λ = 675(10) MeV and G(B) shown in Fig. 5 (left). We find that the coupling constant inferred from lattice baryon masses strongly decreases with increasing magnetic field. This reinforces studies which hand tuned the coupling to a qualitatively similar function in order to achieve the correct T c behavior. As a consistency check, in the right panel of Fig. 5 we show a T = 0 consequence of including the B-dependence in the coupling G(B).
The magnetic field dependence of the quark condensate is compared with lattice QCD results from Ref. [10]. We see that our results are consistent with lattice results in a broad range of magnetic fields, while in a standard PNJL calculation where the coupling constant is a constant value (fixed to our G(B=0)) the two curves seem to diverge for larger magnetic fields.

Phase diagram
We now turn to the results on mapping out the B − T plane by minimizing the potential (3.5) with respect to both M and P using G(B). the standard PNJL we did not even plot the Polyakov-loop expectation value as it remains practically unchanged compared to the B = 0 curve.
We define two pseudo-critical temperatures: the inflection point of the quark condensate, which is identified with the chiral transition temperature, and the inflection point of the Polyakov loop, which in turn is typically associated to the deconfinement transition, however we will only discuss in detail the one obtained from the quark condensate now. In the left panel of Fig. 7 we show the quark condensate curves corresponding to the solutions shown in Fig. 6, which display inverse magnetic catalysis around the transition in the case of the lattice-improved PNJL model, however not in the standard PNJL model, where magnetic catalysis can be seen at all temperatures. In the right panel we compare the T c (B) curves with lattice results from Ref. [9] where we see that after rescaling with the respective T c (B = 0) values the lattice-improved result is consistent with the lattice continuum limit as opposed to the standrad PNJL result. The pseudo critical temperature at vanishing magnetic field in the lattice-improved PNJL model is T c (B = 0) = 204 (3) MeV. The deconfinement temperature defined from the Polyakov-loop seems to be lower (similar behaviour was found in [18]) but the disentanglement of the two transitions may need more in-depth analysis.

Summary
In this paper we performed the first lattice determination of the baryon spectrum in the presence of strong magnetic fields B at the physical point, including a continuum extrapolation. Using the B-dependence of the nucleon and Σ baryon masses, we defined constituent quark masses that were employed as zero-temperature inputs for the Polyakov loop-extended NJL model. The standard variant of this model is known to qualitatively fail in describing the QCD phase diagram in the magnetic field-temperature plane. We demonstrate that our lattice-improved PNJL model reproduces all features of the lattice findings at B > 0, including the inverse magnetic catalysis of the light quark condensate in the transition region as well as the reduction of the chiral crossover temperature by B. This result reveals that the model can be substantially improved if minimal information is fed to it at zero temperature -allowing it to capture the non-trivial dependence ofψψ(B, T ) in a broad range of magnetic fields and temperatures. An obvious extension of our results would be to include the strange quark flavor in the PNJL model or isospin splittings as well as further channels that may emerge at B > 0 [47]. The ideas presented in this work might also be generalized to other low-energy models of QCD.