---
abstract: 'We examine the behaviour of a charged particle in a two dimensional confining potential, in the presence of a magnetic field. The confinement serves to remove the otherwise infinite degeneracy, but additional ingredients are required to produce sensible results. We treat both circular and square geometries, and in the latter we explicitly demonstrate the gauge invariance of the energy levels and wave function amplitudes. Both bulk states and edge states are examined, and in the latter case, with sufficiently high quantum numbers we achieve significant differences in the square and circular geometries. Results are achieved using straightforward matrix mechanics, in a manner that is accessible to novices in the field.'
author:
- Asadullah Bhuiyan and Frank Marsiglio
bibliography:
- 'References.bib'
title: '**Landau Levels and the Issue of Gauge Invariance in Confined Spaces**'
---

[Introduction]{}
================

The motion of a charged particle in the presence of a magnetic field is a subject of considerable interest in many areas of physics. In particular, in condensed matter the presence of the magnetic field alters the particle’s behaviour in a profound way. Ignoring the spin degree of freedom, it is known from classical physics that a vector potential is required to include the effect of a magnetic field, and more than one choice of vector potential is possible. Each of these choices constitutes a particular gauge. For a uniform magnetic field, say in the $z$ direction, common choices in Cartesian coordinates are the Landau Gauge, $\vec{\bf A}_1= xB \hat{\bf y}$ (or $\vec{\bf A}_2 = -yB \hat{\bf x}$), or the so-called Symmetric gauge, $\vec{\bf A}_S = -yB/2 \hat{\bf x} + xB/2 \hat{\bf y}$; all of these choices yield the same field $\vec{\bf B} \equiv \vec{\bf \nabla} \times \vec{\bf A} = B {\hat{\bf z}}$.

Naturally all physical quantities obtained in a calculation should be independent of the gauge choice. Typically a gauge choice is made to take advantage of some symmetry in the problem, so that one can proceed analytically. Just as often the physical system is altered to take advantage of the symmetry afforded by the gauge choice. A good example is the seminal paper by Laughlin, [@laughlin1981quantized] where he wanted to address the quantization of the Hall conductivity in a two-dimensional metal. To this end he chose a ribbon bent into a loop (i.e. open boundary conditions in one direction and periodic boundary conditions in the other). The magnetic field is perpendicular to the plane of the metal, which implies the physically questionable notion of a field either emanating from or converging to an axis at the centre of the loop. That this configuration did not coincide precisely with the experiment was immaterial — the concept of quantization was properly established and what mattered was that the Landau gauge allowed Laughlin to proceed analytically with the argument. Similarly, Halperin [@halperin1982quantized] followed up with a thin film with annular geometry, and adopted the gauge choice $\vec{\bf A} = [B_0r/2 + \Phi/(2 \pi r)] \hat{\bf \theta}$ to produce a uniform field $B_0$ in the $z$-direction along with a central flux of magnitude $\Phi$.

Other, “physical” geometries, suitable to so-called “quantum dots,” were adopted thereafter — see for example, Refs. \[\] and \[\]. In this paper we will review and expand upon some of these calculations, representing charged particles in a confined geometry subjected to a constant magnetic field. A very good reference for the quantum dot work is Chakraborty’s book, [*Quantum Dots*]{},[@chakraborty99] which also contains a number of very useful reprinted articles. As mentioned above, we will focus on the orbital motion of the charged particle (hereafter these will be electrons) and neglect the spin degree of freedom. Spin will cause a breaking of degeneracy in the presence of a magnetic field, owing to the direct coupling with field through the so-called Zeeman term. We begin in the next section with a textbook review of the free particle, from which the idea of Landau levels first arise. In particular we will explore both the Landau and the Symmetric gauges, and seemingly derive very different results for the wave functions. We further note that confinement within a parabolic trap can also be treated analytically, following Rontani.[@rontani99] Next we will review and expand upon Lent’s treatment of the circular quantum dot,[@lent91] and finally, we will show numerical results for the square quantum dot in both aforementioned gauges. This will hopefully make clear the role of degeneracy, and we explicitly illustrate the equivalence of the results in the two gauge formulations by removing the degeneracy. We will present the results of calculations for simple properties only, such as probability density and current density.

[Brief Review of Landau Quantization in Free Space]{}
=====================================================

In the presence of a magnetic field $\mathbf{B}$, the canonical momentum $\mathbf{p}$ is shifted by the magnetic vector potential $\mathbf{A}$ to give us a Hamiltonian of the form: $$H=\frac{\boldsymbol{\pi}^2}{2m}=\frac{(\mathbf{p}-q\mathbf{A})^2}{2m}
    \label{ham_original}$$ where $q$ is the charge and $m$ is the mass of the particle of interest, and $\boldsymbol{\pi}=\mathbf{p}-q\mathbf{A}$ is the *kinetic* momentum operator. Hereafter we adopt $q = -e$ for the electron, where $e>0$ is the magnitude of the charge of the electron. It is important to note that in the presence of a magnetic field, $\boldsymbol{\pi}$ represents the true momentum of the particle rather than $\mathbf{p}$. In this problem, we consider an electron free to move in a 2-dimensional system in a uniform magnetic field, $\mathbf{B}=B\boldsymbol{\hat{z}}$. For such a magnetic field, there are three common gauge choices, as already mentioned in the introduction. The two Landau gauges are handled very similarly mathematically, and so we will focus only on $\mathbf{A}_1$. However, $\mathbf{A}_1$ and $\mathbf{A}_S$ are handled very differently mathematically and we now discuss each in turn. A number of resources are available that provide considerable detail for this problem. These include textbooks like Refs. \[\], \[\], \[\] and \[\] and online lecture notes by Tong[@qheTong] and Murayama.[@Murayama]

[The Symmetric Gauge]{}
-----------------------

Substituting $\mathbf{A}_S$ into the Hamiltonian (\[ham\_original\]) leads to $$H=\frac{\mathbf{p}^2}{2m}-i\frac{\hbar\omega_c}{2}(x\partial_y-y\partial_x)+\frac{1}{8}m\omega_c^2(x^2+y^2),
\label{ham_symm1}$$ where $\partial_x \equiv \partial/\partial x$, and similarly for $y$ and later for radial coordinates as well, We have introduced the classical cyclotron frequency $\omega_c \equiv {eB}/{m}$. Converting to polar coordinates, we get $$H=-\frac{\hbar^2}{2m}(\partial_r^2+\frac{1}{r}\partial_r)+\frac{L_z^2}{2mr^2}+\frac{\omega_c}{2}L_z+\frac{1}{8}m\omega_c^2r^2,
\label{ham_symm2}$$ where $L_z= -i\hbar (x\partial_y - y \partial_x) = -i\hbar \partial_{\phi}$ is the standard operator for the $z$ component of the angular momentum. Note that this operator is *gauge-dependent* and its eigenvalues should not be considered to represent the real, physical angular momentum of the electrons in this system. To consider the *true* angular momentum of the system, we must take into account the effect of the magnetic field. Just like the canonical momentum is shifted to the kinetic momentum in the presence of a magnetic field, we can also define a kinetic angular momentum $\mathbf{\Lambda}$, $$\mathbf{L}=\mathbf{r}\cross\mathbf{p}\rightarrow\mathbf{\Lambda}=\mathbf{r}\cross\boldsymbol{\pi}=\mathbf{r}\cross(\mathbf{p}+e\mathbf{A})=\mathbf{L}+\mathbf{r}\cross e\mathbf{A}.
    \label{kin_ang_mom}$$ Since we are working in 2 dimensions, the kinetic angular momentum is perpendicular to the plane and is given by $$\mathbf{\Lambda}(\mathbf{r})=\left(L_{z}+\frac{1}{2}m\omega_cr^2\right)\boldsymbol{\hat{z}}=\Lambda_z(r)\boldsymbol{\hat{z}},
      \label{kin_ang_mom2}$$ where $\Lambda_z(r)$ is the magnitude of the kinetic angular momentum in the $\boldsymbol{\hat{z}}$ direction. We can also rewrite Eq. (\[ham\_symm2\]) in terms of $\Lambda_z(r)$, $$H=-\frac{\hbar^2}{2m}(\partial_r^2+\frac{1}{r}\partial_r)+\frac{\Lambda_z^2(r)}{2mr^2}.
       \label{ham_symm3}$$ Analogously, this new term can be interpreted as a kinetic centrifugal potential $V_{l}(r)$.

To calculate the eigenstates and eigenvalues of Eq. (\[ham\_symm2\]), we first note that since $[H,L_z]=0$, we can write the eigenstates $\Psi(r,\phi)$ as $$\Psi(r,\phi)=e^{i \ell \phi}\psi(r),
    \label{separable}$$ where $\hbar \ell$ is the eigenvalue of the $L_z$ operator. Using periodic boundary conditions, $\Psi(r,\phi+2\pi)=\Psi(r,\phi)$, requires $\ell$ to be an integer, and the eigenvalues of $L_z$ are thus quantized, with $\ell$ the ‘azimuthal’ or ‘angular momentum’ quantum number. This results in the radial ordinary differential equation, $$\left[-\frac{\hbar^2}{2m}\left(\frac{d^2}{dr^2}+\frac{1}{r}\frac{d}{dr}\right)+\frac{\hbar^2 \ell^2}{2mr^2}+\frac{1}{2}m\left(\frac{\omega_c}{2}\right)^2r^2 +\frac{\hbar \omega_c\ell}{2} \right]\psi(r) = E\psi(r).
      \label{ham_rad}$$ But Eq. (\[ham\_rad\]) is clearly just the 2D radial quantum harmonic oscillator Hamiltonian with frequency ${\omega_c}/{2}$ and with an additional energy shift by the constant term ${\hbar \omega_c\ell}/{2}$; the normalized 2D radial eigenstates are given by:[@jain2007composite] $$\Psi_{n_r,\ell}(r,\phi)=\frac{1}{\ell_B}\sqrt{\frac{n_r!}{\left(n_r+|\ell|\right)!}}\left(\frac{r}{\sqrt{2}\ell_B}\right)^{|\ell|}e^{-\frac{r^2}{4\ell_B^2}}L^{|\ell|}_{n_r}\left[\frac{r^2}{2\ell_B^2}\right]\frac{e^{i\ell\phi}}{\sqrt{2\pi}}
    \label{psi_rad}$$ where $L^{|\ell|}_{n_r}[x]$ are the associated Laguerre polynomials and $\ell_B \equiv \sqrt{{\hbar}/{eB}}$ is the magnetic length. The magnetic length plays an important role in what follows, and represents the characteristic length for the wave function of a charged particle in a magnetic field. For $B \approx 1$ Tesla, $\ell_B \approx 26$ nm. The quantum number $n_r$ is the ‘radial quantum number’ and simply counts the number of nodes that the radial part of the wave function has ($n_r$ is a non-negative integer). The corresponding probability density takes the appearance of concentric rings that are radially localized increasingly outwards from $r=0$ with increasing $|\ell|$. The wave function has $n_r$ nodes and the probability density $|\Psi|^2$ depends on $|\ell|$ and *not the sign* of $\ell$.

The corresponding eigenenergies are given by: $$E_{n_r,\ell} = \frac{\hbar\omega_c}{2}(2n_r+|\ell|+1) + \frac{\hbar \omega_c\ell}{2} =\hbar\omega_c\left(n_r+\frac{|\ell|+\ell}{2}+\frac{1}{2}\right),$$ so for all $\ell \leq 0 \ $, $E_{n_r,\ell} = \hbar\omega_c\left(n_r+\frac{1}{2}\right)$, and the energy spectrum is infinitely degenerate.

