aainstitutetext: Institut für Theoretische Physik, Universität Regensburg,
D-93040 Regensburg, Germany.bbinstitutetext: Tata Institute of Fundamental
Research, Homi Bhabha Road, Mumbai 400005, India.ccinstitutetext: Institut für
Theoretische Physik, Goethe Universität Frankfurt, D-60438 Frankfurt am Main,
Germany.ddinstitutetext: Fakultät für Physik, Universität Bielefeld, D-33615
Bielefeld, Germany.

#  Magnetic susceptibility of QCD matter and its decomposition from the
lattice

Gunnar. S. Bali c,d Gergely Endrődi a and Stefano Piemonte<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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

###### Keywords:

Lattice field theory simulation, Quark Gluon Plasma, QCD Phenomenology,
Lattice QCD, Nonperturbative Effects.

##  1 Introduction

The development of a quantitative and precise understanding of the response of
QCD matter to background (electro)magnetic fields is of vital importance for
furthering our knowledge about a multitude of physical systems. Examples
include the interior of magnetars, neutron star mergers Kiuchi:2015sga ;
Baiotti:2016qnr ; Kawamura:2016nmk , off-central heavy-ion collisions and the
evolution of the universe in its early stages. For general reviews, we refer
the reader to Refs. Kharzeev:2013jha ; Andersen:2014xxa ; Miransky:2015ava . A
characteristic feature of the behavior of strongly interacting quarks and
gluons is rooted in the dependence of the QCD equation of state (EoS) on the
background magnetic field $B$. The EoS enters all the above mentioned
examples: it appears in the gravitational stability conditions of compact
stars, affecting the mass-radius relation Lattimer:2000nx ; it governs the
expansion rate in cosmological models Grasso:2000wj ; Durrer:2013pga  and it
also sets the conditions where freeze-out is reached in heavy-ion collisions,
see, e.g., Ref. Kharzeev:2015znc .

While the equilibration of fireballs produced in heavy-ion collisions is still
a subject of research, at least in astrophysical systems the time and distance
scales over which the magnetic field varies are much larger than those that
govern QCD processes that affect, e.g., the EoS or nucleo-synthesis. With
these applications in mind, solving QCD in a constant background magnetic
field is sufficient and this scenario is amenable to lattice simulations.

The leading dependence of the EoS on $B$ is encoded in the magnetic
susceptibility $\chi$ of QCD matter. Its sign distinguishes between
paramagnets ($\chi>0$), for which the exposure to the background field is
energetically favorable, and diamagnets ($\chi<0$), which repel the external
field. In QCD matter, like in any other material, the origin of the
magnetization and hence of the magnetic susceptibility is related to the spin
and angular momentum of charged particles. At high temperatures the quarks are
the relevant degrees of freedom; at low temperatures the hadrons and in
particular the pions take over, while (valence and sea) quarks contribute just
as their fundamental constituents.

The total angular momentum that gives rise to the magnetization can be
decomposed into contributions from the spins of the quarks of different
flavors and a remainder. The latter contains the quark orbital angular momenta
but also the angular momentum of the gluons, that can split into quark-
antiquark pairs.111This situation is analogous to the decomposition of the
nucleon spin. In particular, vacuum expectation values of the same local
operators appear in the magnetic field background at zero momentum as in the
decomposition Ji:1996ek  of the transversely polarized generalized parton
distribution functions of deep inelastic scattering at leading twist. While
the individual quark spin contributions in both cases are unique and gauge
invariant, the further decomposition of the remainder into quark and gluon
parts is ambiguous: the decomposition of Ref. Ji:1996ek  is based on the
Belinfante-Rosenfeld form 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 Jaffe:1989jz ; Bakker:2004ib .
The quark spin contribution to the magnetic susceptibility is due to the
expectation value
$\left\langle\bar{\psi}_{f}\sigma_{\mu\nu}\psi_{f}\right\rangle$, whose
leading order response is linear in the magnetic field strength tensor
$F_{\mu\nu}$. Depending on the normalization, the slope is proportional to the
magnetic susceptibility of the quark condensate Ioffe:1983ju ; Belyaev:1984ic
; Balitsky:1985aq , or the so-called tensor coefficient, i.e. the
normalization of the photon lightcone distribution amplitude (DA)
Balitsky:1989ry ; Petrov:1998kg ; Ball:2002ps  of finding a quark-antiquark
pair of flavor $f$ in a transverse photon. This in itself appears in a
multitude of applications, e.g., as a correction to the hadronic light-by-
light scattering contribution to the muon anomalous magnetic moment
Nyffeler:2009tw ; Bijnens:2019ghy , within radiative transitions Ball:2003fq ;
Colangelo:2005hv  and in the photo-production of mesons Agaev:2014wna .

Some of the present authors have already addressed several of the above
aspects in Refs. Bali:2012zg ; Bali:2012jv ; Bali:2013esa ; Bali:2014kia .
Here we improve on these studies by employing a novel calculational method
that is based on Ref. Bali:2015msa , by carrying out the QCD renormalization
non-perturbatively with respect to the intermediate RI’-MOM scheme
Martinelli:1994ty ; Chetyrkin:1999pq  and by adding a finer lattice spacing.
In addition we present and exploit new analytical findings. One of the
outcomes will be that in the strongly interacting medium at low to moderately
high temperatures the quark spin-related susceptibility is negative
(diamagnetism) while the part that is due to the orbital angular momentum is
positive. Clearly, this behavior is very different from the response to
magnetic fields of the materials that have so far been accessible to solid
state physics experiments. Therefore, our results offer a glimpse into a
completely new regime of spin physics.

The article is organized as follows. In Sec. 2 we introduce our notations and
the central observables. For conceptual clarity, we carefully address their
divergence structure and renormalization in QED and QCD. In Sec. 3 we then
discuss details of the simulation and, in particular, we introduce our new
method that employs current-current correlators in a mixed coordinate- and
momentum-space representation. This enables us to determine susceptibilities
from lattice simulations at $B=0$. We then present and discuss our results in
Sec. 4, before we summarize. We include several technical appendices: in App.
A we investigate the effects of taste splitting in the staggered formulation
within the hadron resonance gas model. This turns out to be important to avoid
underestimating the systematics of the continuum limit extrapolation. In App.
B we derive the factorization of the susceptibility into quark spin-related
and other contributions, building upon earlier partial results Bali:2012jv .
In the extensive App. C, several derivations are carried out for the free
case, establishing, e.g., the structure of QED divergencies. App. D discusses
the non-perturbative renormalization procedure. In App. E we present more
detail on the derivation of the new current-current method and compare to
numerical results, obtained using conventional background field approaches.
Finally, App. F gives a parametrization of our results for the QCD EoS for a
broad range of temperatures and magnetic field strengths. The corresponding
Python script param_EoS.py is uploaded to the arXiv as ancillary file along
with this paper.

##  2 The response of QCD matter to background fields and the magnetic
susceptibility

Without any loss of generality, below we consider a magnetic field pointing in
the $x_{3}$ direction, with the magnitude $B$. The magnetic susceptibility is
defined via the leading (quadratic) dependence of the QCD free energy density
$f$ on the field strength,

$\chi_{b}=-\left.\frac{\partial^{2}f}{\partial(eB)^{2}}\right|_{B=0},\quad\quad
f=-\frac{T}{V}\log\mathcal{Z}\,,$ (1)

where $\mathcal{Z}$ is the partition function, $T$ the temperature and $V$ the
spatial volume of the system. The product $eB$ of the elementary electric
charge $e$ and the magnetic field is a renormalization group invariant due to
the QED vector Ward identity, see Eq. (3) below. Therefore, $\chi_{b}$ is free
of multiplicative renormalization. Still, the susceptibility undergoes
additive renormalization, which is made explicit by the index $b$, indicating
the bare quantity, as obtained in the lattice scheme at a lattice spacing
$a$.222Below we will indicate quantities that are subject to QED
renormalization with the subscript $b$. This was discussed in Ref.
Bali:2014kia  in depth, but we repeat the argument here for comprehensiveness.

###  2.1 QED renormalization

The total free energy density of the system, that includes both QCD matter and
the classical background field,

$f_{\rm tot}=f+\frac{B_{b}^{2}}{2}\,,$ (2)

is a physical observable and therefore free of divergences. However, the
second term within Eq. (2), involving the bare magnetic field $B_{b}$,
contains a logarithmic divergence in the lattice spacing $a$ due to electric
charge renormalization Schwinger:1951nm ,

$B_{b}^{2}=Z_{e}B^{2},\quad\quad e_{b}^{2}=Z_{e}^{-1}e^{2},\quad\quad
eB=e_{b}B_{b},\quad\quad Z_{e}=1+\beta_{1}(a^{-1})\,e^{2}\,\log\left(\mu_{\rm
QED}^{2}a^{2}\right)\,,$ (3)

where the renormalized quantities $e$ and $B$ depend on the QED
renormalization scheme and on the QED renormalization scale $\mu_{\rm QED}$.
Since the background field is classical, only the leading-order QED
$\beta$-function coefficient $\beta_{1}$ appears here Dunne:2004nc . Note that
$\beta_{1}$ is affected by QCD corrections at the cut-off scale,

