Q:

Bohr-van Leeuwen theorem and quantum mechanics

Preamble:
If one considers an ideal gas of non interacting charged particles of charge $q$ in a uniform magnetic field $\mathbf{B} = \mathbf{\nabla} \wedge \mathbf{A}$, then the classical partition function in the canonical ensemble reads (in SI units):
$Q(\beta,V,N,\mathbf{B}) = \frac{1}{N!}q(\beta,V,\mathbf{B})^N$
where $q(\beta,V,\mathbf{B}) = \int \frac{d\mathbf{p} d \mathbf{r}}{h^3}\:e^{-\frac{\beta}{2m}(\mathbf{p}-q\mathbf{A}(\mathbf{r}))^2}$
If we integrate first with respect to momenta over all possible values from $-\infty$ to $+\infty$ for each component, a simple change of variable leads to 
$q(\beta,V,\mathbf{B})=\frac{V}{\Lambda^3}$ which is the ideal gas result and where $\Lambda$ is the thermal de Brooglie wavelength.
If one then wants to get the magnetization per particle $\mathbf{\mu}$ induced by the field $\mathbf{B}$ it is simply:
$\mathbf{\mu} = -\frac{\partial \langle \epsilon \rangle}{\partial \mathbf{B}} = \frac{\partial }{\partial \mathbf{B}}\left( \frac{\partial \ln(q(\beta,V,\mathbf{B}))}{\partial \beta} \right) = \frac{\partial }{\partial \beta}\left( \frac{\partial \ln(q(\beta,V,\mathbf{B}))}{\partial \mathbf{B}} \right) = \mathbf{0}$
This is one way to state the Bohr-van Leeuwen theorem.
Now, I physically understand this result as coming from some symmetry associated with the momenta (it is as likely to go to the right as it is to go to the left) and the fact that the boundaries of the integral over the momenta are infinite.
If the problem is treated quantum mechanically, the eigenstates of one charge particle are discretized Landau levels with a typical spacing between two neighbouring levels that is $\hbar \omega_c$ where $\omega_c = qB/m$ is the cyclotron frequency and one finds that the sum over these states depends on the magnetic field $\mathbf{B}$.
Question(s):
I am lost in my interpretation of the quantum to classical limit for this system...so far I thought that the quantum -> classical limit for the statistical properties of an individual particle was related to the way of counting the number of states for this particle i.e. whether we consider the set of states as a continuum or as a discrete set. This analogy seems to work in this case as well since the classical limit arises if $k_B T \gg \hbar \omega_c$. However two major points differ from what I am used to:

The quantum treatment of this system yields a non zero magnetic moment (although it vanishes at infinite temperatures) in the limit where $k_B T \gg \hbar \omega_c$ while the classical treatment gives strictly zero.
I do not understand how does the left-right symmetry argument used in the classical partition function disappear in the quantum treatment to yield a partition function that depends on $\mathbf{B}$.
Is there any classical way to assess that quantum corrections will be of order $\mathcal{O}(\Lambda/R_c)$ where $R_c \sim \sqrt{m k_B T}/(qB)$ is the typical size of the radius of the helical paths taken by a charged particle?

Sorry if my questions seem confused, I will try to improve them if they are not clear enough.
EDIT: I realize that one of my points is not very clear and shall explain it with the example of a true harmonic oscillator.
If I consider classical statistical mechanics, I know that $\langle \frac{1}{2}m\omega^2 x^2 \rangle = \frac{1}{2}k_B T$. This tells me that the typical uncertainty on the position of my particle is $\sigma_x = \sqrt{k_B T/(m\omega^2)}$. Incindently this length is also the typical confinement length scale owing to the harmonic potential. One way to semi-classicaly probe the validity of the classical limit is to imagine the particle as a non dispersive wave packet of width $\Lambda = h/\sqrt{2\pi m k_B T}$ and to realize that interferences (ultimately leading to quantization) are not important if $\Lambda \ll \sigma_x$.
This is very appealing because one can then probe the validity of a classical approximation using a $\sigma_x$ that comes from a classical treatment.
My biggest problem with a charged particle in a magnetic field is that the Bohr-van Leewen theorem apparently prevents this typical length scale (that I know for sure is $R_c$) to be found with a classical statistical treatment.

A:

Landau diamagnetism takes place because of the noncommutativity of 
$\mathbf{p}$ and $\mathbf{A}(\mathbf{r})$. In the classical treatment no such noncommutativity exists, thus  the magnetic susceptibility is identically zero. 
One way to appreciate the dependence of the result on $\hbar$ and at the same time perform the full quantum treatment is to perform the computation in the coherent state basis. In this basis, the Landau Hamiltonian has the form of an isotropic two dimensional harmonic oscillator with the Larmor frequency as a natural frequency, (for some detail, please see the following question .
$H = \hbar \omega_c (a^{\dagger}a + \frac{1}{2})$
Due to the nonvanishing commutation relation between the creation and
anihilation operators, the exponential of the Hamiltonian operator is given by:
$\exp(-\beta H) = \exp(-\frac{1}{2}\beta \hbar \omega_c) \exp(-\beta \hbar \omega_ca^{\dagger}a(1-e^{-\beta \hbar \omega_c}))$
In the coherent space basis, the partition function is given by:
$Z = \omega_c \int d^2\alpha \exp(-\frac{1}{2}\beta \hbar \omega_c) \exp(-\beta \hbar \omega_c\bar{\alpha}\alpha(1-e^{-\beta \hbar \omega_c})) = \omega_c (1-e^{-\beta \hbar \omega_c})^{-1}$
(The multiplicative factor proportional to $\omega_c$ is  the Jacobian of the transformation of the phase space volume). (I am being careless in the immaterial constant terms). 
This partition function gives the correct magnetic susceptibility. In addition, In the small \hbar limit the partition function becomes a constant, thus giving the Bohr-van Leewen theorem.

