# Deep inelastic scattering on the quark-gluon plasma

Marco Cè  Theoretical Physics Department, CERN, CH-1211 Geneva 23, Switzerland
Tim Harris  Dipartimento di Fisica, Università di Milano-Bicocca, and INFN,
sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy
Harvey B. Meyer  Helmholtz Institut Mainz, Johannes Gutenberg-Universität
Mainz, D-55099 Mainz, Germany  PRISMA+ Cluster of Excellence & Institut für
Kernphysik, Johannes Gutenberg-Universität Mainz, D-55099 Mainz, Germany
Arianna Toniato  PRISMA+ Cluster of Excellence & Institut für Kernphysik,
Johannes Gutenberg-Universität Mainz, D-55099 Mainz, Germany

###### Abstract

We provide an interpretation of the structure functions of a thermal medium
such as the quark-gluon plasma in terms of the scattering of an incoming
electron on the medium via the exchange of a spacelike photon. We then focus
on the deep-inelastic scattering (DIS) regime, and formulate the corresponding
moment sum rules obeyed by the structure functions. Accordingly, these moments
are given by the thermal expectation value of twist-two operators, which is
computable from first principles in lattice QCD for the first few moments. We
also show how lattice QCD calculations can be used to probe how large the
photon virtuality needs to be in order for the Bjorken scaling of structure
functions to set in. Finally, we provide the parton-model interpretation of
the structure functions in the Bjorken limit and test its consistency. As in
DIS on the proton, the kinematic variable $x$ is proportional to the
longitudinal momentum carried by the partons, however $x$ ranges from zero to
infinity. Choosing the parton momentum parametrization to be $xTu$ where $u$
is the fluid four-velocity and $T$ its temperature in the rest frame, the
parton distribution function for a plasma of non-interacting quarks is
proportional to $x\log(1+e^{-x/2})$.

††preprint: MITP/20-075, CERN-TH-2020-206

##  I Introduction

Strong-interaction matter at high temperature is investigated extensively in
high-energy heavy-ion collisions Braun-Munzinger et al. (2016); Busza et al.
(2018). The time-evolution of the matter produced in these collisions has been
described successfully using relativistic hydrodynamics Ollitrault (2008);
Teaney (2010). Although the reaction occurs extremely rapidly, the system
appears to reach an approximate local thermal equilibrium. The typical
temperatures achieved are around 300 MeV. The value of the shear viscosity to
entropy density ratio (around 0.20 Everett et al. (2020)) extracted from
comparing observed soft particle spectra with the outcome of the hydrodynamics
calculations is the lowest of all known forms of matter, indicating that the
produced matter behaves like an excellent fluid.

The high-temperature phase of QCD is known as the quark-gluon plasma (QGP). At
extremely high temperatures, due to asymptotic freedom, one indeed expects the
thermodynamic and transport properties of the medium to be dominated by
quasiparticles carrying the quantum numbers of either quarks or gluons.
Whether such a simple picture applies even approximately at the temperatures
reached in heavy-ion collisions is questionable in view of the
phenomenological findings described above. If one then accepts the picture of
a strongly correlated fluid for the QGP at $T=300\,$MeV, a new question
arises. Probed on short enough distance and time scales, surely this fluid
must exhibit quark and gluon scatterers; this property is guaranteed by
asymptotic freedom and is analogous to the statement that quarks and gluons
can be ‘seen’ inside the proton when probed in deep inelastic scattering. How
then do the fluid properties of the QGP emerge from these effective degrees of
freedom seen on short time scales? This question provides part of the
motivation behind the present work. It has recently been addressed D’Eramo et
al. (2019) by invoking the large-angle scattering of jets or parts thereof on
the QGP in a heavy-ion collision. Here we invoke the idealised situation of a
light lepton scattering on the QGP via the exchange of a spacelike photon, and
ask the question in reverse: starting from the fluid in thermal equilibrium,
at what resolution scale do the quarks become ‘visible’? Although it appears
unlikely that such scattering could be identified in a heavy-ion collision,
the interaction of a lepton with the QGP is conceptually simpler, and it turns
out that it can be probed in lattice QCD via fixed-virtuality dispersion
relations. At large photon virtuality $Q^{2}$, the cross-section per unit
volume can be understood in the framework of deep inelastic scattering (DIS).
The structure functions of the medium then depend only logarithmically on
$Q^{2}$, and a parton distribution function of the QGP can be defined.

There is a further motivation for studying the structure functions of the QGP.
The latter are connected via a Kubo-Martin-Schwinger relation to the spectral
functions through which the thermal production rate of photons and dileptons
are usually expressed. These rates are of central importance for understanding
the corresponding spectra measured in heavy-ion collisions Braun-Munzinger et
al. (2016); David (2020). As far as the photon emissivity is concerned, it has
recently been shown that it can be probed in lattice QCD at fixed, vanishing
virtuality $Q^{2}=0$ Meyer (2018). Probing the dilepton rate, however, is only
possible via a dispersion relation at fixed spatial momentum $\boldsymbol{q}$
of the lepton pair; see for instance references Ding et al. (2011); Brandt et
al. (2013); Amato et al. (2013); Aarts et al. (2015); Brandt et al. (2016);
Ding et al. (2016). In these calculations, the Euclidean correlator probes
both the spacelike regime of the spectral function and the timelike regime,
since it is represented as an integral over all positive real frequencies
$q^{0}$. For that reason, it turns out that understanding the behaviour of the
spectral functions in the DIS regime can be helpful in probing the low-mass
dilepton production rate with a large spatial momentum relative to the plasma
in lattice QCD. We will return to this point at the end of section II.

We note that certain aspects of DIS on non-Abelian plasmas have been addressed
previously in the references Hatta et al. (2008); Iancu and Mueller (2009). In
the first one, DIS on an ${\cal N}=4$ super Yang-Mills plasma is investigated
using the AdS/CFT correspondence, and the possibility of interpreting the
corresponding structure functions in terms of partons is discussed in detail.
The impact-parameter representation plays an important role in the analysis.
The second reference proposes to evaluate the thermal expectation value of a
particular twist-two operator in QCD, and more concretely in quenched QCD, in
order to judge whether its size is more typical of a weak-coupling or a
strong-coupling ‘scenario’. Its relation to the medium’s structure functions
via a moment sum rule is used to infer a partonic interpretation of these two
scenarios. With respect to these references, our analysis differs by its
emphasis on full QCD and by technical aspects; in particular the range of the
Bjorken-$x$ variable is different, the fixed-virtuality dispersion relations
are new and the connection to the spectral functions discussed in the context
of photon and dilepton production plays a major role. This paper also has
common aspects with the work of Ref. Chambers et al. (2017), which studies the
nucleon structure functions on the lattice via fixed-virtuality dispersion
relations. A difference with the nucleon structure functions and, more
generally, with zero-temperature structure functions, is that in the thermal
case one cannot simply take derivatives of Euclidean correlators with respect
to the variable $q^{0}$ without confronting a numerically ill-posed analytic
continuation problem. This difference stems from the fact that the thermal
correlator can be computed only at the discrete Matsubara frequencies,
$q^{0}=i2\pi Tn$, $n\in\mathbb{Z}$.

The rest of this paper is structured as follows. Section II presents the
derivation of the cross-section per unit volume for a lepton scattering on a
thermal medium in terms of its structure functions, and provides their
connection to the spectral functions. The Bjorken limit corresponding to deep-
inelastic scattering is formulated and the structure functions for a plasma of
non-interacting quarks are given in that limit. In section III, the moment sum
rules are presented for the medium structure functions, allowing the first few
moments to be computed non-perturbatively in lattice QCD. Particular attention
is devoted to the $n=2$ sum rule, corresponding to the momentum sum rule in
the parton model. Opportunities for lattice QCD studies of the medium
structure functions for spacelike momenta are presented in section IV.
Finally, the parton-model interpretation of the structure functions in the DIS
limit is discussed in section V. We summarize and conclude in section VI. A
number of technical aspects are collected in appendix, in particular the
derivation of the moment sum rules and the explicit verification of the first
two sum rules for the plasma of non-interacting quarks.

##  II Lepton scattering on a thermal medium and structure functions of the
latter

In this section, we recall some of the basic aspects of the inelastic
scattering of a lepton on hadronic matter. Our notation is largely taken from
the lecture notes Manohar (1992). Let $k$ be the initial momentum of the
lepton, and $k^{\prime}$ its final momentum. We also use $E$ and $E^{\prime}$
respectively to denote the zeroth components of these vectors. The momentum
transferred to the medium is thus $q=k-k^{\prime}$, and we use
$Q^{2}=-q^{2}\geq 0$ to denote the spacelike virtuality of the exchanged
photon. The lepton mass is neglected throughout.

QGP$e^{-}$$e^{-}$$k$$k^{\prime}$$\gamma$$q$$u$ $X$

$e^{-}$$e^{-}$${\rm q}$${\rm
q}$$k$$k^{\prime}$$p_{\xi}$$p_{\xi}^{\prime}$$\gamma$$q$

Figure 1: Left: Scattering of a lepton on quark-gluon plasma at global thermal
equilibrium moving with a four-velocity $u$ in the one-photon exchange
approximation. The picture shows the interaction occurring with a cubic fluid
cell, producing an unobserved QCD final-state $X$. Right: interpretation of
the process in the deep-inelastic scattering limit as an elastic parton-lepton
collision.

###  II.1 Inelastic cross section and hadronic tensor

First, recall the expression for the inclusive differential cross section for
lepton-proton scattering,

$\displaystyle d\sigma$ $\displaystyle=$
$\displaystyle\sum_{X}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}2E^{\prime}}\;\frac{(2\pi)^{4}\delta^{4}(k+p-k^{\prime}-p_{X})}{(2E)(2M)}\,\frac{e^{4}}{Q^{4}}\,$
$\displaystyle\qquad\qquad\qquad\qquad\times\langle
p,\lambda|j^{\mu}(0)|X\rangle\;\langle
X|j^{\nu}(0)|p,\lambda\rangle\;\ell_{\mu\nu}\;,\phantom{\frac{1}{1}}$
$\displaystyle\ell_{\mu\nu}$ $\displaystyle=$
$\displaystyle\sum_{s_{l}^{\prime}}\langle
k,s_{l}|j_{\ell\mu}(0)|k^{\prime},s_{l}^{\prime}\rangle\,\langle
k^{\prime},s_{l}^{\prime}|j_{\ell\nu}(0)|k,s_{l}\rangle\;,$ (2)

