# The Kepler problem on the lattice

D. Sanjinés, E. Mamani  Carrera de Física, Universidad Mayor de San Andrés, c.
27 Cota - Cota, casilla de correos 8635, La Paz, Bolivia.
email<EMAIL_ADDRESS><EMAIL_ADDRESS>J.
Velasco  Departamento de Física, Universidad Técnica de Oruro, Oruro, Bolivia.
email<EMAIL_ADDRESS>
###### Abstract

We study the motion of a particle in a 3-dimensional lattice in the presence
of a potential $-V_{1}/r$, but we demonstrate semiclassicaly that the
trajectories will always remain in a plane which can be taken as a rectangular
lattice. The Hamiltonian model for this problem is the conservative tight-
binding one with lattice constants $a,b$ and hopping elements $A,B$ in the
$XY$ axes, respectively. We use the semiclassical and quantum formalisms; for
the latter we apply the pseudo-spectral algorithm to integrate the Schrödinger
equation. Since the lattice discrete subspace is not isotropic, the angular
momentum is not conserved, which has interesting consequences as chaotic
trajectories and precession trajectories, similar to the astronomical
precession trajectories due to non-central gravitational forces, notably, the
non-relativistic Mercury’s perihelion precession. Although the elements of the
mass tensor are naturally different in a rectangular lattice, these can be
chosen to be still different in the continuum, which permits to study the
motion with kinetic energies $p_{i}^{2}/2m_{i}$ ($i=x,y$). We calculate also
the contour plots of an initial Gaussian wavepacket as it moves in the lattice
and we propose an “intrinsec angular momentum” $S$ associated to its
asymmetrical deformation, such that the quantum and semiclassical angular
momenta, $L_{q},L_{c}$, respectively, could be related as $L_{q}=L_{c}+\alpha
S$.

###### keywords:

Kepler problem; tight-binding dynamics; semiclassical model; pseudo- spectral
method; physics education.

††journal: Rev. Mex. Phys. E

##  1 Introduction

The Kepler problem on the lattice can be traced back to the anisotropic Kepler
problem about the description of an electron dynamics in the presence of the
Coulomb field of a donor impurity in a Silicon or Germanium semiconductors
[1]. Some authors have revived this issue either in the theoretical or the
experimental contexts, for example, because of more efficient numerical
algorithms and simulation applications, or new and more precise experimental
techniques. We may mention a few authors in the theoretical context: (i) Boris
et al. [2] use quantum analytical and numerical methods to describe the
wavepaket evolution within elliptical and hyperbolic orbits. (ii) Bai and
Zheng [3] invoke classical and quantum mechanics methods for studying the
Kepler problem with a kind of weak anisotropy of the mass tensor. (iii)
Petrova [4] set about solving exactly the full energy spectrum (via Green’s
function methods) of a particle hopping on a Bethe lattice under a Coulomb
potential. (iv) Cortés [5] has particularly caught our attention because he
obtains a chaotic particle’s trajectory by defining an anisotropic mass
tensor, a result that we could also obtain by defining the Kepler problem on a
rectangular lattice, such that in the continuum limit the lattice disappears
but an anisotropic mass tensor remains.

The main motivation for our work presented in this article has an educational
purpose: the Kepler problem is a common topic taught at the undergraduate
level (a usual reference text is Symon [6]) where the conic-section
trajectories are deduced from the conservation of energy and angular momentum.
As pointed out above [5], the definition of an anisotropic mass tensor, and
consequently the non-conservation of angular momentum, leads to a simple
system where chaos was rigorously proved, becoming therefore a common
educational model to teach this subject. Within this same purpose, we realize
that the introduction of a “fake lattice” in the space can be used as a model
to justify the mass asymmetry, not as a consequence of some particle’s
property, but as a property of space itself lacking rotational symmetry.
Moreover, two important and rather different subjects in the undergraduate
physics education, namely, the astronomical Kepler problem (studied in
classical mechanics [6]) and the tight-binding dynamics (studied in solid
state by quantum and semiclassical methods [7]) converge in the “Kepler
problem on the lattice”. This, in turn, is the ocassion to draw other subjects
in physics education, as for example: higher order statistical moments
associated to random distributions [8] and their relation to physical
properties of the quantum probability density (as suggested in this article);
analytical and numerical methods to solve the Schrödinger equation [9] [10].

Our work in this article is organized as follows: In Sections 2 and 3 we set
about developing the semiclassical and quantum formalisms, respectively, at
the common level of the usual textbooks for the undergraduate education in
physics, for example, Symon [6] and Liboff [9]. In Section 4 we present the
central quantum results of our work and compare them to previous semiclassical
results. In Section 5 we give a list of the most relevant results of our work.

##  2 Semiclassical model for the rectangular lattice

The tight-binding Hamiltonian function for an electron moving by the hopping
mechanism in a 3-dimensional orthorhombic lattice with constants $a,b,c$ along
the $X,Y,Z$ axes and hopping elements $A,B,C$, respectively, and in the
presence of a donor impurity potential located at the origin, is:

