Critical point in the QCD phase diagram for extremely strong background magnetic fields

Lattice simulations have demonstrated that a background (electro)magnetic field reduces the chiral/deconfinement transition temperature of quantum chromodynamics for eB<1 GeV^2. On the level of observables, this reduction manifests itself in an enhancement of the Polyakov loop and in a suppression of the light quark condensates (inverse magnetic catalysis) in the transition region. In this paper, we report on lattice simulations of 1+1+1-flavor QCD at an unprecedentedly high value of the magnetic field eB = 3.25 GeV^2. Based on the behavior of various observables, it is shown that even at this extremely strong field, inverse magnetic catalysis prevails and the transition, albeit becoming sharper, remains an analytic crossover. In addition, we develop an algorithm to directly simulate the asymptotically strong magnetic field limit of QCD. We find strong evidence for a first-order deconfinement phase transition in this limiting theory, implying the presence of a critical point in the QCD phase diagram. Based on the available lattice data, we estimate the location of the critical point.

system and the high degeneracy of the lowest Landau-level [9,10]. For low magnetic fields, it can be related to the positivity of the QED -function that fixes the dependence of the condensate on to order 2 [11,12]. Magnetic catalysis even has connections to solid state physics models like the Hofstadter model [13]. In line with these arguments, the catalysis of the condensate at low temperatures was observed in a variety of model settings and effective theories of QCD. However, in most of these models, magnetic catalysis takes place not only for < , but for all temperatures, giving rise to a monotonously increasing dependence of on . Thus, for the phase diagram, these models predict just the opposite of what the above discussed lattice results suggest. For a recent summary on these model approaches and a comparison to the lattice results, see Refs. [3,14].
While the mechanisms behind magnetic catalysis, as mentioned above, are quite transparent, the opposite behavior around -inverse magnetic catalysis -apparently has its origin in the rearrangement of the gluonic configurations that dominate the QCD path integral and is thus highly nontrivial [7]. Several attempts have been made recently to reproduce this behavior in effective approaches to QCD [15], among others, by introducing new, -dependent model parameters. Several of these models exhibit a non-monotonous ( ) dependence, with an initial reduction followed by an enhancement due to the magnetic field. In certain settings, it was even shown that no matter how the existing parameters of the model are tuned as functions of , the transition temperature always tends to rise above a given threshold magnetic field [16].
It is just the apparent universality of magnetic catalysis that has made the lattice results about inverse catalysis and the decreasing ( ) dependence for 0 ≤ < 1 GeV 2 so unexpected. It was speculated that magnetic catalysis should reappear at even stronger magnetic fields, and different hypotheses were recently put forward about the strong regime of the phase diagram. In particular, the strong limit was argued to induce a new critical point [17]. The transition temperature was conjectured to turn around and increase if the magnetic field is sufficiently strong [18][19][20][21]. In other cases, was argued to keep decreasing and to hit zero [22]. The transition temperatures for chiral restoration and for deconfinement were predicted to split [23] in the presence of the magnetic field, and a splitting between the chiral restoration temperature for the up and down quarks was also argued to take place [21,24]. Let us refer the reader to the recent reviews [3,14] for details.
In this paper, we aim to check these conjectures by means of first-principles lattice simulations of 1 + 1 + 1-flavor QCD at an unprecedentedly strong magnetic field = 3.25 GeV 2 . In addition, we also simulate the → ∞ limit directly, by considering the effective theory relevant for this limit [25]. We find strong evidence that this limiting theory has a first-order deconfinement phase transition and, thus, the QCD phase diagram exhibits a critical endpoint at strong magnetic fields, where the analytic crossover terminates. Based on our results, we estimate the location of the endpoint, and sketch the dependence of the deconfinement transition temperature on over a broad range. Besides answering a fundamental question about the QCD phase diagram, we believe that the results will also be useful for building/refining effective theories and models of QCD.
The paper is organized as follows. In Sec. 2 we discuss the details of the simulations and define the observables used to study the phase diagram. Sec. 3 contains the lattice results in ordinary QCD, followed by Sec. 4, where we discuss the simulations of the anisotropic theory in the asymptotic limit. The derivation of this effective theory and the employed simulation algorithms are discussed in the appendices App. A and App. B. Finally, in Sec. 5 we summarize our findings regarding the QCD phase diagram and conclude.

