You would expect the equilibrium distribution in the bulk to be isotropic, to not prefer a particular direction. Thus, likely not the direction of the relative velocity $\vec v=\vec \xi - \vec u$ (with all its three components) will enter the equilibrium distribution $f^{(eq)}$ but instead only the length of this vector, its magnitude $| \vec v|$. This corresponds to a normal distribution in the variable $\vec \xi$ (particle's individual velocity) around the mean $\vec u$ (macroscopic velocity) with the standard deviation corresponding to the isothermal speed of sound $c_s^T = \sqrt{ R_m T}$.
In case you are interested, here a brief non-rigorous derivation of the Maxwell-Boltzmann equilibrium distribution.
Just arguing with symmetries and conservation of moments, one is able to predict a possible distribution function in a state of equilibrium, the so called Maxwell-Boltzmann distribution for mono-atomic gases (simple gases with no internal vibrational degrees of freedom).
Decompositions with symmetries
In a system without boundary conditions with a certain initial perturbation it seems reasonable to assume that after some time the perturbation will even out across the domain and the system will reach an equilibrium distribution that is isotropic in velocity space. This means it should be possible to write an equilibrium distribution for velocity as a joint probability, as the product of three one-dimensional distribution functions in terms of the deviation of the microscopic velocities from the mean value, the macroscopic velocity $\vec v=\vec \xi - \vec u$:
$$ f^{(eq)}(\vec v^2) = f^{(eq)}(v_x^2+v_y^2+v_z^2) = f_{1D}^{(eq)}(v_x^2)\,f_{1D}^{(eq)}(v_y^2)\,f_{1D}^{(eq)}(v_z^2) $$
If the magnitude of the velocity deviation $\vec v^2$ is held constant also $f^{(eq)}(\vec v^2)$ must be a constant and therefore
$$ln[f^{(eq)}(\vec v^2)] = ln[f_{1D}^{(eq)}(v_x^2)]+ln[f_{1D}^{(eq)}(v_y^2)]+ln[f_{1D}^{(eq)}(v_z^2)] = const $$
holds. The simplest non-trivial distribution function fulfilling this criterion is given by $ln[f_{1D}^{(eq)}(v_x^2)] = A-B v_x^2$ due to
$$ ln[f_{1D}^{(eq)}(v_x^2)]+ln[f_{1D}^{(eq)}(v_y^2)]+ln[f_{1D}^{(eq)}(v_z^2)]=3A+B(v_x^2+v_y^2+v_z^2)=3A-B|\vec v|^2. $$
Therefore one possible equilibrium distribution function takes the form of
$$f^{(eq)}(|\vec v|) = e^{3A} e^{-B|\vec v|^2} = C e^{-B|\vec v|^2}.$$
Expected values and macroscopic properties
We have to determine now the constants by enforcing the correlation between expected values of the "microscopic" distribution and the macroscopic flow field.
The zero-th moment of any particle distribution must be equal to the number density $n := \frac{N}{V}$, which is connected to the density of the gas by the mass of a single particle $\rho = n m_p$. The resulting integral converges only for negative exponents $-B < 0$. Again we assume that the integral can be split up into one-dimensional solutions according to
$$\iiint \left( \frac{B}{\pi} \right)^{\frac{3}{2}} e^{-B \vec x^2} d \vec x = \left( \int \sqrt{\frac{B}{\pi}} e^{-B x^2} dx \right) \left( \int \sqrt{\frac{B}{\pi}} e^{-B y^2} dy \right) \left( \int \sqrt{\frac{B}{\pi}} e^{-B z^2} dz \right)$$
This integral can be solved by transformation to polar coordinates
$$ \int\limits_{x=-\infty}^{\infty} e^{-x^2} dx = \sqrt{ \left( \, \int\limits_{y=-\infty}^{\infty} e^{-y^2} dy \right) \left( \, \int\limits_{z=-\infty}^{\infty} e^{-z^2} dz \right) } = \sqrt{ \, \int\limits_{z=-\infty}^{\infty} \, \int\limits_{y=-\infty}^{\infty} e^{-(y^2+z^2)} dy \, dz }
= \sqrt{ \int\limits_{\phi=0}^{2 \pi} \, \int\limits_{r=0}^{\infty} r \, e^{-r^2} dr \, d \phi } = \sqrt{ \int\limits_{\phi=0}^{2 \pi} \left[ - \frac{e^{-r^2}}{2} \right]_{r=0}^{\infty} d \phi } = \sqrt{ \int\limits_{\phi=0}^{2 \pi} \frac{1}{2} d \phi } = \sqrt{ \left[ \frac{\phi}{2} \right]_{\phi=0}^{2 \pi} } = \sqrt{ \pi }
$$
Similarly deriving both sides of the equation
$$\int\limits_{x=-\infty}^{\infty} \sqrt{\frac{C}{\pi}} e^{-C x^2} dx = \left[ \frac{1}{2} erf{(\sqrt{C} x)} \right]_{x=-\infty}^{\infty} = 1$$
with respect to the parameter $C$
$$ \frac{\partial^n}{\partial C^n} \left( \, \int\limits_{x=-\infty}^{\infty} e^{-C x^2} dx \right) = \frac{\partial^n}{\partial C^n} \left( \sqrt{ \frac{\pi}{C} } \right) $$
allows us to determine integration formulas for the all even moments
$$ \int\limits_{x=-\infty}^{\infty} x^{2n} e^{-C x^2} dx = \frac{1 \cdot 3 \cdot 5 \cdot 7 \cdot \ldots \cdot (2n-1) \sqrt{\pi}}{2^n \, C^{n + \frac{1}{2}}} = \frac{(2n)! \sqrt{\pi}}{n! \, 2^{2n} C^{n + \frac{1}{2}}} $$
while due to symmetry the odd-order moments must vanish.
This leads to
$$ C = n \left( \frac{B}{\pi} \right) ^{\frac{3}{2}} $$
meaning this particular equilibrium distribution must take the form
$$ f^{(eq)}(|\vec v|) = n \left( \frac{B}{\pi} \right)^{\frac{3}{2}} e^{-B|\vec v|^2}. $$
Enforcing the second-order momentum for the internal energy
$$ \rho e = \frac{3 \, n \, k_B \, T}{2} = \frac{m_P}{2} \iiint \vec v^2 f^{(eq)}(|\vec v|) \, d \vec v. $$
yields with $\frac{k_B}{m_P} = \frac{R}{M} = R_m$
$$B = \frac{m_P}{2 \, k_B \, T} = \frac{1}{2 \, R_m \, T}.$$
This means the so called Maxwell-Boltzmann distribution in three-dimensional space takes the form of a Gaussian distribution
$$ f^{(eq)}(\vec x, |\vec v|, t) = n \frac{1}{(2 \pi R_m T )^{\frac{3}{2}}} \, e^{-\frac{|\vec v|^2}{2 R_m T}}. $$
Multivariate normal distribution
The equation for a $d$-dimensional Gaussian with the corresponding means $\vec \mu$ and the covariance matrix
$$\Sigma_{ij} = \mathbb{E}[(x_i - \mathbb{E}[x_i])(x_j - \mathbb{E}[x_j])]$$
is given by
$$ \mathcal{N} (\vec x, \vec \mu, \underline{\Sigma}) = \frac{1}{\sqrt{(2 \pi)^d |\underline{\Sigma}|}} e^{- \frac{1}{2} (\vec x - \vec \mu)^T \underline{\Sigma}^{-1} (\vec x - \vec \mu)}.$$
In case of uncorrelated variables and three-dimensional space the covariance matrix degenerates to
$$ \underline{\Sigma} = \begin{pmatrix}
\sigma_x^2 & 0 & 0\\
0 & \sigma_y^2 & 0 \\
0 & 0 & \sigma_z^2
\end{pmatrix}$$
and its inverse can be calculated to
$$ \underline{\Sigma}^{-1} = \begin{pmatrix}
\frac{1}{\sigma_x^2} & 0 & 0\\
0 & \frac{1}{\sigma_y^2} & 0 \\
0 & 0 & \frac{1}{\sigma_z^2}
\end{pmatrix}.$$
If we assume similarity between the three spatial directions
$$ \sigma^2 = \sigma_x^2 = \sigma_y^2 = \sigma_z^2$$
the determinant of the covariance matrix is given by
$$| \underline{\Sigma} | = \det{\underline{\Sigma}} = \sigma_x^2 \, \sigma_y^2 \, \sigma_z^2 = \sigma^6.$$
For this reason the Maxwell-Boltzmann equilibrium distribution can be seen as a three-dimensional ($d = 3$) Gaussian normal distribution of uncorrelated variables $\vec x = \vec \xi$, the particle's individual speeds, around the average macroscopic velocity $\vec \mu = \vec u$ with the standard deviation corresponding to the isothermal speed of sound $\sigma = \sqrt{ R_m T} = c_s^T$
$$ \frac{f^{(eq)}}{n} = \mathcal{N} (\vec x, \vec \mu, \sigma) = \frac{1}{(2 \pi)^\frac{d}{2} \sigma^d} e^{- \frac{1}{2} \left( \frac{\vec x - \vec \mu}{\sigma} \right)^2}.$$