Sometimes a new quantum number $n \equiv n_r+{(|l|+l)}/{2}$ is defined, and then the energies are simply $$E_{n}=\hbar\omega_c\left(n+\frac{1}{2}\right).
    \label{energy_landau_level}$$ Clearly we see that further degeneracies are possible through the use of $\ell > 0$ for $\ell \leq n$. This is shown in the table below, where all energies have an infinite degeneracy due to the negative $\ell$ states that carry on indefinitely to the right for each energy. Additionally, however, as $n$ increases, positive $\ell$ states are also degenerate in growing number as indicated. $$\begin{aligned}
{\rm Energy} (\hbar \omega_c/2) \phantom{aaaaaaaa} n \phantom{aaaaaaaaaaaaaaaaaaaaa} (n_r,\ell) \ \ \phantom{aaaaaaaaaaaaaa} \cr
\vdots \phantom{aaaaaaaaaaaaaaaa} \vdots \phantom{aaaaaaaaaaaaa} \phantom{(0,2) \ \ (1,1) \ \ } \vdots \phantom{\ \ (2,-1) \ \ (2,-2)} \dots \cr
5        \phantom{aaaaaaaaaaaaaaaa} 2         \phantom{aaaaaaaaaa} (0,2) \ \ (1,1) \ \  (2,0) \ \ (2,-1) \ \ (2,-2) \dots \cr
3 \phantom{aaaaaaaaaaaaaaaa} 1        \phantom{aaaaaaaaaa} \phantom{(0,2)} \ \ (0,1) \ \  (1,0) \ \ (1,-1) \ \ (1,-2) \dots \cr 
1 \phantom{aaaaaaaaaaaaaaaa} 0 \phantom{aaaaaaaaaa} \phantom{(0,2) \ \ (1,1)} \ \  (0,0) \ \ (0,-1) \ \ (0,-2) \dots \cr
\end{aligned}$$

[Degeneracy of the Lowest Landau Level (LLL) Wave Function]{}
-------------------------------------------------------------

To better understand the wave functions and issues associated with them we focus on the lowest Landau level (LLL). For the LLL, $n_r = 0$ and $\ell \leq0$. Then $E_0 = \hbar \omega_c/2$ and the associated Laguerre polynomial is identically unity. So an infinite set of degenerate LLL wave functions is given simply by $$\Psi_{LLL}(r,\phi)=\frac{e^{i\ell\phi}}{\sqrt{2\pi}}\left(\frac{r}{\sqrt{2}l_B}\right)^{|\ell|}\frac{e^{-\frac{r^2}{4l_B^2}}}{l_B\sqrt{|\ell|!}}  \ \ \ \ \ \ell \le 0.
    \label{lll_wavefunction}$$ The corresponding probability density $|\Psi_{LLL}|^2$ is a radial Gaussian function centred around $r=0$ but with a maximum at $r_{\ell}=l_B\sqrt{2|\ell|}$ with a spread about this maximum of $\sim\mathcal{O}(l_B)$. As $|\ell|$ increases, the location of the maximum of $|\Psi_{LLL}|^2$, i.e. $r_{\ell}$, moves outwards [*from zero*]{}. On the other hand, increasing the magnetic field $B$ (thus decreasing $l_B$) leads to a decrease in the position of the maximum in $|\Psi_{LLL}|^2$.

It is important to note that these results are all gauge-dependent, even given that we use the Symmetric gauge. In particular, the choice of Symmetric gauge with $A_s$ defined in the first paragraph of the Introduction represents just one possible choice (out of infinitely many) where the vector field is referred to the origin. Since for this problem there is no preferred origin the wave functions will reflect this choice and be gauge-dependent. For this reason we highlighted “[*from zero*]{}” in italics above. An alternative would have been to use $$\mathbf{A}_S=\frac{B}{2}[(x-x_0)\boldsymbol{\hat{y}}-(y-y_0)\boldsymbol{\hat{x}}],
        \label{gauge_as}$$ and then the wave functions would have been different, even though this represents the same problem, i.e. a charged particle in free space in the presence of a uniform magnetic field $B$ (see, for example, Ref. \[\] for a discussion of this point).

Returning to the particular gauge choice given in Eq. (\[gauge\_as\]) with $x_0=y_0=0$, the values of $r_{\ell}$ also reveal that the amount of magnetic flux $\Phi_B$ associated with each wave function in the LLL is quantized in units of the magnetic flux quantum $\Phi_0={h}/{e}$. This is easily seen: consider the magnetic flux $\Phi_B$ between a circle of radius $r_{\ell}$ and that of (smaller) radius $r_{\ell + 1}$ (recall that $\ell \leq 0$, so that $r_{\ell} > r_{\ell +1}$ ): $$\Phi_B=B\pi (r_{\ell}^2-r_{\ell+1}^2) = 2 \pi B \ell_B^2(|\ell| - |\ell+1|) = 2 \pi B\frac{\hbar}{eB}=\Phi_0.
    \label{flux_quantum}$$ Thus, there is exactly a single quantum of magnetic flux between successive (negative) $\ell$ states in the LLL.

There are troublesome issues here. The infinite degeneracy is difficult, but even more so is the dependency of the wave function on gauge choice. There are some gauge invariant properties in such a system, for example, the spectrum of energies. The wave function, on the other hand, is [*not*]{} gauge invariant, and for a change in gauge, $\vec{\bf A}
\rightarrow \vec{\bf A^\prime}$, there is a corresponding expected change in wave function,[@griffiths18] $$\vec{\bf A^\prime} \equiv \vec{\bf A} + \vec{\bf \nabla}\Lambda  \Longrightarrow \Psi^\prime \equiv e^{iq \Lambda/\hbar} \Psi,
\label{psi_transform}$$ where $\Lambda$ is an arbitrary real function of position and time, and $q$ is the charge of the particle. So the probability density, $|\Psi|^2$ is expected to be gauge invariant. And yet, even within the symmetric gauge, Eq. (\[gauge\_as\]), different choices of $(x_0,y_0)$ seemingly result in wave functions that differ by far more than a phase (and therefore probability densities that [*do differ*]{} from one another). The reason for this problem is that a degeneracy exists in the solution. An explicit calculation of the wave function transformation corresponding to that in Eq. (\[psi\_transform\]) between the Landau gauge (to be reviewed below) and the Symmetric gauge has been given recently in Ref. \[\], and this problem was first highlighted in Ref. \[\], where the author presented a more general gauge condition to replace the simple wave function relation in Eq. (\[psi\_transform\]) when degeneracy is present. We will return to this issue later.

For the physical samples considered in this paper, the system is finite, and the electron is confined to some extent, so there must be some finite number of $\ell$’s that constitute the degeneracy number for a given Landau level (since the wave functions grow outward with increasing $|\ell |$). The standard argument (already provided above) is that the expected degeneracy $G$ of the $LLL$, given a sample of area $A$, is the ratio of the total magnetic flux $\Phi$ to the flux quantum, $$G=\frac{\Phi}{\Phi_0}=\frac{BA}{h/e}.
    \label{degeneracy}$$ As demonstrated later, this is not quite true, as sample edges confine the electron and result in higher energies for states near the edges.

Also note that the difference between centres of two successive $\ell$ eigenstates $\Psi_{n_r,\ell}$ and $\Psi_{n_r,\ell+1}$ close to the edge will be much less than the spread $l_B$ \[see Eq. (\[lll\_wavefunction\])\] and will decrease like $\sim\sqrt{|\ell|}-\sqrt{|\ell + 1 |}$, implying that the corresponding eigenstates will increasingly overlap with one another as $|\ell |$ increases.

[The Kinetic Centrifugal Potential]{}
-------------------------------------

Returning to the free electron case (still in the Symmetric gauge) note that the Landau energies clearly depend on the sign of $\ell$, and yet the wave functions (\[psi\_rad\]) clearly depend only on their absolute value, $|\ell|$. The situation where two identical probability densities correspond to completely different energies is not as strange as one might initially think; for example, the difference between a stationary vs. travelling wave packet is a simple phase factor and the energies of these two configurations differ by an arbitrary amount, as the phase factor determines the average velocity of the wave packet, and yet the probability distributions are identical in the two cases.

Here, two eigenstates with oppositely signed $\ell$ values ($|\ell|$ and $-|\ell|$), but with the same value of $n_r$, have energies $E_{n_r,|\ell|}=\left(n_r+|\ell|+\frac{1}{2}\right)\hbar\omega_c$ and $E_{n_r,-|\ell|}=\left(n_r+\frac{1}{2}\right)\hbar\omega_c$, respectively. Yet $|\Psi_{n_r,|\ell|}|^2=|\Psi_{n_r,-|\ell|}|^2$. To understand why, note that the effective potential, given by the kinetic centrifugal potential (\[ham\_symm3\]) is $$V_{\ell}(r)=\frac{\Lambda_z^2(r)}{2mr^2}=\frac{\left(\hbar \ell+\frac{1}{2}m\omega_cr^2\right)^2}{2mr^2},\\[5pt]
    \label{kin_cent_pot}$$ where $L_z$ in the original definition of $\Lambda_z(r)$, Eq. (\[kin\_ang\_mom2\]), is replaced by its eigenvalue $\hbar \ell$. Rewriting this as ($x \equiv r/r_{\ell}$)\
$$V_{\ell}(r)=\frac{\hbar \omega_c |\ell|}{4} \left( x^2 + 2\frac{\ell}{|\ell|} + \frac{1}{x^2} \right) \\[16pt]
    \label{kin_cent_pot2}$$ allows one to easily show that the minimum occurs at $x = 1$ ($r = r_\ell$) and the wave function (not surprisingly) is maximized where the effective potential is a minimum, [*independent of the sign of $\ell$*]{}. The value of the minimum potential is zero for all $\ell \le 0$ (independent of the magnitude of $\ell$), while it is $\hbar \omega_c |\ell|$ for $\ell > 0$ (grows linearly with $\ell$). In fact the difference between positive and negative $\ell$ is more general than this, and one can write that for [*all*]{} $r$, $$V_{|\ell|}(r) = V_{-|\ell|}(r) + \hbar\omega_c |\ell|,
  \label{pot_shift}$$ so there is a constant shift in potential values. The constant shift is consistent with attaining the same probability density for different energies.

![ $V_\ell(r)/(\hbar \omega_c)$ vs. $r/\ell_B$ for various values of $\ell$, both negative and positive. Solid lines correspond to negative $\ell$ while dashed lines correspond to positive $\ell$, and solid and dashed lines of the same colour correspond to the same magnitude, $|\ell|$ and hence have the same $r_\ell$ at which their minima occur. As explained in the text, for negative $\ell$ the minimum value of $V_\ell$ is always at zero while for positive $\ell$, the minimum value of $V_\ell$ increases linearly with $\ell$.[]{data-label="eff_pot"}](Kin-cent-potential.pdf)

Physically, electrons circulating with positive angular momentum in the presence of a magnetic field directed in the positive $z$ direction cost energy, whereas those with negative angular momentum do not. This is consistent with the classical idea that electrons speed up or slow down depending on whether they are circulating with positive or negative angular momentum, respectively.[@griffiths99] This effective potential is illustrated in Fig. (\[eff\_pot\]).

[Fock-Darwin States]{}
----------------------

Fock first addressed the problem of the eigenstates of a charged particle in a uniform magnetic field in 1928.[@fock28] He discovered a generally non-degenerate spectrum, particularly because the charged particle was confined to a parabolic potential. Somewhat later Landau[@landau30] examined a similar (simpler) problem in free space, and found a degenerate spectrum given in Eq. (\[energy\_landau\_level\]), which now bears his name. A year later, Darwin[@darwin31] independently obtained results similar to those of Fock, and figures displaying energy levels as a function of applied field are now called Fock-Darwin Spectra. The presence of a confining potential in the form of a parabolic trap requires very little additional work, so we now describe this work, also featured, for example, in Ref. \[\].

In the Symmetric gauge, with an additional confining potential of the form $m\omega_0^2r^2/2$, the Hamiltonian given in Eq. (\[ham\_symm1\]) becomes $$H=\frac{{\mathbf{p}}^2}{2m}+\frac{\omega_c}{2}L_z+\frac{1}{2}m\Omega^2r^2,
\label{ham_parabolic}$$ where $\Omega \equiv \sqrt{\omega_0^2+\omega_c^2/4}$. With this parabolic confinement, the eigenenergies are $$E_{n_r,l}=2\hbar\Omega \biggl(n_r + \frac{1}{2} [ |\ell | + \frac{\omega_c}{2 \Omega} \ell] + \frac{1}{2} \biggr).
    \label{energies_parabolic}$$ Figure \[fock\_darwin\] shows a typical Fock-Darwin plot. A variation on this plot is shown in Fig. (\[fock\_darwin\_fm\]). As in Fig. \[fock\_darwin\], Eq. (\[energies\_parabolic\]) is used for the energies, but now they are normalized to the energy $\hbar \omega_c$, and are plotted as a function of the (normalized) confinement potential.

![Fock-Darwin spectra obtained from using the symmetric gauge in a parabolic trap, i.e. Eq. (\[energies\_parabolic\]). The crossing of energy levels is the defining feature of the Fock-Darwin spectra. The dashed lines correspond to the free Landau levels. Note however, that there is no condensation of levels as the field increases, in contrast to what happens with more confining traps, to be discussed later. []{data-label="fock_darwin"}](FD-spec-fig-2.pdf)

