# On the Thermodynamics of “Continuous Spin” photons

Philip Schuster<EMAIL_ADDRESS>SLAC National Accelerator
Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA  Gowri Sundaresan
<EMAIL_ADDRESS>SLAC National Accelerator Laboratory, 2575 Sand Hill
Road, Menlo Park, CA 94025, USA  Department of Physics, Stanford University,
Stanford, CA, 94305, USA  Natalia Toro<EMAIL_ADDRESS>SLAC National
Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA

(August 28, 2024)

###### Abstract

Special relativity allows massless particles to have states of different
integer (or half-integer) helicities that mix under boosts, much like the
spin-states of a massive particle. Such massless particles are known as
“continuous spin” particles (CSPs), a term coined by Wigner, and they are
notable for their infinite tower of spin polarizations. The mixing under
boosts is controlled by a spin-scale $\rho$ with units of momentum. Normally,
we assume $\rho=0$. The interactions of CSPs are known to satisfy certain
simple properties, one of which is that the $\rho\rightarrow 0$ limit
_generically_ recovers familiar interactions of massless scalars, photons, or
gravitons, with all other polarizations decoupling in this limit. Thus, one
can ask if the photon of the Standard Model is a CSP at small but non-zero
$\rho$. One concern about this possibility – originally raised by Wigner – is
that the infinite tower of polarizations could pose problems for
thermodynamics. To address this question, we study the thermal evolution of a
CSP photon gas coupled to isothermal matter, across CSP helicity modes and
phase space. We find that the structure of the interactions dictated by
Lorentz symmetry imply well behaved thermodynamics. When the CSP photon’s
interactions to charged matter are turned on, the primary $h=\pm 1$ helicity
modes thermalize quickly, while the other modes require increasingly long
time-scales to thermalize, set by powers of $T/\rho$. In familiar thermal
systems, the CSP photon behaves like the QED photon with small $\rho$\- and
time- dependent corrections to its effective relativistic degrees of freedom.
Sizable departures from familiar thermal behavior arise at energy scales
comparable to $\rho$ and could have testable experimental consequences.

##  I Introduction

At a purely kinematic level, Lorentz symmetry allows massless particles to
have states of different integer (or half-integer) helicities that mix under
boosts much like the spin-states of a massive particle [1]. Such “continuous
spin” particles (CSPs) possess a spin-scale $\rho$ (with units of momentum)
that characterizes the extent of this mixing; in the limit $\rho\rightarrow
0$, the helicity states are exactly Lorentz-invariant. The possibility that
$\rho\neq 0$ has never been experimentally studied.

Until the work of [2, 3], it was usually assumed that theories of CSPs — if
they exist at all — are unrelated to more familiar theories and are physically
irrelevant. However, the soft factor analysis of interacting CSPs in [2, 3]
and the field theory works in [4, 5, 6] provide strong evidence that
consistent interacting CSP theories exist, and their $\rho\rightarrow 0$ limit
_generically_ recovers familiar theories of massless scalars, photons, or
gravitons. In the recent work of [6], a formalism for computing scattering
amplitudes in QED for nonzero $\rho$ was given, allowing the computation and
study of finite $\rho$ corrections to QED. (Also see [7] for a lower
dimensional generalization, [8] for fermionic and [9, 10] for supersymmetric
cases using the “vector superspace” formalism. Other formalisms, including
constrained metric-like, frame-like, and BRST [11, 12] formulations, have also
been used to construct actions for fermionic [13, 14] and supersymmetric [15,
16, 17, 18] continuous spin fields, as well as those in higher-dimensional
[19, 20, 21] and (A)dS [22, 23, 24, 25, 26, 27] spaces. Relations between
these formalisms are discussed in Refs. [28, 20, 9]. For reviews and
discussion, see Refs. [28, 29, 30, 31].)

At least for the case of a photon with finite $\rho$, it is now possible to
more precisely ask the questions: _How much_ does a photon’s helicity
transform under boosts? What constraints can be derived on the spin-scale
$\rho$ of the photon, and how might a non-zero $\rho$ be observed?

To set the stage, we recall a few basic facts about CSP kinematics and
interactions. CSP states of four-momentum $k^{\mu}$ can always be decomposed
into eigenstates of the helicity operator $\hat{k}\cdot\vec{\bf J}$:

$\hat{k}\cdot\vec{\bf J}|k,h\rangle=h|k,h\rangle,$ (1)

with integer or half-integer eigenvalue $h$. The sole change from the familiar
case $\rho=0$ to non-zero $\rho$ is that Lorentz transformations $\Lambda$
transform $|k,h\rangle$ into a superposition of states $|\Lambda
k,h^{\prime}\rangle$ with integer $h^{\prime}-h$. In the limit
$\rho\rightarrow 0$, this superposition becomes simply $|\Lambda k,h\rangle$
times a phase $e^{ih\theta(\Lambda)}$ for a suitable $\theta$. While a theory
of $\rho=0$ particles can consistently include _only_ a single helicity $h$
(or two states with $\pm h$, accounting for CPT symmetry), particles of non-
zero $\rho$ must possess an infinite tower of states of _all_ integer
helicities (or all half-integer helicities, a case we will not consider
further here). This fact suggests, at first glance, that CSP physics must be
quite different from that of the familiar theories of helcity 0, 1 and 2
particles.

However, the application of Weinberg’s famous soft limit of scattering
amplitudes revealed that CSP interactions are necessarily well approximated by
familiar gauge theories in the limit of small $\rho$ [3]. Only three types of
consistent scattering amplitudes can exist in the soft limit, and each
displays a “helicity correspondence” at CSP energies large compared to $\rho$.
In the first, scalar-like amplitude, the helicity-0 mode is emitted most
strongly, with an amplitude that is well approximated by ordinary massless
scalar interactions, while all other helicity $\pm h$ modes have emission
amplitudes suppressed by $(\rho/E)^{h}/h!$. In the other two types of
amplitude, photon-like and graviton-like interactions, respectively, the
helicity $\pm 1$ ($\pm 2$) modes are emitted most strongly, with amplitudes
well approximated by ordinary vector gauge theories (perturbative GR), while
other helicity amplitudes are suppressed by $(\rho/E)^{|h-1|}$
($(\rho/E)^{|h-2|}$.). At energies much larger than $\rho$, the other
helicities must be present in the theory, but induce only small effects. In
the more complete formalism of [6], general scattering amplitudes in QED at
finite $\rho$ retain this behavior, and the sum over all polarization states
is finite, causing no new divergences in cross sections.

The mere presence of the tower of modes raises a thermodynamic puzzle first
mentioned (in a brief remark) by Wigner [32]: in a theory with CSPs, which
have infinitely many polarization states, the vacuum’s heat capacity per unit
volume is formally infinite. This feature, which appears discontinuous between
$\rho=0$ and arbitrarily small non-zero $\rho$, is sometimes taken to preclude
the relevance of CSPs to the physical world (and is sometimes conflated with
possible inconsistencies with black hole thermodynamics in GR). In this paper,
we decisively address Wigner’s original concerns regarding CSP thermodynamics.

The key insight is that any physical system has had only finite time to
(approximately) thermalize. Helicity correspondence implies that at high-
enough energies, only one CSP mode interacts appreciably with matter, while
the others decouple. Hot matter quickly thermalizes the principally
interacting mode of a CSP (helicity $\pm 1$ in the photon-like case), but
other modes take parametrically longer to thermalize. The formal infinity of
polarization modes is never physically relevant, because only a finite number
of modes will reach equilibrium in any finite time. As such, late-time
corrections to $\rho=0$ thermodynamics are calculable and the $\rho\to 0$
limit is smooth. Indeed, for temperatures $T\gg\rho$, the effects of non-zero
$\rho$ on thermodynamics are parametrically small over exponentially long
times.

###  I.1 Summary of results

This paper clarifies and elaborates on the physical conclusion summarized
above, with a focus on QED at finite $\rho_{\gamma}$, as defined in [6], i.e.
a CSP “photon” with vector-like couplings to matter. The qualitative features
of our results would apply equally to scalar- or tensor-like CSPs, but we
choose the CSP photon as our primary example because it is the case of
greatest physical interest.

We will use the terms ‘ordinary photon’ and ‘familiar photon’ to refer to the
QED photon at $\rho_{\gamma}=0$. We will consider basic reactions that produce
and absorb on-shell photons in a bath of charged particles, which is precisely
what the formalism of [6] is appropriate for. However, for our analysis of
elementary thermodynamics, we will use soft CSP momentum limits of the
scattering amplitudes. We do this because (a) it is simpler than using full
amplitudes, and (b) it captures the parametric $\rho$\- and energy- scaling
for arbitrary processes, up to ${\cal O}(1)$ factors.

In this paper, we explore in detail the effects of non-zero spin scale
$\rho_{\gamma}$ on the phase-space evolution of a CSP photon gas coupled to a
thermal gas of charged particles. We show that for energies much greater than
the spin scale $\rho_{\gamma}$ (“correspondence domain” of the CSP gas):

  * •

If brought into contact with a thermal bath at finite temperature, a CSP
photon can be modeled as the familiar photon with small, time- and
$\rho_{\gamma}$-dependent corrections to its effective number of degrees of
freedom.

  * •

There is a strong hierarchy in mode thermalization, where all but the $h=\pm
1$ helicity modes take parametrically long to thermalize in the
$\rho_{\gamma}\rightarrow 0$ limit. The thermalization of the CSP proceeds
increasingly slowly with time.

At energies much smaller than the spin scale (“deviation domain” of the CSP
gas), we find that:

  * •

There is a weak hierarchy in mode thermalization, where a large (but finite,
${\cal O}(\rho_{\gamma}/E)$) number of modes thermalize comparably as the
$h=\pm 1$ modes, but take parametrically longer to do so at progressively
lower energies.

  * •

The radiation power spectrum at ultra-low frequencies is significantly
enhanced relative to the familiar black body spectrum.

At temperatures $T\gg\rho_{\gamma}$, these effects collectively generate
parametrically small corrections to the total energy density in radiation,
even at times exponentially larger than the photon thermalization time.

The remainder of this paper is presented as follows: We will begin our
discussion in section II with the thermodynamic setup. Subsequently, we
discuss CSP behavior in the “correspondence domain” and “deviation domain” in
turn, introducing the key ideas step-by-step in sections III and IV
respectively. Section V will serve as an extended synthesis of results, and we
will highlight a collection of open questions in VI.

##  II Setup

We start by considering an idealized thermal system: a gas of charged
particles held at a constant temperature $T$ and constant number density. The
gas can be relativistic or non-relativistic, we assume it is non-degenerate.
At times $t<0$, the CSP photon phase space density is taken to be $f_{h}=0$
for all $h$.

At time $t=0$, we ‘turn on’ the coupling of the charged particles to CSP
photons. The interactions of charged matter produce CSP photons, which begin
to thermalize. We assume, in keeping with helicity correspondence, that CSP
modes do not have appreciable self-interactions and so thermalize
predominantly via interactions with the charged scatterers. With this
assumption, the Boltzmann equation can be solved exactly to study the
thermodynamic evolution of the photon gas (see Appendix A). The phase space
density for mode $h$ evolves as:

$f_{h}(E,t)=f_{h}^{(eq)}(E)[1-\exp{(-t/\tau_{h}(E))}]$ (2)

Here, $\tau_{h}(E)$ is the ‘characteristic thermalization time’ of mode $h$ at
energy $E$, given by:

$\begin{split}\tau_{h}(E)=f_{h}^{(eq)}(E)\bigg{[}&\int
d\Pi_{in}f_{in}d\Pi_{out}(2\pi)^{4}\\\ &\times\delta^{4}(\Sigma p_{in}-\Sigma
p_{out+\gamma_{h}})|\mathcal{M}|^{2}\frac{1}{E}\bigg{]}^{-1}\end{split}$ (3)

