💼 Case: charge $q$ in a homogeneous magnetic field along $z$-axis $\vec B=(0,0,B)$ with $B=|\vec B|$
We choose the following vector potential which we can verify works
$$ \vec A=(0,Bx,0) \quad \vec B=\vec\nabla \times \vec A =\left ( \frac{\partial}{\partial x}, \,\frac{\partial}{\partial y},\, \frac{\partial}{\partial z}\right )\times (0,Bx,0)=(0,0,B) $$
Since we don’t have an electric field $\Phi= 0$ so $V(\hat{\vec r})=q\Phi=0$ so using this we get
$$ \begin{aligned} \hat H &=\frac{1}{2m} (\hat{\vec p}-q\vec A )^2 \\ &=\frac{1}{2m}[\hat p_x^2 +(\hat p_y - q Bx)^2 +\hat p^2_z] \end{aligned} $$
🗒️ Note: the $z$-component is completely decoupled from $xy$ motion thus $[\hat p_{z},\hat H]=0$
Here since along the $z$ axis we just have $\hat H_z=\hat p^2/2m$ then we use a known wave plane solution
$$ \psi(x,y,z)=\psi(x,y)e^{ik_zz} $$
where $k_z=p_z/\hbar$ with $p_z$ an eigenvalue of the operator $\hat p_z$
Looking at the $xy$ plane we see that the Hamiltonian does not depend on $y$ thus
$$ \begin{aligned} \hat H_{xy}&=\frac{1}{2m}[\hat p^2_x+(\hat p_y -qBx)^2] \\ &=\frac{1}{2m}[\hat p_x^2 + \hat p_y ^2 + q^2 B^2 x^2 - 2qBx\hat p_y] \end{aligned} $$
where we used $[x,\hat p_y]=0$ and we see that $[\hat p_y,\hat H]=0$ since there is no $y$ dependence
We can thus diagonalise $\hat p_y$ to give (set it to an eigenvalue of the Hamiltonian)
$$ \psi_{xy}(x,y)=\psi_{k_y}(x)e^{ik_yy} $$
where $k_y=p_y/\hbar$ with $p_y$ an eigenvalue of the operator $\hat p_y$
The Schrodinger equation for the $xy$ plane is thus
$$ \begin{aligned} \frac{1}{2m} [\hat p_x + \hbar ^2 k_y ^2 +q^2 B ^2 x^2-2qB\hbar k_y x] \psi_{k_y} (x)&=E_{xy} \psi_{k_y}(x) \\ \frac{1}{2m} \left [ \hat p_x^2 + q^2 B^2 \left ( x-\frac{\hbar k_y}{qB} \right )^2 \right ] \psi_{k_y}(x) &=E_{xy} \psi_{k_y}(x) \end{aligned} $$
Subbing in $x_0=(\hbar k_y)/(qB)$ and $\omega_c=(qB)/m$ (cyclotron frequency) we get
$$ \left [ \frac{1}{2m} \hat p_x^2 +\frac{1}{2} m\omega c^2 (x-x_0)^2 \right ]\psi{k_y}(x)=E_{xy} \psi _{k_y} (x) $$
Here we recognise simple harmonic motion thus
$$ \psi_{k_y}(x)=\Phi_n (x-x_0)=\Phi_n \left ( x-\frac{\hbar k_y}{qB} \right ) $$
where $\Phi_n(x)$ is the $n$-th eigenstate of the harmonic oscillator of frequency $\omega_c$
Combining the $xy$ plane we thus have
$$ \psi_{xy}(x,y)=\Phi_n (x-x_0) e^{ik_y y} $$
defines the in-plane eigenstates with energies
$$ E_n=\hbar \omega_c \left ( n+\frac 12 \right ) \qquad \forall n\in \N $$
🗒️ Note: all states with the same $n$ (fix the $x$-dependence) are degenerate (due to $k_y$)
📖 Definition: Landau level, set of degenerate states for a fixed value of $n$
💼 Case: here lets instead chose $\vec A=(-yB,0,0)$ since $\vec B = \vec \nabla \times \vec A =(0,0,B)$ it works
To find the Gauge transform here we can see that the two vector potential differ by
$$ \vec \nabla \lambda =-(yB,Bx,0) \qquad \Rightarrow \qquad \lambda =-Bxy $$
we have the following unitary gauge transformation operator to transform between these states
$$ \hat G_\lambda = \exp \left (-i\frac{qBxy}{\hbar} \right ) $$
The Schrodinger equation in the new gauge becomes
$$ \left [ \frac{1}{2m} \hat p^2 y +\frac 12 m\omega c^2 (y+y_0)^2 \right ]\psi{k_x}(y)=E\psi{k_x} (y) $$
with $y_0=\hbar k_x /qB$
The $xy$-plane eigenstate is now
$$ \psi'_{xy}(x,y)=e^{ik_x x} \Phi_n (y+y_0) $$
with the same energies
$$ E_{n,k_x} =\hbar \omega _c \left ( n+ \frac 12 \right ) \qquad \forall n \in \N $$
💎 Conclusion: due to Landau levels (large degeneracy), transforming the eigenstate in one gauge by $\hat G_\lambda$ will give superpositions(linear combination) of eigenstates in another gauge
💼 Case: going back to $\vec A= (0,Bx,0)$
We rewrite the Hamiltonian as
$$ \begin{aligned} \psi(x,y,z)&=\psi_{k_y} (x) e^{ik_yy} e^{ik_zZ} \\ &=\Phi_n (x-x_0) e^{ik_yY} e^{ik_zZ} \\ &=\mathcal N H_n (x-x_0)\exp \left ( -qB \frac{(x-x_0)^2}{2\hbar } \right ) e^{ik_y y} e^{ik_z z}
\end{aligned} $$
where $H_n (x)$ are Hermite polynomials and $\mathcal N$ is an appropriate normalisation.
We can then sub in to the complete $3$-dimensional Schrodinger equation
$$ \frac{1}{2m} \left [ \hat p_x^2 + q^2 B^2 \left ( x-\frac{\hbar k_y}{qB} \right )^2+\hbar^2 k^2_z \right ]\psi_{k_y}(x)=E\psi_{k_y} (x) $$
We get a continuous energy spectrum ($k_z$) but with discrete Landau levels for fixed $k_z$
$$ E_n=\hbar \omega_c \left ( n+\frac 12 \right )+ \frac{\hbar^2 k_z^2}{2m} \qquad \forall n \in \N $$
🗒️ Note: $E_n$ is independent of $k_y$ meaning that for a fixed $n,k_z$ we can chose $k_y$ (degenerate)