03a2BoltzmnEinst.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Boltzmann-Einstein equation}\label{sec:BoltzmannEinstein}
We now begin a detailed study of the nonequilibrium properties of the neutrino freeze-out\index{neutrino!freeze-out} and its impact on the effective number of neutrinos\index{neutrino!effective number}, an important cosmological observable. We model the dynamics of the neutrino freeze-out using the Boltzmann-Einstein equation\index{Boltzmann-Einstein equation}, also called the general relativistic Boltzmann equation, which describes the dynamics of a gas of particles that travel on geodesics in an general spacetime, with the only interactions being point collisions~\cite{Andreasson:2011ng,cercignani,Choquet-Bruhat:2009xil,ehlers},
\begin{equation}\label{eq:BoltzmannEinstein}
p^\alpha\partial_{x^\alpha}f-\sum_{j=1}^3\Gamma^j_{\mu\nu}p^\mu p^\nu\partial_{p^j}f=C[f]\,.
\end{equation}
Here $ \Gamma^\alpha_{\mu\nu}$ is the affine connection (Christoffel symbols) corresponding to a metric $g_{\alpha\beta}$, and the distribution function $f$ is a function of four-momentum on the mass shell, \ie, that satisfies
 \begin{equation}
g_{\alpha\beta}p^\alpha p^\beta=m^2\,.
\end{equation}
Here and in the following, repeated Greek indices are summed from $0$ to $3$. $C[f]$ is the collision operator\index{collision operator} and encodes all information about point interactions between particles. 

If $C[f]$ vanishes, then the equation is called the Vlasov equation and describes particles that move on geodesics (or free stream). At this point, we are not invoking the assumption that the distribution function has a kinetic equilibrium form\index{kinetic equilibrium}, nor are we assuming a FLRW\index{cosmology!FLRW} Universe; in this section we will discuss general properties of \req{eq:BoltzmannEinstein} before turning to the study of neutrino freeze-out in subsequent sections. We will need the following definitions of entropy current $S^\mu$, stress-energy tensor ${T}^{\mu\nu}$, and number current $N^\mu$,\index{entropy!current}\index{stress-energy tensor}\index{number current}
\begin{align}
\label{smdef} S^\mu&=-\int \left(f\ln(f)\pm(1\mp f)\ln(1\mp f)\right)p^\mu d\pi\,,\\
\label{Tmndef}{T}^{\mu\nu}&=\int p^\mu p^\nu f d\pi\,,\\
\label{nmdef} N^\nu&=\int f p^\nu d\pi\,,\\
d\pi&=\frac{\sqrt{-g}}{p_0}\frac{g_pd^3{\bf p}}{8\pi^3}\,,
\end{align}
where $d\pi$ is the volume element on the future mass shell; here $g$ denotes the determinant of the metric tensor, $p_0=g_{0\alpha} p^\alpha$; non-bold $p$ are four-momenta while bold ${\bf p}$ denotes the spacial components; the upper signs are for fermions and the lower signs for bosons. See \rapp{ch:vol:forms} for the derivation of the form of the volume element.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\para{Collision operator}
We now elaborate on the form of the collision operator. Our presentation is an expanded version of the survey in \cite{ehlers}. Suppose we have a collection of distinct particle and antiparticle types $\mathcal{C}$ with distribution functions $f_{C}$, $C\in\mathcal{C}$, and they partake in some number of reactions or interactions $I=n_{B_1} B_1, n_{B_2}B_2 \ldots \longrightarrow n_{A_1} A_1,n_{A_2}A_2 \ldots $, $A_i\in\mathcal{C}$ distinct and $B_j\in\mathcal{C}$ distinct, where $n_{A_i}$ is the number of particles of type $A_i$ occurring in the interaction (all nonzero) and similarly for $n_{B_i}$. Given an interaction, $I$, we let $r(I)$ be the collection of particle types that are reactants in the interaction, $p(I)$ be the collection of particle types that are products, and we let $\overleftarrow{I}$ denote the reverse reaction, \ie, with reactants and products reversed. We let $int$ denote the set of all interactions and, for any given species $A$, $int(A)$ be the set of all interactions involving $A$ as a product. We will assume that $\overleftarrow{I}\in int$ whenever $I\in int$. With these conventions, the collision operator for particle type $A$ takes the form\index{collision operator}
\begin{align}\label{collision:operator}
C[f_A]=&\sum_{I\in int(A)} \frac{n_A}{\prod_i n_{A_i}!\prod _j n_{B_j}!}\int\left[\left(\prod_j \prod_{l=1}^{n_{B_j}}f_{B_j}(p_{B_j}^l)\right)\left(\prod_i \prod_{k=1}^{n_{A_i}}f^{A_i}(p_{A_i}^k)\right)W^I(p_{B_j}^l,p_{A_i}^k) \right.\notag\\[0.3cm]
 -& \left.\left(\prod_i \prod_{k=1}^{n_{A_i}}f_{A_i}(p_{A_i}^k)\right)\left(\prod_j \prod_{l=1}^{n_{B_j}}f^{B_j}(p_{B_j}^l)\right)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l) \right] \delta(\Delta p)\prod_i \widehat{dV}_{A_i}\prod_j dV_{B_j}\,