Here,‘in’ denotes the incoming scatterers, $f_{in}$ are their phase space
distributions (which can follow any equilibrium statistics), ‘out’ denotes all
outgoing particles except CSP mode $\gamma_{h}$,
$d\Pi\equiv\frac{g}{(2\pi)^{3}}\frac{d^{3}p}{2E}$, and there is an implied
product over all species. Any averages over initial and final spins of
scattering states (besides the photon) as well as symmetry factors are
implicit, and we ignore any effects of Bose enhancement or Fermi suppression
from the outgoing states. $f_{h}^{(eq)}(E)$ is the equilibrium distribution of
mode h, which will follow Bose-Einstein statistics.

$f_{h}^{(eq)}(E)=[\exp(E/T)-1]^{-1}$ (4)

Equation (3) is completely general, but going forward, we will work in the
soft limit, with

$|\mathcal{M}|=|\mathcal{M}|_{0}\;|\Sigma_{i}q_{i}z_{i}F_{h}(\rho_{\gamma}z_{i})|$
(5)

where the subscript ‘0’ denotes the underlying scattering process to which we
attach soft photons in the charged legs denoted by ‘i’, with momenta $p_{i}$
and carrying charges $q_{i}$. The soft factor in (5) follows [2] and uses:

$\displaystyle z_{i}$ $\displaystyle\equiv\frac{\epsilon(k).p_{i}}{k.p_{i}}$
(6) $\displaystyle F_{h}(w)$
$\displaystyle\equiv\frac{(J_{h}(|w|)-c\delta_{h0})\leavevmode\nobreak\
e^{ih\leavevmode\nobreak\ \text{arg}(w)}}{|w|}$ (7)

where $\epsilon(k)$ is the polarization vector, $J_{h}(w)$ is a Bessel
function of the first kind, and $c$ can be any arbitrary constant - due to
charge conservation, the contribution from $c\delta_{h0}$ terms will cancel
when the sum over charged legs is taken. To study the thermodynamics from a
general scattering process, it will be convenient to work with a single
charged leg and leave the sum implicit, and for such calculations it is
convenient to set $c=1$ at higher energies ($E\gg\rho_{\gamma}$) and $c=0$ at
low energies ($E\ll\rho_{\gamma}$). We will use this scheme for calculations
in this paper. Note that in the $\rho_{\gamma}\rightarrow 0$ limit,
$F_{h}(\rho_{\gamma}z_{i})\rightarrow 1$ for $h=\pm 1$,
$F_{h}(\rho_{\gamma}z_{i})\rightarrow 0$ for all other $h$, and we recover the
familiar soft limit factorization of the QED photon amplitude [33] in (5).

We will use the term ‘primary modes’ to refer to $h=\pm 1$ modes of the CSP
photon, ‘partner modes’ to refer to $h=0$ and $\pm 2$ modes, and ‘sub-partner
modes’ to refer to all other helicities. For the discussion going forward, we
allow the charged scatterers to be non-relativistic, with a velocity
$v_{\perp}$ transverse to the soft photon (which was set to unity in I.1). In
the non-relativistic limit,
$\rho_{\gamma}z\approx\frac{\rho_{\gamma}v_{\perp}}{E}$. All numerical
simulations in this paper use a full thermal distribution of $v_{\perp}$,
however we use its average value $\langle v_{\perp}\rangle$ in equations. We
also note here that all numerical simulations in the paper use the full Bessel
function soft factors without any simplifications, but we will present several
equations that use approximate forms of $J_{h}(x)$ in the applicable regimes
as an aid to understanding thermodynamic behavior. Allowing for a non-
relativistic scatterer, the “correspondence domain” will refer to
$E\gg\rho_{\gamma}\langle v_{\perp}\rangle$ and the “deviation domain” will
refer to $E\lesssim\rho_{\gamma}\langle v_{\perp}\rangle$.

From equations (5), (6) and (7), it can be seen that when
$\rho_{\gamma}\langle v_{\perp}\rangle/E\ll 1$, the primary modes of the CSP
photon behave similarly as the QED photon (with corrections discussed in III).
When $\rho_{\gamma}\langle v_{\perp}\rangle/E\gg 1$, the primary modes couple
more weakly to charged matter than the QED photon, and much more
democratically as other helicites. Intuitively, the difference in behavior in
these two energy regimes can be understood as a consequence of $\rho_{\gamma}$
having mass dimension 1 - we can expect deviations from familiar behavior at
lower energies.

In the following two sections, we study these two energy regimes in detail. We
will proceed slowly, introducing key concepts as we proceed, and bring the big
picture together in Section V.

##  III  Behavior at energies above spin scale
$\bigl{(}E\gg\rho_{\gamma}\langle v_{\perp}\rangle\bigr{)}$:
“Correspondence Domain”

In this section, we study the thermal evolution and behavior of the CSP gas at
energies greater than the photon’s spin scale. We demonstrate that in this
energy regime, a CSP photon gas behaves like a thermal gas of the familiar
photon, with sub-dominant, finite, calculable corrections from all other
helicities. Since we expect the spin scale to be small in nature, the
discussion in this section applies to most energy regimes probed in familiar
physical systems (that have $T\gg\rho_{\gamma}$). Even when the thermal system
under study is ultra-cold ($T\ll\rho_{\gamma}$), its dominant behavior could
be in this energy regime if the charged scatterers are highly non-relativistic
($\langle v_{\perp}\rangle\ll 1$).

We will begin by reviewing the characteristic thermalization times in III.1,
and show that mode thermalization in the correspondence domain follows a
strong hierarchy. This sets the stage for our discussion on how the phase
space of the CSP photon is populated in III.2, with most partner and sub-
partner modes only partially thermal, and contributing only small corrections
to the relativistic degrees of freedom from the primary modes. We will
subsequently discuss the energy density of CSP gas overall in III.3. (This
discussion is possible at this juncture because all deviation domain
contributions to CSP number density and energy density are sub-dominant in
familiar thermal systems, and we will justify this in Section IV.) We
calculate the time dependence of evolution in CSP energy density explicitly,
and show that the CSP photon gas increases its internal energy progressively
slower as time progresses - thus, we demonstrate explicitly that despite
having an infinite helicity tower, the CSP photon gas has finite energy at all
finite times. We defer a discussion on the full CSP number density and degrees
of freedom to Section V.

###  III.1 Hierarchy in characteristic thermalization times

Using the equations for characteristic thermalization time introduced in II
((3), (5), (6), (7)), and using the Taylor expansion of Bessel functions of
the first kind [34, (10.2.2)] (which is valid for all helicities throughout
the correspondence domain), we find that the characteristic thermalization
time of mode $h$ of the CSP photon to leading order follows:

$\frac{\tau_{h}(E)}{\tau_{*}(E)}\sim\frac{\bigl{\langle}|z|^{2}\bigr{\rangle}}{\bigl{\langle}|zF_{h}(\rho_{\gamma}z)|^{2}\bigr{\rangle}}\;\approx\;2^{\tilde{h}}(\tilde{h}+1)!^{2}\biggl{(}\frac{E}{\rho_{\gamma}\langle
v_{\perp}\rangle}\biggr{)}^{2\tilde{h}}$ (8)

where we define $\tilde{h}\equiv||h|-1|$ to capture the mode ‘distance from
primary modes’, and we introduce a benchmark thermalization time $\tau_{*}(E)$
\- the characteristic thermalization time of an ordinary photon
($\rho_{\gamma}=0$) at the energy $E$ we are interested in.111Note that a
proper calculation of the thermal averages in (8) in the non-relativistic
limit can reduce the factorial scaling from $(\tilde{h}+1)!^{2}$ to
$(\tilde{h}+1)!$ for some helicities. Such non-relativistic effects are fully
included in all numerical simulations and elaborated in Appendix B but omitted
from the main text of the paper.

$\tau_{*}(E)\sim\tau_{*}(T)\biggl{(}\frac{E}{T}\biggr{)}^{2}$ (9)

For primary modes, $\tilde{h}=0$, and the ratio in (8) reduces to 1 - that is,
the primary modes of the CSP photon behave like the familiar photon at leading
order in this energy range, with corrections only at
$\mathcal{O}(\frac{\rho_{\gamma}\langle v_{\perp}\rangle}{E})^{2}$. The other
modes, all of which have lower production cross sections, have longer
thermalization times. The thermalization of partner modes is suppressed by
$\mathcal{O}(\frac{\rho_{\gamma}\langle v_{\perp}\rangle}{E})^{2}$ relative to
the primary mode. Sub-partner modes thermalize parametrically slower, not only
due to higher order suppressions by $(\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E})$, but also due to suppressions via
$2^{(\tilde{h}+1)}(\tilde{h}+1)!^{2}$, which grows super-exponentially in $h$.
Figure 1 demonstrates the hierarchy in characteristic thermalization times,
and the exponentially longer timescales needed for populating the phase space
of the partner and sub-partner modes.

Figure 1: These figures illustrate that CSP photon modes thermalize with a
strong “double hierarchy” in the correspondence domain
($E\gg\rho_{\gamma}\langle v_{\perp}\rangle$). For these illustrations, we
choose $T=10^{4}\rho_{\gamma}$ and $\langle v_{\perp}\rangle=0.1$.
Left: Shows the characteristic thermalization times $\tau_{h}$ of the CSP
photon modes at 3 chosen energies $E$, per equation (8). The x-axis is the
helicity $|h|$. The y-axis is logarithmic in time, and shown in units of
characteristic thermalization time $\tau_{*}(T)$ of an ordinary photon
undergoing the same thermodynamic process at temperature T.
Right: Shows the rate at which primary, partner and $h=\pm 3$ modes are
populating their phase space at a given energy (chosen to be E = T) as time
evolves, per equations (2) and (8). The x-axis is logarithmic in time, and
shown in units of benchmark thermalization time $\tau_{*}(T)$. The y-axis is
the phase space density $f_{h}$ of the mode at the chosen energy, normalized
to the equilibrium Bose-Einstein distribution.

It is also instructive to review how $\tau_{h}$ of a single mode $h$ depends
on energy. It follows from (8) and (9) that $\tau_{h}(E)\propto
E^{2(\tilde{h}+1)}$ for all $E\gg\rho_{\gamma}\langle v_{\perp}\rangle$.
Therefore the low energy phase space close to $E\sim\rho_{\gamma}\langle
v_{\perp}\rangle$ thermalizes first, whereas the higher energies take longer
to thermalize, with a greater suppression at higher helicities.

Thus, the thermalization of a CSP gas at energies above the spin scale is
“doubly hierarchical” (hierarchical in helicity and energy), with higher
energy phase space of sub-partner modes getting super-exponentially
suppressed. We next present a method to study the effects of this “double
hierarchy”.

###  III.2 Partial thermalization of modes

Due to the strong “doubly hierarchical” thermalization behavior, a CSP gas as
a whole is always partially thermal at any finite time. The phase space
density in a mode $h$ at any time $t$ can be obtained by integrating (2), and
the number density in the mode follows from it [35]:

$\displaystyle f_{h}(t)$
$\displaystyle=\int_{0}^{\infty}dEf_{h}^{(eq)}(E)[1-\exp{(-t/\tau_{h}(E))}]$
(10) $\displaystyle n_{h}(t)$
$\displaystyle=\int_{0}^{\infty}dEn_{h}^{(eq)}(E)[1-\exp{(-t/\tau_{h}(E))}]$
(11) where $\displaystyle
n_{h}^{(eq)}(E)\equiv\frac{1}{2\pi^{2}}f_{h}^{(eq)}(E)E^{2}$ (12)

Note that each mode of the CSP carries a single internal degree of freedom.

At a given time $t$, the low energy phase space for which $\tau_{h}(E)\ll t$
is approximately thermalized. This behavior continues until a maximum energy
$E_{h\wedge}(t)$ that satisfies the condition $\tau_{h}(E)=t$, above which
density is well below the thermal value. Thus, the differential phase space
density per unit energy at time $t$ (integrand of (10)) can be expressed in
terms of its equilibrium value, and the differential number density per unit
energy at time $t$ (integrand of (11)) has the same form. The latter is222Note
that per (11), at $E=E_{h\wedge}(t)$, $n_{h}$ has already fallen well below
its equilibrium value. (13) thus provides a conservative estimate.:

$n_{h}(E,t)\lesssim\begin{cases}n_{h}^{(eq)}(E)&\rho_{\gamma}\langle
v_{\perp}\rangle\leq E\leq E_{h\wedge}(t)\\\
n_{h}^{(eq)}(E)(\frac{E_{h\wedge}(t)}{E})^{2(\tilde{h}+1)}&E\geq
E_{h\wedge}(t)\end{cases}$ (13)

where all variables and parameters are defined as before. Note that a more
complete expression for number density will be discussed in Section V. We can
use (8) and (9) to express $E_{h\wedge}(t)$ as:

$E_{h\wedge}(t)\sim\biggl{[}\,\frac{t}{\tau_{*}(T)}\frac{(\rho_{\gamma}\langle
v_{\perp}\rangle)^{2\tilde{h}}}{(\tilde{h}+1)!^{2}}T^{2}\biggr{]}\,^{\frac{1}{2(\tilde{h}+1)}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for all}\leavevmode\nobreak\ {h}$ (14)

For the equilibrium phase space of a mode $h$ to be fully populated by time
$t$, its $E_{h\wedge}(t)$ has to be $\gg T$. Since $E_{h\wedge}(t)\propto
t^{\frac{1}{2(\tilde{h}+1)}}$, it increases progressively slower for higher
helicities, and remains close to $\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{|h|}$ for many orders of time after $\tau_{*}(T)$ in the
high helicity limit. Differential number density per unit energy at a given
time $t$ is illustrated in Figure 2. As we expected, most modes are only
partially thermal at any given time, and sub-partner modes occupy
progressively smaller fractions of their total available phase space.

Figure 2: This figure illustrates that at any finite time, the CSP photon gas
is only partially thermalized. We show the differential number density per
unit energy of CSP modes at a snapshot in time. The x-axis is the energy $E$,
shown in units of temperature $T$. The y-axis is the mode differential number
density per unit energy, normalized to the equilibrium value at temperature T.
As shown, the helicity $\pm 1$ modes have fully equilibrated, the partner
modes have only thermalized their phase space up to $E\sim T$, and the $h=\pm
3$ modes have thermalized a vanishingly small fraction of their phase space
(up to $E<10^{-2}T$). The parameter choices for this illustration are:
Temperature $T=10^{4}\rho_{\gamma}$, $\langle v_{\perp}\rangle=0.1$, time of
snapshot t = $10^{10}\tau_{*}(T)$.

As mentioned in the section introduction, we will look at CSP energy density
overall now, and will return to CSP number density in Section V.

###  III.3 CSP energy density

We will demonstrate two critical aspects of the CSP energy density:

  1. 1.

The total energy density is finite at all finite times

  2. 2.

The rate of increase in energy density is inverse in time (once the primary
modes are thermal3). That is, the energy density increases more and more
slowly as time evolves

The full derivation including simplifying assumptions made is provided in
Appendix C, but here we present the main results. We will use $\mathcal{E}$ to
denote energy density instead of the more commonly used $\rho$ since we
reserve the latter for spin scale.

The energy density in the CSP gas is:

$\mathcal{E}_{CSP}(t)=\sum_{h}\int_{0}^{\infty}dEn_{h}(E,t)E$ (15)

We start with the time rate of change of this.

$\displaystyle\frac{d}{dt}\mathcal{E}_{CSP}(t)$
$\displaystyle\lesssim\sum_{h}\int_{E_{h\wedge}(t)}^{\infty}dEn_{h}^{(eq)}(E)E\frac{1}{\tau_{h}(E)}$
(16)
$\displaystyle\approx\Bigg{\\{}\sum\limits_{h\,:\,\begin{subarray}{c}\tau_{h}(T)\geq
t\\\ h\neq\pm
1\end{subarray}}\frac{T^{4}}{2\pi^{2}[2\tilde{h}-1]}t^{\left(-1+\frac{3}{2(\tilde{h}+1)}\right)}$
$\displaystyle\times\left(\frac{\rho_{\gamma}^{2}}{\langle
p_{\perp}^{2}\rangle}\right)^{\frac{3\tilde{h}}{2(\tilde{h}+1)}}\left[\frac{1}{\tau_{*}(T)(\tilde{h}+1)!^{2}}\right]^{\frac{3}{2(\tilde{h}+1)}}\Bigg{\\}}$

All variables and parameters in (III.3) are as defined before. The conditional
sum singles out all the helicities that are still populating their phase space
up to $E=T$, since at $E>T$, energy density increase is exponentially
suppressed. We assume the primary modes are thermal and use $\tau_{*}(T)$ as
the benchmark time333When thermalizing, the energy density in primary mode
will not change per the form in (III.3) but instead follows
$\frac{d}{dt}\mathcal{E}_{*}(t)\propto t^{\frac{1}{2}}$. When a mode is
thermalizing, we can generally expect its energy density to change roughly in
line with $TE_{h\wedge}^{3}(t)$. Refer Appendix C for more details. .

Even though (III.3) looks forbidding, its key features are straightforward.
First, the rate of change in energy density in all modes in inverse in time,
since $\frac{d\mathcal{E}_{h}}{dt}\propto
t^{\bigl{(}-1+\frac{3}{2(\tilde{h}+1)}\bigr{)}}$ for all $\tilde{h}\neq 0$.
This means that even though the mode energy density is increasing with time,
it does so more and more slowly. Second, the sum in (III.3) is convergent,
ensured by the $\tilde{h}^{4}$ factor in the denominator. Lastly,
$\rho_{\gamma}$ hierarchically controls which helicities participate in the
conditional sum, via its control on $\tau_{h}(T)$.

Due to these factors, the energy density of the CSP as a whole increases very
slowly and in a well-controlled manner. At any time, the sum over modes in
(III.3) is dominated by the nearest non-thermal mode (smallest $\tilde{h}$
participating in the convergent sum), with thermalized modes dropping out of
the conditional sum. (15) and (III.3) point to a simplified model for the CSP
energy density:

$\mathcal{E}_{CSP}(t)<\sum\limits_{h\,:\,\tau_{h}(T)<t}\frac{\pi^{2}}{30}T^{4}+\sum\limits_{h\,:\,\tau_{h}(T)\gtrsim
t}\int_{0}^{t}dt\frac{d}{dt}\mathcal{E}_{CSP}(t)$ (18)

In (18), both terms are conditional sums: the first sum is taken over modes
that have thermalized by time $t$, whereas the second sum accounts for energy
density in modes that are still thermalizing at time $t$. We note here that
(18) overestimates the energy density at any time, since we assume that all
the available energy density in a mode i.e., $\frac{\pi^{2}}{30}T^{4}$ [35] is
unlocked when it is thermal up to $E=T$. This is a very conservative
assumption since only $\approx 3.5\%$ of total energy density is unlocked by
the time any Bose gas has fully thermalized its phase space up to temperature
$T$. Additionally, specifically for CSP photon modes, unlocking the additional
energy density in the region of phase space above $T$ is parametrically harder
since characteristic thermalization times scale as $E^{2\tilde{h}}$ per (8).
Nevertheless, we adopt this toy model for its tractability and the aid it
provides in building intuition.

From (III.3) and (18), it is easy to see that when a mode is thermalizing,
that is, at times $t\lesssim\tau_{h}(T)$, its energy density is increasing as:

$\mathcal{E}_{h}(t)<\frac{3}{(\tilde{h}+1)^{3}}t^{\frac{3}{2(\tilde{h}+1)}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for}\leavevmode\nobreak\
h\,:\,\begin{subarray}{c}\tau_{h}(T)\geq t\\\ \tilde{h}\neq 0\end{subarray}$
(19)

Table 1 enlists how $\mathcal{E}_{h}(t)$ is varying for select modes, and
compares that with behavior of the primary modes3 given by
$\mathcal{E}_{*}(t)$. (19) implies that the sum over thermalizing modes in
(18) is convergent at any time. Since there are a finite number of thermal
modes at any finite time, the sum over thermalized modes in (18) is also
convergent at any time. Thus, the CSP photon gas always carries finite energy
at all finite times.

Table 1: Upper bounds on increase in mode energy density with time during thermalization (when $\tau_{h}(T)\gtrsim t$) for select modes. $\mathcal{E}_{*}(t)$ is shown in the first row for comparison. Mode | Time dependence of $\mathcal{E}_{h}(t)$ (slower than):
---|---
$\mathcal{E}_{*}(t)$ = $\mathcal{E}_{\pm 1}(t)$ | $t^{3/2}$
$\mathcal{E}_{0,\pm 2}(t)$ | $t^{3/4}$
$\mathcal{E}_{\pm 3}(t)$ | $t^{1/2}$
$\mathcal{E}_{h}(t)\leavevmode\nobreak\ \textrm{as}\leavevmode\nobreak\ h\rightarrow\infty$ | Frozen/ Constant

This is a remarkable result - even when supplied with an isothermal bath that
makes infinite energy available to it, the CSP gas takes infinite time to
access that energy. This is the crux of why Wigner’s infinite heat capacity
argument is not physically relevant. Figure 3 illustrates these aspects of CSP
energy density.

Figure 3: This figure illustrates that the CSP photon gas has finite energy
and increases its energy density progressively slowly with time, as indicated
by (III.3) and (19). The x-axis is logarithmic in time, and shown in units of
benchmark thermalization time $\tau_{*}(T)$ \- the time taken for a QED photon
undergoing the same thermodynamic process to populate its phase space up to
$E=T$. The y-axis is the energy density, normalized to the total energy
density in a fully thermalized Bose-Einstein distribution
$\frac{\pi^{2}}{30}T^{4}$ [35]. We show a log-log plot to make the slowing
growth rate in energy density manifest. It can be seen that in
$10^{100}\tau_{*}(T)$, the CSP photon gas has unlocked $\approx 18.7$ degrees
of freedom of the infinite degrees of freedom it has in principle, with a
decelerating rate of growth. The parameter choices for this illustration are:
Temperature $T=10^{4}\rho_{\gamma}$, and $\langle v_{\perp}\rangle=0.1$.
Simulations account for mode thermalization behaviors over the full energy
range $E=[0,\infty)$, and the dashed lines indicate analytical extensions to
simulations.

##  IV Behavior at energies below spin scale
$\bigl{(}E\lesssim\rho_{\gamma}\langle v_{\perp}\rangle\bigr{)}$: “Deviation
Domain”

As summarized in Section I.1, the CSP photon gas has a fundamentally
different, but still well-controlled behavior at energies
$\lesssim\rho_{\gamma}\langle v_{\perp}\rangle$. Since we expect the spin
scale $\rho_{\gamma}$ to be small, the discussion in this section applies to a
very small volume of the overall phase space of the CSP gas in familiar
thermodynamic systems (with $T\gg\rho_{\gamma}$). Nevertheless, we devote some
attention to this since thermal CSP photons in this energy regime exhibit
interesting deviations vs. thermal QED photons. In ultra-cold thermal systems
(those that have $T\sim\rho_{\gamma}\langle v_{\perp}\rangle$), most of the
phase space will be in the deviation domain, and this could have interesting
phenomenological implications. Additionally, the well-controlled nature of CSP
thermal behavior in the deviation domain is not readily apparent and requires
different arguments vs. in the correspondence domain. As before, we begin by
studying the characteristic thermalization times.

###  IV.1 Weaker hierarchy in
characteristic thermalization times

In this subsection, we show that there is a very weak hierarchy in mode
thermalization at energies much smaller than $\rho_{\gamma}\langle
v_{\perp}\rangle$. We will see that this weak hierarchy is balanced by the
parametrically longer thermalization times as we move down the energy scale
and $E\rightarrow 0$.

Using the equations for characteristic thermalization time introduced in II
((3), (5), (6), (7)), we find that the characteristic thermalization time of
mode $h$ at energy $E\ll\rho_{\gamma}\langle v_{\perp}\rangle$ is444For
$\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{1+\varepsilon}<|h|<\frac{\rho_{\gamma}}{E}$,
Bessel functions cannot be expressed in any simpler form for the entire range
of thermal integration in (3). Nevertheless, it can be numerically verified
that the Bessel amplitudes fall off rapidly as mode numbers increase. For a
relativistic scatterer, the second case in (20) applies for
$|h|\gg(\frac{\rho_{\gamma}}{E})^{1+\varepsilon}$.:

$\displaystyle\frac{\tau_{h}(E)}{\tau_{*}(E)}\sim\begin{cases}\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{3}&|h|\lessapprox\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{1+\varepsilon}\\\
2^{\tilde{h}}(\tilde{h}+1)!^{2}\biggl{(}\frac{E}{\rho_{\gamma}\langle
v_{\perp}\rangle}\biggr{)}^{2\tilde{h}}&|h|>\frac{\rho_{\gamma}}{E}\\\
&\text{and}>\frac{1}{2\langle v_{\perp}\rangle^{2}}\\\ \end{cases}$ (20)

Here, $0<\varepsilon\ll 1$. $\varepsilon$ has been used to indicate that more
modes with $|h|\sim\mathcal{O}\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}$ follow similar behavior as the
$h\lessapprox\bigl{(}\frac{\rho_{\gamma}\langle v_{\perp}\rangle}{E}\bigr{)}$
modes. The benchmark thermalization time $\tau_{*}(E)$ still follows (9). The
low helicity case above uses the asymptotic form [34, (10.17.3)] of Bessel
functions of the first kind. For the high helicity case, its Taylor expansion
[34, (10.2.2)] is valid, just as in the correspondence domain. Equation (20)
implies that in the deviation domain, CSP modes follow a fundamentally
different thermalization behavior vs. in the correspondence domain. We
highlight three key aspects of this behavior.

First, whereas ordinary photons thermalize _more rapidly_ at lower energies
since $\tau_{*}(E)$ scales as $E^{2}$, CSP photon modes with
$|h|\lessapprox\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{1+\varepsilon}$ thermalize _less rapidly_ at low
energies since $\tau_{h}(E)$ scales as $E^{-1}$. So, at any time $t$, we can
define the lowest energy $E_{h\vee}(t)$ down to which a CSP mode is thermal.
We return to this in IV.2 and IV.3.

Second, in the deviation domain the primary mode is no longer ‘special’.
Instead, as we move to lower energies, an increasing number of modes (those
with helicity $\lesssim\frac{\rho_{\gamma}\langle v_{\perp}\rangle}{E}$)
thermalize on nearly identical timescales, but this timescale also increases
as a positive power of $\frac{\rho_{\gamma}\langle v_{\perp}\rangle}{E}$. This
behavior might seem surprising at first - but it is simply a consequence of
lower scattering amplitudes at lower energies, and these scattering amplitudes
having a weaker dependence on helicity (up to a limit). As we will see later,
this beautiful balance between the number of modes thermalizing and the
increasing thermalization time keeps the low energy phase space of the CSP
photon well-controlled at all finite times.

Third, modes with high helicity ($|h|>\frac{\rho_{\gamma}}{E}$) still take
parametrically longer to thermalize - as dicussed above, these modes behave
exactly like they do in the correspondence domain, with their thermalization
times following (8). We provide additional perspectives on this in IV.3.

Thus, we still have a hierarchical thermalization behavior, albeit a much
weaker one when compared to that in the correspondence domain. The hierarchy
in characteristic thermalization times in the deviation domain is illustrated
in Figure 4.

Figure 4: This figure illustrates the controlled thermalization behavior of
the CSP photon modes in the deviation domain. It can be seen that the
thermalization behavior follows a weak hierarchy, as given by (20), with up to
$\mathcal{O}\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{1+\varepsilon}$ modes identically thermalizing
as the primary modes. For the three energies illustrated, it can be seen that
$\varepsilon\ll 1$. A strong hierarchy in thermalization behavior, akin to
that in correspondence domain, sets in at higher helicities. We illustrate
these behaviors at three energies, each separated by an order of magnitude.
The x-axis is the helicity $|h|$. The y-axis is logarithmic in time, and shown
in units of benchmark thermalization time $\tau_{*}(T)$. The parameter choices
for this illustration are: Temperature $T=10^{4}\rho_{\gamma}$ and $\langle
v_{\perp}\rangle=0.1$.

We studied the characteristic thermalization times at
$E\gg\rho_{\gamma}\langle v_{\perp}\rangle$ using (8) and at
$E\ll\rho_{\gamma}\langle v_{\perp}\rangle$ using (20). We now talk about an
intuitive way to understand mode thermalization behavior in the full energy
range, including intermediate energies $E\sim\rho_{\gamma}\langle
v_{\perp}\rangle$.

###  IV.2 Characteristic energy of a mode

We have seen that thermalization time of a given mode increases both at very
high energies and at very low energies, with the former controlled by the
first term in the Taylor expansion of the Bessel function in (7) and the
latter by its large-argument asymptotic scaling. In between these two regimes,
CSP emission amplitudes are not readily approximated but their behavior can be
studied analytically with Bessel functions, and thus can be easily bounded
[34, (10.14.1, 10.14,2)]. An important energy scale in the problem is the
_characteristic energy_
<EMAIL_ADDRESS>2.33333pt}(h)$ of a given mode, at which that mode thermalizes most rapidly.

<EMAIL_ADDRESS>2.33333pt}(h)\approx\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{f^{{}^{\prime}}_{h,1}}=\begin{cases}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{j^{{}^{\prime}}_{h,1}}\sim\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{|h|}&h\neq 0\\\ \frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{1.52}&h=0\end{cases}$ (21)

In (21), $f^{{}^{\prime}}_{h,1}(x)$ refers to the first positive zero of the
derivative of $F_{h}(x)$ as defined in (7), $j^{{}^{\prime}}_{h,1}(x)$ refers
to the first positive zero555The error in approximating
$j^{{}^{\prime}}_{h,1}(x)$ with $|h|$ decreases as $|h|$ increases. See [36]
of the derivative of $J_{h}(x)$ [36], and $f^{{}^{\prime}}_{0,1}(x)$ follows
from using $c=0.5$ in (7). Note that
<EMAIL_ADDRESS>2.33333pt}(h)\leq\rho_{\gamma}\langle v_{\perp}\rangle$ for all $h$.

At energies below its characteristic energy, mode thermalization time
increases inversely with energy, and at $E\ll
<EMAIL_ADDRESS>2.33333pt}(h)$, it follows the equation for the low helicity case in (20). At
energies greater than its characteristic energy, mode thermalization time
increases with energy, and at $E\gg
<EMAIL_ADDRESS>2.33333pt}(h)$, it follows (8). At
<EMAIL_ADDRESS>2.33333pt}(h)$, the mode thermalization time has a global minimum:

<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}}{\tau_{*}(T)}\sim\begin{cases}\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{T}\bigr{)}^{2}|h|&h\neq 0\\\
\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{T}\bigr{)}^{2}1.52&h=0\end{cases}$ (22)

Since a mode thermalizes fastest at its characteristic energy,
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$ gives the lower bound on $\tau_{h}(E)$. Since this
lower bound is monotonic in $|h|$, it is straightforward to see that
$|h|\rightarrow\infty$ modes need infinite time even to populate the phase
space close to their characteristic energy, which is $\approx$ zero.

###  IV.3 Mode thermalization behavior

The notion of characteristic energy (21), and the lower bound on
thermalization time it provides (22) allow us to understand mode
thermalization behavior throughout the deviation domain, including at
intermediate energies $E\sim\rho_{\gamma}\langle v_{\perp}\rangle$. At time
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$, no region of mode phase space has populated
appreciably. At these times,

$n_{h}(E,t)\approx n_{h}^{(eq)}(E)\frac{t}{\tau_{h}(E)}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for all}\leavevmode\nobreak\ \leavevmode\nobreak\
E\leavevmode\nobreak\ \text{at}\leavevmode\nobreak\
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$ (23)

At
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$, the region of phase space close to
<EMAIL_ADDRESS>2.33333pt}(h)$ is mostly populated, but the rest of mode phase space is non-
thermal. As time evolves, phase space at energies higher and lower than the
characteristic energy thermalize.

At energies greater than
<EMAIL_ADDRESS>2.33333pt}(h)$, we can define the maximum energy up to which the mode phase
space is thermal at a given time: $E_{h\wedge}(t)$, introduced and discussed
in Section III.2. To track how the region of phase space at $E\gg
<EMAIL_ADDRESS>2.33333pt}(h)$ is populated as time evolves, (13) and the discussion
associated with its evolution is valid. It can be verified that at
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$, (14) reduces to
<EMAIL_ADDRESS>2.33333pt}(h)$ in the high helicity limit.

At energies less than
<EMAIL_ADDRESS>2.33333pt}(h)$, we can define an analogous quantity $E_{h\vee}(t)$: the
minimum energy down to which the mode phase space is thermal at a given time.
$E_{h\vee}(t)$ tracks how the phase space below
<EMAIL_ADDRESS>2.33333pt}(h)$ is populated as time evolves.

$E_{h\vee}(t)\sim\frac{(\rho_{\gamma}\langle
v_{\perp}\rangle)^{3}}{T^{2}}\frac{\tau_{*}(T)}{t}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$ (24)

So at energies below characteristic energy, the differential number density
per unit energy at time $t$ can be expressed as:

$n_{h}(E,t)\lesssim\begin{cases}n_{h}^{(eq)}(E)&E_{h\vee}(t)\leq E\leq
<EMAIL_ADDRESS>2.33333pt}(h)\\\ n_{h}^{(eq)}(E)\frac{E}{E_{h\vee}(t)}&E\leq
E_{h\vee}(t)\end{cases}$ (25)

We can now gain some intuition for the mode thermalization hierarchies in the
deviation domain (20).

At any given energy $E<\rho_{\gamma}\langle v_{\perp}\rangle$, the modes that
near-identically thermalize are those that have characteristic energies
greater than $E$ \- so they are all evolving in the phase space region lower
than their respective characteristic energies in accordance with (25), with
their $E_{h\vee}(t)$ changing as in (24) - the helicity independence of (24)
explains why their evolution is near-identical. As we go down the energy
scale, we pass the characteristic energies of more modes per (21), and these
new modes also start evolving in accordance with (25) and (24).

At any given energy $E$, the modes that thermalize with parametrically longer
times in (20) are those that have characteristic energies lower than $E$, so
they are evolving in the phase space region higher than their respective
characteristic energies, with thermalization behavior approaching (8) at $E\gg
<EMAIL_ADDRESS>2.33333pt}(h)$. When we get to the correspondence domain, all modes have
characteristic energies smaller than the energy we are interested in, and all
modes thermalize per (8).

###  IV.4 CSP number density, energy density and
spectrum of thermal radiation

We now discuss three key aspects of the CSP number density and energy density
in the deviation domain:

  1. 1.

The total number density and energy density are both finite at all finite
times

  2. 2.

The total energy density in the deviation domain is sub-dominant to the total
energy density in the correspondence domain (for systems with
$T\gg\rho_{\gamma}$)

  3. 3.

The thermal radiation from a CSP photon gas is expected to be stronger vs. the
standard black body spectrum of a thermal QED gas

First, CSP number density and energy density in the deviation domain are
finite since only a finite number of modes contribute to them at all finite
times. Those modes for which
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}>t$ would not have populated any portion of their phase
space to equilibrium densities, and since
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}$ varies as a positive power of $|h|$, infinite time is
needed for the deviation domain of the CSP to fully thermalize. Additionally,
the non-equilibrium contributions summed over all CSP modes is also finite
since the total number of contributing modes is almost exactly offset by the
fractional thermalization at that energy (per (23) and (20)) - so, the low
volume of phase space at these energies fully regulates the number density and
energy density in the deviation domain at all finite times.

