Magnetic susceptibility of QCD matter and its decomposition from the lattice

We determine the magnetic susceptibility of thermal QCD matter by means of first principles lattice simulations using staggered quarks with physical masses. A novel method is employed that only requires simulations at zero background field, thereby circumventing problems related to magnetic flux quantization. After a careful continuum limit extrapolation, diamagnetic behavior (negative susceptibility) is found at low temperatures and strong paramagnetism (positive susceptibility) at high temperatures. We revisit the decomposition of the magnetic susceptibility into spin- and orbital angular momentum-related contributions. The spin term -- related to the normalization of the photon lightcone distribution amplitude at zero temperature -- is calculated non-perturbatively and extrapolated to the continuum limit. Having access to both the full magnetic susceptibility and the spin term, we calculate the orbital angular momentum contribution for the first time. The results reveal the opposite of what might be expected based on a free fermion picture. We provide a simple parametrization of the temperature- and magnetic field-dependence of the QCD equation of state that can be used in phenomenological studies.

1 Introduction The development of a quantitative and precise understanding of the response of QCD matter to background (electro)magnetic fields is of vital importance for furthering our knowledge about a multitude of physical systems. Examples include the interior of magnetars, neutron star mergers [1][2][3], off-central heavy-ion collisions and the evolution of the universe in its early stages. For general reviews, we refer the reader to Refs. [4][5][6]. A characteristic feature of the behavior of strongly interacting quarks and gluons is rooted in the dependence of the QCD equation of state (EoS) on the background magnetic field B. The EoS enters all the above mentioned examples: it appears in the gravitational stability conditions of compact stars, affecting the mass-radius relation [7]; it governs the expansion rate in cosmological models [8,9] and it also sets the conditions where freeze-out is reached in heavy-ion collisions, see, e.g., Ref. [10]. While the equilibration of fireballs produced in heavy-ion collisions is still a subject of research, at least in astrophysical systems the time and distance scales over which the magnetic field varies are much larger than those that govern QCD processes that affect, e.g., the EoS or nucleo-synthesis. With these applications in mind, solving QCD in a constant background magnetic field is sufficient and this scenario is amenable to lattice simulations.
The leading dependence of the EoS on B is encoded in the magnetic susceptibility χ of QCD matter. Its sign distinguishes between paramagnets (χ > 0), for which the exposure to the background field is energetically favorable, and diamagnets (χ < 0), which repel the external field. In QCD matter, like in any other material, the origin of the magnetization and hence of the magnetic susceptibility is related to the spin and angular momentum of charged particles. At high temperatures the quarks are the relevant degrees of freedom; at low temperatures the hadrons and in particular the pions take over, while (valence and sea) quarks contribute just as their fundamental constituents.
The total angular momentum that gives rise to the magnetization can be decomposed into contributions from the spins of the quarks of different flavors and a remainder. The latter contains the quark orbital angular momenta but also the angular momentum of the gluons, that can split into quark-antiquark pairs. 1 The quark spin contribution to the magnetic susceptibility is due to the expectation value ψ f σ µν ψ f , whose leading order response is linear in the magnetic field strength tensor F µν . Depending on the normalization, the slope is proportional to the magnetic susceptibility of the quark condensate [14][15][16], or the so-called tensor coefficient, i.e. the normalization of the photon lightcone distribution amplitude (DA) [17][18][19] of finding a quark-antiquark pair of flavor f in a transverse photon. This in itself appears in a multitude of applications, e.g., as a correction to the hadronic light-by-light scattering contribution to the muon anomalous magnetic moment [20,21], within radiative transitions [22,23] and in the photo-production of mesons [24].
Some of the present authors have already addressed several of the above aspects in Refs. [25][26][27][28]. Here we improve on these studies by employing a novel calculational method that is based on Ref. [29], by carrying out the QCD renormalization non-perturbatively with respect to the intermediate RI'-MOM scheme [30,31] and by adding a finer lattice spacing. In addition we present and exploit new analytical findings. One of the outcomes will be that in the strongly interacting medium at low to moderately high temperatures the quark spin-related susceptibility is negative (diamagnetism) while the part that is due to the orbital angular momentum is positive. Clearly, this behavior is very different from the response to magnetic fields of the materials that have so far been accessible to solid state physics experiments. Therefore, our results offer a glimpse into a completely new regime of spin physics.
The article is organized as follows. In Sec. 2 we introduce our notations and the central observables. For conceptual clarity, we carefully address their divergence structure and renormalization in QED and QCD. In Sec. 3 we then discuss details of the simulation and, in particular, we introduce our new method that employs current-current correlators in a mixed coordinate-and momentum-space representation. This enables us to determine susceptibilities from lattice simulations at B = 0. We then present and discuss our results in Sec. 4, before we summarize. We include several technical appendices: in App. A we investigate the effects of taste splitting in the staggered formulation within the hadron resonance gas model. This turns out to be important to avoid underestimating the systematics of the continuum limit extrapolation. In App. B we derive the factorization of the susceptibility into quark spin-related and other contributions, building upon earlier partial results [26]. In the extensive App. C, several derivations are carried out for the free case, establishing, e.g., the structure of QED divergencies. App. D discusses the non-perturbative renormalization procedure. In App. E we present more detail on the derivation of the new current-current method and compare to numerical results, obtained using conventional background field approaches. Finally, App. F gives a parametrization of our results for the QCD EoS for a broad range of temperatures and magnetic field strengths. The corresponding Python script param_EoS.py is uploaded to the arXiv as ancillary file along with this paper.

The response of QCD matter to background fields and the magnetic susceptibility
Without any loss of generality, below we consider a magnetic field pointing in the x 3 direction, with the magnitude B. The magnetic susceptibility is defined via the leading (quadratic) dependence of the QCD free energy density f on the field strength, where Z is the partition function, T the temperature and V the spatial volume of the system. The product eB of the elementary electric charge e and the magnetic field is a renormalization group invariant due to the QED vector Ward identity, see Eq. (2.3) below. Therefore, χ b is free of multiplicative renormalization. Still, the susceptibility undergoes additive renormalization, which is made explicit by the index b, indicating the bare quantity, as obtained in the lattice scheme at a lattice spacing a. 2 This was discussed in Ref. [28] in depth, but we repeat the argument here for comprehensiveness.

QED renormalization
The total free energy density of the system, that includes both QCD matter and the classical background field, is a physical observable and therefore free of divergences. However, the second term within Eq. (2.2), involving the bare magnetic field B b , contains a logarithmic divergence in the lattice spacing a due to electric charge renormalization [32], where the renormalized quantities e and B depend on the QED renormalization scheme and on the QED renormalization scale µ QED . Since the background field is classical, only the leading-order QED β-function coefficient β 1 appears here [33]. Note that β 1 is affected by QCD corrections at the cut-off scale, where g is the strong coupling and q f denotes the electric charge of the quark flavor f . The coefficients c i of the perturbative series are known up to i = 4 in the MS and MOM schemes [34] and c 1 = 1/(4π 2 ) is universal for massless schemes. These QCD corrections 3 vanish logarithmically with the lattice spacing towards the continuum limit due to the asymptotically free nature of the strong interactions, as is also indicated in Eq. (2.4). Eq. (2. 2) implies that f contains the same additive divergence as B 2 b /2, but with an opposite sign. This propagates into the susceptibility (2.1), resulting in The renormalized susceptibility χ depends on the renormalization scale, which we indicated here explicitly in square brackets. We confirm the presence of the logarithmic divergence in χ b analytically for the free case in App. C.2 and numerically for full QCD in Sec. 4.1. The divergence is independent of the temperature so that it cancels within the difference This definition of the renormalized susceptibility -implying that it vanishes identically at T = 0 -corresponds to a particular choice of the renormalization scale µ QED = µ phys QED . In fact this is the only prescription that adheres to the physical requirement that the magnetic permeability (1 − e 2 χ) −1 should be unity in the vacuum. In the following we suppress the dependence on the QED renormalization scale and simply write χ for the susceptibility, renormalized in this way.