where $M$ is the proton mass and the proton-state normalisation is $\langle
p^{\prime}\lambda^{\prime}|p\lambda\rangle=2E_{\boldsymbol{p}}(2\pi)^{3}\delta(\boldsymbol{p}-\boldsymbol{p}^{\prime})\delta_{\lambda\lambda^{\prime}}$.
The leptonic tensor $\ell_{\mu\nu}$ is easily computed, and its expression can
be found in Manohar (1992). With

$j^{\mu}(x)=\sum_{f=u,d,s,\dots}{\cal
Q}_{f}\;\bar{\psi}_{f}(x)\gamma^{\mu}\psi_{f}(x),$ (3)

the electromagnetic current carried by the quarks111The physical values of the
quark charges are ${\cal Q}_{u}={\textstyle\frac{2}{3}},\leavevmode\nobreak\
{\cal Q}_{d}={\cal Q}_{s}=-{\textstyle\frac{1}{3}},\leavevmode\nobreak\
\dots$, one definition of the hadronic tensor is given by

$\displaystyle W_{>}^{\mu\nu}(p,q)_{\lambda^{\prime}\lambda}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\int d^{4}x\;e^{iq\cdot x}\,\langle
p,\lambda^{\prime}|j^{\mu}(x)\,j^{\nu}(0)|p,\lambda\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\sum_{X}(2\pi)^{4}\delta^{4}(q+p-p_{X})\;\langle
p,\lambda^{\prime}|j^{\mu}(0)|X\rangle\;\langle
X|j^{\nu}(0)|p,\lambda\rangle,$

which allows one to substitute the sum over final states $X$ by
$W_{>}^{\mu\nu}(p,q)$. In fact, the hadronic tensor $W^{\mu\nu}(p,q)$ is
usually defined by replacing the product of currents by their commutator. The
tensors $W_{>}^{\mu\nu}$ and $W^{\mu\nu}$ are equal for the kinematics
accessible in $ep$ scattering. One then obtains for the differential cross
section

$\frac{d^{2}\sigma}{dE^{\prime}d\Omega}=\frac{e^{4}}{8\pi^{2}Q^{4}}\left(\frac{E^{\prime}}{ME}\right)\ell_{\mu\nu}\;\frac{1}{2}\sum_{\lambda}W^{\mu\nu}(p,q)_{\lambda\lambda}.$
(5)

The derivation of the cross-section for a lepton scattering on thermal
hadronic matter is very similar to the case of scattering on a proton. See the
left panel of figure 1 illustrating the process. If $L^{3}$ is the volume of
the thermal system in its rest frame, the differential cross section reads

$\frac{d^{2}\sigma}{dE^{\prime}d\Omega}=\frac{e^{4}\,L^{3}}{8\pi^{2}Q^{4}}\left(\frac{E^{\prime}}{E}\right)\ell_{\mu\nu}W_{>}^{\mu\nu}(u,q),$
(6)

where the hadronic tensor is given by

$W_{>}^{\mu\nu}(u,q)=\frac{1}{4\pi\,Z}\sum_{n}e^{-\beta E_{n}}\int
d^{4}x\;e^{iq\cdot x}\,\langle n|j^{\mu}(x)\,j^{\nu}(0)|n\rangle,$ (7)

$Z=\sum_{n}e^{-\beta E_{n}}$ being the canonical partition function. While the
temperature dependence of $W_{>}^{\mu\nu}$ is not indicated explicitly, the
argument $u^{\mu}$ represents the four-velocity of the fluid. The state
$|n\rangle$ is unit-normalized. In inelastic scattering on the proton,
$q^{0}=E-E^{\prime}\geq 0$, since the energy of the final hadronic state can
only be larger than the proton mass; in scattering on the plasma, it is
possible for the electron to gain energy ($q^{0}<0$), although this is
suppressed as $e^{-\beta|q^{0}|}$. Here, due to our normalisation of the
$|n\rangle$ states, the $(2M)^{-1}$ factor appearing in Eq. (II.1) should not
be included; instead a factor $L^{3}$ should be included. We recall that for a
thermal state, the typical energy of a contributing state $|n\rangle$ is of
order volume, while its spatial momentum is of order $L^{3/2}$ in the
canonical ensemble (see e.g.Giusti and Meyer (2011)). Therefore, working in
the thermodynamic limit and in the rest frame of the fluid, a typical state
$|n\rangle$ can be treated as if it were at rest.

Through the Kubo-Martin-Schwinger relation, $W_{>}^{\mu\nu}$ can be rewritten
in terms of the spectral function $\rho^{\mu\nu}$, which is defined as the
Fourier transform of the commutator of the two currents, analogously to
$W^{\mu\nu}$ in the proton case (see Eq. 58 for the explicit definition).
Using Eq. (38) of Ref. Meyer (2011)222Our present convention for the
normalization of the spectral function is that no factor (of $1/(4\pi)$ or
other) appears in front of the Fourier transform of the current commutator,
which makes the spectral function defined here a factor $(2\pi)$ larger than
in the convention used in Meyer (2011).,

$W_{>}^{\mu\nu}(u,q)=\frac{1}{4\pi(1-e^{-\beta
q^{0}})}\,\rho^{\mu\nu}(q^{0},\boldsymbol{q})\qquad\textrm{(fluid rest
frame).}$ (8)

It is worth pointing out here that in the DIS regime defined below, since
$q^{0}$ is extremely large, the structure function only differs from the
spectral function by the purely conventional factor of $4\pi$. One may
parametrize the tensor $W_{>}^{\mu\nu}$ in a way similar to the proton case,

$W_{>}^{\mu\nu}(u,q)=F_{1}(u\cdot
q,Q^{2})\Big{(}-g^{\mu\nu}+\frac{q^{\mu}q^{\nu}}{q^{2}}\Big{)}+\frac{m_{T}}{u\cdot
q}F_{2}(u\cdot q,Q^{2})\Big{(}u^{\mu}-(u\cdot
q)\frac{q^{\mu}}{q^{2}}\Big{)}\Big{(}u^{\nu}-(u\cdot
q)\frac{q^{\nu}}{q^{2}}\Big{)}.$ (9)

Here $m_{T}$ is some thermal energy scale; we will later choose to set
$m_{T}=T$. The contraction of $W_{>}^{\mu\nu}$ with the symmetric part of the
leptonic tensor

$\ell^{\\{\mu\nu\\}}=2(k^{\mu}k^{\prime}{}^{\nu}+k^{\nu}k^{\prime}{}^{\mu}-g^{\mu\nu}(k\cdot
k^{\prime}))$ (10)

yields

$\ell_{\\{\mu\nu\\}}W_{>}^{\mu\nu}=\frac{4m_{T}(u\cdot
k)}{y}\,\Big{(}xy^{2}F_{1}+(1-y)F_{2}\Big{)}\,,$ (11)

where

$x=\frac{Q^{2}}{2m_{T}(u\cdot q)},\qquad\qquad y=\frac{u\cdot q}{u\cdot k}\leq
1\,.$ (12)

The last inequality comes from the requirement that the lepton cannot lose
more than its initial energy $E$. We note that the definition of $x$ is
analogous to the definition of Bjorken-$x$ in DIS on the proton. One
qualitative difference between scattering on the proton and on the fluid is
that in the latter case, the value $x=1$ does not have any special
significance, while for the proton it corresponds to elastic scattering. Thus
the multiplicative normalization of $x$ has no particular meaning, which is
reflected in the freedom of choosing the precise value of $m_{T}$. For
convenience, we may introduce

$\rho_{1,2}(q^{0},Q^{2})=4\pi(1-e^{-\beta
q^{0}})F_{1,2}(q^{0},Q^{2}),\qquad\textrm{(fluid rest frame)}$ (13)

so that the cross-section may finally be written in terms of the spectral
functions,

$\frac{d^{2}\sigma}{L^{3}dE^{\prime}d\Omega}=\frac{2\alpha^{2}\,E^{\prime}m_{T}}{\pi
Q^{4}y(1-e^{-\beta
q^{0}})}\,\Big{(}xy^{2}\rho_{1}(q^{0},Q^{2})+(1-y)\rho_{2}(q^{0},Q^{2})\Big{)},\qquad\textrm{(fluid
rest frame)}.$ (14)

One easily verifies that $\rho_{1}/q^{0}$ and $\rho_{2}$ are even, non-
negative functions of $q^{0}$ at fixed $|\boldsymbol{q}|$ Meyer (2011), and
hence at fixed $Q^{2}\geq 0$. Therefore the cross-section (14) is guaranteed
to be positive for all values of $y\leq 1$. We recall that the DIS limit
corresponds to $Q^{2}\to+\infty$ with the variable $x$ fixed. In this limit,
due to the asymptotic freedom property of QCD, the structure functions
$F_{1,2}$ and the corresponding spectral functions $\rho_{1,2}$ are expected
to depend mainly on $x$, and only logarithmically on $Q^{2}$. In the parton-
model interpretation of the process, the lepton scattering cross section
arises from the incoherent sum of lepton-parton scatterings; see the right
panel of figure 1. Given that in QCD these partons are massless spin-1/2
particles, this picture implies the Callan-Gross relation, $F_{2}=2xF_{1}$,
which receives corrections logarithmically suppressed in $Q^{2}$.

We set $m_{T}=T$ from here on until section (V), in which the role of this
thermal energy scale is discussed within the parton model.

###  II.2 Structure functions of the non-interacting plasma

Working in the rest frame of the medium, we note that spectral functions are
usually calculated as a function of $q^{0}$ (also often notated $\omega$) and
$\boldsymbol{q}$. We remark that if one chooses $|\boldsymbol{q}|$ and $x$ as
independent variables, the frequency variable takes the form

$q^{0}=\sqrt{\boldsymbol{q}^{2}+T^{2}x^{2}}\;-Tx\stackrel{{\scriptstyle|\boldsymbol{q}|\to\infty}}{{\simeq}}|\boldsymbol{q}|-Tx.$
(15)

The DIS regime can be reached by sending $|\boldsymbol{q}|\to\infty$ at fixed
$x$, since $Q^{2}=2T|\boldsymbol{q}|x+{\rm O}(T^{2})$.

Thermal QGP spectral functions at leading order in perturbation theory are
given for instance in Laine (2013). Taking the leading-order expression of
$\rho_{1}$, which is equal to $\rho^{11}$ in the fluid rest frame for
$q^{\mu}=(q^{0},0,0,q^{3})$, one finds, for $N_{c}$ colors of massless quarks,

$\lim_{Q^{2}\to\infty}F_{1}(x,Q^{2})=\frac{(\sum_{f}\mathcal{Q}_{f}^{2})N_{c}T^{2}}{4\pi^{2}}x\log(1+e^{-x/2}).$
(16)

One also finds that $-\frac{1}{2}W_{>}^{\mu}{}_{\mu}$ has the same limit. It
follows that

