Confining and chiral properties of QCD in extremely strong magnetic fields

We investigate, by numerical lattice simulations, the static quark-antiquark potential, the flux tube properties and the chiral condensate for $N_f = 2+1$ QCD with physical quark masses in the presence of strong magnetic fields, going up to $eB = 9$ GeV$^2$, with continuum extrapolated results. The string tension for quark-antiquark separations longitudinal to the magnetic field is suppressed by one order of magnitude at the largest explored magnetic field with respect to its value at zero magnetic background, but is still non-vanishing; in the transverse direction, instead, the string tension is enhanced but seems to reach a saturation at around 50 % of its value at $B = 0$. The flux tube shows a consistent suppression/enhancement of the overall amplitude, with mild modifications of its profile. Finally, we observe magnetic catalysis in the whole range of explored fields with a behavior compatible with a lowest Landau level approximation, in particular with a linear dependence of the chiral condensate on $B$ which is in agreement, within errors, with that already observed for $eB \sim 1$ GeV$^2$.


I. INTRODUCTION
In the recent past, various analytic and numerical studies have uncovered a plenty of interesting new phenomena regarding the non-perturbative properties of strong interactions in the presence of a magnetic background field . Some of these phenomena might be of direct phenomenological relevance for heavy ion experiments [68][69][70][71][72] or astrophysics [73][74][75], some of them are more speculative but nevertheless interesting. Among these phenomena, a direct impact on the QCD vacuum properties, in particular those regarding the pure gauge sector, is particularly striking, since gluons are not electrically charged, and might be the signal of a stronger impact of the magnetic field on the QCD phase structure.
In particular, the conclusions of Ref. [41] pointed to the possible presence of a critical magnetic field eB 4 GeV 2 , above which the longitudinal string tension * Electronic address: massimo.delia@unipi.it † Electronic address: lorenzo.maio@phd.unipi.it ‡ Electronic address: a.stanzione1@studenti.unipi.it § Electronic address: francesco.sanfilippo@infn.it would vanish, resulting in a different and yet unknown phase of strongly interacting matter. Such conclusions, however, were not based on direct simulations performed at such large values of the magnetic background, but just on the extrapolation of results obtained in a smaller magnetic field range.
Given the new and interesting predicted phenomena, a direct investigation is of utmost importance. As we will better explain in the following Sections, the main difficulty in studying large magnetic backgrounds by lattice simulations is that the ultraviolet (UV) cut-off must be tuned correspondingly in order to keep discretization errors under control and allow for a reliable continuum extrapolation. In this study we will investigate N f = 2 + 1 QCD with physical quark masses and lattice spacings down to a ≃ 0.057 fm, which is around half the finest spacing explored in Ref. [41], with a similar discretization based on stout-improved staggered fermions. That will allow us to obtain continuum extrapolated results for eB up to ∼ 10 GeV 2 , which is enough to confirm or update the prediction of Ref. [41]. In addition to that, we will consider the chiral properties of the theory, in particular the chiral condensate, to investigate if magnetic catalysis is still at work.
The paper is organized as follows. In Section II we provide more details regarding the adopted discretization of N f = 2 + 1 QCD in the presence of a magnetic background, as well as about the lattice observable used to extract the potential and the chromoelectric field between the static quark-antiquark pair. In Section III, after some preliminary details regarding our numerical simulations, we illustrate our results for the chiral condensate, the static potential and the color flux tube. Fi-nally, in Section IV, we summarize our conclusions.