Setup and observables
We consider a spatially symmetric 3 × lattice with spacing so that the temperature is given by = ( ) −1 , the spatial volume by = ( ) 3 and the four-volume by 4 = / . Given this geometry, we simulate 1 + 1 + 1-flavor QCD, described by the partition function, given by the functional integral over the gluonic links . We employ stout smeared rooted staggered quarks described by the fermion matrix . In Eq. (2.1), ≡ 4 is the tree-level Symanzik improved gauge action and = 6/ 2 the inverse gauge coupling. For further details of the simulation setup and algorithm, see Refs. [5,26]. The parameters of the fermion matrix are the quark masses = ̸ = and the electric charges = −2 = −2 = 2 /3 ( > 0 is the elementary charge), which enter in the product with . The quark masses are set to their physical values along the line of constant physics [27]. The magnetic field is oriented along the positive -direction and has the quantized flux where we used that the smallest charge in the system is that of the down quark. Lattice discretization effects are suppressed as long as the flux quantum is much smaller than the period 2 . Previous experience suggests that < 2 /16 is a reasonable choice [5]. In the following, we employ the fixed-approach and change the temperature ( ) = ( ( )) −1 by varying the inverse gauge coupling. This also implies that a given flux quantum corresponds to different magnetic fields at different temperatures, i.e. ∝ 2 ( ). In particular, we choose values where a fixed magnetic field = 3.25 GeV 2 is represented by integer flux quanta. Although this implies that only discrete temperatures are allowed, at this strong magnetic field the temperature differences are small enough in order to map out the transition region (see below).
Next, we define the observables that can be used to pin down the transition temperature. We begin with the quark condensates and susceptibilities, signaling chiral symmetry, and employ the normalization inspired by the Gell-Mann-Oakes-Renner relation, introduced in Ref. [6] for the condensate,

