Steady electric currents in magnetized QCD and their use for the equation of state

: In this paper we study the emergence of steady electric currents in QCD as a response to a non-uniform magnetic background using lattice simulations with 2 + 1 quark flavors at the physical point, as well as leading-order chiral perturbation theory. Using these currents, we develop a novel method to determine the leading-order coefficient of the equation of state in a magnetic field expansion: the magnetic susceptibility of the QCD medium. We decompose the current expectation value into valence-and sea-quark contributions and demonstrate that the dominant contribution to the electric current is captured by the valence term alone, allowing for a comparably cheap determination of the susceptibility. Our continuum extrapolated lattice results for the equation of state confirm the findings of some of the existing studies in the literature, namely that the QCD medium behaves diamagnetically at low and paramagnetically at high temperatures.


Introduction
Background electromagnetic fields are featured by a series of physical systems ranging from magnetized neutron stars through heavy-ion collisions to the evolution of the early universe [1].In particular, magnetic fields play a central role for the phenomenology of non-central heavy-ion collisions (HIC).While a direct measurement of the magnetic field is elusive in such experiments due to its short lifetime, it can be estimated via event-by-event simulations.These simulations predict strong and highly non-uniform electromagnetic fields for peripheral collisions [2], with field strengths ranging from √ eB ∼ 0.1 GeV (at typical collisions at RHIC) up to √ eB ∼ 0.5 GeV (at the LHC).Such fields are comparable in magnitude to the characteristic scale of the underlying theory, QCD, and therefore might impact the initial stages of the quark-gluon plasma (QGP) created at the experiment and potentially the subsequent hydrodynamic evolution [3].In fact, recent results by the STAR Collaboration at RHIC observed charge-dependent particle flow signals consistent with Faraday and Coulomb effects due to the electromagnetic fields in peripheral collisions [4].
On the theoretical side, the impact of electromagnetic fields on QCD matter can be studied by lattice QCD simulations as well as within low-energy QCD models and effective theories (see [5] for a review on the topic).So far, most of these studies have assumed spatially homogeneous fields.A theoretical improvement to this picture is to consider a non-uniform magnetic field that more closely resembles the situation in the early stages of off-central collisions.For a specific, analytically treatable inhomogeneous profile [6], a study using the NJL model [7] predicted non-trivial structures appearing in various observables relevant for thermodynamics.Recent lattice QCD simulations revealed that such fields not only induce an inhomogeneous breaking of chiral symmetry, but also affect the confinement aspects of the theory in a non-uniform manner [8].In addition, the inhomogeneous background fields are expected to generate local electric currents in equilibrium QCD [7].
In general, conserved currents appear as a manifestation of global symmetries.In non-relativistic quantum mechanical systems in equilibrium, the global flow of conserved currents in the ground state is forbidden by Bloch's theorem [9].However, the theorem does not forbid a local current density, as long as its spatial integral vanishes.In particular, this is the case for an electric current due to an inhomogeneous background magnetic field in accordance with Ampère's law.This naturally leads us to the question: could inhomogeneous magnetic fields induce local currents in an equilibrium state of QCD and if yes, what role do they play for QCD thermodynamics?
In this paper, we perform lattice simulations to investigate the emergence of electric current densities as a response of the strongly interacting medium to a non-uniform magnetic field in an equilibrium setup.Our simulations employ 2+1 flavors of stout-improved rooted staggered quarks with physical masses and a continuum extrapolation based on four lattice spacings.Using the currents measured on the lattice and Ampère's law, we present a novel approach to compute the leading-order coefficient of the equation of state in a B-expansion, i.e. the magnetic susceptibility of the thermal QCD medium.This new determination of the susceptibility may be compared with existing lattice techniques in the literature [10][11][12][13][14][15].Our results show that, in the continuum limit, the magnetic susceptibility changes sign across the QCD crossover, characterized by the pseudo-critical temperature T c , indicating that the response of the QCD medium to an external magnetic field comes in at least two phases: 1. a weak diamagnetic phase at temperatures below T c , and 2. a strong paramagnetic phase at temperatures above T c .This confirms the findings of [14].
In addition to computing the electric current expectation value in an ensemble fully incorporating the inhomogeneous magnetic field background, we will introduce what we call the valence approximation, based on the weak-field expansion of fermionic expectation values into sea and valence contributions (see App. A), introduced for the chiral condensate in Ref. [16].We empirically show that, unlike the chiral condensate, the electric current has a negligible sea contribution across a broad range of temperatures and magnetic fields.As a consequence, the electric current at a given magnetic field strength B is well-approximated by the valence term, which only requires gauge configurations at B = 0.This observation decreases the computational cost by a substantial factor.Thereby, we suggest that our method could be applied to compute the magnetic susceptibility using different fermion discretizations using less computational resources.To the best of our knowledge, the cal-culation of the magnetic susceptibility has only been done in the staggered formulation so far.
Our lattice simulations are corroborated by a determination of the electric current induced by the inhomogeneous magnetic field within a low-energy effective theory of QCD, chiral perturbation theory.We work out the details of this approach in App.B.
This article is organized as follows.In Sec. 2, we discuss our method to obtain the magnetic susceptibility from electric currents using Ampère's law.In Sec. 3 we describe our simulation setup and define the electric current density in the staggered formulation.This is followed by Sec. 4, which is devoted to the renormalization of the current operator as well as that of the susceptibility.In Sec. 5 we present our lattice results from the full-current and the valence approximation and finally, in Sec.6 we conclude.

