On electric fields in hot QCD: perturbation theory

We investigate the response of a hot gas of quarks to external electric fields via leading-order perturbation theory. In particular, we discuss how equilibrium is maintained in the presence of the electric field and calculate the electric susceptibility, providing its high-temperature expansion for arbitrary quark mass. Furthermore, we point out that there is a mismatch between this, direct determination of the susceptibility at zero field and the weak-field expansion of the effective action at nonzero electric fields, as obtained using Schwinger's exact propagator. We discuss the origin of this mismatch and elaborate on the generalization of our results to full QCD in electric fields.


Introduction
The quantum vacuum is not empty space but a dynamical state containing virtual particles that can be viewed as a medium. This medium can be most efficiently probed by background electromagnetic fields, both for the case of the strong interactions -virtual quarks and gluons described by Quantum Chromodynamics (QCD), as well as for electromagnetism -virtual electrons and photons described by Quantum Electrodynamics (QED).
For strongly interacting systems, such background electromagnetic fields are known to play a relevant role for a range of physical settings including magnetized neutron stars, the evolution of the early Universe and off-central heavy-ion collisions. In the latter case, the phenomenological impact of magnetic fields B has been recognized primarily in the context of the chiral magnetic effect [1] and collective motion. These collisions feature magnetic fields that by far exceed the QCD scale [2,3], even though their life-time in the plasma is the subject of active debate [3][4][5]. Recently, it was recognized that heavy-ion collisions might also generate strong electric fields E on an event-by-event basis, with components comparable to the magnetic ones in magnitude [6,7]. This has been in the focus of recent interest also due to the isobaric experimental program at RHIC [8].
In turn, the response of the QED vacuum to intense background fields is by now routinely explored by high-intensity lasers [9,10]. These new-generation laser facilities (e.g. Ref. [11]) operate at unprecedented intensities, probing the structure of the quantum vacuum and aiming at recreating astrophysical environments in the laboratory [12]. Current efforts also include collision experiments involving high-energy electron beams and laser pulses, as in the upcoming LUXE experiment [13]. For a recent summary, see e.g. Ref. [14].
Much of our knowledge about the impact of magnetic fields on QCD matter comes from firstprinciples lattice simulations regarding the phase diagram [15][16][17][18][19] and the equation of state [20][21][22], as well as from a wide range of different model and effective theory approaches (for a recent review, see Ref. [23]). In contrast, there are only a few results about the role of electric fields, mainly due to the fact that the QCD action at E = 0 becomes complex, hindering standard simulation algorithms [24][25][26][27]. Owing to asymptotic freedom in QCD, ordinary perturbation theory is expected to work as a high-temperature alternative of the lattice approach. Nevertheless, convergence problems restrict the region of applicability to very high temperatures [28][29][30], unless mended by resumming the perturbative series using, e.g., hard thermal loop techniques [31], bridging the gap towards the lattice approach [32,33]. In stark contrast, the effect of weak magnetic fields is found to be captured very well already at the one-loop order of ordinary perturbation theory down to T ≈ 300 MeV [21,22]. A perturbative treatment at nonzero electric fields is therefore worth pursuing and in this paper we will focus on the one-loop order. In turn, the analogous perturbative calculation is applicable for hot QED as well.
There are two main approaches to study this problem. First, one may use the exact fermion propagator at nonzero E due to Schwinger [34] and construct the free energy density f (E), followed by a weak-field expansion to find the leading coefficient, the electric susceptibility ξ. This method has been employed both with real-time propagators [35][36][37] as well as within the imaginary time formalism [38,39] at arbitrary temperatures. Second, one may consider the expansion of f in the presence of arbitrary external photon fields to express ξ via the photon vacuum polarization diagram directly at E = 0. This approach goes back to Weldon [40], who calculated the leading O(T 2 ) temperature-dependence for massless fermions. The above two ways of treating background fields in hot plasmas have both been instrumental for a whole set of subsequent calculations. On the one hand, Schwinger's approach has been used for constructing Euler-Heisenberg-type effective actions [41] in the context of laser physics, for QED thermodynamics as well as QCD effective models [42]. On the other hand, Weldon's approach, generalized to color interactions, has formed the basis for hard thermal loop perturbation theory [43,44] and was employed for various physical applications within quark-gluon plasma physics, see e.g. Ref. [45].
In this paper we follow Weldon's method and calculate the electric susceptibility for arbitrary fermion masses and present its high-temperature expansion for the first time within this approach. We discuss the importance and effect of the equilibrium charge distribution that forms under the impact of the electric field. Furthermore, we point out that the so defined susceptibility disagrees with the one obtained via Schwinger's approach, revealing the significance of an infrared regularization for this observable and the non-commutativity of the involved limits (zero electric field and infinite volume). Finally, we discuss how our results can be generalized to full QCD.