$\beta_{1}(a^{-1})=\beta_{1}\cdot\Bigg{[}1+\sum_{i\geq
1}c_{i}\,g^{2i}(a^{-1})\Bigg{]}\xrightarrow{a\to
0}\beta_{1},\quad\quad\beta_{1}=\frac{1}{4\pi^{2}}\cdot\sum_{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 $\overline{\mathrm{MS}}$ and $\mathrm{MOM}$ schemes
Baikov:2012zm  and $c_{1}=1/(4\pi^{2})$ is universal for massless schemes.
These QCD corrections333Note that in general disconnected diagrams also start
to contribute for $i\geq 3$, but in the present case these vanish because we
are dealing with the three lightest quark flavors and $\sum_{f}q_{f}/e=0$.
vanish logarithmically with the lattice spacing towards the continuum limit
due to the asymptotically free nature of the strong interactions, as is also
indicated in Eq. (4).

Eq. (2) implies that $f$ contains the same additive divergence as
$B_{b}^{2}/2$, but with an opposite sign. This propagates into the
susceptibility (1), resulting in

$\chi_{b}=\chi[\mu_{\rm QED}]+\beta_{1}(a^{-1})\,\log(\mu_{\rm
QED}^{2}a^{2})\,.$ (5)

The renormalized susceptibility $\chi$ depends on the renormalization scale,
which we indicated here explicitly in square brackets. We confirm the presence
of the logarithmic divergence in $\chi_{b}$ analytically for the free case in
App. C.2 and numerically for full QCD in Sec. 4.1. The divergence is
independent of the temperature so that it cancels within the difference

$\chi\equiv\chi[\mu_{\rm QED}^{\rm phys}]=\chi_{b}(T)-\chi_{b}(T=0)\,.$ (6)

This definition of the renormalized susceptibility – implying that it vanishes
identically at $T=0$ – corresponds to a particular choice of the
renormalization scale $\mu_{\rm QED}=\mu_{\rm QED}^{\rm phys}$. In fact this
is the only prescription that adheres to the physical requirement that the
magnetic permeability $(1-e^{2}\chi)^{-1}$ should be unity in the vacuum. In
the following we suppress the dependence on the QED renormalization scale and
simply write $\chi$ for the susceptibility, renormalized in this way.

###  2.2 The tensor coefficient

Besides the magnetic susceptibility there exist further quantities that
characterize the leading-order response of the QCD medium to the background
magnetic field. For a general background field $F_{\mu\nu}$, the fermion
bilinear involving the relativistic spin operator $\sigma_{\mu\nu}$ develops a
nonzero expectation value Ioffe:1983ju ; Balitsky:1983xk ,

$\left\langle\bar{\psi}_{f}\sigma_{\mu\nu}\psi_{f}\right\rangle=q_{f}F_{\mu\nu}\cdot\tau_{fb}+\mathcal{O}(F^{3})\,,\quad\quad\sigma_{\mu\nu}=\frac{1}{2i}[\gamma_{\mu},\gamma_{\nu}]\,.$
(7)

We will refer to $\tau_{fb}$ as the tensor coefficient for the flavor $f$.
Similarly to $\chi_{b}$, this is also a bare observable that contains additive
logarithmic divergences in the cut-off. For our choice of direction of the
magnetic field the tensor coefficient can be determined as the slope of the
expectation value of the fermion bilinear involving $\sigma_{12}$ at small
values of the magnetic field:

$\tau_{fb}=\lim_{B\to
0}\frac{\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle}{q_{f}B}\,.$
(8)

The tensor coefficient contains a similar logarithmic divergence as
$\chi_{b}$. We demonstrate the reason for this in Sec. 2.3 below. In
particular, the divergence structure takes the form Bali:2012jv ,

$\tau_{fb}=\tau_{f}+\gamma_{1}^{\tau}(a^{-1})\,m_{f}\log(\mu_{\rm
QED}^{2}a^{2})\,,\quad\quad\gamma_{1}^{\tau}(a^{-1})=\gamma^{\tau}_{1}\cdot\left[1+\mathcal{O}(g^{2}(a^{-1}))\right]\,,\quad\;\gamma_{1}^{\tau}=\frac{3}{4\pi^{2}}\,,$
(9)

where $m_{f}$ is the mass of the quark of flavor $f$. Again, due to asymptotic
freedom, QCD corrections to $\gamma_{1}^{\tau}$ vanish in the continuum
limit.444Unlike for $\beta_{1}(a^{-1})$ in Eq. (4), the order $g^{2}$
perturbative coefficient is not known in this case. Therefore, any definition
of a renormalized quark mass $m_{f}$ is valid to this order and we use the
lattice quark mass. Eq. (9) is confirmed in App. C.4 for the free case and
checked numerically in full QCD in Sec. 4.2. Notice that in the chiral limit
the divergent term disappears, so that $\lim_{m_{f}\to 0}\tau_{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_{\gamma}^{\perp}$ of the leading-
twist photon distribution amplitude Balitsky:1989ry ; Petrov:1998kg ;
Ball:2002ps , i.e. of the infinite momentum frame probability amplitude that a
real photon dissociates into a quark-antiquark pair of flavor $f$.

In analogy to Eq. (6), we define the renormalized tensor coefficient by
subtracting its value at zero temperature,

$\tau_{f}=\tau_{fb}(T)-\tau_{fb}(T=0)\,,$ (10)

which again corresponds to a particular choice of the QED renormalization
scale. Unlike for $\chi$, there is no preferred choice in this case, but it is
natural to use the same prescription as in Eq. (6) above.

Besides the (QED-related) additive renormalization detailed above, the tensor
coefficient also undergoes (QCD-related) multiplicative renormalization by the
tensor renormalization constant $Z_{T}$. This introduces a further scheme- and
scale-dependence of this observable. Below we will consider $Z_{T}$ in the
$\overline{\rm MS}$ scheme at the QCD renormalization scale $\mu_{\rm
QCD}=2\textmd{ GeV}$. Unlike in our previous study Bali:2012jv , 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 Martinelli:1994ty ;
Chetyrkin:1999pq  and subsequently translate the result at three-loop order
Gracey:2003yr  into the $\overline{\rm MS}$ scheme. This procedure is detailed
in App. D. We remark that $Z_{T}$ is independent of the temperature, thus the
ordering of the QED renormalization (i.e. the $T=0$ subtraction) and the QCD
renormalization (multiplication by $Z_{T}$) is irrelevant for the
determination of the renormalized tensor coefficient $\tau_{f}$.

###  2.3 Decomposition into spin and orbital angular momentum contributions

One might suspect that $\chi$ and $\tau_{f}$ are not completely unrelated.
Indeed, as we first discussed in Ref. Bali:2012jv , $\tau_{f}$ represents the
contribution of the spin of the quark flavor $f$ to the total magnetic
susceptibility. In particular, $\chi$ can be decomposed into spin-related and
orbital angular momentum-related contributions,555Note that for simplicity we
refer to $\chi^{\rm ang}$ as the orbital angular momentum contribution. In the
interacting case this can be further factorized, separating out the gluon
total angular momentum contribution Bali:2013esa  from those of the quark
angular momenta.

$\chi=\chi^{\rm spin}+\chi^{\rm ang}\,,$ (11)

and the spin term is related to the tensor coefficients as

$\chi^{\rm spin}=\sum_{f}\frac{(q_{f}/e)^{2}}{2m_{f}}\left[\tau_{f}(m_{f}^{\rm
val})-\tau_{f}(m_{f}^{\rm val}\to 0)\right]\cdot Z_{T}Z_{S}\,,$ (12)

where $m_{f}$ denotes the quark mass in the lattice scheme and $Z_{S}$ and
$Z_{T}$ are the scalar and tensor renormalization constants, respectively. The
second term in the square brackets is understood to correspond to the limit of
a vanishing valence quark mass taken at physical, i.e. nonzero, values of the
sea quark masses. We discuss the difference between valence and sea quark
masses in Sec. 3 below. In App. B we prove Eqs. (11)–(12) and illustrate the
origin of the subtraction of the chiral valence quark limit. App. C also
contains an explicit check of Eq. (12) in the free case.

A remark regarding the choice of renormalization scales is in order here. Eq.
(12) is chosen so that $\chi^{\rm spin}$ vanishes at $T=0$. Recall however
that, according to our remark below Eq. (10), we are free to choose an
arbitrary QED renormalization scale for the renormalized tensor coefficient.
This freedom propagates into $\chi^{\rm spin}$ and – through the decomposition
(11) – to $\chi^{\rm ang}$ as well. In contrast, the renormalization scale for
$\chi$ is fixed by the requirement $\chi(T=0)=0$. Thus, in principle both
susceptibility contributions may be shifted by an arbitrary amount, as long as
their sum remains zero at $T=0$. We will follow the choice made in Eq. (12),
which corresponds to setting $\chi^{\rm spin}(T=0)=\chi^{\rm ang}(T=0)=0$. In
the free case this is realized by choosing one and the same QED
renormalization scale for all susceptibility contributions, see App. C.3. Our
numerical results in full QCD below suggest that also in the interacting case
the QED scale that corresponds to the renormalization condition $\chi^{\rm
spin}(T=0)=0$ is consistent with the one obtained from setting $\chi(T=0)=0$.

In Eq. (12) we also carried out the QCD related multiplicative renormalization
by including the tensor renormalization constant $Z_{T}$ required for
$\tau_{f}$ (as mentioned above) as well as the scalar renormalization constant
$Z_{S}=Z_{m}^{-1}$, which multiplies the inverse quark mass.666Like in
massless continuum schemes, in staggered lattice formulations there is no
difference between singlet and non-singlet renormalization factors for these
quark bilinears. This was explicitly demonstrated at order $g^{4}$ in Ref.
Constantinou:2016ieh . Note that these renormalization factors depend on the
QCD regularization scheme and on the QCD renormalization scale. As mentioned
above, we choose the $\overline{\mathrm{MS}}$ scheme and $\mu_{\rm
QCD}=2\textmd{ GeV}$. We remark that, just as for $\tau_{f}$, the ordering of
the QED and the QCD renormalization is irrelevant for $\chi^{\rm spin}$. We
stress again that while the factorization of the total susceptibility $\chi$
into $\chi^{\rm spin}$ and $\chi^{\rm ang}$ depends on the QCD scheme and
scale, $\chi$ itself is a QCD renormalization group invariant.

In the free case the two susceptibility contributions have a constant ratio,
$\chi^{\rm spin}:\chi^{\rm ang}=3:(-1)$, reflecting the well-known response of
a free charged fermion to the magnetic field via its spin and its orbital
angular momentum, dating back to Pauli and Landau Pauli:1927:GPG ; Landau1930
. This ratio, which translates into the rule $\chi^{\rm spin}:\chi=3:2$, holds
identically in the free case, see App. C.3. In contrast, in full QCD it only
applies to the divergence structure, which in the continuum limit – as we have
seen above – is governed by pure QED physics. Below we determine to what
extent the Pauli-Landau decomposition is affected by the strong interactions.

##  3 Lattice methods

We consider spatially symmetric $N_{s}^{3}\times N_{t}$ lattice ensembles,
corresponding to the temperature $T=1/(N_{t}a)$ and the volume
$V=L^{3}=(N_{s}a)^{3}$. The simulations are performed with the tree-level
Symanzik improved gauge action $S_{g}$ and three flavors ($f=u,d,s$) of stout
smeared rooted staggered quarks Aoki:2005vt , described by the Dirac
operator777Due to the electromagnetic charge $q_{f}$, the covariant derivative
depends on the quark flavor $f$.  $\not{D}_{f}+m_{f}$. The quark masses
$m_{f}$ are tuned as a function of the inverse gauge coupling $\beta=6/g^{2}$
along the line of constant physics: $m_{ud}(\beta)\equiv
m_{u}(\beta)=m_{d}(\beta)=m_{s}(\beta)/R$ with $R=28.15$ Borsanyi:2010cj . The
electric charges are set as $q_{d}=q_{s}=-q_{u}/2=-e/3$, where $e>0$ is the
elementary electric charge. The magnetic field enters in $\not{D}_{f}$ via
space-dependent $\mathrm{U}(1)$ phases. Further details of our setup and of
the simulation algorithm are discussed in Ref. Bali:2011qj . The lattice
geometries for our finite temperature lattices are $16^{3}\times 6$,
$24^{3}\times 6$, $24^{3}\times 8$, $28^{3}\times 10$ and $36^{3}\times 12$,
allowing for the investigation of both finite volume and discretization
effects. Our zero-temperature ensembles consist of $24^{3}\times 32$,
$32^{3}\times 48$ and $40^{3}\times 48$ lattices.

In the rooted staggered formulation the partition function and the expectation
value of the tensor bilinear are written as path integrals over the
$\mathrm{SU}(3)$-valued gluonic links $U$ as

$\begin{split}\mathcal{Z}&=\int\mathcal{D}U\,e^{-\beta
S_{g}}\prod_{f^{\prime}=u,d,s}\left[\det(\not{D}_{f^{\prime}}+m_{f^{\prime}}^{\rm
sea})\right]^{1/4}\,,\\\
\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle&=\frac{T}{V}\frac{1}{\mathcal{Z}}\int\mathcal{D}U\,e^{-\beta
S_{g}}\;\textmd{tr}\,\frac{\sigma_{12}}{\not{D}_{f}+m_{f}^{\rm
val}}\prod_{f^{\prime}=u,d,s}\left[\det(\not{D}_{f^{\prime}}+m_{f^{\prime}}^{\rm
sea})\right]^{1/4}\,.\\\ \end{split}$ (13)

Here we distinguished between two different types of masses: the sea quark
masses $m_{f}^{\rm sea}$ which appear in the fermion determinant and thus
affect the generation of gluonic configurations; and the valence quark mass
$m_{f}^{\rm val}$, which enters in the operator and thereby affects the
measurement on a given set of configurations. For usual observables both
masses are equal and set according to the line of constant physics,
$m_{f}^{\rm sea}=m_{f}^{\rm val}=m_{f}$. The spin contribution to the magnetic
susceptibility is exceptional in this sense – as pointed out above in Eq.
(12), it also involves the value of
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$ in the limit
$m_{f}^{\rm val}\to 0$ but keeping $m_{f}^{\rm sea}=m_{f}$. Note that this
does not mean that we are dealing with a non-unitary theory but merely that
$\chi^{\rm spin}$ can be expressed as a difference of two expectation values
involving $\tau_{f}$ at different valence quark mass values, see App. B.

###  3.1 Magnetic flux quantization

In an infinite volume a magnetic field pointing in the $x_{3}$ direction can
be generated by the Landau-gauge electromagnetic potential

$A_{2}=Bx_{1}\,.$ (14)

In a finite volume, in order to comply with periodic boundary conditions for
the electromagnetic parallel transporters $u_{\mu f}=\exp(iq_{f}A_{\mu})$, the
boundary twist term $A_{1}=-Bx_{2}L\,\delta(x_{1}-L)$ needs to be included as
well Martinelli:1982cb . In this setup the flux of the magnetic field is
quantized according to AlHashimi:2008hr

$eB=6\pi N_{B}/L^{2},\quad\quad N_{B}\in\mathbb{Z}\,,$ (15)

so that a differentiation with respect to $eB$ – as required in Eq. (1) – is
not possible in a standard manner. Several methods were developed to overcome
this problem on the lattice, including the anisotropy method Bali:2013esa ;
Bali:2013owa , the finite difference method Bonati:2013lca ; Bonati:2013vba
and the generalized integral method Bali:2014kia . These are all based on
approximating the derivative numerically using finite differences in the
integer variable $N_{B}$. This requires independent simulations using several
values of $N_{B}$, which increases the computational requirements
considerably. Furthermore, an extrapolation $N_{B}\to 0$ becomes necessary,
which inevitably introduces systematic uncertainties. An alternative approach
is the half-half method Levkova:2013qda , which employs a magnetic field
profile that is positive in one half and negative in the other half of the
lattice – this enables taking the derivative with respect to the amplitude of
the field analytically. However, finite volume effects are substantially
enhanced in this case due to the discontinuity of the background field at the
boundaries Bali:2015msa  – even if such effects are expected to cancel in
temperature differences Ludmilla .

###  3.2 Determining the susceptibility via current-current correlators

In view of the above, a method is desirable that only involves measurements at
$B=0$, thereby circumventing the flux quantization problem of the constant
background field profile. In Ref. Bali:2015msa  we demonstrated that at zero
temperature, $\chi_{b}$, as defined in (1), is related to a mixed-
representation two-point function of the electromagnetic current, and can thus
be measured at $B=0$. Below we motivate this method and clarify how it extends
to nonzero temperatures. A detailed derivation can be found in App. E and, for
$T=0$, in Ref. Bali:2015msa .

Before integrating out the fermions in the path integral (13), the vector
potential (14) couples to $i\cdot e$ times the $\mu=2$ component of the
electromagnetic current,

$j_{\mu}=\sum_{f}\frac{q_{f}}{e}\,\bar{\psi}_{f}\gamma_{\mu}\psi_{f}\,.$ (16)

in the action density. Taking derivatives of $\log\mathcal{Z}$ with respect to
$eB$ therefore brings down integrals over the current $j_{2}$ times the
$x_{1}$-dependent term $i\,\partial A_{2}/\partial B$. Thus we can anticipate
the result to take the form of a convolution of the projected correlator,

$G(x_{1})=\int\textmd{d}x_{2}\,\textmd{d}x_{3}\,\textmd{d}x_{4}\,\left\langle
j_{2}(x)j_{2}(0)\right\rangle\,,$ (17)

with an $x_{1}$-dependent kernel.

Instead of directly using the gauge (14), it is instructive to approach the
constant magnetic field background via oscillatory fields that possess nonzero
momentum $p_{1}$ in the $x_{1}$ direction. Using these profiles, we can take
the thermodynamic limit and subsequently the $p_{1}\to 0$ limit. This approach
reveals that the magnetic susceptibility (1) arises as a smooth limit of
susceptibilities with respect to oscillatory fields. As the details of the
derivation are somewhat technical, we delegate them to App. E and only quote
the main results here.

In the thermodynamic limit, where the momentum variable is continuous and the
$p_{1}\to 0$ limit can be taken, the susceptibility is obtained as

$\chi_{b}=-\lim_{p_{1}\to
0}\int\textmd{d}x_{1}\,\frac{\cos(p_{1}x_{1})-1}{p_{1}^{2}}\,G(x_{1})=\int\textmd{d}x_{1}\,\frac{x_{1}^{2}}{2}\,G(x_{1})\,.$
(18)

In finite volumes ($x_{1}\in[0,L]$) the momentum variable $p_{1}$ is discrete,
so that the $p_{1}\to 0$ limit does not exist. Nevertheless, we can safely
employ the formula (18) directly in finite volumes, as long as the linear size
$L$ is much larger than the characteristic length governing the exponential
decay of $G(x_{1})$. Symmetrizing Eq. (18) to comply with periodic boundary
conditions and the symmetry $G(x_{1})=G(L-x_{1})$, we arrive at

$\chi_{b}=\frac{1}{2}\int_{0}^{L}\textmd{d}x_{1}\,G(x_{1})\cdot\begin{cases}x_{1}^{2},&x_{1}\leq
L/2\\\ (x_{1}-L)^{2},&x_{1}>L/2\end{cases}\,.$ (19)

In the representation (19) the current-current correlator is computed in
coordinate space. Only afterwards a Fourier transformation is carried out via
the convolution with the quadratic kernel in order to represent the constant
background field. In this way the problem of flux quantization is avoided.
Notice that there is a remnant of flux quantization in the formula (19),
signaled by the cusp in the kernel at $x_{1}=L/2$. However, this cusp has no
practical relevance, as in the integral it is multiplied by $G(L/2)$, which is
exponentially small. Thus we do not expect to encounter substantial finite
volume effects. This is contrary to the case of the half-half method
Levkova:2013qda , where translational invariance is broken already on the
level of the expectation values, involving a vector potential with kinks.
Nevertheless, we investigate the finite volume effects of the new method
numerically in Sec. 4.

The result (18) can be recast into an alternative form using the vacuum
polarization tensor

$\Pi_{\mu\nu}(p)=\int\textmd{d}^{4}x\,e^{ipx}\left\langle
j_{\mu}(x)j_{\nu}(0)\right\rangle\,,\quad\quad\Pi_{22}(p=\\{p_{1},0,0,0\\})=-p_{1}^{2}\,\Pi(p^{2})\,.$
(20)

The second relation, involving the vacuum polarization form factor $\Pi$, only
holds for this specific choice of spatial indices.888At zero temperature, the
second relation of Eq. (20) follows directly from the decomposition
$\Pi_{\mu\nu}(p)=(p_{\mu}p_{\nu}-\delta_{\mu\nu}p^{2})\,\Pi(p^{2})$. For $T>0$
the Lorentz structure of $\Pi_{\mu\nu}(p)$ is more complicated so that, in
addition to $\Pi$, a form factor $\Pi_{L}$ appears bellac2000thermal ;
kapusta2006finite . However, in the static case ($p_{4}=0$) only $\Pi$
contributes to the spatial components of $\Pi_{\mu\nu}$ so that the second
relation of Eq. (20) continues to hold. Employing these definitions, we can
rewrite

$\chi_{b}=\lim_{p_{1}\to
0}\,\Pi(p^{2})=\Pi(0),\quad\quad\Pi(p^{2})=\int\textmd{d}x_{1}\,\frac{1-\cos(p_{1}x_{1})}{p_{1}^{2}}\,G(x_{1})\,,$
(21)

where we used that the imaginary part of $\Pi(p^{2})$ vanishes.

We note that the vacuum polarization function $\Pi$ has been the subject of
intense research as it appears in the hadronic contribution to the muon
anomalous magnetic moment, see, e.g., the recent review Meyer:2018til . In
that setting the relevant observable is the second moment of the two-point
function of the electromagnetic current (16), projected to zero spatial
momentum Francis:2013fzp . Exchanging the time coordinate $x_{4}$ for the
spatial coordinate $x_{1}$, one can obtain $\chi_{b}$ in an analogous way in
our background field setup Bali:2015msa . The two determinations are
equivalent at zero temperature. For $T>0$ it is important to use spatial
momenta, i.e. a kernel involving spatial coordinates for $\chi_{b}$, since
this encodes the magnetic response.

In App. E we derive a similar representation for $\tau_{fb}$ as well. In this
case the equivalents of Eqs. (17) and (19) become999On general grounds, the
linear response of the expectation value of an $n$-point function with respect
to a background field can always be obtained by computing $(n+1)$-point
functions in the vacuum. Usually, the former method is favorable in terms of
the statistical noise. However, the latter option exempts us from the need of
generating additional gauge ensembles with non-vanishing values of the
background field. As already discussed, in the present context, we also
circumvent the issue of flux quantization.

$\tau_{fb}=\frac{i}{q_{f}/e}\int_{0}^{L}\\!\textmd{d}x_{1}\,H_{f}(x_{1})\cdot\begin{cases}x_{1},&x_{1}\leq
L/2\\\ x_{1}-L,&x_{1}>L/2\end{cases}\,,\quad
H_{f}(x_{1})=\int\\!\textmd{d}x_{2}\,\textmd{d}x_{3}\,\textmd{d}x_{4}\,\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}(x)j_{2}(0)\right\rangle\,.$
(22)

We remark that the second moment of the photon DA is accessible too, replacing
$\sigma_{12}$ by combinations of
$\sigma_{\mu\nu}\left(\overleftarrow{D}_{\rho}\overleftarrow{D}_{\sigma}+\overrightarrow{D}_{\rho}\overrightarrow{D}_{\sigma}-2\overleftarrow{D}_{\rho}\overrightarrow{D}_{\sigma}\right)$
that are antisymmetric in indices equal to $1$ and $2$, symmetrized over all
other non-trivial combinations of indices and with all traces subtracted, see,
e.g., Ref. Braun:2016wnx . This is beyond the scope of the present work.

In summary, via the relation (19) we are able to determine the magnetic
susceptibility using direct measurements at $B=0$. This is certainly
advantageous over calculating the free energy density (which cannot be
obtained as a simple expectation value) at nonzero magnetic fields and
differentiating it numerically. The similar relation for the tensor
coefficient, Eq. (22), might also be used to avoid measuring
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$ at $B>0$.
However, since the latter is a simple one-point function, the gain is not
obvious in this case. In App. E we compare the two methods for this observable
and conclude that the correlator method indeed gives larger statistical
errors. Therefore, we opted to use our earlier results Bali:2012jv  for
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$.

##  4 Results

First we demonstrate that – in accordance with Eq. (21) – $\chi_{b}=\Pi(0)$
arises as a smooth limit of the vacuum polarization function $\Pi(p^{2})$ at
spatial momenta. To this end we calculate the correlator $G(x_{1})$ using
$\mathcal{O}(1000)$ random sources located on three-dimensional $x_{1}$-slices
of our lattices, taking into account both connected and disconnected
contributions. The correlator is shown in the left panel of Fig. 1 for our
$N_{t}=6$ lattices at a high temperature $T\approx 176\textmd{ MeV}$. Here we
compare two different volumes with $N_{s}=24$ and $N_{s}=16$, revealing that
finite size effects in the exponential fall-off are tiny. Subsequently,
$G(x_{1})$ is convoluted with the kernels of Eqs. (19) and (21) to obtain
$\Pi(0)$ and $\Pi(p^{2})$, respectively. In the right panel of Fig. 1 we show
$\Pi(p^{2})$ for low momenta, again at the same temperature $T\approx
176\textmd{ MeV}$. As expected, the zero-momentum limit is approached smoothly
and the two volumes are found to agree perfectly.

Figure 1: Left panel: comparison of the absolute value of the current-current
correlator for two different volumes, $24^{3}\times 6$ (red) and $16^{3}\times
6$ (blue). Filled (open) points indicate positive (negative) values. Right
panel: the vacuum polarization function at spatial momenta for $T\approx
176\textmd{ MeV}$ using two different volumes. The bare magnetic
susceptibility can be read off from the intersect $\Pi(0)$.

Figure 2: Left panel: the susceptibility obtained via a truncation of Eq. (19)
for two different volumes, $24^{3}\times 6$ (red) and $16^{3}\times 6$ (blue).
The inset zooms into the region near $x_{1}^{\rm max}=N_{s}a/2$. Right panel:
relative error of the susceptibility as a function of the number of employed
noisy estimators.

Next we investigate finite volume effects in more detail. In particular, we
truncate the convolution (19) at $x_{1}^{\rm max}\leq L/2$ and plot the so-
obtained truncated susceptibility in the left panel of Fig. 2. This sheds more
light on why volume effects are so small. While for the smaller volume, the
exponential decay is cut off at a lower $x_{1}$, the slight enhancement of $G$
around $L/2$ due to the backward propagating exponential (see the left panel
of Fig. 1) almost completely corrects for this. Finally, we estimated the
deviation of the result from the thermodynamic limit by considering a single-
exponential fit of the correlator at $x_{1}<L/2$ and performing the
convolution (18) for $L/2\leq x_{1}<\infty$. For $N_{s}=16$ this correction is
found to be about half of the statistical error of $\chi_{b}$, while for
$N_{s}=24$ it is found to be two orders of magnitude smaller than that. For
these analyses we considered the results at $T\approx 176\textmd{ MeV}$, where
we have very precise data. For the lower temperatures our data are noisier;
here we find the finite volume errors to be significantly smaller than our
statistical uncertainties already for $N_{s}=16$.

Before turning to the main results, we discuss the statistical accuracy and
the numerical costs of the present method. The right panel of Fig. 2 shows the
relative error of $\chi_{b}$ at $T\approx 113\textmd{ MeV}$ on our
$24^{3}\times 6$ lattices as a function of the employed number $N_{\rm vec}$
of noisy estimators. The evaluation of $G(x_{1})$ requires two inversions for
the light quarks and two for the strange quark for each noisy estimator. The
figure reveals that sub-percent errors can be achieved. For sufficiently high
$N_{\rm vec}$ the connected contributions are found to dominate the error, as
already recognized in Ref. Bali:2015msa . Alternative methods to calculate
$\chi_{b}$ Bonati:2013lca ; Bonati:2013vba ; Bali:2014kia  require several
independent simulations at nonzero $B$ and a reconstruction of
$\log\mathcal{Z}$ for each magnetic field and are therefore much more
expensive than the present approach (of course, in that case the physical
$B>0$ ensembles might also be used for other purposes). To be specific, we
consider our results Bali:2014kia  using the integral method at the same
temperature and lattice spacing as above. In that case we needed to perform
around 40 independent simulations (at nonzero $B$ as well as at different
quark masses), generating several hundred decorrelated configurations and
measuring the quark condensate on each ensemble. We achieved a relative error
of about four percent. More importantly, the present method outperforms
previous alternatives because this determination of $\chi_{b}$ entails no
further systematic uncertainty, unlike approaches Bonati:2013lca ;
Bonati:2013vba ; Bali:2014kia , where a numerical differentiation of
$\log\mathcal{Z}(B)$ is required.

###  4.1 The magnetic susceptibility

Figure 3: Bare magnetic susceptibility at zero temperature versus the
logarithm of the lattice spacing, normalized to $a_{0}=1.46\textmd{
GeV}^{-1}$. Different approaches are compared: the finite difference method
Bonati:2013vba  (red triangles), the generalized integral method Bali:2014kia
(green circles) and the new approach via current-current correlators (blue
squares). The orange band indicates the fit based on perturbation theory, Eq.
(23).

We compare the results of our new method for $\chi_{b}$ to our old data
(generalized integral method) and also to those of Ref. Bonati:2013vba
(finite difference method) in Fig. 3 at zero temperature.101010 Ref.
Bonati:2013vba  employs the same lattice action. For a different action the
bare susceptibilities would not only differ in terms of lattice artifacts but
also by an additive constant, due to the different choice of renormalization
scheme.  Within errors perfect agreement between the three groups of results
is found. The data – plotted in Fig. 3 against $\log(a)$ – clearly reflect the
logarithmic divergence dictated by Eq. (5). Similarly to our fitting strategy
in Ref. Bali:2015msa , here we also include the universal perturbative QCD
corrections to the QED $\beta$-function coefficient $c_{1}$ Baikov:2012zm ,
see Eq. (4), where $g^{2}=6/\beta$ is obtained from the inverse lattice
coupling $\beta$ at the lattice scale $a^{-1}$. We also take into account
$\mathcal{O}(a^{2})$ lattice artifacts so that our fit function reads

$\chi_{b}=2\beta_{1}(a^{-1})\cdot\big{[}\log(a/a_{0})+\log(\mu_{\rm
QED}a_{0})\big{]}\cdot\left[1+z_{1}(a/a_{0})^{2}\right],\quad\quad
a_{0}=1.46\textmd{ GeV}^{-1}.$ (23)

The result of this fit, with the parameter values

$\mu_{\rm QED}=115(3)(5)\textmd{ MeV}\,,\quad\quad z_{1}=-0.05(1),$ (24)

is shown as an error band in Fig. 3. We also considered fits with further
(quartic) lattice artifacts. The impact of this is included in the second
parentheses of Eq. (24) for $\mu_{\rm QED}$ as a systematic error. The
renormalization scale agrees within errors with our earlier determinations
Bali:2014kia ; Bali:2015msa . It also lies near the mass of the lightest
charged hadron (the charged pion) which effectively sets the scale for the
magnetic response of this system. Nevertheless, note that the value of
$\mu_{\rm QED}$ depends on the choice of the regulator.

The formula (23) is used to interpolate $\chi_{b}$ and then employed to
renormalize the susceptibility at nonzero temperatures according to Eq. (6):
$\chi=\chi_{b}(T)-\chi_{b}(0)$. The results are shown in the left panel of
Fig. 4 for a broad range of temperatures and four lattice spacings
$a=1/(N_{t}T)$ with $N_{t}=6,8,10$ and $12$. A continuum extrapolation is
performed by means of a multi-spline fit Endrodi:2010ai  taking into account
$\mathcal{O}(a^{2})$ lattice artifacts. To have acceptable fits of this type,
we needed to discard the coarsest lattices ($N_{t}=6$ points at $T\lesssim
160\textmd{ MeV}$). We also repeated the analysis including
$\mathcal{O}(a^{4})$ discretization errors as well, this time fitting all
available data points. The systematic error was estimated by the difference of
these two extrapolations as well as by varying the spline node points and by
including/excluding $N_{t}=6$ data points at high temperatures for the
$\mathcal{O}(a^{2})$ fit. In addition, we consider a further systematic error
due to lattice artifacts related to the taste splitting of the staggered
spectrum. This effect is particularly relevant at low temperatures and can be
estimated by a generalization of the Hadron Resonance Gas (HRG) model that we
describe in App. A. The light yellow bands in both panels of Fig. 4 indicate
the total systematic uncertainties.

Figure 4: Left panel: renormalized susceptibility at nonzero temperature. The
symbols indicate different lattice spacings and the dark orange band the
continuum limit. The light orange band represents an estimate of systematic
errors of the continuum extrapolation. The dashed gray line is the HRG model
prediction Bali:2014kia . The inset zooms into the low-temperature region to
highlight the diamagnetic response there. Right panel: our results for $\chi$
(orange, labeled “via $\Pi(0)$”) are compared to the results of Ref.
Bonati:2013vba  (green) and those of Ref. Bali:2014kia  (red) as well as to
the HRG model prediction (dashed gray line) and to perturbation theory
Bali:2014kia  (dashed light blue band), see Eqs. (25) and (4).

The results demonstrate strong paramagnetism in the quark-gluon plasma phase,
in agreement with previous lattice studies Bonati:2013vba ; Bonati:2013lca ;
Bali:2013owa ; Bali:2013esa ; Levkova:2013qda ; Bali:2014kia . At high
temperatures the results are well described by the free theory, which predicts
(see App. C)

$\chi(T)=\beta_{1}(\mu_{\rm
therm})\cdot\log\left(\gamma\,\frac{T^{2}}{\mu_{\rm
QED}^{2}}\right)+\mathcal{O}(1/T^{2})\,,$ (25)

where $\gamma=\mathcal{O}(1)$ is a regulator-dependent constant.111111For
example in the free case with cut-off regularization
$\gamma=\pi^{2}e^{-\gamma_{E}}$, see App. C.5. The renormalized susceptibility
$\chi(T)$ and the ratio $\gamma/\mu_{\rm QED}^{2}$ are regulator-independent.
As indicated, QCD corrections at the thermal scale $\mu_{\rm therm}$ affect
the leading behavior. In the right panel of Fig. 4 we also include a
comparison to this perturbation theory formula, with the scheme-independent
$\mathcal{O}(\alpha_{s})$ corrections to $\beta_{1}$ taken into account
Bali:2014kia . Here we set $\gamma=1$ and use the $\overline{\mathrm{MS}}$
scheme definition of the strong coupling $\alpha_{s}=g^{2}/(4\pi)$, running
this at five-loop order Baikov:2016tgj  to the thermal scale $\mu_{\rm
therm}\sim 2\pi T$. We employ the central value $\Lambda_{\rm
QCD}^{\overline{\rm MS}}=0.341\textmd{ GeV}$ from the recent three flavor QCD
determination of $\alpha_{s}$ by the ALPHA Collaboration Bruno:2017gxd . The
band in the figure corresponds to a variation of the thermal scale from $\pi
T$ to $4\pi T$. The perturbative formula agrees surprisingly well with our
results down to temperatures $T\sim 200\textmd{ MeV}$.

In contrast to the paramagnetic behavior at high $T$, towards low temperatures
the continuum extrapolated results become negative, revealing a diamagnetic
response, previously noted in Ref. Bali:2014kia . This behavior is in
agreement, albeit within large errors, with the HRG model prediction
Bali:2014kia  (dashed line in the figures). In the right panel of Fig. 4 we
also include results obtained from other approaches that employed the same
lattice action. Around the pseudo-critical temperature $T_{c}\approx
155\textmd{ MeV}$ a significant difference is visible between earlier
determinations and our present results. Ref. Bonati:2013vba  carried out a
continuum extrapolation using fixed $\beta$ ensembles with lattice spacings
$a\geq 0.125\textmd{ fm}$. In Ref. Bali:2014kia  we used the fixed $N_{t}$
approach with $N_{t}=6,8,10$ ensembles, while in the present study
$N_{t}=6,8,10,12$ lattices are simulated. To highlight the differences between
the continuum extrapolations, in the left panel of Fig. 5 we plot the lattice
spacing-dependence of the susceptibility for all three approaches. We pick one
temperature $T=130\textmd{ MeV}$, where the deviation of the continuum
estimates is substantial.

Figure 5:  Left panel: lattice discretization errors in the renormalized
magnetic susceptibility. The results of Ref. Bonati:2013vba  (green) and those
of Ref. Bali:2014kia  (red) are compared to the present approach (yellow)
including systematic uncertainties (light yellow). The dashed gray and the
light blue bands represent the $T=130\textmd{ MeV}$ slices of the multi-spline
fit involving up to $\mathcal{O}(a^{2})$ and $\mathcal{O}(a^{4})$ lattice
artefacts, respectively. Note that for the green points at $a>0$, $\chi$ was
obtained by temperature-interpolations of the results published in Ref.
Bonati:2013vba . Right panel: parametrization of the relative magnetic
permeability $\mu/\mu_{0}=(1-e^{2}\chi)^{-1}$ via the function (88) of App. F.

The left panel of Fig. 5 reveals the importance of our new $N_{t}=12$
ensemble, showing a significant downward trend as the lattice spacing is
reduced and a negative value in the continuum limit. The downward trend is not
captured by our previous estimate using the integral method on $N_{t}\leq 10$
lattices Bali:2014kia , neither is it visible in the data of Ref.
Bonati:2013vba . We note that in the left panel of Fig. 5, a difference beyond
one standard deviation can only be observed at the smallest lattice spacing.
Nevertheless, our $N_{t}=12$ data lie consistently below the other lattice
spacings for all temperatures (see the left panel of Fig. 4) so that the
downward trend towards $a\to 0$ is statistically significant. We indeed expect
lattice artefacts in $\chi$ to be large and positive in this temperature
region, as predicted by the generalized HRG model of App. A. Finally we remark
that Ref. Bonati:2013vba  performed the continuum extrapolation assuming a
strictly positive function for $\chi(T)$. Excluding the possibility of a
negative susceptibility in the continuum limit might in general underestimate
the systematics of the extrapolation. To clarify this issue, dedicated
simulations should be performed with the same action using all available
methods, preferably at the same temperatures and the same values of the
lattice spacing.

Finally we provide a parametrization that connects all three approaches (HRG,
lattice continuum limit and perturbation theory) and describes $\chi$ for
arbitrary temperatures. The details are discussed in App. F. In the right
panel of Fig. 5 we plot this parametrization, translated to the magnetic
permeability $\mu/\mu_{0}=(1-e^{2}\chi)^{-1}$, expressed in units of the
vacuum permeability $\mu_{0}$. This combination is equal to the ratio of the
magnetic induction and the external field, see, e.g., Refs. Bonati:2013lca ;
Bali:2014kia .

###  4.2 The normalization of the photon distribution amplitude

Here we address the tensor coefficients $\tau_{fb}$ at zero temperature. We
consider a set of independent gauge ensembles, generated at the physical value
of the strange quark mass $m_{s}=m_{s}^{\rm phys}$, but at different values of
the light quark mass: $0.5\,m_{ud}^{\rm phys}\leq m_{ud}\leq m_{s}^{\rm
phys}$. We follow a similar strategy as in Ref. Bali:2012jv , simultaneously
fitting the dependence of $\tau_{ub}\cdot Z_{T}$ on the light quark mass
$m_{ud}$ and on the lattice spacing $a$ according to the ansatz (9). Since
$Z_{T}$ is found to depend very mildly on the lattice spacing within the range
covered (see Fig. 11 of App. D), this does not significantly affect the
functional dependence on $a$. We note that on our coarsest ensembles the
uncertainty of $Z_{T}$ is quite large, which also imprints on the errors of
the renormalized tensor coefficients.

Figure 6: Light quark mass-dependence of the tensor coefficient at $T=0$ for
the up quark using our four finest lattice spacings (green to gray symbols).
The index $b$ indicates that the QED divergence that one encounters at
$m_{ud}>0$ has not been subtracted. The results diverge logarithmically
towards the continuum limit for any $m_{ud}\neq 0$. In contrast, the chiral
limit is free of ultraviolet divergences and a combined chiral and continuum
limit exists (black circle).

Notice that $\tau_{fb}$ diverges for $a\to 0$ for any quark mass, except in
the chiral limit, where it is ultraviolet-finite. This tendency is clearly
visible in Fig. 6, which shows our results for the up 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 $\chi^{\rm spin}$ we will need to take
differences between results obtained at different temperatures (see below).

The ansatz (9) contains the free parameters $\tau_{f}$ and $\mu_{\rm QED}$. In
addition, we include a quadratic mass-dependence and lattice artifacts of
$\mathcal{O}(a^{2})$ to each parameter in the fit. Varying the fit ranges in
$a$ and in $m_{ud}/m_{ud}^{\rm phys}$ as well as the functional form, we carry
out several acceptable fits that are used to build a histogram for the chiral
continuum limit of the tensor coefficient. In this combined limit we obtain in
the $\overline{\mathrm{MS}}$ scheme

$T=0:\quad\quad f_{\gamma}^{\perp}(2\textmd{ GeV})\equiv\lim_{m_{ud}\to
0}\tau_{ub}\cdot Z_{T}(2\textmd{ GeV})=-45.4(1.5)\textmd{ MeV}\,.$ (26)

The central value differs from our previous result Bali:2012jv
$f_{\gamma}^{\perp}=-40.3(1.4)\textmd{ MeV}$, mainly due to the multiplicative
renormalization factor that we determined non-perturbatively here, see Fig. 11
in App. D.

Figure 7: The bare tensor coefficient for the up quark (left panel) and for
the strange quark (right panel) at the physical point and at zero temperature
(blue points), together with an interpolation (orange bands). For the up quark
this interpolation is the $m_{ud}=m_{ud}^{\rm phys}$ slice of a two-
dimensional fit like in Fig. 6. The red dashed lines indicate the leading
logarithmic divergence $\propto m_{f}\log a$ in both fits. The remaining
$a$-dependence is consistent with lattice artifacts. We always use
$m_{s}=m_{s}^{\rm phys}$. The lattice spacing is normalized to
$a_{0}=1.46\textmd{ GeV}^{-1}$.

In the left panel of Fig. 7 we show the $a$-dependence of the $T=0$ light
quark tensor coefficient $\tau_{ub}\cdot Z_{T}$ at the physical point and the
result of the above interpolation, including the systematic error estimated
using the different fits. For demonstration purposes, we also indicate the
leading logarithmic behavior, that we obtain by subtracting the lattice
artifact terms from the central fit. Comparing to the similar plot for
$\chi_{b}$ (Fig. 3), we see that deviations from the continuum behavior are
sizable (and are predominantly due to the fact that we are dealing with a
dimensionful quantity in this case). For this reason, here we cannot reliably
determine the value of $\mu_{\rm QED}$. Nevertheless, we note that fixing the
renormalization scale to its value from Eq. (24) also gives acceptable fits.
This is in agreement with the expectation of Sec. 2.3, as well as with the
results in the free case, see App. C.4. The logarithmic divergence $\propto
m_{f}\log a$ becomes more pronounced for heavy quarks. This is visible in the
right panel of Fig. 7, where we plot the strange quark tensor coefficient
$\tau_{sb}\cdot Z_{T}$ at $m_{ud}=m_{ud}^{\rm phys}$ against the lattice
spacing and again indicate the leading logarithmic term.

As we have discussed above, at non-vanishing values of the quark mass $m_{f}$,
the tensor coefficient diverges logarithmically. In Ref. Bali:2012jv  we
suggested to cancel this by taking the logarithmic derivative with respect to
the quark mass, see also Ref. Ioffe:1983ju :

$f_{\gamma f}^{\perp}=\left(1-m_{f}\frac{\partial}{\partial
m_{f}}\right)\tau_{f}\cdot Z_{T}\,.$ (27)

This renormalization prescription will give identical results for any
regulator, up to the multiplicative factor $Z_{T}$.121212 This construction
will not only cancel the logarithmic divergence but also any finite term
$\propto m_{f}$. Should this be unwanted then one will have to accept a
scheme-dependence and convert between different schemes in a similar way as is
done for the massive chiral condensate, e.g., in Ref. McNeile:2012xh . Using
this prescription, it turns out that $f_{\gamma u}^{\perp}=f_{\gamma
d}^{\perp}=f_{\gamma}^{\perp}$ holds within statistical errors. For the
strange quark we obtain:

$f_{\gamma s}^{\perp}=-68(3)(4)\textmd{ MeV}\,.$ (28)

The first error includes the described variation of the fit while the second
error reflects the uncertainty of the derivative with respect to $m_{s}$ that
we indirectly determine from the dependence of the tensor coefficient on the
light quark mass, following the procedure explained in Ref. Bali:2012jv .

In the literature often the magnetic susceptibility of the quark condensate,

$X_{u}=\frac{\tau_{u}}{\left\langle\bar{\psi}_{u}\psi_{u}\right\rangle}\cdot\frac{Z_{T}}{Z_{S}}\,,$
(29)

is given, rather than $f^{\perp}_{\gamma}=\tau_{u}\cdot Z_{T}$. Since the
latter quantity has a smaller anomalous dimension and its value does not
depend on a separate computation of the chiral condensate, this is the
preferred choice for practical applications. However, for convenience of
comparison, we shall convert it into the other convention. The numerical value
of the quark condensate in the SU(2) chiral limit in the
$\overline{\mathrm{MS}}$ scheme at the scale $\mu_{\rm QCD}=2\textmd{ GeV}$
reads $\langle\bar{\psi}_{u}\psi_{u}\rangle=[272(5)\textmd{ MeV}]^{3}$
Aoki:2019cca . To enable a comparison with other results, below we also list
$X_{u}$ at the scale $\mu_{\rm QCD}=2\textmd{ GeV}$. Since most literature
values refer to a low, sometimes unspecified scale, in addition we run $X_{u}$
as well as $f_{\gamma}^{\perp}$ to the scale $\mu_{\rm QCD}=1\textmd{ GeV}$,
which is used in most sum rule calculations, see, e.g., Refs. Balitsky:1985aq
; Balitsky:1989ry ; Ball:2002ps . This is done, using the five-loop $\beta$\-
and quark mass anomalous dimension $\gamma$-functions Baikov:2016tgj ;
Baikov:2014qja  and the three-loop $\gamma$-function of the tensor current
Broadhurst:1994se ; Gracey:2003yr . The results read

$\displaystyle X_{u}(2\textmd{ GeV})$ $\displaystyle=-\left[665(13)\textmd{
MeV}\right]^{-2}\,,\quad X_{u}(1\textmd{ GeV})=-\left[542(11)\textmd{
MeV}\right]^{-2}\,,$ (30) $\displaystyle f_{\gamma}^{\perp}(1\textmd{ GeV})$
$\displaystyle=-51.1(1.6)\textmd{ MeV}\,,$ (31)

where we have added all errors in quadrature, including the uncertainty of
$f^{\perp}_{\gamma}(2\textmd{ GeV})$, the difference between running with the
two- and three-loop $\gamma$-functions of the tensor current, the uncertainty
of $\langle\bar{\psi}_{u}\psi_{u}\rangle$ and the uncertainty of the strong
coupling parameter Bruno:2017gxd . All the above results are in the
$\overline{\mathrm{MS}}$-scheme.

We summarize earlier results from the literature for comparison. The first sum
rule determination of $X_{u}$ Ioffe:1983ju  suggested a value
$X_{u}(0.5\textmd{ GeV})=-[350(50)\textmd{ MeV}]^{-2}$ while vector meson
dominance yields Balitsky:1983xk  $X_{u}\approx 2/m_{\rho}\approx-(540\textmd{
MeV})^{-2}$. This was improved upon in subsequent sum rule determinations,
see, e.g., Ref. Ball:2002ps  and references therein. The most extensive sum
rule study Ball:2002ps  found $X_{u}(1\textmd{ GeV})\approx-(560\textmd{
MeV})^{-2}$, which agrees reasonably well with our determination. A
comparatively smaller absolute value $f_{\gamma}^{\perp}\approx-38\textmd{
MeV}$ was obtained at a low scale $\mu\sim 600\,\textmd{ MeV}$ in the quark-
soliton model Petrov:1998kg  while the Vainshtein relation Vainshtein:2002nv
suggests an even smaller modulus of the magnetic susceptibility of the quark
condensate $X_{u}=-N_{c}/(4\pi^{2}F_{\pi}^{2})\approx-(335\textmd{
MeV})^{-2}$. This parameter was also considered in holographic studies, with
the result $X_{u}\approx-(295\textmd{ MeV})^{-2}$ Gorsky:2009ma , while NJL-
and quark-meson-model predictions give $X_{u}\approx-(480\textmd{ MeV})^{-2}$
Frasca:2011zn  and $X_{u}\approx-(440\textmd{ MeV})^{-2}$ Frasca:2011zn ,
respectively. Finally, quenched lattice simulations, without renormalization,
gave the values $X_{u}\approx-[804(3)\textmd{ MeV}]^{-2}$ in SU(2) gauge
theory Buividovich:2009ih  and $X_{u}\approx-[486(21)\textmd{ MeV}]^{-2}$ in
SU(3) Braguta:2010ej . Our previous full QCD study Bali:2012jv  resulted in
$X_{u}(2\textmd{ GeV})=-[693(13)\textmd{ MeV}]^{-2}$, however, in that case
the renormalization was only carried out perturbatively.

Our result (28) for the strange quark coefficient translates into

$\displaystyle X_{s}(2\textmd{ GeV})$ $\displaystyle=-\left[565(50)\textmd{
MeV}\right]^{-2}\,,\quad X_{s}(1\textmd{ GeV})=-\left[460(41)\textmd{
MeV}\right]^{-2}\,,$ (32) $\displaystyle f_{\gamma s}^{\perp}(1\textmd{ GeV})$
$\displaystyle=-76.5(5.7)\textmd{ MeV}\,,$ (33)

where we used the ratio
$\langle\bar{\psi}_{s}\psi_{s}\rangle/\langle\bar{\psi}_{u}\psi_{u}\rangle=1.08(17)$
McNeile:2012xh  for the conversion between $f^{\perp}_{\gamma s}$ and $X_{s}$.
The difference between $X_{s}$ and $X_{u}(1\textmd{
GeV})\approx-\left(475\textmd{ MeV}\right)^{-2}$ has been reported to be
negligible in the sum rule calculations Balitsky:1989ry .

###  4.3 The spin contribution at $T>0$

Having interpolated the $T=0$ tensor coefficients, we are now in the position
to perform the additive renormalization (10) by subtracting this contribution
from the finite temperature results. We use our existing $N_{t}=6$, $8$ and
$10$ results from Ref. Bali:2012jv  to approach the continuum limit.
Especially in light of the slow convergence of $\chi$ towards $a\to 0$, see
the right panel of Fig. 5, this extrapolation should be backed up with finer
lattice ensembles in the future. In analogy to the analysis of $\chi$, again
we carry out a multi-spline fit of all data sets, determining a systematic
error by varying the positions of the spline node points. The so-obtained fit
is shown for the up quark and the strange quark in Fig. 8. The results for
$\tau_{d}$ are consistent with $\tau_{u}$ within errors. The large errors of
our $N_{t}=6$ results at low temperatures are due to the uncertainties of the
$T=0$ contributions on our coarse lattices, see Fig. 7.

Figure 8: Tensor coefficients after multiplicative as well as additive
renormalization for the up (left panel) and for the strange quark (right
panel).

After the additive renormalization, the tensor coefficient vanishes by
definition at $T=0$. For the light quarks $\tau_{u}(T)$ grows substantially as
the temperature is increased, before the slope reduces and a plateau is
approached. The inflection point of the continuum curve is found to be at
$T_{c}=158(5)\textmd{ MeV}$. The chiral transition temperature determined from
the inflection point of the quark condensate $T_{c}=155(4)\textmd{ MeV}$
Borsanyi:2010bp  is in agreement with this value. For the strange quark
pseudo-critical thermal effects set in at somewhat higher temperatures
Borsanyi:2010bp . Also in our case $\tau_{s}$ does not appear to exhibit any
inflection point, at least for $T\lesssim 170\textmd{ MeV}$, and below
$T\approx 200\textmd{ MeV}$ no saturation into a plateau is visible. For
sufficiently high temperatures, where the finite quark mass becomes
negligible, we expect the two renormalized tensor coefficients to coincide.

Next, the continuum extrapolated results are inserted into Eq. (12) to
determine the spin contribution $\chi^{\rm spin}$ to the susceptibility. To
this end we need to evaluate the tensor bilinear for massless valence quarks.
Instead of performing measurements at additional valence quark masses, we
estimate this limit using the difference between the results for the strange
quark and for the light quarks. We assume a linear dependence on the valence
quark mass in the range $[0,m_{s}]$, which implies that

$\lim_{m^{\rm val}_{u}\to 0}\tau_{u}=\lim_{m_{s}^{\rm val}\to
0}\tau_{s}\approx\frac{\tau_{u}m_{s}-\tau_{s}m_{ud}}{m_{s}-m_{ud}}=\tau_{u}\frac{R}{R-1}-\tau_{s}\frac{1}{R-1}\,,\quad\quad
R\equiv\frac{m_{s}}{m_{ud}}\,.$ (34)

In this case the contributions of all flavors to $\chi^{\rm spin}$ are
proportional to $\tau_{s}-\tau_{u}$ and the renormalized spin susceptibility
(12) simplifies to

$\chi^{\rm spin}\approx\frac{1}{2m_{ud}}\frac{\tau_{s}-\tau_{u}}{R-1}\cdot
Z_{T}Z_{S}\cdot\sum_{f}(q_{f}/e)^{2}\,.$ (35)

Thus, in this approximation the individual flavors simply contribute in
proportion to their squared electric charges. The scalar renormalization
constants entering this expression are displayed in Fig. 11 of App. D.

The so-obtained estimate of $\chi^{\rm spin}$ is shown in the left panel of
Fig. 9 for three lattice spacings, together with a continuum extrapolation
performed in the same way as for $\tau_{f}$. We observe $\chi^{\rm spin}<0$
for all temperatures, with a minimum somewhat above the pseudo-critical
temperature and an upward trend for high temperatures. The approximation (34)
tends to overestimate the valence chiral limit of the tensor coefficient due
to the presence of logarithmic deviations from a linear behavior in
$m_{f}^{\rm val}$.131313This is also visible in Fig. 6, although the
dependencies on the valence and sea quark masses are not disentangled in that
figure. Consequently, Eq. (35) underestimates $\chi^{\rm spin}$. This is also
the case at high temperatures, as can be checked using the analytic formula
valid for the free case, see App. C. To take this effect into account we
include a systematic error based on the free case formula (70). In particular,
we consider the difference between the approximation and the true value in the
free case and scale it with the typical magnitude of the light quark tensor
coefficient at lower temperatures (see Fig. 8). The so-obtained uncertainty is
also included in the left panel of Fig. 9.

Figure 9: Left panel: spin contribution to the susceptibility using three
lattice spacings (colored symbols) and an extrapolation to the continuum limit
(orange band). A systematic uncertainty, related to the estimation of the
tensor coefficient for massless valence quarks, is indicated by the light
yellow band. Right panel: the total magnetic susceptibility from Fig. 4
(blue), together with the decomposition into spin (orange-yellow) and orbital
angular momentum (green-gray) contributions.

We remark that $\chi^{\rm spin}<0$ for the temperature range covered in our
simulations. This can be understood by noting that Eq. (35) is the
discretization of the mass-derivative of $\tau_{f}$. Increasing the mass
pushes the inflection point of $\tau_{f}$ to higher temperatures (visible in
Fig. 8), thus making the derivative negative around the transition
temperature. Nevertheless, $\chi^{\rm spin}$ will necessarily turn positive
for even higher temperatures. Indeed, for sufficiently high temperatures the
difference $\tau_{f}=\tau_{fb}(T)-\tau_{fb}(T=0)$ will be dominated by the
$T=0$ term, so that Eq. (35) becomes proportional to
$\tau_{ub}(T=0)-\tau_{sb}(T=0)$, which is positive for any lattice spacing
(see Fig. 6). Perturbation theory also predicts $\chi^{\rm spin}>0$ for high
temperatures, see App. C.3.

###  4.4 Pauli and Landau decomposition of the magnetic susceptibility

Finally, we compare the spin contribution to the total susceptibility in order
to learn about the orbital angular momentum-related contribution $\chi^{\rm
ang}=\chi-\chi^{\rm spin}$. All three susceptibilities are included in the
right panel of Fig. 9. While the errors of the two contributions are much
larger than that of the total susceptibility, several qualitative comments can
be made based on this plot. First of all, in the complete temperature range
under study, $\chi^{\rm spin}$ and $\chi^{\rm ang}$ have opposite signs and
$\chi$ emerges as a result of a large cancellation between the two terms. As
we argued above, the spin part will necessarily turn positive for higher
temperatures, eventually approaching $3/2$ times the full susceptibility.
Consequently, $\chi^{\rm ang}$ will turn negative and approach
$-1/2\cdot\chi$. It is intriguing to observe that in the strongly interacting
regime the two contributions have opposite signs than in the usual free
fermion picture according to Pauli and Landau: it is the Landau term that
drives the paramagnetic response of the QCD vacuum up to temperatures
$T\gtrsim 200\textmd{ MeV}$, while the Pauli term reduces the susceptibility
in this region. This unusual behavior becomes possible due to the strong
interaction, which confines quarks into composite hadrons and thereby fixes
the relative orientation of their spins, i.e. their magnetic moments. In
particular, in charged pions one of the constituent quarks is bound to anti-
align its magnetic moment with the background field in order to maintain zero
total spin. Similar effects arise for certain baryons as well. Beyond this
qualitative argument, it is difficult to anticipate the outcome of this
competition between the strong and the electromagnetic forces. Our
quantitative results reveal a peculiar interplay between confinement and spin
physics.

To further our understanding, in principle $\chi$ can also be decomposed into
$\chi_{f}$ for the quark flavors $f$ and a gluonic contribution $\chi_{g}$.
Subtracting this $\chi_{g}$ from $\chi^{\rm ang}$ will isolate the total quark
orbital angular momentum contribution $\sum_{f}(\chi_{f}-\chi_{f}^{\rm
spin})$, in analogy to spin decompositions Ji:1996ek  in deep inelastic
scattering that are based on the Belinfante-Rosenfeld energy-momentum tensor,
in this case of the transverse spin. The unrenormalized qualitative results of
Ref. Bali:2013esa  indicate that $\chi_{g}\sim\chi/3$ at small temperatures.
It may be interesting to address this quantitatively in the future.

##  5 Summary

In this paper we determined the magnetic susceptibility $\chi$ of the thermal
QCD medium via a method introduced originally for $T=0$ Bali:2015msa , which
circumvents the flux quantization problem and allows us to express $\chi$ in
terms of $B=0$ measurements. This considerably reduces the measurement costs
as well as systematic uncertainties compared to previous approaches. The
susceptibility is extrapolated to the continuum limit for a broad range of
temperatures, making contact to the Hadron Resonance Gas (HRG) model at low
$T$ as well as to perturbation theory at high $T$. In the confined phase we
find evidence for a diamagnetic behavior ($\chi<0$), while for $T\gtrsim
150\textmd{ MeV}$ we observe paramagnetism ($\chi>0$). Our continuum
extrapolations are based on four lattice spacings and are guided by a
generalized HRG model taking into account taste splitting (see App. A). A
careful continuum limit is found to be essential to observe diamagnetism at
low $T$ since this is due to light pions – we argue that this behavior was
missed in previous investigations because of large lattice artifacts.

The susceptibility is decomposed into spin- ($\chi^{\rm spin}$) and orbital
angular momentum-related ($\chi^{\rm ang}$) contributions based on our
previous study Bali:2012jv . The spin term is shown to be given in terms of
the mass-dependence of the
$\left\langle\bar{\psi}\sigma_{12}\psi\right\rangle$ fermion bilinear in the
presence of a small magnetic field, see Eq. (12) and App. B. Besides its role
in this decomposition, the tensor bilinear is related to the normalization
$f_{\gamma}^{\perp}$ of the photon distribution amplitude, relevant for a
range of phenomenological applications. We update our previous determination
Bali:2012jv  of the corresponding tensor coefficient in the chiral limit at
$T=0$, by performing the multiplicative renormalization of
$\left\langle\bar{\psi}\sigma_{12}\psi\right\rangle$ non-perturbatively on the
lattice. We obtain the value $f_{\gamma}^{\perp}=-45.4(1.5)\textmd{ MeV}$ for
massless quarks, in the $\overline{\rm MS}$ scheme at a QCD renormalization
scale of $2\textmd{ GeV}$. The values of the tensor coefficient at the
physical light and strange quark masses and at different renormalization
scales are given in Eqs. (28)–(33).

At finite temperatures we performed the continuum extrapolation of $\chi^{\rm
spin}$ and also determined the orbital angular momentum-related susceptibility
$\chi^{\rm ang}$. In the absence of color interactions, the two contributions
exhibit the constant ratio $\chi^{\rm spin}:\chi^{\rm ang}=3:(-1)$ as is well
known since the analysis of the free electron gas by Pauli Pauli:1927:GPG  and
Landau Landau1930 . Around the transition temperature, in full QCD this ratio
is instead found to be close to $(-1):(1.03)$, resulting in a large
cancellation between the two contributions, thereby substantially reducing the
total susceptibility. As the temperature grows the susceptibilities approach
their free-case counterparts, which are discussed in detail in App. C. Still,
it is stunning to observe that in the strongly coupled QCD medium $\chi^{\rm
spin}$ and $\chi^{\rm ang}$ have signs that are opposite to the naive
expectations.

Considering our results at high temperature, it is interesting to make a
comparison to a classical ideal system. In such a setting the Bohr-van Leeuwen
theorem Bohr ; BVL  (for a recent review, see Ref. BVL2 ) holds: the total
magnetization vanishes, since the magnetic field does not transfer any work to
the electric currents in the system. Apparently, the QCD medium does not
become classical in this sense for $T\to\infty$, even if the
$\mathcal{O}(B^{2})$ terms of the free energy density that we have discussed
in this paper are small compared to the dominant $\mathcal{O}(T^{4})$
contributions in that limit. The non-classicality has two different origins.
First, quark spins are of quantum nature and can induce a magnetization by
aligning with the magnetic field. Second, both $\chi^{\rm spin}$ and
$\chi^{\rm ang}$ diverge as $\log T$ for high temperatures. This behavior
stems from the renormalization properties of the bare susceptibilities:
quantum effects give rise to a logarithmic divergence $\propto\log 1/a$ in the
cut-off. In turn, the same behavior shows up in the renormalized
susceptibilities if they are probed by another large dimensionful scale, the
temperature. Note that a similar connection exists between the logarithmic
divergence and the behavior of the renormalized free energy in the
$B\to\infty$ limit Dunne:2004nc .

###### Acknowledgements.

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, Nikolay Kivel and Massimo D’Elia for insightful
comments as well as Ren-Hong Fang for pointing out a mistake in an earlier
version of App. C.5.

##  Appendix A The HRG model and lattice discretization errors

At low temperatures the staggered action suffers from enhanced lattice
artifacts due to taste splitting. Here we attempt to incorporate the effects
of this splitting into the HRG model. The magnetic susceptibility was
calculated in a standard HRG model in Ref. Bali:2014kia . Following Ref.
Huovinen:2009yb  we replace the contribution of pions in the model by a sum
over each taste, weighted by the corresponding degeneracies. The masses of the
individual tastes and their parametrization in the range of our lattice
spacings are taken from Ref. Borsanyi:2010cj . Since pions are dominant for
the susceptibility, the taste splitting for other mesonic and baryonic states
is ignored for simplicity (although the splitting for $\eta$ mesons might also
lead to light mesonic states, see, e.g., Ref. Bazavov:2009bb ). The list of
hadrons taken into account can be found in Ref. Endrodi:2013cs .

Figure 10:  Lattice artifacts in the susceptibility in a generalized HRG model
involving taste splitting.

In Fig. 10 we show the renormalized magnetic susceptibility evaluated at
$T=120\textmd{ MeV}$ as a function of the lattice spacing. The spacings for
our four ensembles $N_{t}=6,8,10$ and $12$ at this temperature are highlighted
in the plot. This reveals slow convergence towards the continuum limit, which
can best be understood by analyzing the mass-dependence of the pionic
contribution $\chi_{\pi}$ to the susceptibility, which takes the form
Bali:2014kia

$\chi_{\pi}(m_{\pi})=-\frac{1}{48\pi^{2}}\int_{0}^{\infty}\frac{\textmd{d}t}{t}e^{-m_{\pi}^{2}t/T^{2}}\left[\Theta_{3}\left(0,e^{-1/(4t)}\right)-1\right]\,,$
(36)

where $\Theta_{3}$ is an elliptic $\Theta$-function. This can be derived by
comparing to the analogous expression for fermions, calculated below in Eq.
(53). The bosonic Matsubara frequencies give rise to the different first
argument in the elliptic function. The prefactor in this case is the scalar
QED $\beta$-function coefficient for one complex scalar field $\beta_{1}^{\rm
scalar}=1/(48\pi^{2})$. The pionic susceptibility diverges logarithmically in
the chiral limit (this can be shown similarly to the calculation below in App.
C.5),

$\chi_{\pi}(m_{\pi})\xrightarrow{m_{\pi}\to 0}-\beta_{1}^{\rm
scalar}\log(T/m_{\pi})^{2}\,,$ (37)

explaining its pronounced dependence on $m_{\pi}$. In turn, nonzero lattice
spacings enhance the masses of most pion tastes, thus, reducing the magnitude
of $\chi_{\pi}$.

Based on the HRG predictions for $\chi(a,T)$ we consider the difference
between a simple $\mathcal{O}(a^{2})$ fit taking into account only $N_{t}\leq
12$ lattices and the true continuum limit. This difference is included as a
lower systematic error of our lattice determination of $\chi(T)$ at low
temperatures, see Fig. 4.

##  Appendix 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. (8)–(12) of the main
text. It is instructive to begin with the first derivative of the free energy
density,

$-\frac{\partial f}{\partial
B}=\frac{T}{V}\sum_{f}\left\langle\textmd{tr}\,\frac{1}{\not{D}_{f}+m_{f}}\frac{\partial\not{D}_{f}}{\partial
B}\right\rangle=\frac{T}{2V}\sum_{f}\left\langle\textmd{tr}\,\frac{1}{(\not{D}_{f}+m_{f})\not{D}_{f}}\frac{\partial\not{D}_{f}^{2}}{\partial
B}\right\rangle\,,$ (38)

where we used the cyclicity of the trace (even though $\not{D}_{f}$ and
$\partial\not{D}_{f}/\partial B$ do not commute, we can symmetrize the
expression in the two operators under the trace). Now we use the relation

$\frac{1}{(\not{D}_{f}+m_{f})\not{D}_{f}}=-\frac{1}{m_{f}}\left[\frac{1}{\not{D}_{f}+m_{f}}-\frac{1}{\not{D}_{f}}\right]\,,$
(39)

and the identities

$\frac{\partial\not{D}_{f}^{2}}{\partial(q_{f}B)}=-\sigma_{12}-L_{12},\quad\quad\sigma_{12}=\frac{1}{2i}[\gamma_{1},\gamma_{2}],\quad\quad
L_{12}=-\frac{\partial D_{f}^{2}}{\partial(q_{f}B)}\,,$ (40)

where $\sigma_{12}$ is the relevant component of the relativistic spin
operator defined in Eq. (7) and $L_{12}$ is a generalized angular momentum
operator, which depends on the electromagnetic as well as the $\mathrm{SU}(3)$
gauge.

Using Eqs. (39) and (40), we can rewrite Eq. (38) as

$-\frac{\partial f}{\partial
B}=\frac{T}{2V}\sum_{f}\frac{q_{f}}{m_{f}}\left\langle\textmd{tr}\,\frac{\sigma_{12}+L_{12}}{\not{D}_{f}+m_{f}}-\textmd{tr}\,\frac{\sigma_{12}+L_{12}}{\not{D}_{f}}\right\rangle=\sum_{f}\frac{q_{f}}{2m_{f}}\bigg{[}1-\lim_{m_{f}^{\rm
val}\to
0}\bigg{]}\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}+\bar{\psi}_{f}L_{12}\psi_{f}\right\rangle\,.$
(41)

