# Polarized Photons from the Early Stages of Relativistic Heavy-Ion Collisions

Sigtryggur Hauksson<EMAIL_ADDRESS>Institut de Physique
Théorique, CEA/Saclay, Université Paris-Saclay, 91191, Gif sur Yvette, France
Charles Gale<EMAIL_ADDRESS>Department of Physics, McGill University,
3600 University Street, Montréal, QC, Canada H3A 2T8

###### Abstract

The polarization of real photons emitted from early-time heavy-ion collisions
is calculated, concentrating on the contribution from bremsstrahlung and
quark-antiquark annihilation processes at leading order in the strong
coupling. The effect of an initial momentum space anisotropy of the parton
distribution is evaluated using a model for the non-equilibrium scattering
kernel for momentum broadening. The effect on the photon polarization is
reported for different degrees of anisotropy. The real photons emitted early
during in-medium interactions will be dominantly polarized along the beam
axis.

##  I Introduction

The theory of the nuclear strong interaction, QCD, features a transition from
a phase where the relevant degrees of freedom are quarks and gluons – at high
temperature – to one where the appropriate basis consists of composite hadrons
at a lower temperature. This transition has been predicted by several
theoretical approaches, including the non-perturbative field-theoretical
framework of lattice QCD [[See, forexample, ][,
andreferencestherein.]Aarts:2023vsf]. Decades of intense theoretical effort
have revealed that the transition from confined hadrons to partons is not a
thermodynamic phase transition in the proper sense but rather an analytic
crossover [[See, forexample, ][, andreferencestherein.]Aoki:2006we]. On the
experimental side, this exotic state of strongly interacting matter – the
quark-gluon plasma (QGP) – has been observed in the relativistic collisions of
nuclei (“heavy-ions”) performed at the Relativistic Heavy-Ion Collider (RHIC),
and its existence has been later confirmed by experiments performed at the
Large Hadron Collider (LHC) [3]. There also is strong evidence supporting the
presence of QGP in smaller systems [4].