II. NUMERICAL METHODS
We consider N f = 2 + 1 QCD in the presence of a uniform and constant, external magnetic field, discretized in terms of the tree-level improved Symanzik action [88,89] for the gauge sector, and of rooted staggered fermions with stout improvement [90,91] for the fermionic sector. The resulting partition function is where [DU ] is the SU (3) group invariant integration measure on the link variables, f the flavor index, is the lattice gauge action, and is the discretized Dirac operator. There, i labels lattice sites and µ the direction, while β is the inverse gauge coupling and a the lattice spacing. The W 1×· i,µν s are the real parts of the trace of the links products along the 1 × 1 and 1 × 2 rectangular closed path, respectively. The η i;ν are the staggered quark phases, and U (2) i;ν is the two times stout smeared link (with isotropic smearing parameter ρ = 0.15).
An external electromagnetic (e.m.) field is added by minimal substitution in the covariant derivative where A a µ (x) are the gluon fields, T a the SU (3) generators, A µ (x) the abelian four-potential, and g 0 and q f are respectively the bare strong coupling constant and the quark electric charge. We consider for simplicity a uniform magnetic field B in theẑ direction, a possible gauge choice is then: That can be discretized on a periodic toroidal lattice in terms of U (1) link variables as follows where L x is the lattice extension along the x direction (in lattice units) and the right hand side (RHS) condition is needed to guarantee smoothness of the magnetic field across the periodic boundaries [92,93]; according to Eq. (5), all other abelian links are set to 1. Notice that, consistently with the required zero net magnetic flux across the lattice torus, the above U (1) links lead to a constant magnetic field but for a single plaquette, which is pierced by an additional Dirac string: invisibility of that string leads to a quantization condition for the magnetic field [92][93][94][95] considering that the smallest quark charge is e/3. The external field is finally added to the discretization of N f = 2 + 1 QCD decribed above by the following substitution in the Dirac operator in Eq. (3): Notice that in this approach the e.m. field is treated as purely external, neglecting the back-reaction of quarks on it, and meaning in practice that no additional integration over the U (1) gauge links is introduced in the partition function.
Bare masses and gauge coupling values have been set in order to move on a line of constant physics, determined in Refs. [96][97][98] to reproduce experimental results for hadronic observables at zero temperature in the continuum limit. The introduction of the external field leads to additional, B-dependent discretization errors. In particular, one should consider that the magnetic field acts in practice through the gauge invariant U (1) phase factors that dynamical quarks pick going through closed loops on the lattice: the smallest non-trivial such loop is the plaquette in the xy plane, for which the phase factor is Systematic errors in the discretization of the magnetic field are under control if such phase is much smaller than 2π: for the up quark, which has the largest electric charge q u = 2e/3, the condition reads: a useful way to visualize such systematics is to think that the existence of this minimal phase pickable by dynamical up quarks is like saying that we are approximating a circle by a regular polygon with ∼ L x L y /(2b z ) sides. All that also sets a natural UV cut-off for the largest magnetic fields which are explorable for a given lattice spacing, which is roughly eB ≤ 2π/a 2 . We have spent a few additional words on these aspects, since this will be essential to properly discuss discretization effects in our investigation, where extremely strong magnetic fields (one order of magnitude larger than the standard QCD scale) are considered.

