The QCD equation of state in background magnetic fields

We determine the equation of state of 2+1-flavor QCD with physical quark masses, in the presence of a constant (electro)magnetic background field on the lattice. To determine the free energy at nonzero magnetic fields we develop a new method, which is based on an integral over the quark masses up to asymptotically large values where the effect of the magnetic field can be neglected. The method is compared to other approaches in the literature and found to be advantageous for the determination of the equation of state up to large magnetic fields. Thermodynamic observables including the longitudinal and transverse pressure, magnetization, energy density, entropy density and interaction measure are presented for a wide range of temperatures and magnetic fields, and provided in ancillary files. The behavior of these observables confirms our previous result that the transition temperature is reduced by the magnetic field. We calculate the magnetic susceptibility and permeability, verifying that the thermal QCD medium is paramagnetic around and above the transition temperature, while we also find evidence for weak diamagnetism at low temperatures.


Introduction
Quantum Chromodynamics (QCD) is the theory of the strong interactions. Its most important properties are the confinement of quarks and gluons at low energies and asymptotic freedom at high scales. Lattice simulations of QCD have unambiguously shown that at zero quark densities these two, fundamentally different regimes are connected by a smooth crossover-type transition [1,2]. This transition -through which the dominant degrees of freedom change from composite objects (hadrons) to colored quarks and gluons -has several characteristics of both theoretical and phenomenological relevance. Besides the nature and the (pseudo)critical temperature of the transition, one such characteristic is the equation of state (EoS), which is the fundamental relation encoding the thermodynamic properties of the system.
In particular, the EoS gives the equilibrium description of QCD matter, in terms of relations between thermodynamic observables like the pressure, the energy density or the entropy density. These observables enter hydrodynamic models that are used to describe the time evolution of the quark-gluon plasma (QGP) produced in heavy-ion collision experiments [3,4]. Besides its role in heavy-ion physics, the EoS affects the mass-radius relation of neutron stars [5] and enters cosmological models of the early universe, with implications, for example, for dark matter candidates [6].
Of particluar relevance is the response of the EoS to changes of the control parameters of the system. These parameters include the temperature, the chemical potentials conjugate to conserved charges and, in the present case, a background (electro)magnetic field = |B|. External magnetic fields play an important role in the evolution of the early universe [7], in strongly magnetized neutron stars [8] and in non-central heavy-ion collisions, see, e.g., the recent review [9]. Magnetic fields induce a variety of exciting effects in the thermodynamics of QCD -for example they significantly affect the phase diagram. The first results in this field were obtained using low-energy models and effective theories of QCD, see the summary in, e.g., Ref. [10]. The QCD transition has also been studied extensively on the lattice; we refer the reader to reviews on the subject in, e.g., Refs. [11][12][13]. A very relevant question in this respect has been the dependence of the transition temperature and of the nature of the transition on the magnetic field. In this paper we also address this issue.
Our main objective is to determine the QCD EoS around the crossover transition at vanishing chemical potentials, for nonzero background magnetic fields. To this end we develop a generalization of the so-called integral method [14], which relies on an integration in the quark masses up to asymptotically large values (a similar integration in the quark masses at = 0 was also considered in Ref. [15]). We calculate thermodynamic observables including the pressure, the energy density, the entropy density, the interaction measure, the magnetization and the susceptibility for magnetic fields of up to = 0.7 GeV 2 for a wide range of temperatures 110 MeV < < 300 MeV, allowing for a comparison with the Hadron Resonance Gas (HRG) model and with perturbation theory, at low and high temperatures, respectively. At high we demonstrate that several aspects of perturbative QED physics are encoded in the EoS observables. Furthermore, our results confirm the observation made in Refs. [16,17] that the transition region is shifted to lower temperatures as grows. Another phenomenological consequence of the magnetic field is that the pressure -if defined as the response against a compression at fixed magnetic flux, see precise definition in Sec. 2 below -becomes anisotropic, and a significant splitting between the components parallel and perpendicular to B is developed.
The change in the EoS due to the magnetic field has further theoretical implications. QCD matter may be thought of as a medium with either para-or diamagnetic properties. We establish that around and above the transition region the magnetization is positive and thus the thermal QCD medium behaves as a paramagnet, with a magnetic permeability larger than unity. A possible implication of this paramagnetism for heavy-ion collisions has been pointed out recently in Ref. [18]. In addition, we present evidence for the emergence of a weakly diamagnetic region at low temperatures due to pions. This paper is organized as follows. First we discuss thermodynamic relations in the presence of the magnetic field from a general point of view in Sec. 2. We proceed by describing the lattice methods that are used to determine the EoS in Sec. 3, with special emphasis on the implications of flux quantization and electric charge renormalization. Sec. 4 contains our main results, followed by the conclusions in Sec. 5.

Thermodynamics in an external magnetic field
The fundamental quantity of thermodynamics is the free energy or thermodynamic potential. In terms of the partition function of the system it reads ℱ = − log . In the presence of an external magnetic field the density = ℱ/ of the free energy in a finite spatial volume can be written as [19] = − = total − − · ℳ, (2.1) where is the energy density of the medium, the entropy density and ℳ the magnetization. Without loss of generality, the magnetic field B = e is taken to point in the direction, and for later convenience, is given in units of the elementary charge > 0. Note that the total energy of the system, total = + field , includes the energy of the medium as well as the work necessary to maintain the constant external field, field = · ℳ [20]. The two expressions in Eq. (2.1) thus correspond to two different conventions for the definition of the energy density. The entropy density and the magnetization can be obtained as 2) The corresponding differential relation for the pressure is somewhat more involved. Since the magnetic field marks a preferred direction, the pressures in the transverse (perpendicular to B) and in the longitudinal (parallel to B) directions may be different. In Ref. [21] we have shown that this possible anisotropy depends on the precise definition of . Writing the volume as the product of linear extents = , the pressure components are related to the response of the system to compressions along the corresponding directions, i.e.
In order to unambiguously define , we have to specify the trajectory in parameter space, along which the partial derivative is evaluated. In Ref. [21], we have distinguished between a setup where the magnetic field is kept fixed during the compression (the " -scheme"), and a setup where the magnetic flux Φ = · is kept fixed (the "Φ-scheme"). The -scheme results in isotropic pressures, whereas the Φ-scheme gives anisotropic pressures: The difference in the transverse components for the Φ-scheme is due to the fact that the compressing force in this case also acts against the magnetic field. Note that the definition of the pressures as spatial diagonal components of the energy-momentum tensor exhibits the Φ-scheme anisotropy [22]. This is due to the fact that the energy-momentum tensor is usually defined through the variation of the action with respect to the metric at fixed Φ, see Ref. [21]. In contrast to , , the longitudinal pressure is independent of the scheme, and in the thermodynamic limit → ∞ simplifies to Note that the appropriate scheme to be used depends on the physical situation that one would like to describe. In particular, it is specified by the trajectory ( ), along which the compression perpendicular to the magnetic field proceeds. As will be explained below, in lattice regularization it is natural to keep the flux fixed, and thus, the lattice measurements correspond directly to the Φ-scheme. However, this does not represent a limitation of the lattice approach, since one can easily translate from one scheme into another. The pressure components for a general ( ) trajectory ("general scheme") can be found by combining our results for the longitudinal pressure and for the magnetization (both are contained in the ancillary files submitted to the arXiv), This relation reproduces the -and Φ-schemes, Eq. (2.4), for the trajectories ( ) = and ( ) = Φ/( ), respectively. Another important observable for the EoS is the interaction measure (trace anomaly), which contains the energy of the medium and the three pressures. Thus, also depends on the scheme 1 : whereas the energy density (like ) is by construction scheme-independent, Eqs. (2.1), (2.5) and (2.8) reveal that the entropy density can also be calculated as Finally, the derivative of the magnetization with respect to at vanishing magnetic field gives the magnetic susceptibility, (2.11)