$\displaystyle H(\textbf{r},\textbf{k})=$ $\displaystyle
2A(1-\cos{ak_{x}})+2B(1-\cos{bk_{y}})$
$\displaystyle+2C(1-\cos{ck_{z}})-V_{1}/\sqrt{x^{2}+y^{2}+z^{2}}.$ (1)

This function has been constructed, in accordance with the semiclassical
method [7], from the spectrum of the tight-binding Hamiltonian operator for a
particle hopping in a periodic potential
$V(\textbf{r}+\textbf{R})=V(\textbf{r})$ with $\textbf{R}=(na,mb,qc)$ ($n,m,q$
integers) and $\hbar\textbf{k}=\hbar(k_{x},k_{y},k_{z})$ the position and
quasimomentum vectors, respectively, in 3D orthorhombic Bravais lattices. In
the limit of the effective mass approximation wherein
$m_{x}=\lim_{A\rightarrow\infty,a\rightarrow 0}\hbar^{2}/(2Aa^{2})$, the terms
$\hbar k_{x}$ and $\hbar^{2}k_{x}^{2}/(2m_{x})$ are the usual momentum and
energy, respectively, of a free particle along the $X$ direction. The term
$V_{1}$ in (1) is the coefficient (in units of joule-meter) of the Coulomb
potential energy due to an external point source (not intrinsic to the
lattice) which is fixed somewhere upon the lattice and referred to in this
article as the impurity. Fig. 1 illustrates a case of the Kepler problem on
the lattice described above.

The corresponding Hamilton equations of motion (with
$\textbf{p}=\hbar\textbf{k}$) are

$\displaystyle\dot{x}=\frac{\partial H}{\partial
p_{x}}=(2Aa/\hbar)\sin{ak_{x}},$ (2) $\displaystyle\dot{y}=\frac{\partial
H}{\partial p_{y}}=(2Bb/\hbar)\sin{bk_{y}},$ (3)
$\displaystyle\dot{z}=\frac{\partial H}{\partial
p_{z}}=(2Cc/\hbar)\sin{ck_{z}},$ (4)
$\displaystyle\hbar\dot{k}_{x}=-\frac{\partial H}{\partial
x}=V_{1}x(x^{2}+y^{2}+z^{2})^{-3/2},$ (5)
$\displaystyle\hbar\dot{k}_{y}=-\frac{\partial H}{\partial
y}=V_{1}y(x^{2}+y^{2}+z^{2})^{-3/2},$ (6)
$\displaystyle\hbar\dot{k}_{z}=-\frac{\partial H}{\partial
z}=V_{1}z(x^{2}+y^{2}+z^{2})^{-3/2}.$ (7)

These equations along with their corresponding initial conditions at $t=0$ for
the position $\textbf{r}(0)$ and velocity $\textbf{v}(0)=\dot{\textbf{r}}(0)$,
constitute the semiclassical model for the tight-binding dynamics of a hopping
particle in a lattice. By choosing the initial conditions such that the
particle moves in the $XY$ plane and is not located at the origin:
$z(0)=\dot{z}(0)=0$, $x^{2}(0)+y^{2}(0)\neq 0$, we deduce from (7) that
$\dot{k}_{z}(0)=0$ and from (4) that $k_{z}(0)=0$; therefore
$\ddot{z}(0)=2Cc^{2}\dot{k}_{z}(0)\cos(ck(0))=0$, i.e., the particle will
continue moving in the $XY$ plane in the abscence of any perpendicular
acceleration that could drag the particle out of the plane. Henceforth, we
will consider the dynamics of the electron in a rectangular lattice described
by Eqs. (2), (3), (5) and (6) above (with $z=0$). Besides, and for simplicity
of notation, we will take a unitary value for the angular Planck constant,
$\hbar=1$. This would imply, of course, that the momentum and energy should be
expressed in units of 1/meter and 1/second, respectively, but the election of
$\hbar=1$ in this article does not bear any specific physical meaning, since
the true numerical value of $\hbar$ would be restored before interpreting
physical measurements.

Figure 1: Case of the Kepler problem on the lattice where the impurity is
located at the origin of the orthorhombic lattice; the electron is initially
located at the position $\textbf{r}_{0}=(0,2b,0)$ subject to the impurity
Coulomb interaction and to the 1st-neighbor interactions with energies
(hopping elements) $A,B,C$ along the $X,Y,Z$ axes, respectively. If the
initial velocity lies on the $XY$ plane, the particles’s trajectory will also
remain in the plane (dotted line).

