# Debye Screening of Non-Abelian Plasmas in Curved Spacetimes

Elba Alonso-Monsalve<EMAIL_ADDRESS>David I. Kaiser<EMAIL_ADDRESS>Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA

###### Abstract

Decades of analytic and computational work have demonstrated that a charge
immersed in a hot plasma is screened. For both Abelian and non-Abelian
interactions, the characteristic screening length $1/m_{D}$ is set by the so-
called Debye mass $m_{D}\sim g_{s}T$, proportional to the plasma temperature
$T$ and the dimensionless gauge coupling $g_{s}$. One of the most interesting
naturally occurring examples is the quark-gluon plasma (QGP) that filled the
early universe prior to the QCD confinement phase transition at $t_{\rm
QCD}\sim 10^{-5}\,{\rm s}$. During this early epoch, regimes of strong
spacetime curvature are of significant cosmological interest, such as near
primordial black holes (PBHs). However, the typical description of Debye
screening only applies within Minkowski spacetime, and is therefore
insufficient to describe the dynamics of charged plasmas near PBHs or other
primordial features. We construct an effective field theory for soft modes of
the gauge field $A_{\mu}^{a}$ to give a full description of Debye screening in
non-Abelian plasmas within arbitrary curved spacetimes, recovering a
temperature-dependent Debye mass that exhibits gravitational redshift. We then
apply our results to some scenarios of cosmological interest: an expanding
FLRW universe and the vicinity of a PBH immersed in a hot QGP.

##  I Introduction

The screening of charges in hot plasmas, for both Abelian and non-Abelian
interactions, has been a topic of significant interest for decades. Upon
resumming hard thermal loops, the leading-order effects yield a characteristic
screening length $\lambda_{D}(T)=1/m_{D}(T)$ set by the Debye mass
$m_{D}(T)\sim g_{s}T$, where $g_{s}$ is the dimensionless gauge coupling
strength and $T$ is the temperature of the plasma. (For reviews, see Refs.
Blaizot and Iancu (2002); Litim and Manuel (2002); Kapusta and Gale (2006);
Bazavov and Weber (2021).)

The Debye mass $m_{D}(T)$ sets the characteristic scale for dynamics of
collective excitations within the plasma. For example, collective oscillations
of “soft” modes, with momenta $k_{\rm soft}\sim m_{D}(T)\ll T$, propagate with
frequency set by the plasma frequency $\omega_{p}\simeq m_{D}(T)$ Matinyan
_et al._ (1981); Klimov (1981, 1982); Weldon (1982); Braaten and Pisarski
(1990a, b, 1992); Blaizot and Iancu (1994a, 1993, b); Blaizot  _et al._
(1995); Mrowczynski  _et al._ (2017). Moreover, spatial regions of uniform
color-charge form within the plasma, with coherence length $1/m_{D}(T)$
Blaizot and Iancu (1994c, b). These effects have been studied extensively in
Minkowski spacetime, using both analytic techniques and lattice simulations
Blaizot and Iancu (2002); Litim and Manuel (2002); Kapusta and Gale (2006);
Bazavov and Weber (2021).

In this paper we consider Debye screening in hot plasmas within curved
spacetimes, with particular interest in cosmological applications. For
example, post-inflation reheating typically yields a universe filled with
Standard Model particles in thermal equilibrium at a temperature as high as
$T\sim{\cal O}(10^{14}\,{\rm GeV})\sim{\cal O}(10^{27}\,K)$ by the time
$t_{\rm therm}\sim{\cal O}(10^{-35}\,{\rm s})$ after the big bang Amin  _et
al._ (2014); Allahverdi  _et al._ (2021). Such temperatures are exponentially
greater than the QCD confinement scale $\Lambda_{\rm QCD}\simeq 0.170\,{\rm
GeV}$. Hence at such early times, the universe was filled with a plasma of
unconfined quarks and gluons, subject to the non-Abelian dynamics of QCD
Mukhanov (2005); Boyanovsky  _et al._ (2006).

As the universe expanded, the temperature of the plasma fell adiabatically:
$T(t)a(t)\simeq{\rm constant}$, where $a(t)$ is the scale factor of the
Friedmann-Lemaître-Robertson-Walker (FLRW) line-element. During the radiation-
dominated phase, $a(t)=(t/t_{\rm therm})^{1/2}$. The QCD confinement
transition, after which quarks and gluons remained bound in color-neutral
hadronic states, occurred at $t_{\rm QCD}\simeq 10^{-5}\,{\rm s}$, when
$T(t_{\rm QCD})\simeq\Lambda_{\rm QCD}$ Mukhanov (2005); Boyanovsky  _et al._
(2006).

Within the range of times $t_{\rm therm}\ll t\ll t_{\rm QCD}$, additional
phenomena of cosmological interest could have occurred. For example, if a
significant population of primordial black holes formed at early times, with
masses in the range $10^{17}\,{\rm g}\leq M(t_{c})\leq 10^{22}\,{\rm g}$, they
could account for the entire abundance of dark matter today Carr and Kühnel
(2020); Green and Kavanagh (2021); Villanueva-Domingo  _et al._ (2021);
Escrivà  _et al._ (2022). Given the dependence of the black holes’ masses on
the time of collapse $t_{c}$, the black-hole mass range of interest to account
for the present dark-matter abundance corresponds to formation times
$10^{-21}\,{\rm s}\leq t_{c}\leq 10^{-16}\,{\rm s}$. At these early times, the
temperature of the plasma within which the primordial black holes formed would
have been $10^{5}\,{\rm GeV}\leq T(t_{c})\leq 10^{7}\,{\rm GeV}$, for which
$T(t_{c})\gg T(t_{\rm QCD})$.

With such cosmological applications in mind, we study Debye screening in hot
non-Abelian plasmas within arbitrary curved spacetimes, including spacetimes
that need not be homogeneous and isotropic (such as in the vicinity of a
primordial black hole). In such general spacetimes, the effective temperature
of the plasma can develop spatial gradients, $T\rightarrow T(x)$, akin to the
familiar Tolman temperature Tolman (1930); Tolman and Ehrenfest (1930);
Santiago and Visser (2018). We identify corrections to the induced current
$j_{\mu}^{a}(x)$ and to the Debye mass $m_{D}(T(x))$ from spacetime curvature.
Much as in Minkowski spacetime, the components $A^{a}_{0}(x)$ of the gauge
field acquire an effective mass in the plasma proportional to $m_{D}(T(x))$,
whereas components $A^{a}_{i}(x)$ remain massless. (See, e.g., Refs. Holcomb
and Tajima (1989); Dettmann  _et al._ (1993); Prokopec and Woodard (2004);
Prokopec and Puchwein (2004); Dodin and Fisch (2010); Popov (2018) on the
behavior of Abelian plasmas within isotropic and homogeneous spacetimes.)

In Section II, we introduce an effective theory for soft modes within the
plasma, generalizing the elegant formalism reviewed in Refs. Blaizot and Iancu
(2002); Litim and Manuel (2002) for application to curved spacetimes. In
Section III we evaluate the induced current $j_{\mu}^{a}(x)$ for the soft
modes that arises from interactions with high-energy quarks and gluons in the
plasma. Section IV considers corrections to the usual expression for the Debye
mass $m_{D}(T)$ arising from spacetime curvature. In Section V we apply these
expressions for $j_{\mu}^{a}(x)$ and $m_{D}(T(x))$ to some scenarios of
cosmological interest.

We restrict attention to $(3+1)$ spacetime dimensions and adopt the metric
signature $(-,+,+,+)$. Greek letters $\mu,\nu=0,1,2,3$ label spacetime
indices, while Latin letters $i,j=1,2,3$ label spatial indices. The color-
charge indices associated with the adjoint representation of the gauge group
${\rm SU}(N_{c})$ range over $a,b,c=1,2,...,N_{c}^{2}-1$ and are raised and
lowered with $\delta_{ab}$. The generators $T^{a}$ of the Lie algebra for
${\rm SU}(N_{c})$ satisfy $[T^{a},T^{b}]=if^{abc}T^{c}$, where $f^{abc}$ are
the totally antisymmetric structure constants, and the generators are
normalized such that $2\,{\rm Tr}(T^{a}T^{b})=\delta^{ab}$. We adopt natural
units in which $c=\hbar=k_{B}=1$, in terms of which the reduced Planck mass is
given by $M_{\rm pl}\equiv 1/\sqrt{8\pi G}=2.43\times 10^{18}$ GeV.

##  II Effective Theory for Soft Modes

At early times $t\ll t_{\rm QCD}$, when the temperature of the plasma filling
the universe satisfied $T\gg\Lambda_{\rm QCD}$, the number density of gluon
degrees of freedom likely dominated those of quarks: gluons can radiate gluons
at tree level, and bosonic statistics allow gluon number densities per mode to
scale as $n_{k}\sim 1/\alpha_{s}$, where $\alpha_{s}\equiv g_{s}^{2}/(4\pi)$
and $g_{s}$ is the QCD coupling constant Gyulassy and McLerran (2005). Given
the running of $\alpha_{s}$ to lower values at higher energies, the number
density of gluons could have greatly exceeded the number density of quarks
within the plasma at early times Blaizot and Iancu (2002); Bazavov and Weber
(2021). Similar behavior has been observed in the plasmas that form
immediately after relativistic heavy-ion collisions Scardina  _et al._ (2013);
Gyulassy and McLerran (2005). Hence we expect gluons to dominate the dynamics,
and consider the effective action for a Yang-Mills model with gauge group
${\rm SU}(N_{c})$ and field strength tensor

$F_{\mu\nu}^{a}=\nabla_{\mu}A_{\nu}^{a}-\nabla_{\nu}A_{\mu}^{a}+g_{s}f^{abc}A_{\mu}^{b}A_{\nu}^{c},$
(1)

where $\nabla_{\mu}$ is the (spacetime) covariant derivative associated with
metric $g_{\mu\nu}(x)$.

Following the approach reviewed in Refs. Blaizot and Iancu (2002); Litim and
Manuel (2002), we construct an effective theory for long-wavelength
excitations in the plasma within the high-temperature limit. In particular, we
consider the effective action for soft modes $\tilde{A}_{\mu}^{a}$ with
typical momenta $k_{\rm soft}\sim g_{s}T$. Thermal and quantum corrections to
the behavior of the soft modes are dominated by quanta $a_{\mu}^{a}$ with
momenta $k_{\rm hard}\sim T$. We therefore write
$A_{\mu}^{a}=\tilde{A}_{\mu}^{a}+a_{\mu}^{a}$, integrate out the high-momentum
modes $a_{\mu}^{a}$, and drop the tilde on the soft modes
$\tilde{A}_{\mu}^{a}\rightarrow A_{\mu}^{a}$.

As described in Refs. Blaizot and Iancu (2002); Litim and Manuel (2002), the
soft modes $A_{\mu}^{a}$ have large occupation numbers per mode and hence
behave as effectively classical fields. The behavior of the soft modes
therefore yields the mean-field dynamics for the system on length-scales
$\lambda\geq 1/k_{\rm soft}\sim 1/(g_{s}T)$. We follow the background-field
method of Ref. Blaizot and Iancu (2002), whereby we choose to decompose the
action of the gauge symmetry on $a_{\mu}^{a}$ and $A_{\mu}^{a}$ in such a way
that leaves the soft modes $A_{\mu}^{a}$ unchanged. One example is to use a
generalized Coulomb-type gauge-fixing term in the full theory,
$\nabla_{i}a_{a}^{i}-g_{s}f^{abc}A_{i}^{b}a_{c}^{i}$ Blaizot and Iancu (2002).
No additional gauge-fixing terms or ghosts are then required in the effective
action for the soft modes.

