Thermal QCD in a non-uniform magnetic background

Off-central heavy-ion collisions are known to feature magnetic fields with magnitudes and characteristic gradients corresponding to the scale of the strong interactions. In this work, we employ equilibrium lattice simulations of the underlying theory, QCD, involving similar inhomogeneous magnetic field profiles to achieve a better understanding of this system. We simulate three flavors of dynamical staggered quarks with physical masses at a range of magnetic fields and temperatures, and extrapolate the results to the continuum limit. Analyzing the impact of the field on the quark condensate and the Polyakov loop, we find non-trivial spatial features that render the QCD medium qualitatively different as in the homogeneous setup, especially at temperatures around the transition. In addition, we construct leading-order chiral perturbation theory for the inhomogeneous background and compare its prediction to our lattice results at low temperature. Our findings will be useful to benchmark effective theories and low-energy models of QCD for a better description of peripheral heavy-ion collisions.


Introduction
Strong magnetic fields are known to appear in several physical systems around us in the Universe.Magnetars, for instance, are strongly magnetized neutron stars with magnetic fields up to 10 15 G ( √ eB ∼ 1 MeV) in their interior [1].During mergers of such compact objects considerably high temperatures are also achieved [2].Experiments of non-central heavy-ion collisions can produce transient fields of strength 10 18 -10 19 G ( √ eB ∼ 0.1 -0.5 GeV) at RHIC and LHC energies, respectively [3].Furthermore, cosmological models predict even stronger fields emerging during the electroweak phase of the early Universe, with magnitudes as high as 10 20 G ( √ eB ∼ 1.5 GeV) [4].In turn, strong magnetic fields might also be viewed as valuable probes to study the intricate structure of the non-perturbative QCD medium.Since the strength of the above fields is comparable to the QCD energy scale, the interplay of strong and electromagnetic forces must be taken into account to reach a better understanding of how quarks and gluons behave in the systems mentioned above.
In this work, we will focus on the impact of magnetic fields in the context of heavy-ion collisions.The case of strong uniform fields has been widely studied both numerically on the lattice and analytically via QCD models.For a recent review on these approaches, see e.g.Ref. [5].However, in heavy-ion collision experiments, the generated fields deviate significantly from the uniform case.In addition, since the magnetic field also changes rapidly with time (see, e.g. the review [6]), it also induces non-uniform time-dependent electric fields.These might play a role in the subsequent dynamics of the byproducts of the collision and the behavior of the quark-gluon plasma.In particular, event-by-event simulations of heavy-ion collisions predicted nontrivial profiles for both the electric and magnetic field components [7,8], with potential impact for particle flow [9].
The explicit time dependence and the presence of (real) electric fields cannot be implemented using standard lattice QCD simulations in equilibrium due to the infamous complex action problem. 1owever, well established simulation algorithms can be employed to include inhomogeneous magnetic fields.This way we can take one step towards a better description of the heavy-ion collision scenario that helps capture the correct physics in peripheral events.In this paper we consider, for the first time in lattice QCD simulations, a non-uniform background magnetic field B(x) modulated in one spatial direction.Our choice cosh −2 (x/ϵ) for the profile is motivated by the inhomogeneities observed in the above mentioned model simulations of heavy-ion collisions, as well as due to the possibility of an analytical treatment of the free Dirac operator in this case [12,13].We note that our approach can be generalized to introduce modulations in the other spatial directions, too.Our preliminary results have already been discussed in [14].
The magnetic field is known to lower the QCD transition temperature, thereby preferring deconfinement and chiral symmetry restoration [15].The presence of an inhomogeneous magnetic field therefore raises the intriguing possibility of an inhomogeneous medium, where regions exposed to strong and weak B resemble different phases of the theory.If the magnetic field is so strong that the homogeneous setup exhibits a first-order phase transition (as recently observed in Ref. [16]), the different regions of the inhomogeneous system might even become separated by sharp boundaries, leading to novel physical consequences.Furthermore, the characteristic features for deconfinement and for chiral symmetry restoration may be affected differently by the magnetic field, making the inhomogeneous field an effective novel probe of the QCD medium.
This paper is organized as follows: in Sec. 2 we discuss the basic aspects of magnetic fields on the lattice, reviewing the flux quantization for homogeneous fields, generalizing it to the inhomogeneous case and discussing valence and sea effects.Sec. 3 summarizes the details of our simulations, followed by Sec. 4, where we present our results for the local observables.Finally, we summarize our findings in Sec. 5.In the appendices we give details on the implementation of the inhomogeneous magnetic field and the eigenvalues of free fermions and bosons on such backgrounds.The latter is relevant for the construction of chiral perturbation theory for this setup, which we perform here for the first time.
2 Magnetic fields in QCD and on the lattice