From the definition of angular momentum,
$\textbf{L}=\textbf{r}\times\textbf{p}$, we have the component of L in the
rectangular lattice as $L_{z}=xk_{y}-yk_{x}$. Its corresponding time-variation
is $\dot{L}_{z}=2\left(aAk_{y}\sin{ak_{x}}-bBk_{x}\sin{bk_{y}}\right)$, which
is different from zero for an arbitrary instant. This means that in the
rectangular lattice the angular momentum $L_{z}$ is not conserved in the
presence of the point impurity located at the origin with a corresponding
Coulomb singular potential $-V_{1}/\sqrt{x^{2}+y^{2}}$, equivalent to the
gravitational potential of a point mass in the two-body Kepler problem.
However, for these problems we know from our background in classical mechanics
that the angular momentum is conserved (along with the total energy) because
the Coulomb and gravitational forces are central and therefore the torque is
zero. What is assumed in these cases is that the space is rotationally
symmetric (isotropic); that is not valid for the rectangular lattice. The
procedure to achieve an isotropic space from the rectangular lattice is,
first, to make a square lattice out of a rectangular one by making $a=b$ and
$A=B$ in the Hamiltonian function $H$ in (1), and second, to take the
continuum limit as $a\rightarrow 0$ while $A\rightarrow\infty$, such that
$1/m=2Aa^{2}$ is defined as the “effective mass”. In fact, since the cosine
approximates as $\cos(ak_{x})\cong 1-a^{2}k_{x}^{2}/2$ for small arguments, we
have that

$\displaystyle
H\left(x,y;k_{x},k_{y}\right)\rightarrow(k_{x}^{2}+k_{y}^{2})/2m-V_{1}/\sqrt{x^{2}+y^{2}},$
(8)

giving thus the usual Hamiltonian function for the 2-dimensional Kepler
problem. We can see that in this case $\dot{L}_{z}\rightarrow 0$, i.e., the
angular momentum is conserved, as expected. Interestingly, if the continuum
limit is obtained from the rectangular lattice by maintaining $a\neq b$ and
$A\neq B$, we obtain

$\displaystyle H\left(x,y;k_{x},k_{y}\right)\rightarrow
k_{x}^{2}/2m_{x}+k_{y}^{2}/2m_{y}-V_{1}/\sqrt{x^{2}+y^{2}},$ (9)

along with $\dot{L}_{z}\neq 0$, which describes the classical anisotropic
Kepler problem reported by Cortés [5] where the angular momentum is not
conserved due to the difference of effective masses $m_{x},m_{y}$ in the $XY$
axes, respectively. As we can see, in this case it is suggested that the
“anisotropy parameter” defined as $\xi\equiv m_{x}/m_{y}$ [5] has its origin
in the quotient $m_{x}/m_{y}=Bb^{2}/(Aa^{2})$ of a subyacent rectangular
lattice in the continuum limit. In order to reduce the free parameters in that
quotient we can invoke an expression reported by Wolf [11] wherein the hopping
element is related to the separation by a lattice constant between 1st.
neighbor sites: $A=2(1+a/a_{0})\exp(-a/a_{0})$ and
$B=2(1+b/a_{0})\exp(-b/a_{0})$, where $a_{0}$ is the Bohr radius.

Figure 2: Trajectory of the particle in the $XY$ plane of a rectangular
lattice with a ratio $a/b\cong 20$ of the lattice constants (solid line). The
initial conditions are $x(0)=0.23a,y(0)=0,k_{x}(0)=0$ and $k_{y}(0)=5.33/b$.
The resulting motion yields an oscillatory trayectory parallel to the $Y$
axis. In the continuum limit $a,b\rightarrow 0$ (maintaining $a/b\cong 20$)
and $A,B\rightarrow\infty$ corresponding to the Hamiltonian (9), and
$k_{y}(0)=400.68/b$, the trajectory becomes the dashed one. This result
reproduces the one reported by Cortés [5] in the continuum, who introduces
different masses $m_{x},m_{y}$ in the $XY$ axes. Our result show that the
rectangular lattice approach to the continuum renders a satisfactory
explanation for the origin of the particle’s mass anisotropy.

In fact, following Cortés [5], we take $\xi=2.94$ as the anisotropy parameter
and solve $\xi=Bb^{2}/(Aa^{2})$ with $a/a_{0}=9.5$ which yields the solution
$b/b_{0}=0.477$. The corresponding trajectory in the $XY$ plane is shown in
Fig. 2 (dashed line); it agrees quite well with the one reported by Cortés.
Moreover, we can preserve the ratio $a/b\cong 20$ for a finite lattice without
taking the continuum limit $a,b\rightarrow 0$; the result for the same initial
conditions is also shown in Fig. 2 (solid line). We can see a radical
different trayectory characterized by an almost 1-dimensional oscillation
along the $Y$ axis and $x=x(0)=0.23$ constant. The explanation of this
phenomenon is the following: in the finite lattice, the ratio $a/b\cong 20$
means that the probability of the particle’s hopping to neighbor sites along
the $X$ axis is neglible as compared to the sites along the $Y$ axis,
resulting thus in an effective motion due to the $Y$ component force of the
Coulombian attraction; this component is a restoring force with a magnitude
$V_{1}(x(0)^{2}+y^{2})^{-3/2}y$. For the values of the simulation in Fig. 2:
$x(0)=0.23$ and $-1<y<1$, the resulting motion is oscillatory; as the
condition $x(0)\gg y$ is approached, the restoring $Y$ force is approximated
as a linear one, resulting thus in a simple harmonic motion. In the continuum
limit the hopping along the $X$ axis sites is no longer neglible and the
particle follows a 2-dimensional trajectory.