\end{align}
where
\begin{align}
f^C=&1\mp f_C, \hspace{2mm} 
\Delta p=\sum_i \sum_{k=1}^{n_{A_i}}p^k_{A_i}-\sum_j \sum_{l=1}^{n_{B_j}}p^l_{B_j}\,,  \hspace{2mm}
\widehat{dV}_{A_i}=\tilde{\pi}_{A_i}\prod_{k=2}^{n_{A_i}}\frac{1}{2}d\pi^k_{A_i}\,, \hspace{2mm} 
dV_{B_j}=(2\pi)^4\prod_{l=1}^{n_{B_j}}\frac{1}{2}d\pi^l_{B_j}\,,\notag\\
 \tilde{\pi}_{A_i}=&\frac{1}{2} \text{ if } A_i=A \text{ and } \tilde{\pi}_{A_i}=\frac{1}{2}d\pi^1_{A_i} \text{ otherwise\,,}\hspace{2mm} 
d\pi_{C}^r=\frac{\sqrt{-g}}{(p_{C}^r)_0}\frac{g_{C}d^3{\bf p}_{C}^r}{8\pi^3}, p_0=g_{0\alpha}p^\alpha.\notag
\end{align}
 The integrations are over the future mass shells of all the particles, so the $p$ are related by $g_{\alpha \beta}p^\alpha p^\beta=m^2$. 
 
 The factorials in the collision term take into account the indistinguishability of the particles and prevent one from over counting the independent ways a reaction can happen when integrating over momentum. The terms $f^A$ are due to quantum statistics and account for Fermi repulsion or Bose attraction (again, upper signs are for fermions and lower signs for bosons). $W^I(p_{B_j}^l,p_{A_i}^k)$, an abbreviation for $W^I(p_{B_1}^1,p_{B_1}^2,\ldots ,p_{B_1}^{n_{B_1}},p_{B_2}^1,\ldots ,p_{A_1}^1,\ldots )$, is the scattering kernel that encodes the probability of $n_{B_j}$ particles of types $B_j$ with momenta $p_{B_j}^l$ interacting to form $n_{A_i}$ particles of types $A_i$ with momenta $p_{A_i}^k$ in the process $I=n_{B_1}B_1,n_{B_2}B_2,\ldots\longrightarrow n_{A_1}A_1,n_{A_1}A_1,\ldots $, and so it is non-negative. 
 
 The delta function in the collision term enforces conservation of four-momentum. The factors of $(2\pi^4)$ and $\frac{1}{2}$ in the definitions of the volume elements come from normalization of the transition functions from quantum scattering calculations. For computational purposes, the expression \eqref{collision:operator} must be further simplified, taking into account the structure of each interaction. For example, see \rapp{ch:coll:simp} for a detailed study of the collision operator in the case of neutrino freeze-out.

 As defined, $C[f_A]$ is a function of $p_{A_i}^1$ where $A=A_i$. The choice not to integrate over $p_{A_i}^1$ rather than any of the other $p_{A_i}^k$ is completely arbitrary, but makes no difference in the result since the interaction does not depend on how we number the participating particles. In terms of the scattering kernels, this means we assume $W^I$ has the property
