from the origin at some later time, it would have to first vanish. Also, since an extreme point has not formed yet,
$(r_0, t_0)$ is an inflection point satisfying $\frac{\partial^2 E}{\partial r^2}(r_0, t_0) = 0$ and
$\frac{\partial^3 E}{\partial r^3}(r_0, t_0) < 0$. Plugging these into our expression above, we see that
$\frac{\partial}{\partial t} \frac{\partial E}{\partial r} (r_0, t_0) = -\beta^{\prime}(E(r_0, t_0)) \frac{\partial^3 E}{\partial r^3}(r_0, t_0) < 0$.
If the radial derivative of $E$ ever gets to zero, it starts decreasing again and never passes zero to become positive.
Therefore the origin is the only local extremum at all times.

§.§.§ Existence of a steady-state

We now have everything we need to discuss the behaviour of solutions in the limit of infinite time. By setting
$\frac{\partial E}{\partial t}$ equal to zero, we can see that the limiting energy distribution satisfies

\begin{equation}
\begin{dcases}
\Delta \beta(E(x)) = 0 & x \in \Omega \\
\nabla E(x) \cdot n = 0 & x \in \partial \Omega
\end{dcases} \; . \nonumber
\end{equation}

The only solutions to the Neumann problem for Laplace's equation are constants. Therefore $\beta(E)$ must be
constant in space if $\frac{\partial E}{\partial t} = 0$. There are two ways for this to happen. One is for $E$ to
solely occupy the Hagedorn regime. If not, $\beta$ is invertible in a neighbourhood of at least one energy in the
steady-state and $E$ must be constant. By energy conservation, the value of this constant is of course

\begin{equation}
\bar{E_0} = \frac{1}{|\Omega|} \int_{\Omega} E_0(x) \textup{d}x \; . \nonumber
\end{equation}

In other words, we know what $\lim_{t \rightarrow \infty} E(\cdot, t)$ must be if it exists but it is not
yet obvious that it exists. Even for an energy distribution whose extrema stay in the same place and smooth out over
time, it is possible for the intermediate regions to constantly oscillate without ever converging to any function. Our
proof that this does not happen for (<ref>) will compare it to the heat equation.

The two energy distributions shown have equal minima, equal maxima and equal masses. They also have the
central peak as their only extreme value. A solution to (<ref>) that alternated between these
functions indefinitely would satisfy all of the properties that we have proven so far but not have a steady-state.

Suppose that $E$ solves

\begin{equation}
\begin{dcases}
\frac{\partial E}{\partial t} (x, t) = - \Delta \beta(E(x, t)) & (x, t) \in \Omega \times (0, \infty) \\
E(x, 0) = E_0(x) & x \in \Omega \\
\nabla E(x, t) \cdot n = 0 & (x, t) \in \partial \Omega \times (0, \infty)
\end{dcases} \nonumber
\end{equation}

for a spherically symmetric $E_0$ with a central maximum as its only local extreme point. If $\beta^{\prime}(E) \leq 0$
becomes strict whenever $E_{\mathrm{min}} < E < \bar{E_0}$, then $E$ converges to its average.

Letting $B^d_R$ denote a ball of radius $R$, the mass satisfies

\begin{eqnarray}
M_R &=& \int_{B^d_R} E \textup{d}x \nonumber \\
\frac{\textup{d} M_R}{\textup{d}t} &=& \int_{B^d_R} \frac{\partial E}{\partial t} \textup{d}x \nonumber \\
&=& -\int_{B^d_R} \Delta \beta(E) \textup{d}x \nonumber \\
&=& -\int_{\mathbb{S}^{d-1}_R} \beta^{\prime}(E) \nabla E \cdot n \textup{d}S_x \nonumber \\
&\leq& 0 \; . \nonumber
\end{eqnarray}

The inequality came from the fact that the integrand is non-negative. The flux is non-positive because the value of $E$
must always decrease as we move away from the origin and $\beta^{\prime}(E)$ is non-positive by assumption. Let
$R(t, \epsilon)$ be the largest radius at which $E(R, t) = \bar{E_0} - \epsilon$. If $R(t)$ ever reaches the radius of
$\Omega$, the solution will have converged to $\bar{E_0}$ by energy conservation and the maximum principle. Therefore
assume the opposite and let $R = \max \{ R(t, \epsilon) : t > 0 \}$. Since $E(R, t) \leq \bar{E_0}$, we can make the
integral larger by replacing $\beta^{\prime}(E)$ with $\beta^{\prime}(\bar{E_0})$.

\begin{equation}
\frac{\textup{d} M_R}{\textup{d}t} \leq -\beta^{\prime}(\bar{E_0}) \int_{\mathbb{S}^{d-1}_R} \nabla E \cdot n \textup{d}S_x < 0 \nonumber
\end{equation}

We now see that $M_R$ for our equation (<ref>) shrinks more quickly than $M_R$ for a heat equation
with diffusion constant $-\beta^{\prime}(\bar{E_0})$. The standard Fourier series method shows that the solution to
said heat equation with initial data $E_0$ will converge to $\bar{E_0}$. Thus $\mathbb{S}^{d-1}_R$, the sphere where
$\bar{E_0}$ is first achieved, eventually encloses a constant function. By energy conservation and the maximum
principle, $E$ must be equal to the same constant outside this sphere as well. Therefore for any $\epsilon > 0$, the
energy at the edge always comes within $\epsilon$ of $\bar{E_0}$ which means that “the outside of this sphere” must
have not existed all along.

§.§ Time scales

So far we have shown that the nonlinear diffusion we are interested in shares a number of intuitive properties with
linear diffusion: energy is conserved, peaks smooth out over time and reasonable configurations of energy will decay to
a distribution that is completely uniform. However, we still have relatively little information on how quickly
these energy distributions decay. Defining and bounding time scales for the diffusion will allow us to compare our
model to what is already known about gravitational bound states.

§.§.§ The concentration comparison theorem

This shows how we must change the function that appears inside the Laplacian in order to satisfy the
hypotheses of the theorem. Bounding the time scales arising from $\beta$ is no more difficult than bounding the
time scales arising from $\Phi$.

If we continued applying the basic properties, we would be able to derive some powerful results that are found in the
literature. One of these is the concentration comparison theorem [64] which applies to equations of the form
(<ref>) called filtration-equations. The only difference is that filtration equations
are typically written $\frac{\partial E}{\partial t} = \Delta \Phi(E)$ where $\Phi$ is weakly increasing. A
fundamental property of diffusion is that the mass contained within a fixed ball at the origin $M_R(t)$ should decrease.
The total mass $M(t)$ which we found to be constant in time is $M_{\infty}(t)$. The concentration comparison theorem,
which we state below, is used to convert an inequality involving $\Phi$ to an inequality involving $M_R$.

Let $\Phi_1, \Phi_2 \in C^1(\mathbb{R})$ be increasing functions sending $0$ to $0$ such that
$\Phi_1^{\prime} \geq \Phi_2^{\prime}$. Suppose that $u_1$ having mass $M_{1, R}(t)$ and $u_2$ having mass
$M_{2, R}(t)$ are solutions to $\frac{\partial u}{\partial t} = \Delta \Phi_1(u)$
and $\frac{\partial u}{\partial t} = \Delta \Phi_2(u)$ respectively with spherically symmetric data.
If $M_{1, R}(0) \leq M_{2, R}(0)$ for all $R$ then $M_{1, R}(t) \leq M_{2, R}(t)$ for all
$R$ and $t$.

To use this theorem, we must convert our function $\beta(E)$ into something satisfying $\Phi(0) = 0$. This is easily
done with $\Phi(E) = \beta(E_{\mathrm{min}}) - \beta(E + E_{\mathrm{min}})$. Figure <ref> shows what
kind of function $\Phi$ is if $\beta$ is the function in Figure <ref>.

We were able to show that our filtration equation had a well defined steady-state by comparing it to a heat
equation. With this linear function as our bounding function, an estimate for the decay time can also be obtained.
However, this would lead to a crude bound because the linear function above is not very similar to the true $\Phi$.

\begin{equation}
\frac{\partial}{\partial t} \left ( E + E_{\mathrm{min}} \right ) = -\Delta \beta \left ( E + E_{\mathrm{min}} \right ) \nonumber
\end{equation}

is now the same as saying

\begin{equation}
\frac{\partial E}{\partial t} = \Delta \Phi(E) \nonumber
\end{equation}

where $E_0$ is now some energy distribution equal to zero on $\partial \Omega$. If we can find an exact solution to
another filtration equation whose mass shrinks more quickly than that of $E$, we will have found a “decay time” that
is shorter than the one we are looking for. Similarly, a mass that shrinks more slowly would be associated with a longer
“decay time”. Time scales can therefore be determined if we choose bounding functions that are “steeper”
than $\Phi$ or “flatter” than $\Phi$. The crudest thing we could do is find a lower bound on the decay time by
drawing a linear function above $\Phi$. The next logical step is to find proper estimates for the
decay time by comparing our filtration equation to something more non-trivial.

§.§.§ Estimates in one dimension

As stated before, we are considering a density of states that has two microcanonical phase transitions. We will call
the energy of the first one $E_{\mathrm{H}}$ for Hagedorn and the energy of the second one $E_{\mathrm{F}}$ for field
theory. A useful definition of decay time for us will be the time required for the maximum of an energy distribution
to descend from $E_{\mathrm{F}}$ to $E_{\mathrm{H}}$. Specifically, we will look at the initial condition

\begin{equation}
E_0(x) =
\begin{cases}
E_{\mathrm{F}} & |x_1|, \dots, |x_d| < a \\
0 & \mathrm{otherwise}
\end{cases}
\label{toy-condition}
\end{equation}

and see how long it takes until $E(x, t) < E_{\mathrm{H}}$ for all values of $x$. The steep vertical jump through the
Hagedorn regime is a feature that (<ref>) has in common with sharply peaked initial energy distributions.
The sharply peaked functions do not need to be constant at $E_{\mathrm{F}}$ but this will happen anyway once they are
allowed to evolve. The distributions in Figure <ref> will flatten out in a relatively short time because
diffusion dominates above $E_{\mathrm{F}}$ but not below. This suggests that the decay time for a general peak is
dominated by the decay time of (<ref>). Further evidence that (<ref>) is representative of
more general initial conditions will appear in the next chapter on numerics.

Except for the piece on top that will quickly diffuse, these energy distributions are qualitatively similar
to the step function we are considering.

Finding the decay time of (<ref>) with respect to a general filtration function is still too hard, so we
will pick a simple filtration function with the goal of using the concentration comparison theorem. The one to pick is

\begin{equation}
\Phi(E) =
\begin{cases}
E & E < E_{\mathrm{H}} \\
E_{\mathrm{H}} & E \geq E_{\mathrm{H}}
\end{cases}
\; . \label{simplified-phi}
\end{equation}

If we take an energy profile solving $\frac{\partial E}{\partial t} = \Delta \Phi(E)$ and look at it at some instant
of time, the parts with energy above $E_{\mathrm{H}}$ should be static and the parts with energy below
$E_{\mathrm{H}}$ should be satisfying the heat equation. Energy crosses $E_{\mathrm{H}}$ at a particular distance
$x^*(t)$. Using the method of [2], we will construct a one-dimensional solution where $x^*$ starts off at $a$
and eventually shrinks to zero indicating that the decay time $T$ has been reached.

\begin{equation}
E(x, t) =
\begin{cases}
E_{\mathrm{F}} & |x| < x^*(t) \\
F(x, t) & |x| > x^*(t)
\end{cases}
\label{1d-solution-ansatz}
\end{equation}

where $F$ solves the regular heat equation. Clearly $F$ cannot be just any solution to the heat equation. To obey the
initial condition (<ref>), we must have $F(x, 0) = 0$ for $|x| > a$. Also, we must impose conservation
of energy. The mass contained between $0$ and $x^*(t)$ is $x^*(t) E_{\mathrm{F}}$. The mass everywhere else is
$\int_{x^*(t)}^{\infty} F(x, t) \textup{d}x$. We want their rates of change to be equal and opposite so

\begin{eqnarray}
E_{\mathrm{F}} \frac{\textup{d}x^*(t)}{\textup{d}t} &=& -\frac{\textup{d}}{\textup{d}t} \int_{x^*(t)}^{\infty} F(x, t) \textup{d}x \nonumber \\
&=& F(x^*(t), t) \frac{\textup{d}x^*(t)}{\textup{d}t} - \int_{x^*(t)}^{\infty} \frac{\partial F}{\partial t}(x, t) \textup{d}x \nonumber \\
&=& E_{\mathrm{H}} \frac{\textup{d}x^*(t)}{\textup{d}t} - \int_{x^*(t)}^{\infty} \Delta F(x, t) \textup{d}x \nonumber \\
&=& E_{\mathrm{H}} \frac{\textup{d}x^*(t)}{\textup{d}t} + \frac{\partial F}{\partial x} (x^*(t), t) \nonumber
\end{eqnarray}

where we have differentiated under the integral sign. Therefore

\begin{equation}
\left ( E_{\mathrm{F}} - E_{\mathrm{H}} \right ) \frac{\textup{d}x^*(t)}{\textup{d}t} = \frac{\partial F}{\partial x} (x^*(t), t) \label{1d-mass-matching}
\end{equation}

is a necessary condition for the ansatz (<ref>) to work. It is also sufficient as we show in the
appendix. The initial condition for $F$ that will give us these necessary and sufficient conditions is
$F(x, 0) = A \Theta(x - a)$ where $A$ is a yet undetermined constant that will turn out to be between $E_{\mathrm{H}}$
and $2E_{\mathrm{H}}$. Figure <ref> shows the basic setup.

[][$t = 0$]
[][Some later time]
The union of the red and purple curves is what we are looking for. It is (<ref>), the
solution to the (<ref>) filtration equation for the (<ref>) Cauchy data. The solution to
the auxillary problem that we use to find it is the union of the blue and purple curves. This is $F$, the solution to
the heat equation for the Heaviside Cauchy data.

To solve for $F$, we first note that the heat kernel is

\begin{equation}
\Phi_{\mathrm{H}}(x,t) = \frac{1}{(4\pi t)^{\frac{d}{2}}} e^{-\frac{x^2}{4t}} \; . \nonumber
\end{equation}

The heat equation's Cauchy problem is solved by taking the convolution of the initial condition with the heat kernel.

\begin{eqnarray}
F(x, t) &=& \frac{A}{\sqrt{4\pi t}} \int_{-\infty}^{\infty} \theta(y - a) e^{-\frac{(x - y)^2}{4t}} \textup{d}y \nonumber \\
&=& \frac{A}{\sqrt{4\pi t}} \int_{-\infty}^a e^{-\frac{(x - y)^2}{4t}} \textup{d}y \nonumber \\
&=& \frac{A}{2} \left ( 1 + \erf \left ( \frac{a - x}{2\sqrt{t}}\right ) \right ) \; . \label{1d-auxillary-solution}
\end{eqnarray}

Differentiating $F$ is now straightforward and the position of the interface can be found by setting
$F(x, t) = E_{\mathrm{H}}$. From this we obtain

\begin{eqnarray}
\frac{\partial F}{\partial x}(x, t) &=& -\frac{A}{\sqrt{4\pi t}} e^{-\frac{(x - a)^2}{4t}} \label{1d-f-derivative} \\
x^*(t) &=& a - 2\sqrt{t} \erf^{-1} \left ( \frac{2E_{\mathrm{H}}}{A} - 1 \right ) \label{1d-f-interface}
\end{eqnarray}

and the ratio between (<ref>) at the interface and the derivative of (<ref>) is constant.

\begin{equation}
\frac{\frac{\partial F}{\partial x}(x^*(t), t)}{\frac{\textup{d}x^*(t)}{\textup{d}t}} = \frac{A}{\sqrt{4\pi}} \frac{e^{-\erf^{-1} \left ( \frac{2E_{\mathrm{H}}}{A} - 1 \right )^2}}{\erf^{-1} \left ( \frac{2E_{\mathrm{H}}}{A} - 1 \right )} \nonumber
\end{equation}

We must set this equal to $E_{\mathrm{F}} - E_{\mathrm{H}}$. The $A$ we obtain from doing so is given by

\begin{eqnarray}
A &=& \frac{2 E_{\mathrm{H}}}{1 + \erf(I)} \in \left [ E_{\mathrm{H}}, 2E_{\mathrm{H}} \right ] \label{a-value} \\
\sqrt{\pi}Ie^{I^2}(1 + \erf(I)) &=& \frac{E_{\mathrm{H}}}{E_{\mathrm{F}} - E_{\mathrm{H}}} \; . \label{transcendental-equation}
\end{eqnarray}

The transcendental equation defining $I$ can be solved because the range of the left hand side includes all positive
real numbers. The last thing we need to do is find a more explicit value of $A$ for the case when
$E_{\mathrm{F}} \gg E_{\mathrm{H}}$. If we were to consider $\sqrt{\pi}Ie^{I^2}(1 + \erf(I)) = 0$, the solution would
simply be $I = 0$. Since $E_{\mathrm{H}}$ is close to zero but not quite, it is appropriate to linearize the left hand
side of (<ref>) giving us $\sqrt{\pi}I = \frac{E_{\mathrm{H}}}{E_{\mathrm{F}} - E_{\mathrm{H}}}$.
Plugging this into (<ref>) gives

\begin{equation}
A \approx \frac{2 E_{\mathrm{H}}}{1 + \erf \left ( \frac{E_{\mathrm{H}} / \sqrt{\pi}}{E_{\mathrm{F}} - E_{\mathrm{H}}} \right )} \; . \nonumber
\end{equation}

Finally, the time scale can be found by setting $x^*(T) = 0$ or $F(0, T) = E_{\mathrm{H}}$.

\begin{eqnarray}
T &=& \frac{1}{4} \left ( \frac{a}{\erf^{-1} \left ( \frac{2E_{\mathrm{H}}}{A} - 1 \right )} \right )^2 \nonumber \\
&\approx& \frac{1}{4} \left ( \frac{a \left ( E_{\mathrm{F}} - E_{\mathrm{H}} \right )}{E_{\mathrm{H}} / \sqrt{\pi}} \right )^2 \nonumber \\
&\approx& \frac{\pi}{4} \left [ \frac{a E_{\mathrm{F}}}{E_{\mathrm{H}}} \right ]^2 \label{1d-timescale}
\end{eqnarray}

§.§.§ Higher dimensional generalization

The procedure above works in an arbitrary number of dimensions if the initial condition has energy $E_{\mathrm{F}}$
inside a hypercube of side length $2a$. We will take $F(x, 0) = A \theta(x_1 - a) \dots \theta(x_d - a)$, solve the
heat equation for $F$ with this Cauchy data and write