Susceptibilities, renormalization and equilibrium
In the presence of weak background electromagnetic fields, QCD matter in thermal equilibrium may be treated as a linear medium. The total free energy density, considered as a function of the thermodynamic variables E and B, reads f tot = B 2 /(2µ) − E 2 /2, where µ is the magnetic permeability and the electric permittivity [46]. These parameters enter the constitutive relations B = µH and D = E. In the following we refer to B and E simply as magnetic and electric fields.
The deviation from the vacuum values ( = µ = 1 in natural units) can be measured by the electric and magnetic susceptibilities as = 1 + e 2 ξ and µ = (1 − e 2 χ) −1 . The factor e 2 -the square of the elementary charge -is included here for later convenience. Below we exclude the energy f γ = (B 2 − E 2 )/2 of the fields themselves and merely consider the matter contribution f = f tot − f γ . In terms of the partition function Z of the system, it is given by f = −T log Z/V , where T is the temperature and V the spatial volume. The susceptibilities are then defined as (2.1) Note that at T = 0, Lorentz invariance ensures that χ b = −ξ b . In Eq. (2.1) the subscript b indicates that both susceptibilities contain ultraviolet divergent terms that must be subtracted via additive renormalization. These divergences originate from the multiplicative divergence in the electric charge e [34,41], with the renormalization constant which up to O(e 2 r ) reads where Λ is an UV cutoff and σ the renormalization scale. Note that the QED Ward identity ensures that e r B r = eB and e r E r = eE are renormalization group invariant combinations. To see how this is related to additive divergences in the bare susceptibilities, consider the total free energy density, which is a finite, physical quantity, For this to be finite, Eqs. (2.2) and (2.3) imply that Being temperature-independent, the additive divergence cancels in which sets ξ = χ = 0 at zero temperature after renormalization, consistent with the requirement that the electric permittivity and magnetic permeability µ of the vacuum both be unity. For a recent summary on the role of electric charge renormalization for susceptibilities and on the choice of the renormalization scale σ, see Ref. [22]. Above we assumed E and B to be homogeneous. While for the magnetic field this assumption is consistent with equilibrium, a constant electric field accelerates charged particles indefinitely and inevitably leads the system out of equilibrium. In other words, homogeneous electric fields necessarily involve some type of infrared regularization, for example a finite volume with boundaries, where charges can accumulate. Alternatively, one can consider spatially oscillating electric fields and their infinite wavelength limit. For such a (time-independent) electric field, the system will find the equilibrium configuration, where electric and diffusion forces cancel each other and no currents flow 1 , illustrated in Fig. 1. In this setup, the impact of the electric field is therefore to displace electric charges and to polarize the medium. Rotational invariance and the homogeneity of the E = 0 medium ensures that for weak fields the polarization oscillates with the same wavelength as the electric field does. We note that this interpretation also holds for time-dependent fields, as long as the corresponding time-like frequency approaches zero before the inverse wavelength does. 2 Without loss of generality, the electric field is assumed to point in the x 1 direction and from now on, E denotes this non-vanishing component of the field unless otherwise noted. The electric field couples to the quark charges q f via the vector potential A ν , for which we choose the static gauge A 0 (x 1 ) and A = 0. A constant electric field is then realized by A 0 = Ex 1 , while an example oscillating field with a well defined equilibrium is E(x 1 ) = E · cos(k 1 x 1 ), generated by A 0 (x 1 ) = E · sin(k 1 x 1 )/k 1 with spatial momentum k 1 , illustrated in Fig 1. For completeness, we will also include results for the magnetic susceptibility χ. For the latter, we choose a gauge potential A 2 (x 1 ), so that B points in the x 3 direction, however we always consider electric and magnetic fields separately. In the following we calculate the susceptibilities for a single fermion with charge e -the result for QCD is obtained via multiplication by the number N c = 3 of colors and by f (q f /e) 2 , where f runs over the quark flavors.