Lattice observables and methods
In what follows we consider a spatially symmetric lattice with isotropic lattice spacing . Here the temperature and the three-volume are given by where and are the number of lattice sites along the spatial and temporal directions, respectively. Using conventional Monte-Carlo methods the free energy ℱ = − log itself is not accessible on the lattice, but only its derivatives with respect to the parameters of the theory. For the case of 2 + 1 flavor QCD coupled to a constant external magnetic field, these parameters are the inverse gauge coupling = 6/ 2 , the lattice quark masses ( = , , labeling the flavors) and the magnetic 1 One may understand the scheme-dependence of as follows. The trace anomaly represents the response to a rescaling of the length scale in the system. To define this rescaling unambiguously, the trajectory ( ) has to be specified, i.e. a scheme has to be chosen. Hence, becomes scheme-dependent. As a simple example, consider the magnetic field in the absence of particles. Taking into account the energy 2 /2 of the magnetic field, one obtains (Φ) = 0, while ( ) = 2 2 . Notice that in the Φ-scheme a dimensionless number characterizes the magnetic field, whereas in thescheme we introduced a dimensionful parameter into the system. This is reflected by the vanishing of the trace anomaly in the former case, and the nonzero value of in the latter.
flux Φ = ( ) 2 . In particular, in the staggered formulation of lattice QCD, is written as where = ( / + ) is the fermion matrix, and the quark charges are set to = = − /2 = − /3. Note that the magnetic field has no dynamics, therefore the charge (or, the elementary charge ) and the magnetic field do not appear separately in the partition function, but always in the combination (or ). The constant Maxwell term 2 /2 is independent of the physical properties of the thermal QCD medium and plays no role in the thermodynamics of the system. It only enters in the renormalization prescription, see Sec. 3.2 below.
We work with the tree-level improved Symanzik gauge action , and stout improved staggered quarks in the fermionic sector. The detailed simulation setup is described in Refs. [16,23]. The quark masses are set to their physical values along the line of constant physics (LCP). This means that are tuned as functions of in a way that "physics remains the same", that is to say, ratios of hadron masses measured on the lattice coincide with their experimental values. This defines the physical quark masses ph for each value of . In particular, our LCP is set by fixing the ratio of the kaon decay constant to the pion mass / and the kaon decay constant to the kaon mass / . This results in the fixed ratio of quark masses = ≡ = /28.15. The lattice spacing ( ) is set using . For additional details on this procedure, see Ref. [15].
The derivatives of log with respect to and are the gauge action density and the quark condensate densities, The interaction measure, Eq. (2.7), can be given in terms of the response of the free energy to an overall change of length scales in the system. On the lattice this amounts to a derivative with respect to the lattice spacing . Employing the -dependence of the lattice parameters and , the densities of Eq. (3.3) enter the Φ-scheme interaction measure in the following way: Note that for the -scheme interaction measure (see Eq. (2.8)), an additional term containing the derivative with respect to the lattice flux 2 appears. For convenience, we also define the change due to for any observable as The renormalization of the above observables will be discussed in Sec. 3.2.

Flux quantization and methods to determine the magnetization
Due to the periodic boundary conditions, the magnetic flux traversing the finite lattice is quantized as where we took into account that the smallest charge in the system (that of the down quark) is = /3. Note that since the flux is quantized, the lattice setup automatically corresponds to the Φ-scheme defined in Sec. 2. Moreover, due to the quantization condition, differentiation with respect to is in principle ill-defined and therefore the magnetization of Eq. (2.2) is not accessible directly. Recently, several methods were developed to circumvent this problem, which we summarize briefly below.
• Anisotropy method. One can make use of the relation (2.4) for the Φ-scheme, and express the magnetization as the difference between the longitudinal and transverse lattice pressures. These can be measured as derivatives of log with respect to anisotropy parameters. This approach was developed and successfully applied in Refs. [18,21]. The advantage of the method is that ℳ is directly obtained as an expectation value for any , while its drawback is that anisotropy renormalization coefficients also need to be determined.
• Half-half method. Instead of the uniform (and, thus, quantized) magnetic field, one can work with an inhomogeneous field which has zero flux, e.g. one that is positive in one half and negative in the other half of the lattice. Since the field strength is now a continuous variable, derivatives of log with respect to are well defined and can be measured on a = 0 lattice ensemble [24]. The second-order derivative directly gives the magnetic susceptibility. However, higher-order terms become increasingly noisy, which limits the applicability of the approach to low fields. Note moreover that the discontinuities in the magnetic field may enhance finite volume effects.
• Finite difference method. The derivative of log with respect to is an unphysical quantity due to the quantization Eq. (3.6). Still, this derivative can be measured for any real value of , and its integral over between two integer values gives the change in log between these two fluxes. In this way, log ( + 1) − log ( ) is constructed as the integral of an oscillatory function. The method is in principle applicable for any magnetic field, but 10-20 independent simulations are necessary to go from one integer flux to the next, making large magnetic fields computationally expensive [25,26].
• This work: generalized integral method. The method we will use in the present paper is based on two observations: that magnetic fields have no effect in pure gauge theory, and that the infinite quark mass limit of QCD (at a fixed magnetic field ≪ 2 ) is pure gauge theory. Based on this, the change in log due to the magnetic field can be expressed as an integral of the quark condensate differences Δ¯over the quark masses, including unphysically heavy quarks. On a finite lattice, this integral is well regulated and can be calculated in a controlled manner by using 10-20 independent simulations for any given value of the magnetic field. Most of these simulations are at large quark masses, where the computation is significantly cheaper. Furthermore, the method automatically gives information on the mass-dependence of log as well. This approach was sketched in Ref. [27] and will be described in detail below.

Renormalization
The free energy density contains additive divergences in the cutoff -i.e. in the inverse lattice spacing. These divergences are independent of , except for one logarithmic divergence of the form − 1 ( ) 2 log( ), where is a renormalization scale. This term is canceled through a redefinition of the energy 2 /2 of the magnetic field itself [28], Eq. (3.7) is equivalent to a simultaneous renormalization of the wave function (magnetic field ) and of the electric charge . The combination is renormalization group invariant and, as such, unaffected by this transformation: (3.8) The purely magnetic contribution 2 /2 is trivial and can be omitted from the Lagrangian. Therefore, the renormalization of the free energy amounts to adding the counter-term 1 ( ) 2 log( ) to Δ . In the following it will be advantageous to consider the extensive quantity Δ log = − 4 Δ at zero temperature, in a box of four-volume 4 . The counter-term then takes the form − 1 Φ 2 log( ) with the flux Φ = 2 . The coefficient 1 of the divergence is related to the QED -function [29][30][31].
Since the magnetic field is external, i.e. there are no U(1) degrees of freedom in the system, only the lowest order QED -function coefficient 1 appears in (however, with a full dependence on the QCD coupling, see Eq. (3.14) below).