(2.4)
Here, = 135 MeV is the pion mass and = 86 MeV the chiral limit of the pion decay constant. Both Σ and Σ are free of additive and of multiplicative divergences [5]. In addition, Σ is normalized to be unity at = = 0 and approaches zero above the transition region [6]. The vacuum values necessary for the additive renormalization were determined in Refs. [5,6].
The approximate order parameter for center symmetry, related to the deconfinement transition, is the Polyakov loop, defined on the lattice as In full QCD, the fermion determinant breaks center symmetry explicitly, so that the spontaneous breaking always occurs towards the real center element and ⟨Re ⟩ is a valid (approximate) order parameter. In pure gauge theory (this will be relevant for the → ∞ limit, see Sec. 4), there is no explicit breaking and the three center sectors are equivalent. In this case, it is convenient to consider the projection of the Polyakov loop to the nearest center element (see, e.g., Ref. [28]), Simulating pure gauge theory on a finite lattice, ⟨ ⟩ always vanishes due to the tunneling between center sectors, while ⟨ pr ⟩ is positive in the deconfined phase. The susceptibility of the projected Polyakov loop is defined as The Polyakov loop renormalizes multiplicatively, with a temperature-dependent renormalization constant which is determined by enforcing ⟨ ( ⋆ , 0)⟩ = ⋆ and we chose ⋆ = 162 MeV and ⋆ = 1. The renormalization was discussed in detail and ( ) was determined in Ref. [7]. Notice that while Re < 3 by construction, the renormalized observable has no upper bound. An observable that strongly correlates with -and, thus, is sensitive to the deconfinement transition -is the strange quark number susceptibility, Note that 2 contains neither additive nor multiplicative divergences. Finally, the trace anomaly can be written as a sum of gluonic, fermionic and magnetic contributions [12], is the lowest-order QED -function coefficient. The magnetic term appears due to electric charge renormalization and stems from the counter-term canceling the -dependent additive divergence of the thermodynamic potential log [12]. Note that this term is finite and independent of the regularization, once the continuum limit is taken, see discussion in Ref. [12]. Notice furthermore that the derivative in the definition of is evaluated at fixed magnetic flux Φ and not at fixed magnetic field [this is indicated by the superscript (Φ)]. The need for distinguishing between the two directional derivatives was first discussed in Ref. [29] and put into practice for the trace anomaly in Ref. [12].
3 Results at = 3.25 GeV 2 We extend the previously published data on the light condensates and susceptibilities [6], on the Polyakov loop [7], on the strange quark number susceptibility [5], and on the trace anomaly [12,29] using our new results at = 3.25 GeV 2 . To achieve this magnetic field strength, a temporal lattice extent = 16 turned out to be necessary. These = 16 lattices are finer than the finite temperature configurations used in Refs. [5][6][7] ( = 6, 8 and 10) to extrapolate to the continuum limit. Thus, our results -although not strictly continuum extrapolated -are expected to lie close to the limit → 0. We use two spatial lattice sizes 32 3 × 16 and 48 3 × 16 to control finite size effects.
Let us start the discussion with the light quark condensates. The average of Σ and Σ is plotted in the left panel of Fig. 1, compared to the = 0 and = 1 GeV 2 continuum extrapolated results [6]. In addition to the data at nonzero temperatures, we also indicate an estimate for the zero-temperature condensate. This is obtained by fitting and extrapolating the available lattice data at = 0 by a free-theory inspired form ∼ log . The systematic uncertainty is taken into account by varying the fit interval.
While the condensate is increased by the magnetic field at low temperatures, reflecting the wellknown magnetic catalysis effect, the results also clearly show the reduction of Σ + Σ in the transition region. Thus, inverse magnetic catalysis is observed to persist in the transition region even for our strong magnetic field = 3.25 GeV 2 , pushing further down. In particular, we employ the inflection point of the average condensate to find {Σ } = 112(3) MeV. As a side-remark, we mention that since the transition region is shifted to considerably lower temperatures, the vacuum values determined  . At the highest magnetic field, the curve is a spline interpolation, while for the lower fields the band is the result of a combined continuum extrapolation and interpolation in [7]. Right panel: the strange quark number susceptibility for the same set of magnetic fields. For the highest magnetic field, a spline interpolation is shown, whereas for the lower fields the bands represent a continuum estimate based on the results of Ref. [5].
for 3.45 < < 3.85 in Refs. [5,6] suffice to perform the additive renormalization of the condensates, and there is no need for additional = 0 simulations on finer lattices.
Due to the different electric charges, Σ is expected to be more sensitive to the magnetic field than Σ . On that account, even a splitting in the transition temperatures might seem plausible, see Refs. [21,24]. To check whether this is the case, in the right panel of Fig MeV. An apparent implication of this finding is that the temperature, at which the transition between the chirally broken and restored phases takes place, is encoded in the gluonic configurations rather than in the operator insertion. In lattice language; seems to be a quantity driven predominantly by sea and not by valence effects. This also suggests that purely gluonic observables would also exhibit similar transition temperatures.
This brings us to the simplest, purely gluonic quantity: the Polyakov loop (2.5). The (real part of the) renormalized observable is plotted in the left panel of Fig. 2, for the same set of magnetic fields, and is observed to be drastically enhanced by the magnetic field for all temperatures. The inflection point of ( ) is much more pronounced as compared to the case at = 1 GeV 2 and is determined to be { } = 109(3) MeV. This value is indeed consistent with the transition temperatures obtained above for the light quark condensates. We conclude that the gluonic configurations are vastly different on the two 'sides' of the transition, and predestine the behavior of the light condensates, independently of the electric charge that appears in the operator. We also observe the strange quark number susceptibility to exhibit an analogous trend, see the right panel of Fig. 2. Performing a similar fit as for , we obtain A further observable of interest for the QCD equation of state is the trace anomaly (2.10). It measures the breaking of conformal symmetry by the gluonic condensate, by the quark condensates and by the magnetic field itself. As grows, the latter effect becomes dominant and (Φ) is increased drastically, as visible in Fig. 3. Since (Φ) containsdependent contributions already at zero temperature, the usual normalization (Φ) / 4 is not useful [12]. This large = 0 contribution also damps the behavior of (Φ) in the transition region. The small kink around , moving towards smaller temperatures as grows, is to some extent still visible. We note that in order to determine further quantities related to the equation of state (e.g. pressure, entropy density etc.), one would need additional simulations at low temperature (see the method developed in Ref. [12]). This is outside the scope of the present paper.
Besides the characteristic temperature, the strength of the transition at high magnetic fields is also of interest. To determine, whether the smooth crossover at < 1 GeV 2 turns into a real phase transition at = 3.25 GeV 2 , we analyze the average of the light quark susceptibilities Σ + Σ . This observable exhibits a peak at the transition temperature, see the left panel of Fig. 4. For real phase transitions, the height ℎ of this peak diverges in the infinite volume limit: ℎ ∝ for first-order transitions and ℎ ∝ with a critical exponent < 1 for second-order transitions. In contrast, for the case of an analytic crossover, ℎ is independent of the volume. To perform this finite size scaling study, we compare the results obtained on the 48 3 × 16 and on the 32 3 × 16 ensembles. The left panel of Fig. 4 shows no sign of a singularity as is increased (note that for a first-order transition, the peak heights for the two volumes would differ by more than three). This leads us to conclude that the transition remains an analytic crossover even at = 3.25 GeV 2 . We mention moreover that finite volume effects are also absent from the other observables discussed above.
Although the transition remains an analytic crossover, it is instructive to analyze the -dependence of the susceptibility peak in more detail. We normalize Σ + Σ such that its peak maximum equals unity, and plot it in the right panel of Fig. 4 against the temperature. Here, is shifted so that the observable equals 0.5 at zero. Then, the peak width ( ) at half maximum can be read off at the rightmost intersection of the observable with 0.5. Clearly, ( ) decreases as grows, signaling that the transition becomes stronger in the presence of the magnetic field. We will return to this observation below in Sec. 5.
Besides being useful for quantifying the strength of the transition, the peak of the susceptibility allows for yet another determination of the transition temperature. Fitting for the peak maximum, we obtain { Σ } = 113(4) MeV, consistent with the results obtained for all other observables. We mention that the peak positions for the up and down quark susceptibilities also agree within errors.
Finally, in Fig. 5 we summarize our determinations of in the QCD phase diagram. We consider the results for < 1 GeV 2 obtained for the light quark condensates and for the strange quark number susceptibility [5]. In addition, we also include the transition temperatures at = 3.25 GeV 2 obtained using the light quark condensates, the strange quark number susceptibility and the Polyakov loop. (Note that the inflection point of at < 1 GeV 2 is not pronounced enough to allow for a stable fit.) To interpolate {Σ } and { 2 } for all magnetic fields, we found the following function sufficient, giving the fit parameters shown in Tab. 1. The resulting fit is also shown in the figure. The QCD phase diagram in the magnetic field-temperature plane. Previous results at weaker magnetic fields [5] are complemented by our findings at high . The points have been slightly shifted horizontally for better visibility. The dotted and the dashed lines show an interpolation of the results for Σ + Σ and for 2 , respectively, according to Eq. (3.1). (0)