Homogeneous magnetic fields
In this section we will review how homogeneous magnetic fields are implemented on the lattice, how key observables are defined and discuss their behaviour to set the stage for the inhomogeneous case.Throughout we will consider lattice geometries with N µ points in the µ direction and a lattice spacing a.The spatial sizes, the volume and the temperature are given by L i = N i a, V = L x L y L z and T = (N t a) −1 .The coordinates run over a symmetric range, i.e. −L x /2 ≤ x < L x /2.Our boundary conditions are periodic in space and (anti-)periodic for bosonic (fermionic) fields in imaginary time.
In order to implement the magnetic field on the lattice, besides the non-Abelian SU(3) links corresponding to the gluon fields in QCD, we must also introduce Abelian links u µ = e iaqAµ ∈ U(1) representing the electromagnetic field A µ for a quark with electric charge q.For a homogeneous field B pointing in the z direction, realized by the gauge A y = Bx, a simple choice of links would be u y = e iaqBx and u x = u z = u t = 1.However, this prescription is not periodic: For lattice simulations it is convenient to have the U(1) links satisfy periodic boundary conditions.Therefore, we perform a gauge transformation and derive an equivalent prescription for the links in the uniform case The periodicity of the lattice imposes that the magnetic flux has to be quantized according to [17] The QCD partition function in the presence of a uniform background magnetic field using the rooted staggered formulation of fermions is where the magnetic field dependence enters through the Dirac operator / D of each flavor.The gluon action S g will be specified below in Sec. 3. The Abelian links entering the Dirac operator of a certain flavor depend on the combination q f B, see (2.1), which we indicated explicitly.We will concentrate on two observables, the quark condensate and the Polyakov loop.The former is defined via the derivative of log Z(B) with respect to m f , with Tr standing for a trace in both Dirac and color spaces as well as a summation over the full lattice.In turn, the Polyakov loop is defined as x,y,z Re tr and the tr stands for a trace only in color space.
Both of these observables need renormalization.The quark condensate, as defined in Eq. (2.4), contains an additive and a multiplicative ultraviolet divergence, which are, however, independent of the temperature and the magnetic field.Following Ref. [18], we renormalize it by defining (2.6) Here, the subtraction of the zero-field and zero-temperature condensate eliminates the additive divergence, whereas the multiplication by m f cancels the multiplicative one.The normalization factor involving the physical pion mass m π = 135 MeV and the chiral limit of the pion decay constant F = 86 MeV is for convenience.The combination (2.6) is unity in the vacuum and, according to the Gell-Mann-Oakes-Renner relation, it approaches zero for high temperatures.We define the condensate for the light quarks as Σ ≡ (Σ u + Σ d )/2, which we use as the approximate order parameter in the chiral sector.The Polyakov loop also contains an ultraviolet divergence, which can be eliminated via [19] where the multiplicative factor W (a, T ) is independent of the magnetic field and has been measured for our action in Ref. [20].

