format: text

text:
# Transport coefficients in high temperature gauge theories: (II) Beyond leading log

## I Introduction

Transport coefficients, such as viscosities, diffusivities, or electric conductivity, characterize the dynamics of long wavelength, low frequency fluctuations in a medium. In a weakly coupled quantum field theory, transport coefficients should, in principle, be calculable purely theoretically. The values of transport coefficients are of interest in cosmological applications such as electroweak baryogenesis baryo1 ; baryo2  and the origin of primordial magnetic fields mag\_fields , as well as for hydrodynamic models of heavy ion collisions (see, for example, Refs. teaney ; heavy-ion1 ; heavy-ion2 ; heavy-ion3 ; heavy-ion4a ; heavy-ion4b ; heavy-ion4c , and references therein). From a purely theoretical perspective, the evaluation of transport coefficients also provides an excellent test of our understanding of dynamic processes in thermal field theory.

In a previous paper AMY1 , we performed leading-log calculations of the shear viscosity, electrical conductivity, and flavor diffusion constants in weakly coupled, high temperature gauge theories111The leading-log result for shear viscosity in pure Yang-Mills theory was first evaluated by Heiselberg Heiselberg  and closely approximated by earlier work of Baym et al. BMPRa . A variety of incorrect leading-log calculations of transport coefficients are also present in the literature; see Ref. AMY1  for a discussion of this earlier work. —that is, neglecting relative corrections suppressed only by powers of the inverse logarithm of the gauge coupling, $`1/\mathrm{log}(1/g)`$. Leading-log calculations may be regarded as improvements over phenomenological estimates based on relaxation time approximations (see, for example, Refs. HosoyaKajantie ; Hosoya\_and\_co ; relax1 ; relax2 ; relax3 ; relax4 ), but no leading-log calculation can be trusted to provide even a factor of two estimate in any real application, because the logarithm of the inverse gauge coupling, even for electromagnetism, is never all that large.

At a minimum, one would like to know several terms in the expansion in inverse powers of $`\mathrm{log}(1/g)`$ in order to assess the utility of this asymptotic series. Even better would be a full leading-order calculation of transport coefficients, by which we mean an evaluation which correctly includes all effects suppressed by powers of $`1/\mathrm{log}(1/g)`$ and only neglects relative corrections suppressed by powers of $`g`$. This is a feasible goal. Transport coefficients are dominantly sensitive222 This is true of the transport coefficients under consideration in this paper. It should be noted that bulk viscosity requires a separate treatment because it vanishes identically in a scale invariant theory, and so requires the correct incorporation of effects which are irrelevant for other transport coefficients. For a discussion of bulk viscosity in the case of scalar field theory, see Refs. Jeon ; JeonYaffe .  to the dynamics of excitations (i.e., quarks and gluons) with typical momenta of order of the temperature $`T`$. And, as we recently discussed in Ref. AMY5 , it is possible to formulate an effective kinetic theory which correctly describes the leading-order dynamics of such excitations. The purpose of this paper is to perform a full leading-order evaluation of shear viscosity, electrical conductivity, and flavor diffusion constants, and also to examine the utility of the asymptotic expansion in inverse powers of $`\mathrm{log}(1/g)`$. We believe this work represents the first reliable analysis of transport coefficients beyond a leading-log approximation in high temperature gauge theories.

We begin in Section II by summarizing the relevant definitions, reviewing the form of the effective kinetic theory, and sketching how the actual calculation of transport coefficients reduces to a linear integral equation which may be accurately solved using a variational approach. The evaluation of matrix elements of the linearized collision operator is discussed in Appendices A and B. Our results for the leading-order fermion diffusion constant and shear viscosity in QCD are presented in Section III, and the expansion in powers of $`1/\mathrm{log}(1/g)`$ is discussed in Section IV. We include here next-to-leading-log results for both Abelian and non-Abelian theories. We find that a next-to-leading-log approximation reproduces the full leading-order results substantially better than one might have expected. Section V presents both leading-order, and next-to-leading-log results for the electrical conductivity in QED (or in quark plus lepton) plasmas. A short conclusion summarizes our results and comments briefly on the relation between our effective kinetic theory treatment and purely diagrammatic approaches. Several more technical appendices follow. The first two present details of the required $`22`$ particle effective scattering rates, and the $`12`$ particle effective splitting rates which characterize nearly collinear bremsstrahlung and pair production processes in the fluctuating soft thermal gauge field background. A final appendix contains a proof that the expansion of leading-order transport coefficients in inverse powers of $`\mathrm{log}(1/g)`$ is only asymptotic, not convergent, and gives a non-rigorous argument suggesting that this expansion is not Borel summable.

## II Ingredients

Throughout this work, we assume that the gauge coupling (defined at the scale of the temperature) is weak, $`g(T)1`$. For QCD, this means that the temperature $`T`$ is asymptotically large compared to $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. We assume that there are no particle masses close to $`T`$. Zero-temperature masses must either be negligible compared to $`gT`$ (so thermal self-energies are large compared to the mass), or large compared to $`T`$ (so the particle decouples and may be ignored completely). In particular this means we will not consider temperatures just below thermal phase transitions (or crossovers). We also assume that any chemical potentials for conserved numbers are small compared to $`T`$.

### II.1 Definitions

Transport coefficients characterize a system’s response to weak, slowly varying inhomogeneities or external forces. If a weak and spatially uniform electric field $`𝑬`$ is applied to a plasma, the resulting induced electric current is

$$𝒋^{\mathrm{EM}}=\sigma 𝑬,$$
(1)
where the electric conductivity $`\sigma `$ is a function of the temperature $`T`$, the electric charge $`e`$, and the particle content of the theory. This constitutive relation is satisfied up to corrections involving gradients of the electric field and higher powers of the field strength; strictly speaking the conductivity is defined by the above relation in the limit of vanishing frequency and wavenumber of an infinitesimal applied electric field.

In a theory with a conserved global charge (such as baryon or lepton number), the associated charge density $`nj^0`$ and current density $`𝒋`$ will satisfy a diffusion equation,

$$𝒋=D\mathbf{}n,$$
(2)
in the local rest frame of the medium.333We work throughout in the local rest frame of some arbitrary point $`x`$ in the system. This frame is defined by the Landau-Lifshitz convention of vanishing momentum density, $`T^{0i}(x)=0`$.  The coefficient $`D`$ is called the diffusion constant.

When one or more diffusing species of excitations carry electric charge, the electric conductivity is connected to the diffusion constants for various species by an Einstein relation (see Ref. AMY1  for more discussion). If the net number of each species of charge carriers is conserved, then

$$\sigma =\underset{a}{}e_a^2D_a\frac{n_a}{\mu _a},$$
(3)
where the sum runs over the different species or flavors of excitations with $`e_a`$ and $`D_a`$ the corresponding electric charge and diffusion constant, respectively.444For a single massless Dirac fermion at vanishing chemical potential, the susceptibility $`n_a/\mu _a=\frac{1}{3}T^2`$ (at leading order in the gauge coupling). For quarks, this should be multiplied by 3 (or the number of colors). However, the electrical conductivity is dominated by the transport of charged leptons, which diffuse more readily than quarks since they only interact electromagnetically.

If the flow velocity of the plasma is not uniform, then the stress tensor (which defines the flux of momentum density) will depart from its perfect fluid form. In the local fluid rest frame at a point $`x`$, the stress tensor, to lowest nontrivial order in the velocity gradient, will have the form

$$T_{ij}(x)=\delta _{ij}𝒫\eta \left[_iu_j+_ju_i\frac{2}{3}\delta _{ij}^lu_l\right]\zeta \delta _{ij}^lu_l,$$
(4)
where $`𝒫`$ is the equilibrium pressure associated with the energy density $`T_{00}(x)`$, and the coefficients $`\eta `$ and $`\zeta `$ are known as the shear and bulk viscosities, respectively. The flow velocity $`𝒖`$ equals the momentum density divided by the enthalpy density (which is the sum of energy density and pressure). We will only be concerned with the shear viscosity in this paper; the bulk viscosity requires a more complicated analysis, which to date has only been performed for a scalar field theory Jeon ; JeonYaffe .

### II.2 Effective kinetic theory

In a weakly coupled plasma, the energy, momentum, electric charge, and other global conserved charges, are predominantly carried by excitations with momenta $`pgT`$. Such momenta will be referred to as “hard.” Weak coupling implies that these excitations are long-lived (and hence well-defined) quasiparticles. Thermal corrections to the dispersion relations of excitations are order $`gT`$ in size \[i.e., $`(p^0)^2=𝒑^2+O(g^2T^2)`$\]. Hence, hard excitations are ultrarelativistic.

The parametrically large separation between quasiparticle lifetimes and the screening length of order $`1/(gT)`$ ensures that the dynamics of these excitations may be reproduced by an effective kinetic theory.555 In particular, the contributions of hard degrees of freedom to the flux of energy, momentum, or other global conserved charges are captured by the two-point Wigner functions for the hard degrees of freedom, up to sub-leading corrections. As discussed in, for instance, Refs. Jeon ; JeonYaffe ; Blaizot&Iancu ; CalzettaHu1 ; CalzettaHu2 ; A&Y , kinetic theory correctly describes the time evolution of the two-point function to within $`O(g)`$ or smaller corrections.  In other words, one may describe the instantaneous state of the system by a set of phase space distribution functions $`f^a(𝒙,𝒑,t)`$. The index $`a`$ will label the relevant species of excitations (i.e., gluon, up-quark, up-antiquark, down-quark, down-antiquark, etc., for QCD). Formally, the distribution functions represent spatially smeared and Wigner transformed non-equilibrium two-point correlation functions. The evolution of phase space distribution functions is governed by effective Boltzmann equations of the form

$$\left[\frac{}{t}+𝒗_𝒑\frac{}{𝒙}+𝑭_{\mathrm{ext}}^a\frac{}{𝒑}\right]f^a(𝒑,𝒙,t)=C_a[f],$$
(5)
where $`𝒗_𝒑p^0/𝒑`$ is the velocity of an excitation with momentum $`𝒑`$, $`𝑭_{\mathrm{ext}}^a(𝒙,t)`$ is an external force acting on excitations of type $`a`$ (due, for example, to an applied electric field), and $`C_a[f]`$ is the collision integral which characterizes the net rate at which species $`a`$ excitations with momentum $`𝒑`$ are lost (or created) due to scattering processes involving other excitations in the plasma. For ultrarelativistic hard excitations, the velocity is a unit vector, $`𝒗_𝒑=\widehat{𝒑}`$, up to irrelevant subleading corrections.

As discussed at length in Ref. AMY5  and proposed earlier by Baier, Mueller, Schiff, and Son BMSS ,666See also Refs. Gelis1 ; AMY2 ; AMY3 ; AMY4 ; AGMZ  for discussion of the closely related processes which contribute to the leading-order photon emission rate.  the collision term needed to reproduce leading-order quasiparticle dynamics must contain both $`22`$ particle scattering terms and effective $`12`$ particle splitting and joining terms. The latter represent nearly-collinear bremsstrahlung and pair production/annihilation processes which take place in the presence of soft ($`gT`$ scale) thermal fluctuations in the background gauge field. Hence, the collision terms take the form

$$C_a[f]=C_a^{22}[f]+C_a^{12}[f],$$
(6)
where

$`C_a^{22}[f](𝒑)`$ $`=`$ $`{\displaystyle \frac{1}{4|𝒑|\nu _a}}{\displaystyle \underset{bcd}{}}{\displaystyle _{𝒌𝒑^{}𝒌^{}}}\left|_{cd}^{ab}(𝒑,𝒌;𝒑^{},𝒌^{})\right|^2(2\pi )^4\delta ^{(4)}(P+KP^{}K^{})`$ (7)
$`\times \{f^a(𝒑)f^b(𝒌)[1\pm f^c(𝒑^{})][1\pm f^d(𝒌^{})]`$
$`f^c(𝒑^{})f^d(𝒌^{})[1\pm f^a(𝒑)][1\pm f^b(𝒌)]\},`$
and
$`C_a^{12}[f](𝒑)`$ $`=`$ $`{\displaystyle \frac{(2\pi )^3}{2|𝒑|^2\nu _a}}{\displaystyle \underset{bc}{}}{\displaystyle _0^{\mathrm{}}}𝑑p^{}𝑑k^{}\delta (|𝒑|p^{}k^{})\gamma _{bc}^a(𝒑;p^{}\widehat{𝒑},k^{}\widehat{𝒑})`$ (8)
$`\times \left\{f^a(𝒑)[1\pm f^b(p^{}\widehat{𝒑})][1\pm f^c(k^{}\widehat{𝒑})]f^b(p^{}\widehat{𝒑})f^c(k^{}\widehat{𝒑})[1\pm f^a(𝒑)]\right\}`$
$`+`$ $`{\displaystyle \frac{(2\pi )^3}{|𝒑|^2\nu _a}}{\displaystyle \underset{bc}{}}{\displaystyle _0^{\mathrm{}}}𝑑k𝑑p^{}\delta (|𝒑|+kp^{})\gamma _{ab}^c(p^{}\widehat{𝒑};𝒑,k\widehat{𝒑})`$
$`\times \left\{f^a(𝒑)f^b(k\widehat{𝒑})[1\pm f^c(p^{}\widehat{𝒑})]f^c(p^{}\widehat{𝒑})[1\pm f^a(𝒑)][1\pm f^b(k\widehat{𝒑})]\right\}.`$
Here, $`P=(p^0,𝒑)`$, $`K=(k^0,𝒌)`$, etc. denote null four-vectors (so that $`p^0|𝒑|`$, etc.), $`\nu _a`$ is the number of spin times color states for species $`a`$ (i.e., 6 for each quark or antiquark and 16 for gluons), and $`_𝒑d^3𝒑/[2|𝒑|(2\pi )^3]`$ denotes Lorentz-invariant momentum integration. As usual, upper signs refer to bosons and lower signs to fermions. $`_{cd}^{ab}(𝒑,𝒌;𝒑^{},𝒌^{})`$ is the effective two body scattering amplitude for the process $`abcd`$, defined with a relativistic normalization for single particle states; its square $`\left|_{cd}^{ab}\right|^2`$ is implicitly understood to be summed (not averaged) over the spins and colors of all four excitations. Similarly, $`\gamma _{cd}^a(𝒑;p^{}\widehat{𝒑},k^{}\widehat{𝒑})`$ is the differential rate for an $`abc`$ effective splitting process (or its time reverse), integrated over the parametrically small transverse momenta of the participants. This rate is likewise understood to be summed over the spins and colors of all three participants. The prefactors of $`1/(4|𝒑|)`$ in the $`22`$ terms are a combination of the $`1/(2|𝒑|)`$ from the relativistic normalization of scattering amplitudes together with a symmetry factor of 1/2 which corrects for double counting of final or initial states (specifically $`𝒑^{},c`$ interchanged with $`𝒌^{},d`$).

The $`22`$ matrix elements $`|_{cd}^{ab}|^2`$ are the usual lowest order vacuum amplitudes, except that thermal self-energies must be included on internal gauge boson or fermion lines in $`t`$ and $`u`$ channel exchange processes AMY5 . These thermal self-energies are evaluated in the hard thermal loop (HTL) approximation, which is adequate for a leading-order analysis. Explicit expressions may be found in Ref. AMY5  and also appear in Appendix A, which discusses the evaluation of the $`22`$ collision terms.