Even though new aspects of the QGP are continuously being discovered, it is
fair to assert that the field is well poised to enter a phase of quantitative
characterization, owing in large part to the large variety of experimental
observables being measured by the experimental collaborations. Measurements
designed to probe the QGP have reported results on soft hadron collective
behavior [5], on QCD jet modification and energy loss [6], on photon [7] and
dilepton [8] production, and on many other aspects as well. The theoretical
tools developed to study the dynamics of nuclear collisions and the formation
of the QGP typically consist of multistage models, rendered necessary by the
complexity of the nuclear reaction. Recent examples of such a composite
theoretical approach are studies [9, 10, 11, 12] where the initial state model
(e.g. TRENTo [13], IP-Glasma [14]) preceded a fluid-dynamical phase ([15, 16,
17, 18]. A hadronic cascade afterburner (e.g. UrQMD [19], SMASH [20]), evolve
the final state hadrons until their measurement. Such composite models can
make statistically significant statements about transport parameters such as
shear and bulk viscosity, as well as quantify the energy loss of energetic QCD
jets [21].

As the multistage modeling of relativistic heavy-ion collisions covers a
variety of dynamical conditions ranging from far from equilibrium initial
states to almost ideal fluid dynamics, it is important to critically examine
its different eras. In searching for observables capable of revealing the
different modeling epochs, penetrating probes such as real and virtual photons
impose themselves. Electromagnetic variables are emitted at all stages of the
collision and as such can report on local conditions at their creation point
[22]. The largest uncertainty in the chain of models currently lie at the
beginning, in the time span preceding “hydrodynamization”. Early in the
history of the collision, photon emission is liable to occur in media far from
equilibrium which necessitate a dedicated theoretical treatment. The emission
of photons from non-equilibrium environments has received some recent
attention [23, 24].

This study will consider photon emission in the early stages of heavy-ion
collisions, and will focus more specifically on the polarization states of
those photons as a probe of the medium at early times. As our calculation only
relies on the medium having a pressure anisotropy, or equivalently a momentum
anisotropy in parton distribution, it holds both before hydrodynamization as
well as in the beginning of the hydrodynamic stage while pressure anisotropy
still persists. Some previous estimates for photon polarization at early times
considered leading order direct photon production channels like those of the
Compton process and $q\bar{q}$ annihilation [25, 26, 27]. It is known that an
equally important contribution as those two – at the same order in
$\alpha_{s}$, the strong coupling constant – is that associated with the
Landau-Pomeranchuk-Migdal effect (LPM) [28, 29, 30]. That contribution,
evaluated for a medium out of equilibrium forms the basis of this work. It is
fair to remind readers that the measurement of real photon polarization states
is challenging, owing to the complications related to the external conversions
into lepton pairs. The angular distribution of this pairs will reflect the
polarization state. A more realistic proposition is a measurement of virtual
photon polarization states, as measured through an internal conversion process
leading to a dilepton final state. Consequently, the goal of our work is to
first set the foundations for subsequent such evaluations and to perform a
first estimate of the polarization signature of an early, non-equilibrium,
strongly interacting medium.

Our paper is organized as follows: Section II lays out the building blocks of
our non-equilibrium formalism. The section following that one discusses the
collision kernel used to model the medium interactions. The numerical methods
used to obtain the production rate of polarized photons are discussed in
Section IV. Results and conclusion constitute Sections V and VI, respectively.
We present analytical and numerical details in the Appendices.

##  II Polarized photon emission

Figure 1: Our choice of coordinate system. The $z$-axis is chosen to be along
the beam axis. For a photon emitted at midrapidity we choose the $x$-axis to
be along its momentum. The photon (in red) can either be transversely
polarized along the beam axis ($z$-axis) or transversely polarized orthogonal
to the beam axis, i.e. along the $y$-axis. Figure 2: Photon radiation through
medium-induced bremsstrahlung off a quark Figure 3: Photon radiation through
medium-induced quark-antiquark pair annihilation

The quark-gluon plasma radiates photons through two different processes at
leading order in perturbation theory. (See [31] for higher order corrections.)
These processes are two-to-two scattering with a photon in the final state
[32, 33], and bremsstrahlung and pair annihilation with a resulting photon
[28, 30]. Two-to-two scattering on one hand and bremsstrahlung and pair
annihilation on the other hand give roughly equal contribution to photon yield
in the plasma [29]. While photons from two-to-two scattering have been studied
extensively in an anisotropic medium and have been shown to be polarized [25],
bremsstrahlung and pair annihilation photons have been studied much less in an
anisotropic plasma due to the more complicated physics involved. We consider
this now. Polarized photon emission has also been studied in other contexts,
including in holography where a background magnetic field is included [34, 35,
36], due to vortical flow in the plasma [37], and due to the chiral magnetic
effect [38, 39]. Dilepton polarization has furthermore been considered in [40,
41, 42]

In bremsstrahlung an on-shell photon is emitted collinearly off a quark or an
antiquark, see Fig. 2. In vacuum this process would be kinematically forbidden
as an on-shell quark cannot emit an on-shell photon. However, in a medium the
process is made possible to due soft gluon kicks from the medium that bring
the quark slightly off-shell. These kicks have momentum $\sim g\Lambda$ where
$\Lambda$ is a hard scale akin to temperature and $g$ is the coupling
constant. The off-shellness of the quark is therefore of order $P^{2}\sim
g^{2}\Lambda^{2}$ meaning that the emission takes time $t\sim p/P^{2}\sim
1/g^{2}\Lambda$ where $p$ is the quark momentum. During this time the quark
can receive arbitrarily many soft gluon kicks since the mean-free time between
two such kicks is also of order $1/g^{2}\Lambda$. All of these kicks need to
be included at leading order in perturbation theory [28, 30, 43]. In an
analogous fashion a quark-antiquark pair can annihilate and radiate a photon
due to medium kicks, see Fig. 3.

We will now show that photons emitted through bremsstrahlung and pair
annihilation are polarized in an anisotropic medium. This polarization can be
described by extending the framework developed in [28, 30, 23], see App. A. To
fix ideas we choose the coordinate system in Fig. 1. The $z$-axis lies along
the beam axis in a heavy-ion collision. We consider a photon at midrapidity
and orient the coordinate system so that the $x$-axis lies along its momentum
$\mathbf{k}$. The momentum of the photon can of course have any direction in
the plane transverse to the beam axis; aligning it with the short axis of the
plasma as in Fig. 1 is simply for illustration. 111In this work, the net
polarization of photons emitted from a fluid cell is independent of the
angular orientation in the transverse plane as we focus on the effect of
longitudinal expansion. This could be generalized to include transverse
expansion which breaks this symmetry. Our formalism could also easily be
extended to photons at finite rapidity. As an on-shell photon is transversely
polarized, the polarization basis can be choosen as $\epsilon_{z}=(0,0,0,1)$
and $\epsilon_{y}=(0,0,1,0)$. The photon is thus polarized along the beam axis
or transverse to the beam axis.

In Appendix A we show that the rate of producing $z$-polarized photons with
momentum $k$ through bremsstrahlung is

$\displaystyle
k\frac{d\Gamma_{z}}{d^{3}\mathbf{k}}=\frac{6\alpha_{EM}\sum_{s}q_{s}^{2}}{(2\pi)^{3}}\,2\int_{0}^{\infty}dp_{x}\;\frac{k}{8p_{x}^{2}(k+p_{x})}$
$\displaystyle\times
n_{f}(k+p_{x})\left[1-n_{f}(p_{x})\right]\left(F_{\mathrm{in}}(\zeta)A_{z}+F_{\mathrm{out}}(\zeta)A_{y}\right)$
(1)

where momenta are defined in Fig. 2 and $\zeta=k/(p_{x}+k)$ is the momentum
fraction of the photon. Similarly, the rate of producing $y$-polarized photons
is

$\displaystyle
k\frac{d\Gamma_{y}}{d^{3}\mathbf{k}}=\frac{6\alpha_{EM}\sum_{s}q_{s}^{2}}{(2\pi)^{3}}2\int_{0}^{\infty}dp_{x}\;\frac{k}{8p_{x}^{2}(k+p_{x})}$
$\displaystyle\times
n_{f}(k+p_{x})\left[1-n_{f}(p_{x})\right]\left(F_{\mathrm{in}}(\zeta)A_{y}+F_{\mathrm{out}}(\zeta)A_{z}\right)$
(2)

where $A_{z}$ and $A_{y}$ have been interchanged relative to Eq. (II). Here
$A_{z}$ and $A_{y}$ quantify the amount of momentum broadening in the $z$\-
and $y-$direction and are defined as

$\displaystyle A_{z}$
$\displaystyle=\mathrm{Re}\,\int\frac{d^{2}\mathbf{p}_{\perp}}{(2\pi)^{2}}\;2p_{z}f_{z}(\mathbf{p}_{\perp})$
$\displaystyle A_{y}$
$\displaystyle=\mathrm{Re}\,\int\frac{d^{2}\mathbf{p}_{\perp}}{(2\pi)^{2}}\;2p_{y}f_{y}(\mathbf{p}_{\perp})$
(3)

where $\mathbf{f}$ solves an integro-differential equation given below.
Furthermore,

$F_{\mathrm{in}}(\zeta)=\frac{(2-\zeta)^{2}}{\zeta}$ (4)

and

$F_{\mathrm{out}}(\zeta)=\zeta,$ (5)

are polarized splitting functions [44]. Finally, there are momentum factors
$n_{f}$ for the incoming and outgoing quarks. 222We have assumed that there is
no chiral imbalance in the medium and that the baryon chemical potential
vanishes, so that quarks and antiquarks of both helicities have the same
momentum distribution $n_{f}$. The analogous expressions for quark-antiquark
pair annihilation are given in Appendix A.

Figure 4: Definition of the quantities $\mathbf{n}$ (the vector between the
outgoing quark and photon) and $\phi_{\mathbf{q}}$ (the angle defining a soft
gluon kick of magnitude $q_{\perp}$). The photon can be polarized in the $y$\-
or the $z$-directions which are transverse to its direction of motion, see
also Fig. 1.

Eqs. (II) and (II) for polarized photon emission can be understood
intuitively. We consider a $z$-polarized photon for concreteness. As the
photon travels in the $x$-direction, the splitting plane of the photon and the
outgoing quark is spanned by $\widehat{\mathbf{e}}_{x}$ and a vector
orthogonal to $\widehat{\mathbf{e}}_{x}$ which we call $\widehat{\mathbf{n}}$,
see Fig. 4. Eqs. (II) and (II) show that we can project the vector
$\widehat{\mathbf{n}}$ on to the $y$\- and the $z$-axes and sum over the
contributions. For $\widehat{\mathbf{n}}$ in the $z$-direction, the
$z$-polarized photon is polarized in the splitting plane and the hard
splitting function is $F_{\mathrm{in}}$. Momentum broadening is quantified by
the $z$-component $A_{z}$. On the other hand, for $\widehat{\mathbf{n}}$ in
the $y$-direction, the $z$-polarized photon is polarized out of the splitting
plane and the hard splitting function is $F_{\mathrm{out}}$. Momentum
broadening is then quantified by $A_{y}$.

The rate equations for photon emission depend on the function
$\mathbf{f}=(f_{z},f_{y})$ which quantifies momentum broadening. It solves the
integro-differential equation

$2\mathbf{p}_{\perp}=i\delta
E\mathbf{f}(\mathbf{p}_{\perp})+\int\frac{d^{2}q_{\perp}}{(2\pi)^{2}}\;\mathcal{C}(\mathbf{q}_{\perp})\left[\mathbf{f}(\mathbf{p}_{\perp})-\mathbf{f}(\mathbf{p}_{\perp}+\mathbf{q}_{\perp})\right]$
(6)

where

$\delta E=\frac{k}{2p(p+k)}\left[p_{\perp}^{2}+m_{\infty}^{2}\right].$ (7)

The central ingredient in this equation is the collision kernel
$\mathcal{C}(\mathbf{q}_{\perp})$ which gives the rate for a quark to receive
soft gluon kicks of transverse momentum $\mathbf{q}_{\perp}$. The collision
kernel gives rise to a gain term and a loss term in Eq. (6).

In an isotropic medium, the collision kernel is by definition isotropic,
$\mathcal{C}(\mathbf{q}_{\perp})=\mathcal{C}(q_{\perp})$ and one can show that
$\mathbf{f}=\mathbf{p}_{\perp}\widehat{f}(p_{\perp})$. This means that
$A_{z}=A_{y}$ and thus there is no net polarization of photons emitted. In an
anisotropic medium, $\mathcal{C}(\mathbf{q}_{\perp})$ depends not only on the
magnitude of the kick $q_{\perp}$ but also on its orientation. In other words,
when writing

$\mathbf{q}_{\perp}=\left(q_{z},q_{y}\right)=q_{\perp}(\cos\phi_{\mathbf{q}},\sin\phi_{\mathbf{q}})$
(8)

the collision kernel depends on both $\phi_{\mathbf{q}}$ and $q_{\perp}$, see
Fig. 4. This leads to $\mathbf{f}$ having more complicated angular dependence
so that $A_{z}\neq A_{y}$. Therefore, photon emission from an anisotropic
medium is polarized.

##  III Model of the collision kernel in an anisotropic plasma

As previously argued, the collision kernel for soft gluon kicks
$\mathcal{C}(\mathbf{q}_{\perp})$ is anisotropic in heavy-ion collisions which
leads to polarization of photons emitted through bremsstrahlung. The ultimate
source of the anisotropy in the kernel is longitudinal expansion of the medium
along the beam axis. Such a longitudinal expansion gives pressure anisotropy
at early and intermediate times with longitudinal pressure $P_{L}$ less than
transverse pressure $P_{T}$ . On a microscopic level, this means that quark
and gluon quasiparticles have an anisotropic momentum distribution with
$\langle p_{z}^{2}\rangle<\langle p^{2}_{x}\rangle,\langle p^{2}_{y}\rangle$.
This is captured by the distribution introduced in Ref. [45]:

$f(\mathbf{p})=\sqrt{1+\xi}f_{\mathrm{iso}}(\sqrt{\mathbf{p}^{2}+\xi
p_{z}^{2}})$ (9)

where $f_{\mathrm{iso}}$ is an isotropic distribution, and $\xi>0$ quantifies
the degree of anisotropy. The prefactor $\sqrt{1+\xi}$ ensures that the number
density of quarks and gluons is the same as in equilibrium. A simple
calculation gives the pressure anisotropy $P_{T}/P_{L}$ in terms of the
momentum space anisotropy $\xi$, linking the macroscopic and microscopic
descriptions.333 Specifically,
$P_{T}/P_{L}=\frac{1}{2}\frac{\sqrt{\xi}+(\xi-1)\arctan\sqrt{\xi}}{\arctan\sqrt{\xi}-\sqrt{\xi}/(1+\xi)}$,
see e.g. [46].

Ideally, one would want to calculate the collision kernel for momentum
broadening, $\mathcal{C}(\mathbf{q}_{\perp})$, directly from Eq. (9). Such a
calculation would use that the hard quasiparticles in Eq. (9) radiate soft
gluons which are then responsible for momentum broadening. This would allow to
quantify the degree of photon polarization in a non-equilibrium medium from
first principles. Unfortunately, going from Eq. (9) to the collision kernel is
difficult in practice, partially due to instabilities that can be present in
an anisotropic plasma [47].

In this study, we will use a simple model of the collision kernel in a
longitudinally expanding medium. We take inspiration from results for the
collision kernel in thermal equilibrium at leading order in perturbation
theory [48],

$\mathcal{C}_{\mathrm{eq}}(\mathbf{q}_{\perp})=g^{2}C_{F}T\left(\frac{1}{q_{\perp}^{2}}-\frac{1}{q_{\perp}^{2}+m_{D0}^{2}}\right)$
(10)

Here $m_{D0}^{2}$ is the equilibrium Debye mass which describes screening of
electric fields. (The equilibrium kernel has furthermore been evaluated at
next-to-leading order [49], as well as on the lattice, see e.g. [50, 51].) In
our anisotropic model, we replace the equilibrium Debye mass by its
anisotropic extension444In [45] this quantity is referred to as $m_{+}^{2}$.
We have expanded in small $\xi$ in Eq. (III) but this is simply for
convenience and not fundamental to the setup we use. found in [45]

$m_{D}^{2}(\phi_{\mathbf{q}})=\left(1-\frac{2\xi}{3}\right)m_{D0}^{2}+\xi
m_{D0}^{2}\cos^{2}\phi_{\mathbf{q}}.$ (11)

so that

$\mathcal{C}(\mathbf{q}_{\perp})=g^{2}C_{F}\Lambda\Big{(}\frac{1}{q_{\perp}^{2}}-\frac{1}{q_{\perp}^{2}+m_{D}^{2}(\phi_{\mathbf{q}})}\Big{)}$
(12)

The anisotropic correction has an angular dependence with more broadening in
the $z$-direction than in the $y$-direction. To simplify calculations, we
expand the collision kernel in $\xi$, writing

$\displaystyle\mathcal{C}(\mathbf{q}_{\perp})\approx g^{2}C_{F}\Lambda$
$\displaystyle\left(\frac{1}{q_{\perp}^{2}}-\frac{1}{q_{\perp}^{2}+m_{D0}^{2}}\right.$
$\displaystyle+\left.\frac{-2\xi m_{D0}^{2}/3+\xi
m_{D0}^{2}\cos^{2}\phi_{\mathbf{q}}}{\left(q_{\perp}^{2}+m_{D0}^{2}\right)^{2}}\right)$
(13)

where $\Lambda$ is a hard scale, akin to temperature.

Eq. (III) is a toy model for the collision kernel, intended to illustrate how
an anisotropic kernel leads to polarized photon emission and to estimate the
magnitude of this effect. This toy model only includes changes to the
screening of electric fields in an anisotropic medium and not the myriad other
non-equilibrium effects that can arise. Nevertheless, this collision kernel
can be motivated by theoretical arguments making us belief that it captures
some of the salient features of the full non-equilibrium kernel.

In general the collision kernel is defined by

$\mathcal{C}(\mathbf{q}_{\perp})=g^{2}C_{F}\int\frac{dq^{0}dq_{x}}{(2\pi)^{2}}\;D_{rr}^{\mu\nu}(Q)v_{\mu}v_{\nu}\,2\pi\delta(v\cdot
Q)$ (14)

where $v^{\mu}=(1,1,0,0)$ is the four-velocity of the quark emitting a photon.
The kernel depends on the statistical correlator for gluons in the medium,

$D_{rr}(Q):=\frac{1}{2}\left\langle\left\\{A,A\right\\}\right\rangle(Q)=D_{\mathrm{ret}}(Q)\Pi(Q)D_{\mathrm{adv}}(Q)$
(15)

which characterizes the occupation density of a pair of soft gluons. We have
omitted Lorentz indices for simplicity. This statistical correlator contains
information on how the soft gluons are emitted by hard quasiparticles with
rate $\Pi(Q)$ and then propagate in the medium according to the retarded
propagator $D_{\mathrm{ret}}(Q)=i/(Q^{2}-\Pi_{\mathrm{ret}})$ and the advanced
propagator $D_{\mathrm{adv}}=-D_{\mathrm{ret}}^{*}$.

Making some heuristic assumptions allows one to employ a sum rule in [48] to
motivate the toy model for the collision kernel in Eq. (III), starting from
the definitions in Eqs. (14) and (15). The goal is to include anisotropic
corrections to the screening of chromoelectric fields, while ignoring
anisotropic corrections to the density of gluons and to change in polarization
that occurs during propagation. We work strictly at small anisotropy $\xi\ll
1$, including only effects of order $\mathcal{O}(\xi)$.

The first heuristic assumption is to employ the identity
$\Pi(Q)=\frac{\Lambda}{q^{0}}2\mathrm{Im}\,\Pi_{\mathrm{\mathrm{ret}}}$ which
is known as the KMS identity and which expresses detailed balance between
production and decay of soft gluons. This is not strictly valid in a non-
equilibrium medium and amounts to ignoring anisotropic corrections to the
density of gluons. Then one can write

$D_{rr}(Q):=\frac{\Lambda}{q^{0}}\left(D_{\mathrm{ret}}-D_{\mathrm{adv}}\right).$
(16)

At small anisotropy the retarded propagator in Eq. (16) is

$D_{\mathrm{ret}}^{\mu\nu}(Q)\approx\frac{P_{T}^{\mu\nu}}{Q^{2}-\Pi_{T}}+\frac{P_{L}^{\mu\nu}}{Q^{2}-\Pi_{L}}.$
(17)

Here we have only included anisotropic corrections to the screening as given
by $\Pi_{T}$ and $\Pi_{L}$, see App. B.

Our second heuristic approximation is to focus on the anisotropic correction
to $\Pi_{L}$ and ignore those in $\Pi_{T}$. Comparing with the equilibrium
calculation [48], this amounts to calculating anisotropic corrections to the
term $1/(q_{\perp}^{2}+m_{D0}^{2})$ in Eq. (10) while leaving the term
$1/q_{\perp}^{2}$ as is. This means that we include anisotropic corrections to
the screening of chromoelectric fields as given by a Debye mass but do not
include anistropic corrections to the screening of chromomagnetic fields.

The reason we use this approximation is that the sum rule we employ does not
work for the transverse screening in $\Pi_{T}$. This is because the term
$1/(Q^{2}-\Pi_{T})$ has a pole in the upper half complex plane of $q_{0}$
corresponding to Weibel instabilities [52]. A formal use of the sum rule would
lead to a contribution of the form $1/(q_{\perp}^{2}-m^{2})$ which is ill-
defined at $q_{\perp}=m$. The solution to this issue is to use a retarded
propagator for soft gluons that includes the mechanism by which non-Abelian
interaction caps the growth of the unstable soft gluon modes. This is beyond
the scope of this project.

Given these two approximations, one can use the sum rule in [48] nearly
directly. The longitudinal retarded self-energy at small anisotropy $\xi$ is

$\displaystyle\Pi_{L}(Q)=\Pi_{L}^{0}(x)+\xi\left[\frac{1}{6}\left(1+3\cos
2\theta_{n}\right)\frac{Q^{2}}{\mathbf{q}^{2}}m_{D0}^{2}\right.$
$\displaystyle\left.+\Pi_{L}^{0}(x)\left(\cos
2\theta_{n}-\frac{x^{2}}{2}\left(1+3\cos 2\theta_{n}\right)\right)\right]$
(18)

where $\theta_{n}$ is the angle between $\mathbf{q}$ and the anisotropy vector
$\mathbf{n}=\mathbf{e}_{z}$ which defines a preferred direction in Eq. (9)
[45]. Here

$\Pi_{L}^{0}(x)=\left(x^{2}-1\right)m_{D0}^{2}\left[\frac{x}{2}\log\frac{x+1+i\epsilon}{x-1+i\epsilon}-1\right]$
(19)

is the equilibrium value. We next do a change of variables in Eq. (14) to
$x=q_{0}/q=q_{x}/\sqrt{q_{x}^{2}+q_{\perp}^{2}}$ so that
$q_{x}^{2}=x^{2}q_{\perp}^{2}/(1-x^{2})$ and $q^{2}=q_{\perp}^{2}/(1-x^{2})$.
Then one should substitute $\cos
2\theta_{n}=2(q_{z}/q)^{2}-1=2(1-x^{2})\cos^{2}\phi_{\mathbf{q}}-1$ to get
dependence only on $x$, $q_{\perp}$ and $\phi_{\mathbf{q}}$. This gives that
the longitudinal contribution to the collision kernel is

$\frac{\Lambda
g^{2}C_{F}}{\pi}\int_{0}^{1}\frac{dx}{x}\frac{\mathrm{Im}\,\widetilde{\Pi}_{L}(x,\phi_{\mathbf{q}})}{\left(\mathbf{q}_{\perp}^{2}+\mathrm{Re}\,\widetilde{\Pi}_{L}(x,\phi_{\mathbf{q}})\right)^{2}+\left(\mathrm{Im}\,\widetilde{\Pi}_{L}(x,\phi_{\mathbf{q}})\right)^{2}}.$
(20)

where

$\displaystyle\widetilde{\Pi}_{L}(x,\phi_{\mathbf{q}})=\Pi_{L}^{0}(x)+\xi\left[\frac{1}{6}\left(1+3\cos
2\theta_{n}\right)(x^{2}-1)m_{D0}^{2}\right.$
$\displaystyle\left.+\Pi_{L}^{0}(x)\left(\cos
2\theta_{n}-\frac{x^{2}}{2}\left(1+3\cos
2\theta_{n}\right)\right)\right]\bigg{\rvert}_{\cos
2\theta_{n}=2(1-x^{2})\cos^{2}\phi_{\mathbf{q}}-1}$ (21)

has no explicit dependence on $q_{\perp}$.

@articleAurenche:2002pd, @articleAurenche:2002pd, An argument nearly
identical555 The retarded propagator
$1/\left[x^{2}q_{\perp}^{2}-q_{\perp}^{2}-(1-x^{2})\Pi_{L}(x,\phi)\right]$ in
Eq. 7 in [48] has an extra pole in the upper half complex plane in the
anisotropic case. This pole can be seen by taking the limit
$x=k_{0}/k\rightarrow\infty$ in which case the propagator becomes $\sim
1/(1-x^{2})/(q_{\perp}^{2}+\left(\frac{1}{3}-\frac{1}{3}\xi\cos^{2}\phi\right)m_{D0}^{2}+\frac{\xi}{3}x^{2}\cos^{2}\phi\,m_{D0}^{2})$
which has a pole which is parametrically of the form $x\sim\pm
iq_{\perp}/m_{D}\sqrt{\xi}$ and thus far from the real axis when $\xi\ll 1$.
This is not in contradiction with the usual properties of the retarded
propagator as we have imposed $q_{0}=q_{x}$ and then search for poles in
$q_{0}$. One can then show that the correction due to this pole to the sum
rule in Eq. 9 in [48] is $\mathcal{O}(\xi^{3/2})$ which is subleading to the
$\mathcal{O}(\xi)$ contributions we consider.  to the one in [48] then shows
that Eq. (20) is

$\displaystyle
g^{2}C_{F}\Lambda\left[\frac{1}{q_{\perp}^{2}+\lim_{x\rightarrow\infty}\widetilde{\Pi}_{L}(x,\phi)}-\frac{1}{q_{\perp}^{2}+\widetilde{\Pi}_{L}(0,\phi)}\right]$
$\displaystyle=-g^{2}C_{F}\Lambda\frac{1}{q_{\perp}^{2}+m_{D}^{2}(\phi_{\mathbf{q}})}$
(22)

since $\widetilde{\Pi}_{L}(0,\phi)=m_{D}^{2}(\phi_{\mathbf{q}})$. Our
conclusion is therefore that given our heuristic approximations, the collision
kernel is given by Eq. (III). We emphasize that this collision kernel is not
intended to capture all of the non-equilibrium physics but to focus on
anisotropic corrections to the screening of chromoelectric fields.

##  IV Numerical method

We wish to evaluate the rate of polarized photon emission in an anisotropic
medium such as that found in a longitudinally expanding quark-gluon plasma.
The starting point is Eqs. (II) and (II) which require solving the integro-
differential equation in Eq. (6) numerically, assuming the model for the
collision kernel given in Eq. (III). To solve Eq. (6) we go to impact
parameter space, i.e. the space Fourier conjugate to $\mathbf{p}_{\perp}$.
Defining

$\mathbf{f}(\mathbf{b})=\int\frac{d^{2}p_{\perp}}{(2\pi)^{2}}\;e^{i\mathbf{p}_{\perp}\cdot\mathbf{b}}\,\mathbf{f}(\mathbf{p}_{\perp}).$
(23)

the equation we wish to solve becomes

$-2i\nabla_{b}\delta^{(2)}(\mathbf{b})=\frac{ik}{2p(p+k)}\left[-\nabla_{b}^{2}+m_{\infty}^{2}\right]\mathbf{f}(\mathbf{b})+\mathcal{C}(\mathbf{b})\mathbf{f}(\mathbf{b}).$
(24)

where

$\mathcal{C}(\mathbf{b})=\int\frac{d^{2}p_{\perp}}{(2\pi)^{2}}\;\left[1-e^{i\mathbf{p}_{\perp}\cdot\mathbf{b}}\,\right]\mathcal{C}(\mathbf{p}_{\perp}).$
(25)

A straightforward calculation shows that the collision kernel from Eq. (III)
is

$\mathcal{C}(\mathbf{b})=\mathcal{C}_{0}(b)+\xi\,\mathcal{C}_{1}^{(a)}(b)+\xi\cos
2\beta\,\mathcal{C}_{1}^{(b)}(b)$ (26)

in impact parameter space where
$\mathbf{b}=(b_{z},b_{y})=\left(\cos\beta,\sin\beta\right)b$. The terms of the
collision kernel are given by

$\mathcal{C}_{0}(\mathbf{b})=\frac{g^{2}C_{F}T}{2\pi}\left[K_{0}(m_{D0}b)+\gamma_{E}+\log\frac{m_{D0}b}{2}\right],$
(27)
$\mathcal{C}_{1}^{(a)}(b)=g^{2}C_{F}T\frac{1}{8\pi}\frac{M^{2}}{m_{D0}^{2}}\left(m_{D0}bK_{1}(m_{D0}b)-1\right)$
(28)

and

$\displaystyle\mathcal{C}_{1}^{(b)}(b)=g^{2}C_{F}T\frac{M^{2}b^{2}}{4\pi}\Bigg{[}\frac{2}{(m_{D0}b)^{4}}$
$\displaystyle-\frac{1}{2m_{D0}b}K_{1}(m_{D0}b)-\frac{1}{\left(m_{D0}b\right)^{2}}K_{2}(m_{D0}b)\Bigg{]}$
(29)

In order to solve Eq. (24) we do an expansion in small $\xi$, giving

$\mathbf{f}(\mathbf{b})=\mathbf{f}_{0}(\mathbf{b})+\xi\mathbf{f}_{1}(\mathbf{b})+\dots$
(30)

The zeroth order solution in $\xi$ satisfies the usual isotropic equation

$\frac{ik}{2p(p+k)}\left[-\nabla_{b}^{2}+m_{\infty}^{2}\right]\mathbf{f}_{0}(\mathbf{b})+\mathcal{C}_{0}(b)\mathbf{f}_{0}(\mathbf{b})=-2i\nabla_{b}\delta^{(2)}(\mathbf{b})$
(31)

and can be shown to have angular dependence
$\mathbf{f}_{0}(\mathbf{b})\sim(\cos\beta,\sin\beta)$. The first order
satisfies

$\displaystyle\frac{ik}{2p(p+k)}$
$\displaystyle\left[-\nabla_{b}^{2}+m_{\infty}^{2}\right]\mathbf{f}_{1}(\mathbf{b})+\mathcal{C}_{0}(b)\mathbf{f}_{1}(\mathbf{b})$
(32) $\displaystyle=-\left[\mathcal{C}_{1}^{(a)}(b)+\cos
2\beta\,\mathcal{C}_{1}^{(b)}(b)\right]\mathbf{f}_{0}(\mathbf{b}).$

Due to the angular dependence of the right hand side we can write in full
generality

$f_{1z}(\mathbf{b})=\cos\beta f_{1}^{(1z)}(b)+\cos 3\beta\,f_{1}^{(3)}(b)$
(33)

and

$f_{1y}(\mathbf{b})=\sin\beta f_{1}^{(1y)}(b)+\sin 3\beta\,f_{1}^{(3)}(b)$
(34)

where these functions solve

$\mathcal{K}\left[f_{1}^{(1z)}(\mathbf{b})\right]+\mathcal{C}_{0}(b)f_{1}^{(1z)}(\mathbf{b})=-\left[\mathcal{C}_{1}^{(a)}(b)+\frac{1}{2}\mathcal{C}_{1}^{(b)}(b)\right]f_{0},$
(35)
$\mathcal{K}\left[f_{1}^{(1y)}(\mathbf{b})\right]+\mathcal{C}_{0}(b)f_{1}^{(1y)}(\mathbf{b})=-\left[\mathcal{C}_{1}^{(a)}(b)-\frac{1}{2}\mathcal{C}_{1}^{(b)}(b)\right]f_{0}$
(36)

with

$\mathcal{K}\left[\mathbf{f}(\mathbf{b})\right]=-\frac{ik}{2p(p+k)}\left[\frac{d^{2}}{db^{2}}+\frac{1}{b}\frac{d}{db}-\frac{1}{b^{2}}-m_{\infty}^{2}\right]\mathbf{f}(\mathbf{b}).$
(37)

(The differential equation for $f_{1}^{(3)}$ is given in App. C.)

Our goal is to evaluate

$A_{z}=2\mathrm{Im}\,\partial_{b_{z}}f_{z}(\mathbf{b})\bigg{\rvert}_{b=0}=2\mathrm{Im}\,\frac{\mathbf{\widehat{b}}\cdot\mathbf{f}_{0}+\xi\,f_{1}^{(1z)}}{b}\bigg{\rvert}_{b=0}$
(38)

and

$A_{y}=2\mathrm{Im}\,\partial_{b_{y}}f_{y}(\mathbf{b})\bigg{\rvert}_{b=0}=2\mathrm{Im}\,\frac{\mathbf{\widehat{b}}\cdot\mathbf{f}_{0}+\xi\,f_{1}^{(1y)}}{b}\bigg{\rvert}_{b=0}$
(39)

which were defined in Eqs. (II) and (II). Thus, we only need to know
$f_{1}^{(1z)}(b)/b$ and $f_{1}^{(1y)}(b)/b$ in the limit $b\rightarrow 0$
where it must be finite. This gives the boundary condition that This is done
by demanding that $f_{1}^{(1z)}$ and $f_{1}^{(1y)}$ vanish at $b=0$. The other
boundary condition is that the functions vanish as $b\rightarrow\infty$ as can
be seen from Eq. (23).

To evaluate $f_{1}^{(1z)}(b)/b$ and $f_{1}^{(1y)}(b)/b$ at $b=0$, we demand
that the functions $f_{1}^{(1z)}$ and $f_{1}^{(1y)}$ vanish at very large $b$
and evolve the functions numerically to small $b$ using Eqs. (35) and (36). In
practice, this means that we start the evolution at a large but finite value
of $b$ where $f_{1}^{(1z)}$ and $f_{1}^{(1y)}$ are initialized to a small
value. A typical numerical solution for $f_{1}^{1z}$ and $f_{1}^{1y}$ then
blows up as evolved towards $b\rightarrow 0$. We must then extract the finite
part of our numerical solution. This is done by matching with known, analytic
solutions of the differential equations in the small $b$ limit.

For instance, focusing on Eq. (35), we call the particular solution $w(b)$ and
the two independent solutions of the homogeneous equation $w_{1}(b)$ and
$w_{2}(b)$. These are known analytically at small $b$, see Appendix C. We can
write our numerical solution in full generality at small $b$ as

$f_{1}^{(1z)}(b)=w(b)+\alpha_{1}w_{1}(b)+\alpha_{2}w_{2}(b)$ (40)

where $\alpha_{1}$ and $\alpha_{2}$ are found numerically. To extract from
this a solution with the right behaviour as $b\rightarrow 0$, one must in
essence subtract a linear combination of $w_{1}$ and $w_{2}$ which satisfies
the boundary condition at $b\rightarrow\infty$. Then one is left with a
solution which satisfies boundary conditions both at $b=0$ and
$b\rightarrow\infty$ and which gives $f_{1}^{(1z)}(b)/b$ at $b=0$. This
procedure is explained in further detail in Appendix C, see also [53, 48, 54]
for earlier work in the isotropic case. A major difference with the isotropic
case is that Eqs. (35) and (36) have a non-trivial right hand side which
complicates the matching procedure. For instance, one must find an analytic
solution $w(b)$ of the full differential equations at small $b$, including the
right hand side. Furthermore, cancellation errors between solutions of the
full differential equation and the homogeneous equation must carefully be
avoided to get reliable results, see Appendix C.

##  V Results

Figure 5: Spectrum of photons coming from bremsstrahlung and pair annihilation
in a plasma at effective temperature $\Lambda$. The anisotropic plasma has
$\xi=1.0$. Figure 6: Degree of polarization of photons emitted from
bremsstrahlung and pair annihilation in an anisotropic plasma with $\xi=1.0$.
The quantity $R$ is defined in Eq. (43).

Figs. 5 and 6 are the main results of this work. They show the rate of photon
production through bremsstrahlung and pair-annihilation in an anisotropic
quark-gluon plasma with fixed anisotropy $\xi$. The collision kernel is given
by Eq. (III) while the momentum distribution of medium quarks is

$f(\mathbf{p})=\frac{\sqrt{1+\xi}}{e^{\sqrt{p^{2}+\xi p_{z}^{2}}/\Lambda}+1}$
(41)

where $\Lambda$ can be seen as an effective temperature. Fig. 5 shows the
total rate for producing photons at mid-rapidity and with momentum $k$, i.e.

$k\frac{d\Gamma}{d^{3}\mathbf{k}}=k\frac{d\Gamma_{z}}{d^{3}\mathbf{k}}+k\frac{d\Gamma_{y}}{d^{3}\mathbf{k}},$
(42)

where $k\frac{d\Gamma_{z}}{d^{3}\mathbf{k}}$ is the rate of producing photons
polarized along the beam axis and $k\frac{d\Gamma_{y}}{d^{3}\mathbf{k}}$ is
the rate for photons polarized orthogonal to the beam axis and to the photon
momentum. Fig. 6 shows the degree of polarization at different momenta defined
as

$R=\frac{k\frac{d\Gamma_{z}}{d^{3}\mathbf{k}}-k\frac{d\Gamma_{y}}{d^{3}\mathbf{k}}}{k\frac{d\Gamma_{z}}{d^{3}\mathbf{k}}+k\frac{d\Gamma_{y}}{d^{3}\mathbf{k}}}.$
(43)

These quantities are shown for three values of anisotropy parameter,
$\xi=0.3$, $\xi=0.6$ and $\xi=0.9$, which correspond to pressure anisotropy of
$P_{L}/P_{T}\approx 0.81$, $P_{L}/P_{T}\approx 0.68$ and $P_{L}/P_{T}\approx
0.57$ respectively. These are rather moderate values of pressure anisotropy
which can be found in the early and intermediate stages of heavy-ion
collisions.

The spectrum in an anisotropic medium is higher than that in an equilibrium
medium at the same effective temperature $\Lambda$, as can be seen in Fig. 5.
This is due to the factor $\sqrt{1+\xi}$ in the momentum distribution in Eq.
(9) which increases the number of quarks with $p_{z}=0$ which can emit photons
at midrapidity. This effect is partially compensated by anisotropic
corrections which reduce the collision kernel
$\mathcal{C}(\mathbf{q}_{\perp})$ meaning that a given quark receives less
momentum broadening, reducing the rate at which it emits photons collinearly.

More interesting is the polarization $R$ as a function of momentum $k$. As
seen in Fig. 6, the polarization has different signs for lower and higher
values of the photon energy $k$: it is along the beam axis for lower values of
$k$ while it is orthogonal to the beam axis at higher values of $k$. This owes
to the interplay between bremsstrahlung and pair annihilation. As shown in
Appendix A, because of the different polarized splitting functions,
bremsstrahlung tends to give photons polarized along the $z$-axis, while pair
annihilation tends to give polarization along the $x$-axis, see Fig. 1. As
bremsstrahlung is suppressed at high $k$ (there are few medium quarks with
energy higher than $k$) $x$-polarization dominates in that regime. On the
contrary, pair annihilation is suppressed at low $k$ since the number of
quarks with energy less than $k/2$ is phase space suppressed. This gives
$z$-polarized photons in that regime.

Despite the complicated dependence of polarization on photon momentum,
polarization along the beam axis dominates. This is simply because there are
many more photons at lower $k$ and thus their polarization is dominant.
Furthermore, in [25] it was shown that photons from two-to-two scattering,
which are equally important as bremsstrahlung and pair annihilation photons,
are also predominantly polarized along the beam axis, with an even greater
magnitude of polarization. Thus a definite prediction of our work is that
medium photons are polarized along the beam axis.

To make contact with potential experiments on photon polarization more work is
needed. Firstly, all photon sources that cannot be subtracted in experiments
need to be included, such as prompt photons and photons from the hadronic
stage. Secondly, calculation of the rate of photon production in an
anisotropic medium need to be folded with hydrodynamic or kinetic theory
simulations of the medium to get a realistic evolution of the medium
anisotropy.

##  VI Conclusion

In this work, we have calculated for the first time the degree of polarization
for photons emitted in bremsstrahlung and off-shell pair annihilation
processes in a hot medium consisting of quarks and gluons. Our evaluation
includes the LPM regime and is at complete leading order in the strong
coupling. The polarization of the real photons originating from bremsstrahlung
and annihilation processes depends on the anisotropy of the original parton
distribution, and therefore the polarization can instruct us on the dynamics
in an environment that is not accessible to the vast majority of probes and
observables measured in relativistic heavy-ion collisions. Specifically, it
gives a measure of the pressure anisotropy at early times.

We trust that the methods and techniques developed and used here will be
useful in the evaluations of polarization signatures of real and virtual
photons, evaluated with scattering kernels for momentum broadening derived
from microscopic theories and using time-evolution models based in QCD.

###### Acknowledgements.

C. G. acknowledges the support of the Natural Sciences and Engineering
Research Council of Canada (NSERC), SAPIN-2020-00048.

##  Appendix A Derivation of rate equations

We wish to show that the total rate for polarized photon production through
bremsstrahlung and pair annihilation can be written as

$\displaystyle k\frac{d\Gamma_{z}}{d^{3}\mathbf{k}}$
$\displaystyle=\frac{6\alpha_{EM}\sum_{s}q_{s}^{2}}{(2\pi)^{3}}\int_{-\infty}^{\infty}dp^{z}\;F(\mathbf{k}+\mathbf{p})\left[1-F(\mathbf{p})\right]$
(44)
$\displaystyle\times\frac{1}{2}\frac{1}{4(p^{x})^{2}(p^{x}+k)^{2}}\left[(2p^{x}+k)^{2}A_{z}+k^{2}A_{y}\right]$

and

$\displaystyle k\frac{d\Gamma_{y}}{d^{3}\mathbf{k}}$
$\displaystyle=\frac{6\alpha_{EM}\sum_{s}q_{s}^{2}}{(2\pi)^{3}}\int_{-\infty}^{\infty}dp^{z}\;F(\mathbf{k}+\mathbf{p})\left[1-F(\mathbf{p})\right]$
(45)
$\displaystyle\times\frac{1}{2}\frac{1}{4(p^{x})^{2}(p^{x}+k)^{2}}\left[k^{2}A_{z}+(2p^{x}+k)^{2}A_{y}\right]$

where $A_{z}$ and $A_{y}$ are defined in Eq. (II). Here the momentum
distributions are contained in

$F(p^{x})=\theta(p^{x})f(\mathbf{p})+\theta(-p^{x})\left(1-f(-\mathbf{p})\right).$
(46)

Looking at the momentum factors in Eqs. (44) and (45), we see that
bremsstrahlung off a quark or an antiquark corresponds to $p^{x}>0$ and
$p^{x}<-k$. These cases can be rewritten to give Eqs. (II) and (II).
Furthermore, $-k<p^{x}<0$ corresponds to quark-antiquark pair annihilation. As
the quark moves in the same direction as the photon we can set e.g.
$f(\mathbf{p})=f(p^{x})$.

The derivation of polarized photon emission in Eqs. (44) and (45) is similar
to that for unpolarized emission found in [28, 30, 23]. The polarized rate
comes from evaluating the diagram in Fig. 8. The gluon ladders which represent
soft kicks from the medium are evaluated in the same way as for unpolarized
photon emission. Physically, this is because the soft kicks do not have enough
energy to flip the helicity of quarks. The gluon ladders are summed to all
orders to give the integral equation in Eq. (6). On the contrary, one must
keep track of polarization at the hard emission vertices to evaluate the
polarized emission rate.

For a bare quark loop without soft gluon kicks, Fig. 7, the hard emission
vertices for a $z-$polarized photon take the form

$H^{zz}:=\epsilon^{z}_{\mu}\epsilon^{z\,*}_{\nu}\mathrm{Tr}\,\left[\gamma^{\mu}(\not{K}+\not{P})\gamma^{\nu}\not{P}\right].$
(47)

Here $\epsilon^{z}_{\mu}$ is the photon polarization tensor for
$z$-polarization. This trace is most easily evaluated by using that e.g.

$\not{P}=\sum_{s}u^{s}(\mathbf{p})\overline{u}^{s}(\mathbf{p})$ (48)

where the sum is over spin states and

$u^{s}(\mathbf{p})=\sqrt{2p}\begin{bmatrix}\frac{1-\mathbf{\sigma}\cdot\widehat{\mathbf{p}}}{2}\xi^{s}\\\
\frac{1+\mathbf{\sigma}\cdot\widehat{\mathbf{p}}}{2}\xi^{s}\end{bmatrix}$ (49)

are helicity states of quarks [55]. Here $\xi^{s}$ form a basis for two-
component spin states. One can then show by an explicit calculation that

$\displaystyle H^{zz}$
$\displaystyle=\sum_{s,t}\left[\gamma^{z}u^{s}(\mathbf{p}+\mathbf{k})\overline{u}^{s}(\mathbf{p}+\mathbf{k})\gamma^{z}u^{t}(\mathbf{p})\overline{u}^{t}(\mathbf{p})\right]$
(50)
$\displaystyle=\frac{8p(p+k)}{p^{2}(p+k)^{2}}\left[(2p+k)^{2}p_{\perp\,z}^{2}+k^{2}p_{\perp\,y}^{2}\right].$

Similarly, the hard emission vertices for emission of $y$-polarized photons
are

$H^{yy}=\frac{8p(p+k)}{p^{2}(p+k)^{2}}\left[(2p+k)^{2}p_{\perp\,y}^{2}+k^{2}p_{\perp\,z}^{2}\right].$
(51)

This is the same as Eq. (50), except that $p_{\perp\,z}^{2}$ and
$p_{\perp\,y}^{2}$ have been interchanged.

To include soft gluon as in Fig. 8 one simply replaces one of the hard
vertices by the dressed vertex $\mathbf{f}(\mathbf{p}_{\perp})$ which includes
gluons rungs and obeys the integral equation in Eq. (6). This means that in
Eqs. (50) and (51) one replaces $p_{\perp\,z}^{2}\longrightarrow
p_{\perp\,z}f_{z}$ and $p_{\perp\,y}^{2}\longrightarrow p_{\perp\,y}f_{y}$.
This reproduces Eqs. (44) and (45).

Figure 7: Definition of momenta in photon emission. Figure 8: Diagram for
medium-induced photon emission. The gluon rungs are responsible for momentum
broadening.

##  Appendix B The retarded gluon propagator at small anisotropy

In a system whose hard quasiparticles have the momentum distribution in Eq.
(9), the retarded propagator for soft gluons is [45, 47]

$\displaystyle D_{\mathrm{ret}}^{\mu\nu}(Q)$
$\displaystyle=\left(P_{T}^{\mu\nu}-C^{\mu\nu}\right)D_{\mathrm{ret}}^{B}$
$\displaystyle+\left[\left(Q^{2}-\Pi_{c}\right)P_{L}^{\mu\nu}+\left(Q^{2}-\Pi_{L}\right)C^{\mu\nu}+\Pi_{d}D^{\mu\nu}\right]D_{\mathrm{ret}}^{A}$
(52)

where

$D_{\mathrm{ret}}^{A}=\frac{1}{\left(Q^{2}-\Pi_{L}\right)\left(Q^{2}-\Pi_{c}\right)-\frac{Q^{2}\mathbf{q}^{2}}{q_{0}^{2}}\Pi_{d}^{2}}$
(53)

and

$D_{\mathrm{ret}}^{B}=\frac{1}{Q^{2}-\Pi_{e}}.$ (54)

The tensors $P_{T}^{\mu\nu}$ and $P_{L}^{\mu\nu}$ are the same as in thermal
equilibrium while $C^{\mu\nu}$ and $D^{\mu\nu}$ are new tensors that depend on
the anisotropy vector $\mathbf{n}$. The self-energy components $\Pi_{e}$,
$\Pi_{L}$, $\Pi_{T}$ and $\Pi_{d}$ are given in [47] which contains further
details.

As $\Pi_{d}=\mathcal{O}(\xi)$ we can approximate

$D_{\mathrm{ret}}^{A}\approx\frac{1}{\left(Q^{2}-\Pi_{L}\right)\left(Q^{2}-\Pi_{c}\right)}$
(55)

up to order $\mathcal{O}(\xi)$. We furthermore ignore $\mathcal{O}(\xi)$
corrections in the numerator as they do not describe anisotropic corrections
to screening which are the non-equilibrium corrections we focus on. Then

$\displaystyle D_{\mathrm{ret}}^{\mu\nu}(Q)$
$\displaystyle\approx\left(P_{T}^{\mu\nu}-C^{\mu\nu}\right)D_{\mathrm{ret}}^{B}$
$\displaystyle\qquad+\left[\left(Q^{2}-\Pi_{c}\right)P_{L}^{\mu\nu}+\left(Q^{2}-\Pi_{L}\right)C^{\mu\nu}\right]D_{\mathrm{ret}}^{A}$
$\displaystyle=\frac{P_{T}^{\mu\nu}-C^{\mu\nu}}{Q^{2}-\Pi_{e}}+\frac{P_{L}^{\mu\nu}}{Q^{2}-\Pi_{L}}+\frac{C^{\mu\nu}}{Q^{2}-\Pi_{c}}$
(56)

The terms with the tensor $C^{\mu\nu}$ go like

$C^{\mu\nu}\left[\frac{-1}{Q^{2}-\Pi_{e}}+\frac{1}{Q^{2}-\Pi_{c}}\right]=C^{\mu\nu}\frac{\Pi_{c}-\Pi_{e}}{(Q^{2}-\Pi_{e})(Q^{2}-\Pi_{c})}$
(57)

which can be ignored as $\Pi_{c}-\Pi_{e}=\mathcal{O}(\xi)$ in the numerator.
We are thus left with

$D_{\mathrm{ret}}^{\mu\nu}(Q)\approx\frac{P_{T}^{\mu\nu}}{Q^{2}-\Pi_{e}}+\frac{P_{L}^{\mu\nu}}{Q^{2}-\Pi_{L}}.$
(58)

at small anisotropy where we only include anisotropic corrections in the
denominators. The tensors $P_{L}$ and $P_{T}$ are the same as in equilibrium
while $\Pi_{e}$ and $\Pi_{L}$ have anisotropic corrections. (We call
$\Pi_{e}=\Pi_{T}$ in the main text of the paper.)

##  Appendix C Details of numerical method

In this Appendix we discuss how to solve Eqs. (35) and (36) numerically.
Unlike the isotropic equation, Eq. (31), the anisotropic equation has a non-
vanishing source term on the right hand side for all values of $b$. One thus
needs a different numerical solution method than that developed in [53, 48,
54] for the equilibrium case. We note that the differential equation for the
function $f_{1}^{(3)}$ is

$\displaystyle\mathcal{K}\left[f_{1}^{(3)}(\mathbf{b})\right]$
$\displaystyle+\frac{ik}{2p(p+k)}\frac{8}{b^{2}}f_{1}^{(3)}(\mathbf{b})+\mathcal{C}_{0}(b)f_{1}^{(3)}(\mathbf{b})$
$\displaystyle=-\frac{1}{2}\mathcal{C}_{1}^{(b)}(b)f_{0}$ (59)

but we will not discuss this function further as it is not needed to evaluate
$A_{z}$ and $A_{y}$ in Eqs. (38) and (39)

For concreteness, we focus on solving Eq. (35). Defining a scaled function

$G=\frac{\pi}{2}\frac{k}{p(p+k)m_{\infty}^{2}}\frac{f_{1}^{(1z)}(b)}{b},$ (60)

as well as scaled collision kernels
$\overline{\mathcal{C}}(t)=\frac{2p(p+k)}{k}\frac{1}{m_{\infty}^{2}}\mathcal{C}(b)$
and variable $t=m_{\infty}b$, this equation becomes

$\displaystyle-i$
$\displaystyle\left[\frac{d^{2}G}{dt^{2}}+\frac{3}{t}\frac{dG}{dt}-G\right]+\overline{\mathcal{C}}_{0}(t)G$
$\displaystyle=-\left[\overline{\mathcal{C}}_{1}^{(a)}(t)+\frac{1}{2}\overline{\mathcal{C}}_{1}^{(b)}(t)\right]\overline{f}_{0}(t)$
(61)

where
$\overline{f}_{0}(t)=\frac{\pi}{2}\frac{k}{p(p+k)m_{\infty}^{2}}\mathbf{b}\cdot\mathbf{f}_{0}/b^{2}$.
As shown in [48], we can write
$\overline{f}_{0}(t)=K_{1}(t)/t+\overline{f}_{0}^{\mathrm{rest}}(t)$ where
$\overline{f}_{0}^{\mathrm{rest}}(t)$ is function that is finite in the limit
$t\rightarrow 0$ and which we know numerically using the methods of [53, 48].

We need to solve Eq. (C), imposing the boundary conditions that
$G(t)\rightarrow 0$ as $t\rightarrow\infty$, as well as that $G(t)$ is finite
as $t\rightarrow 0$. These boundary conditions are difficult to satisfy
simultaneously for a numerical solution. Instead we find a numerical solutions
$g(t)$ of Eq. (C) with $g(t\rightarrow\infty)=0$ and a solution of the
homogeneous equation without the source term that also satisfies
$g_{0}(t\rightarrow\infty)=0$. In general, both $g(t)$ and $g_{0}(t)$ blow up
as $t\rightarrow 0$. However, we know that the solution we are seeking can be
written as

$G(t)=g(t)+Ag_{0}(t)$ (62)

where $A$ is chosen so that $G(0)$ is finite.

We can find an explicit expression of $A$ by using analytic solutions of Eq.
(C) for $t\ll 1$. In that limit the right hand side is

$\displaystyle\left[\overline{\mathcal{C}}_{1}^{(a)}(t)+\frac{1}{2}\overline{\mathcal{C}}_{1}^{(b)}(t)\right]\overline{f}_{0}(t)$
$\displaystyle\approx\left[\overline{\mathcal{C}}_{1}^{(a)}(t)+\frac{1}{2}\overline{\mathcal{C}}_{1}^{(b)}(t)\right]K_{1}(t)/t\approx
a+b\log t$ (63)

where $a$ and $b$ are constants that depend on the momenta $k$ and $p$ as well
as the masses $m_{D}^{2}$ and $m_{\infty}^{2}$. The differential equation
becomes

$-i\left[\frac{d^{2}g}{dt^{2}}+\frac{3}{t}\frac{dg}{dt}-g\right]+\overline{\mathcal{C}}_{0}(t)g=a+b\log
t$ (64)

which has general solution

$g(t)=w(t)+\alpha_{1}w_{1}(t)+\alpha_{2}w_{2}(t)$ (65)

where

$w(t)=-ia-ib(\log t+\frac{2}{t^{2}})$ (66)

is an exact particular solution and $w_{2}(t)=2J_{1}(it)/(it)$ and
$w_{1}(t)=\frac{\pi}{2}Y_{1}(-it)/(-it)-\frac{1}{4}\left(2\gamma_{E}-2\log
2+i\pi\right)w_{2}(t)$ are solutions of the homogenous equation such that
$w_{1}(t)=1/t^{2}+\frac{1}{2}\log t+\mathcal{O}(1)$ and
$w_{2}(t)=\mathcal{O}(1)$ for small $t$. Similarly, the homogeneous equation
can be written as

$g_{0}(t)=\beta_{1}w_{1}(t)+\beta_{2}w_{2}(t).$ (67)

The coefficients $\alpha_{1}$, $\alpha_{2}$, $\beta_{1}$ and $\beta_{2}$ can
be found from Eqs. (65) and (67) by equating the numerical solutions $g(t)$
and $g_{0}(t)$ and their derivatives with the analytic functions at some small
value $t_{\mathrm{min}}\ll 1$.

The small $t$ behaviour of Eqs. (65) and (67) shows that

$A=\frac{2ib-\alpha_{1}}{\beta_{1}}$ (68)

which makes

$G(0)=-ia-\frac{1}{2}ib+\alpha_{2}+\frac{(2ib-\alpha_{1})\beta_{2}}{\beta_{1}}$
(69)

finite. This is the expression that we need in order to evaluate Eq. (38). All
quantities are known for numerical solutions $g(t)$ and $g_{0}(t)$.

Eq. (69) suffers from numerical cancellation errors in the terms
$\alpha_{2}-\alpha_{1}\beta_{2}/\beta_{1}$. It is thus better to rewrite these
terms using Eqs. (65) and (67) which shows that

$\alpha_{2}-\frac{\alpha_{1}\beta_{2}}{\beta_{1}}=\frac{g^{\prime}g_{0}-gg_{0}^{\prime}+g_{0}^{\prime}w-g_{0}w^{\prime}}{g_{0}w_{2}^{\prime}-g_{0}^{\prime}w_{2}}.$
(70)

where all quantities are evaluated at $t_{\mathrm{min}}\ll 1$. The culprit
behind cancellation errors is the term $g^{\prime}g_{0}-gg_{0}^{\prime}$. It
can be evaluated more precisely by noting that

$\mathcal{G}_{W}(t)=g^{\prime}(t)g_{0}(t)-g(t)g_{0}^{\prime}(t).$ (71)

solves the equation

$\frac{d(t^{3}\mathcal{G}_{W})}{dt}=-it^{3}\left[\overline{\mathcal{C}}_{1}^{(a)}(t)+\frac{1}{2}\overline{\mathcal{C}}_{1}^{(b)}(t)\right]\overline{f}_{0}(t)g_{0}$
(72)

which can be integrated to give $\mathcal{G}_{W}(t_{\mathrm{min}})$ and thus a
reliable value of $G(0)$.

## References

  * Aarts  _et al._ [2023] G. Aarts _et al._ , Phase Transitions in Particle Physics – Results and Perspectives from Lattice Quantum Chromo-Dynamics (2023) arXiv:2301.04382 [hep-lat] .
  * Aoki  _et al._ [2006] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, The Order of the quantum chromodynamics transition predicted by the standard model of particle physics, Nature 443, 675 (2006), arXiv:hep-lat/0611014 .
  * Busza  _et al._ [2018] W. Busza, K. Rajagopal, and W. van der Schee, Heavy Ion Collisions: The Big Picture, and the Big Questions, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Strickland [2019] M. Strickland, Small system studies: A theory overview, Nucl. Phys. A 982, 92 (2019), arXiv:1807.07191 [nucl-th] .
  * Snellings [2007] R. Snellings (STAR, ALICE), Anisotropic flow from RHIC to the LHC, Eur. Phys. J. C 49, 87 (2007), arXiv:nucl-ex/0610010 .
  * Connors  _et al._ [2018] M. Connors, C. Nattrass, R. Reed, and S. Salur, Jet measurements in heavy ion physics, Rev. Mod. Phys. 90, 025005 (2018), arXiv:1705.01974 [nucl-ex] .
  * David [2020] G. David, Direct real photons in relativistic heavy ion collisions, Rept. Prog. Phys. 83, 046301 (2020), arXiv:1907.08893 [nucl-ex] .
  * Salabura [2008] P. Salabura, Dilepton production in heavy ion collisions, Acta Phys. Polon. B 39, 307 (2008).
  * Heffernan  _et al._ [2023] M. R. Heffernan, C. Gale, S. Jeon, and J.-F. Paquet, Bayesian quantification of strongly-interacting matter with color glass condensate initial conditions, (2023), arXiv:2302.09478 [nucl-th] .
  * Everett  _et al._ [2021a] D. Everett _et al._ (JETSCAPE), Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging, Phys. Rev. Lett. 126, 242301 (2021a), arXiv:2010.03928 [hep-ph] .
  * Everett  _et al._ [2021b] D. Everett _et al._ (JETSCAPE), Multisystem Bayesian constraints on the transport coefficients of QCD matter, Phys. Rev. C 103, 054904 (2021b), arXiv:2011.01430 [hep-ph] .
  * Bernhard  _et al._ [2019] J. E. Bernhard, J. S. Moreland, and S. A. Bass, Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma, Nature Phys. 15, 1113 (2019).
  * Moreland  _et al._ [2015] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions, Phys. Rev. C92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Schenke  _et al._ [2012] B. Schenke, P. Tribedy, and R. Venugopalan, Fluctuating Glasma initial conditions and flow in heavy ion collisions, Phys. Rev. Lett. 108, 252301 (2012), arXiv:1202.6646 [nucl-th] .
  * Schenke  _et al._ [2010] B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Schenke  _et al._ [2011] B. Schenke, S. Jeon, and C. Gale, Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics, Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Paquet  _et al._ [2016] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Production of photons in relativistic heavy-ion collisions, Phys. Rev. C 93, 044906 (2016), arXiv:1509.06738 [hep-ph] .
  * Shen  _et al._ [2016] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, The iEBE-VISHNU code package for relativistic heavy-ion collisions, Comput. Phys. Commun. 199, 61 (2016), arXiv:1409.8164 [nucl-th] .
  * Bleicher  _et al._ [1999] M. Bleicher _et al._ , Relativistic hadron hadron collisions in the ultrarelativistic quantum molecular dynamics model, J. Phys. G 25, 1859 (1999), arXiv:hep-ph/9909407 .
  * Weil  _et al._ [2016] J. Weil, V. Steinberg, J. Staudenmaier, L. G. Pang, D. Oliinychenko, J. Mohs, M. Kretz, T. Kehrenberg, A. Goldschmidt, B. Bäuchle, J. Auvinen, M. Attems, and H. Petersen, Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions, Phys. Rev. C 94, 054905 (2016).
  * Cao  _et al._ [2021] S. Cao _et al._ (JETSCAPE), Determining the jet transport coefficient q from inclusive hadron suppression measurements using Bayesian parameter estimation, Phys. Rev. C 104, 024905 (2021), arXiv:2102.11337 [nucl-th] .
  * Gale [2010] C. Gale, Photon Production in Hot and Dense Strongly Interacting Matter, Landolt-Bornstein 23, 445 (2010), arXiv:0904.2184 [hep-ph] .
  * Hauksson  _et al._ [2018] S. Hauksson, S. Jeon, and C. Gale, Photon emission from quark-gluon plasma out of equilibrium, Phys. Rev. C 97, 014901 (2018), arXiv:1709.03598 [nucl-th] .
  * Hauksson  _et al._ [2021] S. Hauksson, S. Jeon, and C. Gale, Probes of the quark-gluon plasma and plasma instabilities, Phys. Rev. C 103, 064904 (2021), arXiv:2012.03640 [hep-ph] .
  * Baym and Hatsuda [2015] G. Baym and T. Hatsuda, Polarization of Direct Photons from Gluon Anisotropy in Ultrarelativistic Heavy Ion Collisions, PTEP 2015, 031D01 (2015), arXiv:1405.1376 [nucl-th] .
  * Schenke and Strickland [2007] B. Schenke and M. Strickland, Photon production from an anisotropic quark-gluon plasma, Phys. Rev. D 76, 025023 (2007), arXiv:hep-ph/0611332 .
  * Bhattacharya  _et al._ [2016] L. Bhattacharya, R. Ryblewski, and M. Strickland, Photon production from a nonequilibrium quark-gluon plasma, Phys. Rev. D 93, 065005 (2016), arXiv:1507.06605 [hep-ph] .
  * Arnold  _et al._ [2001a] P. B. Arnold, G. D. Moore, and L. G. Yaffe, Photon emission from ultrarelativistic plasmas, JHEP 11, 057, arXiv:hep-ph/0109064 .
  * Arnold  _et al._ [2001b] P. B. Arnold, G. D. Moore, and L. G. Yaffe, Photon emission from quark gluon plasma: Complete leading order results, JHEP 12, 009, arXiv:hep-ph/0111107 .
  * Arnold  _et al._ [2002] P. B. Arnold, G. D. Moore, and L. G. Yaffe, Photon and gluon emission in relativistic plasmas, JHEP 06, 030, arXiv:hep-ph/0204343 .
  * Ghiglieri  _et al._ [2013] J. Ghiglieri, J. Hong, A. Kurkela, E. Lu, G. D. Moore, and D. Teaney, Next-to-leading order thermal photon production in a weakly coupled quark-gluon plasma, JHEP 05, 010, arXiv:1302.5970 [hep-ph] .
  * Baier  _et al._ [1992] R. Baier, H. Nakkagawa, A. Niegawa, and K. Redlich, Production rate of hard thermal photons and screening of quark mass singularity, Z. Phys. C 53, 433 (1992).
  * Kapusta  _et al._ [1991] J. I. Kapusta, P. Lichard, and D. Seibert, High-energy photons from quark - gluon plasma versus hot hadronic gas, Phys. Rev. D 44, 2774 (1991), [Erratum: Phys.Rev.D 47, 4171 (1993)].
  * Yee [2013] H.-U. Yee, Flows and polarization of early photons with magnetic field at strong coupling, Phys. Rev. D 88, 026001 (2013), arXiv:1303.3571 [nucl-th] .
  * Wu and Yang [2013] S.-Y. Wu and D.-L. Yang, Holographic Photon Production with Magnetic Field in Anisotropic Plasmas, JHEP 08, 032, arXiv:1305.5509 [hep-th] .
  * Ávila  _et al._ [2021] D. Ávila, T. Monroy, F. Nettel, and L. Patiño, Emission of linearly polarized photons in a strongly coupled magnetized plasma from the gauge/gravity correspondence, Phys. Lett. B 817, 136287 (2021), arXiv:2101.08802 [hep-th] .
  * Dong and Lin [2022] L. Dong and S. Lin, Dilepton helical production in a vortical quark-gluon plasma, Eur. Phys. J. A 58, 176 (2022), arXiv:2112.07153 [hep-ph] .
  * Mamo and Yee [2013] K. A. Mamo and H.-U. Yee, Spin polarized photons and dileptons from axially charged plasma, Phys. Rev. D 88, 114029 (2013), arXiv:1307.8099 [nucl-th] .
  * Ipp  _et al._ [2008] A. Ipp, A. Di Piazza, J. Evers, and C. H. Keitel, Photon polarization as a probe for quark-gluon plasma dynamics, Phys. Lett. B 666, 315 (2008), arXiv:0710.5700 [hep-ph] .
  * Baym  _et al._ [2017] G. Baym, T. Hatsuda, and M. Strickland, Virtual photon polarization in ultrarelativistic heavy-ion collisions, Phys. Rev. C 95, 044907 (2017), arXiv:1702.05906 [nucl-th] .
  * Shuryak [2012] E. Shuryak, Monitoring parton equilibration in heavy ion collisions via dilepton polarization, (2012), arXiv:1203.1012 [nucl-th] .
  * Speranza  _et al._ [2018] E. Speranza, A. Jaiswal, and B. Friman, Virtual photon polarization and dilepton anisotropy in relativistic nucleus–nucleus collisions, Phys. Lett. B 782, 395 (2018), arXiv:1802.02479 [hep-ph] .
  * Aurenche  _et al._ [1998] P. Aurenche, F. Gelis, R. Kobes, and H. Zaraket, Bremsstrahlung and photon production in thermal QCD, Phys. Rev. D 58, 085003 (1998), arXiv:hep-ph/9804224 .
  * Ellis  _et al._ [2011] R. K. Ellis, W. J. Stirling, and B. R. Webber,  _QCD and collider physics_, Vol. 8 (Cambridge University Press, 2011).
  * Romatschke and Strickland [2003] P. Romatschke and M. Strickland, Collective modes of an anisotropic quark gluon plasma, Phys. Rev. D 68, 036004 (2003), arXiv:hep-ph/0304092 .
  * Martinez and Strickland [2010] M. Martinez and M. Strickland, Dissipative Dynamics of Highly Anisotropic Systems, Nucl. Phys. A 848, 183 (2010), arXiv:1007.0889 [nucl-th] .
  * Hauksson  _et al._ [2022] S. Hauksson, S. Jeon, and C. Gale, Momentum broadening of energetic partons in an anisotropic plasma, Phys. Rev. C 105, 014914 (2022), arXiv:2109.04575 [hep-ph] .
  * Aurenche  _et al._ [2002a] P. Aurenche, F. Gelis, and H. Zaraket, A Simple sum rule for the thermal gluon spectral function and applications, JHEP 05, 043, arXiv:hep-ph/0204146 .
  * Caron-Huot [2009] S. Caron-Huot, O(g) plasma effects in jet quenching, Phys. Rev. D 79, 065039 (2009), arXiv:0811.1603 [hep-ph] .
  * Panero  _et al._ [2014] M. Panero, K. Rummukainen, and A. Schäfer, Lattice Study of the Jet Quenching Parameter, Phys. Rev. Lett. 112, 162001 (2014), arXiv:1307.5850 [hep-ph] .
  * Moore  _et al._ [2021] G. D. Moore, S. Schlichting, N. Schlusser, and I. Soudi, Non-perturbative determination of collisional broadening and medium induced radiation in QCD plasmas, JHEP 10, 059, arXiv:2105.01679 [hep-ph] .
  * Mrowczynski  _et al._ [2017] S. Mrowczynski, B. Schenke, and M. Strickland, Color instabilities in the quark–gluon plasma, Phys. Rept. 682, 1 (2017), arXiv:1603.08946 [hep-ph] .
  * Aurenche  _et al._ [2002b] P. Aurenche, F. Gelis, G. D. Moore, and H. Zaraket, Landau-Pomeranchuk-Migdal resummation for dilepton production, JHEP 12, 006, arXiv:hep-ph/0211036 .
  * Jeon and Moore [2005] S. Jeon and G. D. Moore, Energy loss of leading partons in a thermal QCD medium, Phys. Rev. C 71, 034901 (2005), arXiv:hep-ph/0309332 .
  * Peskin and Schroeder [1995] M. E. Peskin and D. V. Schroeder,  _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995).