Inhomogeneous magnetic fields
In order to implement inhomogeneous fields on the lattice we need to extend the procedure described in the beginning of Sec.2.1.The details of the generalization are relegated to App.A, here we only summarize the main results concerning our choice of inhomogeneous field profile.We consider a magnetic field pointing in the z direction with the particular form where ϵ is the width of the profile, centered in the middle of the lattice and B is the magnitude of the field at the center.Fig. 1 illustrates the shape of the profile for different values of ϵ.Similarly to the uniform case, the flux of the inhomogeneous field is also quantized, allowing us to parameterize the field with the integer N b and the continuous variable ϵ.As Fig. 9 shows, the homogeneous field is recovered by keeping N b constant and approaching ϵ → ∞.The prescription for the links in this case is (see App. A), It is interesting to note that only the U(1) links in the bulk of the volume are sensitive to the particular spatial dependence of the magnetic field -the u x links originating from the gauge transformation required for periodicity only depend on the flux quantum N b .Consequently, the twisted u x links for the inhomogeneous (2.10) and homogeneous (2.1) profiles are identical.Notice moreover that even though B(x) is periodic, its derivatives differ at x = −L x /2 and x = L x /2.This is, however unproblematic as long as ϵ/L x is sufficiently small, i.e. the volume is large compared to the profile width.We point out here, that from now on, whenever we refer to the value of the field, we mean the amplitude B, unless otherwise explicitly noted.We also extend our list of observables to probe the effects of the inhomogeneous B-field.Specifically, we define local (at least in one spatial dimension) versions of Eqs.(2.4) and (2.5).Then, the local quark condensate for a flavor f at finite B is given by ) where Tr stands for a trace in Dirac and color space and a summation over the y, z and t coordinates, retaining the x-dependence of the local quark condensate.With a similar mindset, the local Polyakov loop is defined as Re tr t U t (x, y, z, t) .
(2.12) The renormalization of the local objects is carried out by the same procedure as for the fully volume averaged quantities of Sec.2.1, as the ultraviolet divergences do not depend on the coordinate.Therefore the renormalized local condensate, generalizing (2.6), is defined as and the average light quark condensate is denoted by Σ(x) ≡ (Σ u (x) + Σ d (x))/2.We can also define the modification of the local chiral condensate relatively to the zero-field one as which will be useful for our analysis later on.Finally, the renormalized local Polyakov loop is (2.15)

Valence and sea effects
As we will see, the magnetic field has a nontrivial impact on the observables that can be best understood in terms of a separation of valence and sea effects.Let us first summarize what we know about these in the uniform magnetic field case.For the condensate, the magnetic field dependence appears in two ways [21].First, in the trace in Eq. (2.4), which represents the direct coupling between the electrically charged quarks and B and is called the valence effect.This coupling tends to enhance the condensate and this enhancement is often referred to as magnetic catalysis.Magnetic catalysis can be understood in various different ways: for strong magnetic fields it emerges due to the dimensional reduction of the system [22], while at low B it is related to the positivity of the QED β-function [23,24].Either way, magnetic catalysis is deeply connected to the proliferation of low Dirac eigenvalues, which has been demonstrated explicitly for full QCD [25].
Second, the magnetic field also has an indirect effect on the distribution of gluonic configurations via sea quark loops, i.e. via the determinant in Eq. (2.4).This is the so-called sea effect.Its main impact for the gluon fields is to enhance the Polyakov loop in the transition region [20], which, in turn, correlates with the suppression of the condensate.The sea effect is found to dominate the valence effect in the transition region, leading to an overall reduction of ⟨ ψψ⟩ for the light quarks due to B, dubbed inverse magnetic catalysis.In turn, away from the transition temperature, magnetic catalysis prevails.Altogether, the effect of the magnetic field is to reduce the transition temperature, defined via the inflection point of the condensate.The picture that we just sketched holds for physical quark masses and the inverse magnetic catalysis phenomenon disappears for heavier quarks.However, the reduction of T c (B) was found to persist for all quark masses [26,27].
In our inhomogenous setup, the condensate (2.11) also has a twofold dependence on B, suggesting an approximate factorization into valence and sea sectors again.Specifically, we define the bare valence and sea condensates by ) where the magnetic field was switched off in the determinant or in the trace, respectively.Using these observables, we define the relative change of the chiral condensate with respect to the zero-field one due to the valence and the sea contributions We can also define the light-quark average of these observables as r val = (r val u + r val d )/2 and r sea = (r sea u + r sea d )/2.For weak fields, the full condensate can be shown to be given by the sum of the sea and valence terms [21].For larger values of B, this factorization is not exact anymore, but still insightful.
As we will see below, the impact of the two effects is quite different on the condensate: the valence effect acts in a predominantly local manner, while the sea effect is more global.Still, r sea (x) can depend on x if the gluon configurations that are preferred for B(x) are themselves inhomogeneous.Exactly this is measured by the local Polyakov loop, which only depends on B via the sea effect.