$\lim_{Q^{2}\to\infty}(F_{2}-2xF_{1})=0,$ (17)

which is the Callan-Gross relation. Thus the leading-order calculation of the
spectral functions confirms that the $Q^{2}$ dependence disappears in the
Bjorken limit, and that the Callan-Gross relation is satisfied.

The graph of expression (16) is displayed in Fig. 2. It is positive and
normalizable as a function of $x$. The approach of $F_{1}$ to the expression
in Eq. (16) is also illustrated the figure. Clearly, $F_{1}$ approaches its
limit sooner at small $x$ than at large $x$. It is not before $Q/T=10$ that
$F_{1}$ qualitatively acquires its asymptotic shape.

We now briefly return to a remark made in the introduction. As noted around
Eq. (15), the DIS regime can be reached by setting $q^{0}=|\boldsymbol{q}|-xT$
for a fixed value of the Bjorken variable $x$, and sending
$|\boldsymbol{q}|\to\infty$. Thus, in a dispersive representation at large and
fixed $|\boldsymbol{q}|$, the Euclidean correlator receives a contribution
from the spectral function in the DIS regime. More precisely, if $Q_{0}^{2}$
denotes the virtuality beyond which the DIS regime applies for $x$ of order
unity, the range of frequencies corresponding to the DIS regime is

$|\boldsymbol{q}|-\pi T\lesssim
q^{0}\leq|\boldsymbol{q}|-\frac{Q_{0}^{2}}{2|\boldsymbol{q}|}.$ (18)

Therefore, understanding the behaviour of the spectral functions in the DIS
regime is helpful for probing the low-mass dilepton production rate with a
large spatial momentum relative to the plasma using lattice QCD correlators in
conjunction with a fixed-$|\boldsymbol{q}|$ dispersion relation.

Figure 2: The structure function $F_{1}(x,Q^{2})$ of the free plasma for
different values of $Q^{2}$, together with the Bjorken limit Eq. (16). Here we
have set $N_{c}=3$ and $\sum_{f}{\cal Q}_{f}^{2}=1$.

##  III Sum rules

A standard tool to analyze the structure functions in the DIS regime is the
operator-product expansion (OPE). Here we review the OPE analysis to formulate
the moment sum rules of the structure functions. We set $N_{c}=3$ in this
section, as appropriate for QCD.

The time-ordered operator-product of two electromagnetic current admits the
asymptotic expansion

$\begin{split}t_{\mu\nu}(q)&\equiv\>i\int\mathrm{d}^{4}x\>e^{iq\cdot x}\>{\rm
T}\\{j_{\mu}(x)j_{\nu}(0)\\}\\\
&\overset{Q\to\infty}{\sim}\sum_{n=2,4,\dots}\sum_{f,j}c_{f;\mu\nu\mu_{1}\dots\mu_{n}}(q)\>M_{fj}(Q,\tilde{\mu})\>\bigl{[}O_{nj}^{\mu_{1}\dots\mu_{n}}\bigr{]}_{\tilde{\mu}}\>,\end{split}$
(19)

where $\tilde{\mu}$ is a reference energy scale at which the local operators
on the right-hand side are renormalized. The twist-two fermionic operator with
flavor $f$ is

$O_{nf}^{\mu_{1}\dots\mu_{n}}=\frac{1}{2}\biggl{(}\frac{i}{2}\biggr{)}^{n-1}\mathcal{S}\\{\bar{\psi}_{f}\gamma^{\mu_{1}}\overleftrightarrow{D}^{\mu_{2}}\dots\overleftrightarrow{D}^{\mu_{n}}\psi_{f}\\}\>,$
(20)

and the gluonic one reads

$O_{ng}^{\mu_{1}\dots\mu_{n}}=-\frac{1}{2}\biggl{(}\frac{i}{2}\biggr{)}^{n-2}\mathcal{S}\\{F^{\mu_{1}\alpha}\overleftrightarrow{D}^{\mu_{2}}\dots\overleftrightarrow{D}^{\mu_{n-1}}F^{\mu_{n}}_{\quad\alpha}\\}\>.$
(21)

Here $\overleftrightarrow{D}=\overrightarrow{D}-\overleftarrow{D}$ and ${\cal
S}$ acts as $1/n!$ times the sum over permutations of the $n$ Lorentz indices
and the subtraction of trace terms, such that contracting any pair of indices
yields zero. The coefficients $c_{f;\mu\nu\mu_{1}\dots\mu_{n}}(q)$ are given
by

$\begin{split}c_{f;\mu\nu\mu_{1}\dots\mu_{n}}(q)=&2\mathcal{Q}_{f}^{2}\bigl{(}-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{q^{2}}\bigr{)}2^{n}\frac{q_{\mu_{1}}\dots
q_{\mu_{n}}}{(Q^{2})^{n}}\>+\\\
&2\mathcal{Q}_{f}^{2}\bigl{(}g_{\mu\mu_{1}}-\frac{q_{\mu}q_{\mu_{1}}}{q^{2}}\bigr{)}\bigl{(}g_{\nu\mu_{2}}-\frac{q_{\nu}q_{\mu_{2}}}{q^{2}}\bigr{)}2^{n}\frac{q_{\mu_{3}}\dots
q_{\mu_{n}}}{(Q^{2})^{n-1}}\>\end{split}$ (22)

and $M_{fj}(Q,\tilde{\mu})$ are the coefficients of the operator mixing. To
leading order in perturbation theory, $M_{fj}^{LO}=\delta_{fj}$, and to next-
to-leading order

$M_{fj}^{NLO}(Q,\tilde{\mu})=\biggl{[}\biggl{(}\frac{\log(Q^{2}/\Lambda^{2})}{\log(\tilde{\mu}^{2}/\Lambda^{2})}\biggr{)}^{a^{(n)}/2b_{0}}\biggr{]}_{fj}\>.$
(23)

In this expression, $\Lambda$ is the energy at which the one-loop renormalized
coupling diverges, $b_{0}=11-\frac{2}{3}n_{f}$ the first coefficient of the
expansion of the beta function in powers of the gauge coupling $g$ and
$\gamma_{ij}^{(n)}=-\frac{g^{2}}{(4\pi)^{2}}a_{ij}^{(n)}$ is the anomalous-
dimension matrix for the twist-two operators of dimension $n+2$. The index $j$
runs over $\\{f,g\\}$.

The strategy we follow to obtain moment sum rules for the structure functions
from the OPE is to first relate the time-ordered product to the corresponding
retarded commutator of currents; and secondly to use the fixed-virtuality
dispersive representation of the retarded correlator. Details of the
derivation are given in appendix A. The result is

$\int_{0}^{\infty}\mathrm{d}x\>x^{n-1}\>[F_{1}(x,Q^{2})]_{\textrm{leading-
twist}}=\sum_{f,j}\frac{\mathcal{Q}_{f}^{2}}{2}M_{fj}(Q,\tilde{\mu})\langle
O_{nj}\rangle\>,\quad n=2,4,\dots\>,$ (24)

where

$\langle O_{nj}^{\mu_{1}\dots\mu_{n}}\rangle=T^{n}[u^{\mu_{1}}\dots
u^{\mu_{n}}-\mathrm{traces}]\langle O_{nj}\rangle\>.$ (25)

Apart from the upper integration limit being infinity rather than unity, which
makes it necessary to isolate the leading-twist part of the structure function
prior to computing its moments, the form of the moment sum rule (24) is the
same as in standard DIS. In appendix B, we verify these moment sum rules at
leading order in the case of non-interacting quarks by separately computing
the left- and right-hand side. Given that $\langle O_{2j}\rangle$ is
proportional to the enthalpy, that the total momentum of the fluid is
generically given by the product of its enthalpy with its four-velocity, and
that $F_{1}$ has a parton-model interpretation as being proportional to the
parton distribution with a known prefactor (Eq. (49) below), the $n=2$ sum
rule expresses in this case the (common sense) fact that the quarks carry the
entire momentum of the fluid. Next, we formulate specifically the $n=2$ sum
rule for Wilson coefficients at next-to-leading order (NLO) accuracy, which is
appropriate for QCD and significantly modifies the OPE prediction for the
$n=2$ moment of the structure function $F_{1}$.

###  III.1 The $n=2$ sum rule to NLO

We now formulate the lowest moment sum rule at NLO accuracy in the Wilson
coefficients. We follow the notation of Ref. Peskin and Schroeder (1995). The
twist-two fermionic local operator renormalized at a scale $Q^{2}$ is, at NLO,

$\begin{split}\bigl{[}O_{2f}^{\mu\nu}\bigr{]}_{Q}&=\frac{1}{2}\frac{1}{16/3+n_{f}}T^{\mu\nu}_{\mathrm{traceless}}\\\
&+\frac{1}{n_{f}(16/3+n_{f})}\biggl{(}\frac{\log(Q^{2}/\Lambda^{2})}{\log(\tilde{\mu}^{2}/\Lambda^{2})}\biggr{)}^{-\frac{4}{3}(\frac{16}{3}+n_{f})/2b_{0}}\biggl{[}\frac{16}{3}\sum_{f^{\prime}}O_{2f^{\prime}}^{\mu\nu}-n_{f}O_{2g}^{\mu\nu}\biggr{]}_{\tilde{\mu}}\\\
&+\biggl{(}\frac{\log(Q^{2}/\Lambda^{2})}{\log(\tilde{\mu}^{2}/\Lambda^{2})}\biggr{)}^{-32/9b_{0}}\biggl{[}O_{2f}^{\mu\nu}-\frac{1}{n_{f}}\sum_{f^{\prime}}O_{2f^{\prime}}^{\mu\nu}\biggr{]}_{\tilde{\mu}}\;,\end{split}$
(26)

where $n_{f}$ is the number of fermion flavors, and
$T^{\mu\nu}_{\mathrm{traceless}}=2\bigl{(}\sum_{f}O_{2f}^{\mu\nu}+O_{2g}^{\mu\nu}\bigr{)}$
is the traceless part of the energy-momentum tensor. In the absence of a
flavor-non-singlet chemical potential, the in-medium expectation value of the
last term vanishes by SU($n_{f}$) flavor symmetry. As a consequence, the
moment sum rule to next-to-leading order in $\alpha_{s}$ is

$\begin{split}&\int_{0}^{\infty}\mathrm{d}x\>x\,[F_{1}(x,Q^{2})]_{\textrm{leading-
twist}}=\sum_{f}\frac{\mathcal{Q}^{2}_{f}}{2}\biggl{[}\frac{1}{2}\frac{1}{16/3+n_{f}}\langle
T\rangle\\\
&\qquad\quad+\frac{1}{n_{f}(16/3+n_{f})}\biggl{(}\frac{\log(Q^{2}/\Lambda^{2})}{\log(\tilde{\mu}^{2}/\Lambda^{2})}\biggr{)}^{-\frac{4}{3}(\frac{16}{3}+n_{f})/2b_{0}}\bigl{(}\frac{16}{3}\sum_{f^{\prime}}\langle
O_{2f^{\prime}}\rangle-n_{f}\langle
O_{2g}\rangle\bigr{)}\biggr{]}\>,\end{split}$ (27)