Calculation of the electric susceptibility
The definition (2.1) of the electric susceptibility can be evaluated in two different ways, as already pointed out in Sec. 1. We will first obtain ξ by evaluating the derivatives of the free energy directly at E = 0 using the expansion of the partition function in terms of arbitrary but weak fields [40], and then we will compare this result to Refs. [35][36][37], where the free energy is calculated at nonzero but homogeneous E and the susceptibility is read off from its weak field expansion. Throughout we are using one-loop perturbation theory. For definiteness, the metric tensor is g µν = diag(1, −1, −1, −1) and four-vectors are denoted as for example K = (k 0 , k) with k = |k| and K 2 = k 2 0 − k 2 .

Expansion at arbitrary small fields
As proposed in Sec. 2, in the presence of the electric field the equilibrium situation involves a fixed charge density profile n(x 1 ) that varies in the x 1 direction in the system (of size L). This follows from the equilibrium requirements of a statistical system in a background field, see e.g. [48]. In this setup the temperature and the total chemical potential are constant, however the pressure and the particle density are inhomogeneous. We will first discuss the implications of such an equilibrium using a homogeneous background field regularized by the finite system size, rather than by oscillating with a well defined wavelength. The generalization for the oscillating fields will follow, however we found it more didactic to present this matter in two steps.

Impact of the charge distribution
For weak fields, the gradient of n(x 1 ) is small so that one can imagine the system as a collection of smaller (but still large) subsystems with each subsystem having a given constant density. This is then described by a canonical thermodynamic potential f that depends on the local density and is related to the grand canonical potential Ω via an trading the dependence on a local chemical potentialμ(x 1 ) for the dependence on n(x 1 ). The fixed density is such that the resulting diffusion force balances the electric field: In other words, the sum of the local chemical potential due to the accumulation of charges and the external potential is constant [48]. A constant shift inμ(x 1 ) would merely change the total charge in the system -we set this to zero assuming a neutral medium. Thus, the canonical potential f , relevant for the equilibrium situation, depends on E in two ways: explicitly as well as implicitly via the charge distribution. This implicit dependence reads Altogether, the susceptibility (2.1) as the second total derivative with respect to E, becomes where we used that Ω is an even function both of µ and of E. Notice that the first term is the curvature of Ωwith respect to the electric field at µ = 0, i.e. without the equilibrated charge distribution. In turn, the second term is related to the free energy of the displaced charges themselves. For homogeneous fieldsμ(x 1 ) = −eEx 1 and this term is quadratically divergent in the system size, ∝ L 2 . We will find that the first term contains the same divergence with an opposite sign so that ξ b is altogether infrared finite.