\begin{equation}
E(x, t) =
\begin{cases}
E_{\mathrm{F}} & x_1 < x_1^*(t), \dots, x_d < x_d^*(t) \\
F(x, t) & \mathrm{otherwise}
\end{cases}
\; . \label{nd-solution-ansatz}
\end{equation}

Performing the convolution with the heat kernel is straightforward because the solution has a product form.

\begin{eqnarray}
F(x, t) &=& \frac{A}{(4\pi t)^{\frac{d}{2}}} \int_{-\infty}^{\infty} \dots \int_{-\infty}^{\infty} \theta(y_1 - a) \dots \theta(y_d - a) e^{-\frac{(x - y)^2}{4t}} \textup{d}y_1 \dots \textup{d}y_d \nonumber \\
&=& \frac{A}{(4\pi t)^{\frac{d}{2}}} \int_{-\infty}^a \dots \int_{-\infty}^a e^{-\frac{(x_1 - y_1)^2 + \dots + (x_d - y_d)^2}{4t}} \textup{d}y_1 \dots \textup{d}y_d \nonumber \\
&=& \frac{A}{2^d} \left ( 1 + \erf \left ( \frac{a - x_1}{2\sqrt{t}}\right ) \right ) \dots \left ( 1 + \erf \left ( \frac{a - x_d}{2\sqrt{t}}\right ) \right ) \label{nd-auxillary-solution}
\end{eqnarray}

From this we see that the surface $x^*(t)$ where $F$ takes on the value $E_{\mathrm{H}}$ is the locus of points

\begin{equation}
x^*(t) = \left [
\begin{tabular}{c}
$a - 2 \sqrt{t} \erf^{-1} \left ( z_1 - 1 \right )$ \\
\vdots \\
$a - 2 \sqrt{t} \erf^{-1} \left ( z_d - 1 \right )$
\end{tabular}
\right ] \label{nd-f-interface}
\end{equation}

where $\prod_i z_i = 2^d \frac{E_{\mathrm{H}}}{A}$. Differentiating $F$ and evaluating it on this surface, we get

\begin{equation}
\frac{\partial F}{\partial x_i}(x^*(t), t) = \frac{2E_{\mathrm{H}}}{z_i} \frac{1}{\sqrt{4\pi}} e^{-\erf^{-1} (z_i - 1)^2} \; . \label{nd-f-derivative}
\end{equation}

The analogue of (<ref>) to the higher dimensional case is

\begin{equation}
\left ( E_{\mathrm{F}} - E_{\mathrm{H}} \right ) \frac{\textup{d}V^*(t)}{\textup{d}t} = \int_{x^*(t)} \nabla F(x, t) \cdot n \textup{d}S_x \label{nd-mass-matching}
\end{equation}

where $V^*(t)$ is $a^d$ minus the volume enclosed by the surface (<ref>) in the positive orthant. This
can be visualized in Figure <ref>. Clearly $V^*(0) = 0$. If we want to calculate the volume
$V^*(t)$, we are looking for $x_i$ values that go from the $x_i^*(t)$ curve to $a$ instead of from
$0$ to the $x_i^*(t)$ curve. This is the same as saying that

\begin{eqnarray}
z_1 &\in& \left [ 1, 2^d \frac{E_{\mathrm{H}}}{A} \right ] \nonumber \\
z_2 &\in& \left [ 1, \frac{2^d}{z_1} \frac{E_{\mathrm{H}}}{A} \right ] \nonumber \\
&\vdots& \nonumber \\
z_d &\in& \left [ 1, \frac{2^d}{z_1 \dots z_{d-1}} \frac{E_{\mathrm{H}}}{A} \right ] \; . \nonumber
\end{eqnarray}

[][Early time]
[][Later time]
The curve where the blue sheet intersects the evolving surface is $x^*(t)$. The area of the blue sheet between
the surface and the co-ordinate planes is $a^2 - V^*(t)$. We are instead calculating the area $V^*(t)$ which has an
equal and opposite rate of change.

Each $x_i$ only depends on one $z_i$ so it is trivial to find the Jacobian determinant.

\begin{eqnarray}
\frac{\textup{d}x_i}{\textup{d}z_i} &=& -\sqrt{\pi t} e^{\erf^{-1}(z_i - 1)^2} \nonumber \\
\left | \frac{\textup{d}x_1}{\textup{d}z_1} \dots \frac{\textup{d}x_d}{\textup{d}z_d} \right | &=& (\pi t)^{\frac{d}{2}} \prod_{i=1}^d e^{\erf^{-1}(z_i - 1)^2} \nonumber
\end{eqnarray}

Also, by multiplying all but one of these eigenvalues together, we can find the area element required to calculate the
surface integral in (<ref>). An expression for the volume is

\begin{eqnarray}
V^*(t) &=& \int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 1}} \frac{E_{\mathrm{H}}}{A}} \left | \frac{\textup{d}x_1}{\textup{d}z_1} \dots \frac{\textup{d}x_d}{\textup{d}z_d} \right | \textup{d}z_d \dots \textup{d}z_1 \nonumber \\
&=& (\pi t)^{\frac{d}{2}} \int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 1}} \frac{E_{\mathrm{H}}}{A}} \prod_{i=1}^d e^{\erf^{-1}(z_i - 1)^2} \textup{d}z_d \dots \textup{d}z_1 \; . \label{nd-denominator}
\end{eqnarray}

We only need to know the length of $\nabla F$ to compute the surface integral. Working this out,

\begin{eqnarray}
\int_{x^*(t)} \nabla F(x, t) \cdot n \textup{d}S_x &=& \int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 2}} \frac{E_{\mathrm{H}}}{A}} \left | \nabla F(x^*(t), t) \right | \left | \frac{\textup{d}x_1}{\textup{d}z_1} \dots \frac{\textup{d}x_{d-1}}{\textup{d}z_{d-1}} \right | \textup{d}z_{d-1} \dots \textup{d}z_1 \nonumber \\
&=& \frac{1}{\sqrt{4\pi t}} (\pi t)^{\frac{d-1}{2}} \int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 1}} \frac{E_{\mathrm{H}}}{A}} 2E_{\mathrm{H}} \sqrt{\sum_{i=1}^d \frac{1}{z_i^2} e^{-2 \erf^{-1}(z_i - 1)^2}} \nonumber \\
&& \prod_{i=1}^{d-1} e^{\erf^{-1}(z_i - 1)^2} \textup{d}z_{d-1} \dots \textup{d}z_1 \label{nd-numerator}
\end{eqnarray}

where we have arbitrarily chosen $z_d$ as the variable to be determined from $z_1, \dots, z_{d-1}$.
The fact that (<ref>) is proportional to one more factor of $t$ than (<ref>) is what allows
(<ref>) to be obeyed. There is a problem with these expressions however that only rears its head when
$d > 1$; not all points on the surface $x^*(t)$ reach the co-ordinate planes at the same time. In two dimensions
for instance, the endpoint

\begin{equation}
\left [
\begin{tabular}{c}
$a - 2 \sqrt{t} \erf^{-1} \left ( \frac{4E_{\mathrm{H}}}{A} - 1 \right )$ \\
\end{tabular}
\right ] \nonumber
\end{equation}

will reach $(0, a)$ before the midpoint

\begin{equation}
\left [
\begin{tabular}{c}
$a - 2 \sqrt{t} \erf^{-1} \left ( \sqrt{\frac{4E_{\mathrm{H}}}{A}} - 1 \right )$ \\
$a - 2 \sqrt{t} \erf^{-1} \left ( \sqrt{\frac{4E_{\mathrm{H}}}{A}} - 1 \right )$
\end{tabular}
\right ] \nonumber
\end{equation}

reaches $(0, 0)$. In other words, for late times, the volume (<ref>) includes some regions outside the
positive orthant. It is not correct to integrate $x_i$ values from $a$ to the $x^*(t)$ curve. We must integrate from
$a$ to $\max \left ( x^*(t), 0 \right )$. The decay time obtained by setting (<ref>) to
$E_{\mathrm{H}}$ is

\begin{equation}
T = \frac{1}{4} \left [ \frac{a}{\erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d}} - 1 \right )} \right ]^2 \; . \nonumber
\end{equation}

The expression (<ref>) is most accurate when $T$ is large, i.e. when
$A \approx 2^d E_{\mathrm{H}}$. If this is true, we can obtain a reasonable approximation to this complicated integral
by setting $V^*(T) \approx a^d$. Doing this, we see that

\begin{equation}
\int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 1}} \frac{E_{\mathrm{H}}}{A}} \prod_{i=1}^d e^{\erf^{-1}(z_i - 1)^2} \textup{d}z_d \dots \textup{d}z_1 \approx \left ( \frac{2}{\sqrt{\pi}} \right )^d \erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d}} - 1 \right )^d \; . \nonumber
\end{equation}

With this, $E_{\mathrm{F}} - E_{\mathrm{H}}$ becomes much easier to express.

\begin{eqnarray}
E_{\mathrm{F}} - E_{\mathrm{H}} &=& \frac{\int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 2}} \frac{E_{\mathrm{H}}}{A}} 2E_{\mathrm{H}} \sqrt{\sum_{i=1}^d \frac{1}{z_i^2} e^{-2 \erf^{-1}(z_i - 1)^2}} \prod_{i=1}^{d-1} e^{\erf^{-1}(z_i - 1)^2} \textup{d}z_{d-1} \dots \textup{d}z_1}{\pi d \int_1^{2^d \frac{E_{\mathrm{H}}}{A}} \dots \int_1^{\frac{2^d}{z_1 \dots z_{d - 1}} \frac{E_{\mathrm{H}}}{A}} \prod_{i=1}^d e^{\erf^{-1}(z_i - 1)^2} \textup{d}z_d \dots \textup{d}z_1} \nonumber \\
&\approx& \frac{2E_{\mathrm{H}} \sqrt{d} \left ( \frac{2}{\sqrt{\pi}} \right )^{d-1} \erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d-1}} - 1 \right )^{d-1}}{\pi d \left ( \frac{2}{\sqrt{\pi}} \right )^2 \erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d}} - 1 \right )^d} \nonumber \\
&\approx& \frac{E_{\mathrm{H}}}{\sqrt{\pi d} \erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d}} - 1 \right )} \left [ \lim_{\frac{2^d E_{\mathrm{H}}}{A} \rightarrow 1} \frac{\erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d-1}} - 1 \right )}{\erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d}} - 1 \right )} \right ]^{d - 1} \nonumber \\
&=& \frac{E_{\mathrm{H}}}{\sqrt{\pi d} \erf^{-1} \left ( 2 \left ( \frac{E_{\mathrm{H}}}{A} \right )^{\frac{1}{d}} - 1 \right )} \left ( \frac{d}{d - 1} \right )^{d - 1} \nonumber
\end{eqnarray}

In the first step, we have used the fact that
$\sqrt{\sum_{i=1}^d \frac{1}{z_i^2} e^{-2 \erf^{-1}(z_i - 1)^2}} \approx \sqrt{d}$
because $z_1, \dots, z_d$ are very close to $1$ while in the last step, we have used l'Hôpital's rule. Linearizing
the inverse error function in the denominator and solving for $A$, we see that

\begin{equation}
\frac{2^d E_{\mathrm{H}}}{A} = \left [ 1 + \frac{E_{\mathrm{H}}}{E_{\mathrm{F}} - E_{\mathrm{H}}} \frac{1}{\sqrt{\pi d}} \left ( \frac{d}{d - 1} \right )^{d - 1} \right ] \approx 1 \; , \nonumber
\end{equation}

which is precisely the condition we needed in the first place to be able to say that $T$ was large and use this
approximation. Therefore when $E_{\mathrm{F}} \gg E_{\mathrm{H}}$, the time scale $T$ in $d$ dimensions becomes