# Leveraging recent advances in Pre-Trained Language Models for Eye-Tracking
Prediction

Varun Madhavan
IIT Kharagpur
varun.m.iitkgp
@gmail.com
Aditya Girish Pawate
IIT Kharagpur
adityagirish
<EMAIL_ADDRESS>
Shraman Pal
IIT Kharagpur
shramanpal12345
@gmail.com
Abhranil Chandra
IIT Kharagpur
abhranil.chandra
@gmail.com

###### Abstract

Cognitively inspired Natural Language Processing uses human-derived behavioral
data like eye-tracking data, which reflect the semantic representations of
language in the human brain to augment the neural nets to solve a range of
tasks spanning syntax and semantics with the aim of teaching machines about
language processing mechanisms. In this paper, we use the ZuCo 1.0 and ZuCo
2.0 dataset containing the eye-gaze features to explore different linguistic
models to directly predict these gaze features for each word with respect to
its sentence. We tried different neural network models with the words as
inputs to predict the targets. And after lots of experimentation and feature
engineering finally devised a novel architecture consisting of RoBERTa Token
Classifier with a dense layer on top for language modeling and a stand-alone
model consisting of dense layers followed by a transformer layer for the extra
features we engineered. Finally, we took the mean of the outputs of both these
models to make the final predictions. We evaluated the models using mean
absolute error (MAE) and the R2 score for each target.