Results in the asymptotic magnetic field limit
Our results in the right panel of Fig. 4 indicate that the transition becomes significantly sharper as the magnetic field increases. This observation raises the question: what happens if is even larger? Does the crossover terminate and turn into a real phase transition? To answer this question, we have to consider the limit ≫ Λ 2 QCD . Asymptotic freedom dictates that in this limit quarks and gluons decouple from each other. Still, the explicit breaking of rotational symmetry by and the corresponding dimensional reduction in the quark sector [9] suggests that this limit is not simply given by a pure gluonic theory plus non-interacting (electrically charged) quarks. Indeed, based on the structure of the gluon propagator in strong magnetic fields, Ref. [25] has shown that the effective action describing this limit is an anisotropic pure gauge theory. The anisotropy amounts to an enhancement of the chromo-dielectric constant in the direction parallel to the magnetic field, characterized by the coefficient , The definition of the gluonic field strength components ℬ and ℰ is given in Eq. (A.2) below. The enhancement of the parallel chromo-dielectric constant implies that the corresponding field strength component ℰ ‖ is suppressed. This tendency is already visible in our full QCD simulations at strong magnetic fields, see Fig. 11 in App. A below. Therefore, as is increased, the QCD effective Lagrangian approaches the anisotropic gauge theory given by Eq. (4.1). Assuming that this theory has a first-order phase transition, Ref. [17] has conjectured that the strong magnetic field region of the QCD phase diagram should exhibit a critical endpoint.
Here we address this question in more detail. First of all, in App. A, we reproduce the results of Ref. [25] for the magnetic field-induced anisotropy using the effective action in the Schwinger propertime formulation. The resulting anisotropic gauge theory can be simulated directly on the lattice. The setup and the simulation algorithm are described in App. B. The main difference to simple pure gauge theory amounts to multiplying the plaquettes lying in the − plane by the anisotropy coefficient . The exact form of the anisotropic action is given in Eq. (B.2).
Before discussing the lattice simulations of the anisotropic theory, let us make one more remark. Besides writing down the effective Lagrangian (4.1), Ref. [25] also predicted that the scale QCD of this theory (generated through dimensional transmutation) is much smaller than the QCD scale at = 0: QCD ≪ Λ QCD . In the absence of further dimensionful scales in the anisotropic theory, this implies that the deconfinement transition temperature is also much smaller, i.e. lim →∞ ≪ ( = 0). Below we will also address this prediction.
Due to the exact Z(3) symmetry of the anisotropic theory, the deconfinement transition is characterized by the projected Polyakov loop (2.6). In the left panel of Fig. 6, this observable is plotted against the inverse gauge coupling for several values of , as measured on the 16 3 × 4 lattices. At = 0, we reproduce the results of Refs. [30,31] -in particular, the deconfinement transition occurs at ≈ 4.07. The results indicate to be strongly reduced as grows 1 . We find that scales approximately with 1/ √ , see the right panel of Fig. 6. Extrapolating to = ∞ we obtain ( = ∞) ≈ 2.42(5). Besides approaching this limit via finite values of the anisotropy coefficient, we also develop an algorithm to simulate directly at = ∞. The corresponding setup is described in App. B. The left panel of Fig. 7 shows pr at infinite anisotropy. We find that the critical inverse coupling on the 16 3 × 4 lattice is comparable with the extrapolation based on finite anisotropies, see the right panel of Fig. 6. In addition, a comparison of the results at different spatial volumes 12 3 . . . 24 3 reveals that the transition becomes sharper as the volume increases, as typical for real phase transitions. We also repeated this analysis on = 8 lattices, see the right panel of Fig. 7. The critical couplings are clearly different, showing that the transition is indeed related to the finite temperature. We also mention that finite volume effects in are observed to be unusually large -above 10% for = 4. (For comparison, the finite volume effects at = 0 on similar lattices are of 0.1% [30].) We suspect that this is due to lattice artefacts -indeed, the effect is considerably smaller for = 8, see the right panel of Fig. 7. To determine the nature of the transition, we calculated the susceptibility of the projected Polyakov loop (2.7). This observable is shown in the left panel of Fig. 8, with a normalization by the spatial volume. Within statistical errors, the height of the normalized susceptibility peak is observed to be independent of . In other words, the peak height scales linearly with , which we take as strong evidence that the transition is of first order. The histogram of pr at = 2.4855 as measured on the 16 3 × 4 lattices is shown in the right panel of Fig. 8, revealing the two-peak structure characteristic for first-order transitions.
Through the equivalence between this anisotropic gauge theory and QCD with asymptotically strong magnetic fields, the above finding implies that the QCD phase diagram exhibits a critical endpoint in the strong magnetic field region, where the crossover turns into a real phase transition. Based on our full QCD results for the light quark susceptibilities, we will estimate the magnetic field corresponding to the endpoint in Sec. 5. The next step is to relate the critical inverse coupling to the critical temperature in physical units. To do so, we must set the lattice scale ( ). In principle, the magnetic field is not expected to change this scaling relation (cf. Ref. [5]). However, to arrive at our anisotropic gauge theory, has been taken to infinity, i.e. it also exceeds the squared lattice cutoff −2 . Clearly, the lattice scale determined at = 0 becomes invalid beyond this point. Thus, in order to determine the lattice spacing, one needs a dimensionful quantity whose value is known in the asymptotic limit -for example a purely gluonic observable, where the -dependence is expected to be only mild. A possible candidate for this role is the parameter 0 defined from the gradient flow of the gauge links [32] that is often used for scale setting in QCD, as suggested in Ref. [33].
We determined 0 on our zero-temperature full QCD ensembles [6] for < 1 GeV 2 and also for = 3.25 GeV 2 at our lowest temperature ≈ 75 MeV. The results for the ratio 0 / 0 ( = 0) are plotted in the left panel of Fig. 9, showing a mild reduction of this parameter as grows. A fit of the form similar to Eq. (3.1) describes the data well and suggests a saturation towards the asymptotically strong magnetic field limit. Nevertheless, we cannot exclude a significant dependence of 0 on for > 3.25 GeV 2 .
In addition, we can also gain some insight by considering the dimensionless combination 0 . How close full QCD at = 3.25 GeV 2 is to the asymptotic limit can then be quantified by matching 0 with the anisotropic theory. Multiplying our full QCD results for 0 by the transition temperature (here we take the definition of employing the inflection point of the strange quark number susceptibility, cf. Fig. 5), 0 is shown in the right panel of Fig. 9. Motivated by the scaling of (cf. the right panel of Fig. 6), the results are plotted against 1/ √ . Employing the result for 0 from Ref. [33], at zero magnetic field we have ( = 0) · 0 ( = 0) = 0.174(3) GeV · 0.1755(19) fm = 0.155(3). To carry out the comparison to the asymptotic limit, we also determined 0 / on symmetric 16 4 and 24 4 anisotropic gauge configurations at the critical couplings corresponding to the 16 3 × 4 ( ≈ 2.47) and to the 24 3 × 8 ( ≈ 2.7) lattices 2 . We observe that the combination 0 = 0 ( )/ · 1/ -similarly to -suffers from large lattice discretization effects and exhibits a downwards trend towards the continuum limit. We take the result for the = 8 data as an upper limit, giving lim →∞ ( 0 ) 0.076. This value is also included in the right panel of Fig. 9. Altogether, the results are compatible with a monotonous dependence of 0 on . To extrapolate 0 reliably to the continuum limit in the anisotropic theory requires further simulations on finer lattices and will be discussed in a forthcoming publication.
To summarize, the lattice results favor a saturation of 0 and a monotonous reduction of 0 as the limit → ∞ is approached. This suggests a monotonous reduction of ( ) towards the asymptotic limit. Nevertheless, based on the available findings, no final statement about lim →∞ can be made. Let us make one more remark about the = ∞ anisotropic theory. Since the parallel chromoelectric component tr ℰ 2 ‖ of the action vanishes, all plaquettes lying in the − plane are unity. This implies that all Wilson loops in this plane are trivial, and the static quark-antiquark potential ∝ log is independent of the distance. Accordingly, there is no force acting on quark-antiquark pairs if they are separated in the direction of the magnetic field, i.e. the string tension ‖ in this direction vanishes. This is in line with recent lattice determinations of the string tension in magnetic fields [34].

