Landau level degeneracy in symmetry gauge, finite system As we know, Landau level degeneracy in a finite rectangular system is $\Phi/\Phi_0$, where $\Phi=BS$ is the total magnetic flux and $\Phi_0=h/q$ is the flux quanta. This can be easily derived using Landau gauge and assuming the periodic boundary condition. 
However, if you choose the symmetric gauge, $l_z=x\hat p_y-y\hat p_x$ commutes with Hamiltonian, the corresponding quantum number $m$ is thus a good one. After some calculation, at last the energy level is written as $$E=\left[n+(m+|m|)/2+1/2\right]\hbar\omega.$$ I want to derive the degeneracy of Landau level for a finite system with radius $R$. 
The thought I got is that because in a finite system the angular momentum is bounded, the maximum value of $l_z$ is $qBR^2$ when you consider the particle doing the circular motion classically. Hence the maximum $m=l_z/\hbar=2\Phi/\Phi_0$, and the degeneracy is thus at least doubled. 
What we have learned is that the gauge choice won't change the observable effect, the finite degeneracy can obviously be observed. So what's wrong with my derivation? Or because they are different just because the system we considered is simply not equivalent?
 A: The problem with my thought is that I mix the "Canonical angular momentum" with the physical one. Since $\hat p_x,\hat p_y$ is canonical momentum, they are explicit gauge dependent. See What is canonical momentum? for a little reference.
The $\hat l_z$ in my question admit no physical meaning, so the bounded condition for it is just invalid. The true physical momentum should be 
$$
\begin{aligned}
\hat L_z=&x\hat P_y-y\hat P_x\\
   =&\hat l_z-\frac{qB}{2}r^2
\end{aligned}
$$
Consider classically of left hand side, $L_z=-qBr^2$, so at last we have: 
$$
l_z^{max}= -qBR^2/2
$$
So we have $|q|BR^2/2\hbar$ different values of $\bar m$ to choose, which give the right degeneracy: $\Phi/\Phi_0$
A: Good to see you solved the paradox yourself! Another way to find the same result is to compute the number of orbits which is possible to stack in a surface equal to the area of the system.
Let's first remind basic things on Landau levels. The hamiltonian of a particle moving in 2D $\{x,y\}\equiv\{r,\theta\}$ plan through a static magnetic field reads :
$$
\mathcal{H}=\frac{(\textbf{p}-q\mathbf{\mathcal{A}})^2}{2M}=\frac{\textbf{p}^2}{2M}+\frac{1}{8}M\,\omega^2_c\,r^2-\frac{\omega_c}{2}\,L_z
$$
where $\mathcal{A}=(-By/2,Bx/2,0)$ is the potential vector in the symmetric gauge, $L_z=xp_y-yp_x=-\mathrm{i}\hbar\,\partial_\theta$ the canonical angular momentum, and $\omega_c=qB/M$ the cyclotron pulsation.
Heisenberg inequalities provide here consistant physical constants of the problem (typical length $\ell$ and velocity $v_m$ of the movement) : 
$$
M\,v_m\,\ell \sim \hbar \quad\text{with}\quad v_m=\ell\,\omega_c
$$
Thus, we find $\ell=\sqrt{\frac{\hbar}{M\omega_c}}$ which is often refered as the magnetic length.
One can then compute the spectrum of $\mathcal{H}$ :
$$
\mathcal{Sp(H)}=\left[ E_{n,m}=\left(n-m+1\right)\frac{\hbar\omega_c}{2};n \geq 0, m=-n,-n+2 \,...\,n-2,n\right]
$$
where $n$ is the quantum number associated to the $\frac{\textbf{p}^2}{2M}+\frac{1}{8}M\,\omega^2_c\,r^2$ part of $\mathcal{H}$, and $m$ is the magnetic quantum number. Note that for a movement in the entire plan $\mathbb{R}^2$, Laundau Levels degeneracy $D$ is infinite.
If we now restrain our discussion to the Lowest Landau Level (LLL) given for $n=m$, it is easy to check that the eigen wave function is something like :
$$
\Psi_{n=m}(r,\theta)\propto \left(r\,e^{\mathrm{i}\theta}\right)^m e^{-(r/2\ell)^2}
$$
The essential features of $\Psi_{n=m}$ is that its RMS width is nothing more than the magnetic length $\ell$ and that it is maximum around a radius $r_{\mathrm{max}}=\sqrt{2m+1}\ell$. 
If we now consider that the movement of the particle is confined in a $R$ radius disk, the degeneracy $D$ of the LLL can be estimated by counting how many orbits one can put in the surface $S=\pi\,R^2$ of the system. Suppose that the system is large enough $S>>\ell^2$ to receive a large number $m>>1$ of orbitals, then according to the expression of $\Psi_{n=m}$, the condition to satisfy to be able to put these orbitals in $S$ is simply :
$$
R^2 > m\,2\ell^2
$$
Thus, the degeneracy is simply :
$$
D\sim\frac{R^2}{2\ell^2}=\frac{S}{2\pi\ell^2}=\frac{\Phi}{\Phi_0}
$$
$2\pi\ell^2$ can easily interpreted as the typical surface of an orbital.