where

$\begin{split}&\langle
O_{2j}^{\mu\nu}\rangle=T^{2}\biggl{[}u^{\mu}u^{\nu}-\frac{1}{4}g^{\mu\nu}\biggr{]}\langle
O_{2j}\rangle\,,\\\ &\langle
T^{\mu\nu}_{\mathrm{traceless}}\rangle=T^{2}\biggl{[}u^{\mu}u^{\nu}-\frac{1}{4}g^{\mu\nu}\biggr{]}\langle
T\rangle\>.\end{split}$ (28)

Given the energy-momentum tensor of a relativistic fluid at thermal
equilibrium

$T^{\mu\nu}=(e+p)u^{\mu}u^{\nu}-pg^{\mu\nu}\>,$ (29)

we find

$\langle T\rangle=\frac{e+p}{T^{2}}\>,$ (30)

where $e$ is the energy density and $p$ the pressure, their sum corresponding
to the enthalpy density. In the extreme $Q\to\infty$ limit, we have

$\lim_{Q^{2}\to\infty}\int_{0}^{\infty}\mathrm{d}x\>x\>[F_{1}(x,Q^{2})]_{\textrm{leading-
twist}}=\bigl{(}\sum_{f}{\mathcal{Q}_{f}^{2}}\bigr{)}\frac{1}{16/3+n_{f}}\frac{e+p}{4T^{2}}\>.$
(31)

Thus the asymptotic $n=2$ moment of the structure function $F_{1}$ is uniquely
determined by the enthalpy density of the medium. By invoking the
interpretation of $F_{1}$ in terms of a parton distribution function (see Eq.
(49) below), and recalling that the total momentum density carried by the
fluid is $(e+p)\boldsymbol{u}$, Eq. (31) also shows that the asymptotic
momentum fraction carried by one out of $n_{f}$ flavors of quarks is
$1/(n_{f}+16/3)$, while the gluons carry the fraction $(16/3)/(n_{f}+16/3)$.
These fractions are exactly the same as in conventional DIS on the nucleon. In
that context, these momentum fractions only apply at extremely large $Q^{2}$.
For $n_{f}=3$, the asymptotic gluon momentum fraction is 0.64. This number is
substantially larger than the enthalpy fraction contributed by the gluons in
the high-temperature limit, which is about 0.34. This difference reflects the
fact that the non-interacting structure function $F_{1}$ satisfies the $n=2$
sum rule with leading-order Wilson coefficients, and not the NLO sum rule (31)
appropriate for QCD. Thus the $n=2$ moment of the structure function $F_{1}$
in the DIS limit is unstable against ‘turning on’ interactions between quarks,
and cannot be expanded in powers of $g^{2}(T)$. We return to this point in the
conclusion.

##  IV Tests of the onset of Bjorken scaling from Euclidean correlation
functions

An interesting question is how large the photon virtuality has to be in order
for the structure functions at fixed $x$ to become $Q^{2}$-independent. It is
tempting to address this question by computing both sides of a moment sum
rule, however the $x$-moments of the spectral functions correspond to
derivatives of Euclidean correlators at vanishing Matsubara frequency, and
these derivatives are only accessible via a numerically ill-posed analytic
continuation. We therefore propose tests that do not require a numerical
analytic continuation.

In the rest frame of the plasma, the spectral function components are related
to the longitudinal and transverse spectral functions $\sigma^{L}$ and
$\sigma^{T}$, as defined in Ref. Meyer (2018), as well as to the spectral
functions $\rho_{1}$ and $\rho_{2}$ introduced in the DIS context according to

$\begin{split}\frac{1}{2}\bigl{(}\delta^{ij}-\frac{q^{i}q^{j}}{\boldsymbol{q}^{2}}\bigr{)}\rho^{ij}(u,q)=-\sigma^{T}(q^{0},Q^{2})&=\rho_{1}(q^{0},Q^{2})\;,\\\
\rho^{00}(u,q)=\frac{\boldsymbol{q}^{2}}{Q^{2}}\sigma^{L}(q^{0},Q^{2})&=-\frac{\boldsymbol{q}^{2}}{Q^{2}}\rho_{1}(q^{0},Q^{2})+\frac{T|\boldsymbol{q}|^{4}}{q^{0}Q^{4}}\rho_{2}(q^{0},Q^{2})\;,\\\
\frac{q^{i}q^{j}}{\boldsymbol{q}^{2}}\rho^{ij}(u,q)=\frac{(q^{0})^{2}}{Q^{2}}\sigma^{L}(q^{0},Q^{2})&=-\frac{(q^{0})^{2}}{Q^{2}}\rho_{1}(q^{0},Q^{2})+\frac{Tq^{0}\boldsymbol{q}^{2}}{Q^{4}}\rho_{2}(q^{0},Q^{2})\>.\end{split}$
(32)

The two longitudinal and transverse components of the retarded polarisation
tensor are denoted as $G_{R}^{L}(q^{0},\boldsymbol{q}^{2})$ and
$G_{R}^{T}(q^{0},\boldsymbol{q}^{2})$; see Eq. (56) for their explicit
definition. The imaginary part, for $q^{0}$ approaching the real axis from
above, yields the corresponding spectral function. At a spacelike or lightlike
point $Q^{2}\geq 0$, the relation to the spectral functions $\sigma^{L,T}$
reads

$2\,{\rm
Im}\,G_{R}^{L,T}(q^{0},\boldsymbol{q}^{2})=\sigma^{L,T}(q^{0},Q^{2}),\qquad
Q^{2}=\boldsymbol{q}^{2}-(q^{0})^{2}.$ (33)

For imaginary $q^{0}=i\omega_{n}$, with $\omega_{n}=2n\pi T>0$ the $n$-th
Matsubara frequency, the retarded correlator coincides with the Euclidean
correlator Meyer (2011),

$G_{R}^{L,T}\bigl{(}i\omega_{n},\boldsymbol{q}^{2}\bigr{)}=G_{E}^{L,T}\bigl{(}\omega_{n},\boldsymbol{q}^{2}\bigr{)}\equiv
H_{E}^{L,T}(\omega_{n};Q^{2})\>,\qquad
Q^{2}=\boldsymbol{q}^{2}+\omega_{n}^{2}.$ (34)

For a given spacelike virtuality $Q^{2}$, the Euclidean correlators can be
obtained in lattice QCD as a function of Matsubara frequency by computing

$\displaystyle H_{E}^{T}(\omega_{n};Q^{2})$ $\displaystyle=$
$\displaystyle\int
d^{4}x\;e^{\sqrt{\omega_{n}^{2}-Q^{2}}\hat{\boldsymbol{q}}\cdot\boldsymbol{x}+i\omega_{n}x_{0}}\;\Big{(}{\textstyle\frac{1}{2}}\bigl{(}\delta_{ik}-\hat{q}_{i}\hat{q}_{k}\bigr{)}\;\langle
j_{i}(x)\,j_{k}(0)\rangle\Big{)},$ (35) $\displaystyle
H_{E}^{L}(\omega_{n};Q^{2})$ $\displaystyle=$ $\displaystyle\int
d^{4}x\;e^{\sqrt{\omega_{n}^{2}-Q^{2}}\hat{\boldsymbol{q}}\cdot\boldsymbol{x}+i\omega_{n}x_{0}}\;\Big{(}\langle
j_{0}(x)\,j_{0}(0)\rangle+\hat{q}_{i}\hat{q}_{k}\langle
j_{i}(x)\,j_{k}(0)\rangle\Big{)}.$ (36)

On the right-hand side of these equations, $\hat{\boldsymbol{q}}$ is the unit
vector representing the direction of the spatial momentum, and we are using
Euclidean notation, implying a sign change in the two-point functions of the
spatial current relative to its Minkowski-space correspondent. We have
anticipated the fact that in order to probe DIS kinematics ($Q^{2}={\rm
O}(Tq^{0})$), the spatial momentum in the Euclidean correlator needs to become
imaginary and close in magnitude to $\omega_{n}$. Thus the weight function of
the coordinate-space correlator is a real exponential, as for the case of
vanishing virtuality Meyer (2018). The dispersion relation (37) below however
also applies to the case where one or more of the spatial momenta are real, as
long a the virtuality is kept fixed. We also remark that Eqs. (35–36) could be
averaged analytically over the direction of $\hat{\boldsymbol{q}}$ using
Legendre polynomials, along the lines of Meyer (2017). This could prove
numerically advantageous in calculations performed on lattices with a large
spatial volume.

The fixed-virtuality dispersion relation Meyer (2018)

$H_{E}^{L,T}(\omega_{n};Q^{2})-H_{E}^{L,T}(\omega_{r};Q^{2})=\int_{0}^{\infty}\frac{\mathrm{d}\omega}{\pi}\>\omega\>\sigma^{L,T}(\omega,Q^{2})\biggl{[}\frac{1}{\omega^{2}+\omega_{n}^{2}}-\frac{1}{\omega^{2}+\omega_{r}^{2}}\biggr{]}\>$
(37)

follows. The subtraction of $H_{E}^{L,T}(\omega_{r};Q^{2})$ is needed to
ensure the convergence of the dispersive integral333The forward Compton
scattering amplitudes of the proton require analogous subtractions. We also
remark that, in the lattice regularization of QCD, an additional subtraction
is in general needed on the right-hand side of Eq. (37) due to an ultraviolet
divergence; see Meyer (2018) for one such suitable subtraction.. We can
rewrite Eq. 37 as an integral over the variable $x=Q^{2}/(2T\omega)$ as
follows

$H_{E}^{L,T}(\omega_{n};Q^{2})-H_{E}^{L,T}(\omega_{r};Q^{2})=\int_{0}^{\infty}\frac{\mathrm{d}x}{\pi}\>x\>\hat{\sigma}^{L,T}(x,Q^{2})\frac{a_{r}^{2}-a_{n}^{2}}{(1+a_{n}^{2}x^{2})(1+a_{r}^{2}x^{2})}\>,$
(38)