Simulation setup
We used N f = 2 + 1 flavors of rooted, stout smeared staggered quarks with physical masses and the tree-level improved Symanzik gauge action [28].The electric charges of the quarks are set to q u = 2e/3, q d = q s = −e/3, where e > 0 is the elementary electric charge and the quark masses in lattice units are tuned along the line of constant physics [28].
We generated gauge configurations with this action on four lattice sizes: 16 3 × 6, 24 3 × 8, 28 3 × 10 and 36 3 × 12 to approach the continuum limit.We varied the coupling β to cover temperatures from 113 MeV to 176 MeV, in order to explore the behavior of QCD matter around the crossover temperature T c ≈ 155 MeV.Similarly, we vary the flux quantum N b , assuming a background field given by Eq. (2.8), to cover values of √ eB from 0 GeV to 1.2 GeV, which is around the phenomenologically relevant strengths.We fixed the width of the profile in physical units to ϵ ≈ 0.6 fm, which produced an appreciable inhomogeneity in the field on the lattice and is consistent with the characteristic width of the field that was found in heavy-ion collision simulations [8].
We remark that our lattices correspond to a fixed physical volume V = L x L y L z .In contrast to standard thermodynamics studies, where the infinite volume limit of homogeneous systems is of interest, here we do not wish to approach that limit but are interested in the local spatial dependence of observables.
For each set of parameters, we measured the local quark condensate and the local Polyakov loop.To calculate the condensate, we need to compute the inverse of the massive Dirac operator M = / D(U, q f B) + m f .We use the standard technique of noisy estimators, i.e.N complex vectors whose components χ k (x, y, z, t, c) are randomly sampled from a normal distribution, such that where the correction terms O(1/ √ N ) vanish in the limit N → ∞.Here, x, y, z and t are coordinates and c a color index.The projected trace, appearing in (2.11), can then be approximated using these random vectors as We found that N = 40 provides a reasonably good estimate for the local condensates.Finally, the lattice data are extrapolated to the continuum limit via a combined fit that is explained in detail in App.D.

Results
Due to the inhomogeneity of B(x), different parts of the system are influenced by magnetic catalysis and inverse catalysis differently, depending on the local value of the magnetic field.As a consequence, non-trivial features appear in the observables -this is what we explore next.

Local magnetic catalysis
We first discuss the B-and T -dependence of the light quark condensate in our inhomogeneous setup at low temperature.In this region, the system is dominated by pions and its response to B(x) can be described by chiral perturbation theory (χPT).We work out the details of the leading-order theory in App.B and App.C. In Fig. 2, we plot the modification of the chiral condensate defined in Eq. (2.14), together with the χPT prediction at T = 113 MeV.Notice the local enhancement of Σ(x) due to B(x) -a behavior that we dub local magnetic catalysis.This local increase in the condensate is described perfectly by χPT at low B (left plot), while the full QCD results start to deviate from its prediction as the magnetic field increases.These deviations are more pronounced at x = 0, where B is the strongest.Altogether, we thus find that the inhomogeneous variant of χPT gives a good description of the lowtemperature and low-magnetic field region, similarly to the conclusions drawn for the homogeneous case [18].
To quantify how much the effect on the condensate follows the local profile of the magnetic field, we define an auxiliary observable, constructed solely from the data for the chiral condensate obtained for homogeneous magnetic fields [18], as  where B(x) is our 1/ cosh 2 (x) profile.We also implicitly define the light-quark average of this observable.The so defined quantity -which we name the point-wise chiral condensate -mimics a completely local response, where the observable at point x only depends on the local strength of B at that point.In Fig. 3, we compare the point-wise observable with the full condensate (2.13) at different temperatures.At low T (left plot), the difference between the peaks of the condensates is only about 6%, which means that the full condensate feels the magnetic field approximately in a local manner.However, in the transition region (right plot), the point-wise condensate would lead to the wrong behavior, as it decreases abruptly at x = 0 compared to the full one.The difference in the peaks is ∼ 60% and it indicates that the global effects of B are much more manifest around the crossover.

