# Boltzmann Equation Field Theory I: Ensemble Averages

Jun Yan Lau1
1UCL Mullard Space Sciences Laboratory, Holmbury St Mary, Surrey RH5 6NT
E-mail<EMAIL_ADDRESS>
###### Abstract

A quasi-nonlinear field theory which describes how to take ensemble averages
that are unique to the Collisionless Boltzmann Equation is described. The
assumption that the ensemble average of the distribution function is equal to
the extremum entropy state, $\langle f\rangle=f_{0}$ is taken apart and shown
to be wrong. An application describes the nonlinear saturation of Jeans’
instability, and the gravitational amplification of Poisson noise.

###### keywords:

Galaxy: kinematics and dynamics – galaxies: kinematics and dynamics – methods:
analytical

##  1 Introduction

The long-standing study of the dynamics of observed stellar systems began with
E. C. Pickering (1891) who put forth the claim that different globular
clusters (GCs) had similar “laws of distributions”, referring to the
distribution of stars inferred from 2D projected images obtained from the
Draper Catalogue.

Plummer (1911) twenty years later proposed to ”pursue the search for a
physical basis on which the distribution of stars in clusters may be
established”. Under the presumption that spherical nebulae that are in
convective equilibrium produce spherical globular clusters of the same mass
distribution, he proposed a physically meaningful and demonstrably suitable
fit ( within $\pm 5\%$) to Pickering’s data; his Plummer model.

Soon after, Eddington (1916) utilised Abel’s theorem to translate 2D projected
observations of GCs into 3D position-space, a result used by Plummer to
produce his formulae connecting $f({\bf x},{\bf v},t)=f({\bf w},t)$ the 6D
distribution of stars in position and velocity space (henceforth phase space)
to their mass distributions in position-space.

And from thereon, the general methodology of fitting $\\{{\bf w}_{i}\\}$, the
observed phase-space coordinates of a set of stars which are believed to form
a gravitationally bound system with $f$, the distribution function (DF) of
these stars, is as follows: First, intuit the symmetries of this observed
system. Secondly, compute a distribution function which obeys these
symmetries. Finally, fit this symmetric distribution function to the observed
data.

Examples of this methodology are found in nearly every system: Galaxies were
first fit by the axisymmetric Mestel (1963) disc profile, which bellied the
flat rotation curve found observationally and the Hernquist (1990) profile
approximates De Vacouleur’s law for the surface brightness of elliptical
galaxies. In the realm of globular clusters, Henon (1959)’s isochrone
distribution is the most general distribution of stars for which the radial
period of orbits is purely a function of their energy. King (1966) championed
the usage of a series of lowered isothermal spheres; systems with finite mass
and extent and with the same velocity dispersion everywhere. Osipkov
(1979)-Merritt (1985) distribution functions describe a simple way to
introduce velocity-space anisotropy to initially isotropic, spherical systems.

These distribution functions are also known as mean-fields, $f({\bf
w},t)=f_{0}({\bf w})$, and share several commonalities:

Firstly, they are static in time. This means that states $f_{0}$ are assumed
to be in equilibrium. Whether this is because equilibrium $f_{0}$ are simply
equal to calculate, or because of the long-standing hypothesis that systems
must be stable to be observed in the present day is up to contention. Just
looking at the extensive literature regarding streams of stars being tidally
torn from globular clusters/dwarf galaxies by the Milky Way, density wave
theory in the context of the development of spirals (see Lin & Shu (1964) for
how it began and Sellwood & Carlberg (2014) for how it is going)/slowing bar
theory (Chiba et al., 2020) is that the assumption that observed systems are
well fit by steady-state DF serves as a zeroth approximation at best.

Secondly, they respect some spatial geometric symmetry. This was dubious even
back in 1911, as Plummer wrote in reference to the fit of his model to an
observed GC, ‘it is clear that the counts are at least as well represented by
the formula we have chosen as they are consistent with the fundamental
hypothesis that we are dealing with a truly spherical distribution’. The
natural asymmetry of astrophysical systems has been known of for over a
century observationally and theoretically. Gravity is an attractive long-
ranged force. Gravitationally stable systems generally exhibit system-scale
fluctuations; observationally we see spiral arms, galactic bars, globular
clusters with well developed dipolar asymmetries, fluctuations which are not
captured by the mean-field symmetries.

I think it is clear that fitting $f_{0}$ to a set of observed particles
$\\{{\bf w}_{i}\\}$ is not the right approach. To put it very simply, it
should be clear that we don’t know what the DF is! We know that there are very
many possible distribution functions which may be fit to any set of observed
particles, and we know we cannot claim certainty in knowing that any single
$f$ fits $\\{{\bf w}_{i}\\}$. So instead of choosing $f=f_{0}$ based on a
physical argument (e.g. maximum entropy, or any of the other reasons conferred
above), we know should be assigning probabilities $P[f]$ to $f$, based on
those physical arguments. This insight—that it is better to accomodate our
ignorance of the DF than to just choose what we think is sensible—is invoked
earlier by Magorrian (2006) who applied it in criticism of maximum likelihood
methods. He showed that the mass of supermassive black holes in toy galaxies
for which mock observations have been obtained and trial DFs have been
inferred via Bayesian inference, are better constrained if one selects for the
mass which fits as many viable trial DFs as possible, and not for the mass
which fits the most viable trial DF possible.

Now it must be noted that Pickering and Plummer conceived of their theories
with only a handful of observations to draw from. Even until the late 1950s,
astrophysicists only had several tens of thousands of stars to draw
measurements from. They had a far simpler task than we do, because their data
was sufficiently coarse to be able to neglect variations in $f$. They had
neither the resolution nor the statistics to make claims regarding the
substructures of astronomical objects, so they did not need to define such a
$P[f]$. We are neither as ignorant (from the statistical mechanical
perspective, not from the intellectual perspective!), nor as unlucky as they
were.

Gaia Collaboration et al. (2022), has revealed a dense network of asymmetries
and large-scale substructure within and without our Galaxy; the Antoja et al.
(2018) spiral, the Gaia-Enceladus sausage, or the newly coined ‘X-shaped,
Peanut, Boxy bar’. The corotation resonances of the bar severing the Hercules
stream (see (Monari et al., 2019)), or the evidence for the slowing bar as
stars trapped within bar resonances (Chiba & Schönrich, 2021). The question of
whether fluctuations matter, not just to the secular evolution of systems but
also to what we see today was pertinent before GAIA, but is crucial to
understand now.

So we need a variational theory in $f$, which allows us to predict the
properties of fluctuations beyond the steady-state distributions: this is a
field theory. In this paper, I will derive such a theory which accounts for
finiteness noise, gravitational correlations, and more, all under the
framework of the CBE!

Section 2 describes essential mathematical tools to understand this paper.
Section 3 describes how fluctuations in $f$ can be sourced from discreteness
noise via an entropy argument. Section 4 describes how we can introduce
desired structures into this noise. Section 5 describes how we can shape
fluctuations into a perturbative field theory. Section 6 describes how we can
reclaim observables using the field theory, Section 7 describes fluctuations
corresponding to a system with a mean energy constraint and Section 8 shows
how we can make predictions. Finally, Sections 9 and 10 summarise and
deliberate over the predictions of this field theory.

##  2 Mathematical background

Here I introduce essential mathematical tools and establish my notation.

###  2.1 Poisson Brackets

The Poisson brackets can be computed as

$[f,g]=\sum_{i}\bigg{(}{\partial f\over\partial q_{i}}{\partial g\over\partial
p_{i}}-{\partial f\over\partial p_{i}}{\partial g\over\partial
q_{i}}\bigg{)}.$ (1)

where the pair of canonical phase-space coordinates $({\bf q},{\bf p})$ obey,

$[p_{i},p_{j}]=0,~{}[q_{i},q_{j}]=0,~{}[q_{i},p_{j}]=\delta_{ij}$ (2)