where $a_{n}\equiv 2T\omega_{n}/Q^{2}$ and for convenience we have defined
spectral functions that take as arguments the Bjorken variables,
$\hat{\sigma}^{L,T}(x,Q^{2})=\sigma^{L,T}(Q^{2}/(2Tx),Q^{2})$. The Euclidean
correlators $H_{E}^{L}$ and $H_{E}^{T}$ can be computed directly in lattice
QCD. Thus Eq. (38) can be used to probe the $Q^{2}$ evolution of the spectral
functions $\hat{\sigma}^{L,T}$. For this purpose, the kinematic variables
$a_{n}$ and $a_{r}$ should be kept fixed as $Q^{2}\to\infty$. In particular,
the extent to which the spectral functions become $Q^{2}$-independent at fixed
$x$ can be probed via the corresponding Euclidean correlators, albeit a global
test involving a weighted integral over $x$. The typical $x$-values
contributing are of order unity, provided $a_{n}$ and $a_{r}$ are chosen of
that order.

The spectral function $\hat{\sigma}^{T}$ coincides up to its sign with
$\rho_{1}$, while using $x=Q^{2}/(2Tq^{0})$, it follows from Eqs. (32) that

$\begin{split}\frac{q^{i}q^{j}}{\boldsymbol{q}^{2}}\rho^{ij}(u,q)-\rho^{00}(u,q)=-\sigma^{L}(q^{0},Q^{2})&=\left(\rho_{1}(q^{0},Q^{2})-\frac{\rho_{2}(q^{0},Q^{2})}{2x}\right)-\frac{2xT^{2}}{Q^{2}}\rho_{2}(q^{0},Q^{2})\>.\end{split}$
(39)

Thus the longitudinal spectral function $\sigma^{L}$, which can be probed in
lattice QCD via Eq. (38), measures the size of the correction to the Callan-
Gross relation (47), up to a term suppressed by the factor $1/Q^{2}$ in the
Bjorken limit. Since in QCD the Callan-Gross relation is asymptotically
violated by terms of order $\alpha_{s}(Q^{2})$ Manohar (1992), where
$\alpha_{s}=g^{2}(Q^{2})/(4\pi)$ is the running coupling, the power-suppressed
term is parametrically irrelevant. Thus, in the Bjorken limit, one expects
$H_{E}^{L}$ to be suppressed (logarithmically in $Q^{2}$) as compared to
$H_{E}^{T}$. It is interesting at this point to note that the longitudinal
channel $H_{E}^{L}$ vanishes for $Q^{2}=0$, while $H_{E}^{T}$ probes the real-
photon emissivity. Thus, while the Bjorken limit is very far from lightlike
kinematics in the sense that $Q^{2}\gg T^{2}$, a similar suppression of the
longitudinal relative to the transverse channel is expected.

In Fig. 3, left panel, we illustrate the approach of $-\hat{\sigma}^{T}$ and
$(\hat{\sigma}^{L}-\hat{\sigma}^{T})$ to their common Bjorken limit in the
free theory. We observe that in both cases the departure from the limit value
increases with $x$, and that $-\hat{\sigma}^{T}$ approaches this limit faster.
The right panel allows one to judge the size of the breaking of the Callan-
Gross relation as a function of $x$ in relation to the magnitude of
$\rho_{1}$.

As noted below Eq. (38), by computing the correlators on the left-hand side of
the equation on the lattice for increasing values of the virtuality $Q^{2}$,
one could observe the approach to Bjorken scaling. In Fig. 4, we show as an
example the corresponding analysis performed in the free theory. By using the
spectral functions given in Ref. Laine (2013), we integrate numerically the
right-hand-side of Eq. 38. In order to keep $a_{n}$ and $a_{r}$ fixed and
$O(1)$ as we increase $Q^{2}$, we set $\omega_{n}=\omega_{r}/2=Q^{2}/(2T)$, a
choice that corresponds to $a_{n}=1$ and $a_{r}=2$. Figure 4 displays the
difference of Euclidean correlators for two choices of polarization which, due
to the Callan-Gross relation, converge to the same value in the Bjorken limit.
The figure illustrates again that the convergence is faster for the transverse
channel $H_{E}^{T}$ than for the difference $(H_{E}^{L}-H_{E}^{T})$. One can
also read off from the figure how suppressed $H_{E}^{L}$ is in comparison to
$H_{E}^{T}$.

Figure 3: Left: Approach to the Bjorken limit of the transverse as well as
transverse-minus-longitudinal spectral function, in the case of a plasma of
free fermions. The free spectral functions from Laine (2013) have been used.
Right: Breaking of the Callan-Gross relation $\rho_{2}=2x\rho_{1}$ at finite
$Q^{2}$ in the free theory. The functions $\rho_{1}$ and $\rho_{2}$ have been
computed from the free spectral functions using Eq. (32). In both panels,
$N_{c}=3$ and $\sum_{f}{\cal Q}_{f}^{2}=1$ has been chosen.  Figure 4: Fixed-
virtuality correlation functions (38) in the free theory, as functions of
$Q^{2}$. As we increase $\omega_{n}$, we set $Q^{2}=2T\omega_{n}$ and
$\omega_{r}=2\omega_{n}$, i.e. $a_{n}=1$ and $a_{r}=2$ are kept constant. Two
choices of polarization are shown, that converge to the same value in the
Bjorken limit. The difference between these two polarizations measures the
breaking of the Callan-Gross relation at the correlator level, up to a
$1/Q^{2}$ suppressed term (see Eq. (39)).

##  V Parton model

In this section, we review the partonic interpretation of the structure
functions of the thermal medium. The steps follow closely the textbook
treatment of unpolarized DIS on the nucleon, a small but significant
difference being the parametrization of the parton momentum; see Eq. (41)
below.

Consider the elastic electron–parton collision via one-photon exchange; see
the right panel of figure 1. The parton is assumed to be a Dirac particle and
to carry an electric charge equal to ${\cal Q}_{f}$ times the electron charge.
It has initial momentum $p_{\xi}$, while the electron has initial momentum
$k$. The differential cross-section reads (see e.g. Ref. Peskin and Schroeder
(1995), Eq. (14.3))

$\frac{d\sigma}{d\hat{t}}=\frac{2\pi\alpha^{2}}{\hat{s}^{2}}\,{\cal
Q}_{f}^{2}\;\frac{\hat{s}^{2}+(\hat{s}+\hat{t})^{2}}{\hat{t}^{2}},$ (40)

where $\hat{s}$ and $\hat{t}$ are the Mandelstam variables of the
electron–parton system. The variable $\hat{t}=q^{2}=-Q^{2}$ is the virtuality
of the photon and coincides with the Mandelstam variable $t$ of the
electron–fluid-cell collision. On the other hand,
$\hat{s}=(p_{\xi}+k)^{2}\approx 2p_{\xi}\cdot k$ is specific to the
electron–parton collision.

In a frame in which the fluid cell has a large $\gamma$ factor, we express the
fact that $p_{\xi}$ is to a good approximation collinear with the fluid
velocity $u$,

$p_{\xi}=\xi\,m_{T}\,u.$ (41)

Let $f_{f}(\xi)\,d\xi$ represent the number of partons of type $f$ in the
fluid cell carrying momentum $\xi\>m_{T}\>u$. One does not distinguish between
quarks and antiquarks, e.g. for the up quark
$f_{u}(\xi)=u_{\uparrow}(\xi)+u_{\downarrow}(\xi)+\bar{u}_{\uparrow}(\xi)+\bar{u}_{\downarrow}(\xi)$.
One expects this number to be of order the volume of the cell. Translating the
electron–parton cross-section into an electron–fluid-cell cross-section, we
obtain (see Peskin and Schroeder (1995) Eq. (14.8))

$\frac{d\sigma}{dQ^{2}}=\sum_{f}{\cal
Q}_{f}^{2}\,f_{f}(\xi)d\xi\;\frac{2\pi\alpha^{2}}{Q^{4}}\,\Big{[}1+\Big{(}1-\frac{Q^{2}}{2\xi
m_{T}(u\cdot k)}\Big{)}^{2}\Big{]}.$ (42)

Recalling the definition (Eq. (12)) of the two dimensionless kinematic
variables $x$ and $y$, Eq. (42) becomes

$\frac{d\sigma}{dQ^{2}}=\sum_{f}{\cal
Q}_{f}^{2}\,f_{f}(\xi)d\xi\;\frac{2\pi\alpha^{2}}{Q^{4}}\,\Big{[}1+\Big{(}1-\frac{xy}{\xi}\Big{)}^{2}\Big{]}.$
(43)

We found Eq. (14) for $\frac{d^{2}\sigma}{dE^{\prime}d\Omega}$ in the fluid
rest frame in terms of its spectral functions. Using the conversion444First,
show that
$\frac{d^{2}\sigma}{dE^{\prime}d\Omega}=\frac{d^{2}\sigma}{dxdyd\phi}\cdot\frac{E^{\prime}}{m_{T}(u\cdot
q)}$ by showing that
$\absolutevalue{\frac{\partial(x,y)}{\partial(E^{\prime},\cos\theta)}}=\frac{E^{\prime}}{m_{T}(u\cdot
q)}$; for that, the relations
$x=\frac{2EE^{\prime}(1-\cos\theta)}{2m_{T}(E-E^{\prime})}$ and
$y=\frac{E-E^{\prime}}{E}$ are useful. Second, show that
$\frac{d^{2}\sigma}{dxdyd\phi}=\frac{d^{2}\sigma}{dxdQ^{2}d\phi}\,\frac{Q^{2}}{y}$
by showing that $\frac{\partial(x,Q^{2})}{\partial(x,y)}=\frac{Q^{2}}{y}$. For
that, the relation $Q^{2}=2xyEm_{T}$ is useful. (see Manohar (1992), Eq.
(2.7))

$\frac{m_{T}(u\cdot
q)}{E^{\prime}}\;\frac{d^{2}\sigma}{dE^{\prime}d\Omega}=\frac{Q^{2}}{y}\,\frac{d^{2}\sigma}{dxdQ^{2}d\phi}\;,$
(44)

the cross-section (14) becomes

$\frac{d^{2}\sigma}{L^{3}dQ^{2}dx}=\frac{2m_{T}}{x}\frac{\alpha^{2}}{Q^{4}(1-e^{-\beta
q^{0}})}\,\Big{[}xy^{2}\rho_{1}(u\cdot q,Q^{2})+(1-y)\rho_{2}(u\cdot
q,Q^{2})\Big{]}.$ (45)

Note that Eq. (45) corresponds to the cross-section for scattering on the full
thermal system, which occupies a volume $L^{3}$ in its rest-frame. In this
section, we are considering the lepton scattering on a particular subvolume of
the full system. Since the cross-section scales with the volume, we now
reinterpret the factor $L^{3}$ as being the rest-frame volume of that
particular fluid cell.