A. Observables
To get the static potential of a qq-pair, similarly to Refs. [40,41], we studied the Wilson loop TrW (a n, an t ) and its dependence on the Euclidean time an t , exploiting the relation TrW (a n, an t ) ∝ e −aV (a n)nt , which holds for large enough an t . In particular, from previous equation one can derive aV (a n) = lim nt→∞ log TrW (a n, an t ) TrW (a n, a(n t + 1)) , (12) so that the potential at fixed n can be obtained by fitting to a constant the log in the RHS of Eq. (12) as a function of n t , at least in a suitable stability range. The color flux tube instead was studied, following Ref. [42], by means of the connected correlator-probe scheme [99][100][101][102][103][104][105]: the observable computed to derive chromoelectric field in-between the static quark-antiquark pair is in this case: where W is the open Wilson loop, P µt is the open plaquette in the µt-plane and L the Schwinger path linking the former two operators, while x t is the distance between the plaquette and Wilson loop plane, i.e. the distance from the quark-antiquark axis. One can easily prove that, in the naive continuum limit, which can be considered as a probe of the color field strength induced by the presence of the quark-antiquark pair, i.e., with some abuse of notation, as a 2 g 0 F µt QQ . The connected correlator ρ conn is pictorially described in Fig. 1  I: Equivalence classes of the relative orientations between the quark-antiquark pair, the magnetic field (fixed along theẑ direction) and the transverse direction xt. Labelŝ µ andρ refer to Fig. 1. Letters T and L stand for transverse or longitudinal with regard to the magnetic field. the midpoint of its temporal extent, it reaches half the distance between the quark-antiquark pair and then it moves x t lattice spacings in one of the directions orthogonal to the plane of the Wilson loop: in this way the flux tube profile is determined at the midpoint of the static color source. In this case, as in Ref. [42], the investigation has been limited to squared Wilson loops.
It is important to note that the presence of the background magnetic field alongẑ breaks the spatial octahedral symmetry, leaving a D 4 symmetry on the xy-plane. In the evaluation of the static potential by means of the Wilson loop, this condition implies that loops in the zdirection are not equivalent to those extending in x-and y-directions, which on the other hand are equivalent to each other: that naturally leads to distinguish between a static potential measured longitudinally (L) to the magnetic field, or transverse (T) to it. Actually, one can consider also generic angles between the magnetic field direction and the quark-antiquark axis: this is best done by considering magnetic fields with a generic orientation relative to the lattice axes. This kind of more general analysis has been performed in Ref. [41], showing however that most of the angular dependence of the static potential can be accounted for by the lowest harmonic. That means that the relevant information is contained in the L and T-potentials, which are therefore the only cases considered in the present investigation.
The study of the correlator ρ µt conn is a bit more involved, due to its three-dimensional shape. Apart from the Tor L-cases characterizing the orientation of the quarkantiquark separation relative to the magnetic field, in the T-case one can further distinguish whether the flux tube profile is studied in the direction parallel or orthogonal to B: the analysis of Ref. [42] shows that some minor anisotropies emerge also in this case, i.e. the flux tube itself loses its axial symmetry. Therefore, as for the flux tube profile, we will consider three different cases: L, TL and TT. We denote byμ the direction of the quarkantiquark axis, so that the possible geometries can be mapped into these three equivalence classes according to Table I. A residual symmetry x t → −x t is preserved but, anyway, it has not been exploited in this work.
For the evaluation of both observables in order to re-duce the UV noise, we applied one step of HYP smearing [106] for temporal links, with the following choice of parameters: α 1 = α 2 = 2α 3 = 1 (as for the HYP2-action defined in Ref. [107]). Moreover, we performed several steps of (spatial) APE smearing [108] on the spatial links, so that the N -times smeared link reads where is the sum of the spatial staples around the link U denotes the projection on the gauge group and the choices for α AP E match those of previous works: 0.25 for the string tension as in [41] and 1/6 for the QCD flux tube as in [42].

III. NUMERICAL RESULTS
Our results are based on simulations corresponding to three different values of eB (0, 4 and 9 GeV 2 ) and three different lattice spacings in each case (a ≃ 0.057, 0.086 and 0.114 fm) in order to allow for a continuum extrapolation; simulations at zero magnetic field have been performed mostly for renormalization purposes. The spatial lattice size has been kept fixed in most cases to aL s ∼ 2.75 fm, with an Euclidean temporal extent twice as large. A summary of all simulation points is reported in Table II. Monte-Carlo sampling of gauge configurations has been performed based on a Rational Hybrid Monte-Carlo (RHMC) algorithm running on GPUs [109,110]. For each simulation we performed O(10 3 ) RHMC steps, taking measures every 10 unit trajectories. The statistical analysis has been based in most cases on a binned bootstrap analysis.
In our analysis we will assume the lattice spacing being independent of eB. This is a reasonable assumption as long as the magnetic field is much smaller than the UV cutoff, hence it is expected to lead to sensible results at least when continuum extrapolations are considered. As a matter of fact, the assumption has been explicitly checked only for smaller values of the magnetic field [21]; however our analysis of the chiral condensate, leading to results in agreement with theoretical expectations, will further support the hypothesis.   [22]. The colored band is the result of a linear fit to the data of Ref. [22] (χ 2 /d.o.f. ≃ 1/3), which after extrapolation turns out to be in nice agreement, within errors, with our present determinations.