\begin{eqnarray}
T &=& \frac{\pi d}{4} \left [ a \frac{E_{\mathrm{F}} - E_{\mathrm{H}}}{E_{\mathrm{H}}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \nonumber \\
&\approx& \frac{\pi d}{4} \left [ \frac{a E_{\mathrm{F}}}{E_{\mathrm{H}}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \; . \label{nd-timescale}
\end{eqnarray}

§.§.§ Upper and lower bounds

It is time to leverage these results to  improve the decay time estimates for the filtration function in Figure
<ref>. To achieve the tightest possible bounds, we will make sure that our more-steep and less-steep
functions flatten out at $E_{\mathrm{H}}$. Figure <ref> shows what this
looks like. The function on the left $\Phi_1$ has a slope $\Phi_1^{\prime}(0)$ equal to $\Phi^{\prime}(0)$.
The function on the left $\Phi_2$ has a slope $\Phi_2^{\prime}(0)$ equal to
$\lim_{E \rightarrow E_{\mathrm{H}}^{-}} \Phi^{\prime}(E)$. When these slopes are
determined, the concentration comparison theorem tells us that our time scale will satisfy

\begin{equation}
\frac{\pi d}{4 \Phi_1^{\prime}(0)} \left [ \frac{a E_{\mathrm{F}}}{E_{\mathrm{H}}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \leq T \leq \frac{\pi d}{4 \Phi_2^{\prime}(0)} \left [ \frac{a E_{\mathrm{F}}}{E_{\mathrm{H}}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \; . \label{prototype-timescale}
\end{equation}

[][Comparing $\Phi$ to $\Phi_1$]
[][Comparing $\Phi$ to $\Phi_2$]
The piecewise-linear bounding functions must relate to our filtration function in this way. On the left,
$\Phi_1$ could cross over to Hagedorn behaviour at some energy $E_1 \geq E_{\mathrm{H}}$ and on the right, $\Phi_2$
could cross over to Hagedorn behaviour at some energy $E_2 \leq E_{\mathrm{H}}$. The tightest possible time scale
bounds are achieved when $E_1 = E_{\mathrm{H}} = E_2$.

This is the main part of the text that will use a specific formula for $\beta(E)$, so first consider the one that is
valid when $C(E) = \rho^{-2}(E)$. Recall that this is
$\beta(E) = cE^{\alpha - 1} = \beta \left ( E_{\mathrm{min}} \right )E_{\mathrm{min}}^{1 - \alpha} E^{\alpha - 1}$
in the lowest energy regime. The slopes that we must insert are simple:

\begin{eqnarray}
\Phi_1^{\prime}(0) &=& -\beta^{\prime} \left ( E_{\mathrm{min}} \right ) = (1 - \alpha)\frac{\beta \left ( E_{\mathrm{min}} \right )}{E_{\mathrm{min}}} \label{slope-1} \\
\Phi_2^{\prime}(0) &=& -\beta^{\prime} \left ( E_{\mathrm{H}} \right ) = (1 - \alpha) \beta \left ( E_{\mathrm{min}} \right ) \frac{E_{\mathrm{H}}^{\alpha - 2}}{E_{\mathrm{min}}^{\alpha - 1}} \; . \label{slope-2}
\end{eqnarray}

Plugging (<ref>) and (<ref>) into (<ref>), this time scale is now solved:

\begin{equation}
\frac{\pi d}{4(1 - \alpha) \beta \left ( E_{\mathrm{min}} \right )} \frac{E_{\mathrm{min}}}{E_{\mathrm{H}}^2} \left [ a E_{\mathrm{F}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \leq T \leq \frac{\pi d}{4(1 - \alpha) \beta \left ( E_{\mathrm{min}} \right )} \frac{E_{\mathrm{min}}^{\alpha - 1}}{E_{\mathrm{H}}^{\alpha}} \left [ a E_{\mathrm{F}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \; . \label{timescale-1}
\end{equation}

This dependence on $\left ( a E_{\mathrm{F}} \right )^2$ is a general feature of long Hagedorn regions as we will see
in our next time scale.

Recall that when we had $C(E) = 1$, the effective inverse temperature function to consider $\tilde{\beta}$ was
defined by

\begin{eqnarray}
\tilde{\beta}^{\prime}(E) &=& \rho^2(E) \beta^{\prime}(E) \; . \nonumber
\end{eqnarray}

We will construct a filtration function from this in the same way as before.
$\tilde{\Phi}(E) = \tilde{\beta} \left ( E_{\mathrm{min}} \right ) - \tilde{\beta} \left ( E + E_{\mathrm{min}} \right )$.
In the low energy regime, we have

\begin{eqnarray}
\tilde{\Phi}^{\prime}(E) &=& -\rho^2 \left ( E + E_{\mathrm{min}} \right ) \beta^{\prime} \left ( E + E_{\mathrm{min}} \right ) = \rho^2 \left ( E + E_{\mathrm{min}} \right ) \Phi^{\prime}(E) \nonumber \\
&=& (1 - \alpha) \beta \left ( E_{\mathrm{min}} \right ) e^{\frac{2}{\alpha} \frac{\beta \left ( E_{\mathrm{min}} \right ) \left ( E + E_{\mathrm{min}} \right )^{\alpha}}{E_{\mathrm{min}}^{\alpha - 1}}} \frac{\left ( E + E_{\mathrm{min}} \right )^{\alpha - 2}}{E_{\mathrm{min}}^{\alpha - 1}} \; . \nonumber
\end{eqnarray}

The first thing to notice is that $E_{\mathrm{H}}$ is the same for $\Phi$ and $\tilde{\Phi}$. One derivative vanishes
if and only if the other does. Also, at an energy of $E_{\mathrm{min}}$, the number of states is very close to $1$
so $\tilde{\Phi}^{\prime}(0) = \Phi^{\prime}(0)$. Finding the slopes of appropriate bounding functions
$\tilde{\Phi}_1$ and $\tilde{\Phi}_2$ comes down to extremizing $\tilde{\Phi}^{\prime}$.

[][$E_{\mathrm{H}} < E^*$]
[][$E_{\mathrm{H}} > E^*$]
The upper bound for the decay time, given by $\tilde{\Phi}_2$ depends on which of the above cases is
realized. The figure on the left would lead to a situation much like the one for $\Phi_2$ that we saw before.
However, on the right, the Hagedorn energy is large enough to make $\tilde{\Phi}$ change concavity before reaching
the Hagedorn phase.

\begin{eqnarray}
\tilde{\Phi}^{\prime\prime}(E^*) &=& 0 \nonumber \\
2 \frac{\beta \left ( E_{\mathrm{min}} \right )}{E_{\mathrm{min}}^{\alpha - 1}} (1 - \alpha) E_{\mathrm{min}}^{\alpha - 1} \beta \left ( E_{\mathrm{min}} \right ) \left ( E^* + E_{\mathrm{min}} \right )^{2\alpha - 3} &=& (2 - \alpha) (1 - \alpha) E_{\mathrm{min}}^{\alpha - 1} \beta \left ( E_{\mathrm{min}} \right ) \left ( E^* + E_{\mathrm{min}} \right )^{\alpha - 3} \nonumber \\
E^* &=& \left [ \frac{2 - \alpha}{2 \beta \left ( E_{\mathrm{min}} \right ) E_{\mathrm{min}}^{1 - \alpha}} \right ]^{\frac{1}{\alpha}} - E_{\mathrm{min}} \nonumber
\end{eqnarray}

Therefore $\tilde{\Phi}$ starts off with a steep derivative at $0$ and becomes less steep until reaching the
inflection point $E^*$. This means that there are two important cases to consider: $E_{\mathrm{H}} < E^*$ and
$E_{\mathrm{H}} > E^*$. The flattest point occurs at $E_{\mathrm{H}}$ in the former case and at $E^*$ in the latter.

\begin{eqnarray}
\tilde{\Phi}_1^{\prime}(0) &=& (1 - \alpha)\frac{\beta \left ( E_{\mathrm{min}} \right )}{E_{\mathrm{min}}} \label{slope-3} \\
\tilde{\Phi}_2^{\prime}(0) &=&
\begin{cases}
(1 - \alpha) \beta \left ( E_{\mathrm{min}} \right ) \frac{E_{\mathrm{H}}^{\alpha - 2}}{E_{\mathrm{min}}^{\alpha - 1}} e^{\frac{2}{\alpha} \frac{\beta \left ( E_{\mathrm{min}} \right ) E_{\mathrm{H}}^{\alpha}}{E_{\mathrm{min}}^{\alpha - 1}}} & E_{\mathrm{H}} < E^* \\
(1 - \alpha) \left ( \frac{2 - \alpha}{2} \right )^{\alpha - \frac{2}{\alpha}} \beta \left ( E_{\mathrm{min}} \right )^{\frac{2}{\alpha}} E_{\mathrm{min}}^{\frac{2}{\alpha} - 2} e^{\frac{2}{\alpha} - 1} & E_{\mathrm{H}} > E^*
\end{cases} \label{slope-4}
\end{eqnarray}

Plugging (<ref>) and (<ref>) into (<ref>), the next time scale is solved as well:

\begin{equation}
\frac{\pi d}{4(1 - \alpha) \beta \left ( E_{\mathrm{min}} \right )} \frac{E_{\mathrm{min}}}{E_{\mathrm{H}}^2} \left [ a E_{\mathrm{F}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \leq T \leq C_T \left ( E_{\mathrm{min}}, E_{\mathrm{H}} \right ) \left [ a E_{\mathrm{F}} \left ( \frac{d - 1}{d} \right )^{d - 1} \right ]^2 \label{timescale-2}
\end{equation}

\begin{equation}
C_T \left ( E_{\mathrm{min}}, E_{\mathrm{H}} \right ) = \begin{cases}
\frac{\pi d}{4(1 - \alpha) \beta \left ( E_{\mathrm{min}} \right )} \frac{E_{\mathrm{min}}^{\alpha - 1}}{E_{\mathrm{H}}^{\alpha}} e^{-\frac{2}{\alpha} \frac{\beta \left ( E_{\mathrm{min}} \right ) E_{\mathrm{H}}^{\alpha}}{E_{\mathrm{min}}^{\alpha - 1}}} & E_{\mathrm{H}} < \left [ \frac{2 - \alpha}{2 \beta \left ( E_{\mathrm{min}} \right ) E_{\mathrm{min}}^{1 - \alpha}} \right ]^{\frac{1}{\alpha}} \\
\frac{\pi d}{4(1 - \alpha)} \left ( \frac{2 - \alpha}{2} \right )^{\frac{2}{\alpha} - \alpha} \beta \left ( E_{\mathrm{min}} \right )^{-\frac{2}{\alpha}} \frac{E_{\mathrm{min}}^{2 - \frac{2}{\alpha}}}{E_{\mathrm{H}}^2} e^{1 - \frac{2}{\alpha}} & E_{\mathrm{H}} > \left [ \frac{2 - \alpha}{2 \beta \left ( E_{\mathrm{min}} \right ) E_{\mathrm{min}}^{1 - \alpha}} \right ]^{\frac{1}{\alpha}}
\end{cases} \; . \nonumber
\end{equation}

There is a third possibility; $E_{\mathrm{H}}$ being so large that $\Phi^{\prime} \left ( E_{\mathrm{H}} \right )$
surpasses $\Phi^{\prime}(0)$. However, $E_{\mathrm{min}}$ can always be made small enough to stop this from happening.

§.§.§ Remaining problems

A useful question to ask is how our decay time differs from a decay time that would arise from linear diffusion.
Solivng the heat equation with the (<ref>) initial data,

\begin{eqnarray}
E(x, t) &=& \frac{E_{\mathrm{F}}}{(4\pi t)^{\frac{d}{2}}} \int_{-a}^a \dots \int_{-a}^a e^{-\frac{(x-y)^2}{4t}} \textup{d}y_1 \dots \textup{d}y_d \nonumber \\
&=& \frac{E_{\mathrm{F}}}{2^d} \prod_{i = 1}^d \left ( \erf \left ( \frac{x_i + a}{2 \sqrt{t}} \right ) - \erf \left ( \frac{x_i - a}{2 \sqrt{t}} \right ) \right ) \; . \nonumber
\end{eqnarray}

The time needed for the central peak to come down to $E_{\mathrm{H}}$ is

\begin{equation}
T = \frac{a^2}{4} \erf^{-1} \left ( \frac{E_{\mathrm{H}}}{E_{\mathrm{F}}} \right )^{-\frac{2}{d}} \approx \frac{a^2}{4} \left ( \frac{2E_{\mathrm{F}}}{\sqrt{\pi}E_{\mathrm{H}}} \right )^{\frac{2}{d}} \; . \nonumber
\end{equation}

For $d > 1$, this decay time is parametrically less than the $O \left ( \left ( a E_{\mathrm{F}} \right )^2 \right )$
result that we found for diffusion with a Hagedorn regime. However, linear diffusion in $d = 1$ is much slower.
To get a sense of why the $d = 1$ times are similiar, we will look at the amount of energy outside $B^1_a$ after a
short amount of time [2]. For diffusion with a Hagedorn regime, we find

\begin{eqnarray}
a E_{\mathrm{F}} - M_a(t) &=& A \int_a^{\infty} 1 + \erf \left ( \frac{a - x}{2 \sqrt{t}} \right ) \textup{d}x \nonumber \\
&=& 2A \sqrt{\frac{t}{\pi}} \nonumber \\
&\approx& 4 E_{\mathrm{H}} \sqrt{\frac{t}{\pi}} \; . \nonumber
\end{eqnarray}

For the heat equation with the same slope, the result is

\begin{eqnarray}
a E_{\mathrm{F}} - M_a(t) &=& E_{\mathrm{F}} \int_a^{\infty} \erf \left ( \frac{x + a}{2 \sqrt{t}} \right ) - \erf \left ( \frac{x - a}{2 \sqrt{t}} \right ) \textup{d}x \nonumber \\
&=& 2 E_{\mathrm{F}} \sqrt{\frac{t}{\pi}} \left ( 1 - e^{-\frac{a^2}{t}} \right ) + 2 E_{\mathrm{F}} a \left ( 1 - \erf \left ( \frac{a}{\sqrt{t}} \right ) \right ) \nonumber \\
&\approx& 2 E_{\mathrm{F}} \sqrt{\frac{t}{\pi}} \nonumber
\end{eqnarray}

where in the last step we have assumed $t \ll a^2$. Therefore, diffusion without a Hagedorn regime is faster at the
beginning, but not over longer time scales. We can explain this effect by noticing that a lump of energy diffusing
linearly becomes remarkably flat towards the “end” of its diffusion. Small spatial derivatives lead to small time
derivatives and a long time constant. Diffusion with a Hagedorn regime does not have this problem. The discontinuous
solutions we saw had steep profiles during all stages of what we called the diffusion. A different effect, the large
static phase, is what lengthened the time constant. Even though diffusion times are similar with and without a Hagedorn
phase, the numerics will show that if the energy distribution begins with a peak well above $E_{\mathrm{F}}$, its
diffusion will slow down significantly upon reaching $E_{\mathrm{F}}$. This does not contradict our analysis because
the diffusion in the high energy phase comes with a much larger constant than the diffusion that crosses over to
Hagedorn behaviour in the other two phases.

Regardless of how many dimensions there are, we should expect that our model overestimates the decay time of a plasma
ball. The energy in our model can only diffuse away to infinity so quickly once it leaves a ball of radius $a$.
In a realstic plasma ball, energy escaping the central region is quickly ejected because it is travelling through a

§.§ Comments on unbounded domains

For unbounded domains such as the whole space, the only sensible energy distributions are the ones that decay to zero
and filtration equations like (<ref>) take on many new properties. Even though the
diffusion of energy governed by a master equation should certainly be well defined for a space that extends infinitely,
we found in the last chapter that our PDE is only a useful model for such behaviour when the energy is high. Therefore,
it is not necessarily true that a result derived for $\Omega = \mathbb{R}^d$ applies to the entropic limit of a field

§.§.§ Barriers to uniqueness

The first thing we note is that

\begin{equation}
\begin{dcases}
\frac{\partial E}{\partial t} (x, t) = - \Delta \beta(E(x, t)) & (x, t) \in \mathbb{R}^d \times (0, \infty) \\
E(x, 0) = E_0(x) & x \in \mathbb{R}^d
\end{dcases} \nonumber
\end{equation}

would have many solutions if we simply looked for them in the set of all differentiable functions. This property is
well known for the heat equation. In fact, for the heat equation [65] and similar equations [66]
in the whole space, nonzero solutions have been found that vanish at $t = 0$. These authors have found a restriction
that must be made in order to recover uniqueness; there must exist a paraboloid in $\mathbb{R}^d$ such that
$\log E$ stays below it at all times. Because our equation is different from the heat equation, it is not yet clear
that the same condition, or a stronger one, will yield unique solutions to (<ref>). Nevertheless,
we will impose the stronger condition that $E$ be uniformly bounded in spacetime.

Requiring uniform boundedness is natural because it forces solutions to still obey the maximum principle. The maximum
principle phrased in terms of subsolutions and supersolutions made no assumption about the domain being bounded,
but in applying it, we set one solution equal to the maximum of the other. Our assumption of uniform boundedness
is the only reason that we can still write $u = \max \{ v(x, 0) : x \in \mathbb{R}^d \}$ because the maximum of a
continuous function is only guaranteed to exist on a compact domain. Having to make assumptions about how solutions
grow is a common theme when extending maximum principles to the unbounded case [61]. We now turn to the
question of uniqueness which is intimately related to energy conservation.

Consider the function $\beta(E) = -\frac{1}{m} E^m$ which is monotonically decreasing for
all $m \in \mathbb{R}$. The equation (<ref>) for this choice

\begin{equation}
\frac{\partial E}{\partial t} = \nabla \cdot \left ( E^{m-1} \nabla E \right ) \label{very-fast}
\end{equation}

has been studied extensively [67, 68, 69, 70]. It is called the porous medium equation for
$m \in (1, \infty)$, the fast diffusion equation for $m \in (0, 1]$ and the very fast diffusion equation for
$m \in (-1, 0]$. We are mainly interested in very fast diffusion because we had $0 < \alpha \leq 1$ and $m$ is
essentially $\alpha - 1$. We are only writing $m$ instead of $\alpha - 1$ here for consistency with the literature.

Trying to naively derive conservation of energy for (<ref>) does not work. This is
related to the fact that $\frac{\partial E}{\partial t}$ does not need to be bounded at a given time, even if $E$ is
uniformly bounded. To illustrate this, let us work with the very fast diffusion equation in one dimension.
For the initial condition $E_0$, it is perfectly valid to choose a function that decays like $|x|^{-k}$ where $k > 1$.
For large $x$, this leads to a time derivative given by:

\begin{eqnarray}
\frac{\partial E}{\partial t} (x, 0) &=& \frac{\partial}{\partial x} \left [ E^{m - 1} (x, 0) \frac{\partial E}{\partial x} (x, 0) \right ] \nonumber \\
&=& \frac{\partial}{\partial x} \left [ E_0^{m - 1} (x) \frac{\partial E_0}{\partial x} (x) \right ] \nonumber \\
&\sim& \frac{\partial}{\partial x} \left ( x^{(1 - m)k} x^{-k - 1} \right ) \nonumber \\
&\sim& x^{-mk - 2} \; . \nonumber
\end{eqnarray}

The exponent tells us that it is possible for $E_0$ to decay “too quickly”. If $k > \frac{2}{-m}$ then
$\frac{\partial E}{\partial t}$ is not bounded at $t = 0$. This affects our ability to compute

\begin{equation}
\frac{\textup{d}M}{\textup{d}t} = \frac{\textup{d}}{\textup{d}t} \int_{-\infty}^{\infty} E \textup{d}x \; . \nonumber
\end{equation}

In order to take the derivative inside the integral, there must exist an $\epsilon > 0$ such that
$\frac{\partial E}{\partial t}$ is bounded on $(-\infty, \infty) \times [0, \epsilon]$. We cannot do this with all
$k$ values, but even if we could, we would still not be out of the woods.

Pick a value of $k$ satisfying $\frac{1}{-m} < k < \frac{2}{-m}$. Since $k < \frac{2}{-m}$
$\frac{\partial E}{\partial t}$ as derived above must be finite at $t = 0$. Computing the time derivative of the
total mass,

\begin{eqnarray}
\frac{\textup{d}M}{\textup{d}t}(0) &=& \int_{-\infty}^{\infty} \frac{\partial E}{\partial t} (x, 0) \textup{d}x \nonumber \\
&=& \int_{-\infty}^{\infty} \frac{\partial}{\partial x} \left [ E^{m - 1} (x, 0) \frac{\partial E}{\partial x} (x, 0) \right ] \textup{d}x \nonumber \\
&=& \left. E_0^{m - 1} (x) \frac{\partial E_0}{\partial x} (x) \right|_{-\infty}^{\infty} \nonumber \\
&=& 2 \lim_{|x| \rightarrow \infty} E_0^{m - 1} (|x|) \frac{\partial E_0}{\partial x} (|x|) \nonumber \\
&\sim& \lim_{|x| \rightarrow \infty} |x|^{-mk - 1} \; . \nonumber
\end{eqnarray}

Our choice $k > \frac{1}{-m}$ tells us that the time derivative of $M$ is initially infinite. Therefore $M$ cannot
possibly be a constant function on $[0, \infty)$. We have demonstrated that there are infinitely many functions
$E_0 \in L^1(\mathbb{R})$ such that there is no energy conserving soluton to the following:

\begin{equation}
\begin{dcases}
\frac{\partial E}{\partial t}(x, t) = \frac{\partial}{\partial x} \left [ E^{m - 1}(x, t) \frac{\partial}{\partial x} E(x, t) \right ] & (x, t) \in \mathbb{R} \times (0, \infty) \\
E(x, 0) = E_0(x) & x \in \mathbb{R}
\end{dcases} \; . \label{naive-problem}
\end{equation}

As it turns out, however, for any $E_0 \in L^1(\mathbb{R})$, there is [67, 69] an energy
conserving solution to the similar problem:

\begin{equation}
\begin{dcases}
\frac{\partial E}{\partial t}(x, t) = \frac{\partial}{\partial x} \left [ E^{m - 1}(x, t) \frac{\partial}{\partial x} E(x, t) \right ] & (x, t) \in \mathbb{R} \times (0, \infty) \\
\lim_{t \rightarrow 0} E(x, t) = E_0(x) & x \in \mathbb{R}
\end{dcases} \; . \label{modified-problem}
\end{equation}

Differential equations can only hold in open sets so the notation $E(x, 0) = E_0(x)$ really means that any derivative
of $E(\cdot, t)$ will converge to the corresponding derivative of $E_0$ as $t \rightarrow 0$. The subtle distinction
between (<ref>) and (<ref>) is that the only thing that needs to converge is the
function itself. The $E(\cdot, t)$ will converge in $L^1(\mathbb{R})$ to $E_0$ but the
$\frac{\partial}{\partial t} E(\cdot, t)$ need not converge to
$\frac{\partial}{\partial x} \left ( E_0^{m - 1} \frac{\partial E_0}{\partial x} \right )$. In other words,
$\frac{\textup{d}M(0)}{\textup{d}t}$ appeared to diverge because we were still calculating it incorrectly.

Esteban, Rodriguez and Vázquez [67] showed that an energy conserving solution to the very fast diffusion
Cauchy problem exists in one dimension. However, there are also solutions that vanish in finite time such as

\begin{equation}
E(x, t) = (T - t) \left [ 1 - \left ( \frac{1 - e^{-\sqrt{2}x}}{1 + e^{-\sqrt{2}x}} \right )^2 \right ] \; . \nonumber
\end{equation}

This solves (<ref>) for $m = 0$ and does not conserve energy. The theorem of [67] still
applies to initial conditions like

\begin{equation}
E_0(x) = T \left [ 1 - \left ( \frac{1 - e^{-\sqrt{2}x}}{1 + e^{-\sqrt{2}x}} \right )^2 \right ] \; , \nonumber
\end{equation}

it just yields a different solution. In fact, for all $-1 < m \leq 0$ and all initial conditions, there are short lived
solutions to (<ref>) in addition to the one that conserves energy.
One of the following equivalent conditions is needed to make sure we are choosing the right solution.

* $\int_{-\infty}^{\infty} E(x, t) \textup{d}x = \int_{-\infty}^{\infty} E_0(x) \textup{d}x$ for $t \in (0, \infty)$.

* $\lim_{x \rightarrow \pm \infty} E^{m - 1}(x, t) \frac{\partial E}{\partial x}(x, t) = 0$ for $t \in (0, \infty)$.

* $E(x, t) > 0$ for $(x, t) \in \mathbb{R} \times (0, \infty)$.

In fact, [68] generalized this to arbitrary flux functions.

Suppose $E_0 \in L^1(\mathbb{R})$ and $f ,g \in L^{\infty} ((0, \infty))$ are non-negative.
Then there exists a unique solution to

\begin{equation}
\begin{dcases}
\frac{\partial E}{\partial t}(x, t) = \frac{\partial}{\partial x} \left ( E^{m - 1}(x, t) \frac{\partial}{\partial x} E(x, t) \right ) & (x, t) \in \mathbb{R} \times (0, \infty) \\
\lim_{t \rightarrow 0} E(x, t) = E_0(x) & x \in \mathbb{R} \\
\lim_{x \rightarrow \infty} E^{m - 1}(x, t) \frac{\partial}{\partial x} E(x, t) = -f(t) & t \in (0, \infty) \\
\lim_{x \rightarrow -\infty} E^{m - 1}(x, t) \frac{\partial}{\partial x} E(x, t) = g(t) & t \in (0, \infty)
\end{dcases} \nonumber
\end{equation}

whose total mass satisfies

\begin{equation}
\int_{-\infty}^{\infty} E(x, t) \textup{d}x = \int_{-\infty}^{\infty} E_0(x) \textup{d}x - \int_0^t f(s) + g(s) \textup{d}s \; . \nonumber
\end{equation}

The energy conserving solution that interests us is precisely the case $f = g = 0$.
It is important to remember that even though the vanishing flux of $E$ is posed as a constraint, we are still free
to consider an initial profile $E_0$ that does not have a vanishing flux. This is an example of a problem with
inconsistent initial and boundary conditions. Some numerical methods have been written specifically to
address this [71]. Infinite propogation speed ensures that the conditions are only inconsistent on a set of
measure zero. Depending on what we specify for the flux, the distribution might “jump” to a function that decays
quickly or to a function that decays slowly. When $m > 0$, energy conservation in (<ref>) is known to
hold without the complication of non-uniqueness [69].

We may summarize this discussion by saying that the maximum principle only holds for a general
(<ref>) if we specify uniform boundedness and that uniqueness only holds for a general
(<ref>) if we specify energy conservation. The very fast diffusion equation in one dimension
provides an example of uniqueness not holding for (<ref>) if we only specify uniform boundedness.
Even when we have a unique energy conserving solution, the time derivative of the total mass will not make sense at
$t = 0$ for the perfectly good initial conditions that are inconsistent with the boundary conditions. For those, the
energy flux that vanishes for the solution at all positive times will not vanish at $t = 0$.

§.§.§ Barenblatt profiles

Solutions to (<ref>) on a bounded domain converge to their average values which means that
solutions to the same equation on the whole space should converge to zero. This will be an $L^{\infty}$ convergence
and not an $L^1$ convergence because of energy conservation. One should check whether the associated decay times
on $\mathbb{R}^d$ have anything to do with our main time scales (<ref>) and (<ref>). Some
known results about the very fast diffusion equation will help us do this. There is an important family of exact
solutions to the porous medium, fast diffusion and very fast diffusion equations given by:

\begin{equation}
U(x, t) = \left [ \frac{\left ( \frac{4}{1 - m} - 2d \right ) t}{|x|^2 + Bt^{\frac{2}{2 - d(1 - m)}}} \right ]^{\frac{1}{1 - m}} \label{barenblatt}
\end{equation}

where $B$ is a positive constant that determines the mass [70]. These self-similar functions called
Barenblatt solutions have a Dirac delta as their initial conditions. Solving for the mass,

\begin{eqnarray}
M(t) &=& d \omega_d \int_0^{\infty} U(r, t) r^{d - 1} \textup{d}r \nonumber \\
&=& d \omega_d \left [ \left ( \frac{4}{1 - m} - 2d \right ) t \right ]^{\frac{1}{1 - m}} \int_0^{\infty} \frac{d^{d - 1}}{\left ( r^2 + B t^{\frac{2}{2 - d (1 - m)}} \right )^{\frac{1}{1 - m}}} \textup{d}r \nonumber \\
&=& d \omega_d \left [ \frac{\left ( \frac{4}{1 - m} - 2d \right ) t}{Bt^{\frac{2}{2 - d (1 - m)}}} \right ]^{\frac{1}{1 - m}} \int_0^{\infty} \frac{B^{\frac{d}{2}} t^{\frac{d}{2 - d (1 - m)}} s^{d - 1}}{\left ( s^2 + 1 \right )^{\frac{1}{1 - m}}} \textup{d}s \nonumber \\
&=& d \omega_d \left ( \frac{4}{1 - m} - 2d \right )^{\frac{1}{1 - m}} B^{\frac{d}{2} - \frac{1}{1 - m}} \int_0^{\infty} \frac{s^{d - 1}}{\left ( s^2 + 1 \right )^{\frac{1}{1 - m}}} \textup{d}s \nonumber \\
&=& \frac{d \omega_d}{2} \left ( \frac{4}{1 - m} - 2d \right )^{\frac{1}{1 - m}} B^{\frac{d}{2} - \frac{1}{1 - m}} \frac{\Gamma \left ( \frac{d}{2} \right ) \Gamma \left ( \frac{1}{1 - m} - \frac{d}{2} \right )}{\Gamma \left ( \frac{1}{1 - m} \right )} \label{barenblatt-mass}
\end{eqnarray}

where we have recognized a beta function. The time that it takes for a Barenblatt profile's peak to reach
$E_{\mathrm{H}}$ is given by:

\begin{equation}
T = \left [ E_{\mathrm{H}} \left ( \frac{B}{\frac{4}{1 - m} - 2d} \right )^{\frac{1}{1 - m}} \right ]^{(1 - m) - \frac{2}{d}} \; . \label{time-with-b}
\end{equation}

If we wanted to call this a Barenblatt profile of mass $a E_{\mathrm{F}}$, we could use the mass relation
(<ref>) to replace $B$:

\begin{equation}
T = \left ( \frac{4}{1 - m} - 2d \right )^{\frac{2m - d(1 - m)}{d(1 - m)}} E_{\mathrm{H}}^{\frac{d(1 - m) - 2}{d}} \left [ \frac{2 a E_{\mathrm{F}}}{d \omega_d \frac{\Gamma \left ( \frac{d}{2} \right ) \Gamma \left ( \frac{1}{1 - m} - \frac{d}{2} \right )}{\Gamma \left ( \frac{1}{1 - m} \right )}} \right ]^{\frac{\frac{4}{d} - 2(1 - m)}{2 - d(1 - m)}} \; . \label{time-without-b}
\end{equation}

Even though this decay time is concerned with a Dirac delta of mass $a E_{\mathrm{F}}$, it upper bounds the decay time
for a box of mass $a E_{\mathrm{F}}$ by the concentration comparison theorem. Things are even nicer than this; the
Barenblatt profiles attract all solutions to the Cauchy problems for these nonlinear diffusion equations. As time goes
on, general solutions $E(\cdot, t)$ will converge to $U(\cdot, t)$ and the rate of convergence has been found

\begin{equation}
\lim_{t \rightarrow \infty} t^{\frac{d}{2 - d(1 - m)}} \left | \left | E(\cdot, t) - U(\cdot, t) \right | \right |_{L^{\infty}} = 0 \; . \nonumber
\end{equation}

This means that (<ref>) still depends on the mass in the same way as the true decay time. Since
$O \left ( \left ( aE_{\mathrm{F}} \right )^{\frac{\frac{4}{d} - 2 (1 - m)}{2 - d (1 - m)}} \right )$ and
$O \left ( \left ( aE_{\mathrm{F}} \right )^2 \right )$ are very different, we should not have used results from an
unbounded domain to derive (<ref>) and (<ref>). Implicit in all of this is a certain relation
between $d$ and $m$. Consider this list of conditions:

* The constant $\frac{4}{1 - m} - 2d$ must be positive for the Barenblatt solutions to be well defined.

* The overall power of $t$ in (<ref>) must be negative in order for the peak value to decay with time.

* In order for the mass (<ref>) to be finite, $\frac{s^{d - 1}}{\left ( s^2 + 1 \right )^{\frac{1}{1 - m}}}$ must approach zero faster than $\frac{1}{s}$.

* An $O \left ( \left ( aE_{\mathrm{F}} \right )^{\frac{\frac{4}{d} - 2 (1 - m)}{2 - d (1 - m)}} \right )$ decay time should decrease with the mass.

Any one of these four requirements will tell us that $m > \frac{d - 2}{d}$. Since diffusion in our model happens
when $0 < \alpha < 1$, our primary interest is $-1 < m < 0$. This is only above the critical exponent in one dimension
providing another reason why it was necessary to demand a bounded domain. We know from [67, 68] that
when $d = 1$, the very fast diffusion equation yields infinitely many solutions and one energy conserving solution.
The latter is indeed asymptotic to the Barenblatt profile having the same mass. However, in $d = 2$ we already see
that there is no energy conserving solution to the very fast diffusion equation. It would be interesting to find
some analogue of (<ref>) for bounded domains that acts as a fundamental solution and stays well defined for
all $m$ and $d$.

§ NUMERICAL ANALYSIS

The expressions for our time scales rely on the fact that $E_{\mathrm{F}} \gg E_{\mathrm{H}}$ in a domain that has a
large but finite size. To see how sensitively the inequalities depend on these factors, we have undertaken a numerical
test of our time scales and (<ref>) in particular. The code used for most simulations can be found in the
appendix. One thing that is immediately visible in the code is our decision to use initial conditions that decay
according to a power law:

\begin{equation}
E_0(x) = E_{\mathrm{max}} \left ( \frac{1}{1 + x^2} \right )^{\frac{k}{2}} \; . \nonumber
\end{equation}

For the very fast diffusion equation on an unbounded domain, we saw that it was possible for the initial condition to
decay too quickly ($k > \frac{1}{-m}$). Since the nicest functions to simulate are the ones with a finite
$\frac{\textup{d}M}{\textup{d}t}(0)$, this rules out functions $E_0$ that decay exponentially or have compact support.
For low energies, the relevant exponent is $m + 1 = \alpha = \frac{9}{10}$. This confines us to using $k$ values smaller
than $10$ and not so much smaller that the decay takes forever.
The numerics are done on a bounded domain so $\frac{\textup{d}M}{\textup{d}t}(0)$ is technically finite for
any $k$. However, if we were to pick $k > 10$, the flux at the edges would experience unbounded growth as
we made the domain larger casting doubt on our ability to trust those results.

§.§ Implementation details

In (<ref>), we have been considering a function $\beta = \frac{\textup{d}S}{\textup{d}E}$ which
has diffusive behaviour for low and high energies with a Hagedorn phase and possibly a clustering phase in between. The
chosen entropy should be of the Super Yang-Mills type (<ref>) where

\begin{equation}
S(E) \sim
\begin{cases}
E^{\frac{9}{10}} & E < E_{\mathrm{H}} \\
E & E_{\mathrm{H}} < E < E^{\prime}_{\mathrm{F}} \\
E^{\frac{8}{7}} & E^{\prime}_{\mathrm{F}} < E < E_{\mathrm{F}} \\
E^{\frac{3}{4}} & E > E_{\mathrm{F}}
\end{cases}
\; . \nonumber
\end{equation}

A smooth approximation to this can be accomplished with the function:

\begin{equation}
S(E) = \left ( \frac{E}{E_{\mathrm{H}}} \right )^{\frac{9}{10}} \left ( 1 + \frac{E}{E_{\mathrm{H}}} \right )^{\frac{1}{10}} \left ( 1 + \frac{E}{E^{\prime}_{\mathrm{F}}} \right )^{\frac{1}{7}} \left ( 1 + \frac{E}{E_{\mathrm{F}}} \right )^{-\frac{11}{28}} \; . \label{numerical-s1}
\end{equation}

However, the clustering phase where $S$ is convex is difficult to simulate numerically and our results including
(<ref>) have neglected it. Assuming that $E^{\prime}_{\mathrm{F}} = E_{\mathrm{F}}$ so that this phase does
not exist, we arrive at:

\begin{equation}
S(E) = \left ( \frac{E}{E_{\mathrm{H}}} \right )^{\frac{9}{10}} \left ( 1 + \frac{E}{E_{\mathrm{H}}} \right )^{\frac{1}{10}} \left ( 1 + \frac{E}{E_{\mathrm{F}}} \right )^{-\frac{1}{4}} \; . \label{numerical-s2}
\end{equation}

[][$E^{\prime}_{\mathrm{F}} \neq E_{\mathrm{F}}$]
[][$E^{\prime}_{\mathrm{F}} = E_{\mathrm{F}}$]
One the left is the derivative of (<ref>) showing the three phases. The simplest choice for
removing the clustering phase leads us to differentiate (<ref>), the plot on the right.

Plotting the $\beta(E)$ derived from these two functions, we see a slight problem. After the flat region of
(<ref>), $\beta(E)$ increases reaching a local maximum. After this it begins to decrease more and more
quickly until it reaches a point of inflection and asymptotically approaches 0. The plot for (<ref>), on
the other hand, does not look very flat in any region and seems to always be asymptotically approaching 0. There is no
point of inflection where some other behaviour “crosses over” to $E^{-\frac{1}{4}}$ behaviour. We will now show by
brute force that there is no inflection point, ruling out the possibility that this only happens for certain
$E_{\mathrm{H}}$ and $E_{\mathrm{F}}$.

For any $N > 0$, $\beta: \mathbb{R}_{> 0} \rightarrow \mathbb{R}$ defined by

\begin{equation}
\beta(x) = \frac{9}{10} x^{-\frac{1}{10}} (1 + x)^{\frac{1}{10}} \left ( 1 + \frac{x}{N} \right )^{-\frac{1}{4}} + \frac{1}{10} x^{\frac{9}{10}} (1 + x)^{-\frac{9}{10}} \left ( 1 + \frac{x}{N} \right )^{-\frac{1}{4}} - \frac{1}{4N} x^{\frac{9}{10}} (1 + x)^{\frac{1}{10}} \left ( 1 + \frac{x}{N} \right )^{-\frac{5}{4}} \nonumber
\end{equation}

has no point of inflection.

First we will write

\begin{eqnarray}
\beta^{\prime}(x) &=& -\frac{9}{100} \left ( 1 + \frac{x}{N} \right )^{-\frac{1}{4}} \left [ x^{-\frac{11}{10}} (1 + x)^{\frac{1}{10}} - 2x^{-\frac{1}{10}} (1 + x)^{-\frac{9}{10}} + x^{\frac{9}{10}} (1 + x)^{-\frac{19}{10}} \right ] \nonumber \\
&& -\frac{1}{20N} \left ( 1 + \frac{x}{N} \right )^{-\frac{5}{4}} \left [ 9 x^{-\frac{1}{10}} (1 + x)^{\frac{1}{10}} + x^{\frac{9}{10}} (1 + x)^{-\frac{9}{10}} \right ] + \frac{5}{16N^2} x^{\frac{9}{10}} (1 + x)^{\frac{1}{10}} \left ( 1 + \frac{x}{N} \right )^{-\frac{9}{4}} \; . \nonumber
\end{eqnarray}

Now the equation $\beta^{\prime\prime}(x) = 0$ becomes:

\begin{eqnarray}
&& \frac{9}{1000} \left ( 1 + \frac{x}{N} \right )^{-\frac{1}{4}} \left [ 11 x^{-\frac{21}{10}} (1 + x)^{\frac{1}{10}} - 3 x^{-\frac{11}{10}} (1 + x)^{-\frac{9}{10}} - 27 x^{-\frac{1}{10}} (1 + x)^{-\frac{19}{10}} + 19 x^{\frac{9}{10}} (1 + x)^{-\frac{29}{10}} \right ] \nonumber \\
&& +\frac{27}{400N} \left ( 1 + \frac{x}{N} \right )^{-\frac{5}{4}} \left [ x^{-\frac{11}{10}} (1 + x)^{\frac{1}{10}} - 2 x^{-\frac{1}{10}} (1 + x)^{-\frac{9}{10}} + x^{\frac{9}{10}} (1 + x)^{-\frac{19}{10}} \right ] \nonumber \\
&& +\frac{3}{32N^2} \left ( 1 + \frac{x}{N} \right )^{-\frac{9}{4}} \left [ 9 x^{-\frac{1}{10}} (1 + x)^{\frac{1}{10}} + x^{\frac{9}{10}} (1 + x)^{-\frac{9}{10}} \right ] - \frac{45}{64N^3} x^{\frac{9}{10}} (1 + x)^{\frac{1}{10}} \left ( 1 + \frac{x}{N} \right )^{-\frac{13}{4}} = 0 \nonumber \\
&& \frac{3}{125} \left ( 1 + \frac{x}{N} \right )^3 \left [ 11 (1 + x)^3 - 3 x (1 + x)^2 -27 x^2 (1 + x) + 19 x^3 \right ] \nonumber \\
&& +\frac{9}{50N} \left ( 1 + \frac{x}{N} \right )^2 \left [ x (1 + x)^3 - 2 x^2 (1 + x)^2 + x^3 (1 + x) \right ] \nonumber \\
&& +\frac{1}{4N^2} \left ( 1 + \frac{x}{N} \right ) \left [ 9 x^2 (1 + x)^3 + x^3 (1 + x)^2 \right ] - \frac{15}{8N^3} x^3 (1 + x)^3 = 0 \nonumber \\
&& 625 x^6 + (2500N + 1625) x^5 + (7250 N + 2275) x^4 + (9250N + 819) x^3 \nonumber \\
&& + (2340N^2 + 3402N) x^2 + (720N^3 + 972N^2) x + 264N^3 = 0 \; . \nonumber
\end{eqnarray}

Even though we cannot factor this sixth degree polynomial, we may conclude that it has no positive real roots via
Descartes' rule of signs.

What we should take away from this is that removing the clustering phase by hand is too crude to give us a good
filtration function to use in the numerics. Instead we will use a function like the one plotted in Figure
<ref> where quadratic interpolation has been used to guarantee that it is differentiable and decreasing
with an almost-flat section and two inflection points.

For this function, the flat section begins at $E = 1.0$. The region $0.5 < E < 1.0$ is used to interpolate
between the flat section and the $E^{-\frac{1}{10}}$ power law. The flat section ends at $E = 0.9 E_{\mathrm{F}}$ and
the $E^{-\frac{1}{4}}$ power law begins at $E = E_{\mathrm{F}}$. In the region
$0.9 E_{\mathrm{F}} < E < E_{\mathrm{F}}$, a parabola again interpolates. The plot is for $E_{\mathrm{F}} = 10.0$ but
this general rule has been followed for all $E_{\mathrm{F}}$.

§.§.§ The Crank-Nicolson method

Now that the form of (<ref>) has been decided, we need a method capable of finding approximate
solutions to it. Our PDE involves time and space, so it is natural to discretize space and reduce it to the following
set of coupled ODEs:

\begin{eqnarray}
\frac{\textup{d}}{\textup{d}t} E_i(t) &=& f_i(E(t), t) \nonumber \\
f_i(E(t), t) &=& -\frac{\beta(E_{i+1}(t)) - 2\beta(E_i(t)) + \beta(E_{i-1}(t))}{(dx)^2} \; . \nonumber
\end{eqnarray}

The simplest method for integrating this type of ODE system is the forward Euler method:

\begin{equation}
E_i(t + dt) = E_i(t) + f_i(E(t), t) dt \; . \label{forward-euler}
\end{equation}

One might expect a method like (<ref>) to be accurate as long as $dx$ and $dt$ are small. However, one
can use the heat equation to show that things are not so simple.

[][$dt = 0.005$, $dx = 0.1$]
[][$dt = 0.005$, $dx = 0.05$]
A standard demonstration that the forward Euler method is numerically unstable for the heat equation.

Figure <ref> shows that making the discretization smaller is not always an improvement. To avoid
numerical instabilities, the restriction that must be obeyed is $dt < \frac{1}{2} (dx)^2$. This is the precise
restriction for the heat equation but a nonlinear equation would suffer from a similar restriction. This is especially
harsh for investigating the long decay times associated with a Hagedorn phase. Although they are more difficult to
apply, there are other methods known to be stable for large timesteps [72]. One of these is the backward Euler

\begin{equation}
E_i(t + dt) = E_i(t) + f_i(E(t + dt), t + dt) dt \; . \label{backward-euler}
\end{equation}

In (<ref>), we evolve $E$ forward in time by adding the derivative that it has now. In
(<ref>) which looks similar, we add the derivative that it will have after we add it. Since
some inversion is clearly necessary, we refer to this type of method as an implicit method. In this work, we
will use the average of the forward and backward Euler methods, known as the Crank-Nicolson method. In order to
isolate $E_i(t + dt)$ in (<ref>) when treating the heat equation, we must solve a system of linear
equations - one for each site in the lattice. For a nonlinear equation like (<ref>), these become
nonlinear algebraic equations, and hence often require numerical methods themselves. Something that might look familiar

[][Linear scale]
[][Modified logarithmic scale]
Typically we will work with an initial condition like that shown in red. The width of the domain is chosen
so that all of the function's mass can fit below the Hagedorn energy in blue. The plot on the left shows that $E_0(x)$
is very flat for most of the $x$ values. An adaptive $dx$ makes sense and will be chosen so that the minimum $dx$ is
much smaller than the distance between the inflection points. The plot on the right shows the graph of the same function
but more conveniently.

\begin{equation}
x_{n + 1} = x_{n} - \frac{f(x_n)}{f^{\prime}(x_n)} \; . \nonumber
\end{equation}

This is Newton's method where the sequence converges to a zero of $f$ unless the initial guess $x_0$ is sufficiently
far away. The derivative becomes the Jacobian when we have more than one variable.

\begin{eqnarray}
g^{\prime}(E^n(t + dt))_{ij} (E_j^{n + 1}(t + dt) - E_j^n(t + dt)) = -g_i(E^n(t + dt)) \nonumber \\
g_i(E(t + dt)) = E_i(t) + f_i(E(t + dt)) dt - E_i(t + dt) \label{newton-method}
\end{eqnarray}

We should solve this (tridiagonal) linear system for an initial guess of $E_i^0(t + dt) = E_i(t)$ since $dt$ is small
enough to not change the energy much with each step. The Crank-Nicolson method for us is then

\begin{eqnarray}
E_i(t + dt) &=& \frac{1}{2} \left [ E_i(t) + f_i(E(t)) dt + E^{\infty}_i(t + dt) \right ] \nonumber \\
f_i(E(t)) &=& -\frac{\beta(E_{i+1}(t)) - 2\beta(E_i(t)) + \beta(E_{i-1}(t))}{(dx)^2} \; , \label{crank-nicolson}
\end{eqnarray}

with $E^{\infty}_i(t + dt)$ generated in the (<ref>) way. One aspect of this which is still undesirable
is having a constant $dx$. Given the energy profiles we wish to evolve, it will be more convenient to have $dx$ depend
on $i$ so that the lattice can be finely grained in regions where the energy profile changes the most. Figure
<ref> shows the need for a variable spatial step and also establishes the conventions we will use to
have logarithmic axes.

§.§.§ Convergence tests

We wish to test this algorithm, but do not have any exact solutions at our disposal. To verify that our finite
difference scheme converges to a solution as $O(h^2) \equiv O ((dx)^2)$, there is a useful test suggested by
[73] based on the Richardson expansion:

\begin{equation}
E^h(x, t) = E(x, t) + e_2(x, t) h^2 + e_4(x, t) h^4 + \dots \; . \label{richardson}
\end{equation}

Here $E^h$ is the approximate solution computed by the Crank-Nicolson method, which differs from the exact
solution according to functions $e_2$, $e_4$, etc. If we choose different discretizations,
$dx \in \{h, \; 2h, \; 4h\}$ for example, (<ref>) states that for a given $(x, t)$,

\begin{equation}
Q(t) = \lim_{h \rightarrow 0} \frac{E^{4h}(x, t) - E^{2h}(x, t)}{E^{2h}(x, t) - E^h(x, t)} = \lim_{h \rightarrow 0} \frac{16h^2 - 4h^2}{4h^2 - h^2} = 4 \; . \nonumber
\end{equation}

As a better check, we will compare the 2-norms of the numerical solutions instead of choosing a single point $x$.

\begin{equation}
Q(t) = \lim_{h \rightarrow 0} \frac{\left | \left | E^{4h}(\cdot, t) - E^{2h}(\cdot, t) \right | \right |_{\ell^2}}{\left | \left | E^{2h}(\cdot, t) - E^h(\cdot, t) \right | \right |_{\ell^2}} = 4 \label{q-factor}
\end{equation}

Plotting this for five different values of $h$, it looks like $Q(t) = 2$. This would be expected for a non-centred
finite difference method whose Richardson expansion looks like:

\begin{equation}
E^h(x, t) = E(x, t) + e_1(x, t) h + e_2(x, t) h^2 + \dots \; . \nonumber
\end{equation}

The factor $Q(t)$ plotted for simulations that lasted $t = 90.0$. The grid was uniform and five step sizes
were used.

Our departure from the expected error $O(h^2)$ is likely a result of the nonlinear algebraic equations. As soon as we
use Newton's method, Crank-Nicolson is not being followed exactly and an $O(h)$ error can be introduced. Linear
combinations such as $2E^h - E^{2h}$ can be used to cancel this error term in a technique called Richardson
extrapolation. This has been successfully used for other diffusion equations with the Crank-Nicolson method

§.§ Simulation results

Going ahead with the Crank-Nicolson method, one can print out slices of the energy every so often to see how it is
diffusing. The numerics make it clear that there are two time scales of interest. What we have predicted
in (<ref>) is the long decay time needed for a distribution of energy to descend below $E_{\mathrm{H}}$.
Before this, the distribution will reach $E_{\mathrm{F}}$ or some number close to it, at what we call the termalization
time. Figure <ref> shows that for a variety of initial conditions, the energy profiles are qualitatively
similar to step functions when this happens.

[][Wider function]
[][Narrower function]
These initial conditions with the same total mass both evolve toward a state that has a flat line
in the high energy phase near $E_{\mathrm{F}} = 10.0$.

§.§.§ Short time dynamics

[][Same peak]
[][Same mass]
Two plots of the thermalization time for field theory energy $E_{\mathrm{F}} = 10.0$ and different values of
$k$ in the initial condition.

A plot of how the peak energy moves down with time. This corresponds to the $k = 9$ decay in Figure

The left plot in Figure <ref> uses the
$E_0(x) = 200 \left ( \frac{1}{1 + x^2} \right )^{\frac{k}{2}}$ form. It is fairly clear from the figure
that the profile is flattening out somewhere above $\log \left (1 + E_{\mathrm{F}} \right ) = \log 11$. After
interpolation, the effective Hagedorn and field theory energies become slightly shifted with respect to the parameters
$E_{\mathrm{H}}$ and $E_{\mathrm{F}}$. Therefore to measure thermalization times as a function of $k$, we have
used $\frac{3}{2} E_{\mathrm{F}}$ as a rule of thumb for where the peak energy should be. After making a plot for
five values of $k$, we see that the thermalization time is much more sensitive to variations in $L E_{\mathrm{F}} / M$
than $E_{\mathrm{F}} / E_{\mathrm{max}}$.

Another trend we may investigate is how $E(0, t)$ behaves as a function of $t$. Plotting this for $k = 9$ will give us
Figure <ref>. After a slow start at early times, the plot becomes steeper before levelling off again. This
basic shape holds for other $k$ values as well. It is different from the decay of a Barenblatt profile which would
always have $E(0, t)$ as a power law. The time $t^*$ when $E(0, t)$ is changing most quickly appears somewhere in
Figure <ref>. If we work with the functions having similar thermalization times (the ones normalized to have
the same mass), we find that $t^* \approx 0.6 T^*$ in all five cases. The ratio between
$\frac{\textup{d}}{\textup{d}t} E(0, t^*)$ and the average rate of change in the interval $(0, T^*)$ is about $1.8$. We
have not found a way to predict these numbers analytically.

§.§.§ Long time dynamics

This shows one of our initial conditions diffusing all the way through the Hagedorn regime. It takes about
a hundred times longer to do this than it does to reach the thermalization time.

A plot of the decay time for $E_{\mathrm{F}} = 10.0$, $E_{\mathrm{H}} = 1.0$ and five values of $k$. Even
though they specify different shapes for the initial condition, the times are all within $4\%$ of eachother.

A plot of the decay time for $k = 9$, $E_{\mathrm{H}} = 1.0$ and five values of $E_{\mathrm{F}}$. Even
though the profiles begin to flatten out at different heights, their final decay times are all within $0.4\%$ of

The next step is to wait until $E(0, t) < E_{\mathrm{H}}$. An example of what this diffusion looks
like is in Figure <ref>. This time, $T \gg T^*$ is predicted in (<ref>) to depend on the
squared mass. After all, a cylinder of radius $a$ and height $E_{\mathrm{F}}$ has a volume of
$2aE_{\mathrm{F}}$ if the base is one-dimensional. When investigating the short time dynamics, initial conditions of
five different shapes were normalized to have the same mass. Since $T^*$ was found to be almost the same for them, it
is no surprise that $T$ shares this property. Figure <ref> shows a plot of these times.
As long as the mass is kept constant, early and late time dynamics are largely insensitive to $k$ so we will use the
initial condition

\begin{equation}
E_0(x) = 200 \left ( \frac{1}{1 + x^2} \right )^{\frac{9}{2}} \label{main-initial-condition}
\end{equation}

from now on. The major quantity we have not changed yet is the $E_{\mathrm{F}} / E_{\mathrm{H}}$ ratio. We have always
had $E_{\mathrm{F}} = 10.0$ and $E_{\mathrm{H}} = 1.0$ so far. As we know in one dimension, $T$ is proportional to
$M^2$ and not a power of $E_{\mathrm{F}}$. Whether or not we can change $E_{\mathrm{F}}$ and leave $T$ invariant will
be the real test of (<ref>). The condition
$1.0 = E_{\mathrm{H}} \ll E_{\mathrm{F}} \ll E_{\mathrm{max}} = 200.0$ needs to be satisfied so we will choose a few
$10.0 \leq E_{\mathrm{F}} \leq 20.0$ values. The decay times for these choices, plotted in Figure
<ref>, are remarkably close. The exact prefactors in (<ref>) still need to be checked
but it turns out that they are not very constraining. We know that $\alpha = \frac{9}{10}$ and taking the length $L$
from the code, we can plug it into (<ref>) to find $E_{\mathrm{min}} = E_0(L)$. Because
$E_{\mathrm{min}}$ is so small, the prefactor multiplying $M^2$ in the upper bound from (<ref>) is very
different from the prefactor multiplying $M^2$ in the lower bound from (<ref>). In this simulation, they
differ by a factor of

\begin{equation}
E_{\mathrm{min}}^{\alpha - 2} = \left ( 8.69 \cdot 10^{-19} \right )^{\frac{9}{10} - 2} = 7.36 \cdot 10^{19} \; , \nonumber
\end{equation}

which is much larger than any of the numbers put into the simulation by hand.

§.§.§ Higher dimensions

[][Zoomed in]
[][Zoomed out]
On the left is a plot of the decay time in two dimensions for $k = 9$, $E_{\mathrm{H}} = 1.0$ and five values
of $E_{\mathrm{F}}$. These are much smaller than the one-dimensional decay times in Figure <ref>.
The zoomed out version on the right shows that even though these dots do not vary linearly, they are sandwiched between
two bounds which do vary linearly. For the sake of the plot, the blue lines have been understated. In reality the
lower bound is much closer to being horizontal and the upper bound is much closer to being vertical.

Good agreement between our prediction and the numerical results has been demonstrated in one dimension. In our case,
extending the numerics to higher dimensions is simple because we are focusing on spherically symmetric data. Our
equation becomes

\begin{eqnarray}
\frac{\partial E}{\partial t} &=& -r^{1 - d} \frac{\partial}{\partial r} \left ( r^{d - 1} \frac{\partial \beta(E)}{\partial r} \right ) \nonumber \\
&=& -\frac{\partial^2 \beta(E)}{\partial r^2} - \frac{d - 1}{r} \frac{\partial \beta(E)}{\partial r} \; . \nonumber
\end{eqnarray}

This adds a term with a single spatial derivative to our one-dimensional equation from before. A plot like Figure
<ref> should no longer be expected because $T$ in (<ref>) is no longer proportional to
$M^2$. In two dimensions,
$\left ( a E_{\mathrm{F}} \right )^2 = \left ( a^2 E_{\mathrm{F}} \right ) E_{\mathrm{F}} \propto M E_{\mathrm{F}}$.
This means we are looking for a decay time that varies linearly with the field theory energy. After modifying the
Crank-Nicolson code in the appendix, two-dimensional simulations testing this have been performed. The results, plotted
in Figure <ref>, do not appear to be proportional to $E_{\mathrm{F}}$. The $E_{\mathrm{F}} = 20.0$ decay
time is much less than twice the $E_{\mathrm{F}} = 10.0$ decay time. Even though the times scale differently from the
bounds in (<ref>), they do not actually violate the bounds because of the extreme prefactors that differ
by a factor of $10^{19}$.

The typical values for $T$ differ greatly between dimensions one and two. More generally, the decay times when
$d \geq 2$ are much smaller than when $d = 1$. One dimension is special because it is only in this case that the
$\alpha < 1$ Barenblatt profile exists. Figure <ref> shows decay times that are only nonzero because of
the bounded domain - a phenomenon known as instantaneous extinction. This states that $T \rightarrow 0$ as
$E_{\mathrm{min}} \rightarrow 0$, or in other words, the process we are simulating is not well defined in infinite

§ ROOM FOR IMPROVEMENT

The dynamics of our model have shown many similarities to plasma balls but there are three main problems. To start,
our solutions show that energy diffuses to infinity rather slowly once it escapes the central region. The phenomenon
of decay via hadron ejection is not realized. A more serious problem is the instantaneous extinction we have seen for
a lattice with two or more dimensions. Finally, we have seen that some punishing prefactors allow the admissible decay
times to span several orders of magnitude.

Though one can imagine several ad hoc changes to our model that might eliminate these problems, an extension
that refers to conserved quantities other than energy is physically well motivated. The remainder of this thesis
will focus on developing an explicit PDE description for the momentum in a Poincaré invariant theory.

§.§ A common approximation

There is already a widely applicable PDE framework called hydrodynamics that puts energy and momentum on equal footing.
Hydrodynamics is concerned with systems that are close to equilibrium. Because of this, a system may deviate from a
translation invariant state only slightly leading to slowly varrying charges $E(\textbf{x}, t)$,
$\textbf{P}(\textbf{x}, t)$ and $Q(\textbf{x}, t)$. These are associated with long distance variations in their sources
$T(\textbf{x}, t)$, $\textbf{v}(\textbf{x}, t)$ and $\mu(\textbf{x}, t)$ respectively.

A theory with a symmetric stress-energy tensor and a $U(1)$ current has $d + 1 + \frac{(d + 1)(d + 2)}{2}$ independent
components. It is therefore not possible to describe every system using the $d + 2$ hydrodynamic variables above.
Nevertheless, hydrodynamics provides a good description of field theories with high occupation numbers and has
demonstrated a particular aptitude for problems relating to plasma balls and black holes.

§.§.§ Basic hydrodynamics

The hydrodynamics equations are nothing but the local conservation laws

\begin{eqnarray}
\partial_{\mu} T^{\mu \nu} &=& 0 \nonumber \\
\partial_{\mu} J^{\mu} &=& 0 \; . \label{hydro-equations}
\end{eqnarray}

The conserved currents are built from $T$, $\mu$ and $u^{\mu}$ making this a system of $d + 2$ equations for $d + 2$
unknown functions (the velocity satisfies $u^2 = -1$). The rules for writing (<ref>) in terms of the
hydrodynamic variables are called the constitutive relations. There is no limit to how complicated the
constitutive relations might be so it is helpful to use the assumption that the functions are slowly varrying. This
allows one to consider different versions of hydrodynamics based on how many derivatives are kept. In zeroth order
hydrodynamics, the temperature, velocity and chemical potential are not differentiated in the expressions for
$T^{\mu\nu}$ and $J^{\mu}$. In first order hydrodynamics, they are differentiated at most once. These orders are often
referred to as “ideal” and “dissipative” hydrodynamics. It is important to note that this is different from the
linearization of hydrodynamics. Nonlinear equations (e.g. the Navier-Stokes equations) can easily arise from
zeroth or first order hydrodynamics, so it is common to introduce a further order-by-order expansion that drops terms
with too many variables multiplied together.

The first step in developing constitutive relations is decomposing the currents into components that are transverse
and longitudinal with respect to $u^{\mu}$. This leads to

\begin{eqnarray}
T^{\mu\nu} &=& \mathcal{E}u^{\mu}u^{\nu} + \mathcal{P}\Delta^{\mu\nu} + \left ( q^{\mu}u^{\nu} + q^{\nu}u^{\mu} \right ) + t^{\mu\nu} \nonumber \\
J^{\mu} &=& \mathcal{N}u^{\mu} + j^{\mu} \label{hydro-decomposition}
\end{eqnarray}

where we have defined the projector $\Delta^{\mu\nu} = \eta^{\mu\nu} + u^{\mu}u^{\nu}$. Transverse quantities like
$j^{\mu}$, $q^{\mu}$ and $t^{\mu\nu}$ cannot be built out of $T$, $\mu$ and $u^{\mu}$ without derivatives. Therefore
ideal hydrodynamics takes the form:

\begin{eqnarray}
T^{\mu\nu} &=& \varepsilon u^{\mu}u^{\nu} + P \Delta^{\mu\nu} \nonumber \\
J^{\mu} &=& n u^{\mu} \; . \label{ideal-constituents}
\end{eqnarray}

The pressure $P$, while not a charge or a source, is usually given in terms of sources by an equation of state
$P(T, \mu)$. Going to the fluid's rest frame makes it clear that $\varepsilon$ is the energy density and $n$ is the
charge density.

Things become more complicated when we include one order of dissipation [75]. Looking at
(<ref>), many quantities such as $u^{\mu}$, $q^{\mu}$ and $j^{\mu}$ depend on position for a
general non-equilibrium configuration. There is nothing that prevents us from changing $u^{\mu}$ to a different
function of spacetime as long as $q^{\mu}$ and $j^{\mu}$ change as well to keep $T^{\mu\nu}$ and $J^{\mu}$ invariant.
This redundancy, similar to a gauge freedom, is called frame invariance in hydrodynamics [76].
Consider a shift $u^{\prime \mu}(\textbf{x}, t) = u^{\mu}(\textbf{x}, t) + \delta u^{\mu}(\textbf{x}, t)$. To preserve
the normalization, $\delta u^{\mu}$ is transverse to $u^{\mu}$. Using the inverse of (<ref>),
corresponding changes in the coefficients can be calculated to first order.

\begin{eqnarray}
\mathcal{E} &=& u_{\mu} u_{\nu} T^{\mu\nu} \nonumber \\
\delta \mathcal{E} &=& 2 u_{\mu} \delta u_{\nu} T^{\mu\nu} = 2 u_{\mu} \delta u_{\nu} t^{\mu\nu} \approx 0 \nonumber \\
\mathcal{P} &=& \frac{1}{d} \Delta_{\mu\nu} T^{\mu\nu} \nonumber \\
\delta \mathcal{P} &=& \frac{2}{d} u_{\mu} \delta u_{\nu} T^{\mu\nu} \approx 0 \nonumber \\
\mathcal{N} &=& -u_{\mu} J^{\mu} \nonumber \\
\delta \mathcal{N} &=& -\delta u_{\mu} J^{\mu} = -\delta u_{\mu} j^{\mu} \approx 0 \; . \nonumber
\end{eqnarray}

Note that purely dissipative quantities like $j^{\mu}$ and $t^{\mu\nu}$ become second order when multiplied by
$\delta u_{\mu}$. These are all zero which explains why we could write down (<ref>) without worrying
about frame invariance. The same analysis would give $\delta t^{\mu\nu} = 0$ as well. Conversely,

\begin{eqnarray}
q_{\mu} &=& -\Delta_{\mu \alpha} u_{\beta} T^{\alpha\beta} \nonumber \\
\delta q_{\mu} &=& -\delta u^{\beta} T_{\mu\beta} - 2 u_{\mu} \delta u_{\alpha} u_{\beta} T^{\alpha\beta} - \delta u_{\mu} u_{\alpha} u_{\beta} T^{\alpha\beta} \nonumber \\
&=& -\mathcal{P} \delta u_{\mu} - q^{\beta} u^{\mu} \delta u_{\beta} - \delta u_{\beta} t^{\mu\beta} + 2 u_{\mu} \delta u_{\alpha} q^{\alpha} - \mathcal{E} \delta u_{\mu} \nonumber \\
&\approx& -(\mathcal{E} + \mathcal{P}) \delta u_{\mu} \nonumber \\
j_{\mu} &=& \Delta_{\mu\nu} J^{\nu} \nonumber \\
\delta j_{\mu} &=& u_{\mu} \delta u_{\nu} J^{\nu} + u_{\nu} \delta u_{\mu} J^{\nu} \nonumber \\
&=& u_{\mu} \delta u_{\nu} j^{\nu} - \mathcal{N} \delta u_{\mu} \nonumber \\
&\approx& -\mathcal{N} \delta u_{\mu} \; . \nonumber
\end{eqnarray}

We see that a suitable definition of local velocity reduces the number of terms in (<ref>). The one
that makes $j^{\mu} = 0$ is called the Eckart frame while the one that makes $q^{\mu} = 0$ is called the Landau frame.
Naturally, $T^{\prime}(\textbf{x}, t) = T(\textbf{x}, t) + \delta T(\textbf{x}, t)$ and $\mu^{\prime}(\textbf{x}, t) = \mu(\textbf{x}, t) + \delta \mu(\textbf{x}, t)$ are
allowed shifts of the other hydro variables. Under these redefinitions, the coefficients from
(<ref>) become

\begin{eqnarray}
\varepsilon(T^{\prime}, \mu^{\prime}) &=& \varepsilon(T, \mu) + \frac{\partial \varepsilon}{\partial T} \delta T + \frac{\partial \varepsilon}{\partial \mu} \delta \mu \nonumber \\
P(T^{\prime}, \mu^{\prime}) &=& P(T, \mu) + \frac{\partial P}{\partial T} \delta T + \frac{\partial P}{\partial \mu} \delta \mu \nonumber \\
n(T^{\prime}, \mu^{\prime}) &=& n(T, \mu) + \frac{\partial n}{\partial T} \delta T + \frac{\partial n}{\partial \mu} \delta \mu \; . \nonumber
\end{eqnarray}

This means we can define temperature and chemical potential such that $\mathcal{E} = \varepsilon$ and
$\mathcal{N} = n$ [75].

Choosing the Landau frame, constitutive relations are expressions for the scalar $\mathcal{P}$, the transverse vector
$j^{\mu}$ and the transverse traceless symmetric tensor $t^{\mu\nu}$. These may contain any combination of $T$,
$\mu$ and $u^{\mu}$ with one derivative. However, the equations of zeroth order hydrodynamics (<ref>)
give relations between many of these terms up to higher order corrections. The Landau frame expressions that follow
from this are:

\begin{eqnarray}
T^{\mu\nu} &=& \varepsilon u^{\mu} u^{\nu} + \left ( P - \zeta \partial_{\lambda}u^{\lambda} \right ) \Delta^{\mu\nu} - \eta \Delta^{\mu\alpha} \Delta^{\nu\beta} \left ( \partial_{\alpha} u_{\beta} + \partial_{\beta} u_{\alpha} - \frac{2}{d} \eta_{\alpha\beta} \partial_{\lambda} u^{\lambda} \right ) \nonumber \\
J^{\mu} &=& nu^{\mu} - \sigma T \Delta^{\mu\nu} \partial_{\nu} \left ( \frac{\mu}{T} \right ) + \chi_T \Delta^{\mu\nu} \partial_{\nu} T \; . \label{dissipative-constituents}
\end{eqnarray}

The frame invariant parameters $\zeta$, $\eta$ and $\sigma$ are functions of $T$ and $\mu$ that must be determined from
experiment or the microscopic theory. They are called the bulk viscosity, shear viscosity and conductivity respectively.
It turns out that $\chi_T$ must be zero for time-reversal invariance to be satisfied. The linearization of
(<ref>) will be important in what follows.

We linearize around the equilibrium solution which has constant hydro variables and zero velocity. This allows us to
write $u^{\mu} = (1, v^i)$. Explicitly, what we seek are evolution equations for

\begin{eqnarray}
J^0(\textbf{x}, t) &=& n(T(\textbf{x}, t), \mu(\textbf{x}, t)) \nonumber \\
&\approx& n(T_0, \mu_0) + \left ( \frac{\partial n}{\partial T} \right )_0 \delta T(\textbf{x}, t) + \left ( \frac{\partial n}{\partial \mu} \right )_0 \delta \mu(\textbf{x}, t) \nonumber \\
&\equiv& n_0 + \tilde{n}(\textbf{x}, t) \nonumber \\
T^{00}(\textbf{x}, t) &=& \varepsilon(T(\textbf{x}, t), \mu(\textbf{x}, t)) \nonumber \\
&\approx& \varepsilon(T_0, \mu_0) + \left ( \frac{\partial \varepsilon}{\partial T} \right )_0 \delta T(\textbf{x}, t) + \left ( \frac{\partial \varepsilon}{\partial \mu} \right )_0 \delta \mu(\textbf{x}, t) \nonumber \\
&\equiv& \varepsilon_0 + \tilde{\varepsilon}(\textbf{x}, t) \nonumber \\
T^{0i}(\textbf{x}, t) &=& (\varepsilon(T(\textbf{x}, t), \mu(\textbf{x}, t)) + P(T(\textbf{x}, t), \mu(\textbf{x}, t))) v^i(\textbf{x}, t) \nonumber \\
&\approx& (\varepsilon_0 + P_0) v^i(\textbf{x}, t) \nonumber \\
&\equiv& \tilde{P}^i(\textbf{x}, t) \; . \nonumber
\end{eqnarray}

It is important to remember that most coefficients in (<ref>) and (<ref>)
are functions of $T$ and $\mu$. Their derivatives evaluated at $T_0$ and $\mu_0$ should therefore appear in the
linearization. These derivatives, known as susceptibilities, do not have a functional dependence on $T$ and $\mu$
anymore; they are simply numbers. For this reason, we should not be worried if the linearized hydrodynamics equations
appear to have more coefficients than the six we have seen so far. One equation of hydrodynamics is always the
continuity equation $\frac{\partial \varepsilon}{\partial t} = \partial_i T^{0i}$. It is easy to see that this becomes

\begin{equation}
\frac{\partial \tilde{\varepsilon}}{\partial t} = \partial_i \tilde{P}^i \label{continuity-equation}
\end{equation}

because it is already linear. The equation involving stresses linearizes to

\begin{eqnarray}
\frac{\partial \tilde{P}^j}{\partial t} &=& \partial_i \left [ \left ( P - \zeta \partial_k v^k \right ) \delta^{ij} - \eta \left ( \partial^i v^j + \partial^j v^i  - \frac{2}{d} \delta^{ij} \partial_k v^k \right ) \right ] \nonumber \\
&=& \partial_i \left [ \left ( P - \frac{\zeta}{\varepsilon_0 + P_0} \partial_k \tilde{P}^k \right ) \delta^{ij} - \frac{\eta}{\varepsilon_0 + P_0} \left ( \partial^i \tilde{P}^j + \partial^j \tilde{P}^i  - \frac{2}{d} \delta^{ij} \partial_k \tilde{P}^k \right ) \right ] \nonumber \\
&=& \partial_i \left [ \delta^{ij} \left ( \frac{\partial P}{\partial \varepsilon} \right )_0  \tilde{\varepsilon} + \delta^{ij} \left ( \frac{\partial P}{\partial n} \right )_0  \tilde{n} + \delta^{ij} \frac{\frac{2}{d} \eta_0 - \zeta_0}{\varepsilon_0 + P_0} \partial_k \tilde{P}^k - \frac{\eta_0}{\varepsilon_0 + P_0} \left ( \partial^i \tilde{P}^j + \partial^j \tilde{P}^i \right ) \right ] \nonumber \\
&=& \left ( \frac{\partial P}{\partial \varepsilon} \right )_0 \partial^j \tilde{\varepsilon} + \left ( \frac{\partial P}{\partial n} \right )_0 \partial^j \tilde{n} + \frac{\frac{2 - d}{d} \eta_0 - \zeta_0}{\varepsilon_0 + P_0} \partial^j \partial_k \tilde{P}^k - \frac{\eta_0}{\varepsilon_0 + P_0} \partial_k \partial^k \tilde{P}^j \; . \label{stress-equation}
\end{eqnarray}

Finally, to handle the $U(1)$ current

\begin{eqnarray}
\frac{\partial \tilde{n}}{\partial t} &=& \partial_i \left [ nv^i - \sigma T \delta^{ij} \partial_j \left ( \frac{\mu}{T} \right ) \right ] \nonumber \\
&=& \partial_i \left [ \frac{n_0}{\varepsilon_0 + P_0} \tilde{P}^i - \sigma T \delta^{ij} \partial_j \left ( \frac{\mu}{T} \right ) \right ] \nonumber \\
&=& \partial_i \left [ \frac{n_0}{\varepsilon_0 + P_0} \tilde{P}^i - \sigma_0 T \delta^{ij} \left ( \frac{\partial \mu / T}{\partial \varepsilon} \right )_0 \partial_j \tilde{\varepsilon} - \sigma_0 T \delta^{ij} \left ( \frac{\partial \mu / T}{\partial n} \right )_0 \partial_j \tilde{n} \right ] \nonumber \\
&=& \frac{n_0}{\varepsilon_0 + P_0} \partial_i \tilde{P}^i - \sigma_0 T \left ( \frac{\partial \mu / T}{\partial \varepsilon} \right )_0 \partial_i \partial^i \tilde{\varepsilon} - \sigma_0 T \left ( \frac{\partial \mu / T}{\partial n} \right )_0 \partial_i \partial^i \tilde{n} \; . \label{current-equation}
\end{eqnarray}

It is interesting to note that the continuity equation (<ref>) effectively splits into two
equations in non-relativistic hydrodynamics. In this limit, one must insert the speed of light back into
(<ref>) and collect powers of $c$. Taking $c \rightarrow \infty$, inverse powers
vanish leaving only $O(c)$ and $O(1)$ terms in the hydro equations. Since these orders can be considered separate,
they yield conservation of mass and conservation of kinetic energy [77].

Matching the predictions of hydrodynamics with those of the interaction picture leads to the Kubo formulas. A simple
Kubo formula can be derived for the diffusion constant in the heat equation
$\frac{\partial n}{\partial t} - D \frac{\partial^2 n}{\partial x^2} = 0$. Performing a two-sided Fourier transform
in space and a one-sided Fourier transform in time, the equation can be solved as

\begin{equation}
n(k, z) = \frac{n_0(k)}{Dk^2 - iz} = \frac{\chi \mu_0(k)}{Dk^2 - iz} \; . \label{kubo-step1}
\end{equation}

We assume that the initial buildup of charge is due to a local chemical potential that existed before $t = 0$. If this
chemical potential increased adiabatically starting at $t = -\infty$, the corresponding perturbation to the Hamiltonian
is given by

\begin{eqnarray}
H &\mapsto& H - \int_{-\infty}^{\infty} \mu(x, t) n(x, t) \textup{d}x \nonumber \\
&=& H - \int_{-\infty}^{\infty} e^{\epsilon t} \mu_0(x) \theta(-t) n(x, t) \textup{d}x \nonumber \\
&\equiv& H - \delta H \; . \nonumber
\end{eqnarray}

Heisenberg's equation now gives

\begin{eqnarray}
\frac{\textup{d}}{\textup{d}t} \left < n(x, t) \right > &=& -i \left < [ n(x, t), \delta H(x, t) ] \right > \nonumber \\
\left < n(x, t) \right > &=& -i \int_{-\infty}^t \int_{-\infty}^{\infty} e^{\epsilon t^{\prime}} \mu_0(x) \theta(-t^{\prime}) \left < [n(x, t), n(x^{\prime}, t^{\prime})] \right > dx^{\prime} dt^{\prime} \nonumber \\
&=& -i \int_{-\infty}^0 \int_{-\infty}^{\infty} e^{\epsilon t^{\prime}} \mu_0(x) \theta(t-t^{\prime}) \left < [n(x, t), n(x^{\prime}, t^{\prime})] \right > dx^{\prime} dt^{\prime} \nonumber \\
&=& -\int_{-\infty}^o \int_{-\infty}^{\infty} e^{\epsilon t^{\prime}} \mu_0(x) G_{nn}^{\mathrm{R}} (t - t^{\prime}, x - x^{\prime}) \textup{d}x^{\prime} \textup{d}t^{\prime} \; . \nonumber
\end{eqnarray}

Here $G_{nn}^{\mathrm{R}} (t - t^{\prime}, x - x^{\prime})$ is the retarded Green's function of $n$ with itself.
The one fact we need is that the zero mode of $G_{nn}^{\mathrm{R}}$ is the susceptibility $-\chi$. It is a matter of
complex analysis to show that

\begin{eqnarray}
\left < n(k, z) \right > &=& -\mu_0(k) \int_{-\infty}^{\infty} G_{nn}^{\mathrm{R}} (\omega, k) \frac{1}{(\epsilon + i\omega)(\epsilon + i (\omega - z))} \frac{\textup{d}\omega}{2\pi} \nonumber \\
&=& \frac{\mu_0(k)}{iz} \left [ G_{nn}^{\mathrm{R}} (0, k) - G_{nn}^{\mathrm{R}} (z, k) \right ] \nonumber \\
&=& \frac{i \mu_0(k)}{z} \left [ \chi + G_{nn}^{\mathrm{R}} (z, k) \right ] \; . \label{kubo-step2}
\end{eqnarray}

Comparing (<ref>) to (<ref>), we see that a quantum theory well described by linear diffusion
should be one whose momentum space retarded Green's function is $\frac{\chi D k^2}{iz - Dk^2}$. A calculation that is
longer but equally straightforward has the linear hydrodynamics equations in place of the heat equation.
This is done in [75]. Additionally a field theory is developed that allows one to compute corrections to the
Green's functions arising from small nonlinearities in the hydrodynamics equations.

§.§.§ The fluid / gravity correspondence

Correlation functions in certain CFTs are related to supergravity amplitudes in certain AdS backgrounds.
The link between Green's functions and linearized hydrodynamics opens up the possibility of describing fluid
phenomena using gravity. This was demonstrated in [78] which found that the shear viscosity in Super
Yang-Mills theory is given by $\eta = \frac{\pi}{8} N^2 T^3$. By now it is known that in arbitrary dimension, the
holographic shear viscosity differs from the entropy density by a factor of $\frac{1}{4\pi}$ [79]. A
major result of 2007 is that all transport coefficients in the fluid stress-energy tensor can be algorithmically found
from the gravity side as well [80]. This expansion, which is valid to all orders, is called the
fluid / gravity correspondence.

In the examples discussed previously, we related the energy of a state on the boundary of AdS to the mass of a black
hole in the bulk. More generally, there is a systematic way to find the boundary stress-energy tensor corresponding to
a given bulk metric [81, 82]. Because every metric yields a conserved stress-energy tensor, it is not
surprising that the fluid / gravity correspondence exists. However, the limit in which temperature and velocity suffice
to describe $T^{\mu\nu}$ could have corresponded to an intractable $g_{\mu\nu}$. Showing otherwise, [80]
started with the metric

\begin{equation}
ds^2 = L^2 \left [ -\rho^2 \left ( 1 - \frac{1}{z_0^d \rho^d} \right ) dv^2 + 2 dv d\rho + \rho^2 dx_i dx^i \right ] \; . \label{ed-fink-poincare}
\end{equation}

This is nothing but the black brane (<ref>) written in ingoing Eddington-Finkelstein
co-ordinates with $\rho = \frac{1}{z}$. Boosting the brane to a particular velocity, this becomes

\begin{equation}
ds^2 = L^2 \left [ -\rho^2 \left ( 1 - \frac{1}{z_0^d \rho^d} \right ) u_{\mu} u_{\nu} dx^{\mu} dx^{\nu} - 2 u_{\mu} dx^{\mu} d\rho + \rho^2 \Delta_{\mu\nu} dx^{\mu} dx^{\nu} \right ] \; . \label{boosted-ed-fink-poincare}
\end{equation}

Greek indices have been used for all co-ordinates except $\rho$ which is not in the field theory. The key step is
to promote $z_0$ (which determines the temperature) and $u^{\mu}$ to slowly varrying functions of spacetime. One
could repeat the calculation for a charged black hole if she wanted chemical potential to vary as well. Clearly, these
functions are heavily constrained for (<ref>) to still solve Einstein's equations. The
constraints turn out to be those of hydrodynamics with infinitely many orders of dissipation. Put another way, we may
say that temperature and velocity functions satisfying SYM hydrodynamics at a given order, cause
(<ref>) to only violate Einstein's equations at a higher order.

Instead of making $z_0$ and $u_{\mu}$ functions of $x$, [80] makes them functions of $\epsilon x$ to keep
track of derivatives before setting $\epsilon = 1$. Using $g^{(0)}$ to denote the metric of
(<ref>), $g^{(0)}$ with unconstrained $z_0$ and $u_{\mu}$ has a stress-energy tensor like
(<ref>) and violates the Einstein equations
at order $\epsilon$. To correct this and make the violation order $\epsilon^2$, one must do two things. The first
is to write $z_0 = z_0^{(0)} + \epsilon z_0^{(1)}$, $u_{\mu} = u_{\mu}^{(0)} + \epsilon u_{\mu}^{(1)}$ so that $g^{(0)}$
picks up terms of order $\epsilon$. When these corrections are explicitly calculated, $z_0^{(1)}$ e.g. will
be some multiple of a derivative of $z_0^{(0)}$. The second is to add a new piece of the metric
$g = g^{(0)} + \epsilon g^{(1)}$. If we plug this new metric into Einstein's equation, a vanishing $\epsilon$ term
will dictate that $z_0^{(0)}$ and $u_{\mu}^{(0)}$ satisfy ideal hydrodynamics (<ref>). It will also
allow us to compute the corrected metric and show that it corresponds to the stress-energy tensor
(<ref>). The result of this is [56]

\begin{eqnarray}
ds^2 &=& L^2 \left [ -\rho^2 \left ( 1 - \frac{1}{z_0^{(0) d} \rho^d} \right ) u^{(0)}_{\mu} u^{(0)}_{\nu} dx^{\mu} dx^{\nu} - 2 u^{(0)}_{\mu} dx^{\mu} d\rho + \rho^2 \Delta^{(0)}_{\mu\nu} dx^{\mu} dx^{\nu} \right. \nonumber \\
&& + 2 \rho^2 z_0^{(0)} F \left ( z_0^{(0)} \rho \right ) \left ( \Delta^{(0)}_{\mu \alpha} \Delta^{(0)}_{\mu \beta} \left ( \partial^{\alpha} u^{(0) \beta} + \partial^{\beta} u^{(0) \alpha} \right ) - \frac{1}{d - 1} \Delta^{(0)}_{\mu \nu} \partial_{\lambda} u^{(0) \lambda} \right ) dx^{\mu} dx^{\nu} \nonumber \\
&& \left. + \frac{2}{d - 1} \rho u^{(0)}_{\mu} u^{(0)}_{\nu} \partial_{\lambda} u^{(0) \lambda} dx^{\mu} dx^{\nu} - \rho u^{(0) \lambda} \partial_{\lambda} \left ( u^{(0)}_{\mu} u^{(0)}_{\nu} \right ) dx^{\mu} dx^{\nu} \right ] \label{fgc-metric}
\end{eqnarray}

\begin{equation}
F(x) = \int_x^{\infty} \frac{y^{d - 1} - 1}{y(y^d - 1)} \textup{d}y \; . \nonumber
\end{equation}

Of course (<ref>) is still in the form of (<ref>). The choice to not collect
all the $dx^{\mu} dx^{\nu}$ terms makes it easy to see the non-derivative part whose components satisfy ideal
hydrodynamics. Continuing the procedure with $z_0^{(2)}$ and $u_{\mu}^{(2)}$ and $g^{(2)}$,
these same components are constrained to satisfy dissipative hydrodynamics. Rather
than (<ref>) associated with the dissipative stress-energy tensor (<ref>),
Einstein's equations then yield a longer metric associated with the stress-energy tensor for a conformal fluid at two
orders of dissipation. This logic continues inductively showing many non-trivial constitutive relations.
Beyond just conformal fluids, similar constitutive relations may be derived for fluids in spaces with compactified
directions. These relations have been used to study plasma balls and other objects that appear in the Witten
model [83, 84, 85]. These studies assume that a black hole is already present because the process of
black hole formation is outside the regime of hydrodynamics [56].

Simply replacing the stochastic model of this thesis with hydrodynamics is not what we plan to do. For one thing, it is
hoped that the initial conditions of Figure <ref> correspond to black holes that have not formed yet.
For hydrodynamics to apply, we must wait for the inhomogeneities to be smoothed out. Therefore, we should only expect
agreement between hydrodynamics and a suitable extension of our model if we linearize both of them. Another reason to
extend our model comes from the importance of the density of states. When we only had energy being transferred, this
function allowed us to read off diffusing or clustering behaviour. It is plausible that something similar can be done
when we include momentum.

§.§ Restricting the density of states

If momentum is viewed as a quantity exchanged between nearest neighbour sites, the rate for each transition naturally
depends on a momentum restricted density of states $\rho(E, \textbf{P})$. While expressions for
$\rho(E, \textbf{P})$ appear in some models of electron structure [86], the restricted density of states for
a field theory is a barely studied quantity. Here, we attempt to rectify this by deriving some properties of the
restricted density of states for simple field theories. A strongly coupled SYM expression analogous to
(<ref>) is most likely beyond our reach.

§.§.§ Some conformal field theory

Some of the most interesting techniques for counting states are applicable to two-dimensional CFTs. Consider a
Minkowskian theory on $\mathbb{S}^1 \times \mathbb{R}$ where the circle has radius 1 by convention. At thermal
equilibrium, this can be thought of as a Euclidean theory on $\mathbb{S}^1 \times \mathbb{S}^1$ since temperature and
time are related by a Wick rotation. The radius of the second circle must be $\beta$. If we rescale each $\mathbb{S}^1$
by the same factor, one such rescaling gives us another Euclidean theory on a circle of radius 1 as shown in
Figure <ref>.

These two theories live on a circle of volume $V$ but their temperatures are different. They are related
through multiplication by $\frac{V}{\beta}$ which is a conformal transformation.

Conformal invariance then tells us that inverse temperature $\beta$ and inverse temperature $\frac{V^2}{\beta}$
are completely equivalent [87]. In terms of the modular parameter we defined for (<ref>),
$\tau = i\frac{\beta}{2\pi}$ becomes $i\frac{4\pi^2 / \beta}{2\pi} = -\frac{1}{\tau}$. For a $\tau$ that is not
imaginary, the factor we trace over to compute the generalized partition function is
$e^{-2\pi \tau_2 H + 2\pi i \tau_1 P}$. Clearly, this does not change if we increase the real part of $\tau$ by 1.
Symmetries of the generalized partition function can therefore be written as

\begin{eqnarray}
\tau \mapsto \tau + 1 \nonumber \\
\tau \mapsto -\frac{1}{\tau} \; . \label{modular-transformations}
\end{eqnarray}

These two transformations generate the modular group which is the discrete but infinite group
$SL(2, \mathbb{Z})$. The utility of modular invariance for studying thermodynamics was first noticed by Cardy in 1986
[88]. Our previous expressions for the partition function of a free theory do not posess modular
invariance. This is because they are only asymptotic expressions for the high temperature limit of $Z$. Low temperature
information is lost when we convert sums to integrals. It is perhaps for this reason that most sources prefer to write
infinite products or special functions when simple expressions like (<ref>) and (<ref>) are all
that are needed.

For an infinite collection of oscillators with frequencies
$1, 2, 3, \dots$, the energy of an arbitrary state $\left | N_1, N_2, N_3, \dots \right >$ is
$N_1 + 2N_2 + 3N_3 + \dots$. The degeneracy of energy $n$ is simply the number of ways to make $n$ by picking a certain
number of ones, a certain number of twos, a certain number of threes and so on. This well-studied object $p(n)$ is
usually called “the partition function” in number theory. We will call this microcanonical quantity
“the partition sum” to distinguish it from the canonical partition function we have been using so far. The first
asymptotic $p(n) \sim \frac{1}{4 \sqrt{3} n} e^{\pi \sqrt{\frac{2}{3} n}}$ due to Hardy and Ramanujan was based on
transformation properties of the generating function

\begin{eqnarray}
F(q) &=& \sum_{n = 1}^{\infty} p(n) q^n \nonumber \\
&=& \left [ 1 + q + q^2 + \dots \right ] \left [ 1 + q^2 + q^4 + \dots \right ] \dots \nonumber \\
&=& \prod_{n = 1}^{\infty} \left ( 1 - q^n \right )^{-1} \; . \label{genfunc}
\end{eqnarray}

If we take $q = e^{2\pi i \tau}$ (and $\tau = i\frac{\beta}{2\pi}$ as before), (<ref>) becomes a sum of
Boltzmann factors. It is not quite the partition function of our oscillator family because we have not yet included
the Casimir effect. Accounting for zero point energies, each term in the product (<ref>) picks up a factor of
$q^{\frac{n}{2}}$ giving us

\begin{eqnarray}
Z_1(\tau) &=& \prod_{n = 1}^{\infty} q^{\frac{n}{2}} \left ( 1 - q^n \right )^{-1} \nonumber \\
&=& q^{-\frac{1}{24}} \prod_{n = 1}^{\infty} \left ( 1 - q^n \right )^{-1} \nonumber \\
&=& \eta^{-1}(\tau) \; . \label{partfunc-half}
\end{eqnarray}

We have analytically continued and written this partition function in terms of the Dedekind eta function

\begin{equation}
\eta(\tau) = e^{\frac{\pi i \tau}{12}} \prod_{n = 1}^{\infty} \left ( 1 - e^{2\pi i n \tau} \right ) \; . \label{eta-function}
\end{equation}

Referring to the eta function, Rademacher gave a convergent series for $p(n)$

\begin{eqnarray}
p(n) &=& \frac{1}{\sqrt{2} \pi} \sum_{k = 1}^{\infty} A_k(n) \sqrt{k} \frac{\textup{d}}{\textup{d}n} \left ( \frac{\sinh \left ( \frac{\pi}{k} \sqrt{\frac{2}{3} \left ( n - \frac{1}{24} \right )} \right )}{\sqrt{n - \frac{1}{24}}} \right ) \label{rademacher} \\
A_k(n) &=& \sum_{m \in \left ( \mathbb{Z} / k\mathbb{Z} \right )^*} e^{\pi i \left ( s(m, k) - \frac{2nm}{k} \right )} \nonumber
\end{eqnarray}

where the coefficients are written in terms of the Dedekind sum

\begin{equation}
s(m, k) = \sum_{l = 1}^{k - 1} \frac{l}{k} \left ( \frac{ml}{k} - \bigg \lfloor \frac{ml}{k} \bigg \rfloor - \frac{1}{2} \right ) \; . \nonumber
\end{equation}

Just as we treated $\rho(E)$ as the inverse Laplace transform of $Z(\beta)$, the partition sum is an integral transform
of its generating function.

\begin{equation}
p(n) = \frac{1}{2\pi i} \int_{\gamma} \frac{F(q)}{q^{n + 1}} \textup{d}q \nonumber
\end{equation}

Instead of a vertical line, the contour $\gamma$ that leads to (<ref>) is a beautiful shape built up
from fractals in the complex plane. A very readable thesis on this derivation is [89].

The main fact about $F(q)$ used in the proof is that it is an imaginary exponential times a function which plays
nicely with the modular group — as it would be even if our theory were not free. Under the modular
generators (<ref>), the Dedekind eta function transforms in the following way:

\begin{eqnarray}
\eta(\tau + 1) &=& e^{\frac{\pi i}{12}} \eta(\tau) \nonumber \\
\eta \left ( -\frac{1}{\tau} \right ) &=& \sqrt{-i \tau} \eta(\tau) \; . \label{modular-eta}
\end{eqnarray}

The first of these equalities is trivial. The second is not but there is a clever proof containing no more
than a page of algebra [90]. A function with transformation properties similar to (<ref>) is
called a modular form. More precisely, a holomorphic function $f$ in the upper half plane is a modular form of
weight $w$ if

\begin{equation}
\left | f \left ( -\frac{1}{\tau} \right ) \right | = \tau^w \left | f(\tau) \right | \; . \nonumber
\end{equation}

According to this definition, $\eta$ is a modular form of weight $\frac{1}{2}$. The Rademacher formula can be
rederived for (multiples of) other modular forms. The result is [91, 92]

\begin{eqnarray}
p_w(n) &=& 2\pi \sum_{m - \frac{c}{24} < 0} \left ( \frac{n - \frac{c}{24}}{\left | m - \frac{c}{24} \right |} \right )^{\frac{w - 1}{2}} p_w(m) \nonumber \\
&& \sum_{k = 1}^{\infty} \frac{1}{k} Kl \left ( n - \frac{c}{24}, m - \frac{c}{24}; k \right ) I_{1 - w} \left ( \frac{4\pi}{k} \sqrt{\left | m - \frac{c}{24} \right | \left ( n - \frac{c}{24} \right )} \right ) \label{rademacher-w}
\end{eqnarray}

where the Kloosterman sum is defined by

\begin{equation}
Kl(n, m; k) = \sum_{d \in \left ( \mathbb{Z} / k\mathbb{Z} \right )^*} e^{\frac{2\pi i}{k} \left ( dn + d^{-1} m \right )} \; . \nonumber
\end{equation}

For our purposes, we will choose a weight $w$ and use $k = m = 0$. Uniqueness of the vacuum gives $p_w(0) = 1$ and the
central charge is 2 when we have left movers and right movers. This formula can be understood as the origin of the
Bessel function in (<ref>).

In order to use the partition sum and its asymptotics to study $\rho(E, P)$, there are some changes that must be made
to our theory of positively indexed oscillators.

The theory we have discussed so far resembles a harmonic potential where arbitrarily many non-interacting
bosons can be piled into each level.

To start, we should introduce negatively indexed oscillators to make it more like a field theory in momentum space.
An excitation of any oscillator increases the energy. Whether it increases or decreases momentum depends on whether
the frequency is positive or negative. This means that in the contribution of the newly introduced oscillators,
$\tau_2$ has the same sign as in $Z_1(\tau)$ and $\tau_1$ has the opposite. Including a sum of
$e^{-2\pi i n \tau^*}$ Boltzmann factors along with our sum of $e^{2\pi i n \tau}$ Boltzmann factors suggests that
the partition function for the doubly infinite family of oscillators is $\left | \eta(\tau) \right |^{-2}$.
This is almost correct, but we need to remember that the Casimir contribution doubles instead of vanishes when we add
more oscillators. Fixing this, the appropriate partition function is given by

\begin{equation}
Z_2(\tau) = \eta^{-1}(\tau) \eta^{-1}(-\tau^*)^* \; . \label{partfunc-full}
\end{equation}

This is still not modular invariant because we have not yet included the oscillator with zero frequency. The zero mode
has a continuous spectrum which is quadratic in the free variable $\xi$ [87]. This can be seen e.g. in
the scalar field mode expansion that was discussed in the context of the string action (<ref>).
The most exact expression for the partition function of scalar field theory is therefore

\begin{eqnarray}
Z_3(\tau) &=& Z_2(\tau) \int_{-\infty}^{\infty} e^{-\pi \Im \tau \xi^2} \frac{\textup{d}\xi}{2\pi} \nonumber \\
&=& \frac{1}{2\pi \sqrt{\Im \tau}} \eta^{-1}(\tau) \eta^{-1}(-\tau^*)^* \; . \label{partfunc-final}
\end{eqnarray}

To summarize these results, the partition function $Z_1$ is a modular form of weight $-\frac{1}{2}$. By
(<ref>), its density of states is asymptotic to the partition sum;
$p_{-\frac{1}{2}}(n) \sim \frac{1}{4 \sqrt{3} n} e^{\pi \sqrt{\frac{2}{3} n}}$. Doubling the number of oscillators,
we made the partition function $Z_2$ which is a modular form of weight $-1$. Its density of states has the
behaviour $p_{-1}(n) \sim \frac{3^{\frac{1}{4}}}{12 n^{\frac{5}{4}}} e^{2\pi \sqrt{\frac{1}{3} n}}$. The partition
function $Z_3$, since it describes a CFT, is a modular form of weight 0. If we take the density of states
$p_0(n) \sim \frac{3^{\frac{1}{4}}}{\sqrt{12} n^{\frac{3}{4}}} e^{2\pi \sqrt{\frac{1}{3} n}}$ and substitute
$n = \frac{VE}{2\pi}$ for the conformal weight, our expression agrees with (<ref>) which gave $\rho(E)$ for a
free CFT from an entirely different setup.

§.§.§ Convoluted functions

Evidently, the Casimir term and the zero mode do not affect $\log p_w (n)$ to leading order. It is therefore valid to
use $F(q)^2$ instead of $Z_3(\tau)$ as the generating function. With this approximation, the convolution formula gives
a very clear picture of how the density of states is comprised.

\begin{eqnarray}
p_0(n) &\approx& \frac{1}{n!} \frac{\textup{d}^n}{\textup{d}q^n} F(q)^2 \nonumber \\
&=& \sum_{k = 0}^n p(k) p(n - k) \nonumber
\end{eqnarray}

The first term $p(0) p(n)$ has $n$ right movers and $0$ left movers and corresponds to energy $n$ and momentum $n$.
The last term $p(n) p(0)$ has $0$ right movers and $n$ left movers and corresponds to energy $n$ and momentum $-n$.
The middle term $p \left ( \frac{n}{2} \right )^2$ has a left mover for every right mover and corresponds to energy $n$
and momentum $0$. This tells us that in a $1 + 1$-dimensional possibly interacting CFT,

\begin{equation}
\rho(E, P) = \frac{1}{2} p \left ( \frac{V}{4\pi} (E - P) \right ) p \left ( \frac{V}{4\pi} (E + P) \right ) \label{rdos-2d}
\end{equation}

is a good approximation to the momentum restricted density of states. This only becomes difficult to evaluate when
$|P| \neq E$ are of the same order.

Essentially we derived (<ref>) by starting with a function whose asymptotic behaviour is
$\log p(n) \sim \pi \sqrt{\frac{2}{3} n}$. The convolution of two of them yields a function that follows
$\log p_0(n) \sim \pi \sqrt{\frac{4}{3} n}$. In other words, convolution doubles the central charge. At least for free
theories, this phenomenon can be readily seen in higher dimensions. The crudest form of (<ref>) is

\begin{equation}
\log \rho(E) \sim \left ( \frac{(d + 1)^{d + 1}}{d^d} AVE^d \right )^{\frac{1}{d + 1}} \; . \nonumber
\end{equation}

It will be helpful to write this density of states explicitly as $\rho(E; A)$ from now on. For a bosonic theory, we had
$A$ proportional to $s$ and for a fermionic theory, we had $A$ proportional to $s^*$. Writing the partition
function $\prod_{\textbf{p}} Z(\textbf{p})^s Z^*(\textbf{p})^{s^*}$, it becomes clear that $\rho(\cdot; A)$ is
$\rho \left ( \cdot; \frac{A}{2} \right )$ convolved with itself. Performing substitutions in the
resulting integral,

\begin{eqnarray}
\rho(E; A) &=& \int_0^E \rho \left (\xi; \frac{A}{2} \right ) \rho \left (E - \xi; \frac{A}{2} \right ) \textup{d}\xi \nonumber \\
&=& \frac{1}{2} \int_{-E}^E \rho \left ( \frac{E - P}{2}; \frac{A}{2} \right ) \rho \left ( \frac{E + P}{2}; \frac{A}{2} \right ) \textup{d}P \nonumber \\
&=& \int_0^E \rho \left ( \frac{E - P}{2}; \frac{A}{2} \right ) \rho \left ( \frac{E + P}{2}; \frac{A}{2} \right ) \textup{d}P \nonumber \\
&=& \int_0^\infty \rho \left ( \frac{E - P}{2}; \frac{A}{2} \right ) \rho \left ( \frac{E + P}{2}; \frac{A}{2} \right ) \textup{d}P \nonumber \\
&=& \frac{1}{d \omega_d} \int_{\mathbb{R}^d} \rho \left ( \frac{E - | \textbf{P} |}{2}; \frac{A}{2} \right ) \rho \left ( \frac{E + | \textbf{P} |}{2}; \frac{A}{2} \right ) \textup{d}\textbf{P} \; . \label{posterior-integral}
\end{eqnarray}

In the second last step we have used the fact that densities of states are zero for negative arguments. Part of the
definition of the restricted density of states is

\begin{equation}
\rho(E; A) = \int_{\mathbb{R}^d} \rho(E, \textbf{P}; A) \textup{d}\textbf{P} \; . \label{prior-integral}
\end{equation}

Comparing (<ref>) to (<ref>), a plausible formula for the restricted density of
states in a free CFT is

\begin{equation}
\log \rho(E, \textbf{P}) \sim \left ( \frac{(d + 1)^{d + 1}}{2(2d)^d} AV(E - | \textbf{P} |)^d \right )^{\frac{1}{d + 1}} + \left ( \frac{(d + 1)^{d + 1}}{2(2d)^d} AV(E + | \textbf{P} |)^d \right )^{\frac{1}{d + 1}} \; . \label{rdos-free}
\end{equation}

In order to use the convolution formula, we factored the partition function in a way that is not always valid. This
means we again have results for a free CFT in arbitrary dimension and an arbitrary CFT in two dimensions. It is
possible that further results could be obtained using an entropy proposed by Erik Verlinde in 2000
[93, 94]. This Cardy-Verlinde formula suggests that higher dimensional CFTs are thermodynamically more
similar to their two-dimensional cousins than previously thought.

§ ENTROPIC DYNAMICS OF MOMENTUM

Retracing our steps, there are nonlinear PDEs associated with changes in the momentum distribution on a lattice.
Each lattice direction, which we will call a large direction, needs to have one evolution equation associated with it.
Additionally, each site is allowed to have compact or small directions which do not show up in the PDEs.
Therefore the $\rho(E, \textbf{P})$ we should use comes from a more restricted function
$\rho(E, \textbf{P}, \textbf{P}^*)$ for the whole field theory and has $\textbf{P}^*$ for the small directions
integrated out. Most of the simulations for our previous model had one large direction. Because we used the density of
states (<ref>) for SYM on $\mathbb{S}^3$, it is tempting to say that these simulations had three small
directions. This is misleading as Figure <ref> shows.

[][Sharing no directions]

[][Sharing one direction]
If one uses a density of states for two spatial dimensions, both of these situations are compatible with the
lattice being a line. One has an overall dimensionality of $2 + 1$, the other $3 + 1$. The task of introducing
momentum to our model forces us to address this ambiguity.

If each site is viewed as a miniature field theory whose spatial directions are all orthogonal to the lattice, there is
nothing left after integrating out $\textbf{P}^*$. The result is simply an unrestricted density of states
$\rho(E)$ precluding our attempts to extend the model. What we really need are sites that include the directions of
the lattice in addition to their compact directions. Therefore what we considered previously was not $\mathbb{S}^3$
with radius $r$ but $\mathbb{S}^2 \times [0, r]$. This is a background with only two small directions because the
non-periodic $[0, r]$ matches up with the lattice. By arguments in [7], the density of states for such a
theory should be similar to (<ref>). Including time, this theory with one large direction and two small
ones lives on $\mathbb{R}^2 \times \mathbb{S}^2$. The Witten model, which is known to host plasma balls, has the
background $\mathbb{R}^4 \times \mathbb{S}^1$ suggesting that its lattice sites look like
$\mathbb{S}^1 \times [0, r]^3$. The equations that follow will have $d$ referring to the number of large spatial
dimensions only.

§.§ Proceeding by analogy

In order to derive our PDE for the energy

\begin{equation}
\frac{\partial E}{\partial t} = \partial_i \left ( C(E) \rho^2(E) \partial_i \frac{\textup{d} \log \rho(E)}{\textup{d}E} \right ) \; , \label{original-equation}
\end{equation}

we took the continuum limit of the equation

\begin{equation}
\frac{\partial n_c}{\partial t} = \sum_{\left < b, c \right >} \sum_{k \neq 0} k W_{(n_c, n_b) \rightarrow (n_c + k, n_b - k)} \; . \label{first-non-continuum}
\end{equation}

The main steps involved were assuming that $k$ was some small amount of energy $\pm \epsilon$ and writing $n_c$ as
$E(\textbf{x})$ where site $b$ was “one lattice vector” away from site $c$. Then (<ref>) became

\begin{equation}
\frac{\partial E(\textbf{x})}{\partial t} = \epsilon \sum_{\textbf{e} \in \{ \pm \textbf{e}_1, \dots, \pm \textbf{e}_d \}} \left [ W_{(E(\textbf{x}), E(\textbf{x} + a\textbf{e})) \rightarrow (E(\textbf{x}) + \epsilon, E(\textbf{x} + a\textbf{e}) - \epsilon)} - W_{(E(\textbf{x}), E(\textbf{x} + a\textbf{e})) \rightarrow (E(\textbf{x}) - \epsilon, E(\textbf{x} + a\textbf{e}) + \epsilon)}\right ] \label{first-continuum}
\end{equation}

and we took $\epsilon$ and $a$ to zero. The sum runs over positive and negative versions of the standard basis
vectors so that in a cubic lattice, all $2d$ nearest neighbours are accessed.

§.§.§ Allowed transitions and rates

In (<ref>) above, $\frac{\partial E(\textbf{x})}{\partial t}$ was proportional to $\epsilon$, the
smallest amount of energy allowed to move in one timestep. Therefore, we should presumably add an equation for
momentum where $\frac{\partial \textbf{P}(\textbf{x})}{\partial t}$ is proportional to some $\textbf{q}$, the smallest
momentum vector that moves in one timestep. This can be written as $q \textbf{e}$, the length times some unit vector,
where $q \rightarrow 0$ just like $\epsilon$ and $a$.

However, $q$ and $\epsilon$ are not independent. A CFT in $d = 1$ for instance has energy and momentum operators
that can be written in terms of left movers and right movers:

\begin{eqnarray}
E &=& N_{\textup{R}} + N_{\textup{L}} - \frac{c}{12} \nonumber \\
P &=& N_{\textup{R}} - N_{\textup{L}} \; . \nonumber
\end{eqnarray}

The only way to add to $E$ is to excite either $N_{\textup{R}}$ or $N_{\textup{L}}$ with $\epsilon$ units of energy
and this necessarily adds $\epsilon$ to or subtracts $\epsilon$ from $P$. In this case we must have $q = \epsilon$. We
will make the assumption that the other theories we wish to model also have the energy transferred equal to the
magnitude of the momentum transferred.

Now that we know $E$ increases or decreases by $\epsilon$ whenever $\textbf{P}$ increases or decreases by
$\epsilon \textbf{e}$, an allowed transition takes the form

\begin{equation}
(E(\textbf{x}), \textbf{P}(\textbf{x}), E(\textbf{x} + a\textbf{e}), \textbf{P}(\textbf{x} + a\textbf{e})) \rightarrow (E(\textbf{x}) \pm \epsilon, \textbf{P}(\textbf{x}) - \epsilon\textbf{e}^{\prime}, E(\textbf{x} + a\textbf{e}) \mp \epsilon, \textbf{P}(\textbf{x} + a\textbf{e}) + \epsilon\textbf{e}^{\prime}) \nonumber
\end{equation}

where $\textbf{e}$ is a unit vector indicating the two sites involved and $\textbf{e}^{\prime}$ is a priori
some other unit vector. However, it does not make sense for momentum in the vertical direction to move between sites
that are horizontally displaced and it does not make sense for momentum pointing in the up direction to move down.
Therefore we require that $\textbf{e}^{\prime} = \textbf{e}$. This is where we see that the vector of large momentum
components $\textbf{P}$ lives in the same number of dimensions as the lattice vector $\textbf{x}$.

When focusing on energy distributions, the form of the transition rates in (<ref>) followed from
thermodynamic arguments as

\begin{equation}
W_{(E(\textbf{x}), E(\textbf{x} + a\textbf{e})) \rightarrow (E(\textbf{x}) + \epsilon, E(\textbf{x} + a\textbf{e}) - \epsilon)} = C\left ( \frac{E(\textbf{x}) + E(\textbf{x} + a\textbf{e})}{2} \right ) \rho(E(\textbf{x}) + \epsilon) \rho(E(\textbf{x} + a\textbf{e}) - \epsilon) \label{first-transition-rates}
\end{equation}

with the number of states having energy $E$ given by $\rho(E)$. The analogous transition rate involving momentum
should depend on $\rho(E, \textbf{P})$ in the same way. We are led to:

\begin{eqnarray}
&& W_{(E(\textbf{x}), \textbf{P}(\textbf{x}), E(\textbf{x} + a\textbf{e}), \textbf{P}(\textbf{x} + a\textbf{e})) \rightarrow (E(\textbf{x}) + \epsilon, \textbf{P}(\textbf{x}) - \epsilon\textbf{e}^{\prime}, E(\textbf{x} + a\textbf{e}) - \epsilon, \textbf{P}(\textbf{x} + a\textbf{e}) + \epsilon\textbf{e}^{\prime})} = \nonumber \\
&& \delta_{\textbf{e}, \textbf{e}^{\prime}} C\left ( \frac{E(\textbf{x}) + E(\textbf{x} + a\textbf{e})}{2}, \frac{\textbf{P}(\textbf{x}) + \textbf{P}(\textbf{x} + a\textbf{e})}{2} \right ) \nonumber \\
&& \rho(E(\textbf{x}) + \epsilon, \textbf{P}(\textbf{x}) - \epsilon\textbf{e}^{\prime}) \rho(E(\textbf{x} + a\textbf{e}) - \epsilon, \textbf{P}(\textbf{x} + a\textbf{e}) + \epsilon\textbf{e}^{\prime}) \; . \label{second-transition-rates}
\end{eqnarray}

§.§.§ A tensor identity

Consider what happens when we differentiate $E(\textbf{x} + a\textbf{e})$ twice with respect to $a$.

\begin{eqnarray}
\frac{\partial}{\partial a} E(\textbf{x} + a\textbf{e}) &=& e_i \partial_i E(\textbf{x} + a\textbf{e}) \nonumber \\
\frac{\partial^2}{\partial a^2} E(\textbf{x} + a\textbf{e}) &=& e_i e_j \partial_i \partial_j E(\textbf{x} + a\textbf{e}) \label{example-derivative}
\end{eqnarray}

Note that $e_i$ and $e_j$ without the boldface are not the $i^{\mathrm{th}}$ and $j^{\mathrm{th}}$ standard basis
vectors, they are the $i^{\mathrm{th}}$ and $j^{\mathrm{th}}$ components of the general unit vector $\textbf{e}$.
All of the terms that we must differentiate with respect to $a$ occur inside a sum over $\textbf{e}$ so we must know
how to deal with sums of components of unit vectors. When we sum $e_i e_j$, there is only one standard basis vector
for which $e_i$ is nonzero and only one standard basis vector for which $e_j$ is nonzero. They have to be the same one
for the resulting sum to be nonzero, so it is clear that:

\begin{equation}
\sum_{\textbf{e} \in \{ \pm \textbf{e}_1, \dots, \pm \textbf{e}_d \}} e_i e_j = 2\delta_{ij} \; . \nonumber
\end{equation}

We used this fact when deriving the PDE for our first model.
For reasons that are not initially clear, we will show that one can get the same answer up to a constant by replacing
the sum over lattice vectors with a surface integral over a sphere.

\begin{equation}
\int_{\mathbb{S}^{d - 1}} e_i e_j \textup{d}S_{\textbf{e}} = \omega_d \delta{ij} \nonumber
\end{equation}

where $\omega_d$ is the volume of the unit ball in $\mathbb{R}^d$.

If $i = j$, we simply apply the divergence theorem to the identity vector field. Let
$\textbf{v}(\textbf{e}) = \textbf{e}$, defined on the unit ball $B^d$. Then,

\begin{equation}
\int_{B^d} \nabla \cdot \textbf{v} \textup{d}\textbf{e} = \int_{\partial B^d} \textbf{v} \cdot \textbf{n} \textup{d}S_{\textbf{e}} \nonumber
\end{equation}

for this vector field reads

\begin{equation}
d \int_{B^d} \textup{d}\textbf{e} = \int_{\mathbb{S}^{d - 1}} (e_1)^2 + \dots + (e_d)^2 \textup{d}S_{\textbf{e}} \nonumber
\end{equation}

because if the unit vector $\textbf{e}$ extends from the origin to a point on $\mathbb{S}^{d - 1}$, the unit normal
to the surface at this point is just $\textbf{e}$. By symmetry, the right hand side is $d$ copies of what we are
trying to compute, so this cancels the $d$ on the left hand side yielding:

\begin{equation}
\omega_d = \int_{\mathbb{S}^{d - 1}} (e_i)^2 \textup{d}S_{\textbf{e}} \; . \nonumber
\end{equation}

To see that this integral vanishes when $i$ and $j$ are different, parametrize the unit sphere as

\begin{eqnarray}
x_1 &=& \cos \theta_1 \nonumber \\
x_2 &=& \sin \theta_1 \cos \theta_2 \nonumber \\
&\dots& \nonumber \\
x_{d - 1} &=& \sin \theta_1 \dots \sin \theta_{d - 2} \cos \theta_{d - 1} \nonumber \\
x_d &=& \sin \theta_1 \dots \sin \theta_{d - 2} \sin \theta_{d - 1} \nonumber
\end{eqnarray}

and without loss of generality, choose $e_i = x_d$ and $e_j = x_{d - 1}$. The desired integral is then

\begin{equation}
\int_0^{2\pi} \int_0^{\pi} \dots \int_0^{\pi} [ \sin^2 \theta_1 \dots \sin^2 \theta_{d - 2} \sin \theta_{d - 1} \cos \theta_{d - 1} ] \sin^{d - 2} \theta_1 \sin^{d - 3} \theta_2 \dots \sin \theta_{d - 2} \textup{d}\theta_1 \dots \textup{d}\theta_{d - 2} \textup{d}\theta_{d - 1} \nonumber
\end{equation}

where the part in brackets is $e_i e_j$. Since the area element only goes up to $\theta_{d - 2}$, the part involving
$\theta_{d - 1}$ is just a sine and a cosine which we know integrates to zero.

Unlike in our first model, equations involving momenta will have several places in which there are four
components of $\textbf{e}$. This will make it important that we integrate rather than sum. A sum of four components
will give zero unless the indices are all the same.

\begin{equation}
\sum_{\textbf{e} \in \{ \pm \textbf{e}_1, \dots, \pm \textbf{e}_d \}} e_i e_j e_k e_l =
\begin{cases}
2 & i = j = k = l \\
0 & \mathrm{otherwise}
\end{cases} \nonumber
\end{equation}

The continuous version, on the other hand, is a much nicer object.

\begin{equation}
\int_{\mathbb{S}^{d - 1}} e_i e_j e_k e_l \textup{d}S_{\textbf{e}} = \frac{\omega_d}{d + 2} \left ( \delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} \right ) \nonumber
\end{equation}

where $\omega_d$ is the volume of the unit ball in $\mathbb{R}^d$.

We will first prove the special case when $i = j = k = l$. Define a vector field
$\textbf{v} : B^d \rightarrow \mathbb{R}^d$ which cubes every component of its argument:

\begin{equation}
v_i(\textbf{e}) = (e_i)^3 \; . \nonumber
\end{equation}

Using the divergence theorem on this, we get:

\begin{eqnarray}
\int_{B^d} \nabla \cdot \textbf{v} \textup{d}\textbf{e} &=& \int_{\partial B^d} \textbf{v} \cdot \textbf{n} \textup{d}S_{\textbf{e}} \nonumber \\
3 \int_{B^d} (e_1)^2 + \dots (e_d)^2 \textup{d}\textbf{e} &=& \int_{\mathbb{S}^{d - 1}} (e_1)^4 + \dots + (e_d)^4 \textup{d}S_{\textbf{e}} \nonumber \\
3d \omega_d \int_0^1 r^{d + 1} \textup{d}r &=& d \int_{\mathbb{S}^{d - 1}} (e_i)^4 \textup{d}S_{\textbf{e}} \nonumber \\
\frac{3 \omega_d}{d + 2} &=& \int_{\mathbb{S}^{d - 1}} (e_i)^4 \textup{d}S_{\textbf{e}} \; . \nonumber
\end{eqnarray}

We will now apply this to the case when $i$ and $k$ are different.

\begin{eqnarray}
\sum_{k \neq i} \int_{\mathbb{S}^{d - 1}} (e_i)^2 (e_k)^2 \textup{d}S_{\textbf{e}} &=& \sum_{k = 1}^d \int_{\mathbb{S}^{d - 1}} (e_i)^2 (e_k)^2 \textup{d}S_{\textbf{e}} - \int_{\mathbb{S}^{d - 1}} (e_i)^4 \textup{d}S_{\textbf{e}} \nonumber \\
&=& \int_{\mathbb{S}^{d - 1}} (e_i)^2 \textup{d}S_{\textbf{e}} - \int_{\mathbb{S}^{d - 1}} (e_i)^4 \textup{d}S_{\textbf{e}} \nonumber \\
&=& \omega_d - \frac{3 \omega_d}{d + 2} \nonumber \\