QCD phase diagram and equation of state in background electric fields

The phase diagram and the equation of state of QCD is investigated in the presence of weak background electric fields by means of continuum extrapolated lattice simulations. The complex action problem at nonzero electric field is circumvented by a novel Taylor expansion, enabling the determination of the linear response of the thermal QCD medium to constant electric fields – in contrast to simulations at imaginary electric fields, which, as we demonstrate, involve an infrared singularity. Besides the electric susceptibility of QCD matter, we determine the dependence of the Polyakov loop on the field strength to leading order. Our results indicate a plasma-type behavior with a negative susceptibility at all temperatures, as well as an increase in the transition temperature as the electric field grows.


I. INTRODUCTION
The phase structure of Quantum Chromodynamics (QCD) in the presence of background electromagnetic fields is an essential attribute of the fundamental theory of quarks and gluons and, accordingly, a subject of active theoretical research.The electromagnetic response of the QCD medium is relevant for a range of physical situations, e.g. the phenomenology of heavy-ion collisions, the description of neutron star interiors or the evolution of our universe in its early stages, see the reviews [1,2].If in these settings the electromagnetic fields are sufficiently long-lived compared to the strong scale, it is appropriate to consider QCD matter in a background magnetic or electric field in equilibrium.
Before equilibration, electric fields E induce a dynamical response via the electrical conductivity of the medium [3].The subsequently emerging equilibrium necessarily involves -in contrast to the case of magnetic fields B -an inhomogeneous charge distribution n(x) in the thermal medium while having constant temperature T everywhere [4].The distribution is uniquely fixed by the requirement that pressure gradients and electric forces cancel each other and thus no currents flow [5].The equilibrium system is therefore described by a local canonical statistical ensemble, where n(x) is held fixed.It differs from the grand canonical ensemble parameterized by chemical potentials, employed usually at E = 0.This aspect renders comparisons between equilibrium systems at E > 0 and E = 0, e.g. by means of lattice simulations, problematic.
Moreover, the proper definition of the equilibrium state at E > 0 requires infrared regularization (e.g. a finite spatial volume V ) that prevents charges to be accelerated to infinity.As we have demonstrated recently within perturbative QED [6], the E → 0 and V → ∞ limits of this setup do not commute at nonzero temperature.This renders approaches based on Schwinger's exact E > 0 infinite-volume propagator [7] and infrared-regularized weak-field expansions in the manner of Weldon [8] inherently different.For a certain physical setting, the boundary conditions determine which is the appropriate limit to consider.The generalization of these ideas to the case of QCD enables one to explore the impact of background electric fields on strongly interacting matter as well as the associated phase diagram: our objectives in the present letter.
The impact of magnetic fields on the QCD crossover [9,10] and the corresponding phase diagram is well understood and has been studied extensively on the lattice [11][12][13][14][15], as well as within models and effective theory approaches (for a recent review, see Ref. [16]).In contrast, electric fields render the QCD action complex, hindering standard lattice simulations.Alternative approaches include Taylor-expansions [17][18][19][20], calculations at imaginary electric fields [21][22][23][24][25][26] and simulations with electric fields that couple to the isospin charge of quarks [27].Still, there are no existing results for the QCD equation of state nor the phase diagram.The latter has only been studied within effective theories like the linear σ model [28], variants of the Nambu-Jona-Lasinio (NJL) model [29][30][31][32] and the Euler-Heisenberg effective action [33].These calculations are all based on the Schwinger propagator.
In this letter, we determine the QCD equation of state and the phase diagram on the lattice for the first time for weak background electric fields.The complex action problem is circumvented via a Taylor-expansion: this corresponds to the Weldon-type regularization of the electrically polarized thermal medium and is the proper description of a finite system, where equilibration takes place in the presence of a weak electric field.The expansion is based on the method we developed in Refs.[6,34], and resembles the analogous approach for background magnetic fields [35][36][37].Besides the leading coefficient -the electric susceptibility of QCD matter -we also determine the leading series of the Polyakov loop.Using this observable, we construct the phase diagram and demonstrate that the transition temperature increases as E grows -contrary to existing model predictions, e.g.[31].Finally, we demonstrate that lattice simulations at nonzero imaginary electric fields cannot be used to directly calculate the electric susceptibility due to the singular change of ensembles between E = 0 and iE ̸ = 0. Some of our preliminary results have already been presented in Ref. [34].

