Debye Screening of Non-Abelian Plasmas in Curved Spacetimes

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 λ D (T ) = 1/m D (T ) set by the Debye mass m D (T ) ∼ g s T , where g s is the dimensionless gauge coupling strength and T is the temperature of the plasma.(For reviews, see Refs.[1][2][3][4].) 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 soft ∼ m D (T ) ≪ T , propagate with frequency set by the plasma frequency ω p ≃ m D (T ) [5][6][7][8][9][10][11][12][13][14][15][16][17].Moreover, nontrivial spatial distributions of color charge among the soft gluons can form, with typical length-scale 1/m D (T ) [18,19].These effects have been studied extensively in Minkowski spacetime, using both analytic techniques and lattice simulations [1][2][3][4].
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 ∼ O(10 14 GeV) ∼ O(10 27 K) by the time t therm ∼ O(10 −35 s) after the big bang [20,21].Such temperatures are exponentially greater than the QCD confinement scale Λ QCD ≃ 0.170 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 [22,23].
As the universe expanded, the temperature of the plasma fell adiabatically: T (t)a(t) ≃ constant, where a(t) is the scale factor of the Friedmann-Lemaître-Robertson-Walker (FLRW) line-element.During the radiationdominated phase, a(t) = (t/t therm ) 1/2 .The QCD confinement transition, after which quarks and gluons remained bound in color-neutral hadronic states, occurred at t QCD ≃ 10 −5 s, when T (t QCD ) ≃ Λ QCD [22,23].
Within the range of times t therm ≪ t ≪ t 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 g ≤ M (t c ) ≤ 10 22 g, they could account for the entire abundance of dark matter today [24][25][26][27].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 s ≤ t c ≤ 10 −16 s.At these early times, the temperature of the plasma within which the primordial black holes formed would have been 10 5 GeV ≤ T (t c ) ≤ 10 7 GeV, for which T (t c ) ≫ T (t 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 → T (x), akin to the familiar Tolman temperature [28][29][30].We identify corrections to the induced current j 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.[31][32][33][34][35][36][37][38][39] on the behavior of Abelian plasmas within curved spacetimes.) In Section II, we introduce an effective theory for soft modes within the plasma, generalizing the elegant formalism reviewed in Refs.[1,2] for application to curved arXiv:2309.15385v2[hep-ph] 19 Nov 2023 spacetimes.In Section III we evaluate the induced current j 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 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 µ, ν = 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 SU(N c ) range over a, b, c = 1, 2, ..., N 2 c − 1 and are raised and lowered with δ ab .The generators T a of the Lie algebra for 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 Tr(T a T b ) = δ ab .We adopt natural units in which c = ℏ = k B = 1, in terms of which the reduced Planck mass is given by M pl ≡ 1/ √ 8πG = 2.43 × 10 18 GeV.

II. EFFECTIVE THEORY FOR SOFT MODES
At early times t ≪ t QCD , when the temperature of the plasma filling the universe satisfied T ≫ Λ 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 ∼ 1/α s , where α s ≡ g 2 s /(4π) and g s is the QCD coupling constant [40].Given the running of α 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 [1,4].Similar behavior has been observed in the plasmas that form immediately after relativistic heavy-ion collisions [40,41].Hence we expect gluons to dominate the dynamics, and consider the effective action for a Yang-Mills model with gauge group SU(N c ) and field strength tensor where ∇ µ is the (spacetime) covariant derivative associated with metric g µν (x).Following the approach reviewed in Refs.[1,2], 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 Ãa µ with typical momenta k soft ∼ g s T .Thermal and quantum corrections to the behavior of the soft modes are dominated by quanta a a µ with momenta k hard ∼ T .We therefore write A a µ = Ãa µ + a a µ , integrate out the highmomentum modes a a µ , and drop the tilde on the soft modes Ãa µ → A a µ .As described in Refs.[1,2], the soft modes A 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 λ ≥ 1/k soft ∼ 1/(g s T ).We follow the background-field method of Ref. [1], whereby we choose to decompose the action of the gauge symmetry on a a µ and A a µ in such a way that leaves the soft modes A a µ unchanged.One example is to use a generalized Coulomb-type gauge-fixing term in the full theory, ∇ i a i a − g s f abc A b i a i c [1].No additional gauge-fixing terms or ghosts are then required in the effective action for the soft modes.
For long length-scales λ ≥ 1/k soft ≫ 1/k hard ∼ 1/T , the dynamics may be described by the effective action where j a µ (x) is the induced current generated by (non-Abelian) self-interactions with the high-momentum modes a a µ (x); the induced current j a µ also includes contributrions from high-momentum quarks.The term L 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 λ ≫ 1/k 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 µν (x) and hence impacts the dynamics of the soft modes A a µ (x).Varying S eff with respect to A a µ yields the equations of motion where the covariant derivative D µ acting on a spacetime tensor that transforms in the adjoint representation of SU(N c ), X a ν1•••νn , is defined as The left-hand side of Eq. ( 3) can be expanded as where R µν is the Ricci tensor for the background spacetime.The equations of motion imply that the induced current j a µ must be covariantly conserved with respect to both the spacetime curvature and the SU(N c ) gauge group: for each color a.
The contribution to the energy-momentum tensor from the Yang-Mills field is given, as usual, by while the additional fluid contributes These terms satisfy the usual covariant conservation equations in the presence of the induced current j a µ .