The tensor coefficient
Besides the magnetic susceptibility there exist further quantities that characterize the leading-order response of the QCD medium to the background magnetic field. For a general background field F µν , the fermion bilinear involving the relativistic spin operator σ µν develops a nonzero expectation value [14,35], We will refer to τ f b as the tensor coefficient for the flavor f . Similarly to χ b , this is also a bare observable that contains additive logarithmic divergences in the cut-off. For our choice of direction of the magnetic field the tensor coefficient can be determined as the slope of the expectation value of the fermion bilinear involving σ 12 at small values of the magnetic field: The tensor coefficient contains a similar logarithmic divergence as χ b . We demonstrate the reason for this in Sec. 2.3 below. In particular, the divergence structure takes the form [26], where m f is the mass of the quark of flavor f . Again, due to asymptotic freedom, QCD corrections to γ τ 1 vanish in the continuum limit. 4 Eq. (2.9) is confirmed in App. C.4 for the free case and checked numerically in full QCD in Sec. 4.2. Notice that in the chiral limit the divergent term disappears, so that lim m f →0 τ f b is ultraviolet-finite. We will carry out this limit at zero temperature in Sec. 4.2 below. In this situation, up to multiplicative renormalization, this object corresponds to the normalization f ⊥ γ of the leading-twist photon distribution amplitude [17][18][19], i.e. of the infinite momentum frame probability amplitude that a real photon dissociates into a quark-antiquark pair of flavor f .
In analogy to Eq. (2.6), we define the renormalized tensor coefficient by subtracting its value at zero temperature, which again corresponds to a particular choice of the QED renormalization scale. Unlike for χ, there is no preferred choice in this case, but it is natural to use the same prescription as in Eq. (2.6) above. Besides the (QED-related) additive renormalization detailed above, the tensor coefficient also undergoes (QCD-related) multiplicative renormalization by the tensor renormalization constant Z T . This introduces a further scheme-and scale-dependence of this observable. Below we will consider Z T in the MS scheme at the QCD renormalization scale µ QCD = 2 GeV. Unlike in our previous study [26], where we calculated Z T perturbatively at the one-loop level, here we carry out a nonperturbative matching to to the RI'-MOM scheme [30,31] and subsequently translate the result at three-loop order [36] into the MS scheme. This procedure is detailed in App. D. We remark that Z T is independent of the temperature, thus the ordering of the QED renormalization (i.e. the T = 0 subtraction) and the QCD renormalization (multiplication by Z T ) is irrelevant for the determination of the renormalized tensor coefficient τ f .

Decomposition into spin and orbital angular momentum contributions
One might suspect that χ and τ f are not completely unrelated. Indeed, as we first discussed in Ref. [26], τ f represents the contribution of the spin of the quark flavor f to the total magnetic susceptibility. In particular, χ can be decomposed into spin-related and orbital angular momentumrelated contributions, 5 χ = χ spin + χ ang , (2.11) and the spin term is related to the tensor coefficients as where m f denotes the quark mass in the lattice scheme and Z S and Z T are the scalar and tensor renormalization constants, respectively. The second term in the square brackets is understood to correspond to the limit of a vanishing valence quark mass taken at physical, i.e. nonzero, values of the sea quark masses. We discuss the difference between valence and sea quark masses in Sec. 3 below.
In App. B we prove Eqs. (2.11)-(2.12) and illustrate the origin of the subtraction of the chiral valence quark limit. App. C also contains an explicit check of Eq. (2.12) in the free case. A remark regarding the choice of renormalization scales is in order here. Eq. (2.12) is chosen so that χ spin vanishes at T = 0. Recall however that, according to our remark below Eq. (2.10), we are free to choose an arbitrary QED renormalization scale for the renormalized tensor coefficient. This freedom propagates into χ spin and -through the decomposition (2.11) -to χ ang as well. In contrast, the renormalization scale for χ is fixed by the requirement χ(T = 0) = 0. Thus, in principle both susceptibility contributions may be shifted by an arbitrary amount, as long as their sum remains zero at T = 0. We will follow the choice made in Eq. (2.12), which corresponds to setting χ spin (T = 0) = χ ang (T = 0) = 0. In the free case this is realized by choosing one and the same QED renormalization scale for all susceptibility contributions, see App. C.3. Our numerical results in full QCD below suggest that also in the interacting case the QED scale that corresponds to the renormalization condition χ spin (T = 0) = 0 is consistent with the one obtained from setting χ(T = 0) = 0.
In Eq. (2.12) we also carried out the QCD related multiplicative renormalization by including the tensor renormalization constant Z T required for τ f (as mentioned above) as well as the scalar renormalization constant Z S = Z −1 m , which multiplies the inverse quark mass. 6 Note that these renormalization factors depend on the QCD regularization scheme and on the QCD renormalization scale. As mentioned above, we choose the MS scheme and µ QCD = 2 GeV. We remark that, just as for τ f , the ordering of the QED and the QCD renormalization is irrelevant for χ spin . We stress again that while the factorization of the total susceptibility χ into χ spin and χ ang depends on the QCD scheme and scale, χ itself is a QCD renormalization group invariant.
In the free case the two susceptibility contributions have a constant ratio, χ spin : χ ang = 3 : (−1), reflecting the well-known response of a free charged fermion to the magnetic field via its spin and its 5 Note that for simplicity we refer to χ ang as the orbital angular momentum contribution. In the interacting case this can be further factorized, separating out the gluon total angular momentum contribution [27] from those of the quark angular momenta. 6 Like in massless continuum schemes, in staggered lattice formulations there is no difference between singlet and non-singlet renormalization factors for these quark bilinears. This was explicitly demonstrated at order g 4 in Ref. [37].
orbital angular momentum, dating back to Pauli and Landau [38,39]. This ratio, which translates into the rule χ spin : χ = 3 : 2, holds identically in the free case, see App. C.3. In contrast, in full QCD it only applies to the divergence structure, which in the continuum limit -as we have seen above -is governed by pure QED physics. Below we determine to what extent the Pauli-Landau decomposition is affected by the strong interactions.

Lattice methods
We consider spatially symmetric N 3 s × N t lattice ensembles, corresponding to the temperature T = 1/(N t a) and the volume V = L 3 = (N s a) 3 . The simulations are performed with the tree-level Symanzik improved gauge action S g and three flavors (f = u, d, s) of stout smeared rooted staggered quarks [40], described by the Dirac operator 7 / D f + m f . The quark masses m f are tuned as a function of the inverse gauge coupling β = 6/g 2 along the line of constant physics: m ud (β) ≡ m u (β) = m d (β) = m s (β)/R with R = 28.15 [41]. The electric charges are set as q d = q s = −q u /2 = −e/3, where e > 0 is the elementary electric charge. The magnetic field enters in / D f via space-dependent U(1) phases. Further details of our setup and of the simulation algorithm are discussed in Ref. [42]. The lattice geometries for our finite temperature lattices are 16 3 × 6, 24 3 × 6, 24 3 × 8, 28 3 × 10 and 36 3 × 12, allowing for the investigation of both finite volume and discretization effects. Our zero-temperature ensembles consist of 24 3 × 32, 32 3 × 48 and 40 3 × 48 lattices.
In the rooted staggered formulation the partition function and the expectation value of the tensor bilinear are written as path integrals over the SU(3)-valued gluonic links U as

Magnetic flux quantization
In an infinite volume a magnetic field pointing in the x 3 direction can be generated by the Landaugauge electromagnetic potential In a finite volume, in order to comply with periodic boundary conditions for the electromagnetic parallel transporters u µf = exp(iq f A µ ), the boundary twist term A 1 = −Bx 2 L δ(x 1 − L) needs to be included as well [43]. In this setup the flux of the magnetic field is quantized according to [44] so that a differentiation with respect to eB -as required in Eq. (2.1) -is not possible in a standard manner. Several methods were developed to overcome this problem on the lattice, including the anisotropy method [27,45], the finite difference method [46,47] and the generalized integral method [28]. These are all based on approximating the derivative numerically using finite differences in the integer variable N B . This requires independent simulations using several values of N B , which increases the computational requirements considerably. Furthermore, an extrapolation N B → 0 becomes necessary, which inevitably introduces systematic uncertainties. An alternative approach is the half-half method [48], which employs a magnetic field profile that is positive in one half and negative in the other half of the lattice -this enables taking the derivative with respect to the amplitude of the field analytically. However, finite volume effects are substantially enhanced in this case due to the discontinuity of the background field at the boundaries, as was demonstrated in Ref. [29].