![The same Fock-Darwin spectra obtained as in the previous figure by confining the electron with a parabolic trap, in the Symmetric gauge, but now $E_{n_r,\ell}/(\hbar \omega_c)$ is plotted vs $\omega_0/\omega_c$ for some selected levels. The numbers in brackets label the $(n_r,\ell)$ quantum numbers. The curves of a given colour all emerge from the horizontal dashed line drawn for that same colour; these denote the infinitely degenerate free Landau level energies ($\omega_0 = 0$), while the labelled curves of the same colour indicate how the degeneracy is broken with increasing confinement (increasing $\omega_0$). It is clear from those levels drawn here that the confinement plays a more important role for negative values of $\ell$ and increasing values of $| \ell |$.[]{data-label="fock_darwin_fm"}](fig3fd-parab.pdf)

Note that we could include all negative $\ell$ states, and these would rise even more than those shown as the confinement potential becomes stiffer. But at the far left of Fig. \[fock\_darwin\_fm\] ($\omega_0 \rightarrow 0$), we approach the free space limit, where the levels all become degenerate (this limit is indicated by the dashed horizontal lines).

[The Landau Gauge]{}
====================

Later on we will illustrate explicitly that correct gauge-invariant properties are obtained for various gauge choices. Here we will not demonstrate this explicitly, but instead point out that differences are apparent in the free space case, and what the root causes for these differences are. We therefore proceed in this section with the Landau gauge discussed in the Introduction, ${\vec{\bf A}}_1 = xB \hat{\bf y}$. Now our Hamiltonian is $$H=\frac{1}{2m}\left[p_x^2+\pi_y^2\right]=\frac{1}{2m}\left[p_x^2+(p_y+eBx)^2\right]$$ where $\pi_y$ is the kinetic momentum operator in the $\boldsymbol{\hat{y}}$ direction. Immediately, we see that $[H,p_y]=0$, so we can separate our eigenstate $\Psi$ as $$\Psi(x,y)\propto e^{ik_yy}\psi(x),$$ where the operator $p_y$ is replaced by its eigenvalue $\hbar k_y$. Substituting this into the Schrödinger Equation leads to the one dimensional differential equation $$\left[\frac{-\hbar^2}{2m}\partial_x^2+\frac{1}{2}m\omega_c^2(x+k_yl_B^2)^2\right]\psi(x)=E\psi(x),$$ which is the usual harmonic oscillator (albeit not centred at zero) equation. Our solution is therefore $$\Psi_{n,k_y}(x,y)=e^{ik_yy}\frac{1}{\sqrt{2^{n}n!
    \sqrt{\pi}l_B}}e^{-\frac{(x+k_yl_B^2)^2}{2l_B^2}}H_{n}\left[\frac{x+k_yl_B^2}{l_B}\right]
    \label{psi_landau}$$ where $H_n(z)$ are the Hermite polynomials. This method of solution has introduced two quantum numbers, $k_y$, a real number, and $n$, a non-negative integer. The corresponding energies are as before, in Eq. (\[energy\_landau\_level\]). Note that the wave functions are one-dimensional-like and certainly don’t resemble those obtained previously. The reason this is possible is because of the degeneracy present (energy does not depend on $k_y$), which we now discuss further.

[Degeneracy and the LLL Wavefunction in the Landau Gauge]{}
-----------------------------------------------------------

In the Landau gauge, the $LLL$ wavefunction must have $n=0$ for any $k_y\in\mathbb{R}$, giving us: $$\Psi_{LLL}(x,y)=e^{ik_yy}\frac{e^{-\frac{(x+k_yl_B^2)^2}{2l_B^2}}}{\sqrt{\sqrt{\pi}l_B}}.$$ Usually periodic boundary conditions are imposed in the $y$-direction, resulting in the requirement that $k_y = 2\pi m_y/a$, where $a$ is the length of the sample (i.e. “ribbon”) in the $y$-direction and now $m_y$ is an integer. We can estimate the degeneracy here as well. The magnetic flux $\Phi_B$ between the centres of two successive $m_y$ eigenstates $\Psi_{n,m_y+1}$ and $\Psi_{n,m_y}$ is $$\Phi_B=Ba(x_{k+1}-x_k)=Ba\left(\frac{2\pi (m_y+1)}{a}l^2_B-\frac{2\pi m_y}{a}l^2_B\right)=B\frac{2\pi\hbar}{eB}=\frac{h}{e}=\Phi_0,
   \label{degeneracy_landau}$$ so, as in the Symmetric gauge, each Landau gauge eigenstate contains roughly a single quantum of flux, $\Phi_0$. Thus, the expected degeneracy $G$ in a finite system of area $A$ is $G = BA/(h/e) = \Phi/\Phi_0$, as before. The infinite degeneracy in both cases allows the wave functions in the two gauges to differ by more than the simple “textbook” phase, as discussed by Swenson[@swenson89] and Wakamatsu et al.[@wakamatsu18]

In summary, the two gauges discussed both represent the same uniform magnetic field $\mathbf{B}=B\boldsymbol{\hat{z}}$. They result in the same energy spectrum, and yet they produce two very different-looking eigenstates (see Fig. \[fig: fs symm vs landau\]). There are gauge-invariant similarities as well; both gauges produce infinitely degenerate Landau levels in free space and both gauges have degeneracies that are expected to be approximately $G \equiv BA/\Phi_0$ in a finite space with area $A$. Eq. (\[psi\_transform\]) suggests that the two eigenstates should differ by a simple phase factor. That does not happen here, because of the degeneracy associated with the eigenstates. The dilemma is solved by taking an appropriate linear combination of the degenerate eigenstates (which is still an eigenstate), and these can be chosen to result in a state that differs from one in the other gauge by a simple phase factor, as illustrated in Refs. \[\] and \[\]. We will encounter this dilemma once more in a situation where the degeneracy has (nominally) been lifted, and outline a resolution by means different than that used in Refs. \[\] and \[\].

![Contour plots of the free space probability density vs. $x$ and $y$, in length units of $\ell_B$, for (a) the Symmetric gauge centred about the origin, with $\ell=5$ and (b) the Landau gauge, with $k_y\ell_B=5$. The probability densities in each gauge have strikingly different spatial structure. It is clear that these eigenstates are not related by just a simple phase factor.[]{data-label="fig: fs symm vs landau"}](free-space-landau-symm.pdf)