For long length-scales $\lambda\geq 1/k_{\rm soft}\gg 1/k_{\rm hard}\sim 1/T$,
the dynamics may be described by the effective action

$S_{\rm eff}=\int d^{4}x\sqrt{-g}\left[\frac{M_{\rm
pl}^{2}}{2}R-\frac{1}{4}F_{\mu\nu}^{a}F_{a}^{\mu\nu}-j_{\mu}^{a}A_{a}^{\mu}+{\cal
L}_{\rm fluid}\right],$ (2)

where $j_{\mu}^{a}(x)$ is the induced current generated by (non-Abelian) self-
interactions with the high-momentum modes $a_{\mu}^{a}(x)$; the induced
current $j_{\mu}^{a}$ also includes contributrions from high-momentum quarks.
The term ${\cal L}_{\rm fluid}$ represents the contributions to the evolution
of the spacetime curvature from constituents other than the soft modes. For
example, the high-momentum modes in the plasma (coarse-grained over a length-
scale $\lambda\gg 1/k_{\rm hard}$) behave as a neutral fluid with a radiation-
dominated equation of state, and thereby contribute to the time evolution of
the scale factor $a(t)$ in an FLRW background, while masses $M$ of compact
objects, such as black holes, influence the spacetime curvature in their
vicinity. Each of these contributions affects $g_{\mu\nu}(x)$ and hence
impacts the dynamics of the soft modes $A_{\mu}^{a}(x)$.

Varying $S_{\rm eff}$ with respect to $A_{\mu}^{a}$ yields the equations of
motion

$\begin{split}D^{\mu}F^{a}_{\mu\nu}=j_{\nu}^{a},\end{split}$ (3)

where the covariant derivative $D_{\mu}$ acting on a spacetime tensor that
transforms in the adjoint representation of ${\rm SU}(N_{c})$,
$X^{a}_{\nu_{1}\cdots\nu_{n}}$, is defined as

$D_{\mu}X^{a}_{\nu_{1}\cdots\nu_{n}}\equiv\nabla_{\mu}X^{a}_{\nu_{1}\cdots\nu_{n}}+g_{s}f^{abc}A^{b}_{\mu}X^{c}_{\nu_{1}\cdots\nu_{n}}.$
(4)

The left-hand side of Eq. (3) can be expanded as

$\begin{split}D^{\mu}F^{a}_{\mu\nu}&=\nabla^{\mu}F^{a}_{\mu\nu}+g_{s}f^{abc}A^{\mu}_{b}F^{c}_{\mu\nu}\\\
&=D^{\mu}D_{\mu}A_{\nu}^{a}-g_{s}f^{abc}A^{\mu}_{b}\nabla_{\nu}A_{\mu}^{c}-R_{\mu\nu}A^{\mu}_{a},\end{split}$
(5)

where $R_{\mu\nu}$ is the Ricci tensor for the background spacetime. The
equations of motion imply that the induced current $j_{\mu}^{a}$ must be
covariantly conserved with respect to both the spacetime curvature and the
${\rm SU}(N_{c})$ gauge group:

$D^{\mu}j_{\mu}^{a}(x)=0,$ (6)

for each color $a$.

The contribution to the energy-momentum tensor from the Yang-Mills field is
given, as usual, by

$T^{(A)}_{\mu\nu}=F^{a}_{\mu\lambda}F_{\nu}^{a\>\lambda}-\frac{1}{4}g_{\mu\nu}F_{a}^{\lambda\sigma}F^{a}_{\lambda\sigma},$
(7)

while the additional fluid contributes

$T^{\rm fluid}_{\mu\nu}=-\frac{2}{\sqrt{-g}}\functionalderivative{\cal L_{\rm
fluid}}{g^{\mu\nu}}.$ (8)

These terms satisfy the usual covariant conservation equations

$\begin{split}\nabla_{\mu}T^{\mu\nu}_{(A)}&=j_{\mu}^{a}F^{\mu\nu}_{a},\\\
\nabla_{\mu}T^{\mu\nu}_{\rm fluid}&=0\end{split}$ (9)

in the presence of the induced current $j_{\mu}^{a}$.

##  III Induced Current

The induced current $j_{\mu}^{a}(x)$ for the soft modes $A_{\mu}^{a}(x)$
arises from interactions with high-momentum quarks and gluons in the plasma.
Lattice simulations have confirmed the analytic expectation that for
temperatures $T\gg T_{\rm QCD}$, high-energy quarks and gluons attain an
equilibrium equation of state akin to that of a gas of non-interacting,
massless particles. In particular, the so-called “trace anomaly,”
$(\rho-3P)/T^{4}$, where $\rho$ is the energy density and $P$ the pressure of
particles in the plasma, tends rapidly toward zero for $T\gg T_{\rm QCD}$, and
therefore $\rho\approx 3P$. (See, e.g., Refs. Borsanyi  _et al._ (2012, 2014);
Bazavov  _et al._ (2014); Pasechnik and Šumbera (2017).)

At high temperatures, the soft quantum fields of the plasma are well
approximated by classical fields, due to Bose-Einstein enhancement from high
occupation numbers. Furthermore, the hard (high-momentum) modes can be
approximated by an ensemble of classical point particles to leading order in
the coupling $g_{s}$, since they are weakly interacting. Specifically, the
effects of the high-momentum particles on the soft modes within the plasma are
dominated by processes involving “hard thermal loops”: one-loop diagrams with
arbitrary numbers of low-momentum external legs (with $k_{\rm soft}\sim
g_{s}T$) and hard internal momenta (with $k_{\rm hard}\sim T$) Pisarski
(1989); Braaten and Pisarski (1990b, c, d, 1992). These effects can be
analyzed within a mean-field approximation of the (truncated) Schwinger-Dyson
equations Blaizot and Iancu (1994a, 1993, b); Blaizot  _et al._ (1995);
Blaizot and Iancu (2002); Mrowczynski  _et al._ (2017), or by simply adopting
classical transport equations for high-momentum particles Litim and Manuel
(2002); Kelly  _et al._ (1994a, b); Pisarski (1998); Mrowczynski  _et al._
(2017). In this section we adapt the latter approach to arbitrary curved
spacetimes. (See also Ref. Brandt  _et al._ (1995).)

The equations of motion for a classical point particle with color charge in a
soft gauge-field background are well known Wong (1970); Linden  _et al._
(1996); Kelly  _et al._ (1994b); Brandt  _et al._ (1995):

$\begin{split}\derivative{x^{\mu}}{\lambda}&=P^{\mu},\\\
\derivative{P_{\mu}}{\lambda}&=-\frac{1}{2}\partialderivative{g^{\alpha\beta}}{x^{\mu}}P_{\alpha}P_{\beta}-g_{s}Q^{a}F_{\mu\nu}^{a}P^{\nu},\\\
\derivative{Q^{a}}{\lambda}&=-g_{s}f^{abc}A_{\mu}^{b}Q^{c}P^{\mu}.\end{split}$
(10)

Here $x^{\mu}$ is the position of the particle,
$P^{\mu}=\derivative*{x^{\mu}}{\lambda}$ its kinetic 4-momentum, $\lambda$ an
affine parameter, and $Q^{a}$ its ${\rm SU}(N_{c})$ charge. We take advantage
of the freedom to linearly rescale the affine parameter $\lambda$ so that
$\derivative*{x^{\mu}}{\lambda}$ has units of energy. (We do not explicitly
rescale by the particle mass $m$ so that our formalism can be applied to
massless particles.) Notice that, due to non-Abelian self-interactions, the
charges $Q^{a}$ are dynamical, unlike in electromagnetism.

The induced current $j^{\mu}_{a}(x)$, summed over all hard particles, is given
by Blaizot and Iancu (2002); Litim and Manuel (2002)

$j^{\mu}_{a}(x)=g_{s}\int
dQ\,d\omega\,\derivative{x^{\mu}}{\lambda}Q_{a}\delta f(x,p,Q),$ (11)

where $\delta f(x,p,Q)$ represents the deviation from equilibrium of the
distribution function for the high-momentum charge-carrying particles. The
integration measure includes the momentum volume form $d\omega$ and the
measure for the space of color charges $dQ$, subject to the physical
constraints of on-shell mass condition, positivity of energy, and conservation
of the $N_{c}-1$ group Casimirs. The momentum measure is given below, in Eq.
(36), while the color-space measure is

$dQ=d^{3}Q\,c_{R}\,\delta(Q^{a}Q^{a}-q_{2})$ (12)

for the gauge group SU(2) ($N_{c}=2$) and

$dQ=d^{8}Q\,c_{R}\,\delta(Q^{a}Q^{a}-q_{2})\delta(f^{abc}Q^{a}Q^{b}Q^{c}-q_{3})$
(13)

for the gauge group SU(3) ($N_{c}=3$). The representation-dependent constant
$c_{R}$ is fixed by the normalization $\int dQ=1$, while the constants
$q_{2},q_{3}$ are further fixed by the first and second Casimirs,
respectively. Thus, the integration over color charges serves to enforce the
conservation of the Casimir invariants Litim and Manuel (2002); Kelly  _et
al._ (1994b, a).

We consider small departures from equilibrium, and therefore work
perturbatively in powers of the coupling $g_{s}$, so we may write the full
distribution function as

$f(x,p,Q)=f^{(0)}(x,p)+g_{s}f^{(1)}(x,p,Q)+{\cal O}(g_{s}^{2}).$ (14)

The relevant dynamics for the induced current $j^{\mu}_{a}(x)$ in Eq. (11) are
captured to leading order by identifying $\delta
f(x,p,Q)=g_{s}f^{(1)}(x,p,Q)$, as demonstrated by Refs. Kelly  _et al._
(1994b, a). The dynamical evolution of $\delta f(x,p,Q)$ is governed by the
collisionless Boltzmann equation for the full distribution function
$f(x,p,Q)$, which in the framework of Hamiltonian mechanics corresponds to the
conservation of $f(x,p,Q)$ along dynamical trajectories (Hamiltonian flow) in
phase space.

To solve the Boltzmann equation, it will be convenient to have an explicit
Hamiltonian formulation of the dynamics of the high-momentum particles. To
this end, we first construct an action which yields the correct equations of
motion, Eq. (10). This is nontrivial because there is no such action involving
the charges $Q^{a}$ directly as dynamical variables. To bypass this problem,
following Refs. Balachandran  _et al._ (1977); Barducci  _et al._ (1977);
Linden  _et al._ (1996), we introduce new dynamical variables $q^{a}$ which,
like $Q^{a}$, transform under the adjoint representation of SU($N_{c}$), but
which are anticommuting (Grassmann-valued). These new variables $q^{a}$ are
useful as an intermediate step, in terms of which we may construct the charges
$Q^{a}$ as

$Q^{a}=-\frac{i}{2}f^{abc}q^{b}q^{c}.$ (15)

An action for a single high-momentum particle that yields the appropriate
equations of motion may then be written in terms of the dynamical variables
$x^{\mu}$ and $q^{a}$ as

$\begin{split}S_{\rm 1p}=\int
d\lambda&\bigg{\\{}\frac{1}{2}g_{\mu\nu}\derivative{x^{\mu}}{\lambda}\derivative{x^{\nu}}{\lambda}+\frac{i}{2}q^{a}\derivative{q^{a}}{\lambda}\\\
&\quad-\frac{i}{2}g_{s}f^{abc}A_{\mu}^{a}q^{b}q^{c}\derivative{x^{\mu}}{\lambda}\bigg{\\}}.\end{split}$
(16)