Magnetic susceptibility and Ampère's law
In classical electrodynamics, the total magnetic field in a magnetized medium is given by where H is the external field that one would apply to the system in an experiment and M is the magnetization, i.e. the response of the system to the applied field.Therefore, B represents the total magnetic field that would be measured by an external observer.In Eq. (2.1), we used natural units for the vacuum magnetic permeability µ 0 = 1.Moreover, we expressed the magnetization in units of the elementary electric charge e > 0 -this will be convenient later when we discuss renormalization.To each field in Eq. (2.1) there are associated electric current densities, which follow a similar relation where j H is the free current density that generates H, j is the magnetization current density, and j B is the total current density.Since we are interested in the response of the QCD medium to the background field, we will focus on j.
In a static medium, the magnetization at weak magnetic fields is related to the total field by where V is the volume, χ the local magnetic susceptibility and we neglected terms of O(B 3 ).In general, χ may be a function of the coordinates, but the translational invariance of the medium constrains this dependence.We made this fact explicit by writing χ(r, r ′ ) = χ(r − r ′ ).The magnetization current can be obtained from M via Ampère's law for the magnetization, where we integrated by parts and relabeled r ′ → r − r ′ , making use of translation invariance again.
Next, we specialize to a magnetic field of the form B(r) = B(x) e z , where x is the first component of the vector r.In this case, Eq. (2.4) implies that the current j points in the y direction and is given by The Fourier transform of Eq. (2.5) gives where p denotes the momentum in the x direction.This result is still obtained within a weak-field expansion, i.e. we neglected O(B 3 ) terms on the right hand side.In this work, we will be mostly interested in the susceptibility at zero momentum.This encodes the response of the medium to homogeneous magnetic fields, and we will denote the corresponding, leading weak-field coefficient as χ m = χ(p = 0).Alternatively to the above derivation, the magnetic susceptibility χ m can be obtained directly from the partition function Z of the system in a homogeneous background magnetic field, oriented along the z direction without loss of generality, where T denotes the temperature.This definition of χ m is consistent with the p = 0 limit of χ(p) considered in Eq. (2.6).The factors of e appear in a manner so that the derivative of log Z in Eq. (2.7) is taken with respect to the renormalization-group-invariant combination eB.