A. Magnetic catalysis at extremely large magnetic fields
Before focusing on the confining aspects of the theory, let us discuss its chiral properties, for which predictions are well established. In particular one expects, at least for T = 0, the magnetic catalysis phenomenon, with an enhancement of chiral symmetry breaking induced by the magnetic background field and detectable as an increase of the chiral condensate.
In order to compare with previous results in the literature, we will consider the change of the light quark condensate due to the magnetic field, renormalized as in Ref. [22]: with q = u, d, where ψ ψ q is determined as usual in terms of the volume normalized trace of the inverse fermion matrix (computed by noisy estimators), m π = 135 MeV is the pion mass and F π = 86 MeV is the pion decay constant in the chiral limit.
The average quantity (∆Σ u + ∆Σ d )/2 is displayed in Fig. 2 for the two explored values of eB as a function of the squared lattice spacing a 2 , together with a continuum extrapolation obtained assuming O(a 2 ) corrections. It is interesting to notice that continuum corrections are significantly larger for eB = 9 GeV 2 : that can be easily understood in terms of what discussed above regarding the B-dependent discretization errors (see Eq. (10) and comments thereafter). For eB = 4 GeV 2 , corresponding to b z = 41, these kind of discretization errors, at the three different lattice spacings, can be put in analogy with those that one has by approximating a circle by a regular polygon with respectively (from coarsest to finest) 7, 12, 28 sides, which is reasonable right from the beginning; for eB = 9 GeV 2 instead, corresponding to b z = 93, the approximation starts with a triangle (which is far from good) and ends with a dodecagon (which is reasonable). These considerations make it clearer why, having in mind to perform a reliable continuum extrapolation, it is not reasonable to consider larger values of eB, unless smaller lattice spacings are computationally affordable.
In Fig. 3 continuum extrapolated results are compared to the analogous ones obtained, for eB ≤ 1 GeV 2 , in Ref. [22]. The large field behavior of the magnetic catalysis phenomenon can be thoretically predicted in terms of a lowest Landau level (LLL) approximation. The higher energy levels increase proportionally to √ eB, thus they become practically irrelevant to the dynamics of the system. The LLL is instead independent of eB, while its degeneracy linearly increases with it: that leads to predict a linear behavior in the density of near-zero modes, hence in the chiral condensate by the Banks-Casher relation. This linear behavior is nicely reproduced in Fig. 3: in particular, in the figure we display the result of a linear fit to data from Ref. [22] which, when extrapolated to the large magnetic fields explored in this study, is perfectly compatible with our results within errors.
Notice that the lattice spacing enters with a fourth power in fixing the renormalization group invariant quantity in Eq. (16): hence, the nice consistency with data from Ref. [22] and with the LLL prediction supports the assumption that the lattice spacing is indeed independent of B.