\begin{equation}\label{reorderProperty}
W^I(p_{A_1}^{\sigma_1},p^{\sigma 2}_{A_1},\ldots )=W^I(p_{A_1}^1,p_{A_1}^2,\ldots )\,,
\end{equation}
for any permutation $\sigma$, and similarly for any other permutation with one of the collections $p_{A_i}^k$ or $p_{B_j}^l$ for any choice of $i$ or $j$. For economy of notation in these derivations, we will employ the additional abbreviations for a given interaction $I=n_{B_i}B_i\longrightarrow n_{A_i}A_i$:
\begin{align}
f_{p,I}(p^k_{A_i})&\equiv f_{p,I}(p^1_{A_i},p^2_{A_i},\ldots ,p^{n_{A_i}}_{A_i})\equiv \prod_i \prod_{k=1}^{n_{A_i}}f_{A_i}(p_{A_i}^k)\,,\\
f^{p,I}(p^k_{A_i})&=f^{p,I}(p^1_{A_i},p^2_{A_i},\ldots ,p^{n_{A_i}}_{A_i})=\prod_i \prod_{k=1}^{n_{A_i}}f^{A_i}(p_{A_i}^k)\,,\notag\\
f_{r,I}(p^l_{B_j})&\equiv f_{r,I}(p^1_{B_j},p^2_{B_j},\ldots ,p^{n_{B_j}}_{B_j})\equiv \prod_j \prod_{l=1}^{n_{B_j}}f_{B_j}(p_{B_j}^l)\,,\notag\\
f^{r,I}(p^l_{B_j})&=f^{r,I}(p^1_{B_j},p^2_{B_j},\ldots ,p^{n_{B_j}}_{B_j})=\prod_j \prod_{l=1}^{n_{B_j}}f^{B_j}(p_{B_j}^l)\,,\notag\\
n_I&=\prod_i n_{A_i}!\prod _j n_{B_j}!\,,\notag\\
\widehat{dV}_I&=\delta(\Delta p)\prod_i\widehat{dV}_{A_i}\prod_jdV_{B_j}\,,\notag\\
dV_I&=\delta(\Delta p)\prod_idV_{A_i}\prod_jdV_{B_j}\,,\notag
\end{align}
where the $r$ and $p$ sub and superscripts stand for reactants and products respectively. See \rapp{ch:vol:forms} for more information on the precise meaning and properties of the delta function factors.

In the following subsections we derive several important properties of the equation \eqref{eq:BoltzmannEinstein}. While in principle these properties are well known~\cite{Andreasson:2011ng,cercignani,Choquet-Bruhat:2009xil,ehlers}, here we prove them at a level of generality that, to the authors knowledge, is not available in other references, \ie, for a general collection of interactions as encapsulated in \req{collision:operator}. We note that Riemannian normal coordinates will a key tool in these derivations. These are coordinates centered at a chosen point, $x$, in spacetime wherein the geodesics through $x$ are straight lines in the coordinate system and the derivatives of the metric in the coordinate system vanish at $x$. In particular, the Christoffel symbols vanish at $x$; see, e.g., page 42 in \cite{Wald:1984rg} or pages 72-73 of \cite{o1983semi}. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\para{Conserved currents}
Suppose all the interactions of interest conserve some charge $b_A$, \ie,
 \begin{align}\label{eq:conservedCharge}
\sum_{A\in p(I)} n_Ab_A=\sum_{A\in r(I)} n_Ab_A
\end{align}
for all $I\in int$. We can construct a $4$-vector current corresponding to this charge as follows:
\begin{equation}
B^\mu=\sum_A b_A N_A^\mu\,,
\end{equation}
where $N^\mu_A$ are the number currents of the particle species \req{nmdef}. In this section we show that $B^\mu$ has vanishing divergence, \ie, a $B^\mu$ satisfies a conservation law.