[Circular Confinement: Symmetric Gauge]{} {#circ}
=========================================

We will revisit the issue of degeneracy in different gauges in Section \[square\] when we discuss an electron confined to a square dot. First, however, partly to review results due to Lent,[@lent91] and partly to introduce some technical aspects and more general results to be used later, we present a discussion of results for circular confinement. We will tend to use dimensionless units wherever possible. However, it is important to convert our results to physical values. As we have already seen, the flux quantum plays an important role in these problems; its value in SI units is $h/e \approx 4140$ T (nm)$^2$. We will use a dimensionless measure of flux, $G = BA/(h/e)$, where $A$ is the area of the sample through which a (constant) magnetic field penetrates. For circular geometry, we use $\gamma \equiv B\pi a^2/(h/e)$, where $a$ is the radius. So, for example, applying a $10$ T magnetic field through a sample with radius 100 nm gives $\gamma \approx 76$.

[The Infinite Circular Well]{}
------------------------------

Let us first briefly review the problem of a particle confined to an infinite circular well. The infinite circular well potential $V(r)$ is simply given by: $$V(r) = \left\{
        \begin{array}{ll}
            0 & \quad r \leq a \\
            \infty & \quad r > a
        \end{array}
    \right. 
    \label{circ_pot}$$ where $a$ is the radius of the well. With this potential, the Hamiltonian in polar coordinates becomes: $$H=-\frac{\hbar^2}{2m}(\partial_r^2+\frac{1}{r}\partial_r)+\frac{L_z^2}{2mr^2}.
    \label{ham_polar}$$ Since $[H,L_z]=0$, we again write the eigenstates $\Psi(r,\phi) $ of the Hamiltonian as a separable product, $\Psi=e^{i\ell\phi}\psi(r)$, where $\hbar \ell$ is the same $L_z$ eigenvalue mentioned earlier and $\ell$ must be an integer. Substituting this into the Schrödinger Equation and introducing the wave vector $k=\frac{\sqrt{2mE}}{\hbar}$ results in $$\frac{d^2\psi}{dr^2}(r)+\frac{1}{r}\frac{d\psi}{dr}(r)+\left(k^2-\frac{l^2}{r^2}\right)\psi(r)=0.
\label{bess_eqn}$$ Equation (\[bess\_eqn\]) is the Bessel equation, with solutions given by $$\psi(r)=AJ_\ell(kr)+BN_\ell(kr)
    \label{bess_fn}$$ where $A$ and $B$ are arbitrary constants, and $J_{\ell}(r)$ and $N_{\ell}(r)$ are the $\ell^{\rm th}$ order Bessel and Neumman functions, respectively. Because the Neumman functions are non-normalizable at the origin, we discard this solution, and apply the boundary condition: $AJ_{\ell}(ka)=0$. The solution is given by $ka=\beta_{n,\ell}$, where $\beta_{n,\ell}$ is the $n$th zero of the $\ell$th order Bessel function. We thus obtain the eigenvalues, $$E_{n,\ell}=\frac{\hbar^2\beta_{n,\ell}^2}{2ma^2},
    \label{bess_eigen}$$ where $n$ is a positive integer that enumerates the zeros for the $\ell^{\rm th}$ order Bessel function. Notice the lack of degeneracy for all $n,\ell$. Bessel functions of the same order are both normalizable and orthogonal for different zeros over $r \in [0,a]$ for all integer $\ell$,[@Weisstein1] $$\int_{0}^aJ_\ell\left(\beta_{n,\ell}\frac{r}{a}\right)J_\ell\left(\beta_{n',\ell}\frac{r}{a}\right)rdr = \delta_{n,n'} \frac{a^2}{2}J_{\ell+1}^2\left(\beta_{n,\ell}\right).
\label{ortho}$$ Because of the orthogonality of $e^{i\ell \phi}$, we need not worry about the orthogonality relations for Bessel functions of different order. These solutions imply that a suitable basis set for the “empty” infinite circular well are the functions $$\Psi_{n,\ell}(r,\phi)=\frac{e^{i\ell\phi}}{\sqrt{2\pi}}\frac{\sqrt{2}}{a}\frac{J_\ell\left(\beta_{n,\ell}\frac{r}{a}\right)}{J_{\ell+1}(\beta_{n,\ell})}.
    \label{circ_basis}$$ Since $J_{-\ell}(r)=(-1)^\ell J_\ell (r)$ then wave functions with opposite angular momentum are degenerate, as expected in the absence of a magnetic field.

[Matrix Mechanics with Bessel Functions]{}
------------------------------------------

Matrix mechanics techniques have been utilized for problems like this before, albeit for a study of disordered systems, with a different gauge and confining potential.[@ohtsuki1988numerical] Here, when we include a magnetic field with the Symmetric gauge, we obtain $$H = \frac{\mathbf{p}^2}{2m} + \frac{\omega_c}{2}L_z+\frac{1}{8}m\omega_c^2r^2.
        \label{ham_circ}$$ To solve the Schrödinger equation with this Hamiltonian we follow the usual procedure of expanding in a basis (see, for example, Refs. \[,\]). We use $$\ket{\Psi}=\sum_{n^\prime=1}^{\infty}\sum_{\ell^\prime= -\infty}^{\infty}c_{n^\prime,\ell^\prime} \ket{n^\prime,\ell^\prime}$$ where $\ket{n',\ell'}$ is the ket corresponding to the circular well basis state, Eq. (\[circ\_basis\]). Acting on this state with the Hamiltonian, Eq. (\[ham\_circ\]), and taking the inner product with an arbitrary basis state $\ket{n,\ell}$ results in the matrix equation: $$\sum_{n^\prime=1}^{\infty}\sum_{\ell^\prime=-\infty}^{\infty} H_{n \ell, n^\prime \ell^\prime} c_{n^\prime,\ell^\prime} =E_{n,\ell}c_{n,\ell}
    \label{matrix1}$$ where $$H_{n \ell, n^\prime,\ell^\prime}=\bra{n,\ell}H\ket{n^\prime,\ell^\prime},
\label{hmatrix}$$ and, using Eq. (\[ham\_circ\]), $$H_{n \ell, n^\prime\ell^\prime}= \delta_{n,n^\prime}\delta_{\ell,\ell^\prime}E^{(0)}_{n,l}+\frac{1}{2}\omega_c\bra{n,\ell}L_z\ket{n^\prime,\ell^\prime}+\frac{1}{8}m\omega_c^2\bra{n,\ell}r^2\ket{n^\prime,\ell^\prime}.
  \label{hmatrix2}$$ Here, $E_{n,l}^{(0)}$ is simply the infinite circular well eigenenergy Eq. (\[bess\_eigen\]) and $\delta_{ij}$ is the standard Kronecker delta function. It should now be apparent that the entire matrix is diagonal in $\ell$, i.e. $ H_{n \ell,n^\prime \ell^\prime} = \delta_{\ell, \ell^\prime} H_{n,n^\prime}$, and the eigenvalue problem [*for a given $\ell$ only*]{} needs be solved. While a diagonalization is now required for each $\ell$ separately, this greatly reduces the computational cost. Evaluating the last two inner products of Eq. (\[hmatrix2\]), we obtain all the matrix elements for the Hamiltonian, $$H_{n,n^\prime}=\delta_{n,n^\prime}\left[\frac{\hbar^2\beta_{n,\ell}^2}{2ma^2}+\frac{\hbar\omega_c\ell}{ {2}}\right]+\frac{1}{8}m\omega_c^2a^2\rho^2_{n,n^\prime}.
\label{hmatrix3}$$ Using energy units of $\hbar^2/(2ma^2)$, we can write this in dimensionless form $$\frac{H_{n,n^\prime}}{\hbar^2/(2ma^2)} \equiv h_{n,n^\prime} = \delta_{n,n^\prime}\left[\beta_{n,\ell}^2+2{\gamma} \ell \right]+
    {\gamma^2}\rho^2_{n,n^\prime}
    \label{hmatrix4}$$ with $$\rho^2_{n,n^\prime}=\frac{\bra{n}r^2\ket{n^\prime}}{a^2}=\frac{2\int_{0}^1J_\ell(\beta_{n^\prime,\ell}\rho)\rho^2J_\ell(\beta_{n,\ell}\rho)\rho d\rho}{J_{\ell+1}(\beta_{n^\prime,\ell})J_{\ell+1}(\beta_{n,\ell})}.
    \label{rho_defn}$$\
Note that we have suppressed the dependency on $\ell$ in the matrix labels, but it is carried as a parameter in all of these matrices. Recall that the constant $G$ is related to the ratio of flux to the quantum of flux (given by $\gamma$), through $$\gamma = \pi G = \frac{B\pi a^2}{h/e}=\frac{\Phi}{\Phi_0}.
    \label{gamma_defn}$$ As was argued earlier, we expect $\gamma$ to roughly determine the degeneracy in this system. The integral in Eq. (\[hmatrix\]) can be evaluated analytically; we get $$\rho^2_{n,n^\prime} = \delta_{n,n^\prime} \Bigl[ \frac{\beta^2_{n,\ell} + 2\ell^2 -2} {3 \beta^2_{n,\ell}} \Bigr] + (1 - \delta_{n,n^\prime}) \Bigl[\frac{8\beta_{n^\prime,\ell}\beta_{n,\ell}}{(\beta_{n^\prime,\ell}^2-\beta_{n,\ell}^2)^2} \frac{J_{\ell-1}(\beta_{n^\prime,\ell})J_{\ell-1}(\beta_{n,\ell})}{J_{\ell+1}(\beta_{n^\prime,\ell})J_{\ell+1}(\beta_{n,\ell})} \Bigr]
\label{bess_int1}$$ We now have all the matrix elements analytically, to insert into Eq. (\[matrix1\]) for diagonalization. The big advantage of having adopted the symmetric gauge with this geometry is that the resulting Hamiltonian is diagonal in one of the quantum numbers, $\ell$. This represents a tremendous computational gain. We mention it here because we will [*not*]{} have this possibility when we study the square geometry in the next section. There we will have to diagonalize $N^2 \times N^2$ matrices, where $N$ is a large number (here typically something like $50$ or $200$. Also, for circular geometry we can monitor the quantum number $\ell$ precisely, which means we can track the radial quantum number $n_r$ as well. In the Symmetric gauge the Landau level quantum number $n$ is given simply by $n = n_r + (\ell + |\ell|)/2$. We will not have this possibility with the square geometry.

[Results and Discussion]{}
--------------------------

To be clear, in practice in this section we diagonalize $N \times N$ matrices, and only show results that have converged as $N$ increases (recall that $N$ should be infinite, but is chosen to be finite and represents the number of radial basis states (i.e. number of $n$) for a given $\ell$ in Eq. (\[circ\_basis\]) to be used in the matrix diagonalization). We will do this for a number of different $\ell$’s (typically $N$ of these), again so that all presented properties are converged as a function of basis size.

### [EigenEnergies]{}

Having carried out this programme, we can now order the eigenvalues in increasing value. We do this to mimic the result we would have achieved had we not recognized that the matrix is diagonal in $\ell$ (as will be the case with the square geometry). The eigenvalue results for a few values of $\gamma$ are shown in Fig. \[circular\_eigen\], and show in all cases plateaus (i.e. degenerate levels) separated by regions in which the energies increase to the next plateau. The number of points from the start of one plateau to the start of the next plateau is essentially $\gamma$.

![Dimensionless eigenvalues ($E_n/(\hbar^2/(2ma^2))$) as a function of total quantum number $n$ for various values of $\gamma$. Here, $n$ is used as a label for the horizontal axis to indicate the energies are plotted in ascending order and does **not** refer to any previously defined quantum numbers also labelled with $n$. The truncation sizes of the matrices produced are given by $N = 200$ in all cases. Note that these plots were produced as a result of generating and diagonalizing matrices for $\ell\in [-200,200]$ followed by aggregating the resulting eigenenergies acquired from all matrices and organizing them in ascending order. This results in over 80000 eigenvalues, and the lowest 200 of these are shown here. Note that the degeneracy increases with the value of $\gamma$ according to our expectations, and in fact there are precisely $\gamma$ eigenvalues occurring before the next plateau begins. Also note that the plateaus occur at dimensionless energies given by those expected from the free space Landau levels: $E_{n}/(\hbar^2/(2ma^2)) = 2\gamma(2n+1)$.[]{data-label="circular_eigen"}](circ-confine-eigenenergies.pdf)

The so-called “plateaus” are not degenerate to infinite precision, but they are degenerate to the precision of our computer (16 digit accuracy). These plateaus appear to be the Landau levels, originally infinitely degenerate, but now with a degeneracy of order $\gamma$, and in practice less. The increase in eigenvalues towards the end of the plateau regions is indicative that these eigenfunctions can feel the edge of the confining potential, and these are usually referred to as “edge states”, because, as we shall see, their probability density is concentrated there. One might ask, where have the rest of the (originally) infinite states in a given Landau level gone? Fig. \[circular\_eigen\] gives the impression that there were $\gamma$ such states, and about three quarters of them have remained degenerate, while the remaining one quarter have had their energies elevated due to the edge.

A clue that this interpretation is incorrect arises from something barely discernible on this plot, which is that the increase in eigenvalues at the ends of the plateaus become increasingly less smooth for successive plateaus. The reason for this is made clear in Fig. (\[circular\_eigen\_ell\]a), where now the energy levels are plotted as a function of $\ell$. Again we emphasize that we are fortunate that this is possible here, because different matrices were diagonalized for individual values of $\ell$. In Fig. (\[circular\_eigen\_ell\]b) we show the same energy levels vs. $r_\ell/a \equiv \sqrt{|\ell|/\gamma}$, to show how the maximum in the wave function moves towards the edge as $\ell$ increases (with the cautionary proviso that $r_\ell$ represents the maximum only for cases well away from the wall – which is why $r_\ell/a$ ends up exceeding unity in this plot — but the qualitative trend is correctly portrayed). Figure (\[circular\_eigen\_ell\]) illustrates that the eigenenergies are smoothly increasing when plotted vs. $\ell$. The raggedness that exists in Fig. (\[circular\_eigen\]) is due to the toggling back and forth between edges states of similar energy from the various Landau levels, uniquely identified in Fig. (\[circular\_eigen\_ell\]) \[(and not in Fig. (\[circular\_eigen\])\]. The Landau levels are distinct and clear, except now they carry on indefinitely, to the left in Fig. (\[circular\_eigen\_ell\]a), and to the right in Fig. (\[circular\_eigen\_ell\]b).

![Lowest 1000 Eigenenergies plotted vs. (a) $\ell$ and (b) $r_\ell/a=\sqrt{|\ell|/\gamma}$, for $\gamma=150$ and matrix truncation size, $\text{N}=200$. $\ell$ is the angular momentum quantum number and $r_\ell$ is the centre of the probability density of the corresponding eigenstate as long as the eigenstate is not too close to the walls. This means that $r_\ell$ gives us a good idea about where the eigenstates corresponding to the eigenergies are located in real space, as long as $r_\ell/a { { \atop <} \atop {\sim \atop }} 0.9$. Values of $r_\ell/a>1$ should be interpreted as corresponding to eigenstates piling up near the boundary. This plot (as opposed to Fig. \[circular\_eigen\]), makes clear that the original Landau level with infinite degeneracy ([*all*]{} negative values of $\ell$) has the very large $|\ell|$ states pile up near the boundary, with increasingly higher energies, thus breaking the original free Landau (infinite) degeneracy. As is clear from the figure, though, a quasi-degeneracy of order $\gamma$ (proportional to the applied magnetic field) remains.[]{data-label="circular_eigen_ell"}](E-vs-rl-l.pdf)

![The first 100 levels of the Fock-Darwin spectrum for the circular well. Also shown (with red dashed lines) are the free Landau levels. In this case, in contrast to that shown in Fig. (\[fock\_darwin\]) with parabolic confinement, the levels “condense” to these free Landau levels as the magnetic field increases, as illustrated by Lent (his Fig. 1).[@lent91][]{data-label="fock_darwin1"}](FD-spec-circ-confine.pdf)

As already stated and further corroborated below, the increase in energy is due to the confinement of the edge. The ‘toggling’ is evident if one examines the states at a dimensionless energy value like $\approx 1400$ in Fig. (\[circular\_eigen\_ell\]a). The next energy level that would be placed in Fig. (\[circular\_eigen\]) would come from one of the two branches either emanating from the first or the second Landau level, and these would go back and forth. At a higher value, say $\approx 2000$, toggling would occur between three branches, and so on, leading to increasing “raggedness” in Fig. (\[circular\_eigen\]) as the energy goes up, which is precisely what we observed.

Finally, we also include a Fock-Darwin spectrum, as in Lent,[@lent91] shown in Fig. (\[fock\_darwin1\]), where the so-called “condensation” of levels occurs.[@robnik86] As $\gamma$ increases, more and more levels coalesce as the free space Landau levels become more applicable since the wave functions become more compact and away from the edge. The thick dashed red lines indicate the expected Landau levels indicative of free space. To the far left, with the exception of $\ell = 0$ the eigenvalues come in pairs, corresponding to positive and negative $\ell$ which are degenerate in the absence of a magnetic field. Such condensation is absent with parabolic confinement, which is presumably why these are called ‘Landau’ levels and not ‘Fock’ levels.

### [Probability Density]{}

We have been discussing edge states, and using the parameter $r_\ell \equiv \ell_B\sqrt{2|\ell|}$ (recall $\ell_B  \equiv \sqrt{\hbar/(eB)}$ is the magnetic length) to indicate the location of the wave function. Of course $r_\ell$ is the radius of the maximum of the wave function (from some specified origin) in free space only, so here we examine the actual probability densities. Recall that the probability density is gauge invariant, except in cases of degenerate wave functions. Here the degeneracy is removed because of the confining cylinder, but as already discussed, a ‘practical’ degeneracy remains, so the statement about gauge invariance is no longer true,[@swenson89] as we shall see below and in the next section.

![Contour plots of $|\Psi_{n_r,\ell}(r,\phi)|^2$ for various values of $n_r$ and $\ell$, for $\gamma = 150$. The number of nodes increases with increasing $n_r$, more clearly seen for low values of $|\ell|$, and the radius of the maximal amplitude increases with increasing $|\ell|$. The notion that large $|\ell|$ states are “edge” states is clear from this plot, as the entire probability density resides on the circumference for these states. []{data-label="prob_density1"}](circ-confine-contours-l-nr.pdf)

![Contour plots of $|\Psi_{n}(r,\phi)|^2$ as a function of total quantum number $n$ for various values of $\gamma$, as indicated at the top of each column. Trends similar to those shown in Fig. (\[prob\_density1\]) are evident, but in addition one can note changes in $\ell$ and $n_r$ quantum numbers. For example, $|\Psi_{50}|^2$ for $\gamma = 50$ has $n_r = 0$ and $\ell = +1$, and is located at the start of the 2nd plateau visible in Fig. (\[circular\_eigen\]). It is identical to $|\Psi_2(r,\phi)|^2$ for $\gamma = 50$ displayed here in the 2nd row. The wave function $\Psi_2(r,\phi)$ has quantum numbers $n_r=0$ and $\ell = -1$ (recall that $|\Psi|^2$ was insensitive to the sign of $\ell$). Note that this plot is orderly and symmetric because we have artificially broken the remaining quasi-degeneracy apparent in Figs. (\[circular\_eigen\]) and (\[circular\_eigen\_ell\]) with a [*very*]{} shallow parabolic trap, centred at the middle of the cylinder. This shallow trap changes nothing that is visible to the eye, but it does remove the remaining quasi-degeneracy so that a proper ordering is established. See the text and the next section for further explanation of this additional potential, whose sole purpose is to remove the degeneracy.[]{data-label="prob_density1b"}](circ-confine-contours-G-state.pdf)

In Fig. (\[prob\_density1\]) we show contours of $|\Psi_{n_r,\ell}(r,\phi)|^2$ for $\gamma=150$ for various values of $n_r$ and $\ell$. These indicate that the wave function amplitude moves out from the centre as $|\ell|$ increases, as noted analytically for the free space result. We should note that Fig. (\[prob\_density1\]) looks very “orderly” when it comes to the progression of the number of nodes (as one moves to the right, increasing $n_r$) and the radius of maximum amplitude (as one moves down, increasing $|\ell|$). This is because we have complete control of the quantum numbers, as stressed with respect to Fig. (\[circular\_eigen\_ell\]). This scenario is different when we simply order all the states according to their energies \[this was done to produce Fig. (\[circular\_eigen\])\]. This will make the progression of states somewhat disordered, [*for a reason entirely different*]{} than the disorder already noted in connection with Fig. (\[circular\_eigen\]). [*That*]{} disorder was visible to the eye. Now we are referring to the quasi-degenerate states in the plateau regions, where the disorder [*is not visible, but nonetheless exists at a level lower than than $10^{-16}$*]{} (so it is invisible even to the computer, using double precision accuracy). For this reason, even in the plateau regions in Fig. (\[circular\_eigen\]), the states may be “out of order,” and a state with a higher $|\ell|$ may be ranked lower (according to energy) compared with a lower $|\ell|$ state. To avoid this we artificially include a [*very shallow*]{} parabolic trap, whose sole purpose is to break the remaining quasi-degeneracy, so that the progression of states is orderly. The resulting contour plots are shown in Fig. (\[prob\_density1b\]). We will discuss in more detail the introduction of this “degeneracy-breaking” potential in the next section. Including it as we do here does not change any of the results concerning energy that are visible to the eye; it does order the quasi-degenerate energy levels at the $10^{-10}$ level, so the computer “knows” which states are supposed to come first, and does not resort to a somewhat random linear combination, as happens when there is a degenerate subspace of solutions.

![Dimensionless radial probability density $(r/a)a^2|\psi_{n_r,l}(r/a)|^2$ vs $r/a$, for $\gamma = 150$, for various $|\ell|$ with $n_r=0$. We have used a matrix size to ensure that the results are converged ($N=200$). Note that only positive $\ell$ (i.e. $|\ell|$) values are used since the probability density does not depend on the sign of $\ell$. Also note that the area under each curve is unity, as required by normalization. Finally, note that the edge states are those with large values of $|\ell |$, as expected.[]{data-label="prob_density2"}](radial-prob-density-plots-circ-confine.pdf)

Because of the circular symmetry, it is enlightening to examine the probability density as a function of the radial coordinate, $r$. Defining a purely radial wave function through $\Psi_{n_r,\ell}(r,\phi) \equiv \frac{e^{i\ell\phi}}{\sqrt{2\pi}}\psi_{n_r,\ell}(r)$, we plot in Fig. (\[prob\_density2\]) the radial probability density, $({r/a}) a^2|\psi_{n_r,\ell}(r)|^2$ vs. $r/a$ for a variety of values of $\ell$, with $n_r = 0$, for $\gamma = 150$. These are properly normalized, and show clearly that the large $|\ell|$ states are edge states.

### [Probability Current]{}

In the Integer Quantum Hall Effect, edge states (eigenstates localized along a boundary) are often cited as a key ingredient to explain the quantization of hall conductance. [@qheTong] Classically speaking, the idea is that if an electron encounters a boundary while undergoing a cyclotron orbit, it will reflect off the wall and once again try to complete its orbit. This causes the electron to undergo ’skipping’ orbits along the boundary, resulting in chiral edge currents along the boundary of the sample. This was the motivation for examining the (probability) current associated with these edge states. One would expect that edge states would contribute non-zero current to the system while states in the bulk would not contribute any current. Lent[@lent91] illustrated such currents, and since the matrix mechanics technique allows for easy numerically exact calculations of this property, we will illustrate them here as well.

The probability current density $\mathbf{J}$ of an electron in an eigenstate $\Psi$ immersed in a magnetic field represented by a gauge choice $\mathbf{A}$ is given by [@landau77; @wysin] $$\mathbf{J}=\frac{1}{2m}\left[\Psi^*\mathbf{p}\Psi-\Psi\mathbf{p}\Psi^*+2e|\Psi|^2\mathbf{A}\right].
    \label{current1}$$ Since we are working in two dimensions, ${\mathbf{J}}$ is a *surface* current density with dimensions $IL^{-1}$. Using the symmetric gauge $\vec{\bf A}_s$, the probability current density associated with an eigenstate $\Psi_{n_r,\ell}(r,\phi)=\frac{e^{i\ell\phi}}{\sqrt{2\pi}}\psi_{n_r,\ell}(r)$ is, in polar coordinates, $$\mathbf{J}_{n_r,\ell}(r)=\frac{|\psi_{n_r,\ell}(r)|^2}{2\pi}\left[\frac{\hbar \ell}{mr}+\frac{1}{2}\omega_cr\right]\boldsymbol{\hat{\phi}}=\frac{|\psi_{n_r,\ell}(r)|^2}{2\pi}\frac{\Lambda_z(r)}{mr}
    \boldsymbol{\hat{\phi}}
    \label{current2}$$ where $\Lambda_z(r)=\hbar \ell+\frac{1}{2}m\omega_cr^2$ is the kinetic angular momentum in the $\hat{\boldsymbol{z}}$ direction. Note that this can also be expressed as $$\mathbf{J}_{n_r,\ell}(r)=\frac{|\psi_{n_r,\ell}(r)|^2}{2\pi}\sqrt{\frac{2V_\ell(r)}{m}}\boldsymbol{\hat{\phi}}
    \label{current3}$$ where $V_\ell(r)=\frac{\Lambda_z^2}{2mr^2}$ is the kinetic centrifugal potential. This form is particularly interesting because we see that ${\mathbf{J}}$ is simply the product of the radial probability density and the classical velocity of a particle in the radial potential $V_l(r)$. As can be inferred from Eq. (\[current2\]), and will be evident in the figures shown below, ${\mathbf{J}}$ has non-zero magnitude only in the $\boldsymbol{\hat{\phi}}$ direction, vanishes and has an inflection point at $r_\ell$ for $\ell<0$, is maximum at $r_\ell$ for $\ell>0$, and only vanishes at $r_\ell=0$ for $\ell=0$. To obtain the current from a surface current density, we integrate it with respect to a scalar differential ‘strip’ orthogonal to the direction of flow. The state current ${\mathbf{I}}_{n_r,\ell}$ of a given eigenstate is thus given by $${\mathbf{I}}_{n_r,\ell}=\int_{0}^{a} {\mathbf{J}}_{n_r,\ell}(r)dr.
    \label{current4}$$

In dimensionless length ($a$) and current density units ($\hbar/(2\pi ma^3)$, we have $$\mathbf{j}_{n_r,\ell}(x)=|\tilde{\psi}_{n_r,\ell}(x)|^2\left[\frac{\ell}{x}+\gamma x\right]\boldsymbol{\hat{\phi}}
    \label{current6}$$ where $x \equiv r/a$, $\tilde{\psi}_{n_r,\ell} \equiv a\psi_{n_r,\ell}$, and $\gamma \equiv B\pi a^2/(h/e)$ as before. The state current is now given by $${\mathbf{I}}_{n_r,\ell}/\frac{\hbar}{2\pi ma^2}={\mathbf{i}}_{n_r,\ell}=\int_{0}^{1} {\mathbf{j}}_{n_r,\ell}(x)dx.
\label{current7}$$ We can estimate roughly what to expect for a state far from the edge; for example, for $n_r=0$ in free space, we have $$|\tilde{\psi}^{(FS)}_{0,\ell}|^2=\frac{2\gamma^{|\ell|+1}x^{2|\ell|}}{|\ell|!}e^{-\gamma x^2},
    \label{current_fs}$$ so the free space (FS) state current magnitude $i^{(FS)}_{0,\ell}$ is $$\begin{aligned}
    i^{(FS)}_{0,\ell}&=\int^{\infty}_0\left[\frac{\ell}{x}+\gamma x\right]|\tilde{\psi}^{(FS)}_{0,\ell}|^2dx\\[5pt]
    &=\frac{2\gamma^{|\ell|+1}\ell}{|\ell|!}\int^{\infty}_0x^{2|\ell|-1}e^{-\gamma x^2}dx+g\underbrace{\int_0^{\infty}|\tilde{\psi}^{(FS)}_{0,\ell}|^2xdx}_{=1}\\
    &=\gamma \frac{\ell}{|\ell|!}\underbrace{\left(\int^{\infty}_0x^{|\ell|-1}e^{-x}dx\right)}_{=\Gamma(|\ell|)=(|\ell|-1)!}+\gamma \\
    &=\gamma [\text{sgn}(\ell)+1].
    \label{current9}\end{aligned}$$ Therefore, the state current in free space is given by: $$i^{(FS)}_{0,\ell}
=\begin{cases}
    2\gamma \ &\text{if} \ \ell>0\\
    \gamma \ &\text{if} \ \ell=0\\
    0 \ &\text{if} \ \ell<0.\\
    \end{cases}
    \label{current10}$$ Thus, the dimensionless state current is constant in free space with a value of $2\gamma$ for $\ell>0$ and with a value of $0$ for $\ell<0$. The case $\ell=0$ in free space actually holds for all $n_r$ due to the orthonormality of the eigenstate, $$i^{(FS)}_{n_r,0}=\gamma \int_0^{\infty}|\tilde{\psi}^{(FS)}_{n_r,0}|^2xdx=\gamma.
    \label{current11}$$ This case is also exactly calculable in our confined system as well, since we can expand $\tilde{\psi}_{n_r,0}$ as a series of Bessel functions: $$i_{n_r,0}=\gamma \int_0^{1}|\tilde{\psi}_{n_r,0}|^2xdx=\gamma.
    \label{current12}$$ To better understand how the bulk current of our confined system will behave in comparison, we obtain $$\begin{aligned}
    i_{0,\ell}&=\int^{1}_0\left[\frac{\ell}{x}+\gamma x\right]|\tilde{\psi}_{0,\ell}|^2dx\\[5pt]
    &=\int^{1}_0\frac{\ell}{x}|\tilde{\psi}_{0,\ell}|^2dx+\int^{1}_0\gamma x|\tilde{\psi}_{0,\ell}|^2dx\\[5pt]
    &=\int^{1}_0\frac{\ell}{x}|\tilde{\psi}_{0,\ell}|^2dx+\gamma
    \label{current13}\end{aligned}$$ The integral term left in the last line of Eq. (\[current13\]) is a bit tricky to calculate with a series expansion; since we are interested in the current in the bulk of the well (i.e. $|\ell|<<\gamma$) and the free space solutions are fairly localized, we approximate $|\tilde{\psi}_{0,\ell}|^2$ with its free space solution Eq. (\[current\_fs\]), $$\begin{aligned}
|\tilde{\psi}_{0,\ell}|^2&\approx|\tilde{\psi}^{(FS)}_{0,\ell}|^2 \ \text{for} \ |\ell|<<\gamma \\[5pt]
    \Rightarrow i_{0,\ell}&\approx\frac{2\gamma^{|\ell|+1}\ell}{|\ell|!}\int^{1}_0x^{2|\ell|-1}e^{-\gamma x^2}dx+\gamma \\[5pt]
    &=\gamma \frac{\ell}{|\ell|!}\left(\int^{\gamma}_0x^{|\ell|-1}e^{-x}dx\right)+\gamma\\[5pt]
    &=\gamma \left[\frac{\ell}{|\ell |!}\bar{\gamma}(|\ell|,\gamma)+1\right]
    \label{current14}\end{aligned}$$ where $\bar{\gamma}(|\ell|,\gamma)$ is the lower incomplete gamma function given by [@Weisstein2] $$\begin{aligned}
    &\bar{\gamma}(n,z)=\int^z_0e^{-t}t^{n-1}dt\\
    &\bar{\gamma}(n,z){\underset{z\rightarrow\infty}{\sim}}\Gamma(n)=(n-1)!
    \label{gamma_fn}\end{aligned}$$ where we have also shown the asymptotic form for large $z$. Thus, in our confined system, we expect the state current magnitude far from the boundary to approach the free space results Eq. (\[current10\]) for large $\gamma$: $$\begin{aligned}
    i_{0,\ell}&{\underset{\gamma \rightarrow\infty}{\sim}}\gamma \left[\frac{\ell}{|\ell|!}\Gamma(|\ell|)+1\right]=\gamma[\text{sgn}(\ell)+1]
    =\begin{cases}
    2\gamma \ &\text{if} \ \ell>0\\
    \gamma \ &\text{if} \ \ell=0\\
    0 \ &\text{if} \ \ell<0\\
    \end{cases}\nonumber \\
    \Rightarrow i_{0,\ell}&{\underset{\gamma \rightarrow\infty}{\sim}}i^{(FS)}_{0,\ell}.
    \label{current15}\end{aligned}$$ So, the current contribution from the $l=0$ state is given by $i_{0,0}=\gamma$.   This surprising result gives us both an expectation for our numerical results and implies that near the centre of the well, the probability current is nonzero. Moreover, it seems for very large $\gamma$, the current is expected to be a constant as a function of $\ell$ of a given sign. For $\ell<0$, the current contribution is approximately $0$, which is what would be expected in the bulk. However, for $\ell>0$ it is clear that the current contribution will be approximately $2\gamma$, which is unexpected because it implies that there exists non-zero current in the bulk of the sample.

![Probability current density $J$ (recall $\mathbf{J}=J\hat{\mathbf{\phi}}$) (in units of $\frac{\hbar}{2\pi ma^3}$) vs $r/a$ for (a) various negative and (b) various positive values of $\ell$, with $n_r=0$, and $\gamma=150$. Convergence was achieved with $N= 200$. The $\ell=0$ state is also shown for reference. In contrast to the negative $\ell$ states, the positive $\ell$ states all circulate in the direction expected classically (through the Lorentz force), but are energetically unfavourable. As with the negative $\ell$ states, for larger $|\ell|$ states, the current density is concentrated near the edge, but now positive. Note how the $\ell=0$ state has current density concentrated near the origin, while other (negative) $\ell$ states have current density located (naturally) where their probability density dominates, but with an inflection point \[as predicted from Eq. (\[current2\])\] right at $r = r_\ell$, [*as long as the probability density is well away from the edge.*]{}[]{data-label="prob_current2"}](radial-prob-curr-circ-confine.pdf)

![A vector field plot of the probability current for a “bulk” state in (a) ($\ell = -25$, $n_r = 0$) and an edge state in (b) ($\ell=-300$, $n_r = 0$). In both cases we used $\gamma=100$, and convergence was attained with a matrix size of $N=200$. The expanded portions show what should already be clear from Fig. (\[prob\_current2\]).[]{data-label="circ_vec_field"}](circ_confine_vec_field.pdf)

Numerical results are shown in Fig. (\[prob\_current2\]) for (a) negative $\ell$ states and for (b) positive $\ell$ states. The expectations from our earlier analytical analysis are borne out by these figures. In Fig. (\[circ\_vec\_field\]) we show a vector plot of the current for a “bulk” state in (a) ($\ell = -25$) and an edge state in (b) ($\ell = -300$), with the plots on the right showing the same result in more detail.

The calculation of current is shown in Figure (\[current\]); these are numerical results for the current in the $\hat{\phi}$ direction, and the agreement with the analytical estimate, Eq. (\[current15\]) is remarkable. Our numerical results are also in qualitative agreement with those shown by Lent,[@lent91] although he utilized smaller values of $\gamma$ (his $\beta \equiv 2\gamma$), and he showed pictorial vector plots for fairly small magnitudes of $\ell$ (which he calls $m$). He also provides a nice description in terms of classical orbits, which we will not repeat here; the reader is referred to Ref. \[\] for this description.

![Individual state current $I_{n_r=0,\ell}$ vs (a) $\ell$ and vs (b) $r_\ell/a$. We include (b) to illustrate how the constant current begins to be impacted by the edge, once $r_\ell$ is within a distance $\ell_B$ of the edge. Here $\ell_B \approx 0.06a$, and the deviation from the constant value predicted in Eq. (\[current15\]) starts shortly before this value is reached. Note the symmetry between negative and positive $\ell$ states. Also note how these numerical results agree very well with our analytical estimates, Eq. (\[current15\]).[]{data-label="current"}](curr-vs-rl-l.pdf)

[Square Confinement]{} {#square}
======================

Relatively little work has been done to date for a square geometry. In this section we present results in both the Symmetric and Landau gauges, and illustrate the difficulties encountered for large fields (or samples) due to the practical degeneracy that remains in these cases. We also note some peculiarities with the edge states for large quantum numbers.

[Symmetric Gauge]{}
-------------------

### [Matrix Element Calculation]{}

The confining potential is defined as $$V(x,y)=\begin{cases}
    0 \ &\text{if} \ 0< x < a \ \  {\rm and} \ \ 0<y<a\\
    \infty \ &\text{otherwise}
    \label{pot_confining_sq}
    \end{cases}$$ so $a$ is the length of a side of the two-dimensional infinite square well representing the square quantum dot. Given this potential, a convenient set of normalized basis states is $$\braket{x,y}{n_x,n_y}=\frac{2}{a}\sin{\left(\frac{n_x\pi x}{a}\right)}\sin{\left(\frac{n_y\pi y}{a}\right)},
    \label{basis_sq}$$ where $n_x$ and $n_y$ are positive integers. For convenience, we modify the symmetric gauge ${\mathbf{A}}_S$ such that it is centred in the potential well: $${\mathbf{A}}'_{S}=\frac{B}{2}\left[\left(x-\frac{a}{2}\right)\hat{{\mathbf{y}}}-\left(y-\frac{a}{2}\right)\hat{{\mathbf{x}}}\right].
   \label{symm_shifted}$$ With this gauge ${\mathbf{A}}'_S$, our Hamiltonian takes the form: $$H=\frac{\textbf{p}^2}{2 m_e}-i\frac{\hbar\omega_c}{2}\left[\left(x-\frac{a}{2}\right)\partial_y-\left(y-\frac{a}{2}\right)\partial_x\right]+\frac{1}{8}{m_e}\omega_c^2\left[\left(x-\frac{a}{2}\right)^2+\left(y-\frac{a}{2}\right)^2\right]
    \label{ham_sq}$$ where we use $m_e$ as the mass of the electron to avoid confusion with matrix indices. Following the steps leading to Eq. (\[matrix1\]) we then expand the eigenstate of (\[ham\_sq\]), $\ket{\Psi}$, as: $$\ket{\Psi}=\sum_{n_x,n_y}c_{n_x,n_y}\ket{n_x,n_y}
    \label{expand}$$ Thus, in energy units of $E_0 \equiv \hbar^2\pi^2/(2{m_e}a^2)$, and using $n\equiv (n_x,n_y)$ and $m \equiv (m_x,m_y)$ for short, our dimensionless matrix elements, $h_{n,m} \equiv H_{n,m}/E_0$ for the Hamiltonian for this problem are given by\
$$\begin{split}
    h_{n,m}&=\delta_{n_x,m_x}\delta_{n_y,m_y} \left[n_x^2+n_y^2+\frac{G^2}{2}\left(\frac{1}{3}-\frac{1}{(\pi n_x)^2}-\frac{1}{(\pi n_y)^2}\right)\right]\\[5pt]
        &+2 \left(\frac{G}{\pi}\right)^2 \Bigl( (1-\delta_{n_x,m_x})\delta_{n_y,m_y}  g_e(n_x,m_x)    + (1-\delta_{n_y,m_y})\delta_{n_x,m_x}  g_e(n_y,m_y) \Bigr)\\[5pt]
    &-i\frac{16G} {\pi^3}(1-\delta_{n_y,m_y})(1-\delta_{n_x,m_x})\left[\left(g_o(n_y,m_y)f(n_x,m_x)-g_o(n_x,m_x)f(n_y,m_y)\right)\right]
\end{split}
\label{ham_sq_dim}$$ with $$g_e(n,m)=\frac{1 + (-1)^{n+m}}{ 2}\left(\frac{1}{(n-m)^2}-\frac{1}{(n+m)^2}\right)=\begin{cases}
    0 \ &\text{if} \ n+m=\text{odd}\\
    \frac{1}{(n-m)^2}-\frac{1}{(n+m)^2} \ &\text{if} \ n+m=\text{even}
    \end{cases}
    \label{ge_def}$$ $$g_o(n,m)=\frac{1 - (-1)^{n+m}} {2}\left(\frac{1}{(n-m)^2}-\frac{1}{(n+m)^2}\right)=\begin{cases}
    \frac{1}{(n-m)^2}-\frac{1}{(n+m)^2} \ &\text{if} \ n+m=\text{odd}\\
    0 \ &\text{if} \ n+m=\text{even}
    \end{cases}
        \label{go_def}$$ $$f_o(n,m)=\frac{1-(-1)^{n+m}} {2}\frac{nm}{n^2-m^2}=\begin{cases}
     \frac{nm}{n^2-m^2} \ &\text{if} \ n+m=\text{odd}\\
    0 \ &\text{if} \ n+m=\text{even}
    \end{cases}
        \label{fo_def}$$ where $$G=\frac{Ba^2}{h/e} \equiv \frac{\Phi} {\Phi_0}$$ as defined earlier, and the subscripts ‘e’ and ‘o’ in the definitions, Eqs. (\[ge\_def\],\[go\_def\],\[fo\_def\]) serve to remind us that these are respectively non-zero for even or odd sums of the indices only.

### [Results and Discussion]{}

Note that a lower truncation is [*necessarily*]{} used here since matrices with block matrices as elements have to be constructed and diagonalized, which is much more computationally intensive than in the problem with circular confinement. Thus a truncation of $N=100$ for one of the indices, say $n_x$, implies a total truncation size of $N^2$ for both $n_x$ and $n_y$. We will refer only to the block matrix truncation size ($N$) in the following results (so $N=100$ values of $n_x$ requires diagonalization of a $10000 \times 10000$ matrix).

[*Eigenenergies*]{}

The computed eigenvalues are shown in Fig. (\[square\_eigen\]) for various values of $G$ as indicated. This figure shares many things in common with its circular counterpart, Fig. (\[circular\_eigen\]), and so it should, since for ‘bulk-like’ states we have argued that they resemble the free space results and hence, do not “feel” the edge, i.e. the electrons don’t even know if they are confined in a circular or square geometry.

![Eigenvalues shown in ascending order vs quantum number for $G=50$, $100$, and $150$. This figure is the ‘square’ version of its circular counterpart, Fig. (\[circular\_eigen\]), and we [*do not*]{} have the counterpart to Fig. (\[circular\_eigen\_ell\]), which would allow us to sort out here the two quantum numbers that are present in that case. Just like its circular well counterpart, the energies plateau at free space Landau level energies. Two noticeable features are the plateaus, as in the circular case, especially for large values of $G$, and the number of states present between plateaus, precisely $G$ here, just as it was $\gamma$ for the circular case. Also noticeable is a similar “raggedness” for energy values beyond the higher plateaus, presumably due to the same “toggling” phenomenon that we could explicitly identify in the circular case. Note that the free space Landau level energies are given by: $E_{n_L}/\left(\frac{\hbar^2\pi^2}{2ma^2}\right)=\frac{4G}{\pi}\left(n_L+\frac{1}{2}\right)$ (here we use $n_L$ for the Landau level number to avoid confusion; this plot has energy levels corresponding to the lowest 2 or 3 Landau levels). We used $N=150$ here.[]{data-label="square_eigen"}](square-confine-symm-eigenenergies.pdf)

However, the concept of a good quantum number $\ell$ does not exist here, at least not explicitly. So we do not have the benefit of a figure like Fig. (\[circular\_eigen\_ell\]) for the square geometry. One key similarity is the role of the ratio of the total flux to the flux quantum, which gives precisely the number of states between plateaus, regardless of geometry. In Fig. (\[square\_fock\_darwin\]) we show the Fock-Darwin spectrum for the case of square confinement. As with circular confinement the energy levels eventually “condense” and become coalesced into the degenerate Landau levels, as indicated. Both the circular and square confinement show this “condensation” phenomenon, whereas the parabolic confinement does not.\

![The Fock-Darwin spectrum for the two dimensional infinite square well. Here $N=100$. This plot shares with its counterpart for the circular well, Fig. (\[fock\_darwin1\]) the idea that for large enough field the levels “condense” to a set of degenerate Landau levels (indicated by thicker dashed red lines for the first 3 levels. For weaker fields the levels are fairly disordered, reflecting the more complicated geometry of the square. Note that the zero field energies no longer consistently come in pairs, as the degeneracy pattern for an infinite square well is more complicated than that for an infinite circular well. []{data-label="square_fock_darwin"}](FD-spec-square-confine.pdf)

[*Probability Densities*]{}

As we just remarked, for the so-called “bulk-like” states we expect the results for the square well to be very similar to those of the circular well, simply because the magnetic field keeps the electron sufficiently contained in the central region of the well so that the electron does not “know” the geometry of the confining potential. In Fig. (\[prob\_density\_sq1\]) we show contours of the probability density for a variety of values of G. This plot is to be compared with Fig. (\[prob\_density1\]), which showed probability densities for the circular well. In that plot, however, we could specify both $n_r$ and $\ell$, whereas here we can only order the plots according to the value of the energy. Some features of this plot are immediately apparent. First, results with small quantum number and small $G$ are consistent with what we found for circular geometry. This is expected. More excited states become fairly diffuse throughout the square, and retain the symmetry of the square; this is also expected. However, as $G$ increases, the results become somewhat irregular. The problem here was already alluded to in the case of the circular well: a high degree of “practical” degeneracy remains as $G$ increases. We fixed the problem with little explanation in that case, so now we discuss this issue in more detail.

![Probability density, $|\Psi_{n}(x,y)|^2$ contour plots for various values of total quantum $n$, as ordered by the diagonalization subroutine. Here the confining region is a square, outlined in black. The results only look sensible in the first two columns; the low quantum number results resemble those obtained from the circular geometry, as expected, and then the probability densities become more “square-like”, reflecting the geometry of the confinement, also as expected. As one increases $G$, however, the results become less understandable, and appear to be wrong. The difficulty, as discussed in the text, is the high degree of degeneracy that remains when large values of $G$ are used. We used $N=150$ for all the square well results.[]{data-label="prob_density_sq1"}](square-confine-symm-contour-eps0.pdf)

![Probability density, $|\Psi_{n}(x,y)|^2$ contour plots for the same parameters as in Fig. (\[prob\_density\_sq1\]), but now with an additional confining potential given by Eq. (\[hamp\]) with $\omega_0/\omega_c  =1.3\times 10^{-6}$. Now there is a more systematic progression, for example, for the ground state as one increases $G$; the wave function becomes more localized at the centre of the square as $G$ increases. In fact this could be happening in a circular geometry, and the result would be the same. As we move down, i.e. to more excited states for a given $G$, the wave function becomes more delocalized, and takes the shape of the edges in the square; these are veritable edge states. One notable exception is the result for $|\psi_{50}|^2$ with $G=50$, where the wave function becomes more concentrated; this is likely because we have suddenly moved on to the next Landau level (recall how we could track this explicitly in the case of the circular well). Note that in general the results compare quite well with those of the circular well \[see Fig. (\[prob\_density1b\])\], except fo the obvious “squareness” of the contours of the higher excited states.[]{data-label="prob_density_sq2"}](square-confine-symm-contour-eps1e-3.pdf)

Because of the remaining “practical” degeneracy, the diagonalization routine ends up picking some linear combination of these “degenerate” states and ordering them in some fashion. The arbitrary linear combination easily leads to a wave function that does not have the symmetry of the confining potential, as is evident in Fig. (\[prob\_density\_sq1\]). We will see the same problem emerge in the Landau gauge, which [*does not have*]{} the symmetry of the confinement potential, so it is prudent to emphasize that this aspect (having or not having a gauge with the symmetry of the confining potential) is not so important. Here, where the gauge choice is the symmetric one, we still encounter this problem, and it is because of the “practical” degeneracy (i.e. energies that are within $<10^{-16}$ of one another) that remains in the solutions. The ‘practical’’ degeneracy remains because, for large values of $G$, there is enough space in the confined potential for the “centre” of the wave function to be arbitrarily located.

We have resolved this problem by adding an additional confining potential with the symmetry of the square. Our choice is gauge invariant and is $$H' = \frac{1}{2} m_e \omega_0^2 \left[\left(x-\frac{a}{2}\right)^2+\left(y-\frac{a}{2}\right)^2\right],
 \label{hamp}$$ where $\omega_0$ is the characteristic frequency to characterize the perturbation potential. It is important to emphasize that this additional potential is minute: typically $(\omega_0/\omega_c)^2 \approx 10^{-6}$ for $G=50$, so that the energies are not really affected at any level the eye can detect, but the degeneracy is lifted sufficiently to allow the diagonalization algorithm to properly and unambiguously order the otherwise quasi-degenerate eigenstates. Most importantly, the centre of all wave functions is now correctly placed at the centre of the square.

When this tiny potential is included we obtain the results for the probability densities shown in Fig. (\[prob\_density\_sq2\]), where now the results look correct and more sensible \[the same tiny potential was actually used in Fig. (\[prob\_density1b\])\]. For the sake of completeness, we include the required matrix elements for our current problem, $$\begin{split}
    {(h')}_{n,m}= \frac{H'_{n,m}}{E_0} = \frac{\pi^2 \epsilon^2}{4} &\Biggl[ \frac{1}{12} \delta_{n_x,m_x}\delta_{n_y,m_y} \left( 1 - \frac{6}{(\pi n_x)^2} + 1 - \frac{6}{(\pi n_y)^2} \right)\\[5pt]
    &+\frac{2}{\pi^2} \Bigl( \delta_{n_y,m_y} [1-\delta_{n_x,m_x}] g_e(n_x,m_x)  +  \delta_{n_x,m_x}  [1-\delta_{n_y,m_y}] g_e(n_y,m_y) \Bigr) \Biggr]
\end{split}
\label{ham_sq_pert}$$ where $g_e(n,m)$ is as defined in Eq. (\[ge\_def\]) and $\epsilon \equiv \hbar \omega_0/E_0 = (4G/\pi)(\omega_0/\omega_c)$; with this definition, $\epsilon^2$ provides the relative energy scale of the perturbing potential.

[Landau Gauge]{}
----------------

Since the square geometry is rather difficult in either gauge, we thought it worthwhile to provide results for the same geometry in the Landau gauge. That this will work illustrates the remarkable machinery of gauge invariance, as now the gauge choice results in a very non-symmetric looking potential. In fact as we will shortly see it doesn’t work for large values of $G$, for the reasons just described above, but our very slight perturbing confining potential will fix the problem (i.e. lift the degeneracy) so that results identical to the symmetric gauge are achieved.

### [Matrix Element Calculation]{}

Now we use the (square well centred) Landau gauge, given by $${\mathbf{A}}'_L=B\left(x-\frac{a}{2}\right)\hat{{\mathbf{y}}}
    \label{landau_shifted}$$ so the Hamiltonian is $$H=\frac{\textbf{p}^2}{2{m_e}}-i\hbar\omega_c\left(x-\frac{a}{2}\right)\partial_y+\frac{1}{2}{m_e}\omega_c^2\left(x-\frac{a}{2}\right)^2.
    \label{ham_landau}$$ The Hamiltonian \[units of $\hbar^2\pi^2/(2{m_e}a^2)$\] becomes $$\begin{split}
    h_{n,m}&=\delta_{n_x,m_x}\delta_{n_y,m_y} \left[ n_x^2+n_y^2 + \frac{G^2}{3} \left( 1 - \frac{6}{(\pi n_x)^2} \right) \right]\\[5pt]
    &+8 \left(\frac{G}{\pi}\right)^2  (1-\delta_{n_x,m_x})\delta_{n_y,m_y}  g_e(n_x,m_x)  \\[5pt]
    &+i\frac{32G}{\pi^3}(1-\delta_{n_x,m_x})(1-\delta_{n_y,m_y}) g_o(n_x,m_x)f_o(n_y,m_y)
    \label{ham_landau_nm}
\end{split}$$ where the functions $g_e$, $g_o$, and $f_o$ were all defined earlier. This is a very different matrix than that generated for the symmetric gauge, and moreover is very asymmetric in $x$ and $y$. And yet we expect to obtain the same results as in that gauge.

### [Results and Discussion]{}

[*Eigenenergies*]{}

Here we could plot the eigenvalues generated in the Landau gauge, as in Fig. (\[square\_eigen\]), but we will not, because they are identical. Similarly, the results for the Fock-Darwin spectrum, Fig. (\[square\_fock\_darwin\]), are identical as well.

[*Probability Densities*]{}

Without using the perturbing potential, however, the probability density does not agree. These are shown in Fig. (\[prob\_density\_sq\_lan1\]), where significant discrepancies with Fig. (\[prob\_density\_sq1\]) are apparent, particularly for larger values of $G$.

![Contour plots of the probability density, $|\Psi_{n}(x,y)|^2$, for various values of $n$, calculated with the Landau gauge, and as ordered by the diagonalization subroutine. The results agree with those in the symmetric gauge, Fig. (\[prob\_density\_sq1\]), only in the first two columns. Once again, the reason for the discrepancy in the remaining columns is the degeneracy that is still present and this is remedied in the same manner as with the symmetric gauge results.[]{data-label="prob_density_sq_lan1"}](square-confine-landau-contour-eps0.pdf)

The results of Fig. (\[prob\_density\_sq\_lan1\]) show somewhat random character, particularly for larger values of $G$, in a manner similar to that seen earlier \[Fig. (\[prob\_density\_sq1\])\] in the symmetric gauge, and for the same reason — persistent numerical degeneracy. Nonetheless the addition of the perturbing potential, Eq. (\[hamp\]) with $\epsilon = 10^{-3}$ solves the problem as expected, and results identical (i.e. numerically, to about 5 digit accuracy) to those of Fig. (\[prob\_density\_sq2\]) are attained. These are shown in Fig. (\[prob\_density\_sq\_lan2\]).

![Contour plots of the probability density in the Landau Gauge with the perturbative trap included. While these results are identical to those obtained in the symmetric gauge in Fig. (\[prob\_density\_sq2\]), we have repeated them here to emphasize that this has been achieved with different gauges, and matrix diagonalization on completely different-looking matrices.[]{data-label="prob_density_sq_lan2"}](square-confine-landau-contour-eps1e-3.pdf)

[Other Gauge Possibilities]{}
-----------------------------

In the following section we present briefly our results on the probability current density. Before doing so, however, we should note that we explored other gauge choices, namely $${\mathbf{A}}^{\prime\prime}_{S}=\frac{B}{2}\left[\left(x-\frac{b_x}{2}\right)\hat{{\mathbf{y}}}-\left(y-\frac{b_y}{2}\right)\hat{{\mathbf{x}}}\right],
   \label{symm_shiftedb}$$ for the symmetric gauge, and $${\mathbf{A}}^{\prime \prime}_L=B\left(x-\frac{b_x}{2}\right)\hat{{\mathbf{y}}}
    \label{landau_shiftedb}$$ for the Landau gauge, where $b_x$ and $b_y$ are [*arbitrary*]{}. In particular we formulated the problem for $b_x = b_y = 0$, knowing that this was putting ourselves at a disadvantage. Nonetheless, even this choice works, but at a cost that much larger matrices are required. This is expected, since now the gauge potential is very asymmetric throughout the square — in either case the parabolic confinement arising from the applied magnetic field is centred at one corner of the square. While we refrain from showing results here, we found identical results as earlier for weak magnetic fields (where we could pursue convergence as a function of matrix size). It is noteworthy that in this regime there is no practical degeneracy and so the perturbation potential, Eq. (\[hamp\]), is not required.

[Probability Current]{}
-----------------------

The expression for the probability current density is given by Eq. (\[current1\]). Adopting the expansion of the wave function in the square well (\[expand\]) with basis states (\[basis\_sq\]) and arbitrary gauge choice $\mathbf{A}$, the probability current density of the $n$th total quantum state $\mathbf{J}_n$ can be expressed as: $$\begin{aligned}
    \mathbf{J}_{n}(x,y)&= \frac{\hbar}{m_e}\left[{\rm Im}({\Psi}_n^{\ast} {\mathbf{\Theta}}_n) +\frac{e\mathbf{A}}{\hbar}|\Psi_n|^2, \right] \ \ \ \ \text{where}
    \label{sq_prob_curr_dens}\\[5pt]
    \Psi_n(x,y)&={2 \over a} \sum_{m_x,m_y=1}^{\infty}c^{(n)}_{m_x,m_y} \sin\left(\frac{m_x\pi x}{a}\right)\sin\left(\frac{m_y\pi y}{a}\right), \ \ \ \ \text{and}\\[5pt]
    \begin{split}
     \mathbf{\Theta}_n(x,y)&\equiv \mathbf{\nabla}\Psi_n(x,y)={2 \over a}\sum_{m_x,m_y=1}^{\infty}c^{(n)}_{m_x,m_y}\frac{\pi m_x}{a}\cos\left(\frac{m_x\pi x}{a}\right)\sin\left(\frac{m_y\pi y}{a}\right) \ \mathbf{\hat{x}}\\[5pt]
     &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +{2 \over a}\sum_{m_x,m_y=1}^{\infty}c^{(n)}_{m_x,m_y}\frac{\pi m_y}{a}\sin\left(\frac{m_x\pi x}{a}\right)\cos\left(\frac{m_y\pi y}{a}\right) \ \mathbf{\hat{y}}.
     \end{split}\end{aligned}$$ Note that the eigenvector coefficient corresponding to the $n^{\rm th}$ eigenvector, $c^{(n)}_{m_x,m_y}$, is in general a complex number.

For the Symmetric \[Eq. (\[symm\_shifted\])\] and Landau \[Eq. (\[landau\_shifted\])\] gauge choices, the probability current density is given by $$\mathbf{J}^{(S)}_n(x,y)=\frac{\hbar}{m_e a}\left[a \ {\rm Im}({\Psi}_n^{\ast}{\mathbf{\Theta}}_n)+ \pi G|{\Psi}_n|^2\left(\left(\frac{x}{a}-\frac{1}{2}\right)\mathbf{\hat{y}}-\left(\frac{y}{a}-\frac{1}{2}\right)\mathbf{\hat{x}}\right)\right]
    \label{shifted_symm_prob_curr}$$ for the shifted-symmetric gauge, and $$\mathbf{J}^{(L)}_n(x,y)=\frac{\hbar}{m_e a}\left[a \ {\rm Im}({\Psi}_n^{\ast}{\mathbf{\Theta}}_n)+ 2\pi G|{\Psi}_n|^2\left(\frac{x}{a}-\frac{1}{2}\right)\mathbf{\hat{y}} \right]
    \label{shifted_landau_prob_curr}$$ for the shifted-Landau gauge. While the structure of the equations for the two gauge choices Eqs. (\[shifted\_symm\_prob\_curr\]) and (\[shifted\_landau\_prob\_curr\]) is quite different, the probability current density is a gauge-invariant quantity, and we have confirmed that both results are identical.

![A vector field plot of the probability current density with $G=100$ for a “bulk” state in (a) ($n = 10$) and an edge state in (b) ($n=95$). Further details of the region enclosed by the red square are provided in each case to the right. Convergence was attained with a matrix size of $N=150$ and $\epsilon=10^{-3}$ was used for the perturbing parabolic trap. Although this result was generated with the symmetric gauge, we confirmed that an identical result is found for the Landau gauge. In the first case an inner shell of clockwise-circulating current is followed by a concentric shell of counter-clockwise-circulating current; this is very similar to what was found with circular geometry \[see Fig. (\[circ\_vec\_field\]a)\]. However, the edge state current density is very different. As in the circular case it is [*primarily*]{} a clockwise-circulating current, but the corners cause a vortex-anti-vortex pair to be created, as is clear from the blow-up on the right.[]{data-label="fig: square_well_vec_field"}](square_confine_vec_field.pdf)

Equipped with equations (\[sq\_prob\_curr\_dens\]-\[shifted\_landau\_prob\_curr\]) it is straightforward to evaluate the probability current density for a given eigenstate. Two vector field plots of typical results are shown in Fig. (\[fig: square\_well\_vec\_field\]) (in units of $h/(m_e a^3)$). We see that (a) is a “bulk” state while (b) is an edge state. In fact, the former result is qualitatively indistinguishable from that attained for a circular geometry in Fig. (\[circ\_vec\_field\]), as is apparent from the result. In all cases we have also used a perturbing potential as described earlier, with $\epsilon = 10^{-3}$. In (b) there is a clear difference from the circular case, and the current tends to follow the geometry of the square boundary. While one may suspect from this that edge state currents simply follow the geometry of the boundary, a closer look at the expanded portion in (b) shows very interesting structure. Rather than smoothly following the geometry of the boundary at the corners of the square well, the probability current density appears to form stationary vortices rotating in a direction opposite to that of the nearby principal current density flow. This is a very different scenario than that encountered with the circular geometry, and is reminiscent of the ‘corner modes’ noted recently for topological insulators.[@yan18] Here, these are not zero energy modes, but appear as a simple consequence of the square geometry, and occur at large quantum number. Thus they may in fact be the precursor of semi-classical behaviour (e.g. the skipping orbits described in Ref. \[\]). Further investigation is clearly required.

One can also characterize the probability current density at slices through the sample. Defining $J_{x,n}\equiv\mathbf{J}_n\cdot\mathbf{\hat{x}}$, and $J_{y,n}\equiv\mathbf{J}_n\cdot\mathbf{\hat{y}}$, we show $J_{x,n}$ as a function of $y$ (a vertical slice at $x=a/2$) in Fig. (\[fig: x=a/2\_prob\_curr\_slice\]) through the centre of the square. For low quantum numbers the result will resemble that of the circle, shown in Fig. (\[prob\_current2\]). Here we cannot classify the states according to their $\ell$ quantum number, as we did in that case. Also note that Fig. (\[prob\_current2\]) displays results across a radius, i.e. [*half*]{} the sample, whereas for the square, Fig. (\[fig: x=a/2\_prob\_curr\_slice\]) shows results across [*the entire*]{} sample, and therefore has an inherent asymmetry, This is because, in spite of the [*square*]{} geometry, the current is [*circulating*]{}. For the same reason, $J_{y,n}$ along this vertical slice was found to be zero (within numerical noise). Slices related through symmetry operations yield entirely equivalent results. For example, the $J_{y,n}$ plotted across a horizontal slice at $y=a/2$ looks identical to $J_{x,n}$ at $x=a/2$.

![The $x$-component of the current density $J_{x,n}$ as a function of $y$ taken through $x=a/2$. This represents a vertical slice through the centre of the sample. We show results for various total quantum numbers $n$, as labeled. Results obtained in the Landau (Symmetric) gauge are shown as curves (points), and obviously agree with one another. Here we used $G=100$ and $N=150$. Also note that $J_{y,n}\sim0$ over this slice.[]{data-label="fig: x=a/2_prob_curr_slice"}](square-confine-prob-curr-gauge-compare-x-a-2-xcomponent-G100-N150-eps1e-3.pdf)

As seen before in the circular well, the current contribution for a particular state becomes non-zero as we move from the bulk region to the edge region in Fig. (\[fig: x=a/2\_prob\_curr\_slice\]). Moving from $y=a/2$ to $y=a$, we see that the total area under the lower states (with the exception of the ground state) vanishes due to a sign change across a node. States near the edge ($J_{x,95}$, for example) do not experience this nodal sign change, resulting in non-zero area, and thus non-zero current contribution. Also note that the edge state current has opposite sign on opposite edges of the confinement, confirming its inherent chirality.\

![$J_{y,n}$ as a function of $x$, taken along the diagonal from the upper left of the square down to the lower right. As $x$ varies from $0$ to $a$, $y = a-x$ varies also, from $a$ to $0$. As in Fig. (\[fig: x=a/2\_prob\_curr\_slice\]), results are shown for various quantum numbers and in both gauges. We used $G=100$ and $N=150$. $J_{x,n}$ mirrors this result about the vertical axis. []{data-label="fig: diagonal_prob_curr_slice"}](square-confine-prob-curr-gauge-compare-y-x-ycomponent-G100-N150-eps1e-3.pdf)

In Fig. (\[fig: diagonal\_prob\_curr\_slice\]) we show results across another slice, a diagonal across the square from the upper left to the lower right. We show only $J_{y,n}$ but in this case $J_{x,n}$ is simply the mirror of $J_{y,n}$ about the vertical axis. This plot is not so different from the previous one, but can give us a glimpse of the vorticies seen only at the corners of the square. The occurrence of these vortices clearly merits further study.

[Summary]{}
===========

We have presented an elementary discussion of the issues concerning gauge invariance, first in free space, and, more pertinently, in confined geometries. Practical degeneracies still exist in confined geometries, or quantum dots, in particular where a level “floor” in the potential well exists; this is not the case with a parabolic trap, for example. We used a quantum dot with circular symmetry to illustrate some of these degeneracies. However, in this case it is easy to overlook some of the subtleties, as we conveniently have a good quantum number, $\ell$, that can both simplify the calculation, and can be used to organize the results in an unambiguous manner. In this way, even states whose energy difference cannot be distinguished to numerical precision will nonetheless be ordered properly through our knowledge of these quantum numbers.

The circular potential also “begs” for the Symmetric Gauge to be utilized. Doing so simplifies the problem immensely, so only a one-dimensional equation requires solution. We could have used the Landau gauge, but then the problem would have been significantly more difficult. For the case of a two-dimensional square, however, the problem is more difficult right from the start, regardless of which gauge is chosen. Partly for this reason, it became a good testing ground for comparing gauge choices, especially given our method of solution, matrix mechanics. To our knowledge this problem has not been previously solved in either gauge.

Choosing a geometry also highlights that a gauge choice also involves a choice of origin for the gauge. Whether we use the Symmetric or Landau gauge, the natural origin is the centre of the sample, but we showed that this is not required. Non-optimal choices (i.e. [*not*]{} the centre of the sample) generally require more basis states with our method, so there is an additional numerical cost for a non-optimal choice. One can compute any property desired (we showed probability and current densities in this study), since we have the eigenvalues and eigenvectors. A further study will explore the susceptibility and other properties that can readily be measured.

Finally, we also suggested a simple way of removing the degeneracy at the $10^{-10}$ level (so the computer could tell the difference). As plotted, the eigenvalues will still appear to be degenerate, and this is the correct physics. However, for purposes of organization it is necessary to have a method to distinguish these from one another, and the (very) shallow gauge-invariant parabolic potential that we proposed does the trick.

The results for a square geometry are new; we expected that, for the properties studied in this work, the square geometry would not produce anything qualitatively new beyond the results for the circular quantum dot. Instead, as we saw, the current density for the edge states contains vortex-anti-vortex pairs at the corners, which are completely absent for the circular dot for the same range of quantum numbers. We plan to carry out a more in-depth investigation of the conditions under which such modes appear.

[Acknowledgments]{}
===================

This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC). These studies were initiated through a number of former students; in particular Sophie Taylor carried out some initial studies that were very helpful to get us started. We also thank Ania Michalik for helpful calculations and discussions. Mason Protter and Joel Hutchinson were instrumental in the initial stages with assistance with numerically technical matters. We also thank Joseph Maciejko and Jorge Hirsch for very helpful discussions on various aspects of this problem.
