text:
# Square billiard with a magnetic flux

## I Introduction

The behavior of mobile charged particles confined to some region and also subjected to a magnetic flux has interested physicists since the discovery of the Hall effect over a century ago. Quantum effects turned out to be subtle and surprising, as attested by Landau diamagnetism and the Aharonov-Bohm effect. About twenty years ago, with the advent of the quantum Hall effect and mesoscopic systems, the two dimensional case became prominent. At about the same time, the development of the subject of “quantum chaos” also focussed interest on such systems as being among the simplest of the “Gaussian unitary ensemble” universality class.

There is thus a long history of work on the confined quantum motion of charged particles in a magnetic flux. It is remarkable that all previous workers overlooked the fact that many fundamental cases, some of which have been extensively studied numerically, can be solved and classified analytically to good approximation. Moreover, the results have an interesting and suggestive complexity.

In this paper we obtain good approximate solutions to a couple of simply posed and well studied problems of this type. A preliminary version has appeared electronically. Rather than stress the generality of our method, we focus on a typical problem: a two dimensional charged particle confined to a square billiard in a perpendicular magnetic flux. Certain conditions on the flux are required to justify the approximations, and we also require the energy of the particle to be large. These approximate solutions are compared to numerical solutions. The two flux configurations considered are a uniform flux and an Aharonov-Bohm flux line. It is crucial that the square billiard is integrable. Systems other than the square to which our methods apply will be mentioned at the end.

There are many recent research papers in which the basic system studied is a square or rectangular billiard with a magnetic flux. A number of these are inspired by the experiments of Lévy et al., which measure the magnetic susceptibility of a collection of a considerable number of mesoscopic, two dimensional metallic systems, each approximately a square. The field is nearly uniform over the square in this case.

In the presence of a magnetic flux, there is the possibility that persistent currents exist. In other words, it is possible that the equilibrium state has a nontrivial current. Indeed, we find that eigenstates of the quantum system have interesting current densities. Although the wavefunction has a fairly simple representation, the currents can be quite complex. Some states are predominantly paramagnetic, others are predominantly diamagnetic. Still others may support both strong paramagnetic and strong diamagnetic currents which in total nearly cancel. However, these states strongly affected by the field are rather rare and most states have weak persistent currents.

Ideas from the field of quantum chaos have also motivated much work. Since the magnetic field intuitively has circular symmetry, which ‘conflicts’ with the symmetry of the square, one might expect chaos to ensue, as is the case with the Sinai billiard. Another theme of quantum chaos is that of energy level and wavefunction statistics. These statistics depend on whether time reversal symmetry (or other antiunitary symmetry) is in force. A natural way to break time reversal is by a magnetic flux. \[The square with a uniform flux still has an antiunitary symmetry, however.\]

Diffraction effects, in which a classical length shorter than the wavelength becomes important, are much studied in this context. This is obviously the case for the zero radius Aharonov-Bohm flux line. The sharp corners of the square also cause diffractive effects in the presence of a uniform magnetic field. These effects are much smaller than for the flux line, of course. We give estimates for the parameter range in which such diffraction becomes important, although we defer study of these effects.

Our main motivation however, is that we add to the store of solvable problems, and perhaps suggest some experiments. In the textbooks, there are relatively few such integrable problems, basically, only those which reduce to one dimension, or separate into several one-dimensional problems. The square without a magnetic flux is such a case in which the $`x`$ and $`y`$ motion separates.

The traditional way to widen the class of approximately solvable problems is perturbation theory and indeed, our approach is a form of quantum perturbation theory. We are able to study systems which are integrable except for a ‘classically weak’ perturbation. Of course, treating weak perturbations classically is challenging, because the long time behavior may be chaotic. However, quantum perturbation theory is better behaved and depends on the short time rather than the long time classical behavior.

On the other hand, perturbations small in this sense can give rise to very large quantum effects, especially on the wave functions. Moreover, even perturbations which change some length scale by an amount $`\delta L<<\lambda ,`$ where $`\lambda `$ is the wavelength, can have big effects. In terms of the standard nomenclature, a degenerate perturbation theory is required, and a number of unperturbed states are strongly mixed together to give the final result. Usually, this is done by diagonalization of a small matrix, but in our case the ‘matrix’ can be quite large. However, rather than just diagonalizing some matrix by a computer calculation, we obtain the result by an intuitively appealing Schrödinger differential equation.

Having solved, for the first time, this interesting class of problems, we discovered some other methods of obtaining the solution at the same level of approximation. We shall present these methods elsewhere.

## II Square in a uniform field

We begin with the case of uniform field. For given velocity $`v`$, the cyclotron radius is $`R_c=v/\omega _c=cp/eB`$ where $`\omega _c=eB/mc,`$ and $`p=mv.`$ The momentum $`p=\mathrm{}k=h/\lambda `$ is quantally related to wavenumber $`k`$ and wavelength $`\lambda `$. We define the classical small parameter $`ϵ=L/R_c=eBL/\mathrm{}ck=2\pi \varphi /\varphi _0kL`$. Here $`L`$ is the length of the side of the square, $`\varphi `$ is the magnetic flux $`BL^2`$ and $`\varphi _0`$is the flux quantum $`hc/e.`$ Small $`ϵ`$ allows us to approximate orbits within the square as straight lines, to first approximation. This is sometimes known as the Aharonov-Bohm regime, since the leading quantum effects come from the phase interference effects associated with the vector potential, and do not depend on the change of classical orbital motion caused by the Lorentz force. Many potential experiments are in this parameter range.

We choose units such that the dimensionless field is $`B2\pi \varphi /\varphi _0,`$ i.e., $`2\pi `$ times the number of flux quanta in the square. We take $`L,\mathrm{}`$ and $`2m`$to be unity so that

$$ϵB/k<<1.$$
(1)
The dimensionless wavenumber $`k`$ is the number of wavelengths in a side of the square, up to a factor of $`2\pi .`$ It satisfies

$$k>>1,$$
(2)
which is the basis for the quasiclassical approximation.

We shall see that the condition for standard quantum perturbation theory to work is

$$k\sqrt{ϵ}<<1,$$
(3)
or, in other words, $`\sqrt{kB}<<1.`$ This is not completely obvious, and in other contexts, it has been guessed incorrectly that the quantum ‘perturbation border’ is $`kϵ<<1,`$ i.e. $`B<<1,`$ which has the simple meaning that the number of flux quanta in the square is small. We, however, find that nothing much changes at the border $`B1.`$

There may also be a ‘high energy’ condition in the form of a requirement that $`kϵ^b<<1.`$ The exponent $`b`$ depends on the smoothness of the perturbation. We find that for the uniform field, $`b=2`$, while for the ideal flux line $`b=1.`$ Note that for fixed $`B`$ the energy can be arbitrarily high, but if instead $`ϵ`$ is kept fixed, there is a limitation on the energy.

This high energy condition is basically the requirement that diffraction effects not be too important. From a semiclassical perspective, diffraction effects occur where the classical system has a length scale as short or shorter than $`\lambda .`$ An ideal flux line obviously gives rises to diffraction effects. Billiards also have such a length scale of course, namely the distance it takes the confining potential to change from zero inside the billiard, to infinity outside the billiard. This can be taken into account by a ‘Maslov phase’ of $`\pi `$ at the boundary, however. There are also the sharp corners of the square. The square corners whose angular opening is $`\pi /N,`$ where $`N=2`$ is an integer, are a special case at which no diffraction occurs. With such an angle, the billiard can be extended by reflection, and the corner in effect disappears. However, in the presence of a magnetic field, this reflection technique does not work, and with sufficiently large field, orbits which hit directly into the corner eventually become important to the semiclassics.

## III Bogomolny’s Quasiclassical Surface of Section Method

Our approach utilizes the quasiclassical surface of section \[SS\] method of Bogomolny. Poincaré’s surface of section is a surface in classical phase space through which all interesting orbits repeatedly pass. For two dimensional systems, the surface of section is a two dimensional phase space. For a billiard, the Birkhoff surface of section is often chosen. Namely, the space part of the surface of section represents a point on the boundary at which the orbit bounces and is usually measured by the distance along the boundary of the billiard. The variable conjugate to this is the component of momentum parallel to the boundary at the moment of contact. However, many possible surfaces of section can be considered, and some are more convenient than others.