Second, the total energy density in the region of phase space below
$\rho_{\gamma}\langle v_{\perp}\rangle$ remains sub-dominant (assuming
$T\gg\rho_{\gamma}\langle v_{\perp}\rangle$), despite many more modes
thermalizing simultaneously. It can be shown that666The total energy density
in any mode is most sensitive to its $E_{h\wedge}(t)$, and relatively
insensitive to its $E_{h\vee}(t)$, so we can get an order of magnitude
estimate of the total energy in any mode by tracking only the behavior of its
$E_{h\wedge}(t)$. This simplifying assumption can be further justified by
noting that the total energy contribution to the CSP from thermalization of
every mode’s phase space region $E_{h\vee}(t)\leq E\leq
<EMAIL_ADDRESS>2.33333pt}(h)$ is at all times bounded by $\mathcal{O}((\rho_{\gamma}\langle
v_{\perp}\rangle)^{3}T)$. The total contribution to CSP deviation domain
energy from every mode’s
<EMAIL_ADDRESS>2.33333pt}(h)\leq E\leq E_{h\wedge}(t)$ is dominated by the contribution from
the highest $E_{h\wedge}(t)$ at any time $t$ i.e., the energy at the ‘boundary
region’ $E\sim\rho_{\gamma}\langle v_{\perp}\rangle$. we can get an order of
magnitude estimate of the total energy in the deviation domain of the CSP gas
using only the energy in the ‘boundary region’ $E\sim\rho_{\gamma}\langle
v_{\perp}\rangle$, which can be expressed in the suggestive form:

$\mathcal{E}_{deviation\leavevmode\nobreak\
domain}\sim\mathcal{O}\biggl{(}\leavevmode\nobreak\
\biggl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{T}\biggr{)}^{3}h_{max}T^{4}\biggr{)}$ (26)

where $h_{max}$ is the highest helicity that has populated its phase space at
$E=\rho_{\gamma}\langle v_{\perp}\rangle$, given by the mode that saturates
the condition $E_{h\wedge}(t)=\rho_{\gamma}\langle v_{\perp}\rangle$ in (14):

$\frac{t}{\tau_{*}(T)}\biggl{(}\frac{T}{\rho_{\gamma}\langle
v_{\perp}\rangle}\biggr{)}^{2}\approx 2^{h_{max}}\leavevmode\nobreak\
h_{max}!^{2}$ (27)

From (26), it is clear that for systems with $\rho_{\gamma}\langle
v_{\perp}\rangle\ll T$, the deviation domain of the CSP gas contains a
negligible fraction of the total energy in even a single fully thermalized
mode of a Bose gas, which is $\mathcal{O}(T^{4})$. The only way this can be
circumvented is if at a given time, $h_{max}$ can be
$>\mathcal{O}\biggl{(}\frac{T}{\rho_{\gamma}\langle
v_{\perp}\rangle}\biggr{)}^{3}$ _and_ those modes which have
$E_{h\wedge}(t)>\rho_{\gamma}\langle v_{\perp}\rangle$ at that time (which are
necessarily smaller helicities than $h_{max}$) have not thermalized much of
their correspondence domain. This condition can never be met since the energy
scaling of characteristic thermalization time at $E\gg
<EMAIL_ADDRESS>2.33333pt}(h)$ is $\tau_{h}\propto E^{2(\tilde{h}+1)}$ whereas the helicity
scaling is the much stronger $\tau_{h}\propto(\tilde{h}+1)^{2(\tilde{h}+1)}$.
Simply put, it is _much_ easier for smaller helicities to populate their
higher energy phase space than it is for larger helicities to even populate
their phase space at energies close to their characteristic energies - at any
time, it is impossible for $h_{max}$ to grow large enough to compensate for
the energy contained in the correspondence domain of the smaller helicities.

Third, despite the sub-dominance of the total deviation domain energy of the
thermal CSP gas, its black body spectrum shows substantial deviations from
that of the familiar thermal QED gas, both in radiated power and spectrum
shape. Due to near-identical thermalization of up to
$\mathcal{O}\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{1+\varepsilon}$ modes at energies
$E_{h\vee}(t)\leq E\leq\rho_{\gamma}\langle v_{\perp}\rangle$, the CSP
differential energy density stays nearly linear777There will be deviations
from linearity because: a) $0<\varepsilon\ll 1$, and b) If we wait long
enough, we can expect contributions from the modes that are populating the
region of their phase space above their characteristic energies. down to a
time-dependent low energy cutoff. Below this low energy cutoff, no mode has
populated its phase space to the equilibrium value, but non-equilibrium
contributions of up to $\mathcal{O}\bigl{(}\frac{\rho_{\gamma}\langle
v_{\perp}\rangle}{E}\bigr{)}^{1+\varepsilon}$ modes can still add up per (23),
resulting in a quadratic spectrum shape similar to QED gas in the deep IR,
albeit with stronger radiated power888Small deviations from quadratic form can
exist in principle because $0<\varepsilon\ll 1$. However, as can be seen in
Figure 4, $\varepsilon$ decreases as we move to lower energies, so these
deviations are likely to be negligible.. Thus the CSP gas radiates stronger
than the QED gas at all frequencies in the deviation domain, with exact
spectrum shape dependent on allowed thermalization time, the spin scale and
temperature. This is in contrast to the QED photon thermal spectrum that
remains quadratic in frequency at all $\omega\ll T$. Since the power radiated
by the CSP photon at these very low frequencies is spread across modes with
suppressed matter interactions, further study is needed to evaluate the
detectability of this excess.

Figure 5 illustrates all 3 aspects of CSP deviation domain behavior discussed
above.

Figure 5: This figure illustrates the differences in the spectrum of thermal
radiation of a CSP photon gas vs. a QED photon gas, at time
$t=10^{10}\tau_{*}(T)$, with discernible deviations from the familiar photon
blackbody spectrum prominent at all deviation domain frequencies. We show a
log-log plot to make the behavior across frequencies manifest. The x-axis is
the radiated frequency $\omega$ in units of temperature $T$. The y-axis shows
the power radiated at that frequency, in units of power radiated by a fully
thermalized QED photon gas at the frequency $\omega=T$. The lowest frequency
that has been populated to equilibrium density by any mode at this time is
$\omega\approx 10^{-25}T$. The dashed portion of the green CSP line indicates
an analytical extension to lower frequencies, where all energy radiated is
from non-equilibrium contributions summed over all modes. Weaker deviations
from the standard spectrum are also apparent at frequencies close to $T$ since
the partner modes have populated some of the correspondence domain phase space
at this time. Despite these strong deviations from standard spectrum in the
deviation domain, the total energy radiated by the CSP photon gas is strongly
dominated by frequencies higher than $\rho_{\gamma}\langle v_{\perp}\rangle$
\- in this figure, the total energy radiated at frequencies above
$\rho_{\gamma}\langle v_{\perp}\rangle$ is $\approx 10^{15}\times$ total
energy radiated at frequencies below $\rho_{\gamma}\langle v_{\perp}\rangle$,
a fact which might be obscured by the log-log scale. The parameter choices for
this illustration are: Temperature $T=10^{4}\rho_{\gamma}$, $\langle
v_{\perp}\rangle=0.1$ and time of snapshot $t=10^{10}\tau_{*}(T)$.

##  V Synthesis: The effect of a non-zero spin scale on thermodynamic behavior

So far, we have discussed the thermalization behavior of the CSP photon in the
two energy regimes ($E\gg\rho_{\gamma}\langle v_{\perp}\rangle$ and
$E\leq\rho_{\gamma}\langle v_{\perp}\rangle$), examining the key aspects of
the behavior and building intuition step by step. In this section, we bring
together a synthesis of the key results already discussed and present the
complete picture of CSP thermalization behavior across all energies and times.
We supplement this synthesis with some additional salient aspects of CSP
thermodynamics whose discussion required the complete picture.

###  V.1 Overall mode thermalization behavior

This subsection first brings together the key aspects of thermalization
behavior of a single mode across phase space. Subsequently, we synthesize all
the helicity-based thermalization hierarchies.

####  V.1.1 The complete picture of mode thermalization

Working in the soft limit, we saw that the thermalization behavior of all CSP
modes follows the same pattern: Each mode has a characteristic energy, given
by (21), at which it its thermalization time is the shortest. As we move along
the energy scale in either direction and away from the characteristic energy,
the characteristic thermalization time of the mode increases, with
$\tau_{h}\rightarrow\infty$ as $E\rightarrow 0$ and $E\rightarrow\infty$. At
energies much greater than its characteristic energy, a mode’s characteristic
thermalization time grows with energy, following (8). At energies much less
than its characteristic energy, the mode’s characteristic thermalization time
is inverse in energy, following the low helicity case in (20). Notably, this
behavior is also followed by the primary modes, whose characteristic energy is
$\sim\rho_{\gamma}\langle v_{\perp}\rangle$.

Figure 6 illustrates the behavior of the characteristic thermalization time
across phase space for select modes.

Figure 6: This log-log plot illustrates that for every thermalizing CSP photon
mode, the characteristic thermalization times diverge in the UV and IR (but
approach these divergences at different rates), keeping the total energy in
the CSP photon gas well controlled at all times. The x-axis is shown in units
of $\rho_{\gamma}\langle v_{\perp}\rangle$, and the y-axis is shown in units
of benchmark thermalization time $\tau_{*}(T)$. We illustrate the behavior for
primary modes, partner modes, and select high helicities chosen for
illustrative purposes. It can be seen that: a) All modes have a characteristic
energy
<EMAIL_ADDRESS>2.33333pt}(h)$ given by (21), at which it thermalizes fastest (given by (22))
- these can be read off at the minima for each mode in the figure b) Higher
helicities have lower characteristic energies, and take increasingly longer to
populate their phase space even at these characteristic energies. Note that
$h=\pm 10,000$ modes have already increased their $\tau_{h}(E)/\tau_{*}(T)$ to
$\mathcal{O}(10^{104})$ (beyond the y-axis cutoff in figure) at less than one
order of magnitude in energy above their characteristic energy. The dashed
lines denote the specific energies we used to illustrate mode thermalization
hierarchies in Figure 1 and Figure 4. In this figure, we include the hierarchy
across energies to complete the picture of “double hierarchy” in UV
thermalization and “weak hierarchy” in IR thermalization. The parameter
choices for this illustration are: Temperature $T=10^{4}\rho_{\gamma}$, and
$\langle v_{\perp}\rangle=0.1$.

This thermalization behavior means that each mode achieves equilibrium number
density at its characteristic energy first. At any given time $t$, we can
identify the lowest energy $E_{h\vee}(t)$, and the highest energy
$E_{h\wedge}(t)$, between which a mode has populated its equilibrium phase
space. At times
<EMAIL_ADDRESS>2.33333pt}(h))$,
<EMAIL_ADDRESS>2.33333pt}(h)$, which is as yet not fully populated. As time evolves above
<EMAIL_ADDRESS>2.33333pt}(h))$, $E_{h\wedge}(t)$ is given by (14) and $E_{h\vee}(t)$ by (24).
Putting these together, the complete equation for the differential mode number
density at any time is given by:

$\displaystyle n_{h}(E,t)\lessapprox\begin{cases}\text{At \leavevmode\nobreak\
time\leavevmode\nobreak\
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}:\\\ n_{h}^{(eq)}(E)\frac{t}{\tau_{h}(E)}&\text{for
all}\leavevmode\nobreak\ \leavevmode\nobreak\ E\\\ \\\ \text{At
\leavevmode\nobreak\ time\leavevmode\nobreak\
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}:\\\ n_{h}^{(eq)}(E)\frac{E}{E_{h\vee}(t)}&0\leq E\leq
E_{h\vee}(t)\\\ n_{h}^{(eq)}(E)&E_{h\vee}(t)\leq E\leq E_{h\wedge}(t)\\\
n_{h}^{(eq)}(E)(\frac{E_{h\wedge}(t)}{E})^{2(\tilde{h}+1)}&E_{h\wedge}(t)\leq
E<\infty\\\ \end{cases}$ (28)

####  V.1.2 The complete picture of mode thermalization hierarchies

While each mode follows this behavior, there are multiple hierarchies in the
thermalization across modes that govern how the CSP gas as a whole behaves.
The following hierarchies were discussed in III and IV:

  1. (i)

