What's the coordinate of the Landau tube?

The Hamiltonian of a three-dimensional electron gas in a static magnetic field $$\vec{B}$$ is : $$\hat{H}=\frac{(\hat{\vec{p}}+e\hat{\vec{A}})^2}{2m}$$

If choosing $$\vec{B}=B\vec{e_z}$$ and using Landau gauge : $$\vec{A}=(0,Bx,0)$$, we can rewrite the Hamiltonian in a quasi harmonic oscillator form:

$$\hat{H}=\frac{\hat{p}_x^2}{2m}+\frac{1}{2}m\omega_c^2\left(\hat{x}+\frac{\hat{p}_y}{m\omega_c}\right)^2+\frac{\hat{p}_z^2}{2m}$$ where $$\omega_c$$ is the cyclotron frequency.

For the first two terms, they are just in harmonic oscillator form, so the eigen-energy corresponding to this part is $$(n+\frac{1}{2})\hbar\omega_c$$. For the last term, since the magnetic field doesn't break the translation symmetry along z-direction (i.e., $$[\hat{H},\hat{p}_z]=0$$), the eigen-energy corresponding to this part is $$(\hbar k_z)^2/2m$$.

Hence we achieve the eigen-energy of the Hamiltonian with two quantum number:

$$E_{nk_z}=\left(n+\frac{1}{2}\right)\hbar\omega_c+\frac{(\hbar k_z)^2}{2m}$$and use landau tube to describe the degeneracy for each eigen-energy:

We can find that for a fixed $$n$$ and $$k_z$$, the states of the same energy are just denoted by a circle with centre axis along $$k_z$$ axis and satisfying $$\frac{\hbar^2 (k_x^2+ k_y^2)}{2m}=(n+\frac{1}{2})\hbar\omega_c$$.

Here is my question : Are the coordinates of Landau tube ($$k_x$$, $$k_y$$ and $$k_z$$) corresponding to three components of mechanical momentums ($$\hat{\vec{\pi}}=\hat{\vec{p}}+e\hat{\vec{A}}$$) or canonical momentums($$\hat{\vec{p}}$$)?

If they are mechanical, why can we use both $$k_x$$ and $$k_y$$ to label a state when $$[\hat{\vec{\pi_x}},\hat{\vec{\pi_y}}]\neq0$$?

If they are canonical, why can we easily find the states of same energy by satisfying $$\frac{\hbar^2 (k_x^2+ k_y^2)}{2m}=\left(n+\frac{1}{2}\right)\hbar\omega_c$$ when the form of Hamiltonian is not simple $$\frac{\hat{p}_x^2+ \hat{p}_y^2+\hat{p}_z^2}{2m}$$?

The coordinates of the Landau tube should be the canonical momentum, since you have already replaced the vector potential in the kinetic momentum by $$Bx$$ which led to the displaced harmonic oscillator in the $$x$$-direction. And note that your mechanical and canonical momentum are essentially the same in the $$y$$- and $$z$$-direction.
With respect to your follow up question, I am not quite sure if I understand it correctly, but, the simple form $$\frac{p_x^2+p_y^2+p_z^2}{2m}$$ is the Hamiltonian for a free particle in all dimensions which is apparently not the case here. Thus the energy has the constraint that $$\frac{\hbar^2 (k_x^2+k_y^2)}{2m} = \hbar \omega_c(n+\frac{1}{2})$$ and is not simply given by $$E = \frac{\hbar^2 (k_x^2+k_y^2+k_z^2)}{2m}$$.