Introduction: thermodynamics and Landau argument
We introduce the main ideas of the Fermi liquid and motivate many-body perturbation theory. A reference for the Fermi liquid ideas is the first chapter of Landau and Lifshitz volume 9.
The assumption of Fermi liquid theory is that, whatever the interactions may be, we can identify "elementary excitations" that are like particles (and hence are called "quasi-particles"): they have well-defined momentum, spin-$\frac{1}{2}$, and charge $e$. Furthermore, they are spatially "local" relative to macroscopic sizes like a sample size, and are long-lived if their excitation energy is small. This last point is quite important: quasiparticles of finite excitation energy $E>0$ do interact with each other and can decay, but the lifetime becomes larger as $E \rightarrow 0$.
As an intuitive picture of the quasiparticles, think of them as a single electron dressed by a cloud of electron-hole pairs: as a result, the effective mass may be modified by the "screening cloud". But, since the electron-hole pairs are neutral and bosonic (actually we'll assume them usually to be spin-zero), the charge and fermionic statistics are unmodified.
The Fermi momentum is determined through
$$\begin{equation}
N/V = {2 \cdot 4 \pi p_F^3 / 3 (2 \pi \hbar)^3}.
\end{equation}$$
We will assume that this continues to hold because of the adiabatic connection between quasiparticles and original fermions. Define the normalized quasiparticle distribution through
$$\begin{equation}
\sum_\alpha \int n_{\alpha \alpha}\,d\tau = {\rm tr}\ \int {\hat n}\,d\tau = {N \over V}, d\tau = {d^3 p \over (2 \pi \hbar)^2}.
\end{equation}$$
In the following we will use the Einstein summation convention that repeated indices are implicity summed over. The meaning of the indices on the Hermitian density matrix ${\hat n}$ is that diagonal elements correspond to the density of spin-up or spin-down electrons.
In case you haven't seen a density matrix in a while, here is a lightning review. Recall that the off-diagonal term reflects the fact that spin is a quantum variable. For instance, for a single particle in a mixed state with half spin-up and half spin-down, we would have
$$\begin{equation}
n=\left(\matrix{\frac{1}{2} & 0\cr 0&\frac{1}{2}}\right)
\end{equation}$$
while for a pure state with spin aligned along some axis in the $x$-$y$ plane,
$$\begin{equation}
n=\left(\matrix{\frac{1}{2} & \frac{1}{2}\cr \frac{1}{2}&\frac{1}{2}}\right).
\end{equation}$$
The spin operator along the $x$-axis is given by $\sigma_x$ in this two-component space, so
$$\begin{equation}
\langle s_x \rangle = {\rm tr}\ n \sigma_x
\end{equation}$$
where here this is a quantum statistical expectation value, which for a pure state coincides with the ordinary quantum expectation value. Pure states are essentially projection operators, so $n^2 = n$; any state, pure or mixed, has a trace related to the total number of particles.
The change in energy due to a change $\delta n$ in the quasiparticle occupancies should be written, for a small change from the equilibrium distribution obtained below, as
$$\begin{equation}
{\delta E \over V} = \int \epsilon({\bf p}) \delta n d\tau.
\end{equation}$$
Actually, let us make a slightly more general form to account for situations where the distribution of up-spin particles is different from that of down-spin particles.
$$\begin{equation}
{\delta E \over V} = \int \epsilon_{\alpha \beta} ({\bf p}) \delta n_{\beta \alpha} d\tau =
{\rm tr}\ \int {\hat \epsilon}({\bf p}) \delta {\hat n} d\tau.
\label{echange}
\end{equation}$$
For the spin-symmetric case, the tensors ${\hat \epsilon}$ and $\delta {\hat n}$ are both diagonal: $\epsilon_{\alpha \beta} = \epsilon \delta_{\alpha \beta}$, $n_{\alpha \beta} = n \delta_{\alpha \beta}$.
The assumption of Fermi liquid theory can now be stated: we will allow the effective energy ${\hat \epsilon}({\bf p})$ in equation (\ref{echange}) to depend on the occupancy of other quasiparticle states, in a simple way. Denote by $\delta \epsilon_{\alpha \beta}({\bf p})$ the change in the effective energy of quasiparticles of momentum ${\bf p}$ induced by a given nonequilibrium distribution of the other quasiparticles:
$$\begin{equation}
\delta \epsilon_{\alpha \beta}({\bf p}) = \int f_{\alpha \gamma, \beta \delta}({\bf p},{\bf p}^\prime) \delta n({\bf p^\prime})_{\gamma \delta}\,d\tau^\prime.
\end{equation}$$
Here $f$ is some effective interaction about which we will have much more to say. The important thing is that we kept only the linear term $\delta n$ on the right-hand side and neglected higher powers in the deviation from equilibrium. So, even though the interaction $f$ in the above may be quite strong, any three-body interactions or other interactions not of the above form are neglected.
Now let me explain why the equilibrium distribution of quasiparticles has formally the same form as in the noninteracting case, even though the quasiparticles are interacting and the meaning of this distribution is quite different. The simplest derivation of the equilibrium distribution comes from maximizing the entropy
$$\begin{equation}
{S \over V} = - {\rm tr}\ \int \left( {\hat n} \log {\hat n} + (1- {\hat n}) \log(1 - {\hat n})\right)d\tau
\end{equation}$$
subject to the constraints of constant total particle number and energy,
$$\begin{equation}
{\delta N \over V} = {\rm tr}\ \int \delta {\hat n}\,d\tau = 0,\quad
{\delta E \over V} = {\rm tr}\ \int {\hat \epsilon} \delta {\hat n}\,d\tau = 0.
\end{equation}$$
Note that this expression for the entropy depends only on the fermionic statistics of the quasiparticles.
Introducing Lagrange multipliers $\lambda_1$ and $\lambda_2$, the total variation is
$$\begin{equation}
\log {\hat n \over 1 - \hat n} - \lambda_1 - \lambda_2 {\hat \epsilon} = 0.
\end{equation}$$
We leave it to the reader to show that $\lambda_1 = \beta \mu$ and $\lambda_2 = -\beta$, so a solution of the above is
$$\begin{equation}
{\hat n} = {1 \over e^{\beta({\hat \epsilon} - \mu)} + 1}
\end{equation}$$
since this gives
$$\begin{equation}
\log e^{\beta(\mu - {\hat \epsilon})} - \beta \mu + \beta {\hat \epsilon} = 0.
\end{equation}$$
For spin-symmetric quasiparticle energies, this becomes a single-component equation
$$\begin{equation}
n = {1 \over e^{\beta(\epsilon - \mu)} + 1}
\end{equation}$$
This, of course, looks the same as what we had for free particles. However, note that now it is actually a complicated self-consistency equation, since $\epsilon$ on the right side is effectively a function of $n$.
There is a simple self-consistent calculation that partly illustrates why the Fermi liquid is so stable to repulsive interactions in three dimensions. Consider two excitations of small energy above the Fermi level, and ask what is the cross-section for their scattering. Suppose the two initial particles have total energy $E$ and total momentum ${\bf p}$.
The scattering linewidth, from Fermi's golden rule, is
$$\begin{equation}
\Gamma = {2 \pi \over \hbar} \int |M|^2 \delta(E_1 + E_2 - E) \delta({\bf p}_1 + {\bf p}_2 - {\bf p}).
\end{equation}$$
In the above the integral is over all final states, and the matrix element $M$ is something dependent on the details of the interparticle interaction. Suppose that $M$'s variation is unimportant, and just look at the available phase space for scattering as $E$ is small. At zero temperature, the only final states available are outside the Fermi sea. Noting that there are six integrals and four constraints, and linearizing around the Fermi surface as usual: $E_k - \mu = v_F (k-k_F)$, we guess that $\Gamma \sim (E-\mu)^2$, which is borne out by a detailed calculation.
Hence in 3D $\Gamma \ll E-\mu$, so the lifetime of quasiparticles increases rapidly as the
quasiparticle moves toward the Fermi level. Another way to look at this, which will be used later, is that the linewidth of the quasiparticle becomes very narrow as $E \rightarrow \mu$. This gives some support for the contention that the Fermi liquid is a stable picture in 3D.