Susceptibilities for oscillatory fields
Now we generalize (3.4) for the case of oscillating fields. In the spirit of the end of Sec. 2, we consider the general linear response of a translationally and rotationally symmetric system thus redefining the electric susceptibility as the limit where the time-like frequency k 0 approaches zero first, as discussed in Sec. 2.
To find ξ b in this setup, we need Ω for an arbitrary background gauge field A α in order to express the derivatives appearing in Eq. (3.4). To this end, we write down the leading order expansion of Ω for general background gauge fields [40], where A α (K) is the momentum-space gauge field and Π µν (K) the photon vacuum polarization tensor, see Fig. 2. Based on Eq. (3.7), we can evaluate the partial derivatives in Eq. (3.4) by choosing specific gauge fields representing either the oscillating electric field E(x 1 ) = E · cos(k 1 x 1 ) or the chemical potential. Explicitly, for the ∂ 2 Ω/∂(eE) 2 term appearing in (3.4), we have to consider A α (K) = δ α0 E(K)/(ik 1 ). In turn, for ∂ 2 Ω/∂µ 2 we need A α (K) = δ α0 µ δ (4) (K), and finally (∂μ(x 1 )/∂(eE)) 2 becomes δA 0 (x 1 ) Combining all this we find where in the second term it is understood that k 0 approaches zero before k does. Since nothing depends on x 1 anymore, the averaging is trivial. Finally, we subtract the zero temperature part to renormalize as prescribed in (2.6), µν is the thermal part of the vacuum polarization diagram. Here we replaced k 2 1 by k 2 , restoring the rotational symmetry of the system, which we explicitly broke by choosing E to point in a particular direction. Moreover, we rewrote the k → 0 limit of the finite difference as half of the second derivative with respect to k. This is legitimate, as Π T =0 00 has no linear term in k around k = 0 due to parity symmetry. The first term in the square brackets of Eq. (3.9) agrees with Weldon's formula [40], while the second term stems from the local charge distribution, which was not taken into account in Ref. [40].
The analogous expression for the magnetic susceptibility does not have the second term, since a modified equilibrium charge distribution is not formed for purely magnetic backgrounds. For the gauge A 2 (x 1 ) we chose for the magnetic field, the component Π 22 shows up in the magnetic susceptibility based on (3.7). One can restore rotational symmetry by including the other spatial components (see e.g. [40]), and the magnetic susceptibility becomes Unlike the electric component of the gauge field, the magnetic ones do not acquire thermal masses and therefore Π T =0 S (0, k) ∼ k 2 in the low-momentum region. Thus we could replace division by k 2 by half of the second derivative with respect to k, just as above. Eqs. (3.9) and (3.10) are our final formulas for the susceptibilities. Now we proceed by calculating the vacuum polarization tensor at finite temperature. To calculate the thermal part Π T =0 µν (k 0 , k), we follow and extend the calculation of Weldon [40] to arbitrary fermion masses m. The general formula is

Photon vacuum polarization diagram
where the fermion propagators running in the loop are the finite temperature Feynman propagators with n F (|p 0 |) = (e β|p 0 | + 1) −1 being the Fermi-Dirac distribution. We relegate the details of the calculation to App. A; here we only summarize the results. Notice that in the limit k → 0 when evaluating ξ, the leading contribution Π T =0 00 (0, 0), i.e. the thermal photon mass, drops out and ξ b is infrared finite. The exact value of this term we merely use to check against previous results (see (A.4) and (A.5)). Nevertheless, since the thermal mass is non-vanishing, should one neglect the Edependence of the equilibrium charge distribution, the susceptibility would be infrared divergent, as we already mentioned above in Sec. 3.1.1. This infrared divergence was already observed by Weldon [40] (in the massless case), but its interpretation in terms of the equilibrium charge distribution was not discussed.
The second order term in the expansion of Π T =0 00 (0, k) around k = 0 gives then the susceptibility, which reads (see the details again in App. A) and hence (3.14) One can expand in terms of m/T to obtain a high-temperature expansion (HTE) expression where ζ (z) is the derivative of the Riemann ζ function. The leading term in (3.15) agrees with the results of Refs. [40,49] for massless fermions -in that case m is replaced by the renormalization scale under the logarithm. For the magnetic case using (A.13) from App. A we find which can be expanded in m/T to obtain (3.17)