Bogolmony’s method is a generalization of the “boundary integral method”, applicable for billiards, and based on Birkhoff’s surface of section, to much more general systems and surfaces of section. The boundary integral method introduces an operator $`K(x,x^{},E)`$ and an integral equation $`\psi (x)=𝑑x^{}K(x,x^{};E)\psi (x^{})`$. This exact equation has nontrivial solutions only when $`E`$ is on the spectrum. The SS wavefunction $`\psi (x)`$ is the normal derivative of the full wavefunction, $`\psi (x)=\mathrm{\Psi }(𝐫)/n,`$ when $`𝐫`$ is at the boundary point $`x.`$ We should mention that only recently has the boundary integral method been extended to uniform magnetic fields in the case that $`ϵ`$ is of order unity.

Bogolmony’s operator $`T(x,x^{};E)`$ is basically the quasiclassical approximation to $`K.`$ It thus takes the particle crossing the SS at position $`x^{}`$ to its next crossing at position $`x,`$ all at energy $`E=k^2.`$ The quasiclassical approximation to the spectrum is determined by the existence of solutions of $`T\psi =\psi .`$ If only the spectrum is of interest, as it has been for many authors, the condition may be expressed as $`det\left[1T(E)\right]=0.`$

The operator $`T`$ is given quite generally by

$$T(x,x^{};E)=\left[\frac{1}{2\pi i}\left|\frac{^2S(x,x^{};E)}{xx^{}}\right|\right]^{\frac{1}{2}}\mathrm{exp}\left[iS(x,x^{};E)\right]$$
(4)
where $`S=_x^{}^x𝐩𝑑𝐫`$ is the action integral along the classical path from $`x^{}`$ to $`x`$. Note that, rather than giving position and momentum on the SS, positions at two sequential crossings of the SS are given. It is assumed, for notational convenience, that there is a unique orbit from $`x^{}`$ to $`x.`$ We also suppress the Maslov phase. Note that $`T`$ is semiclassically unitary.

The classical action $`S(x,x^{})`$ generates the surface of section map. Namely, the momenta conjugate to $`x,x^{}`$ are given by

$$p=\frac{S(x,x^{})}{x};p^{}=\frac{S(x,x^{})}{x^{}}.$$
(5)
Eq. (5) implicitly gives the surface of section map $`\{p,x\}=\{p^{},x^{}\}.`$

We attempt to simplify $`T`$ by astute choice of the surface of section. It seems simpler to use just one side of the square, rather than all four sides. It is even easier to use a method of images. Namely, we consider, instead of a unit square, $`x,y[\frac{1}{2},\frac{1}{2}][\frac{1}{2},\frac{1}{2}],`$ an infinite channel of width 2 obtained by reflecting the original square first about $`x=\frac{1}{2}`$ and then about $`y=\frac{1}{2},`$ and finally repeating the resulting $`2\times 2`$ square periodically to $`x=\pm \mathrm{}`$. The flux changes sign in neighboring squares. This geometry is shown in Fig. 1. There are a continuum of channel solutions. The solutions to the original square are a subset of these, which exist only at certain quantized energies. This quantization can be carried out in several ways, one of which is shown below.

The SS is taken as the axis $`y=\frac{1}{2}`$ which is identified with $`y=\frac{3}{2}.`$ Because the field is classically weak, the path used to calculate the action is approximated by a straight line. We immediately find

$$S(x,x^{})=k\sqrt{4+\left(xx^{}\right)^2}+\mathrm{\Phi }(x,x^{}),$$
(6)
i. e. the flux free result plus $`\mathrm{\Phi }=(e/c)𝐀𝑑𝐫`$, where the integral is done along the straight line path.

Our scheme finds solutions of $`T\psi =\psi `$ by a perturbation theory. Taking advantage of the fact that $`T`$ is unitary, we first solve

$$𝑑x^{}T(x,x^{};E)\psi (x^{})=e^{i\omega (k)}\psi (x),$$
(7)
treating $`k=\sqrt{E}`$ as a parameter, and then find the energies by solving $`\omega (k)=2\pi n.`$ Given $`\psi (x),`$ a quadrature which can be carried out quasiclassically yields the full wave function $`\mathrm{\Psi }(x,y).`$ Details are given in Appendix A.

## IV Resonances

Classical and quantum perturbation expansions in powers of $`ϵ`$ fail near resonances, necessitating modifications which introduce $`\sqrt{ϵ}`$. Classical resonances correspond to periodic orbits of the unperturbed system. Periodic orbits on the square correspond to straight line orbits in the channel from ($`x^{},\frac{1}{2}`$) to ($`x=x^{}+2p/q,`$ $`\frac{3}{2}`$). Here $`q`$ is a positive integer and $`p`$ is a positive or negative integer relatively prime to $`q.`$ Negative and positive $`p`$ are not equivalent if there is a magnetic flux. Between resonances, or near resonances with large $`p`$ and $`q,`$ ordinary perturbation theory works. See Ref.  for a fuller discussion.

We now specialize to the $`(\pm 1,1)`$ resonances. These are the simplest resonances depending strongly on the field and, as we shall see, in some sense dominate the magnetic response.

We look for a solution of Eq. (7) of the form $`\psi (x)=e^{i\kappa x}u_m(x)`$ where $`\kappa =k\mathrm{cos}45^{}=k/\sqrt{2}`$ and $`u_m`$ varies much more slowly than the exponential. The reason for this choice is that the phase factor $`e^{i\kappa x}`$ makes the rapidly varying phases in the integral $`𝑑x^{}T\psi `$ to be stationary at $`x^{}=x2.`$ This corresponds to rectangular shaped periodic orbits of the original square whose sides make angles of $`45^{}`$ with the $`x`$ axis. Such orbits are shown in Fig. 1.

Because $`ϵ`$ is small, the phase $`\mathrm{\Phi }`$does not greatly influence the position of the stationary phase and it suffices to evaluate $`\mathrm{\Phi }(x,x^{})`$ at $`\mathrm{\Phi }(x,x2)=\mathrm{\Phi }(x+2,x).`$ \[The accuracy of this approximation depends on the smoothness of $`\mathrm{\Phi },`$ and failure of the approximation is related to the onset of diffraction effects mentioned earlier.\] $`\mathrm{\Phi }(x+2,x)`$ is obtained by integrating the vector potential about the closed rectangular loop, and somewhat remarkably is independent of gauge. The result is that $`u_m`$ satisfies the Schrödinger equation

$$u_m^{\prime \prime }+V(x)u_m=E_mu_m$$
(8)
where $`V(x)=k\mathrm{\Phi }(x+2,x)/`$ and $`=\sqrt{8}`$ is the length of the periodic orbits. Thus we convert the phase $`\mathrm{\Phi }`$ to a ‘potential’ $`V.`$ The transverse energy $`E_m`$ enters into the total energy of the eigenstate, according to Eq. (17) below, and $`m`$ is one of two quantum numbers classifying the states.

## V Transverse potential and periodic orbits

For a given resonance and surface of section, there is an effective potential $`V`$ which determines the functions $`u_m`$ and the energy $`E_m.`$ The resonance classically corresponds to a continuous set of nonisolated periodic orbits of the integrable problem. Before perturbation, each of these orbits has the same action. The potential $`V`$, to leading order, is proportional to the change of the action under perturbation, calculated along the unperturbed path. Each such path is labelled by the parameter $`x`$, where it crosses the surface of section. Higher order corrections may also be found.

Knowledge of the potential gives much qualitative insight into the problem. Its minima, \[if smooth\], are at stable periodic orbits, as a rule, and its maxima are at the unstable orbits. In that sense, it represents a classical island chain. Of course, it was known how to quantize states near the stable periodic orbits, if a harmonic expansion is allowed. However, the states that can be found with the aid of $`V`$ are much more general and in particular the states with energies $`E_m`$ near or even above the maxima of the potential can also be found.

In general, isolated unstable periodic orbits do not support wavefunctions, but rather ‘scar’ them. In other words, there appears some excess weight on the wavefunction near the unstable orbit. In the sense of Feynman’s path integral formulation, there are not enough classical paths ‘near’ the unstable orbit, to build a complete wavefunction. Here ‘near’ means that the paths are close to the periodic orbit in the sense of being well approximated by a quadratic expansion about the periodic orbit.

The same is true in the present case, and wavefunctions cannot be built just from orbits near an unstable periodic orbit. However, because of the small parameter $`ϵ`$, we can approximate well an entire shell of orbits in the Feynman integral, and express the result in terms of the potential $`V(x).`$ This shell can support many states which we find. There are states whose energies are near the maxima of $`V`$ and thus have extra weight near the unstable periodic orbits.