The corresponding Hamiltonian is simply

$H=\frac{1}{2}g^{\mu\nu}P_{\mu}P_{\nu}.$ (17)

The kinetic momentum $P^{\mu}\equiv\derivative*{x^{\mu}}{\lambda}$ is related
to the canonical momentum $p_{\mu}$, conjugate to $x^{\mu}$, by

$P_{\mu}=p_{\mu}-g_{s}Q^{a}A_{\mu}^{a}.$ (18)

Then Hamilton’s equations are Eq. (10), as desired. Dynamical evolution in
phase space is generated by the Liouville vector field $X_{H}$:

$\begin{split}X_{H}\equiv\derivative{\lambda}&=\derivative{x^{\mu}}{\lambda}\partialderivative{x^{\mu}}+\derivative{p_{\mu}}{\lambda}\partialderivative{p_{\mu}}+\derivative{Q^{a}}{\lambda}\partialderivative{Q^{a}}\\\
&=P^{\mu}\partialderivative{x^{\mu}}\bigg{|}_{P}-\frac{1}{2}\partialderivative{g^{\alpha\beta}}{x^{\mu}}\bigg{|}_{P}P_{\alpha}P_{\beta}\partialderivative{P_{\mu}}\\\
&\quad\quad-g_{s}Q^{a}P^{\mu}F_{\mu\nu}^{a}\partialderivative{P_{\nu}}\\\
&\quad\quad-
g_{s}f^{abc}A^{b}_{\mu}Q^{c}P^{\mu}\partialderivative{Q^{a}}\bigg{|}_{P},\end{split}$
(19)

where, in the second line, we have written the vector field in noncanonical
coordinates $\\{x^{\mu},P_{\mu},Q^{a}\\}$, and we have noted explicitly which
derivatives are taken keeping $P_{\mu}$ fixed, as opposed to $p_{\mu}$.

The collisionless Boltzmann equation for the distribution function $f(x,p,Q)$
can be written in terms of the Liouville vector field as $X_{H}[f]=0$, and is
equivalently known as the Liouville equation. Given the expansion in Eq. (14),
we require that the Liouville equation be satisfied order by order in $g_{s}$:

$P^{\mu}\partialderivative{f^{(0)}}{x^{\mu}}\bigg{|}_{P}-\frac{1}{2}\partialderivative{g^{\alpha\beta}}{x^{\mu}}\bigg{|}_{P}P_{\alpha}P_{\beta}\partialderivative{f^{(0)}}{P_{\mu}}=0$
(20)

and

$\begin{split}P^{\mu}\partialderivative{f^{(1)}}{x^{\mu}}\bigg{|}_{P}-\frac{1}{2}\partialderivative{g^{\alpha\beta}}{x^{\mu}}\bigg{|}_{P}P_{\alpha}P_{\beta}\partialderivative{f^{(1)}}{P_{\mu}}=Q^{a}P^{\mu}F_{\mu\nu}^{a}\partialderivative{f^{(0)}}{P_{\nu}}.\end{split}$
(21)

We assume that all dependence of $f(x,p,Q)$ on the charges $Q^{a}$ can only
appear at ${\cal O}(g_{s})$ or above, and therefore $\partial f^{(0)}/\partial
Q^{a}=0$. Notice that, as emphasized in Ref. Pisarski (1998), the non-Abelian
terms in $X_{H}$, proportional to $g_{s}f^{abc}$, make no contribution to the
evolution of $f^{(1)}$, given our perturbative expansion in $g_{s}$.

We are free to choose any distribution function for our fluid, as long as it
solves the Boltzmann equation and its physical meaning is compatible with
small deviations from thermal equilibrium. For this reason, we aim to find a
solution of the collisionless Boltzmann equation $X_{H}[f]=0$ such that the
zeroth order in $g_{s}$ is of the usual form for a canonical ensemble,
$f^{(0)}\sim\exp[-\beta E]$, for a constant $\beta$ and a quantity $E$ that
may be interpreted as an energy. The equation $X_{H}[f]=0$ is typically solved
to ${\cal O}(g_{s})$ by employing Green’s function techniques Pisarski (1998);
Blaizot and Iancu (2002); Litim and Manuel (2002). We will instead show that
one can further exploit the Hamiltonian formalism to solve efficiently for the
distribution function to the same order. The two approaches are equivalent up
to ${\cal O}(g_{s})$, because quantum corrections to the classical equations
of motion in Eq. (10) only arise at ${\cal O}(g_{s}^{2})$ and above.

It is a classic result that Hamiltonian mechanics in the fully covariant phase
space $\\{x^{\mu},p_{\mu},Q^{a}\\}$ with evolution parameterized by an affine
$\lambda$ and generated by the Hamiltonian in Eq. (17) is equivalent to
Hamiltonian mechanics in the reduced phase space $\\{x^{i},p_{i},Q^{a}\\}$
with evolution parameterized by coordinate time $t\equiv x^{0}$ and generated
by the reduced Hamiltonian

$\bar{H}\equiv-p_{0}(t,x^{i},p_{i},Q^{a};h).$ (22)

This means that both formalisms yield the same equations of motion. This is
possible thanks to the redundancy in the choice of affine parameter $\lambda$,
as well as the conservation of $H(x^{\mu},p_{\mu},Q)$ along phase space
trajectories, in the covariant formalism. In Eq. (22) we have solved for
$p_{0}$ in terms of $\\{t,x^{i},p_{i},Q^{a},h\\}$ using
$H(x^{\mu},p_{\mu},Q)\equiv h$, with $h$ a constant. In the covariant
formalism, the conservation of $H(x^{\mu},p_{\mu},Q)$ follows from Eq. (10)
and the fact that $H$ cannot depend explicitly on $\lambda$. The reduced
phase-space formalism does not yield the conservation of
$H(x^{\mu},p_{\mu},Q^{a})$, so we must impose it by hand. For details of the
proof, see Sections 44 and 45 in Chapter 9 of Ref. Arnold (1978).

If the reduced Hamiltonian $\bar{H}$ from Eq. (22) does not depend explicitly
on $t$, then it is also conserved along phase-space trajectories,
$X_{H}[\bar{H}]=0$, so we may exploit $\bar{H}$ to construct a distribution
function. That is, we may use the fact that any scalar function of $\bar{H}$
will solve the collisionless Boltzmann equation to devise a valid distribution
function for our fluid. We therefore introduce a quasi-stationary
approximation: we consider only scenarios in which $\partial_{t}g_{\mu\nu}$
and $\partial_{t}A_{\mu}^{a}$ remain subdominant. This approximation is
appropriate, since we are interested in the behavior of the high-momentum
particles, whose dynamics are governed by the time-scale $1/k_{\rm hard}\sim
1/T$, whereas we expect the soft modes $A_{\mu}^{a}(x)$ to evolve on time-
scales set by $1/k_{\rm soft}\gg 1/k_{\rm hard}$. Likewise, our effective
description can only resolve dynamics up to scales set by $1/k_{\rm soft}$,
which bounds how sharply the background spacetime can evolve within our self-
consistent expansion as well. In particular, if we allowed $\bar{H}$ to have
an explicit dependence on $t$, then

$X_{H}[\bar{H}]=\frac{1}{2}\partial_{t}g^{\mu\nu}P_{\mu}P_{\nu}-g_{s}Q^{a}P^{\mu}\partial_{t}A_{\mu}^{a},$
(23)

and thus a distribution function constructed from $\bar{H}$ would yield self-
consistent dynamics as long as

$|\partial_{t}g^{\mu\nu}|,\frac{|\partial_{t}A_{\mu}^{a}|}{|A_{\mu}^{a}|}\ll
k_{\rm soft}.$ (24)

As we will see in Section V, this quasi-stationary approximation is easily
satisfied in many cosmological applications of interest.

We choose a distribution function $f=\exp[-\beta_{T}\bar{H}]$, where
$\beta_{T}$ is a constant Sato (2021). As explained just above, this satisfies
the collisionless Boltzmann equation by construction. We will see in the next
section that the zeroth-order term $f^{(0)}$ in the $g_{s}$ expansion
corresponds to a canonical ensemble (for fluid undergoing normal flow with
respect to coordinate time $t$), so it is a physically reasonable choice
because it represents thermal equilibrium. This term is

$f^{(0)}=e^{\beta_{T}P_{0}}.$ (25)

Given $f^{(0)}$, we may evaluate the equilibrium occupation numbers per mode
in the usual way,

$n_{\rm B,F}^{(0)}\equiv\frac{\sum n_{i}e^{n_{i}\beta_{T}P_{0}}}{\sum
e^{n_{i}\beta_{T}P_{0}}},$ (26)

where the sums run from $n_{i}=0$ to $n_{i}=\infty$ for bosons and to
$n_{i}=1$ for fermions. Thus we obtain the Bose-Einstein and Fermi-Dirac
distributions at zeroth-order:

$n_{\rm B,F}^{(0)}=\frac{1}{e^{-\beta_{T}P_{0}}\mp 1}.$ (27)

Substituting into Eq. (21) yields

$\begin{split}n_{\rm B,F}^{(1)}&=Q^{a}A_{0}^{a}\partialderivative{n_{\rm
B,F}^{(0)}}{P_{0}}\\\
&=-\beta_{T}\frac{e^{-\beta_{T}P_{0}}}{(e^{-\beta_{T}P_{0}}\mp
1)^{2}}Q^{a}A_{0}^{a}\end{split}$ (28)

for bosons and fermions, respectively. Summing over species and helicities, we
have $\delta f=g_{S}(2n_{\rm B}^{(1)}+4N_{f}n_{\rm F}^{(1)})$ in Eq. (11).

To evaluate $j^{\mu}_{a}(x)$ we first perform the $Q$ integral with the
measure in Eq. (13), which yields factors proportional to the index of the
representation: $\int dQ\,Q_{a}Q_{b}=C_{\rm B,F}\,\delta^{a}_{\>b}$, with
$C_{\rm B}=N_{c}$ for bosons and $C_{\rm F}=1/2$ for fermions. Then the
current becomes

$\begin{split}j^{\mu}_{a}(x)&=g_{s}^{2}\int dQ\,d\omega\left(2n_{\rm
B}^{(1)}+4N_{f}n_{\rm F}^{(1)}\right)Q^{a}P^{\mu}\\\
&=g_{s}^{2}A_{0}^{a}(x)\int d\omega\,{\cal N}(x,P)\,P^{\mu},\end{split}$ (29)

where we have defined

$\begin{split}{\cal N}(x,P)&\equiv 2N_{c}\partialderivative{n_{\rm
B}^{(0)}}{P_{0}}+2N_{f}\partialderivative{n_{\rm F}^{(0)}}{P_{0}}\\\
&=-2\beta_{T}\left[N_{c}\frac{e^{-\beta_{T}P_{0}}}{(e^{-\beta_{T}P_{0}}-1)^{2}}+N_{f}\frac{e^{-\beta_{T}P_{0}}}{(e^{-\beta_{T}P_{0}}+1)^{2}}\right].\end{split}$
(30)

Eq. (29) includes contributions from two polarization states for each
effectively massless particle, and $2N_{f}$ fermion species ($N_{f}$ each for
quarks and antiquarks). Note that, as in Minkowski spacetime Blaizot and Iancu
(2002), $j^{\mu}_{a}(x)$ is proportional to $A_{0}^{a}(x)$.

##  IV Debye Mass

