Magnetic susceptibility of QCD matter and its decomposition from the lattice
PPrepared for submission to JHEP
Magnetic susceptibility of QCD matter and itsdecomposition from the lattice
Gunnar. S. Bali, a,b
Gergely Endrődi, c,d and Stefano Piemonte a a Institut für Theoretische Physik, Universität Regensburg, D-93040 Regensburg, Germany. b Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India. c Institut für Theoretische Physik, Goethe Universität Frankfurt, D-60438 Frankfurt am Main, Germany. d Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany.
E-mail: [email protected] , [email protected] , [email protected] Abstract:
We determine the magnetic susceptibility of thermal QCD matter by means of firstprinciples lattice simulations using staggered quarks with physical masses. A novel method is em-ployed that only requires simulations at zero background field, thereby circumventing problems relatedto magnetic flux quantization. After a careful continuum limit extrapolation, diamagnetic behavior(negative susceptibility) is found at low temperatures and strong paramagnetism (positive suscepti-bility) at high temperatures. We revisit the decomposition of the magnetic susceptibility into spin-and orbital angular momentum-related contributions. The spin term – related to the normalizationof the photon lightcone distribution amplitude at zero temperature – is calculated non-perturbativelyand extrapolated to the continuum limit. Having access to both the full magnetic susceptibility andthe spin term, we calculate the orbital angular momentum contribution for the first time. The resultsreveal the opposite of what might be expected based on a free fermion picture. We provide a simpleparametrization of the temperature- and magnetic field-dependence of the QCD equation of statethat can be used in phenomenological studies.
Keywords:
Lattice field theory simulation, Quark Gluon Plasma, QCD Phenomenology, LatticeQCD, Nonperturbative Effects. a r X i v : . [ h e p - l a t ] A p r ontents T > C.1 Magnetic susceptibility 23C.2 Ultraviolet divergences and QED renormalization 24C.3 Spin contribution 24C.4 Tensor bilinear 25C.5 High-temperature expansion 26
D Multiplicative QCD renormalization 27E Susceptibilities via current-current correlators 29
E.1 Magnetic susceptibility from correlators 29E.2 Tensor coefficient from correlators 31
F Parametrization of the equation of state 32 – 1 –
Introduction
The development of a quantitative and precise understanding of the response of QCD matter to back-ground (electro)magnetic fields is of vital importance for furthering our knowledge about a multitudeof physical systems. Examples include the interior of magnetars, neutron star mergers [1–3], off-centralheavy-ion collisions and the evolution of the universe in its early stages. For general reviews, we referthe reader to Refs. [4–6]. A characteristic feature of the behavior of strongly interacting quarks andgluons is rooted in the dependence of the QCD equation of state (EoS) on the background magneticfield B . The EoS enters all the above mentioned examples: it appears in the gravitational stabilityconditions of compact stars, affecting the mass-radius relation [7]; it governs the expansion rate incosmological models [8, 9] and it also sets the conditions where freeze-out is reached in heavy-ioncollisions, see, e.g., Ref. [10].While the equilibration of fireballs produced in heavy-ion collisions is still a subject of research,at least in astrophysical systems the time and distance scales over which the magnetic field variesare much larger than those that govern QCD processes that affect, e.g., the EoS or nucleo-synthesis.With these applications in mind, solving QCD in a constant background magnetic field is sufficientand this scenario is amenable to lattice simulations.The leading dependence of the EoS on B is encoded in the magnetic susceptibility χ of QCDmatter. Its sign distinguishes between paramagnets ( χ > ), for which the exposure to the backgroundfield is energetically favorable, and diamagnets ( χ < ), which repel the external field. In QCD matter,like in any other material, the origin of the magnetization and hence of the magnetic susceptibilityis related to the spin and angular momentum of charged particles. At high temperatures the quarksare the relevant degrees of freedom; at low temperatures the hadrons and in particular the pions takeover, while (valence and sea) quarks contribute just as their fundamental constituents.The total angular momentum that gives rise to the magnetization can be decomposed into con-tributions from the spins of the quarks of different flavors and a remainder. The latter contains thequark orbital angular momenta but also the angular momentum of the gluons, that can split intoquark-antiquark pairs. The quark spin contribution to the magnetic susceptibility is due to theexpectation value (cid:10) ¯ ψ f σ µν ψ f (cid:11) , whose leading order response is linear in the magnetic field strengthtensor F µν . Depending on the normalization, the slope is proportional to the magnetic susceptibil-ity of the quark condensate [14–16], or the so-called tensor coefficient, i.e. the normalization of thephoton lightcone distribution amplitude (DA) [17–19] of finding a quark-antiquark pair of flavor f in a transverse photon. This in itself appears in a multitude of applications, e.g., as a correction tothe hadronic light-by-light scattering contribution to the muon anomalous magnetic moment [20, 21],within radiative transitions [22, 23] and in the photo-production of mesons [24].Some of the present authors have already addressed several of the above aspects in Refs. [25–28].Here we improve on these studies by employing a novel calculational method that is based on Ref. [29], This situation is analogous to the decomposition of the nucleon spin. In particular, vacuum expectation values ofthe same local operators appear in the magnetic field background at zero momentum as in the decomposition [11] ofthe transversely polarized generalized parton distribution functions of deep inelastic scattering at leading twist. Whilethe individual quark spin contributions in both cases are unique and gauge invariant, the further decomposition of theremainder into quark and gluon parts is ambiguous: the decomposition of Ref. [11] is based on the Belinfante-Rosenfeldform of the energy momentum tensor, which is also the natural starting point for lattice QCD, but one may also, e.g.,resort to the canonical definition of the angular momentum [12, 13]. – 2 –y carrying out the QCD renormalization non-perturbatively with respect to the intermediate RI’-MOM scheme [30, 31] and by adding a finer lattice spacing. In addition we present and exploit newanalytical findings. One of the outcomes will be that in the strongly interacting medium at low tomoderately high temperatures the quark spin-related susceptibility is negative (diamagnetism) whilethe part that is due to the orbital angular momentum is positive. Clearly, this behavior is verydifferent from the response to magnetic fields of the materials that have so far been accessible to solidstate physics experiments. Therefore, our results offer a glimpse into a completely new regime of spinphysics.The article is organized as follows. In Sec. 2 we introduce our notations and the central observ-ables. For conceptual clarity, we carefully address their divergence structure and renormalization inQED and QCD. In Sec. 3 we then discuss details of the simulation and, in particular, we introduce ournew method that employs current-current correlators in a mixed coordinate- and momentum-spacerepresentation. This enables us to determine susceptibilities from lattice simulations at B = 0 . Wethen present and discuss our results in Sec. 4, before we summarize. We include several technical ap-pendices: in App. A we investigate the effects of taste splitting in the staggered formulation within thehadron resonance gas model. This turns out to be important to avoid underestimating the systematicsof the continuum limit extrapolation. In App. B we derive the factorization of the susceptibility intoquark spin-related and other contributions, building upon earlier partial results [26]. In the extensiveApp. C, several derivations are carried out for the free case, establishing, e.g., the structure of QEDdivergencies. App. D discusses the non-perturbative renormalization procedure. In App. E we presentmore detail on the derivation of the new current-current method and compare to numerical results,obtained using conventional background field approaches. Finally, App. F gives a parametrization ofour results for the QCD EoS for a broad range of temperatures and magnetic field strengths. Thecorresponding Python script param_EoS.py is uploaded to the arXiv as ancillary file along with thispaper. Without any loss of generality, below we consider a magnetic field pointing in the x direction, withthe magnitude B . The magnetic susceptibility is defined via the leading (quadratic) dependence ofthe QCD free energy density f on the field strength, χ b = − ∂ f∂ ( eB ) (cid:12)(cid:12)(cid:12)(cid:12) B =0 , f = − TV log Z , (2.1)where Z is the partition function, T the temperature and V the spatial volume of the system. Theproduct eB of the elementary electric charge e and the magnetic field is a renormalization group invari-ant due to the QED vector Ward identity, see Eq. (2.3) below. Therefore, χ b is free of multiplicativerenormalization. Still, the susceptibility undergoes additive renormalization, which is made explicitby the index b , indicating the bare quantity, as obtained in the lattice scheme at a lattice spacing a . This was discussed in Ref. [28] in depth, but we repeat the argument here for comprehensiveness. Below we will indicate quantities that are subject to QED renormalization with the subscript b . – 3 – .1 QED renormalization The total free energy density of the system, that includes both QCD matter and the classical back-ground field, f tot = f + B b , (2.2)is a physical observable and therefore free of divergences. However, the second term within Eq. (2.2),involving the bare magnetic field B b , contains a logarithmic divergence in the lattice spacing a dueto electric charge renormalization [32], B b = Z e B , e b = Z − e e , eB = e b B b , Z e = 1 + β ( a − ) e log (cid:0) µ a (cid:1) , (2.3)where the renormalized quantities e and B depend on the QED renormalization scheme and on theQED renormalization scale µ QED . Since the background field is classical, only the leading-order QED β -function coefficient β appears here [33]. Note that β is affected by QCD corrections at the cut-offscale, β ( a − ) = β · (cid:34) (cid:88) i ≥ c i g i ( a − ) (cid:35) a → −−−→ β , β = 14 π · (cid:88) f ( q f /e ) , (2.4)where g is the strong coupling and q f denotes the electric charge of the quark flavor f . The coefficients c i of the perturbative series are known up to i = 4 in the MS and MOM schemes [34] and c = 1 / (4 π ) is universal for massless schemes. These QCD corrections vanish logarithmically with the latticespacing towards the continuum limit due to the asymptotically free nature of the strong interactions,as is also indicated in Eq. (2.4).Eq. (2.2) implies that f contains the same additive divergence as B b / , but with an oppositesign. This propagates into the susceptibility (2.1), resulting in χ b = χ [ µ QED ] + β ( a − ) log( µ a ) . (2.5)The renormalized susceptibility χ depends on the renormalization scale, which we indicated hereexplicitly in square brackets. We confirm the presence of the logarithmic divergence in χ b analyticallyfor the free case in App. C.2 and numerically for full QCD in Sec. 4.1. The divergence is independentof the temperature so that it cancels within the difference χ ≡ χ [ µ physQED ] = χ b ( T ) − χ b ( T = 0) . (2.6)This definition of the renormalized susceptibility – implying that it vanishes identically at T =0 – corresponds to a particular choice of the renormalization scale µ QED = µ physQED . In fact thisis the only prescription that adheres to the physical requirement that the magnetic permeability (1 − e χ ) − should be unity in the vacuum. In the following we suppress the dependence on the QEDrenormalization scale and simply write χ for the susceptibility, renormalized in this way. Note that in general disconnected diagrams also start to contribute for i ≥ , but in the present case these vanishbecause we are dealing with the three lightest quark flavors and (cid:80) f q f /e = 0 . – 4 – .2 The tensor coefficient Besides the magnetic susceptibility there exist further quantities that characterize the leading-orderresponse of the QCD medium to the background magnetic field. For a general background field F µν , the fermion bilinear involving the relativistic spin operator σ µν develops a nonzero expectationvalue [14, 35], (cid:10) ¯ ψ f σ µν ψ f (cid:11) = q f F µν · τ fb + O ( F ) , σ µν = 12 i [ γ µ , γ ν ] . (2.7)We will refer to τ fb as the tensor coefficient for the flavor f . Similarly to χ b , this is also a bareobservable that contains additive logarithmic divergences in the cut-off. For our choice of directionof the magnetic field the tensor coefficient can be determined as the slope of the expectation value ofthe fermion bilinear involving σ at small values of the magnetic field: τ fb = lim B → (cid:10) ¯ ψ f σ ψ f (cid:11) q f B . (2.8)The tensor coefficient contains a similar logarithmic divergence as χ b . We demonstrate the reasonfor this in Sec. 2.3 below. In particular, the divergence structure takes the form [26], τ fb = τ f + γ τ ( a − ) m f log( µ a ) , γ τ ( a − ) = γ τ · (cid:2) O ( g ( a − )) (cid:3) , γ τ = 34 π , (2.9)where m f is the mass of the quark of flavor f . Again, due to asymptotic freedom, QCD correctionsto γ τ vanish in the continuum limit. Eq. (2.9) is confirmed in App. C.4 for the free case and checkednumerically in full QCD in Sec. 4.2. Notice that in the chiral limit the divergent term disappears, sothat lim m f → τ fb is ultraviolet-finite. We will carry out this limit at zero temperature in Sec. 4.2 below.In this situation, up to multiplicative renormalization, this object corresponds to the normalization f ⊥ γ of the leading-twist photon distribution amplitude [17–19], i.e. of the infinite momentum frameprobability amplitude that a real photon dissociates into a quark-antiquark pair of flavor f .In analogy to Eq. (2.6), we define the renormalized tensor coefficient by subtracting its value atzero temperature, τ f = τ fb ( T ) − τ fb ( T = 0) , (2.10)which again corresponds to a particular choice of the QED renormalization scale. Unlike for χ , thereis no preferred choice in this case, but it is natural to use the same prescription as in Eq. (2.6) above.Besides the (QED-related) additive renormalization detailed above, the tensor coefficient alsoundergoes (QCD-related) multiplicative renormalization by the tensor renormalization constant Z T .This introduces a further scheme- and scale-dependence of this observable. Below we will consider Z T in the MS scheme at the QCD renormalization scale µ QCD = 2
GeV. Unlike in our previousstudy [26], where we calculated Z T perturbatively at the one-loop level, here we carry out a non-perturbative matching to to the RI’-MOM scheme [30, 31] and subsequently translate the result atthree-loop order [36] into the MS scheme. This procedure is detailed in App. D. We remark that Z T is independent of the temperature, thus the ordering of the QED renormalization (i.e. the T = 0 subtraction) and the QCD renormalization (multiplication by Z T ) is irrelevant for the determinationof the renormalized tensor coefficient τ f . Unlike for β ( a − ) in Eq. (2.4), the order g perturbative coefficient is not known in this case. Therefore, anydefinition of a renormalized quark mass m f is valid to this order and we use the lattice quark mass. – 5 – .3 Decomposition into spin and orbital angular momentum contributions One might suspect that χ and τ f are not completely unrelated. Indeed, as we first discussed inRef. [26], τ f represents the contribution of the spin of the quark flavor f to the total magneticsusceptibility. In particular, χ can be decomposed into spin-related and orbital angular momentum-related contributions, χ = χ spin + χ ang , (2.11)and the spin term is related to the tensor coefficients as χ spin = (cid:88) f ( q f /e ) m f (cid:104) τ f ( m val f ) − τ f ( m val f → (cid:105) · Z T Z S , (2.12)where m f denotes the quark mass in the lattice scheme and Z S and Z T are the scalar and tensorrenormalization constants, respectively. The second term in the square brackets is understood tocorrespond to the limit of a vanishing valence quark mass taken at physical, i.e. nonzero, values of thesea quark masses. We discuss the difference between valence and sea quark masses in Sec. 3 below.In App. B we prove Eqs. (2.11)–(2.12) and illustrate the origin of the subtraction of the chiral valencequark limit. App. C also contains an explicit check of Eq. (2.12) in the free case.A remark regarding the choice of renormalization scales is in order here. Eq. (2.12) is chosen sothat χ spin vanishes at T = 0 . Recall however that, according to our remark below Eq. (2.10), we arefree to choose an arbitrary QED renormalization scale for the renormalized tensor coefficient. Thisfreedom propagates into χ spin and – through the decomposition (2.11) – to χ ang as well. In contrast,the renormalization scale for χ is fixed by the requirement χ ( T = 0) = 0 . Thus, in principle bothsusceptibility contributions may be shifted by an arbitrary amount, as long as their sum remains zeroat T = 0 . We will follow the choice made in Eq. (2.12), which corresponds to setting χ spin ( T = 0) = χ ang ( T = 0) = 0 . In the free case this is realized by choosing one and the same QED renormalizationscale for all susceptibility contributions, see App. C.3. Our numerical results in full QCD belowsuggest that also in the interacting case the QED scale that corresponds to the renormalizationcondition χ spin ( T = 0) = 0 is consistent with the one obtained from setting χ ( T = 0) = 0 .In Eq. (2.12) we also carried out the QCD related multiplicative renormalization by includingthe tensor renormalization constant Z T required for τ f (as mentioned above) as well as the scalarrenormalization constant Z S = Z − m , which multiplies the inverse quark mass. Note that theserenormalization factors depend on the QCD regularization scheme and on the QCD renormalizationscale. As mentioned above, we choose the MS scheme and µ QCD = 2
GeV. We remark that, just asfor τ f , the ordering of the QED and the QCD renormalization is irrelevant for χ spin . We stress againthat while the factorization of the total susceptibility χ into χ spin and χ ang depends on the QCDscheme and scale, χ itself is a QCD renormalization group invariant.In the free case the two susceptibility contributions have a constant ratio, χ spin : χ ang = 3 : ( − ,reflecting the well-known response of a free charged fermion to the magnetic field via its spin and its Note that for simplicity we refer to χ ang as the orbital angular momentum contribution. In the interacting casethis can be further factorized, separating out the gluon total angular momentum contribution [27] from those of thequark angular momenta. Like in massless continuum schemes, in staggered lattice formulations there is no difference between singlet andnon-singlet renormalization factors for these quark bilinears. This was explicitly demonstrated at order g in Ref. [37]. – 6 –rbital angular momentum, dating back to Pauli and Landau [38, 39]. This ratio, which translatesinto the rule χ spin : χ = 3 : 2 , holds identically in the free case, see App. C.3. In contrast, in full QCDit only applies to the divergence structure, which in the continuum limit – as we have seen above – isgoverned by pure QED physics. Below we determine to what extent the Pauli-Landau decompositionis affected by the strong interactions. We consider spatially symmetric N s × N t lattice ensembles, corresponding to the temperature T =1 / ( N t a ) and the volume V = L = ( N s a ) . The simulations are performed with the tree-levelSymanzik improved gauge action S g and three flavors ( f = u, d, s ) of stout smeared rooted staggeredquarks [40], described by the Dirac operator /D f + m f . The quark masses m f are tuned as a functionof the inverse gauge coupling β = 6 /g along the line of constant physics: m ud ( β ) ≡ m u ( β ) = m d ( β ) = m s ( β ) /R with R = 28 . [41]. The electric charges are set as q d = q s = − q u / − e/ , where e > is the elementary electric charge. The magnetic field enters in /D f via space-dependent U(1) phases.Further details of our setup and of the simulation algorithm are discussed in Ref. [42]. The latticegeometries for our finite temperature lattices are × , × , × , × and × ,allowing for the investigation of both finite volume and discretization effects. Our zero-temperatureensembles consist of × , × and × lattices.In the rooted staggered formulation the partition function and the expectation value of the tensorbilinear are written as path integrals over the SU(3) -valued gluonic links U as Z = (cid:90) D U e − βS g (cid:89) f (cid:48) = u,d,s (cid:2) det( /D f (cid:48) + m sea f (cid:48) ) (cid:3) / , (cid:10) ¯ ψ f σ ψ f (cid:11) = TV Z (cid:90) D U e − βS g tr σ /D f + m val f (cid:89) f (cid:48) = u,d,s (cid:2) det( /D f (cid:48) + m sea f (cid:48) ) (cid:3) / . (3.1)Here we distinguished between two different types of masses: the sea quark masses m sea f which appearin the fermion determinant and thus affect the generation of gluonic configurations; and the valencequark mass m val f , which enters in the operator and thereby affects the measurement on a given set ofconfigurations. For usual observables both masses are equal and set according to the line of constantphysics, m sea f = m val f = m f . The spin contribution to the magnetic susceptibility is exceptionalin this sense – as pointed out above in Eq. (2.12), it also involves the value of (cid:10) ¯ ψ f σ ψ f (cid:11) in thelimit m val f → but keeping m sea f = m f . Note that this does not mean that we are dealing with anon-unitary theory but merely that χ spin can be expressed as a difference of two expectation valuesinvolving τ f at different valence quark mass values, see App. B. In an infinite volume a magnetic field pointing in the x direction can be generated by the Landau-gauge electromagnetic potential A = Bx . (3.2) Due to the electromagnetic charge q f , the covariant derivative depends on the quark flavor f . – 7 –n a finite volume, in order to comply with periodic boundary conditions for the electromagneticparallel transporters u µf = exp( iq f A µ ) , the boundary twist term A = − Bx L δ ( x − L ) needs to beincluded as well [43]. In this setup the flux of the magnetic field is quantized according to [44] eB = 6 πN B /L , N B ∈ Z , (3.3)so that a differentiation with respect to eB – as required in Eq. (2.1) – is not possible in a stan-dard manner. Several methods were developed to overcome this problem on the lattice, includingthe anisotropy method [27, 45], the finite difference method [46, 47] and the generalized integralmethod [28]. These are all based on approximating the derivative numerically using finite differencesin the integer variable N B . This requires independent simulations using several values of N B , whichincreases the computational requirements considerably. Furthermore, an extrapolation N B → be-comes necessary, which inevitably introduces systematic uncertainties. An alternative approach is thehalf-half method [48], which employs a magnetic field profile that is positive in one half and negativein the other half of the lattice – this enables taking the derivative with respect to the amplitude ofthe field analytically. However, finite volume effects are substantially enhanced in this case due tothe discontinuity of the background field at the boundaries, as was demonstrated in Ref. [29]. In view of the above, methods are desirable that only involve measurements at B = 0 , therebycircumventing the flux quantization problem of the constant background field profile. In Ref. [29] weshowed that χ b can be related to correlators of the electromagnetic current: j µ = (cid:88) f q f e ¯ ψ f γ µ ψ f . (3.4)Specifically, for our electromagnetic gauge (3.2) the relevant role is played by the x -dependence ofthe correlator of the µ = 2 component, G ( x ) = (cid:90) d x d x d x (cid:104) j ( x ) j (0) (cid:105) , (3.5)where the zero-momentum projection p = p = p = 0 was performed. Using this projected correla-tor, the magnetic susceptibility reads χ b = 12 (cid:90) L d x G ( x ) · (cid:40) x , x ≤ L/ x − L ) , x > L/ . (3.6)In Ref. [29], this relation was derived for χ b in momentum space. In App. E we repeat the argumentin coordinate space.In the representation (3.6) the current-current correlator is computed for a general electromagneticfield in coordinate space. Only afterwards a Fourier transformation is carried out via the convolutionwith the quadratic kernel in order to enforce the constancy of the background field. In this way theproblem of flux quantization is avoided. Notice that there is a remnant of flux quantization in theformula (3.6), signaled by the cusp in the kernel at x = L/ . However, this cusp has no practical– 8 –elevance, as in the integral it is multiplied by G ( L/ , which is exponentially small. We checked thatfinite volume effects are indeed small, see Fig. 10 in App. E.At zero temperature, the result (3.6) is very similar to the zero-momentum limit of the formfactor Π appearing in the vacuum polarization tensor Π µν . In fact our result can be expressed as χ b = − ∂ Π ∂p (cid:12)(cid:12)(cid:12)(cid:12) p =0 = Π(0) , Π µν ( p ) = (cid:90) d x e ipx (cid:104) j µ ( x ) j ν (0) (cid:105) = (cid:2) p µ p ν − δ µν p (cid:3) · Π( p ) . (3.7)The vacuum polarization function Π has been the subject of intense research as it appears in thehadronic contribution to the muon anomalous magnetic moment, see, e.g., the recent review [49]. Inthat setting the relevant observable is the second moment of the two-point function of the electro-magnetic current (3.4), projected to zero spatial momentum [50]. Exchanging the time coordinate x for the spatial coordinate x , one can obtain χ b in an analogous way in our background field setup.Since this change is irrelevant at T = 0 , at zero temperature χ b = Π(0) holds [29]. For nonzerotemperatures the Lorentz structure of Π µν ( p ) is more complicated than in Eq. (3.7) and, in additionto Π , a form factor Π L appears [51, 52]. However, in the static case ( p = 0 ) only Π contributes tothe spatial components of Π µν so that the first relation of Eq. (3.7) still holds at T > .In App. E we derive a similar representation for τ fb as well. In this case the equivalents ofEqs. (3.5) and (3.6) become τ fb = iq f /e (cid:90) L d x H f ( x ) · (cid:40) x , x ≤ L/ x − L, x > L/ , H f ( x ) = (cid:90) d x d x d x (cid:10) ¯ ψ f σ ψ f ( x ) j (0) (cid:11) . (3.8)We remark that the second moment of the photon DA is accessible too, replacing σ by combinationsof σ µν (cid:16) ←− D ρ ←− D σ + −→ D ρ −→ D σ − ←− D ρ −→ D σ (cid:17) that are antisymmetric in indices equal to and , symmetrizedover all other non-trivial combinations of indices and with all traces subtracted, see, e.g., Ref. [53].This is beyond the scope of the present work.In summary, via the relation (3.6) we are able to determine the magnetic susceptibility using directmeasurements at B = 0 . This is certainly advantageous over calculating the free energy density (whichcannot be obtained as a simple expectation value) at nonzero magnetic fields and differentiating itnumerically. The similar relation for the tensor coefficient, Eq. (3.8), might also be used to avoidmeasuring (cid:10) ¯ ψ f σ ψ f (cid:11) at B > . However, since the latter is a simple one-point function, the gain isnot obvious in this case. In App. E we compare the two methods for this observable and concludethat the correlator method indeed gives larger statistical errors. Therefore, we opted to use our earlierresults [26] for (cid:10) ¯ ψ f σ ψ f (cid:11) . On general grounds, the linear response of the expectation value of an n -point function with respect to a backgroundfield can always be obtained by computing ( n +1) -point functions in the vacuum. Usually, the former method is favorablein terms of the statistical noise. However, the latter option exempts us from the need of generating additional gaugeensembles with non-vanishing values of the background field. As already discussed, in the present context, we alsocircumvent the issue of flux quantization. – 9 – igure 1 . Bare magnetic susceptibility at zero temperature versus the logarithm of the lattice spacing, normal-ized to a = 1 . GeV − . Different approaches are compared: the finite difference method [47] (red triangles),the generalized integral method [28] (green circles) and the new approach via current-current correlators (bluesquares). The orange band indicates the fit based on perturbation theory, Eq. (4.1). The correlator G ( x ) is calculated using O (1000) random sources located on three-dimensional x -slices of our lattices. Subsequently, G ( x ) is convoluted with the quadratic kernel to calculate χ b viaEq. (3.6). The results of this method, compared to our old data (generalized integral method) andalso to those of Ref. [47] (finite difference method) are shown in Fig. 1 at zero temperature. Within errors perfect agreement between the three groups of results is found. The data – plottedin Fig. 1 against log( a ) – clearly reflect the logarithmic divergence dictated by Eq. (2.5). Similarlyto our fitting strategy in Ref. [29], here we also include the universal perturbative QCD correctionsto the QED β -function coefficient c [34], see Eq. (2.4), where g = 6 /β is obtained from the inverselattice coupling β at the lattice scale a − . We also take into account O ( a ) lattice artifacts so thatour fit function reads χ b = 2 β ( a − ) · (cid:2) log( a/a ) + log( µ QED a ) (cid:3) · (cid:2) z ( a/a ) (cid:3) , a = 1 . GeV − . (4.1)The result of this fit, with the parameter values µ QED = 115(3)(5)
MeV , z = − . , (4.2)is shown as an error band in Fig. 1. We also considered fits with further (quartic) lattice artifacts.The impact of this is included in the second parentheses of Eq. (4.2) for µ QED as a systematic error.The renormalization scale agrees within errors with our earlier determinations [28, 29]. It also liesnear the mass of the lightest charged hadron (the charged pion) which effectively sets the scale for the Ref. [47] employs the same lattice action. For a different action the bare susceptibilities would not only differ interms of lattice artifacts but also by an additive constant, due to the different choice of renormalization scheme. – 10 – igure 2 . Left panel: renormalized susceptibility at nonzero temperature. The symbols indicate differentlattice spacings and the dark orange band the continuum limit. The light orange band represents an estimateof systematic errors of the continuum extrapolation. The dashed gray line is the HRG model prediction [28].The inset zooms into the low-temperature region to highlight the diamagnetic response there. Right panel:our results for χ (orange, labeled Π(0) ) are compared to the results of Ref. [47] (green) and those of Ref. [28](red) as well as to the HRG model prediction (dashed gray line) and to perturbation theory [28] (dashed lightblue band), see Eqs. (4.3) and (2.4). magnetic response of this system. Nevertheless, note that the value of µ QED depends on the choiceof the regulator.The formula (4.1) is used to interpolate χ b and then employed to renormalize the susceptibilityat nonzero temperatures according to Eq. (2.6): χ = χ b ( T ) − χ b (0) . The results are shown in theleft panel of Fig. 2 for a broad range of temperatures and four lattice spacings a = 1 / ( N t T ) with N t = 6 , , and . A continuum extrapolation is performed by means of a multi-spline fit [54] takinginto account O ( a ) lattice artifacts. The coarsest lattices ( N t = 6 points at T (cid:46) MeV) were foundto lie outside the scaling region and were therefore excluded from the analysis. The systematic errorwas estimated by varying the spline node points and by including/excluding N t = 6 data pointsat high temperatures. In addition, we consider a further systematic error due to lattice artifactsrelated to the taste splitting of the staggered spectrum. This effect is particularly relevant at lowtemperatures and can be estimated by a generalization of the Hadron Resonance Gas (HRG) modelthat we describe in App. A. The light yellow bands in both panels of Fig. 2 include this uncertainty.The results demonstrate strong paramagnetism in the quark-gluon plasma phase, in agreementwith previous lattice studies [27, 28, 45–48]. At high temperatures the results are well described bythe free theory, which predicts (see App. C) χ ( T ) = β ( µ therm ) · log (cid:32) γ T µ (cid:33) + O (1 /T ) , (4.3)where γ = O (1) is a regulator-dependent constant. As indicated, QCD corrections at the thermal For example in the free case with cut-off regularization γ = π e − γ E , see App. C.5. The renormalized susceptibility χ ( T ) and the ratio γ/µ are regulator-independent. – 11 – igure 3 . Left panel: lattice discretization errors in the renormalized magnetic susceptibility. The resultsof Ref. [47] (green) and those of Ref. [28] (red) are compared to the present approach (yellow) includingsystematic uncertainties (light yellow). Note that for the green points at a > , χ was obtained by temperature-interpolations of the results published in Ref. [47]. Right panel: parametrization of the relative magneticpermeability µ/µ = (1 − e χ ) − via the function (F.7) of App. F. scale µ therm affect the leading behavior. In the right panel of Fig. 2 we also include a comparisonto this perturbation theory formula, with the scheme-independent O ( α s ) corrections to β takeninto account [28]. Here we set γ = 1 and use the MS scheme definition of the strong coupling α s = g / (4 π ) , running this at five-loop order [55] to the thermal scale µ therm ∼ πT . We employ thecentral value Λ MSQCD = 0 . GeV from the recent three flavor QCD determination of α s by the ALPHACollaboration [56]. The band in the figure corresponds to a variation of the thermal scale from πT to πT . The perturbative formula agrees surprisingly well with our results down to temperatures T ∼ MeV.In contrast to the paramagnetic behavior at high T , towards low temperatures the continuumextrapolated results become negative, revealing a diamagnetic response, previously noted in Ref. [28].This behavior is in agreement, albeit within large errors, with the HRG model prediction [28] (dashedline in the figures). In the right panel of Fig. 2 we also include results obtained from other approachesthat employed the same lattice action. Around the pseudo-critical temperature T c ≈ MeV asignificant difference is visible between earlier determinations and our present results. Ref. [47] carriedout a continuum extrapolation using fixed β ensembles with lattice spacings a ≥ . fm. In Ref. [28]we used the fixed N t approach with N t = 6 , , ensembles, while in the present study N t = 6 , , , lattices are simulated. To highlight the differences between the continuum extrapolations, in the leftpanel of Fig. 3 we plot the lattice spacing-dependence of the susceptibility for all three approaches. Wepick one temperature T = 130 MeV, where the deviation of the continuum estimates is substantial.The left panel of Fig. 3 reveals the importance of our new N t = 12 ensemble, showing a significantdownward trend as the lattice spacing is reduced and a negative value in the continuum limit. Thedownward trend is not captured by our previous extrapolation using the integral method on N t ≤ lattices [28], neither is it visible in the data of Ref. [47]. We remark furthermore that Ref. [47]performed the continuum extrapolation assuming a strictly positive function for χ ( T ) . Excluding the– 12 –ossibility of a negative susceptibility in the continuum limit might also bias the extrapolation. Insummary, we suspect that the qualitative deviation between the continuum limits below the pseudo-critical temperature is due to the absence of N t = 12 ensembles in Ref. [28] and to systematicuncertainties in the a → extrapolation procedure of Ref. [47]. To clarify this issue, dedicatedsimulations should be performed, preferably at the same temperature and the same values of thelattice spacing.Finally we provide a parametrization that connects all three approaches (HRG, lattice continuumlimit and perturbation theory) and describes χ for arbitrary temperatures. The details are discussedin App. F. In the right panel of Fig. 3 we plot this parametrization, translated to the magneticpermeability µ/µ = (1 − e χ ) − , expressed in units of the vacuum permeability µ . This combinationis equal to the ratio of the magnetic induction and the external field, see, e.g., Refs. [28, 46]. Here we address the tensor coefficients τ fb at zero temperature. We consider a set of independentgauge ensembles, generated at the physical value of the strange quark mass m s = m phys s , but atdifferent values of the light quark mass: . m phys ud ≤ m ud ≤ m phys s . We follow a similar strategy asin Ref. [26], simultaneously fitting the dependence of τ ub · Z T on the light quark mass m ud and onthe lattice spacing a according to the ansatz (2.9). Since Z T is found to depend very mildly on thelattice spacing within the range covered (see Fig. 9 of App. D), this does not significantly affect thefunctional dependence on a . We note that on our coarsest ensembles the uncertainty of Z T is quitelarge, which also imprints on the errors of the renormalized tensor coefficients. Figure 4 . Light quark mass-dependence of the tensor coefficient at T = 0 for the up quark using our four finestlattice spacings (green to gray symbols). The index b indicates that the QED divergence that one encountersat m ud > has not been subtracted. The results diverge logarithmically towards the continuum limit for any m ud (cid:54) = 0 . In contrast, the chiral limit is free of ultraviolet divergences and a combined chiral and continuumlimit exists (black circle). Notice that τ fb diverges for a → for any quark mass, except in the chiral limit, where itis ultraviolet-finite. This tendency is clearly visible in Fig. 4, which shows our results for the up– 13 – igure 5 . The bare tensor coefficient for the up quark (left panel) and for the strange quark (right panel) at thephysical point and at zero temperature (blue points), together with an interpolation (orange bands). For theup quark this interpolation is the m ud = m phys ud slice of a two-dimensional fit like in Fig. 4. The red dashed linesindicate the leading logarithmic divergence ∝ m f log a in both fits. The remaining a -dependence is consistentwith lattice artifacts. We always use m s = m phys s . The lattice spacing is normalized to a = 1 . GeV − . quark. Therefore, we can define an ultraviolet-finite observable for the light quarks, without any zero-temperature subtraction, namely the chiral limit of the tensor coefficient. In contrast, to calculate χ spin we will need to take differences between results obtained at different temperatures (see below).The ansatz (2.9) contains the free parameters τ f and µ QED . In addition, we include a quadraticmass-dependence and lattice artifacts of O ( a ) to each parameter in the fit. Varying the fit ranges in a and in m ud /m phys ud as well as the functional form, we carry out several acceptable fits that are usedto build a histogram for the chiral continuum limit of the tensor coefficient. In this combined limitwe obtain in the MS scheme T = 0 : f ⊥ γ (2 GeV ) ≡ lim m ud → τ ub · Z T (2 GeV ) = − . . MeV . (4.4)The central value differs from our previous result [26] f ⊥ γ = − . . MeV, mainly due to themultiplicative renormalization factor that we determined non-perturbatively here, see Fig. 9 in App. D.In the left panel of Fig. 5 we show the a -dependence of the T = 0 light quark tensor coefficient τ ub · Z T at the physical point and the result of the above interpolation, including the systematicerror estimated using the different fits. For demonstration purposes, we also indicate the leadinglogarithmic behavior, that we obtain by subtracting the lattice artifact terms from the central fit.Comparing to the similar plot for χ b (Fig. 1), we see that deviations from the continuum behaviorare sizable (and are predominantly due to the fact that we are dealing with a dimensionful quantityin this case). For this reason, here we cannot reliably determine the value of µ QED . Nevertheless,we note that fixing the renormalization scale to its value from Eq. (4.2) also gives acceptable fits.This is in agreement with the expectation of Sec. 2.3, as well as with the results in the free case, seeApp. C.4. The logarithmic divergence ∝ m f log a becomes more pronounced for heavy quarks. Thisis visible in the right panel of Fig. 5, where we plot the strange quark tensor coefficient τ sb · Z T at m ud = m phys ud against the lattice spacing and again indicate the leading logarithmic term.– 14 –s we have discussed above, at non-vanishing values of the quark mass m f , the tensor coefficientdiverges logarithmically. In Ref. [26] we suggested to cancel this by taking the logarithmic derivativewith respect to the quark mass, see also Ref. [14]: f ⊥ γf = (cid:18) − m f ∂∂m f (cid:19) τ f · Z T . (4.5)This renormalization prescription will give identical results for any regulator, up to the multiplicativefactor Z T . Using this prescription, it turns out that f ⊥ γu = f ⊥ γd = f ⊥ γ holds within statistical errors.For the strange quark we obtain: f ⊥ γs = − MeV . (4.6)The first error includes the described variation of the fit while the second error reflects the uncertaintyof the derivative with respect to m s that we indirectly determine from the dependence of the tensorcoefficient on the light quark mass, following the procedure explained in Ref. [26].In the literature often the magnetic susceptibility of the quark condensate, X u = τ u (cid:10) ¯ ψ u ψ u (cid:11) · Z T Z S , (4.7)is given, rather than f ⊥ γ = τ u · Z T . Since the latter quantity has a smaller anomalous dimension andits value does not depend on a separate computation of the chiral condensate, this is the preferredchoice for practical applications. However, for convenience of comparison, we shall convert it into theother convention. The numerical value of the quark condensate in the SU(2) chiral limit in the MS scheme at the scale µ QCD = 2
GeV reads (cid:104) ¯ ψ u ψ u (cid:105) = [272(5) MeV ] [58]. To enable a comparison withother results, below we also list X u at the scale µ QCD = 2
GeV. Since most literature values refer toa low, sometimes unspecified scale, in addition we run X u as well as f ⊥ γ to the scale µ QCD = 1
GeV,which is used in most sum rule calculations, see, e.g., Refs. [16, 17, 19]. This is done, using thefive-loop β - and quark mass anomalous dimension γ -functions [55, 59] and the three-loop γ -functionof the tensor current [36, 60]. The results read X u (2 GeV ) = − [665(13) MeV ] − , X u (1 GeV ) = − [542(11) MeV ] − , (4.8) f ⊥ γ (1 GeV ) = − . . MeV , (4.9)where we have added all errors in quadrature, including the uncertainty of f ⊥ γ (2 GeV ) , the differencebetween running with the two- and three-loop γ -functions of the tensor current, the uncertainty of (cid:104) ¯ ψ u ψ u (cid:105) and the uncertainty of the strong coupling parameter [56]. All the above results are in the MS -scheme.We summarize earlier results from the literature for comparison. The first sum rule determina-tion of X u [14] suggested a value X u (0 . GeV ) = − [350(50) MeV ] − while vector meson dominanceyields [35] X u ≈ /m ρ ≈ − (540 MeV ) − . This was improved upon in subsequent sum rule deter-minations, see, e.g., Ref. [19] and references therein. The most extensive sum rule study [19] found X u (1 GeV ) ≈ − (560 MeV ) − , which agrees reasonably well with our determination. A comparatively This construction will not only cancel the logarithmic divergence but also any finite term ∝ m f . Should this beunwanted then one will have to accept a scheme-dependence and convert between different schemes in a similar way asis done for the massive chiral condensate, e.g., in Ref. [57]. – 15 –maller absolute value f ⊥ γ ≈ − MeV was obtained at a low scale µ ∼ MeV in the quark-solitonmodel [18] while the Vainshtein relation [61] suggests an even smaller modulus of the magnetic suscep-tibility of the quark condensate X u = − N c / (4 π F π ) ≈ − (335 MeV ) − . This parameter was also con-sidered in holographic studies, with the result X u ≈ − (295 MeV ) − [62], while NJL- and quark-meson-model predictions give X u ≈ − (480 MeV ) − [63] and X u ≈ − (440 MeV ) − [63], respectively. Finally,quenched lattice simulations, without renormalization, gave the values X u ≈ − [804(3) MeV ] − inSU(2) gauge theory [64] and X u ≈ − [486(21) MeV ] − in SU(3) [65]. Our previous full QCD study [26]resulted in X u (2 GeV ) = − [693(13) MeV ] − , however, in that case the renormalization was only car-ried out perturbatively.Our result (4.6) for the strange quark coefficient translates into X s (2 GeV ) = − [565(50) MeV ] − , X s (1 GeV ) = − [460(41) MeV ] − , (4.10) f ⊥ γs (1 GeV ) = − . . MeV , (4.11)where we used the ratio (cid:104) ¯ ψ s ψ s (cid:105) / (cid:104) ¯ ψ u ψ u (cid:105) = 1 . [57] for the conversion between f ⊥ γs and X s . Thedifference between X s and X u (1 GeV ) ≈ − (475 MeV ) − has been reported to be negligible in thesum rule calculations [17]. T > Having interpolated the T = 0 tensor coefficients, we are now in the position to perform the additiverenormalization (2.10) by subtracting this contribution from the finite temperature results. We useour existing N t = 6 , and results from Ref. [26] to approach the continuum limit. Especially inlight of the slow convergence of χ towards a → , see the right panel of Fig. 3, this extrapolationshould be backed up with finer lattice ensembles in the future. In analogy to the analysis of χ , againwe carry out a multi-spline fit of all data sets, determining a systematic error by varying the positionsof the spline node points. The so-obtained fit is shown for the up quark and the strange quark inFig. 6. The results for τ d are consistent with τ u within errors. The large errors of our N t = 6 resultsat low temperatures are due to the uncertainties of the T = 0 contributions on our coarse lattices,see Fig. 5.After the additive renormalization, the tensor coefficient vanishes by definition at T = 0 . For thelight quarks τ u ( T ) grows substantially as the temperature is increased, before the slope reduces and aplateau is approached. The inflection point of the continuum curve is found to be at T c = 158(5) MeV.The chiral transition temperature determined from the inflection point of the quark condensate T c =155(4) MeV [66] is in agreement with this value. For the strange quark pseudo-critical thermal effectsset in at somewhat higher temperatures [66]. Also in our case τ s does not appear to exhibit anyinflection point, at least for T (cid:46) MeV, and below T ≈ MeV no saturation into a plateau isvisible. For sufficiently high temperatures, where the finite quark mass becomes negligible, we expectthe two renormalized tensor coefficients to coincide.Next, the continuum extrapolated results are inserted into Eq. (2.12) to determine the spincontribution χ spin to the susceptibility. To this end we need to evaluate the tensor bilinear formassless valence quarks. Instead of performing measurements at additional valence quark masses, weestimate this limit using the difference between the results for the strange quark and for the lightquarks. We assume a linear dependence on the valence quark mass in the range [0 , m s ] , which implies– 16 – igure 6 . Tensor coefficients after multiplicative as well as additive renormalization for the up (left panel)and for the strange quark (right panel). that lim m val u → τ u = lim m val s → τ s ≈ τ u m s − τ s m ud m s − m ud = τ u RR − − τ s R − , R ≡ m s m ud . (4.12)In this case the contributions of all flavors to χ spin are proportional to τ s − τ u and the renormalizedspin susceptibility (2.12) simplifies to χ spin ≈ m ud τ s − τ u R − · Z T Z S · (cid:88) f ( q f /e ) . (4.13)Thus, in this approximation the individual flavors simply contribute in proportion to their squaredelectric charges. The scalar renormalization constants entering this expression are displayed in Fig. 9of App. D.The so-obtained estimate of χ spin is shown in the left panel of Fig. 7 for three lattice spacings,together with a continuum extrapolation performed in the same way as for τ f . We observe χ spin < for all temperatures, with a minimum somewhat above the pseudo-critical temperature and an upwardtrend for high temperatures. The approximation (4.12) tends to overestimate the valence chiral limitof the tensor coefficient due to the presence of logarithmic deviations from a linear behavior in m val f . Consequently, Eq. (4.13) underestimates χ spin . This is also the case at high temperatures, as can bechecked using the analytic formula valid for the free case, see App. C. To take this effect into accountwe include a systematic error based on the free case formula (C.25). In particular, we consider thedifference between the approximation and the true value in the free case and scale it with the typicalmagnitude of the light quark tensor coefficient at lower temperatures (see Fig. 6). The so-obtaineduncertainty is also included in the left panel of Fig. 7.We remark that χ spin < for the temperature range covered in our simulations. This can beunderstood by noting that Eq. (4.13) is the discretization of the mass-derivative of τ f . Increasing This is also visible in Fig. 4, although the dependencies on the valence and sea quark masses are not disentangledin that figure. – 17 – igure 7 . Left panel: spin contribution to the susceptibility using three lattice spacings (colored symbols) andan extrapolation to the continuum limit (orange band). A systematic uncertainty, related to the estimationof the tensor coefficient for massless valence quarks, is indicated by the light yellow band. Right panel: thetotal magnetic susceptibility from Fig. 2 (blue), together with the decomposition into spin (orange-yellow) andorbital angular momentum (green-gray) contributions. the mass pushes the inflection point of τ f to higher temperatures (visible in Fig. 6), thus makingthe derivative negative around the transition temperature. Nevertheless, χ spin will necessarily turnpositive for even higher temperatures. Indeed, for sufficiently high temperatures the difference τ f = τ fb ( T ) − τ fb ( T = 0) will be dominated by the T = 0 term, so that Eq. (4.13) becomes proportionalto τ ub ( T = 0) − τ sb ( T = 0) , which is positive for any lattice spacing (see Fig. 4). Perturbation theoryalso predicts χ spin > for high temperatures, see App. C.3. Finally, we compare the spin contribution to the total susceptibility in order to learn about the orbitalangular momentum-related contribution χ ang = χ − χ spin . All three susceptibilities are included inthe right panel of Fig. 7. While the errors of the two contributions are much larger than that ofthe total susceptibility, several qualitative comments can be made based on this plot. First of all,in the complete temperature range under study, χ spin and χ ang have opposite signs and χ emergesas a result of a large cancellation between the two terms. As we argued above, the spin part willnecessarily turn positive for higher temperatures, eventually approaching / times the full suscepti-bility. Consequently, χ ang will turn negative and approach − / · χ . It is intriguing to observe thatin the strongly interacting regime the two contributions have opposite signs than in the usual freefermion picture according to Pauli and Landau: it is the Landau term that drives the paramagneticresponse of the QCD vacuum up to temperatures T (cid:38) MeV, while the Pauli term reduces thesusceptibility in this region. This unusual behavior becomes possible due to the strong interaction,which confines quarks into composite hadrons and thereby fixes the relative orientation of their spins,i.e. their magnetic moments. In particular, in charged pions one of the constituent quarks is boundto anti-align its magnetic moment with the background field in order to maintain zero total spin.Similar effects arise for certain baryons as well. Beyond this qualitative argument, it is difficult to– 18 –nticipate the outcome of this competition between the strong and the electromagnetic forces. Ourquantitative results reveal a peculiar interplay between confinement and spin physics.To further our understanding, in principle χ can also be decomposed into χ f for the quark flavors f and a gluonic contribution χ g . Subtracting this χ g from χ ang will isolate the total quark orbitalangular momentum contribution (cid:80) f ( χ f − χ spin f ) , in analogy to spin decompositions [11] in deepinelastic scattering that are based on the Belinfante-Rosenfeld energy-momentum tensor, in this caseof the transverse spin. The unrenormalized qualitative results of Ref. [27] indicate that χ g ∼ χ/ atsmall temperatures. It may be interesting to address this quantitatively in the future. In this paper we determined the magnetic susceptibility χ of the thermal QCD medium via a methodintroduced originally for T = 0 [29], which circumvents the flux quantization problem and allows usto express χ in terms of B = 0 measurements. This considerably reduces the measurement costs aswell as systematic uncertainties compared to previous approaches. The susceptibility is extrapolatedto the continuum limit for a broad range of temperatures, making contact to the Hadron ResonanceGas (HRG) model at low T as well as to perturbation theory at high T . In the confined phase wefind evidence for a diamagnetic behavior ( χ < ), while for T (cid:38) MeV we observe paramagnetism( χ > ). Our continuum extrapolations are based on four lattice spacings and are guided by ageneralized HRG model taking into account taste splitting (see App. A). A careful continuum limitis found to be essential to observe diamagnetism at low T since this is due to light pions – we arguethat this behavior was missed in previous investigations because of large lattice artifacts.The susceptibility is decomposed into spin- ( χ spin ) and orbital angular momentum-related ( χ ang )contributions based on our previous study [26]. The spin term is shown to be given in terms ofthe mass-dependence of the (cid:10) ¯ ψσ ψ (cid:11) fermion bilinear in the presence of a small magnetic field, seeEq. (2.12) and App. B. Besides its role in this decomposition, the tensor bilinear is related to thenormalization f ⊥ γ of the photon distribution amplitude, relevant for a range of phenomenologicalapplications. We update our previous determination [26] of the corresponding tensor coefficient in thechiral limit at T = 0 , by performing the multiplicative renormalization of (cid:10) ¯ ψσ ψ (cid:11) non-perturbativelyon the lattice. We obtain the value f ⊥ γ = − . . MeV for massless quarks, in the MS scheme ata QCD renormalization scale of GeV. The values of the tensor coefficient at the physical light andstrange quark masses and at different renormalization scales are given in Eqs. (4.6)–(4.11).At finite temperatures we performed the continuum extrapolation of χ spin and also determinedthe orbital angular momentum-related susceptibility χ ang . In the absence of color interactions, thetwo contributions exhibit the constant ratio χ spin : χ ang = 3 : ( − as is well known since the analysisof the free electron gas by Pauli [38] and Landau [39]. Around the transition temperature, in full QCDthis ratio is instead found to be close to ( −
1) : (1 . , resulting in a large cancellation between thetwo contributions, thereby substantially reducing the total susceptibility. As the temperature growsthe susceptibilities approach their free-case counterparts, which are discussed in detail in App. C.Still, it is stunning to observe that in the strongly coupled QCD medium χ spin and χ ang have signsthat are opposite to the naive expectations.Considering our results at high temperature, it is interesting to make a comparison to a classicalideal system. In such a setting the Bohr-van Leeuwen theorem [67, 68] (for a recent review, see– 19 –ef. [69]) holds: the total magnetization vanishes, since the magnetic field does not transfer any workto the electric currents in the system. Apparently, the QCD medium does not become classical in thissense for T → ∞ , even if the O ( B ) terms of the free energy density that we have discussed in thispaper are small compared to the dominant O ( T ) contributions in that limit. The non-classicalityhas two different origins. First, quark spins are of quantum nature and can induce a magnetization byaligning with the magnetic field. Second, both χ spin and χ ang diverge as log T for high temperatures.This behavior stems from the renormalization properties of the bare susceptibilities: quantum effectsgive rise to a logarithmic divergence ∝ log 1 /a in the cut-off. In turn, the same behavior showsup in the renormalized susceptibilities if they are probed by another large dimensionful scale, thetemperature. Note that a similar connection exists between the logarithmic divergence and thebehavior of the renormalized free energy in the B → ∞ limit [33]. Acknowledgments
This research was funded by the DFG (Emmy Noether Program EN 1064/2-1 and SFB/TRR 55).The authors would like to thank Falk Bruckmann and V. M. Braun for enlightening discussions aswell as Nikolay Kivel and Massimo D’Elia for insightful comments.
A The HRG model and lattice discretization errors
At low temperatures the staggered action suffers from enhanced lattice artifacts due to taste splitting.Here we attempt to incorporate the effects of this splitting into the HRG model. The magneticsusceptibility was calculated in a standard HRG model in Ref. [28]. Following Ref. [70] we replacethe contribution of pions in the model by a sum over each taste, weighted by the correspondingdegeneracies. The masses of the individual tastes and their parametrization in the range of our latticespacings are taken from Ref. [41]. Since pions are dominant for the susceptibility, the taste splittingfor other mesonic and baryonic states is ignored for simplicity (although the splitting for η mesonsmight also lead to light mesonic states, see, e.g., Ref. [71]). The list of hadrons taken into accountcan be found in Ref. [72].In Fig. 8 we show the renormalized magnetic susceptibility evaluated at T = 120 MeV as afunction of the lattice spacing. The spacings for our four ensembles N t = 6 , , and at thistemperature are highlighted in the plot. This reveals slow convergence towards the continuum limit,which can best be understood by analyzing the mass-dependence of the pionic contribution χ π to thesusceptibility, which takes the form [28] χ π ( m π ) = − π (cid:90) ∞ d tt e − m π t/T (cid:104) Θ (cid:16) , e − / (4 t ) (cid:17) − (cid:105) , (A.1)where Θ is an elliptic Θ -function. This can be derived by comparing to the analogous expressionfor fermions, calculated below in Eq. (C.8). The bosonic Matsubara frequencies give rise to thedifferent first argument in the elliptic function. The prefactor in this case is the scalar QED β -function coefficient for one complex scalar field β scalar1 = 1 / (48 π ) . The pionic susceptibility divergeslogarithmically in the chiral limit (this can be shown similarly to the calculation below in App. C.5), χ π ( m π ) m π → −−−−→ − β scalar1 log( T /m π ) , (A.2)– 20 – igure 8 . Lattice artifacts in the susceptibility in a generalized HRG model involving taste splitting. explaining its pronounced dependence on m π . In turn, nonzero lattice spacings enhance the massesof most pion tastes, thus, reducing the magnitude of χ π .Based on the HRG predictions for χ ( a, T ) we consider the difference between a simple O ( a ) fittaking into account only N t ≤ lattices and the true continuum limit. This difference is includedas a lower systematic error of our lattice determination of χ ( T ) at low temperatures, see Fig. 2. B Separation into quark spin and other angular momentum contributions
Here we derive the relation between the spin contribution to the susceptibility and the tensor bilinear,as shown in Eqs. (2.8)–(2.12) of the main text. It is instructive to begin with the first derivative ofthe free energy density, − ∂f∂B = TV (cid:88) f (cid:28) tr /D f + m f ∂ /D f ∂B (cid:29) = T V (cid:88) f (cid:42) tr /D f + m f ) /D f ∂ /D f ∂B (cid:43) , (B.1)where we used the cyclicity of the trace (even though /D f and ∂ /D f /∂B do not commute, we cansymmetrize the expression in the two operators under the trace). Now we use the relation /D f + m f ) /D f = − m f (cid:20) /D f + m f − /D f (cid:21) , (B.2)and the identities ∂ /D f ∂ ( q f B ) = − σ − L , σ = 12 i [ γ , γ ] , L = − ∂D f ∂ ( q f B ) , (B.3)where σ is the relevant component of the relativistic spin operator defined in Eq. (2.7) and L is ageneralized angular momentum operator, which depends on the electromagnetic as well as the SU(3) gauge. – 21 –sing Eqs. (B.2) and (B.3), we can rewrite Eq. (B.1) as − ∂f∂B = T V (cid:88) f q f m f (cid:28) tr σ + L /D f + m f − tr σ + L /D f (cid:29) = (cid:88) f q f m f (cid:20) − lim m val f → (cid:21) (cid:10) ¯ ψ f σ ψ f + ¯ ψ f L ψ f (cid:11) . (B.4)Thus, in the language of Eq. (3.1), we need the difference of two terms: one with valence quark mass m val f = m f and one with m val f → . The sea quark mass is kept fixed in both cases: m sea f = m f . Weremark that the vanishing valence quark mass needs to be defined as a limit in finite volumes (seebelow). Also note that the fermion bilinears are defined to include the volume factor T /V .Differentiating Eq. (B.4) once more with respect to B at B = 0 and dividing by e , we recoverthe bare magnetic susceptibility (2.1) on the left hand side, χ b = (cid:88) f ( q f /e ) m f (cid:20) − lim m val f → (cid:21) lim B → (cid:10) ¯ ψ f σ ψ f + ¯ ψ f L ψ f (cid:11) q f B . (B.5)The slope of the tensor bilinear (cid:10) ¯ ψ f σ ψ f (cid:11) for small values of B gives the tensor coefficient τ fb asdefined in Eq. (2.8). After subtracting its value at T = 0 and multiplying by the relevant QCDrenormalization factors, this term gives the spin contribution to the susceptibility χ spin , as we wrotein the main text, Eq. (2.12). In turn, the magnetic field-dependence of the bilinear involving the gen-eralized angular momentum operator L is related to χ ang . The latter term cannot be implementedstraightforwardly due to its gauge-dependence and magnetic flux quantization.In Ref. [26] we already discussed the separation of the magnetic susceptibility into quark spin- andother angular momentum-related contributions. There, the m val f = 0 term of Eq. (B.4) was arguednot to contribute – indeed, in a finite volume the massless limit of fermion bilinears always vanishes.However, in the thermodynamic limit this is not the case if chiral symmetry is broken spontaneously.To elucidate this point in more detail, let us rewrite the trace in Eq. (B.4) using the eigenmodes ofthe Dirac operator, /D f χ fλ = iλχ fλ , (B.6)so that, exploiting chiral symmetry { γ , /D f } = 0 , TV (cid:42) tr σ /D f + m val f (cid:43) = T m val f V (cid:42) tr σ − /D f + ( m val f ) (cid:43) V →∞ −−−−→ (cid:90) ∞ d λ m val f λ + ( m val f ) (cid:68) ρ f ( λ ; m sea f ) χ † fλ σ χ fλ (cid:69) , (B.7)where ρ f ( λ ; m sea f ) is the spectral density of /D f in the infinite volume, determined in an ensemble gen-erated with sea quark masses m sea f . Towards the valence chiral limit the kernel becomes proportionalto the δ -distribution, so that we have a Banks-Casher-type [73] relation, (cid:10) ¯ ψ f σ ψ f (cid:11) V →∞ , m val f → −−−−−−−−−→ π (cid:90) ∞ d λ δ ( λ ) (cid:68) ρ f ( λ ; m sea f ) χ † fλ σ χ fλ (cid:69) = π (cid:68) ρ f (0; m sea f ) χ † f σ χ f (cid:69) . (B.8)On the one hand, this limit is zero if chiral symmetry is intact and the spectral density vanishes at theorigin. On the other hand, a nonzero chiral condensate (cid:68) ρ f (0; m sea f ) (cid:69) , together with the polarization– 22 – χ f = χ f of the low modes will turn the chiral limit of the tensor bilinear nonzero. Our latticeresults reveal a nonzero value for (cid:10) ¯ ψ f σ ψ f (cid:11) in the full chiral limit at low temperatures, see Fig. 4.Clearly, the fermion bilinear remains nonzero also if only m val f is sent to zero. This is in accordancewith the recent findings of Ref. [74] about the Dirac spectrum at B > , where the low modes wereindeed found to exhibit almost perfect spin-polarization. C Susceptibilities in the free case
Here we consider the free case (i.e. we set the color charges of quarks to zero) to exemplify the mostimportant relations of the main text. These include the proportionality between the tensor bilinearand the spin contribution to the susceptibility, the ultraviolet divergences of the susceptibilities atzero temperature as well as the high-temperature behavior of the renormalized susceptibilities. Thesecalculations include our previous results [26, 28], which we also show here for completeness.Below we will extensively use Schwinger’s proper time formulation [32]. This is based on theMellin transform E − z = 1Γ( z/ (cid:90) ∞ d t t z/ − e − E t , (C.1)and its inverse e − lE/T = 12 πi (cid:90) c + i ∞ c − i ∞ d z Γ( z ) l − z E − z T z , (C.2)which are valid for Re z > , c > and E > . Moreover, taking the derivative of Eq. (C.1) withrespect to z at z = 0 gives the standard ζ -function regularization result [75], log E = − ∂ ( E ) − z/ ∂z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z =0 = − ∂∂z (cid:12)(cid:12)(cid:12)(cid:12) z =0 z/ (cid:90) ∞ d t t z/ − e − E t . (C.3) C.1 Magnetic susceptibility
We consider one quark flavor ψ with electric charge q and mass m in a volume V = L at temperature T , exposed to a background magnetic field B . For convenience we assume that the magnetic fieldis oriented in the x direction and qB > . The free energy density in this setting reads (see, e.g.,Ref. [76]): f ( B, T ) = − N c qB π ∞ (cid:88) k =0 (cid:88) s = ± / T ∞ (cid:88) n = −∞ (cid:90) d p π log ω n + E p,s,k T , (C.4)where ω n = (2 n + 1) πT is the n -th fermionic Matsubara frequency. Moreover, p , s and k are themomentum, spin and angular momentum in the direction of the magnetic field, N c = 3 is the numberof colors and the energies are given by the Landau levels, E p,s,k = (cid:112) p + m + (2 k + 1 − s ) qB . (C.5)Rewriting the logarithm using Eq. (C.3), the integral over p becomes Gaussian and can be solved.Furthermore, the sums over n , k and s are T ∞ (cid:88) n = −∞ e − ω n t = 12 √ πt Θ (cid:16) π , e − / (4 tT ) (cid:17) , ∞ (cid:88) k =0 e − (2 k +1) qB t = 12 sinh( qBt ) , (cid:88) s = ± / e − sqB t = 2 cosh( qBt ) , (C.6)– 23 –here Θ is an elliptic function. Inserting these in Eq. (C.4) and performing the derivative withrespect to z , we obtain f ( B, T ) = N c qB π (cid:90) ∞ d tt e − m t coth( qBt ) Θ (cid:16) π , e − / (4 tT ) (cid:17) . (C.7)Taking the second derivative with respect to eB to obtain the bare magnetic susceptibility (2.1) resultsin χ b ( T ) = − N c π ( q/e ) (cid:90) ∞ d tt e − m t Θ (cid:16) π , e − / (4 tT ) (cid:17) . (C.8) C.2 Ultraviolet divergences and QED renormalization
To determine the ultraviolet structure of the magnetic susceptibility, we consider Eq. (C.8) at zerotemperature. For T = 0 the elliptic function Θ approaches unity. The resulting expression needs tobe regularized, for example by setting an ultraviolet cut-off / Λ as the lower limit of the proper timeintegral. Performing the integral and expanding for large Λ we obtain, χ b ( T = 0) = N c π ( q/e ) (cid:20) log Λ m − γ E (cid:21) + O (Λ − ) , (C.9)where γ E is the Euler-Mascheroni constant. Thus, the coefficient of the logarithmic divergence indeedequals the lowest-order QED β -function coefficient β (for one quark flavor with electric charge q ),demonstrating the validity of Eq. (2.5). In fact, this relation continues to hold in full QCD as well,owing to the fact that towards the continuum limit QCD corrections to β at the scale /a approachzero due to asymptotic freedom (see Eq. (4.1)). We note moreover that in the proper time formulationthe renormalization scale is set by the mass – in fact µ QED = m e γ E / for our choice of the regulator Λ – explaining the appearance of m in the argument of the logarithm in Eq. (C.9).The additive renormalization can be performed by subtracting χ b ( T = 0) from Eq. (C.8): χ = χ b ( T ) − χ b ( T = 0) = − N c π ( q/e ) (cid:90) ∞ d tt e − m t (cid:104) Θ (cid:16) π , e − / (4 tT ) (cid:17) − (cid:105) . (C.10)As we mentioned after Eq. (2.6), this corresponds to the choice of a physical, albeit scheme-dependent,QED renormalization scale. C.3 Spin contribution
The contribution of orbital angular momentum to the total susceptibility can be calculated by simplyreplacing the fermion with two ghost particles (spin-zero but antiperiodic in Euclidean time) in theabove calculation. This removes the − sqB from the energies (C.5) and excludes the spin sum in thefree energy density (C.4). Consequently, the magnetic field-dependent part in Eq. (C.7) changes as coth( qBt ) (cid:55)→ / sinh( qBt ) . This merely changes the second derivative of the free energy density at B = 0 by a factor − / . Thus, for the renormalized susceptibility we arrive at χ ang ( T ) = − · χ ( T ) , (C.11)which also implies χ spin ( T ) = 32 · χ ( T ) , (C.12)– 24 –onfirming the − ratio of the two contributions to the total susceptibility. We mention thata similar argument has been used in perturbative QCD (with chromomagnetic background fields) torelate asymptotic freedom to spin effects [77]. C.4 Tensor bilinear
For the tensor bilinear we begin with the result of the fermionic path integral, (cid:10) ¯ ψσ ψ (cid:11) = TV tr σ /D + m = T mV tr σ − /D + m , (C.13)where we used chiral symmetry { γ , /D } = 0 . The trace is represented using the eigenbasis of − /D ,giving the eigenvalues ω n + E p,s,k . Since [ /D , σ ] = 0 , the spin operator can also be diagonalizedin this basis and its eigenvalues are minus two times the spin: σ → − s . Taking into account the N c · ( qBL ) / (2 π ) -fold degeneracy of the eigenvalues, we obtain (cid:10) ¯ ψσ ψ (cid:11) = N c qB mπ ∞ (cid:88) k =0 (cid:88) s = ± / T ∞ (cid:88) n = −∞ (cid:90) d p π − sω n + p + m + (2 k + 1 − s ) qB . (C.14)In the sum the contributions { k, s = 1 / } and { k + 1 , s = − / } cancel, leaving only the unpairedlowest Landau level { k = 0 , s = 1 / } . Hence we get (cid:10) ¯ ψσ ψ (cid:11) = − N c qB mπ T ∞ (cid:88) n = −∞ (cid:90) d p π ω n + p + m . (C.15)Note that, unlike in full QCD, here the tensor bilinear is exactly linear in the magnetic field. Thus,the tensor coefficient τ b of Eq. (2.8) is obtained by simply dividing Eq. (C.15) by qB .Using Eq. (C.1) with E = (cid:112) ω n + p + m , performing the Gaussian integral over p and theMatsubara sum (C.6) over ω n , we arrive at τ b ( T ) = − N c m π (cid:90) ∞ d tt e − m t Θ (cid:16) π , e − / (4 tT ) (cid:17) . (C.16)A comparison to Eq. (C.8) reveals that this quantity contains the same logarithmic divergence as χ b ,just with a different coefficient. Using a cut-off regulator as in Eq. (C.9), we obtain at T = 0 , τ b ( T = 0) = N c π m (cid:20) log Λ m − γ E (cid:21) + O (Λ − ) , (C.17)confirming Eq. (2.9). The same considerations regarding QCD corrections to the coefficient and therenormalization scale µ QED apply as in Sec. C.2 for χ b .The difference τ = τ b ( T ) − τ b ( T = 0) is ultraviolet-finite. We can compare this with Eqs. (C.10)and (C.12) to conclude that ( q/e ) m [ τ ( m ) − τ ( m → χ spin , (C.18)confirming the relation (2.12) and Eq. (C.12). Notice that τ vanishes for m → , so the subtractionof the massless limit is irrelevant in the free case (but it is relevant for the interacting system withspontaneous chiral symmetry breaking, see App. B).– 25 – .5 High-temperature expansion The temperature-dependent part of the free energy density (C.4) can be simplified using the well-known trick [52] of differentiating and subsequently integrating the integrand with respect to E p,s,k .The result is f ( B, T ) − f ( B,
0) = − N c qB π ∞ (cid:88) k =0 (cid:88) s = ± / (cid:90) d p π T log (cid:104) e − E p,s,k /T (cid:105) . (C.19)The energy levels are given in Eq. (C.5) above. To obtain the high-temperature expansion in a closedform, we need to replace the logarithm by its series expansion log(1 + x ) = − ∞ (cid:88) l =1 ( − x ) l l . (C.20)This approach was used, e.g., in Ref. [78] for scalars at nonzero chemical potential.Inserting the expansion (C.20) into (C.19) and rewriting the exponentials using Eq. (C.1) resultsin f ( B, T ) − f ( B,
0) = N c qB T π ∞ (cid:88) l =1 ( − l l ∞ (cid:88) k =0 (cid:88) s = ± / (cid:90) d p πi (cid:90) c + i ∞ c − i ∞ d z Γ( z ) l − z E − zp,s,k T z . (C.21)Inserting the Mellin transform (C.1) for E − zp,s,k renders the integral over p Gaussian. We can reuse theangular momentum and spin-sums from Eq. (C.6), giving f ( B, T ) − f ( B,
0) = N c qB π / πi (cid:90) c + i ∞ c − i ∞ d z Γ( z )Γ( z/ T z +1 ∞ (cid:88) l =1 ( − l l z (cid:90) ∞ d t t ( z − / e − m t coth( qBt ) . (C.22)Differentiating the above expression twice with respect to eB at B = 0 gives (minus) the renor-malized magnetic susceptibility χ . The integral over t can be solved via the Mellin transform (C.1)and gives a Γ -function, while the sum over l results in a ζ -function: ∞ (cid:88) l =1 ( − l l z = ζ (1 + z ) · (2 − z − . (C.23)Using the duplication formula [79] for the ratio of Γ -functions, we arrive at χ = − N c π m ( q/e ) πi (cid:90) c + i ∞ c − i ∞ d z Γ (cid:18) z + 12 (cid:19) Γ (cid:18) z + 12 (cid:19) ζ (1 + z ) (1 − z ) m − z T z +1 . (C.24)For the validity of the Mellin transforms we needed to assume c > (as well as m > ). Thefinal integral over z can be solved using Cauchy’s theorem, closing the integral towards the leftand calculating the residue at the poles. There is a double pole at z = − and simple poles at z = − , − , . . . . These z -values set the powers of T that appear in the high-temperature expansion.Keeping the leading terms (i.e. z = − and z = − ), we finally obtain, χ = N c π ( q/e ) (cid:20) log T π m − γ E (cid:21) + N c π ( q/e ) m T (cid:20) C log T m + D (cid:21) + O ( m /T ) , (C.25)– 26 –ith the numerical constants C ≈ . and D ≈ . . Notice that the coefficient of the leadinglogarithmic term is equal to β (for one flavor with electric charge q ), confirming Eq. (4.3). As wehave seen below Eq. (C.9), in the proper time formulation the renormalization scale is intrinsicallyset by the quark mass, µ QED = m e γ E / . We may express the square bracket in the leading termas log( γ T /µ ) with γ = π e − γ E . Clearly, γ = O (1) depends on the definition of the regulator.The general form is again expected to hold in full QCD [28]: in this case QCD corrections at scales T (cid:29) µ QED are small due to asymptotic freedom.
D Multiplicative QCD renormalization
Since lattice perturbation theory is slowly convergent and high-loop results are unavailable, we firstmatch the local lattice QCD operators of interest non-perturbatively to the regulator independentRI’-MOM scheme [30, 31] and subsequently translate the result at three-loop order [36] to the MS scheme.The quark bilinear operators are renormalized by computing the corresponding amputated flavornon-singlet vertex functions for different momenta on Landau gauge-fixed ensembles. We wish torenormalize light- and strange-quark bilinears, which can be written as linear combinations of thediagonal SU(3) flavor-octet and -singlet currents. In continuum schemes, with the exception of theaxial current that we do not discuss here, the renormalization of flavor singlet and non-singlet operatorsof dimension three is the same. This also appears to hold for the staggered action [37, 80]. We remarkthat, instead of extrapolating to the N f = 3 massless case, we use physical quark masses, whichmay be problematic, in particular regarding the strange quark mass. However, in Ref. [81] it wasdemonstrated that the effect of the mass-dependence is tiny for the perturbative momentum transfersthat we are interested in. Moreover, the difference is expected to vanish after a continuum limitextrapolation of a renormalized matrix element is carried out because our quark masses are tuned toa line of constant physics.Since the spin degrees of freedom are spread over hypercubes for staggered fermions, the determi-nation of the vertex function in momentum space is more challenging than for Wilson fermions. Wefollow the approach described in Ref. [82]: the taste and spin degrees of freedom are reconstructedfrom different momentum combinations. The quark propagator for a given momentum, as any vertexfunction, will be a matrix of size × , after averaging over the color degrees of freedom. Our choiceof the scalar and tensor currents, where, in the latter case, we employ a two-link operator, is detailedin Ref. [26]. β a /fm Z T Z S Table 1 . Conversion factors to the MS scheme at µ QCD = 2
GeV. – 27 – igure 9 . Multiplicative renormalization constants as a function of β . The symbols have been slightly shiftedhorizontally for better visibility and connected by lines to guide the eye. Also shown as dashed lines are theone-loop perturbative expectations [26] that will be approached as β → ∞ . The error of the final renormalization constants is dominated by systematics. On the one hand,the conversion factors from the RI’-MOM to the MS-scheme are only known up to a fixed order inperturbation theory (three loops in our case). Hence high momenta are preferable. On the otherhand, at momentum scales close to the lattice cut-off the intermediate matching to the RI’-MOMscheme will significantly be affected by lattice artifacts. Therefore, we are restricted to a “window” ofintermediate momentum values. We employ combinations along symmetric lattice directions, wherethe lattice corrections are smallest. Another complication is that due to the choice of the staggeredaction, the maximum momentum scale that can be achieved on a four-dimensional lattice is π/a ,rather than π/a . As a compromise, on the finest three lattices we interpolate the RI’-MOM resultto the fixed scale µ QCD = 2
GeV. Subsequently, this is perturbatively converted to the MS -scheme.We estimate the uncertainty by adding the difference between the scheme conversion at two- and atthree-loop order and the (statistical and systematic) interpolation uncertainty in quadrature. Thelatter contribution is negligible. At the coarsest two lattice spacings, µ QCD = 2
GeV is too close tothe cut-off scale to obtain reliable results. Therefore, at a ≈ . fm and at a ≈ . fm, we convertthe RI’-MOM results at µ QCD = 1 . GeV and at µ QCD = 1 . GeV, respectively, to the MS -schemeand evolve the result to µ QCD = 2
GeV. We replicate the same procedure at a ≈ . fm and addthe difference that we obtain at this lattice spacing between the matching at µ QCD = 2
GeV andthe matching at these lower scales in quadrature to the systematic error at µ QCD = 1 . GeV and µ QCD = 1 . GeV.The results are listed in Table 1 and shown in Fig. 9. We also include the lattice perturbativetheory one-loop expectations [26] in the figure. The comparatively larger value of Z T results in alarger modulus of the renormalized tensor coefficient.– 28 – Susceptibilities via current-current correlators
In order to circumvent the issue of flux quantization, we consider a background field that possessesnonzero momentum p in the x direction. The constant field setup will be approached via the p → limit. In this section we are working in the infinite volume, where p is a continuous variable and thislimit exists. We will discuss finite volume corrections below. This approach has been described indetail in Ref. [29] for χ b in momentum space. Here we repeat the argument in coordinate space andalso generalize it for τ fb . E.1 Magnetic susceptibility from correlators
We consider an oscillatory magnetic field and the corresponding Landau-gauge vector potential, B ( x ) = B · cos( p x ) , A ( x ) = B · sin( p x ) p . (E.1)The latter couples to i · e times the current (3.4) in the action density. We can define the associatedsusceptibility just like in Eq. (2.1), χ p , cos b = − ∂ f∂ ( eB ) (cid:12)(cid:12)(cid:12)(cid:12) B =0 = − TV (cid:90) d y d z sin( p y ) p sin( p z ) p (cid:104) j ( y ) j ( z ) (cid:105) , (E.2)where each derivative brought down an integral over the current j times the coordinate-dependenceof A and we used (cid:104) j (cid:105) = 0 . Changing the integration variable from z to x = z − y and exploitingthe translational invariance of the current-current correlator, the integrals over y , y and y can becarried out, χ p , cos b = − L (cid:90) d y d x sin( p y ) sin( p ( y + x )) p G ( x ) , (E.3)where the projected correlator G ( x ) , defined in Eq. (3.5), appears. In the p → limit, B ( x ) becomes homogeneous and χ p , cos b equals the ordinary susceptibility χ b .For reasons that will become clear in a moment, let us consider a different background field, B ( x ) = B · sin( p x ) , A ( x ) = − B · cos( p x ) p , (E.4)for which the associated oscillatory susceptibility, similarly to Eq. (E.3), reads χ p , sin b = − L (cid:90) d y d x cos( p y ) cos( p ( y + x )) p G ( x ) . (E.5)In this case the p → limit does not reproduce χ b . Instead, A ( x ) becomes homogeneous: it actsas if we had introduced a constant imaginary ‘chemical potential’ in the x direction, with magnitude µ = − eB/p . Therefore the oscillatory susceptibility becomes proportional to the leading responseto this spatial chemical potential, χ p , sin b p → −−−→ c p , c = − L (cid:90) d y d x G ( x ) . (E.6)– 29 –his detour was necessary to simplify the p → limit of the oscillatory susceptibilities. Specifi-cally, let us examine the following combination: χ p , cos b + χ p , sin b − c p = − L (cid:90) d y d x sin( p y ) sin( p ( y + x )) + cos( p y ) cos( p ( y + x )) − p G ( x ) . (E.7)This approaches χ b for p → . Using the trigonometric identity for the cosine of the difference ofangles in the numerator of the kernel reveals that the integrand is independent of y . Integrating overthis variable cancels the prefactor /L , resulting in χ b = − lim p → (cid:90) d x cos( p x ) − p G ( x ) = (cid:90) d x x G ( x ) , (E.8)where we finally performed the p → limit.In finite volumes ( x ∈ [0 , L ] ) the momentum variable p is discrete, so that the p → limit doesnot exist. Nevertheless, we can safely employ the final formula Eq. (E.8) directly in finite volumes, aslong as the linear size L is much larger than the characteristic length governing the exponential decayof G ( x ) . Symmetrizing Eq. (E.8) to comply with periodic boundary conditions and the symmetry G ( x ) = G ( L − x ) , we arrive at Eq. (3.6) of the main text.We investigated finite volume effects by comparing the bare susceptibilities obtained on × and × ensembles. The results are found to agree within statistical errors. In Fig. 10 we showthe vacuum polarization function Π( p ) of Eq. (3.7) for spatial momenta p = ( p , , , , as obtainedvia Eq. (E.8) before taking the limit p → . Here we considered the results at T ≈ MeV but theoverall picture is very similar for the other temperatures.
Figure 10 . The volume dependence of the determination of the vacuum polarization function at T ≈ MeV.We compare the × (red) and × (blue) results. The bare magnetic susceptibility can be read off fromthe intersect Π(0) . – 30 – .2 Tensor coefficient from correlators We generalize the above derivation for τ fb , which can be written as τ fb = 1 q f /e ∂∂ ( eB ) (cid:12)(cid:12)(cid:12)(cid:12) B =0 TV (cid:90) d x (cid:10) ¯ ψ f σ ψ f ( x ) (cid:11) . (E.9)Again we consider oscillatory magnetic fields of the types (E.1) and (E.4). These give rise to modulatedtensor bilinears of the forms ¯ ψ f σ ψ f ( x ) cos( p x ) and ¯ ψ f σ ψ f ( x ) sin( p x ) , respectively, which enterthe corresponding oscillatory tensor coefficients τ p , cos fb and τ p , sin fb : τ p , cos fb = iq f /e L (cid:90) d y d x cos( p y ) sin( p ( y + x )) p H f ( x ) ,τ p , sin fb = − iq f /e L (cid:90) d y d x sin( p y ) cos( p ( y + x )) p H f ( x ) , (E.10)where the projected tensor-vector correlator H f ( x ) , defined in Eq. (3.8), appears. Here we performedthe same variable substitution as in Eq. (E.3) above.In this case, τ p , cos fb approaches τ fb for p → , while τ p , sin fb vanishes in that limit. Thus we needto consider the sum of the two coefficients. Employing the trigonometric identity for the sine of thedifference of angles and carrying out the integral over y gives τ fb = lim p → (cid:104) τ p , cos fb + τ p , sin fb (cid:105) = lim p → iq f /e (cid:90) d x sin( p x ) p H f ( x ) = iq f /e (cid:90) d x x H f ( x ) . (E.11)In finite volumes we carry out the same symmetrization as in Eq. (3.6), this time taking into accountthat H f ( x ) = − H f ( L − x ) to finally arrive at Eq. (3.8) of the main text.This method to calculate τ fb is compared to the results for (cid:10) ¯ ψ f σ ψ f (cid:11) measured at B > on × lattices at T = 113 MeV in Fig. 11. For the light quarks we obtain consistent results, however,
Figure 11 . Comparison of different methods to calculate τ fb for all three flavors. Simulations at nonzero(quantized) values of the magnetic field (points) are compared with a direct determination of the slope at B = 0 (colored bands). – 31 –or τ sb the correlator tends to give values that slightly differ from the slope of a linear fit to the lowestfew available points. Since lattice artifacts and finite volume effects might be different in the twocases, such slight differences are not unexpected.In addition, we find that a linear fit to results from simulations at B > has smaller uncertaintiesthan extracting the slope at B = 0 using the correlator method. In the main text we therefore useour earlier results for (cid:10) ¯ ψ f σ ψ f (cid:11) from Ref. [26].We note that the tensor-vector correlators at nonzero spatial momenta might also be useful forextracting further features of the photon distribution amplitude. F Parametrization of the equation of state
Up to O ( B ) , the magnetic field-dependence of the complete EoS can be calculated from the magneticsusceptibility χ ( T ) . Here we provide a parametrization for this observable and also collect the relevantthermodynamical relations, which were also summarized in Ref. [28].First of all, we remind the reader that in the presence of a background magnetic field, the differentcomponents of the pressure – defined by considering an infinitesimal compression of the system in therespective direction – might become anisotropic [27]. In particular, one should distinguish between the Φ -scheme , where the flux of the magnetic field is kept constant during the compression (superscript (Φ) below), and the B -scheme , where the magnetic field strength is kept constant (superscript ( B ) ).On the one hand, the B -scheme pressure is isotropic and equals the negative of the free energy densityin the thermodynamic limit, p ( B )1 , = p = − f . (F.1)On the other hand, in the Φ -scheme the pressure components are related by the magnetization M , p (Φ)1 , = p − eB · M , M = − ∂f∂ ( eB ) . (F.2)The entropy density s and the energy density (cid:15) are scheme-independent, s = − ∂f∂T , (cid:15) = f + T s , (F.3)whereas also the interaction measure (trace anomaly) I differs between the two schemes, I ( B ) = (cid:15) − p , I (Φ) = (cid:15) − p (Φ)1 , − p = I ( B ) + 2 eB · M . (F.4)Using Eqs. (2.1) and (F.1), the leading-order expansion in the magnetic field takes the form p ( T, B ) = p ( T,
0) + χ ( T ) ( eB ) , M ( T, B ) = χ ( T ) eB . (F.5)Together with Eqs. (F.1)–(F.4) these specify the B -dependence of all relevant observables up to O ( B ) .At B = 0 the pressure is isotropic, and can be obtained from the interaction measure as p ( T, B = 0) T = (cid:90) T (cid:48) d T (cid:48) I ( T (cid:48) , B = 0) T . (F.6) We note that Eq. (F.6)) remains valid also for
B > in the B -scheme but not in the Φ -scheme. – 32 –hus, to calculate the complete EoS including B and B effects, altogether it suffices to parameterize I ( T, and χ ( T ) . For the latter we consider a parametrization of the continuum extrapolated latticeresults that smoothly approach the HRG model prediction (see Fig. 2) at low and the perturbationtheory formula (4.3)) at high temperatures. We found the following parametric form to be sufficientfor this, χ ( T ) = exp( − h /t ) · g /t + g /t + g /t g /t + g /t + g /t · β log tq , t = T GeV . (F.7)Eq. (F.7) incorporates the non-perturbative temperature-dependence predicted by the HRG model(see App. A) at low T and the logarithmic rise at high temperatures. The β coefficient is fixed toits perturbative value (2.4), while the scale q inside the logarithm is allowed to be a free parameter.The rational function involving the g i parameters interpolates between the two limiting behaviors.The so-obtained parametrization is shown in the left panel of Fig. 3 in the main text.For the interaction measure we take the parametrization of Ref. [83], I ( T, T = exp( − h /t − h /t ) · (cid:18) h + f · [tanh( f · t + f ) + 1]1 + k · t + k · t (cid:19) , t = T . GeV . (F.8)The parameters of both functions are included in Table. 2. The two parametrizations, together withthe implementations of the formulae (F.1)–(F.6) are included in the Python script param_EoS.py thatis submitted to arXiv.org together with this manuscript. β h g g g g g g q / (6 π ) h h h f f f k k Table 2 . Parameters of the functions (F.7) and (F.8).
This parametrization is valid for low magnetic fields. To be more quantitative, we compare our O ( B ) truncated results for the longitudinal pressure to the complete magnetic field-dependence fromRef. [28] for T (cid:38) MeV. We find agreement within errors in the range B/ ( πT ) (cid:46) . This upperlimit is hard-coded in the Python script as well. One final remark about the parametrization is inorder. All truncated thermodynamic observables approach zero for T → , such that a normalizationby the corresponding powers of the temperature (i.e. p /T , s/T and so on) produces sensible plots.This is not the case if O ( B ) terms are also included: at this order vacuum contributions arise andthe equation of state depends on B already at T = 0 , rendering a normalization like p /T ill-definedin the T → limit [28]. References [1] K. Kiuchi, P. Cerdá-Durán, K. Kyutoku, Y. Sekiguchi and M. Shibata,
Efficient magnetic-fieldamplification due to the Kelvin-Helmholtz instability in binary neutron star mergers , Phys. Rev.
D92 (2015) 124034, [ ].[2] L. Baiotti and L. Rezzolla,
Binary neutron star mergers: a review of Einstein’s richest laboratory , Rept.Prog. Phys. (2017) 096901, [ ]. – 33 –
3] T. Kawamura, B. Giacomazzo, W. Kastaun, R. Ciolfi, A. Endrizzi, L. Baiotti et al.,
Binary neutron starmergers and short gamma-ray bursts: effects of magnetic field orientation, equation of state, and massratio , Phys. Rev.
D94 (2016) 064012, [ ].[4] D. Kharzeev, K. Landsteiner, A. Schmitt and H.-U. Yee,
Strongly interacting matter in magnetic fields , Lect. Notes Phys. (2013) 1–624.[5] J. O. Andersen, W. R. Naylor and A. Tranberg,
Phase diagram of QCD in a magnetic field: A review , Rev. Mod. Phys. (2016) 025001, [ ].[6] V. A. Miransky and I. A. Shovkovy, Quantum field theory in a magnetic field: From quantumchromodynamics to graphene and Dirac semimetals , Phys. Rept. (2015) 1–209, [ ].[7] J. Lattimer and M. Prakash,
Neutron star structure and the equation of state , Astrophys. J. (2001)426, [ astro-ph/0002232 ].[8] D. Grasso and H. R. Rubinstein,
Magnetic fields in the early universe , Phys. Rept. (2001) 163–266,[ astro-ph/0009061 ].[9] R. Durrer and A. Neronov,
Cosmological magnetic fields: their generation, evolution and observation , Astron. Astrophys. Rev. (2013) 62, [ ].[10] D. E. Kharzeev, J. Liao, S. A. Voloshin and G. Wang, Chiral magnetic and vortical effects inhigh-energy nuclear collisions – A status report , Prog. Part. Nucl. Phys. (2016) 1–28, [ ].[11] X.-D. Ji, Gauge-invariant decomposition of nucleon spin , Phys. Rev. Lett. (1997) 610–613,[ hep-ph/9603249 ].[12] R. L. Jaffe and A. Manohar, The G problem: Fact and fantasy on the spin of the proton , Nucl. Phys.
B337 (1990) 509–546.[13] B. L. G. Bakker, E. Leader and T. L. Trueman,
A critique of the angular momentum sum rules and anew angular momentum sum rule , Phys. Rev.
D70 (2004) 114001, [ hep-ph/0406139 ].[14] B. L. Ioffe and A. V. Smilga,
Nucleon magnetic moments and magnetic properties of vacuum in QCD , Nucl. Phys.
B232 (1984) 109.[15] V. M. Belyaev and Y. I. Kogan,
Calculation of quark condensate magnetic susceptibility by QCD sumrule method , Yad. Fiz. (1984) 1035–1038.[16] I. I. Balitsky, A. V. Kolesnichenko and A. V. Yung, On vector dominance in sum rules forelectromagnetic hadron characteristics (in Russian) , Sov. J. Nucl. Phys. (1985) 178.[17] I. I. Balitsky, V. M. Braun and A. V. Kolesnichenko, Radiative decay Σ + → pγ in QuantumChromodynamics , Nucl. Phys.
B312 (1989) 509–550.[18] V. Y. Petrov, M. V. Polyakov, R. Ruskov, C. J. Weiss and K. Goeke,
Pion and photon light cone wavefunctions from the instanton vacuum , Phys. Rev. D (1999) 114018, [ hep-ph/9807229 ].[19] P. Ball, V. M. Braun and N. Kivel, Photon distribution amplitudes in QCD , Nucl. Phys.
B649 (2003)263–296, [ hep-ph/0207307 ].[20] A. Nyffeler,
Hadronic light-by-light scattering in the muon g-2: A New short-distance constraint onpion-exchange , Phys. Rev.
D79 (2009) 073012, [ ].[21] J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez,
Short-distance constraints for theHLbL contribution to the muon anomalous magnetic moment , Phys. Lett.
B798 (2019) 134994,[ ]. – 34 –
22] P. Ball and E. Kou, B → γeν transitions from QCD sum rules on the light cone , JHEP (2003) 029,[ hep-ph/0301135 ].[23] P. Colangelo, F. De Fazio and A. Ozpineci, Radiative transitions of D ∗ sJ (2317) and D sJ (2460) , Phys.Rev.
D72 (2005) 074004, [ hep-ph/0505195 ].[24] S. S. Agaev, V. M. Braun, N. Offen, F. A. Porkert and A. Schäfer,
Transition form factors γ ∗ γ → η and γ ∗ γ → η (cid:48) in QCD , Phys. Rev.
D90 (2014) 074019, [ ].[25] G. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. Katz et al.,
QCD quark condensate in externalmagnetic fields , Phys. Rev.
D86 (2012) 071502, [ ].[26] G. S. Bali, F. Bruckmann, M. Constantinou, M. Costa, G. Endrődi, S. D. Katz et al.,
Magneticsusceptibility of QCD at zero and at finite temperature from the lattice , Phys. Rev.
D86 (2012) 094512,[ ].[27] G. Bali, F. Bruckmann, G. Endrődi, F. Gruber and A. Schäfer,
Magnetic field-induced gluonic (inverse)catalysis and pressure (an)isotropy in QCD , JHEP (2013) 130, [ ].[28] G. Bali, F. Bruckmann, G. Endrődi, S. Katz and A. Schäfer,
The QCD equation of state in backgroundmagnetic fields , JHEP (2014) 177, [ ].[29] G. Bali and G. Endrődi,
Hadronic vacuum polarization and muon g − from magnetic susceptibilitieson the lattice , Phys. Rev.
D92 (2015) 054506, [ ].[30] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas,
A general method fornonperturbative renormalization of lattice operators , Nucl. Phys.
B445 (1995) 81–105,[ hep-lat/9411010 ].[31] K. G. Chetyrkin and A. Retey,
Renormalization and running of quark mass and field in theregularization invariant and MS schemes at three loops and four loops , Nucl. Phys.
B583 (2000) 3–34,[ hep-ph/9910332 ].[32] J. S. Schwinger,
On gauge invariance and vacuum polarization , Phys. Rev. (1951) 664–679.[33] G. V. Dunne, Heisenberg-Euler effective Lagrangians: Basics and extensions , in
From fields to strings:Circumnavigating theoretical physics. Ian Kogan memorial collection (3 volume set) (M. Shifman,A. Vainshtein and J. Wheater, eds.), pp. 445–522. World Scientific, Singapore, 2004. hep-th/0406216 .DOI.[34] P. Baikov, K. Chetyrkin, J. Kuhn and J. Rittinger,
Vector correlator in massless QCD at order α s andthe QED β -function at five loop , JHEP (2012) 017, [ ].[35] I. I. Balitsky and A. V. Yung,
Proton and neutron magnetic moments from QCD sum rules , Phys. Lett. (1983) 328–334.[36] J. A. Gracey,
Three loop anomalous dimension of nonsinglet quark currents in the RI’ scheme , Nucl.Phys.
B662 (2003) 247–278, [ hep-ph/0304113 ].[37] M. Constantinou, M. Hadjiantonis, H. Panagopoulos and G. Spanoudes,
Singlet versus nonsingletperturbative renormalization of fermion bilinears , Phys. Rev.
D94 (2016) 114513, [ ].[38] W. Pauli,
Über Gasentartung und Paramagnetismus , Zeitschrift für Physik (June, 1927) 81–102.[39] L. Landau, Diamagnetismus der Metalle , Zeitschrift für Physik (Sep, 1930) 629–637.[40] Y. Aoki, Z. Fodor, S. Katz and K. Szabó, The Equation of state in lattice QCD: With physical quarkmasses towards the continuum limit , JHEP (2006) 089, [ hep-lat/0510084 ]. – 35 –
41] S. Borsányi, G. Endrődi, Z. Fodor, A. Jakovác, S. D. Katz et al.,
The QCD equation of state withdynamical quarks , JHEP (2010) 077, [ ].[42] G. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. Katz et al.,
The QCD phase diagram for externalmagnetic fields , JHEP (2012) 044, [ ].[43] G. Martinelli, G. Parisi, R. Petronzio and F. Rapuano,
The proton and neutron magnetic moments inlattice QCD , Phys. Lett.
B116 (1982) 434–436.[44] M. H. Al-Hashimi and U.-J. Wiese,
Discrete accidental symmetry for a particle in a constant magneticfield on a torus , Annals Phys. (2009) 343–360, [ ].[45] G. S. Bali, F. Bruckmann, G. Endrődi and A. Schäfer,
Paramagnetic squeezing of QCD matter , Phys.Rev. Lett. (2014) 042301, [ ].[46] C. Bonati, M. D’Elia, M. Mariti, F. Negro and F. Sanfilippo,
Magnetic susceptibility of stronglyinteracting matter across the deconfinement transition , Phys. Rev. Lett. (2013) 182001, [ ].[47] C. Bonati, M. D’Elia, M. Mariti, F. Negro and F. Sanfilippo,
Magnetic susceptibility and equation ofstate of N f = 2 + 1 QCD with physical quark masses , Phys. Rev.
D89 (2014) 054506, [ ].[48] L. Levkova and C. DeTar,
Quark-gluon plasma in an external magnetic field , Phys. Rev. Lett. (2014) 012002, [ ].[49] H. B. Meyer and H. Wittig,
Lattice QCD and the anomalous magnetic moment of the muon , Prog. Part.Nucl. Phys. (2019) 46–96, [ ].[50] A. Francis, B. Jäger, H. B. Meyer and H. Wittig,
A new representation of the Adler function for latticeQCD , Phys. Rev.
D88 (2013) 054502, [ ].[51] M. Bellac,
Thermal Field Theory . Cambridge Monographs on Mathematical Physics. CambridgeUniversity Press, 2000.[52] J. Kapusta and C. Gale,
Finite-Temperature Field Theory: Principles and Applications . Cambridgemonographs on mechanics and applied mathematics. Cambridge University Press, 2006.[53] V. M. Braun et al.,
The ρ -meson light-cone distribution amplitudes from lattice QCD , JHEP (2017)082, [ ].[54] G. Endrődi, Multidimensional spline integration of scattered data , Comput. Phys. Commun. (2011)1307–1314, [ ].[55] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn,
Five-loop running of the QCD coupling constant , Phys.Rev. Lett. (2017) 082002, [ ].[56]
ALPHA collaboration, M. Bruno, M. Dalla Brida, P. Fritzsch, T. Korzec, A. Ramos, S. Schaefer et al.,
QCD Coupling from a nonperturbative determination of the three-flavor Λ parameter , Phys. Rev. Lett. (2017) 102001, [ ].[57] C. McNeile, A. Bazavov, C. T. H. Davies, R. J. Dowdall, K. Hornbostel, G. P. Lepage et al.,
Directdetermination of the strange and light quark condensates from full lattice QCD , Phys. Rev.
D87 (2013)034503, [ ].[58]
Flavour Lattice Averaging Group collaboration, S. Aoki et al.,
FLAG Review 2019 , Eur. Phys. J.
C80 (2020) 113, [ ].[59] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn,
Quark mass and field anomalous Dimensions to O ( α s ) , JHEP (2014) 076, [ ]. – 36 –
60] D. J. Broadhurst and A. G. Grozin,
Matching QCD and HQET heavy-light currents at two loops andbeyond , Phys. Rev.
D52 (1995) 4082–4098, [ hep-ph/9410240 ].[61] A. Vainshtein,
Perturbative and nonperturbative renormalization of anomalous quark triangles , Phys.Lett.
B569 (2003) 187–193, [ hep-ph/0212231 ].[62] A. Gorsky and A. Krikun,
Magnetic susceptibility of the quark condensate via holography , Phys. Rev.
D79 (2009) 086015, [ ].[63] M. Frasca and M. Ruggieri,
Magnetic susceptibility of the quark condensate and polarization from chiralmodels , Phys. Rev.
D83 (2011) 094024, [ ].[64] P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya and M. I. Polikarpov,
Chiral magnetization ofnon-Abelian vacuum: A lattice study , Nucl. Phys.
B826 (2010) 313–327, [ ].[65] V. V. Braguta, P. V. Buividovich, T. Kalaydzhyan, S. V. Kuznetsov and M. I. Polikarpov,
The ChiralMagnetic Effect and chiral symmetry breaking in SU(3) quenched lattice gauge theory , PoS
LATTICE2010 (2010) 190, [ ].[66] S. Borsányi et al.,
Is there still any T c mystery in lattice QCD? Results with physical masses in thecontinuum limit III , JHEP (2010) 073, [ ].[67] N. Bohr,
Studier over metallernes elektrontheori , Thesis, K obenhavns Universitet (translated andreprinted) (1911) .[68] Van Leeuwen, Hendrika J.,
Problèmes de la théorie électronique du magnétisme , J. Phys. Radium (1921) 361–377.[69] B. Savoie, A rigorous proof of the Bohr-van Leeuwen theorem in the semiclassical limit , Reviews inMathematical Physics (Oct, 2015) 1550019–143, [ ].[70] P. Huovinen and P. Petreczky, QCD equation of state and Hadron Resonance Gas , Nucl. Phys.
A837 (2010) 26–53, [ ].[71]
MILC collaboration, A. Bazavov et al.,
Nonperturbative QCD simulations with 2+1 flavors of improvedstaggered quarks , Rev. Mod. Phys. (2010) 1349–1417, [ ].[72] G. Endrődi, QCD equation of state at nonzero magnetic fields in the Hadron Resonance Gas model , JHEP (2013) 023, [ ].[73] T. Banks and A. Casher,
Chiral symmetry breaking in confining theories , Nucl. Phys.
B169 (1980)103–125.[74] F. Bruckmann, G. Endrődi, M. Giordano, S. D. Katz, T. G. Kovács, F. Pittler et al.,
Landau levels inQCD , Phys. Rev.
D96 (2017) 074506, [ ].[75] E. Elizalde, S. D. Odintsov, A. Romeo, A. A. Bytsenko and S. Zerbini,
Zeta regularization techniqueswith applications . World Scientific, Singapore, 1994.[76] E. S. Fraga,
Thermal chiral and deconfining transitions in the presence of a magnetic background , Lect.Notes Phys. (2013) 121–141, [ ].[77] N. Nielsen,
Asymptotic freedom as a spin effect , Am. J. Phys. (1981) 1171.[78] D. J. Toms, The Effective action at finite temperature and density with application to Bose-Einsteincondensation , cond-mat/9612003 .[79] “ NIST Digital Library of Mathematical Functions .” http://dlmf.nist.gov/ , Release 1.0.23 of2019-06-15. For the duplication formula, see https://dlmf.nist.gov/5.5.E5 . – 37 –
80] W.-J. Lee and S. R. Sharpe,
Partial flavor symmetry restoration for chiral staggered fermions , Phys.Rev.
D60 (1999) 114503, [ hep-lat/9905023 ].[81] M. Göckeler et al.,
Perturbative and nonperturbative renormalization in lattice QCD , Phys. Rev.
D82 (2010) 114511, [ ].[82] A. T. Lytle and S. R. Sharpe,
Nonperturbative renormalization for improved staggered bilinears , Phys.Rev.
D88 (2013) 054506, [ ].[83] S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabó,
Full result for the QCDequation of state with 2+1 flavors , Phys. Lett.
B730 (2014) 99–104, [ ].].