The interpretation of $`u_m`$ is that it gives the structure of the wavefunction ‘transverse’ to the resonant periodic orbits. Along the periodic orbits, the wavefunction varies rapidly, but transversely, it varies relatively slowly. The ‘longitudinal’ and transverse motions are weakly coupled, because $`V`$ and thus $`u_m`$ and $`E_m`$ depend on $`k,`$ but this is easy to take into account.

The concept of ‘transverse’ is a little murky in the quantum case, although there are cases, including the one under study where it can be made more precise. We shall not dwell on this further in this paper, however.

We also remark that the Schrödinger equation (8) requires boundary conditions, in order to pick out the physically interesting solutions. These boundary conditions come from the properties imposed on the solution by the physics of the problem, and are usually simplified by symmetries of the problem.

## VI Uniform field solution

### A Effective potential

For the uniform field, the potential is

$`V(x)`$ $`=`$ $`Bk\left(\frac{1}{2}2x^2\right)/;x[\frac{1}{2},\frac{1}{2}],`$ (9)
$`V(x)`$ $`=`$ $`+Bk\left[\frac{1}{2}2(x+1)^2\right]/;x[\frac{3}{2},\frac{1}{2}],`$ (10)
$`V(x)`$ $`=`$ $`V(x+2).`$ (11)
The factor $`\frac{1}{2}2x^2`$ is simply the area enclosed by a periodic $`(1,1)`$ resonant orbit in the shape of a rectangle which bounces from the bottom of the square at $`x.`$ In Appendix B we give the corresponding potential for other resonances.

This periodic potential consists of alternating positive and negative harmonic potential wells of depth $`Bk/2.`$ At the boundaries $`x=\pm \frac{1}{2},`$ the second derivative of the potential is discontinuous, a fact which leads to the mentioned diffraction effects at sufficiently large $`B^2/k.`$

For the $`(1,1)`$ resonance, whose orbits are time reversed $`(1,1)`$ orbits, $`V(x)`$ changes sign. This would not be true if $`V`$ had its origin in a time reversal invariant perturbation of the square, for example, a small change of shape. We can include the $`(1,1)`$ resonance in the present scheme by attributing the region $`1/2<x<3/2`$ to that resonance. This extension of the $`x`$ coordinate is thus similar to use of a ‘angle’ variable, with positive $`x`$-velocity $`v_x`$ for $`x[\frac{1}{2},\frac{1}{2}],`$ and negative $`v_x`$ for $`x[\frac{1}{2},\frac{3}{2}].`$

If $`\sqrt{Bk}=k\sqrt{ϵ\text{ }}`$ is small, the potential $`V(x)`$ can be treated perturbatively. On the other hand, for sufficiently large $`Bk/`$, Eq. (8) will have low lying tight binding harmonic oscillator type solutions centered at $`x=0,`$ (if $`B>0`$), with energies approximately given by

$$E_m=\frac{1}{2}Bk/+(m+\frac{1}{2})\sqrt{8Bk/}.$$
(12)
This formula holds for $`m<<\sqrt{Bk/}`$. The lowest wavefunction is approximately $`u_0(x)=e^{\sqrt{Bk/2}x^2}`$ which is arbitrarily narrow at large energy. These states are paramagnetic, as follows from the fact that $`E_m/B<0.`$ This will be seen more clearly below.

Eq. (8) is valid for larger $`m.`$ Although very simple analytic answers are not available, the problem is the well known one of a particle in a one dimensional periodic potential. We shall see below that we need only consider the boundary conditions $`u(x+2)=\pm u(x).`$ This simplification is a consequence of the symmetry of the square, and something slightly more complicated would be needed for the rectangle.

The solutions to Eq. (8) may be put into four classes, $`A,B,C,D`$. Class $`A`$ states are those with ‘low’ energies near the bottom of the well, $`E_m\frac{1}{2}Bk/`$. For these cases $`u(x)`$ has support only near $`x=0`$, $`\pm 2,\pm 4,\mathrm{}`$ These localized states are strongly paramagnetic, that is, the current circulates in the opposite direction from that of the particle in free space. In this case $`dE_m/dB<0.`$ \[We shall see that the transverse energy $`E_m`$ carries nearly all the field dependence of the total energy of the corresponding two dimensional eigenstates.\]

Class $`D`$ states have energies much greater than the maximum potential energy, i.e. $`E_m>>\frac{1}{2}Bk/`$. In this case, the magnetic field is a small perturbation, since the ‘potential energy’ $`V(x)`$ in Eq. (8) is small compared with the ‘kinetic energy’ given by $`u^{\prime \prime }.`$ These states are weakly diamagnetic. We shall not consider further this case. Of course, the approximation of expanding about the $`(1,1)`$ resonance eventually breaks down, and higher order resonances are eventually involved.

Class $`C`$ has total transverse energy near the top of the potential $`V(\pm 1),`$ that is, $`E_m\frac{1}{2}Bk/`$ and the states are strongly affected by the magnetic field. Very crudely, they are somewhat localized or ‘scarred’ near $`x=\pm 1,`$ since they spend more time in that region. This means that they are strongly influenced by the $`(1,1)`$ resonance. They are diamagnetic and $`dE_m/dB>0.`$

States of class $`B`$ form a transition region between the low-energy paramagnetic states, and the higher energy diamagnetic ones, i.e. near $`E_m0`$. These are states strongly affected by the field, but are such that $`dE_m/dB0.`$

### B Quantization

There are two states with identical energy in the repeated square scheme. These are $`\psi _I=e^{i\kappa x}u_m(x),`$ and $`\psi _{II}=e^{i\kappa x}u_m(x1).`$ \[Changing the sign of $`\kappa `$ is equivalent to changing the sign of the field, which in turn can be accomplished by replacing $`V(x)`$ by $`V(x+1).`$\]

Rather than finding the eigenvalues by imposing conditions directly on the $`\psi `$’s, as in the Appendix A, we produce the four two-dimensional solutions $`\mathrm{\Psi }`$ corresponding to $`\psi _{I,II}.`$ \[Each $`\psi (x)`$ gives two $`\mathrm{\Psi }(x,y)`$’s because $`y`$ and $`1y`$ in the strip represent the same point in the original square.\] We show elsewhere that these states can also be found directly by a Born-Oppenheimer approximation. One of these states may be written

$$\mathrm{\Psi }_0(x,y)=e^{i\kappa (x+y)}u_m(xy\frac{1}{2})$$
(13)
The remaining three states, $`1,2,3,`$ can be obtained by rotations, e.g. $`\mathrm{\Psi }_1(x,y)=\mathrm{\Psi }_0(x,y)=\mathrm{\Psi }_0(y,x),`$ etc. Here $`:(x,y)(y,x)`$ is the rotation by $`90^{}.`$ The gauge can be chosen so that the Hamiltonian is invariant under $``$. Therefore, the symmetry of an eigenstate can be labelled by $`r=0,1,2,3`$, where the eigenvalue of $``$ is $`i^r.`$

Thus, an eigenfunction with symmetry $`r`$ is given by

$$\mathrm{\Psi }_{(r)}(x,y)=\left(\underset{s=0}{\overset{3}{}}i^{rs}^s\right)e^{i\kappa (x+y)}u_m(xy\frac{1}{2}).$$
(14)
\[The sequence of rapidly varying phase factors is $`\{e^{i\kappa (x+y)},e^{i\kappa (x+y)},e^{i\kappa (x+y)},e^{i\kappa (xy)}\}.`$ These in turn are rapidly varying in the $`45^{}`$ directions of the sides of the periodic orbits.\] In general, a solution of Eq. (8) satisfies the boundary condition $`u_m(x+2)=e^{i\beta }u_m(x).`$ We need to find the allowed values for $`\beta `$ and $`\kappa `$ which will give the quantized energies. These conditions are obtained by requiring $`\mathrm{\Psi }_{(r)}(x,\frac{1}{2})`$ to vanish, corresponding to Dirichlet conditions in the original problem. If the wavefunction vanishes on the bottom, it will by symmetry vanish on the boundary of the square.

Clearly, the sum of the two terms \[$`s=0,3`$\] in Eq. (14) which are proportional to $`e^{+i\kappa x}`$ must vanish. This implies $`u_m(x)=i^{3r}e^{i\kappa }u_m(x).`$ The reflection symmetry $`V(x)=V(x)`$ allows us to take $`u_m(x)=(1)^mu_m(x)`$. In turn, this allows quantization of $`\kappa `$ in the form $`\kappa =n\pi /2,`$where $`n`$ is an integer satisfying certain conditions depending on $`r`$ and $`m.`$ Similarly the two terms proportional to $`e^{i\kappa x}`$ in Eq. (14) give the condition $`e^{i\beta }=(1)^r.`$ The relationship of $`n`$ to $`r`$ and $`m`$ is

