I am trying to integrate over $v_y$ and $v_z$ the following function in Mathematica 13.1 \begin{equation} n(v_x, v_y, v_z)= \frac{1}{\exp\left[\frac{\frac{m}{2}\left(v_x^2 + v_y^2 + v_z^2 \right) - \mu}{k_B T}\right]+ 1}. \end{equation} All parameters are real. Moreover $T$, $\mu$, $k_B$, $m$ and $h$ are positive constants.
I know that the result must be something like \begin{equation} n(v_x) =\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}n(v_x, v_y, v_z)\,\mathrm{d}v_y\,\mathrm{d}v_z=\frac{2\pi k_B T}{m} \log\left[1 + \exp\left(-\frac{\frac{1}{2}m v_x^2 - \mu}{k_B T }\right)\right] \end{equation}
My code is
Clear["Global`*"]
$Assumptions =
Vx \[Element] Reals && Vy \[Element] Reals && Vz \[Element] Reals &&
kB > 0 && m > 0 && T > 0 && \[Mu] > 0; Integrate[
1/(Exp[(m/2*(Vx^2 + Vy^2 + Vz^2) - \[Mu])/(k*T)] +
1), {Vy, -\[Infinity], +\[Infinity]}, {Vz, -\[Infinity], +\
\[Infinity]}]
Mathematica 13.1 does not find the solution.
The thing that puzzles me is that a previous version of Mathematica gave me the correct result.