Figure 3: “Precession of the perihelion” trajectory of a particle in the $XY$
plane of a square lattice with adimensional parameters $a=1$, $A=125$,
$V_{1}=20000$ in the 2-dimensional Hamiltonian (1) and conditions at
$t_{0}\cong 0.17\;T$ in the “perihelion” position: $(x_{0},y_{0})=(0,20)$,
$(k_{x0},k_{y0})=(-1,0)$. The impurity is located at $(0,-140)$. In the
continuum limit the precession effect disappears and the quasi-elliptical
orbits collapse into a unique elliptical orbit (dashed line) with period $T$.

Interestingly, if the rectangular lattice is transformed into a square one
($a=b,A=B$) and we take the continuum limit, the resulting trajectory is an
ellipse (Fig. 3, dashed line), as expected in an isotropic continuum space
with a central force field, while in the lattice the trajectory is a succesion
of quasi-elliptical orbits (Fig. 3, solid line) giving rise to the familiar
trajectory showing the precession of the perihelion due to non-central forces
in Newtonian celestial mechanics [6]. The behaviour of the angular momentum
$L_{z}=xk_{y}-yk_{x}$ corresponding to the cases in Fig. 3 is shown in Fig. 4.
The solid line shows the abrupt change of $L_{z}$ in the perihelion, as
expected, while in the rest of the quasi-elliptical trajectory passing through
the aphelion $L_{z}\cong L_{z}^{\textrm{aph}}$ remains almost constant and
equal to their values for other orbits. In the continuum limit, $L_{z}$ tends
to a constant value $L_{z}^{\textrm{cont}}$ (dashed line), again, as expected.
The evident difference between $L_{z}^{\textrm{aph}}$ and
$L_{z}^{\textrm{cont}}$ is due to the fact that in order to maintain constant
the value of the energy $E=2A(1-\cos ak_{x}(0))-V_{1}/y_{0}$, the initial
condition for $k_{x}(0)$ at $(0,y_{0})$ should change according to

$\displaystyle k_{x}(0)$
$\displaystyle=\frac{1}{a}\cos^{-1}\left[1-\frac{a^{2}}{2}(k^{\textrm{cont}}_{x}(0))^{2}\right],$
(10)

where $k_{x}(0)\rightarrow k^{\textrm{cont}}_{x}(0)$ as $a\rightarrow 0$.

Figure 4: Angular momentum $L_{z}(t)=xk_{y}-yk_{x}$ corresponding to the cases
in Fig. 3. $L_{z}$ changes abruptly at the perihelion positions of the
precession orbits (solid line) while in the continuum limit $L_{z}$ tends to a
constant value $L_{z}^{\textrm{cont}}$ (dashed line). The difference of
$L_{z}^{\textrm{cont}}$ and the plateau value of $L_{z}\cong
L_{z}^{\textrm{aph}}$ near the aphelion position is due to the difference of
values for $k_{x}(0)$ at the position $(0,y_{0})$ that guarantee the
conservation of energy.

The results reported in this section concern the application of the
semiclassical model to the Kepler problem on the lattice, taking as the
starting point the Hamiltonian function (1) which was constructed from the
tight-binding energy of the lattice and the Coulomb potential energy due to
the impurity. Since such a tight-binding energy is deduced within the quantum
formalism (see, for example, ch.10 in [7]), the semiclassical results can only
be considered as approximate, being the quantum formalism the exact
theoretical framework to study the Kepler problem on the lattice. It is known
that the semiclassical position r and momentum $\hbar\textbf{k}$ defined in
(2)–(7) correspond to their quantum analogues as the mean values of the
position $\langle\Psi|\textbf{r}|\Psi\rangle$ and momentum
$\langle\Psi|-i\nabla|\Psi\rangle$ operators for a normalized quantum state
$\Psi(\textbf{r},t)$. For such an approximation to hold, the semiclassical
model is restricted by the conditions that the quantum wavepacket width should
be much larger than the lattice constants $a,b,c$ in (1) and much smaller than
the length of an appreciable variation of the external Coulomb field. Besides,
it is required that the energy $V_{1}/r$ of the external field in (1) should
not be too large so that the wavepacket’s energy remains within a single band
(see, for example, ch.12 in [7]). It is also known that the wavepacket spreads
as it propagates in the lattice, being in general such spreading non-
symmetrical with respect to the wavepacket’s center of mass. Therefore, if a
such a non-symmetrical spreading could be related to a kind of “internal
degree of freedom” of the wavepacket, then it would be possible to establish a
significant relation between other dynamical variables than the position and
momentum. Those dynamical variables are the semiclassical angular momentum and
its quantum equivalent. This is done in the next two sections.

##  3 Quantum formalism in the lattice

The formal solution of the time-dependent Schrödinger equation
$i\partial\Psi/\partial t=\widehat{H}\Psi$ (recall that $\hbar=1$) for a
general time-independent one-particle Hamiltonian operator $\widehat{H}$, is
given as [9]