$$nmod4=\left[2(1mmod2)+r\right]mod4.$$
(15)
It is straightforward to find for the eigenwavenumber

$$k_{n,m}=2\pi n/+E_m/k.$$
(16)
Eq. (16) should be solved iteratively. For example, the first approximation replaces the $`k`$ dependence of the term $`E_m/k`$ by $`2\pi n/.`$ Equivalently, the energy

$$E_{n,m}=4\pi ^2n^2/^2+2E_m.$$
(17)
Note that $`E_m`$ depends, relatively weakly, on $`n`$, since the $`k`$ in formula Eq.(12) should be replaced by $`2\pi n/`$. The dependence of the total energy on $`B`$comes through the term $`E_m.`$ Eqs. (16) and (17) hold for all symmetries and successive values of $`n`$ at fixed $`m`$ cycle through the representations of $``$. Note that, since $`E_m/k<<k`$ the wavelength is given approximately by $`/n,`$ i.e. the length of the classical orbits is an integer number of wavelengths.

Thus we have an expression for the energies of a class of states, namely the $`(\pm 1,1)`$ resonant states. They are labelled by integer $`n`$ which effectively give the number of wavelengths measured along the $`(1,1)`$ periodic orbits, and by a second integer $`m`$ which gives the number of ‘nodes’ ‘perpendicular’ to this orbit. The very low $`m`$ states could very well have been found by earlier methods, since they can be obtained by expansions about the stable periodic orbits. However, these remarkable states do not seem to have been noticed heretofore.

### C Orbital magnetism

The $`(1,1)`$ states just obtained dominate the magnetic orbital susceptibility in a parameter range appropriate to experiments. The susceptibility for the square is on a scale rather larger than the Landau diamagnetism. It is of course not necessary to find the states, or for that matter, their energies, to calculate the susceptibility. That is because the susceptibility depends only on the density of states smoothed over an energy width proportional to the temperature. The Gutzwiller or better, the perturbed Berry-Tabor trace formula is designed to give exactly that quantity in quasiclassical approximation. Nevertheless, it’s interesting and previously unremarked, that a small subset of states accounts for most of the magnetism.

We start by finding the orbital susceptibility $`\chi `$ of a system of noninteracting electrons in a grand canonical ensemble. This is given by $`\chi =/B`$ where the magnetization $`=\mathrm{\Omega }(T,\mu ,B)/B.`$ Here the grand potential is

$$\mathrm{\Omega }(T,\mu ,B)=k_BT\underset{a}{}\mathrm{ln}\left[1+e^{(E_a\mu )/k_BT}\right].$$
(18)
The temperature is $`T`$, $`k_B`$ is Boltzmann’s constant, and $`\mu =k_F^2`$ is the chemical potential. The dependence of $`\mathrm{\Omega }`$ on $`B`$ comes only because the eigenenergies $`E_a`$ depend on $`B.`$ The sum is over all eigenstates labelled by $`a`$.

We divide the states $`a`$ into those relatively few whose energies depend appreciably on the field and the rest. These field dependent states are exactly the $`(1,1)`$ states found above, plus possibly states classically associated with a few other low resonances, e.g. $`(1,3)`$. The reason for this is that the $`(1,1)`$ states enclose the maximum directed area. They also have the shortest length $``$ which we will see plays a role. The even shorter $`(0,1)`$ periodic orbits do not enclose any flux in the approximation of neglecting the curvature of the orbits although at higher fields they eventually become important.

Thus, we replace $`_a=_b+_{n,m}`$ and we can neglect the sum $`b`$ over field independent states. The second sum, over $`(1,1)`$ resonance states, has many fewer terms than the first in a given range of energy. Since $`\mu `$ is related to the number of particles, it is nearly independent of $`B.`$ It is possible to find $`\mu =k_F^2+\delta \mu (B)`$ and make a consistent expansion, and that is indeed necessary if an average over a large number of squares with canonical statistics is done. However, we just want to illustrate how the $`(1,1)`$ states dominate the susceptibility, and we will not consider this further average. Then, we may approximate

$$(T,\mu ,B)=\underset{n,m}{}\frac{E_{n,m}}{B}f_D\left[E_{n,m}(B)\right].$$
(19)
and $`f_D`$ is the Fermi-Dirac distribution function.

Using the Poisson sum formula, replace the sum on $`n`$ in Eq. (19) by an integral over $`k`$, and do the integral to obtain

$``$ $`=`$ $`{\displaystyle \frac{k_BT}{k_F}}{\displaystyle \underset{r,m,s=0}{\overset{\mathrm{}}{}}}\alpha _m\mathrm{exp}\left({\displaystyle \frac{\omega _rs}{2k_F}}\right)`$ (21)
$`\times \mathrm{sin}\left[s\left(k_F{\displaystyle \frac{E_m}{k_F}}\right)\right].`$
Here, $`\omega _r=\pi (2r+1)k_BT`$ is the Matsubara frequency and $`\alpha _m=E_m/B.`$ \[We have dropped the ‘leading’ term in the Poisson formula which totally neglects the discrete quantum nature of the states and which therefore cannot produce a magnetization.\] As an example, take $`k_BT`$ ten times the level spacing of all levels, i.e. $`k_BT=20\pi `$ in our units. Then, $`\omega _0/2300.`$ If $`k_F300600,`$ so that the square contains about 2-6$`\times `$10<sup>4</sup> electrons, the exponential suppression will not be too serious for $`r=0,s=1`$. However, larger $`r`$ or $`s`$ do not contribute much. \[In the trace formula approach, $`s`$ gives the number of repetitions of the primitive periodic orbit and the sum over $`r`$ is explicitly carried out.\]

Eq. (21) shows that states with larger $`,`$ i.e. smaller spacing, are suppressed, exactly as seen from the trace formula in terms of periodic orbits. It also shows that relatively large field dependence of the levels, $`\alpha _m,`$ is important. For the square, the $`(1,1)`$ states have the smallest $``$ and also the largest $`\alpha _m`$. The $`(2,1)`$ resonance does not couple to a small constant field.

It is also seen that if Eq. (21) is averaged over many squares of somewhat different sizes, because of the oscillations of the sine, the result is much reduced and it is necessary to go to higher order in $`\delta \mu .`$

In Fig. 2 we show $`E_m/\widehat{B}`$ as a function of $`E_m/\widehat{B}`$ for the (1,1) resonant states, where $`\widehat{B}=Bk/2`$. According to Eq. (21), if the sign of $`E_m/\widehat{B}`$ is negative the contribution of the corresponding state is paramagnetic. For an exact, two dimensional wavefunction $`\mathrm{\Psi }_{n,m}`$ with energy $`E_{n,m}`$, it is known that $`E_{n,m}/\widehat{B}`$ $`=(/k)d^2r[𝐫\times 𝐣(x,y)]_z`$ , i.e. the expectation value of the $`z`$-component of the magnetization density. Here the current is

$$𝐣(x,y)=2Re\mathrm{\Psi }^{}(x,y)\left(\frac{1}{i}𝐀(x,y)\right)\mathrm{\Psi }(x,y).$$
(22)
In our approximation, according to Eq. (17) above, $`E_{n,m}/\widehat{B}`$ is equal to $`2E_m/\widehat{B}.`$ It follows from Eqs. (8) and (11) above, that $`E_m/\widehat{B}=\widehat{V}(x)_m=u_m\left|\widehat{V}(x)\right|u_m`$ where $`\widehat{V}=V/\widehat{B}.`$ Since $`\widehat{V}<0`$ for $`x[\frac{1}{2},\frac{1}{2}],`$ and $`\widehat{V}>0`$ for $`x[\frac{1}{2},\frac{3}{2}],`$ etc. we see that the sign of the magnetic response of a given wavefunction depends on which region of $`\widehat{V}`$ dominates the expectation value. Of course, the classical periodic orbits in these two regions have the expected sense.

Finally, we may express $`E_m/\widehat{B}`$ quasiclassically as