For any point $x$ in spacetime, by transforming to Riemannian normal coordinates at $x$ and using \eqref{eq:BoltzmannEinstein} along with the fact that the first derivatives of the metric vanish at $x$, one can compute
\begin{equation}\label{useNormalCoords}
\nabla_\mu N_A^\mu=\int p^\mu \partial_{x^\mu} f d\pi_A=\int C[f_A] d\pi_A
\end{equation}
at $x$. The left and right hand sides are scalars and therefore they are equal in any coordinate system. Noting this, we can then calculate
\begin{align}\label{eq:charge_current_div}
\nabla_\mu B^\mu=&\sum_A b_A\int C[f_A]d\pi_A\notag\\
=&\sum_A\sum_{I\in int(A)} \frac{n_Ab_A}{n_I}\int\int\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) %\right.\notag\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right)\widehat{dV}_I d\pi_A\notag\\
=&\sum_A\sum_{I\in int(A)} \frac{n_Ab_A}{n_I}\int\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) %\right.\notag\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\;. 
\end{align}
Now observe that, for any collection of finite sets $D_j$ indexed by a finite set $J$ with $\bigcup_{j\in J}D_j=D$ and any function $h:J\times D\rightarrow \mathbb{R}^m$, we have
\begin{equation}\label{sumLemma}
\sum_{j\in J}\sum_{x\in D_j} h(j,x)=\sum_{x\in D}\sum_{\{j:x\in D_j\}}h(j,x)\,.
\end{equation}
Using this fact, we can switch the order of the sums to obtain
\begin{align}\label{delB}
&\nabla_\mu B^\mu=\sum_{I\in int}\sum_{A\in p(I)} n_Ab_A R_I \,,\\
&R_I\equiv \frac{1}{n_I}\int\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) -f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\,.\notag
\end{align}
The sum over all interactions splits over a sum over symmetric interactions, $int_{s}$, and a sum over asymmetric interactions. For each asymmetric interaction, pair it up with its reverse and arbitrarily choose one of them to call the forward direction. Let the set of these forward interactions be denoted $\overrightarrow{int}$. Then the sum in \req{delB} splits as follows
\begin{equation}
\nabla_\mu B^\mu=\sum_{I\in int_s}R_I\sum_{A\in p(I)} n_A b_A+\sum_{I\in\overrightarrow{int}}R_I\sum_{A\in p(I)} n_Ab_A+\sum_{I\in\overrightarrow{int}}R_{\overleftarrow{I}}\sum_{A\in p(\overleftarrow{I})} n_Ab_A\,.
\end{equation}
For every $I\in int_s$ we have $W^I=W^{\overleftarrow{I}}$, $f_{A_i}=f_{B_i}$, and $f^{A_i}=f^{B_i}$, and therefore
\begin{align}
R_I=&\frac{1}{n_I}\left(\int f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) dV_I 
%\right.\\&\left. 
-\int f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l) dV_I\right)\notag\\
=&\frac{1}{n_I}\left(\int f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) dV_I %\right.\notag\\&\left. 
-\int f_{r,I}(p_{A_i}^k)f^{p,I}(p_{B_j}^l)W^{I}(p_{A_i}^k,p_{B_j}^l) dV_I\right)
%\notag\\=&0\,,\notag
=0\,,
\end{align}
as the two integrals differ only by a relabeling of integration variables. Asymmetric interactions satisfy
\begin{align}
R_{\overleftarrow{I}}=&\frac{1}{n_I}\int\left(f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l) 
-f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k)\right) dV_I\notag\\
=&-R_I.
\end{align}
Combining this with \req{eq:conservedCharge}, we find
\begin{align}
\nabla_\mu B^\mu=\sum_{I\in\overrightarrow{int}}R_I\left(\sum_{A\in p(I)} n_Ab_A-\sum_{A\in p(\overleftarrow{I})} n_Ab_A\right)
=\sum_{I\in\overrightarrow{int}}R_I\left(\sum_{A\in p(I)} n_Ab_A-\sum_{A\in r(I)} n_Ab_A\right)=0\,.
\end{align}
Therefore $B^\mu$ is a conserved current, as claimed.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\para{Divergence freedom of stress energy tensor}
The Einstein equation implies that the total stress energy tensor\index{stress-energy tensor} of all matter coupled to gravity is divergence free. Here we show that the relativistic Boltzmann stress energy tensor \req{Tmndef} has this property, and is therefore a natural candidate matter model for coupling to gravity. 