The cross-section (45) is in a form to which the parton-model calculation can
easily be matched. By analogy with DIS on the proton, we anticipate that $\xi$
is proportional to $x$. Noting that the cross-section has the form of a
quadratic polynomial in $y$, the terms in the square bracket proportional to
$-y$ and $1$ are equal in Eq. (45), and they are also equal in Eq. (43)
precisely when

$x/\xi=1\,.$ (46)

Continuing from here, matching the ratio of the term proportional to $y^{2}$
to the term proportional to $1$ between Eqs. (45) and (43) yields

$\rho_{2}(u\cdot q,Q^{2})=2x\rho_{1}(u\cdot q,Q^{2}).$ (47)

This is the Callan-Gross relation. Matching the absolute normalization of the
cross-sections for $y=1$, one obtains

$\sum_{f}{\cal
Q}_{f}^{2}f_{f}(\xi)=L^{3}\frac{m_{T}}{\pi}\frac{\rho_{1}(u\cdot
q,Q^{2})}{1-e^{-\beta q^{0}}}.$ (48)

Recalling the KMS relation (13) between spectral and structure functions, the
choice $m_{T}=T$ allows us to rewrite

$F_{1}(u\cdot q,Q^{2})=\frac{1}{4L^{3}T}\sum_{f}{\cal Q}_{f}^{2}\;f_{f}(x).$
(49)

In words, $4\,F_{1}\cdot dx$ is the ${\cal Q}_{f}^{2}$-weighted number of
partons carrying a momentum $xT$ times the fluid four-velocity $u^{\mu}$ per
unit transverse area in a slab of fluid which in its rest frame has thickness
$1/T$ in the longitudinal direction. This makes the choice $m_{T}=T$ natural.

The interpretation of the structure function as a parton distribution function
is confirmed by the following observation. The first moment
$\int_{0}^{\infty}dxf(x)\;xTu^{k}$ should yield the spatial momentum of the
fluid cell carried by quarks and antiquarks. The quantity playing the role of
momentum density in ideal hydrodynamics is $(e+p)u^{k}$, therefore we would
expect to find $(e+p)L^{3}u^{k}$ for the first moment in the case of free
quarks. The connection (49) of the parton distribution function $f$ to the
structure function $F_{1}$ together with the $n=2$ moment sum rule obeyed by
$F_{1}$ (Eq. (24)), which is derived from the operator-product expansion,
independently confirm this expectation; see the remark below Eq. (79).

##  VI Conclusions and outlook

In this paper we have provided a physics interpretation of the two thermal
spectral functions characterizing the electric-current correlator in the
spacelike regime. They are related via the KMS relation to the structure
functions which describe the total cross-section for a light lepton scattering
on the thermal medium. What makes this observation interesting is that a lot
is known about the physics of DIS, in particular, we have shown that the same
moment sum rules apply to the thermal QCD medium as to the nucleon. Therefore,
the first few moments of the structure functions in the DIS limit can be
computed in lattice QCD via the expectation value of local twist-two
operators. Furthermore, Euclidean correlation functions can be computed in
lattice QCD probing the transverse and the longitudinal spectral function at a
fixed spacelike virtuality $Q^{2}$, thus allowing the onset of the DIS regime
to be probed. This approach offers a rare opportunity to answer the question,
“On what length scales do the quarks and gluons in the quark-gluon plasma
start to become ‘visible’?” in a theoretically clean way. Importantly, none of
these lattice-QCD based studies requires addressing the numerically ill-posed
problem of resolving the $x$-dependence of the structure functions from the
Euclidean correlators. On the other hand, the DIS kinematic regime is
numerically challenging to reach, since it involves using an imaginary spatial
momentum close in magnitude to the Matsubara frequency, requiring one to
achieve excellent control over non-static screening correlators at long
distances. For that purpose, a dedicated study of the corresponding non-static
correlation lengths Brandt et al. (2014) is a good preparatory step.

It is well-known from the OPE analysis of standard DIS on the nucleon that the
non-trivial one-loop anomalous dimensions of the twist-two operators and their
mixing lead to a modified prediction of the infinite-$Q^{2}$ limit of the
structure functions, as compared to a leading-order analysis. The same
observation applies to the structure functions of the thermal medium. As we
have verified explicitly for the first two non-trivial moments, the moment sum
rules (24) using leading-order Wilson coefficients are consistent with the
structure functions of thermal, non-interacting quarks. There is no analogue
of this calculation in standard DIS, because a nucleon would simply not exist
in the absence of the ${\rm SU}(3)_{\rm color}$ gauge field, whereas a plasma
of non-interacting quarks makes perfect sense. In QCD however, even at a very
high temperature where one might have expected interactions to play a
subleading role, the infinite-$Q^{2}$ limit of the structure functions differs
from that of non-interacting quarks by an amount that is not suppressed by the
strong coupling constant. Indeed, the $Q^{2}\to\infty$ asymptotic momentum
fractions carried by the quarks and the gluons according to the OPE with NLO
Wilson coefficients is the same as in DIS on the nucleon – see the discussion
below Eq. (31), whereas explicitly taking the $n=2$ moment of the leading-
order spectral function (16) leads to a different result. Thus it is clear
that the structure functions of free quarks in the DIS regime are unstable
against ‘turning on’ their interactions. Clearly, this issue deserves further
investigation, which could be carried out with the help of the next-to-leading
results Jackson (2019); Jackson and Laine (2019) for the vector spectral
functions. We also conclude from these remarks that the content of the moment
sum rules is far from trivial, and could be used to test future weak-coupling
calculations of the thermal spectral functions, including the sophisticated
resummations performed in this context; see references Arnold et al. (2001a,
b); Ghiglieri et al. (2013) for the case of vanishing virtuality. In
particular, calculations of thermal spectral functions in the DIS regime can
benefit from the existing three-loop results for the splitting functions and
anomalous dimensions of the twist-two operators Vogt et al. (2004).

###### Acknowledgements.

We thank Csaba Török for discussions. This work was supported by the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program through Grant Agreement No. 771971-SIMDAMA. The work of
M.C. is supported by the European Union’s Horizon 2020 research and innovation
program under the Marie Skłodowska-Curie Grant Agreement No. 843134-multiQCD.

##  Appendix A Derivation of the moment sum rules

In this appendix we provide a derivation of Eq. (24).

###  A.1 OPE prediction for the invariant amplitudes of the time-ordered
correlator

Let $G_{T}^{\mu\nu}(q)$ be the thermal expectation value of the Fourier-
transformed, time-ordered product of electromagnetic currents,

$G_{T}^{\mu\nu}(q)=i\int\mathrm{d}^{4}x\>e^{iq\cdot x}\langle{\rm
T}\\{j^{\mu}(x)j^{\nu}(0)\\}\rangle=\langle t^{\mu\nu}(q)\rangle.$ (50)

In the limit $Q^{2}=-q^{2}\to\infty$, the OPE of $t^{\mu\nu}(q)$ given in Eq.
(19) holds. Recalling the definition (25) of the reduced thermal expectation
values $\langle O_{nj}\rangle$ in terms of their corresponding twist-two
operator,

$\langle O_{nj}^{\mu_{1}\dots\mu_{n}}\rangle=T^{n}[u^{\mu_{1}}\dots
u^{\mu_{n}}-\mathrm{traces}]\langle O_{nj}\rangle\>,$ (51)

we have

$\begin{split}G_{T}^{\mu\nu}(q)&\overset{Q\to\infty}{\sim}\sum_{n=2,4,\dots}\sum_{f,j}2\mathcal{Q}_{f}^{2}M_{fj}(Q,\tilde{\mu})\biggl{[}\bigl{(}-g^{\mu\nu}+\frac{q^{\mu}q^{\nu}}{q^{2}}\bigr{)}\bigl{(}\frac{2T(u\cdot
q)}{Q^{2}}\bigr{)}^{n}\langle O_{nj}\rangle\\\ &+\bigl{(}u^{\mu}-(u\cdot
q)\frac{q^{\mu}}{q^{2}}\bigr{)}\bigl{(}u^{\nu}-(u\cdot
q)\frac{q^{\nu}}{q^{2}}\bigr{)}\bigl{(}2T\bigr{)}^{n}\frac{(u\cdot
q)^{n-2}}{(Q^{2})^{n-1}}\langle O_{nj}\rangle\biggr{]}\>.\end{split}$ (52)

We consider the decomposition of the time-ordered correlation function
$G_{T}^{\mu\nu}$ in invariant amplitudes,

$G_{T}^{\mu\nu}(q)=\tilde{F}_{1}(u\cdot
q,Q^{2})\bigl{(}-g^{\mu\nu}+\frac{q^{\mu}q^{\nu}}{q^{2}}\bigr{)}+\tilde{F}_{2}(u\cdot
q,Q^{2})\frac{T}{u\cdot q}\bigl{(}u^{\mu}-(u\cdot
q)\frac{q^{\mu}}{q^{2}}\bigr{)}\bigl{(}u^{\nu}-(u\cdot
q)\frac{q^{\nu}}{q^{2}}\bigr{)}\>,$ (53)

and find the following expansions for the invariant amplitudes,

$\displaystyle\tilde{F}_{1}(u\cdot q,Q^{2})$
$\displaystyle\overset{Q\to\infty}{\sim}$
$\displaystyle\sum_{n=2,4,\dots}\sum_{f,j}2\>\mathcal{Q}_{f}^{2}\>M_{fj}(Q,\tilde{\mu})\>\bigl{(}\frac{2T(u\cdot
q)}{Q^{2}}\bigr{)}^{n}\langle O_{nj}\rangle$ (54)
$\displaystyle\tilde{F}_{2}(u\cdot q,Q^{2})$
$\displaystyle\overset{Q\to\infty}{\sim}$
$\displaystyle\sum_{n=2,4,\dots}\sum_{f,j}4\>\mathcal{Q}_{f}^{2}\>M_{fj}(Q,\tilde{\mu})\>\bigl{(}\frac{2T(u\cdot
q)}{Q^{2}}\bigr{)}^{n-1}\langle O_{nj}\rangle\>.$ (55)

###  A.2 Dispersive integral for the retarded correlator

In the rest frame of the plasma, we introduce the retarded and advanced
correlators, as well as the spectral function,

$G_{R}^{\mu\nu}(q)=i\int\mathrm{d}^{4}x\>e^{iq\cdot
x}\>\theta(x^{0})\>\langle[j^{\mu}(x),j^{\nu}(0)]\rangle$ (56)
$G_{A}^{\mu\nu}(q)=-i\int\mathrm{d}^{4}x\>e^{iq\cdot
x}\>\theta(-x^{0})\>\langle[j^{\mu}(x),j^{\nu}(0)]\rangle$ (57)
$\rho^{\mu\nu}(q)=\int\mathrm{d}^{4}x\>e^{iq\cdot
x}\>\langle[j^{\mu}(x),j^{\nu}(0)]\rangle=-i(G_{R}^{\mu\nu}(q)-G_{A}^{\mu\nu}(q))$
(58)