Thus, in the language of Eq. (13), we need the difference of two terms: one
with valence quark mass $m_{f}^{\rm val}=m_{f}$ and one with $m_{f}^{\rm
val}\to 0$. The sea quark mass is kept fixed in both cases: $m_{f}^{\rm
sea}=m_{f}$. We remark that the vanishing valence quark mass needs to be
defined as a limit in finite volumes (see below). Also note that the fermion
bilinears are defined to include the volume factor $T/V$.

Differentiating Eq. (41) once more with respect to $B$ at $B=0$ and dividing
by $e^{2}$, we recover the bare magnetic susceptibility (1) on the left hand
side,

$\chi_{b}=\sum_{f}\frac{(q_{f}/e)^{2}}{2m_{f}}\bigg{[}1-\lim_{m_{f}^{\rm
val}\to 0}\bigg{]}\lim_{B\to
0}\frac{\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}+\bar{\psi}_{f}L_{12}\psi_{f}\right\rangle}{q_{f}B}\,.$
(42)

The slope of the tensor bilinear
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$ for small values
of $B$ gives the tensor coefficient $\tau_{fb}$ as defined in Eq. (8). After
subtracting its value at $T=0$ and multiplying by the relevant QCD
renormalization factors, this term gives the spin contribution to the
susceptibility $\chi^{\rm spin}$, as we wrote in the main text, Eq. (12). In
turn, the magnetic field-dependence of the bilinear involving the generalized
angular momentum operator $L_{12}$ is related to $\chi^{\rm ang}$. The latter
term cannot be implemented straightforwardly due to its gauge-dependence and
magnetic flux quantization.