Susceptibility on the torus
The discussion so far corresponds to the infinite volume system.Next we will consider a finite volume with periodic boundary conditions, in order to comply with the standard setup in lattice simulations.Specifically, we take our volume to be a symmetric box with side length L, and label the coordinates so that −L/2 ≤ x ≤ L/2.Furthermore, we will work with an inhomogeneous magnetic field profile, where ϵ is the inhomogeneity parameter.In the discussion below we do not exploit the specific form of B(x), but only the fact that it is an even function, B(−x) = B(x).The latter also implies j y (−x) = −j y (x).
Let us now come back to Eq. (2.6), and focus on the zero-momentum limit of χ(p).In the infinite volume, when p is continuous, this can be achieved via L'Hôpital's theorem, On the torus, however, p is discrete, and we are formally not allowed to take the p → 0 limit.Therefore, we consider χ(p) at finite momentum and rewrite it using the parity symmetry of the current, j y (−x) = −j y (x), and the magnetic field, B(−x) = B(x), where in the last equality we integrated by parts.For low p, we can approximate the sine in the integral above using its Taylor series.For small x, we can expand the sine around 0, whereas for x ∼ L/2, we can expand it around L/2. Therefore, in the p → 0 limit, we approximate the sine in Eq. (2.10) by a wedge-like function in the interval [0, L/2], which has the same parity symmetry as the sine in that interval, and it satisfies the boundary conditions.Hence, the zero-momentum susceptibility becomes (2.12) Notice that the magnetic field (2.8) decays exponentially towards the boundaries of the volume and so does the electric current.Therefore, the contribution from 0 ≤ x ≤ L/4 dominates the integral as L → ∞ and Eq.(2.12) approaches Eq. (2.9) in the infinite-volume limit.In Sec. 5, we will apply Eq. (2.12) to calculate χ(0) on the lattice.Just like Eq. (2.6), all equations above are formulated as a weak-field expansion, with O(B 3 ) terms neglected.To recover the magnetic susceptibility (2.7) from the observables determined at B > 0, we need to extrapolate to B = 0, At B = 0, both j y (x) and B(x) vanish but the ratio in the right-hand side of Eq. (2.13) remains finite.Furthermore, as we will discuss in Sec. 3, in a finite box the magnetic field is quantized and we only have access to discrete values of B. Therefore we extrapolate to vanishing fields using the ansatz c 1 + c 2 (eB) 2 , valid for low B, based on the symmetry of the partition function Z(B) = Z(−B), which only allows for even powers of B in the susceptibility.We obtained the parameters c 1 and c 2 by fitting our finite B lattice data.The so obtained χ m contains an additive logarithmic divergence which has to be cancelled via renormalization.We will discuss its origin and our renormalization scheme for χ m and j y in Sec. 4.

Non-uniform magnetic fields on the lattice
In a finite box, the flux of the magnetic field is quantized.In order to couple a quark flavor f , with charge q f , to a magnetic field background given by Eq. (2.8) with amplitude B in the z-direction, where −L x /2 ≤ x < L x /2, we need the following U(1) factors (see [8] for more details), in addition to a quantization condition for the magnetic flux in a finite box, Here, q is a reference charge.For a many-flavor system, the reference charge is given by the modulus of the greatest common divisor of the charges.In the case of u, d, and s quarks: The parameter ϵ controls the width of the profile.The value of ϵ was chosen such that B(−L x /2) = B(L x /2) ≈ 0 to avoid substantial effects from the boundaries.Rewriting Eq. (3.2) in terms of the aspect ratio N s /N t , where N s and N t are the number spatial and temporal sites, respectively where T = 1/aN t .This equation implies that at fixed T , ϵ, and N b the physical value of the magnetic field B is the same between different ensembles if they have the same aspect ratio.This observation will be useful in Sec.5.1 when we discuss the results.In Fig. 1 we illustrate the magnetic field profile implemented by the links in (3.1) and its first spatial derivative, which will be relevant for the following discussion regarding the electric current.

Full-current setup
We simulated 2 + 1 flavors of rooted staggered fermions with physical masses.The Dirac operator / D is a one-link operator including two steps of stout smearing.In the rooted staggered formalism, the µ-component of the electric current density is given by ) where M(q f B) = / D(U, q f B) + m f is the massive Dirac operator, and Tr x represents the trace over y, z, t, and color but keeping the x coordinate fixed.Γ µ are the staggered equivalents of the Euclidean Dirac matrices, defined as in [17], where η µ (x) are the staggered phases, and u f µ (x) are the U(1) gauge links, defined in Eq. (3.1).To avoid confusion with the x coordinate, we denoted the spacetime dependence of the Γ-matrices as x ≡ (x, y, z, t).Note that in Eq. (3.5) we made the dependence of Γ µ on the quark flavor implicit.The total electric current density due to all quark flavors is given by Note that this definition of the current is consistent with Eq. (2.4).
For our magnetic field profile, the only non-zero component is for µ = y.In Sec.5.2, we will also discuss the correlator between the electric current and the Polyakov loop and between the chiral condensate and the Polyakov loop.These operators can be defined in a similar fashion [8] In this setup, we computed the electric current operator in Eq. (3.6) using four different lattice spacings in order to study its continuum limit.To resolve the x-dependence in the fermionic traces, we employed random sources in the same way as in Ref. [8].Our parameter set covers a range of temperatures from 113 MeV to 300 MeV, and various values for the magnetic flux quantum N b .Moreover, we fixed the physical value of the inhomogeneity parameter to ϵ ≈ 0.6 fm, which is compatible with what one would expect from model simulations of off-central heavy-ion collisions [2].For the renormalization of some observables, we also used T ≈ 0 gauge configurations.In Table 1, we show the simulation parameters used in the full-current setup.
T > 0 ensemble Table 1: Simulation parameters for the full-current setup.To keep the physical value of the magnetic profile width fixed at ϵ ≈ 0.6 fm, we tuned the ratio ϵ/a for every temperature.