III. INDUCED CURRENT
The induced current j a µ (x) for the soft modes A 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 ≫ T 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 socalled "trace anomaly," (ρ − 3P )/T 4 , where ρ is the energy density and P the pressure of particles in the plasma, tends rapidly toward zero for T ≫ T QCD , and therefore ρ ≈ 3P .(See, e.g., Refs.[42][43][44][45].) At high temperatures, the soft quantum fields of the plasma are well approximated by classical fields, due to their large occupation numbers per mode.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 soft ∼ g s T ) and hard internal momenta (with k hard ∼ T ) [10,11,[46][47][48].These effects can be analyzed within a mean-field approximation of the (truncated) Schwinger-Dyson equations [1,12,13,[15][16][17], or by simply adopting classical transport equations for highmomentum particles [2,17,[49][50][51].In this section we adapt the latter approach to arbitrary curved spacetimes.(See also Ref. [52].) The equations of motion for a classical point particle with color charge in a soft gauge-field background are well known [50,[52][53][54]: Here x µ is the position of the particle, P µ = dx µ /dλ its kinetic 4-momentum, λ an affine parameter, and Q a its SU(N c ) charge.We take advantage of the freedom to linearly rescale the affine parameter λ so that dx µ /dλ 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 µ a (x), summed over all hard particles, is given by [1,2] where δ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ω 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).The phase-space structure for the color-charge degrees of freedom is unaffected by the spacetime structure, so we may use the now-standard parameterization developed for previous studies within Minkowski spacetime, as reviewed, e.g., in Ref. [2]; an explicit construction for the cases SU(2) and SU( 3) is available in Section III of Ref. [50].In particular, for the gauge group SU(2) the color-space measure is and for the gauge group SU(3) we have For the SU(3) case, the d abc are the totally symmetric constants given by d abc = 2 Tr({T a , T b }, T c ).The representation-dependent constant c R is fixed by the normalization 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 [2,49,50].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 The relevant dynamics for the induced current j µ a (x) in Eq. ( 11) are captured to leading order by identifying δf (x, p, Q) = g s f (1) (x, p, Q), as demonstrated by Refs.[49,50].The dynamical evolution of δ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.[54][55][56], 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 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 µ and q a as The corresponding Hamiltonian is simply The kinetic momentum P µ ≡ dx µ /dλ is related to the canonical momentum p µ , conjugate to x µ , by Then Hamilton's equations are Eq.( 10), as desired.Dynamical evolution in phase space is generated by the Liouville vector field X H : where, in the second line, we have written the vector field in noncanonical coordinates {x µ , P µ , Q a }, and we have noted explicitly which derivatives are taken keeping P µ fixed, as opposed to p µ .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 : and We assume that all dependence of f (x, p, Q) on the charges Q a can only appear at O(g s ) or above, and therefore ∂f (0) /∂Q a = 0. Notice that, as emphasized in Ref. [51], 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) ∼ exp[−βE], for a constant β and a quantity E that may be interpreted as an energy.The equation X H [f ] = 0 is typically solved to O(g s ) by employing Green's function techniques [1,2,51].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 O(g s ), because quantum corrections to the classical equations of motion in Eq. ( 10) only arise at O(g 2 s ) and above.It is a classic result that Hamiltonian mechanics in the fully covariant phase space {x µ , p µ , Q a } with evolution parameterized by an affine λ 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 ≡ x 0 and generated by the reduced Hamiltonian This means that both formalisms yield the same equations of motion.This is possible thanks to the redundancy in the choice of affine parameter λ, as well as the conservation of H(x µ , p µ , 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 µ , p µ , Q) ≡ h, with h a constant.In the covariant formalism, the conservation of H(x µ , p µ , Q) follows from Eq. ( 10) and the fact that H cannot depend explicitly on λ.The reduced phase-space formalism does not yield the conservation of H(x µ , p µ , Q a ), so we must impose it by hand.For details of the proof, see Sections 44 and 45 in Chapter 9 of Ref. [57].
If the reduced Hamiltonian H from Eq. ( 22) does not depend explicitly on t, then it is also conserved along phase-space trajectories, X H [ H] = 0, so we may exploit H to construct a distribution function.That is, we may use the fact that any scalar function of 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 ∂ t g µν and ∂ t A 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 hard ∼ 1/T , whereas we expect the soft modes A a µ (x) to evolve on time-scales set by 1/k soft ≫ 1/k hard .Likewise, our effective description can only resolve dynamics up to scales set by 1/k soft , which bounds how sharply the background spacetime can evolve within our self-consistent expansion as well.In particular, if we allowed H to have an explicit dependence on t, then and thus a distribution function constructed from H would yield self-consistent dynamics as long as Eq. ( 24) involves coordinate-dependent quantities.However, any change in x µ would be accompanied by changes in the conjugate momenta p µ , such that Eq. ( 24) remains meaningful, given that k soft is a momentum scale.This follows from the invariance of the phase-space volume under coordinate transformations.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 −β T H , where β T is a constant [58].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 Given f (0) , we may evaluate the equilibrium occupation numbers per mode in the usual way, where the sums run from n i = 0 to n i = ∞ for bosons and to n i = 1 for fermions.Thus we obtain the Bose-Einstein and Fermi-Dirac distributions at zeroth-order: Substituting into Eq.( 21) yields for bosons and fermions, respectively.Summing over species and helicities, we have δf = g S (2n F ) in Eq. (11).
To evaluate j µ a (x) we first perform the Q integral with the measure in Eq. ( 13), which yields factors proportional to the index of the representation: dQ Q a Q b = C B,F δ a b , with C B = N c for bosons and C F = 1/2 for fermions.Then the current becomes where we have defined 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 [1], j µ a (x) is proportional to A a 0 (x).

IV. DEBYE MASS
In Minkowski spacetime, the induced current reduces (in the static limit) to j a µ (x) = m 2 D A a 0 δ 0 µ , effectively giving a constant mass to the soft gluon components A a 0 , which is responsible for Debye screening [1][2][3][4].In this section, we evaluate the induced current j µ a (x) of Eq. ( 29) and find that in spacetimes of interest, it is proportional to the square of an effective mass, m 2 D (x), which has spatial dependence.
For any (3 + 1)-dimensional, globally hyperbolic spacetime, we can choose coordinates x µ to put the metric in the Arnowitt-Deser-Misner (ADM) form (see, e.g., Ref. [59]) where x 0 = t is a global time function, N (x) is the lapse funtion, and β i (x) the shift vector, whose indices are raised and lowered with γ ij : β i = γ ij β j .Hypersurfaces of constant time t are Cauchy surfaces by construction, with induced metric γ µν = g µν + n µ n ν , where We normalize t such that N 2 → 1 on the spatial (possibly asymptotic) boundary.The components of the inverse metric are given by In these coordinates, the kinetic momentum takes the form where k ≡ γ µν P µ P ν = γ ij P i P j (35) is the magnitude of the momentum projected on constant-t hypersurfaces.
The spacetime metric g µν induces a metric in all of phase space, known as the Sasaki metric [60][61][62], such that the invariant momentum volume form is We restrict to the mass-shell by integrating over 2δ(P 2 )Θ(P 0 ).To lowest order in g s , the kinetic and canonical momenta coincide, P µ → p µ , so where the last step follows upon noting that √ −g = N √ γ.The current j µ a (x) in Eq. ( 29) may then be written We next consider a change in momentum coordinates ki ≡ (γ 1/2 ) ij p j , where (γ 1/2 ) ij is the square matrix that satisfies (γ 1/2 ) mi γ ij (γ 1/2 ) jℓ = δ mℓ .Then d 3 p (i) = √ γ d 3 k(i) , and k 2 = k2 = δ ij kikj , which is even in the components ki .One could perform an additional coordinate transformation as in Ref. [52], p ′ i = p i − β 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, β i = 0. Then the integration in Eq. ( 38) may be performed exactly.In such cases, we find The quantity N (x, k) is even in the components ki , so that j i a (x) vanishes identically.We further note that the momentum coordinates ki are Cartesian, so j 0 a (x) is equivalent to with Eq. ( 41) corresponds to The induced current in Eq. ( 43) is equivalent to inserting an effective mass m D (x) for the A a 0 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) → 1 and identifying 1/β T = T with the temperature of the plasma [1-4, 8, 16, 17, 63-66].In our case, the Boltzmann factor for the high-momentum particles is exp ], since we are considering spacetimes in which β i = 0.This is exactly what one would expect for particles in a fluid that is undergoing normal flow, with four-velocity n µ (x) and therefore local energy per particle E ≡ −n µ p µ = k, if we identify T (x) ≡ 1/(β T N (x)) with the local temperature.We may therefore identify the constant β T ≡ 1/T 0 and write Given that we have normalized t such that N (x) → 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 [28][29][30].The Debye mass then takes the form 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.

A. Debye screening in Minkowski spacetime
We first consider the behavior of soft gluon modes A 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, β i = 0.In the static limit, ∂ 0 A a µ (x) = 0, the exact equations of motion of Eq. ( 5) reduce to where repeated indices are summed and we have adopted Cartesian coordinates for the spatial sections.The induced current of Eq. ( 43) reduces to A generic feature of non-Abelian field theories, which has been well-studied for the case of self-interacting gauge fields in Minkowski spacetime, is the existence of monopole-like solutions among the spatial components A a i (x).Such solutions can be found in simplest form for SU(2) [1,[67][68][69], and have been generalized to monopoles charged under various 3-dimensional subgroups of SU(N c ) [70][71][72].Moreover, it is well-known that self-consistent, static solutions to the non-Abelian equations of motion in Minkowski spacetime generically include a Yukawa-like screened behavior for the component A a 0 (x) combined with the Wu-Yang monopole solution for the components A a i (x) [1,69].Such solutions underscore the important physical point that only the components A a 0 (x) acquire a nonzero mass within a medium in the state described in Section IV-as indicated by the form of the induced current j a µ (x) in Eq. ( 47)-and hence only those components undergo screening, with amplitude proportional to exp[−m D r].
Our aim in this section is to examine how well-known field configurations in Minkowski spacetime generalize to various curved spacetimes of interest, in which additional length-scales become relevant.We therefore begin by considering an exact solution to Eq. ( 46) that consists of a superposition of a screened component A a 0 with a Wu-Yang monopole solution for the components A a i .For the case of SU(2), the solution takes the form [1,[67][68][69]] where r is a unit vector, and ϵ ijk is the usual Levi-Civita symbol.Note that the solution mixes spatial indices and color-space indices, which is straightforward for SU(2), since both i, j = 1, 2, 3 and a, b = 1, 2, 3. (To confirm that A a µ (x) in Eq. ( 48) satisfies Eq. ( 46), it is helpful to construct the projection operator P a b = δ a b − xa xb , noting that xa P a b = 0 = xb P a b and ∂ j xi = r −1 P i j .)For SU(3), exact monopole solutions have been found for subgroups such as U(1) × U(1) and U(2), for which the components A a i have comparable asymptotic behavior to the solution in Eq. ( 48) [70][71][72].
The solution for A a µ (x) in Eq. ( 48) carries both chromoelectric charge Q and chromomagnetic charge P. We consider the quasi-local charges defined in terms of the chromoelectric and chromomagnetic fields, E a i ≡ F a 0i and B a i ≡ − 1 2 ϵ ijk F jk a , respectively.Specifically, we identify charges via the relations where the spatial and color indices are summed over.Note that the fields E a i and B a i transform covariantly in color space under gauge transformations, and therefore the bilinear combinations in Eq. ( 49) are gauge-invariant.For the solution given in Eq. ( 48), the charges are where is the chromoelectric charge measured at the origin.Because of screening, the value of Q(r) measured at a finite distance from the origin will be reduced compared to Q 0 .The chromomagnetic charge, on the other hand, is not screened, and hence P = 1/g s everywhere.
Whereas the magnitude of the chromomagnetic charge P is fixed by 1/g s , the chromoelectric charge 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 ≤ r ≤ r abel , within which |A a 0 (x)| ≫ |A a i (x)|, with r abel given by For 0 ≤ r ≤ r abel , the system is quasi-Abelian, with A a µ (x) ≃ A a 0 (x)δ 0 µ , and hence non-Abelian contributions of the form f abc A a µ A b ν remain subdominant.Within this regime, one may solve the exact equations of motion of Eqs. ( 46)- (47) with the ansatz A a µ (x) = A a 0 (x)δ 0 µ .The solution for A a µ (x) is then of the form in Eq. ( 48), but with A a i (x) = 0 and Q a 0 an arbitrary, constant vector in color space.

B. Debye screening in an FLRW spacetime
We next consider Debye screening within a hot quarkgluon plasma in an expanding FLRW spacetime at early times, prior to the QCD confinement transition (t < t QCD ∼ 10 −5 s).The line-element for a spatially flat FLRW universe may be written where a(t) is the scale factor and r is a comoving radial coordinate.The Hubble parameter is given by H(t) ≡ ȧ/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 hard and k soft .Note also that in these coordinates the lapse function is N (x) = 1, so the temperature has no spatial gradients.The term L 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 ρ and pressure P associated with the high-momentum particles dominate T µν , so that the Friedmann equation takes the form 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 topquark mass (m t = 173 GeV), g * = 106.75[22,23].During the radiation-dominated phase, a(t) ∝ t 1/2 , so we have H(t) = 1/(2t) and consistent with adiabatic expansion, T (t) a(t) = constant.
The equations of motion from Eq. ( 5) are where repeated indices are summed.The Hubble parameter H(t) = ȧ/a sets the 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(t)) are either proportional to H(t) or include a time derivative Ȧc µ .We consider an ansatz for A c µ (x) in which the only time dependence arises from the scale factor a(t).By the chain rule, dA c µ dt = ȧ ∂A c µ ∂a = H ∂A c µ ∂ ln a .Therefore any term with a time derivative Ȧc µ will be proportional to H.If H ≪ k 2 soft /k hard ∼ g 2 s T , then all terms proportional to H, which include those with time derivatives Ȧc µ , will be subleading compared to s T is satisfied at early times, it will be even more easily satisfied at later times, since H ∼ T 2 /M pl from Eq. ( 54) and T /M pl decreases as the universe cools down, while g s gradually runs to larger values.Now we proceed to show that H ≪ g 2 s T holds at early times of interest.Consider, as an example, dynamics at the time t c ∼ 10 −21 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 [24][25][26][27].At such early times, the Hubble scale H(t c ) = 3.0×10 −4 GeV and the fluid filling the FLRW spacetime had a temperature T (t c ) = 10 7 GeV.At such high energy scales, we must take into account the running of the strong coupling α s ≡ g 2 s /(4π).To two-loop order, the running with energy scale µ is given by [4] dα s (µ) with Upon normalizing α s (m Z ) = 0.118 at the Z boson mass (m Z = 91.2GeV), we find α s = 0.046 at the energy scale µ = T (t c ) = 10 7 GeV.This yields a ratio g 2 s T /H ∼ 10 9 .Hence we may self-consistently neglect H ≪ g 2 s T when considering the dynamics of the gluon modes on lengthscales λ 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 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 → 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 a 0 with a Wu-Yang monopole for the components A a i , which in the case of SU(2) takes a similar form to Eq. ( 48), where we have rescaled the spatial coordinates back so that r is the comoving radial coordinate from Eq. ( 52).The solution for A a µ (x) in Eq. ( 58) remains consistent with the quasi-stationary approximation of Eq. ( 24) for λ D ≤ r ≪ r H .

C. 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πM 2  pl /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 [24][25][26][27].
Any PBHs that formed at times t < t QCD = 10 −5 s would undergo collapse amid a hot plasma of unconfined quarks and gluons.On long length-scales at such early times, λ ≫ 1/k soft , we expect that the plasma would have attained a charge-neutral equilibrium distribution.Yet on shorter length-scales, set by λ D = 1/m D ∼ 1/k soft ∼ 1/(g s T ), spatial regions of nonvanishing net color charge can form, within which the charges for most soft gluons align in color space [18,19].In such scenarios, the PBHs would form by absorbing one or more net-color regions.Depending on the ratio of the PBH radius (∼ GM PBH ) and the Debye length (∼ 1/m D ) at the time of formation, PBHs therefore could form with net color charge Q 0 .Such a scenario is distinct from the examples reviewed in Ref. [73,74], 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 Q 0 ; the case of larger Q 0 is treated in Ref. [75].In particular, we consider the regime Within this regime, backreaction on spacetime from the energy associated with the enclosed charge Q 0 remains subdominant, and we may consider the dynamics of soft gluon modes A a µ (x) within a fixed background geometry.(For Q 2 0 ∼ GM 2 PBH , one must solve the Einstein field equations as well as the equations of motion for A a µ (x), which we consider in separate work [76].) We consider a spherically symmetric spacetime with a black hole of mass M 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 PBH not change appreciably over timescales set by 1/k soft .Two competing effects could change M PBH : evaporation due to Hawking radiation (which would reduce M PBH over time) and accretion from the surrounding medium (which would increase M PBH over time).For the regime of interest, we find that both of these effects are negligible over the relevant dynamical time-scale, and hence we may neglect ṀPBH .
Consider first evaporation from Hawking radiation.Because we are considering PBHs with modest charge, subject to the inequality of Eq. ( 59), we may approximate the Hawking temperature based on that for an (uncharged) Schwarzschild black hole of mass M PBH .(Any net enclosed charge would decrease the black hole's surface gravity, and hence its Hawking temperature, compared to the zero-charge case, thus rendering evaporation even less efficient.)The Hawking temperature for a Schwarzschild black hole is given by [59] The typical mass for a PBH is set by the mass enclosed within a Hubble volume M H (t c ) at the time of collapse t c [77][78][79][80]: where γ ≃ 0.2 and T c is the temperature of the plasma at the time t c .The values of M PBH relevant to account for dark matter lie within the range 10 17 g ≤ M PBH ≤ 10 22 g [24-27]; from Eq. ( 61), these correspond to plasma temperatures 10 5 GeV ≤ T c ≤ 10 7 GeV.Meanwhile, for masses within the dark-matter range, Eq. ( 60) yields 10 −9 GeV ≤ T H ≤ 10 −4 GeV, exponentially lower than the temperature of the surrounding plasma.Such PBHs will therefore be net absorbers (rather than emitters) around the time of their formation.
Next we may consider accretion.One might suppose that if the fluid were moving with respect to the black hole, M 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 PBH would only increase by a fraction ∆M/M PBH ∼ O(1) R T /M pl over an entire Hubble time, where R is the ratio of the black hole radius to the Debye screening length and T the temperature of the plasma.As in the previous section, we are interested in early times at which the plasma temperature could have been as high as T ∼ 10 7 GeV, in which case T /M pl ∼ 10 −11 .We consider black holes for which the ratio R is not exponentially larger than one, so that the PBH forms with some small, residual color charge, and therefore ∆M/M remains negligible.A more formal calculation of the Bondi accretion rate for PBHs in this scenario yields the same conclusion [81][82][83][84].
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 [85][86][87].A convenient parameterization may be written [87] with and r s ≡ 2GM PBH the usual Schwarzschild radius.
We found in Section V B that in cosmological scenarios of interest, H ≪ m D , even at very early times.For the remainder of this section, we will therefore set H ∼ 0, for which a(t) ∼ constant (which we will scale to 1).In that limit, the lapse function reduces to N (x) → f (r), with f (r) = 1 − r s /r, and the shift vector vanishes, β i → 0.
Next we consider an appropriate range for the enclosed charge Q 0 , subject to the constraint in Eq. ( 59).The shortest scales that can be resolved within our EFT are set by λ D = 1/m D , so we restrict r s ≥ λ D , which in turn requires M PBH ≥ 1/(2Gm D ).Then Eq. ( 59) corresponds to the upper bound around t c ∼ 10 −21 s.Meanwhile, as discussed around Eq. ( 51), if Q 0 ≫ 1/g s , then the soft modes A a µ will assume a quasi-Abelian form, with |A a 0 | ≫ |A a i |, for r < r abel .The estimate of r 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) ∼ g s T 0 / f (r).Following the same steps that led to Eq. ( 51), we find where λD is the Debye length associated with the temperature T 0 .For r s / λD = 3.5 and Q 0 ∼ 10 4 , this yields r abel ∼ 3 r s , while Q 0 ∼ 10 10 yields r abel ∼ 7 r s .The last departure to consider from the previous examples concerns the effect of the local temperature T (x) on the coupling α 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 α s with energy scale µ → T (r).To resolve the behavior of α s near r s , we adopt the "tortoise" radial coordinate, r * ≡ r + r s ln[(r/r s ) − 1], for which r = r s corresponds to r * → −∞ [59].See Fig. 1.
In the vicinity of a primordial black hole, the evolution of the soft modes A a µ (x) is therefore characterized by a hierarchy of length-scales: For r s ≤ r ≤ r abel , the soft modes will evolve as a quasistatic, quasi-Abelian system within a fixed background spacetime.In that regime, A a µ (x) ≃ A a 0 (r)δ 0 µ +O(r/r abel ), and the equations of motion in Eq. ( 5) reduce to where mD = 1/ λD is the Debye mass associated with the temperature T 0 .Including the running of g s with r, we The QCD coupling strength αs ≡ g 2 s /(4π) 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 * ≡ r + rsln[(r/rs) − 1]; r * /rs = −36 corresponds to r/rs = 1 + 10 −16 .We have set T0 = 10 7 GeV on the asymptotic boundary.
2. The soft mode A a 0 (r) (solid line) for rs/ λD = 3.5 and T0 = 10 7 GeV.Given the greater density of plasma near the black hole, Debye screening is particularly effective for r ≳ rs.For r rs, the behavior of A a 0 (r) asymptotes to the expected slope, exp[−mDr]/r (dashed line).may solve Eq. ( 67) numerically, with a typical example shown in Fig. 2. The component A a 0 (r) undergoes strong screening for r ≳ r s , since more plasma gathers near the event horizon, yielding a higher density and temperature than at locations r ≫ r s .Far from the black hole, A a 0 (r) asymptotes to similar behavior as found in Sections V A and V B, with A a 0 ∼ exp[−m D r]/r.Note that an observer at r ≳ r s would measure an effective charge Q(r) much smaller than the charge Q 0 contained within the black hole.
Within our coordinate system, the Tolman temperature gradient drives T (r) → ∞ as r → 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.[88][89][90] and evaluated the field with a boundary condition at a "stretched horizon," r s + ϵ, rather than at r s , with ϵ/r s = 10 −6 .In forthcoming work [76], we address this limitation by solving the coupled Einstein field equations and equation of motion for A a 0 (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 λ D (x) = 1/m D (x) is set by a Debye mass m D (x) ∼ 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 [28][29][30].
To characterize Debye screening, we constructed an effective theory for long-wavelength excitations in a hot Yang-Mills plasma.We analyzed the dynamics of highmomentum particles in the plasma (with momentum k hard ∼ 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 a 0 components of the soft modes (with momentum k soft ∼ g s T ).
We applied our results to solve for the gauge potential A a µ in a few cases of interest.We recovered the wellknown Wu-Yang monopole solution in the Minkowski limit and generalized it FLRW spacetimes, demonstrating the self-consistency of the quasi-static approximation in regimes of interest, for which cosmological time-scales may be neglected compared to 1/k soft .
Lastly, we analyzed Debye screening in the vicinity of a primordial black hole immersed in a hot quark-gluon plasma.Such black holes can form by absorbing regions of the plasma with nonvanishing net color charge, with characteristic size set by λ D ; hence the resulting dial 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 a µ (x) ≈ A a 0 (r)δ 0 µ .Given the Tolman temperature gradients, the interaction strength α 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 a 0 (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 [76].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 [75].