$\displaystyle\Psi(\textbf{r},t)$
$\displaystyle=e^{-i\widehat{H}t}\Psi(\textbf{r},0),$ (11)
$\displaystyle\Psi(\textbf{r},t+\Delta t)$
$\displaystyle=e^{-i\widehat{H}\Delta t}\Psi(\textbf{r},t)$ (12)
$\displaystyle=e^{-i\widehat{V}\Delta t/2}e^{-i\widehat{T}\Delta
t}e^{-i\widehat{V}\Delta t/2}\Psi(\textbf{r},t)+O(\Delta t^{3}),$ (13)

where the Hamiltonian was separated, within the single-band approximation
[12], into a kinetic energy operator $\widehat{T}(\textbf{p})$ (with
$\textbf{p}=-i\nabla$ in the position representation) and a external potential
energy operator $\widehat{V}(\textbf{r})$; the periodicity of the lattice is
incorporated within $\widehat{T}(\textbf{p})$ while the potential energy
corresponds to externally applied forces not intrinsic to the lattice (such as
the electric field). The further splitting of the Hamiltonian as
$\widehat{H}=\widehat{V}/2+\widehat{T}+\widehat{V}/2$ follows from the
application of Baker-Campbell-Hausdorff formula [13] for the non-commuting
operators $\widehat{T}$ and $\widehat{V}$, giving thus an approximate time-
iterative formula whose integration yields the result in (11) within a quite
good approximation for a finite $\Delta t$.

The next step is to transform the algorithm (13) into a practical useful tool
for integrating the Schrödinger equation by expressing (13) in terms of
Fourier ($\mathcal{F}$) and inverse Fourier ($\mathcal{F}^{-1}$) transforms as

$\displaystyle\Psi(\textbf{r},t+\Delta t)=$ $\displaystyle
e^{-iV(\textbf{r})\Delta t/2}\mathcal{F}^{-1}\left[e^{-iT(\textbf{k})\Delta
t}\right.$ $\displaystyle\left.\times\mathcal{F}\left[e^{-iV(\textbf{r})\Delta
t/2}\Psi(\textbf{r},t)\right]\right]+O(\Delta t^{3}).$ (14)

In this expression, the operator $\widehat{T}(\textbf{p})$ is transformed into
the function $T(\textbf{p})$ (with $\textbf{p}=\textbf{k}$) in the reciprocal
space where the inverse Fourier transform ($\mathcal{F}^{-1}$) acts upon [10].
The operator $\widehat{V}(\textbf{r})$ is just the function $V(\textbf{r})$ in
the direct space where the Fourier transform ($\mathcal{F}$) acts upon.
Specifically, for the Kepler problem on the lattice defined by the Hamiltonian
function $H$ in (1), $T(\textbf{k})=2A(1-\cos ak_{x})+2B(1-\cos bk_{y})$ and
$V(\textbf{r})=-V_{1}(x^{2}+y^{2})^{-1/2}$, already applied within the $XY$
plane, as explained below (7).

To complete the description of the quantum formalism in the lattice, we use
the representation of the state $\Psi(\textbf{r},t)$ in the basis of the
complete set of orthonormal Wannier functions $\Phi(\textbf{r}-\textbf{R})$
(the Wannier basis) as
$\Psi(\textbf{r},t)=\sum_{n,m}C_{nm}(t)\Phi(\textbf{r}-\textbf{R})$, where
these functions are localized about the lattice site $\textbf{R}=(na,mb)$
($n,m$ integers) with an extent of the order of the lattice constants $a$ and
$b$ in the rectangular lattice [19]. Therefore, the algorithm (14) is to be
applied to the coefficients $C_{nm}(t)$ (which represent the state
$\Psi(\textbf{r},t)$) to calculate the $p$-th order statistical moments.

The case of $p=0$ corresponds to the normalization condition of the quantum
wavepacket, $\sum_{n}|C_{nm}|^{2}=1$. For $p=1$ we have the mean values of the
position along the $X$ and $Y$ axes: $z_{x}=\sum_{n,m}an|C_{nm}|^{2}$ and
$z_{y}=\sum_{n,m}bm|C_{nm}|^{2}$, respectively. The 2nd moment about the mean
(along the $X$ axis) is the variance: $\sigma_{x}^{2}=\sum_{n,m}(an-
z_{x})^{2}|C_{nm}|^{2}$, whose square root is the standard deviation
$\sigma_{x}$. Interestingly, the mean values of the position $z_{x},z_{y}$
correspond to their classical analogues $x,y$ in the semiclassical model under
the specific physical conditions referred to at the end of Section 2. The
variances $\sigma_{x}^{2},\sigma_{y}^{2}$ represent the dispersion of the
quantum wavepackets relative to their mean positions but these variances do
not have a dynamical variable counterpart. The 3th moments will be considered
in this work as dynamical significant in the sense we will explore below. They
will be defined here as $s_{x}^{3}\equiv\sum_{n,m}(an-z_{x})^{3}|C_{nm}|^{2}$
and the corresponding expression for $s_{y}^{3}$, instead of the so-called
“skewness” which is defined as $s_{x}^{3}/\sigma_{x}^{3}$ and
$s_{y}^{3}/\sigma_{y}^{3}$ [8]. The skewness is a measure of the asymmetry of
the wavepacket distribution about the mean. Notice that while the skewness has
no physical dimensions, $s_{x},s_{y}$ have dimensions of length, and this will
suit our purposes hereafter.