Valence-current setup
Above we have seen that the magnetic field appears in the current expectation value (3.4) in two different ways: in the current operator (valence contribution) as well as in the fermion determinant (sea contribution).In Sec.5.2 we will demonstrate that the former contribution is by far the dominant one and therefore it is a good approximation to only include the magnetic field in the operator.For this valence approximation, we also employed rooted staggered fermions at the physical point with two stout smearing steps.We start by decomposing the same staggered current into operators that are relevant to quantify the impact of sea and valence effects in the weak-field regime, ) where the superscripts indicate the valence and sea terms, respectively.We can also define the valence and sea currents due to all flavors as and (3.11) respectively.For a weak magnetic field, we can show that these currents follow an additivity relation (see Appendix A for more details) From now on, we refer to the terms on the right-hand side as valence and sea currents.
On the lattice, the sea current corresponds to a B-independent operator computed on configurations at B ̸ = 0, whereas the valence one corresponds to a B-dependent operator computed on gauge configurations at B = 0.This makes the valence term computationally much cheaper than the sea, since we only need a single ensemble at B = 0 to resolve the B-dependence of j val µ .Notice that both currents contain all orders in B. The O(B 3 ) notation in Eq. (3.13) simply emphasizes that the decomposition of j µ is only exact at leading order in the magnetic field.In Sec.5.2, we empirically show that j sea µ is negligible for a broad range of temperatures and up to rather strong magnetic fields.This result justifies what we call the valence approximation: j µ ≈ j val µ , where the dominant contribution to the full current is captured by the valence term.
Within the valence approximation, we explored two cases: • Constant ϵ: to validate the valence approximation, we computed the valence current using the same set of parameters from Table 1, i.e. setting ϵ ≈ 0.6 fm, and compared it to the full current.We will discuss the corresponding results in Sec.5.2.
• Constant ϵ/L: to have a better control over finite volume effects at hight T , we computed the valence current using the parameters from Table 2, i.e. with fixed ϵ/L = 0.12 for all T and lattices.This ensures that the inhomogeneity of the magnetic field fits into the box regardless of the parameters.Note that the susceptibility defined via Eq.(2.13) does not depend on ϵ/L in the limit of B → 0. In this case, we used existing ensembles from [14].We will discuss the corresponding results in Sec.5.3.
The reason for studying two cases was to understand whether boundary effects become important at high T .Since the physical value of ϵ is fixed in the first case of the valence setup, as well as in the full-current setup, the ratio ϵ/L increases at high T , where due to the smaller lattice spacings, the physical volume decreases.This could mean large contributions to the susceptibility coming from the edges.Therefore, the second case of the valence setup was designed to check if the high-T part of the susceptibility is affected by the boundary effects.
At this point, we make two important remarks on the computation of the susceptibility in the valence setup.First, using the sea-valence expansion in Eq. (2.12) yields Therefore, the decomposition of the susceptibility is also of the form χ m = χ val m + χ sea m .In other words, the accuracy of using the valence approximation to compute the susceptibility is only limited by the magnitude of the sea term and is independent of the higher-order terms, for they vanish in the B = 0 limit.Second, since the valence currents for different values of B are computed on the same gauge configurations, we employed a correlated fit when carrying out the B = 0 extrapolation to account for the correlation of the data.
T > 0 ensemble Table 2: Simulation parameters for the valence-current setup, where we kept the ratio ϵ/L fixed at 0.12.In this case, ϵ/a = 0.12N s only depends on the number of sites in the spatial direction.