In Ref. Bali:2012jv  we already discussed the separation of the magnetic
susceptibility into quark spin- and other angular momentum-related
contributions. There, the $m_{f}^{\rm val}=0$ term of Eq. (41) was argued not
to contribute – indeed, in a finite volume the massless limit of fermion
bilinears always vanishes. However, in the thermodynamic limit this is not the
case if chiral symmetry is broken spontaneously. To elucidate this point in
more detail, let us rewrite the trace in Eq. (41) using the eigenmodes of the
Dirac operator,

$\not{D}_{f}\chi_{f\lambda}=i\lambda\chi_{f\lambda}\,,$ (43)

so that, exploiting chiral symmetry $\\{\gamma_{5},\not{D}_{f}\\}=0$,

$\displaystyle\frac{T}{V}\left\langle\\!\textmd{tr}\,\frac{\sigma_{12}}{\not{D}_{f}+m_{f}^{\rm
val}}\\!\right\rangle$ $\displaystyle=\frac{T\,m_{f}^{\rm
val}}{V}\left\langle\\!\textmd{tr}\,\frac{\sigma_{12}}{-\not{D}_{f}^{2}+(m_{f}^{\rm
val})^{2}}\\!\\!\right\rangle$
$\displaystyle\xrightarrow{V\to\infty}\int_{0}^{\infty}\\!\\!\\!\textmd{d}\lambda\,\frac{2m_{f}^{\rm
val}}{\lambda^{2}+(m_{f}^{\rm
val})^{2}}\left\langle\\!\rho_{f}(\lambda;m_{f}^{\rm
sea})\,\chi_{f\lambda}^{\dagger}\sigma_{12}\chi_{f\lambda}\right\rangle\,,$
(44)