In Minkowski spacetime, the induced current reduces (in the static limit) to
$j_{\mu}^{a}(x)=m_{D}^{2}A_{0}^{a}\delta_{\mu}^{0}$, effectively giving a
constant mass to the soft gluon components $A_{0}^{a}$, which is responsible
for Debye screening Blaizot and Iancu (2002); Kapusta and Gale (2006); Bazavov
and Weber (2021); Litim and Manuel (2002). In this section, we evaluate the
induced current $j^{\mu}_{a}(x)$ of Eq. (29) and find that in spacetimes of
interest, it is proportional to the square of an effective mass,
$m_{D}^{2}(x)$, which has spatial dependence.

For any $(3+1)$-dimensional, globally hyperbolic spacetime, we can choose
coordinates $x^{\mu}$ to put the metric in the Arnowitt-Deser-Misner (ADM)
form (see, e.g., Ref. Wald (1984))

$\begin{split}ds^{2}&=-N^{2}dt^{2}+\gamma_{ij}\left(dx^{i}+\beta^{i}dt\right)\left(dx^{j}+\beta^{j}dt\right)\\\
&=-\left(N^{2}-\beta_{i}\beta^{i}\right)dt^{2}+2\beta_{i}\,dtdx^{i}+\gamma_{ij}dx^{i}dx^{j},\end{split}$
(31)

where $x^{0}=t$ is a global time function, $N(x)$ is the lapse funtion, and
$\beta^{i}(x)$ the shift vector, whose indices are raised and lowered with
$\gamma_{ij}$: $\beta_{i}=\gamma_{ij}\,\beta^{j}$. Hypersurfaces of constant
time $t$ are Cauchy surfaces by construction, with induced metric
$\gamma_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$, where

$n_{\mu}=-(dt)_{\mu}=\left(-N,{\bf 0}\right).$ (32)

We normalize $t$ such that $N^{2}\to 1$ on the spatial (possibly asymptotic)
boundary. The components of the inverse metric are given by

$g^{00}=-\frac{1}{N^{2}},\>g^{0i}=\frac{\beta^{i}}{N^{2}},\>g^{ij}=\gamma^{ij}-\frac{\beta^{i}\beta^{j}}{N^{2}}.$
(33)

In these coordinates, the kinetic momentum takes the form

$P_{\mu}=\left(-Nk+\beta^{i}P_{i},P_{i}\right),\>\>P^{\mu}=\left(\frac{k}{N},\gamma^{ij}P_{j}-\frac{k}{N}\beta^{i}\right),$
(34)

where

$k\equiv\sqrt{\gamma^{\mu\nu}\,P_{\mu}P_{\nu}}=\sqrt{\gamma^{ij}\,P_{i}P_{j}}$
(35)

is the magnitude of the momentum projected on constant-$t$ hypersurfaces.

The spacetime metric $g_{\mu\nu}$ induces a metric in all of phase space,
known as the Sasaki metric Sasaki (1958); Israel (1989); Acuña Cárdenas  _et
al._ (2022), such that the invariant momentum volume form is

$d\omega=\frac{d^{4}p}{(2\pi)^{3}\sqrt{-g}}.$ (36)

We restrict to the mass-shell by integrating over
$2\delta(P^{2})\Theta(P^{0})$. To lowest order in $g_{s}$, the kinetic and
canonical momenta coincide, $P_{\mu}\rightarrow p_{\mu}$, so

$\frac{d^{4}p}{(2\pi)^{3}\sqrt{-g}}\rightarrow\frac{d^{3}p_{(i)}}{(2\pi)^{3}p^{0}\,N\sqrt{\gamma}},$
(37)

where the last step follows upon noting that $\sqrt{-g}=N\sqrt{\gamma}$. The
current $j^{\mu}_{a}(x)$ in Eq. (29) may then be written

$\begin{split}j^{\mu}_{a}(x)&=\frac{g_{s}^{2}}{(2\pi)^{3}}\frac{A_{0}^{a}}{N}\\\
&\times\int\frac{d^{3}p_{(i)}}{\sqrt{\gamma}}{\cal
N}(x,p)\left[\delta^{\mu}_{0}+\delta^{\mu}_{i}\left(\beta^{i}-\frac{N}{k}\gamma^{ij}p_{j}\right)\right].\end{split}$
(38)

We next consider a change in momentum coordinates
$\tilde{k}^{i}\equiv(\gamma^{1/2})^{ij}p_{j}$, where $(\gamma^{1/2})^{ij}$ is
the square matrix that satisfies
$(\gamma^{1/2})_{mi}\gamma^{ij}(\gamma^{1/2})_{j\ell}=\delta_{m\ell}$. Then
$d^{3}p_{(i)}=\sqrt{\gamma}\,d^{3}\tilde{k}^{(i)}$, and
$k^{2}=\tilde{k}^{2}=\delta_{ij}\tilde{k}^{i}\tilde{k}^{j}$, which is even in
the components $\tilde{k}^{i}$. One could perform an additional coordinate
transformation as in Ref. Brandt  _et al._ (1995),
$p_{i}^{\prime}=p_{i}-\beta_{i}p_{0}/g_{00}$, to aid in evaluating the
integral in Eq. (38). For the applications of interest to us, we will instead
restrict attention to spacetimes in which the shift vector vanishes,
$\beta^{i}=0$. Then the integration in Eq. (38) may be performed exactly. In
such cases, we find

$\displaystyle j^{0}_{a}(x)$
$\displaystyle=\frac{g_{s}^{2}}{(2\pi)^{3}}\frac{A_{0}^{a}}{N}\int
d^{3}\tilde{k}^{(i)}\,{\cal N}(x,\tilde{k}),$ (39) $\displaystyle
j^{i}_{a}(x)$
$\displaystyle=-\frac{g_{s}^{2}}{(2\pi)^{3}}A_{0}^{a}(\gamma^{1/2})^{ij}\delta_{j\ell}\int
d^{3}\tilde{k}^{(i)}\,{\cal N}(x,\tilde{k})\frac{\tilde{k}^{\ell}}{k}.$ (40)

The quantity ${\cal N}(x,\tilde{k})$ is even in the components
$\tilde{k}^{i}$, so that $j^{i}_{a}(x)$ vanishes identically. We further note
that the momentum coordinates $\tilde{k}^{i}$ are Cartesian, so $j^{0}_{a}(x)$
is equivalent to

$\begin{split}j^{0}_{a}(x)&=\frac{g_{s}^{2}}{2\pi^{2}}\frac{A_{0}^{a}}{N}\int_{0}^{\infty}d\tilde{k}\,\tilde{k}^{2}\,{\cal
N}(x,\tilde{k})\\\ &=-\frac{m_{D}^{2}(x)}{N^{2}(x)}A_{0}^{a}(x),\end{split}$
(41)

with

$m_{D}^{2}(x)\equiv\frac{1}{6}\left(2N_{c}+N_{f}\right)\frac{g_{s}^{2}}{(\beta_{T}N(x))^{2}}.$
(42)

Eq. (41) corresponds to

$j_{\mu}^{a}(x)=m_{D}^{2}(x)A_{0}^{a}\,\delta_{\mu}^{0}.$ (43)

The induced current in Eq. (43) is equivalent to inserting an effective mass
$m_{D}(x)$ for the $A_{0}^{a}$ soft gluon components within the effective
action of Eq. (2).

The expression for $m_{D}(x)$ in Eq. (42) reduces to the usual expression in
Minkowski spacetime upon setting $N(x)\rightarrow 1$ and identifying
$1/\beta_{T}=T$ with the temperature of the plasma Blaizot and Iancu (2002);
Kapusta and Gale (2006); Bazavov and Weber (2021); Weldon (1982); Rebhan
(1993, 1994); Braaten and Nieto (1996); Boyanovsky  _et al._ (2000); Blaizot
_et al._ (1995); Litim and Manuel (2002); Mrowczynski  _et al._ (2017). In our
case, the Boltzmann factor for the high-momentum particles is
$\exp[-\beta_{T}\bar{H}^{(0)}]=\exp[\beta_{T}N(x)k(x)]$, since we are
considering spacetimes in which $\beta^{i}=0$. This is exactly what one would
expect for particles in a fluid that is undergoing normal flow, with four-
velocity $n^{\mu}(x)$ and therefore local energy per particle
$E\equiv-n^{\mu}p_{\mu}=k$, if we identify $T(x)\equiv 1/(\beta_{T}N(x))$ with
the local temperature. We may therefore identify the constant $\beta_{T}\equiv
1/T_{0}$ and write

$T(x)=\frac{T_{0}}{N(x)}.$ (44)

Given that we have normalized $t$ such that $N(x)\rightarrow 1$ on the spatial
(possibly asymptotic) boundary, $T_{0}$ is the temperature associated with
time $t$ on the spatial boundary. We have thus recovered the familiar Tolman
temperature gradient in a curved spacetime Tolman and Ehrenfest (1930); Tolman
(1930); Santiago and Visser (2018). The Debye mass then takes the form

$m_{D}^{2}(x)=\frac{1}{6}\left(2N_{c}+N_{f}\right)g_{s}^{2}\,T^{2}(x),$ (45)

with the local temperature $T(x)$ given by Eq. (44).

##  V Cosmological Applications

In this section we consider a few specific examples. We begin with the
familiar case of Debye screening of a non-Abelian plasma in Minkowski
spacetime, to identify several regimes of interest. Next we generalize to the
case of a spatially flat Friedmann-Lemâitre-Robertson-Walker (FLRW) spacetime,
which introduces a new scale (compared to the Minkowski case), set by the
Hubble radius. In the last subsection, we examine Debye screening in the
vicinity of a primordial black hole.

###  V.1 Debye screening in Minkowski spacetime

We first consider the behavior of soft gluon modes $A_{\mu}^{a}(x)$ within a
hot plasma in Minkowski spacetime, to clarify notation and identify several
physical regimes of interest. In that case, the lapse function becomes
trivial, $N(x)=1$, and the shift vector vanishes, $\beta^{i}=0$. In the static
limit, $\partial_{0}A_{\mu}^{a}(x)=0$, the exact equations of motion of Eq.
(5) reduce to

$\begin{split}\delta_{\nu}^{0}&\bigg{\\{}\delta^{ik}\partial_{k}F^{a}_{i0}+g_{s}f^{abc}\delta^{ik}A_{k}^{b}F_{i0}^{c}\bigg{\\}}\\\
&+\delta_{\nu}^{i}\bigg{\\{}\delta^{\ell k}\partial_{k}F_{\ell
i}^{a}+g_{s}f^{abc}\left(A_{0}^{b}F_{i0}^{c}+\delta^{\ell k}A_{k}^{b}F_{\ell
i}^{c}\right)\bigg{\\}}=j_{\nu}^{a},\end{split}$ (46)

where we have adopted Cartesian coordinates for the spatial sections. The
induced current of Eq. (43) reduces to

$j_{\nu}^{a}({\bf x})=m_{D}^{2}A_{0}^{a}({\bf x})\delta_{\nu}^{\>\>0}+{\cal
O}(g_{s}^{3}).$ (47)

A well-known exact solution to Eq. (46) is a superposition of a screened
component $A_{0}^{a}$ with a Wu-Yang monopole solution for the components
$A_{i}^{a}$. For the case of ${\rm SU}(2)$, the solution takes the form Wu and
Yang (1976a, b); Jackiw  _et al._ (1994); Blaizot and Iancu (2002)

$A_{0}^{a}({\bf x})=\frac{{\cal A}^{a}\,e^{-m_{D}r}}{r},\>\>A_{i}^{a}({\bf
x})=\frac{\epsilon_{aij}\hat{x}^{j}}{g_{s}r},$ (48)