Renormalization scheme
In this section, we will discuss the renormalization of our observables.It is well known that the magnetic susceptibility contains an additive divergence proportional to the logarithm of the lattice spacing.This divergence stems from the multiplicative renormalization of the bare electric coupling [13,18].Since the divergence in χ m is temperature-independent, to obtain a finite result we can define the renormalized susceptibility as As pointed out in Ref. [14], this scheme is compatible with the requirement that the magnetic permeability µ equals unity in the vacuum.From Eq. (2.4), it is clear that the induced electric current inherits the divergence from the susceptibility.This is not what one would a priori expect from a conserved current in a quantum field theory since these, in general, avoid renormalization, with only few exceptions (see, for instance, Ref. [19]).Nevertheless, this is indeed the case for the electric current operator in consideration.To clarify why, let us go back to the separation (2.2) of the total current into free and magnetization currents.Suppose that our system is subject to an external magnetic field H, induced by the free current j H .The medium responds to the external field with a magnetization M, carried by the electric current j.An external observer, however, does not detect j H or j individually, but only the total current j B , which indeed does not renormalize.By separating the medium response j from the free current j H , one finds ultraviolet divergent terms in both, nevertheless, they mutually cancel and give a finite total current.This means that although j is a conserved current, it is infinite in the ultraviolet due to the fact that we separated it from the free current.
Having clarified this aspect, we proceed to find the renormalization prescription for the induced current.Using Eq. (2.12) and the definition of the renormalized susceptibility, we find that from which we conclude that a suitable definition of the renormalized current is Note that this form of the renormalized current is consistent with the fact that the local susceptibility in (2.5) renormalizes with a contact term, χ(x ′ ) − Lχ m (T = 0) • δ(x ′ ) or, that the divergent part of the susceptibility in Fourier-space (2.6) is momentum-independent.We stress that Eq. (4.3) is not equivalent to a subtraction of the T = 0 current.The latter would imply a subtraction to all orders in the magnetic field, while the correct renormalization only affects the term proportional to ∂B(x)/∂x.

Results
In this section, we discuss our results for the electric current density, and the magnetic susceptibility derived from it.To clearly distinguish the two approaches used, we start by showing the results from the full-current in Sec.5.1 and later we discuss the valence-current approach in Secs.5.2 and 5.3.

Full-current results
In Fig. 2, we show the bare electric current densities on the lattice computed using Eq.(3.4) at a given temperature for two magnetic field strengths.Note that, from Eq. (3.3), the physical magnetic field eB is different for the ensembles with slightly different aspect ratios.In particular, for the 28 3 × 10 lattice, which has aspect ratio 2.8 instead of 3, the physical magnetic field is slightly higher than the other two lattices.
As discussed in Sec. 2, these profiles can be combined with Ampère's law to obtain the magnetic susceptibility.The magnetic flux quantization only allows us to access discrete values of N b .Therefore, to carry out the zero-B limit, we fitted a quadratic function to the ratio in Eq. (2.12).In Fig. 3, we illustrate the extrapolation of the bare magnetic susceptibility to B = 0 for two lattice spacings, as described in Sec. 2, and the corresponding temperature dependence.
For some temperatures, we found that the data points at N b = 1 show a slight deviation from the expected quadratic behavior.This issue can be explained in the following way: it has been shown that finite-volume effects induce sinusoidal modulations in some thermodynamic quantities even at uniform background magnetic field [20], where the frequency of the oscillation increases with N b .In particular, we found this to also hold for the electric current.At N b = 1, the finite-volume contamination overlaps with the wedge function in Eq. (2.12), and its signal is therefore enhanced by it.This, however, was found to be unproblematic for the B → 0 extrapolations and did not have a substantial impact on the renormalized susceptibility.
To renormalize the current and the susceptibility, we needed to carry out the extrapolation to B = 0 for the T = 0 ensemble and subtract the corresponding divergent contribution.In Fig. 4, we show the result of Eq. (2.13) at T = 0 and the linear fit in log(a/a 0 ), where a is the lattice spacing and a 0 ≡ 1.46 GeV −1 is a normalization factor.
Using χ m (T = 0) obtained from the linear fit, we renormalized the current as described in Sec. 4. Our continuum limit extrapolation was done by fitting the lattice data in x, B, and the lattice spacing simultaneously using a Monte Carlo approach.We carried out this procedure for the chiral condensate and Polyakov loop in Ref. [8].For more details on  the method, see [21].In Fig. 5, we show the continuum limit surface fitted to the data to illustrate the continuum extrapolation procedure.In Fig. 6, we show slices of the continuum limit surface as a function of x at specific magnetic field strengths.In Appendix B, we also compute the induced electric current in chiral perturbation theory from the momentum-  Notice that the peak of the electric current appears around x ∼ ±0.28 fm, where the curl of the magnetic field is maximal.In Fig. 7, we show the behavior of this peak as a function of the magnetic field for various temperatures.In Fig. 8, we show the renormalized susceptibilities at nonzero lattice spacings and their continuum limit extrapolation.At high T , we see that χ R m > 0, indicating a strong paramagnetic phase, whereas at T < T c , the sign of χ R m changes and we observe a weak diamagnetic behavior.This is in agreement with the lattice results obtained with a different method [14].