$$\frac{E_m}{\widehat{B}}=\frac{𝑑x\widehat{V}(x)\left[\widehat{E}\widehat{V}(x)\right]^{\frac{1}{2}}}{𝑑x\left[\widehat{E}\widehat{V}(x)\right]^{\frac{1}{2}}}.$$
(23)
Here $`\widehat{E}=E_m/\widehat{B}`$ can be treated as a continuous variable, so that $`E_m/\widehat{B}`$ as a function of $`\widehat{E}`$ falls on a continuous curve which in this approximation is independent of $`\widehat{B}.`$ The integrals are between the turning points. For $`\widehat{E}`$ near the minimum $`\widehat{V}_{\mathrm{min}}`$ of $`\widehat{V}`$ , $`E_m/\widehat{B}\widehat{V}_{\mathrm{min}}.`$ For $`\widehat{E}`$ at the maximum of $`\widehat{V}`$, the integrals diverge at $`x=1`$, thus making $`E_m/\widehat{B}`$ positive, and also giving the cusp in Fig. 2. Of course, our simple quasiclassical approximation needs corrections in this case.

The quantization of the transverse motion relates $`\widehat{E}=E_m/\widehat{B}`$and $`m`$ according to the standard formula

$$𝑑x\sqrt{\widehat{E}\widehat{V}(x)}=\pi (m+\frac{1}{2})/\sqrt{\widehat{B}}.$$
(24)
The spacing of the quantum states therefore depends on $`\widehat{B}`$. The level $`m_{top}`$ such that $`E_{m_{top}}V(1)`$ satisfies, approximately, $`m_{top}0.6\sqrt{\widehat{B}}.`$ The states with $`k=142`$ and $`B=25`$ were chosen in part because with these parameters there is a state very near the diamagnetic maximum $`m_{top}=14,`$ and another state with $`m=10`$ for which $`E_m/\widehat{B}`$ is very small. On the other hand, with these parameters, while our approximate states are very good for small $`m,`$ the transverse momenta are becoming large enough (as compared with the longitudinal momenta) that our approximation deteriorates significantly at $`m10,`$ especially in regions of the square where the wavefunction is small.

### D Visual representations of eigenstates

Eigenstates of Hamiltonians not satisfying time reversal invariance are rarely shown graphically in the literature, except for trivial cases. In the invariant case, numerically produced picture galleries of eigenstates for systems such as the Bunimovich stadium, have led to a great deal of interest and insight, both theoretical and experimental.

‘Magnetic’ states are not real but are inherently complex. A complete graphical picture of a complex state would seem to require twice the number of pictures as that necessary for a real state satisfying time reversal invariance. In addition, the states themselves are gauge dependent, and only gauge invariant quantities are physically meaningful.

A difference of the magnetic case as compared with the time reversal invariant case is that generally $`\mathrm{\Psi }`$ or $`𝐣`$ do not vanish along nodal lines, but rather only at isolated nodal points. There may of course be symmetries or boundary conditions requiring these quantities to vanish along a line, but in general the representation of wavefunctions by their nodal patterns is not available in the absence of time reversal invariance.

Two gauge invariant quantities we choose to display are the absolute value squared of the wavefunction $`\left|\mathrm{\Psi }(x,y)\right|^2`$ and the current density. The current density is a two dimensional vector field that is divergence free, $`𝐣=0.`$ The one-dimensional surface of section states are also of interest.

We first show a picture of a class $`A`$ state, which is quite simple to represent. Fig. 3 shows $`\left|\mathrm{\Psi }(x,y)\right|^2,`$ for $`n=62,k2\pi 62/140,B=31.4,`$ $`\sqrt{Bk}66,`$ and $`m=0,`$ which implies the symmetry $`r=0.`$ For such a well localized $`u_m`$, each term in Eq. (14) dominates one side of the rectangular periodic orbit. For example, near $`x=y=\frac{1}{4}`$ only the first term, $`s=0,`$ in the sum (14) is appreciable. In this region $`\left|\mathrm{\Psi }(x,y)\right|`$ $`u_0(xy\frac{1}{2})`$ which is well approximated by a Gaussian.

There is interference near the square edges (e.g. near ($`0,\frac{1}{2}`$)), and two terms contribute appreciably. Near this point then

$`\mathrm{\Psi }(x,y)`$ (25)
$``$ $`\left|e^{\frac{1}{2}n\pi iy}u_0(xy\frac{1}{2})+e^{\frac{1}{2}n\pi iy}u_0(xy\frac{1}{2})\right|`$ (26)
$``$ $`\left|2u_0(x)\mathrm{cos}\frac{1}{2}n\pi y\right|`$ (27)
so that $`\left|\mathrm{\Psi }(0,\frac{1}{2}+\frac{1}{n})\right|`$ at its first maximum near $`(0,\frac{1}{2})`$ is about twice as large as $`\left|\mathrm{\Psi }(\frac{1}{4},\frac{1}{4})\right|`$. \[Note that $`nmod4=2.`$\] The current density is thus largest close to the middle of the square edges and it is of course nearly parallel to the edge there. In this case, the shape of $`\mathrm{\Psi }`$ close to an edge is given by $`u_m(x),`$ as can be seen in Fig. 3.

Fig. 4 shows the streamlines of the current in this state. The direction of the streamline gives the direction of the current flow, while the density of streamlines is proportional to the magnitude of the current density. That is, between any two neighboring streamlines, the same total current flows. The state in Fig. 4 is paramagnetic, that is, the current circulates in the opposite sense from that of a free particle in the field. The choice of which streamline to display is easily obtained in this case, since each line crosses a symmetry line like $`x=0,y[\frac{1}{2},0],`$ once and only once. The current has but one interior zero, at the center of the square.

Fig. 5. shows the state $`n=60,`$ $`B=31.4,`$ $`m=1`$, which has one ‘transverse node’. This node does not give rise to a nodal line in the total wavefunction, of course, although the wavefunction is small in regions corresponding to the node. The states have rotational symmetry, and we show a different representation in each corner of the unit square.

One remarkable feature of these states, $`m=0,1`$ is that they are localized very near the central paramagnetic diamond orbit. Theoretically, for large $`Bk`$ this localization can be as tight as one pleases. Although this is a result of our theory which uses a ‘potential’ function $`V(x),`$ it is clear that $`V(x)`$ is really something that arises from Aharonov-Bohm phases, and there is no classical localization based on energy considerations. Although a detailed analysis of this can be a little tricky, we find the same kind of effect in the AB flux line case, for large $`Bk`$ and small $`B.`$ This localization is thus like Anderson localization in the sense that, absent phase interference, localization would not exist. Of course, the random disorder aspects of Anderson localization are absent.

The next figures, Figs. 6-11, have $`n=62,B=25.`$ Figs. 6,7 show a state of class $`C`$ corresponding to $`m=14`$ which is energetically at the top of the periodic potential corresponding to the $`(1,1)`$ resonance. This state is diamagnetic with the current circulating in the opposite sense from the states with $`m=0,1.`$ The theoretical and numerical wavefunctions are shown.

The theoretical predictions in this case are considerably less good than for the states with small $`m.`$ First, the transverse wavenumber is no longer quite so small compared with the wavenumber along the path. Second, since according to Eq. (14), as many as four approximately found component states are added, there may be relatively large errors, especially in those parts of the square where destructive interference is important and the final wavefunction is small. Nevertheless, Eq. (14) is a quite good representation of the more accurate numerical results.

In Fig. 8 we show the state $`m=10`$, for which $`E_{10}/B`$ is very small. The streamlines of this state are rather striking. Note that there are large current loops which have opposite magnetic polarity. Again, the approximate wavefunction captures many features of the exact one, although it does not reproduce the finer details. In this case we show also, in Fig. 9, the transverse state $`u_{10}(x)`$ and the normal derivative of $`\mathrm{\Psi }_{62,10}`$ on the surface of section. Although for small $`m,`$ the normal derivative on the surface of section bears an understandable relation to $`u_m,`$ it is quite complicated for transverse energies this large.

Finally, in Fig. 10 we show streamlines for the sequence of states $`m=6,8,10,12.`$ Although there are systematic changes of pattern, we have not tried to rationalize these changes. We conclude that even though the wavefunction of Eq. (14) is fairly simple, it is difficult to foresee interference patterns when four terms are important.

These relatively complicated states are harder to represent adequately. In both theory and numerics the $`𝐣=0`$ character of the current is not exact. Numerically following a streamline, particularly in the neighborhood of a zero of the current is difficult. We therefore imposed the divergence free character by representing $`𝐣(x,y)=\times \widehat{𝐳}\chi (x,y).`$ We calculated $`\chi `$ as a symmetrized integral of $`j_x,`$ where $`j_x`$ was obtained either numerically or theoretically. The streamlines are then contour lines of $`\chi `$.