where $\rho_{f}(\lambda;m_{f}^{\rm sea})$ is the spectral density of
$\not{D}_{f}$ in the infinite volume, determined in an ensemble generated with
sea quark masses $m_{f}^{\rm sea}$. Towards the valence chiral limit the
kernel becomes proportional to the $\delta$-distribution, so that we have a
Banks-Casher-type Banks:1979yr  relation,

$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle\xrightarrow{V\to\infty,\,m_{f}^{\rm
val}\to
0}\pi\int_{0}^{\infty}\\!\\!\textmd{d}\lambda\,\delta(\lambda)\left\langle\rho_{f}(\lambda;m_{f}^{\rm
sea})\,\chi_{f\lambda}^{\dagger}\sigma_{12}\chi_{f\lambda}\right\rangle=\pi\left\langle\rho_{f}(0;m_{f}^{\rm
sea})\,\chi^{\dagger}_{f0}\sigma_{12}\chi_{f0}\right\rangle\,.$ (45)

On the one hand, this limit is zero if chiral symmetry is intact and the
spectral density vanishes at the origin. On the other hand, a nonzero chiral
condensate $\left\langle\rho_{f}(0;m_{f}^{\rm sea})\right\rangle$, together
with the polarization $\sigma_{12}\chi_{f0}=\chi_{f0}$ of the low modes will
turn the chiral limit of the tensor bilinear nonzero. Our lattice results
reveal a nonzero value for
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$ in the full
chiral limit at low temperatures, see Fig. 6. Clearly, the fermion bilinear
remains nonzero also if only $m_{f}^{\rm val}$ is sent to zero. This is in
accordance with the recent findings of Ref. Bruckmann:2017pft  about the Dirac
spectrum at $B>0$, where the low modes were indeed found to exhibit almost
perfect spin-polarization.

##  Appendix C Susceptibilities in the free case

Here we consider the free case (i.e. we set the color charges of quarks to
zero) to exemplify the most important relations of the main text. These
include the proportionality between the tensor bilinear and the spin
contribution to the susceptibility, the ultraviolet divergences of the
susceptibilities at zero temperature as well as the high-temperature behavior
of the renormalized susceptibilities. These calculations include our previous
results Bali:2012jv ; Bali:2014kia , which we also show here for completeness.

Below we will extensively use Schwinger’s proper time formulation
Schwinger:1951nm . This is based on the Mellin transform

$E^{-z}=\frac{1}{\Gamma(z/2)}\int_{0}^{\infty}\textmd{d}t\,t^{z/2-1}\,e^{-E^{2}t}\,,$
(46)

and its inverse

$e^{-lE/T}=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\textmd{d}z\,\Gamma(z)\,l^{-z}E^{-z}\,T^{z}\,,$
(47)

which are valid for $\textmd{Re}\,z>0$, $c>0$ and $E>0$. Moreover, taking the
derivative of Eq. (46) with respect to $z$ at $z=0$ gives the standard
$\zeta$-function regularization result Elizalde:1994gf ,