Charge renormalization -free case
It is instructive to first discuss the renormalization procedure in the free case -i.e. for electrically charged quarks in the absence of strong interactions. In this case the free energy can be calculated analytically (see, e.g., Refs. [30][31][32]). We consider quark flavors of charges and, for simplicity, we assume degenerate masses = for all . The discussion is easily generalized to unequal masses. For = 3 colors, the QED -function coefficient reads At zero temperature, the expansion of Δ log in the magnetic field is given by where we also included the counter-term. Taking the derivative with respect to the mass of the quark flavor , we obtain the corresponding quark condensate at = 0, showing that the condensate contains no -dependent divergences, and that it is also independent of the renormalization scale . Note also that the sign of the magnetic field-induced change in the condensate is, to leading order, determined by the sign of free 1 . Since QED is not asymptotically free, free 1 is positive and the condensate undergoes magnetic catalysis at = 0 to quadratic order in (we have already presented this argument in Refs. [27,32]). Note that approaching the chiral limit (i.e. / 2 → ∞), the magnetic field-expansion in Eq. (3.11) becomes ill-defined. In fact, in the → 0 limit Δ¯f ree vanishes (see, e.g., Ref. [33]) for any magnetic field. However, the condensate difference is expected to be positive if any weak attractive interaction is turned on [34].
Let us now calculate the interaction measure. We resort to the Φ-scheme of Sec. 2, as this is the natural one in the lattice setup. Contrary to the case of the condensate, here the counter-term also contributes a finite term free 1 ( ) 2 . The remainder of Eq. (3.10) depends only on the combination , thus the derivative with respect to log is equivalent to that with respect to log . Using Eqs. (3.10) and (3.11), we therefore obtain which is again finite and -independent. The trace anomaly difference contains the two well-known sources of scale violation [35]: the classical breaking through the condensates 2 and the anomalous one through the running of the electric charge. In our case, the two contributions cancel each other to (( ) 2 ), since at this order drops out of Eq. (3.10). We shall return to this observation below.
We have seen that the condensate difference and the trace anomaly difference are independent of the renormalization scale. However, in order to define the renormalized free energy and pressures, we need to specify in Eq. (3.10). Setting the renormalization scale equal to the mass, = , means that all terms quadratic in the magnetic field are canceled. This scheme 3 is intrinsic to the Schwinger proper time representation [28], and coincides with the one used in Ref. [32]. Since in this scheme the expansion of log starts as ( ) 4 at = 0, so do the expansions of the pressures, of the energy density and of · ℳ. Therefore -being a linear combination of the former -has an expansion starting with a quartic term as well, consistent with Eq. (3.12).
We remind the reader that we excluded the renormalized pure magnetic energy 2 /2 above, which depends explicitly on the renormalization scale . To restore this term in log one needs to add where em ( ) = 2 ( )/(4 ) is the running QED coupling defined at the scale .

Charge renormalization -full QCD
We proceed by applying the renormalization prescription discussed for free quarks above to the case of full QCD. With the strong interactions taken into account, 1 will contain QCD corrections, which, in a perturbative expansion in the strong coupling take the form where the coefficients are independent of the quark masses and have been calculated in the MS scheme up to = 4 in Ref. [38]. Note that the running of the QCD coupling -governed by the QCD 2 Note that the usual definition of the condensate (with¯< 0) differs from our convention by a minus sign. 3 We remark that renormalization schemes with different choices for the scale have also been used in the literature. For example, is taken to be proportional to √ in the schemes employed in Refs. [36,37], which are connected to our choice by a finite (albeit mass-dependent) renormalization. Our scheme has the advantage that the leading magnetic field-dependence of the total free energy is simply 2 /2. Moreover, the → ∞ limit of the magnetization vanishes, in accordance with the expectation that magnetic fields should have no effect on static non-relativistic particles (see the discussion in Ref. [32]).
-function -induces a dependence of 1 on the regulator, which on the lattice amounts to a dependence on the lattice spacing . Thus, due to the asymptotically free nature of the strong interactions, QCD corrections vanish in the continuum limit, and 1 ( ) approaches its free value, as indicated in Eq. (3.14). We will see that for the lattice spacings we employ, these corrections are already tiny, see the right panel of Fig. 3 below.
The consistency with charge renormalization ensures that the free energy is again of the form Eq. (3.10). Contrary to the free case, inside the logarithm of the divergent term, an additional dimensionful hadronic scale Λ H appears, which may depend on . The expansion of Δ log at zero temperature then reads From this -in analogy to the free case -we can extract the leading dependence of the condensate difference and of the interaction measure difference on the magnetic field: We again conclude that the -dependent divergence is absent from the condensate [16,21]. Moreover, in the renormalization group invariant combination Δ¯, multiplicative divergences cancel as well. To quadratic order in the sign of the change of the condensate is related to the sign of 1 and to that of Λ H / , which we will revisit in Sec. 4.1. The interaction measure difference is also explicitly finite, as noted in Ref. [21], where we determined the gluonic and fermionic contributions to Δ (Φ) separately. Similarly to the free case, Eq. (3.12), Δ (Φ) receives a finite contribution from the counter-term in Δ log , which, due to Eq. (3.14), equals its free-case equivalent in the continuum limit (the QCD -function damps the second term in the square brackets as → 0). In the continuum limit, the contribution from the counter-term results in a cancellation to (( ) 2 ) in the total interaction measure, as already stated in Eq. (3.16). For later reference, the renormalization of Δ (Φ) thus reads To discuss the renormalization of log itself, we have to specify the renormalization scale. We may again choose such that the quadratic term in log at = 0 is completely subtracted in the renormalization process: = Λ H . This is the equivalent of the on-shell renormalization scheme in the free case. The renormalization prescription for the free energy (and, similarly, for the longitudinal pressure) at = 0 in this scheme reads where we defined as the operator that projects out the (( ) 2 ) term from an observable : (3.20) We remark that at finite temperature, thermal contributions induce additional finite terms that are quadratic in . Thus, the subtraction of [ ] is to be performed at = 0, as indicated in Eq. (3.20).