In astrophysics, there are two important choices of such coordinates, the 3D
Cartesian positions/velocities, ${\bf w}=({\bf x},{\bf v})$, and their
generalisations, the angle/actions, ${\bf
w}=({\mathchoice{\hbox{{\gkvecten\char 18\relax}}}{\hbox{{\gkvecten\char
18\relax}}}{\hbox{{\gkvecseven\char 18\relax}}}{\hbox{{\gkvecseven\char
18\relax}}}},{\bf J})$. The angle-action coordinates generalise positions and
velocities for inhomogeneous systems. The actions ${\bf J}$ are constants of
motion for orbits defined by the flow of the globally integrable Hamiltonian
$H=H({\bf J})$, whereas the angles ` are the conjugates of ${\bf J}$, and thus
obey ${\mathchoice{\hbox{{\gkvecten\char 18\relax}}}{\hbox{{\gkvecten\char
18\relax}}}{\hbox{{\gkvecseven\char 18\relax}}}{\hbox{{\gkvecseven\char
18\relax}}}}(t)={\mathchoice{\hbox{{\gkvecten\char
18\relax}}}{\hbox{{\gkvecten\char 18\relax}}}{\hbox{{\gkvecseven\char
18\relax}}}{\hbox{{\gkvecseven\char 18\relax}}}}_{0}+\mathbf{\Omega}t$,
$\Omega=\partial H/\partial{\bf J}$ by Hamilton’s equation. Due to their
canonicity, the phase-space element may be expressed ${\rm d}^{6}{\bf w}={\rm
d}^{3}{\bf x}{\rm d}^{3}{\bf v}={\rm d}^{3}{\mathchoice{\hbox{{\gkvecten\char
18\relax}}}{\hbox{{\gkvecten\char 18\relax}}}{\hbox{{\gkvecseven\char
18\relax}}}{\hbox{{\gkvecseven\char 18\relax}}}}{\rm d}^{3}{\bf J}$.

While action-angle coordinates provides for a simple way to track how a star
moves along its trajectory, their usage presumes that the Hamiltonian is
globally integrable. A field theory integrates over all distribution
functions, and general Hamiltonians only admit local integrability, (i.e.
different parts of a system have different conserved quantities). Thus, I
choose to use ${\bf w}=({\bf x},{\bf v})$ unless otherwise stated.

###  2.2 CBE

The collisionless Boltzmann equation governs the evolution of the one-particle
distribution function $f({\bf w},t)$;

${{\rm d}f\over{\rm d}t}={\partial f\over\partial t}+[f,H]=0$ (3)

where $H$ is the Hamiltonian of the system, and $f$ is the 1-particle
distribution function. In a self-consistent system (that is a system which
evolves under a force law which mediates inter-particle interactions) the
Hamiltonian takes the form:

$H[f]({\bf w},t)={{\textstyle{1\over 2}}}m{\bf v}^{2}+\Phi[f]({\bf x},t)$ (4)

for self-gravitating systems, the potential $\Phi$ is defined via Poisson’s
equation,

$\Phi({\bf x},t)=-GMm\int{\rm d}^{3}{\bf w}^{\prime}~{}{1\over|{\bf x}-{\bf
x}^{\prime}|}f({\bf w}^{\prime},t)$ (5)

Boltzmann’s equation can be derived from the BBGKY hierarchy (Binney &
Tremaine, 2008) which begins with an N-particle distribution function
$f^{(N)}({\bf w}_{1},{\bf w}_{2},...,{\bf w}_{N},t)$, which is reduced to
1-particle distribution functions $f({\bf w},t)$ under the assumption that no
particle is correlated with any other particle.

In general, $f=f({\bf w},t)$ is a function of phase-space coordinates and
time. I have chosen to exclude the additional dependence on time from my
notation for the sake of brevity, since it is only rarely involved in the
context of time derivatives for which the time-dependence is obvious. Where it
appears, the maximum entropy state $f_{0}({\bf w})$ has no explicit time-
dependence.

In this writing, the usual choice of normalisation for the DF $f$ is

$\int{\rm d}^{6}{\bf w}~{}f({\bf w})=1$ (6)

. Astrophysical systems are usually of finite mass, and thus bear that natural
normalisation. However, exceptions exist in the treatment of systems with
infinite mass: then, the normalisation is usually done with respect to the
density at some scale radii/the spatial density, for homogeneous systems.

###  2.3 Functional Analysis

A functional $G[f]({\bf w})$ is in general, an integral over a function of the
function $f$:

$G[f]({\bf w})=\int{\rm d}^{6}{\bf w}_{a}~{}g(f({\bf w}_{a}),{\bf w}).$ (7)

The functional integral of a functional $G[f]$ with respect to the function
$f({\bf w})$ is formally defined as

$\int^{a}_{b}\mathcal{D}f~{}G[f]=\int^{a}_{b}...\int^{a}_{b}\Pi_{{\bf w}}{\rm
d}f({\bf w})~{}G[f]$ (8)

where the product goes over all points in phase-space. In physics, this
integral is often used to indicate integration over all possible paths or all
possible states, where these paths and/or states are defined by $f$.

A simpler, more intuitive way of arriving at the non-rigorous definition
provided above is to consider first a world in which $f({\bf w})=f_{i}$ is
discretised over some phase-space volume $\Delta$. These phase-space pixels
can individually be varied, and all possible $f$ can thus be represented by
combinations of all possible $f_{i}$.

Then the functional integral is defined as,

$\int^{a}_{b}\mathcal{D}f~{}G[f]=\Bigg{(}\prod_{i}\int^{a}_{b}{\rm
d}f_{i}\Bigg{)}G(f_{1},f_{2},...)$ (9)

and in the limit $\Delta\rightarrow 0$ we reclaim the original result.

The functional derivative obeys the usual rules derivatives do, and is defined
by the following relation:

${\delta f({\bf w}^{\prime})\over\delta f({\bf w})}=\delta^{6}({\bf w}-{\bf
w}^{\prime}).$ (10)

###  2.4 Generalised Equipartition Theorem

I will now derive a more general variant of the equipartition theorem (Tolman,
1938). Consider an abstract Gibbs ensemble in which the probability associated
with a state described by ${\bf s}$ is;

$P({\bf s})={1\over\mathcal{Z}}\exp\Bigg{(}-\frac{1}{2}{\bf s}^{T}M{\bf
s}\Bigg{)}$ (11)

where $M$ is an invertible matrix, summation is implied in the scalar product,
and

$\mathcal{Z}=\int\mathcal{D}{\bf s}~{}\exp(-{\textstyle{1\over 2}}{\bf
s}^{T}M{\bf s})$ (12)

is a normalisation factor where the functional integral over $({\bf
s})_{i}=s_{i}$ goes between $(-\infty,\infty)$. Then ensemble averages with
respect to this probability take the form

$\langle A\rangle=\int\mathcal{D}{\bf s}~{}P({\bf s})A({\bf s}).$ (13)

In this paper, we are concerned with ensemble averages of functionals $A$ that
take the form,

$\langle s_{1}s_{1}s_{3}...s_{n(\mathcal{N})}\rangle=\int\mathcal{D}{\bf
s}~{}P({\bf s})\prod_{i=1}^{\mathcal{N}}s_{n(i)}$ (14)

Solving for these $\mathcal{N}$-point correlation functions can be done more
simply by considering a related quantity:

$P({\bf s},{\bf u})=\exp\Bigg{(}-{1\over 2}{\bf s}^{T}M{\bf s}+{\bf s}^{T}{\bf
u}\Bigg{)}$ (15)

which enables us to write $s_{i}P({\bf s},{\bf u})={\partial/\partial
u_{i}}P[{\bf s},{\bf u}]$. Completing the square, we find:

$\displaystyle P({\bf s},{\bf u})$ $\displaystyle=\exp\Bigg{(}-{1\over 2}({\bf
s}-M^{-1}{\bf u})^{T}M({\bf s}-M^{-1}{\bf u})\Bigg{)}$ (16)
$\displaystyle\cdot\exp\Bigg{(}{1\over 2}{\bf u}^{T}M^{-1}{\bf u}\Bigg{)}$

and then we find that the normalisation factor of this quantity is:

$\mathcal{Z}({\bf u})=\mathcal{Z}(0)\exp({\textstyle{1\over 2}}{\bf
u}^{T}M^{-1}{\bf u})$ (17)

Utilising equation (15) we observe that equation (14) can be rewritten as,

$\left\langle\prod_{i}s_{i}\right\rangle={1\over\mathcal{Z}[0]}\Bigg{(}\prod_{i}{\partial\over\partial
u_{n(i)}}\mathcal{Z}[{\bf u}]\Bigg{)}\Bigg{|}_{{\bf u}=0}$ (18)

This reveals the structure of these correlation functions as products of the
$2$-point correlation functions when $\mathcal{N}$ is even.

$\left\langle\prod_{i}s_{n(i)}\right\rangle=\prod_{i,j,i\neq
j}M^{-1}_{n(i)n(j)}.$ (19)

This result also proves that correlation functions for odd $\mathcal{N}$ go to
zero. Funnily enough, it also masquerades as the traditional equipartition
theorem when we set the state vectors to be positions and velocities in 1D,
${\bf s}=(x,v)$ and choose $M$ such that ${\textstyle{1\over 2}}{\bf
s}^{T}M{\bf s}$ is the product of the Hamiltonian of a harmonic oscillator and
its Lagrange multiplier, $\beta H=\beta({\textstyle{1\over
2}}mv^{2}+{\textstyle{1\over 2}}kx^{2})$.

These choices produce $\langle x^{2}\rangle={1/k\beta}$ and $\langle
v^{2}\rangle={1/m\beta}$, precisely as one should expect!

Later in this paper, the matrix multiplications will be replaced with phase-
space integrals, and the partial derivatives with functional derivatives, but
the essence of this calculation remains unchanged.

##  3 Probability of probabilities

An $N$-body system can be summarised as a set of phase-space coordinates,
${\bf w}_{i}=({\bf x}_{i},{\bf v}_{i}),i=1,2,....,N$. The equations of motion
for such a system are trivial to conceive of–they are $6N$ coupled
differential equations, which integrate the particles forward in the
(gravitational) potential of their peers. The complete Liouville equation is
precisely this,

${{\rm d}f^{(N)}\over{\rm d}t}={\partial f^{(N)}\over\partial
t}+\Big{[}f^{(N)},H^{(N)}\Big{]}=0$ (20)

describing the evolution of a distribution of $N$ particles, $f^{(N)}$ under
the dynamics prescribed by a corresponding $N$-particle Hamiltonian $H^{(N)}$.
Such a set of equations, however, are devoid of insight, and it is obviously
too soon to leave everything to the computers.

We want to predict important features in our $N$-body system; collective
effects which are insensitive to exactly the choice of particles comprising a
spiral arm, or a Langmuir wave. Such an assertion is equivalent to saying that
we don’t care about correlations between particles - or that the particles are
uncorrelated. Then, we have a no-correlations ansatz:

$f^{(N)}(\\{{\bf w}_{i}\\},t)=\prod_{i=1}^{N}f({\bf w}_{i},t).$ (21)

which states that the N-particle distribution function $f^{(N)}$ is the
$N$-fold product of the one-particle distribution function $f$ at the
uncorrelated phase-space coordinates $\\{{\bf w}_{i}\\}$ and gives us, for
large $N$,

${{\rm d}f\over{\rm d}t}={\partial f\over\partial t}+[f,H]=0$ (22)

which is the Collisionless Boltzmann Equation (CBE).

The central assumption we must embrace when we use the CBE is therefore that
we must randomly sample our distribution $f$ to obtain coordinates $\\{{\bf
w}_{i}\\}$.

We know that we can obtain almost any set of coordinates $\\{{\bf w}_{i}\\}$
from the distribution function $f$ insofar as $f$ evaluated at each coordinate
is non-zero,

$f\rightarrow\\{{\bf w}_{i}\\}\textrm{ if }0\notin\\{f({\bf w}_{i})\\}$ (23)

but this also means that almost any distribution function $f$ is consistent
with sampling the particles $\\{{\bf w}_{i}\\}$,

$\\{{\bf w}_{i}\\}\rightarrow f\textrm{ if }0\notin\\{f({\bf w}_{i})\\}$ (24)

and we know the distribution function $f$ is only an artifice of the mind,
while the particles $\\{{\bf w}_{i}\\}$ are real! A theory of fluctuations
should include finiteness noise. Discarding the preconception that we should
fit finite systems with mean field distributions allow us to ask a more
fundamental question: What is the optimal distribution function for a sampling
of particles, $\\{{\bf w}_{i}\\}$?

Now, say the optimal distribution function for a sampling of particles
$\\{{\bf w}_{i}\\}$ exists. Pursuing this program, we might first guess that
the optimal distribution function is the one which maximises the probability
of obtaining $\\{{\bf w}_{i}\\}$, i.e. $f^{(N)}(\\{{\bf w}_{i}\\},t)$, but we
soon see that the valid, discrete probability distribution

$f_{d}({\bf w})={1\over N}\sum_{i=1}^{N}\delta^{6}({\bf w}-{\bf w}_{i})$ (25)

is a profoundly unhelpful result which singularly maximises the product sum.
This result should come at no surprise, as having absolute confidence that the
coordinates $\\{{\bf w}_{i}\\}$ completely defines $f$ relegates us to having
to solve the initial problem of $6N$ coupled force equations.

Now, how confident should we be, really? We know that random sampling a
distribution to get a sample is an irreversible process, so we obviously
shouldn’t be confident at all that any one realisation of $\\{{\bf w}_{i}\\}$
constrains $f$ in the slightest. I suggest that much in the same way we can
sample one choice of $f$ to obtain many $\\{{\bf w}_{i}\\}$, we should
associate with each distribution function, $f$, a joint probability $P_{J}$
that a distribution function $f$ was chosen and sample $\\{{\bf w}_{i}\\}$ was
drawn from it, $P_{J}=P_{J}[f,\\{{\bf w}_{i}\\}]$.

$P_{J}[f,\\{{\bf w}_{i}\\}]=P[f]\prod^{N}_{i=1}f({\bf w}_{i}).$ (26)

We know that once a choice of distribution function $f$ is made, we can only
randomly sample it to obtain particles $\\{{\bf w}_{i}\\}$. This implies that
$P_{J}$ can only depend on $\\{{\bf w}_{i}\\}$ through $\prod_{i}f({\bf
w}_{i})$, or that $P=P[f]$ cannot depend on the sample.

###  3.1 Constraining $P[f]$

Constraining $P[f]$ can be done through a maximum entropy argument.

$P_{J}$ admits a Shannon entropy; which we will denote the joint entropy
$S_{J}$. Integrating over all distribution functions and all samplings of
particles,

$S_{J}=\int\mathcal{D}f\int{\rm d}^{6N}{\bf w}~{}\left(-P_{J}\ln
P_{J}\right).$ (27)

Substituting equation $\eqref{eq:PJnf}$ and simplifying leads us to express
$S_{J}$ in terms of $P[f]$ and $f$,

$S_{J}=\int\mathcal{D}f~{}(-P[f]\ln P[f]+P[f]S_{g}[f])$ (28)

where $S_{g}$ is the Gibbs entropy

$S_{g}=-N\int{\rm d}^{6}{\bf w}f({\bf w})\ln f({\bf w}).$ (29)

Consequently, $S_{g}$ may be interpreted as the number of ways in which one
might randomly sample $N$ particles from the DF $f$.

With this foundation set in stone, we can then ask a perfectly well-defined
and valid question: “Given $f$, how can I choose $P[f]$ so as to maximise my
ignorance with regards to the connection between $f$ and $\\{{\bf w}_{i}\\}$?”

The answer, of course, is to calculate $P[f]$ by extremising the combined
entropy $S_{J}$ with respect to $P[f]$ for some known $f$,

${\delta S_{J}\over\delta P[f]}\Bigg{|}_{f}=0$ (30)

. This grants us the unique probability of $f$ being the ignorance-maximising
distribution function corresponding to a system which evolves via the CBE:

$P[f]\propto\exp(S_{g}[f])$ (31)

which upon proper normalisation such that $\int\mathcal{D}f~{}P[f]=1$ via
dividing by the partition function, $\mathcal{Z}$, we get:

$P[f]={1\over\mathcal{Z}}\exp(S_{g}[f]).$ (32)

The significance of this result is three-fold. Firstly, it confers physical
meaning to the Gibbs entropy $S_{g}$, concretising its dependence on $N$ oft
left neglected. Secondly, it reveals that the states $f$ with the largest
$P[f]$ are also those with the largest $S_{g}[f]$, appealing to our intuition
regarding ignorance: when we know nothing, the best guesses which maximise
$P[f]$ also maximises our ignorance regarding how the particles $\\{{\bf
w}_{i}\\}$ are arranged i.e. maximises $S_{g}$. Finally, it also shows quite
explicitly that choices of $f$ near the state of maximum Gibbs entropy which
we will denote $f_{0}=f|_{\mathrm{max}(S_{g})}$, have $P[f]$ comparable to
$P[f_{0}]$.

This has the unfortunate side-effect of ensuring that the ensemble average of
$f$ (hereon out, the mean-field), is not $f_{0}$ in general,

$\langle f\rangle=\int\mathcal{D}f~{}P[f]f\neq f_{0}$ (33)

since contributions from near-extremum entropy states also influence the
integral, and $P[f]$ is not even about $f=f_{0}$.

Returning to the original problem of finding the optimal $f$ for some
unconstrained sampling of particles $\\{{\bf w}_{i}\\}$, we can substitute for
$P[f]$ from equation (26) to obtain $P_{J}$ which reflects equation (30)’s
entropy-maximisation condition:

$P_{J}[f,\\{{\bf w}_{i}\\}]={1\over\mathcal{Z}}\exp({S_{g}[f]})\prod_{i}f({\bf
w}_{i}).$ (34)

We notice that there is no one optimal choice of $f$, only a distribution of
$f$ which all contribute to the sampling.

It is also evident that $S_{J}$, the Shannon entropy of $P_{J}$ is not wholly
maximised. If it were, then we would instead find that $P_{J}$ is a constant,
as that is the unconstrained prior for any probability subject to entropy
maximisation.

###  3.2 Luck is in the Eye

This then raises the question: What does this partial maximisation of the
entropy mean, and why does it matter?

Let us consider a world where we only have the sample $\\{{\bf w}_{i}\\}$ and
the knowledge that the sample was randomly sampled. We may think we know
nothing about the distribution from which the sample was drawn, but given that
we know we are random sampling an unknown distribution function $f$, we must
at the very least be able to discern between different $f$ on the basis of not
wanting $\\{{\bf w}_{i}\\}$ to be an outlier of our guess of $f$.

Intuitively, this feature of ‘outlierliness’ in random sampling is captured by
comparing the size of the $N$-particle distribution function (equation (21)):
$f^{(n)}=\prod_{i}f({\bf w}_{i})$ with the size of the $N$-particle
distribution function of an average sampling of $f$. If they are of comparable
size, then we would say that $\\{{\bf w}_{i}\\}$ is not an outlier, that we
expected $\\{{\bf w}_{i}\\}$ to be sampled from $f$.

In this subsection, I will use Shannon’s definition of typicality to encode
this intuition. More specifically, I will show that $P_{J}[f,\\{{\bf
w}_{i}\\}]$ as defined in equation (34) encodes an assumption: that choices of
$f$ for which sampling $\\{{\bf w}_{i}\\}$ is perfectly typical, are equally
to likely to have produced the sample. The joint probability $P_{J}[f,\\{{\bf
w}_{i}\\}]$ is well-positioned to describe this, since it describes the
probability that one samples $f$, and then samples $\\{{\bf w}_{i}\\}$ from
$f$.

The $N$-particle distribution function can be expressed in terms of the sample
entropy $S_{s}[f,\\{{\bf w}_{i}\\}]$ of a sample

$S_{s}=-\sum_{i}\ln f({\bf w}_{i}).$ (35)

The sample entropy is a measure of how surprising it is that a distribution
$f$ produces a sample $\\{{\bf w}_{i}\\}$. This can be seen from how $-\ln
f({\bf w}_{i})$ is monotonically increasing as $f({\bf w}_{i})$ diminishes. A
comparatively large sample entropy therefore indicates a very lucky sampling
(or a wrong guess of $f$!). Observe:

$\displaystyle f^{(N)}({\bf w}_{i})$ $\displaystyle=\prod_{i}f({\bf w}_{i})$
(36) $\displaystyle=\exp\Bigg{(}\sum_{i}\ln f({\bf w}_{i})\Bigg{)}$
$\displaystyle=\exp(-S_{s}[f,\\{{\bf w}_{i}\\}]).$

This result is promising as a measure of outlierliness in random sampling;
however it falls short because only a comparatively large sample entropy
indicates an outlierly sampling, and we have no measure of what the average
sample entropy is.

The answer was found by Shannon (1948), who applied Monte Carlo integration to
the sampling entropy. Assuming that $\\{{\bf w}_{i}\\}$ are sampled from $f$,

$\lim_{N\rightarrow\infty}{1\over N}S_{s}[f,\\{{\bf w}_{i}\\}]=-\int{\rm
d}^{6}{\bf w}~{}f\ln f={1\over N}S_{g}.$ (37)

where $\\{{\bf w}_{i}\\}$ are $N$-length sequences. He then went on to intuit
that for finite length samples there must exist some ‘normal-looking’-typical
samples $\\{{\bf w}_{i}\\}_{T}$ that obeyed a similar, but weaker condition
owing to the law of large numbers. He showed that for any $\epsilon>0$ it was
possible to find a sample size $N$ such that the entropy of the typical sample
is close to the Gibbs entropy in the sense that:

${1\over N}S_{g}[f]-\epsilon<{1\over N}S_{s}[f,\\{{\bf w}_{i}\\}_{T}]<{1\over
N}S_{g}[f]+\epsilon.$ (38)

Shannon showed that these typical samples $\\{{\bf w}_{i}\\}_{T}$ bore other
desirable qualities; for sufficiently large $N$ (or sufficiently small
$\epsilon$) the set of typical samples occupies $1-\epsilon$, nearly all of
probability-space, while comprising only a vanishing fraction of sample-space.
In simpler terms, this means that almost every random sampling of $f$ of a
sufficiently large sample size is typical of $f$.

It seems reasonable to say that the perfectly typical sample $\\{{\bf
w}_{i}\\}_{\rm{PT}}$, which defines the perfectly typical sample entropy,
satisfies:

$S_{s}[f,\\{{\bf w}_{i}\\}_{\rm{PT}}]=S_{g}[f].$ (39)

This result defines a sample which is absolutely non-outlierly no matter the
choice of $\epsilon$, is a sample which we should expect and a perfectly
typical sample. We find that we cannot distinguish between different $f$ for
which the sample $\\{{\bf w}_{i}\\}=\\{{\bf w}_{i}\\}_{\rm{PT}}$ is perfectly
typical. Thus we must assign to each $P_{J}[f,\\{{\bf
w}_{i}\\}_{\rm{PT}}]={1/\mathcal{Z}}$ an equal probability of contributing to
the sample; reflecting our ignorance. Then we have

$\displaystyle P_{J}[f,\\{{\bf w}_{i}\\}_{PT}]$
$\displaystyle=P[f]\exp(-S_{s}[f,\\{{\bf w}_{i}\\}_{PT}])$ (40)
$\displaystyle=P[f]\exp(-S_{g}[f])$ $\displaystyle=1/\mathcal{Z}$

that defines $P[f]=\exp(S_{g}[f])/\mathcal{Z}$, e.g. we reclaim equations (32)
and (34).

###  3.3 No Best Coarse-Grained Distribution

This notion of typicality in random sampling can be used to illustrate many
results previously only reachable by intuition. Amongst these results is
perhaps one of the most annoying ones: that there’s no way to smooth a
sampling of particles $\\{{\bf w}_{i}\\}$ to obtain a coarse-grained
distribution (think a histogram!) function $f_{\rm{cg}}$ without introducing
some kind of constraint/belief of what the grain size should be.

If the volume of each coarse-grained phase-space element is $\Delta$, and the
number of particles found within the volume indexed by $\ell$ is $n_{\ell}$,
then the coarse-grained DF $f_{\rm{cg}}$ is

$f_{\rm{cg}\ell}={n_{\ell}\over N\Delta}.$ (41)

By definition, all coarse-grained DFs are perfectly typical of their samples
(and vice versa):

$\displaystyle\prod_{i}f_{\rm{cg}}({\bf w}_{i})$
$\displaystyle=\prod_{\ell}f_{\rm{cg}\ell}^{n_{\ell}}$ (42)
$\displaystyle=\exp\Bigg{(}\sum_{\ell}n_{\ell}\ln f_{\rm{cg}\ell}\Bigg{)}$
$\displaystyle=\exp\Bigg{(}N\sum_{\ell}\Delta f_{\rm{cg}\ell}\ln
f_{\rm{cg}\ell}\Bigg{)}$ $\displaystyle=\exp(-S_{g}[f_{\rm{cg}}]).$

This means $P_{J}[f_{\rm{cg}},\\{{\bf w}_{i}\\}]={1/\mathcal{Z}}$, coarse-
grained distribution functions $f_{\rm{cg}}$ of any $\Delta$ have the same
joint probability with the sample $\\{{\bf w}_{i}\\}$.

The space of coarse-grained DFs includes both the discrete distribution
function (equation (25)) in the limit $\Delta\rightarrow 0$ and the constant
DF in the opposite limit $\Delta\rightarrow\infty$. Observe that each member
of this space is assigned equal probabilities of being the correct DF; and are
equally likely to be the source of the sample.

Let us think of the $\Delta\rightarrow\infty$ limit, as a DF which allows
anything and everything to happen. Sampling this DF can thus be modelled by a
monkey, banging away on a typewriter. On the converse side of this problem,
let us think of the $\Delta\rightarrow 0$ coarse-grained DF as a divinely
exact mandate. Some deity has defined a distribution function for which
particles can only be placed at specific locations; letting us choose
(randomly) where we want to place them from amongst $N$ locations.

Inspecting $\\{{\bf w}_{i}\\}$, how can we possibly know if it was generated
by a monkey banging on a typewriter, or by exercising what little freedom we
have under a divine mandate? We cannot: unless we have the ability to gain
additional samples of $f$, so as to gain a better understanding of $f$, we
cannot investigate the veracity of either possibility.

And so we can only assume that all $f_{\rm{cg}}$, everything between monkeys
on typewriters and divine law are equally probable in having generated the
sample $\\{{\bf w}_{i}\\}$, which is the only observable. This result is the
cornerstone of this field theory, and establishes the primacy of the Shannon
entropy of the DF $f$ (i.e. the Gibbs entropy) as the correct entropy to use
in this random-sampling problem.

To put it more simply, this theory begins by acknowledging that samples
$\\{{\bf w}_{i}\\}$ cannot be binned into histograms of the form $f_{\rm{cg}}$
without choosing a bin-size. This is blindingly obvious.

###  3.4 Applying $P_{J}$?

$P_{J}$ is the probability that first, the DF $f$ is chosen, and then a series
of samples are sampled. Given that we don’t know which $f$ is chosen, but we
do know $\\{{\bf w}_{i}\\}$ it is only sensible to define the sample
probability $\mathcal{P}$ by integrating $P_{J}$ over $f$:

$\displaystyle\mathcal{P}[\\{{\bf w}_{i}\\}]$
$\displaystyle=\int\mathcal{D}f~{}P_{J}[f,\\{{\bf w}_{i}\\}]$ (43)
$\displaystyle=\int\mathcal{D}f~{}P[f]\prod_{i}f({\bf w}_{i})$
$\displaystyle=\left\langle\prod_{i}f({\bf w}_{i})\right\rangle.$

$\mathcal{P}[\\{{\bf w}_{i}\\}]$ is a DF independent measure of the
probability of sampling the sample $\\{{\bf w}_{i}\\}$ that pushes the onus of
defining the sample probability onto the constraints imposed on the space of
DFs. Invoking a binomial expansion by substituting $f=f_{0}+\delta f$, we see
we can express the sample probability as an expansion in $\mathcal{N}$-point
correlation functions $C_{\mathcal{N}}$ (that measure the correlations in
$\delta f$ at $\mathcal{N}$ points in phase-space):

$\displaystyle\mathcal{P}[\\{{\bf w}_{i}\\}]=\prod_{i}f_{0}({\bf
w}_{i})\Bigg{(}1+\sum_{j}\left\langle{\delta f({\bf w}_{j})\over f_{0}({\bf
w}_{j})}\right\rangle$ (44) $\displaystyle+{1\over 2!}\sum_{j,k,j\neq
k}\left\langle{\delta f({\bf w}_{j})\delta f({\bf w}_{k})\over f_{0}({\bf
w}_{j})f_{0}({\bf w}_{k})}\right\rangle+...\Bigg{)}$
$\displaystyle=\prod_{i}f_{0}({\bf
w}_{i})\Bigg{(}\sum_{\mathcal{N}=0}^{N}{1\over\mathcal{N}!}\sum_{\\{{\bf
w}_{i}\\}_{P}\in P(\mathcal{N})}{C_{\mathcal{N}}[\\{{\bf
w}_{i}\\}_{P}]\over\prod f_{0}[\\{{\bf w}_{i}\\}_{P}]}\Bigg{)}$

where $P(\mathcal{N})$ represents the set of all permutations of selecting
$\mathcal{N}$ particles (via their phase-space coordinates) from $N$
particles, non-repeating.

This result is great fun. In the construction of $P_{J}$, we have assumed that
each sampling in $\\{{\bf w}_{i}\\}$ is sampled independently from each other
from $f$, $\mathcal{P}$ need not reflect random sampling: we manifestly see
$\mathcal{N}$-point correlations appearing in the calculation.

Since each sampling of $f$ is independent of the previous one, we must realise
that these departures from ‘random sampling’ arise because $f$ itself has
natural internal correlations, which when sampled are imparted onto $\\{{\bf
w}_{i}\\}$.

Later in this paper, we will see that these $\mathcal{N}$-point correlations
are sourced from physical interactions. In the next section, I describe how we
can insert constraints (physical or not) into this statistical theory to this
effect.

##  4 Entropy Maximisation with constraints

Maximising $S_{J}$ with respect to multiple constraints

$\langle G_{n}[f]\rangle=G_{n0}$ (45)

on the space of distribution $f$ can be done via the method of Lagrange
multipliers. Let us denote $P_{\\{G_{n}\\}}[f]$ as the $P[f]$ which is
obtained by maximising

$S_{J}-\int\mathcal{D}f~{}P_{\\{G_{n}\\}}[f]\sum_{n}\beta_{n}G_{n}[f]=0$ (46)

where $\beta_{n}$ are the Lagrange multipliers, which are defined to maintain
the constraints. Then we find:

$P_{\\{G_{n}\\}}[f]={1\over\mathcal{Z}}\exp\left(S_{g}[f]-\sum_{n}\beta_{n}G_{n}[f]\right).$
(47)

In studies of collisionless systems, we will desire constraints to represent
quantities which are globally conserved by the time-evolution of the
collisionless Boltzmann equation, i.e. total energy $E[f]$, total angular
momentum ${\bf L}[f]$, etc. $S_{g}[f]$ is already conserved as one of the
Casimir invariants of the CBE, which are integrals of functions of the DF
$g(f)$ which go to zero at infinity rapidly enough to satisfy the last
equality.

${{\rm d}\over{\rm d}t}G[f]=\int{\rm d}^{6}{\bf w}~{}{\partial\over\partial
t}g(f)=\int{\rm d}^{6}{\bf w}~{}[g(f),-H]=0$ (48)

Notice that when $\\{G_{n}[f]\\}$ are conserved by the CBE, that ${\rm
d}G_{n}[f]/{\rm d}t=0$ then ${\rm d}P_{\\{G_{n}\\}}[f]/{\rm d}t=0$, which
means this method of assigning a probability to each $f$ is strictly time-
independent (even in the nonlinear regime!). This ensures that the ensemble
averages obey the CBE if their arguments also obey the CBE, i.e. ${{\rm
d}\langle f\rangle/{\rm d}t}=0$ since ${\rm d}f/{\rm d}t=0$.

###  4.1 Nonlinear Equilibria

It is obvious that the above statement is not terribly impressive, only
delineating that $\langle f\rangle$ is a solution to Boltzmann’s equation. A
stronger version of this result can be proven with a little more work. Here I
will prove that when $\\{G_{n}[f]\\}$ are conserved by the CBE, their mean
fields $\langle f\rangle$ are in equilibria; that is that they have no
explicit time dependence.

Consider first the system with a mean energy imposed $\langle E\rangle=E_{0}$
so that $P[f]=P_{E}[f]$ takes the form

$P_{E}[f]={1\over\mathcal{Z}}\exp(S_{g}[f]-\beta E[f]).$ (49)

$E=E[f]$ is the energy of a system; which is always conserved if the system is
isolated and lacks dissipation

$E[f]=N\int{\rm d}^{6}{\bf w}~{}f(H-{{\textstyle{1\over 2}}}m\Phi[f])$ (50)

and

$H({\bf w})={1\over N}{\delta E[f]\over\delta f({\bf w})}=m({\textstyle{1\over
2}}{\bf v}^{2}+\Phi({\bf x})[f]+\Phi_{e}({\bf x}))$ (51)

is the Hamiltonian of the system. I choose a Hamiltonian with self-consistent
gravitational potential $\Phi[f]$,

$\Phi[f]=-\int{\rm d}^{6}{\bf w}^{\prime}~{}{\mathrm{G}M\over|{\bf x}-{\bf
x}^{\prime}|}f({\bf w}^{\prime})$ (52)

and external potential $\Phi_{e}$. Note then,

$\displaystyle{\partial\langle f\rangle_{E}\over\partial t}$
$\displaystyle=\int\mathcal{D}f~{}P_{E}[f]{\partial f\over\partial t}$ (53)
$\displaystyle=\int\mathcal{D}f~{}P_{E}[f]([f,-H])$
$\displaystyle={1\over\beta}\int\mathcal{D}f~{}P_{E}[f]([f,-\beta H])$
$\displaystyle={1\over\beta}\int\mathcal{D}f~{}P_{E}[f]\Bigg{(}\left[f,{\delta
S_{g}\over\delta f}-\beta{\delta E\over\delta f}\right]\Bigg{)}$
$\displaystyle={1\over\beta}\int\mathcal{D}f~{}\left[f,{\delta
P_{E}[f]\over\delta f}\right]$ $\displaystyle=0.$

Between the third and fourth equalities, I have used the fact that
$[f,g(f)]=0$ for any differentiable function $g$ of $f$, and between the
fourth and the fifth equalities, I have used the chain rule. The final
equality is found via integration by parts, and the fact that $\delta f({\bf
w})/\delta f({\bf w}^{\prime})=\delta^{6}({\bf w}-{\bf w}^{\prime})$ is
symmetric, but the Poisson brackets are asymmetrical. The boundary terms of
the functional integral go to zero because they are constants, and are
eliminated by the Poisson brackets.

This result can easily be extended to show that an ensemble defined by a mean
energy and a mean angular momentum is in equilibrium, but only in the
corotating frame, where the Hamiltonian in the rotating frame is
$\tilde{H}=H+\Omega_{p}J_{z}$ owing to the canonical transformation to said
rotating frame.

Here, $P[f]=P_{E,L_{z}}[f]$ takes the form,

$P_{E,L_{z}}[f]={1\over\mathcal{Z}}\exp(S_{g}[f]-\beta
E[f]-\beta_{L_{z}}L_{z}[f])$ (54)

If the angular momentum of a system is $L_{z}$,

$L_{z}[f]=N\int{\rm d}^{6}{\bf w}~{}f\left({\bf x}-\int{\rm d}^{6}{\bf
w}^{\prime}~{}f({\bf w}^{\prime}){\bf x}^{\prime}\right)\times{\bf v}.$ (55)

then the angular momentum of a particle in the system is:

$\displaystyle J_{z}({\bf w})$ $\displaystyle={1\over N}{\delta
L_{z}[f]\over\delta f({\bf w})}$ (56) $\displaystyle=\left({\bf x}-\int{\rm
d}^{6}{\bf w}^{\prime}~{}f({\bf w}^{\prime}){\bf x}^{\prime}\right)\times{\bf
v}$ $\displaystyle-{\bf x}\times\int{\rm d}^{6}{\bf w}^{\prime}f({\bf
w}^{\prime}){\bf v}^{\prime}.$

where the first term of $J_{z}$ is the angular momentum we know and love, and
the second term is the correction term to the angular momentum arising from
the centre of mass drift.

I identify the angular speed of the corotating frame as
$\Omega_{p}={\beta_{L_{z}}/\beta}$, and find that in the corotating frame:

$\displaystyle{\partial\langle f\rangle_{E,L_{z}}\over\partial t}$
$\displaystyle=\int\mathcal{D}f~{}P_{E,L_{z}}{\partial f\over\partial t}$ (57)
$\displaystyle=\int\mathcal{D}f~{}P_{E,L_{z}}[f]([f,-\tilde{H}])$
$\displaystyle={1\over\beta}\int\mathcal{D}f~{}P_{E,L_{z}}[f]\left(\left[f,-\beta
H-\beta_{L_{z}}J_{z}\right]\right)$
$\displaystyle={1\over\beta}\int\mathcal{D}f~{}P_{E,L_{z}}[f]$
$\displaystyle\cdot\left(\left[f,{\delta S_{g}\over\delta f}-\beta{\delta
E\over\delta f}-\beta_{L_{z}}{\delta L_{z}\over\delta f}\right]\right)$
$\displaystyle={1\over\beta}\int\mathcal{D}f~{}\left[f,{\delta
P_{E,L_{z}}[f]\over\delta f}\right]$ $\displaystyle=0.$

It is an elementary exercise to prove that ensembles with a mean drift in the
position of the barycenter as well as its velocity, ${\bf x},{\bf v}$, in
addition to means imposed for $E,L_{z}$, are in equilibrium in the corotating,
co-drifting frame. A similar extension to Casimir invariants can be inserted
between the third and fourth equalities.

Thus we have a non-perturbative result which defines a family of nonlinear
equilibria. These equilibria need not be stable to all perturbations; only
those induced by the structure of $P[f]$, which are natural to the system.
Actually calculating the distributions of these nonlinear equilibria is far
more involved, and can be done perturbatively. We investigate how such
calculations can be made in the next section.

##  5 Field theory of distributions

The crowning jewel of Section 3 is the representation of correlation functions
$C_{\mathcal{N}}$ as contributors to the sampling probability
$\mathcal{P}[\\{{\bf w}_{i}\\}]$ (equation (44)). This result describes how
correlations raise or lower the probability of observing a constrained
$N$-particle system in a certain configuration.

I will now explain (with the simplest possible toy model) how we can proceed
to calculate these correlation functions, and in doing so open CBE
calculations up to the powerful field theory formalism.

###  5.1 Field Theory Basics

A perturbation field theory is very much like any other perturbation theory:
it begins by finding a nice, physically helpful point to Taylor expand around.
In this theory, that point is the state of extremised Gibbs entropy, $f_{0}$,
and this expansion is well-known as the Saddle Point Approximation (of
Lagrangian mechanics fame).

I propose we study a toy model in which I choose where I expand the
perturbation theory around. This equivalent to choosing the distribution
function of maximum/extremum Gibbs entropy $S_{g}$ (and in doing so
constraining the underlying ensemble of sampled particles). Do note that this
choice of constraint is incredibly artificial!

$P_{f_{0}}[f]={1\over\mathcal{Z}}~{}\exp\Bigg{(}\int{\rm d}^{6}{\bf
w}~{}-Nf\ln f-\beta_{f_{0}}({\bf w})f\Bigg{)}$ (58)

We select $\beta_{f_{0}}$ such that,

$\beta_{f_{0}}({\bf w})={\delta S_{g}\over\delta f}\Bigg{|}_{f=f_{0}}$ (59)

to eliminate the first order in the Taylor expansion of equation (58), where
we expand $f=f_{0}+\delta f$. What remains then after we have fixed our choice
of $f_{0}$, and neglected constants is:

$P_{f_{0}}[\delta f]={1\over\mathcal{Z}}\exp\Bigg{(}\int{\rm d}^{6}{\bf
w}~{}-{1\over 2!}{N\delta f^{2}\over f_{0}}+{1\over 3!}{N\delta f^{3}\over
f_{0}^{2}}+...\Bigg{)}.$ (60)

Notice that we have teased a Gaussian form out of $P_{f_{0}}[f]$. The
integration that has to occur over all $\delta f$ however is not between
$(-\infty,\infty)$ but rather between $(-f_{0},\infty)$, which means that we
are not quite able to solve for $\langle A\rangle$ just yet.

Intuition into this issue can be obtained by observing that $P[\delta
f=-f_{0}]={1/\mathcal{Z}}\exp(-N)$, which is vanishingly small compared to the
maximum of $P_{f_{0}}[\delta f]={1/\mathcal{Z}}$ with even $N\approx 50$ or
so, indicating that the correction to probability-space where the Gaussian
integrand is significant is negligible. $P[f]$ does become large and highly
oscillatory at large, negative $\delta f$, however this only increases our
motivation in using a truncated series in $\delta f$ to represent
$P_{f_{0}}[\delta f]$, since we know that including more terms in the
expansion brings us further away from the underlying model where the negative
$f$ domain is excluded.

So now the ensemble average of some system-wide quantity represented by the
functional $A[f]$ is:

$\langle A\rangle=\int^{\infty}_{\infty}\mathcal{D}\delta
f~{}{1\over\mathcal{Z}}\exp\Bigg{(}\int{\rm d}^{6}{\bf w}~{}-{1\over
2!}{N\delta f^{2}\over f_{0}}+...\Bigg{)}A[\delta f]$ (61)

Now consider a theory in which only the first non-trivial term in the argument
of the exponent exists; i.e. the quadratic term which then turns the integral
into a Gaussian integral, as presented in equation (61). This is known as a
free theory (Peskin & Schroeder, 1995), simply because there are no
’interactions’, which are the subsequent cubic/quartic/… terms.

Denoting correlator brackets taken with respect to the free theory with the
hyperscript $0$, we note that the first non-trivial correlator evaluates to
zero, $\langle\delta f\rangle^{0}_{f_{0}}=0$, due to the even nature of the
Gaussian. The second non-trivial correlator may be computed via usage of the
equipartition theorem, and takes the form:

$C^{0}_{2,f_{0}}({\bf w},{\bf w}^{\prime})=\langle\delta f({\bf w})\delta
f({\bf w}^{\prime})\rangle^{0}_{f_{0}}={f_{0}\over N}\delta^{6}({\bf w}-{\bf
w}^{\prime}).$ (62)

Higher order correlation functions for the free theory are only products of
two-point correlation functions, which are removed via division by the
partition function $\mathcal{Z}$. (e.g. Feynman’s result that only connected
diagrams contribute to the correlations.) Thus the theory is ‘free’.

Observe that the two-point correlation function in equation (62) is inversely
proportional to $N$, that is that couplings between different points in phase-
space are controlled by the size of $N$—or that $N$ plays the role of a
coupling parameter.

The effect of a coupling parameter is more easily examined if the free theory
(which usually is defined without any inherent couplings, and is thus ’free’),
has no explicit dependence on it. Thus we now change the normalisation of $f$,
defining:

$\delta\mu=\sqrt{N}\delta f$ (63)

and so the newly normalised free theory two-point correlator is

$\tilde{C}_{f_{0}}({\bf w},{\bf w}^{\prime})=\langle\delta\mu({\bf
w})\delta\mu({\bf w}^{\prime})\rangle=f_{0}\delta^{6}({\bf w}-{\bf
w}^{\prime}),$ (64)

and the equivalent $P_{f_{0}}[\delta\mu]$ is now:

$P_{f_{0}}[\delta\mu]={1\over\mathcal{Z}}\exp\Bigg{(}\int{\rm d}^{6}{\bf
w}-{1\over 2}{\delta\mu^{2}\over f_{0}}+{1\over\sqrt{N}}{1\over
3!}{\delta\mu^{3}\over f_{0}^{2}}+...\Bigg{)}$ (65)

where the new normalisation of the correlator has made explicit the dependence
of the interaction terms on the coupling parameters, which are negative
integer powers of $\sqrt{N}$.

The existence of couplings beyond the free theory allows us to calculate
(small) corrections to it. If the coupling parameter $1/\sqrt{N}$ is
sufficiently small, we can apply the expansion,

$\exp\Bigg{(}\int{\rm d}^{6}{\bf w}{1\over\sqrt{N}}{1\over
3!}{\delta\mu^{3}\over f_{0}^{2}}\Bigg{)}\approx 1+\int{\rm d}^{6}{\bf
w}{1\over\sqrt{N}}{1\over 3!}{\delta\mu^{3}\over f_{0}^{2}}+...$ (66)

which produces a series of interaction terms which diminish with increasing
order as the expansion parameter $1/\sqrt{N}$.

###  5.2 Feynman Diagrams

Feynman (1949) describes a way to calculate correlators which respect Wick
(1950)’s theorem via a diagrammatic approach to calculating combinatorial
factors. We use his elegant theory like a club, so as to compute the simplest
correlation function imaginable.

I want to compute the first non-trivial correction to the $1$-point
correlation function $\langle\delta f({\bf w})\rangle$ because it is also the
first non-trivial correction to the sample probability $\mathcal{P}[\\{{\bf
w}_{i}\\}]$ (equation (43)). This can be done by starting with one factor of
$\delta f$, and then looking for the first interaction term which bumps the
number of factors of $\delta f$ up to an even number, so as to satisfy the
generalised equipartition theorem.

In this theory, that is the first order in the cubic expansion we have just
seen. Thus, we may write:

$\displaystyle\langle\delta f({\bf w})\rangle\approx$
$\displaystyle\int\mathcal{D}\delta\mu~{}\delta\mu({\bf w}){1\over
3!N}\Bigg{(}\int{\rm d}^{6}{\bf w}_{a}{\delta\mu_{a}^{3}\over
f_{0a}^{2}}\Bigg{)}$ (67) $\displaystyle\cdot\exp\Bigg{(}-\int{\rm d}^{6}{\bf
w}{1\over 2}{\delta\mu^{2}\over f_{0}}\Bigg{)}$

This calculation is simplified by the generalised equipartition theorem, which
reduces it to calculating:

$\displaystyle\langle\delta f({\bf w})\rangle\approx\int{\rm d}^{6}{\bf
w}_{a}{1\over
3!Nf_{0a}^{2}}\prod_{\forall\textrm{Perm.}}\tilde{C}_{f_{0}}({\bf w},{\bf
w}_{a})\tilde{C}_{f_{0}}({\bf w}_{a},{\bf w}_{a})$ (68)

where the product goes over all permutations of the two-point correlators.
While this combinatorial calculation is simple enough to do by inspection
(there are three ways to connect ${\bf w}$ to three different ${\bf w}_{a}$s),
it generally becomes difficult really quickly.

Feynman describes how to use his diagrams to solve for these combinatorial
constants:

$\displaystyle\langle\delta f({\bf w})\rangle$
$\displaystyle=\Bigg{(}{1\over\sqrt{N}}\delta f({\bf w})~{}\vbox{\hbox{
\leavevmode\hbox to13.9pt{\vbox to4.24pt{\pgfpicture\makeatletter\hbox{\hskip
27.72864pt\lower-2.12134pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-14.22638pt}{0.0pt}\pgfsys@lineto{-25.6073pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@moveto{-23.48596pt}{0.0pt}\pgfsys@curveto{-23.48596pt}{1.17159pt}{-24.43571pt}{2.12134pt}{-25.6073pt}{2.12134pt}\pgfsys@curveto{-26.77888pt}{2.12134pt}{-27.72864pt}{1.17159pt}{-27.72864pt}{0.0pt}\pgfsys@curveto{-27.72864pt}{-1.17159pt}{-26.77888pt}{-2.12134pt}{-25.6073pt}{-2.12134pt}\pgfsys@curveto{-24.43571pt}{-2.12134pt}{-23.48596pt}{-1.17159pt}{-23.48596pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{-25.6073pt}{0.0pt}\pgfsys@fill\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-25.6073pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\times{1\over
3!\sqrt{N}}~{}\vbox{\hbox{ \leavevmode\hbox to21.31pt{\vbox
to23.96pt{\pgfpicture\makeatletter\hbox{\hskip
13.50226pt\lower-11.97751pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{}{{}}{}
{}{}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{-11.38092pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@moveto{-9.25958pt}{0.0pt}\pgfsys@curveto{-9.25958pt}{1.17159pt}{-10.20934pt}{2.12134pt}{-11.38092pt}{2.12134pt}\pgfsys@curveto{-12.5525pt}{2.12134pt}{-13.50226pt}{1.17159pt}{-13.50226pt}{0.0pt}\pgfsys@curveto{-13.50226pt}{-1.17159pt}{-12.5525pt}{-2.12134pt}{-11.38092pt}{-2.12134pt}\pgfsys@curveto{-10.20934pt}{-2.12134pt}{-9.25958pt}{-1.17159pt}{-9.25958pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{-11.38092pt}{0.0pt}\pgfsys@fill\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-11.38092pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{5.69046pt}{9.85617pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@moveto{7.8118pt}{9.85617pt}\pgfsys@curveto{7.8118pt}{11.02776pt}{6.86205pt}{11.97751pt}{5.69046pt}{11.97751pt}\pgfsys@curveto{4.51888pt}{11.97751pt}{3.56912pt}{11.02776pt}{3.56912pt}{9.85617pt}\pgfsys@curveto{3.56912pt}{8.68459pt}{4.51888pt}{7.73483pt}{5.69046pt}{7.73483pt}\pgfsys@curveto{6.86205pt}{7.73483pt}{7.8118pt}{8.68459pt}{7.8118pt}{9.85617pt}\pgfsys@closepath\pgfsys@moveto{5.69046pt}{9.85617pt}\pgfsys@fill\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{5.69046pt}{9.85617pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}{{}}{}
{}{}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@lineto{5.69046pt}{-9.85617pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@moveto{7.8118pt}{-9.85617pt}\pgfsys@curveto{7.8118pt}{-8.68459pt}{6.86205pt}{-7.73483pt}{5.69046pt}{-7.73483pt}\pgfsys@curveto{4.51888pt}{-7.73483pt}{3.56912pt}{-8.68459pt}{3.56912pt}{-9.85617pt}\pgfsys@curveto{3.56912pt}{-11.02776pt}{4.51888pt}{-11.97751pt}{5.69046pt}{-11.97751pt}\pgfsys@curveto{6.86205pt}{-11.97751pt}{7.8118pt}{-11.02776pt}{7.8118pt}{-9.85617pt}\pgfsys@closepath\pgfsys@moveto{5.69046pt}{-9.85617pt}\pgfsys@fill\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{5.69046pt}{-9.85617pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\Bigg{)}$ (69)
$\displaystyle={3\over 3!N}\Bigg{(}\delta f({\bf w})~{}\vbox{\hbox{
\leavevmode\hbox to42.27pt{\vbox
to29.09pt{\pgfpicture\makeatletter\hbox{\hskip
27.72864pt\lower-14.54355pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{14.14355pt}{0.0pt}\pgfsys@curveto{14.14355pt}{7.81136pt}{7.81136pt}{14.14355pt}{0.0pt}{14.14355pt}\pgfsys@curveto{-7.81136pt}{14.14355pt}{-14.14355pt}{7.81136pt}{-14.14355pt}{0.0pt}\pgfsys@curveto{-14.14355pt}{-7.81136pt}{-7.81136pt}{-14.14355pt}{0.0pt}{-14.14355pt}\pgfsys@curveto{7.81136pt}{-14.14355pt}{14.14355pt}{-7.81136pt}{14.14355pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{0.0pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {}{{}}{}
{}{}{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}{}\pgfsys@moveto{-14.22638pt}{0.0pt}\pgfsys@lineto{-25.6073pt}{0.0pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{
}{{}{{{}}}{{}}{}{}{}{}{}{}{}{}{}{\pgfsys@moveto{-23.48596pt}{0.0pt}\pgfsys@curveto{-23.48596pt}{1.17159pt}{-24.43571pt}{2.12134pt}{-25.6073pt}{2.12134pt}\pgfsys@curveto{-26.77888pt}{2.12134pt}{-27.72864pt}{1.17159pt}{-27.72864pt}{0.0pt}\pgfsys@curveto{-27.72864pt}{-1.17159pt}{-26.77888pt}{-2.12134pt}{-25.6073pt}{-2.12134pt}\pgfsys@curveto{-24.43571pt}{-2.12134pt}{-23.48596pt}{-1.17159pt}{-23.48596pt}{0.0pt}\pgfsys@closepath\pgfsys@moveto{-25.6073pt}{0.0pt}\pgfsys@fill\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-25.6073pt}{0.0pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}}}\Bigg{)}$
$\displaystyle={3\over 3!N}\delta^{6}(0).$

Where between the first and the second equalities, we have matched up the
external vertex $\delta f({\bf w})$ to all three internal vertices (internal
being within an integral), producing a combinatorial factor of $3$. Then there
are no other ways to pair up the remaining two vertices in each case. Finally,
the loop integral amounts to $\int{\rm d}^{6}{\bf
w}_{a}\delta^{6}(0)\delta^{6}({\bf w}_{a}-{\bf w})=\delta^{6}(0)$, completing
the calculation.

Now there are two results worth taking note of here. Note $C_{1}=\langle\delta
f\rangle$ and $C_{2}=\langle\delta f\delta f^{\prime}\rangle$ are of the same
leading order in the coupling parameter. This is a ‘lucky’ feature: $C_{1}$’s
leading order is $0$, and its second order is $O(1/N)$. $C_{2}$’s first order
is $O(1/N)$, and so to leading order, they are of the same size in the
coupling parameter. Generally however, $C_{n}\sim O(1/N^{n-1})$.

##  6 Coarse Graining; Discretisation

Perhaps more pertinently, there is a huge, glaring $\delta^{6}(0)$ sitting
front and centre in the expression which arises as a feature of our
calculation. This divergence arises because there are uncountably infinite
degrees of freedom in the defining of $f$ while we only have a finite number
of particles to constrain $f$, and can be mastered by discretising the space
of distribution functions.

Discretising by setting $f({\bf w}_{i})=f_{i}$,

$f_{i}={1\over\Delta}\int_{\Delta}{\rm d}^{6}{\bf w}~{}f({\bf w})$ (70)

as the mean of $f$ within a discrete phase-space element (henceforth, pixel)
of volume $\Delta$ and $f({\bf w}_{i})$ is a constant within that element.
This choice of coarse-graining causes $f({\bf w}_{i})$ to only change due to
fluxes of probability at the boundaries of each pixel, thus decoupling
microscopic dynamics happening within a pixel, from the macroscopic
observables taken by studying a set of pixels. The correlation functions then
map accordingly;

$\displaystyle\langle\delta f({\bf w})\delta f({\bf
w}^{\prime})\rangle={f_{0}\over N}\delta^{6}({\bf w}-{\bf
w}^{\prime})\rightarrow\langle\delta f_{i}\delta f_{j}\rangle={f_{0i}\over
N\Delta}\delta_{ij}.$ (71)

From hereon out I have chosen to retain the continuum formalism, since
$\Delta$ is the smallest scale we are interested in.

The regularised first order correction to the expectation of $f({\bf w})$ can
be expressed as,

$\langle\delta f({\bf w})\rangle\approx{\delta^{6}(0)\over
2N}\rightarrow{1\over 2N\Delta}$ (72)

and the general expansion to order $n$ for the corrections to $\langle\delta
f\rangle_{f}$ is,

$\langle f({\bf w})\rangle=f_{0}({\bf w})\sum_{i=0}^{n}c_{i}\Bigg{(}{1\over
N\Delta f_{0}}\Bigg{)}^{i}$ (73)

with combinatorial coefficients $c_{i}$. Discretisation reveals that the
expansion parameter is the inverse of $N_{f}=N\Delta f_{0}$, which is the
number of particles expected to lie within a discrete phase-space element from
an $N$-fold sampling of the mode-field.

If $N_{f}$ is large, then the expansion (for some proper choice of $n$)
approaches asymptotically the real value of $\langle\delta f({\bf
w})\rangle_{f}$; the statistics we might extract from a gridded ’average
N-particle cloud’ becomes more refined the more particles lie within each
grid, becoming ill-defined when the expected number of particles is less than
$1$. This arises because the contribution from a distribution function $f$ and
sampling of particles $\\{{\bf w}_{i}\\}$ to the ensemble average (equation
(32)), $\mathcal{P}[f](\\{{\bf w}_{i}\\})$ only constrains $P[f]$ if $f({\bf
w}_{i})$ is non-zero. This causes a massive degeneracy as grids with no
particles do not contribute to the ensemble average and do not constrain
$P[f]$ whatsoever.

While the Boltzmann equation itself is a continuous probability flow down to
and even beyond the level of individual particles, taking the ensemble average
introduces divergences which must be regularised via coarse-graining. This
necessitates a ’smoothing’-like procedure in which the distribution function
is gridded over distances larger than the inter-particle separation. I
therefore introduce a statistically motivated method of coarse-graining which
is independent of any underlying physics.

Under this model, the distribution functions strictly evolve under the CBE;
however if and when we decide to take an observable, we must coarse-grain to
obtain ensemble averages. This is one solution to the long-standing issue that
no coarse-grained equivalent to the CBE exists, though perhaps not a
satisfying one—I claim there is no need for a coarse-grained CBE!

Further discussion of this coarse-graining regularisation scheme and its
impact on the Boltzmann theory will be reserved for the discussion section.

##  7 Mean Energy Theory

Now of chief interest is the imposition of a mean energy constraint, which I
denote with the $E$ subscript:

$\langle E\rangle_{f_{0},E}=E_{0}$ (74)

Setting a mean energy leads us to write,

$\displaystyle P_{f_{0},E}[f]={1\over\mathcal{Z}}\exp\Bigg{(}$ $\displaystyle
S_{g}-\beta E[f]$ (75) $\displaystyle-N\int{\rm d}^{6}{\bf w}~{}f(-\ln
f_{0}+1+\beta H_{0})\Bigg{)}$

where we have redefined $\beta_{f}({\bf w})$ so as to maintain our choice of
the mean field $f_{0}$, and $\beta>0$ means $P[f]$ favours bound states over
unbound states.

Expansion in the parameter $\delta\mu$ about
$f=f_{0}+{1\over\sqrt{N}}\delta\mu$ from equation (63) gives us;

$\displaystyle P_{f_{0},E}[\delta\mu]={1\over\mathcal{Z}}\exp$
$\displaystyle\Bigg{(}-{1\over 2}\int{\rm d}^{6}{\bf w}{\rm d}^{6}{\bf
w}^{\prime}~{}\delta\mu\delta\mu^{\prime}$ (76)
$\displaystyle\times\Bigg{[}{1\over f_{0}}\delta^{6}({\bf w}-{\bf
w}^{\prime})-{\mathrm{G}Mm\beta\over|{\bf x}-{\bf x}^{\prime}|}\Bigg{]}$
$\displaystyle+\int{\rm d}^{6}{\bf w}{1\over\sqrt{N}}{1\over
3!}{\delta\mu^{3}\over f_{0}^{2}}+...\Bigg{)}\ $

and the corresponding normalised two-point correlation function
$\tilde{C}_{2,f_{0},E}=\langle\delta\mu({\bf w})\delta\mu({\bf
w}^{\prime})\rangle_{f_{0},E}$ is

$\displaystyle\tilde{C}_{f_{0},E}({\bf w},{\bf w}^{\prime})$
$\displaystyle=f_{0}\delta^{6}({\bf w}-{\bf w}^{\prime})$ (77)
$\displaystyle+f_{0}({\bf w})f_{0}({\bf w}^{\prime})X({\bf x},{\bf
x}^{\prime})$

where the spatial correlation function $X({\bf x},{\bf x}^{\prime})$ satisfies
the differential equation,

${1\over 4\pi A}\nabla^{2}X+\int f_{0}({\bf w}){\rm d}^{3}{\bf v}\times
X=-\delta^{3}({\bf x}-{\bf x}^{\prime})$ (78)

and $A=\beta GMm$. One can check via substitution that $\langle\delta\mu({\bf
w})\delta\mu({\bf w}^{\prime})\rangle_{f,E}$ is indeed the inverse of the
expression in the square brackets in equation (76).

Then the first order correction to $\langle f\rangle_{f_{0},E}$ is;

$\displaystyle\langle\delta f({\bf w})\rangle_{f_{0},E}$
$\displaystyle\approx{1\over 2}\int{\rm d}^{6}{\bf w}_{a}~{}{1\over
Nf_{0}^{2}}$ (79) $\displaystyle\times\tilde{C}_{f,E}({\bf w}_{a},{\bf
w}_{a})\tilde{C}_{f,E}({\bf w}_{a},{\bf w})$

Where $\tilde{C}_{f,E}({\bf w},{\bf w})$

$\tilde{C}_{f,E}({\bf w},{\bf w})={f_{0}\over\Delta}$ (80)

is finite due to our coarse-graining in phase-space, and the gravitational
correlation is set to zero. In the idealised picture of the Boltzmann
equation, phase-space is incompressible; what this means is that a high
density phase-space element cannot grow in density by gravitationally
depleting neighbouring elements of particles; it can only attract other
elements of high density together, displacing low density elements in the
process. This inability to enrich oneself is what justifies setting the
spatial self-correlation to zero, even though the spatial correlations are the
strongest for neighbouring position-space elements.

$\displaystyle\langle$ $\displaystyle\delta f({\bf w})\rangle_{f_{0},E}$ (81)
$\displaystyle\approx{1\over 2}\int{\rm d}^{6}{\bf w}_{a}~{}{1\over
Nf_{0}^{2}}\left({f_{0}({\bf w}_{a})\over\Delta}\right)$
$\displaystyle\times\left(f_{0}\delta^{6}({\bf w}-{\bf w}_{a})+f_{0}({\bf
w})f_{0}({\bf w}_{a})X({\bf x},{\bf x}_{a})\right)$ $\displaystyle={1\over
2}f_{0}({\bf w})\Bigg{(}{1\over N_{f}({\bf w})}+\int{\rm d}^{6}{\bf
w}_{a}~{}{f_{0}({\bf w}_{a})\over N_{f}({\bf w}_{a})}X({\bf x}_{a},{\bf
x})\Bigg{)}$ $\displaystyle={1\over 2N\Delta}\Bigg{(}1+f_{0}({\bf w})\int{\rm
d}^{6}{\bf w}_{a}~{}X({\bf x}_{a},{\bf x})\Bigg{)}$

where the integral,

$\mathcal{X}({\bf x})=\int{\rm d}^{3}{\bf x}_{a}~{}X({\bf x}_{a},{\bf x})$
(82)

is a solution of the ODE,

$\Bigg{(}{1\over 4\pi A}\nabla^{2}+\int{\rm d}^{3}{\bf
v}~{}f_{0}\Bigg{)}\mathcal{X}({\bf x})=-1.$ (83)

with boundary conditions (for the homogeneous solution) determined by the
boundary conditions of the problem. I will illustrate this with an example in
the next section.

##  8 Non-Zero Non-Fluctuations

Having calculated these corrections for the mean energy theory, we can now ask
a most crucial question: are these corrections which define the mean-field
significant?

If they are, then how we linearise and study the CBE must change, because the
ensemble average is modified by terms which are themselves a function of
phase-space.

Let us consider what might be the simplest non-trivial system in which these
corrections might matter: the self-gravitating Maxwellian. The Maxwellian is
well-known as a stable equilibria: Indeed, we see that there is no need to
impose a constraint to ‘choose’ $f_{0}$: Isothermal distributions arise
naturally when we expand (e.g. equation (49))

$P_{E}[f]=\exp(S_{g}-\beta E)$ (84)

about its saddle-point, the isothermal distribution defined by,

$\Bigg{(}{\delta S_{g}\over\delta f}-\beta{\delta E\over\delta
f}\Bigg{)}\Bigg{|}_{f=f_{0}}=0\implies f_{0}\propto\exp(-\beta H_{0}).$ (85)

In particular, I will consider a self-gravitating Maxwellian which is self-
consistent within a maximum radius $r_{m}$ from the origin, and is otherwise
‘frozen’ outside that radius, much in the same way that Zang (1976) ‘froze’
the Mestel disc. This can be done by considering the potential imposed on the
system from outside $r\leq r_{m}$ as an external potential.

The maximum entropy state is the Maxwellian, which we have normalised to
$\int{\rm d}^{3}{\bf v}~{}f_{0}=\rho_{0}$, which differs from our usual
normalisation of $\int{\rm d}^{6}{\bf w}~{}f_{0}=1$. This maps $A=\beta
GMm\rightarrow{A}=\beta Gm$.

$f_{0}={\rho_{0}\over(2\pi/\beta)^{3/2}}\exp(-\beta{{\textstyle{1\over
2}}}m{\bf v}^{2})$ (86)

For this system, we find that the spatial correlation function obtained by
substituting equation (86) into equation $\eqref{eq:Spatial}$ is:

$X(|{\bf x}-{\bf x}^{\prime}|)=A{\cos(k_{J}r)\over r};~{}k_{J}=\sqrt{4\pi
A\rho_{0}}$ (87)

where $k_{J}=2\pi/\lambda_{J}$ is the Jeans wavenumber and $\lambda_{J}$ is
the Jeans length.

Then $\mathcal{X}(\mathbf{0})$ is simply,

$\displaystyle\mathcal{X}(\mathbf{0})$ $\displaystyle=\int{\rm d}^{3}{\bf
x}^{\prime}~{}X(|{\bf x}^{\prime}|)$ (88)
$\displaystyle=4\pi\int_{0}^{r_{m}}{\rm d}r~{}Ar\cos(k_{J}r)$
$\displaystyle={4\pi A\over
k_{J}^{2}}\Bigg{(}k_{J}r_{m}\sin(k_{J}r_{m})+\cos(k_{J}r_{m})-1\Bigg{)}$
$\displaystyle={1\over\rho_{0}}\Bigg{(}k_{J}r_{m}\sin(k_{J}r_{m})+\cos(k_{J}r_{m})-1\Bigg{)}$

This lets us constrain $\mathcal{X}({\bf x})$, which obeys two boundary
conditions. The first is that $\mathcal{X}=\mathcal{X}(r)$, due to the
rotational symmetry of the problem. The second is that $\mathcal{X}(0)$ is
finite and takes on the calculated value. Applying these constraints gives us,

$\mathcal{X}(r)=C_{1}{\exp({\rm i}k_{J}r)\over r}+C_{2}{\exp(-ik_{J}r)\over
r}-{1\over\rho_{0}},$ (89)

$C_{1}=-C_{2}$ and

$\mathcal{X}(0)=C_{1}(2{\rm i}k_{J})-{1\over\rho_{0}}.$ (90)

Then,

$\mathcal{X}(r)=\Bigg{(}k_{J}r_{m}\sin(k_{J}r_{m})+\cos(k_{J}r_{m})\Bigg{)}{\sin(k_{J}r)\over\rho_{0}k_{J}r}-{1\over\rho_{0}}$
(91)

and the first order correction to $\langle\delta f\rangle_{f,E}$ is,

$\langle\delta f\rangle_{f_{0},E}={1\over
2N\Delta}\Bigg{(}1+{V_{v}}f_{0}\mathcal{X}({\bf x})\Bigg{)}+....$ (92)

where $V_{v}$ is the volume of velocity-space, $V_{v}=\int{\rm d}^{3}{\bf v}$.
Now define $N_{E}$,

$N_{E}={N\Delta\rho_{0}\over V_{v}}$ (93)

and we find that,

$\langle\delta f\rangle({\bf w})={1\over 2N_{f}}f_{0}+{1\over
2N_{E}}f_{0}(\mu(k_{J}r_{m}){\textrm{sinc}(k_{J}r)}-1)+...,$ (94)

where $N_{E}$ plays the role of an expansion parameter much like $N_{f}$ does,
and $\mu(x)=x\sin(x)+\cos(x)$ controls the size of the correction. Observe
that

$N_{E}/N_{f}={\rho_{0}}/V_{v}f_{0}\ll 1$ (95)

, implying that the gravitational interaction has blown up the ’small’
$1/N_{f}$ noise. This is roughly analogous to the dressing of Poisson noise.

This calculation induces a central-peak in the once spatially homogeneous
Maxwellian via the sinc function with wavelength equal to the Jeans length,
and appears to represent a system which has undergone growth via the linear
Jeans instability, and then nonlinearly come to saturation.

That $|\mu(k_{J}r_{m})|$ for most values of $r_{m}$ grows with $r_{m}$
supports this hypothesis, since more mass would be drawn under the influence
of the Jeans instability. The zeroes of $\mu(k_{J}r_{m})$ are nonlinearly
spaced for low radii, but approach a linear spacing for large $r_{m}$,
$k_{J}r_{m}=n\pi$ where is $n$ large, or simply that $2r_{m}=n\lambda_{J}$.
This appears to indicate that suitably sized self-consistent spheres are
capable of holding fluctuating standing waves. These fluctuations evade the
ensemble average that can only capture static features which are not removed
in the process of averaging.

We can infer that the ensemble average is only sensitive to corrections which
are non-zero on average (as we might expect), and this can lead to us not
seeing features when they are in fact dynamic and not static.

##  9 Discussion

###  9.1 Typicality and Shannon’s Entropy

Perhaps the most questionable (and most interesting) part of this theory lies
in its conception: in the construction of the typicality-based calculation
with which we justified the maximum entropy approach to obtain
$P[f]=\exp(S_{g}[f])$.

More specifically, in Section 3 I conjecture that all $f$ for which an
observed sample, $\\{{\bf w}_{i}\\}$ is perfectly typical are equally likely
to be the sampled $f$. This underpins the field theory.

Entropy maximisation with respect to $P[f]$ does not require that the Gibbs
entropy $S_{g}$ is maximised, and thus does not require the particles to
achieve equipartition. I instead describe how our lack of understanding with
respect to the ‘true’ $f$ necessitates that we assume the widest plausible
group of $f$ contributes equiprobably to the sample $\\{{\bf w}_{i}\\}$: that
is, the group of $f$ for which $\\{{\bf w}_{i}\\}$ is always not an outlier
are equally likely to be the ‘true’ $f$.

This is the only sensible choice to make in the absence of knowledge
differentiating between different $f$, and is what puts the Gibbs entropy
$S_{g}$ above the menagerie of generalised entropies (for the random sampling
problem we have here!).

The Gibbs entropy (but in the hands of Shannon) finds its place most naturally
in the $N$-particle distribution function—what Shannon would call the
probability of randomly sampling a sequence—where it arises as a natural
byproduct of simultaneously taking the logarithm and exponentiating a product
of individual sampling probabilities, and then using the law of large numbers.
(e.g. compose equations (36) and (38)).

There there are no axioms of what the Shannon entropy is, or what it should
do. It is merely a quantity which measures the expected probability of
sampling a sequence, given that each member of the sequence is sampled from
$f$. It defines what is typical, and for sufficiently large $N$ we learn that
probability space is dominated by typical samples. It seems that the Shannon
entropy is prevalent in physics because real systems which are of large $N$
tend to appear typical, and measurements reflect this.

As mentioned above, we do not need to maximise the Gibbs entropy here. Non-
perturbatively, it is integrated over as a part of $P[f]$. Perturbatively,
while we expand about $f_{0}$, we only do so to facilitate the Saddle Point
Approximation. However, one should note that the partition function
$\mathcal{Z}$ does not converge if $P[f]$ does not go to zero at
$f\rightarrow\infty$. This is usually an issue because it signifies that there
are states $f$ which do exist that occupy infinite probability $P[f]$. Not so
in our case, where $f$ is normalised to $1$.

###  9.2 $\langle f\rangle\neq f_{0}$

The standard program of fitting $f_{0}$ to a distribution is flawed when it
comes to the study of dynamical systems. Such an assumption is made whenever
we try to fit a distribution of stars, or a plasma density, with isothermal
$f_{0}$, Hernquist disks, Plummer spheres, and even exponential disks! The
error is made when we make claims regarding what the distribution function is,
as opposed to what it might be.

Instead of making statements regarding the DF, which we have no direct handle
on, we should instead be making statements regarding the state of an ensemble
of DFs, amongst which we have beliefs regarding which to favour; whether that
is a belief regarding a mean energy, or a mean angular momentum, or anything
else we would like.

Embracing this inherent uncertainty regarding our knowledge of the DF has to
be done through the space of distribution functions, $f$, which necessitates
the use of a field theory: a perturbative field theory in which we have the
freedom to vary $f({\bf w})$ at each point in phase-space. Owing to the
stochastic sourcing of this deviation, we find that the corrections to $f_{0}$
on $\langle f\rangle$ are represented by the expansion parameter $1/N_{f}$,
which describes the Poissonian fractional uncertainty of the number of
particles sampled within an element of phase-space, and is also the parameter
controlling the size of corrections to the CBE field theory. To leading order,
however, these corrections have no phase-space dependence and thus are
dynamically irrelevant.

When a mean energy is imposed, another such expansion parameter is produced:
$1/N_{E}>>1/N_{f}$, which controls the size of the corrections based on
correlations introduced by the Hamiltonian of the system. The relative sizes
of these expansion parameters show that the dressing of Poisson noise arises
because discreteness introduces fluctuations in phase-space, which are
compounded because force calculations between particles separated in position-
space sum the force contributions from fluctuations at all velocities within
each spatial element (see equation (95)).

###  9.3 Nonlinearity

The field theory does not rely on linearising the CBE; and though we split
$f=f_{0}+\delta f$ it must be understood that each $\delta f$ is a unique and
different $f$ from every other $\delta f$; that is we are studying ensembles
of independently realised distribution functions, not fluctuations on the same
distribution function. The difference here is that different fluctuations on
the same system must share the same momentum, energy, et cetera, while
different distribution functions can obviously differ in these regards.

This freedom is necessary because of the nature of random sampling. We know
that it is almost certain that randomly sampling a stationary, spherical and
isotropic distribution function like the isochrone distribution will produce a
discrete $N$-body system which is neither stationary, spherical nor isotropic.
Thus we must accept that the reverse is true too: a globular cluster with zero
velocity drift could be sampled from a distribution function which does drift
in velocities. We must allow ourselves to consider all distribution functions
$f$, of which the vast majority are out of equilibrium.

And yet we see in Section 4 that despite this, ensemble averaged distribution
functions for systems with physically meaningful constraints are in
equilibria. This result is hilarious and eye-opening: formalising the belief
that mean-fields should be in equilibria, while embracing the fact that almost
all distributions are out of equilibrium.

This is also a proof that isothermal distributions and exponential disks (that
is, isothermal and iso-angular-momentum), once dressed in stochastic
gravitational fluctuations, are nonlinear equilibria of the CBE, but only in
the sense of the mean field being static: meaning there are still time-varying
fluctuations ‘dancing’ around these systems. This is exciting because it
illustrates that well-structured noise can preserve the mean-field they
surround, and not induce secular evolution.

A perturbative calculation of the self-gravitating Maxwellian reveals that the
nonlinearly stable mean field which correponds to it bears signs of having
undergone Jeans collapse. I hypothesise that gravitational dressing is made
manifestly calculable through this apparatus! Further confirmation must be
done by comparison with simulation/data.

The field theory presents a quasi-nonlinear theory; it is quasi-nonlinear in
that we can only expand to some finite order in the expansion, but we know
that we are close to the true nonlinear solution.

###  9.4 Collisionality and Coarse-Graining

Regularising this field theory involves coarse-graining phase-space; i.e.
taking $\delta^{6}(0)\rightarrow 1/\Delta$.

This represents how a finite number of particles cannot possibly constrain the
uncountably infinite degrees of freedom captured by the continuous
distribution function. In Section 3 I described how it is not possible to
choose a coarse-graining grain size without an explicit belief regarding what
it should be. This belief was then asserted in Section 6 in the form of a
minimum grain-size based on the finiteness of $N$.

The imposition of a minimum grain-size was based on our need to truncate the
asymptotic series of $\langle f\rangle$, so as to produce finite, well-defined
outputs with the perturbative field theory. This result also makes statistical
sense: there is little reason to believe a single particle in a finite phase-
space element can constrain the value of $f$ at that point. Note, however,
that such a cut is unncessary: taking infinitesimal phase-space elements (i.e.
$\Delta\rightarrow 0$) is equivalent to taking the strong coupling limit of a
Quantum Field Theory, for which perturbation field theory does not apply,
since the expansions do not converge.

The concept of coarse-graining phase space given a sample $\\{{\bf w}_{i}\\}$
so as to reclaim $f$ is not new: see e Silva et al. (2019) for a salient
analysis on how we can use the Shannon-Nyquist (i.e. anti-aliasing) theorem to
set bounds on the fineness of distributions which can be uniquely attributed
to a sample. Their analysis concurs with ours: in our notation, roughly they
conclude that $1/\Delta\leq N$, or $1/N\Delta\leq 1$ is the limit, where
$\Delta$ is the phase-space volume of the smallest significant feature they
capture. Plasma physicists have also used the Debye sphere, a natural physical
scale at which the electrostatic interaction falls off due to electron
shielding, e.g. see Ewart et al. (2023) for a more careful treatment in the
context of Lynden-Bell statistics. That considerably larger $\Delta$ allows
for $1/N\Delta<<1$. But here we should ask ourselves: do we choose the theory
which admits as many choices of $\Delta$ as possible, or do we select a theory
with one, physically motivated $\Delta$?

Perhaps less well-motivated is a presumption I made when I introduced the mean
energy constraint in Section 8, and then specified that there are no
gravitational correlations within the same discrete phase-space element. While
this is most certainly true in the continuum limit, since single particles
occupying infinitesimal packets of phase-space cannot possibly gravitationally
enhance themselves due to the incompressibility of phase-space, discrete
phase-space elements contain a finite volume of phase-space, and so
infinitesimal packets of phase-space can enter and leave the volume. The
discrete phase-space element can thus harbour gravitational enhancements.

What I neglect are the local gravitational collisions which occur between
particles within a discrete phase-space element. Such local scatterings are
well described by Chandrasekhar (1949). Choosing a size for the discrete
phase-space element therefore requires deciding for oneself, at what point
does local scattering become less important than the Poisson noise one suffers
if there are only a few stars within the element? It can be seen that the more
we suppress Poisson noise through discretising over larger distances, the more
we end up neglecting scattering.

We might understand the relationship between the CBE and local scattering as
being one in which the CBE handles all long-distance relationships between the
finite phase-space elements, while local scattering handles the local
relationships between particles within the phase-space elements. One cannot
look at the Coulomb logarithm, which describes how local scattering
logarithmically diverges as one considers scatterings between a source
particle and particles which lie successively further from the origin, and
Fouvry et al. (2021) which features a logarithm which approaches the
divergence from global scattering, and not feel some hope that somewhere in
between they might meet to form a complete theory of collisions.

This is an avenue for future work, which I am quite eager to pursue.

###  9.5 The Sample Probability

Last, but certainly not least, we have the sample probability, $\mathcal{P}$
(equation (44)). Connecting the probability required to sample the sequence of
particles, $\\{{\bf w}_{i}\\}$ with the correlation functions that describe
how $f$ is correlated with itself between different points in phase-space,
proves that correlations within $f$ at the points $\\{{\bf w}_{i}\\}$ are
equal to the correlations between indistinguishable particles $\\{{\bf
w}_{i}\\}$.

This makes sense: if an ensemble of $f$ bears internal correlations, then an
ensemble of samples drawn from the ensemble of $f$ should share the same
internal correlations. This is a result which cannot be described by the
Liouville equation, which acts on individual $f$, and cannot incorporate
information regarding an ensemble of $f$.

##  10 Conclusions

I have presented a perturbative field theory which allows us to calculate the
mean distribution function, $\langle f\rangle$ and higher moments, with
respect to the nonlinear Boltzmann equation.

I show that there is a unique method to assigning probabilities to
distribution functions, $P[f]$ for the CBE, conditioned on Shannon’s
typicality condition. This allows us to calculate two-point correlations,
discovering how a distribution function is correlated with itself across
phase-space, via an ensemble average.

Defining the ensemble average is important because it allows us to employ the
field theory formalism. We are able to explore nonlinear equilibria, calculate
the saturated states of certain linear instabilities, and (I think) most
importantly, express long-term interactions like global gravitational dressing
in a sensible manner.

The next papers in this series will describe how we can use this theory to
describe secular evolution, and present results for rotating systems (think
discs!).

Further work should also include an application to the periodic cube: such an
application is just waiting to be completed, and the periodic cube is
promising as an instructional testbed for nonlinear effects. It might be a
long shot, but what interests me the most is the connection between this
theory and the extensive (but with non-additive entropies) statistical
mechanics of Tsallis’ (see Tsallis (2011) for a review). In parallel, I think
it is worthwhile to attempt introducing collision operators into this
currently collisionless theory.

I heartily welcome collaborators to aid in ironing out the creases in this
theory, and look forward to cooperating on more applications within and beyond
the astrophysical context!

## Acknowledgements

Jun Yan Lau gratefully acknowledges support from University College London’s
Overseas and Graduate Research Scholarships. I also gratefully acknowledge
support from my parents. I also thank Minghao Li, Robert Ewart, Thormund Tay,
Ralph Schoenrich, Kinwah Wu, James Binney, Walter Dehnen, Rimpei Chiba, John
Magorrian and Alexander Schekochihin for many insightful conversations, and
their heartening interest in my babbling.

## Data Availability

No new data was generated or analysed in support of this research.

## References

  * Antoja et al. (2018) Antoja T. et al., 2018, Nat, 561, 360
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Chandrasekhar (1949) Chandrasekhar S., 1949, Reviews of Modern Physics, 21, 383
  * Chiba et al. (2020) Chiba R., Friske J. K. S., SchÃ¶nrich R., 2020, Monthly Notices of the Royal Astronomical Society, 500, 4710
  * Chiba & Schönrich (2021) Chiba R., Schönrich R., 2021, Monthly Notices of the Royal Astronomical Society, 505, 2412
  * e Silva et al. (2019) e Silva L. B., de Siqueira Pedra W., Valluri M., 2019, The Astrophysical Journal, 872, 20
  * Eddington (1916) Eddington A. S., 1916, MNRAS, 76, 572
  * Ewart et al. (2023) Ewart R. J., Nastac M. L., Schekochihin A. A., 2023, Non-thermal particle acceleration and power-law tails via relaxation to universal lynden-bell equilibria
  * Feynman (1949) Feynman R. P., 1949, Phys. Rev., 76, 769
  * Fouvry et al. (2021) Fouvry J.-B., Hamilton C., Rozier S., Pichon C., 2021, Monthly Notices of the Royal Astronomical Society, 508, 2210
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, Astronomy and Astrophysics, arXiv:2208.00211
  * Henon (1959) Henon M., 1959, Annales d’Astrophysique, 22, 126
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * King (1966) King I. R., 1966, AJ, 71, 64
  * Lin & Shu (1964) Lin C. C., Shu F. H., 1964, ApJ, 140, 646
  * Magorrian (2006) Magorrian J., 2006, Monthly Notices of the Royal Astronomical Society, 373, 425
  * Merritt (1985) Merritt D., 1985, MNRAS, 214, 25P
  * Mestel (1963) Mestel L., 1963, Monthly Notices of the Royal Astronomical Society, 126, 553
  * Monari et al. (2019) Monari G., Famaey B., Siebert A., Wegg C., Gerhard O., 2019, A&A, 626, A41
  * Osipkov (1979) Osipkov L. P., 1979, Pis ma Astronomicheskii Zhurnal, 5, 77
  * Peskin & Schroeder (1995) Peskin M. E., Schroeder D. V., 1995
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Sellwood & Carlberg (2014) Sellwood J. A., Carlberg R. G., 2014, ApJ, 785, 137
  * Shannon (1948) Shannon C. E., 1948, The Bell System Technical Journal, 27, 379
  * Tolman (1938) Tolman R., 1938, The Principles of Statistical Mechanics. Oxford University Press
  * Tsallis (2011) Tsallis C., 2011, Entropy, 13, 1765
  * Wick (1950) Wick G. C., 1950, Phys. Rev., 80, 268
  * Zang (1976) Zang T. A., 1976, PhD thesis, Massachusetts Institute of Technology, Cambridge, Massachusetts