Valence-current results at fixed ϵ
In this section, we present our results employing the sea-valence decomposition, as defined in Sec.3.3.In Fig. 9, we compare the magnitude of the unrenormalized valence current with the full one for different temperatures and magnetic fields.In Fig. 10, we also show the bare sea contribution to the current for various magnetic field values on one particular lattice.
Figs. 9 and 10 provide empirical evidence that the valence term captures most of the physics of the electric current in the weak-field regime.Since the sea term quantifies the effect of quark loops in gluonic interactions at low B, its tiny magnitude suggests that the electric current operator is very weakly sensitive to the gluon fluctuations.To further quantify this effect, we define a correlator between the bare electric current and the bare Polyakov loop.The latter operator encodes the most relevant properties of the gluon fields in the thermal medium.The correlator reads where we normalized by the cubed pion mass to have a dimensionless combination.The integration is carried out on half of the x-direction so that C j is non-zero.This correlator is shown in Fig. 11 to demonstrate that it vanishes within errors.
For comparison, we can define the analogous correlator between the Polyakov loop and the chiral condensate (see [8]).
where ψψ is the average over the light-quark condensates.Fig. 11 shows that the magnitude of C ψψ is substantially larger, peaking around the transition temperature where correlations are most prominent.We note moreover that the low sensitivity of the current to gluons also manifests itself in a reduction of statistical uncertainties associated with this operator on a given ensemble.For instance, we found that the magnitude of the errors σ are typically much lower in the case of the current as compared to the chiral condensate, which has a non-vanishing sea contribution.In Fig. 11, we also include a comparison between C j and C ψψ and the ratio σ ψψ /σ j as functions of T .For concreteness, we note that for the same ensemble, the condensate can have up to ∼ 10 times larger statistical errors than the current.When analyzed on a single gauge configuration, σ j and σ ψψ behave rather similarly as a function of the number of random vectors, indicating that, indeed, the factor of 10 difference between them is better explained by the gauge noise.

Valence-current results at fixed ϵ/L
The procedure for the valence-current method is entirely analogous to the one of the fullcurrent, except that now we take j y as the valence term in Eq. (3.9).In Fig. 12, we compare the susceptibilities estimated from the full and valence currents, as well as an alternative method in the literature employing current-current correlators [14].Here we also include the prediction of the hadron resonance gas (HRG) model [13,22].

Discussion
In the previous section, we investigated the appearance of an electric current profile at non-uniform magnetic fields using three approaches: 1. Employing the full electric current operator at fixed ϵ, 2. Making use of the valence approximation at fixed ϵ, and 3. Making use of the valence approximation at fixed ϵ/L x .The precise agreement between approaches 1 and 2 on the electric current throughout a broad range of temperatures confirms the validity of this new approximation scheme developed in this work.It also paved the way for our third approach, which was used to compute the magnetic susceptibility in the valence approximation scheme.Interestingly, the valence approximation suggests that in a magnetized medium, valence quarks are more efficient carriers of the electric current than sea quarks.As we have shown, this physical picture is supported by two of our results from the previous section: 1. the tiny magnitude of the sea contribution to the total current, as seen in Fig. 10, even at strong magnetic fields and 2. the fact that the correlator between the current and the Polyakov loop vanishes within errors.This finding also reveals the low sensitivity of the electric current to the gluonic fluctuations in the medium, which explains the smaller statistical errors on the computation of ψΓ µ ψ compared to the chiral condensate on a given ensemble.in the valence setup as a function of temperature.The points represent the lattice data.Right plot: comparison between the magnetic susceptibility estimated using the full-current and the valencecurrent methods.In addition, the continuum extrapolated results of [14] using an alternative method (yellow band), as well as the result of the HRG model [13] (black line) are also included.
Another interesting behavior of the chiral condensate in the presence of magnetic fields that is completely absent in the electric current is the so-called inverse magnetic catalysis [23].Around T c , sea effects dominate over valence effects and the chiral condensate is known to be supressed by the field.However, as shown in Fig. 7, the peak of the electric current is a monotonic function of B for all temperatures, where the dependence on the temperature was also found to be mild.This is directly related to the fact that sea effects on the electric current are negligible at low B and correlations with the Polyakov loop are practically absent.
At stronger B, however, the valence approximation breaks down.For instance, in Figure 9, at T = 200 MeV, we start to see a small deviation from the full current at eB ≳ 0.9 GeV 2 .Nevertheless, the sea term alone cannot explain the difference between the full current and the valence current, as one can see in Fig. 10.Therefore, the origin of the slight discrepancy between full current and valence current is attributed to higher-order contributions in the sea-valence expansion (see Eq. (A.7)).
A question that arises is whether one can extend this approximation to other B-odd operators, for instance those related to anomalous transport phenomena [24].Moreover, due to the low computational cost of the valence current computation, the approach we developed in this work might be easily implemented within other fermion discretizations to study the weak-field regime of operators.In addition, the method may also be extended to arbitrary profiles of the magnetic background.