The integral method at nonzero magnetic fields
To determine the free energy -or, equivalently, the longitudinal pressure , see Eq. (2.5) -on the lattice, we employ a variation of the so-called integral method [14]. The basic idea is to construct by integrating its partial derivatives in Eq. (3.3) along a particular path in the parameter space spanned by the parameters { , , Φ}. Since the magnetization is not accessible as a derivative (see Sec. 3.1), one is only allowed to integrate along a constant-Φ trajectory in this parameter space. For a lattice of fixed size 3 × , the magnetic field thus changes as ∼ −2 ∼ 2 along such a path. Specifically, we consider a trajectory at constant Φ, from 1 to 2 with the quark masses tuned along the LCP ph . Then, the integral method is written down for the change Δ in the pressure: the difference of Δ at the two endpoints equals the integral of the gradient of Δ along this trajectory. Using the definitions Eq. (3.3) of the subtracted lattice observables Δ and Δ¯, we obtain Here, the endpoints of the integral correspond to the temperatures , and tuning the quark masses along the LCP resulted in the factor ( ph )/ .
The expression (3.21) gives the difference between the dimensionless pressure differences Δ / 4 at two distinct temperatures for a given Φ. To determine the change in the pressure at one temperature, for each such Φ additional information is necessary, which corresponds to fixing an integration constant. In the conventional integral method [14] at = 0, one exploits the fact that / 4 vanishes at zero temperature, therefore the integration constant at = 0 is zero. Here, this method is not applicable, since in the presence of a magnetic field, the zero-temperature pressure is no longer zero, see Eq. (3.15). Instead, we propose to use a different region of the parameter space to fix the integration constant, namely the = ∞ line, which corresponds to pure gauge theory plus free static quarks. Since the external magnetic field couples only to quarks, in pure gauge theory has by definition no effect, and Δ (Φ, ) is given solely by the contribution Δ free (Φ, ) of free heavy quarks, which is naively expected to vanish for any and any finite Φ in the limit 2 ≫ . However, in the continuum theory the bare Δ free for static quarks contains the ultraviolet divergent term ∝ free 1 Φ 2 (higher orders in Φ vanish in the static limit), see Eq. (3.10). Therefore Δ free only vanishes in the infinite mass limit after the renormalization has been carried out. Nevertheless, in the lattice regularization Δ free is suppressed as 1/( ) 4 once the quark mass exceeds the lattice scale 1/ , see App. A and the discussion in Sec. 4.1 below. Thus we conclude that at finite lattice spacings, Δ vanishes in the asymptotic quark mass limit. Therefore, integrating down to the physical quark masses ph at fixed , we obtain for an arbitrary temperature Thus, the pressure difference is expressed as an integral of Δ¯over all higher-than-physical quark masses. In practice we first integrate over the two light quark masses up to the point where all three masses coincide (the = 3 theory with different quark charges). Second we integrate over the quark masses simultaneously 4 up to = ∞. The integrand for the up quark is shown in Fig. 1 as a function of the light lattice quark mass for three values of the magnetic field, as measured on the = 6 lattices. At = 113 MeV (left panel of the figure), the difference Δ¯is positive, reflecting the well-known magnetic catalysis of the condensate at low temperatures, see, e.g., Refs. [34,39]. As the mass is increased and the quark decouples, this difference eventually approaches zero. In the right panel of Fig. 1 the same observable is shown, but in the high-temperature phase, at = 189 MeV. At the physical light quark mass the difference is close to zero (see also the results presented in Ref. [39]), and as the mass is increased a peak-like structure is revealed. This structure is a consequence of the strong dependence of the transition temperature on the light quark mass: Around chiral symmetry is restored, the condensate is strongly suppressed and magnetic catalysis is not effective anymore. While at the physical point ≈ 150 MeV [40,41], in pure gauge theory ≈ 260 MeV [42]. At = 189 MeV we start in the chirally restored phase, but as the masses are increased, at some point the transition line is crossed and we enter the chirally broken phase where magnetic catalysis is dominant and Δ¯is large. Eventually, for → ∞ the difference again approaches zero. We note that the down quark condensate shows a very similar behavior for both temperatures = 113 MeV and = 189 MeV.
Through Eq. (3.22) we have determined the integration constant. Now, complemented by Eq. (3.21), the change of the pressure Δ (Φ, ) can be determined at any temperature and at any magnetic flux Φ. Next, the renormalization is performed according to Eq. (3.19) to obtain the renormalized pressure Δ , (Φ, ). The resulting curves are interpolated to compute Δ , ( , ) for any and . This is finally shifted by the zero-field pressure, which we take from Ref. [43], to obtain the full pressure for a range of magnetic fields and temperatures. From the longitudinal pressure, all other thermodynamic observables can be calculated using the relations of Sec. 2.

Lattice ensembles
Before presenting the results, we briefly describe the lattice ensembles we used. These consist of two sets of lattice configurations: one at high temperatures, necessary for the determination of thedependence of the EoS, and one at effectively zero temperature, necessary for the renormalization. The high-ensemble contains = 6, 8 and 10 lattices with various values of the inverse gauge coupling , such that the temperature range 113 MeV < < 300 MeV can be scanned and a continuum estimate can be given (note that at a fixed temperature, the lattice spacing is proportional to 1/ such that the continuum limit corresponds to → ∞, see Eq. (3.1)). These configurations correspond to physical quark masses, tuned along the LCP as discussed at the beginning of Sec. 3. This ensemble was mainly generated in Ref. [16] for the study of the phase diagram and is supplemented in the present analysis by configurations at = 250 MeV and = 300 MeV.
Based on detailed comparisons to our zero-temperature 24 3 × 32 ensembles (see Sec. 4.2), it turned out that the renormalization factors can be determined reliably at our lowest 'finite-temperature' point, = 113 MeV. At this temperature we included two additional lattice spacings with = 12 and 16, allowing for a determination of the renormalization factors down to small lattice spacings, and a matching with perturbation theory. For each we generated configurations ranging from = ph up to = 1200 · ph . The simulation parameters are listed in Table 1.

Condensates, the -function and a comment on magnetic catalysis
We start the presentation of our results with additional details on and implications of the generalized integral method. Let us consider the integrand on the right hand side of Eq. (3.22), and expand it in powers of . For asymptotically large quark masses, quarks and gluons decouple, and Δā pproaches its free theory value. In this limit, we obtain from Eq. (3.11), where is the projector defined in Eq. (3.20). In Fig. 2 we plot [ Δ¯] normalized by ( ) 2 as a function of the light quark mass at = 113 MeV. We perform a combined continuum extrapolation of all five lattice spacings ( = 6, 8, 10, 12 and 16), assuming ( 2 ) discretization errors. The plotted combination has a finite continuum limit (see the discussion after Eq. (3.16)), but discretization errors become large when the lattice quark mass approaches unity. For asymptotically large masses the change of the condensate is proportional to ( ) −5 , as shown in App. A. The finer the lattice, the later this lattice artefact sets in, and the larger are the quark masses that can be reached. With our present lattice spacings we can control the continuum extrapolation up to / ph ≈ 100 − 200. Here the extrapolated values are consistent with the free theory prediction, Eq. (4.1), showing that in the continuum limit Δ¯∝ 1/ for large masses ≫ Λ QCD , to quadratic order in . This implies that the (( ) 2 ) contribution to the integral on the right hand side of Eq. (3.22) diverges logarithmically. On the lattice this divergence is regulated by the inverse lattice spacing such that the cutoff ≈ 1 plays the role of the upper limit of the mass-integral, see the sharp drop in Fig. 2 for asymptotically large masses −1 . In the continuum limit the logarithmic divergence reappears (cf. Eq. (3.15)), and has to be subtracted via charge renormalization. Let us now determine the chiral limit of the combination shown in Fig. 2 using chiral perturbation theory ( PT). The lowest excitation in this case is the charged pion, implying that to leading order Λ H = in Eqs. (3.15) and (3.16). Moreover, due to the spin-zero nature of the pion, the scalar QED -function appears instead of the spinor -function. Using the Gell-Mann-Oakes-Renner relation 2 2 =¯(0)( + ) we obtain for the light flavors = , , as was already pointed out within the Hadron Resonance Gas (HRG) model [32]. To derive Eq. (4.2), we considered equal masses for the light flavors = . Note that the first relation in Eq. (4.2) is understood to hold at = 0, where the charged and neutral pion masses are equal. We indicate the PT prediction in Fig. 2 by the dashed-dotted line, showing a good agreement with the lattice data at physical quark masses, see also the comparison in Ref. [39].
Altogether we observe that the QCD quark condensate, as determined in a fully non-perturbative treatment, interpolates between the PT prediction at small masses and the free-theory limit at large masses. We have seen in Eq. (3.16) that the sign of the condensate Δ¯equals the product of the sign of the factor Λ H / and that of 1 . While for large quark masses ≫ Λ QCD one has Λ H = , towards the chiral limit Λ H = such that the first factor is in both cases positive. It would be quite unexpected to have an intermediate mass where this factor turned negative, nevertheless, we did not find any strict proof of this positivity. In any case we conclude that the (( ) 2 ) magnetic catalysis of the quark condensate and the positivity of the scalar/spinor QED -function are intimately related phenomena. This picture was first described in Ref. [32] and also discussed in Ref. [27]. Note that the above argument concerns the (( ) 2 ) behavior of the condensate and does not address what happens in the large-limit, where the dimensional reduction is expected to be the dominant driving force of magnetic catalysis [34].