First use Riemannian normal coordinates to compute
\begin{align}
\nabla_\mu T^{\mu\nu}=&\sum_A \int p_A^\nu C[f_A]d\pi_A \\
=&\sum_A\sum_{I\in int(A)}\frac{n_A}{n_I}\int (p^1_{A_\ell})^\nu\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) 
%\right.\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\,,\notag
\end{align}
where $\ell$ is the unique index such that $A_\ell=A$ ($\ell$ depends on $A$ and $I$, but we suppress this dependence for simplicity of notation). Using \req{sumLemma} we can switch the summation order to get
\begin{align}
\nabla_\mu T^{\mu\nu}=&\sum_{I\in int}\sum_{A\in p(I)} \frac{n_A}{n_I}\int (p^1_{A_\ell})^{\nu}\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k)
%\right.\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\,.
\end{align}
By \req{reorderProperty} and the surrounding remarks, we can rewrite this as
\begin{align}\label{delTSum}
\nabla_\mu T^{\mu\nu}=&\sum_{I\in int}\sum_{A\in p(I)} \frac{1}{n_I}\sum_{a=1}^{n_A}\int (p^a_{A_\ell})^{\nu}\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) 
%\right.\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\notag\\
=&\sum_{I\in int} \frac{1}{n_I}\sum_{\ell}\sum_{a=1}^{n_{A_\ell}}\int (p_{A_\ell}^a)^{\nu}\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) %\right.\notag\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\,.
\end{align}
As before, we can break the sum over $I$ into a sum over symmetric processes and two other sums over forward and backward asymmetric processes, respectively. For a symmetric interaction $I=\overleftarrow{I}$ and $f_{A_i}=f_{B_i}$ for all $i$, hence
\begin{align}
&\sum_\ell\sum_{a=1}^{n_{A_\ell}}\int (p_{A_\ell}^a)^{\nu}\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) 
%\right.\\ &\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\notag\\
&\hspace{2cm}=\int\sum_\ell\sum_{a=1}^{n_{A_\ell}}\left((p_{A_\ell}^a)^\nu- (p_{B_\ell}^a)^{\nu}\right) f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) dV_I
%\notag\\=& 0\,,\notag
=0\,,
\end{align}
due to the delta function $\delta(\Delta p)$ in the volume form $dV_I$. Therefore the terms in the sum \req{delTSum} corresponding to symmetric interactions vanish. For every pair of forward and backward asymmetric interactions we obtain
\begin{align}
&\sum_\ell\sum_{a=1}^{n_{A_\ell}}\int (p_{A_\ell}^a)^{\nu}\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) 
%\right.\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\notag\\
&+\sum_{\tilde\ell}\sum_{c=1}^{n_{B_{\tilde\ell}}} \int (p_{B_{\tilde\ell}}^c)^{\nu}\left(f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l) 
%\right.\notag\\&\left. 
-f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k)\right) dV_I\notag\\
=&\int\left( \sum_\ell\sum_{a=1}^{n_{A_\ell}}(p_{A_\ell}^a)^{\nu} -\sum_{\tilde \ell}\sum_{c=1}^{n_{B_{\tilde\ell}}} (p_{B_{\tilde\ell}}^c)^{\nu}\right)f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) dV_I\notag\\
&+\int \left(\sum_{\tilde\ell}\sum_{c=1}^{n_{B_{\tilde\ell}}}(p_{B_{\tilde\ell}}^c)^{\nu}-\sum_\ell\sum_{a=1}^{n_{A_\ell}}(p_{A_\ell}^a)^{\nu}\right)f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l) dV_I
%\notag\\=&0\,,\notag
=0\,,
\end{align}
again because of $\delta(\Delta p)$ in the volume forms. This shows $\nabla_\mu T^{\mu\nu}=0$, as claimed.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\index{Boltzmann!H-theorem}
\para{Entropy and Boltzmann's H-Theorem}
Finally, we prove that the entropy four-current satisfies $\nabla_\mu S^\mu\geq 0$, known as Boltzmann's H-theorem. This result requires the additional assumption that the interactions are time-reversal symmetric, \ie,
\begin{equation}\label{timeSymmetry}
W^I(p_{B_j}^l,p_{A_i}^k)=W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)
\end{equation}
for all $I$. Working in Riemannian normal coordinates once again, we can compute
\begin{align}
\nabla_\mu S^\mu=&-\sum_A\int p^\mu \partial_{x^\mu}\left(f_A\ln\left(f_A\right)\pm\left(1\mp f_A\right)\ln\left(1\mp f_A\right)\right)d\pi_A
%\\=&
=\sum_A\int\ln\left(1/f_A\mp 1\right)C[f_A]d\pi_A\,. 
\end{align}
Similar reasoning to the above two subsections then gives
\begin{align}
\nabla_\mu S^\mu=&\sum_{I\in int}\frac{1}{n_I}\sum_\ell\sum_{a=1}^{n_{A_\ell}}\int\ln\left(1/f_{A_\ell}(p_{A_\ell}^a)\mp1\right)\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) \right.\\
&\left. -f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\,.\notag
\end{align}
Once again, we break the summation into a sum over symmetric processes and two other sums over forward and backward asymmetric processes respectively. Each symmetric process contributes a term of the form
\begin{align}
&\int\sum_\ell\sum_{a=1}^{n_{A_\ell}}f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)\left(\ln\left(1/f_{A_\ell}(p_{A_\ell}^a)\mp 1\right)
%\right.\\&\left.\qquad\qquad\qquad\qquad
-\ln\left(1/f_{B_\ell}(p_{B_\ell}^a)\mp 1\right)\right) W^I(p_{B_j}^l,p_{A_i}^k) dV_I\notag\\
=& \int\ln\left(\frac{f^{p,I}(p_{A_i}^k)f_{r,I}(p_{B_j}^l)}{f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)}\right) f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) dV_I\notag\\
=& \frac{1}{2}\int\ln\left(\frac{f^{p,I}(p_{A_i}^k)f_{r,I}(p_{B_j}^l)}{f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)}\right)\left( f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)
%\right.\notag\\&\left.\qquad\qquad\qquad\qquad
- f_{p,I}(p_{A_j}^l)f^{r,I}(p_{B_i}^k)\right)W^I(p_{B_j}^l,p_{A_i}^k) dV_I\,,
\end{align}
where to obtain the last line we used the time-reversal property \eqref{timeSymmetry}.