The mode characteristic energy, given by (21), is monotonic in $\tilde{h}$,
with
<EMAIL_ADDRESS>2.33333pt}(h)\rightarrow 0$ for $\tilde{h}\rightarrow\infty$ and
<EMAIL_ADDRESS>2.33333pt}(h)\rightarrow\rho_{\gamma}\langle v_{\perp}\rangle$ for
$\tilde{h}\rightarrow 0$. This means modes with higher helicities first
populate their phase space at smaller energies.

  2. (ii)

The mode thermalization time at its characteristic energy, given by (22),
which gives the shortest time taken by a mode to achieve equilibrium density
at any region in its phase space, is also monotonic in $\tilde{h}$, with
<EMAIL_ADDRESS>2.33333pt}(h)\bigr{)}\rightarrow\infty$ as $\tilde{h}\rightarrow\infty$. This
means modes with higher helicities take parametrically longer to even populate
the region of phase space close to their (already smaller) characteristic
energies.

  3. (iii)

In the correspondence domain, all modes have a “double hierarchy” in their
characteristic thermalization times relative to the primary mode per (8), due
to cross sections suppressed via $\biggl{(}\frac{E}{\rho_{\gamma}\langle
v_{\perp}\rangle}\biggr{)}^{2\tilde{h}}$ as well as
$2^{\tilde{h}}(\tilde{h}+1)!^{2}$. That is, higher helicities find it super-
exponentially difficult to populate any region of their phase space, with the
difficulty also increasing polynomially with energy. This hierarchy is
reflected in the behavior of $E_{h\wedge}(t)$, which grows as
$t^{\frac{1}{2(\tilde{h}+1)}}$ per (14).

  4. (iv)

In the deviation domain, characteristic thermalization times follow a much
weaker hierarchy, given by (20).

The mode thermalization hierarchies at energies much greater than and much
less than the spin scale ((iii) and (iv) above) are illustrated in Figures 1
and 4 respectively. The hierarchies in mode characteristic energies and in the
mode thermalization times at characteristic energy ((i) and (ii) above) are
evident in Figure 6.

###  V.2 Overall CSP behavior

We first present the analyses estimating the effective relativistic degrees of
freedom $g_{CSP}(t)$. We then synthesize the salient aspects of CSP
thermalization behavior, energy density and spectrum of thermal radiation
discussed previously across sections III.2, III.3 and IV.4.

####  V.2.1 Effective relativistic degrees of freedom

The total number density of the CSP gas is given by the sum over mode
densities:

$\displaystyle n_{CSP}(t)$
$\displaystyle=\sum_{h}\int_{0}^{\infty}dEn_{h}(E,t)$ (29)
$\displaystyle\equiv g_{CSP}(t)\frac{\zeta(3)}{\pi^{2}}T^{3}$ (30)

where we separate the familiar form for a thermalized Bose gas [35] to define
the effective internal degrees of freedom for the CSP gas. $g_{CSP}(t)$
accounts for the fractional thermalization of every CSP mode and is a useful
state variable of this system. Since the CSP gas is always partially thermal,
it should be immediately apparent that $g_{CSP}(t)\lll$ the naive guess
$\infty$, at any finite time. Since $g_{CSP}(t)$ depends on how much of its
phase space each mode has populated by the time $t$, it implicitly depends on
$E_{h\wedge}(t)$ and $E_{h\vee}(t)$ of every CSP mode.

Given the hierarchical, controlled thermalization behavior of the CSP at all
energies, we can model the CSP as the familiar photon, but receiving time-
dependent corrections to its relativistic degrees of freedom:

$g_{CSP}(t)\approx[2-\Delta(t)]+\delta g(t)$ (31)

where the $2-\Delta(t)$ degrees of freedom correspond to the primary modes,
and $\delta g(t)$ accounts for the correction from partial thermalization of
all other modes. We estimate these effects using (28):

$\begin{split}g_{CSP}(t)\approx\sum_{h}\frac{1}{\zeta(3)}&\biggl{[}\mathrm{Li}_{3}(e^{-\frac{E_{h\vee}(t)}{T}})\\\
&-\mathrm{Li}_{3}(e^{-\frac{E_{h\wedge}(t)}{T}})\\\
&+\frac{E_{h\vee}(t)}{T}\mathrm{Li}_{2}(e^{-\frac{E_{h\vee}(t)}{T}})\\\
&-\frac{E_{h\wedge}(t)}{T}\mathrm{Li}_{2}(e^{-\frac{E_{h\wedge}(t)}{T}})\\\
&+\frac{E_{h\vee}(t)^{2}}{2T^{2}}\mathrm{Li}_{1}(e^{-\frac{E_{h\vee}(t)}{T}})\\\
&-\frac{E_{h\wedge}(t)^{2}}{2T^{2}}\mathrm{Li}_{1}(e^{-\frac{E_{h\wedge}(t)}{T}})\biggr{]}\\\
\end{split}$ (32)

Comparing with (31) and noting that the primary modes will thermalize fastest,
we recognize that at times $t\gg\tau_{*}(T)$, the two primary modes will
contribute 2 full degrees of freedom, but with very small deviations due to
their $E_{h\vee}(t)>0$ at all finite times:

$\begin{split}g_{\pm
1}(t)\equiv[2-\Delta(t)]\approx\frac{2}{\zeta(3)}\biggl{[}&\mathrm{Li}_{3}(e^{-\frac{E_{h\vee}(t)}{T}})\\\
&+\frac{E_{h\vee}(t)}{T}\mathrm{Li}_{2}(e^{-\frac{E_{h\vee}(t)}{T}})\\\
&+\frac{E_{h\vee}(t)^{2}}{2T^{2}}\mathrm{Li}_{1}(e^{-\frac{E_{h\vee}(t)}{T}})\biggr{]}\\\
\end{split}$ (33)

Note that time- and $\rho_{\gamma}$\- dependence in (32) and (33) can be made
explicit by using the equations for $E_{h\vee}(t)$ in (24) and
$E_{h\wedge}(t)$ in (14). As $E_{h\vee}(t)\rightarrow 0$ and
$E_{h\wedge}(t)\rightarrow\infty$, a full degree of freedom is unlocked from a
mode. We will consider these conditions as met when $E_{h\vee}(t)\ll T$ and
$E_{h\wedge}(t)\gg T$.

For familiar thermodynamic systems that have $T\gg\rho_{\gamma}$,
$E_{h\vee}(t)$ is always $\ll T$, so the behavior of the relativistic degrees
of freedom tracks the evolution of $E_{h\wedge}(t)$ only. Specifically, it is
driven by the strong mode thermalization hierarchies at energies much greater
than the spin scale. For such thermal systems, the primary modes will be
contributing 2 full degrees of freedom at all $t\gg\tau_{*}(T)$, and any
deviations from this are negligible until the partner modes start appreciably
thermalizing several time orders of magnitude after $\tau_{*}(T)$. Ultra cold
thermal systems (that have $T\sim\rho_{\gamma}\langle v_{\perp}\rangle$) are
likely to see detectable deviations in $g_{CSP}(t)$ and energy density in
short times, and require further investigation.

####  V.2.2 Synopsis of CSP thermalization behavior,
energy density and spectrum of thermal radiation

The essence of the effects of a small non-zero $\rho_{\gamma}$ on photon
thermalization is two-fold:

  1. (i)

At energies above the spin scale (“correspondence domain”), we recover
familiar thermodynamic behavior with $\rho_{\gamma}$-dependent corrections to
all thermodynamic quantities that become manifest as the system is allowed to
thermalize over long periods of time. Specifically, in an isothermal system,
we get corrections to internal energy, thermal power spectrum and relativistic
degrees of freedom that are tell-tale signs of a CSP, but such deviations are
observable only with exponentially long thermalization time scales. All
thermodynamically relevant quantities remain finite at all finite times, and
the rate of increase in CSP energy density is inverse in time.

  2. (ii)

We unlock an entirely new region of phase space at energies less than the spin
scale (“deviation domain”), with novel behavior from all helicities. The CSP
gas is populated nearly identically by increasing number of modes at
progressively lower energies, but with thermalization timescales that also
increase in tandem, keeping the IR phase space of the CSP gas well-behaved at
all finite times. The total energy density in the deviation domain remains
sub-dominant, but the power radiated at these very low frequencies shows large
fractional deviations from QED, with calculable deviations in spectrum shape
as well.

Since spin scale has mass dimensions, a thermodynamic system with non-zero
$\rho_{\gamma}$ has a new natural scale (in addition to temperature). The
interplay of these two energy scales sets the overall thermal behavior of the
CSP gas. Specifically, for a system with $T\gg\rho_{\gamma}$, the dominant
behavior of the CSP gas follows (i) above, whereas for a system with
$T\sim\rho_{\gamma}\langle v_{\perp}\rangle$, we get entirely novel behavior
overall, dominated by (ii) above.

Since the spin scale is expected to be small in nature, the thermalization
behavior of the CSP photon is identical to that of the ordinary photon for
most familiar thermodynamic systems (which will have $T\gg\rho_{\gamma}$). In
such systems, any deviations in thermodynamic behavior due to the familiar
photon potentially having a non-zero spin scale will be manifest only with
very long thermalization timescales (assuming $\rho_{\gamma}$ non-zero but not
so small that it evades detectability in the age of the universe) and/or in
the deviation domain behavior (assuming $\rho_{\gamma}$ large enough to be
detectable with available low energy technologies).

##  VI Open Questions

In this section, we briefly discuss the open questions that arise from the
thermodynamic study of CSP photons presented in this paper. Resolving these is
beyond the scope of this work, but could inspire aspects of our future study.

We expect the spin scale of the CSP photon $\rho_{\gamma}$ to be small in
nature. If this weren’t the case, the CSP photon gas would rapidly thermalize
its partner and sub-partner modes, with rapid increase in its internal energy
and degrees of freedom - essentially acting as a system with a high heat
capacity that grows discernibly with time. In an isothermal bath, such as the
one in early universe shortly before recombination, a photon gas with
$\rho_{\gamma}$ much larger than $meV$-scale would have exhibited a
fundamentally different behavior from what has been observed. For
$\rho_{\gamma}$ smaller than $meV$-scale, we expect only small departures from
standard thermal behavior, but it would be interesting to study how best to
use early universe thermal signatures to constrain the spin-scale of the
photon.

In this paper, we investigated an isothermal system, with unbounded energy
available to thermalize the CSP gas. Even with such a set up, we saw that the
CSP gas takes infinite time to access that energy. Despite possessing infinite
internal degrees of freedom in principle, a thermodynamic system with a CSP
photon gas cannot access those degrees of freedom in any finite time - not
even when supplied with infinite energy to do so. Now, let us consider a more
physically realistic thermodynamic system. If we relax the isothermal
assumption and supply fixed energy to the CSP gas, for instance. The primary
modes of the CSP photon gas will still thermalize rapidly. The other modes
still thermalize on exponentially longer time scales, but we can expect that
some of the energy increase in modes that thermalize later will come from
leakage of energy from the already thermalized modes. In such thermodynamic
systems, the total energy in the CSP gas will be finite and bounded by the
available energy, even with infinite time. We can expect the gas to slowly
lower its overall temperature and increase its entropy as more modes
thermalize. CSP thermalization behavior in this and other thermodynamic
situations requires additional investigation.

Finally, ultra-cold thermal systems (those that have
$T\lesssim\rho_{\gamma}\langle v_{\perp}\rangle$) would provide a potentially
interesting regime in which to study signals of non-zero $\rho_{\gamma}$. To
study such systems completely, we need to include Bose condensation effects
and work with full scattering amplitudes, not just soft limits. This motivates
further study of CSP physics (in QED) at low temperatures in future work.

###### Acknowledgements.

We thank Javier Acevedo, Lance Dixon, Saarik Kalia, Aidan Reilly and Kevin
Zhou for useful discussions over the course of this work. The authors were
supported by the U.S. Department of Energy under contract number DE-
AC02-76SF00515 at SLAC.

##  Appendix A Phase space evolution of the thermalizing CSP photon gas

The microscopic evolution of the phase space distribution of every mode $h$ of
the photon gas is governed by the Boltzmann equation [35]:

$\hat{\textbf{L}}[f_{h}]=\textbf{C}[f_{h}]$ (34)

where C is the collision operator and $\hat{\textbf{L}}$ is the Liouville
operator. The covariant, relativistic Liouville operator is [35]:

$\hat{\textbf{L}}=p^{\mu}\frac{\partial}{\partial
x^{\mu}}-\Gamma^{\mu}_{\nu\sigma}p^{\nu}p^{\sigma}\frac{\partial}{\partial
p^{\mu}}$ (35)

In a Minkowski background, the Liouville operator simplifies to:

$\hat{\textbf{L}}[f_{h}(E,t)]=E\frac{\partial}{\partial t}f_{h}(E,t)$ (36)

The collision term for the process $a+b+...\longleftrightarrow
i+j+\gamma_{h}+...$ is given by [35]:

$\displaystyle\hat{\textbf{C}}[$ $\displaystyle f_{h}(E,t)]=\int
d\Pi_{a}d\Pi_{b}....\leavevmode\nobreak\ d\Pi_{i}d\Pi_{j}....$
$\displaystyle\times(2\pi)^{4}\delta^{4}(\Sigma p_{in}-\Sigma
p_{out+\gamma_{h}})$
$\displaystyle\times\bigg{[}|\mathcal{M}|^{2}_{a+b+..\rightarrow
i+j+\gamma_{h}+..}f_{a}f_{b}...(1\pm f_{i})(1\pm f_{j})(1\pm f_{h})..$
$\displaystyle-|\mathcal{M}|^{2}_{i+j+\gamma_{h}+..\rightarrow
a+b+..}f_{i}f_{j}f_{h}..(1\pm f_{a})(1\pm f_{b})..\bigg{]}$ (37)

where ‘in’ denotes the incoming scatterers, ‘out’ denotes all outgoing
particles except the CSP mode $\gamma_{h}$ that we are interested in,
$f_{\psi}$ refers to the phase space distribution function of particle $\psi$,
$d\Pi\equiv\frac{g}{(2\pi)^{3}}\frac{d^{3}p}{2E}$, and $(1\pm f_{\psi})$
factors capture Bose enhancement/ Fermi suppression respectively. We invoke
time invariance to assume $|\mathcal{M}|^{2}_{a+b+..\rightarrow
i+j+\gamma_{h}+..}=|\mathcal{M}|^{2}_{i+j+\gamma_{h}+..\rightarrow
a+b+..}\equiv|\mathcal{M}|^{2}$. Using the energy conservation part of the
delta function, and assuming chemical and kinetic equilibrium of all other
species, (A) simplifies to:

$\begin{split}\hat{\textbf{C}}[f_{h}&(E,t)]=\int
d\Pi_{a}d\Pi_{b}\leavevmode\nobreak\ ...\leavevmode\nobreak\
d\Pi_{i}d\Pi_{j}\leavevmode\nobreak\ ...\\\ &\times(2\pi)^{4}\delta^{4}(\Sigma
p_{in}-\Sigma p_{out+\gamma_{h}})|\mathcal{M}|^{2}\\\ &\times
f_{a}f_{b}...(1\pm f_{i})(1\pm
f_{j})...\bigg{[}1+f_{h}[1-\exp{\frac{E}{T}}]\bigg{]}\end{split}$ (38)

where T is the temperature. Note that while all equations upto (A) are
explicitly Lorentz covariant, (38) is not, since we single out the frame of
reference in which the temperature T is a monopole to do thermodynamics. This
is the frame in which we will specify all particle distribution functions.

Recognizing that in equilibrium, $f_{h}$ will follow Bose-Einstein statistics,
we can use $f_{h}^{(eq)}(E)=\frac{1}{\exp{\frac{E}{T}}-1}$ to rewrite the
$\bigg{[}1+f_{h}[1-\exp{\frac{E}{T}}]\bigg{]}$ factor in (38) as
$\bigg{[}1-\frac{f_{h}(E,t)}{f_{h}^{(eq)}(E)}\bigg{]}$. Using (38) and (36),
the Boltzmann equation (34) reduces to a differential equation that can be
solved under the assumptions already stated, to give:

$f_{h}(E,t)=f_{h}^{(eq)}(E)[1-\exp{(-t/\tau_{h}(E))}]$ (39)

where $\tau_{h}(E)$ can be recognized as:

$\begin{split}\tau_{h}(E)=&f_{h}^{(eq)}(E)\bigg{[}\int d\Pi_{in}f_{in}(1\pm
f_{i})(1\pm f_{j})...d\Pi_{out}\\\ &\times(2\pi)^{4}\delta^{4}(\Sigma
p_{in}-\Sigma
p_{out+\gamma_{h}})|\mathcal{M}|^{2}\frac{1}{E}\bigg{]}^{-1}\end{split}$ (40)

Ignoring the Bose enhancement/ Fermi suppression factors from the outgoing
states, i.e., taking $(1\pm f_{\psi})\approx 1$, we get (3). When we consider
a multi-photon scattering process, this assumption means that we ignore any
Bose enhancement to the production of partner and sub-partner modes from the
faster thermalization of primary modes. We make the reasonable assumption that
multi-photon CSP processes are sub-dominant to single photon scattering
processes, as in familiar QED. Additionally, familiar QED has the same IR
divergence from Bose-Einstein statistics - so this aspect is not unique to CSP
photons. Hence, we find it best to keep the Bose enhancement aside and focus
only on aspects of thermal behavior arising from a non-zero spin scale.

##  Appendix B Modifications for a non-relativistic thermal scatterer

In the main paper, we used the average velocity $\langle v_{\perp}\rangle$ in
equations, and used Taylor expansion and/or asymptotic forms of
$J_{h}(\frac{\rho_{\gamma}v_{\perp}}{E})$ where it was valid to do so. These
simplifications were made mainly to aid intuitive understanding of CSP
behavior. This appendix provides some details on modifications needed to some
of these simplified equations when the scatterer is non-relativistic. Note
that all numerical simulations presented in the paper used the full thermal
distribution of $v_{\perp}$ as well as the full Bessel function form of the
soft limit scattering amplitudes (without any simplifications).

Since scattering cross sections have a velocity dependence in the soft limit
per (5), (6) and (7), the velocity distribution of a non-relativistic
scatterer needs to be accounted for when calculating the mode thermalization
times. Equation (3) includes calculations with the following form, using a
Boltzmann distribution for the non-relativistic scatterer:

$\tau_{h}(E)\supset\int_{0}^{1}dv_{\perp}v_{\perp}\exp\biggl{(}{-\frac{v_{\perp}^{2}}{2\langle
v_{\perp}\rangle^{2}}}\biggr{)}\leavevmode\nobreak\
\bigg{|}J_{h}(\frac{\rho_{\gamma}v_{\perp}}{E})\bigg{|}^{2}$ (41)

We can approximate $J_{\alpha}(x)$ with the first term in its Taylor expansion
[34, (10.2.2)] when:

$\displaystyle J_{\alpha}(x)$
$\displaystyle=\sum_{m=0}^{\infty}\frac{(-1)^{m}(\frac{x}{2})^{2m+\alpha}}{m!\leavevmode\nobreak\
\Gamma(m+\alpha+1)}$ (42)
$\displaystyle\approx\frac{(\frac{x}{2})^{\alpha}}{\Gamma(\alpha+1)}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{when}\leavevmode\nobreak\ {x\ll 2\sqrt{\alpha+1}}$
(43)

When the simplifying condition in (43) is valid for the entire range of the
integration in (41), the thermal averaging will give a lower incomplete gamma
function $\gamma(\tilde{h}+1,\frac{1}{2\langle v_{\perp}\rangle^{2}})$, which
is $\approx(\tilde{h}+1)!$ only if $\tilde{h}+1<\frac{1}{2\langle
v_{\perp}\rangle^{2}}$. This means that when working with the simplified forms
of the Bessel scattering cross sections, we need to be careful about the range
of validity of the simplifications. The rest of this appendix discusses the
modifications that need to be made for certain equations in the main paper to
account for the effect of thermal distribution of non-relativistic scatterer
velocities.

Equation (8) gets modified as:

$\frac{\tau_{h}(E)}{\tau_{*}(E)}\sim\frac{\bigl{\langle}|z|^{2}\bigr{\rangle}}{\bigl{\langle}|zF_{h}(\rho_{\gamma}z)|^{2}\bigr{\rangle}}\;\approx\;2^{\tilde{h}}(\tilde{h}+1)!^{(2-\delta)}\biggl{(}\frac{E}{\rho_{\gamma}\langle
v_{\perp}\rangle}\biggr{)}^{2\tilde{h}}$ (44)

where the continuous parameter $\delta$ accounts for the weaker suppression of
smaller helicities due to the thermal distribution in scatterer velocities
when we consider a non-relativistic scatterer. $\delta=0$ for all helicities
when the scatterer is relativistic. When we consider non-relativistic
scatterers, $0\leq\delta\leq 1$. Smaller helicities will see a shorter time
vs. that obtained using $\langle v_{\perp}\rangle$ in lieu of a full thermal
distribution of $v_{\perp}$ i.e., $\delta\rightarrow 1$ for only those modes
with $\tilde{h}<\frac{1}{2\langle v_{\perp}\rangle^{2}}$.

Equation (III.3) requires a minor modification with $(\tilde{h}+1)!^{2}$ in
the denominator replaced with $\tilde{h}+1)!^{(2-\delta)}$. This manifests as
a change in the helicity scaling of $\frac{d}{dt}\mathcal{E}_{CSP}(t)$, which
we discussed in the main paper as following $\tilde{h}^{-4}$, where 3 powers
of $\tilde{h}$ came from $(\tilde{h}+1)!^{\frac{3}{(\tilde{h}+1)}}$ and 1 came
directly from the $2\tilde{h}$ in the denominator. Including the $\delta$
parameter, the helicity scaling varies as $\tilde{h}^{2.5-4}$ for a non-
relativistic scatterer. Low helicities in the non-relativistic limit have
$\delta\approx 1$ as above and $(\tilde{h}+1)!^{\frac{1.5}{(\tilde{h}+1)}}$
grows as $(\tilde{h}+1)^{1.5}$. High helicities have $\delta=0$ as explained
above and $(\tilde{h}+1)!^{\frac{3}{(\tilde{h}+1)}}$ grows as
$(\tilde{h}+1)^{3}$. In the relativistic limit, all modes get suppressed with
$(\tilde{h}+1)^{3}$. These modifications do not change anything significant
about $\mathcal{E}_{CSP}(t)$ since its behavior is controlled by other aspects
of (III.3). The only key property of $\mathcal{E}_{CSP}(t)$ that depended on
the helicity scaling is the sum convergence, and this continues to hold with
the modified scaling for low helicities.

Equation (19), which directly followed from (III.3), gets the same
modifications discussed in the previous paragraph. For a relativistic
scatterer, the power law dependence for all modes has $(\tilde{h}+1)^{3}$ in
the denominator. For a non-relativistic scatterer, this gets modified to be
$(\tilde{h}+1)$ for low helicities with $\tilde{h}<\frac{1}{2\langle
v_{\perp}\rangle^{2}}$ and $(\tilde{h}+1)^{3}$ for higher helicities.

Note that equation (20), which is valid for $E\ll\rho_{\gamma}\langle
v_{\perp}\rangle$ already implicitly accounted for the ranges in $\delta$.

##  Appendix C Time evolution of CSP energy density

The energy density in the CSP gas can be obtained using [35]:

$\mathcal{E}_{CSP}(t)=\sum_{h}\int_{0}^{\infty}dEn_{h}(E,t).E$ (45)

Using (28), we can write (45) as:

$\begin{split}\mathcal{E}_{CSP}(t)\lessapprox\sum_{h}\bigg{[}&\int_{0}^{E_{h\wedge}(t)}dEn_{h}^{(eq)}(E)E\\\
+&\int_{E_{h\wedge}(t)}^{\infty}dEn_{h}^{(eq)}(E)E\frac{t}{\tau_{h}(E)}\bigg{]}\end{split}$
(46)