Conclusions
In this paper, we determined the nature and the characteristic temperature of the chiral/deconfinement transition of QCD at an extremely strong background magnetic field = 3.25 GeV 2 . The results for various observables consistently show that the transition temperature is further decreased compared to its value at lower magnetic fields. For the light quark condensates, the reduction of is due to the so-called inverse magnetic catalysis: between = 1 GeV 2 and = 3.25 GeV 2 , Σ and of Σ are significantly reduced in the transition region. At the same time, the condensates are enhanced by both for ≪ and for ≫ (the latter effect is small, since the condensate is suppressed at high temperatures).
Comparing the behavior of the up and down quark condensates and that of the Polyakov loop also revealed that there is no splitting between the transition temperatures for the individual flavors, neither is there significant difference between the chiral and the deconfinement transition temperatures. On the contrary, the different definitions of tend to approach each other as grows and at = 3.25 GeV 2 all observables exhibit a single transition temperature of around 109−112 MeV, see Fig. 5. Furthermore, we performed a finite size scaling analysis of the light quark susceptibilities, which has revealed that there is no singularity in the infinite volume limit and, thus, the transition remains an analytic crossover even at = 3.25 GeV 2 . In addition, we considered the asymptotically strong magnetic field limit, and simulated the corresponding effective theory on the lattice. This limiting effective theory -an anisotropic pure gauge theory -was found to exhibit a first-order deconfinement phase transition. Together with our findings above, this implies the existence of a critical endpoint in the QCD phase diagram. To provide a first estimate for the magnetic field CEP corresponding to the endpoint, let us return to our results about the width ( ) of the light quark susceptibilities. We have seen that the width is reduced as the magnetic field grows, see the right panel of Fig. 4. Assuming a linear dependence of on and extrapolating in the magnetic field we find that vanishes at CEP ≈ 10(2) GeV 2 . (5.1) In light of the fact that the -dependence of some of our observables (e.g. of and of 0 ) tends to flatten out as grows, this first estimate should rather be taken as a lower bound for CEP . We mention that in order to simulate with magnetic fields of strengths comparable to that in Eq. (5.1), lattices with 28 are required, out of reach for current computational resources. In the absence of a priori known dimensionful scales in the → ∞ system, we could not determine lim →∞ in physical units. Nevertheless, the deconfinement transition temperature of the anisotropic theory is expected to be much smaller than ( = 0) [25]. Our results for the combination 0 are compatible with this prediction and suggest a gradual reduction of the deconfinement transition temperature as is increased 3 . Taking these aspects into account, Fig. 10 represents a sketch of the deconfinement transition line in the QCD phase diagram for a broad range of magnetic fields.
The reader might wonder whether it is possible that the crossover at ≤ 3.25 GeV 2 and the first-order transition in the asymptotic limit are not connected by a single line. To see that this is not the case, note that by varying the anisotropy parameter , one can continuously deform the anisotropic theory to usual pure gauge theory, as was demonstrated in Fig. 6. Furthermore, the isotropic pure gauge theory can be thought of as QCD with infinitely heavy quarks and thus can be continuously transformed into full QCD by increasing the inverse quark masses from zero to their physical values. Thus, the transition we identified at → ∞ is indeed the same deconfinement transition that occurs Figure 10. The deconfinement transition temperature against the background magnetic field. The results of our full lattice QCD simulations (white background) are complemented by the prediction (gray background) based on the results corresponding to the → ∞ limit and on the extrapolation of the light quark susceptibility peak to high magnetic fields (see the text).
at low magnetic fields.
Let us highlight that according to this discussion, having a decreasing deconfinement transition temperature is actually natural to QCD. Furthermore, since the → ∞ limit is independent of the quark masses 4 , a similar reduction of by the magnetic field should also take place in QCD with heavier-than-physical quarks. However, in the latter case this reduction most probably follows an initial increase in the transition temperature, cf. Refs. [5,35]. Indeed, recent lattice results employing overlap fermions and pion masses of about 500 MeV indicate inverse catalysis to occur around the transition temperature at the magnetic field ≈ 1.3 GeV 2 [8]. Finally, we note that magnetic fields well above the strength (5.1) are predicted to be generated during the electroweak phase transition in the early universe [36]. If these fields remain strong enough until the QCD epoch, the emerging first-order phase transition might have several exciting consequences. Via supercooling, bubbles of the confined phase can be formed as the temperature drops below , leading to large inhomogeneities, important for nucleosynthesis [37]. Collisions between the bubbles can also lead to the emission of gravitational waves and, thus, leave an imprint on the primordial gravitational spectrum [38]. An absence of such signals, in turn, would imply an upper limit for the strength of the primordial magnetic fields.