Conclusions
In this article, we studied the emergence of local electric currents in QCD in the presence of magnetic fields with non-zero curl using lattice simulations.We extrapolated our results to the continuum limit and determined the spatial structure of these currents, as well as their magnetic field and temperature dependence.As discussed in Sec. 1, such currents are allowed by Bloch's theorem, as long as the their integral over the full volume vanishes, which is the case at hand.Furthermore, using the sea-valence expansion of the electric current, we empirically established that sea quark effects are negligible for its weak-field behavior.This lead us to the valence approximation as a computationally comparably cheap technique.Furthermore, we determined that the sea contribution to the electric current remains very close to zero even at strong fields.Therefore, we attributed the slight discrepancy between the full current and the valence current to higher-order terms in the sea-valence expansion.
Employing the full current operator, and also the valence approximation, we showed two novel ways of computing the magnetic susceptibility of QCD applying Ampère's law.This method draws a direct connection between the susceptibility (a two-point function) and the electric current (a one-point function), namely Eq. (2.12), which is only possible due to the inhomogeneity of the field.In this sense, the current-current correlator method of [14,25] may be viewed as the exact weak-field expansion of the local currents in our approach.
The magnetic susceptibility obtained from the full-current method and from the valence approximation agree with the results obtained in [14] on the quantitative level.We thereby confirmed that the susceptibility changes sign across the QCD crossover, going from diamagnetic to paramagnetic matter, which pinpoints the transition from a system dominated by hadrons to a system with quasi-free quarks.
In the context of off-central heavy-ion collision phenomenology, the benchmarking of the initial magnetic field strength has recently shifted to the focus of interest, with various suggestions for observables that may be measured both experimentally as well as in lattice simulations, see e.g.Ref. [26].The induced currents that we explored in this paper are paramagnetic at high temperatures, i.e. these reinforce the background magnetic field and might impact on observables tied to the early stages of the collisions.They may be taken into account in a similar manner as currents originating from Faraday and Hall effects [3].
In accordance with the definitions (3.9) and (3.10) in the main text, the valence and the sea terms for the operator O correspond to the following path integrals The valence term corresponds to measuring the B-dependent operator O(B) on B = 0 configurations, while the sea term corresponds to measuring the B = 0 operator O(0) on B-dependent configurations.
To apply this procedure to the electric current operator, we first make contact with the expectation value of j µ by defining O ≡ Tr Γ µ M −1 inside the path integral.Then, the sea-valence decomposition reads However, the electric current expectation value is an odd function of B, hence ⟨ j µ (0) ⟩ 0 = 0.Moreover, both j val µ and j sea µ are also odd in B due to parity symmetry, therefore which is Eq.(3.13) of the main text.Notice that for parity-even observables, like the quark condensate, all three terms in the right hand side of Eq. (A.7) are O(B 2 ) and there is no order in B where the third term may be neglected.

B Induced currents in chiral perturbation theory
The low-energy effective theory of QCD is chiral perturbation theory (χPT).In this appendix, we calculate the induced currents due to the inhomogeneous magnetic field within this framework.In χPT, the effective degrees of freedom are the pions, among which the charged ones are responsible for the electric current.These are described, to leading order, via scalar QED coupled to the background eletromagnetic field.The theory is thus defined by the Minkowski action where D µ = ∂ µ − ieA µ is the covariant derivative in the presence of the background field A µ and m π the pion mass.We turn the argumentation of Sec. 2 around: we calculate the momentum-dependent magnetic susceptibility of the charged pions, and use that to build the inhomogenous current with the magnetic field profile used in the main text.In our choice of gauge (A 2 (x) to represent the magnetic field), the magnetic susceptibility can be written as where Π 22 is the photon polarization tensor of scalar QED and p denotes, just like in Sec. 2 of the main text, the momentum in the x direction.The calculations to be made are text-book, and can be found in e.g.Chapter 16 of [27].Two one-loop order diagrams contribute: a momentum-dependent bubble-type, and a momentum-independent tadpoletype.However, when putting all external momentum components, except for p, to zero, the final result is considerably simplified, even though both diagrams contribute.In this case the bare polarization tensor in dimensional regularization reads where γ E is the Euler-Mascheroni constant and µ the renormalization scale.As well known [22,28], the prefactor of the divergent term is proportional to the leading-order β-function coefficient β 1 = 1/(48π 2 ) of scalar QED.Notice that this singular term leads to a momentum independent divergence for the susceptibility itself, which we can take care of in a similar manner as in the case of full QCD, by a zero-momentum subtraction.This also sets the renormalization scale implicitly, such that Calculating the above current using the physical pion mass at weak magnetic fields, one should find reasonable agreement with the continuum extrapolation of our lattice simulation results.Here, we follow a slightly different route instead and incorporate staggered lattice artefacts into the χPT description.Following [14,29], we average over the 16, partially degenerate taste copies of pions appearing in staggered lattice simulations, of which only the lightest one has the physical value.The taste splitting of the masses for our action is taken from [30].
Our result for the χPT prediction of the current is shown in Fig. 13, where we plot j R y (x) for both the continuum setup with the physical pion mass, and for the above described tasteaveraged setup.The curves are compared to our lattice data obtained on a 24 3 × 32 lattice at the lowest possible magnetic field, eB = 0.88 GeV 2 with a lattice spacing of a = 0.29 fm.We observe a remarkable agreement between the lattice data and the taste-averaged χPT current.In turn, the continuum curve with physical pion mass lies above the lattice results by about a factor of two.This shows that lattice artefacts on this coarse lattice are substantial, but still in a region where a continuum extrapolation is under control.