##  1 Introduction

When reading, we humans process language “automatically” without reflecting on
each step — we string words together into sentences, understand the meaning of
spoken and written ideas, and process language without thinking too much about
how the underlying cognitive process happens. This process generates cognitive
signals that could potentially facilitate natural language processing tasks.
Our understanding of the field of language processing is highly dependent on
accurately modeling eye-tracking features (example of cognitive data) as they
provide millisecond-accurate records of where humans look while reading, thus
shedding some light on human attention mechanisms and language comprehension
skills. Eye movement data is used extensively in various fields including but
not limited to Natural Language Processing and Computer Vision. The recent
introduction of standardized datasets has enabled the comparison of the
capabilities of different cognitive modeling and linguistically motivated
approaches to model and analyze human patterns of reading.

##  2 Problem Description

Feature Name |  Feature Description
---|---
nFix |  Number of fixations on the word
FFD |  Duration of the first fixation
TRT |  Total fixation duration, including regressions
GPT |  Sum of all prior fixations
fixProp |  Proportion of participants that fixated the current word
Table 2.1: Target Variables

Our brain processes language and generates cognitive processing data such as
gaze patterns and brain activity. These signals can be recorded while reading.
The problem is formulated as a regression task to predict the five eye-
tracking features for each token of a sentence in the test set. Each training
sample consists of a token in a sentence and the corresponding features. We
are required to predict the five features: nFix, FFD, TRT, GPT, and fixProp.
We were given a CSV file with sentences to train our model. Tokens in the
sentences were split in the same manner as they were presented to the
participants during the reading experiments. Hence, this did not necessarily
follow linguistically correct tokenization. Sentence endings were marked with
an <EOS> symbol. A prediction had to be made for each token for all the five
target variables. The evaluation metric used was the mean absolute error.