Finally, we calculate the mean value $\langle\textbf{L}\rangle$ of the angular
momentum operator $\textbf{L}=\textbf{r}\times\textbf{p}$
($\textbf{p}=-i\nabla$) in the quantum state approximated as
$\phi(\textbf{r})=D\exp\left({i\int^{\textbf{r}}\textbf{k}(\textbf{r}^{\prime})\cdot
d\textbf{r}^{\prime}}\right)$ according to the zeroth-order term of the WKB
expansion [9]. The result is $\langle\textbf{L}\rangle\equiv
L_{q}=|\textbf{z}\times\textbf{k}|=z_{x}k_{y}-z_{y}k_{x}$, where $z_{x},z_{y}$
are the mean values of the $x,y$ coordinates calculated in the lattice as in
(15) and (16) below.

##  4 Quantum numerical simulations

In this Section we apply the algorithm (14) in terms of direct and inverse
Fourier transforms to the coefficients $C_{nm}(t)$ which represent the quantum
state $\Psi(\textbf{r},t)$. The recursive application of (14) over a time
interval is thus equivalent to solve numerically the Schrödinger equation.
Then, with the values of $C_{nm}(t)$, we calculate specifically the first
order and third order statistical moments as defined in Section 3, which
correspond to the wavepaket’s center of mass and “skewness”, respectively.

For 3-dimensional quantum simulations in a cubic lattice ($a=b=c$) we use

$\displaystyle z_{x}$ $\displaystyle=\sum_{n,m,p}an|C_{nmp}|^{2},$ (15)
$\displaystyle z_{y}$ $\displaystyle=\sum_{n,m,p}am|C_{nmp}|^{2},$ (16)
$\displaystyle z_{z}$ $\displaystyle=\sum_{n,m,p}ap|C_{nmp}|^{2},$ (17)

for the mean-value positions of the quantum wavepacket.

Figure 5: Quasi-Bloch oscillations: values of $z_{y}(t)$ for the initial
Gaussian wavepacket with conditions $(x_{0},y_{0})=(0,32a)$ and
$(k_{x0},k_{y0})=(0,0)$; the donor impurity acting as center of the Coulomb
attraction is located at $(0,-120a)$ and the intensity of the Coulomb
potential in (1) is $V_{1}=307200$ (adimensional units). Figure (a)
corresponds to the 3-dimensional simulation with negligible and equal values
of $z_{x}$ and $z_{z}$; figure (b) corresponds to the 2-dimensional simulation
in the $XY$ plane where $z_{x}$ (not shown) is the same as in case (a).

In Fig. 5(a) we show the numerical simulations of the mean-value positions
corresponding to an initial Gaussian wavepacket with conditions
$(x_{0},y_{0},z_{0})=(0,32,0)$ (in units of the lattice constant $a$) and
$(k_{x0},k_{y0},k_{z0})=(0,0,0)$ (in units of $1/a$); the donnor impurity
acting as center of the Coulomb attraction is located at
$(x_{0},y_{0},z_{0})=(0,-120,0)$ and the intensity of the field is
$V_{1}=307200$ (adimensional units) in (1). The motion of the particle is an
approximate Bloch oscillation along the $Y$ axis with negligible values of
$z_{x}(t)$ and $z_{z}(t)$ as compared to $z_{y}(t)$. The corresponding
2-dimensional simulation for $z_{y}(t)=\sum_{n,m}am|C_{nm}|^{2}$ is shown in
Fig. 5(b). Thus, we confirm that the symmetrical dispersion of the wavepacket
along the $Z$ axis does not affect the quasi-Bloch oscillation effect in the
$XY$ plane (with dominant motion along the $Y$ axis), consistently with the
semiclassical prediction about the 2-dimensional dynamics for the Kepler
problem, being the symmetric dispersion of the wavepacket along the $Z$ axis
an independent phenomenon which does not affect the evolution of $z_{x}$ and
$z_{y}$. Further numerical confirmation that the quasi Bloch oscillation is
well represented in the $XY$ plane is shown in Fig. 6 for the countour plot of
the probability density $|C_{nm}|^{2}$. In fact, the value of $z_{x}(t)$ at
100 % of its time evolution is negative, as can be seen in Fig. 5(a).

Figure 6: Contour plots in the $XY$ plane corresponding to the case of the
quasi-Bloch oscillations shown in Fig. 5.