Local inverse magnetic catalysis
Next, we discuss our full results for Σ(x) in more detail.we observe an overall decrease of the condensate for increasing B. Second, Σ(x) is also deformed and we observe the formation of dips in the tails of the curves, i.e., localized regions where Σ(x) is smaller even compared to the edges of the lattice, where the magnetic field practically vanishes.This nontrivial behavior is due to the interplay of the sea effect (reducing the condensate globally) and the valence effect (increasing it locally) that we will explain in more detail below.
The appearance of the dips, and the nontrivial behavior of Σ(x) is a new type of indication of the inverse magnetic catalysis phenomenon.
There are two avenues to understand this behavior better: via a separate analysis of the valence and sea contributions to Σ(x), and via an investigation of the Polyakov loop P R (x).First we discuss the contributions from the sea and the valence effects, as defined in Eq. (2.18).These are plotted in Fig. 5 and compared to the full condensate.At low temperature (left plot), the dominant contribution is clearly given by the valence part, whereas the sea contribution is only slightly above zero.On the other hand, at T ≈ T c the situation changes and the sea contribution becomes not only relevant but also negative, therefore acting to decrease the condensate.For even higher B, the negative contribution of the sea effect tends to dominate over the valence effect, although in this regime the additive property of the two effects no longer holds.
Notice that the sea contribution in both cases tends to be more broadly distributed, whereas the valence one is more sharply peaked.That difference is essential for the appearance of the dips in the chiral condensate and can be explained via the Polyakov loop, which we discuss in the next section.

Local Polyakov loop
The sea effect is known to be rooted in the behavior of the Polyakov loop; therefore we turn to this observable next.Its continuum extrapolation is also discussed in App.D. In Fig. 6, we show the continuum extrapolated, renormalized local Polyakov loop as a function of x for two different magnetic fields at the temperature (T = 155 MeV), where we found the dips to develop in the light quark condensate.We first notice that the Polyakov loop is also enhanced where the magnetic field is stronger.As we discussed above, the Polyakov loop is a purely gluonic quantity that does not couple to the magnetic field directly, and is only affected by the sea effect through the fermion determinant.Since the latter depends on the magnetic field at all lattice points, this implies that B has a smeared impact on P .Therefore, the Polyakov loop feels the magnetic field in a broader region of the lattice compared to the quark condensate.This is manifested by the broader peaks in P R (x) in the figure.What is the impact of such a Polyakov loop profile on the condensate?In the homogeneous case, it is known that P and ψψ anticorrelate with each other in the transition region [20].To answer this question for our inhomogeneous setup, we need to quantify the extent of local correlations between the two observables.To this end, we define the correlator of the bare quark condensate and the bare Polyakov loop where we employed a convenient normalization to make C dimensionless. 2 In Fig. 7, we show a heat plot of C(x, x ′ ) for different temperatures at approximately the same magnetic field.We can see that there are no sizeable correlations at low or high temperatures.In contrast, there is a strong local anticorrelation between the Polyakov loop and the condensate in the transition region.This anticorrelation is most prominent along the diagonal of the figure, i.e. when both quantities are taken at the same coordinate, but the negative correlation persists even away from the diagonal.Thus, the suppressing effect of the Polyakov loop can impact the condensate from a distance of about 0.5 fm.
We can now put together the valence and sea effects for the condensate.The valence effect enhances Σ val (x) in roughly the same shape as that of the magnetic field B(x).The sea effect acts via the Polyakov loop, which feels the magnetic field in a smeared (global) way and is enhanced in a broader region.The anticorrelation with the local condensate thus suppresses Σ sea (x) in a broader interval, too.The total effect arises as a combination of both, explaining the formation of the dips in Σ(x) in the transition region.This combined effect is illustrated again in Fig. 8, where we show a comparison between the renormalized quantities, namely the Polyakov loop and the chiral condensate, normalized by their sum over the whole lattice for T < T c and T ∼ T c .To highlight the local variation of the observables, here we subtract their values at the edge of the volume x = ±L x /2.

