The Liouville equation and the BBGKY hierarchy. The Liouville equation of motion is written in terms of an $N$ particle distribution $f_N$. 
\begin{equation}
\frac{\partial f_N}{\partial t}=\{H,f_N\}
\end{equation}
Where $\{\cdot ,\cdot \}$ is the Poisson bracket and $f_N=f_N(q_1,\dots ,q_N,p_1,\dots ,p_N)$. Let us now define an $n$ particle probability distribution function $f_n$ with $n< N$. 
\begin{equation}
f _n(q_1,\dots,q_n,p_1,\dots ,p_n,t)=\frac{N!}{(N-n)!}\int \prod ^N_{i=n+1}d q^id p_if_{N+1}( q_1,\dots ,q_{N+1},p_1 ,\dots,p_{N+1},t)
\end{equation}
Now $f_n$ satisfies,
\begin{equation}
\frac{\partial f _n}{\partial t}=\{H_n,f_n\}+\sum ^n_{i=1}\int dq_{n+1}dp_{n+1}\frac{\partial U(q_i-q_{n+1})}{\partial q_i}\frac{\partial f _{n+1}}{\partial p_i}\ \ \ \ \ \ \ \ \ \ \ \ \ \ (*)
\end{equation}
With the $n$-body Hamiltonian $H_n$,
\begin{equation}
H_n=\sum ^n_{i=1}\bigg(\frac{p_i^2}{2m}+U(q_i)\bigg)+\sum _{i<j\leq n}U(q_{ij})
\end{equation}
And $q_{ij}=q_i-q_j$. Here $(*)$ is the BBGKY hierarchy. I am reading out of the following notes.
Despite reading the linked notes and wikipedia pages etc I am struggling to understand how the BBGKY hierarchy is related to the Liouville equation. In particular taking $n=N$ does not (to my naive understanding) regenerate the Liouville equation. Why do we not require $f_{N+1}$ in the Liouville equation by the logic of the $n$ particle distribution function? Lastly is the Boltzmann equation defined for $f_N$ or $f_1$ (or is it irrelevant, the equation holding in any case?).  
Any help on the BBGKY formalism is appreciated so much!  
 A: The Liouville equation for the $N$ particle system, describes the time evolution of the phase space N-particle probability density, which you can also neatly rewrite with the Liouville operator: $f^{N}(t)= e^{-iLt}f^{N}(0).$ Now almost always we're interested in a smaller subset of only $n$ particles, for which then we have to define a reduced phase space distribution function, obtained by integrating out the remaining degrees of freedom ($6(N-n)).$ So the reduced $f^n$ has the form:
$$
f^n (\mathbf{r}^n,\mathbf{p}^n,t) = \frac{N!}{(N-n)!}\int f^N(\mathbf{r}^N,\mathbf{p}^N,t)d\mathbf{r}^{N-n}d\mathbf{p}^{N-n}
$$ 
For our purposes here, to derive an equation of motion for $f^n,$ let's consider a simpler general Hamiltonian:
$$
\frac{\partial H}{\partial \mathbf{r}_i}=-\mathbf{X}_i-\sum_{j=1}^{N}\mathbf{F}_{ij}
$$
With $\mathbf{X}$ denotes the external forces due to an external field and $\mathbf{F}_{ij}$ a pairwise interparticle interaction. Inserting into the Liouville equation, for $f^N,$ with a slight rearrangement:
$$
\left(\frac{\partial}{\partial t}+\sum_{i=1}^{N}\frac{\mathbf{p}_i}{m}\frac{\partial}{\partial\mathbf{r}_i}+\sum_{i=1}^N \mathbf{X}_i\frac{\partial}{\partial \mathbf{p}_i}\right)f^N=-\sum_{i=1}^{N}\sum_{j=1}^{N}\mathbf{F}_{ij}\frac{\partial f^N}{\partial \mathbf{p}_i}
$$
Now by integrating out the unwanted $6(N-n)$ degrees of freedome, we obtain:
$$
\left(\frac{\partial}{\partial t}+\sum_{i=1}^{n}\frac{\mathbf{p}_i}{m}\frac{\partial}{\partial\mathbf{r}_i}+\sum_{i=1}^n \mathbf{X}_i\frac{\partial}{\partial \mathbf{p}_i}\right)f^n=-\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbf{F}_{ij}\frac{\partial f^n}{\partial \mathbf{p}_i}- \underbrace{\frac{N!}{(N-n)!}\sum_{i=1}^n \sum_{j=n+1}^N \iint \mathbf{F}_{ij} \frac{\partial f^N}{\partial \mathbf{p}_i} d\mathbf{r}^{N-n}}_{(*)}
d\mathbf{p}^{N-n}
$$
Notice that all the terms with $i>n$ have vanished, to ensure a valid density function $f^n$ for the subspace. Working with identical particles here, we know $f^N$ is symmetric w.r.t. interchange of individual particle labels and the second sum in $(*)$ from $n+1$ to $N$ can be replaced by $N-n$ times the value of any one term (identical particles). With this simplification we can re-write the above equation as:
$$
\left(\frac{\partial}{\partial t}+\sum_{i=1}^{n}\frac{\mathbf{p}_i}{m}\frac{\partial}{\partial\mathbf{r}_i}+\sum_{i=1}^n \left(\mathbf{X}_i +\sum_{j=1}^n \mathbf{F}_{ij}\right)\frac{\partial}{\partial \mathbf{p}_i}\right)f^n=-\sum_{i=1}^n \iint \mathbf{F}_{i,n+1} \frac{\partial f^{n+1}}{\partial \mathbf{p}_i} d\mathbf{r}_{n+1}
d\mathbf{p}_{n+1}
$$
We have derived the BBGKY hierarchy, for a relatively less general Hamiltonian. But note that it is an exact expression as it stands, linking the one particle phase space density to the two particle density which itself linked the three particle...all the way up to the N particle phase space density. Although such expression doesn't lend itself to much use in practice (for one thing, we don't know any of the $f^n$'s), it can become useful by providing some form of a closure relation, in other words without a closure relation, all BBGKY does is: expresses one unknown density function $f^n$ in terms of another unknown $f^{n+1}.$ It does not however represent an equivalent formulation of the Liouville equation for $f^N$ when you set $n=N$, it only gives you a link with the density function of subspaces of N, which you still don't know. If no closure is provided, BBGKY is not really useful! Have a look at the book of JP Hansen and IR McDonald for a more in depth coverage. 