$\log E^{2}=-2\left.\frac{\partial\,(E^{2})^{-z/2}}{\partial
z}\right|_{z=0}=-2\left.\frac{\partial}{\partial
z}\right|_{z=0}\frac{1}{\Gamma(z/2)}\int_{0}^{\infty}\textmd{d}t\,t^{z/2-1}\,e^{-E^{2}t}\,.$
(48)

###  C.1 Magnetic susceptibility

We consider one quark flavor $\psi$ with electric charge $q$ and mass $m$ in a
volume $V=L^{3}$ at temperature $T$, exposed to a background magnetic field
$B$. For convenience we assume that the magnetic field is oriented in the
$x_{3}$ direction and $qB>0$. The free energy density in this setting reads
(see, e.g., Ref. Fraga:2012rr ):

$f(B,T)=-N_{c}\,\frac{qB}{2\pi}\sum_{k=0}^{\infty}\sum_{s=\pm
1/2}T\sum_{n=-\infty}^{\infty}\int\frac{\textmd{d}p}{2\pi}\log\frac{\omega_{n}^{2}+E_{p,s,k}^{2}}{T^{2}}\,,$
(49)

where $\omega_{n}=(2n+1)\pi T$ is the $n$-th fermionic Matsubara frequency.
Moreover, $p$, $s$ and $k$ are the momentum, spin and angular momentum in the
direction of the magnetic field, $N_{c}=3$ is the number of colors and the
energies are given by the Landau levels,

$E_{p,s,k}=\sqrt{p^{2}+m^{2}+(2k+1-2s)qB}\,.$ (50)

Rewriting the logarithm using Eq. (48), the integral over $p$ becomes Gaussian
and can be solved. Furthermore, the sums over $n$, $k$ and $s$ are

$\displaystyle
T\\!\sum_{n=-\infty}^{\infty}\\!\\!e^{-\omega_{n}^{2}t}=\frac{1}{2\sqrt{\pi
t}}\,\Theta_{3}\\!\left(\frac{\pi}{2},e^{-1/(4tT^{2})}\right)$
$\displaystyle,\;\,\sum_{k=0}^{\infty}e^{-(2k+1)qB\,t}=\frac{1}{2\sinh(qBt)},$
$\displaystyle\sum_{s=\pm 1/2}\\!\\!e^{-2sqB\,t}$
$\displaystyle=2\cosh(qBt)\,,$ (51)

where $\Theta_{3}$ is an elliptic function. Inserting these in Eq. (49) and
performing the derivative with respect to $z$, we obtain

$f(B,T)=N_{c}\,\frac{qB}{8\pi^{2}}\int_{0}^{\infty}\frac{\textmd{d}t}{t^{2}}\,e^{-m^{2}t}\,\coth(qBt)\;\Theta_{3}\\!\left(\frac{\pi}{2},e^{-1/(4tT^{2})}\right)\,.$
(52)

Taking the second derivative with respect to $eB$ to obtain the bare magnetic
susceptibility (1) results in

$\chi_{b}(T)=-\frac{N_{c}}{12\pi^{2}}(q/e)^{2}\int_{0}^{\infty}\frac{\textmd{d}t}{t}e^{-m^{2}t}\,\;\Theta_{3}\\!\left(\frac{\pi}{2},e^{-1/(4tT^{2})}\right)\,.$
(53)

###  C.2 Ultraviolet divergences and QED renormalization

To determine the ultraviolet structure of the magnetic susceptibility, we
consider Eq. (53) at zero temperature. For $T=0$ the elliptic function
$\Theta_{3}$ approaches unity. The resulting expression needs to be
regularized, for example by setting an ultraviolet cut-off $1/\Lambda^{2}$ as
the lower limit of the proper time integral. Performing the integral and
expanding for large $\Lambda$ we obtain,

$\chi_{b}(T=0)=\frac{N_{c}}{12\pi^{2}}(q/e)^{2}\left[\log\frac{\Lambda^{2}}{m^{2}}-\gamma_{E}\right]+\mathcal{O}(\Lambda^{-2})\,,$
(54)

where $\gamma_{E}$ is the Euler-Mascheroni constant. Thus, the coefficient of
the logarithmic divergence indeed equals the lowest-order QED $\beta$-function
coefficient $\beta_{1}$ (for one quark flavor with electric charge $q$),
demonstrating the validity of Eq. (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 $\beta_{1}$ at the scale $1/a$ approach zero due to
asymptotic freedom (see Eq. (23)). We note moreover that in the proper time
formulation the renormalization scale is set by the mass – in fact $\mu_{\rm
QED}=m\,e^{\gamma_{E}/2}$ for our choice of the regulator $\Lambda$ –
explaining the appearance of $m$ in the argument of the logarithm in Eq. (54).

The additive renormalization can be performed by subtracting $\chi_{b}(T=0)$
from Eq. (53):

$\chi=\chi_{b}(T)-\chi_{b}(T=0)=-\frac{N_{c}}{12\pi^{2}}(q/e)^{2}\int_{0}^{\infty}\frac{\textmd{d}t}{t}e^{-m^{2}t}\,\;\left[\Theta_{3}\\!\left(\frac{\pi}{2},e^{-1/(4tT^{2})}\right)-1\right]\,.$
(55)

As we mentioned after Eq. (6), this corresponds to the choice of a physical,
albeit scheme-dependent, QED renormalization scale.

###  C.3 Spin contribution

The contribution of orbital angular momentum to the total susceptibility can
be calculated by simply replacing the fermion with two ghost particles (spin-
zero but antiperiodic in Euclidean time) in the above calculation. This
removes the $-2sqB$ from the energies (50) and excludes the spin sum in the
free energy density (49). Consequently, the magnetic field-dependent part in
Eq. (52) changes as $\coth(qBt)\mapsto 1/\sinh(qBt)$. This merely changes the
second derivative of the free energy density at $B=0$ by a factor $-1/2$.
Thus, for the renormalized susceptibility we arrive at

$\chi^{\rm ang}(T)=-\frac{1}{2}\cdot\chi(T)\,,$ (56)

which also implies

$\chi^{\rm spin}(T)=\frac{3}{2}\cdot\chi(T)\,,$ (57)

confirming the $3:(-1)$ ratio of the two contributions to the total
susceptibility. We mention that a similar argument has been used in
perturbative QCD (with chromomagnetic background fields) to relate asymptotic
freedom to spin effects Nielsen:1980sx .

###  C.4 Tensor bilinear

For the tensor bilinear we begin with the result of the fermionic path
integral,

$\left\langle\bar{\psi}\sigma_{12}\psi\right\rangle=\frac{T}{V}\,\textmd{tr}\,\frac{\sigma_{12}}{\not{D}+m}=\frac{T\,m}{V}\,\textmd{tr}\,\frac{\sigma_{12}}{-\not{D}^{2}+m^{2}}\,,$
(58)

where we used chiral symmetry $\\{\gamma_{5},\not{D}\\}=0$. The trace is
represented using the eigenbasis of $-\not{D}^{2}$, giving the eigenvalues
$\omega_{n}^{2}+E_{p,s,k}^{2}$. Since $[\not{D}^{2},\sigma_{12}]=0$, the spin
operator can also be diagonalized in this basis and its eigenvalues are minus
two times the spin: $\sigma_{12}\to-2s$. Taking into account the
$2N_{c}\cdot(qBL^{2})/(2\pi)$-fold degeneracy of the eigenvalues, we obtain

$\left\langle\bar{\psi}\sigma_{12}\psi\right\rangle=N_{c}\frac{qB\,m}{\pi}\sum_{k=0}^{\infty}\sum_{s=\pm
1/2}T\sum_{n=-\infty}^{\infty}\int\frac{\textmd{d}p}{2\pi}\frac{-2s}{\omega_{n}^{2}+p^{2}+m^{2}+(2k+1-2s)qB}\,.$
(59)

In the sum the contributions $\\{k,s=1/2\\}$ and $\\{k+1,s=-1/2\\}$ cancel,
leaving only the unpaired lowest Landau level $\\{k=0,s=1/2\\}$. Hence we get

$\left\langle\bar{\psi}\sigma_{12}\psi\right\rangle=-N_{c}\frac{qB\,m}{\pi}\,T\sum_{n=-\infty}^{\infty}\int\frac{\textmd{d}p}{2\pi}\frac{1}{\omega_{n}^{2}+p^{2}+m^{2}}\,.$
(60)

Note that, unlike in full QCD, here the tensor bilinear is exactly linear in
the magnetic field. Thus, the tensor coefficient $\tau_{b}$ of Eq. (8) is
obtained by simply dividing Eq. (60) by $qB$.

Using Eq. (46) with $E=\sqrt{\omega_{n}^{2}+p^{2}+m^{2}}$, performing the
Gaussian integral over $p$ and the Matsubara sum (51) over $\omega_{n}$, we
arrive at

$\tau_{b}(T)=-N_{c}\frac{m}{4\pi^{2}}\int_{0}^{\infty}\frac{\textmd{d}t}{t}\,e^{-m^{2}t}\;\Theta_{3}\\!\left(\frac{\pi}{2},e^{-1/(4tT^{2})}\right)\,.$
(61)

A comparison to Eq. (53) reveals that this quantity contains the same
logarithmic divergence as $\chi_{b}$, just with a different coefficient. Using
a cut-off regulator as in Eq. (54), we obtain at $T=0$,

$\tau_{b}(T=0)=\frac{N_{c}}{4\pi^{2}}\,m\left[\log\frac{\Lambda^{2}}{m^{2}}-\gamma_{E}\right]+\mathcal{O}(\Lambda^{-2})\,,$
(62)

confirming Eq. (9). The same considerations regarding QCD corrections to the
coefficient and the renormalization scale $\mu_{\rm QED}$ apply as in Sec. C.2
for $\chi_{b}$.

The difference $\tau=\tau_{b}(T)-\tau_{b}(T=0)$ is ultraviolet-finite. We can
compare this with Eqs. (55) and (57) to conclude that

$\frac{(q/e)^{2}}{2m}\left[\tau(m)-\tau(m\to 0)\right]=\chi^{\rm spin}\,,$
(63)

confirming the relation (12) and Eq. (57). Notice that $\tau$ vanishes for
$m\to 0$, so the subtraction of the massless limit is irrelevant in the free
case (but it is relevant for the interacting system with spontaneous chiral
symmetry breaking, see App. B).

###  C.5 High-temperature expansion

The temperature-dependent part of the free energy density (49) can be
simplified using the well-known trick kapusta2006finite  of differentiating
and subsequently integrating the integrand with respect to $E_{p,s,k}$. The
result is

$f(B,T)-f(B,0)=-2N_{c}\frac{qB}{2\pi}\sum_{k=0}^{\infty}\sum_{s=\pm
1/2}\int\frac{\textmd{d}p}{2\pi}\,T\log\left[1+e^{-E_{p,s,k}/T}\right]\,.$
(64)

The energy levels are given in Eq. (50) above. To obtain the high-temperature
expansion in a closed form, we need to replace the logarithm by its series
expansion

$\log(1+x)=-\sum_{l=1}^{\infty}\frac{(-x)^{l}}{l}\,.$ (65)

This approach was used, e.g., in Ref. Toms:1996dg  for scalars at nonzero
chemical potential.

Inserting the expansion (65) into (64) and rewriting the exponentials using
Eq. (46) results in

$f(B,T)-f(B,0)=N_{c}\frac{qB\,T}{2\pi^{2}}\sum_{l=1}^{\infty}\frac{(-1)^{l}}{l}\sum_{k=0}^{\infty}\sum_{s=\pm
1/2}\int\textmd{d}p\,\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\\!\\!\textmd{d}z\,\Gamma(z)\,l^{-z}E_{p,s,k}^{-z}\,T^{z}\,.$
(66)

Inserting the Mellin transform (46) for $E_{p,s,k}^{-z}$ renders the integral
over $p$ Gaussian. We can reuse the angular momentum and spin-sums from Eq.
(51), giving

$f(B,T)-f(B,0)=N_{c}\frac{qB}{2\pi^{3/2}}\,\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\\!\\!\\!\textmd{d}z\,\frac{\Gamma(z)}{\Gamma(z/2)}\,T^{z+1}\sum_{l=1}^{\infty}\frac{(-1)^{l}}{l^{1+z}}\int_{0}^{\infty}\\!\\!\textmd{d}t\,t^{(z-3)/2}\,e^{-m^{2}t}\,\coth(qBt)\,.$
(67)

Differentiating the above expression twice with respect to $eB$ at $B=0$ gives
(minus) the renormalized magnetic susceptibility $\chi$. The integral over $t$
can be solved via the Mellin transform (46) and gives a $\Gamma$-function,
while the sum over $l$ results in a $\zeta$-function:

$\sum_{l=1}^{\infty}\frac{(-1)^{l}}{l^{1+z}}=\zeta(1+z)\cdot(2^{-z}-1)\,.$
(68)

Using the duplication formula NIST:DLMF  for the ratio of $\Gamma$-functions,
we arrive at

$\chi=-\frac{N_{c}}{6\pi^{2}m}(q/e)^{2}\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\textmd{d}z\,\Gamma\left(\frac{z+1}{2}\right)\Gamma\left(\frac{z+1}{2}\right)\zeta(1+z)\,(1-2^{z})\,m^{-z}\,T^{z+1}\,.$
(69)

For the validity of the Mellin transforms we needed to assume $c>0$ (as well
as $m>0$). The final integral over $z$ can be solved using Cauchy’s theorem,
closing the integral towards the left and calculating the residue at the
poles. There is a double pole at $z=-1$ and simple poles at $z=-3,-5,\ldots$.
These $z$-values set the powers of $T$ that appear in the high-temperature
expansion. Keeping the leading terms (i.e. $z=-1$ and $z=-3$), we finally
obtain,

$\chi=\frac{N_{c}}{12\pi^{2}}(q/e)^{2}\left[\log\frac{T^{2}\pi^{2}}{m^{2}}-2\gamma_{E}+\frac{7\,\zeta(3)}{4\pi^{2}}\,\frac{m^{2}}{T^{2}}\right]+\mathcal{O}(m^{4}/T^{4})\,,$
(70)

reproducing the results of Ref. Cangemi:1996tp . Notice that the coefficient
of the leading logarithmic term is equal to $\beta_{1}$ (for one flavor with
electric charge $q$), confirming Eq. (25), in agreement with Refs.
Loewe:1991mn ; Elmfors:1993bm ; Gies:1999xn . As we have seen below Eq. (54),
in the proper time formulation the renormalization scale is intrinsically set
by the quark mass, $\mu_{\rm QED}=m\,e^{\gamma_{E}/2}$. We may express the
square bracket in the leading term as $\log(\gamma\,T^{2}/\mu_{\rm QED}^{2})$
with $\gamma=\pi^{2}\,e^{-\gamma_{E}}$. Clearly, $\gamma=\mathcal{O}(1)$
depends on the definition of the regulator. The general form is again expected
to hold in full QCD Bali:2014kia : in this case QCD corrections at scales
$T\gg\mu_{\rm QED}$ are small due to asymptotic freedom.

##  Appendix D Multiplicative QCD renormalization

Since lattice perturbation theory is slowly convergent and high-loop results
are unavailable, we first match the local lattice QCD operators of interest
non-perturbatively to the regulator independent RI’-MOM scheme
Martinelli:1994ty ; Chetyrkin:1999pq  and subsequently translate the result at
three-loop order Gracey:2003yr  to the $\overline{\rm MS}$ scheme.

The quark bilinear operators are renormalized by computing the corresponding
amputated flavor non-singlet vertex functions for different momenta on Landau
gauge-fixed ensembles. We wish to renormalize light- and strange-quark
bilinears, which can be written as linear combinations of the diagonal SU(3)
flavor-octet and -singlet currents. In continuum schemes, with the exception
of the axial current that we do not discuss here, the renormalization of
flavor singlet and non-singlet operators of dimension three is the same. This
also appears to hold for the staggered action Constantinou:2016ieh ;
Lee:1999zxa . We remark that, instead of extrapolating to the $N_{f}=3$
massless case, we use physical quark masses, which may be problematic, in
particular regarding the strange quark mass. However, in Ref. Gockeler:2010yr
it was demonstrated that the effect of the mass-dependence is tiny for the
perturbative momentum transfers that we are interested in. Moreover, the
difference is expected to vanish after a continuum limit extrapolation of a
renormalized matrix element is carried out because our quark masses are tuned
to a line of constant physics.

Since the spin degrees of freedom are spread over hypercubes for staggered
fermions, the determination of the vertex function in momentum space is more
challenging than for Wilson fermions. We follow the approach described in Ref.
Lytle:2013qoa : the taste and spin degrees of freedom are reconstructed from
different momentum combinations. The quark propagator for a given momentum, as
any vertex function, will be a matrix of size $16\times 16$, after averaging
over the color degrees of freedom. Our choice of the scalar and tensor
currents, where, in the latter case, we employ a two-link operator, is
detailed in Ref. Bali:2012jv .

$\beta$ |  $a$/fm | $Z_{T}$ | $Z_{S}$
---|---|---|---
3.45 | 0.282 | 1.07(12) | 1.14(17)
3.55 | 0.217 | 1.114(45) | 0.829(12)
3.67 | 0.153 | 1.125(19) | 0.788(41)
3.75 | 0.125 | 1.123(19) | 0.723(38)
3.85 | 0.099 | 1.100(18) | 0.660(34)
Table 1: Conversion factors to the $\overline{\rm MS}$ scheme at $\mu_{\rm
QCD}=2\textmd{ GeV}$. Figure 11: Multiplicative renormalization constants as a
function of $\beta$. The symbols have been slightly shifted horizontally for
better visibility and connected by lines to guide the eye. Also shown as
dashed lines are the one-loop perturbative expectations Bali:2012jv  that will
be approached as $\beta\rightarrow\infty$.

The error of the final renormalization constants is dominated by systematics.
On the one hand, the conversion factors from the RI’-MOM to the
$\overline{\textrm{MS}}$-scheme are only known up to a fixed order in
perturbation theory (three loops in our case). Hence high momenta are
preferable. On the other hand, at momentum scales close to the lattice cut-off
the intermediate matching to the RI’-MOM scheme will significantly be affected
by lattice artifacts. Therefore, we are restricted to a “window” of
intermediate momentum values. We employ combinations along symmetric lattice
directions, where the lattice corrections are smallest. Another complication
is that due to the choice of the staggered action, the maximum momentum scale
that can be achieved on a four-dimensional lattice is $\pi/a$, rather than
$2\pi/a$. As a compromise, on the finest three lattices we interpolate the
RI’-MOM result to the fixed scale $\mu_{\rm QCD}=2\textmd{ GeV}$.
Subsequently, this is perturbatively converted to the $\overline{\rm
MS}$-scheme. We estimate the uncertainty by adding the difference between the
scheme conversion at two- and at three-loop order and the (statistical and
systematic) interpolation uncertainty in quadrature. The latter contribution
is negligible. At the coarsest two lattice spacings, $\mu_{\rm QCD}=2\textmd{
GeV}$ is too close to the cut-off scale to obtain reliable results. Therefore,
at $a\approx 0.28\textmd{ fm}$ and at $a\approx 0.22\textmd{ fm}$, we convert
the RI’-MOM results at $\mu_{\rm QCD}=1.1\textmd{ GeV}$ and at $\mu_{\rm
QCD}=1.5\textmd{ GeV}$, respectively, to the $\overline{\rm MS}$-scheme and
evolve the result to $\mu_{\rm QCD}=2\textmd{ GeV}$. We replicate the same
procedure at $a\approx 0.15\textmd{ fm}$ and add the difference that we obtain
at this lattice spacing between the matching at $\mu_{\rm QCD}=2\textmd{ GeV}$
and the matching at these lower scales in quadrature to the systematic error
at $\mu_{\rm QCD}=1.1\textmd{ GeV}$ and $\mu_{\rm QCD}=1.5\textmd{ GeV}$.

The results are listed in Table 1 and shown in Fig. 11. We also include the
lattice perturbative theory one-loop expectations Bali:2012jv  in the figure.
The comparatively larger value of $Z_{T}$ results in a larger modulus of the
renormalized tensor coefficient.

##  Appendix E Susceptibilities via current-current correlators

Here we derive Eqs. (18) and (22) of the main text. To this end we consider a
background field that possesses nonzero momentum $p_{1}$ in the $x_{1}$
direction. The constant field setup will be approached via the $p_{1}\to 0$
limit. This approach has been described in detail in Ref. Bali:2015msa  for
$\chi_{b}$ in momentum space. Here we repeat the argument in coordinate space
and also generalize it for $\tau_{fb}$.

###  E.1 Magnetic susceptibility from correlators

We consider an oscillatory magnetic field and the corresponding Landau-gauge
vector potential,

$B(x_{1})=B\cdot\cos(p_{1}x_{1}),\quad\quad
A_{2}(x_{1})=B\cdot\frac{\sin(p_{1}x_{1})}{p_{1}}\,.$ (71)

The latter couples to $i\cdot e$ times the current (16) in the action density.
We can define the associated susceptibility just like in Eq. (1),

$\chi_{b}^{p_{1},{\rm
cos}}=-\left.\frac{\partial^{2}f}{\partial(eB)^{2}}\right|_{B=0}=-\frac{T}{V}\int\textmd{d}^{4}y\,\textmd{d}^{4}z\,\frac{\sin(p_{1}y_{1})}{p_{1}}\frac{\sin(p_{1}z_{1})}{p_{1}}\left\langle
j_{2}(y)j_{2}(z)\right\rangle\,,$ (72)

where each derivative brought down an integral over the current $j_{2}$ times
the coordinate-dependence of $A_{2}$ and we used $\left\langle
j_{2}\right\rangle=0$. Changing the integration variable from $z$ to $x=z-y$
and exploiting the translational invariance of the current-current correlator,
the integrals over $y_{2}$, $y_{3}$ and $y_{4}$ can be carried out,

$\chi_{b}^{p_{1},{\rm
cos}}=-\frac{1}{L}\int\textmd{d}y_{1}\,\textmd{d}x_{1}\,\frac{\sin(p_{1}y_{1})\sin(p_{1}(y_{1}+x_{1}))}{p_{1}^{2}}\,G(x_{1})\,,$
(73)

where the projected correlator $G(x_{1})$, defined in Eq. (17), appears. In
the $p_{1}\to 0$ limit, $B(x_{1})$ becomes homogeneous and
$\chi_{b}^{p_{1},{\rm cos}}$ equals the ordinary susceptibility $\chi_{b}$.

For reasons that will become clear in a moment, let us consider a different
background field,

$B(x_{1})=B\cdot\sin(p_{1}x_{1}),\quad\quad
A_{2}(x_{1})=-B\cdot\frac{\cos(p_{1}x_{1})}{p_{1}}\,,$ (74)

for which the associated oscillatory susceptibility, similarly to Eq. (73),
reads

$\chi_{b}^{p_{1},{\rm
sin}}=-\frac{1}{L}\int\textmd{d}y_{1}\,\textmd{d}x_{1}\,\frac{\cos(p_{1}y_{1})\cos(p_{1}(y_{1}+x_{1}))}{p_{1}^{2}}\,G(x_{1})\,.$
(75)

In this case the $p_{1}\to 0$ limit does not reproduce $\chi_{b}$. Instead,
$A_{2}(x_{1})$ becomes homogeneous: it acts as if we had introduced a constant
imaginary ‘chemical potential’ in the $x_{2}$ direction, with magnitude
$\mu_{2}=-eB/p_{1}$. Therefore the oscillatory susceptibility becomes
proportional to the leading response to this spatial chemical potential,

$\chi_{b}^{p_{1},{\rm sin}}\xrightarrow{p_{1}\to
0}\frac{c_{2}}{p_{1}^{2}},\quad\quad
c_{2}=-\frac{1}{L}\int\textmd{d}y_{1}\,\textmd{d}x_{1}\,G(x_{1})\,.$ (76)

This detour was necessary to simplify the $p_{1}\to 0$ limit of the
oscillatory susceptibilities. Specifically, let us examine the following
combination:

$\displaystyle\chi_{b}^{p_{1},{\rm cos}}+\chi_{b}^{p_{1},{\rm
sin}}-\frac{c_{2}}{p_{1}^{2}}$
$\displaystyle\qquad=-\frac{1}{L}\int\textmd{d}y_{1}\,\textmd{d}x_{1}\,\frac{\sin(p_{1}y_{1})\sin(p_{1}(y_{1}+x_{1}))+\cos(p_{1}y_{1})\cos(p_{1}(y_{1}+x_{1}))-1}{p_{1}^{2}}\,G(x_{1})\,.$
(77)

This approaches $\chi_{b}$ for $p_{1}\to 0$. Using the trigonometric identity
for the cosine of the difference of angles in the numerator of the kernel
reveals that the integrand is independent of $y_{1}$. (This is why we needed
to consider both the $\cos$\- and $\sin$-type fields.) Integrating over
$y_{1}$ cancels the prefactor $1/L$, resulting in

$\chi_{b}=-\lim_{p_{1}\to
0}\int\textmd{d}x_{1}\,\frac{\cos(p_{1}x_{1})-1}{p_{1}^{2}}\,G(x_{1})=\int\textmd{d}x_{1}\,\frac{x_{1}^{2}}{2}\,G(x_{1})\,,$
(78)

where we finally performed the $p_{1}\to 0$ limit. This proves Eq. (18) of the
main text. We note that the crucial point of the derivation was Eq. (77),
where the kernel was shown to only depend on the distance $x_{1}$ between the
two current insertions. This was done conveniently using the combination of
the oscillatory fields – in contrast, it would have been more tedious if we
started directly with a constant background.

###  E.2 Tensor coefficient from correlators

We generalize the above derivation for $\tau_{fb}$, which can be written as

$\tau_{fb}=\frac{1}{q_{f}/e}\left.\frac{\partial}{\partial(eB)}\right|_{B=0}\frac{T}{V}\int\textmd{d}^{4}x\,\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}(x)\right\rangle\,.$
(79)