Quadratic contribution to the EoS
Next we calculate the (( ) 2 ) contribution to the longitudinal pressure at effectively zero temperature. This term will be subtracted through charge renormalization, according to Eq. (3.19). We perform the integral Eq. (3.22) to determine Δ at various values of the magnetic flux (see the left panel of Fig. 3). To extract the (( ) 2 ) contribution, we fit the data to a quadratic function in ( ) 2 and take the → 0 limit of Δ /( ) 2 , represented by the bars at = 0 in the figure. We perform this analysis for the 24 3 ×6 ensemble at our lowest finite-temperature point, = 113 MeV, and on the 24 3 ×32 ensemble, which corresponds to = 0. The → 0 limits at these two temperatures are found to coincide within our statistical errors. This implies that thermal contributions to the (( ) 2 ) pressure are still strongly suppressed at = 113 MeV, in agreement with our previous findings from the anisotropy method [18], with the results of Ref. [24] using the half-half method and with those of Ref. [26] using the finite difference method (for a brief description of these approaches, see Sec. 3.1). Therefore, we conclude that within our present statistics, it is safe to use = 113 MeV as the reference temperature for the quadratic subtraction. This observation allows us to omit expensive lattice simulations at = 0, and to substitute these by cheaper runs at a finite but low temperature.
In addition, we also check the → 0 extrapolation of our results by directly determining the second derivative of Δ with respect to the magnetic field at = 0 using the half-half method at = 0. For this measurement we employ the setup of Ref. [24] and use 400 noisy estimators to calculate the trace of the necessary operators (for details see Ref. [24]). The result is indicated by the black point in the left panel of Fig. 4, showing that the two methods agree perfectly. Note that in the generalized integral method we use results at all > 0 to extract the quadratic part, resulting in smaller errors. Increasing instead the statistics at = 0 for the half-half method, we could also improve the signalto-noise ratio of the latter, however finite volume effects should also be studied carefully in this case. Since the results at > 0 will be in any case necessary for calculating the higher order contributions to the renormalized pressure, it is advantageous to use the generalized integral method to extract the quadratic term as well.
We proceed by discussing the dependence of [Δ ] on . First we perform the → 0 extrapolation separately for each of our five lattice spacings. The resulting values are expected to lie on the curve 1 ( )·log(Λ H ), see Eq. (3.15). We consider the universal one-loop QCD corrections to 1 ( ) -i.e. terms up to = 1 in Eq. (3.14). The strong coupling in the lattice scheme is defined as 2 (1/ ) = 6/ ( ). The so obtained function 1 ( ) · log(Λ H ) is fitted to the data with Λ H considered as a free parameter. The result is indicated by the orange band in the right panel of Fig. 3. For comparison we also carry out a similar fit in the free case, which corresponds to a simple linear fit with fixed slope free 1 = 0.0169 (which is obtained using Eq. (3.9) for three flavors). The two fits agree within errors for the whole range, indicating that for our lattice spacings, the QCD corrections to the QED -function in Eq. (3.14) are smaller than our statistical errors.
Exploiting the fact that the free scaling describes our data to a very good accuracy, we consider an alternative strategy to determine the quadratic contribution to Δ . Namely, we fit the results for Δ /( ) 2 at all magnetic fields and lattice spacings together using a combined inter/extrapolation in and . We consider the fit function which takes into account the logarithmic divergence of the constant term (with the free scaling coefficient). Here 0 = 1.47 GeV −1 is our largest lattice spacing (corresponding to = 6). We found it necessary to include ( 4 ) lattice discretization effects in the ( ) 2 part and ( 2 ) terms in the ( ) 4 part of the fit function. The results of this combined fit are shown in the left panel of Fig. 4. Considering higher orders in the magnetic field or more lattice artefacts in Eq. (4.3) did not improve the quality of   Table 2. As a cross-check we carried out a similar fit with 1 as a free parameter, resulting in a value consistent with the expected continuum value free 1 . Matching Eq. The scale Λ H depends on the regularization scheme (i.e. on the lattice action). However, towards the chiral limit it is expected to approach a hadronic scale, Λ PT H = (see the discussion in Sec. 4.1), so that this scheme-dependence should only be mild. We stress that Λ H is no free parameter but is automatically determined by the lattice implementation of the renormalization prescription Eq. (3.19). Λ H will appear below in the perturbative description of the pressure as an input from the lattice side.  Subtracting terms quadratic in gives the renormalized pressure Δ , = (1 − )[Δ ] at = 113 MeV. Its continuum limit is given by the → 0 limit of our fit function Eq. (4.3) and is shown in the right panel of Fig. 4. Note that the positivity of Δ , indicates the response to be paramagnetic. Besides the curve for the physical quark masses, we also include here the results obtained for heavierthan-physical quark masses. The analysis is the same in this case except for the lower endpoint of the integral in Eq. (3.22) which we set to = 15 · ph . These light quark masses correspond to a pion mass of about 500 MeV, and the fit according to Eq. (4.3) gives Λ H (15 ph ) = 159(10) MeV. The plot shows that the pressure is clearly less sensitive on as quarks become heavier. In addition we also indicate the Hadron Resonance Gas (HRG) model prediction [32] for the physical pion mass. We will get back to the visible discrepancy between the lattice results and the model curve in Sec. 4.3 below.