space-time is
and the quark determinant will be regularized using Schwinger's proper time formulation [39]. Since the electromagnetic field exceeds all scales in the system and in particular, ( ) 2 ≫ tr 2 , we may approximate the the chromo-fields in the fermionic action to be weak. In addition, we assume the chromo-fields to be covariantly constant, = 0 to enable a fully analytical treatment of the problem. Given this condition, the field strength can be gauge transformed to be constant in spacetime and diagonal in color space [40], = diag( ) with the color index = 1, 2, 3. Let us decompose the chromo-fields to chromomagnetic/chromoelectric components, parallel or perpendicular to the electromagnetic field . The leading terms in the strong -expansion are quadratic in the chromo-fields and thus, to find the coefficients of the respective components, it suffices to consider separately the effect of and ℬ ‖ , and ℬ ⊥ , and ℰ ‖ and and ℰ ⊥ . The effective Lagrangian for these components for small background magnetic fields was considered in Ref. [29]. Let us first take the case of and ℬ ‖ . For each flavor we may choose our coordinate system such that is positive. Then, each color component experiences a total (positive) magnetic field +ℬ ‖ so that Since ℬ ‖ only appears in the sum with , the effective Lagrangian becomes independent of the chromomagnetic field in the limit ≫ ℬ ‖ . This implies that quarks become insensitive to ℬ ‖ , i.e. decouple from this gluonic component.
Next we take the case with and ℬ ⊥ . Rotating our coordinate axes for each color component such that the axis points in the direction of the total magnetic field we get (A.4) In the strong limit, this becomes independent of ℬ ⊥ , signaling that quarks decouple from the perpendicular chromomagnetic component of the gluons as well.
For a perpendicular chromoelectric field, for each color component we can perform the Lorentz transformation that eliminates the electric field and, in turn, gives a total magnetic field √︁ ( ) 2 + ℰ 2 ⊥ . The corresponding effective Lagrangian equals Eq. (A.4), but with ℬ ⊥ replaced by ℰ ⊥ . This implies the decoupling of quarks from the perpendicular chromoelectric fields.
Finally, for a parallel chromoelectric field ℰ ‖ , we have a Landau problem in the − as well as in the − planes, giving Taking the limit ≫ ℰ ‖ , 2 , we see that -unlike for the other components above -a non-trivial dependence on ℰ ‖ remains. We are interested in the quadratic term, proportional to tr ℰ 2 ‖ , which contributes 5 to the gluonic Lagrangian tr 2 of Eq. (A.1), Altogether, the asymptotically strong magnetic field limit of the QCD effective Lagrangian indeed equals Eq. (4.1). Thus we find that the chromo-dielectric constant is enhanced in the direction of the background magnetic field, and the coefficient ( ) coincides with the result of Ref. [25]. Having ≫ 1 in the action implies that the corresponding gluonic field strength component tr ℰ 2 ‖ is strongly suppressed. In other words, the anisotropy in the chromoelectric part of the action density, is enhanced by . To back up this prediction, in Fig. 11 we plot (ℰ) as a function of the magnetic field, based on our zero-temperature results at < 1 GeV 2 [29] and the measurements at = 3.25 GeV 2 at our lowest temperature ≈ 75 MeV. The results clearly indicate that the anisotropy is positive and strongly increased as grows. We note that due to the coupling between the gluonic field strength components, a similar anisotropy in the chromomagnetic sector also appears, altogether giving rise to the hierarchy tr ℬ 2 ‖ > tr ℬ 2 ⊥ = tr ℰ 2 ⊥ > tr ℰ 2 ‖ at low temperatures. The same hierarchy is also observed in the anisotropic gauge theory 6 .
We note that the calculation leading to Eq. (4.1) can also be generalized to incorporate nonzero temperatures. At > 0 an additional factor appears in the proper time integral due to the sum over Matsubara frequencies, containing an elliptic Θ-function. This factor decouples from the -dependence, implying that even for > 0, only the chromo-dielectric constant is affected. The coefficient is, 5 Here we omitted a divergent term of the form tr ℰ 2 ‖ log Λ/ , where Λ is a cutoff entering as the lower endpoint of the proper time integration 0 ∝ 1/Λ 2 . This divergence can be eliminated by the multiplicative renormalization of the wave function ℰ ‖ and of the gauge coupling [39]. Closer inspection of Eqs. (A.3) and (A.4) shows that the same type of divergence is present for the other components as well. Thus, these -independent terms merely represent an isotropic redefinition of the gauge coupling , which does not alter the form of the effective Lagrangian for strong magnetic fields. Another divergence, independent of the gluonic field strengths, takes the form ( ) 2 log Λ/ and is canceled by the renormalization of and of [39]. Thus, the necessary renormalizations at → ∞ are of the same type as for the theory at small magnetic fields. 6 To see how the anisotropic dielectric constant affects the chromomagnetic components, it is instructive to consider the gauge potential . A large value of implies a suppression of tr ℰ 2 ‖ and a corresponding suppression of the fluctuations in and in . This suppression propagates into the magnetic sector and creates the anisotropy between tr ℬ 2 ⊥ and tr ℬ 2 ‖ . Indeed, while the former contains , the latter does not.
however, altered as The integral over equals unity at = 0 and is reduced monotonously (and smoothly) as the temperature grows. Simulating the anisotropic gauge theory according to the Lagrangian (4.1) on the lattice, we found that the theory exhibits a first-order phase transition. Thus, since the smooth ( ) dependence does not affect the discontinuous transition, in order to locate the critical temperature it suffices to simulate the theory at fixed (large) values.