II. LATTICE SETUP
QCD matter in thermal equilibrium is a medium that can be polarized by weak background electromagnetic fields.The associated static linear response is characterized by the electric and magnetic susceptibilities (we employ the same notation as in Ref. [6]).These are defined via the matter free energy density f , Here, the subscript b indicates that both susceptibilities contain ultraviolet divergent terms that must be subtracted via additive renormalization, see below.The elementary charge e is included so that we can work with the renormalization group invariants eE and eB.
The matter free energy density can be rewritten using the partition function Z of the system.Using the rooted staggered formalism of lattice QCD, it is given by the Euclidean path integral over the gluon links U , where β = 6/g 2 is the inverse gauge coupling and m f denotes the quark masses with f = u, d, s running over the quark flavors.The simulations are done in a periodic spatial volume V = L 3 with linear size L. Note that Z corresponds to the grand canonical ensemble; its relation to the canonical one at E > 0 is discussed in App. A. In Eq. ( 2), S g is the gluon action (in our discretization, the tree-level improved Symanzik action) and / D f is the staggered Dirac operator (including a twofold stout smearing of the links) that contains the quark charges q u /2 = −q d = −q s = e/3.The quark masses are set to their physical values as a function of the lattice spacing a [38].Further details of the action and of our simulation algorithm are given in Refs.[12,39].
The electromagnetic vector potential A ν enters the Dirac operator in the form of temporal parallel transporters u ν,f = exp(iaq f A ν ) multiplying the gluon links U ν .We choose a gauge where A 0 (x 1 ) represents the electric field and A 2 (x 1 ) the magnetic field (both pointing in the x 3 direction).While magnetic fields are identical in Minkowski and Euclidean space-times, the vector potential relevant for the electric field undergoes a Wick rotation so that A 4 = iA 0 , similarly to the case of a chemical potential µ.Finally we mention that in our setup, quarks do not couple to dynamical photons but only to the external gauge field.The independent thermodynamic variable is the field E that enters the Dirac operator, analogously to the situation for magnetic fields [40].

III. OBSERVABLES
As we demonstrated in Ref. [6], the susceptibilities of Eq. ( 1) are related to derivatives of the electromag-netic vacuum polarization tensor with respect to spatial momenta.For our gauge choice, these relations read in terms the Euclidean polarization tensor Π µν , with a spatial momentum k = (k 1 , 0, 0, 0).In other words, the zero momentum limit is considered at vanishing time-like frequency, reflecting the static nature of the susceptibilities.The negative sign for ξ b in (3) appears due to the Wick rotation of the electric field.We highlight that the equilibrium systems at different values of E exhibit different charge profiles n(x 1 ), and this implicit E-dependence is taken into account properly in Eq. ( 3) for the calculation of ξ b [6].In fact, without this contribution, ξ b would diverge in the k 1 → 0 limit.The vacuum polarization tensor is defined as the correlator of the electromagnetic current j µ = f q f e ψf γ µ ψ f , for which we use the conserved (one-link) staggered vector current.It is convenient to evaluate (3) in coordinate space, where the bare susceptibilities become [6,34] containing the second moment of a partially zeromomentum projected two-point function The Grassmann integral over quark fields is understood to be implicitly carried out on the right hand side of the last equation.Both susceptibilities undergo additive renormalization.This originates from the multiplicative divergence in the electric charge e [36,41,42].Being temperatureindependent, the divergence cancels in (8) which sets ξ = χ = 0 at zero temperature.In fact, at 44 , implying that the bare magnetic and electric susceptibilities coincide up to a minus sign.To renormalize the electric susceptibility, we can therefore employ the existing results for χ b (T = 0) from Ref. [36].
Next we consider the bare Polyakov loop operator, Its expectation value is related to the free energy of a static, electrically neutral color charge and is often taken as a measure of deconfinement.Just as for ξ b , the contribution of the equilibrium charge profile needs to be taken into account for the E-dependence of the Polyakov loop as well.As we show in App.A, the proper second-order expansion of ⟨P b ⟩ is given by the correlator where the superscript n on the left denotes that the derivative is evaluated along the equilibrium condition specified by the local charge profiles.Analogously, the magnetic derivative of ⟨P b ⟩ can be obtained by replacing 22 in Eq. (10), although in that case nontrivial charge distributions do not appear.
The bare Polyakov loop is subject to multiplicative renormalization [43], where the renormalization factor is independent of the background field and has been determined for our lattice spacings in Ref. [44].This renormalization fixes P = P ⋆ at T = T ⋆ and E = B = 0.In our renormalization scheme we choose T ⋆ = 162 MeV and P ⋆ = 1.