Complete magnetic field dependence of the EoS
With the quadratic contribution determined as a function of the lattice spacing, we can carry out the renormalization of the pressure according to Eq. (3.19), and that of the interaction measure using Eq. (3.18) for arbitrary temperatures. Using these two renormalized observables, all other EoS-related quantities can be calculated via the thermodynamic relations of Sec. 2. At = 0 the usual normalization of, e.g., the pressure is / 4 . In our case this may not be the optimal choice, since contains terms of (( ) 4 ) at zero temperature, which give rise to a ∼ 1/ 4 divergence towards = 0. We demonstrate this in Fig. 5 for the case of the longitudinal pressure (left panel) and the interaction measure in the Φ-scheme (right panel). It is therefore instructive to plot the observables without this normalization with respect to 4 . First we show the change in the EoS induced by the magnetic field, Δ , and Δ (Φ) for two values of , see Fig. 6. We find that in the low-temperature region our three lattice spacings do not suffice to perform a controlled continuum   We proceed by discussing the full EoS at nonzero magnetic fields and concentrate on the lowtemperature region where the HRG model is expected to be a valid description. In Fig. 7 we show the longitudinal pressure and the ratio of pressure over total energy density. Comparing , with the HRG model reveals that the model overestimates the pressure at nonzero magnetic fields. This mismatch between the model and the lattice results was not yet visible in our previous comparison [27], where only the = 6 lattice data were available. The ratio , / total (see Eq. (2.1) for our definition of total ) exhibits a shallow dip in the transition region, which moves towards the left, signaling the reduction of the transition temperature as grows, in accordance with our earlier findings with other thermodynamic observables [16,17]. This dip is, however, not pronounced enough to enable us to reliably determine the position of its minimum. We will study the dependence of the transition temperature on in Sec. 4.6.
Up to now we only discussed the longitudinal pressure. Depending on which scheme is used, the magnetic field can induce an anisotropy between the pressure components, see Eq. (2.4). The Φ-scheme describes a situation, in which the transverse compression of the system proceeds at fixed magnetic flux, thereby inducing an anisotropy that is proportional to the magnetization. This scheme is adequate for systems where the magnetic field lines are frozen in the medium -i.e. where the electric conductivity is infinite. The parallel and perpendicular components of the Φ-scheme pressure are shown in the left panel of Fig. 8. The splitting between the components grows as increases, due to the logarithmic rise of the magnetization, cf. Sec. 4.4 below. Note that for large magnetic fields the transverse pressure components become negative. This is due to the positivity of ℳ, which implies that the free energy decreases with growing , i.e. the system prefers large magnetic fields. For fixed Φ, this preference leads to a collapse of the medium in the transverse directions and is signalled by the negative transverse pressure. This unphysical instability invalidates the Φ-scheme for large magnetic fields, and is avoided if a finite (physical) electric conductivity is considered, such that magnetic field lines are not completely frozen and magnetic flux is not completely conserved. We emphasize again that the notion of transverse pressure depends on its precise definition (i.e. the scheme) and should be specified for each problem in question. Our lattice results for and for ℳ are reliable for the whole magnetic field range under study, and can be combined to obtain the transversal pressures in an arbitrary "general" scheme according to Eq. (2.6).

Figure 8.
Left panel: splitting of the Φ-scheme pressure components due to the magnetic field for a few temperatures. The upper branches correspond to (Φ) , , whereas the lower ones to , . Right panel: consistency check at = 0 (see the text).
Since the above described analysis to obtain the equation of state involves several interpolations, we perform one additional consistency check at = 0. Here the longitudinal pressure and the energy density coincide (the -and -directions are indistinguishable even in the presence of the magnetic field, i.e. = − , = ), which gives a relation between the trace of the energy-momentum tensor (the interaction measure in the Φ-scheme) and the pressures. For the renormalized quantities this reads The left hand side of this relation can be obtained without new inter/extrapolations. The renormalization involves subtracting 4[ 0 + free 1 log( / 0 )] · ( ) 2 for the pressure part and free 1 · ( ) 2 for the interaction measure part, respectively (the necessary parameter values are listed in Table 2). The right hand side is obtained by interpolating and differentiating the renormalized longitudinal pressure (as was already done to obtain ℳ of Fig. 6 at > 0), cf. Eqs. (2.2) and (2.5). The two sides are compared in the right panel of Fig. 8 for a zero-temperature lattice at = 3.55 (corresponding to = 1.09 GeV −1 ), showing nice agreement. Note that the relation Eq. (4.5) is only valid at vanishing temperature and is subject to corrections as grows.

Magnetic susceptibility and permeability
The low-behavior of the magnetization provides the magnetic susceptibility , as defined in Eq. (2.11). In the left panel of Fig. 9 we show our results for , compared to those obtained with various other methods (see the discussion in Sec. 3.1). The new results agree within errors with our previous results employing the anisotropy method [18] and with the results of Ref. [26] using the finite difference method. In these studies the same lattice action (stout improved staggered quarks with physical quark masses, / = 28.15) were employed, and the data in each case correspond to a similar continuum estimate as we discussed above. We stress that the three approaches are completely different, but nevertheless show excellent agreement. We also include the susceptibility obtained in Ref. [24], where the HISQ action with nearly physical quark masses ( / = 20) was used and the half-half method was employed on = 8 lattices. The susceptibility is observed to be somewhat smaller than in the other approaches, which may be related to the larger value of the quark mass. In summary, all lattice results indicate that the susceptibility is positive and increases as grows, thus signalling the paramagnetic nature of the thermal QCD medium for temperatures around and above the transition region. The susceptibility can also be calculated using the HRG model. Here the EoS is represented by a sum over contributions of non-interacting hadrons and resonances. At low temperatures and magnetic fields, the dominant term in the sum is given by the lightest hadrons, i.e. pions. The pionic contribution to is negative (see the discussion in App. B), suggesting the presence of a (weakly) diamagnetic region at low . Contrary to pions, higher-spin hadrons give positive contributions to the susceptibility, such that is expected to bend back towards positive values as grows. The HRG prediction is plotted in the left panel of Fig. 9, indeed exhibiting a negative region at low temperatures. Note that at (( ) 4 ), the HRG magnetization receives contributions from the = 0 vacuum term as well, which are positive both for pions and for spin-1/2 and spin-1 hadrons. Eventually, this results in a positive magnetization for magnetic fields exceeding ≈ 0.05 GeV 2 for all temperatures [32], in agreement with our lattice results in the lower left panel of Fig. 6.
In order to extend the temperature range of the lattice results, we include two more sets of 28 3 × 10 lattice ensembles at = 90 MeV ( = 3.55) and at = 105 MeV ( = 3.6). We make use of Eq. (3.21) to obtain the magnetic field-dependence of the pressure, which directly gives the susceptibility. The result is indeed found to be consistent with the HRG prediction, see the red squares in the left panel of Fig. 9. At = 90 MeV we also measure using the half-half method, and use the 24 3 × 32 ensemble for renormalization. Consistent with the above lattice determination, we find ( = 90 MeV) = −0.002 (2). Altogether, the lattice results are compatible with the existence of a weakly diamagnetic region at low temperatures. Our conclusion from this comparison is that the HRG model correctly captures the diamagnetic effect of pions, but overestimates their role for 120 MeV. Note that a negative magnetization has also been obtained within the Parton-Hadron-String Dynamics approach [44]. This model, however, predicts that the diamagnetic response persists to higher magnetic fields ( ≈ 0.2 − 0.7 GeV 2 ) and is even enhanced as the temperature grows. This is in conflict with both the lattice results and the HRG model description presented above. The magnetic susceptibility has also been calculated in a model with free quarks coupled to the Polyakov loop, giving results consistent with the lattice data at high temperatures [45].
The susceptibility can also be translated to the magnetic permeability of the thermal QCD medium. To write down the relation between and , we need to distinguish between the magnetic induction ind and the external field ext that would be present in the absence of the medium. The two fields are connected by the magnetization [19], (4.6) In the present study the magnetic field corresponds to ind , since it is the field that traverses the lattice and that quarks couple to, such that ℳ( ind ≈ 0) = · ind . Using this, the external field ext can be found from Eq. (4.6). Reinserting the factors of = √ 4 em in the definitions of and of ℳ, the magnetic permeability reads In the SI system this is the magnetic permeability in units of the vacuum permeability 0 , cf. Ref. [25]. In the right panel of Fig. 9 we plot / 0 , and compare it to perturbation theory. We discuss some details of this perturbative expansion in the following. It turns out that -even for the lowest-order perturbative expansion of the susceptibility -a non-perturbative parameter is necessary to carry out the comparison with the lattice results in a consistent manner. This parameter is the scale Λ H , which plays the role of the renormalization scale, = Λ H , see the discussion of Sec. 3.2.2. The necessity of using Λ H in this comparison is related to the entanglement of ultraviolet and infrared divergences in the presence of the magnetic field: the behavior of the bare free energy at → 0 is identical 6 to that of the renormalized free energy at → ∞. Thus, the high-temperature susceptibility is again governed by the QED -function [18,46,47], ( ) = 2 · 1 · log( /Λ H ). (4.8) and inserting the value Λ H = 0.12 GeV that we determined in Sec. 4.2 for physical quark masses, we obtain the blue dashed curve in the right panel of Fig. 9. To improve this perturbative expression, we also take into account QCD corrections to free 1 using Eq. (3.14), but this time at the thermal scale th ∼ instead of at the lattice regulator 1/ . To calculate 2 ( th ) we use the four-loop running coupling and Λ MS QCD = 0.34 GeV for three-flavor QCD [48]. Considering QCD corrections up to various orders in = 2 /(4 ) we obtain the dashed bands in the figure. The width of the bands correspond to the uncertainty of the thermal scale, which we allow to vary between and 4 . This range is physically motivated based on convergence arguments [49]. The perturbative expansion 7 seems to show a fast convergence even for reasonably low temperatures, and agrees nicely with the lattice data in the temperature region 200 MeV < < 300 MeV.
It is important to stress that Λ H appeared in the perturbative description due to its role as the 6 This may be understood as follows. As the temperature increases, gradually adopts the role of the largest scale in the system and replaces the scale ΛH in Eq. (3.15). The renormalization prescription, however, remains unchanged and still amounts to the subtraction at a renormalization scale = ΛH, as contained in Eq. (3.15). Altogether one obtains, to quadratic order, a logarithmic term log( /ΛH), from which Eq. (4.8) follows directly. 7 We mention that to obtain the perturbative improvement of Eq. (4.8) we only considered the QCD corrections in the coefficient of the leading logarithm and ignored possible corrections that arise in the sub-leading constant terms (i.e. those that could modify the dependence on the renormalization scale ΛH). These could be accounted for by considering the scale th (proportional to ) instead of in Eq. (4.8). Note moreover that both the higher-order corrections to the -function and ΛH depend on the renormalization scheme, but this dependence must cancel in the susceptibility, being a physical observable. renormalization scale. Since the susceptibility contains an ultraviolet divergence, its renormalization inevitably introduces an ambiguity, expressed as a dependence on the renormalization scale. To derive Eq. (4.8) we relied on two observations: that the zero-temperature free energy to (( ) 2 ) is determined exclusively by the QED -function and that at high temperatures replaces the regulator 1/ in this expression. Eq. (4.8) can be explicitly checked in the free case [18,46,47]. Note that our arguments about charge renormalization only relate to the (( ) 2 ) contributions, whereas the full, -dependent free energy contains much more information and is also considerably more complicated. Its calculation to ( ) was performed recently using the lowest-Landau-level approximation, valid for large magnetic fields [50].

