The following is a self answer concerning the rigorous derivation of the relative velocity in a Maxwellian gas. It is here for comparison against other answers.
Assumptions are
- Two particle collisions are the most probable ones, and collisions involving more than this, contribute to the mean relative velocity by a negligible amount.
- The mean relative velocity is strictly positive.
I think that this derivation differs from the one by @Thorondor because within it, $\mathbf{v_1}^2+\mathbf{v_2}^2\neq\mathbf{v_r}^2$. Which is to say that the dot product of $\mathbf{v_1}$, and $\mathbf{v_2}$ is not assumed to be zero.
This forces us to integrate over all possible angles between the velocities, which may account for the ~0.07 difference. It is also why the derivation is significantly longer.
The mean free path of an atom/molecule in a Maxwellian gas, depends upon the average relative velocity of each particle to one another. In order to obtain it, we first find the magnitude of velocity $\mathbf{v}_1$ relative to all other particles moving with $\mathbf{v}_2$. This expression is then averaged for all values of $\mathbf{v}_2$, from zero to infinity, producing a mean relative speed given $\mathbf{v}_1$ exists.
Following modus ponens, the mean is multiplied by the probability that an atom/molecule has velocity $\mathbf{v}_1$, and then averaged for all $\mathbf{v}_1$ from zero to infinity.
Consider two velocities $\mathbf{v}_1$ and $\mathbf{v}_2$ inclined at an angle $\theta$, there relative velocity will be described by
\begin{align}
(v_r)_{12}&=\sqrt{v_1^2+v_2^2-2v_1v_2cos(\theta)}.
\end{align}
All directions are equally probable for $\mathbf{v}_2$. To calculate the average of $(v_r)_{12}$, we multiply it by the probability that it lies within some solid angle $d\Omega$.
\begin{align}
\langle(v_r)_{12}\rangle&=\int_{\Omega}\frac{d\Omega}{4\pi}(v_r)_{12}\\
&=\int_{\Omega}\frac{2\pi sin(\theta)d\theta}{4\pi}(v_r)_{12}\\
&=\frac{1}{2}\int_{0}^{\pi}sin(\theta)d\theta(v_1^2+v_2^2-2v_1v_2cos(\theta))^\frac{1}{2}
\end{align}
Consider a change of basis, where $cos(\theta)=x$, so $sin(\theta)d\theta=-dx$. The limits of integration will change to $cos(0)=1$ and $cos(\pi)=-1$.
\begin{align}
\therefore~\langle(v_r)_{12}\rangle&=\frac{1}{2}\int_{1}^{-1}dx(v_1^2+v_2^2-2v_1v_2x)^\frac{1}{2}\\
&=\frac{1}{3}\frac{(v_1^2+v_2^2-2v_1v_2x)^\frac{3}{2}}{2v_1v_2}\Biggr|_{1}^{-1}\\
&=\frac{1}{6v_1v_2}\left((v_1^2+v_2^2+2v_1v_2)^\frac{3}{2}-(v_1^2+v_2^2-2v_1v_2)^\frac{3}{2}\right)\\
&=\frac{1}{6v_1v_2}\left((v_1+v_2)^\frac{3}{2}(v_1+v_2)^\frac{3}{2}+(v_1-v_2)^\frac{3}{2}(v_1-v_2)^\frac{3}{2}\right)\\
&=\frac{1}{6v_1v_2}\left((v_1+v_2)^3+|v_1-v_2|^3\right)\\
\end{align}
We take the magnitude of $v_1-v_2$ because $\langle(v_r)_{12}\rangle$ should always be positive, and therefore $\left((v_1-v_2)^2\right)^\frac{3}{2}$ must also be. This splits the average into two parts.
\begin{align}
\langle(v_r)_{12}\rangle&=\frac{1}{6v_1v_2}\left((v_1+v_2)^3-(v_1-v_2)^3\right),~when~v_1\geq v_2,\\
&=\frac{1}{6v_1v_2}\left[2v_2^3+6v_1^2v_2\right]\\
&=\frac{3v_1^2+v_2^2}{3v_1}.\\
\langle(v_r)_{12}\rangle&=\frac{1}{6v_1v_2}\left((v_1+v_2)^3+(v_1-v_2)^3\right),~when~v_2>v_1,\\
&=\frac{3v_2^2+v_1^2}{3v_2}.
\end{align}
The average relative velocity of an atom/molecule possessing magnitude and direction $\mathbf{v_1}$, with respect to all other particles moving with $\mathbf{v_2}$, lying within speeds of zero to infinity is then
\begin{align}
\langle(v_r)_1\rangle&=\int_{0}^{\infty}P(v_2)\langle(v_r)_{12}\rangle dv_2,\\
P(v_2)dv_2&=4\pi\left(\frac{m}{2\pi k_BT}\right)^\frac{3}{2}e^{-\frac{mv_2^2}{2k_BT}}v_2^2dv_2.
\end{align}
Note that $\langle(v_r)_1\rangle$ has two distinct forms depending on the difference between speeds. Because of this we break the integral into two parts.
\begin{align}
\langle(v_r)_1\rangle=4\pi\left(\frac{m}{2\pi k_BT}\right)^\frac{3}{2}&\left[\int_{0}^{v_1}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2\right.\\
&+\left.\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_2^2+v_1^2}{3v_2}\right)v_2^2\right]
\end{align}
To arrive at the relative velocity of any atom with another, we form a product of the above, and the probability a particle has velocity $v_1+dv_1$.
\begin{align}
\langle v_r \rangle&=\int_{0}^{\infty}P(v_1)\langle(v_r)_1\rangle dv_1
\end{align}
\begin{align}
\langle v_r \rangle=\left(4\pi\left[\frac{m}{2\pi k_BT}\right]^\frac{3}{2}\right)^2\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2&\left[\int_{0}^{v_1}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2\right.\\
&+\left.\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_2^2+v_1^2}{3v_2}\right)v_2^2\right]
\end{align}
which may be re-written as
\begin{align}
\langle v_r\rangle=\left(\frac{2m}{k_BT}\right)^3&\left[\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\int_{0}^{v_1}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2\right.\\
&+\left.\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_2^2+v_1^2}{3v_2}\right)v_2^2\right],
\end{align}
or
\begin{align}
\langle v_r\rangle&=\left(\frac{2m}{k_BT}\right)^3\left[I_1+I_2\right],
\end{align}
where $I_1$ stands for the first integral in the square brackets of the above, and $I_2$ the second one.
\begin{align}
I_1&=\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\int_{0}^{v_1}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2\\
I_2&=\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_2^2+v_1^2}{3v_2}\right)v_2^2
\end{align}
We wish to show that these two definite integrals are equivalent. To do so, a substitution, then change in the order of integration following Fubini's theorem is employed.
Consider the piecewise function,
\begin{align}
\mathbb{I}_{\{v_2\leq v_1\}}&=\begin{cases}
1&v_1\leq v_2\\
0&v_1>v_2\\
\end{cases}
\end{align}
which is zero outside the limit of the second integral of $I_1$. Using it, we can rewrite equation $I_1$ as
\begin{align}
I_1&=\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\int_{0}^{\infty}dv_2\mathbb{I}_{\{v_2\leq v_1\}}e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2.\\
\end{align}
We can combine these two integrals with a substitution $a=v_1$, allowing us to write the following
\begin{align}
I_1&=\int_{0}^{\infty}\int_{0}^{\infty}\mathbb{I}_{\{v_2\leq v_1\}}e^{-\frac{ma^2}{2k_BT}}a^2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2dv_2da.\\
\end{align}
Since the integrals limits are constant, we can apply Fubini's theorem exchange there order.
\begin{align}
I_1&=\int_{0}^{\infty}\int_{0}^{\infty}\mathbb{I}_{\{v_2\leq v_1\}}e^{-\frac{ma^2}{2k_BT}}a^2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_2^2dadv_2\\
&=\int_{0}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}v_2^2\int_{0}^{\infty}dv_1\mathbb{I}_{\{v_2\leq v_1\}}e^{-\frac{mv_1^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_1^2\\
&=\int_{0}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}v_2^2\int_{v_2}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}\left(\frac{3v_1^2+v_2^2}{3v_1}\right)v_1^2\\
\end{align}
Then, because the integral is definite, we can interchange $v_1$ with $v_2$, and obtain an equivalent volume.
\begin{align}
I_1&=\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_2^2+v_1^2}{3v_2}\right)v_2^2\\
&=I_2
\end{align}
Hence, the average relative velocity can be rewritten as
\begin{align}
\langle v_r\rangle&=2\left(\frac{2m}{k_BT}\right)^3I_2.
\end{align}
To evaluate $I_2$, first consider the integral
\begin{align}
I_0&=\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}\left(\frac{3v_2^2+v_1^2}{3v_2}\right)v_2^2.
\end{align}
We can break this into parts,
\begin{align}
I_0&=\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}v_2^3+\frac{v_1^2}{3}\int_{v_1}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}v_2,
\end{align}
then consider the substitution $\frac{m}{2k_BT}v_2^2=x$, so $dv_2=\frac{k_BT}{mv_2}$, and the lower limit of integration changes to $\frac{mv_1^2}{2k_BT}$.
\begin{align}
I_0&=\int_{\frac{mv_1^2}{2k_BT}}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}v_2^3+\frac{v_1^2}{3}\int_{\frac{mv_1^2}{2k_BT}}^{\infty}dv_2e^{-\frac{mv_2^2}{2k_BT}}v_2\\
&=2\left(\frac{k_BT}{m}\right)^2\int_{\frac{mv_1^2}{2k_BT}}^{\infty}e^{-x}xdx+\frac{v_1^2}{3}\frac{k_BT}{m}\int_{\frac{mv_1^2}{2k_BT}}^{\infty}e^{-x}dx
\end{align}
Since
\begin{align}
\int_{a}^{b}xe^{-x}dx&=-(x+1)e^{-x}\Biggr|_{a}^{b},
\end{align}
$I_0$ simplifies to
\begin{align}
I_0&=2\left(\frac{k_BT}{m}\right)^2\left[\frac{m}{2k_BT}v_1^2+1\right]e^{\frac{-mv_1^2}{2k_BT}}\\
&=e^{-\frac{mv_1^2}{2k_BT}}\left[\frac{4}{3}\left(\frac{k_BT}{m}\right)v_1^2+2\left(\frac{k_BT}{m}\right)^2\right]
\end{align}
Substituting this into $I_2$ yields
\begin{align}
I_2=&\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\left[\frac{4}{3}\left(\frac{k_BT}{m}\right)v_1^2+2\left(\frac{k_BT}{m}\right)^2\right]
\\=&2\left(\frac{k_BT}{m}\right)^2\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^2\\
&+\frac{4}{3}\frac{k_BT}{m}\int_{0}^{\infty}dv_1e^{-\frac{mv_1^2}{2k_BT}}v_1^4\\
\end{align}
$I_2$ is actually a standard integral of the form,
\begin{align}
\int_{0}^{\infty}e^{-ax^2}x^n&=\frac{1}{2a^{\frac{n+1}{2}}}\Gamma\left(\frac{n+1}{2}\right),
\end{align}
where $\Gamma$ is the gamma function, a generalization of the factorial. $I_2$ now simplifies to
\begin{align}
I_2&=\left(\frac{k_BT}{m}\right)^{\frac{7}{2}}\left[\Gamma\left(\frac{3}{2}\right)+\frac{2}{3}\Gamma\left(\frac{5}{2}\right)\right]\\
&=\left(\frac{k_BT}{m}\right)^{\frac{7}{2}}\left[\frac{\sqrt{\pi}}{2}+\frac{2}{3}\frac{3\sqrt{\pi}}{4}\right]\\
&=\left(\frac{k_BT}{m}\right)^{\frac{7}{2}}\sqrt{\pi}
\end{align}
Hence, the average relative speed of any Maxwellian gas particle with respect to any other is
\begin{align}
\langle v_r\rangle&=2\left[4\pi\left(\frac{m}{2\pi k_BT}\right)^\frac{3}{2}\right]^2\left(\frac{k_BT}{m}\right)^\frac{7}{2}\sqrt{\pi},\\
&=\frac{2\sqrt{\pi}\cdot16\pi^2}{8\pi^3}\left(\frac{k_BT}{m}\right)^{-3}\left(\frac{k_BT}{m}\right)^\frac{7}{2},\\
&=\frac{4}{\sqrt{\pi}}\sqrt{\frac{k_BT}{m}},\\
&=\sqrt{\frac{16}{\pi}\frac{k_BT}{m}}.
\end{align}
And the ratio of mean velocity $\langle v_r\rangle$ to mean relative velocity $\langle v\rangle$ is
\begin{align}
\frac{\langle v\rangle}{\langle v_r\rangle}&=\frac{\sqrt{\frac{8}{\pi}\frac{k_BT}{m}}}{\sqrt{\frac{16}{\pi}\frac{k_BT}{m}}}\\
&=\frac{1}{\sqrt{2}}
\end{align}
And thus the mean free path $\ell$ is
\begin{align}
\ell=\frac{1}{n\sigma\sqrt{2}}
\end{align}