Figure 1 :
Figure 1: Magnetic field profile (upper plot) and its first derivative (lower plot) in lattice units for N b = 6 as a function of the x coordinate.The connecting line are only included to guide the eye.

Figure 2 :
Figure 2: Left plot: bare electric current profiles as functions of x at T = 113 MeV and N b = 1.Right plot: similar plot at N b = 6.

Figure 3 :
Figure 3: Right column plots: χ m as a function of eB obtained from Eq. (2.12) with quadratic extrapolation (gray bands) for different temperatures and two lattice spacings.The width of the bands depict the statistical error of the fits to the data on the right-hand-side.Left column plots: projection of the zero-magnetic-field χ m on the temperature axis for two lattice sizes.The gray bands, representing the statistical error of the fits, were extended to the T -axis to guide the eye.

Figure 4 :
Figure4: Bare magnetic susceptibility at zero-temperature and zero-field obtained in the fullcurrent setup with linear perturbative fit as a function of log(a/a 0 ). a is the lattice spacing and a 0 ≡ 1.46 GeV −1 is a reference lattice spacing, as defined in Ref.[14].

Figure 5 :
Figure 5: Continuum limit surface and renormalized electric current data at T = 113 MeV for N b = 1, 2, 3.The higher N b values were included in the global fit but omitted from the plot for clarity.

Figure 6 :
Figure 6: Continuum limit of the renormalized electric current profiles projected on the x direction for two magnetic fields at T = 124 MeV.The points represent two lattices used in the multidimensional continuum extrapolation.The dashed lines are the results of the splines at the corresponding lattice spacing.

Figure 7 :
Figure 7: Renormalized electric currents evaluated at x = 0.28 fm (at the peak position) in the continuum limit as a function of the magnetic field for three temperatures.

Figure 8 :
Figure 8: Continuum limit extrapolation of the magnetic susceptibility (green band) in the fullcurrent setup as a function of temperature.The colored points represent the data for different lattice spacings.

Figure 9 :
Figure 9: Comparison between valence and full currents for various magnetic fields at temperatures T < T c and T > T c .In each plot, we indicate the magnetic field strength, in physical units, corresponding to the highest N b shown in the plot, i.e.N b = 4.

Figure 10 :
Figure 10: Sea contributions to the electric current density at T = 113 MeV (left) and at T = 155 MeV (right) as functions of the x-coordinate for different magnetic field strengths on a 28 3 × 10 lattice.The sea term fluctuates closely to zero even in the strong-field regime.

Figure 11 :
Figure 11: Upper left panel: correlators between the Polyakov loop and the electric current (blue) and between the Polyakov loop and the chiral condensate (red) as functions of T .Lower left panel: ratio between the average statistical errors for the chiral condensate (σ ψψ ) and the electric current (σ j ) on this ensemble.The dashed line represents σ ψψ /σ j = 1.Right panels: average error coming from the random vectors on one gauge configuration for the current and chiral condensate at T = 124 MeV and T = 155 MeV as a function of the number of random vectors.

Figure 12 :
Figure 12: Left plot: continuum limit extrapolation of the magnetic susceptibility (purple band)

Figure 13 :
Figure 13: Comparison of our lattice data for the renormalized current with chiral perturbation theory estimates.A correct description at finite lattice spacings needs a correction for the breaking of chiral symmetry on the taste level.