Background field approach
In this subsection we will summarize the results described in detail in Ref. [37], with similar calculations carried out in e.g. Refs. [35,38,39]. The method relies on the calculations of Schwinger [34], where the exact fermion propagator in a constant electromagnetic background was obtained and used to evaluate the one loop effective potential in the vacuum. The extension of this method to finite temperature is carried out in [37], providing a separation of the effective potential at nonzero E (and B) into a vacuum and a temperature-dependent part. We concentrate here on the latter. The one loop effective action is then obtained as a series where f 2n (ω) is of the order 2n in the electromagnetic fields. Therefore we only need n = 1 to read off the susceptibility 3 leading to The expression for χ agrees with (3.16), demonstrating the equivalence of the expansion approach and the background field approach for the magnetic case. In contrast, the same is not true for the electric response: in ξ the term involving the derivative of the Fermi-distribution comes with the opposite sign as in the expansion approach, (3.14). The formula (3.20) leads to the high-temperature expansion differing from the result (3.15) of the expansion approach. Most prominently, the sign of the O(m 0 ) term next to the logarithmic term is the opposite here 4 . Notice that the two approaches differ by the order of two distinct limits -while in the small field expansion approach the infrared regulator (the external momentum) is sent to zero after E → 0, here the opposite is done, since no explicit long wavelength regulator had to be introduced.

Numerical comparison
Next we check the analytical one-loop results within the small field expansion by means of numerical calculation of the free fermion determinant on the lattice. To this end, we need to perform a Wick rotation to imaginary times. Using the particular gauge choices A 4 = iEx 1 and A 2 = Bx 1 , the sum of the electric and magnetic susceptibilities, Eqs. (3.9) and (3.10), becomes It is advantageous to work with this sum, because it is free of additive divergences, since at strictly zero temperature χ b (T = 0) = −ξ b (T = 0) due to Lorentz invariance. Therefore the sum of the renormalized susceptibilities can be calculated without explicit T = 0 subtractions, cf. Eq. (2.6). Going over to coordinate space, we obtain The two methods therefore lead to opposite signs for the sum χ + ξ. While it is positive for the background field approach, it is negative for the small field expansion. As well known [50], this sum of susceptibilities enters the speed of light in the medium, v 2 /c 2 = 1/( µ) = 1 − e 2 (ξ + χ), see also Ref. [39]. As a light wave with momentum K propagates, the electric field that it involves attempts to polarize the medium. The equilibrium charge distribution that we discussed in this paper is not reached in this case. Thus, one should exclude the term due to the nontrivial charge profile from ξ, leading to ξ = lim k→0 m 2 γ /k 2 with m 2 γ = Π T =0 00 /e 2 = T 2 /3, see App. A. This turns the sign of χ + ξ positive. In fact, ξ diverges in this case for k → 0, signalling that such waves cannot enter the thermal medium due to the nonzero photon mass mγ. For a more detailed discussion on light speeds in media, see also Ref. [51].  involving the mixed-representation current-current correlator combination G(x 1 ). Here we also highlighted that our numerical calculations are carried out in a finite volume of linear spatial extent L. Placing the origin of the coordinate system in the middle makes the formula (4.2) consistent with periodic boundary conditions. Truncating the integral at ±L/2 introduces exponentially small errors, because both components of the coordinate-space vacuum polarizations (i.e. current-current correlators) decay exponentially with the distance. This point is discussed in detail for example in Ref. [22].
We work with the staggered fermion discretization for the Dirac operator on high-temperature lattices N 3 s × N t with spacing a at nonzero fermion mass m. The spatial size and the temperature are L = N s a and T = 1/(N t a), respectively. We calculate G(x 1 ) using a few thousand random wall sources [52]. Since Π 44 and Π 22 correlate with each other, it is advantageous to take their difference for each random source to cancel some of the statistical noise. The continuum limit of the result is defined by keeping LT and T /m fixed and sending N s , N t → ∞. In turn, the thermodynamic limit is obtained by fixing ma and T a = 1/N t and sending only N s → ∞.
We consider a high temperature with T /m = 12.5. The continuum extrapolation on LT = 4 lattices is shown in the left panel of Fig. 3, revealing perfect agreement with the analytical prediction. Note that only the three smallest lattice spacings were used in the extrapolation, indicated by the grey dotted line on the plot. We also compare different spatial sizes and do not observe significant finite volume effects. The ratio of the continuum extrapolation and the N t = 6 result is used as an improvement factor. The right panel of Fig. 3 shows the temperature-dependence of the 24 3 × 6 data, after a multiplication by this improvement factor. The small deviations from the analytical result are explained by temperature-dependent lattice artefacts. In summary, the numerical data fully confirm our analytical findings.
Finally, we comment on lattice simulations at nonzero (imaginary) electric fields and nonzero temperature. As the analytic continuation of the setup with real fields, the equilibrium charge −∆Ω/T 4 ieE/T 2 oscillatory k=1 200 3 ×20 oscillatory k=2 200 3 ×20 oscillatory k=3 200 3 ×20 oscillatory k=1 300 3 ×20 periodic 200 3 ×20 periodic 300 3 ×20 Figure 4. The imaginary electric field dependent part of the grand canonical free energy density as a function of the amplitude of the field strength. The lines represent oscillatory electric fields with different wavelengths, while the dots mark the results for homogeneous fields. Both setups show that towards the thermodynamic limit the grand canonical free energy density develops a discontinuous behaviour at E = 0.
distribution in this case is also inhomogeneous 5 . However, such simulations are always naturally performed in the grand canonical ensemble, characterized by a chemical potential and not in the canonical one, which would fix the equilibrium charge distribution, see Sec. 3.1.1. This implies that there is a change of relevant ensembles between E = 0 and any E = 0, signalled by an abrupt discontinuity in physical observables in the infinite volume limit.
We demonstrate this effect for the grand canonical free energy density Ω = −T /V log det( / D + m) in the free case for two different setups. First, homogeneous imaginary electric fields satisfying periodic boundary conditions, where the field is quantized as eE = 2πN E T /L with N E ∈ Z. Second, oscillatory imaginary fields with profile E(x 1 ) = E cos(kx 1 ), again with periodic boundary conditions so that k = 2πN k /L and N k ∈ Z. The determinant is calculated exactly by diagonalizing the Dirac operator on two-dimensional subspaces. Fig. 4 shows the negative of ∆Ω = Ω(E) − Ω(0) for both setups on 200 3 × 20 and 300 3 × 20 lattices at T /m = 12.5. In the quantized case, the data at discrete electric field values collapse on a line at E > 0 that clearly has a different intersect as the E = 0 value. In the oscillatory case, the curves also approach a discontinuity at zero field as k decreases (either by reducing N k or by increasing L). To eliminate this discontinuity, the equilibrium construction at every nonzero E would have to be performed according to Eq. (3.1). The effects of the discontinuity will also be inherited by other observables derived from Ω.