where ${\cal A}^{a}={\cal A}\hat{x}^{a}$, $\hat{x}^{j}\equiv x^{j}/r$ is a
unit vector, and $\epsilon_{ijk}$ is the usual Levi-Civita symbol. Note that
the solution mixes spatial indices and color-space indices, which is
straightforward for ${\rm SU}(2)$, since both $i,j=1,2,3$ and $a,b=1,2,3$. (To
confirm that $A_{\mu}^{a}({\bf x})$ in Eq. (48) satisfies Eq. (46), it is
helpful to construct the projection operator
$P^{a}_{\>\>b}=\delta^{a}_{\>\>b}-\hat{x}^{a}\hat{x}_{b}$, noting that
$\hat{x}_{a}P^{a}_{\>\>b}=0=\hat{x}^{b}P^{a}_{\>\>b}$ and
$\partial_{j}\hat{x}^{i}=r^{-1}P^{i}_{\>\>j}$.) For SU(3), exact monopole
solutions have been found for subgroups such as ${\rm U}(1)\times{\rm U}(1)$
and U(2), for which the components $A_{i}^{a}$ have comparable asymptotic
behavior to the solution in Eq. (48) Sinha (1976); Preskill (1984); Kunz and
Masak (1987).

The solution for $A_{\mu}^{a}({\bf x)}$ in Eq. (48) carries both
chromoelectric charges ${\cal Q}^{a}$ and chromomagnetic charges ${\cal
P}^{a}$. We consider the quasi-local charges defined as

$\begin{split}{\cal Q}^{a}&\equiv\frac{1}{4\pi}\int_{\partial V}\left(\star
F^{a}\right)_{\mu\nu}dx^{\mu}\wedge dx^{\nu},\\\ {\cal
P}^{a}&\equiv\frac{1}{4\pi}\int_{\partial V}F^{a}_{\mu\nu}\,dx^{\mu}\wedge
dx^{\nu},\end{split}$ (49)

where the integrals are taken over the boundary of a subset $V\subseteq\Sigma$
of a Cauchy surface $\Sigma$ of constant time and, as usual, $(\star
F^{a})_{\mu\nu}\equiv\sqrt{-g}\,\epsilon_{\mu\nu\alpha\beta}F^{\alpha\beta}_{a}$
with $\epsilon_{\mu\nu\alpha\beta}$ the totally antisymmetric Levi-Civita
symbol. Although ${\cal Q}^{a}$ and ${\cal P}^{a}$ are each gauge-dependent
quantities, the bilinear combinations ${\cal Q}^{2}=\delta_{ab}{\cal
Q}^{a}{\cal Q}^{b}$ and ${\cal P}^{2}=\delta_{ab}{\cal P}^{a}{\cal P}^{b}$ are
each gauge-invariant. For the solution in Eq. (48), we find the charges at the
spatial origin to be

${\cal Q}^{a}|_{r\rightarrow 0}={\cal A}^{a},\>\>{\cal
P}^{a}=\frac{\hat{x}^{a}}{g_{s}}.$ (50)

Because of screening, the value of ${\cal Q}^{2}$ measured by an observer at
finite distance from the origin will be reduced compared to ${\cal
A}^{2}=\delta_{ab}{\cal A}^{a}\,{\cal A}^{b}$. We will write ${\cal
Q}_{0}=\sqrt{{\cal A}^{2}}$ to indicate the magnitude of the charge located at
the origin. The chromomagnetic charge, on the other hand, is not screened, and
hence ${\cal P}^{2}=1/g_{s}^{2}$ everywhere.

Whereas the magnitude of the chromomagnetic charge ${\cal P}$ is fixed by
$1/g_{s}$, the chromoelectric charge ${\cal Q}_{0}$ at the origin is a free
parameter and can be arbitrarily larger than $1/g_{s}$. Such a scenario would
be compatible with a collection of test charges at the origin such that the
energy associated with the charges does not backreact on the spacetime itself.
In that case, there will exist a self-consistent regime, $0\leq r\leq r_{\rm
abel}$, within which $|A_{0}^{a}({\bf x})|\gg|A_{i}^{a}({\bf x})|$, with
$r_{\rm abel}$ given by

$r_{\rm abel}\equiv\frac{1}{m_{D}}{\rm ln}\left[g_{s}{\cal Q}_{0}\right].$
(51)

For $0\leq r\leq r_{\rm abel}$, the system is quasi-Abelian, with
$A_{\mu}^{a}({\bf x})\simeq A_{0}^{a}({\bf x})\delta_{\mu}^{\>\>0}$, and hence
non-Abelian contributions of the form $f^{abc}A_{\mu}^{a}A_{\nu}^{b}$ remain
subdominant. Within this regime, the solution for $A_{0}^{a}({\bf x})$ is of
the form in Eq. (48), with ${\cal A}^{a}$ arbitrary and constant,
$\partial_{i}{\cal A}^{a}=0$.

###  V.2 Debye screening in an FLRW spacetime

We next consider Debye screening within a hot quark-gluon plasma in an
expanding FLRW spacetime at early times, prior to the QCD confinement
transition ($t<t_{\rm QCD}\sim 10^{-5}\,{\rm s}$). The line-element for a
spatially flat FLRW universe may be written

$ds^{2}=-dt^{2}+a^{2}(t)\,\left[dr^{2}+r^{2}d\Omega_{(2)}^{2}\right],$ (52)

where $a(t)$ is the scale factor and $r$ is a comoving radial coordinate. The
Hubble parameter is given by $H(t)\equiv\dot{a}/a$, where overdots denote
derivatives with respect to cosmic time $t$, and the Hubble radius is
$r_{H}=1/H$. The presence of $r_{H}$ introduces a new scale compared to
Minkowski spacetime, so now we must consider $H$ as well as $k_{\rm hard}$ and
$k_{\rm soft}$. Note also that in these coordinates the lapse function is
$N(x)=1$, so the temperature has no spatial gradients.

The term ${\cal L}_{\rm fluid}$ in the effective action of Eq. (2) includes
contributions from the high-momentum quarks and gluons in the QGP, which, as
noted in Section III, behave to leading order as a gas of non-interacting
particles with a radiation-dominated equation of state. We assume that the
(coarse-grained) energy density $\rho$ and pressure $P$ associated with the
high-momentum particles dominate $T_{\mu\nu}$, so that the Friedmann equation
takes the form

$H^{2}=\frac{1}{3M_{\rm pl}^{2}}\left(\frac{\pi^{2}}{30}g_{*}T^{4}\right)$
(53)

corresponding to a fluid of $g_{*}$ effectively massless degrees of freedom in
equilibrium at temperature $T$; for the Standard Model at temperatures much
greater than the top-quark mass ($m_{t}=173\,{\rm GeV}$), $g_{*}=106.75$
Mukhanov (2005); Boyanovsky  _et al._ (2006). During the radiation-dominated
phase, $a(t)\propto t^{1/2}$, so we have $H(t)=1/(2t)$ and

$T(t)=\left(\frac{90}{\pi^{2}g_{*}}\right)^{1/4}\left(\frac{M_{\rm
pl}}{2t}\right)^{1/2},$ (54)

consistent with adiabatic expansion, $T(t)\,a(t)={\rm constant}$.

The equations of motion from Eq. (5) are

$\begin{split}j_{\nu}^{a}&=\frac{1}{a^{2}}\,\delta_{\nu}^{0}\bigg{\\{}\partial_{i}F^{a}_{i0}+g_{s}f^{abc}A_{i}^{b}F_{i0}^{c}\bigg{\\}}\\\
&+\delta_{\nu}^{i}\bigg{\\{}HF_{i0}^{a}+\partial_{0}F_{i0}^{a}+\frac{1}{a^{2}}\partial_{j}F_{ji}^{a}\\\
&\qquad\qquad+g_{s}f^{abc}\left(A_{0}^{b}F_{i0}^{c}+\frac{1}{a^{2}}A_{j}^{b}F_{ji}^{c}\right)\bigg{\\}},\end{split}$
(55)

where repeated indices are summed. The Hubble parameter $H$ sets a time-scale
over which we expect cosmological dynamics to be relevant. Notice that the
terms that make Eq. (55) differ from Eq. (46) (up to factors of $a$) are
either proportional to $H$ or include a time derivative $\dot{A}_{\mu}^{a}$.
We consider an ansatz for $A_{\mu}^{a}(x)$ in which the only time dependence
arises from the scale factor $a(t)$. By the chain rule, any term with a time
derivative $\dot{A}_{\mu}^{a}$ will then also be proportional to
$H=\dot{a}/a$. If $H\ll k_{\rm soft}^{2}/k_{\rm hard}\sim g_{s}^{2}T$, then
all such terms will be subleading compared to $j_{\nu}^{a}\sim
m_{D}^{2}A_{0}^{a}\sim k_{\rm soft}^{2}k_{\rm hard}$.

If the hierarchy $H\ll g_{s}^{2}T$ is satisfied at early times, it will be
even more easily satisfied at later times, since $H\sim T^{2}/M_{\rm pl}$ from
Eq. (54) and $T/M_{\rm pl}$ decreases as the universe cools down, while
$g_{s}$ gradually runs to larger values. Now we proceed to show that $H\ll
g_{s}^{2}T$ holds at early times of interest.

Consider, as an example, dynamics at the time $t_{c}\sim 10^{-21}\,{\rm s}$,
the earliest time that a population of primordial black holes (PBHs) could
have formed, if the PBHs are to constitute all of dark matter today Carr and
Kühnel (2020); Green and Kavanagh (2021); Villanueva-Domingo  _et al._ (2021);
Escrivà  _et al._ (2022). At such early times, the Hubble scale
$H(t_{c})=3.0\times 10^{-4}\,{\rm GeV}$ and the fluid filling the FLRW
spacetime had a temperature $T(t_{c})=1.5\times 10^{7}\,{\rm GeV}$. At such
high energy scales, we must take into account the running of the strong
coupling $\alpha_{s}\equiv g_{s}^{2}/(4\pi)$. To two-loop order, the running
with energy scale $\mu$ is given by Bazavov and Weber (2021)

$\derivative{\alpha_{s}(\mu)}{\,{\rm
ln}\mu}=-2\alpha_{s}\left[b_{0}\frac{\alpha_{s}}{4\pi}+b_{1}\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\right],$
(56)

with

$\begin{split}b_{0}&=\frac{11}{3}N_{c}-\frac{2}{3}N_{f},\\\
b_{1}&=\frac{32}{3}N_{c}^{2}-\frac{10}{3}N_{c}N_{f}-\left(\frac{N_{c}^{2}-1}{N_{c}}\right)N_{f}.\end{split}$
(57)

Upon normalizing $\alpha_{s}(m_{Z})=0.118$ at the $Z$ boson mass
($m_{Z}=91.2\,{\rm GeV}$), we find $\alpha_{s}=0.046$ at the energy scale
$\mu=T(t_{c})=1.5\times 10^{7}$ GeV. This yields a ratio $g_{s}^{2}T/H\sim
10^{9}$. Hence we may self-consistently neglect $H\ll g_{s}^{2}T$ when
considering the dynamics of the gluon modes on length-scales $\lambda_{D}$,
even at very early times in cosmic history. In other words, the FLRW spacetime
is consistent with our quasi-stationary approximation from Eq. (24), since
$|\partial_{t}g^{\mu\nu}|\sim H\ll g_{s}k_{\rm soft}<k_{\rm soft}$.

Having shown that we may neglect cosmological dynamics, the scale factor
$a(t)$ will be approximately constant in regimes of interest. We may therefore
rescale the spatial coordinates as $x^{i}\rightarrow x^{i}/a$, and thus
recover identical equations of motion to those in the Minkowski background,
Eq. (46). A solution is the superposition of a screened component $A_{0}^{a}$
with a Wu-Yang monopole for the components $A_{i}^{a}$, which in the case of
${\rm SU}(2)$ takes a similar form to Eq. (48),