Entropy density and the Adler function
As emphasized in Sec. 4.3, most observables contain magnetic field-induced contributions at = 0, making the usual normalization (e.g., with respect to 4 for the pressure) disadvantageous. In this respect, the entropy density is special: being the derivative of , with respect to the temperature, it vanishes identically at = 0. This may be understood from the fact that the vacuum contribution is a pure quantum effect (it emerges from the interaction of virtual quarks with the external field) and thus it cannot produce entropy. Note that at > 0, the magnetic field changes the thermal distribution and is expected to modify . The vanishing of ( = 0) allows for a study of the usual normalization / 3 , plotted in the left panel of Fig. 10. Note that the errors become large at low temperatures as increases, due to the cancellation of the ( −3 ) magnetic field-induced contributions in . For the lowest three magnetic fields at low temperatures, we also show the HRG prediction in the plot. Moreover, we indicate the Stefan-Boltzmann limits for each value of . These can be calculated from the dependence of the free pressure on and , where we considered three massless quark flavors and used the high-temperature limit of the magnetic susceptibility from Eq. (4.8). Note that the renormalization scale Λ H cancels in the entropy density. In the right panel of Fig. 10 we show the magnetic field-induced part Δ / 3 , compared to the expected perturbative behavior. Note that while ( = 0) is by 20−30% below its Stefan-Boltzmann limit at our highest temperature, Δ almost perfectly follows the free-case prediction, already for 170 MeV. A similar behavior is observed for the pressure and, thus, for all other observables as well. In other words, the -dependence of the EoS in the deconfined phase is predominantly dictated by the ( 2 ) free-case behavior -in sharp contrast to the = 0 EoS, which shows strong deviations from the perturbative predictions within this range of temperatures. Similar conclusions were drawn using a perturbative treatment of QCD in magnetic fields, where the ( 2 ) term was shown to be suppressed with respect to the free-case contribution [50].
The correspondence between the entropy density and perturbative QED physics can be pushed even further: we find that the second derivative of the entropy density with respect to the magnetic field at = 0 is related to the Adler function (for its definition, see, e.g., Ref. [51]) at high temperatures. To understand the origin of this relation, it is advantageous to consider a (perturbative) diagrammatic representation of the free energy ( ): it consists of all closed loop diagrams containing virtual quarks and Figure 11. Second derivative of the entropy density with respect to at high temperatures, compared to perturbation theory.
gluons. On the one hand, taking the second derivative with respect to gives the susceptibility . On the other hand, this derivative effectively pulls out two photon legs (these photons correspond to the background magnetic field), giving the photon vacuum polarization diagram Π. This diagram is usually considered with an inflowing momentum , and the Adler function is defined from it via ( ) = 12 2 · Π/ log 2 . The equivalent of the momentum in our setup is expected to be the largest scale in the system: at high temperatures, for example, ∼ (remember that to define , the magnetic field has already been set to zero). Thus, in this region we expect that the Adler function is approached as where we interchanged the derivatives with respect to and to . Eq. (4.10) reveals a relation between the magnetic field-dependence of the entropy and the Adler function at a thermal scale th ∼ . In Fig. 11 we show our continuum estimate for the right-hand-side of Eq. (4.10) and a comparison to the perturbative expansion [51] of ( th ) (where we used th = 2 ). Note that the above correspondence is only expected to be valid for ≫ Λ H , where the relevant scale is uniquely defined by the temperature. Note also that according to this argument, Eq. (4.10) fixes the asymptotic dependence of 2 / ( ) 2 on any external parameter (e.g., on chemical potentials), as long as this parameter represents the largest scale in the system. The correspondence between the susceptibility and the Adler function clearly deserves a more detailed investigation, which we plan to conduct in the near future.