##  3 Dataset

We use the eye-tracking data of the Zurich Cognitive Language Processing
Corpus (ZuCo 1.0 and ZuCo 2.0) Hollenstein et al. (2020) recorded during
normal reading. The training and test data contains 800 and 191 sentences
respectively.In recent years, collecting these signals has become increasingly
easy and less expensive Papoutsaki et al. (2016); as a result, using cognitive
features to improve NLP tasks has become more popular. The data contains
features averaged over all readers and scaled in the range of 0 to 100,
facilitating evaluation via mean absolute average (MAE). The 5 eye-tracking
features included in the dataset are shown in 2.1. The heatmap of the
correlation matrix of the 5 target variables have been shown in Figure 1. It
is observed that apart from GPT all the features are correlated with each
other to a large extent.

Figure 3.1: Correlation Heatmap of Targets

##  4 Preprocessing

We first removed the <EOS> token at the end of each sentence. We looked at
similar word Clifton et al. (2007), Demberg and Keller (2008) and Just and
Carpenter (1980) to find features that we could add. As GPT consistently
performed the worst in all our initial experiments, we tried replacing GPT
with a new feature (TRT - GPT) during training. The results improved slightly
and we decided to keep this change for the remaining experiments. To
incorporate word level characteristics (i.e characteristics that belong to a
word; eg: word length) and syntactic information we introduced the following
features -

  * •