B. Static quark-antiquark potential
The static quark-antiquark potential has been derived as described above, by looking for a plateau, as a function of the temporal extent n t , for the logarithm of Wilson loop ratios reported in Eq. (12). An example is showed in Fig. 4, where we report data obtained for both the trasverse and longitudinal direction at eB = 9 GeV 2 and R = 5 for the finest lattice spacing. The two bands show our final determination of the potential for the two cases and have been obtained considering Wilson loops after 30 spatial APE smearing steps, however we report in the figure also data obtained after 10 and 20 smearing steps, which are practically indistinguishable. A similar stability under APE smearing is observed for all explored values of eB, lattice spacing and quark-antiquark separation R.
In Fig. 5 we show the final determination of the static potential obtained for the finest lattice spacing and all the explored values of eB. The anisotropy which is present when eB = 0 is clearly evident, even if in the transverse direction it shows a non-monotonic behavior with eB, with a tendency for a slight decrease of the slope when going from 4 to 9 GeV 2 , at least for this value of the lattice spacing.
As a preliminary analysis, we have considered results for the potential at eB = 0 and compared them with previous results in the literature, in order to check consistency. In particular, in Fig. 6 we compare results obtained for the string tension in this work with those obtained in Ref. [41] using the same lattice discretization but different lattice spacings. The two sets of results are perfectly compatible with each other and a combined continuum extrapolation assuming O(a 2 ) corrections returns a continuum value χ 2 /d.o.f. = 1.3/3 (the point on the coarsest lattice being discarded), which is perfectly compatible with phenomenological predictions and lattice determinations for σ [111]. We would like to stress the importance of this consistency check: since the physical spatial lattice sizes adopted in this work and in Ref. [41] are, for computational reasons, quite different (∼ 3 fm vs ∼ 5 fm), the agreement we find shows that, at least for what concerns the static quark-antiquark potential, finite size effects are not significant. Next, in order to assess to the fate of the static potential anisotropy in the explored range of magnetic fields, we consider, as in Ref. [41], the dimensionless ratios σ(eB)/σ(0), which are reported in Fig. 7 as a function of a 2 for both magnetic fields and for both the longitudinal and the transverse directions, together with continuum extrapolations assuming O(a 2 ) corrections. Regarding the string tension in the longitudinal direction, we con-  The dashed gray regions correspond to the continuum extrapolations of Ref. [41] firm the findings of Ref. [41]: it is a consistently decreasing function of eB, both for at finite lattice spacing and in the continuum limit; however, contrary to the hypothesis put forward in Ref. [41], we find a non-zero string tension, within two standard deviations, even at the largest explored value of eB. Regarding the transverse direction, instead, we can appreciate from Fig. 7 that, at least for finite lattice spacing, the trend for an increasing string tension observed in Ref. [41] seems inverted. However, when considering continuum extrapolations, one realizes that σ(eB)/σ(0) actually reaches a saturation at large eB.
Such results are better appreciable in Fig. 8, where the continuum extrapolated values for σ(B)/σ(0) in the T-and L-directions obtained in this study are compared with the continuum extrapolation of Ref. [41], which is plotted only in the relevant range of eB where simulations  Table I for details). The physical distance between the color charges is d = 0.68 fm.
of Ref. [41] were performed. Present results are not inconsistent with those of Ref. [41], however they clarify the perspective for the large-eB limit of the string tension. In the trasverse direction, the string tension seems to reach a saturation at a value which is around 50 % larger than the zero field value. In the longitudinal direction the string tension keeps decreasing as a function of eB, but is still significantly different from zero for eB ∼ 4 GeV 2 , contrary to what the continuum extrapolation of Ref. [41] could have suggested: the large field behavior is difficult to predict precisely, it could be either an exponential-like decreasing behavior for which the string tension never vanishes, or a more straight decrease where σ L finally vanishing at some critical value for eB 10 GeV 2 . Such possibility should be further explored by future studies, capable of approaching even smaller values of the lattice spacing, which has been the main constraint limiting our simulations to eB 10 GeV 2 in order to access properly extrapolated continuum results.