$A_{0}^{a}(x)=\frac{{\cal
A}^{a}\,e^{-m_{D}ar}}{ar},\>\>A_{i}^{a}(x)=\frac{\epsilon_{aij}\hat{x}^{j}}{g_{s}ar},$
(58)

where we have rescaled the spatial coordinates back so that $r$ is the
comoving radial coordinate from Eq. (52). The solution for $A_{\mu}^{a}(x)$ in
Eq. (58) remains consistent with the quasi-stationary approximation of Eq.
(24) for $\lambda_{D}\leq r\ll r_{H}$.

###  V.3 Debye screening near a primordial black hole

Our final example concerns Debye screening within the hot plasma surrounding a
primordial black hole (PBH) that formed early in cosmic history. PBHs form via
direct gravitational collapse of primordial overdensities, and their masses
are typically proportional to the Hubble mass $M_{H}$ at the time of formation
$t_{c}$, where $M_{H}(t_{c})=4\pi M_{\rm pl}^{2}/H(t_{c})$ is the mass
contained within a Hubble radius $H^{-1}(t_{c})$. Such black holes may
therefore form with a huge range of masses, depending on their time of
formation Carr and Kühnel (2020); Green and Kavanagh (2021); Villanueva-
Domingo  _et al._ (2021); Escrivà  _et al._ (2022).

Any PBHs that formed at times $t<t_{\rm QCD}=10^{-5}\,{\rm s}$ would undergo
collapse amid a hot plasma of unconfined quarks and gluons. On long length-
scales at such early times, $\lambda\gg 1/k_{\rm soft}$, we expect that the
plasma would have attained a charge-neutral equilibrium distribution. Yet on
shorter length-scales, set by $\lambda_{D}=1/m_{D}\sim 1/k_{\rm soft}\sim
1/(g_{s}T)$, color-coherent regions would persist, within which the color-
space vectors $Q^{a}$ for nearly all (soft) gluons would align Blaizot and
Iancu (1994c, b). The PBHs would therefore form by absorbing one or more
color-coherent regions. Depending on the ratio of the PBH radius ($\sim
GM_{\rm PBH}$) and the Debye length ($\sim 1/m_{D}$) at the time of formation,
PBHs could form with net color charge ${\cal Q}_{0}$. Such a scenario is
distinct from the examples reviewed in Ref. Volkov and Gal’tsov (1999); Volkov
(2017), which concern black hole solutions to the Einstein-Yang-Mills
equations in vacuum. In our case, the PBHs are immersed in a hot, active
medium.

We consider a scenario in which a PBH forms with a small net charge ${\cal
Q}_{0}$; the case of larger ${\cal Q}_{0}$ is treated in Ref. Alonso-Monsalve
and Kaiser (2023). In particular, we consider the regime

${\cal Q}_{0}^{2}\ll GM_{\rm PBH}^{2}.$ (59)

Within this regime, backreaction on spacetime from the energy associated with
the enclosed charge ${\cal Q}_{0}$ remains subdominant, and we may consider
the dynamics of soft gluon modes $A_{\mu}^{a}(x)$ within a fixed background
geometry. (For ${\cal Q}_{0}^{2}\sim GM_{\rm PBH}^{2}$, one must solve the
Einstein field equations as well as the equations of motion for
$A_{\mu}^{a}(x)$, which we consider in separate work Alonso Monsalve and
Kaiser .)

We consider a spherically symmetric spacetime with a black hole of mass
$M_{\rm PBH}$ at the origin surrounded by hot plasma with a radiation-
dominated equation of state. Our formalism holds in a quasi-static limit,
which requires that $M_{\rm PBH}$ not change appreciably over time-scales set
by $1/k_{\rm soft}$. One might suppose that if the fluid were moving with
respect to the black hole, $M_{\rm PBH}$ would increase by absorbing some of
the nearby fluid. Nonetheless, even if the relative speed between the fluid
and the black hole approached the speed of light, and the black hole absorbed
all the fluid in its path, then $M_{\rm PBH}$ would only increase by a
fraction $\delta M/M_{\rm PBH}\sim{\cal O}(1)\,{\cal R}\,(T/M_{\rm pl})$ over
an entire Hubble time, where ${\cal R}$ is the ratio of the black hole radius
to the Debye screening length. As in the previous section, we are interested
in early times at which the plasma temperature could have been as high as
$T\sim 10^{7}$ GeV, in which case $T/M_{\rm pl}\sim 10^{-11}$. We consider
black holes for which the ratio ${\cal R}$ is not exponentially larger than
one, so that the PBH forms with some small, residual color charge, and
therefore $\delta M/M$ remains negligible. A more formal calculation of the
Bondi accretion rate for PBHs in this scenario yields the same conclusion
Shapiro and Teukolsky (1983); Rice and Zhang (2017); De Luca  _et al._ (2020);
Richards  _et al._ (2021).

Within the regime indicated in Eq. (59), the spacetime can be described by the
McVittie line-element, which reduces to an ordinary Schwarzschild spacetime
near the origin and asymptotes to a spatially flat FLRW spacetime at large
distances Kaloper  _et al._ (2010); Faraoni  _et al._ (2014); Gaur and Visser
(2022). A convenient parameterization may be written Gaur and Visser (2022)

$\begin{split}ds^{2}&=-f(t,r)\,dt^{2}\\\
&\quad+\left[\frac{a(t)\left\\{dr+H(t)\,rdt\right\\}}{\sqrt{f(t,r)}}-H(t)\,a(t)\,rdt\right]^{2}\\\
&\quad+a^{2}(t)\,r^{2}d\Omega_{(2)}^{2},\end{split}$ (60)

with

$f(t,r)\equiv 1-\frac{r_{s}}{a(t)r}$ (61)

and $r_{s}\equiv 2GM_{\rm PBH}$ the usual Schwarzschild radius.

We found in Section V.2 that in cosmological scenarios of interest, $H\ll
m_{D}$, even at very early times. For the remainder of this section, we will
therefore set $H\sim 0$, for which $a(t)\sim{\rm constant}$ (which we will
scale to 1). In that limit, the lapse function reduces to
$N(x)\rightarrow\sqrt{f(r)}$, with $f(r)=1-r_{s}/r$, and the shift vector
vanishes, $\beta^{i}\rightarrow 0$.

Next we consider an appropriate range for the enclosed charge ${\cal Q}_{0}$,
subject to the constraint in Eq. (59). The shortest scales that can be
resolved within our EFT are set by $\lambda_{D}=1/m_{D}$, so we restrict
$r_{s}\geq\lambda_{D}$, which in turn requires $M_{\rm PBH}\geq 1/(2Gm_{D})$.
Then Eq. (59) corresponds to the upper bound

${\cal Q}_{0}\ll\frac{1}{\sqrt{2\alpha_{s}}}\,\frac{M_{\rm pl}}{T}\sim
10^{12}$ (62)

around $t_{c}\sim 10^{-21}\,{\rm s}$. Meanwhile, as discussed around Eq. (51),
if ${\cal Q}_{0}\gg 1/g_{s}$, then the soft modes $A_{\mu}^{a}$ will assume a
quasi-Abelian form, with $|A_{0}^{a}|\gg|A_{i}^{a}|$, for $r<r_{\rm abel}$.
The estimate of $r_{\rm abel}$ in Eq. (51) holds in homogeneous spacetimes,
for which the temperature has no spatial gradients. In the present case, given
Eq. (44), $m_{D}(x)\sim g_{s}T_{0}/\sqrt{f(r)}$. Following the same steps that
led to Eq. (51), we find

$\frac{r_{\rm
abel}}{r_{s}}=\frac{1}{2}\left[1+\sqrt{1+\left(\frac{2\tilde{\lambda}_{D}}{r_{s}}{\rm
ln}[g_{s}{\cal Q}_{0}]\right)^{2}}\right],$ (63)

where $\tilde{\lambda}_{D}$ is the Debye length associated with the
temperature $T_{0}$. For $r_{s}/\tilde{\lambda}_{D}=3.5$ and ${\cal Q}_{0}\sim
10^{4}$, this yields $r_{\rm abel}\sim 3\,r_{s}$, while ${\cal Q}_{0}\sim
10^{10}$ yields $r_{\rm abel}\sim 7\,r_{s}$.

The last departure to consider from the previous examples concerns the effect
of the local temperature $T(x)$ on the coupling $\alpha_{s}$. Since the local
temperature $T(r)$ increases as $r$ approaches $r_{s}$, the local QCD strength
runs toward zero as one approaches the black hole event horizon: an example of
asymptotic freedom within an inhomogeneous spacetime. To quantify the effect,
we use Eq. (56) for the running of $\alpha_{s}$ with energy scale
$\mu\rightarrow T(r)$. To resolve the behavior of $\alpha_{s}$ near $r_{s}$,
we adopt the “tortoise” radial coordinate, $r_{*}\equiv r+r_{s}{\rm
ln}[(r/r_{s})-1]$, for which $r=r_{s}$ corresponds to
$r_{*}\rightarrow-\infty$ Wald (1984). See Fig. 1.

Figure 1: The QCD coupling strength $\alpha_{s}\equiv g_{s}^{2}/(4\pi)$ runs
to lower values near the event horizon of a black hole, as the temperature of
the surrounding plasma increases. Here $r_{*}$ is the “tortoise” radial
coordinate, $r_{*}\equiv r+r_{s}{\rm ln}[(r/r_{s})-1]$; $r_{*}/r_{s}=-36$
corresponds to $r/r_{s}=1+10^{-16}$. We have set $T_{0}=10^{7}\,{\rm GeV}$ on
the asymptotic boundary.

In the vicinity of a primordial black hole, the evolution of the soft modes
$A_{\mu}^{a}(x)$ is therefore characterized by a hierarchy of length-scales:

$\lambda_{D}<r_{s}<r_{\rm abel}\ll r_{H}.$ (64)

For $r_{s}\leq r\leq r_{\rm abel}$, the soft modes will evolve as a quasi-
static, quasi-Abelian system within a fixed background spacetime. In that
regime, $A_{\mu}^{a}(x)\simeq A_{0}^{a}(r)\delta_{\mu}^{0}+{\cal O}(r/r_{\rm
abel})$, and the equations of motion in Eq. (5) reduce to

$\partial_{r}^{2}A_{0}^{a}(r)+\frac{2}{r}\partial_{r}A_{0}^{a}(r)-\frac{\tilde{m}_{D}^{2}}{f^{2}(r)}A_{0}^{a}(r)=0,$
(65)

where $\tilde{m}_{D}=1/\tilde{\lambda}_{D}$ is the Debye mass associated with
the temperature $T_{0}$. Including the running of $g_{s}$ with $r$, we may
solve Eq. (65) numerically, with a typical example shown in Fig. 2. The
component $A_{0}^{a}(r)$ undergoes strong screening for $r\gtrsim r_{s}$,
since more plasma gathers near the event horizon, yielding a higher density
and temperature than at locations $r\gg r_{s}$. Far from the black hole,
$A_{0}^{a}(r)$ asymptotes to similar behavior as found in Sections V.1 and
V.2, with $A_{0}^{a}\sim\exp[-m_{D}r]/r$. Note that an observer at $r\gtrsim
r_{s}$ would measure an effective charge ${\cal Q}(r)$ much smaller than the
charge ${\cal Q}_{0}$ contained within the black hole.

