Integral of the Fermi function

Question

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.

Answer 1

It seems that in 13.2, MMA will do the first integral wrt Vy and returns a polylog function, then get stuck there.

We can use the new function IntegrateChangeVariables in 13.1,

integral=Inactive[Integrate][1/(Exp[x^2+y^2+z^2-1]+1),{y,-Infinity,Infinity},{z,-Infinity,Infinity}]

to transform the integral into the radial coordinates,

IntegrateChangeVariables[integral,{r,θ},{y==r Cos[θ],z==r Sin[θ]},Assumptions->0<θ<2π&&r>0]//Activate
(*π Log[1+E^(1-x^2)]*)