B Simulating anisotropic pure gauge theory on the lattice
In this appendix we discuss the simulation algorithm for the anisotropic pure gauge theory described by the Lagrangian (4.1). The corresponding path integral can be simulated directly on the lattice. Here, denotes the gauge links, = 6/ 2 is the inverse gauge coupling and the anisotropic gauge action reads where are linear combinations of closed loops lying in the − plane. We take the tree-level Symanzik improved gauge action such that these loops include the 1 × 1 plaquettes 1×1 and the 2 × 1 rectangles 2×1 with appropriately tuned coefficients [41], The correspondence between the continuum and lattice expressions reads To simulate this theory, we use an overrelaxation/heatbath algorithm, based on the isotropic pure gauge implementation by the MILC collaboration [42]. One trajectory consists of one overrelaxation step followed by four heatbath steps. The simulation at finite anisotropy coefficient simply involves multiplying the plaquettes and rectangles lying in the − plane by . We observe that autocorrelation times grow large as increases, similarly to the issue of critical slowing down of the isotropic theory at large . This prohibits approaching → ∞, necessary for the asymptotically strong magnetic field limit. However, it is possible to modify the algorithm to simulate directly at = ∞. In this limit, the Indeed, any fluctuation in the link variables that leads off of this subspace makes the action infinitely large and is thus forbidden. (Note that if the 1 × 1 plaquette equals the unit matrix, then so does 2×1 .) We thus have to parameterize the subspace Ω[ ] in terms of the gauge links . Let us label the lattice sites by = ( , , , ) with 0 ≤ < . To find the parameterization of Eq. (B.5), it is advantageous to fix the links to 1 on a so-called maximal tree. The specific choice for the tree is shown in the left panel of Fig. 12. (Note that Faddeev-Popov fields are absent for such a gauge fixing [43].) In order to have unit plaquettes for < − 1 and < − 1, all -links must be set to unity at these sites. To have unit plaquettes on the last -slice, all -links at − 1 must be set equal, denoted by . Similarly, all the -links at − 1 must be set equal, denoted by , see the visualization in the right panel of Fig. 12. These remaining links correspond to the local Polyakov loops in the -and in the -direction, lying in the − plane at a given and . Finally, to ensure that the plaquette at the corner = − 1, = − 1 is unity, we need † † = 1, i.e. and must commute. Altogether, the subspace in question reads Notice that the 'degenerate' timelike Polyakov loop (represented by the blue arrows in the right panel of Fig. 12) appears multiple times in the action -in fact, in 4 plaquettes and in 12 rectangles. The corresponding 'staples' are all taken into account in the update of (and similarly for ). There are several ways to fulfill the commutativity relation [ , ] = 0. One possibility (setup A) is to simply set = 1 for all and . Another approach (setup B) is to constrain to be a center element, ( , ) ∈ Z 3 . The two setups only differ on a set whose measure vanishes in the limit, where all lattice extents are taken to infinity. Note that the expectation value of the average -Polyakov loop ( ) [defined similarly as the usual Polyakov loop , Eq. (2.5)] is three for setup A, whereas it is zero for setup B, if the spatial size of the system is large enough. Nevertheless, we checked that observables sensitive to the finite temperature transition ( , the gauge action, etc.) all have vanishing correlators with ( ) . In fact, we found that the setups A and B give identical results for and for the gauge action for all values of the inverse gauge coupling on the 16 3 × 4 lattices. In other words, center symmetry breaking in the direction appears to be completely irrelevant for the deconfinement phase transition.
In addition, we also tried allowing both and to be general SU(3) matrices (which violates the commutativity relation). This approach (setup C) turned out to introduce negligible differences in the results 7 . In the following, we consider setup A and set = 1 throughout the lattice. As a consistency check, besides the overrelaxation/heatbath algorithm, we also considered a hybrid Monte-Carlo update and found that the two give fully consistent results.