C. Color flux tubes
The analysis of color flux tubes is expected, in general, to confirm results obtained by the analysis of the static quark-antiquark potential: this is indeed the outcome of Ref. [42], showing that the main effect of the magnetic field is an overall suppression/enhancement of the flux tube in the L/T directions, with a slight modification of its profile, which however can be still nicely described by models inspired to dual superconductivity of the QCD vacuum.
We measured the longitudinal component E l of the chromo-electric field (directed along the quark-antiquark axis) since previous studies showed that it is by far the dominant one (see [105] and references therein). We de-note byμ the direction where the color charges lay, so that the longitudinal chromo-electric field is given by where d is the separation distance between the quarkantiquark pair and x t is the transverse distance at which the field is probed (see Fig. 1). The use of smearing techniques introduces a non trivial dependence on the amount of smearing adopted. On the other hand, the analysis of Ref. [42] showed that ratios of observables with and without the magnetic field, such as E l (B)/E l (B = 0), are insensitive to the number of smearing steps. We verified that this feature holds true for the extreme magnetic fields investigated in this study, as illustrated in Fig. 9, where we show the ratios of the chromo-electric fields obtained at a = 0.0572 fm for each magnetic field choice and inequivalent class of Table I. We display the flux tubes at x t = 0 but, anyway, the independence is observed for each value of the transverse distance. Since we are interested in similar ratios of observables, we unambigously decided to fix N AP E = 30 for the analysis carried out for each lattice spacing and value of x t , so that the dependence on N AP E will be dropped in the following.
In Fig. 10 and Fig. 11 we show some results for the flux tube extracted from simulations performed at lattice spacings a = 0.0858, 0.0572 fm, using Wilson loops of spatial size d = 0.68 fm. The influence of the background field on the color flux tube is compatible with the findings of the previous section: strong anisotropies are induced depending on the magnitude and orientation of the external field. In detail, in the L case the flux tube monotonically decreases as the magnetic field grows; in transverse cases (TT-TL), the chromo-electric field is enhanced at eB = 4 GeV 2 while a non-trivial dependence on the lattice spacing is exhibited for eB = 9 GeV 2 : the flux tube is suppressed by the magnetic field at a = 0.0858 fm and compatible with the eB = 0 case at a = 0.0572 fm. We stress that this trend is consistent with the scaling dependence on a observed for the string tension extracted in transverse cases at eB = 9 GeV 2 (see Fig. 7).
In Fig. 12 we show the ratio E l (B, x t )/E l (0, x t ) for each magnetic field value and geometry class at a = 0.0858 fm. Results clearly point out the loss of the cylindrical symmetry in transverse configurations, since the TT an TL cases are not equivalent. Furthermore, both at eB = 4 and 9 GeV 2 the ratio E l (B, x t )/E l (0, x t ) is a decreasing function of x t in the L case, meaning that the color flux tube gets squeezed by the background field. This behaviour had already been outlined at weaker fields in Ref. [42], whose results are displayed together with our findings in Fig. 13. The comparison points out two effects: • the squeezing phenomenon seems to reach a saturation for large fields, since the decreasing dependence on x t is quite similar for eB ≥ 3.12 GeV 2 ; • flux tubes are monotonically suppressed by the increasing background magnetic field. On the other hand, a qualitative weakening of the dependece on eB can be noticed, in agreement with the behaviour of the string tension outlined in Fig. 8 and discussions thereafter.
Despite the observed deformations, the functional dependence of the flux tube profile does not seem to change significantly. We employ a parametrization inspired by the form of magnetic fields inside vortices in type II superconductors to fit the data for the longitudinal chromoelectric field. In particular, we follow the parametrization proposed in [112], where an expression for the magnetic flux tube which solves the Ginzburg-Landau equations is obtained by a variational model for the normalized order parameter of an isolated vortex. This expression, often called Clem ansatz, reads where K n are the modified Bessel functions of the second kind of order n while α, µ and φ are fit parameters. Previous studies showed that the flux tube profile is well described by the Clem function, also in the presence of an external field [42]. Remarkably, we find that the ex-  pression in Eq. (18) is a suitable model even at the large magnetic fields explored in this work. As an example, in Fig. 14 we show the chromo-electric fields obtained at the finest lattice spacing in L configurations together with the best fit functions. We checked that the fit works reasonably well for all the choices of magnetic field, geometry class and lattice spacing.
Flux tubes can be used to compute a more significant parameter: the linear energy density ǫ(B). Since the transverse components of the chromo-electric field are negligible, the energy density just reads The integration is performed over the section orthogonal to the quark-antiquark pair axis, requiring to explicitly know the angular dependence of the color flux tubes. However, the chromo-electric field possess a cylindrical symmetry over the plane when B is directed along the charges (L configurations), so that it is independent of the azimuthal angle. In this case, the integration procedure can be pursued based on data extracted along just one direction on the plane. Furthermore, a numerical integration is not needed: assuming the expression in Eq. (18), then the known integrals of the modified Bessel functions can be exploited (see, e.g., Eq. (5.52.1) and Eq. (5.54.2) in Ref. [113]), leading to so that the linear energy density can be expressed in terms of best fit parameters. This allows to avoid problems regarding the numerical integration instability and the systematic uncertainties which would arise due to the sharp peaks of the flux tube profile. So, we show in Fig. 15 the ratio ǫ(B)/ǫ(0) extracted from L configurations for both the values of the background field and each lattice spacing, together with the continuum extrapolations performed assuming O(a 2 ) corrections. In a classical picture, the energy density per unit length is strictly related to the string tension, since the latter is nothing but the slope of the linear term in the potential. A direct comparison is not possible due to the strong dependence of ǫ(B) on the smearing procedure. However, this issue is overcome by taking into account the ratio ǫ(B)/ǫ(0), where the dependence on N AP E is expected to disappear. Actually, the fact that the ratio E l (B)/E l (0) is independent of the smearing procedure, as seen in Fig. 9, does not imply a priori that ratios of fit parameters (and so ǫ(B)) are independent too. Nevertheless, this independence is numerically observed in the energy density for all the values of eB. The comparison is hence possible and it is performed in Fig. 16, where the continuum limit of the ratio ǫ(B)/ǫ(0) is shown together with the ratio of the string tension σ(B)/σ(0) extracted in the longitudinal case: results are in nice agreement, whithin errors.