Number of characters in a word ($wordlen$)

  * •

Number of characters in the lemmatized word subtracted from the number of
characters in the word. Obtained using NLTK package Bird and Loper (2004)
($lemwordlen$)

  * •

If a word was a stopword or not, we appended 1 if stopword and -1
otherwise($stopword$)

  * •

If a word was a number, we appended 1 if it was a number and -1
otherwise($number$)

  * •

If a word was the at end of the sentence or not, we appended 1 if last word
and -1 otherwise ($endword$). Our hypothesis for adding the endword was that
it had a strong correlation with GPT as endwords had significantly higher GPT
values compared to the other words.

  * •

Part of Speech Tags for each word in a sentence ($postag$). We one hot encoded
the POS TAGS which gave us a sparse matrix of 20 features. The POS tags have a
high statistical correlation with all the target variables as they represent
the semantic structure of a sentence. Obtained using NLTK package.

  * •

TF-IDF value of each word by considering each sentence as a separate document
and each word as a separate term($tfidf$). For the tfidf score, words
unrecognized by the in-built tokenizer due to being single digit numbers,
having “-” etc were replaced by 0.01 if it was among “a”, “A” or “I” and 0.1
otherwise. This provides a hyperparameter(TF-IDF Error) that could be tuned
for better results. Obtained using sklearn package. Pedregosa et al. (2011)