We also introduce functions that exhibit the virtuality dependence explicitly,

$H_{R,A}^{\mu\nu}(q^{0},Q^{2};\;\hat{q})\equiv
G_{R,A}^{\mu\nu}\bigl{(}q^{0},\boldsymbol{q}=\sqrt{(q^{0})^{2}+Q^{2}}\>\hat{q}\bigr{)}\>,$
(59)

where $\hat{q}\equiv\boldsymbol{q}/\absolutevalue{\boldsymbol{q}}$ is a unit
vector. We now assume that, at fixed $Q^{2}$, $H_{R}^{\mu\nu}$ is analytic for
$\imaginary q^{0}>0$ and $H_{A}$ for $\imaginary q^{0}<0$, as was shown in
Meyer (2018) in the lightlike case $Q^{2}=0$. Note that the same property
holds in the case of DIS on the nucleon Peskin and Schroeder (1995). We define

$H^{\mu\nu}(q^{0},Q^{2};\;\hat{q})\equiv\begin{cases}H_{R}^{\mu\nu}(q^{0},Q^{2};\;\hat{q})&\imaginary
q^{0}>0\\\ H_{A}^{\mu\nu}(q^{0},Q^{2};\;\hat{q})&\imaginary
q^{0}<0\end{cases}\>.$ (60)

We now integrate along the contour $C$ in the complex-$q^{0}$ plane consisting
of a large half-circle in the upper plane with its diameter running along the
real axis from left to right, and a second half-circle in the lower plane with
its diameter running along the real axis from right to left; see Fig. 1 in
Ref. Meyer (2018). We obtain

$\begin{split}\frac{1}{2\pi
i}\oint_{C}\frac{H^{\mu\nu}(\tilde{q}^{0},Q^{2};\;\hat{q})}{(\tilde{q}^{0}-q^{0})^{m}}\>\mathrm{d}\tilde{q}^{0}&=\frac{1}{(m-1)!}\biggl{[}\frac{\mathrm{d}^{m-1}}{\mathrm{d}(\tilde{q}^{0})^{m-1}}\>H^{\mu\nu}(\tilde{q}^{0},Q^{2};\;\hat{q})\biggr{]}_{\tilde{q}^{0}=q^{0}}\\\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{\rho^{\mu\nu}(\tilde{q}^{0},Q^{2};\;\hat{q})}{(\tilde{q}^{0}-q^{0})^{m}}\>\mathrm{d}\tilde{q}^{0}\>,\quad
m=1,2,\dots\>.\end{split}$ (61)

Here
$\rho^{\mu\nu}(q^{0},Q^{2};\;\hat{q})\equiv\rho^{\mu\nu}(q^{0},\boldsymbol{q}=\sqrt{(q^{0})^{2}+Q^{2}}\>\hat{q})$
is the fixed-virtuality spectral function. Setting $q^{0}=i\epsilon$, we
obtain

$\frac{1}{(m-1)!}\biggl{[}\frac{\mathrm{d}^{m-1}}{\mathrm{d}(\tilde{q}^{0})^{m-1}}\>H_{R}^{\mu\nu}(\tilde{q}^{0},Q^{2};\;\hat{q})\biggr{]}_{\tilde{q}^{0}=i\epsilon}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{\rho^{\mu\nu}(\tilde{q}^{0},Q^{2};\;\hat{q})}{(\tilde{q}^{0}-i\epsilon)^{m}}\>\mathrm{d}\tilde{q}^{0}\>.$
(62)

###  A.3 Relation between the retarded and the time-ordered correlator

The relation of $W_{>}^{\mu\nu}$ to the time-ordered correlator can be worked
out starting from the identity

${\rm
T}\left\\{j^{\mu}(x)\,j^{\nu}(0)\right\\}=\theta(x^{0})\,[j^{\mu}(x),\,j^{\nu}(0)]+j^{\nu}(0)j^{\mu}(x).$
(63)

Second, using translational invariance, we have

$\int d^{4}x\;e^{iq\cdot x}\,\langle n|j^{\nu}(0)\,j^{\mu}(x)|n\rangle=\int
d^{4}x\;e^{-iq\cdot x}\,\langle n|j^{\nu}(x)\,j^{\mu}(0)|n\rangle.$ (64)

Therefore

$\displaystyle i\int d^{4}x\;e^{iq\cdot x}\,\langle n|{\rm
T}\left\\{j^{\nu}(0)\,j^{\mu}(x)\right\\}|n\rangle$ $\displaystyle=$
$\displaystyle i\int d^{4}x\;e^{iq\cdot x}\,\theta(x^{0})\,\langle
n|[j^{\mu}(x),\,j^{\nu}(0)]|n\rangle$ $\displaystyle+i\int d^{4}x\;e^{-iq\cdot
x}\,\langle n|j^{\nu}(x)\,j^{\mu}(0)|n\rangle.$

Recalling the definitions

$\displaystyle G_{T}^{\mu\nu}(q)$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{n}e^{-\beta E_{n}}\,i\int d^{4}x\;e^{iq\cdot
x}\,\langle n|{\rm T}\left\\{j^{\nu}(0)\,j^{\mu}(x)\right\\}|n\rangle,$ (66)
$\displaystyle G_{R}^{\mu\nu}(q)$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{n}e^{-\beta E_{n}}\,i\int d^{4}x\;e^{iq\cdot
x}\,\theta(x^{0})\,\langle n|[j^{\mu}(x),\,j^{\nu}(0)]|n\rangle,$ (67)

we have the identity

$G_{T}^{\mu\nu}(q)=G_{R}^{\mu\nu}(q)+i\,\frac{\rho^{\nu\mu}(-q)}{1-e^{\beta
q^{0}}}.$ (68)

However, for $\mu=\nu$, the spectral function $\rho^{\nu\mu}$ is real,
therefore

${\rm Re}(G_{T}^{\mu\nu}(q))={\rm
Re}(G_{R}^{\mu\nu}(q))\qquad(\mu=\nu,\;q^{0}\in\mathbb{R}).$ (69)

###  A.4 The sum rules

It follows from Eq. (69) that

$\real(G^{11}_{R}(q))=\real(G^{11}_{T}(q))\quad q^{0}\in\mathbb{R}\>,$ (70)

and the equality also holds for $q^{0}=i\epsilon$. In the rest frame of the
fluid, and choosing the orientation $q=(q^{0},0,0,\sqrt{(q^{0})^{2}+Q^{2}})$,
we find that

$\real(H_{R}^{11}(q^{0},Q^{2};\;\hat{z}))=\real(\tilde{F}_{1}(q^{0},Q^{2}))\>,$
(71)

and Eq. (62) reads in this case

$\frac{1}{(m-1)!}\biggl{[}\frac{\mathrm{d}^{m-1}}{\mathrm{d}(\tilde{q}^{0})^{m-1}}\>\real(\tilde{F}_{1}(\tilde{q}^{0},Q^{2}))\biggr{]}_{\tilde{q}^{0}=i\epsilon}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{\rho_{1}(\tilde{q}^{0},Q^{2})}{(\tilde{q}^{0}-i\epsilon)^{m}}\>\mathrm{d}\tilde{q}^{0}\>,$
(72)

where the functions $\rho_{1,2}$ are defined via a decomposition of the
spectral function $\rho^{\mu\nu}$ analogous to Eq. (53), see Eqs. (13) and
(9). In the limit $Q\to\infty$, we can insert the expansion of Eq. (54). Due
to the fact that in the limit $\epsilon\to 0$

$\frac{1}{(m-1)!}\biggl{[}\frac{\mathrm{d}^{m-1}}{\mathrm{d}(q^{0})^{m-1}}\>\bigl{(}\frac{2Tq^{0}}{Q^{2}}\bigr{)}^{n}\biggr{]}_{q^{0}=i\epsilon}=\biggl{(}\frac{2T}{Q^{2}}\biggr{)}^{n}\delta_{n,m-1}\>,$
(73)

we find

$\biggl{(}\frac{2T}{Q^{2}}\biggr{)}^{n}\sum_{f,j}2\>\mathcal{Q}_{j}^{2}\>M_{fj}(Q,\tilde{\mu})\langle
O_{nj}\rangle=\frac{1}{\pi}\int_{0}^{\infty}\rho_{1}(\tilde{q}^{0},Q^{2})\;\frac{(\tilde{q}^{0})^{n+1}}{((\tilde{q}^{0})^{2}+\epsilon^{2})^{n+1}}\>\mathrm{d}\tilde{q}^{0}\>,\quad
n=2,4,\dots\>.$ (74)

To rewrite the right-hand-side we have used the fact that
$\rho_{1}(q^{0},Q^{2})$ is odd in its first argument and that $n+1$ is always
an odd number in this context. We now remark that if we first truncate
$\rho_{1}$ to its leading-twist contribution, the parameter $\epsilon$, which
plays the role of an infrared regulator, can be set to zero. This expectation
is mathematically confirmed by the free-quark thermal spectral function; and
once only the leading-twist part is kept, partons with $x\gg 1$ are expected
to be exponentially rare on physical grounds, so that
$(\rho_{1})_{\textrm{leading-twist}}$ should fall off exponentially in $x$.
Proceeding in this way and changing the integration variable to
$x=Q^{2}/(2Tq^{0})$, we obtain the moment sum rule

$\frac{1}{2\pi}\int_{0}^{\infty}\mathrm{d}x\>x^{n-1}[\rho_{1}(x,Q^{2})]_{\textrm{leading-
twist}}=\sum_{f,j}\mathcal{Q}_{f}^{2}M_{fj}(Q,\tilde{\mu})\langle
O_{nj}\rangle\>,\quad n=2,4,\dots\>.$ (75)

Taking into account the relation (13), and considering that $e^{-\beta q^{0}}$
is zero in the Bjorken limit, Eq. (24) is recovered from the sum rule Eq.
(75). We remark that having to isolate the leading-twist part of the spectral
function in the dispersion integral distinguishes the thermal system from the
case of the proton, for which this operation is not necessary. This difference
is related to the physical range of the variable $x$ being unbounded from
above in the thermal case.

##  Appendix B Verifying leading-twist predictions for the structure functions
of the plasma of non-interacting massless Dirac fermions

In this appendix, we explicitly verify the validity of the moment sum rules
for non-interacting quarks using the leading-order Wilson coefficients, i.e.
Eq. (24) with $M_{fj}\doteq\delta_{fj}$.

We note that the Euclidean-notation operator which coincides with
$O_{nf}^{\mu_{1}\dots\mu_{n}}$ defined in Eq. (20) when all indices are
temporal is