Summary and conclusions
In this paper we presented the first lattice simulations of QCD with physical quark masses in the presence of inhomogeneous external fields.In particular, we considered a localized background magnetic field B(x) ∝ cosh −2 (x/ϵ).Choosing a profile width ϵ = 0.6 fm, this setup resembles the conditions relevant for the early stages of off-central heavy-ion collisions.We calculated observables relevant for the thermodynamics of the system: the local light quark condensates as well as the local Polyakov loop.Our results are provided in ancillary files submitted to arxiv.orgalongside with this publication and will serve as useful benchmarks for effective theories and low-energy models of QCD, in a similar way to the results for homogeneous fields [15,18].This is in particular the case because the profile of the field allows for an analytical treatment in NJL-type models, as in Ref. [13].Our findings may also provide useful input for hydrodynamic models aiming to describe off-central heavy-ion collisions.
At low temperatures, we established the local analogon of the well-known magnetic catalysis phenomenon: the enhancement of the local condensate due to the local magnetic field.We worked out how to treat this setup within leading-order chiral perturbation theory (for details, see App.C).The corresponding prediction for the condensate was found to exhibit good agreement with the lattice results for weak magnetic fields at low temperature.Moreover, we found that this effect mainly stems from the valence sector and acts in a local manner on the condensate.
In contrast, the behavior of the condensate is much more complicated in the transition region.Here, we performed a systematic study of the dependence of the local condensate and the local Polyakov loop on the magnetic field, the temperature as well as the coordinate x.Our results show an enhancement of the condensate Σ(x) similar to the shape of B(x), together with a reduction in a broader interval, leading to the development of the 'dips' away from the center.This is a novel, nontrivial manifestation of the inverse magnetic catalysis phenomenon, arising from a delicate interplay between the local Polyakov loop and the local condensate via their local anticorrelations.
In summary, our results reveal that the inhomogeneous background magnetic field can realize situations where different parts of the system feature substantially different values for the relevant thermodynamic observables.In terms of the local Polyakov loop for example, one might describe the center of the system to behave more as deconfined, whereas the outer regions more as confined matter.However, unlike in the case of homogeneous fields, this picture is more complicated due to the quark condensate being also more enhanced in the center and reducing towards the edges.In other words, the extent to which deconfinement and chiral symmetry restoration occur simultaneously, is different point by point.Due to the crossover nature of the transition, a smooth dependence of all observables on x is natural.An intriguing question is whether the local variations in these observables become abrupt when magnetic fields are investigated for which the QCD transition becomes first order [16], even though these are not relevant for heavy-ion phenomenology anymore.

