Rigorous derivation of the mean free path in a gas 
Can anyone supply me with a derivation of the mean free path, of
  particles in a Maxwell Boltzmann Gas?

Cited in various literature is the formula,
\begin{align}
\begin{split}
\ell&=\frac{1}{\sqrt{2}n\sigma},
\end{split}
\end{align}
which describes the mean free path of an atom/molecule in a Maxwell Boltzmann gas. Within it, $n$ is the gas density (assumed to be homogeneous) and $\sigma$ is the cross section of colliding particles.
Derivations of the mean free path tend to include a comment about the relative velocity of particles. However, most literature will not treat this as a distribution, and will instead follow a two part approximate method.
First, magnitude of the relative velocity is expressed using equation (1). Within it, $\theta$ is the angle between relative velocities $2$ and $1$.
\begin{align}
\tag{1}
|\mathbf{v}_{r}|&=\sqrt{\mathbf{v}_{r}\cdot\mathbf{v}_{r}}\\
&=\sqrt{(\mathbf{v_2}-\mathbf{v_1})\cdot(\mathbf{v_2}-\mathbf{v_1})}\\
&=\sqrt{\mathbf{v_1}^2+\mathbf{v_2}^2-2\mathbf{v_1}\cdot\mathbf{v_2}}\\
v_{r}&=\sqrt{v_1^2+v_2^2-2v_1v_2cos(\theta)}\\
\end{align}
Once equation (1) has been found, most authors take an average of both sides. Instead of evaluating this by treating each velocity as a distribution however, they instead take the average under the square root and onto each summed term (equation 2)
\begin{align}
\tag{2}
\langle v_{r}\rangle&=\big\langle\sqrt{v_1^2+v_2^2-2v_1v_2cos(\theta)}\big\rangle\\
&=\sqrt{\big\langle v_1^2+v_2^2-2v_1v_2cos(\theta)\big\rangle}\\
&=\sqrt{\big\langle v_1^2\big\rangle+\big\langle v_2^2\big\rangle-\big\langle2v_1v_2cos(\theta)\big\rangle}\\
&=\sqrt{\langle v_1\rangle^2+\langle v_2\rangle^2-0}\\
&=\sqrt{\langle v_1\rangle^2+\langle v_2\rangle^2}\\
&=\sqrt{2\langle v_1\rangle^2}\\
&=\sqrt{2}~\langle v_1\rangle\\
\end{align}
A large amount of (2) is wrong. Exceptions are the forth and fifth steps, where an average of the dot product between two arbitrary velocities is zero, justified by the velocity in a Maxwell Boltzmann distribution having no preferred direction. The remaining step involves the average velocity of any particle in a Boltzmann gas being the same ($\langle v_1\rangle=\langle v_2\rangle$).
One might think to replace the average of the square of velocity with an rms velocity, however doing so yields the rms relative velocity distribution, which is out by $8\%$ from the mean relative velocity.
I would like to see a derivation of the same result without these approximate methods. Failing that, reference to a textbook that provides them would be helpful as well. 
References (all use the approximation):


*

*Reif, F. (1965) Fundamentals of Statistical and Thermal Physics. McGraw Hill, New York, 273-278.

*formula for mean free path in two dimensions

*Why is the mean free path divided by $\sqrt{2}$?

*RMS Free Path vs Mean Free Path
 A: Another way:
$$\langle v_r \rangle = \mathcal{N}^2 \iint \vert \vec{v}_1 - \vec{v}_2\vert \ e^{-\frac{m(\vec{v}_1^2+\vec{v}_2^2)}{2k_B T}}\ \mathrm{d}\vec{v}_1  \mathrm{d}\vec{v}_2,$$
wherein the normalized Boltzmann distribution function is $\mathcal{N} e^{-\frac{m \vec{v}^2}{2k_B T}}$, $\ \mathcal{N}$ being the numerical factor that normalizes the Gaussian distribution.
Rotating velocity coordinates by $45^{\circ}$, define new coordinates
$$\vec{U}= \frac{\vec{v}_1+\vec{v}_2}{\sqrt{2}},$$
$$\vec{V}= \frac{-\vec{v}_1+\vec{v}_2}{\sqrt{2}}.$$
Since the above transformation is a rotation, the Jacobian is 1.
Rewriting $\langle v_r \rangle$ in terms of new coordinates,
\begin{align}
\langle v_r \rangle &= \mathcal{N}^2 \iint \vert \sqrt{2} \vec{V}\vert e^{-\frac{m(\vec{V}^2+\vec{U}^2)}{2k_B T}}\ \mathrm{d}\vec{V} \mathrm{d}\vec{U} \\
&=  \sqrt{2} \ \underbrace{\int  \vert \vec{V}\vert \ \mathcal{N} e^{-\frac{m\vec{V}^2}{2k_B T}}\ \mathrm{d}\vec{V}}_{\text{equals } \langle v \rangle} \ \ \underbrace{\int \mathcal{N} e^{-\frac{m\vec{U}^2}{2k_B T}} \ \mathrm{d}\vec{U}}_{\text{equals }1}  \\
& = \sqrt{2} \ \langle v \rangle. 
\end{align}
The relation, $\langle v_r \rangle = \sqrt{2} \ \langle v \rangle $, seems to hold independent of the spatial dimensionality.
A: Edit: This answer is wrong.  I made an error that is equivalent to implicitly assuming that the particles' velocities are always perpendicular (see below).  The question author has asked me to leave it up as it is interesting that the final result is so nearly correct despite the mistake.
The Maxwell-Boltzmann distribution is
$$B(\mathbf{v}) = \left( \frac{m}{2\pi k T} \right)^{3/2} e^{-\frac{m|\mathbf{v}|^2}{2kT}}$$
This function has a few nice features that make the following calculations easier.  First, it's normalized, so integrals such as $\int_{\mathbb{R}^3} B(\mathbf{v}) \, d^3\mathbf{v}$ will always evaluate to 1.  Second, since it is a radial function, we'll be able to take advantage of the trick
$$\int_{\mathbb{R}^n} f(|\mathbf{x}|) \, d^n\mathbf{x} = \omega_{n-1} \int_0^{\infty}f(r) r^{n-1} \,dr$$
where $\omega_{n-1}$ is the area of an $n$-sphere of radius 1.
All right, there's no way to fully answer your question without doing a few ugly integrals, so let's dive in.  The mean speed of a particle is
\begin{align}
\langle v \rangle &= \frac{\int_{\mathbb{R}^3} v \, B(\mathbf{v}) \, d^3\mathbf{v}}{\int_{\mathbb{R}^3} B(\mathbf{v}) \, d^3\mathbf{v}} \\
&= \int_{\mathbb{R}^3} v \, B(\mathbf{v}) \, d^3\mathbf{v} \\
&= \left( \frac{m}{2\pi k T} \right)^{3/2} \int_{\mathbb{R}^3} v \, e^{-\frac{m|\mathbf{v}|^2}{2kT}} \, d^3\mathbf{v} \\
&= \left( \frac{m}{2\pi k T} \right)^{3/2} \int_0^{\infty} v \, e^{-\frac{mv^2}{2kT}} \, 4\pi v^2 dv \\
&= \left( \frac{m}{2kT} \right)^{3/2} \frac{4}{\sqrt{\pi}} \int_0^{\infty} v^3 \, e^{-\frac{mv^2}{2kT}} \, dv \\
&= \left( \frac{m}{2kT} \right)^{3/2} 2\pi \left( \frac{m}{2kT} \right)^{-2} \\
&= \left( \frac{m}{2kT} \right)^{3/2} \frac{2}{\sqrt{\pi}} \left( \frac{m}{2kT} \right)^{-2} \\
&= \frac{2}{\sqrt{\pi}} \left( \frac{2kT}{m} \right)^{1/2}
\end{align}
[Edit: if we make the approximation that the particles' velocities are perpendicular,] the mean relative speed of two particles is
\begin{align}
\langle v_r \rangle &\approx \frac{\int_{\mathbb{R}^3} \int_{\mathbb{R}^3} v_r \, B(\mathbf{v})B(\mathbf{v}') \, d^3\mathbf{v} \, d^3\mathbf{v}'}{\int_{\mathbb{R}^3}\int_{\mathbb{R}^3} B(\mathbf{v})B(\mathbf{v}') \, d^3\mathbf{v} \, d^3\mathbf{v}'} \\
&= \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} v_r \, B(\mathbf{v})B(\mathbf{v}') \, d^3\mathbf{v} \, d^3\mathbf{v}' \\
&= \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} v_r \, \left( \frac{m}{2\pi k T} \right)^{3/2} e^{-\frac{m|\mathbf{v}|^2}{2kT}} \left( \frac{m}{2\pi k T} \right)^{3/2} e^{-\frac{m|\mathbf{v'}|^2}{2kT}} \, d^3\mathbf{v} \, d^3\mathbf{v}' \\
&= \left( \frac{m}{2\pi k T} \right)^3 \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} v_r \,  e^{-\left(\frac{m|\mathbf{v}|^2}{2kT} + \frac{m|\mathbf{v'}|^2}{2kT}\right)} \, d^3\mathbf{v} \, d^3\mathbf{v}' \\
&= \left( \frac{m}{2\pi k T} \right)^3 \int_{\mathbb{R}^3}\int_{\mathbb{R}^3} v_r \, e^{-\frac{mv_r^2}{2kT}} \, d^3\mathbf{v} \, d^3\mathbf{v}' \\
&= \left( \frac{m}{2\pi k T} \right)^3 \int_0^{\infty} v_r \, e^{-\frac{mv_r^2}{2kT}} \, \pi^3 v_r^5 \, dv_r \\
&= \left( \frac{m}{2 kT} \right)^3 \int_0^{\infty} v_r^6 \, e^{-\frac{mv_r^2}{2kT}} \, dv_r \\
&= \left( \frac{m}{2 kT} \right)^3 \frac{15}{16}\sqrt{\pi} \left( \frac{m}{2kT} \right)^{-7/2} \\
&= \frac{15}{16}\sqrt{\pi} \, \left( \frac{2kT}{m} \right)^{1/2}
\end{align}
where on line 6, $\pi^3 v_r^5$ is the surface area of a 6-sphere of radius $v_r$, and the integral on line 8 was calculated using Wolfram Alpha.
Thus we find that $\langle v_r \rangle / \langle v \rangle = 15\pi/32 \approx 1.4726$ which is almost, but not exactly equal to the textbook result $\sqrt{2} \approx 1.4142$.
[Edit: the textbook result $\langle v_r \rangle / \langle v \rangle = \sqrt{2}$ is correct.  For the full, correct proof without any simplifying assumptions, please see the other answer.]