IV. RESULTS: SUSCEPTIBILITY
We have measured the zero-momentum projected correlator (7) for a broad range of temperatures on N t = 6, 8, 10 and 12 lattice ensembles.Finite volume effects were checked using 16 3 × 6 and 24 3 × 6 lattices.We report on the details of the measurements and the analysis in App.B. The correlator is convolved with the quadratic kernel according to Eq. ( 6) to find the bare electric susceptibility ξ b , and its renormalization ( 8) is carried out by subtracting the zero-temperature contribution.
The negative of the so obtained ξ is plotted in the upper panel of Fig. 1.A continuum extrapolation is performed via a multi-spline fit of all data points, taking into account O(a 2 ) lattice artefacts.The systematic error of the fit is estimated by varying the spline node points and including O(a 4 ) discretization errors in the fit at low temperatures.For all temperatures we observe ξ < 0, translating to an electric permittivity below unity -a characteristic feature of plasmas [45].At high T , our results may be compared to the high-temperature limit calculated for non-interacting quarks of mass m [6], where N c = 3 is the number of colors [46].In full QCD, the quark mass is replaced by a QED renormalization FIG. 1. Upper panel: the negative of the renormalized electric susceptibility as a function of the temperature for four lattice spacings (colored symbols) and a continuum extrapolation (orange band).Lower panel: continuum extrapolated magnetic (green) and electric (orange) susceptibilities (solid) compared to leading-order perturbation theory (dashed).
scale µ QED that can be determined at T = 0 and is found to be µ QED = 115(6) MeV [36], close to the mass of the lightest charged hadron i.e. the pion.Moreover, QCD corrections are included by taking into account O(α s ) effects in the prefactor, the QED β-function [36,47].The associated thermal scale is varied between πT and 4πT for error estimation.The so obtained curve lies very close to our results at high temperature, as visible in the lower panel of Fig. 1, where we also show the corresponding results for χ from Ref. [36].

V. RESULTS: PHASE DIAGRAM
Next we turn to the Polyakov loop.Its leading expansion is given by Eq. ( 10), containing the correlator of the bare observable with −G peratures, i.e. a reduction of the Polyakov loop by the electric field.Finite volume effects are found to be small, although the results at low temperature have large statistical uncertainties.Using the results for the Polyakov loop at E = 0 [44] and the multiplicative renormalization factor from Eq. ( 11), we construct the E-dependence of ⟨P ⟩, see the lower panel of Fig. 2. The Polyakov loop is known to exhibit a smooth temperature-dependence, so that a precise determination of its inflection point is cumbersome already at E = 0.As an alternative, we associate the transition temperature T c with the point where ⟨P ⟩ = 1 holds.Defined in this manner, the lower panel of Fig. 2 clearly shows that T c is increased by E.
We repeat this analysis for all four lattice spacings.Our results for the transition temperature are shown in Fig. 3, confirming the significant enhancement of T c as the electric field grows.We perform the continuum extrapolation by a quadratic fit of T c (E) taking into account O(a 2 ) lattice artefacts.To estimate the systematic error, we vary the fit range and also allow a quartic term in the fit.The fits are found to be stable for the region E ≲ 0.3 GeV 2 [48].The curvature of the transition line is found to be Furthermore, we find the transition to get stronger as E grows, revealed by an enhancement of the slope of the Polyakov loop as a function of T , see Fig. 2.However, due to the large uncertainties at low temperatures, we FIG. 3. Transition temperature as a function of the electric field for different lattice spacings (colored symbols) and a continuum extrapolation (yellow band).Higher-order effects in eE become non-negligible for eE ≳ 0.3 GeV 2 , indicated by the dashed section of the fits.
cannot make a quantitative statement about this aspect.