The diamagnetic states are not so spatially localized as the low $`m`$ states, although they continue to have a sort of localization in ‘momentum’ space, as we show below. More generally, as $`m`$ increases, the states become more delocalized, and eventually become independent of $`B.`$ This means that for larger $`m`$ all four terms in Eq. (14) make comparable contributions at an arbitrary typical point $`x,y`$, whereas in the localized case, only one or two terms contribute. This gives interference oscillations in $`\left|\psi _m\right|`$ near $`\left|x\right|\frac{1}{2}`$ as shown in Fig. 9.

### E Momentum localization

We have assumed that $`n>>m,`$ and that $`u_m`$ is slowly varying compared with $`e^{i\pi nx/2}.`$ We may expand $`u_m=\widehat{u}_{m,l}e^{i\pi lx}`$, where $`l`$ is an integer for $`r`$ even and half odd integer for $`r`$ odd. Also $`\widehat{u}_{m,l}=(1)^m\widehat{u}_{m,l}.`$ The unperturbed states \[$`\mathrm{sin}\pi p(x+\frac{1}{2})\mathrm{sin}\pi q(y+\frac{1}{2})`$\] can be labelled by integers $`p,q`$ with unperturbed energies $`(p^2+q^2),`$ dropping a factor $`\pi ^2.`$

Eq. (14) is a superposition of unperturbed states with quantum numbers $`p=\frac{1}{2}n+l,`$ $`q=\frac{1}{2}nl.`$ In particular, the energies $`p^2+q^2=\frac{1}{2}n^2+2l^2`$ are closer to the base energy $`\frac{1}{2}n^2`$ than to the base energy of the next representation, $`ϵ_{n+1}\frac{1}{2}n^2+n.`$ Of course, if the perturbation is symmetric under rotation, the next base energies coupled are $`ϵ_{n\pm 4}ϵ_n\pm 4n.`$

There are, however, other unperturbed states with $`p^2+q^2ϵ_n.`$ For example, $`7^2+49^2=ϵ_{70}.`$ However, the matrix elements of a smoothly perturbed Hamiltonian, $`_{pq,p^{}q^{}}`$, in the unperturbed basis, are small if $`\left|pp^{}\right|`$ or $`\left|qq^{}\right|`$ is large. Because the perturbation due to a uniform field has a singular third derivative, these perturbations drop off as a power law, and this relatively long range effect in momentum space is the source of the diffraction corrections.

In Fig. 11 we show the magnitude of the amplitudes of the unperturbed states combining to make the state $`n=62,`$ $`B=25,`$ $`m=14.`$ The area of a circle is proportional to the square of the amplitude. Above the main diagonal we show the theoretical result. All amplitudes lie on the line $`31l,`$ $`31+l.`$ Below, we show the result for the numerical wavefunction. The circle $`p^2+q^2=231^2`$ is also shown. Our theory is starting to need corrections for $`m`$ this large.

Thus, an interpretation of our method which yields the $`(1,1)`$ resonance states of Eq. (14) is that we effectively diagonalize the Hamiltonian in a basis restricted to the unperturbed states nearly ‘degenerate’ with $`ϵ_n`$ and close to $`\frac{1}{2}n,\frac{1}{2}n.`$ This is the case for the uniform field, and indeed, we achieve agreement between full numerical diagonalization, diagonalization restricted to ‘degenerate’ states, and the procedure using the solution of the differential equation Eq. (8).

## VII Aharonov-Bohm Flux Line

The above approach can be generalized to deal with nonuniform flux configurations. To get the potential $`V`$ associated with the $`(1,1)`$ resonance, all that is needed is to be able to calculate the flux contained within a $`(1,1)`$ periodic orbit. We consider here the case of the Aharonov-Bohm flux line \[ABFL\]. Some further results are published elsewhere.

While much of the above discussion can be carried over to the ABFL, in the ideal case of a zero radius line there are strong diffraction effects which limit the applicability of our theory. We therefore begin by defining the alternative problem of a finite size flux line or tube. This is more realistic if actual experiments are contemplated.

Let $`\rho `$ give the linear scale of the flux tube. We may think of the flux as uniform inside a tube of this radius, or as being distributed in some way, say as a gaussian, with $`\rho `$ giving the scale of the distribution. \[We actually used a square tube of side $`\rho .`$ Even more accurately, we used four quarter strength tubes symmetrically located which allows a symmetry reduction in the numerics. In the approximation of our theory, this gives essentially the same result as a single circular tube. The field inside a single flux tube is $`B_0=\varphi /4\rho ^2`$ and $`\varphi `$ is the total flux.\] The typical angular deflection suffered by a classical particle traversing this field is

$$\delta \theta (\varphi /\varphi _0)/k\rho .$$
(28)
To avoid diffraction we require $`\delta \theta `$ to be small. This can be achieved, of course, if $`k\rho `$ is large, but that is not necessary. In the numerical work shown, we take $`\varphi /\varphi _0=0.1,`$ and $`\rho =0.01,`$ while $`k140.`$ An alternative and equivalent condition is to insist that, on the appropriately defined average, the terms in the Hamiltonian satisfy $`(eA/c)^2/2m`$ $`<<e𝐩𝐀/mc`$.

The results depend on where the flux line is located. We consider first the case that it is located at $`x=0,`$ $`y=\frac{1}{2}+a`$, where $`0a\frac{1}{2}.`$ Again, we consider states related to the (1,1) resonance. This leads as before to Eq. (8) but now we find a potential

$`V_{AB}(x)`$ $`=`$ $`Bk/;x[a,a],`$ (29)
$`V_{AB}(x)`$ $`=`$ $`+Bk/;x[1a,1+a],`$ (30)
$`V_{AB}(x)`$ $`=`$ $`0;x[a,a][1a,1+a]`$ (31)
for points in $`[\frac{3}{2},\frac{1}{2}],`$ and $`V_{AB}`$ is extended periodically by $`V_{AB}(x+2)=V_{AB}(x).`$ For the flux line, $`B=2\pi \varphi /\varphi _0.`$

Actually, Eq. (31) is for the ideal flux line. The ideal case has step function jumps in the potential which are smoothed out at finite $`\rho `$. In other words, the sharp jump at $`x=a`$, is replaced by a smooth rise beginning at $`x=a\rho `$ and ending at $`x=a+\rho .`$ The exact shape of the rise depends on the distribution of flux in the line. As long as the transverse wavelengths in the solution of Eq. (8) are long compared with $`\rho ,`$ the finiteness of $`\rho `$ does not play a significant role in our theory at this level of approximation.

For the ideal, zero radius ABFL, there are significant deviations from this scenario. Indeed, most matrix elements of the ABFL perturbed Hamiltonian in the unperturbed basis are infinite. However, it is a weak, logarithmic infinity, and our theory seems to capture the main shape of the wavefunction, although at the relatively low energies for which numerical results are available there are significant corrections. We consider these to be diffractive corrections, arising from a characteristic length shorter than the wavelength.

It is clear that $`Bk/`$ is the important parameter in the UF case, while for the ABFL, both $`Bk/`$ and $`a`$ are important. We begin with the case $`a=\frac{1}{4},`$ which has a ‘square well potential’ of width $`\frac{1}{2}`$ near $`x=0.`$ For sufficiently large $`Bk/,`$ there will be ‘tight binding’ solutions approximately $`u_m(x)=\mathrm{cos}(m+1)\pi x/2a,\left|x\right|<a,`$ and zero elsewhere. This expression holds for sufficiently small even $`m,`$ and for odd $`m`$ the cosine is replaced by the sine. The energy $`E_mBk/+\pi ^2(m+1)^2/4a^2.`$

Fig. 12 is for the ABFL case, with $`n=86,`$ $`m=0,`$ $`r=0`$ and $`a=\frac{1}{4}.`$ Fig. 12(a) shows $`u_0(x)`$ and $`u_0(x1)`$, its extension into $`x<\frac{1}{2},`$ reflected. For these parameters, $`u_0`$ is not extremely localized, and extends significantly outside $`[\frac{1}{2},\frac{1}{2}].`$ The remaining plots give $`\left|\psi _m\right|=\left|\mathrm{\Psi }/n\right|.`$ Fig. 12(b) plots Eq. (14) and 12(d) is from diagonalization in the limited basis of the ‘degenerate’ states. Clearly, these two approximations are nearly the same. Fig. 12(c) and 12(e) are obtained by numerical diagonalization in the complete unperturbed basis, for the finite size flux tube and the ideal flux line respectively. Clearly, for these parameters, there is significant diffraction from the flux tube. It is somewhat less than for the zero radius line but the two patterns have some resemblance. Diffraction evidently modifies the interference between different parts of the wavefunction in an irregular way. The overall shape of $`\left|\mathrm{\Psi }/n\right|`$ is well predicted by the theoretical $`u_0`$. To give a sense of the irregularity of diffraction effects on the wave function, we show in Fig. 12(f) the function for $`n=70,`$ a somewhat lower energy. In this case, according to theory, the $`u_0`$ is essentially the same, and the interference fringes have a somewhat longer wavelength, on the average. The diffractive effects are quite different in detail, however.