By changing the initial momentum from $k_{x0}=0$ (quasi-Bloch oscillation) to
$k_{x0}=-1$, the resulting trajectory of the mean-value position
$\textbf{z}=(z_{x},z_{y})$ is shown in Fig. 7 (solid line); the contour plots
show the probability density $|C_{nm}|^{2}$ of the initial Gaussian wavepacket
and the final wavepaket with a skewness vector $\textbf{s}=(s_{x},s_{y})$,
which indicates the asymmetry due to its deformation over time. The dashed
line corresponds to the semiclassical trajectory of the position
$\textbf{r}=(x,y)$ according to the Hamilton equations (2), (3), (5), (6).

Figure 7: Contour plots in the $XY$ plane of the probability density
$|C_{nm}|^{2}$ of the initial (Gaussian) and final wavepackets. The initial
conditions are $(x_{0},y_{0})=(0,32)$ (in units of the lattice constant $a$)
and $(k_{x0},k_{y0})=(-1,0)$ (in units of $1/a$); the donnor impurity acting
as center of the Coulomb attraction is located at $(0,-120)$ and the intensity
of the field is $V_{1}=307200$ (adimensional units) in (1). The solid-line
trajectory of the mean-value position $\textbf{z}=(z_{x},z_{y})$ ends up at
the position $(40.42,10.23)$. The arrow indicates the skewness vector
$\textbf{s}=(-10.50,9.39)$ of the final wavepacket relative to its mean
position. The dashed line corresponds to the semiclassical trajectory of the
position $\textbf{r}=(x,y)$. Figure 8: Numerical evidence of
$s_{x}(t)\rightarrow 0$ in the continuum limit for a particle moving under a
Coulomb potential in a 1-dimensional lattice.

Next, we calculate the time evolution of the quantum angular momentum of the
wavepacket approximated as
$L_{q}=|\textbf{z}\times\textbf{k}|=z_{x}k_{y}-z_{y}k_{x}$, where the mean
value of the position z corresponds to the solid line in Fig. 7, and we
compare it to the semiclassical angular momentum
$L_{c}=|\textbf{r}\times\textbf{k}|=xk_{y}-yk_{x}$ where the position r
corresponds to the dashed line in Fig. 7. In both cases, z and r are measured
relative to the impurity position which is in this case $\textbf{r}=(0,-120)$.
The approximation referred to above for the quantum $L_{q}$ is that we use the
same semiclassical momentum $\textbf{k}=(k_{x},k_{y})$ as for $L_{c}$. Such an
approximation could be justified upon the following basis: it is known that in
the continuum limit and under specific conditions [7], the quantum and
semiclassical positions approximate very well for any time,
$\textbf{z}(t)\cong\textbf{r}(t)$. In the lattice, however, a deviation should
be introduced as $\textbf{z}=\textbf{r}+\alpha\textbf{s}$ such that
$\textbf{s}\rightarrow 0$ in the continuum limit and $\alpha$ being a
dimensionless scale factor. It happens that we already know of such a
quantity: the modified “skewness” defined in this work as
$\textbf{s}=(s_{x},s_{y})$ with $s_{x}^{3}\equiv\sum_{n,m}(an-
z_{x})^{3}|C_{nm}|^{2}$ and $s_{y}^{3}\equiv\sum_{n,m}(an-
z_{y})^{3}|C_{nm}|^{2}$. In Fig. 8 a simulation of $s_{x}$ as $a\rightarrow 0$
was carried out for a 0hopping electron in a 1-dimensional lattice under the
Coulomb field of a donor impurity located at the origin [15]; we can verify
that in fact $s_{x}\rightarrow 0$. We may reasonably suppose that in the
2-dimensional generalization of this problem, the Kepler problem in the
lattice, it also occurs that $s_{y}\rightarrow 0$[16]. Now, we propose in this
work that $L_{q}=L_{c}+\alpha S$ with $L_{q}=|\textbf{z}\times\textbf{k}|$,
$L_{c}=|\textbf{r}\times\textbf{k}|$ and $S=|\textbf{s}\times\textbf{k}|$ (z,
r, s and k are coplanar vectors). The quantity $S$ represents an intrinsic
angular momentum of the wavepacket due to its deformation in the lattice
during its time evolution. In the limit of the continuum, $S\rightarrow 0$ and
the quantum and semiclassical angular momenta coincide, $L_{q}=L_{c}$, as
expected.

Figure 9: Time evolution of the angular momenta corresponding to the case of
Fig. 7: quantum $L_{q}(t)$ (solid line) and semiclassical $L_{c}(t)$ (dashed
line). Figure 10: Time evolution of the momenta difference $L_{q}(t)-L_{c}(t)$
corresponding to Fig. 9 and the intrinsec angular momentum $S(t)$. The
qualitative features of $L_{q}(t)-L_{c}(t)$ are roughly represented by $S(t)$,
specifically: their change of sign at $t\cong 0.1$ and in the interval
$0.4<t<0.6$.

In Fig. 9 we show the numerical values of $L_{q}(t)$ (solid line) and
$L_{c}(t)$ (dashed line). We can see that there is a quite good approximation
between them. In order to quantify such an approximation we calculate their
difference $L_{q}-L_{c}$ and compare it to the value of the intrinsic angular
momentum $S$ (Fig. 10), where an arbitrary scale factor $\alpha$ has been
chosen for the purpose of a better visual comparison. It is to be mentioned
that other odd higher order statistical moments beyond the 3th (skewness) and
usually refered to as “hyperskewness” could also be used to refine our
results. This task is left as an interesting issue for future work.