VI. RESULTS: IMAGINARY ELECTRIC FIELDS
Finally, we consider lattice simulations at constant imaginary electric fields iE.In a finite periodic volume at nonzero temperature, the allowed field values are quantized as ieE = 6πT /L • N e with the 'flux' quantum N e ∈ Z [49].This setup does not correspond to the analytic continuation of the local canonical ensemble as described in the introduction.Nevertheless, it involves a global constraint: the total electric charge in the periodic volume vanishes.As a consequence, this setup is independent of the global imaginary chemical potential iµ.Indeed, including any iµ ̸ = 0 can be canceled in the gauge field by a mere coordinate translation x 1 → x 1 − µ/E.This is in stark contrast to the situation at E = 0, where a dependence on iµ is naturally present.
To discuss this issue, we neglect gluonic interactions in the following.In this simplified setting, we can calculate the free energy density directly via exact diagonalization of the Dirac operator [50].In the right side of Fig. 4 we show the results for ∆f = f − f (E = µ = 0) obtained on a 200 3 × 20 lattice with quark mass m/T = 0.08.As expected, f is found to be independent of the imaginary chemical potential in the whole range 0 ≤ iµ ≤ πT at any iE ̸ = 0.The comparison to a larger volume 300 3 × 20 shows that the smallest allowed electric field value approaches zero in the thermodynamic limit, but lim iE→0 f (iE) ̸ = f (iE = 0).Instead, the data points rather accumulate towards the average of f over all possible imaginary chemical potential values -i.e. a canonical setup where the total charge is constrained to zero.Altogether, we conclude that the dependence of f on iE is singular at E = 0 in the thermodynamic limit, ren- dering simulations with homogeneous imaginary electric fields unsuited for the evaluation of ξ.
In addition, the left side of Fig. 4 shows ∆f for oscillatory imaginary electric fields with the profile E(x 1 ) = E √ 2 cos(2πnx/L).In this case, the role of the infrared regulator is played by the wave number n and not by the volume.Moreover, here E is a continuous variable but n ∈ Z is discrete.This setup does not fix the overall charge and, therefore, maintains the dependence of f on iµ.Indeed, the results reveal a continuous behavior as a function of iE and iµ.However, as visible in the plot, the results again approach a singular behavior as the infrared regulator is removed: the curves collapse to a set of iµ-independent nodepoints approaching the iE = 0 axis.In particular, the curvature of f with respect to iE diverges for n → 0. Thus, the homogeneous limit of the setup with oscillatory imaginary fields reproduces what we have already seen for the homogeneous case.

VII. DISCUSSION
In this letter we studied the thermodynamics of QCD at nonzero background electric fields E via lattice simulations with physical quark masses.To avoid the complex action problem at E > 0, we employed a leading-order Taylor-expansion.This approach is more complicated than the analogous expansion in a chemical potential, because the impact of E on the equilibrium charge distribution needs to be taken into account [6].Our results, measured on four different lattice spacings and extrapolated to the continuum limit, demonstrate two main effects.First, that QCD matter is described by a negative electric susceptibility at all temperatures.Second, that the QCD transition, as defined in terms of the Polyakov loop, is shifted to higher temperatures as the electric field grows, leading to the phase diagram in Fig. 3. Furthermore, we showed that lattice simulations employing imaginary electric fields cannot be used to directly assess these aspects due to a singular behavior around E = 0.
We mention that the susceptibility and the phase diagram are both encoded by the thermal contributions to the real part of the free energy density.These are therefore not impacted by Schwinger pair creation, which is related to the imaginary part of f and is known to be independent of the temperature [51,52].In other words, the equilibrium charge profile and the polarization of the medium are related to the distribution of thermal charges and not of those created from the vacuum via the Schwinger effect.
Finally we point out that calculations within the PNJL model [31], employing the Schwinger propagator, predict the opposite picture for the phase diagram as compared to our findings.Whether the same tendency holds for the Weldon-type regularization within this model, is an open question calling for further study.Besides this aspect, the PNJL model is known to miss important gluonic effects in the presence of electromagnetic fields and fails to correctly describe the phase diagram at B > 0 [16].It would be interesting to see whether improvements that were found to correct these shortcomings of the model in the magnetic setting [53] also work in the E > 0 case.transformation [6], The local chemical potential is fixed by the requirement that diffusion and electric forces cancel, i.e. μ(x 1 ) = −eEx 1 .This choice corresponds to a globally neutral system, where the volume average of the chemical potential vanishes.Including the bare Polyakov loop in the action with a coefficient α and taking the derivative of (A1) with respect to α at α = 0 results in giving the expectation value of the Polyakov loop in the local canonical ensemble.Taking the second total derivative of (A2) with respect to eE, and evaluating it at E = 0 (implying μ = 0), we obtain for (10) with The Polyakov loop operator P b does not depend explicitly on the electric field nor on the chemical potential.The derivatives of ⟨P b ⟩ therefore merely involve the derivative of the weight in the path integral (2).Let us first discuss φ µ .The chemical potential multiplies the volume integral of j 4 in the Euclidean action (before integrating out fermions), therefore (A5) where we used that ⟨j 4 (y)⟩ = 0 due to parity symmetry.Substituting the integration variable z by u = z − y, exploiting the translational invariance of the correlators and using the definition (7) of the projected correlator, we arrive at Next, we turn to φ E .This time, the Euclidean action contains the four-volume integral of ieA 4 (y 1 ) • j 4 (y 1 ) with A 4 (y 1 ) = −iE y 1 .The second derivative therefore becomes 1 )/2 and use that the second term can be replaced by y 2 1 as it multiplies a factor that is symmetric under the exchange of y 1 and z 1 under the integrals.With the same variable substitution as above, the use of translational invariance of the correlators this time gives The second term in (A8) is clearly divergent in the thermodynamic limit.Coming back to (A3), we see that this infrared singular term exactly cancels in φ n E , rendering the curvature of the Polyakov loop expectation value finite when evaluated along the equilibrium condition involving the inhomogeneous charge profile.Finally, employing the u 1 ↔ −u 1 symmetry of the E = 0 system, we end up with Eq. ( 10) of the main text, involving the second moment G (2) 44 defined in Eq. ( 6).There is one more aspect regarding the dependence of the Polyakov loop on E that deserves mentionting.In lattice simulations with constant imaginary electric fields iE at nonzero temperature, the Polyakov loop was observed to develop a local phase proportional to the local vector potential arg P b (x 1 ) ∝ ieEx 1 /T [26] (see also the analogous study [54]).This results from the preference of local Polyakov loops towards different center sectors for different x 1 .Together with the quantization condition for the imaginary electric flux, this corresponds to a topological behavior of the Polyakov loop angle winding around the lattice.Thus, the volume-averaged P b vanishes in these iE ̸ = 0 simulations, as opposed to its nonzero value at E = 0, showing the singular change of relevant ensembles as the electric field is switched on.Again, we conclude that simulations with homogeneous imaginary electric fields cannot be used for a direct comparison to the E = 0 system.