B Classical Dirac equation with an inhomogeneous magnetic field
Our choice (2.8) is special in the sense that (in the absence of color interactions) the Dirac equation can be solved analytically for this profile [29], just like for the standard Landau problem with homogeneous magnetic fields.The lowest Landau levels are the zero modes of the two-dimensional Dirac operator (in the x − y plane).Protected by the index theorem, these remain zero modes even in the inhomogeneous case.In contrast, higher Landau levels split up into N b branches, where 2πN b is the flux of the magnetic field in a finite volume.However, besides bound states, the spectrum also contains a continuum of scattering states that are difficult to treat analytically.
We have calculated the Dirac eigenvalues λ n and eigenvectors χ n numerically on a 16 2 lattice using the staggered discretization [30].The lattice eigenvalue spectrum λ n (N b ), plotted in Fig. 10, constitutes a non-trivial generalization of Hofstadter's butterfly [31,32], which is recovered in the ϵ → ∞ limit.Most of the recursive pattern of the homogeneous butterfly is now lost, since the magnetic field breaks translational symmetry and the spectrum is not a periodic function of the flux quantum N b anymore. 3The characteristic gaps in the spectrum are also filled with eigenvalues, except for the largest gap separating lowest Landau levels from other eigenstates -this remains a pronounced feature of the inhomogeneous butterfly due to the topological nature of the zero modes.It is interesting to point out that the same holds also when the homogeneous system is perturbed by color interactions [25].
Using the so calculated eigenvalues, together with the chiral symmetry of the Dirac operator (for the staggered case, reconstructed as ψψ where p z a = sin(2πn z /N z ) and p t a = sin(2π(n t + 1/2)/N t ) with integers 0 ≤ n z < N z and 0 ≤ n t < N t are the lattice momenta and fermionic Matsubara frequencies for the plane waves in the z and t direction, respectively, and the prefactor 1/4 emerges due to rooting.Owing to the normalization of the modes, x |χ n (x)| 2 = 1, this reproduces the usual expression for the average quark condensate x ψψ(x)/L.Fig. 11 shows the condensate for a quark mass of ma = 0.025 on 16 3 × 4 lattices.(Note that m π = 2m simply in this non-interacting case but we use the physical value for m π /f π .)The left panel Figure 13: Numerical results for the average light quark condensate obtained using the free pion propagator for homogeneous magnetic fields on N 4 s lattices, together with a comparison to standard, continuum chiral perturbation theory [35].
where the lattice momenta and Matsubara frequencies for scalar particles, p z a = 2 sin(πn z /N z ) and p t a = 2 sin(πn t /N t ) with 0 ≤ n z < N z and 0 ≤ n t < N t enter [33].Using the Gell-Mann-Oakes-Renner relation (m u + m d ) ψψ = m 2 π F 2 , we can trade the derivative with respect to m 2 π for a derivative with respect to the quark masses.For degenerate quark masses, m u = m d = m ud , we can express the average light quark condensate in our normalization (2.6) as The results are checked by comparing against the well-known chiral perturbation theory formula in the homogeneous magnetic field case at zero temperature [34,35].This comparison is shown in Fig. 13 for a set of T ≈ 0 lattices using m π L = 6.The continuum limit is approached by keeping m π L fixed and sending N s = N t → ∞, revealing perfect agreement with the analytical formula.
Having cross-checked our approach for homogeneous fields, we proceed to calculate ∆Σ(x) for inhomogeneous profiles.Fig. 14 shows a scan on 16 4 lattices for a range of widths and magnetic fluxes.Our results demonstrate clearly what we call 'local magnetic catalysis' -the inhomogeneous enhancement of the condensate.Just as in the fermionic case, we again find the spatial average of the condensate to be practically independent of ϵ.We compare the so obtained profiles to our full QCD results in the main body of the text.

D Lattice data and continuum extrapolation
In order to compute the continuum limit of our observables in the full range of magnetic fields, we need to interpolate over the discrete N b values, given by the quantization of the magnetic flux in a finite box.In this appendix, we discuss in detail our method of interpolation for the chiral condensate and the Polyakov loop.In order to have a reliable extrapolation to the continuum and to suppress local fluctuations of the data due to statistical noise, we applied an interpolation procedure that combines a multidimensional spline fit -defined upon a set of spline nodepoints -and the continuum extrapolation.This is done by promoting the coefficients of the spline to be a-dependent and finding the best fitting surface that minimizes the 'action' -containing the χ 2 /dof as well as a term suppressing oscillatory solutions.Middle plots: chiral condensate at x = 0 as a function of eB obtained on different lattices (data points) and the continuum limit extracted from the spline fit (green band).Bottom plots: a-scaling of the spline fits at x = 0 for different eB with continuum limit projection (green band).The error bands include statistical and systematic errors added in quadrature.
To assess systematic errors, a Monte-Carlo procedure is set up to generate fits having different numbers and positions of spline nodepoints, according to this action.The details of the spline fit and the Monte-Carlo algorithm are explained in Refs.[36] and [37,38], respectively.
In particular, at each T , we considered the two-dimensional space spanned by x and eB including a lattice spacing dependence, and performed a series of Monte Carlo updates for the corresponding spline fits.Owing to the scaling properties of our lattice action, the lattice spacing-dependence of the spline was fixed to be quadratic and we included the leading O(a 2 ) lattice artefact term.To constrain our spline fits, we exploited parity symmetry, i.e. that the observables are symmetric in x around x = 0 and in B around B = 0.Thus, we imposed the conditions ∂Σ(x, B; a) ∂x and similarly for the Polyakov loop.We chose the initial position of the spline nodepoints individually for each T , given that the observables, for instance the quark condensate, have various degrees of complexity with x for different T .The systematic errors are accounted for by averaging over different possible spline solutions, generated from the Monte-Carlo algorithm.Note that the B = 0 slice of the splines is not constrained to give an x-independent constant -they are observed to fluctuate around the B = 0 value of the observables.In Fig. 15, we show the results of our spline fit for the chiral condensate at a selected set of parameters.To avoid enhanced lattice discretization effects, in the case of the chiral condensate we excluded our 16 3 × 6 data at T ≳ 140 MeV from the global fit.

Figure 1 :
Figure 1: The profile (2.8) of the magnetic field in lattice units for different values of the parameter ϵ on an N x = 16 lattice.For ϵ → ∞ the homogeneous profile is recovered.

Figure 3 :
Figure 3: Point-wise condensate compared to the full one for low T (left plot) and T ≈ T c (right plot) on a 16 3 × 6 lattice.The peaks of the point-wise and the full condensates differ by ∼ 6% at T = 113 MeV and by ∼ 60% at T = 155 MeV.

Figure 4 :
Figure 4: Continuum extrapolated results for the renormalized quark condensate as a function of x for different magnetic fields and temperatures.

Figure 5 :
Figure 5: Sea and valence contributions to the chiral condensate compared to the full one for T = 113 MeV and T = 155 MeV.

Figure 7 :
Figure 7: Heat plot of the correlation C(x, x ′ ) between the bare quark condensate (at point x) and the bare Polyakov loop (at point x ′ ) for temperatures T < T c (upper left panel), T ≈ T c (upper right and lower left panels) and T > T c (lower right panel) on a 24 3 × 8 lattice.The rough size ϵ of the magnetic field is indicated by the gray circle.The scale of the plot is kept fixed in physical units, implying that the heat map area slightly shrinks as a reduces towards higher temperatures.

Figure 8 :
Figure 8: Top panels: Local profiles of the quark condensate (red) and the Polyakov loop (blue) on a 28 3 × 10 lattice at T = 124 MeV (left) and at T = 155 MeV (right).For the sake of visualization, we subtracted the values of both observables at the edge, x = ±L x /2 (denoted by ∆) and normalized by the area.Notice the broader peak of the Polyakov loop and the formation of the dips in the condensate in the transition region.Lower panels: renormalized chiral condensate at the point x minus the one at the edge (±L x /2) showing the formation of the dips for increasing magnetic fields on the N t = 8 and N t = 12 lattices.

Figure 9 :
Figure 9: A gauge transformation at the last x-slice makes the u y links (red arrows) exactly periodic, and also affects the u x links (green arrows).

Figure 10 :Figure 11 :
Figure 10: Eigenvalue spectrum of the staggered Dirac operator in the presence of an inhomogeneous magnetic field characterized by a flux N b and a width ϵ.The three panels correspond to the homogeneous case (ϵ → ∞, left), ϵ = 10a (middle) and ϵ = 4a (right).The first N b eigenvalues (the equivalents of lowest Landau levels) are marked by blue, the next 2N b eigenvalues (first Landau levels) by red, and the rest of the spectrum by orange dots.Note that the eigenvalues come in complex conjugate pairs and the plots only show the positive (imaginary) eigenvalues.

Figure 14 :
Figure 14: The coordinate dependence of the subtracted condensate (2.13) for free charged pions for various values of N b and ϵ.The lattice spacing is set by the pion mass, m π a = 0.375.

Figure 15 :
Figure 15: Top plots: chiral condensate as a function of x for different T .The lattice results (data points) are compared to the one-dimensional slice of the spline fit at the corresponding values of a and eB.Middle plots: chiral condensate at x = 0 as a function of eB obtained on different lattices(data points) and the continuum limit extracted from the spline fit (green band).Bottom plots: a-scaling of the spline fits at x = 0 for different eB with continuum limit projection (green band).The error bands include statistical and systematic errors added in quadrature.