The $`12`$ splitting rates $`\gamma _{cd}^a`$ encapsulate the effect of near-collinear processes in which any number of soft scatterings with other excitations in the plasma occur during the emission. The importance of including such multiple soft scatterings is known as the Landau-Pomeranchuk-Migdal (LPM) effect. It leads to a two-dimensional linear integral equation which must be solved, for each value of $`|𝒑|`$ and $`k`$, to determine the $`12`$ transition rate $`\gamma _{bc}^a`$. This is discussed further in Ref. AMY5  and in Appendix B, which describes the evaluation of these $`12`$ splitting rates.

### II.3 Linearization

Given an effective kinetic theory, transport coefficients are calculated by linearizing the Boltzmann equation about local equilibrium. The required procedure is discussed in some detail in our earlier paper AMY1 . In brief, one writes each distribution function as a slowly varying local equilibrium part $`f_0^a(𝒙,𝒑)`$ plus a small departure from (local) equilibrium $`f_1^a(𝒙,𝒑)`$. For electric conductivity, $`f_0^a`$ may equal a homogeneous equilibrium distribution. For diffusion, $`f_0^a`$ is a local equilibrium distribution in which the chemical potential coupled to the conserved charge of interest varies slowly in space. For shear viscosity, $`f_0^a`$ is a local equilibrium distribution in which the flow velocity has a non-zero shear. In each case, the leading contribution on the left side of the Boltzmann equation (5) comes from either the convective derivative or external force term acting on $`f_0^a`$. Because the collision terms vanish identically for any local equilibrium distribution, the leading contribution on the right side of the equation is linear in the departure $`f_1^a`$ from equilibrium. One solves the resulting linearized Boltzmann equation for $`f_1^a`$, uses the result to evaluate the stress tensor or current density of interest, and then reads off the appropriate transport coefficient.

The convective derivative acting on $`f_0^a`$ is proportional to the relevant driving term, which we will denote by