Scatter plots of word length(number of characters) and the 5 target variables
have been shown. It is observed that there is quite a positive correlation
between them. we expected the last word of the sentence to have an impact on
the GPT score which is evident from the barplot as shown in figure 2f. The GPT
score for a endword is usually quite high compared to a non endword. We also
created barplots of stopword and the target variables. It is evident that for
stopwords the values of all the targets variables are comparitively lower
which points to stopword being a significant feature. The effect of adding
these features is further explored as an ablation study in section 6.3.

Figure 4.1: nFix vs
wordlen

Figure 4.2: FFD vs
wordlen

Figure 4.3: GPT vs
wordlen

Figure 4.4: TRT vs
wordlen

Figure 4.5: GPT vs
endword

Figure 4.6: nFix vs
stopword

Figure 4.7: FFD vs
stopword

Figure 4.8: GPT vs
stopword

Figure 4.9: TRT vs
stopword

Figure 4.10: fixProp vs
stopword

##  5 Model Description

The model has two major components, i.e. the Feature Model and the Language
Model, which encode the features and the sentence respectively. The resulting
encoded representations are mean pooled and passed to a fully connected layer
and sigmoid activation to output a (5x1) score vector for each token in the
sentence. The score vector represents the values of nFix, FFD, TRT, GPT and
fixProp respectively.