Figure 2: The soft mode $A_{0}^{a}(r)$ (solid line) for
$r_{s}/\tilde{\lambda}_{D}=3.5$ and $T_{0}=10^{7}\,{\rm GeV}$. Given the
greater density of plasma near the black hole, Debye screening is particularly
effective for $r\gtrsim r_{s}$. For $r\gg r_{s}$, the behavior of
$A_{0}^{a}(r)$ asymptotes to the expected slope, $\exp[-m_{D}r]/r$ (dashed
line).

Within our coordinate system, the Tolman temperature gradient drives
$T(r)\rightarrow\infty$ as $r\rightarrow r_{s}$. The gradient is physical, but
the divergence is an artifact of our fixed-background approximation; we have
not allowed the spacetime to backreact. To produce Fig. 2, we followed the
example of Refs. Thorne and MacDonald (1982); Thorne  _et al._ (1986);
Macdonald and Suen (1985) and evaluated the field with a boundary condition at
a “stretched horizon,” $r_{s}+\epsilon$, rather than at $r_{s}$, with
$\epsilon/r_{s}=10^{-6}$. In forthcoming work Alonso Monsalve and Kaiser , we
address this limitation by solving the coupled Einstein field equations and
equation of motion for $A_{0}^{a}(x)$ within the quasi-Abelian regime.

##  VI Discussion

We have generalized the description of Debye screening of charges in hot
plasmas to curved spacetime backgrounds. For fluids undergoing normal flow in
approximately static spacetimes, we found that the characteristic screening
length $\lambda_{D}(x)=1/m_{D}(x)$ is set by a Debye mass $m_{D}(x)\sim
g_{s}T(x)$ given in Eq. (45), where $g_{s}$ is the dimensionless gauge
coupling. This Debye mass is the natural generalization of the Minkowski
result, upon allowing the temperature to have spatial gradients due to
gravitational redshift, thereby reproducing Tolman’s classic result Tolman and
Ehrenfest (1930); Tolman (1930); Santiago and Visser (2018).

To characterize Debye screening, we constructed an effective theory for long-
wavelength excitations in a hot Yang-Mills plasma. We analyzed the dynamics of
high-momentum particles in the plasma (with momentum $k_{\rm hard}\sim T$) by
considering classical transport equations, and exploited the structure of
Hamiltonian mechanics to show that non-Abelian self-interactions induce an
effective local mass $m_{D}(x)$ for the $A_{0}^{a}$ components of the soft
modes (with momentum $k_{\rm soft}\sim g_{s}T$).

We applied our results to solve for the gauge potential $A_{\mu}^{a}$ in a few
cases of interest. We recovered the well-known Wu-Yang monopole solution in
the Minkowski limit and generalized it to FLRW spacetimes, demonstrating the
self-consistency of the static approximation in regimes of interest, for which
cosmological time-scales may be neglected compared to $1/k_{\rm soft}$.

Lastly, we analyzed Debye screening in the vicinity of a primordial black hole
immersed in a hot quark-gluon plasma. Such black holes form by absorbing
color-coherent regions of the plasma, with characteristic size set by
$\lambda_{D}$; hence the resulting primordial black holes can have a residual,
net color charge. Building on the examples involving homogeneous spacetimes,
we identified a regime in which the soft modes of the gauge field outside the
black hole exhibit quasi-Abelian behavior, $A_{\mu}^{a}(x)\approx
A_{0}^{a}(r)\delta_{\mu}^{0}$. Given the Tolman temperature gradients, the
interaction strength $\alpha_{s}$ runs to smaller values near the event
horizon. Incorporating this unusual example of asymptotic freedom, we solved
numerically for the soft-mode gauge potential $A_{0}^{a}(r)$, and found
enhanced screening of the charge enclosed within the black hole, due to an
increased density of the plasma near the event horizon.

Our examples have been restricted so far to fixed background spacetimes.
Future work will focus on exploiting our EFT to study realistic cosmological
scenarios involving primordial black holes. This will require solving the
coupled Einstein-Yang-Mills equations for a black hole in a hot, non-Abelian
plasma Alonso Monsalve and Kaiser . This formalism can also be used to
consider primordial black holes with substantial QCD color charge, which could
have formed at very early times in cosmic history, well before the QCD
confinement transition Alonso-Monsalve and Kaiser (2023).

## Acknowledgements

It is a pleasure to thank Chris Akers, Thomas Baumgarte, Jolyon Bloomfield,
Daniel Harlow, Scott Hughes, Mikhail Ivanov, Patrick Jefferson, Jérôme Martin,
Govert Nijs, Krishna Rajagopal, Phiala Shanahan, Vincent Vennin, and Xiaojun
Yao for helpful discussions. Portions of this work were conducted in MIT’s
Center for Theoretical Physics and supported in part by the U. S. Department
of Energy under Contract No. DE-SC0012567. EAM is also supported by a
fellowship from the MIT Department of Physics.