##  5 Conclusions

We have studied the motion of a charged particle (electron) in a rectangular
lattice in the presence of the Coulomb electric field of another fixed charged
particle (donor impurity). Although the field is 3-dimensional, we
demonstrated that the motion will occur in a plane. The methods used in this
study were the semiclassical model (Hamilton equations) and the quantum
formalism; for the latter we invoked the tight-binding model for the
Hamiltonian and the pseudo-spectral algorithm for the time integration of the
Schrödinger equation. The total energy is conserved but the angular momentum
is not since the lattice discrete subspace has not rotational symmetry. We
propose that these subjects are ideally suited for their pedagogical
adaptation in undergraduate courses as classical mechanics, quantum mechanics
and solid state.

Specifically, we report the following results:

i) We have reproduced a relevant result reported by Cortés [5] about the
trajectory of a particle in an attractive central field. In [5] the author
invokes different masses in the $X$ and $Y$ axes for the same particle. We
have shown that this can be understood in terms of the continuum limit of a
rectangular lattice.

ii) We have obtained the precession effect of a trajectory in a square lattice
as a result of the non-conservation of the angular momentum. In the limit of
the continuum we obtained a closed elliptical orbit (the precession
disappears) and the angular momentum is conserved.

iii) We have applied the quantum formalism in 3 dimensions and integrated
numerically the Schrödinger equation. We have shown that the mean value of the
particle’s position remains in the lattice plane while the dispersion of the
wavepacket in the normal direction does not affect the mean value of the
position in the plane.

iv) We have used the quantum WKB approximation consistent with the
semiclassical method to obtain the mean value of the angular momentum $L_{q}$.
We have compared it with its equivalent semiclassical angular momentum
$L_{c}$. We have proposed that their numerical difference can be attributed to
a “intrinsec angular momentum” $S$ of the wavepacket due to its asymmetrical
deformation in the lattice: $L_{q}=L_{c}+\alpha S$, where $S$ is proportional
to the skewness of a distribution (3th statistical moment). In the continuum
limit ($a\rightarrow 0$), $S\rightarrow 0$ and $L_{q}\rightarrow L_{c}$.

## References

## References

  * [1] J. M. Luttinger, W. Kohn, Phys. Rev. 96 (1954) 802.
https://doi.org/10.1103/PhysRev.96.802.

  * [2] S. D. Boris, S. Brandt, H. D. Dahmen, T. Stroh and M. L. Larsen, Phys. Rev. A 48, 4 (1993) 2574.
https://doi.org/10.1103/PhysRevA.48.2574.

  * [3] Z.-Q. Bai, W.-M. Zheng, Phys. Lett. A 300 (2002) 259-264.
https://doi.org/10.1016/S0375-9601(02)00778-8.

  * [4] O. Petrova, R. Moessner, Phys. Rev. E 93 (2016) 012115.
https://doi.org/10.1103/PhysRevE.93.012115.

  * [5] E. Cortés, Lat. Am. J. Phys. Educ. 2, 2 (2008) 165.
  * [6] K. R. Symon, Mechanics (Addison-Wesley, Cambridge MA, 1953).
  * [7] N. W. Ashcroft and Mermin N. D., Solid State Physics (Philadelphia, PA: Saunders College, 1976).
  * [8] R. E. Walpole and J. E. Freund, Mathematical Statistics (Prentice-Hall, 1980).
  * [9] R. L. Liboff, Introductory Quantum Mechanics (Oakland, CA: Holden-Day, 1980).
  * [10] DeVries P. L., A First Course in Computational Physics (New York: Wiley, 1994).
  * [11] E. L. Wolf, Graphene, a new paradigm in condensed matter and device physics, (1st ed., Oxford University Press, New York, 2014).
  * [12] D. C. Mattis, Rev. Mod. Phys. 58 (1986) 361.
https://doi.org/10.1103/RevModPhys.58.361.

  * [13] Swanson M., Path Integrals and Quantum Processes (Boston, MA: Academic, 1992).
  * [14] J. A. Velasco and D. Sanjinés, Rev. Bol. Fis. 39 (2021) 17.
https://doi.org/10.53287/rkvz2772qg20p.

  * [15] This calculation is a natural extension of the study entitled The dynamics of a particle in a tight-binding lattice under a Coulomb potential [14] wherein we have used both quantum and semiclassical approaches, showing that the results of the former can be quite well approximated by the semiclassical approach.
  * [16] The numerical procedure to verify this assertion on the basis of the algorithm given in (14) is too demanding for a computer since the limit of the continuum $a\rightarrow 0$ implies the discretization of the lattice in a greater number of cells. For the numerical results presented in this work we have taken a region of $127\times 127$ cells where a typical wavepacket evolution lasts around 30 minutes. The further segmentation $a\rightarrow a/2$ would imply a 4-fold time increment, and so on.