IV. CONCLUSIONS
The motivation for the present work lies in Refs. [40][41][42], and in particular Ref. [41], where a prediction was made for a possible vanishing of the string tension for quark-antiquark separations in the direction longitudinal to a magnetic background field and for field values eB 4 GeV 2 , based however on the extrapolation of results obtained from simulations at smaller field values. Which kind of new QCD phase could emerge, if any, where the string tension vanishes in just one direction, is an intriguing question which deserves an answer.
In order to make progress in this direction, in this study we have pushed the range of magnetic background fields explorable by lattice simulations, with a control over the continuum extrapolation, by considering a set of three different lattice spacings, going down to a ≃ 0.057 fm, and a discretization of N f = 2 + 1 QCD similar to that of Refs. [40][41][42], i.e. based on stout improved rooted staggered fermions. In this way, we have been able to reach eB ≃ 9 GeV 2 .
The main result is that, contrary to the expectations of Ref. [41], the string tension in the longitudinal direction is clearly non-vanishing for eB ≃ 4 GeV 2 and still at two standard deviations from zero even at eB ≃ 9 GeV 2 , where however it is suppressed by one order of magnitude with respect to its value at zero magnetic background. On the other hand, the enhancement of the string tension, as a function of eB, in the transverse direction seems to reach a saturation at around 50 % of the string tension value at B = 0. The analysis of the color flux tube shows a consistent suppression/enhancement of its overall amplitude, with mild modifications of its profile, consistent with those already observed in Ref. [42]. In particular, one observes a mild squeezing of the flux tube of quark-antiquark separations parallel to the magnetic field, and a loss of cylindrical symmetry for transverse separations. Notwithstanding such deformations, the flux tube profile is still describable by models inspired to dual superconductivity of the QCD vacuum in all the explored cases.
Finally, the analysis of the chiral condensate shows a persistence of magnetic catalysis in the whole range of explored fields, with a behavior compatible with a lowest Landau level approximation, in particular with a linear dependence of the chiral condensate on B which is in agreement, within errors, with that already observed for eB 1 GeV 2 in Ref. [22].
To summarize, present results postpone to even larger magnetic fields the possibile emergence of a new phase of strong interactions, characterized by the vanishing of the string tension for quark-antiquark separation parallel to the magnetic field, and by other possible associated new phenomena which have not been observed so far. The critical field could be not far from where we are now, since the longitudinal string tension is at just two standard deviations from zero at the largest explored field, however a careful investigation will require simulations on finer lattices: in the future we plan to put further efforts along this direction. A different direction is to investigate QCD at finite temperature for the same lattice spacings and magnetic background fields explored in the present study, since that could give indications about the phase structure from a different perspective: work is in progress along this line [114].