A pair of forward and backward asymmetric interactions combine to give a term of the form
\begin{align}
&\sum_\ell\sum_{a=1}^{n_{A_b}}\int\ln\left(1/f_{A_\ell}(p_{A_\ell}^a)\mp 1\right)\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k)
%\right.\\&\left. 
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l)\right) dV_I\notag\\
&+\sum_{\tilde\ell}\sum_{c=1}^{n_{B_{\tilde\ell}}} \int\ln\left(1/f_{B_{\tilde\ell}}(p_{B_{\tilde\ell}}^c)\mp 1\right)\left(f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{\overleftarrow{I}}(p_{A_i}^k,p_{B_j}^l) 
%\right.\notag\\&\left. 
-f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k)\right) dV_I\notag\\
=&\int\left( \sum_\ell\sum_{a=1}^{n_{A_\ell}}\ln\left(1/f_{A_\ell}(p_{A_\ell}^a)\mp 1\right) 
%\right.\notag\\&\left.\qquad\qquad\qquad
-\sum_{\tilde\ell}\sum_{c=1}^{n_{B_{\tilde\ell}}} \ln\left(1/f_{B_{\tilde\ell}}(p_{B_{\tilde\ell}}^c)\mp 1\right)\right)f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)W^I(p_{B_j}^l,p_{A_i}^k) dV_I\notag\\
&-\int \left(\sum_\ell\sum_{a=1}^{n_{A_\ell}}\ln\left(1/f_{A_\ell}(p_{A_\ell}^a)\mp 1\right) %\right.\notag\\&\left.\qquad\qquad\qquad
-\sum_{\tilde\ell}\sum_{c=1}^{n_{B_{\tilde\ell}}} \ln\left(1/f_{B_{\tilde\ell}}(p_{B_{\tilde\ell}}^c)\mp 1\right)\right)f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)W^{I}(p_{B_j}^l,p_{A_i}^k) dV_I\,,
\end{align}
where to obtain the first equality we used the time-reversal property \eqref{timeSymmetry}. Combining the symmetric and asymmetric cases we find
\begin{align}\label{eq:entropy_current}
\nabla_\mu S^\mu=&\sum_{I\in int_s} \frac{1}{2n_I}\int\ln\left(\frac{f^{p,I}(p_{A_i}^k)f_{r,I}(p_{B_j}^l)}{f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)}\right)\left( f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)
%\right.\\&\left.
- f_{p,I}(p_{A_j}^l)f^{r,I}(p_{B_i}^k)\right)W^I(p_{B_j}^l,p_{A_i}^k) dV_I\notag\\
&+\sum_{I\in \overrightarrow{int}}\frac{1}{n_I}\int\ln\left(\frac{f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)}{f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)}\right)\left(f_{r,I}(p_{B_j}^l)f^{p,I}(p_{A_i}^k)
%\right.\notag\\&\left.
-f_{p,I}(p_{A_i}^k)f^{r,I}(p_{B_j}^l)\right)W^I(p_{B_j}^l,p_{A_i}^k) dV_I\,.
\end{align}
Each term in either sum is the integral of a non-negative quantity $W^{I}$ times a quantity of the form $(a-b)\ln(a/b)$, $a,b>0$, which is easily seen to be non-negative. Therefore we obtain the claimed result $\nabla_\mu S^\mu\geq 0$. The entropy four-current is future directed, due to the volume element being supported on the future mass shell. Therefore, given any splitting of spacetime into space and time $M=S\times T$, Boltzmann's H-theorem implies that the total entropy is non-decreasing on $T$.