It is clear from this result that the ABFL localizes the particle, and to do so, it must of course exert a force. That a flux line can exert a transverse force is now well established.

In Fig. 13 we show currents from a state with $`a=\frac{1}{2},`$ i.e. the much studied case with the flux line at the center of the square. Again, $`m=0,`$ $`r=0,`$ and $`n=`$ $`82`$. The upper part of the figure shows numerical streamlines for a corner of the square. Although the state is spatially not well localized, it does have a strongly ‘paramagnetic’ current structure. By paramagnetic we mean $`dE_m/d\varphi <0`$.

The lower part shows the theoretical current density $`j_y(x)`$ for $`y=0,x[\frac{1}{2},0],`$ that is, the current density along the upper edge of the upper figure. A simple approximate formula for $`j_y`$ is $`j_y(x)\left(\mathrm{cos}n\pi x/2\right)^2\left[u_m\left(x\frac{1}{2}\right)^2u_m\left(x\frac{1}{2}\right)^2\right]`$ which gives double zeroes of the current at equally spaced $`x=(2l+1)/n.`$ In this case $`j_y`$ is negative and the factor depending on $`u_m`$ has no zeroes except at $`x=0,`$ for $`m=0.`$

The maxima are at $`j_y=0`$; $`j_y0.`$ Theory and numerics closely agree on the period and shape of the oscillations, but the lower envelope of the two differ. We mark with horizontal bars the lower envelope of the numerical calculation, which is considerably more irregular than given by our theory. This structure of zeroes of $`j_x`$ is a consequence of a symmetry, namely $`yy`$, together with complex conjugation. This means that $`\mathrm{\Psi }(x,0)`$ is real and therefore generically will have zeroes as a function of $`x.`$

## VIII Summary

### A Extensions of the results

We have shown how to classify and find eigenstates of a charged particle in a square billiard subjected to a magnetic flux which is classically weak. It is not necessarily weak quantally, however, and large remarkable changes in the wave functions are found which sometimes significantly localize the wave function. We used two basic flux configurations: a uniform field and an Aharonov-Bohm flux line. The placement of the ABFL is important. Of lesser importance is any finite radius to the line, at least in the range of parameters we use. We exhibited some of the wavefunctions in several forms and sequences, concentrating on wavefunctions connected with the $`(1,1)`$ resonance, that is, connected with periodic orbits of the flux free square whose velocities make angles of 45 with the coordinate axes.

These results can be readily generalized to higher order resonances. For a given energy the effects diminish quite rapidly as the order increases, but in principle, at sufficiently high energy any given resonance can show strong magnetic effects. We give a few results in Appendix B.

One can extend these results to integrable systems other than the square billiard, and to other flux configurations. There are other billiard shapes, such as the rectangle, certain triangles, the circle and the ellipse. One may also study ‘soft billiards’, e.g. a confining potential of the form $`U(x,y)=U_1(x)+U_2(y)`$ where $`U_i`$ is some sort of anharmonic potential. Soft billiards are more accurate representations of mesoscopic systems than are hard wall billiards, but the additional theoretical effort they require is not usually made. A practical difficulty, although not one of principle, is that it may be necessary to resort to action-angle variables, and it may be tedious to find the periodic orbits of the flux-free system, if some tiling trick cannot be used.

In any case, the first steps of the program are to choose a convenient surface of section, choose the resonance of interest, find the AB phase $`\mathrm{\Phi }(x,x^{}),`$ and the effective potential $`V(x).`$ Much insight can be gained at this level.

A much larger class of systems solvable by this technique are nearly integrable systems subjected to a flux. For example, one could start from a nearly square trapezoid. Then one would have a double perturbation of the square, one from the flux, the other from the change of shape. Each of these perturbations can have big quantum effects, and the combination of the two can be quite different from each one separately, especially since one breaks time reversal invariance and the other doesn’t. Again, finding the effective potential is key to understanding the qualitative results.

### B Possible experiments

There are possible experiments and even interesting devices which might be made, if time reversal can be broken. Having unusual wavefunctions suggests some of the possibilities. Probes of the system will depend on whether the wavefunction is large or small at the position of the probe.

Mesoscopic systems are of great current interest. In this case, electrons are the particles and there is a magnetic flux. The system of Lévy et al. consisting of isolated metallic ‘squares’ is a case of this type, but one can imagine leads weakly connected to a metallic square whose coupling depends on the shape of the wavefunction at the Fermi level and the position of the lead.

Another type of mesoscopic system is formed of surface electrons in ‘corrals’ on a metallic surface. The wave functions of such electrons can be probed with atomic accuracy with a scanning tunnelling microscope. Achieving an appropriate parameter range will be difficult, but perhaps not impossible. The corrals are leaky and do not really confine the electrons to their interior, but in many cases that idea seems to work, at least qualitatively.

Another kind of system is the shallow square or rectangular container of liquid, which is vibrated to set up standing waves. Instead of a flux, the system can be rotated with a uniform angular velocity. This differs from the uniform magnetic field in that the $`A^2`$ is absent from the Hamiltonian, but since under our assumptions that term can be neglected anyway, the same kind of results are expected. One can also introduce a nonideal ABFL. This is done by making a small hole in the tank and allowing the water to flow out with a certain vorticity. Experiments on such a system have been performed, although scattering rather than eigenstates was studied. Such an experiment would certainly have pedagogical value, and if the experiment is done, it may in fact be the first detailed observation of a nontrivial persistent current state.

Still another system is the thin square microwave cavity. Here the ‘quantum’ waves are microwaves. There is no quantity directly equivalent to the magnetic field. However, one can replace part of the boundary by a ferrite strip, say between $`[a,a]`$ on one side of the cavity. The phase shift of the microwave upon reflection from the ferrite depends on the direction of magnetization of the ferrite and the direction of incidence of the microwave. This gives as an effective potential exactly that of Eq. (31). There is then a localized eigenstate circulating the cavity in one direction, but not the other. A similar situation from the point of view of diffraction and effective potential, but maintaining time reversal invariance, is the ‘step’ billiard. This is a square billiard of side $`L`$ say, with one side moved down by $`\delta L=ϵL`$ for $`\left|x\right|<a.`$ The perturbed $`(1,0)`$ resonance states ‘see’ a square well effective potential to first approximation. If $`k\delta L`$ is small, but $`k\sqrt{L\delta L}1,`$ our theory works, gives nontrivial localization, and is similar to the flux line with small $`\varphi /\varphi _0.`$ The range $`k\delta L=\pi ,`$ i.e. $`\delta L`$ a half wavelength, gives no phase shift for normal incidence and is similar to an ABFL with $`\varphi /\varphi _0=1.`$ Experiments on this system could be carried out in several contexts.

### C Conclusions

We have solved a characteristic example of quantum states in a weakly perturbed integrable system. The states and their energies can be classified and found to good approximation. The wavefunctions are quite nontrivial and interesting, as compared with the states of the integrable system. The states of a hard chaotic system are also not individually interesting, except for some weak ‘scars’, requiring resort to statistical studies and averages over many states in such systems.

Our technique applies to a very large class of systems, which includes experimental systems. The particular case we emphasize breaks time reversal invariance, and the eigenstates have persistent currents. As far as we know, these are the first published nontrivial examples such states.

## IX Acknowledgments

Supported in part by the United States NSF grant DMR-9625549 and United States-Israel Binational Science Foundation, grant 99800319. R.N. was partially supported by the NSF grant DMR98-70681 and the University of Kentucky. We thank Prof. Director Peter Fulde for hospitality at the Max-Planck-Institut für Physik komplexer Systeme in Dresden, where some of this work was done. R.N. and O.Z. thank Dr. R. Seiler for hospitality at the SFB 288 “Differentialgeometrie und Quantenphysik”, TU Berlin. O.Z. thanks Dr. F. Haake for hospitality at the Universität GH Essen.

## A Two dimensional wavefunctions

In this Appendix we find the two dimensional wavefunction. If a solution to $`\psi =T\psi `$ is known, the wavefunction in the square can be found by evaluating the integral over the surface of section.