Conclusion
We calculated the leading order electric response of a non-interacting charged thermal medium using the small-E expansion. Our calculation extends the results of Weldon [40] to massive fermions, but more importantly, we take into account another source of E-dependence in the free energy of the system. Since the electric field exerts a force on the charges in the medium, an equilibrium can only be reached if there is an equal force of opposite sign to cancel this effect. This can be any force, e.g. the strong force in a QCD medium or the electric repulsion in the case of QED. In our case of non-interacting fermions it solely stems from diffusion forces, or in other words the degeneracy pressure. This means that for any electric field profile, a non-trivial charge density distribution develops and keeps the medium in equilibrium. The electric susceptibility has to take this into account, by considering the free energy change added by the redistribution of the charge density. On a more pragmatic level, we find that this contribution is needed to make the susceptibility infrared finite, otherwise one would find a divergence in the homogeneous electric field limit. Our method of calculation can also be applied to magnetic fields, however in this case such contributions do not arise as magnetic fields do not lead to non-trivial equilibrium charge distributions.
Comparing our method to results obtained by using the Schwinger fermion propagator in a homogeneous electromagnetic background, the details of which are described in Ref. [37], we find that the two differ by an ordering of limits. Within the small-E expansion, the electric field approaches zero first, followed by the infinite volume limit (or, infinite wavelength limit for oscillatory fields). In contrast, the Schwinger propagator method involves a direct calculation in the infinite volume for homogeneous fields, and a zero field limit at the end. The two methods give rise to electric susceptibilities that differ by a finite term, see the corresponding high-temperature expansions in Eqs. (3.15) and (3.22). In turn, since magnetic fields do not involve such an infrared divergence, the two limits can be exchanged for the magnetic susceptibility and, accordingly, both approaches give the same result.
Our method can be generalized to include QCD interactions. In that case, the diffusion forces are accompanied by the strong force, affecting the equilibrium charge distribution and, thus, leading to corrections in the electric susceptibility. These QCD effects may be taken into account via perturbation theory (at high temperatures) or in full lattice QCD simulations. One advantage of our approach is that it allows for a straightforward implementation on the lattice, similarly to the magnetic field case, by means of current-current correlators [22]. The analogous density-density correlators for the electric response have also been explored in Ref. [54]. In contrast, QCD simulations at nonzero electric fields would not be possible using standard Monte Carlo methods due to the complex action problem.
Finally, we stress that the calculations presented in this paper (for both approaches) only concern the real part of the thermal free energy, and we do not discuss the imaginary part at T = 0. It is important to distinguish the two. The thermal effect is the polarization of charges that are generated due to thermal fluctuations, while the T = 0 effect is the decay of the vacuum due to virtual charges via the Schwinger mechanism [34].