Current-current correlators
In view of the above, methods are desirable that only involve measurements at B = 0, thereby circumventing the flux quantization problem of the constant background field profile. In Ref. [29] we showed that χ b can be related to correlators of the electromagnetic current: Specifically, for our electromagnetic gauge (3.2) the relevant role is played by the x 1 -dependence of the correlator of the µ = 2 component, where the zero-momentum projection p 2 = p 3 = p 4 = 0 was performed. Using this projected correlator, the magnetic susceptibility reads In Ref. [29], this relation was derived for χ b in momentum space. In App. E we repeat the argument in coordinate space. In the representation (3.6) the current-current correlator is computed for a general electromagnetic field in coordinate space. Only afterwards a Fourier transformation is carried out via the convolution with the quadratic kernel in order to enforce the constancy of the background field. In this way the problem of flux quantization is avoided. Notice that there is a remnant of flux quantization in the formula (3.6), signaled by the cusp in the kernel at x 1 = L/2. However, this cusp has no practical relevance, as in the integral it is multiplied by G(L/2), which is exponentially small. We checked that finite volume effects are indeed small, see Fig. 10 in App. E.
At zero temperature, the result (3.6) is very similar to the zero-momentum limit of the form factor Π appearing in the vacuum polarization tensor Π µν . In fact our result can be expressed as The vacuum polarization function Π has been the subject of intense research as it appears in the hadronic contribution to the muon anomalous magnetic moment, see, e.g., the recent review [49]. In that setting the relevant observable is the second moment of the two-point function of the electromagnetic current (3.4), projected to zero spatial momentum [50]. Exchanging the time coordinate x 4 for the spatial coordinate x 1 , one can obtain χ b in an analogous way in our background field setup. Since this change is irrelevant at T = 0, at zero temperature χ b = Π(0) holds [29]. For nonzero temperatures the Lorentz structure of Π µν (p) is more complicated than in Eq. (3.7) and, in addition to Π, a form factor Π L appears [51,52]. However, in the static case (p 4 = 0) only Π contributes to the spatial components of Π µν so that the first relation of Eq. (3.7) still holds at T > 0.
In App. E we derive a similar representation for τ f b as well. In this case the equivalents of Eqs. (3.5) and (3.6) become 8 (3.8) We remark that the second moment of the photon DA is accessible too, replacing σ 12 by combinations are antisymmetric in indices equal to 1 and 2, symmetrized over all other non-trivial combinations of indices and with all traces subtracted, see, e.g., Ref. [53]. This is beyond the scope of the present work.
In summary, via the relation (3.6) we are able to determine the magnetic susceptibility using direct measurements at B = 0. This is certainly advantageous over calculating the free energy density (which cannot be obtained as a simple expectation value) at nonzero magnetic fields and differentiating it numerically. The similar relation for the tensor coefficient, Eq. (3.8), might also be used to avoid measuring ψ f σ 12 ψ f at B > 0. However, since the latter is a simple one-point function, the gain is not obvious in this case. In App. E we compare the two methods for this observable and conclude that the correlator method indeed gives larger statistical errors. Therefore, we opted to use our earlier results [26] for ψ f σ 12 ψ f .