$$\mathrm{\Psi }(𝐫)=G(𝐫,x^{};E)\psi (x^{})𝑑x^{}.$$
(A1)
In this case

$`G(𝐫,x^{};E)`$ $`=`$ $`{\displaystyle \underset{\text{cl. tr.}}{}}{\displaystyle \frac{1}{2}}\left|{\displaystyle \frac{1}{2\pi k}}{\displaystyle \frac{^2S}{x^{}\mathrm{}_{}}}\right|^{1/2}`$ (A3)
$`\times \mathrm{exp}\left[iS(𝐫,x^{};E)i{\displaystyle \frac{\pi }{2}}\nu \right]`$
in our units. Here the sum is over the (unperturbed) classical trajectories going from a point $`x^{}`$ on the surface of section to a point $`𝐫`$ inside the square, $`S(𝐫,x^{};E)`$ is the reduced action for this trajectory, $`\mathrm{}_{}`$ is the direction perpendicular to the trajectory at point $`𝐫`$, and $`\nu `$ is the Maslov index. In a single square scheme we take the surface of section to be the lower side of the square $`1/2x^{}1/2`$, $`y^{}=1/2`$. Then $`\psi (x^{})`$ is a linear combination of functions $`\psi _I`$ and $`\psi _{II}`$ described in the text such that $`\psi (\pm 1/2)=0`$ (there should be no flow through the ends). It follows then that $`u_m(x^{}+2)=(1)^ru_m(x^{})`$ and

$$\psi (x^{})=e^{i\kappa x^{}}u_m(x^{})+i^re^{i\kappa x^{}}u_m(x^{}1)$$
(A4)
with $`r`$ defined by Eq. (15). Here we have also used a zero field quantization condition $`\kappa =\pi n/2`$ and the property $`u_m(x^{})=(1)^mu_m(x^{})`$. Note that in the surface of section picture $`r`$ does not have a direct interpretation as a representation label.

We evaluate the integral (A1) in the stationary phase approximation. As we shall see, only the 45 orbits survive. For a given point $`𝐫=(x,y)`$ there are four such orbits: two starting at $`x_1^{}`$ and two at $`x_2^{}`$ (Fig. 14). Each orbit gives one term in Eq. (14). Note that the orbits with the positive (negative) $`x`$-projection of the momentum at $`x^{}`$ are generated by the first (second) term in Eq. (A4).

Although the following calculation depends in which sector inside the square the point lies, the final result is independent of sector. So for definitiveness we assume $`y<x<y`$, that is the point lies below the diagonals $`y=x`$ and $`y=x`$. We have numbered the orbits by $`s=0,\mathrm{},3`$ according to the terms in Eq. (14) they represent. For $`s=0`$ orbit, expressing the action in terms of the distance, $`S(𝐫,x^{};E)=kL(𝐫,x^{})k\sqrt{(xx^{})^2+(y+1/2)^2}`$, we find the prefactor in Eq. (A3)

$$\left|\frac{1}{8\pi k}\frac{^2S}{x^{}\mathrm{}_{}}\right|^{1/2}=\left(\frac{1}{8\sqrt{2}\pi L}\right)^{1/2}.$$
(A5)
The stationary point in the integral (A1) is determined by the exponents in Eq. (A3) and in the first term of Eq. (A4). It is $`x_1^{}=xy1/2`$, as in a 45 orbit. The $`s=0`$ wavefunction is then

$$\mathrm{\Psi }_0(x,y)=e^{i\kappa (x+y)}u_m\left(xy\frac{1}{2}\right)$$
(A6)
dropping a constant factor $`\kappa ^{1/2}e^{i(\kappa /2\pi /4)}/2`$. The terms with $`s=1,2,3`$ can be obtained in the same manner with an appropriate addition of the Maslov phase for each bounce.

If $`𝐫`$ is located in a different sector of the square, the results remain the same if proper account is made of the Maslov phases. The labeling of the orbits is invariant if done by the rule: $`s=0`$ trajectory arrives to $`𝐫`$ from the South-West, $`s=1`$ from SE, $`s=2`$ from NE, and $`s=3`$ from NW.

It is worth noticing that Eqs. (A6) or (14) are valid only up to the square root order of the small parameter $`ϵ=B/k`$ of the perturbation theory behind our work. If we wish to obtain the results valid to the first order of $`ϵ`$ we should (a) use a better approximation for $`u_m(x^{})`$, (b) add the vector potential term to the action $`S(𝐫,x^{};E)`$, and (c) find a correction to the stationary point due to $`u_m(x^{})`$.

Alternatively, the above calculation can be carried out in the repeated square scheme defined in Sec. III. The surface of section is then the line $`y=1/2`$ (identified with $`y=3/2`$). We can restrict ourselves to the trajectories with the positive $`x`$\- and $`y`$-projections of the momentum in this extended manifold. Then only $`\psi _I(x^{})`$ is needed to generate a complete two-dimensional wavefunction. In order for $`\psi _I`$ to have a period 2, the condition $`u_m(x^{}+2)=(1)^ru_m(x^{})`$ must be satisfied. A point $`𝐫`$ in the original square will be represented by its images in four (up to a translation by $`\mathrm{\Delta }x=2`$) domains in the manifold. Each of these images is connected by one 45 orbit to the surface of section. These four orbits generate the four terms $`s=0,\mathrm{},3`$ in Eq. (14) as follows: $`s=0`$ is produced by the point $`𝐫`$ in the original domain, $`s=1,3,2`$ by its reflection about the right side of the square, upper side of the square, and the composition of both, respectively. There are no Maslov phases in Eq. (A3) in this case. However the terms $`s=1,3`$, which are obtained by the odd number of reflections, should be summed with an additional minus sign. Indeed, the Jacobian of the coordinate transformation from the original square to the extended manifold is singular on the boundary. So, when the wavefunction on the manifold is folded back to the physical domain a phase difference $`\pi `$ will be accumulated between each pair of domains related by one reflection. Of course, this phase is analogous to the Maslov phase for an impenetrable wall.

## B Higher resonances

For completeness, we give some results for higher resonances, which can be labelled with relatively prime integers, $`p,q.`$ These resonances are closely related to states of a rectangle of side $`1/p,`$ $`1/q.`$ Obviously, for the square, the energies for $`p,q`$ are the same as those for $`q,p.`$ The corresponding states are also the same after a $`90^{}`$ rotation. However, for $`\left|pq\right|/\left|p+q\right|`$ not too small, the $`p,q`$ states and the $`q,p`$ states are different and are nearly uncoupled and there is a nearly doubly degenerate set of energy levels. The splitting of the exact levels, which are combinations of $`pq`$ and $`qp`$ that are eigenstates of the rotation operator $``$, can be estimated by using an analog of ‘chaos assisted tunnelling’.

The classical periodic orbits make $`q`$bounces from the $`x`$ sides and $`p`$ bounces from the $`y`$ sides. These resonant orbits all have the same length $`_{pq}=2\sqrt{p^2+q^2}.`$ The maximal directed area enclosed by an orbit is $`\pm 1/2pq.`$ The energy is $`k_{nm}^2`$ where

$$k_{nm}^2\left(\frac{2\pi n}{_{pq}}\right)^2+\left[1+\frac{p^2}{q^2}\right]E_m^{(q)}$$
(B1)
and $`n`$ is an integer. The first term is much larger than the second. Again, we have the approximate quantization of $`n`$ wavelengths in an orbit length.

The potential $`V_{pq}(x)`$ depends on $`q`$ being associated with the $`x`$-direction, thus the notation $`E_m^{(q)}`$. This ‘transverse’ contribution to the energy is of course symmetric in $`p`$ and $`q,`$ that is, $`q^2E_m^{(p)}=p^2E_m^{(q)}.`$ The potential is

$$V_{pq}(x)=\frac{q}{p}\left[\frac{_{11}}{_{pq}}\right]^3V\left[q\left(x+\frac{1}{2}\right)\frac{1}{2}\right]$$
(B2)
where $`V`$ is given by Eq. (11). $`V_{pq}`$ has period $`2/q`$, rather than $`2.`$ The boundary condition on the eigenstates is $`u(x2)=e^{i\beta }u(x)`$ where $`\beta =2\pi frac\left[pn/\left(p^2+q^2\right)\right]`$ and $`frac`$ indicates the fractional part. The potential can be regarded as weak if $`kB<<pq_{pq}^3.`$ If this condition is satisfied the $`pq`$ resonance can be ignored.

Similar but somewhat more complex results can be obtained for rectangles and certain triangles.