Address after May 1, 2022: ]Institute of Physics AS CR, ELI-Beamlines project,
Za Radnicí 835, 25241 Dolní Břežany, Czech Republic.

# Dynamic Magnetic Response of Quark-Gluon Plasma to Electromagnetic Fields

Christopher Grayson<EMAIL_ADDRESS>Martin Formanek
<EMAIL_ADDRESS>[  Johann Rafelski<EMAIL_ADDRESS>Department of Physics, The University of Arizona, Tucson, AZ, 85721, USA
Berndt Müller<EMAIL_ADDRESS>Department of Physics, Duke University,
Durham, NC 27708-0305, USA

(4/29/2022)

###### Abstract

We investigate the electromagnetic response of a viscous quark-gluon plasma in
the framework of the relativistic Boltzmann equation with current conserving
collision term. Our formalism incorporates dissipative effects at all orders
in linear response to the electromagnetic field while accounting for the full
space- and time-dependence of the perturbing fields. As an example, we
consider the collision of two nuclei in a stationary, homogeneous quark-gluon
plasma. We show that for large collision energies the induced magnetic fields
are governed by the response of quark-gluon plasma along the light-cone. In
this limit we derive an analytic expression for the magnetic field along the
beam axis between the receding nuclei and show that its strength varies only
weakly with collision energy for $\sqrt{s_{\mathrm{NN}}}\geq 30$ GeV.

##  I Introduction

The electromagnetic fields generated by colliding relativistic heavy ions at
RHIC and LHC are some of the largest in nature, on the order of
$m_{\pi}^{2}/e\approx 10^{14}$ T, but they are very short-lived with a time
constant $t_{\text{coll}}=2R/\gamma\sim 10^{-25}\,\textrm{s}$, where $R$ is
the nuclear radius and $\gamma$ is the Lorentz factor. These fields occur in
the presence of a quark-gluon plasma (QGP), the high temperature matter that
forms in the region of overlap between the colliding ions. Within the plasma
these strong electromagnetic fields generate a response that scales with
$Z\alpha$, where $Z$ is the charge of the ions and $\alpha$ is the fine-
structure constant. The polarization current is subject to the strong
dissipative effects in the QGP that render the QGP an over-damped system which
sustains slowly decaying magnetic wake fields.

In the past, most analytic calculations [1, 2, 4, 3, 5, 7, 6] assumed that the
response of the QGP is static and obtained the late time dependence of the
induced magnetic field by solving Maxwell’s equations in the presence of
static electric conductivity in a hydrodynamically evolving QGP. These
approaches ignore the fact that the time scale of variation of the
electromagnetic fields, $t_{\text{coll}}$, is shorter than the characteristic
time scales of the plasma response, given by the plasma frequency $\omega_{p}$
and the collisional damping rate $\kappa$, as shown in Table 1. Some numerical
calculations [8, 9] have attempted to incorporate the dynamical response of
QGP by numerically solving the coupled magnet-hydrodynamic equations for a
conducting quark-gluon plasma in the presence of the colliding nuclear
charges. More recent calculations [10, 11] also incorporate the frequency and
wave-vector dependence of QGP response to electromagnetic fields by solving
the coupled Vlasov-Boltzmann–Maxwell equations. The disadvantage of these
approaches, apart from the fact that they are prone to numerical imprecision,
is that they do not easily allow for parametric predictions of how the medium
response depends on the properties of the QGP and energy of the colliding
nuclei.

Here we address this issue by calculating the response of a QGP to heavy ion
fields using the full frequency and wave-vector dependent relativistic
polarization tensor previously obtained in [12] by using relativistic kinetic
theory to derive a gauge invariant polarization tensor within the linear
response formalism. As we are interested in high collision energies, we
consider the case of a baryon number-symmetric plasma. The result of [12]
accounts for both, collective plasma response and viscous effects, and allows
us to obtain analytical expressions for the induced fields and currents that
depend parametrically on the electromagnetic Debye mass $m_{D}$ and the
average collision rate $\kappa$. We expect our results to have validity as an
effective description of the plasma response that is independent of the
detailed assumptions for the QGP transport model.

The average collision rate $\kappa$ is inversely related to the mean-free time
of partons in the plasma, $\kappa=1/\tau_{\text{rel}}$, also called relaxation
time. We will show in Sect. III.3 that, for a particle-antiparticle plasma,
current conservation implies energy momentum conservation to first order in
the plasma response. The collision term generates damping in the evolution of
the one-particle distribution function $f(x,p)$ in response of the
electromagnetic currents and velocity gradients. This gives the plasma
dissipative properties parametrized by the shear viscosity and the electric
conductivity, which are related to the collision rate $\kappa$ by (see e. g.
[13]):

$\eta=\frac{sT}{5\kappa}\,,\qquad\sigma_{0}=\frac{e^{2}T^{2}}{3\kappa}\,,$ (1)

where $s$ denotes the entropy density and $T$ is the temperature. The
treatment of conservation laws in transport coefficients in the relaxation
time approximation is discussed in [14], where it is shown that in the case of
electrical conductivity and shear viscosity the BGK (Bhatnagar-Gross-Krook)
term [15] is the only required correction to the collision term.

Here we focus on the damping effect of the collision rate $\kappa$ on the
induced magnetic fields in nuclear collisions. The collision term generates a
nonvanishing conductivity via the imaginary part of the polarization tensor.
This conductivity manifests itself in poles of the resummed propagator in the
lower complex $\omega$-plane that generate long-range tails or wake fields
that extend far beyond the collision time. In Table 1 we collect the relevant
timescales of the problem in ascending order. The collision time
$t_{\text{coll}}$ is much shorter than all other relevant time scales. The
collective oscillations of the plasma are highly damped by the large value of
the relaxation time giving rise to over-damped behavior.

Table 1:  Approximate time scales relevant to the electromagnetic response of QGP for an Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\,$MeV with QGP temperature $T=300$ MeV. Time scales are shown in ascending order. Time Scale | Formula | Time (fm/c)
---|---|---
Collision Time | $t_{\text{coll}}=2R/\gamma$ | 0.086111Calculated using the Gaussian radius $R=4.33$ fm defined in Eq. (61).
Relaxation Time | $\tau_{\text{rel}}=1/\kappa$ | 0.36
Freeze-out time | $t_{f}$ |  $5$222Estimated using 2+1 dimensional hydrodynamic evolution [16].
Decay Time | $t_{\sigma}=1/\sigma_{0}=\kappa/\omega_{p}^{2}$ | 59333The decay time is the large damping $\kappa/\omega_{p}$ expansion of the plasma oscillation frequency Eq. (36).

The layout of this work is as follows. We begin by solving Maxwell’s equations
in a homogeneous, stationary, polarizable medium in the limit of linear
response in Sect. II. We go on to summarize the calculation of the
electromagnetic polarization tensor of QGP following [12] and briefly discuss
energy-momentum conservation in Sect. III, where we also present our choice of
the QGP parameters. We then explore the magnetic fields generated in a nuclear
collision in Sect. IV. We compare the magnetic field calculated with full
space-time dependence in the polarization tensor to previously used
approximations and find an analytic expression for the late time magnetic
field at the collision center. We summarize our results in Sect. V.

##  II Maxwell equations in QGP

We begin by solving Maxwell’s equations in an infinite homogeneous and
stationary polarizable medium. In this medium Maxwell’s equations take on the
usual form

$\partial^{[\mu}F^{\nu\rho]}(x)=0,\quad\partial_{\mu}F^{\mu\nu}(x)=\mu_{0}J^{\nu}(x)\,,$
(2)

or in Fourier space using the replacement $\partial^{\mu}\rightarrow-ik^{\mu}$

$-ik^{[\mu}\widetilde{F}^{\nu\rho]}(k)=0,\quad-
ik_{\mu}\widetilde{F}^{\mu\nu}(k)=\mu_{0}\widetilde{J}^{\nu}(k)\,.$ (3)

The properties of the medium are introduced by writing the 4-current
$\widetilde{J}^{\mu}$ in terms of its induced and external parts

$\widetilde{J}^{\mu}(k)=\widetilde{j}_{\mathrm{ext}}^{\mu}(k)+\widetilde{j}_{\mathrm{ind}}^{\mu}(k)\,.$
(4)

The induced current $\widetilde{j}_{\mathrm{ind}}^{\mu}$ to leading order is
given by the polarization tensor through the covariant formulation of Ohm’s
law [17]

$\widetilde{j}_{\mathrm{ind}}^{\mu}(k)={\Pi^{\mu}}_{\nu}(k)\widetilde{A}^{\nu}(k)\,,$
(5)

where $\Pi^{\mu}_{\nu}$ is the polarization tensor and $\widetilde{A}^{\mu}$
is the 4-vector potential. The electromagnetic tensor $\widetilde{F}^{\mu\nu}$
can be written in terms of the 4-vector potential

$\widetilde{F}^{\mu\nu}(k)=-i\left(k^{\mu}\widetilde{A}^{\nu}(k)-k^{\nu}\widetilde{A}^{\mu}(k)\right)\,.$
(6)

Solving Maxwell’s equations in the Lorentz gauge $k\cdot\widetilde{A}=0$ one
finds the usual expression

$\widetilde{A}^{\mu}(k)=-\frac{\mu_{0}}{k^{2}}\left(\widetilde{j}_{\mathrm{ext}}^{\mu}(k)+\widetilde{j}_{\mathrm{ind}}^{\mu}(k)\right)\,$
(7)

where $\mu_{0}$ denotes the magnetic permittivity of the vacuum.

We proceed by algebraically solving for the self-consistent potential that
contains the back-reaction of the induced current on the field to all orders
in the Debye screening mass $m_{D}$. To do this, we first note that in a
homogeneous medium the response depends only on two independent scalar
polarization functions $\Pi_{\parallel}$ and $\Pi_{\perp}$ describing
polarization in the parallel and transverse directions relative to the wave-
vector $\boldsymbol{k}$ [18]. The polarization tensor may be written in terms
of these polarization functions as

$\Pi^{\mu\nu}(k,u)=\Pi_{\parallel}(k)L^{\mu\nu}(k,u)+\Pi_{\perp}(k)S^{\mu\nu}(k,u)\,,$
(8)

where $k^{\mu}$ is the 4-momentum of the field and $u^{\mu}$ is the 4-velocity
of the medium. Conventions for the longitudinal and transverse projection
tensors, $L^{\mu\nu}$ and $S^{\mu\nu}$ respectively, may be found in [19].
These tensors are reproduced here for convenience

$L^{\mu\nu}\equiv\frac{k^{2}}{(k\cdot
u)^{2}-k^{2}}\bigg{[}\frac{k^{\mu}u^{\nu}}{(k\cdot
u)}+\frac{k^{\nu}u^{\mu}}{(k\cdot u)}\\\ -\frac{k^{2}u^{\mu}u^{\nu}}{(k\cdot
u)^{2}}-\frac{k^{\mu}k^{\nu}}{k^{2}}\bigg{]}\,,$ (9) $S^{\mu\nu}\equiv
g^{\mu\nu}+\frac{1}{(k\cdot u)^{2}-k^{2}}\bigg{[}k^{\mu}k^{\nu}\\\ -(k\cdot
u)(k^{\mu}u^{\nu}+k^{\nu}u^{\mu})+k^{2}u^{\mu}u^{\nu}\bigg{]}\,.$ (10)

These projections are equivalent to ones defined in [18] up to an overall
normalization. To simplify the calculation the wave-vector $\boldsymbol{k}$ is
chosen, without loss of generality, to point along the third spatial direction
($\mu=3$):

$\Pi^{\mu}_{\nu}(\omega,\boldsymbol{k})=\left[\begin{array}[]{cccc}-\frac{|\mathbf{k}|^{2}}{\omega^{2}}\Pi_{\parallel}&0&0&\frac{|\mathbf{k}|}{\omega}\Pi_{\parallel}\\\
0&\Pi_{\perp}&0&0\\\ 0&0&\Pi_{\perp}&0\\\
-\frac{|\mathbf{k}|}{\omega}\Pi_{\parallel}&0&0&\Pi_{\parallel}\\\
\end{array}\right]\,.$ (11)

Utilizing this decomposition the spatial component of the potential
$\widetilde{\boldsymbol{A}}$ is expressed as

$\widetilde{\boldsymbol{A}}=\widetilde{A}_{\parallel}\hat{\boldsymbol{k}}+\widetilde{\boldsymbol{A}}_{\perp}\,,$
(12)

which implies

$\widetilde{A}_{\parallel}=\frac{\boldsymbol{k}\cdot\widetilde{\boldsymbol{A}}}{|\boldsymbol{k}|},\quad\widetilde{\boldsymbol{A}}_{\perp}=\widetilde{\boldsymbol{A}}-\widetilde{A}_{\parallel}\hat{\boldsymbol{k}}\,,$
(13)

with analogous definitions for the current, $\widetilde{j}_{\parallel}$ and
$\widetilde{j}_{\perp}$. Note that the Lorentz gauge condition
$\partial_{\mu}A^{\mu}=0$ implies

$\widetilde{A}_{\parallel}=\frac{\omega}{|\boldsymbol{k}|}\widetilde{\phi}\,$
(14)

with $\phi=A^{0}$. The resulting induced charge can be calculated using the
projected polarization tensor Eq. (11):

$\widetilde{\rho}_{\text{ind}}(\omega,\boldsymbol{k})=\Pi^{0}_{\nu}\widetilde{A}^{\nu}=-\frac{|\boldsymbol{k}|^{2}}{\omega^{2}}\Pi_{\parallel}\widetilde{\phi}+\frac{|\boldsymbol{k}|}{\omega}\Pi_{\parallel}\widetilde{A}_{\parallel}\,.$
(15)

For the Lorentz gauge condition Eq. (14) one finds

$\widetilde{\rho}_{\text{ind}}(\omega,\boldsymbol{k})=\Pi_{\parallel}\widetilde{\phi}\left(1-\frac{|\boldsymbol{k}|^{2}}{\omega^{2}}\right)\,.$
(16)

Similarly,

$\widetilde{j}_{\parallel\text{ind}}(\omega,\boldsymbol{k})=\Pi^{z}_{\nu}\widetilde{A}^{\nu}=\Pi_{\parallel}\frac{\omega}{|\boldsymbol{k}|}\widetilde{\phi}\left(1-\frac{|\boldsymbol{k}|^{2}}{\omega^{2}}\right)\,,$
(17)

expressing current conservation, $\partial_{\mu}j^{\mu}=0$. The induced
transverse current is

$\boldsymbol{j}_{\perp\text{ind}}(\omega,\boldsymbol{k})=\Pi_{\perp}\widetilde{A}_{\perp}\,.$
(18)

Solving for the potential on both sides of Eq. (7) with the help of Eqs.
(16-18) gives the self-consistent solutions,

$\displaystyle\widetilde{\phi}(\omega,\boldsymbol{k})=\frac{\widetilde{\rho}_{\text{ext}}(\omega,\boldsymbol{k})}{\varepsilon_{0}(\boldsymbol{k}^{2}-\omega^{2})\left(\Pi_{\parallel}/(\omega^{2}\varepsilon_{0})+1\right)}\,,$
(19)
$\displaystyle\widetilde{\boldsymbol{A}}_{\perp}(\omega,\boldsymbol{k})=\frac{\mu_{0}\widetilde{\boldsymbol{j}}_{\perp\text{ext}}(\omega,\boldsymbol{k})}{\boldsymbol{k}^{2}-\omega^{2}-\mu_{0}\Pi_{\perp}}\,.$
(20)

The gauge condition Eq. (14) gives the self-consistent potential
$\widetilde{A}_{\parallel}$. These self-consistent potentials determine the
electric and magnetic fields via the usual relations

$\widetilde{\boldsymbol{B}}(\omega,\boldsymbol{k})=i\boldsymbol{k}\times\widetilde{\boldsymbol{A}}_{\perp}\,,\quad\widetilde{\boldsymbol{E}}(\omega,\boldsymbol{k})=-i\boldsymbol{k}\widetilde{\phi}+i\omega\widetilde{\boldsymbol{A}}\,.$
(21)

The magnetic field of interest to us is obtained by Fourier transforming these
expressions back to position space (for details see Appendix B). We note here
that it is important to calculate the fields with the resummed propagator in
Eqs. (19-20), rather than by using its series representation, in order to
correctly capture the pole structure of the propagator that governs the late-
time dependence of the magnetic field. For ease of calculation, we specify the
external 4-current using two colliding Gaussians charge distributions
normalized to the nuclear rms radius $R$ and charge $Z$:

$\rho_{\text{ext}\pm}(t,\boldsymbol{x})=\frac{Zq\gamma}{\pi^{3/2}R^{3}}e^{-\frac{1}{R^{2}}(x\mp
b/2)^{2}}e^{-\frac{1}{R^{2}}y^{2}}\\\ \times
e^{-\frac{\gamma^{2}}{R^{2}}(z\mp\beta t)^{2}}\,,$ (22)

where $\gamma$ and $\beta$ are the Lorentz factor and speed, respectively, and
$b$ is the impact parameter of the collision. The plus and minus signs
indicate motion in the $\pm\hat{z}$-direction (beam-axis). This charge
distribution corresponds to the vector current

$\boldsymbol{j}_{\text{ext}\pm}(t,\boldsymbol{x})=\pm\beta\hat{\boldsymbol{z}}\rho_{\text{ext}\pm}(t,\boldsymbol{x})\,.$
(23)

Further details of the external charge distribution for two colliding nuclei
are presented in Appendix A..

##  III The QGP Polarization tensor

###  III.1 Derivation of $\Pi^{\mu}_{\nu}$

In this Section we derive the polarization tensor, including damping, for the
idealized case where the QGP is homogeneous and stationary. We follow the
derivation presented in [12] for the damped polarization tensor of an
electron-positron plasma. The calculation differs slightly from [12], since in
QGP we consider three quark species: up, down, and strange. We start from the
Vlasov-Boltzmann equation for each quark flavor:

$(p\cdot\partial)f_{f}(x,p)+q_{f}F^{\mu\nu}p_{\nu}\frac{\partial
f_{f}(x,p)}{\partial p^{\mu}}=(p\cdot u)C_{f}(x,p)\,,$ (24)

The collision term $C_{f}(x,p)$ in the BGK form is given by

$C_{f}(x,p)=\kappa_{f}\left(f^{(\mathrm{eq})}_{f}(p)\frac{n_{f}(x)}{{n^{(\mathrm{eq})}_{f}}}-f_{f}(x,p)\right)\,,$
(25)

where plasma constituents collide on a momentum-averaged time scale
$\tau_{\text{rel}}=\kappa^{-1}$. The collision term is constructed such that
Eq. (24) retains current conservation [15]. We show in Sect. III.3 that energy
is also conserved for the case of a neutral particle-antiparticle plasma at
linear order in the external field.

The induced current $j_{\mathrm{ind}}^{\mu}$ can be written in terms of the
phase-space distribution of quarks and anti-quarks as

$j_{\mathrm{ind}}^{\mu}(x)=2N_{c}\int(dp)p^{\mu}\\\
\times\sum_{u,d,s}q_{f}(f_{f}(x,p)-f_{\bar{f}}(x,p))\,,$ (26)

where $N_{c}$ is the number of colors, and we sum over the quark flavors with
charges $q_{f}$. One can calculate the induced current for small perturbations
away from equilibrium for each quark flavor

$f_{f}(x,p)={f^{(\mathrm{eq})}_{f}}(p)+\delta f_{f}(x,p)\,,$ (27)

Note that the equilibrium contributions ${f^{(\mathrm{eq})}_{f}}(p)$ do not
contribute to Eq. (26) because of the opposite sign of the charges of
particles and antiparticles, but the perturbations $\delta f$ add up due the
change in sign of the external force $qF^{\mu\nu}p_{\nu}$:

$j_{\mathrm{ind}}^{\mu}(x)=2N_{c}\int(dp)p^{\mu}\sum_{u,d,s}q_{f}(\delta
f_{f}(x,p)-\delta f_{\bar{f}}(x,p))\\\
=4N_{c}\int(dp)p^{\mu}\sum_{u,d,s}q_{f}^{2}\delta f(x,p)\\\
=4N_{Q}e^{2}\int(dp)p^{\mu}\delta f(x,p)\,.$ (28)

In the second line we pulled out a factor of electric charge $\delta
f_{f}=q_{f}\delta f$. The perturbations $\delta f$ are identical for all quark
species in the ultrarelativistic limit. The result Eq. (28) differs from that
found in the case of an electron-positron plasma by the factor

$N_{Q}\equiv N_{c}\sum_{f}(q_{f}/e)^{2}=2\,,$ (29)

where the numerical value holds for three light quarks flavors ($u,d,s$). We
refer to [12] for the derivation of the polarization tensor in terms of
integrals over the phase-space distribution $\delta f$, because the only
difference is the overall factor $N_{Q}$.

As noted in the previous Section, the polarization tensor in Eq. (8) may be
written in terms of two independent components: the longitudinal polarization
function $\Pi_{\parallel}$, which describes response parallel to wave-vector
$\boldsymbol{k}$, and the transverse polarization function $\Pi_{\perp}$,
which describes response in the plane perpendicular to wave-vector
$\boldsymbol{k}$. When the $\mu=3$ ($z$) axis is chosen along the wave-vector
$\boldsymbol{k}$, the longitudinal and transverse polarization functions
relate to the components of the polarization tensor Eq. (11) along the
coordinate axes as

$\Pi_{\parallel}=\Pi^{3}_{3},\quad\Pi_{\perp}=\Pi^{1}_{1}=\Pi^{2}_{2}\,.$ (30)

In the ultrarelativistic limit, neglecting quark masses, one finds [12]:

$\displaystyle\Pi_{\parallel}(\omega,|\boldsymbol{k}|)=m_{D}^{2}\frac{\omega^{2}}{\boldsymbol{k}^{2}}\left(1-\frac{\omega\Lambda}{2|\boldsymbol{k}|-i\kappa\Lambda}\right)\,,$
(31)
$\displaystyle\Pi_{\perp}(\omega,|\boldsymbol{k}|)=\frac{m_{D}^{2}\,\omega}{4|\boldsymbol{k}|}\left(\Lambda\left(\frac{\omega^{\prime
2}}{\boldsymbol{k}^{2}}-1\right)-\frac{2\omega^{\prime}}{|\boldsymbol{k}|}\right)\,,$
(32)

where $\Lambda(\omega,\boldsymbol{k})$ is defined as

$\displaystyle\Lambda\equiv\ln\frac{\omega^{\prime}+|\boldsymbol{k}|}{\omega^{\prime}-|\boldsymbol{k}|}\,,\quad\text{with}\quad\omega^{\prime}=\omega+i\kappa.$
(33)

The natural logarithm leads to branch cut in the complex $\omega$ plane
running from $-|\boldsymbol{k}|-i\kappa$ to $|\boldsymbol{k}|-i\kappa$ as
noted in [20]. The parallel and transverse polarization functions have the
same form as in [12] except for an overall factor $N_{Q}$ that is contained in
the leading order electromagnetic Debye mass for the QGP plasma [21]:

${m_{D}}^{2}_{(\text{EM})}=\sum_{u,d,s}q^{2}_{f}T^{2}\frac{N_{c}}{3}=N_{Q}\frac{e^{2}T^{2}}{3}\equiv
C_{\text{em}}T^{2}\,,$ (34)

where $C_{\text{em}}=2e^{2}/3$. In the following we will use $m_{D}$ as short-
hand notation for the electromagnetic screening mass since we do not discuss
color screening here.

The polarization tensor may be written in any general frame by using Eq. (8),
but for our purposes it will be simpler to carry out calculations in the
coordinate system where $\boldsymbol{k}$ aligns with the $z$-axis so that the
polarization tensor takes the form shown in Eq. (11).

###  III.2 QGP parameters

The strength of the plasma response to an external magnetic field depends on
the values of two physical parameters: the quark damping rate $\kappa$, and
the electromagnetic screening mass $m_{D}$. In this Section we provide
estimates for these parameters.

We adopt the perturbative result Eq. (34) to estimate $m_{D}$. Higher-order
corrections to this expression can been derived from higher-order calculation
of the vector spectral function in thermal perturbation theory (see [22] and
references cited therein).

The scale of the collisional quark damping $\kappa$ is much larger than the
electromagnetic Debye mass $m_{D}$ because it depends on the strong coupling
constant $\alpha_{s}$, not the electromagnetic coupling $\alpha$. Solving the
dispersion relation

$\frac{1}{(k\cdot
u)^{2}}(k^{2}+\mu_{0}\Pi_{\parallel}(\omega,k))(k^{2}+\mu_{0}\Pi_{\perp}(\omega,k))^{2}=0\,,$
(35)

see [19], in the limit $\boldsymbol{k}\rightarrow 0$ one finds for the plasma
oscillation frequency [12]

$\omega_{p}^{\pm}=-\frac{i\kappa}{2}\pm\sqrt{\frac{m_{D}^{2}}{3}-\frac{\kappa}{4}^{2}}\,.$
(36)

We see that if $\kappa>\tfrac{2}{\sqrt{3}}m_{D}$, the plasma oscillations are
over-damped.

Figure 1: Plot of the QED Debye mass and the QCD dampening rate $\kappa$ as a
function of temperature. At temperature $T=300\,$MeV used in the plots below,
$\kappa=4.86\,m_{D}$.

The collision rate $\kappa$ is related to the inverse of the mean-free time of
quarks in QGP. In kinetic theory the mean-free time is given by the product of
the parton density in the QGP and the quark-parton transport cross section,
leading to the expression [23]

$\kappa(T)=\frac{10}{17\pi}(9N_{f}+16)\zeta(3)\alpha_{s}^{2}\ln\left(\frac{1}{\alpha_{s}}\right)T\,,$
(37)

where $N_{f}$ is the number of flavors, $\zeta(x)$ denotes the Riemann zeta
function, and $\alpha_{s}(T)$ is the running QCD coupling. We model the
running of the QCD coupling constant as a function of temperature in the range
$T<5T_{c}$ using a fit provided in [24]:

$\alpha_{s}(T)\approx\frac{\alpha_{s}(T_{c})}{1+C\ln(T/T_{c})}\,,$ (38)

where $C=0.760\pm 0.002$. For the QCD (pseudo-)critical temperature we use
$T_{c}=160\,$MeV. The QED Debye mass is compared to $\kappa(T)$ in Fig. 1.

From $\kappa(T)$ in Eq. (37) and the running of the coupling in Eq. (38), we
calculate the static conductivity using the leading order electromagnetic
Debye mass $m_{D}$. The momentum dependent transverse conductivity
$\sigma_{\perp}$, which controls the response of the plasma to magnetic
fields, is related to the imaginary part of the transverse polarization
function $\Pi_{\perp}$ as follows [19]:

$\sigma_{\perp}(\omega,\boldsymbol{k})=-i\frac{\Pi_{\perp}(\omega,\boldsymbol{k})}{\omega}\,.$
(39)

In the long wavelength limit $\boldsymbol{k}\rightarrow 0$, the branch cut in
Eq. (33) shrinks to a single pole at $\omega=-i\kappa$, and the conductivity
has the simple form

$\sigma_{\perp}(\omega,0)=\sigma_{\parallel}(\omega,0)=\frac{\sigma_{0}}{1-i\omega/\kappa}\,.$
(40)

We will refer to $\sigma_{\perp}(\omega,0)$ as the Drude model [25]. In the
static limit $\omega\rightarrow 0$ the parallel and perpendicular
conductivities are the same, and the static conductivity $\sigma_{0}$ is given
by

$\sigma_{0}=\frac{m_{D}^{2}}{3\kappa}\,.$ (41)

The static conductivity determines the late time behavior of the magnetic
field.

###  III.3 Energy-momentum conservation

In general, the modified BGK collision term Eq. (25) violates energy and
momentum conservation. Rocha et al. [14] recently showed how energy-momentum
conservation can be restored by introducing a linearized collision operator
that is projected on eigenfunctions of the conserved quantities with
eigenvalue zero. Here we show explicitly that for a symmetric particle-
antiparticle plasma the energy momentum violations cancel at linear order in
the external field.

Recall that the energy momentum tensor $T^{\mu\nu}$ of the plasma is given by

$T^{\mu\nu}=2\int(dp)p^{\mu}p^{\nu}(f_{-}(x,p)+f_{+}(x,p))\,,$ (42)

where the factor of two accounts for spin and $f_{\pm}(x,p)$ represent the
distributions of particles ($+$) and antiparticles ($-$), respectively. We
recall that for $T^{\mu\nu}$ to be conserved the covariant divergence

$\partial_{\mu}T^{\mu\nu}=2\,\partial_{\mu}\int(dp)p^{\mu}p^{\nu}\left(f_{-}(x,p)+f_{+}(x,p)\right)$
(43)

must vanish. In linear response the distribution functions $f_{\pm}(x,p)$ are
given by

$f_{\pm}(x,p)=\delta f_{\pm}(x,p)+f^{(\mathrm{eq})}(p)\,.$ (44)

Equation Eq. (28) indicates that the perturbation $\delta f_{\pm}$ is linear
in the quark charge

$\delta f_{\pm}=\pm q\,\delta f\,.$ (45)

This leads to a cancellation of the particle and antiparticle perturbations in
the energy-momentum tensor at linear order:

$\partial_{\mu}T^{\mu\nu}=4\,\partial_{\mu}\left(\int(dp)p^{\mu}p^{\nu}f^{(\mathrm{eq})}(p)\right)=0\,.$
(46)

Thus for a symmetric particle-antiparticle plasma corrections to the energy-
momentum tensor appear only at second order in external field. This is a
general consequence of CPT symmetry of the medium.

##  IV Magnetic field in a nuclear collision

In this Section we calculate the magnetic field at the center of the heavy ion
collision by Fourier transforming the momentum space magnetic field Eq. (21)
to position space. We calculate the self-consistent magnetic field using the
potentials given in Eqs. (19-20) and model the response of QGP using the
idealized case of a homogeneous, stationary plasma detailed in Sects. II and
III. The external fields are specified by the moving Gaussian charge
distributions defined in Eqs. (22-23).

The magnetic field is of particular interest due to its role in the separation
of electric charge in the QGP through the chiral magnetic effect (CME) [26].
In the large magnetic fields that occur in heavy ion collisions the electric
current generated by the CME could lead to a charge separation along the
direction of the magnetic field. Whether this effect is observable depends
strongly on the size of the magnetic field. If a magnetic field of meaningful
strength survives until the time of hadronization of the QGP, it could also
lead to a difference in the global polarization of $\Lambda$ hyperons and
antihyperons [28].

We chose the collision center as the origin of our spatial coordinate system
and align the spatial $z$-axis with the beam direction. We calculate the
magnetic field along the $z$-axis by Fourier transforming the momentum space
expressions given in Eqs. (20-21):

$\boldsymbol{B}(t,z)=\int\frac{d^{4}k}{(2\pi)^{4}}e^{-i\omega
t+ik_{z}z}\frac{\mu_{0}i\boldsymbol{k}\times\widetilde{\boldsymbol{j}}_{\perp\text{ext}}(\omega,\boldsymbol{k})}{\boldsymbol{k}^{2}-\omega^{2}-\mu_{0}\Pi_{\perp}(\omega,\boldsymbol{k})}$
(47)

to position space. It is convenient to perform the Fourier integrals in
cylindrical coordinates $(\boldsymbol{x}_{\perp},z)$. The angular integral
$d\theta$ and the integral over momentum along the beam axis $dk_{z}$ can be
performed exactly. The $dk_{z}$ integral is trivial due to the delta function
in the external charge distribution Eq. (66). The frequency integral $d\omega$
and the transverse momentum integral $dk_{\rho}$ must, in general, be done
numerically. We present the details of this calculation in Appendix B. Due to
the symmetry of the colliding ions, the only nonzero component of the magnetic
field along the $z$-axis points out of the collision plane ($x-y$ plane). In
our coordinate system, described in Appendix A, this corresponds to the
$y$-component of the magnetic field. The numerical results for the position-
space magnetic field are shown in Fig. 3 and compared with earlier results.

To connect to these previous studies, we compute the magnetic field in
position space at the origin in various levels of approximation defined in
Eqs. (39-41) and Eq. (51).

Table 2:  Conductivity models used to calculate the resulting magnetic field. Each conductivity represents the response of QGP with a different spacetime dependence.  Conductivity | Dependence | Formula
---|---|---
Full | $\sigma_{\perp}(\omega,\boldsymbol{k})$ | $-i\Pi_{\perp}(\omega,\boldsymbol{k})/\omega$
Light-cone | $\sigma_{\perp}(\omega=|\boldsymbol{k}|)$ | Eq. (51)
Drude | $\sigma_{\perp}(\omega,0)$ | $\sigma_{0}/(1-i\omega/\kappa)$
Static | $\sigma_{0}$ | $m_{D}^{2}/(3\kappa)$

These conductivities, collected in Table 2, refer to different treatments of
the frequency $\omega$ and wave-vector $\boldsymbol{k}$ dependence of the
conductivity $\sigma_{\perp}(\omega,\boldsymbol{k})$. For instance, solving
for the magnetic field in the limit $\boldsymbol{k}\rightarrow 0$ assumes that
the spatial dependence of the external field can be neglected, not
superficially a good approximation because at any given time $t$ the field
varies rapidly with $z$. The levels of approximation we consider include: the
full space- and time-dependence of the conductivity
$\sigma_{\perp}(\omega,\boldsymbol{k})$, the Drude model Eq. (40), and the
static response $\sigma_{\perp}(0,0)$. We list these limits in Eqs. (39-41),
respectively.

The fourth limit we are considering is the conductivity along the light-cone
$\sigma_{\perp}(\omega,k_{z}=\pm\omega,k_{\rho}=0)$. We now show that the
light-cone limit closely resembles the Drude model. We first recall that the
frequency dependence of the free charge distribution in cylindrical
coordinates Eq. (66) has the form

$\widetilde{\rho}_{\text{ext}\pm}(\omega,\boldsymbol{k})=2\pi
Zq\,e^{-(k_{\rho}^{2}+k_{z}^{2}/\gamma^{2})\frac{R^{2}}{4}}\\\ \times
e^{\mp\frac{ik_{\rho}b\cos\theta}{2}}\delta(\omega\mp k_{z}\beta)\,.$ (48)

After performing the Fourier transform over the parallel component of the
wave-vector $k_{z}$ using the delta function, the magnitude of the wave-vector
$|\boldsymbol{k}|$ is effectively set to the light-cone
$\omega\approx|\boldsymbol{k}|$, with a small deviation due to the transverse
dependence of the field,

$|\boldsymbol{k}|^{2}=k_{z}^{2}+k_{\rho}^{2}\rightarrow(\omega/\beta)^{2}+k_{\rho}^{2}\,.$
(49)

Inspecting the external charge distribution after this replacement

$\widetilde{\rho}_{\text{ext}\pm}(\omega,k_{\rho})=2\pi
Zq\,e^{-(k_{\rho}^{2}+\omega^{2}/(\beta\gamma)^{2})\frac{R^{2}}{4}}\\\ \times
e^{\mp\frac{ik_{\rho}b\cos\theta}{2}}\,,$ (50)

we can see that the size of the deviation from the light-cone due to
$k_{\rho}$ is or order $O(1/R)$, while the width of the current distribution
in frequency space is of order $O(\beta\gamma/R)$.

Figure 2: The magnitude of the polarization tensor is plotted in momentum
space showing deviations in $k_{\rho}$ from the light-cone
$\omega=|\boldsymbol{k}|$ on the horizontal axis. The contours show lines of
constant magnitude of $\Pi_{\perp}(\omega,|\boldsymbol{k}|)$; lighter shading
indicates increasing magnitude. The dashed line encapsulates the $2\sigma$
support of the external charge distribution. The width of the external charge
distribution in momentum space is $\sqrt{2}/R$ in the transverse direction and
$\beta\gamma\sqrt{2}/R$ along the light-cone. One can see that in the region
sampled by the external charge distribution the polarization tensor is
effectively constant as a function of $k_{\rho}$.

The region of two-sigma support of the Gaussian charge distribution is shown
as the region enclosed by the dashed line in Fig. 2. The polarization tensor
is approximately constant as a function of $k_{\rho}$ in this region. This
implies that one can approximate the integral in Eq. (47) using the
polarization function at $k_{\rho}=0$, i. e. on the light-cone. This means
that the fields of the ions, traveling near the speed of light, probe the
polarization tensor along the light-cone. In this limit, the transverse
conductivity near the light-cone is

$\sigma_{\perp}(\omega=|\boldsymbol{k}|)=i\frac{m_{D}^{2}}{4\omega}\left(\frac{\kappa^{2}}{\omega^{2}}\xi\ln\xi+\frac{i\kappa}{\omega}\left(\xi+1\right)\right)\,,$
(51)

where $\xi$ is defined as

$\xi\equiv 1-2i\frac{\omega}{\kappa}\,.$ (52)

Since the light-cone conductivity only depends on a single variable
($\omega=|\boldsymbol{k}|$) it simplifies integrals involved in the Fourier
transform of fields back into position space.

Our results for the magnetic field at the collision center $B_{y}(t,0)$ are
shown in Fig. 3. The right panel of the figure shows the field at early times
($|t|<0.25~{}\text{fm}/c$) on a linear scale, the left panel shows the field
over a wider time range on a logarithmic scale. The most general case
$\sigma_{\perp}(\omega,\boldsymbol{k})$, shown as the dashed red curve in Fig.
3, includes the full time- and space-dependent response of the medium to the
fields of the colliding ions. The blue dashed curve shows the magnetic field
in the Drude model approximation Eq. (40), where the response depends only on
time. The magnetic field using the light-cone conductivity is seen as the gray
line overlapping the red dashed line in Fig. 3, where $\sigma_{0}$ is defined
in Eq. (41). The result of Fourier transforming this expression is shown as
the brown dotted curve in Fig. 3. Our results differ slightly from those of
[3] because here we account for the finite size of the ions and use a slightly
different conductivity value. Phantom text

Figure 3: The magnetic field at the collision center as a function of time,
with $T=300$ MeV for a Au-Au collisions ($Z=79$) at $\sqrt{s_{\text{NN}}}=200$
GeV and impact parameter $b=6.4\,$fm. The left panel shows the magnetic field
on a semi-logarithmic scale up to $ct=5$ fm. The right panel shows the early-
time magnetic field on a linear scale. At the chosen temperature the
electromagnetic Debye mass is $m_{D}=74\,$MeV and the quark damping rate is
$\kappa=4.86\,m_{D}$. This gives a static conductivity of
$\sigma_{0}=5.01\,$MeV. Comparing the different approximations we see that all
of them have similar asymptotic behavior. Only the Drude conductivity, the
light-cone limit of the conductivity, and the full conductivity
$\sigma_{\perp}(\omega,\boldsymbol{k})$ describe the field at early times.
Note that here that the plasma is considered homogeneous and stationary. In a
more realistic situation the field would become screened only after the QGP is
formed in the collision.

The magnetic field in the presence of a QGP was previously calculated using a
static conductivity in [3]. In this case, the magnetic field in Fourier space
has the form

$\widetilde{\boldsymbol{B}}(\omega,\boldsymbol{k})=\frac{\mu_{0}i\boldsymbol{k}\times\widetilde{\boldsymbol{j}}_{\perp\text{ext}}}{\boldsymbol{k}^{2}-\omega^{2}-i\omega\sigma_{0}}\,,$
(53)

Looking at the left panel of Fig. 3, one can see that every model of the
response function predicts similar values for the magnetic field approaching
the freeze-out time $t_{f}$. This is because the static conductivity
determines the late-time dependence of the magnetic field. As we discuss in
Appendix B, we can expect the static solution to match the full solution when
$t>1/\kappa$. The static conductivity initially overestimates the magnetic
field after the external field begins to fall, since the effect of dynamic
screening is not captured. This matches the qualitative picture given by the
detailed numerical transport calculation done in [11]. The full space-time
dependent model and the Drude model model behave similarly for most times, and
are almost identical for $t>1/\kappa\approx 0.36$ fm/c. The approximation of
the polarization tensor on the light-cone agrees very closely tracks the full
solution at all times.

We can use the light-cone conductivity in Eq. (51) to understand why the Drude
model $\sigma_{\perp}(\omega,0)$ matches the full solution for times
$t>1/\kappa$. Late times probe the small frequency limit of the conductivity.
An expansion of Eq. (51) in $\omega/\kappa$ yields

$\sigma_{\perp}(\omega=|\boldsymbol{k}|)=\sigma_{0}\left(1+i\omega/\kappa\right)\\\
-\frac{6\sigma_{0}}{5}\frac{\omega^{2}}{\kappa^{2}}+O\left(\frac{\omega^{3}}{\kappa^{3}}\right)\,.$
(54)

We then compare to the same expansion for the Drude conductivity

$\sigma_{\perp}(\omega,0)=\frac{\sigma_{0}}{1-i\omega/\kappa}\approx\sigma_{0}\left(1+i\omega/\kappa\right)\\\
-\sigma_{0}\frac{\omega^{2}}{\kappa^{2}}+O\left(\frac{\omega^{3}}{\kappa^{3}}\right)\,.$
(55)

The lowest-order term, which coincides with the expression for the Drude
model, closely approximates the full solution when $\kappa\gg\omega$ as shown
in Fig 4. Since $\kappa t_{f}\gg 1$ for the QGP, the series converges rapidly
for times of the order of the freeze-out time $t_{f}$.

Figure 4: Comparison of the conductivity on the light-cone to
$\sigma_{\perp}(\omega,\boldsymbol{k}\rightarrow 0)$, scaled with the static
conductivity. We see that at small $\omega/\kappa$, i.e. times much larger
than $1/\kappa$, both approximations converge to the static case, while they
diverge $\omega/\kappa>1$. This predicts that the Drude model will
underestimate screening at small times, which is exactly what we observe in
Fig. 3.

The simple form of the Drude approximation Eq. (40) allows one to find the
poles of the denominator in Eq. (47), analytically. The frequency integral can
then be done using the residue theorem, allowing for an approximate analytical
expression for the late-time magnetic field. This is done in Appendix B. In
the ultrarelativistic limit $\gamma\gg 1$ and large times $t\gg 1/\kappa$
gives

$B_{y}(t)=-\mu_{0}\frac{Zq\beta}{(2\pi)}\frac{b\sigma_{0}}{4t^{2}}e^{\frac{-b^{2}\sigma_{0}}{16t}}\,,$
(56)

This result differs from the “diffusive” solution of Tuchin [3] by a factor
1/4 in the exponent, due to their convention for impact parameter
$b\rightarrow 2b$. The reason why $\kappa$ does not appear in the expression
Eq. (56) for the late-time magnetic field lies in the hierarchy of time scales
$t_{\text{coll}}\ll 1/\kappa\ll t_{f}$, which makes plasma damping irrelevant
during the spike of the external field as well as at freeze-out.

Interestingly, this solution has a finite limit when $\gamma\rightarrow\infty$
as it depends only on $\beta$, but not on $\gamma$. This property, which was
first observed by Tuchin [3], can be understood as follows: For late times the
Fourier integral of Eq. (66) is dominated by contributions from small
frequencies $\omega$, and it is sufficient to consider the $\omega\rightarrow
0$ limit of the Fourier spectrum of the external charge distributions
$\widetilde{\rho}_{f\pm}$ given in Eq. (72). In this limit Eq. (66) takes the
form

$\widetilde{\rho}_{\text{ext}\pm}(0,\boldsymbol{k})\rightarrow 2\pi
Ze\,e^{-k_{\rho}^{2}R^{2}/4}e^{\mp\frac{ik_{\rho}b\cos\theta}{2}}\delta(k_{z}\beta)\,,$
(57)

which is independent of $\gamma$. This occurs because

$\widetilde{\rho}_{\text{ext}\pm}(0,\boldsymbol{k})=\int dt\int
d^{3}xe^{-\boldsymbol{k}\cdot\boldsymbol{x}}\rho_{\text{ext}\pm}(0,\boldsymbol{x})$
(58)

integrates over the passage of the entire nucleus at a given location
$\boldsymbol{x}$ and thus is independent of $\gamma$ as the total charge is
Lorentz invariant. We conclude that, quite generally, for high collision
energies the remnant magnetic field at late times is determined by the time-
integrated action of the external electromagnetic pulse on the QGP. In a more
realistic calculation, where the QGP is not present for the entire duration of
the pulse, because it is created during the collision, the remnant magnetic
field will be diminished as only a fraction of the pulse acts on the QGP. We
therefore expect our result to represent an upper bound to the late-time
magnetic field in a realistic collision scenario.

Figure 5: Plot of the freeze-out magnetic field for $T=150$ MeV. We expect
that around this temperature QGP will hadronize into a mixed phase [29]. The
approximate late time solution Eq. (56) shown as a black line is compared to
numerical calculations using the full polarization tensor Eq. (47). the
approximate solution does not fully match the ultrarelativistic limit until
times $t>t_{\sigma}\approx 59$ fm/c. The magnetic field is independent of the
beam energy over a wide range of $\gamma$ but begins to diverge slowly from
the ultrarelativistic case at around $\gamma\leq 15$ for the time window shown
in the figure. Lower beam energies result in a somewhat larger field at late
time.

The approximation used to derive this solution holds for
$\gamma\beta\gg\sqrt{\kappa/\sigma_{0}}\approx 12$. In Fig. 5 we compare Eq.
(83) to the full numerical result to explore its dependence on $\gamma$. One
can see that the ultrarelativistic case (black solid line) begins to diverge
from the numerical solution at around $\gamma\approx 15$ for the times shown.
The early time magnetic field is not shown because the initial temperature of
QGP will depend strongly on the collision energy. The times are chosen such
that they cover the range of freeze-out times predicted for QGP for the range
of experimental collision energies used [27]. We do not show curves for
$\gamma<10$ because we expect the effects of chemical potential will become
important, yet here chemical potential $\mu$ is set to zero. In Fig. 5 one can
see that the late-time magnetic field has a very soft dependence on collision
energy. The time at which the magnetic field freezes out, which varies with
collision energy, has a much stronger effect on the magnitude of the freeze-
out field.

As the QGP begins to hadronize at time $t_{f}$, one may expect hadrons to be
statistically polarized with respect to the magnetic field. In [28] the
measured difference in global polarization of hyperons and antihyperons is
used to give an upper bound on the magnetic field at QGP freeze-out, $B\sim
2.7\times 10^{-3}\,m_{\pi}^{2}$ for Au+Au collisions at
$\sqrt{s_{\text{NN}}}=200$ GeV. Looking at Fig. 5 the magnetic field for
$\gamma=100$ at QGP freeze-out $t_{f}\approx 5$ fm/c is predicted to be
$B\approx 1.2\times 10^{-3}\,m_{\pi}^{2}$, somewhat below this upper bound.
Note that this assumes the polarization rapidly equilibriates in the plasma.
It also neglects any interactions during the hadron gas phase of the
collision.

Figure 6: Comparison of the magnetic field for different values of quark
damping rate or, equivalently, electric conductivity. Larger values of the
damping rate $\kappa$ represent smaller conductivities and vice versa as
indicated by Eq. (41). The black dashed line and the solid black line
represent the limits of zero and infinite conductivity, respectively. One can
see that as $\kappa$ increases the asymptotic value of the magnetic field
decreases.

In Fig. 6 we look at the magnetic field at the origin for different values of
$\kappa$. Increasing $\kappa$ reduces the static conductivity $\sigma_{0}$
which decreases the asymptotic value of the magnetic field as indicated by Eq.
(56). As $\kappa$ goes to zero the results converge to the case of ideal
conductivity $\sigma_{0}\rightarrow\infty$ where the magnetic field quickly
approaches a constant value. This case was studied in [2] where the authors
considered a magnetic field that falls to a constant value and then decreases
with $1/t$ due to Bjorken flow. More recent calculations [10, 11] solve the
Vlasov-Boltzmann equation numerically with parton-parton scattering. The
magnetic field predicted by [10] is around $\sim 10^{-4}m_{\pi}^{2}$ after
$t\approx 2$ fm/c, which is two orders of magnitude lower than the value found
here (see Fig. 5). However, the magnetic field predicted by [11] is around
$\sim 10^{-2}m_{\pi}^{2}$ after $t\approx 2$ fm/c, which is in agreement with
our model.

In Fig. 7 we show a space-time contour plot of the magnetic field. The field
at the higher collision energy (on the left) has a higher peak magnetic field.
For lower collision energy (on the right) the field is less Lorentz
contracted, and leads to a magnetic field at late times that is a factor of
$\sim 1.1$ larger. The freeze-out magnetic field will increase at lower
collision energy mainly due to the decreasing freeze-out time.

Phantom text

Figure 7: Space-time plot of the magnetic field on the beam axis ($x=y=0$) in
eternal (pre-existent) QGP with $T=300$ MeV for a Au-Au collision at impact
parameter $b=6.4\,$fm. Left panel: collision energy $\sqrt{s_{\text{NN}}}=200$
GeV; right panel: collision energy $\sqrt{s_{\text{NN}}}=17$ GeV. The same
value of $\kappa$ is used as in Fig. 3. In a more realistic scenario, where
the QGP is formed during the collision, the field would only create induced
currents in the upper light-cone.

##  V Summary and Conclusion

We have studied the linear response of a stationary, homogeneous QGP to the
electromagnetic field of two colliding nuclei in relativistic kinetic theory
with collisional damping. As an application, we calculated the magnetic field
between the two nuclei. We showed that the response to the external field is
controlled by the polarization function along the light-cone,
$\Pi^{\mu}_{\nu}(\omega,|\boldsymbol{k}|\approx\omega)$. This allowed us to
derive an approximate analytic solution for the magnetic field that takes into
account the dynamical medium response. We showed that the late-time magnetic
field is mainly determined by the static electric conductivity of the QGP.
Since the remnant magnetic field at hadronization does not depend strongly on
the collision energy an experimental measurement of the magnetic field at
different collision energies could permit a determination the electrical
conductivity of the QGP.

Our calculation can be improved in various ways to make it more realistic.
Such improvements would include a realistic space-time dependence of the
medium (formation and hydrodynamical evolution), nonzero net baryon density,
and viscous corrections to the unperturbed phase-space distribution used to
calculate the polarization tensor. These improvements will require numerical
solutions of the linear response equations. It would also be of interest to
study the energy-momentum deposition into the QGP by the external
electromagnetic field.

Beyond heavy-ion collisions, the homogeneous polarization tensor with
collisional damping studied in our work has applicability in cosmology, where
a QGP existed during the first $10~{}\mu$s of the early universe, and possibly
in astrophysics where quark matter can exist at the core of collapsed stars.
In these situations the assumption of homogeneity and stationarity of the
medium on the scale of the relevant parameters, $m_{D}$ and $\kappa$, is well
justified.

###### Acknowledgements.

An allocation of computer time from the UA Research Computing High Performance
Computing (HPC) at the University of Arizona is gratefully acknowledged. BM
acknowledges support by the U.S. Department of Energy Office of Science under
Grant DE-FG02-05ER41367.

##  Appendix A Electric current of two colliding nuclei

Here we define the free charge and current density used to describe heavy ion
collisions. We wish to model two nuclei moving at constant velocity $\pm\beta$
along the collision axis ($\hat{z}$ direction) that are offset by $\pm b/2$
within the collision plane ($\hat{x}$ direction). For simplicity we model the
charge distribution as a gaussian in all directions

$\rho_{\text{ext}\pm}(t,\boldsymbol{x})=\frac{Ze\gamma}{\pi^{3/2}R^{3}}e^{-\frac{1}{R^{2}}(x\mp
b/2)^{2}}e^{-\frac{1}{R^{2}}y^{2}}\\\ \times
e^{-\frac{\gamma^{2}}{R^{2}}(z\mp\beta t)^{2}}\,,$ (59)

where the normalization is chosen in such a way that

$\int\rho_{\text{ext}\pm}(t,\boldsymbol{x})d^{3}\boldsymbol{x}=Ze\,,$ (60)

is the total charge of the heavy ion nucleus and $\gamma$ is the usual
relativistic factor. The Gaussian radius parameter $R$ is related to the mean
squared radius of the nucleus, $\langle r^{2}\rangle$ at rest, ($\gamma=1$) by

$\langle r^{2}\rangle=\frac{1}{Ze}\int
r^{2}\rho_{\text{ext}\pm}(\boldsymbol{x})d^{3}\boldsymbol{x}=\frac{3}{2}R^{2}\,,$
(61)

which is measured experimentally for a gold nucleus to be $\sqrt{\langle
r^{2}\rangle}\approx 5.30\,$fm [30] .

At time $t=0$ both nuclei are localized at the $z=0$ plane and we assume that
before and after the collision they continue moving on a straight line along
the $z$-axis. The gaussian form of the charge distributions allows us to
evaluate the Fourier transformations easily. The transforms in the transverse
directions are

$\displaystyle\int_{-\infty}^{\infty}dy\ e^{-ik_{y}y}e^{-y^{2}/R^{2}}$
$\displaystyle=R\sqrt{\pi}e^{-k_{y}^{2}R^{2}/4}\,,$ (62)
$\displaystyle\int_{-\infty}^{\infty}\ dxe^{-ik_{x}x}e^{-(x\mp
b/2)^{2}/R^{2}}$ $\displaystyle=R\sqrt{\pi}e^{-k_{x}^{2}R^{2}/4}e^{\pm
ik_{x}b/2}\,.$ (63)

The last two integrals are a bit more complicated because they are coupled

$\int_{-\infty}^{\infty}e^{i\omega
t}\left(\int_{-\infty}^{\infty}e^{-ik_{z}z}e^{-\frac{\gamma^{2}}{R^{2}}(z^{2}\pm
2z\beta t)}dz\right)\\\ \times
e^{-\frac{\gamma^{2}}{R^{2}}\beta^{2}t^{2}}dct=\frac{R\sqrt{\pi}}{\gamma}e^{-\frac{k_{z}^{2}R^{2}}{4\gamma^{2}}}\int_{-\infty}^{\infty}e^{i(\omega\pm
k_{z}\beta)t}dt\\\
=\frac{2R\pi^{3/2}}{\gamma}e^{-\frac{k_{z}^{2}R^{2}}{4\gamma^{2}}}\delta(\omega\pm
k_{z}\beta)\,,$ (64)

where delta function appears because both nuclei move at a constant velocity.
Altogether the Fourier transformed charge distributions are

$\widetilde{\rho}_{\text{ext}\pm}(\omega,\boldsymbol{k})=2\pi
Ze\,e^{-(k_{x}^{2}+k_{y}^{2}+k_{z}^{2}/\gamma^{2})\frac{R^{2}}{4}}\\\ \times
e^{\mp\frac{ik_{x}b}{2}}\delta(\omega\mp k_{z}\beta)\,,$ (65)

which may be written in cylindrical coordinates

$\widetilde{\rho}_{\text{ext}\pm}(\omega,\boldsymbol{k})=2\pi
Ze\,e^{-(k_{\rho}^{2}+k_{z}^{2}/\gamma^{2})\frac{R^{2}}{4}}\\\ \times
e^{\mp\frac{ik_{\rho}b\cos\theta}{2}}\delta(\omega\mp k_{z}\beta)\,.$ (66)

The current densities are obtained from

$\widetilde{\boldsymbol{j}}_{\text{ext}\pm}(\omega,\boldsymbol{k})=\pm\beta\hat{\boldsymbol{z}}\widetilde{\rho}_{\text{ext}\pm}(\omega,\boldsymbol{k})\,.$
(67)

The transverse component of the current given by

$\widetilde{\boldsymbol{j}}_{\perp,\text{ext}}=\widetilde{\boldsymbol{j}}_{\text{ext}}-(\hat{\boldsymbol{k}}\cdot\widetilde{\boldsymbol{j}}_{\text{ext}})\hat{\boldsymbol{k}}\\\
=(\hat{\boldsymbol{z}}-\hat{k}_{z}\hat{\boldsymbol{k}})\beta(\widetilde{\rho}_{\text{ext}+}-\widetilde{\rho}_{\text{ext}-})\,.$
(68)

##  Appendix B Magnetic field at the collision center

The magnetic field inside the plasma is given in Fourier space by

$\widetilde{\boldsymbol{B}}=i\boldsymbol{k}\times\widetilde{\boldsymbol{A}}=i\boldsymbol{k}\times\widetilde{\boldsymbol{A}}_{\perp}\,,$
(69)

where the potential $\widetilde{\boldsymbol{A}}$ has been projected into
components transverse and longitudinal to $\boldsymbol{k}$. In the following
we represent the wave-vector in cylindrical coordinates
$\boldsymbol{k}=(k_{\rho}\cos\theta,k_{\rho}\sin\theta,k_{z})$. Using the
expression for the self-consistent vector potential Eq. (20) we find for the
magnetic field,

$\widetilde{\boldsymbol{B}}=\frac{\mu_{0}i\boldsymbol{k}\times\widetilde{\boldsymbol{j}}_{\perp\text{ext}}}{\boldsymbol{k}^{2}-\omega^{2}-\mu_{0}\Pi_{\perp}}\,.$
(70)

Given the definition of $\widetilde{\boldsymbol{j}}_{\perp,\text{ext}}$ in Eq.
(68) we can replace the perpendicular component of the current by its full
form, adding the $\pm$ components of Eq. (67):

$\widetilde{\boldsymbol{B}}=\frac{\mu_{0}i\boldsymbol{k}\times\widetilde{\boldsymbol{j}}_{\text{ext}}}{\boldsymbol{k}^{2}-\omega^{2}-\mu_{0}\Pi_{\perp}}\,.$
(71)

We now Fourier transform this quantity back to position space in order to
calculate the magnetic field at the collision center as a function of time.
Due to symmetry the only nonvanishing component of the magnetic field at this
location will be the $y$-component:

$\widetilde{B}_{y}=\mu_{0}\frac{ik_{x}\beta(\widetilde{\rho}_{\text{ext}-}-\widetilde{\rho}_{\text{ext}+})}{\boldsymbol{k}^{2}-\omega^{2}-\mu_{0}\Pi_{\perp}}\,.$
(72)

The Fourier transform at any point along the collision axis ($x=y=0$) is given
by

$B_{y}(z,t)=\int\frac{d^{4}k}{(2\pi)^{4}}e^{-i\omega
t+ik_{z}z}\widetilde{B}_{y}(\omega,\boldsymbol{k})\,.$ (73)

In cylindrical coordinates the integral can be written as

$B_{y}(z,t)=\frac{1}{(2\pi)^{4}}\int k_{\rho}dk_{\rho}d\omega dk_{z}d\theta
e^{-i\omega t+ik_{z}z}\\\ \mu_{0}\frac{i\beta
k_{\rho}\cos\theta(\widetilde{\rho}_{f-}-\widetilde{\rho}_{f+})}{\boldsymbol{k}^{2}-\omega^{2}-\mu_{0}\Pi_{\perp}(\omega,|\boldsymbol{k}|)}\,.$
(74)

We can use the the delta function in the Fourier transformed current Eq. (66)
to trivially perform the $k_{z}$ integral:

$B_{y}(z,t)=-\mu_{0}\frac{Ze\beta}{(2\pi)^{3}}\int dk_{\rho}d\omega d\theta\
e^{-i\omega t}\\\
\frac{2k_{\rho}^{2}\cos\theta\sin\left(k_{\rho}\cos\theta\frac{b}{2}-\frac{\omega
z}{\beta}\right)e^{-(k_{\rho}^{2}+\omega^{2}/(\beta\gamma)^{2})\frac{R^{2}}{4}}}{\omega^{2}/(\gamma\beta)^{2}+k_{\rho}^{2}-\mu_{0}\Pi_{\perp}\left(\omega,\sqrt{k_{\rho}^{2}+\omega^{2}/\beta^{2}}\right)}\,.$
(75)

We next perform the angular integration

$B_{y}(z,t)=-\mu_{0}\frac{Ze\beta}{(2\pi)^{2}}\int dk_{\rho}d\omega\
e^{-i\omega t}\\\
\frac{2k_{\rho}^{2}J_{1}\left(\frac{k_{\rho}b}{2}\right)\cos\left(\frac{\omega
z}{\beta}\right)e^{-(k_{\rho}^{2}+\omega^{2}/(\beta\gamma)^{2})\frac{R^{2}}{4}}}{\omega^{2}/(\gamma\beta)^{2}+k_{\rho}^{2}-\mu_{0}\Pi_{\perp}\left(\omega,\sqrt{k_{\rho}^{2}+\omega^{2}/\beta^{2}}\right)}\,,$
(76)

where $J_{1}$ is a Bessel function of the first kind. The remaining integrals
have to be performed numerically. From here forward we pull the factor of
$\mu_{0}$ into the Debye mass $m_{D}^{2}$, such that the factor of $e^{2}$
goes to $4\pi\alpha$ in Eq. (34).

We can obtain an analytical expression for the Drude approximation Eq. (40) in
the limit $\gamma\beta\gg\sqrt{\kappa/\sigma_{0}}$, which is valid when
$\gamma\gg 12$ for the values of $\sigma_{0}$ and $\kappa$ adopted here. In
this limit we can neglect the first term in the denominator of Eq. (76), which
now takes the simple form

$k_{\rho}^{2}-i\omega\frac{\omega_{p}^{2}}{\kappa-i\omega}\,.$ (77)

Note that using Eq. (41) we can see
$\omega_{p}^{2}=m_{D}^{2}/3=\kappa\sigma_{0}$. The integrand of Eq. (76) then
has a single pole at

$\omega=-i\frac{k_{\rho}^{2}\kappa}{k_{\rho}^{2}+\omega_{p}^{2}}\,,$ (78)

and the frequency integral can be performed by contour integration in the
lower complex plane. Consistently neglecting the term proportional to
$\omega^{2}/(\beta\gamma)^{2}$ in the exponent, the integration yields:

$B_{y}(z,t)\approx-\mu_{0}\frac{Ze\beta}{2\pi}\int dk_{\rho}\,2\kappa
k_{\rho}^{2}\\\
\frac{\omega_{p}^{2}}{(k_{\rho}^{2}+\omega_{p}^{2})^{2}}J_{1}\left(\frac{k_{\rho}b}{2}\right)e^{-k_{\rho}^{2}R^{2}/4}\\\
\cosh\left(\frac{k_{\rho}^{2}\kappa}{k_{\rho}^{2}+\omega_{p}^{2}}\frac{z}{\beta}\right)\exp\left(-\frac{k_{\rho}^{2}\kappa
t}{k_{\rho}^{2}+\omega_{p}^{2}}\right)\,.$ (79)

For late times $t$ the exponential factor only samples the small $k_{\rho}$
region of the integrand ($k_{\rho}^{2}<\sigma_{0}/t$). We can then neglect
$k_{\rho}$ with respect to $\omega_{p}$ in the integrand provided that
$\sigma_{0}/t\ll\omega_{p}^{2}$, which is satisfied when $t\gg
1/\kappa=t_{\text{rel}}\approx 1\,\text{fm}/c$. The expression then takes the
simplified form: yielding

$B_{y}(z,t)\approx-\mu_{0}\frac{Ze\beta}{2\pi}\int
dk_{\rho}\,\frac{2k_{\rho}^{2}}{\sigma_{0}}J_{1}\left(\frac{k_{\rho}b}{2}\right)\\\
e^{-k_{\rho}^{2}R^{2}/4}\cosh\left(\frac{k_{\rho}^{2}}{\sigma_{0}}\frac{z}{\beta}\right)e^{-k_{\rho}^{2}t/\sigma_{0}}\,.$
(80)

The integral over $k_{\rho}$ can now be performed analytically resulting in:

$B_{y}(z,t)\approx-\mu_{0}\frac{Ze\beta}{2\pi}\frac{b}{8\sigma_{0}}\left[\frac{e^{-\frac{b^{2}}{16L_{+}}}}{L_{+}^{2}}+\frac{e^{-\frac{b^{2}}{16L_{-}}}}{L_{-}^{2}}\right]$
(81)

with

$L_{\pm}=\frac{R^{2}}{4}+\frac{t\pm z/\beta}{\sigma_{0}}\,.$ (82)

At the collision center $(z=0)$ and for $t\gg\sigma_{0}R^{2}/4$ our result
simplifies to

$B_{y}(0,t)\approx-\mu_{0}\frac{Ze\beta}{2\pi}\frac{b\sigma_{0}}{4t^{2}}e^{-\frac{\sigma_{0}b^{2}}{16t}}\,.$
(83)

This result differs from Tuchin’s [3] by a factor 1/4 in the exponent due to
their convention for impact parameter $b\rightarrow 2b$.

## References

  * [1] K. Tuchin, Phys. Rev. C 82, 034904 (2010) [erratum: Phys. Rev. C 83 (2011), 039903] [arXiv:1006.3051 [nucl-th]].
  * [2] W. T. Deng and X. G. Huang, Phys. Rev. C 85, 044907 (2012) [arXiv:1201.5108 [nucl-th]].
  * [3] K. Tuchin, Phys. Rev. C 88, 024911 (2013) [arXiv:1305.5806 [hep-ph]].
  * [4] L. McLerran and V. Skokov, Nucl. Phys. A 929, 184 (2014) [arXiv:1305.0774 [hep-ph]].
  * [5] U. Gursoy, D. Kharzeev and K. Rajagopal, Phys. Rev. C 89, 054905 (2014) [arXiv:1401.3805 [hep-ph]].
  * [6] V. Roy, S. Pu, L. Rezzolla and D. Rischke, Phys. Lett. B 750, 45 (2015) [arXiv:1506.06620 [nucl-th]].
  * [7] H. Li, X. l. Sheng and Q. Wang, Phys. Rev. C 94, 044903 (2016) [arXiv:1602.02223 [nucl-th]].
  * [8] G. Inghirami, L. Del Zanna, A. Beraudo, M. H. Moghaddam, F. Becattini and M. Bleicher, Eur. Phys. J. C 76, 659 (2016) [arXiv:1609.03042 [hep-ph]].
  * [9] G. Inghirami, M. Mace, Y. Hirono, L. Del Zanna, D. E. Kharzeev and M. Bleicher, Eur. Phys. J. C 80, 293 (2020) [arXiv:1908.07605 [hep-ph]].
  * [10] L. Yan and X. G. Huang, [arXiv:2104.00831 [nucl-th]].
  * [11] Z. Wang, J. Zhao, C. Greiner, Z. Xu and P. Zhuang, [arXiv:2110.14302 [hep-ph]].
  * [12] M. Formanek, C. Grayson, J. Rafelski and B. Müller, Annals Phys. 434, 168605 (2021) [arXiv:2105.07897 [physics.plasm-ph]].
  * [13] W. Florkowski, M. P. Heller and M. Spalinski, Rept. Prog. Phys. 81, 046001 (2018) [arXiv:1707.02282 [hep-ph]].
  * [14] G. S. Rocha, G. S. Denicol and J. Noronha, Phys. Rev. Lett. 127, 042301 (2021) [arXiv:2103.07489 [nucl-th]].
  * [15] P. L. Bhatnagar, E. P. Gross and M. Krook, Phys. Rev. 94, 511 (1954).
  * [16] H. Song and U. W. Heinz, Phys. Rev. C 77, 064901 (2008) [arXiv:0712.3715 [nucl-th]].
  * [17] R. Starke and G. A. H. Schober, Int. J. Mod. Phys. D 25, 1640010 (2016) [arXiv:1409.3723 [math-ph]].
  * [18] H. A. Weldon, Phys. Rev. D 26, 1394 (1982).
  * [19] D. Melrose, Quantum Plasmadynamics: Unmagnetized Plasmas, Lect. Notes Phys. 735 (Springer, New York, 2008).
  * [20] P. Romatschke, Eur. Phys. J. C 76, 352 (2016) [arXiv:1512.02641 [hep-th]].
  * [21] J. I. Kapusta, Phys. Rev. D 46, 4749 (1992)
  * [22] G. Jackson, Phys. Rev. D 100, 116019 (2019) [arXiv:1910.07552 [hep-ph]].
  * [23] S. Mrówczyński, Acta Phys. Polon. B 19, 91 (1988).
  * [24] J. Letessier, and J. Rafelski, Hadrons and Quark-Gluon Plasma ( Cambridge University Press, 2002)
  * [25] P. Drude, Ann. Phys. (Leipzig) 1, 566 (1900).
  * [26] D. E. Kharzeev, L. D. McLerran and H. J. Warringa, Nucl. Phys. A 803, 227 (2008) [arXiv:0711.0950 [hep-ph]].
  * [27] S. A. Bass and A. Dumitru, Phys. Rev. C 61, 064909 (2000) [arXiv:nucl-th/0001033 [nucl-th]].
  * [28] B. Müller and A. Schäfer, Phys. Rev. D 98, 071902 (2018) [arXiv:1806.10907 [hep-ph]].
  * [29] J. Letessier, A. Tounsi, U. W. Heinz, J. Sollfrank and J. Rafelski, Phys. Rev. Lett. 70, 3530 (1993) [arXiv:hep-ph/9711349 [hep-ph]].
  * [30] H. De Vries, C. W. De Jager and C. De Vries, Atom. Data Nucl. Data Tabl. 36, 495 (1987)