Why the minimum point for this Hamiltonian is not a maximum for the probability density function

$\newcommand{\diff}{\operatorname{d}}$ $\newcommand{\deriv}{\frac{\diff #1}{\diff #2}}$ $\newcommand{\pderiv}{\frac{\partial #1}{\partial #2}}$ $\renewcommand{\vec}{\boldsymbol{#1}}$I have a system of $N$ identical independent non-interacting particles with single mass $m$ forced to move on a plane with single particle Hamiltonian $$H^{(1)}(\vec{p}, \vec{q}) \equiv H(\vec{p}, \vec{q}) = \frac{\vec{p}^2}{2m} + A\vec{q}^4$$ with $\vec{p}, \vec{q}\in\mathbb{R}^2$. Now I want to find the most probable value of the quantity $q = |\vec{q}|$, say $q^*$, so the point is to find the ensemble density distribution $\rho(\vec{p}, \vec{q})$ for the system and maximise it in the variable $q$.

It's known that the canonical ensemble density function fon a single particle is given by $$\rho(\vec{p}, \vec{q}) = \frac{e^{-\beta H(\vec{p}, \vec{q})}}{Z_1},$$ where $\beta \equiv T^{-1}$ ($k_B=1$) and $Z_1$ is the single particle partition function defined by $$Z_1 = \iint\limits_{\vec{p},\vec{q}\in\mathbb{R}^2}\frac{\diff\vec{p}\diff\vec{q}}{h^2}e^{-\beta H(\vec{p}, \vec{q})}.$$ Integrating $\rho(\vec{p}, \vec{q})$ on $\diff\vec{p}$ and $\diff\theta_q$ the $\rho(q)$ is easily find to be (by marginalising) $$\rho(q) = \int_{\mathbb{R}^2}\frac{\diff\vec{p}}{h^2}\int_0^{2\pi}\diff\theta_q\rho(\vec{p}, \vec{q}) = \frac{\int_{\mathbb{R}^2}\diff\vec{p}\int_0^{2\pi}\diff\theta_q e^{-\beta H(\vec{p}, \vec{q})}}{\iint\limits_{\vec{p},\vec{q}\in\mathbb{R}^2}\diff\vec{p}\diff\vec{q}e^{-\beta H(\vec{p}, \vec{q})}} = \frac{\int_0^{2\pi}\diff\theta_q \cdot q\cdot e^{-A\beta\vec{q}^4}}{\int_0^{2\pi}\diff\theta_q \int_0^{\infty} \diff{q}\ q\cdot e^{-A\beta\vec{q}^4}} = \frac{4\sqrt{A\beta}}{\sqrt{\pi}}q\cdot e^{-A\beta\vec{q}^4}.$$ The maximum of $\rho(q)$ is then obtained by impose the first derivative equal to 0 $$\deriv{\rho(q)}{q} = 0 \implies q^* = \left(\frac{T}{4A}\right)^{\frac14}.$$ As you can see the maximum of the $q$ distribution is not $q* = 0$, but a value greater than zero, that if is plugged in the hamiltonian does not minimise the energy. Is there anything wrong in these reasoning? Because I naively thought that the lowest energy state should also be the most probable, but in this case the value found is bigger than zero, so represents an higher energy state.

And also how it is possible that if i put $q \equiv |\vec{q}| = 0$ inside the $\rho(\vec{q}) \equiv \int_{\mathbb{R}^2}\diff\vec{p}\rho(\vec{p}, \vec{q})$ I obtain the maximum (which is compatible with the above reasoning), while if I put it inside the $\rho(q)$ I got $0$? I guess it must be some misunderstanding with the change of coordinates (from cartesian ones to polar ones) but I can't figure it out.

• Indeed, it is to do with the change in variables. $\rho(q)$ and $\rho(\mathbf{q})$ are two different (but related) probability densities. Your variable $q$, for instance, can never be negative, so its average is bound to be positive. There's a fairly lengthy explanation here. A similar question comes up when one compares the most probable velocity with the most probable speed in the Maxwell-Boltzmann distribution.
– user197851
Aug 26 '18 at 11:02

Indeed, it is to do with the change in variables, and you have presented all the relevant mathematical expressions in your question. $\rho(q)$ and $\rho(\mathbf{q})$ are two different (but related) probability densities.
The maximum of $\rho(\mathbf{q})$ is at the minimum of the configurational part of $H$, $A|\mathbf{q}|^4=Aq^4$, i.e. at $\mathbf{q}=(q_x,q_y)=(0,0)$, where the magnitude $q$ is zero. We are discussing minima/maxima with respect to the two components of $\mathbf{q}$. So this is the most probable (vector) $\mathbf{q}$.
However, the maximum of $\rho(q)$, with respect to the single variable $q$, is at a nonzero value of $q$. Mathematically, there is a Jacobian involved in changing variables to polar coordinates, and an extra factor $q$ ends up in the expression for $\rho(q)$. You derived this in your question. Physically, points in the two-dimensional plane $(q_x,q_y)$ having values of $|\mathbf{q}|$ in the range $q\ldots q+dq$ lie in an annulus of area $2\pi \, q \, dq$, and this area increases as $q$ becomes larger, given constant $dq$. The probability density is "probability per unit $dq$".
Questions of this kind often arise when one meets the Maxwell-Boltzmann distribution of velocities, i.e. "why is the most probable speed $v=|\mathbf{v}|$ different from the most probable velocity $\mathbf{v}$?". It's exactly the same argument as in this case. For 3D systems, of course, the Jacobian introduces an extra weighting of $v^2$, since the relevant velocities lie in a spherical shell of volume $4\pi \, v^2 \,dv$.