Again we consider oscillatory magnetic fields of the types (71) and (74).
These give rise to modulated tensor bilinears of the forms
$\bar{\psi}_{f}\sigma_{12}\psi_{f}(x)\cos(p_{1}x_{1})$ and
$\bar{\psi}_{f}\sigma_{12}\psi_{f}(x)\sin(p_{1}x_{1})$, respectively, which
enter the corresponding oscillatory tensor coefficients $\tau_{fb}^{p_{1},{\rm
cos}}$ and $\tau_{fb}^{p_{1},{\rm sin}}$:

$\begin{split}\tau_{fb}^{p_{1},{\rm
cos}}&=\frac{i}{q_{f}/e}\,\frac{1}{L}\int\textmd{d}y_{1}\,\textmd{d}x_{1}\,\cos(p_{1}y_{1})\,\frac{\sin(p_{1}(y_{1}+x_{1}))}{p_{1}}\,H_{f}(x_{1})\,,\\\
\tau_{fb}^{p_{1},{\rm
sin}}&=\frac{-i}{q_{f}/e}\,\frac{1}{L}\int\textmd{d}y_{1}\,\textmd{d}x_{1}\,\sin(p_{1}y_{1})\,\frac{\cos(p_{1}(y_{1}+x_{1}))}{p_{1}}\,H_{f}(x_{1})\,,\end{split}$
(80)

where the projected tensor-vector correlator $H_{f}(x_{1})$, defined in Eq.
(22), appears. Here we performed the same variable substitution as in Eq. (73)
above.

Figure 12:  Comparison of different methods to calculate $\tau_{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).

In this case, $\tau_{fb}^{p_{1},{\rm cos}}$ approaches $\tau_{fb}$ for
$p_{1}\to 0$, while $\tau_{fb}^{p_{1},{\rm sin}}$ vanishes in that limit. Thus
we need to consider the sum of the two coefficients. Employing the
trigonometric identity for the sine of the difference of angles and carrying
out the integral over $y_{1}$ gives

$\tau_{fb}=\lim_{p_{1}\to 0}\left[\tau_{fb}^{p_{1},{\rm
cos}}+\tau_{fb}^{p_{1},{\rm sin}}\right]=\lim_{p_{1}\to
0}\frac{i}{q_{f}/e}\,\int\textmd{d}x_{1}\,\frac{\sin(p_{1}x_{1})}{p_{1}}\,H_{f}(x_{1})=\frac{i}{q_{f}/e}\int\textmd{d}x_{1}\,x_{1}\,H_{f}(x_{1})\,.$
(81)

In finite volumes we carry out the same symmetrization as in Eq. (19), this
time taking into account that $H_{f}(x_{1})=-H_{f}(L-x_{1})$ to finally arrive
at Eq. (22) of the main text.

This method to calculate $\tau_{fb}$ is compared to the results for
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$ measured at $B>0$
on $24^{3}\times 6$ lattices at $T=113\textmd{ MeV}$ in Fig. 12. For the light
quarks we obtain consistent results, however, for $\tau_{sb}$ the correlator
tends to give values that slightly differ from the slope of a linear fit to
the lowest few available points. Since lattice artifacts and finite volume
effects might be different in the two cases, such slight differences are not
unexpected.

In addition, we find that a linear fit to results from simulations at $B>0$
has smaller uncertainties than extracting the slope at $B=0$ using the
correlator method. In the main text we therefore use our earlier results for
$\left\langle\bar{\psi}_{f}\sigma_{12}\psi_{f}\right\rangle$ from Ref.
Bali:2012jv .

We note that the tensor-vector correlators at nonzero spatial momenta might
also be useful for extracting further features of the photon distribution
amplitude.

##  Appendix F Parametrization of the equation of state

Up to $\mathcal{O}(B^{2})$, the magnetic field-dependence of the complete EoS
can be calculated from the magnetic susceptibility $\chi(T)$. Here we provide
a parametrization for this observable and also collect the relevant
thermodynamical relations, which were also summarized in Ref. Bali:2014kia .

First of all, we remind the reader that in the presence of a background
magnetic field, the different components of the pressure – defined by
considering an infinitesimal compression of the system in the respective
direction – might become anisotropic Bali:2013esa . In particular, one should
distinguish between the _$\Phi$ -scheme_, where the flux of the magnetic field
is kept constant during the compression (superscript $(\Phi)$ below), and the
_$B$ -scheme_, where the magnetic field strength is kept constant (superscript
$(B)$). On the one hand, the $B$-scheme pressure is isotropic and equals the
negative of the free energy density in the thermodynamic limit,

$p_{1,2}^{(B)}=p_{3}=-f\,.$ (82)

On the other hand, in the $\Phi$-scheme the pressure components are related by
the magnetization $\mathcal{M}$,

$p_{1,2}^{(\Phi)}=p_{3}-eB\cdot\mathcal{M},\quad\quad\mathcal{M}=-\frac{\partial
f}{\partial(eB)}\,.$ (83)

The entropy density $s$ and the energy density $\epsilon$ are scheme-
independent,

$s=-\frac{\partial f}{\partial T},\quad\quad\epsilon=f+Ts\,,$ (84)

whereas also the interaction measure (trace anomaly) $I$ differs between the
two schemes,

$I^{(B)}=\epsilon-3p_{3},\quad\quad
I^{(\Phi)}=\epsilon-2p_{1,2}^{(\Phi)}-p_{3}=I^{(B)}+2eB\cdot\mathcal{M}\,.$
(85)

Using Eqs. (1) and (82), the leading-order expansion in the magnetic field
takes the form

$p_{3}(T,B)=p_{3}(T,0)+\chi(T)\,\frac{(eB)^{2}}{2},\quad\quad\mathcal{M}(T,B)=\chi(T)\,eB\,.$
(86)

Together with Eqs. (82)–(85) these specify the $B$-dependence of all relevant
observables up to $\mathcal{O}(B^{2})$.

At $B=0$ the pressure is isotropic, and can be obtained from the interaction
measure as141414We note that Eq. (87)) remains valid also for $B>0$ in the
$B$-scheme but not in the $\Phi$-scheme.

$\frac{p(T,B=0)}{T^{4}}=\int_{0}^{T^{\prime}}\\!\\!\textmd{d}T^{\prime}\,\frac{I(T^{\prime},B=0)}{T^{5}}\,.$
(87)

Thus, to calculate the complete EoS including $B^{0}$ and $B^{2}$ effects,
altogether it suffices to parameterize $I(T,0)$ and $\chi(T)$. For the latter
we consider a parametrization of the continuum extrapolated lattice results
that smoothly approach the HRG model prediction (see Fig. 4) at low and the
perturbation theory formula (25)) at high temperatures. We found the following
parametric form to be sufficient for this,

$\chi(T)=\exp(-h_{3}/t)\cdot\frac{1+g_{0}/t+g_{1}/t^{2}+g_{2}/t^{3}}{1+g_{3}/t+g_{4}/t^{2}+g_{5}/t^{3}}\cdot
2\beta_{1}\log\frac{t}{q_{0}}\,,\quad\quad t=\frac{T}{1\textmd{ GeV}}\,.$ (88)

Eq. (88) 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 $\beta_{1}$ coefficient is fixed to its perturbative value
(4), while the scale $q_{0}$ inside the logarithm is allowed to be a free
parameter. The rational function involving the $g_{i}$ parameters interpolates
between the two limiting behaviors. The so-obtained parametrization is shown
in the left panel of Fig. 5 in the main text.

For the interaction measure we take the parametrization of Ref.
Borsanyi:2013bia ,

$\frac{I(T,0)}{T^{4}}=\exp(-h_{1}/t-h_{2}/t^{2})\cdot\left(h_{0}+\frac{f_{0}\cdot\left[\tanh(f_{1}\cdot
t+f_{2})+1\right]}{1+k_{1}\cdot t+k_{2}\cdot t^{2}}\right)\,,\quad\quad
t=\frac{T}{0.2\textmd{ GeV}}\,.$ (89)