Appendix B: Correlators
Here we discuss the determination of the correlator G 44 (x 1 ) and the bare electric susceptibility in more detail.The density-density correlator is calculated using O(1000) random sources located on three-dimensional x 1 -slices of our lattices.We take into account both connected and disconnected contributions in the two-point function.More details regarding the implementation can be found in [35].We note that the same two-point function is required, at zero temperature, for the calculation of the hadronic contribution to the muon anomalous magnetic moment, see e.g.Ref. [55].
In Fig. 5 we show the zero-momentum projected density-density correlator ⟨G 44 ⟩ as a function of the coordinate at T ≈ 176 MeV.For comparison, the currentcurrent correlator ⟨G 22 ⟩, relevant for the magnetic response, is also included.A substantial difference is visible, reflecting the absence of Lorentz symmetry at this high temperature.It is interesting to note the systematic oscillation of ⟨G 44 (x 1 )⟩ between even and odd distances -related to the use of staggered fermions -which is however absent for ⟨G 22 (x 1 )⟩.
To assess finite volume effects, we consider the convolution (6) and truncate it at a distance x max 1 .The so obtained truncated susceptibility approaches the full susceptibility at x max 1 = L/2 and is plotted in Fig. 6 for two different volumes.On the 24 3 × 6 lattice, the plot shows that contributions coming from the middle of the lattice volume are exponentially small, as expected.Moreover, at x max 1 = L/2, the results obtained on the two volumes agree with each other within errors.The dominant systematic error for the determination of our final result is found to come from the continuum extrapolation, which is discussed in the main text.

( 2 ) 44 .
FIG. 2. Upper panel: the leading expansion coefficient of the bare Polyakov loop as a function of the temperature as obtained on our 24 3 × 6 (red) and 16 3 × 6 (green) ensembles.Lower panel: renormalized Polyakov loop at nonzero electric fields, constructed from the leading Taylor series.The crossing point with the dashed yellow line is identified with Tc.

FIG. 4 .
FIG.4.Free energy density as a function of homogeneous (right side) and oscillatory (left side) imaginary electric fields.The results for different imaginary chemical potentials correspond to the set of curves in the left (iµ grows from 0 to πT from the bottom to the top), while they lie on top of each other on the right.