$O_{{\rm
E},nf}^{\mu_{1}\dots\mu_{n}}=\frac{1}{2}\left(\frac{-1}{2}\right)^{n-1}{\cal
S}\\{\bar{\psi}_{f}\gamma^{\mu_{1}}\overleftrightarrow{D}^{\mu_{2}}\dots\overleftrightarrow{D}^{\mu_{n}}\psi_{f}\\},$
(76)

since $\frac{\partial}{\partial x^{0}}=i\frac{\partial}{\partial x_{0}^{\rm
E}}$.

###  B.1 The case $n=2$

The twist-two, dimension-four operator reads

$O_{{\rm
E},2f}^{\mu_{1}\mu_{2}}=\frac{-1}{4}\,\Big{(}\bar{\psi}_{f}\gamma^{\\{\mu_{1}}\overleftrightarrow{D}^{\mu_{2}\\}}\psi_{f}-\frac{1}{4}g^{\mu_{1}\mu_{2}}\bar{\psi}_{f}\overleftrightarrow{D}\\!\\!\\!\\!\\!/\;\psi_{f}\Big{)}.$
(77)

Now

$\langle O_{2f}^{00}\rangle_{\rm rest\leavevmode\nobreak\
frame}=-\frac{6T^{4}N_{c}}{\pi^{2}}\sum_{n=1}^{\infty}\frac{(-1)^{n}}{n^{4}}=\frac{7\pi^{2}T^{4}N_{c}}{120}.$
(78)

Since $\langle
O_{2f}^{\mu_{1}\mu_{2}}\rangle=T^{2}[u^{\mu_{1}}u^{\mu_{2}}-\frac{1}{4}g^{\mu_{1}\mu_{2}}]\langle
O_{2f}\rangle$, we have for the RHS of Eq. (24)

$\frac{1}{2}\langle O_{2f}\rangle=\frac{2}{3T^{2}}\langle
O_{2f}^{00}\rangle_{\rm rest\leavevmode\nobreak\
frame}=\frac{7\pi^{2}T^{2}N_{c}}{180}.$ (79)

Using Eq. (16), one finds that the LHS of Eq. (24) agrees with Eq. (79). We
remark that the last member of Eq. (79) is equal to $\frac{1}{4T^{2}}$ times
the enthalpy density of $N_{c}$ massless Dirac fermions. This observation,
given the interpretation (49) of the structure function $F_{1}$ in terms of a
parton distribution function, confirms that the non-interacting quarks carry
the entire momentum $L^{3}(e+p)u$ of the fluid.

###  B.2 The case $n=4$

At tree-level, it suffices to write four permutations for the principal term
(not containing the metric),

$\displaystyle O_{E,4f}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ $\displaystyle=$
$\displaystyle\frac{-1}{2}\frac{6}{4!\cdot
8}\Big{(}\bar{\psi}_{f}\gamma^{\mu_{1}}\overleftrightarrow{D}^{\mu_{2}}\overleftrightarrow{D}^{\mu_{3}}\overleftrightarrow{D}^{\mu_{4}}\psi_{f}+{\rm
3\leavevmode\nobreak\ perms}$
$\displaystyle-\frac{1}{4}\Big{[}g^{\mu_{1}\mu_{2}}\bar{\psi}_{f}\overleftrightarrow{D}\\!\\!\\!\\!\\!/\;\overleftrightarrow{D}^{\mu_{3}}\overleftrightarrow{D}^{\mu_{4}}\psi_{f}+{\rm
perms}\Big{]}$
$\displaystyle-\frac{1}{8}\Big{[}g^{\mu_{2}\mu_{3}}\bar{\psi}_{f}\gamma^{\mu_{1}}\overleftrightarrow{D}^{2}D^{\mu_{4}}\psi_{f}+{\rm
perms}\Big{]}$
$\displaystyle+\frac{1}{12}\bar{\psi}_{f}\overleftrightarrow{D}\\!\\!\\!\\!\\!/\;\overleftrightarrow{D}^{2}\psi_{f}\Big{[}g^{\mu_{1}\mu_{2}}g^{\mu_{3}\mu_{4}}+g^{\mu_{1}\mu_{3}}g^{\mu_{2}\mu_{4}}+g^{\mu_{1}\mu_{4}}g^{\mu_{2}\mu_{3}}\Big{]}\Big{)}.$

The coefficients of the trace terms have been determined in such a way that
contracting any two indices of the expression annihilates it.

The thermal expectation value of the zero-components is given by

$\langle
O_{4j}^{0000}\rangle=\frac{-1}{2}\langle\bar{\psi}_{f}\gamma^{0}(D^{0})^{3}\psi_{f}\rangle^{T}_{\rm
vacuum}=\frac{-120T^{6}N_{c}}{\pi^{2}}\sum_{n=1}^{\infty}\frac{(-1)^{n}}{n^{6}}=\frac{31\pi^{4}T^{6}N_{c}}{252}.$
(81)

The reduced matrix element is determined by

$\displaystyle\langle O_{4j}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\rangle$
$\displaystyle=$ $\displaystyle
T^{4}\Big{(}u^{\mu_{1}}u^{\mu_{2}}u^{\mu_{3}}u^{\mu_{4}}-\frac{1}{8}\Big{[}g^{\mu_{1}\mu_{2}}u^{\mu_{3}}u^{\mu_{4}}+g^{\mu_{1}\mu_{3}}u^{\mu_{2}}u^{\mu_{4}}+\dots\Big{]}$
(82)
$\displaystyle+\frac{1}{48}\Big{[}g^{\mu_{1}\mu_{2}}g^{\mu_{3}\mu_{4}}+g^{\mu_{1}\mu_{3}}g^{\mu_{2}\mu_{4}}+g^{\mu_{1}\mu_{4}}g^{\mu_{2}\mu_{3}}\Big{]}\Big{)}\langle
O_{4j}\rangle.$

Thus

$\frac{1}{2}\langle O_{4j}\rangle=\frac{16}{5T^{4}}\langle
O_{4j}^{0000}\rangle=\frac{62\pi^{4}T^{2}N_{c}}{315}.$ (83)

One finds that this matches the LHS of Eq. (24) for $n=4$.

## References

  * Braun-Munzinger et al. (2016) P. Braun-Munzinger, V. Koch, T. Schaefer, and J. Stachel, Phys. Rept. 621, 76 (2016), eprint 1510.00442.
  * Busza et al. (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), eprint 1802.04801.
  * Ollitrault (2008) J.-Y. Ollitrault, Eur. J. Phys. 29, 275 (2008), eprint 0708.2433.
  * Teaney (2010) D. A. Teaney, _Viscous Hydrodynamics and the Quark Gluon Plasma_ (2010), pp. 207–266, eprint 0905.2433.
  * Everett et al. (2020) D. Everett et al. (JETSCAPE) (2020), eprint 2010.03928.
  * D’Eramo et al. (2019) F. D’Eramo, K. Rajagopal, and Y. Yin, JHEP 01, 172 (2019), eprint 1808.03250.
  * David (2020) G. David, Rept. Prog. Phys. 83, 046301 (2020), eprint 1907.08893.
  * Meyer (2018) H. B. Meyer, Eur. Phys. J. A 54, 192 (2018), eprint 1807.00781.
  * Ding et al. (2011) H.-T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann, et al., Phys.Rev. D83, 034504 (2011), eprint 1012.4963.
  * Brandt et al. (2013) B. B. Brandt, A. Francis, H. B. Meyer, and H. Wittig, JHEP 1303, 100 (2013), eprint 1212.4200.
  * Amato et al. (2013) A. Amato, G. Aarts, C. Allton, P. Giudice, S. Hands, et al., Phys.Rev.Lett. 111, 172001 (2013), eprint 1307.6763.
  * Aarts et al. (2015) G. Aarts, C. Allton, A. Amato, P. Giudice, S. Hands, and J.-I. Skullerud, JHEP 02, 186 (2015), eprint 1412.6411.
  * Brandt et al. (2016) B. B. Brandt, A. Francis, B. Jaeger, and H. B. Meyer, Phys. Rev. D93, 054510 (2016), eprint 1512.07249.
  * Ding et al. (2016) H.-T. Ding, O. Kaczmarek, and F. Meyer, Phys. Rev. D94, 034504 (2016), eprint 1604.06712.
  * Hatta et al. (2008) Y. Hatta, E. Iancu, and A. Mueller, JHEP 01, 063 (2008), eprint 0710.5297.
  * Iancu and Mueller (2009) E. Iancu and A. Mueller, Phys. Lett. B 681, 247 (2009), eprint 0906.3175.
  * Chambers et al. (2017) A. Chambers, R. Horsley, Y. Nakamura, H. Perlt, P. Rakow, G. Schierholz, A. Schiller, K. Somfleth, R. Young, and J. Zanotti, Phys. Rev. Lett. 118, 242001 (2017), eprint 1703.01153.
  * Manohar (1992) A. V. Manohar, in _Lake Louise Winter Institute: Symmetry and Spin in the Standard Model_ (1992), pp. 1–46, eprint hep-ph/9204208.
  * Giusti and Meyer (2011) L. Giusti and H. B. Meyer, Phys. Rev. Lett. 106, 131601 (2011), eprint 1011.2727.
  * Meyer (2011) H. B. Meyer, Eur. Phys. J. A 47, 86 (2011), eprint 1104.3708.
  * Laine (2013) M. Laine, JHEP 11, 120 (2013), eprint 1310.0164.
  * Peskin and Schroeder (1995) M. E. Peskin and D. V. Schroeder, _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995), ISBN 978-0-201-50397-5.
  * Meyer (2017) H. B. Meyer, Eur. Phys. J. C77, 616 (2017), eprint 1706.01139.
  * Brandt et al. (2014) B. Brandt, A. Francis, M. Laine, and H. Meyer, JHEP 1405, 117 (2014), eprint 1404.2404.
  * Jackson (2019) G. Jackson, Phys. Rev. D 100, 116019 (2019), eprint 1910.07552.
  * Jackson and Laine (2019) G. Jackson and M. Laine, JHEP 11, 144 (2019), eprint 1910.09567.
  * Arnold et al. (2001a) P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 0111, 057 (2001a), eprint hep-ph/0109064.
  * Arnold et al. (2001b) P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 12, 009 (2001b), eprint hep-ph/0111107.
  * Ghiglieri et al. (2013) J. Ghiglieri, J. Hong, A. Kurkela, E. Lu, G. D. Moore, et al., JHEP 1305, 010 (2013), eprint 1302.5970.
  * Vogt et al. (2004) A. Vogt, S. Moch, and J. Vermaseren, Nucl. Phys. B 691, 129 (2004), eprint hep-ph/0404111.