In writing (46), we have used the following: i) Overall energy density in any
mode is most sensitive to the highest occupied energy, and comparatively
insensitive to the lowest occupied energy, allowing us to take
$E_{h\vee}(t)\rightarrow 0$ ii) $E_{h\wedge}(t)$ tracks (14), with
characteristic thermalization time tracking (8) since $E_{h\wedge}(t)\geq
<EMAIL_ADDRESS>2.33333pt}(h)$ always.

Taking the time derivative of (46), and noting that $n_{h}(E)$ is continuous
at all energies including $E_{h\wedge}(t)$, we get:

$\frac{d}{dt}\mathcal{E}_{CSP}(t)\lessapprox\sum_{h}\int_{E_{h\wedge}(t)}^{\infty}dEn_{h}^{(eq)}(E)\frac{E}{\tau_{h}(E)}$
(47)

To make (47) tractable for further study, we approximate $f_{h}^{(eq)}(E)$ as
$\frac{T}{E}\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for
all}\leavevmode\nobreak\ E\leq T$ and as
$\exp(-\frac{E}{T})\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for
all}\leavevmode\nobreak\ E>T$. We study the two regimes separately, and write:

$\begin{split}\frac{d}{dt}\mathcal{E}_{CSP}(t)\lesssim\frac{1}{2\pi^{2}}\bigg{[}&\sum\limits_{h\,:\,\begin{subarray}{c}\tau_{h}(T)\geq
t\\\ h\neq\pm
1\end{subarray}}\int_{E_{h\wedge}(t)}^{\infty}dE\frac{TE^{2}}{\tau_{h}(E)}\\\
+&\sum\limits_{h\,:\,\begin{subarray}{c}\tau_{h}(T)<t\\\ h\neq\pm
1\end{subarray}}\int_{E_{h\wedge}(t)}^{\infty}dEe^{-\frac{E}{T}}\frac{E^{3}}{\tau_{h}(E)}\bigg{]}\end{split}$
(48)

In (48), the former sum is taken over modes that are still thermalizing, given
by the condition $\tau_{h}(T)\geq t$. The latter sum is taken over modes that
have already thermalized their phase space upto energies $E>T$, given by the
criterion $\tau_{h}(T)<t$. We exclude the primary modes from these sums since
we assume them to be thermal in this analysis, and will use $\tau_{*}(T)$ as
the benchmark time in the next step. Substituting (8) in (48), the integration
in both sums can be exactly done.

$\begin{split}\frac{d}{dt}\mathcal{E}_{CSP}(t)\lesssim&\sum\limits_{h\,:\,\begin{subarray}{c}\tau_{h}(T)\geq
t\\\ h\neq\pm
1\end{subarray}}\frac{T^{4}}{2\pi^{2}[2\tilde{h}-1]}t^{-1}\left(\frac{E_{h\wedge}(t)}{T}\right)^{3}\\\
+&\sum\limits_{h\,:\,\begin{subarray}{c}\tau_{h}(T)<t\\\ h\neq\pm
1\end{subarray}}\bigg{[}\frac{T^{4}}{2\pi^{2}}t^{-1}\left(\frac{E_{h\wedge}(t)}{T}\right)^{2(\tilde{h}+1)}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\times\Gamma\big{(}2-2\tilde{h},\frac{E_{h\wedge}(t)}{T}\big{)}\bigg{]}\end{split}$
(49)

In (49), $\Gamma(x,y)$ is the upper incomplete gamma function. Since the
second sum is taken over modes that have $E_{h\wedge}(t)>T$, the gamma
function falls exponentially or faster with time for $\tilde{h}\geq 1$. So, as
long as there are at least some modes with $\tau_{h}(T)\geq t$ i.e., at all
finite times, the dual sum in (49) is always dominated by the sum over
thermalizing modes. The time dependence in (49) can be made explicit with
(14). Doing so, and dropping the sum over thermal modes that is sub-dominant
gives the expression in (III.3).

## References

  * Wigner [1939] E. P. Wigner, On Unitary Representations of the Inhomogeneous Lorentz Group, Annals of Mathematics 40, 149 (1939).
  * Schuster and Toro [2013a] P. Schuster and N. Toro, On the theory of continuous-spin particles: wavefunctions and soft-factor scattering amplitudes, JHEP 09, 104, arXiv:1302.1198 [hep-th] .
  * Schuster and Toro [2013b] P. Schuster and N. Toro, On the theory of continuous-spin particles: helicity correspondence in radiation and forces, JHEP 09, 105, arXiv:1302.1577 [hep-th] .
  * Schuster and Toro [2015a] P. Schuster and N. Toro, Continuous-spin particle field theory with helicity correspondence, Phys. Rev. D 91, 025023 (2015a), arXiv:1404.0675 [hep-th] .
  * Schuster  _et al._ [2023] P. Schuster, N. Toro, and K. Zhou, Interactions of Particles with ”Continuous Spin” Fields, JHEP 04, 010, arXiv:2303.04816 [hep-th] .
  * Schuster and Toro [2024] P. Schuster and N. Toro, Quantum Electrodynamics Mediated by a Photon with Generalized (Continuous) Spin, Accepted to Phys. Rev. D (2024), arXiv:2308.16218 [hep-th] .
  * Schuster and Toro [2015b] P. Schuster and N. Toro, A new class of particle in 2 + 1 dimensions, Phys. Lett. B 743, 224 (2015b), arXiv:1404.1076 [hep-th] .
  * Bekaert  _et al._ [2016] X. Bekaert, M. Najafizadeh, and M. R. Setare, A gauge field theory of fermionic continuous-spin particles, Phys. Lett. B 760, 320 (2016), arXiv:1506.00973 [hep-th] .
  * Najafizadeh [2020] M. Najafizadeh, Supersymmetric continuous spin gauge theory, JHEP 03, 027, arXiv:1912.12310 [hep-th] .
  * Najafizadeh [2022] M. Najafizadeh, Off-shell supersymmetric continuous spin gauge theory, JHEP 02, 038, arXiv:2112.10178 [hep-th] .
  * Metsaev [2018a] R. R. Metsaev, BRST-BV approach to continuous-spin field, Phys. Lett. B 781, 568 (2018a), arXiv:1803.08421 [hep-th] .
  * Buchbinder  _et al._ [2018] I. L. Buchbinder, V. A. Krykhtin, and H. Takata, BRST approach to Lagrangian construction for bosonic continuous spin field, Phys. Lett. B 785, 315 (2018), arXiv:1806.01640 [hep-th] .
  * Alkalaev  _et al._ [2018] K. Alkalaev, A. Chekmenev, and M. Grigoriev, Unified formulation for helicity and continuous spin fermionic fields, JHEP 11, 050, arXiv:1808.09385 [hep-th] .
  * Buchbinder  _et al._ [2020a] I. L. Buchbinder, S. Fedoruk, A. P. Isaev, and V. A. Krykhtin, Towards Lagrangian construction for infinite half-integer spin field, Nucl. Phys. B 958, 115114 (2020a), arXiv:2005.07085 [hep-th] .
  * Buchbinder  _et al._ [2019a] I. L. Buchbinder, S. J. Gates, and K. Koutrolikos, Superfield continuous spin equations of motion, Phys. Lett. B 793, 445 (2019a), arXiv:1903.08631 [hep-th] .
  * Buchbinder  _et al._ [2019b] I. L. Buchbinder, M. V. Khabarov, T. V. Snegirev, and Y. M. Zinoviev, Lagrangian formulation for the infinite spin $N$=1 supermultiplets in $d$=4, Nucl. Phys. B 946, 114717 (2019b), arXiv:1904.05580 [hep-th] .
  * Buchbinder  _et al._ [2020b] I. L. Buchbinder, S. Fedoruk, and A. P. Isaev, Massless infinite spin (super)particles and fields, Proc. Steklov Inst. Math. 309, 46 (2020b), arXiv:1911.00362 [hep-th] .
  * Buchbinder  _et al._ [2022] I. L. Buchbinder, S. A. Fedoruk, A. P. Isaev, and V. A. Krykhtin, On the off-shell superfield Lagrangian formulation of 4$D$, $N$=1 supersymmetric infinite spin theory, Phys. Lett. B 829, 137139 (2022), arXiv:2203.12904 [hep-th] .
  * Zinoviev [2017] Y. M. Zinoviev, Infinite spin fields in $d$ = 3 and beyond, Universe 3, 63 (2017), arXiv:1707.08832 [hep-th] .
  * Alkalaev and Grigoriev [2018] K. B. Alkalaev and M. A. Grigoriev, Continuous spin fields of mixed-symmetry type, JHEP 03, 030, arXiv:1712.02317 [hep-th] .
  * Burdík  _et al._ [2020] C. Burdík, V. K. Pandey, and A. Reshetnyak, BRST–BFV and BRST–BV descriptions for bosonic fields with continuous spin on $R^{1,d-1}$, Int. J. Mod. Phys. A 35, 2050154 (2020), arXiv:1906.02585 [hep-th] .
  * Metsaev [2017a] R. R. Metsaev, Continuous spin gauge field in (A)dS space, Phys. Lett. B 767, 458 (2017a), arXiv:1610.00657 [hep-th] .
  * Metsaev [2017b] R. R. Metsaev, Fermionic continuous spin gauge field in (A)dS space, Phys. Lett. B 773, 135 (2017b), arXiv:1703.05780 [hep-th] .
  * Metsaev [2018b] R. R. Metsaev, Continuous-spin mixed-symmetry fields in AdS(5), J. Phys. A 51, 215401 (2018b), arXiv:1711.11007 [hep-th] .
  * Khabarov and Zinoviev [2018] M. V. Khabarov and Y. M. Zinoviev, Infinite (continuous) spin fields in the frame-like formalism, Nucl. Phys. B 928, 182 (2018), arXiv:1711.08223 [hep-th] .
  * Metsaev [2019] R. R. Metsaev, Light-cone continuous-spin field in AdS space, Phys. Lett. B 793, 134 (2019), arXiv:1903.10495 [hep-th] .
  * Metsaev [2021] R. R. Metsaev, Mixed-symmetry continuous-spin fields in flat and AdS spaces, Phys. Lett. B 820, 136497 (2021), arXiv:2105.11281 [hep-th] .
  * Rivelles [2015] V. O. Rivelles, Gauge theory formulations for continuous and higher spin fields, Phys. Rev. D 91, 125035 (2015), arXiv:1408.3576 [hep-th] .
  * Rivelles [2017] V. O. Rivelles, Remarks on a gauge theory for continuous spin particles, Eur. Phys. J. C 77, 433 (2017), arXiv:1607.01316 [hep-th] .
  * Bekaert and Skvortsov [2017] X. Bekaert and E. D. Skvortsov, Elementary particles with continuous spin, Int. J. Mod. Phys. A 32, 1730019 (2017), arXiv:1708.01030 [hep-th] .
  * Najafizadeh [2018] M. Najafizadeh, Modified Wigner equations and continuous spin gauge field, Phys. Rev. D 97, 065009 (2018), arXiv:1708.00827 [hep-th] .
  * Wigner [1963] E. P. Wigner, Invariant Quantum Mechanical Equations of Motion, Theoretical Physics, International Atomic Energy Agency, Vienna (1963).
  * Weinberg [1964] S. Weinberg, Photons and Gravitons in $S$-Matrix Theory: Derivation of Charge Conservation and Equality of Gravitational and Inertial Mass, Phys. Rev. 135, B1049 (1964).
  * [34] DLMF, NIST Digital Library of Mathematical Functions, https://dlmf.nist.gov/, Release 1.2.0 of 2024-03-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Kolb and Turner [1994] E. Kolb and M. Turner,  _The Early Universe_ (Westview Press, 1994).
  * Mecholsky  _et al._ [2021] N. A. Mecholsky, S. Akhbarifa, W. Lutze, M. Brandys, and I. L. Pegg, Bessel function $j_{n}$ maxima and minima, Data in Brief 39, https://doi.org/10.1016/j.dib.2021.107508 (2021).