The magnetic susceptibility
The correlator G(x 1 ) is calculated using O(1000) random sources located on three-dimensional x 1slices of our lattices. Subsequently, G(x 1 ) is convoluted with the quadratic kernel to calculate χ b via Eq. (3.6). The results of this method, compared to our old data (generalized integral method) and also to those of Ref. [47] (finite difference method) are shown in Fig. 1 at zero temperature. 9 Within errors perfect agreement between the three groups of results is found. The data -plotted in Fig. 1 against log(a) -clearly reflect the logarithmic divergence dictated by Eq. (2.5). Similarly to our fitting strategy in Ref. [29], here we also include the universal perturbative QCD corrections to the QED β-function coefficient c 1 [34], see Eq. (2.4), where g 2 = 6/β is obtained from the inverse lattice coupling β at the lattice scale a −1 . We also take into account O(a 2 ) lattice artifacts so that our fit function reads The result of this fit, with the parameter values is shown as an error band in Fig. 1. We also considered fits with further (quartic) lattice artifacts. The impact of this is included in the second parentheses of Eq. (4.2) for µ QED as a systematic error. The renormalization scale agrees within errors with our earlier determinations [28,29]. It also lies near the mass of the lightest charged hadron (the charged pion) which effectively sets the scale for the 9 Ref. [47] employs the same lattice action. For a different action the bare susceptibilities would not only differ in terms of lattice artifacts but also by an additive constant, due to the different choice of renormalization scheme. The symbols indicate different lattice spacings and the dark orange band the continuum limit. The light orange band represents an estimate of systematic errors of the continuum extrapolation. The dashed gray line is the HRG model prediction [28]. The inset zooms into the low-temperature region to highlight the diamagnetic response there. Right panel: our results for χ (orange, labeled Π(0)) are compared to the results of Ref. [47] (green) and those of Ref. [28] (red) as well as to the HRG model prediction (dashed gray line) and to perturbation theory [28] (dashed light blue band), see Eqs. magnetic response of this system. Nevertheless, note that the value of µ QED depends on the choice of the regulator. The formula (4.1) is used to interpolate χ b and then employed to renormalize the susceptibility at nonzero temperatures according to Eq. (2.6): The results are shown in the left panel of Fig. 2 for a broad range of temperatures and four lattice spacings a = 1/(N t T ) with N t = 6, 8, 10 and 12. A continuum extrapolation is performed by means of a multi-spline fit [54] taking into account O(a 2 ) lattice artifacts. The coarsest lattices (N t = 6 points at T 160 MeV) were found to lie outside the scaling region and were therefore excluded from the analysis. The systematic error was estimated by varying the spline node points and by including/excluding N t = 6 data points at high temperatures. In addition, we consider a further systematic error due to lattice artifacts related to the taste splitting of the staggered spectrum. This effect is particularly relevant at low temperatures and can be estimated by a generalization of the Hadron Resonance Gas (HRG) model that we describe in App. A. The light yellow bands in both panels of Fig. 2 include this uncertainty.
The results demonstrate strong paramagnetism in the quark-gluon plasma phase, in agreement with previous lattice studies [27,28,[45][46][47][48]. At high temperatures the results are well described by the free theory, which predicts (see App. C) where γ = O(1) is a regulator-dependent constant. 10 As indicated, QCD corrections at the thermal scale µ therm affect the leading behavior. In the right panel of Fig. 2 we also include a comparison to this perturbation theory formula, with the scheme-independent O(α s ) corrections to β 1 taken into account [28]. Here we set γ = 1 and use the MS scheme definition of the strong coupling α s = g 2 /(4π), running this at five-loop order [55] to the thermal scale µ therm ∼ 2πT . We employ the central value Λ MS QCD = 0.341 GeV from the recent three flavor QCD determination of α s by the ALPHA Collaboration [56]. The band in the figure corresponds to a variation of the thermal scale from πT to 4πT . The perturbative formula agrees surprisingly well with our results down to temperatures T ∼ 200 MeV.
In contrast to the paramagnetic behavior at high T , towards low temperatures the continuum extrapolated results become negative, revealing a diamagnetic response, previously noted in Ref. [28]. This behavior is in agreement, albeit within large errors, with the HRG model prediction [28] (dashed line in the figures). In the right panel of Fig. 2 we also include results obtained from other approaches that employed the same lattice action. Around the pseudo-critical temperature T c ≈ 155 MeV a significant difference is visible between earlier determinations and our present results. Ref. [47] carried out a continuum extrapolation using fixed β ensembles with lattice spacings a ≥ 0.125 fm. In Ref. [28] we used the fixed N t approach with N t = 6, 8, 10 ensembles, while in the present study N t = 6, 8, 10, 12 lattices are simulated. To highlight the differences between the continuum extrapolations, in the left panel of Fig. 3 we plot the lattice spacing-dependence of the susceptibility for all three approaches. We pick one temperature T = 130 MeV, where the deviation of the continuum estimates is substantial.
The left panel of Fig. 3 reveals the importance of our new N t = 12 ensemble, showing a significant downward trend as the lattice spacing is reduced and a negative value in the continuum limit. The downward trend is not captured by our previous extrapolation using the integral method on N t ≤ 10 lattices [28], neither is it visible in the data of Ref. [47]. We remark furthermore that Ref. [47] performed the continuum extrapolation assuming a strictly positive function for χ(T ). Excluding the possibility of a negative susceptibility in the continuum limit might also bias the extrapolation. In summary, we suspect that the qualitative deviation between the continuum limits below the pseudocritical temperature is due to the absence of N t = 12 ensembles in Ref. [28] and to systematic uncertainties in the a → 0 extrapolation procedure of Ref. [47]. To clarify this issue, dedicated simulations should be performed, preferably at the same temperature and the same values of the lattice spacing.
Finally we provide a parametrization that connects all three approaches (HRG, lattice continuum limit and perturbation theory) and describes χ for arbitrary temperatures. The details are discussed in App. F. In the right panel of Fig. 3 we plot this parametrization, translated to the magnetic permeability µ/µ 0 = (1−e 2 χ) −1 , expressed in units of the vacuum permeability µ 0 . This combination is equal to the ratio of the magnetic induction and the external field, see, e.g., Refs. [28,46].

The normalization of the photon distribution amplitude
Here we address the tensor coefficients τ f b at zero temperature. We consider a set of independent gauge ensembles, generated at the physical value of the strange quark mass m s = m phys s , but at different values of the light quark mass: 0.5 m phys ud ≤ m ud ≤ m phys s . We follow a similar strategy as in Ref. [26], simultaneously fitting the dependence of τ ub · Z T on the light quark mass m ud and on the lattice spacing a according to the ansatz (2.9). Since Z T is found to depend very mildly on the lattice spacing within the range covered (see Fig. 9 of App. D), this does not significantly affect the functional dependence on a. We note that on our coarsest ensembles the uncertainty of Z T is quite large, which also imprints on the errors of the renormalized tensor coefficients. . Light quark mass-dependence of the tensor coefficient at T = 0 for the up quark using our four finest lattice spacings (green to gray symbols). The index b indicates that the QED divergence that one encounters at m ud > 0 has not been subtracted. The results diverge logarithmically towards the continuum limit for any m ud = 0. In contrast, the chiral limit is free of ultraviolet divergences and a combined chiral and continuum limit exists (black circle).
Notice that τ f b diverges for a → 0 for any quark mass, except in the chiral limit, where it is ultraviolet-finite. This tendency is clearly visible in Fig. 4, which shows our results for the up quark. Therefore, we can define an ultraviolet-finite observable for the light quarks, without any zerotemperature subtraction, namely the chiral limit of the tensor coefficient. In contrast, to calculate χ spin we will need to take differences between results obtained at different temperatures (see below).
The ansatz (2.9) contains the free parameters τ f and µ QED . In addition, we include a quadratic mass-dependence and lattice artifacts of O(a 2 ) to each parameter in the fit. Varying the fit ranges in a and in m ud /m phys ud as well as the functional form, we carry out several acceptable fits that are used to build a histogram for the chiral continuum limit of the tensor coefficient. In this combined limit we obtain in the MS scheme The central value differs from our previous result [26] f ⊥ γ = −40.3(1.4) MeV, mainly due to the multiplicative renormalization factor that we determined non-perturbatively here, see Fig. 9 in App. D.
In the left panel of Fig. 5 we show the a-dependence of the T = 0 light quark tensor coefficient τ ub · Z T at the physical point and the result of the above interpolation, including the systematic error estimated using the different fits. For demonstration purposes, we also indicate the leading logarithmic behavior, that we obtain by subtracting the lattice artifact terms from the central fit.
Comparing to the similar plot for χ b (Fig. 1), we see that deviations from the continuum behavior are sizable (and are predominantly due to the fact that we are dealing with a dimensionful quantity in this case). For this reason, here we cannot reliably determine the value of µ QED . Nevertheless, we note that fixing the renormalization scale to its value from Eq. (4.2) also gives acceptable fits. This is in agreement with the expectation of Sec. 2.3, as well as with the results in the free case, see App. C.4. The logarithmic divergence ∝ m f log a becomes more pronounced for heavy quarks. This is visible in the right panel of Fig. 5, where we plot the strange quark tensor coefficient τ sb · Z T at m ud = m phys ud against the lattice spacing and again indicate the leading logarithmic term.
As we have discussed above, at non-vanishing values of the quark mass m f , the tensor coefficient diverges logarithmically. In Ref. [26] we suggested to cancel this by taking the logarithmic derivative with respect to the quark mass, see also Ref. [14]: This renormalization prescription will give identical results for any regulator, up to the multiplicative factor Z T . 11 Using this prescription, it turns out that f ⊥ γu = f ⊥ γd = f ⊥ γ holds within statistical errors. For the strange quark we obtain: The first error includes the described variation of the fit while the second error reflects the uncertainty of the derivative with respect to m s that we indirectly determine from the dependence of the tensor coefficient on the light quark mass, following the procedure explained in Ref. [26].
In the literature often the magnetic susceptibility of the quark condensate, is given, rather than f ⊥ γ = τ u · Z T . Since the latter quantity has a smaller anomalous dimension and its value does not depend on a separate computation of the chiral condensate, this is the preferred choice for practical applications. However, for convenience of comparison, we shall convert it into the other convention. The numerical value of the quark condensate in the SU(2) chiral limit in the MS scheme at the scale µ QCD = 2 GeV reads ψ u ψ u = [272(5) MeV] 3 [58]. To enable a comparison with other results, below we also list X u at the scale µ QCD = 2 GeV. Since most literature values refer to a low, sometimes unspecified scale, in addition we run X u as well as f ⊥ γ to the scale µ QCD = 1 GeV, which is used in most sum rule calculations, see, e.g., Refs. [16,17,19]. This is done, using the five-loop β-and quark mass anomalous dimension γ-functions [55,59] and the three-loop γ-function of the tensor current [36,60]. The results read where we have added all errors in quadrature, including the uncertainty of f ⊥ γ (2 GeV), the difference between running with the two-and three-loop γ-functions of the tensor current, the uncertainty of ψ u ψ u and the uncertainty of the strong coupling parameter [56]. All the above results are in the
Our result (4.6) for the strange quark coefficient translates into where we used the ratio ψ s ψ s / ψ u ψ u = 1.08(17) [57] for the conversion between f ⊥ γs and X s . The difference between X s and X u (1 GeV) ≈ − (475 MeV) −2 has been reported to be negligible in the sum rule calculations [17].

The spin contribution at T > 0
Having interpolated the T = 0 tensor coefficients, we are now in the position to perform the additive renormalization (2.10) by subtracting this contribution from the finite temperature results. We use our existing N t = 6, 8 and 10 results from Ref. [26] to approach the continuum limit. Especially in light of the slow convergence of χ towards a → 0, see the right panel of Fig. 3, this extrapolation should be backed up with finer lattice ensembles in the future. In analogy to the analysis of χ, again we carry out a multi-spline fit of all data sets, determining a systematic error by varying the positions of the spline node points. The so-obtained fit is shown for the up quark and the strange quark in Fig. 6. The results for τ d are consistent with τ u within errors. The large errors of our N t = 6 results at low temperatures are due to the uncertainties of the T = 0 contributions on our coarse lattices, see Fig. 5.
After the additive renormalization, the tensor coefficient vanishes by definition at T = 0. For the light quarks τ u (T ) grows substantially as the temperature is increased, before the slope reduces and a plateau is approached. The inflection point of the continuum curve is found to be at T c = 158(5) MeV. The chiral transition temperature determined from the inflection point of the quark condensate T c = 155(4) MeV [66] is in agreement with this value. For the strange quark pseudo-critical thermal effects set in at somewhat higher temperatures [66]. Also in our case τ s does not appear to exhibit any inflection point, at least for T 170 MeV, and below T ≈ 200 MeV no saturation into a plateau is visible. For sufficiently high temperatures, where the finite quark mass becomes negligible, we expect the two renormalized tensor coefficients to coincide.
Next, the continuum extrapolated results are inserted into Eq. (2.12) to determine the spin contribution χ spin to the susceptibility. To this end we need to evaluate the tensor bilinear for massless valence quarks. Instead of performing measurements at additional valence quark masses, we estimate this limit using the difference between the results for the strange quark and for the light quarks. We assume a linear dependence on the valence quark mass in the range [0, m s ], which implies In this case the contributions of all flavors to χ spin are proportional to τ s − τ u and the renormalized spin susceptibility (2.12) simplifies to (4.13) Thus, in this approximation the individual flavors simply contribute in proportion to their squared electric charges. The scalar renormalization constants entering this expression are displayed in Fig. 9 of App. D. The so-obtained estimate of χ spin is shown in the left panel of Fig. 7 for three lattice spacings, together with a continuum extrapolation performed in the same way as for τ f . We observe χ spin < 0 for all temperatures, with a minimum somewhat above the pseudo-critical temperature and an upward trend for high temperatures. The approximation (4.12) tends to overestimate the valence chiral limit of the tensor coefficient due to the presence of logarithmic deviations from a linear behavior in m val f . 12 Consequently, Eq. (4.13) underestimates χ spin . This is also the case at high temperatures, as can be checked using the analytic formula valid for the free case, see App. C. To take this effect into account we include a systematic error based on the free case formula (C.25). In particular, we consider the difference between the approximation and the true value in the free case and scale it with the typical magnitude of the light quark tensor coefficient at lower temperatures (see Fig. 6). The so-obtained uncertainty is also included in the left panel of Fig. 7.
We remark that χ spin < 0 for the temperature range covered in our simulations. This can be understood by noting that Eq. (4.13) is the discretization of the mass-derivative of τ f . Increasing the mass pushes the inflection point of τ f to higher temperatures (visible in Fig. 6), thus making the derivative negative around the transition temperature. Nevertheless, χ spin will necessarily turn positive for even higher temperatures. Indeed, for sufficiently high temperatures the difference τ f = τ f b (T ) − τ f b (T = 0) will be dominated by the T = 0 term, so that Eq. (4.13) becomes proportional to τ ub (T = 0) − τ sb (T = 0), which is positive for any lattice spacing (see Fig. 4). Perturbation theory also predicts χ spin > 0 for high temperatures, see App. C.3.

Pauli and Landau decomposition of the magnetic susceptibility
Finally, we compare the spin contribution to the total susceptibility in order to learn about the orbital angular momentum-related contribution χ ang = χ − χ spin . All three susceptibilities are included in the right panel of Fig. 7. While the errors of the two contributions are much larger than that of the total susceptibility, several qualitative comments can be made based on this plot. First of all, in the complete temperature range under study, χ spin and χ ang have opposite signs and χ emerges as a result of a large cancellation between the two terms. As we argued above, the spin part will necessarily turn positive for higher temperatures, eventually approaching 3/2 times the full susceptibility. Consequently, χ ang will turn negative and approach −1/2 · χ. It is intriguing to observe that in the strongly interacting regime the two contributions have opposite signs than in the usual free fermion picture according to Pauli and Landau: it is the Landau term that drives the paramagnetic response of the QCD vacuum up to temperatures T 200 MeV, while the Pauli term reduces the susceptibility in this region. This unusual behavior becomes possible due to the strong interaction, which confines quarks into composite hadrons and thereby fixes the relative orientation of their spins, i.e. their magnetic moments. In particular, in charged pions one of the constituent quarks is bound to anti-align its magnetic moment with the background field in order to maintain zero total spin. Similar effects arise for certain baryons as well. Beyond this qualitative argument, it is difficult to anticipate the outcome of this competition between the strong and the electromagnetic forces. Our quantitative results reveal a peculiar interplay between confinement and spin physics.
To further our understanding, in principle χ can also be decomposed into χ f for the quark flavors f and a gluonic contribution χ g . Subtracting this χ g from χ ang will isolate the total quark orbital angular momentum contribution f (χ f − χ spin f ), in analogy to spin decompositions [11] in deep inelastic scattering that are based on the Belinfante-Rosenfeld energy-momentum tensor, in this case of the transverse spin. The unrenormalized qualitative results of Ref. [27] indicate that χ g ∼ χ/3 at small temperatures. It may be interesting to address this quantitatively in the future.

Summary
In this paper we determined the magnetic susceptibility χ of the thermal QCD medium via a method introduced originally for T = 0 [29], which circumvents the flux quantization problem and allows us to express χ in terms of B = 0 measurements. This considerably reduces the measurement costs as well as systematic uncertainties compared to previous approaches. The susceptibility is extrapolated to the continuum limit for a broad range of temperatures, making contact to the Hadron Resonance Gas (HRG) model at low T as well as to perturbation theory at high T . In the confined phase we find evidence for a diamagnetic behavior (χ < 0), while for T 150 MeV we observe paramagnetism (χ > 0). Our continuum extrapolations are based on four lattice spacings and are guided by a generalized HRG model taking into account taste splitting (see App. A). A careful continuum limit is found to be essential to observe diamagnetism at low T since this is due to light pions -we argue that this behavior was missed in previous investigations because of large lattice artifacts.
The susceptibility is decomposed into spin-(χ spin ) and orbital angular momentum-related (χ ang ) contributions based on our previous study [26]. The spin term is shown to be given in terms of the mass-dependence of the ψ σ 12 ψ fermion bilinear in the presence of a small magnetic field, see Eq. (2.12) and App. B. Besides its role in this decomposition, the tensor bilinear is related to the normalization f ⊥ γ of the photon distribution amplitude, relevant for a range of phenomenological applications. We update our previous determination [26] of the corresponding tensor coefficient in the chiral limit at T = 0, by performing the multiplicative renormalization of ψ σ 12 ψ non-perturbatively on the lattice. We obtain the value f ⊥ γ = −45.4(1.5) MeV for massless quarks, in the MS scheme at a QCD renormalization scale of 2 GeV. The values of the tensor coefficient at the physical light and strange quark masses and at different renormalization scales are given in Eqs. (4.6)-(4.11).
At finite temperatures we performed the continuum extrapolation of χ spin and also determined the orbital angular momentum-related susceptibility χ ang . In the absence of color interactions, the two contributions exhibit the constant ratio χ spin : χ ang = 3 : (−1) as is well known since the analysis of the free electron gas by Pauli [38] and Landau [39]. Around the transition temperature, in full QCD this ratio is instead found to be close to (−1) : (1.03), resulting in a large cancellation between the two contributions, thereby substantially reducing the total susceptibility. As the temperature grows the susceptibilities approach their free-case counterparts, which are discussed in detail in App. C. Still, it is stunning to observe that in the strongly coupled QCD medium χ spin and χ ang have signs that are opposite to the naive expectations.
Considering our results at high temperature, it is interesting to make a comparison to a classical ideal system. In such a setting the Bohr-van Leeuwen theorem [67,68] (for a recent review, see Ref. [69]) holds: the total magnetization vanishes, since the magnetic field does not transfer any work to the electric currents in the system. Apparently, the QCD medium does not become classical in this sense for T → ∞, even if the O(B 2 ) terms of the free energy density that we have discussed in this paper are small compared to the dominant O(T 4 ) contributions in that limit. The non-classicality has two different origins. First, quark spins are of quantum nature and can induce a magnetization by aligning with the magnetic field. Second, both χ spin and χ ang diverge as log T for high temperatures. This behavior stems from the renormalization properties of the bare susceptibilities: quantum effects give rise to a logarithmic divergence ∝ log 1/a in the cut-off. In turn, the same behavior shows up in the renormalized susceptibilities if they are probed by another large dimensionful scale, the temperature. Note that a similar connection exists between the logarithmic divergence and the behavior of the renormalized free energy in the B → ∞ limit [33].

A The HRG model and lattice discretization errors
At low temperatures the staggered action suffers from enhanced lattice artifacts due to taste splitting.
Here we attempt to incorporate the effects of this splitting into the HRG model. The magnetic susceptibility was calculated in a standard HRG model in Ref. [28]. Following Ref. [70] we replace the contribution of pions in the model by a sum over each taste, weighted by the corresponding degeneracies. The masses of the individual tastes and their parametrization in the range of our lattice spacings are taken from Ref. [41]. Since pions are dominant for the susceptibility, the taste splitting for other mesonic and baryonic states is ignored for simplicity (although the splitting for η mesons might also lead to light mesonic states, see, e.g., Ref. [71]). The list of hadrons taken into account can be found in Ref. [72].
In Fig. 8 we show the renormalized magnetic susceptibility evaluated at T = 120 MeV as a function of the lattice spacing. The spacings for our four ensembles N t = 6, 8, 10 and 12 at this temperature are highlighted in the plot. This reveals slow convergence towards the continuum limit, which can best be understood by analyzing the mass-dependence of the pionic contribution χ π to the susceptibility, which takes the form [28] where Θ 3 is an elliptic Θ-function. This can be derived by comparing to the analogous expression for fermions, calculated below in Eq. (C.8). The bosonic Matsubara frequencies give rise to the different first argument in the elliptic function. The prefactor in this case is the scalar QED βfunction coefficient for one complex scalar field β scalar 1 = 1/(48π 2 ). The pionic susceptibility diverges logarithmically in the chiral limit (this can be shown similarly to the calculation below in App. C.5), explaining its pronounced dependence on m π . In turn, nonzero lattice spacings enhance the masses of most pion tastes, thus, reducing the magnitude of χ π . Based on the HRG predictions for χ(a, T ) we consider the difference between a simple O(a 2 ) fit taking into account only N t ≤ 12 lattices and the true continuum limit. This difference is included as a lower systematic error of our lattice determination of χ(T ) at low temperatures, see Fig. 2.

B Separation into quark spin and other angular momentum contributions
Here we derive the relation between the spin contribution to the susceptibility and the tensor bilinear, as shown in Eqs. (2.8)-(2.12) of the main text. It is instructive to begin with the first derivative of the free energy density, where we used the cyclicity of the trace (even though / D f and ∂ / D f /∂B do not commute, we can symmetrize the expression in the two operators under the trace). Now we use the relation and the identities where σ 12 is the relevant component of the relativistic spin operator defined in Eq. (2.7) and L 12 is a generalized angular momentum operator, which depends on the electromagnetic as well as the SU (3) gauge.

Using Eqs. (B.2) and (B.3), we can rewrite Eq. (B.1) as
(B.4) Thus, in the language of Eq. (3.1), we need the difference of two terms: one with valence quark mass m val f = m f and one with m val f → 0. The sea quark mass is kept fixed in both cases: m sea f = m f . We remark that the vanishing valence quark mass needs to be defined as a limit in finite volumes (see below). Also note that the fermion bilinears are defined to include the volume factor T /V . Differentiating Eq. (B.4) once more with respect to B at B = 0 and dividing by e 2 , we recover the bare magnetic susceptibility (2.1) on the left hand side, The slope of the tensor bilinear ψ f σ 12 ψ f for small values of B gives the tensor coefficient τ f b as defined in Eq. (2.8). After subtracting its value at T = 0 and multiplying by the relevant QCD renormalization factors, this term gives the spin contribution to the susceptibility χ spin , as we wrote in the main text, Eq. (2.12). In turn, the magnetic field-dependence of the bilinear involving the generalized angular momentum operator L 12 is related to χ ang . The latter term cannot be implemented straightforwardly due to its gauge-dependence and magnetic flux quantization. In Ref. [26] we already discussed the separation of the magnetic susceptibility into quark spin-and other angular momentum-related contributions. There, the m val f = 0 term of Eq. (B.4) was argued not to contribute -indeed, in a finite volume the massless limit of fermion bilinears always vanishes. However, in the thermodynamic limit this is not the case if chiral symmetry is broken spontaneously. To elucidate this point in more detail, let us rewrite the trace in Eq. (B.4) using the eigenmodes of the Dirac operator, so that, exploiting chiral symmetry {γ 5 , / D f } = 0, where ρ f (λ; m sea f ) is the spectral density of / D f in the infinite volume, determined in an ensemble generated with sea quark masses m sea f . Towards the valence chiral limit the kernel becomes proportional to the δ-distribution, so that we have a Banks-Casher-type [73] relation, (B.8) On the one hand, this limit is zero if chiral symmetry is intact and the spectral density vanishes at the origin. On the other hand, a nonzero chiral condensate ρ f (0; m sea f ) , together with the polarization σ 12 χ f 0 = χ f 0 of the low modes will turn the chiral limit of the tensor bilinear nonzero. Our lattice results reveal a nonzero value for ψ f σ 12 ψ f in the full chiral limit at low temperatures, see Fig. 4. Clearly, the fermion bilinear remains nonzero also if only m val f is sent to zero. This is in accordance with the recent findings of Ref. [74] about the Dirac spectrum at B > 0, where the low modes were indeed found to exhibit almost perfect spin-polarization.

C Susceptibilities in the free case
Here we consider the free case (i.e. we set the color charges of quarks to zero) to exemplify the most important relations of the main text. These include the proportionality between the tensor bilinear and the spin contribution to the susceptibility, the ultraviolet divergences of the susceptibilities at zero temperature as well as the high-temperature behavior of the renormalized susceptibilities. These calculations include our previous results [26,28], which we also show here for completeness.
Below we will extensively use Schwinger's proper time formulation [32]. This is based on the Mellin transform and its inverse which are valid for Re z > 0, c > 0 and E > 0. Moreover, taking the derivative of Eq. (C.1) with respect to z at z = 0 gives the standard ζ-function regularization result [75],

C.1 Magnetic susceptibility
We consider one quark flavor ψ with electric charge q and mass m in a volume V = L 3 at temperature T , exposed to a background magnetic field B. For convenience we assume that the magnetic field is oriented in the x 3 direction and qB > 0. The free energy density in this setting reads (see, e.g., Ref. [76]): where ω n = (2n + 1)πT is the n-th fermionic Matsubara frequency. Moreover, p, s and k are the momentum, spin and angular momentum in the direction of the magnetic field, N c = 3 is the number of colors and the energies are given by the Landau levels, Rewriting the logarithm using Eq. (C.3), the integral over p becomes Gaussian and can be solved. Furthermore, the sums over n, k and s are where Θ 3 is an elliptic function. Inserting these in Eq. (C.4) and performing the derivative with respect to z, we obtain Taking the second derivative with respect to eB to obtain the bare magnetic susceptibility (2.1) results in

C.2 Ultraviolet divergences and QED renormalization
To determine the ultraviolet structure of the magnetic susceptibility, we consider Eq. (C.8) at zero temperature. For T = 0 the elliptic function Θ 3 approaches unity. The resulting expression needs to be regularized, for example by setting an ultraviolet cut-off 1/Λ 2 as the lower limit of the proper time integral. Performing the integral and expanding for large Λ we obtain, where γ E is the Euler-Mascheroni constant. Thus, the coefficient of the logarithmic divergence indeed equals the lowest-order QED β-function coefficient β 1 (for one quark flavor with electric charge q), demonstrating the validity of Eq. (2.5). In fact, this relation continues to hold in full QCD as well, owing to the fact that towards the continuum limit QCD corrections to β 1 at the scale 1/a approach zero due to asymptotic freedom (see Eq. (4.1)). We note moreover that in the proper time formulation the renormalization scale is set by the mass -in fact µ QED = m e γ E /2 for our choice of the regulator Λ -explaining the appearance of m in the argument of the logarithm in Eq. (C.9). The additive renormalization can be performed by subtracting χ b (T = 0) from Eq. (C.8): As we mentioned after Eq. (2.6), this corresponds to the choice of a physical, albeit scheme-dependent, QED renormalization scale.

C.3 Spin contribution
The contribution of orbital angular momentum to the total susceptibility can be calculated by simply replacing the fermion with two ghost particles (spin-zero but antiperiodic in Euclidean time) in the above calculation. This removes the −2sqB from the energies (C.5) and excludes the spin sum in the free energy density (C.4). Consequently, the magnetic field-dependent part in Eq. (C.7) changes as coth(qBt) → 1/ sinh(qBt). This merely changes the second derivative of the free energy density at B = 0 by a factor −1/2. Thus, for the renormalized susceptibility we arrive at which also implies confirming the 3 : (−1) ratio of the two contributions to the total susceptibility. We mention that a similar argument has been used in perturbative QCD (with chromomagnetic background fields) to relate asymptotic freedom to spin effects [77].

C.4 Tensor bilinear
For the tensor bilinear we begin with the result of the fermionic path integral, where we used chiral symmetry {γ 5 , / D} = 0. The trace is represented using the eigenbasis of − / D 2 , giving the eigenvalues ω 2 n + E 2 p,s,k . Since [ / D 2 , σ 12 ] = 0, the spin operator can also be diagonalized in this basis and its eigenvalues are minus two times the spin: σ 12 → −2s. Taking into account the 2N c · (qBL 2 )/(2π)-fold degeneracy of the eigenvalues, we obtain (C.14) In the sum the contributions {k, s = 1/2} and {k + 1, s = −1/2} cancel, leaving only the unpaired lowest Landau level {k = 0, s = 1/2}. Hence we get Note that, unlike in full QCD, here the tensor bilinear is exactly linear in the magnetic field. Thus, the tensor coefficient τ b of Eq. (2.8) is obtained by simply dividing Eq. (C.15) by qB. Using Eq. (C.1) with E = ω 2 n + p 2 + m 2 , performing the Gaussian integral over p and the Matsubara sum (C.6) over ω n , we arrive at A comparison to Eq. (C.8) reveals that this quantity contains the same logarithmic divergence as χ b , just with a different coefficient. Using a cut-off regulator as in Eq. (C.9), we obtain at T = 0, confirming Eq. (2.9). The same considerations regarding QCD corrections to the coefficient and the renormalization scale µ QED apply as in Sec. C.2 for χ b .
We can compare this with Eqs. (C.10) and (C.12) to conclude that confirming the relation (2.12) and Eq. (C.12). Notice that τ vanishes for m → 0, so the subtraction of the massless limit is irrelevant in the free case (but it is relevant for the interacting system with spontaneous chiral symmetry breaking, see App. B).

C.5 High-temperature expansion
The temperature-dependent part of the free energy density (C.4) can be simplified using the wellknown trick [52] of differentiating and subsequently integrating the integrand with respect to E p,s,k .
The result is The energy levels are given in Eq. (C.5) above. To obtain the high-temperature expansion in a closed form, we need to replace the logarithm by its series expansion This approach was used, e.g., in Ref. [78] for scalars at nonzero chemical potential. Inserting the expansion (C.20) into (C. 19) and rewriting the exponentials using Eq. (C.1) results in Inserting the Mellin transform (C.1) for E −z p,s,k renders the integral over p Gaussian. We can reuse the angular momentum and spin-sums from Eq. (C.6), giving (C.22) Differentiating the above expression twice with respect to eB at B = 0 gives (minus) the renormalized magnetic susceptibility χ. The integral over t can be solved via the Mellin transform (C.1) and gives a Γ-function, while the sum over l results in a ζ-function: Using the duplication formula [79] for the ratio of Γ-functions, we arrive at For the validity of the Mellin transforms we needed to assume c > 0 (as well as m > 0). The final integral over z can be solved using Cauchy's theorem, closing the integral towards the left and calculating the residue at the poles. There is a double pole at z = −1 and simple poles at z = −3, −5, . . .. These z-values set the powers of T that appear in the high-temperature expansion. Keeping the leading terms (i.e. z = −1 and z = −3), we finally obtain, with the numerical constants C ≈ 0.21314 and D ≈ 0.59836. Notice that the coefficient of the leading logarithmic term is equal to β 1 (for one flavor with electric charge q), confirming Eq. (4.3). As we have seen below Eq. (C.9), in the proper time formulation the renormalization scale is intrinsically set by the quark mass, µ QED = m e γ E /2 . We may express the square bracket in the leading term as log(γ T 2 /µ 2 QED ) with γ = π 2 e −γ E . Clearly, γ = O(1) depends on the definition of the regulator. The general form is again expected to hold in full QCD [28]: in this case QCD corrections at scales T µ QED are small due to asymptotic freedom.

D Multiplicative QCD renormalization
Since lattice perturbation theory is slowly convergent and high-loop results are unavailable, we first match the local lattice QCD operators of interest non-perturbatively to the regulator independent RI'-MOM scheme [30,31] and subsequently translate the result at three-loop order [36] to the MS scheme.
The quark bilinear operators are renormalized by computing the corresponding amputated flavor non-singlet vertex functions for different momenta on Landau gauge-fixed ensembles. We wish to renormalize light-and strange-quark bilinears, which can be written as linear combinations of the diagonal SU(3) flavor-octet and -singlet currents. In continuum schemes, with the exception of the axial current that we do not discuss here, the renormalization of flavor singlet and non-singlet operators of dimension three is the same. This also appears to hold for the staggered action [37,80]. We remark that, instead of extrapolating to the N f = 3 massless case, we use physical quark masses, which may be problematic, in particular regarding the strange quark mass. However, in Ref. [81] it was demonstrated that the effect of the mass-dependence is tiny for the perturbative momentum transfers that we are interested in. Moreover, the difference is expected to vanish after a continuum limit extrapolation of a renormalized matrix element is carried out because our quark masses are tuned to a line of constant physics.
Since the spin degrees of freedom are spread over hypercubes for staggered fermions, the determination of the vertex function in momentum space is more challenging than for Wilson fermions. We follow the approach described in Ref. [82]: the taste and spin degrees of freedom are reconstructed from different momentum combinations. The quark propagator for a given momentum, as any vertex function, will be a matrix of size 16 × 16, after averaging over the color degrees of freedom. Our choice of the scalar and tensor currents, where, in the latter case, we employ a two-link operator, is detailed in Ref. [26].  The error of the final renormalization constants is dominated by systematics. On the one hand, the conversion factors from the RI'-MOM to the MS-scheme are only known up to a fixed order in perturbation theory (three loops in our case). Hence high momenta are preferable. On the other hand, at momentum scales close to the lattice cut-off the intermediate matching to the RI'-MOM scheme will significantly be affected by lattice artifacts. Therefore, we are restricted to a "window" of intermediate momentum values. We employ combinations along symmetric lattice directions, where the lattice corrections are smallest. Another complication is that due to the choice of the staggered action, the maximum momentum scale that can be achieved on a four-dimensional lattice is π/a, rather than 2π/a. As a compromise, on the finest three lattices we interpolate the RI'-MOM result to the fixed scale µ QCD = 2 GeV. Subsequently, this is perturbatively converted to the MS-scheme. We estimate the uncertainty by adding the difference between the scheme conversion at two-and at three-loop order and the (statistical and systematic) interpolation uncertainty in quadrature. The latter contribution is negligible. At the coarsest two lattice spacings, µ QCD = 2 GeV is too close to the cut-off scale to obtain reliable results. Therefore, at a ≈ 0.28 fm and at a ≈ 0.22 fm, we convert the RI'-MOM results at µ QCD = 1.1 GeV and at µ QCD = 1.5 GeV, respectively, to the MS-scheme and evolve the result to µ QCD = 2 GeV. We replicate the same procedure at a ≈ 0.15 fm and add the difference that we obtain at this lattice spacing between the matching at µ QCD = 2 GeV and the matching at these lower scales in quadrature to the systematic error at µ QCD = 1.1 GeV and µ QCD = 1.5 GeV.
The results are listed in Table 1 and shown in Fig. 9. We also include the lattice perturbative theory one-loop expectations [26] in the figure. The comparatively larger value of Z T results in a larger modulus of the renormalized tensor coefficient.

E Susceptibilities via current-current correlators
In order to circumvent the issue of flux quantization, we consider a background field that possesses nonzero momentum p 1 in the x 1 direction. The constant field setup will be approached via the p 1 → 0 limit. In this section we are working in the infinite volume, where p 1 is a continuous variable and this limit exists. We will discuss finite volume corrections below. This approach has been described in detail in Ref. [29] for χ b in momentum space. Here we repeat the argument in coordinate space and also generalize it for τ f b .

E.1 Magnetic susceptibility from correlators
We consider an oscillatory magnetic field and the corresponding Landau-gauge vector potential, The latter couples to i · e times the current (3.4) in the action density. We can define the associated susceptibility just like in Eq. (2.1), where each derivative brought down an integral over the current j 2 times the coordinate-dependence of A 2 and we used j 2 = 0. Changing the integration variable from z to x = z − y and exploiting the translational invariance of the current-current correlator, the integrals over y 2 , y 3 and y 4 can be carried out, where the projected correlator G(x 1 ), defined in Eq. (3.5), appears. In the p 1 → 0 limit, B(x 1 ) becomes homogeneous and χ p 1 ,cos b equals the ordinary susceptibility χ b . For reasons that will become clear in a moment, let us consider a different background field, for which the associated oscillatory susceptibility, similarly to Eq. (E.3), reads In this case the p 1 → 0 limit does not reproduce χ b . Instead, A 2 (x 1 ) becomes homogeneous: it acts as if we had introduced a constant imaginary 'chemical potential' in the x 2 direction, with magnitude µ 2 = −eB/p 1 . Therefore the oscillatory susceptibility becomes proportional to the leading response to this spatial chemical potential, This detour was necessary to simplify the p 1 → 0 limit of the oscillatory susceptibilities. Specifically, let us examine the following combination: This approaches χ b for p 1 → 0. Using the trigonometric identity for the cosine of the difference of angles in the numerator of the kernel reveals that the integrand is independent of y 1 . Integrating over this variable cancels the prefactor 1/L, resulting in where we finally performed the p 1 → 0 limit. In finite volumes (x 1 ∈ [0, L]) the momentum variable p 1 is discrete, so that the p 1 → 0 limit does not exist. Nevertheless, we can safely employ the final formula Eq. (E.8) directly in finite volumes, as long as the linear size L is much larger than the characteristic length governing the exponential decay of G(x 1 ). Symmetrizing Eq. (E.8) to comply with periodic boundary conditions and the symmetry G(x 1 ) = G(L − x 1 ), we arrive at Eq. (3.6) of the main text.
We investigated finite volume effects by comparing the bare susceptibilities obtained on 16 3 × 6 and 24 3 × 6 ensembles. The results are found to agree within statistical errors. In Fig. 10 we show the vacuum polarization function Π(p 2 ) of Eq. (3.7) for spatial momenta p = (p 1 , 0, 0, 0), as obtained via Eq. (E.8) before taking the limit p 1 → 0. Here we considered the results at T ≈ 124 MeV but the overall picture is very similar for the other temperatures. Figure 10. The volume dependence of the determination of the vacuum polarization function at T ≈ 124 MeV. We compare the 24 3 × 6 (red) and 16 3 × 6 (blue) results. The bare magnetic susceptibility can be read off from the intersect Π(0).

E.2 Tensor coefficient from correlators
We generalize the above derivation for τ f b , which can be written as Again we consider oscillatory magnetic fields of the types (E.1) and (E.4). These give rise to modulated tensor bilinears of the formsψ f σ 12 ψ f (x) cos(p 1 x 1 ) andψ f σ 12 ψ f (x) sin(p 1 x 1 ), respectively, which enter the corresponding oscillatory tensor coefficients τ p vanishes in that limit. Thus we need to consider the sum of the two coefficients. Employing the trigonometric identity for the sine of the difference of angles and carrying out the integral over y 1 gives In finite volumes we carry out the same symmetrization as in Eq. (3.6), this time taking into account that H f (x 1 ) = −H f (L − x 1 ) to finally arrive at Eq. (3.8) of the main text. This method to calculate τ f b is compared to the results for ψ f σ 12 ψ f measured at B > 0 on 24 3 × 6 lattices at T = 113 MeV in Fig. 11. For the light quarks we obtain consistent results, however, Figure 11. Comparison of different methods to calculate τ f b for all three flavors. Simulations at nonzero (quantized) values of the magnetic field (points) are compared with a direct determination of the slope at B = 0 (colored bands).
for τ sb the correlator tends to give values that slightly differ from the slope of a linear fit to the lowest few available points. Since lattice artifacts and finite volume effects might be different in the two cases, such slight differences are not unexpected.
In addition, we find that a linear fit to results from simulations at B > 0 has smaller uncertainties than extracting the slope at B = 0 using the correlator method. In the main text we therefore use our earlier results for ψ f σ 12 ψ f from Ref. [26].
We note that the tensor-vector correlators at nonzero spatial momenta might also be useful for extracting further features of the photon distribution amplitude.

F Parametrization of the equation of state
Up to O(B 2 ), the magnetic field-dependence of the complete EoS can be calculated from the magnetic susceptibility χ(T ). Here we provide a parametrization for this observable and also collect the relevant thermodynamical relations, which were also summarized in Ref. [28].
First of all, we remind the reader that in the presence of a background magnetic field, the different components of the pressure -defined by considering an infinitesimal compression of the system in the respective direction -might become anisotropic [27]. In particular, one should distinguish between the Φ-scheme, where the flux of the magnetic field is kept constant during the compression (superscript (Φ) below), and the B-scheme, where the magnetic field strength is kept constant (superscript (B)). On the one hand, the B-scheme pressure is isotropic and equals the negative of the free energy density in the thermodynamic limit, p Thus, to calculate the complete EoS including B 0 and B 2 effects, altogether it suffices to parameterize I(T, 0) and χ(T ). For the latter we consider a parametrization of the continuum extrapolated lattice results that smoothly approach the HRG model prediction (see Fig. 2) at low and the perturbation theory formula (4.3)) at high temperatures. We found the following parametric form to be sufficient for this, χ(T ) = exp(−h 3 /t) · 1 + g 0 /t + g 1 /t 2 + g 2 /t 3 1 + g 3 /t + g 4 /t 2 + g 5 /t 3 · 2β 1 log Eq. (F.7) incorporates the non-perturbative temperature-dependence predicted by the HRG model (see App. A) at low T and the logarithmic rise at high temperatures. The β 1 coefficient is fixed to its perturbative value (2.4), while the scale q 0 inside the logarithm is allowed to be a free parameter. The rational function involving the g i parameters interpolates between the two limiting behaviors. The so-obtained parametrization is shown in the left panel of Fig. 3 in the main text.
For the interaction measure we take the parametrization of Ref. [83], The parameters of both functions are included in Table. 2. The two parametrizations, together with the implementations of the formulae (F.1)-(F.6) are included in the Python script param_EoS.py that is submitted to arXiv.org together with this manuscript. This parametrization is valid for low magnetic fields. To be more quantitative, we compare our O(B 2 ) truncated results for the longitudinal pressure to the complete magnetic field-dependence from Ref. [28] for T 180 MeV. We find agreement within errors in the range B/(πT ) 2 1. This upper limit is hard-coded in the Python script as well. One final remark about the parametrization is in order. All truncated thermodynamic observables approach zero for T → 0, such that a normalization by the corresponding powers of the temperature (i.e. p 3 /T 4 , s/T 3 and so on) produces sensible plots. This is not the case if O(B 4 ) terms are also included: at this order vacuum contributions arise and the equation of state depends on B already at T = 0, rendering a normalization like p 3 /T 4 ill-defined in the T → 0 limit [28].