## References

  * Blaizot and Iancu (2002) Jean-Paul Blaizot and Edmond Iancu, “The quark gluon plasma: Collective dynamics and hard thermal loops,” Phys. Rept. 359, 355–528 (2002), arXiv:hep-ph/0101103 .
  * Litim and Manuel (2002) Daniel F. Litim and Cristina Manuel, “Semiclassical transport theory for nonAbelian plasmas,” Phys. Rept. 364, 451–539 (2002), arXiv:hep-ph/0110104 .
  * Kapusta and Gale (2006) J. I. Kapusta and Charles Gale,  _Finite-Temperature Field Theory: Principles and Applications_ (Cambridge University Press, Cambridge, UK, 2006).
  * Bazavov and Weber (2021) Alexei Bazavov and Johannes Heinrich Weber, “Color screening in quantum chromodynamics,” Prog. Part. Nucl. Phys. 116, 103823 (2021), arXiv:2010.01873 [hep-lat] .
  * Matinyan  _et al._ (1981) Sergei G. Matinyan, G. K. Savvidy, and N. G. Ter-Arutunian Savvidy, “Classical Yang-Mills mechanics: Nonlinear color oscillations,” Sov. Phys. JETP 53, 421–425 (1981).
  * Klimov (1981) V. V. Klimov, “Spectrum of Elementary Fermi Excitations in Quark Gluon Plasma. (In Russian),” Sov. J. Nucl. Phys. 33, 934–935 (1981).
  * Klimov (1982) V. V. Klimov, “Collective Excitations in a Hot Quark Gluon Plasma,” Sov. Phys. JETP 55, 199–204 (1982).
  * Weldon (1982) H. Arthur Weldon, “Covariant calculations at finite temperature: The relativistic plasma,” Phys. Rev. D 26, 1394 (1982).
  * Braaten and Pisarski (1990a) Eric Braaten and Robert D. Pisarski, “Resummation and Gauge Invariance of the Gluon Damping Rate in Hot QCD,” Phys. Rev. Lett. 64, 1338 (1990a).
  * Braaten and Pisarski (1990b) Eric Braaten and Robert D. Pisarski, “Soft Amplitudes in Hot Gauge Theories: A General Analysis,” Nucl. Phys. B 337, 569–634 (1990b).
  * Braaten and Pisarski (1992) Eric Braaten and Robert D. Pisarski, “Simple effective Lagrangian for hard thermal loops,” Phys. Rev. D 45, R1827 (1992).
  * Blaizot and Iancu (1994a) Jean Paul Blaizot and Edmond Iancu, “Soft collective excitations in hot gauge theories,” Nucl. Phys. B 417, 608–673 (1994a), arXiv:hep-ph/9306294 .
  * Blaizot and Iancu (1993) Jean Paul Blaizot and Edmond Iancu, “Kinetic equations for long wavelength excitations of the quark - gluon plasma,” Phys. Rev. Lett. 70, 3376–3379 (1993), arXiv:hep-ph/9301236 .
  * Blaizot and Iancu (1994b) Jean-Paul Blaizot and Edmond Iancu, “NonAbelian plane waves in the quark - gluon plasma,” Phys. Lett. B 326, 138–144 (1994b), arXiv:hep-ph/9401323 .
  * Blaizot  _et al._ (1995) Jean-Paul Blaizot, Edmond Iancu, and Rajesh R. Parwani, “On the screening of static electromagnetic fields in hot QED plasmas,” Phys. Rev. D 52, 2543–2562 (1995), arXiv:hep-ph/9504408 .
  * Mrowczynski  _et al._ (2017) Stanislaw Mrowczynski, Bjoern Schenke, and Michael Strickland, “Color instabilities in the quark–gluon plasma,” Phys. Rept. 682, 1–97 (2017), arXiv:1603.08946 [hep-ph] .
  * Blaizot and Iancu (1994c) Jean-Paul Blaizot and Edmond Iancu, “NonAbelian excitations of the quark - gluon plasma,” Phys. Rev. Lett. 72, 3317–3320 (1994c), arXiv:hep-ph/9401210 .
  * Amin  _et al._ (2014) Mustafa A. Amin, Mark P. Hertzberg, David I. Kaiser, and Johanna Karouby, “Nonperturbative dynamics of reheating after inflation: A review,” Int. J. Mod. Phys. D 24, 1530003 (2014), arXiv:1410.3808 [hep-ph] .
  * Allahverdi  _et al._ (2021) Rouzbeh Allahverdi _et al._ , “The first three seconds: A review of possible expansion histories of the early universe,” Open J. Astrophys. 4 (2021), 10.21105/astro.2006.16182, arXiv:2006.16182 [astro-ph.CO] .
  * Mukhanov (2005) V. Mukhanov,  _Physical Foundations of Cosmology_ (Cambridge University Press, New York, 2005).
  * Boyanovsky  _et al._ (2006) D. Boyanovsky, H. J. de Vega, and D. J. Schwarz, “Phase transitions in the early and the present universe,” Ann. Rev. Nucl. Part. Sci. 56, 441–500 (2006), arXiv:hep-ph/0602002 .
  * Carr and Kühnel (2020) Bernard Carr and Florian Kühnel, “Primordial Black Holes as Dark Matter: Recent Developments,” Ann. Rev. Nucl. Part. Sci. 70, 355–394 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Green and Kavanagh (2021) Anne M. Green and Bradley J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Villanueva-Domingo  _et al._ (2021) Pablo Villanueva-Domingo, Olga Mena, and Sergio Palomares-Ruiz, “A brief review on primordial black holes as dark matter,” Front. Astron. Space Sci. 8, 87 (2021), arXiv:2103.12087 [astro-ph.CO] .
  * Escrivà  _et al._ (2022) Albert Escrivà, Florian Kuhnel, and Yuichiro Tada, “Primordial Black Holes,” (2022), arXiv:2211.05767 [astro-ph.CO] .
  * Tolman (1930) Richard C. Tolman, “On the Weight of Heat and Thermal Equilibrium in General Relativity,” Phys. Rev. 35, 904–924 (1930).
  * Tolman and Ehrenfest (1930) Richard Tolman and Paul Ehrenfest, “Temperature Equilibrium in a Static Gravitational Field,” Phys. Rev. 36, 1791–1798 (1930).
  * Santiago and Visser (2018) Jessica Santiago and Matt Visser, “Tolman-like temperature gradients in stationary spacetimes,” Phys. Rev. D 98, 064001 (2018), arXiv:1807.02915 [gr-qc] .
  * Holcomb and Tajima (1989) K. A. Holcomb and T. Tajima, “General relativistic plasma physics in the early universe,” Phys. Rev. D 40, 3809–3818 (1989).
  * Dettmann  _et al._ (1993) C. P. Dettmann, N. E. Frankel, and V. Kowalenko, “Plasma electrodynamics in the expanding Universe,” Phys. Rev. D 48, 5655–5667 (1993).
  * Prokopec and Woodard (2004) Tomislav Prokopec and Richard P. Woodard, “Vacuum polarization and photon mass in inflation,” Am. J. Phys. 72, 60–72 (2004), arXiv:astro-ph/0303358 .
  * Prokopec and Puchwein (2004) Tomislav Prokopec and Ewald Puchwein, “Photon mass generation during inflation: de Sitter invariant case,” JCAP 04, 007 (2004), arXiv:astro-ph/0312274 .
  * Dodin and Fisch (2010) I. Y. Dodin and N. J. Fisch, “On the evolution of linear waves in cosmological plasmas,” Phys. Rev. D 82, 044044 (2010).
  * Popov (2018) Fedor K. Popov, “Debye mass in de Sitter space,” JHEP 06, 033 (2018), arXiv:1711.11010 [hep-th] .
  * Gyulassy and McLerran (2005) Miklos Gyulassy and Larry McLerran, “New forms of QCD matter discovered at RHIC,” Nucl. Phys. A 750, 30–63 (2005), arXiv:nucl-th/0405013 .
  * Scardina  _et al._ (2013) F. Scardina, M. Colonna, S. Plumari, and V. Greco, “Quark-to-gluon composition of the quark-gluon plasma in relativistic heavy-ion collisions,” Phys. Lett. B 724, 296–300 (2013), arXiv:1202.2262 [nucl-th] .
  * Borsanyi  _et al._ (2012) Sz. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, “Precision SU(3) lattice thermodynamics for a large temperature range,” JHEP 07, 056 (2012), arXiv:1204.6184 [hep-lat] .
  * Borsanyi  _et al._ (2014) Szabocls Borsanyi, Zoltan Fodor, Christian Hoelbling, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo, “Full result for the QCD equation of state with 2+1 flavors,” Phys. Lett. B 730, 99–104 (2014), arXiv:1309.5258 [hep-lat] .
  * Bazavov  _et al._ (2014) A. Bazavov _et al._ (HotQCD), “Equation of state in (2+1)-flavor QCD,” Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Pasechnik and Šumbera (2017) Roman Pasechnik and Michal Šumbera, “Phenomenological Review on Quark–Gluon Plasma: Concepts vs. Observations,” Universe 3, 7 (2017), arXiv:1611.01533 [hep-ph] .
  * Pisarski (1989) Robert D. Pisarski, “Scattering Amplitudes in Hot Gauge Theories,” Phys. Rev. Lett. 63, 1129 (1989).
  * Braaten and Pisarski (1990c) Eric Braaten and Robert D. Pisarski, “Deducing Hard Thermal Loops From Ward Identities,” Nucl. Phys. B 339, 310–324 (1990c).
  * Braaten and Pisarski (1990d) Eric Braaten and Robert D. Pisarski, “Calculation of the gluon damping rate in hot QCD,” Phys. Rev. D 42, 2156–2160 (1990d).
  * Kelly  _et al._ (1994a) P. F. Kelly, Q. Liu, C. Lucchesi, and C. Manuel, “Deriving the hard thermal loops of QCD from classical transport theory,” Phys. Rev. Lett. 72, 3461–3463 (1994a), arXiv:hep-ph/9403403 .
  * Kelly  _et al._ (1994b) P. F. Kelly, Q. Liu, C. Lucchesi, and C. Manuel, “Classical transport theory and hard thermal loops in the quark - gluon plasma,” Phys. Rev. D 50, 4209–4218 (1994b), arXiv:hep-ph/9406285 .
  * Pisarski (1998) Robert D. Pisarski, “NonAbelian Debye screening, tsunami waves, and worldline fermions,” NATO Sci. Ser. C 511, 195–210 (1998), arXiv:hep-ph/9710370 .
  * Brandt  _et al._ (1995) Fernando T. Brandt, J. Frenkel, and J. C. Taylor, “High temperature QCD and the classical Boltzmann equation in curved space-time,” Nucl. Phys. B 437, 433–446 (1995), arXiv:hep-th/9411130 .
  * Wong (1970) S. K. Wong, “Field and particle equations for the classical Yang-Mills field and particles with isotopic spin,” Nuovo Cim. A 65, 689–694 (1970).
  * Linden  _et al._ (1996) N. Linden, A. J. Macfarlane, and J. W. van Holten, “Particle motion in a Yang-Mills field: Wong’s equations and spin 1/2 analogs,” Czech. J. Phys. 46, 209–215 (1996), arXiv:hep-th/9512071 .
  * Balachandran  _et al._ (1977) A. P. Balachandran, Per Salomonson, Bo-Sture Skagerstam, and Jan-Olof Winnberg, “Classical Description of Particle Interacting with Nonabelian Gauge Field,” Phys. Rev. D 15, 2308–2317 (1977).
  * Barducci  _et al._ (1977) A. Barducci, R. Casalbuoni, and L. Lusanna, “Classical Scalar and Spinning Particles Interacting with External Yang-Mills Fields,” Nucl. Phys. B 124, 93–108 (1977).
  * Arnold (1978) V. I. Arnold,  _Mathematical Methods of Classical Mechanics_ (Springer, New York, 1978).
  * Sato (2021) Naoki Sato, “The effect of spacetime curvature on statistical distributions,” Class. Quant. Grav. 38, 165003 (2021), arXiv:2104.09019 [gr-qc] .
  * Wald (1984) Robert M. Wald,  _General Relativity_ (University of Chicago Press, Chicago, 1984).
  * Sasaki (1958) Shigeo Sasaki, “On the differential geometry of tangent bundles of Riemannian manifolds,” Tohoku Mathematical Journal 10, 338 – 354 (1958).
  * Israel (1989) Werner Israel, “Covariant fluid mechanics and thermodynamics: An introduction,” in  _Relativistic Fluid Dynamics_, edited by A. M. Anile and Y. Choquet-Bruhat (Springer, Berlin, 1989) pp. 152–210.
  * Acuña Cárdenas  _et al._ (2022) Rubén O. Acuña Cárdenas, Carlos Gabarrete, and Olivier Sarbach, “An introduction to the relativistic kinetic theory on curved spacetimes,” Gen. Rel. Grav. 54, 23 (2022), arXiv:2106.09235 [gr-qc] .
  * Rebhan (1993) A. K. Rebhan, “The NonAbelian Debye mass at next-to-leading order,” Phys. Rev. D 48, R3967–R3970 (1993), arXiv:hep-ph/9308232 .
  * Rebhan (1994) Anton K. Rebhan, “NonAbelian Debye screening in one loop resummed perturbation theory,” Nucl. Phys. B 430, 319–344 (1994), arXiv:hep-ph/9408262 .
  * Braaten and Nieto (1996) Eric Braaten and Agustin Nieto, “Free energy of QCD at high temperature,” Phys. Rev. D 53, 3421–3437 (1996), arXiv:hep-ph/9510408 .
  * Boyanovsky  _et al._ (2000) D. Boyanovsky, H. J. de Vega, and M. Simionato, “Nonequilibrium quantum plasmas in scalar QED: Photon production, magnetic and Debye masses and conductivity,” Phys. Rev. D 61, 085007 (2000), arXiv:hep-ph/9909259 .
  * Wu and Yang (1976a) Tai Tsun Wu and Chen Ning Yang, “Dirac’s Monopole Without Strings: Classical Lagrangian Theory,” Phys. Rev. D 14, 437–445 (1976a).
  * Wu and Yang (1976b) Tai Tsun Wu and Chen Ning Yang, “Dirac Monopole Without Strings: Monopole Harmonics,” Nucl. Phys. B 107, 365 (1976b).
  * Jackiw  _et al._ (1994) R. Jackiw, Q. Liu, and C. Lucchesi, “Hard thermal loops, static response and the composite effective action,” Phys. Rev. D 49, 6787–6793 (1994), arXiv:hep-th/9401002 .
  * Sinha (1976) A. Sinha, “SU(3) Magnetic Monopoles,” Phys. Rev. D 14, 2016–2022 (1976).
  * Preskill (1984) John Preskill, “Magnetic monopoles,” Ann. Rev. Nucl. Part. Sci. 34, 461–530 (1984).
  * Kunz and Masak (1987) J. Kunz and D. Masak, “Finite energy SU(3) monopoles,” Phys. Lett. B 196, 513–518 (1987).
  * Volkov and Gal’tsov (1999) Mikhail S. Volkov and Dmitri V. Gal’tsov, “Gravitating nonAbelian solitons and black holes with Yang-Mills fields,” Phys. Rept. 319, 1–83 (1999), arXiv:hep-th/9810070 .
  * Volkov (2017) Mikhail S. Volkov, “Hairy black holes in the XX-th and XXI-st centuries,” in  _14th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories_, Vol. 2 (2017) pp. 1779–1798, arXiv:1601.08230 [gr-qc] .
  * Alonso-Monsalve and Kaiser (2023) Elba Alonso-Monsalve and David I. Kaiser, “Primordial black holes with QCD color charge,” in prep. (2023).
  * (71) Elba Alonso-Monsalve and David I. Kaiser, in prep .
  * Shapiro and Teukolsky (1983) Stuart L. Shapiro and Saul A. Teukolsky,  _Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects_ (Wiley, New York, 1983).
  * Rice and Zhang (2017) Jared R. Rice and Bing Zhang, “Cosmological evolution of primordial black holes,” JHEAp 13-14, 22–31 (2017), arXiv:1702.08069 [astro-ph.HE] .
  * De Luca  _et al._ (2020) V. De Luca, G. Franciolini, P. Pani, and A. Riotto, “Constraints on Primordial Black Holes: the Importance of Accretion,” Phys. Rev. D 102, 043505 (2020), arXiv:2003.12589 [astro-ph.CO] .
  * Richards  _et al._ (2021) Chloe B. Richards, Thomas W. Baumgarte, and Stuart L. Shapiro, “Relativistic Bondi accretion for stiff equations of state,” Mon. Not. Roy. Astron. Soc. 502, 3003–3011 (2021), [Erratum: Mon.Not.Roy.Astron.Soc. 506, 3935 (2021)], arXiv:2101.08797 [astro-ph.HE] .
  * Kaloper  _et al._ (2010) Nemanja Kaloper, Matthew Kleban, and Damien Martin, “McVittie’s Legacy: Black Holes in an Expanding Universe,” Phys. Rev. D 81, 104044 (2010), arXiv:1003.4777 [hep-th] .
  * Faraoni  _et al._ (2014) Valerio Faraoni, Andres F. Zambrano Moreno, and Angus Prain, “The charged McVittie spacetime,” Phys. Rev. D 89, 103514 (2014), arXiv:1404.3929 [gr-qc] .
  * Gaur and Visser (2022) Rudeep Gaur and Matt Visser, “Cosmology in Painlevé-Gullstrand coordinates,” JCAP 09, 030 (2022), arXiv:2207.08375 [gr-qc] .
  * Thorne and MacDonald (1982) K. S. Thorne and D. MacDonald, “Electrodynamics in Curved Spacetime - 3+1 Formulation,” Mon. Not. Roy. Astron. Soc. 198, 339 (1982).
  * Thorne  _et al._ (1986) Kip S. Thorne, Richard H. Price, and Douglas A. MacDonald,  _Black holes: The membrane paradigm_ (Yale University Press, New Haven, CT, 1986).
  * Macdonald and Suen (1985) Douglas A. Macdonald and Wai-Mo Suen, “Membrane viewpoint on black holes: Dynamical electromagnetic fields near the horizon,” Phys. Rev. D 32, 848–871 (1985).