Acknowledgments
This research was funded by the DFG (Emmy Noether Programme EN 1064/2-1 and SFB TRR 211 -project number 315477589). The authors are grateful to Holger Gies and Bo-Sture Skagerstam for their enlightening insight on calculations of the effective action and to Andrei Alexandru, Bastian Brandt, Sándor Katz, Tamás Kovács and Zsolt Szép for helpful discussions.

A Calculation of the photon polarization tensor
We are looking at the real and non-zero temperature part of Π T =0 00 first, therefore we need to evaluate the integral where we introduced E p = p 2 + m 2 and the trace is in Dirac space. Evaluating the latter together with the trivial p 0 integral allows us to take the k 0 → 0 limit to find a manifest real expression as the sum of a complex number and its conjugate. Then Carrying out the angle integrals and the trivial → 0 limit yields In order to carry out the limit in (3.9), Π T =0 00 (k 0 = 0, k) needs to be expanded around k = 0 to quadratic order. The leading order can simply be obtained by taking the limit k → 0 and yields Π T =0 00 (0, k = 0) = 2 e 2 π 2 ∞ 0 dp n F (E p ) As a check we expand this result around infinite temperature using standard HTE tools and recover the well known result for the photon mass in a plasma (see, e.g., Ref. [40]), The O(k 2 ) term in the expansion of the polarization tensor needs a more careful analysis. Taking the second derivative with respect to k can be replaced in this case by a k 2 derivative and a multiplication by 2 (the difference is higher order in k), leading to (A.6) where P denotes principal value integration. Naively taking k → 0 leads to an infrared divergent integral, arising from the non-commutativity of the principal value integral and the k → 0 limit, since for vanishing external momentum the pole is shifted onto the lower border of the integration. One particular way of handling this problem is by adding 0 = − e 2 π 2 P ∞ 0 dp 4 3 (A.7) Then the k → 0 limit of the sum of (A.6) and (A.7) is which is well behaved around p = 0. One can do an integration by parts on the term proportional to 1/3p 2 and change the integration variable to ω = E p to cast the final result into the simpler form We also need the sum of the three spatial components Π T =0 S , for the magnetic susceptibility. Before carrying out the Dirac trace, it reads Tr γ i ( / P + / K + m)γ i ( / P + m) π E p n F (|p 0 |)δ(P 2 − m 2 ) (P + K) 2 + m 2 + i .
(A.10) Evaluating the trace, the sum for the spatial components as well as the p 0 integral yields for k 0 = 0 After carrying out the angle integral, the → 0 limit of the real part reads ∞ 0 dp p 2 n F (E p ) E p 2 + k 2 + 4p 2 4kp log (k − 2p) 2 (k + 2p) 2 . (A.12) Then taking the k → 0 limit of the second derivative with respect to k is straightforward, and we obtain