$$X_{i\mathrm{}j}\{\begin{array}{cc}E_i,\hfill & \text{(conductivity),}\hfill \\ _i\mu _\alpha ,\hfill & \text{(diffusion),}\hfill \\ \frac{1}{\sqrt{6}}\left(_iu_j+_ju_i\frac{2}{3}\delta _{ij}𝒖\right),\hfill & \text{(shear viscosity).}\hfill \end{array}$$
(9)
This tensor has $`\mathrm{}=1`$ angular dependence for conductivity and diffusion, and $`\mathrm{}=2`$ angular dependence for shear viscosity. The departure from equilibrium which solves the linearized Boltzmann equation must be proportional to $`X_{i\mathrm{}j}`$. It is convenient to express the departure in the local rest frame (where $`f_0^a`$ is isotropic) as

$$f_1^a(𝒑)=\beta ^2f_0^a(p)[1\pm f_0^a(p)]X_{i\mathrm{}j}\chi _{i\mathrm{}j}^a(𝒑),$$
(10)
which defines $`\chi _{i\mathrm{}j}^a(𝒑)`$. The resulting linearized Boltzmann equation for the functions $`\chi _{i\mathrm{}j}^a(𝒑)`$ may be written compactly in the form

$$𝒮_{i\mathrm{}j}^a(𝒑)=(𝒞\chi _{i\mathrm{}j})^a(𝒑),$$
(11)
where $`𝒞`$ is a linearized collision operator defined below. The source term is

$$𝒮_{i\mathrm{}j}^a(𝒑)Tq^af_0^a(p)[1\pm f_0^a(p)]I_{i\mathrm{}j}(\widehat{𝒑}).$$
(12)
Here $`q^a`$ denotes the relevant conserved charge carried by species $`a`$ associated with the transport coefficient of interest,

$$q^a\{\begin{array}{cc}e_{\mathrm{EM}}^a,\hfill & \text{(conductivity),}\hfill \\ q_\alpha ^a,\hfill & \text{(diffusion),}\hfill \\ |𝒑|,\hfill & \text{(shear viscosity),}\hfill \end{array}$$
(13)
where in the case of diffusion, $`\alpha `$ is simply a label for the flavor symmetry of interest (e.g. quark number or lepton number). $`I_{i\mathrm{}j}`$ is the unique $`\mathrm{}=1`$ or $`\mathrm{}=2`$ rotationally covariant tensor depending only on $`\widehat{𝒑}`$,

$$I_{i\mathrm{}j}(\widehat{𝒑})\{\begin{array}{cc}\widehat{p}_i,\hfill & \mathrm{}=1\text{ (conductivity/diffusion),}\hfill \\ \sqrt{\frac{3}{2}}(\widehat{p}_i\widehat{p}_j\frac{1}{3}\delta _{ij}),\hfill & \mathrm{}=2\text{ (shear viscosity).}\hfill \end{array}$$
(14)
The normalization on $`I_{i\mathrm{}j}`$ was chosen so that

$$I_{i\mathrm{}j}(\widehat{𝒑})I_{i\mathrm{}j}(\widehat{𝒑})=1$$
(15)
(which is the reason for the peculiar normalization of $`X_{ij}`$ in the shear viscosity case). More generally,

$$I_{i\mathrm{}j}(\widehat{𝒑})I_{i\mathrm{}j}(\widehat{𝒌})=P_{\mathrm{}}(\widehat{𝒑}\widehat{𝒌}),$$
(16)
where $`P_{\mathrm{}}(x)`$ is the $`\mathrm{}`$’th Legendre polynomial. Rotational invariance of the collision operator (in the local rest frame) implies that

$$\chi _{i\mathrm{}j}^a(𝒑)=I_{i\mathrm{}j}(\widehat{𝒑})\chi ^a(p)$$
(17)
for some scalar function $`\chi ^a(p)`$. Here and throughout, $`p|𝒑|`$. Solving the linearized Boltzmann equation means determining $`\chi ^a(p)`$.

Let $`N_s`$ be the number of relevant species, and define an inner product on $`N_s`$-component functions of momenta,

$$(f,g)\beta ^3\underset{a}{}\nu _a\frac{d^3𝒑}{(2\pi )^3}f^a(𝒑)g^a(𝒑).$$
(18)
One may show that the linearized collision operator $`𝒞`$ is symmetric with respect to this inner product. It is a positive semi-definite operator, which is strictly positive definite in the $`\mathrm{}=1`$ and $`\mathrm{}=2`$ channels relevant for diffusion or shear viscosity. Consequently, the linearized Boltzmann equation (11) is precisely the condition for maximizing the functional

$$Q[\chi ](\chi _{i\mathrm{}j},𝒮_{i\mathrm{}j})\frac{1}{2}(\chi _{i\mathrm{}j},𝒞\chi _{i\mathrm{}j}).$$
(19)
The explicit forms of the source and collision parts of this quadratic functional are

$`(\chi _{i\mathrm{}j},𝒮_{i\mathrm{}j})`$ $`=\beta ^2{\displaystyle \underset{a}{}}{\displaystyle \frac{d^3𝒑}{(2\pi )^3}f_0(p)[1\pm f_0(p)]q^a\chi ^a(p)},`$ (20)
and
$`(\chi _{i\mathrm{}j},𝒞\chi _{i\mathrm{}j})`$ $`=(\chi _{i\mathrm{}j},𝒞^{22}\chi _{i\mathrm{}j})+(\chi _{i\mathrm{}j},𝒞^{12}\chi _{i\mathrm{}j}),`$ (21)
with
$`(\chi _{i\mathrm{}j},𝒞^{22}\chi _{i\mathrm{}j})`$ $`{\displaystyle \frac{\beta ^3}{8}}{\displaystyle \underset{abcd}{}}{\displaystyle _{𝒑𝒌𝒑^{}𝒌^{}}}\left|_{cd}^{ab}(𝒑,𝒌;𝒑^{},𝒌^{})\right|^2(2\pi )^4\delta ^{(4)}(P+KP^{}K^{})`$
$`\times f_0^a(p)f_0^b(k)[1\pm f_0^c(p^{})][1\pm f_0^d(k^{})]`$
$`\times \left[\chi _{i\mathrm{}j}^a(𝒑)+\chi _{i\mathrm{}j}^b(𝒌)\chi _{i\mathrm{}j}^c(𝒑^{})\chi _{i\mathrm{}j}^d(𝒌^{})\right]^2,`$ (22)
$`(\chi _{i\mathrm{}j},𝒞^{12}\chi _{i\mathrm{}j})`$ $`{\displaystyle \frac{\beta ^3}{2}}{\displaystyle \underset{abc}{}}4\pi {\displaystyle _0^{\mathrm{}}}𝑑p^{}𝑑p𝑑k\gamma _{bc}^a(p^{};p,k)\delta (p^{}pk)`$
$`\times f_0^a(p^{})[1\pm f_0^b(p)][1\pm f_0^c(k)]\left[\chi ^a(p^{})\chi ^b(p)\chi ^c(k)\right]^2.`$ (23)
We have used the normalization condition (15) to simplify the expression (23). \[One may also use the identity (16) to contract the spatial indices in the $`22`$ piece (22).\] The overall coefficients of $`1/8`$ and $`1/2`$ in the collision parts (22) and (23) are symmetry factors which compensate for multiple counting of the same physical process. We have dropped the directional information $`\widehat{𝒑}`$ in the arguments of the collinear splitting functions $`\gamma _{bc}^a`$ because, in the linearized theory, these are to be evaluated in the background of the isotropic distribution $`f_0(𝒑)`$ and so do not depend on direction. See Appendix B for further details.

Maximizing the functional $`Q[\chi ]`$ is equivalent to computing an expectation value of the inverse of the linearized collision operator,

$$Q_{\mathrm{max}}=\frac{1}{2}(𝒮_{i\mathrm{}j},𝒞^1𝒮_{i\mathrm{}j}).$$
(24)
The value of the extremum directly determines the actual transport coefficients,

$`\sigma `$ $`=`$ $`\frac{2}{3}Q_{\mathrm{max}}|_{\mathrm{}=1,q=q_{\mathrm{EM}}},`$ (25a)
$`D_\alpha `$ $`=`$ $`\frac{2}{3}Q_{\mathrm{max}}|_{\mathrm{}=1,q=q_\alpha }\left({\displaystyle \frac{n_\alpha }{\mu _\alpha }}\right)^1,`$ (25b)
$`\eta `$ $`=`$ $`\frac{2}{15}Q_{\mathrm{max}}|_{\mathrm{}=2,q=|𝒑|}.`$ (25c)
The charge susceptibility appearing in expression (25b) for $`D_\alpha `$ is

$$\frac{n_\alpha }{\mu _\alpha }=\frac{T^2}{12}\underset{a}{}\nu _a\lambda _a(q_\alpha ^a)^2,$$
(26)
where $`\lambda _a`$ is 1 for fermionic and 2 for bosonic species. If $`n_\alpha `$ is the net number density of a single fermion flavor, then this susceptibility is $`\frac{1}{3}N_\mathrm{c}T^2`$ for a massless Dirac fermion in the fundamental representation of SU($`N_\mathrm{c}`$).

In numerical computations, we will focus in this paper on transport coefficients at zero chemical potential, so that the local equilibrium distributions $`f_0^a(p)`$ appearing in the pieces (20) and (21) of $`Q[\chi ]`$ are simply

$$f_0^a(p)=\frac{1}{e^{\beta p}1}.$$
(27)
We will primarily discuss the quark number diffusion constant $`D_\mathrm{q}`$ and the shear viscosity $`\eta `$ in QCD, although we will also report results for the electrical conductivity in a QED plasma (with or without quarks). Because we treat all relevant quarks as massless, the diffusion constants for individual flavors of quarks in QCD are all the same and equal to $`D_\mathrm{q}`$, which is also the diffusion constant for baryon number.

### II.4 Variational solution

To maximize the functional $`Q[\chi ]`$ exactly, one must work in the infinite dimensional space of arbitrary functions $`\chi (p)`$. However, as with many other variational problems, one can obtain highly accurate approximate results by performing a restricted extremization within a well chosen finite dimensional subspace. We will select a finite set of basis functions, $`\{\varphi ^{(m)}(p)\}`$, $`m=1\mathrm{}K`$, and only consider functions $`\chi ^a(p)`$ which are linear combinations of these basis functions,

$$\chi ^a(p)=\underset{m=1}{\overset{K}{}}\stackrel{~}{\chi }_m^a\varphi ^{(m)}(p).$$
(28)
(More precisely, we will choose increasingly large sets of basis functions, in order to study the convergence of results with the size of the basis.) Restricted to this subspace, the source and collision parts of the functional $`Q[\chi ]`$ become linear and quadratic combinations of the arbitrary coefficients $`\{\stackrel{~}{\chi }_m^a\}`$,

$`(\chi _{i\mathrm{}j},𝒮_{i\mathrm{}j})`$ $`=`$ $`{\displaystyle \underset{a,m}{}}\stackrel{~}{\chi }_m^a\stackrel{~}{S}_m^a,`$ (29)
$`(\chi _{i\mathrm{}j},𝒞\chi _{i\mathrm{}j})`$ $`=`$ $`{\displaystyle \underset{a,m}{}}{\displaystyle \underset{b,n}{}}\stackrel{~}{\chi }_m^a\stackrel{~}{C}_{mn}^{ab}\stackrel{~}{\chi }_n^b.`$ (30)
The matrix elements $`\stackrel{~}{C}_{mn}^{ab}`$ may be regarded as forming a square matrix $`\stackrel{~}{C}`$ of dimension $`N_sK`$, while the components $`\{\stackrel{~}{S}_m^a\}`$ and $`\{\stackrel{~}{\chi }_m^a\}`$ comprise vectors $`\stackrel{~}{S}`$ and $`\stackrel{~}{\chi }`$ of the same size.777As discussed in Ref. AMY1 , when computing transport coefficients in a plasma with vanishing chemical potentials, charge conjugation symmetry relates and constrains the departures from equilibrium for different species. For a QCD-like theory where all matter fields are the same type (e.g., fundamental representation fermions), the net effect is that there is only a single independent departure from equilibrium in the C-odd, $`\mathrm{}=1`$ channel relevant for flavor diffusion, while there are two independent functions (fermion and gauge boson) in the C-even, $`\mathrm{}=2`$ channel relevant for shear viscosity. Hence the actual size of the matrix $`\stackrel{~}{C}`$ one needs to deal with is just the basis size $`K`$ for diffusion, or $`2K`$ for shear viscosity.  The functional $`Q[\chi ]`$ restricted to this subspace is $`\stackrel{~}{Q}[\stackrel{~}{\chi }]=\stackrel{~}{\chi }^{}\stackrel{~}{S}\frac{1}{2}\stackrel{~}{\chi }^{}\stackrel{~}{C}\stackrel{~}{\chi }`$. The extremum of $`\stackrel{~}{Q}[\stackrel{~}{\chi }]`$ occurs at $`\stackrel{~}{\chi }=\stackrel{~}{C}^1\stackrel{~}{S}`$, and

$$\stackrel{~}{Q}_{\mathrm{max}}=\frac{1}{2}\stackrel{~}{\chi }^{}\stackrel{~}{S}=\frac{1}{2}\stackrel{~}{S}^{}\stackrel{~}{C}^1\stackrel{~}{S}.$$
(31)
Note that the variational estimate $`\stackrel{~}{Q}_{\mathrm{max}}`$ gives a lower bound on the true extremum $`Q_{\mathrm{max}}`$, and a nested sequence of variational estimates must converge monotonically upward to the true result.

The functional $`Q[\chi ]`$ is most sensitive to the behavior of $`\chi ^a(p)`$ for momenta near $`T`$. Having factored $`f_0^a[1\pm f_0^a]`$ out of the first order correction to local equilibrium, as done in Eq. (10), one may show888In the absence of $`12`$ processes, the asymptotic large and small $`p`$ behaviors are both proportional to $`p^{\mathrm{}}`$. Because of LPM suppressed $`12`$ processes, the large $`p`$ behavior is modified to $`p^{\mathrm{}{\scriptscriptstyle \frac{1}{2}}}`$. However, $`22`$ processes are much more efficient so this behavior only sets in at very large momenta which are irrelevant for transport coefficients. The soft region is also modified by $`12`$ processes, leading to $`p^\mathrm{}1`$ behavior. This only occurs for $`p<m_\mathrm{D}`$, and is an accident of our including finite Debye screening masses in $`22`$ processes but treating $`12`$ processes as perfectly collinear.  that the resulting functions $`\chi ^a(p)`$ grow no faster than $`p^{\mathrm{}}`$ as $`p\mathrm{}`$. Also, $`\chi ^a(p)`$ vanishes at least as fast as $`p^\mathrm{}1`$ at small $`p`$. Consequently, one reasonable choice of basis functions is

$$\varphi ^{(m)}(p)=\frac{p^\mathrm{}1(p/T)^{m1}}{(1+p/T)^{K2}},m=1,\mathrm{},K.$$
(32)
These functions are not orthogonal; they do not need to be. As noted below, choosing strictly positive basis functions improves the accuracy of numerical integrations. These basis functions span nested subspaces and, as $`K\mathrm{}`$, the basis set becomes complete.999For each $`K`$, our functions $`\varphi ^{(m)}`$ are linear combinations of the first $`K`$ functions in the sequence $`\chi _n(p)=p^\mathrm{}1(1+p/T)^{1n}`$, $`n=0,1,\mathrm{}`$. Therefore, as we increase $`K`$, the span of our functions strictly increases. One may easily show that the functions $`\{\chi _n\}`$ form a complete basis in $`L^2(_+,d\mu )`$ with measure $`d\mu =w(p)p^2dp`$ and $`w(p)`$ any weight function which is real, positive, smooth, bounded, and falls faster than $`p^{2l4}`$ as $`p\mathrm{}`$. The most natural weight function for our application is $`w(p)=f_0(p)[1\pm f_0(p)]`$.  that we expect the large $`K`$ limit of the finite $`K`$ variational extremum to give the extremal value in the full space of allowed $`\chi `$.

Evaluating matrix elements of the linearized $`22`$ collision operator (22) requires performing an eight dimensional integral (after accounting for energy and momentum conservation). By suitably choosing variables, as discussed in Appendix A, three of these integrals represent overall rotations in momentum space and are trivial, and one remaining angular integral may be performed analytically. This leaves a non-trivial four dimensional integral which must be evaluated numerically for each pair of basis functions. It is computationally challenging to perform these numerical integrations both efficiently and accurately. Use of a basis of strictly positive, non-orthogonal functions, such as the set (32), prevents cancellations between different regions of the integration which would otherwise degrade the accuracy of the numerically integration. Of course, one does not want to use a basis which causes the resulting matrix $`\stackrel{~}{C}`$ to be so ill-conditioned that its inversion in formula (31) becomes a problem. The choice of the functions (32) was motivated by the need to balance these two conflicting goals.

Due to the collinearity of the $`12`$ processes, computing matrix elements of this part of the linearized collision operator (23) involves only a two dimensional integral. As noted earlier, however, evaluating the integrand itself requires the solution of a two-dimensional linear integral equation. This evaluation is discussed in Appendix B; one can either convert the integral equation to a differential equation/boundary value problem which can be solved with a shooting method Gelis2 ; AGMZ , or one may again use a finite basis variational approach AMY3 .

### II.5 Thermal masses

In the next two sections we will present results for SU($`N`$) or U(1) gauge theories with $`N_\mathrm{f}`$ Dirac fermions in the fundamental representation.101010In a theory with chiral (Weyl) fermions, such as the standard electroweak theory (ignoring Yukawa interactions), each Weyl fermion contributes half as much as a Dirac fermion; so for instance, SU(2) theory with 12 left handed doublets behaves the same as SU(2) theory with 6 Dirac doublets.  Because the weak coupling behavior of transport coefficients in hot gauge theories is not given by a simple power series in $`g`$, leading order results cannot be presented just by reporting a coupling-independent value of the leading coefficient. Rather, leading order results contain non-trivial dependence on the ratio of the effective thermal mass for hard gauge bosons to the temperature, $`m_{\mathrm{eff},\mathrm{g}}/T`$, and the analogous ratio involving the effective thermal mass for hard fermions, $`m_{\mathrm{eff},\mathrm{f}}/T`$. These masses depend on the gauge coupling and on the particle content of the theory; at leading order,

$`m_{\mathrm{eff},\mathrm{g}}^2`$ $`=`$ $`\frac{1}{6}\left(C_\mathrm{A}+N_\mathrm{f}C_\mathrm{F}{\displaystyle \frac{d_\mathrm{F}}{d_\mathrm{A}}}\right)g^2T^2\frac{1}{2}m_\mathrm{D}^2,`$ (33)
$`m_{\mathrm{eff},\mathrm{f}}^2`$ $`=`$ $`\frac{1}{4}C_\mathrm{F}g^2T^22m_\mathrm{F}^2.`$ (34)
As indicated here, the thermal masses for hard excitations, $`m_{\mathrm{eff},\mathrm{g}}`$ and $`m_{\mathrm{eff},\mathrm{f}}`$, differ by factors of $`\sqrt{2}`$ from the more commonly used Debye mass $`m_\mathrm{D}`$ or thermal quark mass $`m_\mathrm{F}`$ (which is the thermal energy of a fermion with zero momentum). We will generally use the latter variables below. Here, $`d_\mathrm{F}`$ and $`d_\mathrm{A}`$ denote the dimensions of the fundamental and adjoint representations, respectively, while $`C_\mathrm{F}`$ and $`C_\mathrm{A}`$ are the corresponding quadratic Casimirs.111111For U(1), $`d_\mathrm{F}=d_\mathrm{A}=C_\mathrm{F}=1`$ and $`C_\mathrm{A}=0`$.
For SU(2), $`d_\mathrm{F}=C_\mathrm{A}=2`$, $`C_\mathrm{F}=3/4`$, and $`d_\mathrm{A}=3`$.
For SU(3), $`d_\mathrm{F}=C_\mathrm{A}=3`$, $`C_\mathrm{F}=4/3`$, and $`d_\mathrm{A}=8`$.
For SU($`N`$), $`d_\mathrm{F}=C_\mathrm{A}=N`$, $`C_\mathrm{F}=(N^21)/(2N)`$, and $`d_\mathrm{A}=N^21`$.  Note that the ratio of (leading-order) thermal masses, $`m_{\mathrm{eff},\mathrm{g}}/m_{\mathrm{eff},\mathrm{f}}`$, is independent of the gauge coupling and only depends on the particle content of the theory.

## III Leading order diffusion and shear viscosity in QCD

In Figure 1 we plot the fermion flavor diffusion constant $`D_\mathrm{q}`$ (multiplied by $`g^4T`$) as a function of $`m_\mathrm{D}/T`$ for QCD with several different values of $`N_\mathrm{f}`$. Figure 2 shows the shear viscosity, multiplied by $`g^4/T^3`$, as a function of $`m_\mathrm{D}/T`$ for the same set of theories. The results in these figures were computed using four basis functions of the form (32). Truncation error due to the use of a finite basis set is smaller than one part in $`10^4`$, or smaller than the width of the lines in Figs. 1 and 2. Errors dues to the finite numerical integration precision are also smaller than the widths of the lines. If only two basis functions are used, then errors of about 1% result.

### III.1 Quality of various simplifying approximations

To illustrate the relative importance of different parts of the collision term, and to examine the sensitivity of results to the correct inclusion of thermal self-energies, we plot in Figure 3 the relative change in the quark diffusion constant of three flavor QCD resulting from various (over)simplifications. The different levels of damage are as follows.

1. Neglect LPM corrections to the rate of $`12`$ processes.
2. Neglect collinear $`12`$ processes altogether.
3. Neglect $`12`$ processes and those $`22`$ terms which do not contribute at leading-log order. In other words, include only the underlined terms in the $`22`$ matrix elements shown in Table 3 in Appendix A.
4. Drop all terms which do not contribute at leading-log order, and then replace the correct momentum-dependent HTL self-energies in the exchanged propagator by just the Debye mass or thermal quark mass. That is, replace $`1/t^2`$ by $`1/(tm_\mathrm{D}^2)^2`$ in gauge boson exchange diagrams, and replace $`1/t`$ by $`t/(tm_\mathrm{F}^2)^2`$ in fermion exchange diagrams.

None of these approximations are correct at leading order in $`g`$, but some are far more numerically significant than others. As one may see from Fig. 3, neglecting the LPM effect makes quite a small change in the results, as does dropping $`22`$ terms which do not contribute at leading-log order.121212Both of these approximations actually overestimate the collision term. The LPM effect suppresses scatterings, so its neglect increases the collinear splitting rate. And, in SU(3) gauge theory, dropping the subdominant $`22`$ terms increases the collision rate because the most important neglected term is an interference term which is negative in a non-Abelian theory.  However, the collinear $`12`$ processes are important at the 10% level, while the approximation sometimes used in the literature of replacing the bosonic self-energy by just the Debye mass is quite poor and results in errors at the 50% level.

### III.2 Sensitivity to higher-order corrections

The leading-order results have relative corrections of order $`g`$, and hence are only reliable when $`m_\mathrm{D}/T`$ is sufficiently small. Of course, for a given level of precision, specifying just how small is “small enough” is not possible without knowing more about the actual size of sub-leading corrections.

It should be emphasized that our leading order results do depend on specific choices which were made in defining the linearized collision operator $`𝒞`$. We would obtain somewhat different numerical results if we had included thermal self-energies on internal lines of diagrams even when this was not strictly necessary, if we had used full one-loop self-energies instead of their HTL approximation (valid for momenta small compared to $`T`$), if we had included thermal corrections to the on-shell dispersion relations in $`22`$ processes, or if we had not approximated the nearly collinear $`12`$ processes as strictly collinear. All of these effects produce relative changes suppressed by at least one power of $`g`$, so we are justified in neglecting them. However, handling any of these issues differently can lead to different, but equally valid, formally leading-order results.

Examining the sensitivity of results to changes in the precise definition of the linearized collision operator $`𝒞`$ (which are equally valid at leading-order) is one way to get a “hint” as to the likely size of some actual $`O(g)`$ sub-leading effects. Fig. 4 shows results for the diffusion constant in three flavor SU(3) gauge theory, computed using several different, but equally valid at leading-order, definitions of the linearized collision operator.

We have only considered modifications of the way in which the thermal gauge field self-energy is introduced in $`t`$ or $`u`$ channel gauge boson exchange diagrams. We focus on these modifications because these diagrams numerically dominate the collision term, and because such modifications are easy to study. Our “standard” choice is discussed in Appendix A; it consists of writing $`t`$ (or $`u`$) channel exchange diagrams as the analogous result for scalar quarks, computed with the HTL self-energy included in the gauge boson propagator, plus a spin-dependent infrared safe remainder in which the HTL self-energy may be neglected. Our first alternative is to instead write $`t`$ or $`u`$ channel gauge boson exchanges as the result for fermion-from-fermion scatterings, plus a (different) IR safe piece. Expressions for fermion-from-fermion scattering (with the HTL gauge boson self-energy included) may be found in Appendix B of Ref. largeN .

The inclusion of hard thermal loop self-energies on $`t`$ or $`u`$ channel exchange lines is only justified for soft exchange momentum, since otherwise the HTL approximation to the full one-loop self-energy is invalid. For large exchange momentum, inserting the HTL self-energy is no more correct than neglecting the self-energy altogether, but in this regime self-energy corrections are $`O(g^2)`$ effects. So in a leading-order treatment, one valid approximation is to multiply the HTL self-energy by the step function $`\mathrm{\Theta }(T^2Q^2)`$, with $`Q`$ the exchanged 4-momentum.131313We use a $`(+++)`$ spacetime metric, so $`Q^2>0`$ for a spacelike momentum exchange.  Alternatively, one may multiply the HTL self-energy by $`\mathrm{\Theta }(T^2q^2)`$, with $`q=|𝒒|`$ the spatial 3-momentum transfer.

These various possibilities are compared in Fig. 4. The differences between these curves are all formally at most $`O(g)`$. The width of the band of results, for a particular value of $`m_\mathrm{D}/T`$, may be viewed as a guess as to the size of certain types of actual sub-leading corrections. The figure also includes the next-to-leading logarithmic approximation from the next section, for comparison. One sees that the different, but equally valid, leading-order implementations agree to within 15% provided $`m_\mathrm{D}/T0.8`$. This suggests that the leading-order results might have about as large a range of utility as one could have reasonably hoped for; certainly one should not expect a leading-order weak-coupling analysis to be accurate when $`m_\mathrm{D}T`$.

An obvious question is whether there is some natural way to define a unique leading-order result which would systematically drop all $`O(g)`$ corrections while retaining the full dependence on $`1/\mathrm{ln}g^1`$. This is directly related to the summability of the asymptotic expansion in $`1/\mathrm{ln}g^1`$, and will be discussed in Appendix C.

## IV Expansion in inverse logs

As discussed in Refs. BMPRa ; Heiselberg ; AMY1 , part of the linearized collision operator $`𝒞`$ contains a logarithmic enhancement proportional to $`\mathrm{ln}(T/m_\mathrm{D})\mathrm{ln}(1/g)`$ arising from $`t`$ (or $`u`$) channel $`22`$ exchange processes with momentum transfer between the Debye screening scale $`m_\mathrm{D}`$ and $`T`$. One may separate the full leading-order collision operator $`𝒞`$ into a logarithmically enhanced piece, which we will denote by $`𝒞_{\mathrm{LL}}`$, and a remainder $`\delta 𝒞`$,

$$𝒞=𝒞_{\mathrm{LL}}+\delta 𝒞.$$
(35)
Expanding the inverse collision operator appearing in $`Q_{\mathrm{max}}`$ \[c.f., Eq. (24)\] in powers of $`\delta 𝒞`$ will generate an asymptotic expansion of transport coefficients in powers of $`1/\mathrm{ln}g^1`$.

To make this separation precise, let $`\overline{𝒞}𝒞/(g^4T)`$ and note that this rescaled operator is dimensionless and depends on $`g`$ only through the ratios $`m_\mathrm{D}/T`$ and $`m_\mathrm{F}/T`$. Equivalently, since the (leading-order) ratio $`m_\mathrm{F}/m_\mathrm{D}`$ has a fixed value in a given theory, the $`g`$ dependence of $`\overline{𝒞}`$ may be regarded as arising only through dependence on $`m_\mathrm{D}/T`$. This dependence may be isolated by introducing a separation scale $`q_{}`$ satisfying $`m_\mathrm{D}q_{}T`$ and splitting the relevant $`t`$ and $`u`$ channel exchange parts of matrix elements of $`\overline{𝒞}`$ into contributions from exchange momentum less or greater than $`q_{}`$,

$$(\chi _{i\mathrm{}j},\overline{𝒞}\chi _{i\mathrm{}j})=(\chi _{i\mathrm{}j},\overline{𝒞}_{[q<q_{}]}\chi _{i\mathrm{}j})+(\chi _{i\mathrm{}j},\overline{𝒞}_{[q>q_{}]}\chi _{i\mathrm{}j}).$$
(36)
When $`q>q_{}`$, one may safely set $`m_\mathrm{D}=m_\mathrm{F}=0`$. This contribution is therefore independent of $`m_\mathrm{D}`$ and $`m_\mathrm{F}`$, up to corrections subleading in $`g`$. When $`q<q_{}`$, one may expand in $`q`$ and keep only the first nontrivial order. (This is safe because in all but a $`g^2`$ suppressed part of the integration domain, all other momenta are large compared to $`q_{}`$.) As far as the $`q`$ and $`m_\mathrm{D}`$ dependence is concerned, the result has the form Heiselberg ; AMY1

$$_0^q_{}\frac{dq}{q}f\left(\frac{m_\mathrm{D}}{q}\right).$$
(37)
Rescaling $`qq/m_\mathrm{D}`$, the integrand becomes $`g`$ independent and all remaining $`g`$ dependence is isolated in the upper limit of $`q_{}/m_\mathrm{D}`$. Simplifying the integrand at $`q=q_{}`$, using the assumed scale separation $`m_\mathrm{D}q_{}T`$, leads to an explicit form for

$$𝒜\underset{m_\mathrm{D}/T0}{lim}m_\mathrm{D}\frac{\overline{𝒞}(m_\mathrm{D}/T)}{m_\mathrm{D}},$$
(38)
given in Ref. AMY1 . For a gauge theory with $`N_\mathrm{f}`$ Dirac fermions all in the same representation,

$`(\chi _{i\mathrm{}j},𝒜\chi _{i\mathrm{}j})=`$ $`{\displaystyle \frac{N_\mathrm{f}d_\mathrm{F}C_\mathrm{F}^2\beta ^2}{32\pi ^3}}{\displaystyle _0^{\mathrm{}}}𝑑ppf_0^f(p)[1+f_0^g(p)]\left([\chi ^f(p)\chi ^g(p)]^2+[\chi ^{\overline{f}}(p)\chi ^g(p)]^2\right)`$
$`+`$ $`{\displaystyle \frac{d_\mathrm{A}C_\mathrm{A}\beta ^3m_\mathrm{D}^2}{(2\pi )^3g^2}}{\displaystyle _0^{\mathrm{}}}𝑑pf_0^g(p)[1+f_0^g(p)]\left(p^2[\chi ^g(p)^{}]^2+\mathrm{}(\mathrm{}+1)[\chi ^g(p)]^2\right)`$
$`+`$ $`{\displaystyle \frac{N_\mathrm{f}d_\mathrm{F}C_\mathrm{F}\beta ^3m_\mathrm{D}^2}{(2\pi )^3g^2}}{\displaystyle _0^{\mathrm{}}}dpf_0^f(p)[1f_0^f(p)](p^2[\chi ^f(p)^{}]^2+\mathrm{}(\mathrm{}+1)[\chi ^f(p)]^2`$
$`+p^2[\chi ^{\overline{f}}(p)^{}]^2+\mathrm{}(\mathrm{}+1)[\chi ^{\overline{f}}(p)]^2).`$ (39)
Here $`\chi ^g`$, $`\chi ^f`$, and $`\chi ^{\overline{f}}`$ are the departures from equilibrium for gauge bosons, fermions, and anti-fermions, respectively, and primes denote derivatives. The above form assumes that departures from equilibrium are fermion flavor independent, as is appropriate for computing shear viscosity or baryon number diffusion. The prefactor of $`\beta ^2m_\mathrm{D}^2/g^2`$ in two of the terms is just a compact way of writing $`\frac{1}{3}(C_\mathrm{A}+N_\mathrm{f}C_\mathrm{F}d_\mathrm{F}/d_\mathrm{A})`$. The first integral comes from fermion exchange diagrams whose infrared behavior is regulated by $`m_\mathrm{F}`$, while the latter two integrals arise from gauge boson exchange diagrams whose IR behavior is regulated by $`m_\mathrm{D}`$.

Since the operator $`𝒜`$ is itself $`m_\mathrm{D}`$ independent, the definition (38) implies that the limit

$$\delta \overline{𝒞}(\mu /T)\underset{m_\mathrm{D}0}{lim}\left[\overline{𝒞}(m_\mathrm{D}/T)𝒜\mathrm{ln}(\mu /m_\mathrm{D})\right]$$
(40)
exists. Therefore the original linearized collision operator, up to $`O(g)`$ corrections, may be written as

$`𝒞`$ $`=`$ $`𝒞_{\mathrm{LL}}(\mu )+\delta 𝒞(\mu ),`$ (41)
where
$`𝒞_{\mathrm{LL}}(\mu )`$ $``$ $`g^4T𝒜\mathrm{ln}(\mu /m_\mathrm{D}),`$ (42)
and
$`\delta 𝒞(\mu )`$ $``$ $`g^4T\delta \overline{𝒞}(\mu /T).`$ (43)
Here $`\mu `$ is an arbitrary scale which should be chosen to be $`O(T)`$ so that $`\delta 𝒞(\mu )`$ does not contain large logs.

Expanding the inverse collision operator in a geometric series,

$$𝒞^1=\left[𝒞_{\mathrm{LL}}(\mu )+\delta 𝒞(\mu )\right]^1\underset{n=0}{\overset{\mathrm{}}{}}(1)^n𝒞_{\mathrm{LL}}(\mu )^1\left[\delta 𝒞(\mu )𝒞_{\mathrm{LL}}(\mu )^1\right]^n,$$
(44)
directly gives its asymptotic expansion in powers of $`[\mathrm{ln}(\mu /m_\mathrm{D})]^11/\mathrm{ln}g^1`$. Inserting this expansion into $`Q_{\mathrm{max}}=\frac{1}{2}(𝒮_{i\mathrm{}j},𝒞^1𝒮_{i\mathrm{}j})`$ and using the definitions (42) and (43) yields the inverse log expansion of $`Q_{\mathrm{max}}`$ (and hence of transport coefficients),

$$Q_{\mathrm{max}}\frac{1}{g^4T}\underset{n=1}{\overset{\mathrm{}}{}}Q_n(\mu )\left(\mathrm{ln}\frac{\mu }{m_\mathrm{D}}\right)^n,$$
(45)
with

$$Q_n(\mu )\frac{1}{2}(1)^{n1}(𝒮_{i\mathrm{}j},𝒜^1\left[\delta \overline{𝒞}(\mu )𝒜^1\right]^{n1}𝒮_{i\mathrm{}j}).$$
(46)
Except for $`Q_1`$, the coefficients $`Q_n(\mu )`$ are $`\mu `$-dependent. However, the series is formally $`\mu `$ independent in exactly the same way that perturbative series in QCD are formally renormalization point independent even though individual $`n`$-loop contributions do depend on the renormalization point.

To avoid the presence of large logarithms in the coefficients $`Q_n(\mu )`$, one must choose $`\mu `$ to be $`O(T)`$ parametrically, but the exact coefficient is not uniquely prescribed. One somewhat natural choice is to select the value $`\mu _{}`$ for which $`Q_2(\mu _{})`$ vanishes. This can be termed the fastest apparent convergence (FAC) choice of $`\mu `$ at next-to-leading log order. Given results for $`Q_1`$ and $`Q_2(\mu )`$ at some other value of $`\mu `$,

$$\mu _{}=\mu \mathrm{exp}[Q_2(\mu )/Q_1].$$
(47)
To evaluate the coefficients $`Q_n(\mu )`$ of the expansion (45), we use exactly the same finite basis set approach described in section II.4. The linear operators $`𝒜`$ and $`\delta \overline{𝒞}`$ are replaced by their matrix representations $`\stackrel{~}{A}`$ and $`\stackrel{~}{\delta C}`$ in the finite basis set, and

$$\stackrel{~}{Q}_{\mathrm{max}}\frac{1}{g^4T}\underset{n=1}{\overset{\mathrm{}}{}}\stackrel{~}{Q}_n(\mu )\left(\mathrm{ln}\frac{\mu }{m_\mathrm{D}}\right)^n,$$
(48)
with

$$\stackrel{~}{Q}_n(\mu )=\frac{1}{2}(1)^{n1}\stackrel{~}{S}^{}\stackrel{~}{A}^1\left[\stackrel{~}{\delta 𝒞}(\mu )\stackrel{~}{A}^1\right]^{n1}\stackrel{~}{S}.$$
(49)
The limit (40) defining $`\delta \overline{𝒞}(\mu )`$ is performed numerically by evaluating each matrix element of $`\overline{𝒞}`$ at several small values of $`m_\mathrm{D}`$, subtracting off the leading log piece (whose matrix elements are easy to evaluate), and then extrapolating to vanishing Debye mass. We find that this extrapolation is quite well behaved, although the numerical integrals at small $`m_\mathrm{D}`$ or $`m_\mathrm{F}`$ become rather demanding.

Using the basis functions (32), we previously found AMY1  that the fractional difference between $`Q_1`$ and its finite basis approximation $`\stackrel{~}{Q}_1`$ is less than $`10^5`$ with 4 basis elements and less than $`10^6`$ using six. The higher $`Q_n`$ are more sensitive both to basis size and to numerical integration errors: $`Q_2`$ can be reliably determined with 4 basis functions, and $`Q_3`$ and $`Q_4`$ can be found with reasonably small errors using 6 to 8 basis functions, but higher moments become rapidly more difficult to evaluate, showing poor convergence with basis set size and high numerical integration error sensitivity when the basis sets become very large. Consequently we have been unable to go very deep in the $`Q_n`$ series.

Table 1 shows results for $`\mu _{}`$ and the first coefficient of the inverse log expansion, for the case of fermion flavor diffusion and shear viscosity in theories with various gauge groups and the indicated number $`N_\mathrm{f}`$ of Dirac fermion flavors in the fundamental representation. Specifically, we show the first coefficients $`D_1`$ and $`\eta _1`$ of the series141414 The diffusion constant $`D`$ and shear viscosity $`\eta `$ are related to their respective $`Q_{\mathrm{max}}`$’s as shown in Eq. (25). Explicitly, $`\eta =\frac{2}{15}Q_{\mathrm{max}}^{(\eta )}`$ and $`D=\frac{2}{3}Q_{\mathrm{max}}^{(D)}/(N_\mathrm{f}T^2)`$ with $`Q^{(D)}`$, in this normalization, being associated with total quark number. Hence, the inverse log expansion coefficients in Eqs. (50) and (51) are related to the previous coefficients $`Q_n(\mu )`$ via $`D_n(\mu )=\frac{2}{3}Q_n^{(D)}(\mu )/(N_\mathrm{f}T^2)`$ and $`\eta _n(\mu )=\frac{2}{15}Q_n^{(\eta )}(\mu )/T^4`$.

$$D\frac{1}{g^4T}\underset{n=1}{\overset{\mathrm{}}{}}D_n(\mu )\left(\mathrm{ln}\frac{\mu }{m_\mathrm{D}}\right)^n$$
(50)
and

$$\eta \frac{T^3}{g^4}\underset{n=1}{\overset{\mathrm{}}{}}\eta _n(\mu )\left(\mathrm{ln}\frac{\mu }{m_\mathrm{D}}\right)^n.$$
(51)
Together with the corresponding values of $`\mu _{}/T`$, these numbers determine the next-to-leading-log (NLL) approximation to the respective transport coefficients,

$`D_{\mathrm{NLL}}`$ $`=`$ $`{\displaystyle \frac{1}{g^4T}}\left[{\displaystyle \frac{D_1}{\mathrm{ln}(\mu _{}/m_\mathrm{D})}}\right]`$ (52)
and
$`\eta _{\mathrm{NLL}}`$ $`=`$ $`{\displaystyle \frac{T^3}{g^4}}\left[{\displaystyle \frac{\eta _1}{\mathrm{ln}(\mu _{}/m_\mathrm{D})}}\right].`$ (53)
Note that in SU(3) gauge theory, $`\mu _{}/T`$ is quite close to 3 for both transport coefficients, regardless of the number of fermion flavors.

We have computed further terms in the inverse log expansion in the case of SU(3) gauge theory with 3 fermion flavors. For flavor diffusion, we find

$$D_3(\mu _{})=2.436(2),D_4(\mu _{})=0.11(1),D_5(\mu _{})=1.7(1),$$
(54)
while for shear viscosity

$$\eta _3(\mu _{})=27(1),\eta _4(\mu _{})=6(5),\eta _5(\mu _{})=100(100).$$
(55)
\[And $`D_2(\mu _{})=0=\eta _2(\mu _{})`$, by our definition of $`\mu _{}`$.\] The third-order coefficients $`D_3(\mu _{})`$ and $`\eta _3(\mu _{})`$ are roughly one quarter the size of $`D_1`$ and $`\eta _1`$, respectively. The next order coefficients $`D_4(\mu _{})`$ and $`\eta _4(\mu _{})`$ are yet smaller, but subsequent coefficients appear to grow.

The behavior of the expansion in inverse powers of $`\mathrm{ln}(\mu _{}/m_\mathrm{D})`$, truncated at second or third order, is compared to the full leading-order result in Figure 5, for the case of flavor diffusion in three flavor QCD. As the figure makes clear, the next-to-leading-log result is remarkably close to the full leading order result out to $`m_\mathrm{D}/T=1`$, but going beyond second order in the inverse log expansion has very little practical utility.

Appendix C discusses the asymptotics of the inverse log expansion and proves that this expansion is, in fact, a typical asymptotic expansion with zero radius of convergence. An argument is also given suggesting that this asymptotic expansion is not Borel summable. As is well known, the presence of singularities in the Borel transform on the positive real axis generates ambiguities in the inverse Borel transform. In the typical case of a power series in $`g^2`$, this irreducible ambiguity is non-perturbative, behaving as $`\mathrm{exp}(c/g^2)`$ where $`c`$ is the location of the singularity of the Borel transform nearest to the origin on the positive real axis. (See, for example, Refs. Borel1 ; Borel2 .) In the present case where the expansion parameter is an inverse log of the coupling, a singularity in the Borel transform instead indicates an inherent ambiguity in the Borel sum of the asymptotic series which is a power of coupling. The estimate of appendix C suggests that this is an $`O(g^5)`$ ambiguity.

## V Electrical conductivity

We have also computed the electrical conductivity for a high temperature plasma of leptons, or leptons plus quarks. As mentioned earlier, the electrical conductivity is related to the diffusion constants of charged species via the Einstein relation (3). In a plasma of leptons plus quarks, we make the same $`\alpha _{\mathrm{EM}}^2\alpha _\mathrm{s}^2`$ approximation used in previous work BaymHeiselberg ; AMY1 : we neglect the electric current directly carried by quarks, and only compute the charged lepton diffusion constant. Because quarks undergo efficient QCD scattering (as compared to QED scatterings), their departure from equilibrium in the presence of an electric field is negligible compared to that of charged leptons, and hence so is their contribution to the electric current. This approximation amounts to the neglect of relative $`O(\alpha _{\mathrm{EM}}^2/\alpha _\mathrm{s}^2)`$ corrections to the conductivity. Quarks remain relevant, however, as excitations off of which charged leptons can scatter. Of course, for plasmas containing quarks we still require $`\alpha _\mathrm{s}1`$, so that quarks may be treated as nearly free massless excitations.

For simplicity we have only analyzed physically relevant combinations of leptons and quarks. We have evaluated the electrical conductivity $`\sigma `$ at next-to-leading-log order (NLL FAC) as well as at full leading order. In the next-to-leading-log form,

$$\sigma _{\mathrm{NLL}}=\frac{T}{e^2}\left[\frac{\sigma _1}{\mathrm{ln}(\mu _{}/m_\mathrm{D})}\right],$$
(56)
the Debye mass now refers to the inverse QED screening length given by

$$m_\mathrm{D}^2=\frac{1}{3}e^2T^2\left(_iq_i^2\right).$$
(57)
$`q_i`$ is the charge assignment of a given species, and the sum runs over all Dirac fermions. Instead of presenting plots showing the conductivity as a function of $`\alpha _{\mathrm{EM}}`$, we have simply set $`\alpha _{\mathrm{EM}}=1/137.04`$.151515 We use this many significant digits, ignoring the running of the coupling, for no reason other than to show the precision of the numerics and to compare different results.  Our results are presented in Table 2. For this quite small value of coupling, one sees that the NLL FAC treatment agrees with the full leading order result to better than 1%.

## VI Conclusion

We have performed complete leading-order calculations of shear viscosity, electrical conductivity, and fermion diffusivity in QCD and QED. “Leading-order” means that all neglected effects are suppressed by one or more powers of the gauge coupling $`g(T)`$. To our knowledge, this is the first time any transport coefficient has been evaluated with leading-order accuracy in a high temperature gauge theory. Due to the presence of Coulomb logarithms arising from small angle scattering, the coefficient of the leading power of $`g(T)`$ is not a simple number, as in scalar theories, but rather is a non-trivial function of $`\mathrm{ln}(g^1)`$.

Leading-order results for transport coefficients may themselves be expanded in powers of $`1/\mathrm{ln}(g^1)`$. We have explicitly computed both leading and first sub-leading terms for shear viscosity and quark diffusivity in U(1), SU(2), and SU(3) gauge theories with various numbers of fermion fields (as well as several more terms for three flavor QCD). For QCD, the next-to-leading log result (with the sub-leading term absorbed by suitably shifting the scale inside the leading log) was found to be remarkably close to the full leading-order result as long as $`m_\mathrm{D}/T1`$. This is a much larger domain of utility than one might have expected. For these transport coefficients, we also find that only roughly 10% errors are made if one neglects near-collinear $`12`$ particle splitting processes, which are considerably more difficult to analyze than $`22`$ particle scattering processes. (However, it should be noted that some transport coefficients which we have not analyzed, such as bulk viscosity, depend primarily on particle number-changing processes and so may be expected to depend essentially on $`12`$ processes.)

Because the expansion in inverse powers of $`\mathrm{ln}(g^1)`$ is only asymptotic, not convergent, as demonstrated in Appendix C, we are not able to give a unique, unambiguous prescription for separating leading-order contributions from higher-order effects. As discussed in Appendix C, it appears that the inverse log expansion is not Borel summable, which would imply that no clean separation is possible. In practice, this means that any specific calculation yielding results of leading-order accuracy will necessarily include some contributions from higher-order effects. However, our examination of several different prescriptions for computing leading-order results suggests that this is not a significant issue for $`m_\mathrm{D}\text{ }<\text{ }\text{ }0.8T`$.

Our tool for studying transport coefficients has been kinetic theory, specifically the effective kinetic theory presented in our previous paper AMY5 . As originally shown by Jeon Jeon , in the context of weakly-coupled relativistic scalar theories, it is also possible to compute transport coefficients diagrammatically starting from the appropriate Kubo formulae involving current-current or stress-stress correlators. Such a diagrammatic approach amounts to a complicated way to derive the appropriate linearized Boltzmann equation specialized to the particular symmetry channel of interest. For gauge theories, this diagrammatic approach has been applied, only at leading logarithmic order, to the electrical conductivity by Valle Basagoiti Basagoiti . Trying to use a diagrammatic approach for a complete leading order calculation would be an enormously more difficult task. For instance, a typical diagram which we believe contributes at leading order to the $`T_{ij}T_{kl}`$ correlator, needed for the shear viscosity, is depicted in Figure 6. Note that the hard (nearly) on-shell propagators require self-energy resummations, illustrated at the bottom of the figure, which go far beyond the HTL approximation. The one and two loop self-energy contributions shown account for scattering via $`22`$ processes, while the very complicated self-energy diagram represents one contribution to the effective $`12`$ splitting process; the reason it needs so many loops is that this process can involve any number of soft scatterings off of other particles in the plasma \[the number of such scatterings is summed over by Eq. (88)\]. The complicated “cross-rungs” in the upper diagram are the result of opening up one of the lines in any one of the self-energy contributions. For more discussion of this point, see Ref. GDMSEWM02 .

An interesting problem for the future is to understand the accuracy of leading-order calculations of transport coefficients by calculating higher order effects explicitly. In the case of the QCD free energy it is known that, beyond the ideal gas result, the perturbative expansion in powers of $`g(T)`$ is quite poorly behaved hot-free1 ; hot-free2 ; hot-free3 , except for unrealistically small values of $`g(T)`$. (Specifically, $`T`$ must substantially exceed the Planck scale.) Does this same unpleasant behavior apply to transport coefficients? At the moment, the only known test case is a many flavor limit of QCD largeN , where the leading-order result (as well as the next-to-leading-log approximation thereto) is quite successful — its accuracy is comparable to the renormalization point sensitivity. It would be useful to know if this holds more generally.

###### Acknowledgements.

This work was supported, in part, by the U.S. Department of Energy under Grant Nos. DE-FG03-96ER40956 and DE-FG02-97ER41027.

## Appendix A $`\mathrm{𝟐}\mathbf{}\mathrm{𝟐}`$ Matrix Elements

The matrix elements for all $`22`$ particle processes in a QCD-like theory, neglecting thermal self-energy corrections, are listed in Table 3. These matrix elements arise from the diagrams shown in Fig. 7. Terms in Table 3 with underlined denominators are sufficiently infrared sensitive that thermal self-energy corrections must be included, as discussed in Ref. AMY5 . Singly-underlined denominators indicate IR sensitive contributions arising from soft gauge boson exchange, while double-underlined denominators indicate IR sensitive contributions from a soft exchanged fermion.

### A.1 Self-energy corrections to matrix elements

#### Fermion self-energy

When $`22`$ particle processes involving $`t`$ (or $`u`$) channel fermion exchange are computed using free propagators, the resulting squared matrix elements (shown with double-underlined denominators in Table 3) generate logarithmic infrared divergences in the collision term (22). This logarithmic infrared sensitivity is cut off by the inclusion of the retarded thermal self-energy $`\mathrm{\Sigma }(Q)`$, so that the internal fermion propagator appearing inside the matrix element $``$ is $`[\text{ / }Q\text{ / }\mathrm{\Sigma }(Q)]^1`$. Its conjugate, the advanced self energy $`\mathrm{\Sigma }^{}(Q)`$, appears in $`^{}`$. Since the exchange four-momentum $`Q^\mu (\omega ,𝒒)`$ is spacelike, and the thermal self-energy is only relevant when $`Q`$ is soft, we only need the self-energy for spacelike momenta, $`q|𝒒|>|\omega |`$, in the hard thermal loop limit, $`qT`$. In this regime the self-energy was originally evaluated by Klimov Klimov3  (and independently by Weldon Weldon2 ), and is given by

$`\mathrm{\Sigma }^0(Q)`$ $`=`$ $`{\displaystyle \frac{m_\mathrm{F}^2}{2q}}\left[\mathrm{ln}\left({\displaystyle \frac{q+\omega }{q\omega }}\right)i\pi \right],`$ (58)
$`𝚺(Q)`$ $`=`$ $`\widehat{𝒒}{\displaystyle \frac{m_\mathrm{F}^2}{q}}\left\{1{\displaystyle \frac{\omega }{2q}}\left[\mathrm{ln}\left({\displaystyle \frac{q+\omega }{q\omega }}\right)i\pi \right]\right\}.`$ (59)
Here $`m_\mathrm{F}=\sqrt{C_\mathrm{F}/8}gT`$ is the (leading-order) “fermion thermal mass,” equal to the thermal energy of a fermion at zero momentum.

For the process $`fgfg`$, the net effect of the inclusion of the fermion self-energy is to make the replacement

$`{\displaystyle \frac{s}{u}}`$ $``$ $`{\displaystyle \frac{\frac{1}{8}\mathrm{tr}\left[\text{ / }P\gamma ^\mu \text{ / }𝒬\gamma ^\nu \text{ / }P^{}\gamma _\nu \text{ / }𝒬^{}\gamma _\mu \right]}{\left|𝒬𝒬\right|^2}}={\displaystyle \frac{4\mathrm{Re}[(P𝒬)(P^{}𝒬^{})]+t𝒬𝒬^{}}{\left|𝒬𝒬\right|^2}},`$ (60)
\[with $`(+++)`$ metric convention\] in the contribution shown in Table 3. Here $`P^\mu =(|𝒑|,𝒑)`$ and $`P^\mu =(|𝒑^{}|,𝒑^{})`$ are the incoming and outgoing fermion 4-momenta, respectively, and $`𝒬^\mu P^\mu K^\mu \mathrm{\Sigma }^\mu (PK^{})`$ with $`K^{}`$ the outgoing gauge boson 4-momentum.

The analogous replacements needed in the $`ffgg`$ squared matrix elements are

$`{\displaystyle \frac{u}{t}}`$ $``$ $`{\displaystyle \frac{4\mathrm{Re}[(P𝒬)(K𝒬^{})]+s𝒬𝒬^{}}{\left|𝒬𝒬\right|^2}}|_{𝒬^\mu =P^\mu P^\mu \mathrm{\Sigma }^\mu (PP^{})},`$ (61)
and
$`{\displaystyle \frac{t}{u}}`$ $``$ $`{\displaystyle \frac{4\mathrm{Re}[(P𝒬)(K𝒬^{})]+s𝒬𝒬^{}}{\left|𝒬𝒬\right|^2}}|_{𝒬^\mu =P^\mu K^\mu \mathrm{\Sigma }^\mu (PK^{})},`$ (62)
where $`P`$ and $`K`$ are the incoming fermion momenta, and $`P^{}`$ and $`K^{}`$ the outgoing gauge boson momenta.

#### Gauge boson self-energy

Processes involving $`t`$ or $`u`$ channel gauge boson exchange require inclusion of the thermal gauge boson self-energy on the internal propagator to cut off the infrared sensitivity of these processes. Because the self-energy only matters when the exchange momentum is soft, one may exploit the fact that soft gluon exchange between hard particles is spin-independent (to leading order) AMY4 . If one separates the IR sensitive matrix elements (those with singly-underlined denominators in Table 3) into (i) the result one would have with fictitious scalar quarks plus (ii) a spin-dependent remainder, then all the IR sensitivity resides in the first spin-independent piece. This is the only piece which must be recomputed with the thermal self-energy included. For the $`t`$-channel exchange terms, this amounts to using the exact identities

$$\frac{s^2+u^2}{t^2}=\frac{1}{2}+\frac{1}{2}\frac{(su)^2}{t^2},\frac{su}{t^2}=\frac{1}{4}\frac{1}{4}\frac{(su)^2}{t^2},$$
(63)
and then replacing

$$\frac{(su)^2}{t^2}\left|D(PP^{})_{\mu \nu }(P+P^{})^\mu (K+K^{})^\nu \right|^2,$$
(64)
where $`D(Q)_{\mu \nu }`$ is the retarded thermal equilibrium gauge field propagator, evaluated in the HTL approximation.

The HTL result for the above replacement (64) does not depend on gauge choice. One convenient choice is Coulomb gauge, where BraatenPisarski

$`D_{00}(\omega ,𝒒)`$ $`=`$ $`{\displaystyle \frac{1}{q^2+\mathrm{\Pi }_{00}(\omega ,q)}},`$ (65)
$`D_{ij}(\omega ,𝒒)`$ $`=`$ $`{\displaystyle \frac{\delta _{ij}\widehat{𝒒}_i\widehat{𝒒}_j}{q^2\omega ^2+\mathrm{\Pi }_\mathrm{T}(\omega ,q)}},`$ (66)
$`D_{0i}(\omega ,𝒒)`$ $`=`$ $`D_{i0}(\omega ,𝒒)=\mathrm{\hspace{0.25em}0}.`$ (67)
The equilibrium transverse and longitudinal gauge boson self-energies are Klimov1 ; Weldon1

$`\mathrm{\Pi }_\mathrm{T}(\omega ,𝒒)`$ $`=`$ $`m_\mathrm{D}^2\left\{{\displaystyle \frac{\omega ^2}{2q^2}}+{\displaystyle \frac{\omega (q^2\omega ^2)}{4q^3}}\left[\mathrm{ln}\left({\displaystyle \frac{q+\omega }{q\omega }}\right)i\pi \right]\right\},`$ (68)
$`\mathrm{\Pi }_{00}(\omega ,𝒒)`$ $`=`$ $`m_\mathrm{D}^2\left\{1{\displaystyle \frac{\omega }{2q}}\left[\mathrm{ln}\left({\displaystyle \frac{q+\omega }{q\omega }}\right)i\pi \right]\right\},`$ (69)
where we have assumed $`|\omega |<q`$, which is the only case of relevance.

### A.2 Integration variables

Since all external particles are to be treated as massless, the domain of the phase space integrations appearing in matrix elements of the linearized collision operator (22) are the same for all $`22`$ processes. One rather straightforward method for doing the multi-dimensional numerical integration is to use an adaptive Monte Carlo integrator. This can give reasonable accuracy at a tolerable investment of computational effort, but for the highest accuracy it is preferable to use nested one-dimensional adaptive Gaussian integration.

In order to handle efficiently the infrared-sensitive terms in $`t`$ and $`u`$ channel processes which give rise to leading-log contributions, it is useful to pick the exchange momentum $`q`$ and energy $`\omega `$ as two of the integration variables. When doing nested adaptive quadrature integrations, it is especially advantageous to choose integration variables in a manner which allows one to perform analytically as many of the integrations as possible. In particular, it is convenient to use different parameterizations for terms in Table 3 having different denominators. For terms having denominators of $`t`$ (or $`t^2`$), the $`t`$ channel parameterization described below allows all but four integrations to be done analytically. And similarly, the $`u`$ and $`s`$ channel parameterizations described below allow all but four integrations to be done analytically for terms with denominators of $`u`$ (or $`u^2`$), or $`s`$ (or $`s^2`$), respectively. The constant term in the $`gggg`$ matrix element can be handled using any of these parameterizations. The only other terms in Table 3 are those involving $`s^2/tu`$ or $`u^2/st`$. These can be reduced to the previous cases by rewriting

$`{\displaystyle \frac{s^2}{tu}}`$ $`=`$ $`{\displaystyle \frac{s}{t}}{\displaystyle \frac{s}{u}},{\displaystyle \frac{u^2}{st}}={\displaystyle \frac{u}{s}}{\displaystyle \frac{u}{t}}.`$ (70)
In what follows, our convention for labeling momenta in $`22`$ processes is that $`P,KP^{},K^{}`$.

#### $`t`$ channel parameterization

For terms containing $`t=(P^{}P)^2`$ in the denominator, it is convenient to use the spatial delta function in (22) to perform the $`𝒌^{}`$ integration, and to shift the $`𝒑^{}`$ integration into an integration over $`𝒑^{}𝒑𝒒`$. The angular integrals may be written in spherical coordinates defined such that the $`z`$ axis is in the direction of $`𝒒`$ while $`𝒑`$ lies in the $`xz`$ plane. This yields

$`(\chi _{i\mathrm{}j},𝒞^{22}\chi _{i\mathrm{}j})`$ $`=`$ $`{\displaystyle \frac{\beta ^3}{(4\pi )^6}}{\displaystyle \underset{abcd}{}}{\displaystyle _0^{\mathrm{}}}q^2𝑑qp^2𝑑pk^2𝑑k{\displaystyle _1^1}d\mathrm{cos}\theta _{pq}d\mathrm{cos}\theta _{kq}{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle \frac{1}{pkp^{}k^{}}}`$ (71)
$`\times \left|_{cd}^{ab}\right|^2\delta (p+kp^{}k^{})f_0^a(p)f_0^b(k)[1\pm f_0^c(p^{})][1\pm f_0^d(k^{})]`$
$`\times \left[\chi _{i\mathrm{}j}^a(𝒑)+\chi _{i\mathrm{}j}^b(𝒌)\chi _{i\mathrm{}j}^c(𝒑^{})\chi _{i\mathrm{}j}^d(𝒌^{})\right]^2,`$
where $`p`$, $`k`$, and $`q`$ denote to the magnitudes of the corresponding three-momenta (not the associated 4-momenta), $`p^{}|𝒒+𝒑|`$ and $`k^{}|𝒌𝒒|`$ are the magnitudes of the outgoing momenta, $`\varphi `$ is the azimuthal angle of $`𝒌`$ (and $`𝒌^{}`$) \[i.e., the angle between the $`𝒑`$,$`𝒒`$ plane and the $`𝒌`$,$`𝒒`$ plane\], and $`\theta _{pq}`$ is the angle between $`𝒑`$ and $`𝒒`$ (so $`\mathrm{cos}\theta _{pq}\widehat{𝒑}\widehat{𝒒}`$), etc.

Following Baym et al. BMPRa , it is convenient to introduce a dummy integration variable $`\omega `$, defined to equal the energy transfer $`p^{}p`$, so that

$$\delta (p+kp^{}k^{})=_{\mathrm{}}^{\mathrm{}}𝑑\omega \delta (\omega +pp^{})\delta (\omega k+k^{}).$$
(72)
Evaluating $`p^{}=|𝒑+𝒒|`$ in terms of $`p`$, $`q`$, and $`\mathrm{cos}\theta _{pq}`$, and defining $`t=\omega ^2q^2`$ (which is the usual Mandelstam variable), one finds

$`\delta (\omega +pp^{})`$ $`=`$ $`{\displaystyle \frac{p^{}}{pq}}\delta \left(\mathrm{cos}\theta _{pq}{\displaystyle \frac{\omega }{q}}{\displaystyle \frac{t}{2pq}}\right)\mathrm{\Theta }(\omega +p),`$ (73)
$`\delta (\omega k+k^{})`$ $`=`$ $`{\displaystyle \frac{k^{}}{kq}}\delta \left(\mathrm{cos}\theta _{kq}{\displaystyle \frac{\omega }{q}}+{\displaystyle \frac{t}{2kq}}\right)\mathrm{\Theta }(k\omega ),`$ (74)
with $`\mathrm{\Theta }`$ the unit step function. The $`\mathrm{cos}\theta `$ integrals may now be trivially performed and yield one provided $`p>\frac{1}{2}(q\omega )`$, $`k>\frac{1}{2}(q+\omega )`$, and $`|\omega |<q`$; otherwise the delta functions cannot both be satisfied. The remaining integrals are

$`(\chi _{i\mathrm{}j},𝒞^{22}\chi _{i\mathrm{}j})`$ $`=`$ $`{\displaystyle \frac{\beta ^3}{(4\pi )^6}}{\displaystyle \underset{abcd}{}}{\displaystyle _0^{\mathrm{}}}𝑑q{\displaystyle _q^q}𝑑\omega {\displaystyle _{\frac{q\omega }{2}}^{\mathrm{}}}𝑑p{\displaystyle _{\frac{q+\omega }{2}}^{\mathrm{}}}𝑑k{\displaystyle _0^{2\pi }}𝑑\varphi `$ (75)
$`\times \left|_{cd}^{ab}\right|^2f_0^a(p)f_0^b(k)[1\pm f_0^c(p^{})][1\pm f_0^d(k^{})]`$
$`\times \left[\chi _{i\mathrm{}j}^a(𝒑)+\chi _{i\mathrm{}j}^b(𝒌)\chi _{i\mathrm{}j}^c(𝒑^{})\chi _{i\mathrm{}j}^d(𝒌^{})\right]^2,`$
with $`p^{}=p+\omega `$ and $`k^{}=k\omega `$. For evaluating the final factor of (75), we use the relationship (16) that

$$I_{i\mathrm{}j}(\widehat{𝒑})I_{i\mathrm{}j}(\widehat{𝒌})=P_{\mathrm{}}(\mathrm{cos}\theta _{pk}).$$
(76)
One therefore needs expressions for the angles between all species, as well as the remaining Mandelstam variables $`s`$ and $`u`$, which may appear in $`^2`$. They are

$`s`$ $`=`$ $`{\displaystyle \frac{t}{2q^2}}\left\{\left[(p+p^{})(k+k^{})+q^2\right]\mathrm{cos}\varphi \sqrt{\left(4pp^{}+t\right)\left(4kk^{}+t\right)}\right\},`$ (77a)
$`u`$ $`=`$ $`ts,`$ (77b)
and

$`\mathrm{cos}\theta _{pq}`$ $`=`$ $`{\displaystyle \frac{\omega }{q}}+{\displaystyle \frac{t}{2pq}},\mathrm{cos}\theta _{p^{}q}={\displaystyle \frac{\omega }{q}}{\displaystyle \frac{t}{2p^{}q}},`$ (78a)
$`\mathrm{cos}\theta _{kq}`$ $`=`$ $`{\displaystyle \frac{\omega }{q}}{\displaystyle \frac{t}{2kq}},\mathrm{cos}\theta _{k^{}q}={\displaystyle \frac{\omega }{q}}+{\displaystyle \frac{t}{2k^{}q}},`$ (78b)
$`\mathrm{cos}\theta _{pp^{}}`$ $`=`$ $`1+{\displaystyle \frac{t}{2pp^{}}},\mathrm{cos}\theta _{kk^{}}=1+{\displaystyle \frac{t}{2kk^{}}},`$ (78c)
$`\mathrm{cos}\theta _{pk^{}}`$ $`=`$ $`1+{\displaystyle \frac{u}{2pk^{}}},\mathrm{cos}\theta _{p^{}k}=1+{\displaystyle \frac{u}{2p^{}k}},`$ (78d)
$`\mathrm{cos}\theta _{pk}`$ $`=`$ $`1{\displaystyle \frac{s}{2pk}},\mathrm{cos}\theta _{p^{}k^{}}=1{\displaystyle \frac{s}{2p^{}k^{}}}.`$ (78e)
For terms in which only $`t`$ appears in the denominator of the matrix element, the $`\varphi `$ integration can be easily performed analytically, leaving four integrals which must be evaluated numerically.

#### $`u`$ channel parameterization

For terms in which $`u=(K^{}P)^2`$ appears in the denominator, exchanging $`𝒑^{}`$ and $`𝒌^{}`$ in the $`t`$ channel parameterization provides the natural choice of variables.

#### $`s`$ channel parameterization

For terms in which $`s=(P+K)^2`$ appears in the denominator, one may use the spatial delta function in Eq. (22) to perform the $`𝒌^{}`$ integration, and then shift the $`𝒌`$ integration to an integral over $`𝒒=𝒑+𝒌`$, the total incoming spatial momentum (and the momentum on the internal propagator in $`s`$ channel exchange processes). Again choosing spherical coordinates so that $`𝒒`$ lies on the $`z`$ axis and $`𝒑`$ lies in the $`xz`$ plane, the $`22`$ contribution (22) becomes

$`(\chi _{i\mathrm{}j},𝒞^{22}\chi _{i\mathrm{}j})`$ $`=`$ $`{\displaystyle \frac{\beta ^3}{(4\pi )^6}}{\displaystyle \underset{abcd}{}}{\displaystyle _0^{\mathrm{}}}q^2𝑑qp^2𝑑pp_{}^{}{}_{}{}^{2}𝑑p^{}{\displaystyle _1^1}d\mathrm{cos}\theta _{pq}d\mathrm{cos}\theta _{p^{}q}{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle \frac{1}{pkp^{}k^{}}}`$ (79)
$`\times \left|_{cd}^{ab}\right|^2\delta (p+kp^{}k^{})f_0^a(p)f_0^b(k)[1\pm f_0^c(p^{})][1\pm f_0^d(k^{})]`$
$`\times \left[\chi _{i\mathrm{}j}^a(𝒑)+\chi _{i\mathrm{}j}^b(𝒌)\chi _{i\mathrm{}j}^c(𝒑^{})\chi _{i\mathrm{}j}^d(𝒌^{})\right]^2,`$
where now $`k=|𝒒𝒑|`$, $`k^{}=|𝒒𝒑^{}|`$, and $`\varphi `$ is the azimuthal angle of $`𝒌`$ (and $`𝒌^{}`$). Introducing the total energy $`\omega `$ via

$$\delta (p+kp^{}k^{})=_0^{\mathrm{}}𝑑\omega \delta (\omega pk)\delta (\omega p^{}k^{}),$$
(80)
and defining $`s=\omega ^2q^2`$ (which is the usual Mandelstam variable), one finds

$`\delta (\omega pk)`$ $`=`$ $`{\displaystyle \frac{k}{pq}}\delta \left(\mathrm{cos}\theta _{pq}{\displaystyle \frac{\omega }{q}}+{\displaystyle \frac{s}{2pq}}\right)\mathrm{\Theta }(\omega p),`$ (81)
$`\delta (\omega p^{}k^{})`$ $`=`$ $`{\displaystyle \frac{k^{}}{p^{}q}}\delta \left(\mathrm{cos}\theta _{p^{}q}{\displaystyle \frac{\omega }{q}}+{\displaystyle \frac{s}{2p^{}q}}\right)\mathrm{\Theta }(\omega p^{}).`$ (82)
Integration over $`\mathrm{cos}\theta _{pq}`$ and $`\mathrm{cos}\theta _{p^{}q}`$ yields unity provided $`q<\omega `$, $`|2p\omega |<q`$, and $`|2p^{}\omega |<q`$ (and zero otherwise). Therefore,

$`(\chi _{i\mathrm{}j},𝒞^{22}\chi _{i\mathrm{}j})`$ $`=`$ $`{\displaystyle \frac{\beta ^3}{(4\pi )^6}}{\displaystyle \underset{abcd}{}}{\displaystyle _0^{\mathrm{}}}𝑑\omega {\displaystyle _0^\omega }𝑑q{\displaystyle _{\frac{\omega q}{2}}^{\frac{\omega +q}{2}}}𝑑p{\displaystyle _{\frac{\omega q}{2}}^{\frac{\omega +q}{2}}}𝑑p^{}{\displaystyle _0^{2\pi }}𝑑\varphi `$ (83)
$`\times \left|_{cd}^{ab}\right|^2f_0^a(p)f_0^b(k)[1\pm f_0^c(p^{})][1\pm f_0^d(k^{})]`$
$`\times \left[\chi _{i\mathrm{}j}^a(𝒑)+\chi _{i\mathrm{}j}^b(𝒌)\chi _{i\mathrm{}j}^c(𝒑^{})\chi _{i\mathrm{}j}^d(𝒌^{})\right]^2,`$
with $`k=\omega p`$ and $`k^{}=\omega p^{}`$. The other Mandelstam variables are

$`t`$ $`=`$ $`{\displaystyle \frac{s}{2q^2}}\left\{\left[(pk)(p^{}k^{})q^2\right]+\mathrm{cos}\varphi \sqrt{(4pks)(4p^{}k^{}s)}\right\},`$ (84a)
$`u`$ $`=`$ $`st,`$ (84b)
and the angles between $`𝒒`$ and the external momenta are

$`\mathrm{cos}\theta _{pq}`$ $`=`$ $`{\displaystyle \frac{\omega }{q}}{\displaystyle \frac{s}{2pq}},\mathrm{cos}\theta _{p^{}q}={\displaystyle \frac{\omega }{q}}{\displaystyle \frac{s}{2p^{}q}},`$ (85a)
$`\mathrm{cos}\theta _{kq}`$ $`=`$ $`{\displaystyle \frac{\omega }{q}}{\displaystyle \frac{s}{2kq}},\mathrm{cos}\theta _{k^{}q}={\displaystyle \frac{\omega }{q}}{\displaystyle \frac{s}{2k^{}q}}.`$ (85b)
Eqs. (78c)–(78e) for the angles between external momenta still hold. For terms in which only $`s`$ appears in a denominator, the $`\varphi `$ integration can be easily performed analytically, leaving four integrals to do numerically.

## Appendix B Matrix elements for $`12`$ processes

In this appendix, we will review the integral equations that determine the splitting/joining functions $`\gamma _{bc}^a`$ for effective $`12`$ processes, appearing in our various formulations (8) and (23) of the Boltzmann collision term. These equations are summarized in Ref. AMY5 , but here we will make a number of simplifications applicable to the problem at hand.

The only $`12`$ processes with $`O(g^4T^4)`$ total rates per unit volume involve hard, collinear creation or destruction of a gauge boson. In QCD, the relevant processes are $`qqg`$, $`\overline{q}\overline{q}g`$, $`gq\overline{q}`$, and $`ggg`$. By “hard” we mean that the outgoing states each have $`O(T)`$ energy, and by “collinear” we mean that the angles between the three external momentum vectors are all $`O(g)`$. \[These parametrically small opening angles are ignored when evaluating distribution functions in the $`12`$ piece of our collision term (8), which only makes a sub-leading error in the evaluation of transport coefficients.\]

The duration of such processes (also known as the formation time of the gauge boson) is parametrically $`O(1/g^2T)`$, which is the same as the mean free path for hard particles to undergo soft scattering (with momentum exchange of order $`gT`$). For this reason, one must sum all possible number of soft scatterings with other excitations during the emission process. (In realistic theories, at least one such soft scattering is required to allow for energy-momentum conservation in the $`12`$ process.) This summation can be implemented by an integral equation which must be solved numerically. For a complete discussion, and a derivation of the initial integral equation presented below, see Refs. AMY5 ; AMY2 ; AMY4 .

### B.1 The integral equation for $`𝜸_{𝒃𝒄}^𝒂`$

For non-Abelian gauge theories such as QCD, the splitting/joining functions $`\gamma _{bc}^a(p^{};p,k)`$ for particle types $`abc`$ with momenta $`p^{}pk`$ are given in equilibrium by AMY5

$`\gamma _{qg}^q(p^{};p,k)`$ $`=`$ $`\gamma _{\overline{q}g}^{\overline{q}}(p^{};p,k)={\displaystyle \frac{p_{}^{}{}_{}{}^{2}+p^2}{p_{}^{}{}_{}{}^{2}p^2k^3}}_\mathrm{q}(p^{},p,k),`$ (86a)
$`\gamma _{q\overline{q}}^g(p^{};p,k)`$ $`=`$ $`{\displaystyle \frac{p^2+k^2}{k^2p^2p_{}^{}{}_{}{}^{3}}}_\mathrm{q}(k,p,p^{}),`$ (86b)
$`\gamma _{gg}^g(p^{};p,k)`$ $`=`$ $`{\displaystyle \frac{p_{}^{}{}_{}{}^{4}+p^4+k^4}{p_{}^{}{}_{}{}^{3}p^3k^3}}_\mathrm{g}(p^{},p,k),`$ (86c)
where

$$_s(p^{},p,k)\frac{d_sC_s\alpha }{2(2\pi )^3}\frac{d^2h}{(2\pi )^2}\mathrm{\hspace{0.33em}2}𝒉\mathrm{Re}𝑭_s(𝒉;p^{},p,k)$$
(87)
and $`\alpha g^2/(4\pi )`$. The function $`𝑭_s(𝒉;p^{},p,k)`$, for fixed given values of $`p^{}`$, $`p`$, and $`k`$, depends on a two-dimensional vector $`𝒉`$ which is related to transverse momentum during the splitting process. $`𝑭_s`$ is the solution to the linear integral equation

$`2𝒉=i\delta E(𝒉;p^{},p,k)𝑭_s(𝒉;p^{},p,k)`$
$`+g^2{\displaystyle }{\displaystyle \frac{d^2𝒒^{}}{(2\pi )^2}}𝒜(𝒒^{})\{(C_s\frac{1}{2}C_\mathrm{A})`$ $`\left[𝑭_s(𝒉;p^{},p,k)𝑭_s(𝒉k𝒒^{};p^{},p,k)\right]`$
$`+\frac{1}{2}C_\mathrm{A}`$ $`\left[𝑭_s(𝒉;p^{},p,k)𝑭_s(𝒉+p^{}𝒒^{};p^{},p,k)\right]`$
$`+\frac{1}{2}C_\mathrm{A}`$ $`[𝑭_s(𝒉;p^{},p,k)𝑭_s(𝒉p𝒒^{};p^{},p,k)]\},`$ (88)
where

$$\delta E(𝒉;p^{},p,k)=\frac{m_{\mathrm{eff},\mathrm{g}}^2}{2k}+\frac{m_{\mathrm{eff},s}^2}{2p}\frac{m_{\mathrm{eff},s}^2}{2p^{}}+\frac{𝒉^2}{2pkp^{}}$$
(89)
and

$$𝒜(𝒒^{})\frac{dq^z}{2\pi }A^{}(Q)[A^{}(Q)]^{}|_{q^0=q^z}.$$
(90)
Here $`Q=(q^0,𝒒^{},q^z)`$, $`A^{}A^0A^z`$, and $`A^{}(Q)[A^{}(Q)]^{}`$ is the thermal Wightman correlator evaluated in the hard-thermal-loop approximation. Explicit formulas for this correlator in equilibrium may be found, for example, in Ref. AMY4 . However, we will make use here of the wonderful simplification found by Aurenche et al. Gelis2  showing that the integral (90) has a remarkably simple form,

$$𝒜(𝒒^{})=T\left(\frac{1}{|𝒒^{}|^2}\frac{1}{|𝒒^{}|^2+m_\mathrm{D}^2}\right).$$
(91)
### B.2 Variational solution

One way to solve the integral equation for $`𝑭_s`$ is to use a variational method, similar to the method used in the main text for the Boltzmann equation and implemented in the Abelian case for $`12`$ processes in Ref. AMY3 . Some further simplifications are then possible, as we shall describe below. The variational formulation is

$$_s(p^{},p,k)\frac{d_sC_s\alpha }{(2\pi )^3}(𝒬_s)_{\mathrm{extremum}},$$
(92)
with

$$𝒬_s[𝑭]=\mathrm{Re}\left[(2𝒉,𝑭)\frac{1}{2}(𝑭,(i\delta E+𝒦_s)𝑭)\right],$$
(93)
$$(𝒬_s)_{\mathrm{extremum}}=\frac{1}{2}\mathrm{Re}(2𝒉,(i\delta E+𝒦_s)^12𝒉),$$
(94)
where in this context

$$(𝒇,𝒈)\frac{d^2𝒉}{(2\pi )^2}𝒇(𝒉)𝒈(𝒉).$$
(95)
The analog $`𝒦_s`$ of the collision operator is given by

$$(𝑭,𝒦_s𝑭)=g^2\left\{(C_s\frac{1}{2}C_\mathrm{A})(𝑭,K(k)𝑭)+\frac{1}{2}C_\mathrm{A}(𝑭,K(p^{})𝑭)+\frac{1}{2}C_\mathrm{A}(𝑭,K(p)𝑭)\right\},$$
(96)
where

$$(𝑭,K(r)𝑭)\frac{1}{2}\frac{d^2𝒉}{(2\pi )^2}\frac{d^2𝒒^{}}{(2\pi )^2}𝒜(𝒒^{})\left[𝑭(𝒉)𝑭(𝒉r𝒒^{})\right]^2.$$
(97)
Rotational invariance implies that the extremum must have the rotationally covariant form

$$𝑭_s(𝒉;p^{},p,k)=𝒉\stackrel{~}{F}_s(h;p^{},p,k).$$
(98)
The number of integrations necessary to evaluate $`(𝑭,K(r)𝑭)`$ can then be reduced by expanding the square in (97), switching integration variables to the dimensionless variables $`𝒖=𝒉/(rm_\mathrm{D})`$ and $`𝒘=(𝒉r𝒒^{})/(rm_\mathrm{D})`$, and then performing the angular integrations using the explicit form (90) for $`𝒜(𝒒^{})`$:

$`(𝒉\stackrel{~}{F},K(r)𝒉\stackrel{~}{F})=T{\displaystyle \frac{(rm_\mathrm{D})^4}{(4\pi )^2}}{\displaystyle _0^{\mathrm{}}}𝑑uu𝑑ww`$ $`\left\{{\displaystyle \frac{u^2+w^2}{|u^2w^2|}}{\displaystyle \frac{(u^2+w^2+1)}{\left[\left((u+w)^2+1\right)\left((uw)^2+1\right)\right]^{1/2}}}\right\}`$
$`\times \left[\stackrel{~}{F}(rm_\mathrm{D}u)\stackrel{~}{F}(rm_\mathrm{D}w)\right]^2.`$ (99)
Making a finite basis expansion of the form,

$$\stackrel{~}{F}(h)=\underset{m=1}{\overset{N_\mathrm{r}}{}}\stackrel{~}{F}_\mathrm{r}^{(m)}\mathrm{\Phi }_\mathrm{r}^{(m)}(h)+i\underset{m=1}{\overset{N_\mathrm{i}}{}}\stackrel{~}{F}_\mathrm{i}^{(m)}\mathrm{\Phi }_\mathrm{i}^{(m)}(h),$$
(100)
an efficient set of basis elements is AMY3

$`\mathrm{\Phi }_\mathrm{r}^{(m)}(h)`$ $`=`$ $`{\displaystyle \frac{(h^2/A)^{m1}}{\left(1+h^2/A\right)^{N_\mathrm{r}+2}}},m=1,\mathrm{},N_\mathrm{r},`$ (101)
$`\mathrm{\Phi }_\mathrm{i}^{(m)}(h)`$ $`=`$ $`{\displaystyle \frac{(h^2/A)^{m1}}{\left(1+h^2/A\right)^{N_\mathrm{i}}}},m=1,\mathrm{},N_\mathrm{i},`$ (102)
where the scale $`A`$ is chosen to optimize convergence as the basis set is increased.161616 An efficient choice of $`A`$ can be found by numerical search as described in Ref. AMY3 . However, in the present context, we have also found that the simple choice $`A=pp^{}m_{\mathrm{eff},\mathrm{g}}^2+k(p^{}p)m_{\mathrm{eff},s}^2`$ works reasonably well. This is the value of $`h^2`$ for which the last term in Eq. (89) equals the preceding terms.  Since the basis elements $`\mathrm{\Phi }`$ are only functions of $`h^2/A`$, one can see from (99) that matrix elements such as $`(\mathrm{\Phi }_\mathrm{r}^{(m)},K(r)\mathrm{\Phi }_\mathrm{r}^{(n)})`$ will equal $`T(m_\mathrm{D}r)^4`$ times a dimensionless function depending only on $`m`$, $`n`$, and $`z(m_\mathrm{D}r)^2/A`$. That means that one can numerically generate these functions just once for a given basis size (by evaluating on a fine mesh of points in $`z`$ and spline interpolating), and then repeatedly re-use their tabulated values for calculations in different theories with different values of $`m_\mathrm{D}`$ and $`A`$.

For further details relating to setting up the framework for numerical evaluation, see the treatment of photo-emission in Ref. AMY3 . The solution of the integral equation is somewhat involved but can be performed with very good numerical accuracy ($`<10^4`$ relative errors). Once the values of the $`K(r)`$ matrix elements have been tabulated, the solution of the splitting functions $`\gamma _{bc}^a`$ is quite fast. Consequently, even though these splitting functions appear inside the integrals of the Boltzmann collision term (23), the numerical cost of nested integral equations is ultimately not very large.

## Appendix C Asymptotic behavior of inverse log expansion

Here we prove that the expansion (45) in inverse powers of $`\mathrm{ln}(\mu /m_\mathrm{D})`$ is an asymptotic series with vanishing radius of convergence. We also give a (non-rigorous) argument that the expansion coefficients have non-alternating factorial growth, implying the presence of singularities in the Borel transform on the positive real axis.

### C.1 What we need to show

Within the $`\mathrm{}=2`$ or C-odd $`\mathrm{}=1`$ symmetry channels of interest, both the linearized collision operator $`𝒞`$, and its leading-log piece $`𝒜`$ defined in Eq. (38), are real, symmetric, positive definite linear operators.171717 More generally, $`𝒞`$ is positive semi-definite with zero modes associated with conserved charges: one C-even $`\mathrm{}=0`$ zero mode associated with energy conservation, one C-even $`\mathrm{}=1`$ zero mode associated with momentum conservation, and various C-odd $`\mathrm{}=0`$ zero modes associated with conservation of quark flavors. The leading-log operator $`𝒜`$ has the same structure except for one additional $`\mathrm{}=0`$ zero mode associated with total particle number (gluon plus quark plus anti-quark) conservation.  This can be seen by looking at the explicit forms (21) and (39). In other words, within the symmetry channels of interest, the operators $`𝒜`$ and $`𝒞`$ are both invertible.

However, we are not assured that the difference operator $`\overline{\delta 𝒞}(\mu /T)`$ which generates the inverse log expansion is positive definite, because its definition (40) involves a limiting procedure in which we subtract larger and larger multiples of $`𝒜`$ as we take the $`m_\mathrm{D}/T0`$ limit. There is no guarantee that adding a finite multiple of $`𝒜`$ to $`\overline{\delta 𝒞}(\mu /T)`$ will yield a positive operator, and in fact it does not. In other words, even though $`𝒞`$ is a positive definite operator, removing all subleading (in $`g`$) contributions by defining

$$𝒞|_{\mathrm{pure}\mathrm{leading}\mathrm{order}}𝒞_{\mathrm{LL}}(\mu )+\delta 𝒞(\mu )=g^4T\left[\mathrm{ln}(\mu /m_\mathrm{D})𝒜+\overline{\delta 𝒞}(\mu /T)\right],$$
(103)
yields an operator which is not positive definite. To show this, we will demonstrate that there is a family $`\varphi ^{(m)}`$ of functions such that

$$\underset{m\mathrm{}}{lim}\frac{(\varphi _{i\mathrm{}j}^{(m)},\overline{\delta 𝒞}(\mu /T)\varphi _{i\mathrm{}j}^{(m)})}{(\varphi _{i\mathrm{}j}^{(m)},𝒜\varphi _{i\mathrm{}j}^{(m)})}=\mathrm{},$$
(104)
which is enough to ensure that $`[\mathrm{ln}(\mu /m_\mathrm{D})𝒜+\overline{\delta 𝒞}(\mu /T)]`$ cannot be free of negative eigenvalues for any finite $`\mu /m_\mathrm{D}`$. Since $`𝒜`$ has positive spectrum, this implies that the expansion (44) in inverse powers of $`\mathrm{ln}(\mu /m_\mathrm{D})`$ must have a vanishing radius of convergence.181818 Here’s a general proof. Let $`C(t)=A+tB`$, with $`A`$ an Hermitian positive definite and hence invertible operator, $`B`$ Hermitian, and $`t`$ real. A necessary (but not sufficient) condition for the convergence of the Taylor expansion of $`C(t)^1`$ in powers of $`t`$ at $`t=\overline{t}`$, is the existence of $`C(t)^1`$ for all $`t`$ between 0 and $`\overline{t}`$. As $`t`$ increases from zero, $`C(t)`$ first fails to be invertible when some eigenvalue crosses zero. If $`(v,C(\overline{t})v)<0`$ for some vector $`v`$, then one or more eigenvalues must have crossed zero for $`t<\overline{t}`$, implying that $`\overline{t}`$ is outside the radius of convergence.

### C.2 The test function sequence $`\mathbf{\{}\mathit{\varphi }^{\mathbf{(}𝒎\mathbf{)}}\mathbf{\}}`$

In the limit $`m_\mathrm{D}/T0`$, which is used to define $`𝒜`$ and $`\overline{\delta 𝒞}(\mu /T)`$, the leading-log contribution dominates transport coefficients. The leading-log result is associated with $`t`$-channel (or $`u`$-channel) $`22`$ processes with momentum transfers $`𝒒`$ in the range $`m_\mathrm{D}qT`$. For the hard particle momenta $`𝒑𝒌T`$ that dominate transport, one may then approximate

$$\chi _{i\mathrm{}j}(𝒑^{})=\chi _{i\mathrm{}j}(𝒑+𝒒)\chi _{i\mathrm{}j}(𝒑)+𝒒\mathbf{}_𝒑\chi _{i\mathrm{}j}(𝒑)$$
(105a)
and
$$\chi _{i\mathrm{}j}(𝒌^{})=\chi _{i\mathrm{}j}(𝒌𝒒)\chi _{i\mathrm{}j}(𝒌)𝒒\mathbf{}_𝒌\chi _{i\mathrm{}j}(𝒌)$$
(105b)
in the collision term (7). This approximation, plus the neglect of screening, reduces $`𝒞^{22}/(g^4T)`$ to $`𝒜`$ times a logarithmically divergent integral $`𝑑q/q`$. One or the other of these simplifications break down outside the region $`m_\mathrm{D}qT`$, and so this integral more properly gives a factor of $`\mathrm{ln}(T/m_\mathrm{D})`$, sometimes called a Coulomb logarithm. For details, see Ref. AMY1  (or the somewhat analogous discussion of the non-relativistic case in Ref. LifshitzPitaevskii ).

Now suppose that we artificially consider functions $`\chi _{i\mathrm{}j}(𝒑)`$ that are similar to those one finds in a leading-log calculation of transport coefficients but which abruptly cut off above some momentum scale $`p_{}`$ with $`m_\mathrm{D}p_{}T`$. Specifically, consider replacing $`\chi (p)`$ by the test function

$$\varphi ^{()}(p)=\frac{p^{\mathrm{}}}{T}e^{p/p_{}}.$$
(106)
The small $`𝒒`$ approximations (105) now break down for $`pp_{}`$ rather than $`pT`$. Therefore, after making the same approximations as before, one might expect to find the same leading-log result for the collision matrix element except for a replacement of $`\mathrm{ln}(T/m_\mathrm{D})`$ by $`\mathrm{ln}(p_{}/m_\mathrm{D})`$. In other words,

$$(\varphi _{i\mathrm{}j}^{()},\overline{𝒞}\varphi _{i\mathrm{}j}^{()})=(\varphi _{i\mathrm{}j}^{()},𝒜\varphi _{i\mathrm{}j}^{()})\left[\mathrm{ln}(p_{}/m_\mathrm{D})+O(1)\right].$$
(107)
We will sketch a more detailed argument momentarily. Subtracting the same matrix element of $`𝒜\mathrm{ln}(\mu /m_\mathrm{D})`$ from both sides, and then taking the limit $`m_\mathrm{D}/T0`$ (as dictated in the definition (40) of $`\overline{\delta 𝒞}`$) with $`p_{}/T`$ held fixed, yields

$$\frac{(\varphi _{i\mathrm{}j}^{()},\overline{\delta 𝒞}(\mu /T)\varphi _{i\mathrm{}j}^{()})}{(\varphi _{i\mathrm{}j}^{()},𝒜\varphi _{i\mathrm{}j}^{()})}=\mathrm{ln}(p_{}/\mu )+O(1).$$
(108)
Consequently, to generate a sequence of functions demonstrating (104) one may simply take, for example, $`p_{}=e^mT`$ with $`m=1,2,3,\mathrm{}`$.

### C.3 Details

To justify the estimate (107) adequately, we must show that there are no other contributions to the matrix element which are as large as the small exchange momentum contribution already considered. For simplicity of presentation, we will choose to discuss explicitly the case where all departure from equilibrium is carried by fermions.

As a benchmark for further discussion, it is convenient first to note the parametric size of matrix elements of $`𝒜`$. Direct evaluation of Eq. (39) using $`\chi (p)=\varphi ^{()}(p)`$ shows that the expectation value scales like $`p_{}^{2\mathrm{}+1}`$, so that

$$(\varphi _{i\mathrm{}j}^{()},𝒜\varphi _{i\mathrm{}j}^{()})\frac{p_{}^{2\mathrm{}+1}}{T^3}$$
(109)
with comparable contributions coming from gauge boson exchange and fermion exchange diagrams. Now consider the contribution of $`22`$ processes to the same expectation value of $`\overline{𝒞}`$. For definiteness, start by considering $`t`$-channel processes and work in the $`t`$-channel representation of (75). The integrand of the $`22`$ contribution (22) will be exponentially suppressed by the test function factor $`[\varphi _{i\mathrm{}j}^{()}(𝒑)+\mathrm{}\varphi _{i\mathrm{}j}^{()}(𝒌^{})]^2`$ unless at least one of the momenta $`p`$, $`p^{}`$, $`k`$, or $`k^{}`$ are $`O(p_{})`$. This is only possible for $`q|\omega |\text{ }<\text{ }\text{ }p_{}`$, and either $`p`$ or $`k`$ within $`O(p_{})`$ of its lower integration limit. This gives an $`O(p_{}^2)`$ phase space suppression. For $`qT`$, we will have $`tp_{}T`$ and dominantly $`sup_{}T`$ as well \[because each Mandelstam variable can be written in the form $`(P_{}\pm P)^2=\pm p_{}p(1\mathrm{cos}\theta )`$ for some one of the external momenta $`P`$ with $`pT`$\]. Therefore,191919 We can ignore the possibility of a Bose enhancement $`f_0(p_{})T/p_{}`$ from the external particle with $`O(p_{})`$ momentum that is associated with $`\chi `$, because we are only considering the case where the fluctuation $`\chi `$ is in the fermionic sector.  the contribution of $`t`$-channel processes from the region $`qT`$ to the matrix element $`(\varphi ^{()},\overline{𝒞}\varphi ^{()})`$ is $`O(p_{}^{2\mathrm{}+2}/T^4)`$. This is smaller, by one power of $`p_{}/T`$, than the small-$`q`$ contribution which leads to the matrix element (109) for $`𝒜`$ (which also multiplies a logarithm). Contributions from smaller exchange momenta, $`qT`$, will contain greater phase space suppression which can only be compensated in those terms whose matrix elements contain small $`q`$ divergences.

For $`s`$-channel contributions, and interference terms, one finds the same result. To see this explicitly, take the matrix elements of Table 3 with the substitution (70) applied to the $`s^2/tu`$ and $`u^2/st`$ terms, so that every denominator is either a power of just $`s`$, $`t`$, or $`u`$ (or simply constant). Terms with just $`s`$ or $`s^2`$ in the denominator do not have small $`s`$ divergences because $`t`$ and $`u`$ become small as well. For terms with small $`t`$ or $`u`$ divergences, we have seen that only the small $`q`$ region can be non-negligible. Consider matrix elements that diverge as $`t^n`$ or $`u^n`$ with $`n=1`$ or 2, and focus on the case of $`t`$. For $`qp_{}`$, we have $`tp_{}^2`$, which is suppressed by a single power of $`p_{}`$ from the size $`tp_{}T`$ relevant for $`qT`$. The contribution of these terms from $`qp_{}`$ is then enhanced by $`(T/p_{})^{n1}`$ compared to the $`qT`$ region discussed previously, with one power of $`p_{}/T`$ reflecting the greater phase space suppression on $`q`$. Additionally, for the soft fermion exchange contributions to the $`qggq`$ and $`\overline{q}gg\overline{q}`$ processes, there will be an additional Bose enhancement factor of order $`f_0^\mathrm{g}(p_{})T/p_{}`$ because, for such small momentum transfers, one of the final or initial gluons must have $`O(p_{})`$ momentum if one of the initial or final quarks does.

For terms which do not involve soft fermion exchange, the resulting contribution to the matrix element of $`\overline{C}`$ is therefore $`O(p_{}^{2\mathrm{}n+3}/T^{5n})`$, which is $`(p_{}/T)^{2n}`$ times the size (109) of the matrix element of $`𝒜`$. For $`qggq`$ and $`\overline{q}g\overline{q}`$ soft fermion exchange contributions, there is an additional enhancement of $`T/p_{}`$, giving a result which is $`(p_{}/T)^{1n}`$ times (109). Since the degree of divergence $`n`$ is at most 2 for gluon exchange, and 1 for fermion exchange, this shows that even for test functions which are peaked at $`p_{}T`$, the dominant $`22`$ matrix elements are the underlined ones of Table 3 — precisely those matrix elements that contribute to the leading-log result. The resulting contributions to $`(\varphi ^{()},\overline{𝒞}\varphi ^{()})`$ are indeed dominated by $`qp_{}`$, as claimed earlier, and the contribution from $`qp_{}`$ is of the same order as the result (109) for $`𝒜`$. For $`qp_{}`$, the small $`q`$ divergences are softened by the cancellation of $`\chi `$ factors in (7), which was the basis for the leading-log approximation, reviewed earlier in discussing (105). As discussed previously, the contribution from $`m_\mathrm{D}qp_{}`$ then exceeds, by a logarithm, the contribution from $`qp_{}`$.

Analyzing $`12`$ processes is substantially more complicated, but eventually leads to the result that the expectation value $`(\varphi ^{()},\overline{𝒞}^{12}\varphi ^{()})`$ is $`O(p_{}^{2\mathrm{}+1}/T^3)`$, without any logarithmic dependence on $`p_{}/T`$ or $`m_\mathrm{D}/T`$. Hence, $`12`$ contributions remain sub-leading by a logarithm even when $`p_{}T`$. We will only outline the analysis. The factors of the test function $`\varphi ^{()}`$ in Eq. (23) cause the integrand to be exponentially small unless either $`p`$ or $`k`$ are $`\text{ }<\text{ }\text{ }p_{}`$. Symmetry of the integrand allows us to focus on $`p\text{ }<\text{ }\text{ }p_{}`$. The dominant contribution comes from $`kT`$ (and so $`p^{}T`$), as can be separately verified. Estimating the collision operator (96) as $`𝒦_s𝑭(𝒉)g^2m_\mathrm{D}^2Tk^2_𝒉^2𝑭`$, and noting that $`\delta Em_\mathrm{D}^2/p+𝒉^2/(pk^2)`$, we find that the $`\delta E`$ term dominates, so the integral equation can be solved by expanding in $`\delta E𝒦_sg^2m_\mathrm{D}^2k^2T/𝒉^2`$ with the result that $`\mathrm{Re}𝑭𝒦_s𝒉/(\delta E)^2`$. Taking $`kT`$ and substituting into all terms in Eq. (23), we find that $`_qp^2g^4T^5`$, $`\gamma _{qg}^q\gamma _{q\overline{q}}^g(g^2T)^2`$, and finally obtain $`(\varphi ^{()},\overline{𝒞}^{12}\varphi ^{()})p^{2\mathrm{}+1}/T^3`$. The contribution from $`kp_{}`$ is suppressed by a further power of $`p_{}/T`$.

### C.4 Lack of Borel summability

Having seen that $`\overline{\delta 𝒞}`$ is not relatively bounded by $`𝒜`$ \[c.f. Eq. (104)\] due to a relative logarithmic enhancement at small momentum, we can obtain a (non-rigorous) estimate of the large order behavior of the coefficients $`Q_n`$ of the inverse log expansion by restricting the domain of $`\overline{\delta 𝒞}(\mu /T)`$ to functions whose support lies in the momentum region $`p<p_{\mathrm{max}}`$ for some $`p_{\mathrm{max}}T`$. In other words, multiply $`\overline{\delta 𝒞}`$ on both sides by a projection onto this subspace. Within this subspace,202020 In more detail, let $`\varphi ^{(p_1)}`$ represent any function roughly similar to (106) with $`p_{}p_1`$ — that is, a function which is cut off for momenta large compared to $`p_1T`$. Now consider matrix elements of $`\overline{\delta 𝒞}`$ of the form $`(\varphi ^{(p_1)},\overline{\delta 𝒞}\varphi ^{(p_2)})`$. One can repeat the analysis of Appendix C.2 with the only new ingredient being that the matrix element need not be diagonal. However $`\overline{\delta 𝒞}`$ is diagonally-dominant: if $`p_1p_2`$, then the normalized off-diagonal matrix element $`(\varphi ^{(p_1)},\overline{\delta 𝒞}\varphi ^{(p_2)})/[(\varphi ^{(p_1)},\varphi ^{(p_1)})(\varphi ^{(p_2)},\varphi ^{(p_2)})]^{1/2}`$ is small compared to the normalized diagonal elements $`(\varphi ^{(p_i)},\overline{\delta 𝒞}\varphi ^{(p_i)})/(\varphi ^{(p_i)},\varphi ^{(p_i)})`$ for $`i=1,2`$ \[as can be seen from the explicit form (39)\]. When $`p_1p_2`$, the momenta $`p`$ that dominate the integrals representing the matrix element will be $`pp_1p_2`$. The logarithm factor in the discussion of Appendix C.2 will become $`\mathrm{ln}(p_1/\mu )\mathrm{ln}(p_2/\mu )`$ and we have simply replaced this by $`\mathrm{ln}(p/\mu )`$, given which momenta $`p`$ dominate. The particular ordering of factors in the expression (110) is not important at this level of precision; the given form is Hermitian (as it must be) and convenient for what follows.

$$\overline{\delta 𝒞}(\mu /T)𝒜^{1/2}\left[\mathrm{ln}(p/\mu )+K\right]𝒜^{1/2},$$
(110)
for some constant $`K`$ which depends on just how the cutoff is implemented. Inserting this approximation into the result (46) for the expansion coefficients $`Q_n`$, one has

$$Q_{n+1}(\mu )Q_{n+1}^{\mathrm{approx}}(\mu )\frac{1}{2}(1)^n(𝒮_{i\mathrm{}j},𝒜^{1/2}𝒫(p)\left[\mathrm{ln}(p/\mu )+K\right]^n𝒜^{1/2}𝒮_{i\mathrm{}j}),$$
(111)
where $`𝒫(p)=\mathrm{\Theta }(p_{\mathrm{max}}p)`$ denotes the projection operator restricting the domain of $`\overline{\delta 𝒞}`$. Assuming that $`p_{\mathrm{max}}<e^K\mu `$, the approximation (110) to $`\overline{\delta 𝒞}`$ is negative-definite, implying that the expression (111) for $`Q_{n+1}(\mu )`$ is strictly positive. Given this positivity, we can apply the same variational approach discussed in section II.4. Specifically, using the test functions $`\varphi ^{()}`$, one has the lower bound

$$Q_{n+1}^{\mathrm{approx}}(\mu )\frac{\frac{1}{2}(𝒮_{i\mathrm{}j},\varphi _{i\mathrm{}j}^{()})(\varphi _{i\mathrm{}j}^{()},𝒮_{i\mathrm{}j})}{(\varphi _{i\mathrm{}j}^{()},𝒜^{1/2}\left[\mathrm{ln}(p/\mu )K\right]^n𝒜^{1/2}\varphi _{i\mathrm{}j}^{()})}.$$
(112)
Consider now, for simplicity, the case of $`\mathrm{}=1`$ with all departure from equilibrium carried by fermions, as is relevant for flavor diffusion. An easy estimate shows that

$$(𝒮_i,\varphi _i^{()})p_{}^4,$$
(113)
while from Eq. (109) one has

$$(\varphi _i^{()},𝒜^{1/2}\left[\mathrm{ln}(p/\mu )K\right]^n𝒜^{1/2}\varphi _i^{()})p_{}^3\left[\mathrm{ln}(p_{}/\mu )K_1\right]^n,$$
(114)
for some new constant $`K_1`$. Hence

$$Q_{n+1}^{\mathrm{approx}}(\mu )K_2p_{}^5\left[\mathrm{ln}(p_{}/\mu )K_1\right]^n,$$
(115)
for some positive constant $`K_2`$. (The $`O(1)`$ constants $`K_1`$ and $`K_2`$ do not depend on $`p_{}`$ or $`n`$.) Since this is a lower bound, we may maximize over all values of $`p_{}`$ (subject to $`p_{}<e^{K_1}\mu `$), which gives $`p_{}|_{\mathrm{opt}}=e^{n/5}e^{K_1}\mu `$ and results in the estimate

$$Q_{n+1}^{\mathrm{approx}}(\mu )K_2(e^{K_1}\mu )^5e^n\left(\frac{n}{5}\right)^n.$$
(116)
The case of shear viscosity, where $`\mathrm{}=2`$ and bosons are also out of equilibrium, is similar. The non-alternating factorial growth of (116) with $`n`$ immediately implies that the inverse log expansion (45) has zero radius of convergence, and moreover that its Borel transform has a singularity on the positive real axis (at 5). Although the estimate (116) is based on the approximation (110), we do not see how inclusion of the non-logarithmically enhanced parts of $`\overline{\delta 𝒞}`$ can affect the $`n!`$ growth of the above estimate. Assuming that this estimate does give the correct large order behavior of the inverse log expansion, it is worth noting that the resulting renormalon ambiguity associated with the singularity of the Borel transform would only be of order $`(m_\mathrm{D}/T)^5g^5`$.