The parameters of both functions are included in Table. 2. The two
parametrizations, together with the implementations of the formulae (82)–(87)
are included in the Python script param_EoS.py that is submitted to arXiv.org
together with this manuscript.

$\beta_{1}$ | $h_{3}$ | $g_{0}$ | $g_{1}$ | $g_{2}$ | $g_{3}$ | $g_{4}$ | $g_{5}$ | $q_{0}$
---|---|---|---|---|---|---|---|---
$1/(6\pi^{2})$ | 0.1544 | 23.99 | -2.085 | 0.1290 | 21.35 | -6.201 | 0.5766 | 0.1497

$h_{0}$ | $h_{1}$ | $h_{2}$ | $f_{0}$ | $f_{1}$ | $f_{2}$ | $k_{1}$ | $k_{2}$
---|---|---|---|---|---|---|---
0.1396 | -0.1800 | 0.0350 | 1.05 | 6.39 | -4.72 | -0.92 | 0.57

Table 2: Parameters of the functions (88) and (89).

This parametrization is valid for low magnetic fields. To be more
quantitative, we compare our $\mathcal{O}(B^{2})$ truncated results for the
longitudinal pressure to the complete magnetic field-dependence from Ref.
Bali:2014kia  for $T\gtrsim 180\textmd{ MeV}$. We find agreement within errors
in the range $B/(\pi T)^{2}\lesssim 1$. This upper limit is hard-coded in the
Python script as well. One final remark about the parametrization is in order.
All truncated thermodynamic observables approach zero for $T\to 0$, such that
a normalization by the corresponding powers of the temperature (i.e.
$p_{3}/T^{4}$, $s/T^{3}$ and so on) produces sensible plots. This is not the
case if $\mathcal{O}(B^{4})$ terms are also included: at this order vacuum
contributions arise and the equation of state depends on $B$ already at $T=0$,
rendering a normalization like $p_{3}/T^{4}$ ill-defined in the $T\to 0$ limit
Bali:2014kia .

## References

  * (1) K. Kiuchi, P. Cerdá-Durán, K. Kyutoku, Y. Sekiguchi and M. Shibata, _Efficient magnetic-field amplification due to the Kelvin-Helmholtz instability in binary neutron star mergers_ , _Phys. Rev._ D92 (2015) 124034, [1509.09205].
  * (2) L. Baiotti and L. Rezzolla, _Binary neutron star mergers: a review of Einstein’s richest laboratory_ , _Rept. Prog. Phys._ 80 (2017) 096901, [1607.03540].
  * (3) T. Kawamura, B. Giacomazzo, W. Kastaun, R. Ciolfi, A. Endrizzi, L. Baiotti et al., _Binary neutron star mergers and short gamma-ray bursts: effects of magnetic field orientation, equation of state, and mass ratio_ , _Phys. Rev._ D94 (2016) 064012, [1607.01791].
  * (4) D. Kharzeev, K. Landsteiner, A. Schmitt and H.-U. Yee, _Strongly interacting matter in magnetic fields_ , _Lect. Notes Phys._ 871 (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._ 88 (2016) 025001, [1411.7176].
  * (6) V. A. Miransky and I. A. Shovkovy, _Quantum field theory in a magnetic field: From quantum chromodynamics to graphene and Dirac semimetals_ , _Phys. Rept._ 576 (2015) 1–209, [1503.00732].
  * (7) J. Lattimer and M. Prakash, _Neutron star structure and the equation of state_ , _Astrophys. J._ 550 (2001) 426, [astro-ph/0002232].
  * (8) D. Grasso and H. R. Rubinstein, _Magnetic fields in the early universe_ , _Phys. Rept._ 348 (2001) 163–266, [astro-ph/0009061].
  * (9) R. Durrer and A. Neronov, _Cosmological magnetic fields: their generation, evolution and observation_ , _Astron. Astrophys. Rev._ 21 (2013) 62, [1303.7121].
  * (10) D. E. Kharzeev, J. Liao, S. A. Voloshin and G. Wang, _Chiral magnetic and vortical effects in high-energy nuclear collisions – A status report_ , _Prog. Part. Nucl. Phys._ 88 (2016) 1–28, [1511.04050].
  * (11) X.-D. Ji, _Gauge-invariant decomposition of nucleon spin_ , _Phys. Rev. Lett._ 78 (1997) 610–613, [hep-ph/9603249].
  * (12) R. L. Jaffe and A. Manohar, _The $G_{1}$ 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 a new 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 sum rule method_ , _Yad. Fiz._ 40 (1984) 1035–1038.
  * (16) I. I. Balitsky, A. V. Kolesnichenko and A. V. Yung, _On vector dominance in sum rules for electromagnetic hadron characteristics (in Russian)_ , _Sov. J. Nucl. Phys._ 41 (1985) 178.
  * (17) I. I. Balitsky, V. M. Braun and A. V. Kolesnichenko, _Radiative decay $\Sigma^{+}\to p\gamma$ in Quantum Chromodynamics_, _Nucl. Phys._ B312 (1989) 509–550.
  * (18) V. Yu. Petrov, M. V. Polyakov, R. Ruskov, C. Weiss and K. Goeke, _Pion and photon light cone wave functions from the instanton vacuum_ , _Phys. Rev._ D59 (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 on pion-exchange_ , _Phys. Rev._ D79 (2009) 073012, [0901.1172].
  * (21) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, _Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment_ , _Phys. Lett._ B798 (2019) 134994, [1908.03331].
  * (22) P. Ball and E. Kou, _$B\to\gamma e\nu$ transitions from QCD sum rules on the light cone_, _JHEP_ 04 (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 $\gamma^{*}\gamma\to\eta$ and $\gamma^{*}\gamma\to\eta^{\prime}$ in QCD_, _Phys. Rev._ D90 (2014) 074019, [1409.4311].
  * (25) G. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. Katz et al., _QCD quark condensate in external magnetic fields_ , _Phys. Rev._ D86 (2012) 071502, [1206.4205].
  * (26) G. S. Bali, F. Bruckmann, M. Constantinou, M. Costa, G. Endrődi, S. D. Katz et al., _Magnetic susceptibility of QCD at zero and at finite temperature from the lattice_ , _Phys. Rev._ D86 (2012) 094512, [1209.6015].
  * (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_ 1304 (2013) 130, [1303.1328].
  * (28) G. Bali, F. Bruckmann, G. Endrődi, S. Katz and A. Schäfer, _The QCD equation of state in background magnetic fields_ , _JHEP_ 1408 (2014) 177, [1406.0269].
  * (29) G. Bali and G. Endrődi, _Hadronic vacuum polarization and muon $g-2$ from magnetic susceptibilities on the lattice_, _Phys. Rev._ D92 (2015) 054506, [1506.08638].
  * (30) G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, _A general method for nonperturbative 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 the regularization invariant and $\overline{\mathrm{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._ 82 (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 $\alpha_{s}^{4}$ and the QED $\beta$-function at five loop_, _JHEP_ 1207 (2012) 017, [1206.1284].
  * (35) I. I. Balitsky and A. V. Yung, _Proton and neutron magnetic moments from QCD sum rules_ , _Phys. Lett._ 129B (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 nonsinglet perturbative renormalization of fermion bilinears_ , _Phys. Rev._ D94 (2016) 114513, [1610.06744].
  * (38) W. Pauli, _Über Gasentartung und Paramagnetismus_ , _Zeitschrift für Physik_ 41 (June, 1927) 81–102.
  * (39) L. Landau, _Diamagnetismus der Metalle_ , _Zeitschrift für Physik_ 64 (Sep, 1930) 629–637.
  * (40) Y. Aoki, Z. Fodor, S. Katz and K. Szabó, _The Equation of state in lattice QCD: With physical quark masses towards the continuum limit_ , _JHEP_ 0601 (2006) 089, [hep-lat/0510084].
  * (41) S. Borsányi, G. Endrődi, Z. Fodor, A. Jakovác, S. D. Katz et al., _The QCD equation of state with dynamical quarks_ , _JHEP_ 1011 (2010) 077, [1007.2580].
  * (42) G. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. Katz et al., _The QCD phase diagram for external magnetic fields_ , _JHEP_ 1202 (2012) 044, [1111.4956].
  * (43) G. Martinelli, G. Parisi, R. Petronzio and F. Rapuano, _The proton and neutron magnetic moments in lattice 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 magnetic field on a torus_ , _Annals Phys._ 324 (2009) 343–360, [0807.0630].
  * (45) G. S. Bali, F. Bruckmann, G. Endrődi and A. Schäfer, _Paramagnetic squeezing of QCD matter_ , _Phys. Rev. Lett._ 112 (2014) 042301, [1311.2559].
  * (46) C. Bonati, M. D’Elia, M. Mariti, F. Negro and F. Sanfilippo, _Magnetic susceptibility of strongly interacting matter across the deconfinement transition_ , _Phys. Rev. Lett._ 111 (2013) 182001, [1307.8063].
  * (47) C. Bonati, M. D’Elia, M. Mariti, F. Negro and F. Sanfilippo, _Magnetic susceptibility and equation of state of $N_{f}=2+1$ QCD with physical quark masses_, _Phys. Rev._ D89 (2014) 054506, [1310.8656].
  * (48) L. Levkova and C. DeTar, _Quark-gluon plasma in an external magnetic field_ , _Phys. Rev. Lett._ 112 (2014) 012002, [1309.1142].
  * (49) L. Levkova. private communication.
  * (50) M. Bellac, _Thermal Field Theory_.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 2000.
  * (51) J. Kapusta and C. Gale, _Finite-Temperature Field Theory: Principles and Applications_.  Cambridge monographs on mechanics and applied mathematics. Cambridge University Press, 2006.
  * (52) H. B. Meyer and H. Wittig, _Lattice QCD and the anomalous magnetic moment of the muon_ , _Prog. Part. Nucl. Phys._ 104 (2019) 46–96, [1807.09370].
  * (53) A. Francis, B. Jäger, H. B. Meyer and H. Wittig, _A new representation of the Adler function for lattice QCD_ , _Phys. Rev._ D88 (2013) 054502, [1306.2532].
  * (54) V. M. Braun et al., _The $\rho$-meson light-cone distribution amplitudes from lattice QCD_, _JHEP_ 04 (2017) 082, [1612.02955].
  * (55) G. Endrődi, _Multidimensional spline integration of scattered data_ , _Comput. Phys. Commun._ 182 (2011) 1307–1314, [1010.2952].
  * (56) P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Five-loop running of the QCD coupling constant_ , _Phys. Rev. Lett._ 118 (2017) 082002, [1606.08659].
  * (57) 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 $\Lambda$ parameter_, _Phys. Rev. Lett._ 119 (2017) 102001, [1706.03821].
  * (58) C. McNeile, A. Bazavov, C. T. H. Davies, R. J. Dowdall, K. Hornbostel, G. P. Lepage et al., _Direct determination of the strange and light quark condensates from full lattice QCD_ , _Phys. Rev._ D87 (2013) 034503, [1211.6577].
  * (59) Flavour Lattice Averaging Group collaboration, S. Aoki et al., _FLAG Review 2019_ , _Eur. Phys. J._ C80 (2020) 113, [1902.08191].
  * (60) P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Quark mass and field anomalous Dimensions to ${\cal O}(\alpha_{s}^{5})$_, _JHEP_ 10 (2014) 076, [1402.6611].
  * (61) D. J. Broadhurst and A. G. Grozin, _Matching QCD and HQET heavy-light currents at two loops and beyond_ , _Phys. Rev._ D52 (1995) 4082–4098, [hep-ph/9410240].
  * (62) A. Vainshtein, _Perturbative and nonperturbative renormalization of anomalous quark triangles_ , _Phys. Lett._ B569 (2003) 187–193, [hep-ph/0212231].
  * (63) A. Gorsky and A. Krikun, _Magnetic susceptibility of the quark condensate via holography_ , _Phys. Rev._ D79 (2009) 086015, [0902.1832].
  * (64) M. Frasca and M. Ruggieri, _Magnetic susceptibility of the quark condensate and polarization from chiral models_ , _Phys. Rev._ D83 (2011) 094024, [1103.1194].
  * (65) P. V. Buividovich, M. N. Chernodub, E. V. Luschevskaya and M. I. Polikarpov, _Chiral magnetization of non-Abelian vacuum: A lattice study_ , _Nucl. Phys._ B826 (2010) 313–327, [0906.0488].
  * (66) V. V. Braguta, P. V. Buividovich, T. Kalaydzhyan, S. V. Kuznetsov and M. I. Polikarpov, _The Chiral Magnetic Effect and chiral symmetry breaking in SU(3) quenched lattice gauge theory_ , _PoS_ LATTICE2010 (2010) 190, [1011.3795].
  * (67) S. Borsányi et al., _Is there still any $T_{c}$ mystery in lattice QCD? Results with physical masses in the continuum limit III_, _JHEP_ 1009 (2010) 073, [1005.3508].
  * (68) N. Bohr, _Studier over metallernes elektrontheori_ , _Thesis, Kobenhavns Universitet (translated and reprinted)_ (1911) .
  * (69) Van Leeuwen, Hendrika J., _Problèmes de la théorie électronique du magnétisme_ , _J. Phys. Radium_ 2 (1921) 361–377.
  * (70) B. Savoie, _A rigorous proof of the Bohr-van Leeuwen theorem in the semiclassical limit_ , _Reviews in Mathematical Physics_ 27 (Oct, 2015) 1550019–143, [1403.2834].
  * (71) P. Huovinen and P. Petreczky, _QCD equation of state and Hadron Resonance Gas_ , _Nucl. Phys._ A837 (2010) 26–53, [0912.2541].
  * (72) MILC collaboration, A. Bazavov et al., _Nonperturbative QCD simulations with 2+1 flavors of improved staggered quarks_ , _Rev. Mod. Phys._ 82 (2010) 1349–1417, [0903.3598].
  * (73) G. Endrődi, _QCD equation of state at nonzero magnetic fields in the Hadron Resonance Gas model_ , _JHEP_ 1304 (2013) 023, [1301.1307].
  * (74) T. Banks and A. Casher, _Chiral symmetry breaking in confining theories_ , _Nucl. Phys._ B169 (1980) 103–125.
  * (75) F. Bruckmann, G. Endrődi, M. Giordano, S. D. Katz, T. G. Kovács, F. Pittler et al., _Landau levels in QCD_ , _Phys. Rev._ D96 (2017) 074506, [1705.10210].
  * (76) E. Elizalde, S. D. Odintsov, A. Romeo, A. A. Bytsenko and S. Zerbini, _Zeta regularization techniques with applications_.  World Scientific, Singapore, 1994.
  * (77) E. S. Fraga, _Thermal chiral and deconfining transitions in the presence of a magnetic background_ , _Lect. Notes Phys._ 871 (2013) 121–141, [1208.0917].
  * (78) N. Nielsen, _Asymptotic freedom as a spin effect_ , _Am. J. Phys._ 49 (1981) 1171.
  * (79) D. J. Toms, _The Effective action at finite temperature and density with application to Bose-Einstein condensation_ , cond-mat/9612003.
  * (80) “NIST Digital Library of Mathematical Functions.” http://dlmf.nist.gov/, Release 1.0.23 of 2019-06-15. For the duplication formula, see https://dlmf.nist.gov/5.5.E5.
  * (81) D. Cangemi and G. V. Dunne, _Temperature expansions for magnetic systems_ , _Annals Phys._ 249 (1996) 582–602, [hep-th/9601048].
  * (82) M. Loewe and J. C. Rojas, _Thermal effects and the effective action of quantum electrodynamics_ , _Phys. Rev._ D46 (1992) 2689–2694.
  * (83) P. Elmfors, D. Persson and B.-S. Skagerstam, _Real time thermal propagators and the QED effective action for an external magnetic field_ , _Astropart.Phys._ 2 (1994) 299–326, [hep-ph/9312226].
  * (84) H. Gies, _Light cone condition for a thermalized QED vacuum_ , _Phys. Rev._ D60 (1999) 105033, [hep-ph/9906303].
  * (85) W.-J. Lee and S. R. Sharpe, _Partial flavor symmetry restoration for chiral staggered fermions_ , _Phys. Rev._ D60 (1999) 114503, [hep-lat/9905023].
  * (86) M. Göckeler et al., _Perturbative and nonperturbative renormalization in lattice QCD_ , _Phys. Rev._ D82 (2010) 114511, [1003.5756].
  * (87) A. T. Lytle and S. R. Sharpe, _Nonperturbative renormalization for improved staggered bilinears_ , _Phys. Rev._ D88 (2013) 054506, [1306.3881].
  * (88) S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabó, _Full result for the QCD equation of state with 2+1 flavors_ , _Phys. Lett._ B730 (2014) 99–104, [1309.5258].