Phase diagram
Let us now use the dependence of the EoS on and on to discuss the QCD phase diagram in the − plane. To this end we need to define ( ) through characteristic points of some observables.
We have seen that most observables are nonzero at = 0, making the usual normalization by 4 disadvantageous. The shift in , at zero temperature is of (( ) 4 ) and positive, as we have seen in Sec. 4.3. It is simple to show using this observation and the thermodynamic relations of Sec. 2 that · ℳ, (Φ) and total are also of (( ) 4 ) and positive at = 0. Conversely, and ( ) are of (( ) 4 ) and negative at zero temperature. A special combination is total − 3 , = (Φ) − · ℳ for which quartic terms also cancel and which is of (( ) 6 ). Let us consider this combination in somewhat more detail. It is plotted in the left panel of Fig. 12. This combination coincides with the renormalized interaction measure at = 0. Accordingly, it exhibits a pronounced peak, which is observed to move towards lower temperatures as increases. Eventually, the (( ) 6 ) terms start to contribute to this combination, making it diverge at low temperatures and washing out the peak-like structure. Still, up to = 0.4 GeV 2 we can use the peak of this observable to characterize the transition region. For a first-order phase transition, the inflection point of this observable would turn into a discontinuity, thus this point marks the (pseudo-critical) transition temperature ( ) for the crossover. Similarly, the inflection point of / 3 (see the left panel of Fig. 10) also represents a candidate for defining the transition temperature. Together with the maximum of ( total − 3 , )/ 4 (which, however, does not correspond to a pseudo-critical temperature, but is merely a characteristic point), we show the -dependence of these definitions in the right panel of Fig. 12, and compare them to our earlier determinations of the phase diagram using the strange quark number susceptibility and the light quark condensates. The results consistently show a reduction of the transition temperature as grows. Note that the difference between different determinations of reflects the crossover nature of the transition. The variance between the four definitions is found to remain constant within the errors.  [16] for the inflection point of the strange quark number susceptibility (light blue shaded band) and that of the light quark condensates (light red shaded band).

Summary
Using a novel 'generalized integral method', we determined the QCD equation of state for a wide range of temperatures and magnetic fields up to = 0.7 GeV 2 . Results were presented for a variety of thermodynamic observables, indicating that the EoS is significantly affected by the magnetic field, even at moderate values of . The tabulated data is available in the ancillary files table_EoS_B.dat and table_EoS_B_derivatives.dat, submitted to the arXiv along with this paper.
The thermodynamic structure of QCD is altered by the magnetic field in several aspects: • Vacuum term. The magnetic field induces a vacuum contribution to most observables, such that e.g. the pressure does not vanish at = 0. This vacuum term makes the usual normalization of the affected observables (e.g. / 4 ) ill-defined at low temperatures.
• Pressure anisotropy. The magnetic field creates an anisotropy between the pressure components if the transverse pressure is defined at constant magnetic flux (Φ-scheme). This anisotropy is sizeable and becomes comparable to the longitudinal pressure as or increase.
• Para-/diamagnetism. The leading-order response of the system to is characterized by the magnetic susceptibility . On the one hand, the thermal QCD medium is a strong paramagnet ( > 0) around and above the transition region. On the other hand, there appears to be a weakly diamagnetic region ( < 0) at low temperatures 100 MeV, where pions dominate.
• Validity of HRG and of PT. For the susceptibility, the Hadron Resonance Gas model breaks down already at ≈ 120 MeV. However, perturbation theory successfully describes the lattice data at suprisingly low temperatures ≈ 200 − 300 MeV.
The presence of the background magnetic field necessitates the renormalization of the electric charge. This has several implications: • Renormalization. The pressure undergoes additive renormalization at = 0. The divergent term that needs to be subtracted is logarithmic in the lattice spacing and its coefficient equals the lowest-order QED -function (with QCD corrections at the scale 1/ ).
• Magnetic catalysis. At zero temperature, the phenomenon of magnetic catalysis (the enhancement of the quark condensate by the magnetic field) to quadratic order in is related to the positivity of the QED -function.
• Susceptibility. For high temperatures, the magnetic susceptibility increases logarithmically with , at a rate given by the QED -function (with QCD corrections at the scale ).
• Adler function. The second derivative of the entropy density with respect to the magnetic field is related to the perturbative Adler function at high temperatures.
In addition, we considered characteristic points of a few observables to explore the QCD phase diagram in the − plane. This analysis indicates that the transition temperature decreases as grows, in agreement with our previous results where other thermodynamic observables (light quark condensate, strange quark number susceptibility and Polyakov loop) were used [16,17]. Lattice results indicating this tendency have also been obtained using overlap fermions in = 2 QCD [52] and in two-color QCD with four equally charged staggered quark flavors [53]. The reduction of ( ) has been reproduced within the bag model [37] and is also supported by large arguments [54]. Nevertheless, this feature remains a property of the chiral/deconfinement transition that most low-energy effective theories or models cannot reproduce, or only for a limited range of magnetic fields, see, e.g., Refs. [55][56][57]. Recent studies of the Nambu-Jona-Lasinio model, however, indicate that taking into account a -dependent Polyakov loop scale parameter [58], or the magnetic field-induced running of the strong coupling [59][60][61] might resolve this discrepancy.

B Magnetic susceptibility in the HRG model
In this appendix we calculate within the Hadron Resonance Gas model. The EoS at nonzero magnetic fields was determined in Ref. [32] by writing the free energy density as an integral over the longitudinal momentum and a sum over Landau-levels, both of which can be performed numerically. To obtain , it is instead advantageous to use the proper-time representation [28], where the expansion of in the magnetic field can be written down directly. For hadrons with electric charge and spin = 0, 1/2 or 1, the energy levels in the magnetic field read ( , , ) = √︀ 2 + 2 + (2 + 1 − 2 ), where is the projection of the spin on the magnetic field, and we approximated the gyromagnetic factor of the hadron as = 2. To calculate the renormalized susceptibility, it suffices to determine the difference of free energies at and at = 0, where the prefactor (−1) 2 +1 reflects the fermionic/bosonic nature of the hadron, the elliptic Θ-function results from summing over Matsubara-frequencies and the factor 2 sinh( ) from summing over the angular momenta in Eq. (B.1). The first argument of the Θ-function is = 0 for = 0, 1 and = /2 for = 1/2, according to the lowest Matsubara-mode. The −1 in the square brackets corresponds to the subtraction of the = 0 term. Note that Eq. (B.2) gives the contribution of a particle and its antiparticle to the free energy density.
Keeping quadratic terms in the magnetic field gives the susceptibility, and made a change in the integration variable. The remaining integral over can be performed numerically. Inspecting the behavior of the Θ-function we see that ( ) is negative for = 0 and positive for = 1/2, 1: charged pions contribute to diamagnetism (this was also recognized in Ref. [25]), whereas protons and charged -mesons to paramagnetism. This tendency can be qualitatively understood invoking the following argument (see the discussion in Ref. [32]). The susceptibility is determined by the thermal part of the free energy, which contains exp(− ( )/ ) where ( ) is the effective mass of the hadron at non-zero magnetic fields. According to the structure of the lowest Landau-level (Eq. (B.1) with = = 0 and = ), for scalar (vector) hadrons this effective mass increases (decreases) as grows. Therefore, the magnitude of the thermal free energy is suppressed by the magnetic field for pions, whereas it is enhanced for -mesons, responsible for the different signs of in the two cases. In the = 1/2 case, the lowest level is independent of , thus the sign of the susceptibility cannot be anticipated from this argument. As already mentioned in Sec. 4.4, at ( 4 ) the vacuum part starts to contribute to the free energy as well and turns the magnetization positive for each spin channel, given that 0.05 GeV 2 . The total susceptibility is obtained as the sum over all hadrons ( ) = ∑︁ ℎ ℎ · ℎ ( ), (B.5) with multiplicities ℎ . The list of hadrons taken into account can be found in Ref. [32].