Figure 5.1: Feature Model

Figure 5.2: Language Head

Figure 5.3: Complete Model Flow

##  6 Experiments

###  6.1 Evaluation Metric

We utilised the R2 score metric instead of the Mean Absolute Error score
because R2 provides a better understanding of the model’s performance in
absence of a baseline score. As the R2 score has an upper bound of 1 it is
possible to understand how good the model is performing for a particular value
of R2 score.

$R^{2}(y,\hat{y})=1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}$

###  6.2 Models

We experimented with a variety of models, in particular recent advances in
attention-basedVaswani et al. (2017) pre-trained language models like RoBERTa
Liu et al. (2019), ELECTRA Clark et al. (2020), BERT Devlin et al. (2019),
DistilBERTSanh et al. (2020). Each language model was used to encode text as
described in section 5 and the rest of the pipeline was kept the same. We
compare this method with Hollenstein et al. (2019), which employed GloVe(200d)
Pennington et al. (2014) and a Bidirectional LSTM Hochreiter and Schmidhuber
(1997) to encode sentences. In addition to GloVe, we also tried the more
recent ConceptNet Numberbatch embeddings (300D) Speer et al. (2018) that have
been shown to outperform GloVe on a variety of tasks. The results for all
experiments are tabulated in Table 6.3.

###  6.3 Features

In addition to various architectures, we tried a number of feature engineering
techniques with inferred/transformed features as described in section 4. We
performed an ablation study with each feature to judge the importance/benefit
of each adding each of the new features.

We note that increasing the number of features does not always help, so
analysing the data we saw that most of the POS tags rarely occurred and there
were 6-7 frequently occurring tags. We decided to club the remaining
significant tags under their respective larger categories (eg: noun
phrases(NNP) were clubbed into noun category(NN)) and marked the remaining
small number of categories as <UNK> while also one-hot encoding all of them.
Lastly we standardized the targets for better training stability and it was
applied after transforming GPT to TRT-GPT.

Feature Name |  Feature Description |  Average Score and Std Deviation
---|---|---
word len |  Number of characters in the word |  0.6145
0.07363
is stopword |  Whether a word is a stopword |  0.6475 0.06043
is number |  Whether a word is a number |  0.6194 0.06731
lem word len |  Difference in the number of characters in the word and its lemmatized form |  0.6125 0.06649
is endword |  Whether the word is the last word of its sentence |  0.5918 0.10606
tfidf score |  TF-IDF score of the word in the corresponding sentence |  0.6265 0.06367
pos tag |  Part of Speech Tag of the word |  0.6284 0.059786
Table 6.1: Features and their Description and the Mean and Std Deviation of R2 scores of target variables Hyperparameter | Value
---|---
Batch Size | 4
Num Warm Up Steps | 4
Learning Rate | 3e-5
Max epochs | 120
Beta 1 | 0.91
Delta | 0.0001
Beta 2 | 0.998
Early stopping patience | 8
Weight Decay | 1e-5
TF-IDF Error | 0.1
Validation Split Ratio | 0.2
Training Split Ratio | 0.8
Table 6.2: Hyperparameter Values

  * •

Effect of scaling:  Due to a large number of features with a large difference
in their means, it is helpful to scale the features to the same range/same
mean and standard deviation. Hence we compared two scaling methods -

    * –

Min-max scaling:
$X_{norm}=(X-X_{min})/(X_{max}-X_{min})$

    * –

Standard scaling:
$X_{norm}=(X-X_{mean})/(X_{std})$

Model | Feature Norm | Extra Features |  nFix |  FFD |  GPT |  TRT |  fixProp
---|---|---|---|---|---|---|---
RoBERTa | Std Scl | STT |  0.8842 |  0.9246 |  0.7343 |  0.8823 |  0.9509
RoBERTa | Min Max | All |  0.8835 |  0.9132 |  0.7383 |  0.8542 |  0.9461
RoBERTa | Min Max | STT |  0.8417 |  0.9090 |  0.6456 |  0.8152 |  0.9492
BERT | Min Max | All |  0.7433 |  0.8528 |  0.4574 |  0.6925 |  0.9173
BiLSTM Minibatch | Min Max | All + GloVE |  0.6595 |  0.7275 |  0.6281 |  0.5988 |  0.7583
BiLSTM Minibatch | Min Max | All + NB111https://github.com/commonsense/conceptnet-numberbatch |  0.6379 |  0.7203 |  0.5992 |  0.5388 |  0.7592
RoBERTa Token | Std Scl | All |  0.7038 |  0.6512 |  0.4697 |  0.6877 |  0.7146
BERT Token | Std Scl | All |  0.7231 |  0.7107 |  0.3440 |  0.6376 |  0.7635
BERT | Std Scl | All |  0.6497 |  0.6456 |  0.5409 |  0.6088 |  0.7171
BiLSTM MiniBatch | Std Scl | All+GloVE |  0.6850 |  0.5976 |  0.5267 |  0.5940 |  0.7158
RoBERTa Token | Std Scl | All |  0.6483 |  0.5861 |  0.5475 |  0.6280 |  0.6566
RoBERTa | Std Scl | All |  0.6360 |  0.5869 |  0.4798 |  0.6061 |  0.6623
BiLSTM Stochastic | Std Scl | All+GloVE |  0.6282 |  0.6247 |  0.4394 |  0.5860 |  0.6893
RoBERTa Token | Std Scl | All |  0.6349 |  0.5914 |  0.4767 |  0.6123 |  0.6452
RoBERTa Token | Std Scl | All |  0.5991 |  0.5847 |  0.5417 |  0.5686 |  0.6524
BERT Token | Std Scl | All |  0.6182 |  0.5677 |  0.4495 |  0.6026 |  0.6960
BERT | Std Scl | All |  0.6019 |  0.5847 |  0.4159 |  0.5771 |  0.6404
BERT | Std Scl | All |  0.5681 |  0.5437 |  0.4891 |  0.5170 |  0.6371
DistilBert | No Scaling | All |  0.4804 |  0.5292 |  0.0983 |  0.4183 |  0.6324
Table 6.3: Compiled results of all the top performing models. STT means Single
Target Training (one model to predict one feature). Numberbatch refers to the
ConceptNet Numberbatch word embeddings used Speer et al. (2017). The
architecture was kept same (except for input dimensions), the learning rate
was reduced for smaller number of features to achieve optimal values which we
observed through trial and error.

With the same model (multi-target with RoBERTa encoder), using standard
scaling gave an R2 score of 0.62 while min-max scaling gave a score of 0.81.
The reason for min-max scaling giving significantly better performance might
be because the targets in this dataset are always positive.

  * •

Single Target vs. Multi-Target models:  With a large number of input features
and a relatively small datasets, we hypothesized that large pre-trained models
may fail to converge to a suitable minima in the limited fine-tuning. Hence,
we trained separate instances of the same model on each of nFix, FFD, TRT, GPT
and fixProp. This further increased the R2 score on the validation set to
0.88.

##  7 Results

In 6.3, we tabulate the R2 scores of each of our models on the validation set.

## References

  * Bird and Loper (2004) Steven Bird and Edward Loper. 2004.  NLTK: The natural language toolkit.  In _Proceedings of the ACL Interactive Poster and Demonstration Sessions_ , pages 214–217, Barcelona, Spain. Association for Computational Linguistics.
  * Clark et al. (2020) Kevin Clark, Minh-Thang Luong, Quoc V. Le, and Christopher D. Manning. 2020.  Electra: Pre-training text encoders as discriminators rather than generators.
  * Clifton et al. (2007) Charles Clifton, Adrian Staub, and Keith Rayner. 2007.  Chapter 15 - eye movements in reading words and sentences.  In Roger P.G. Van Gompel, Martin H. Fischer, Wayne S. Murray, and Robin L. Hill, editors, _Eye Movements_ , pages 341–371. Elsevier, Oxford.
  * Demberg and Keller (2008) Vera Demberg and Frank Keller. 2008.  Data from eye-tracking corpora as evidence for theories of syntactic processing complexity.  _Cognition_ , 109(2):193–210.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  Bert: Pre-training of deep bidirectional transformers for language understanding.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780.
  * Hollenstein et al. (2019) Nora Hollenstein, Maria Barrett, Marius Troendle, Francesco Bigiolli, Nicolas Langer, and Ce Zhang. 2019.  Advancing nlp with cognitive language processing signals.
  * Hollenstein et al. (2020) Nora Hollenstein, Marius Troendle, Ce Zhang, and Nicolas Langer. 2020.  Zuco 2.0: A dataset of physiological recordings during natural reading and annotation.
  * Just and Carpenter (1980) Marcel Adam Just and Patricia A. Carpenter. 1980.  A theory of reading: from eye fixations to comprehension.  _PSYCHOLOGICAL REVIEW_ , 87(4).
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized BERT pretraining approach.  _CoRR_ , abs/1907.11692.
  * Papoutsaki et al. (2016) Alexandra Papoutsaki, Patsorn Sangkloy, James Laskey, Nediyana Daskalova, Jeff Huang, and James Hays. 2016.  Webgazer: Scalable webcam eye tracking using user interactions.  In _Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, IJCAI 2016, New York, NY, USA, 9-15 July 2016_ , pages 3839–3845. IJCAI/AAAI Press.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. 2011.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher D. Manning. 2014.  Glove: Global vectors for word representation.  In _Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1532–1543.
  * Sanh et al. (2020) Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf. 2020.  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.
  * Speer et al. (2017) Robyn Speer, Joshua Chin, and Catherine Havasi. 2017.  ConceptNet 5.5: An open multilingual graph of general knowledge.  pages 4444–4451.
  * Speer et al. (2018) Robyn Speer, Joshua Chin, and Catherine Havasi. 2018.  Conceptnet 5.5: An open multilingual graph of general knowledge.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.