Quantum Mechanics on Curves & Surfaces

Printed from: Ian-Eaves/bentsurfaces2.tex on May 15, 2013

Robert Gilmore
Physics Department, Drexel University, Philadelphia, PA 19104
May 15, 2013

Abstract
In order to describe the quantum mechanics of a particle confined to a
curve or a surface in R3 Schrödinger’s equation must be modified in two
ways. The kinetic energy operator must be written in terms of the metric
g ij induced on the surface from the flat space metric Gij = δij . The
potential energy term must be modified by the inclusion of a geometric
potential, necessary to constrain the particle to the curve or surface.

1 Preliminaries
Schrödinger’s Action Integral on a manifold M embedded in R3 has the form
Z  2 
~ ij 2
I= g ∂i ψ∂j ψ + V (x)ψ dS = 0 (1)
2m
This integral over the manifold is minimized (made stationary) subject to the
condition that the solution ψ is nonzero. This condition is usually imposed
through use of a Lagrange multiplier.

1.1 Geometric Parts of the Action Integral

The metric tensor g ij measures distances in the manifold. It is induced from
the Euclidean metric on R3 . The metric tensor on the manifold is determined
by equating distances as measured in the manifold with distances in R3 . For
example, if dui are the coordinate differences between two points in M and dxα
are the distances between these points in a Euclidean coordinate system in R3 ,

δαβ dxα dxβ = ds2 = gij dui duj (2)

and the two (covariant) metrics are related to each other by taking the appro-
priate partial derivatives.
The potential is the sum of two terms: V (x) = Vphysical (x) + Vgeometric (x).
The physical potential describes purely physical (electric, magnetic) interac-
tions. The geometric or confining potential is present in order to confine the
particle to the D-dimensional manifold M in R3 . It is expressed in terms of the
principal curvatures Ri at each point

1
 
D D
!2 
~2  1 X 1 1X 1
Vgeom = − +  (3)
2m 2 i=1 Ri2 2 i=1 Ri

We will be principally interested in cases in which the physical potential is

absent.
We consider two distinct cases. The first case is that of closed curves (knots)
in R3 . The second case is of two-dimensional surfaces in R3 .

2 Curves
2.1 Parameterization
A closed curve is usefully parameterized by trigonometric functions cos(θ + φ).
There is a theorem that every knot is a “Chebyschev” knot: K = C(a, b, c, φ),
where a, b, c are relatively prime integers and C(a, b, c, φ) describes the (x, y, z)
components as (x, y, z) = (Ta (cos(θ)), Tb (cos(θ)), Tc (cos(θ) + φ)) = (cos(aθ),
cos(bθ), cos(cθ + φ)) For now we simply use the parameterization (x(θ), y(θ), z(θ)).

2.2 Metric
We need a metric tensor to describe the kinetic part of the Schrödinger action.
It is determined from
dx dx
ds2 = gij dxi dxj = gθθ dθ2 ⇒ gθθ = · (4)
dθ dθ
Here dx
dθ is proportional to the tangent vector t̂, in fact t̂ =
dx dx
dθ /| dθ |. The kinetic
energy term in the Schrödinger action is
2 −1  2
~2 θθ ~2
 
dψ dx dx dψ
g = · (5)
2m dθ 2m dθ dθ dθ

2.3 Geometric Potential

~2 1


The geometric potential needed to confine a particle to the curve is 2m − 4R 2 .

Here R is the radius of curvature at a point, related to the curvature κ at a

point by R(θ)κ(θ) = 1. The curvature is related to the derivative of the unit
tangent vector t̂ by

dt̂ dt̂ dt̂

= κn̂ ⇒ κ2 = · (6)
dθ dθ dθ
Here n̂ is the unit normal to the curve and t̂ · n̂ = 0. Then

dt̂ d ẋ ẍ ẋ(ẋ · ẍ) ẍ(ẋ · ẋ) − ẋ(ẋ · ẍ)

κn̂ = = √ =√ − 3/2
= 3/2
(7)
dθ dθ ẋ · ẋ ẋ · ẋ (ẋ · ẋ) (ẋ · ẋ)

2
and

1 (ẋ · ẋ)(ẍ · ẍ) − (ẋ · ẍ)(ẋ · ẍ)

κ2 = 2
= √ (8)
R ẋ · ẋ
Here ẋ = dx/dθ, etc.

2.4 Schrödinger Action and Equation

The Schrodinger Action Integral is

~2 1 dψ ∗ dψ 1 (ẋ · ẋ)(ẍ · ẍ) − (ẋ · ẍ)(ẋ · ẍ) ∗

I  
I= − √ ψ ψ dθ (9)
2m ẋ · ẋ dθ dθ 4 ẋ · ẋ
We searchH for solutions that make I stationary, subject to the normalization
condition ψ ∗ ψdθ = 1.
The arguments within the action integral are all functions of the angle θ
that the parameterization of the knot depends on. For Chebyshev knots (a, b, c)
with φ = 0, ±π/2, π the following symmetries hold

θ → −θ
ẋ · ẋ → +ẋ · ẋ
(10)
ẋ · ẍ → −ẋ · ẍ
ẍ · ẍ → +ẍ · ẍ
As a result the argument of the action integral has a two-fold symmetry and the
set of solutions breaks into two subsets: one with even symmetry under θ → −θ
and the other with odd symmetry. The even and odd sets can be expressed as
sines and cosines:

a0 X X
ψeven (θ) = cos(0θ) + aj cos(jθ) ψodd (θ) = bj sin(jθ) (11)
2 j>0 j>0

The variational equation for the amplitudes aj is

X
I→ aj 0 (KE j 0 j + Vj 0 j ) aj (12)
j 0 j≥0

and the matrix eigenvalue equation is

X πE

KE j 0 j + Vj 0 j aj − δj 0 j =0 (13)
~2 /2m
j≥0

The matrix elements are

I
KE j 0 j = dθ j 0 sin(j 0 θ)g θθ (a, b, c)j sin(jθ)
I 
1 2
 (14)
0
Vj 0 j = dθ cos(j θ) − κ (a, b, c) cos(jθ)
4

3
The functions g θθ (a, b, c) and κ2 (a, b, c) are well-defined by Eqs. (4) and (8)
once the knot has been chosen and its Chebyshev expression is introduced to
parameterize the knot.
The variational expression for odd functions is similar, with the obvious
replacements.
Schrödinger’s second order differential equation is obtained in the usual fash-
ion by integrating by parts:
d θθ dψ 1 2 E
− g (a, b, c) − κ (a, b, c)ψ − 2 ψ=0 (15)
dθ dθ 4 ~ /2m

2.5 Propagation Along a Curve

Propagation along the curve is governed by Eq. (15). Since this is an ordinary
differential equation it is useful to write this as a pair of couple ordinary differ-
ential equations. We do this by defining A(θ) = ψ(θ) and B(θ) = g θθ dA(θ)/dθ.
The equations of motion are, in matrix form

1/g θθ
       
d A 0 A 0 gθθ A
= = (16)
dθ B −(E − V ) 0 B −(E − V ) 0 B
This equation must hold for an arbitrary initial condition, so it holds for a 2 × 2
transfer matrix M :
    
d M11 M12 0 gθθ M11 M12
= (17)
dθ M21 M22 −(E − V ) 0 M21 M22
This equation can be integrated from an initial value θi to some final value θf
with M (θi , θi ) = I2 . The matrix M (θf , θi ) relates the amplitudes at the initial
position θi with those at the final position θf .

2.6 Transfer and Scattering Matrix

At any point along the curve the wave function can be written in terms of
the coefficients A = ψ and B = g ∗∗ ψ 0 or in terms of waves travelling to the
right (W+ e+ikθ ) and to the left (W− e−ikθ ). Equating the wavefunction and its
derivative in these two representations leads to the relation
     +ikθ  
ψ=A 1 1 e 0 W+
= (18)
ψ 0 = g∗∗ B ik −ik 0 e−ikθ W−
The two amplitudes W+ , W− are

−1  −1  −1 

e+ikθ
   
W+ 0 1 1 1 0 A
= (19)
W− 0 e−ikθ ik −ik 0 g ∗∗ B
This result can be combined with the expression for the relation between the
amplitudes at θi and θf

4
    
A A
= M (θf , θi ) (20)
B f
B i

to give an expression for the transfer matrix for the amplitudes between these
points:
   
W+ W+
= (GKΘ)−1 f M (θ ,
f i θ )(GKΘ) i (21)
W− f W− i
     +ikθ 
1 0 1 1 e 0
where G = , K = , and Θ = .
0 g ∗∗ ik −ik 0 e−ikθ
The scattering matrix S relates the incoming amplitudes W+,i , W−,f to the
outgoing amplitudes W+,f , W−,i :
   
W+,f W+,i
=S (22)
W−,i W−,f
The matrix elements of the S matrix are obtained from those of the transfer
matrix T in Eq. (21) by simple algebraic manipulations.

3 Surfaces
Although it is often possible to construct analytic representations of surfaces
in R3 , we will often represent them by specifying a set of vertices xi (i =
1, 2, · · · NV ) that are distributed over the surface. This representation is called a
tessellation. We describe how to construct the geometric inputs to the Schrödinger
Action Integral for both the analytic rerpesentation and the tessellation of the
surface.

3.1 Gauss Fundamental Forms

Gauss introduced two fundamental forms in order to describe the properties of
surfaces embedded in R3 :
I Distance ds2 = Edx2 + 2F dxdy + Gdy 2
II Curvature 2dz = Ldx2 + 2M dxdy + N dy 2
Here ds is the distance between the origin and a point with coordinates
(dx, dy) in the surface and 2dz is twice the distance between a point with co-
ordinates (dx, dy) in the surface and with the same coordinates in the plane
tangent to the surface at the origin. In order to construct the Schrödinger ac-
tion integral we will have to construct: the first fundamental form for the surface
at a point in order to construct the kinetic energy part of the integral; the sec-
ond fundamental form for the surface at the same point in order to construct
the geometric potential energy part of the integral.

5
3.2 Construction of the Metric
Choose a vertex x0 and a sufficient number (n ≥ 3) of neighbors of this point.
Set ∆xi = xi − x0 and treat it as a three-component row vector. Nearby points
lie approximately in a plane. The plane can be determined by carrying out a
singular value decomposition on the 3 × 3 covariance matrix constructed from
coordinate differences:
n
X n
X
CM = ∆xti ∆xi = (xi − x0 )t (xi − x0 ) (23)
i=1 i=1
The two large eigenvalues λ1 ≥ λ2  λ3 ≥ 0 are associated with eigenvectors e1
and e2 that span the plane tangent to the surface at x0 . The third eigendirection
e3 is perpendicular to the surface at x0 : it is in the direction of the gradient to
the surface at x0 . The three unit vectors ei are orthonormal.
In the tangent plane Gauss’ First Fundamental Form is gij = ei ·ej so the
kinetic energy term g ij ∂i ψ∂j ψ is the straightforward Euclidean operator. The
kinetic energy term is created without the complication of computing the matrix
elements (E, F, G) = (1, 0, 1). To construct the kinetic energy matrix for the
triangle with three vertices x0 , xi , xj it is sufficient to project the differences
(xi − x0 ) onto the plane: ∆xi = (xi − x0 )·e1 , ∆yi = (xi − x0 )·e2 , and similarly
for xj , and then use the results for the kinetic energy operator on a triangle in
a plane.

3.3 Construction of the Geometric Potential

The geometric potential at a point on a curved surface is a function of the
principal radii of curvature at that point:
2
~2 1 1

1
Vgeometric (x) = − − (24)
2m 4 R1 R2
Here R1 , R2 are the principal radii of curvature of the surface atx0 . These

L M
in turn are the reciprocals of the eigenvalues of the 2 × 2 matrix
M N
of Gauss’ Second Fundamental Form when the First Fundamental Form is the
identity matrix I2 .
These cofficients are determined as follows. Construct the matrices

   
dx1 dx1 2dx1 dy1 dy1 dy1   2dz1  
 dx2 dx2 2dx2 dy2 dy2 dy2  L  2dz2  L  
M =  or T  M  = 2dz
  
 .. .. ..   ..
 . . .  N  .  N
dxn dxn 2dxn dyn dyn dyn 2dzn
(25)
In this expression (dxi , dyi , dzi ) = (e1 , e2 , e3 )·(xi − x0 ). Multiply both sides by
the transpose matrix T t . The matrix T t T is nonsingular. Multiply both sides
by its inverse to obtain

6
  
L
 M  = (T t T )−1 T t 2dz
 
(26)
N
 
L M
Now construct the eigenvalues λ1 = 1/R1 , λ2 = 1/R2 of and use
M N
these values in Eq. (24) for the geometric potential.

3.4 Parameterization Patches

There are many ways to parameterize a surface. One useful method involves
mapping a planar patch with coordinates (u, v) into R3 by means of (u, v) → x ∈
R3 with x(u) = (x(u, v), y(u, v), z(u, v)). Then tangent vectors to the surface
at a point x0 are
∂x ∂x
tu = = xu tv = = xv N = tu ×tv (27)
∂u ∂v
These vectors are not normalized. The unit normal vector is N̂ = xu ×xv /|xu ×xv |.
In this parameterization the First and Second Fundamental forms are

       
E F xu ·xu xu ·xv xuv ·N̂
L M xuu ·N̂
I= = II = =
F G xv ·xu xv ·xv M
xvv ·N̂N xvu ·N̂
(28)
The principal curvatures κ1 = 1/R1 , κ2 = 1/R2 are obtained by solving the
eigenvalue equation:
 
du
([II] − κ [I]) =0 (29)
dv
The corresponding eigenvectors are the principal directions.
If only the principla curvatures are required, then it is useful to observe that
the eigenvalues can be obtained by rewriting Eq. (28) as
−1 −1
K − κI2 = [II] [I] − κI2 = [I] [II] − κI2 (30)
The geometric potential is then
2 !
~2

1 1 2
Vgeom = tr(K) − tr(K ) (31)
2m 2 2
 
−1 a b
This eliminates the need for diagonalizing a matrix. If [II] [I] =
c d
2 2 2
then (κ1 − κ2 ) = (1/R1 − 1/r2 ) = (a + d) − 4(ad − bc).
The kinetic energy can also be determined in this parameterization. Assume
we have three vertices in the parameterizing patch: u1 , u2 , u3 with coordinates
(u, v)i , i = 1, 2, 3. Then the basis function at vertex 1 is

7
    
 u v 1   u1 v1 1 
   
φ1 (u, v) =  u2 v2 1  /  u2
  v2 1 
 (32)
 u3 v3 1   u3 v3 1 

with φ2 , φ3 defined similarly. The gradient of ψ(u, v) = c1 φ1 + c2 φ2 + c3 φ3 is

constructed from

   
v2 − v3 u2 − u3
∂ψ ∗

 v3 − v1  /D ∂ψ ∗  u3 − u1  /D

∗
= c1 c2 c3 c1 c2 c3 =−
∂u ∂v
v1 − v2 u1 − u2
(33)
where
R R αβ D is the
√ denominator in Eq. (32). The integral over this simplex is
g ψα ψβ detg∗∗ du ∧ dv. The derivatives ψα are independent of (u, v) and
the metric tensor can be approximated as the value of the metric at any of
the three vertices, or even the average of the values at the three vertices. The
intermediate result is

 
v2 − v3 −(u2 − u3 )  11
g 12

A

 v3 − v1 g
K= 2 c1 c2 c3 −(u3 − u1 ) ×
g 21 g 22

D
v1 − v2 −(u1 − u2 )

 
  c1
v2 − v3 v3 − v1 v1 − v2 p
 c2  × detg∗∗ (34)
−(u2 − u3 ) −(u3 − u1 ) −(u1 − u2 )
c3

Here A is 12 D. We have explicit expressions for the other terms in this expression:
 
p p
2 ∗∗ 1 G −F
detg∗∗ = EG − F g = (35)
EG − F 2 −F E
Putting this all together, we find
 
  c1
1 1 E F
M t  c2 


K= √ c1 c2 c3 M (36)
2D EG − F 2 F G
c3
with
 
u2 − u3 v2 − v3
M =  u3 − u1 v3 − v1  (37)
u1 − u2 v1 − v2
For the flat metric g ∗∗ = δ ij the result reduces to the standard result shown in
Eq. (40).

8
3.5 Constraint Functions
Another useful way to describe a surface is through a function f (x, y, z) =
const. A single constraint on three variables typically produces a smooth two-
dimensional surface. In the event a function is not prescribed, the surface can
be locally parameterized in this way by fitting the vertices in the neighborhood
of a vertex of interest (x0 ) to a smooth function using a best fit or singular
value procedure (Appendix). For our purposes a fit containing only the first
and second degree terms in the Taylor series expansion is sufficient.
The normal at a point x0 is N = ∇f |x0 . Tangent vectors to the surface are
conveniently obtained as follows. Determine the magnitude of the components
of N. If |fy = f2 | is smaller than the other two partial derivatives, then take
the first tangent vector as t1 ' (f3 , 0, −f1 ) and the other tangent vector as
t2 ' N × t1 .
In the tangent plane when t1 and t2 are normalized to one, the metric tensor
is flat: gij = g ij = δij and construction of the kinetic energy matrix elements
proceeds as in flat space.
The curvature form is obtained from the first and second derivatives as
follows. Construct the 3 × 3 matrix K as follows:
 
1 (fir fr )(fs fsj )
K= fij − (38)
|N| (ft ftu fu )
The matrix K has one vanishing eigenvalue, corresponding to the normal direc-
tion. The other two eigenvalues are the curvatures in the principal directions.
As in subsection (3.4), the geometric potential is
2 !
~2

1 1 2
Vgeom = tr (K) − tr (K ) (39)
2m 2 2
This again eliminates the need for diagonalizing a matrix.

3.6 Constructing the FEM Matrices

The kinetic energy matrix can be constructed using the vector positions x1 , x2 , x3
of the three vertices of a triangle. The matrix elements are expressed in terms
of inner products as follows
 
(x2 − x3 )·(x2 − x3 ) −(x1 − x3 )·(x2 − x3 ) −(x1 − x2 )·(x3 − x2 )
1 
K4(x1 x2 x3 ) = −(x2 − x3 )·(x1 − x3 ) (x3 − x1 )·(x3 − x1 ) −(x2 − x1 )·(x3 − x1 ) 
4A
−(x3 − x2 )·(x1 − x2 ) −(x3 − x1 )·(x2 − x1 ) (x1 − x2 )·(x1 − x2 )
(40)
where A is the area of the triangle, 2A = |(x2 − x1 )×(x3 − x1 )|.
The potential energy matrix over the same triangle is
 
(6, 2, 2) (2, 2, 1) (2, 1, 2)
A 
V4(x1 x2 x3 ) = (2, 2, 1) (2, 6, 2) (1, 2, 2)  (41)
60
(2, 1, 2) (1, 2, 2) (2, 2, 6)

9
where the meaning of (2, 1, 2) in the (1, 3) matrix element is (2, 1, 2) → 2 ×
V (x1 ) + 1 × V (x2 ) + 2 × V (x3 ).
The overlap matrix is obtained from the potential energy matrix by inserting
the constant potential V = 1 everywhere:
 
2 1 1
A 
O4(x1 x2 x3 ) = 1 2 1  (42)
12
1 1 2

3.7 Periodic Boundary Conditions

Application of Dirichlet and von Neumann boundary conditions within the Fi-
nite Element Method is well-understood. Here is how to impose periodic bound-
ary conditions. Assume the matrix elements have been computed for two ver-
tices, and the vertices are labelled i = (1, 2, 3) for one and r = (4, 5, 6) for the
second. Suppose also that after all the quadratic forms have been constructed,
we impose the condition that vertex 1 is the same as vertex 6. Then c1 = c6 = cj
(j for “joint”). In each quadratic form we replace c1 by cj wherever it occurs,
and we also replace c6 by cj wherever it occurs. The quadratic form marix will
have one fewer row and column. It is constructed by adding rows 1 and 6 and
also adding columns 1 and 6 together. If whole edges are to be identified, as in
cylinder or torus boundary conditions, the rows for the pairs (triples, quartets)
of vertices in the tesselaltion that are identified are to be combined by addition.
Example: A square [0, 2π] × [0, 2π] is to be mapped onto a torus. Divide
the θ into 100 intervals with 101 points. Similarly for the φ direction. There
will be a total of 1012 vertices. The 99 interior vertices along both the θ and
φ directions are pairwise identified and added, eliminating 2 × 99 rows and
columns. The four vertices at [0, 0] , [0, 2π] , [2π, 0] , [2π, 2π] are identifed with
a single point on the torus, squeezing four rows/columns down to one. The
number of rows/columns in the matrices representing the quadratic forms for
the Schrödinger equation on the torus is 1012 − 2 × 99 − (4 − 1) = 1002 , as one
would hope. The generalized eigenvalue equation then provides approximations
to the solution on the boundary-free torus.

3.8 Scaling
Suppose we map a tessellated unit square into a torus. Assume the tessellation
has 101 points per edge, and the torus has radii of the major and minor circles
a > b. How do the energies/wavefunctions scale if we mapped the same square
into a larger (smaller) torus whose length parameters are λa > λb?
Since length in the embedding scales like λ2 and ∆x → λ∆x, the metric
tensor g∗∗ doesn’t scale, nor does its inverse g ∗∗ . That means the kinetic energy
∂ψ ∂ψ
KE ' g ∗∗ λ∂x 2
∗ λ∂x∗ scales like 1/λ . Similarl;y, the curvatures scale like λ so

the geometric potential scales like 1/λ2 . As a result, the energies scale like 1/λ2 .
Since the surfaace area scales like λ2 and the output column vectors from
the diagonaliztion are unchanges, these column vectors are mapped onto spacial

10
wavefunctions by multiplying by the vactor 1/λ2 .

4 Examples
4.1 Ellipse
An ellipse in a plane can be parameterized in terms of an angle θ:

(x, y) = (a cos θ, b sin θ) ẋ = (−a sin θ, b cos θ) ẍ = (−a cos θ, b sin θ) (43)

The metric and curvature are

ds2 = dx2 + dy 2 = (a2 sin2 θ + b2 cos2 θ)(dθ)2 ⇒ gθθ = (a2 sin2 θ + b2 cos2 θ)
(a2 sin2 θ + b2 cos2 θ)(a2 cos2 θ + b2 sin2 θ) − (a2 − b2 ) sin θ cos θ
κ2 =
(a2 sin2 θ + b2 cos2 θ)3
(44)
In particular,
a b
κ(±a,0) = κ(0,±b) = (45)
b2 a2

4.2 Ellipsoid
A triaxial ellipsoid is determined by the function

x2 y2 z2
2
+ 2 + 2 = cst. (46)
a b c
At any point (x, y, z) on this surface the normal is

N ' (2x/a2 , 2y/b2 , 2z/c2 ) (47)

Unit vectors in the tangent plane at this point can be taken as

xy x2 z2
z  
x yz
e1 = n1 2
, 0, − 2 e2 = n2 − 2 2 , 4 + 4 , − 2 2 (48)
c a a b a c b c
where n1 , n2 are appropriate normalization coefficients. The curvatures at
(x, y, z) are obtained from Eq. (38):
  x  
 a 4 

  y  h
x y z
i 

  1    b4  a4

b 4 c 4

 a2 0 0
 
z
 

1 c4
 0
b 2 0  − y
(( a3 )2 +( b3 )2 +( c3 )2 ) 
x z
0 0 c12

 

 


 

 
K= p (49)
(x/a2 )2 + (y/b2 )2 + (z/c2 )2

11
The curvatures at the endpoints of the principal axes are
a a b b c c
κ±a,0,0 = , κ0,±b,0 = , κ0,0,±c = , (50)
b2 c2 c2 a2 a2 b2

4.3 Tori
A recent theorem guarantees that every knot can be expressed as a “Chebyshev”
knot: K(θ) = (cos aθ, cos bθ, cos cθ + ψ) where a, b, c are relatively prime positive
integers and ψ is some angle (0 or π/2 are the usual culprits). The first and
second derivatives are

K̇(θ) = (−a sin aθ, −b sin bθ, −c sin cθ + ψ)

(51)
K̈(θ) = (−a2 cos aθ, −b2 cos bθ, −c2 cos cθ + ψ)
These are related to the tangent and normal vectors by t̂ = K̇(θ)/|K̇(θ)| and
 
n̂ = α K̈(θ) − (K̈(θ)·t̂)·t̂ (52)

with α the normalization constant. Define the binormal b̂ = t̂ × n̂. Then the
torus T (θ, φ) surrounding the knot K(θ) is defined by

T (θ, φ) = K(θ) + a(n̂ cos φ + b̂ sin φ) (53)

for a sufficiently small.

4.4 Surfaces Without Boundaries: Torus

The parameterization of a torus as described above involves mapping a square
0 ≤ θ, φ ≤ 2π onto the torus surface using the mapping described in Eq. (53).
A tessellation of the square is easily made in the plane. This tessellation can be
lifted to the torus with ease. However, there is a technical detail. Vertices at
opposite sides of the square map to the same vertices on the torus.
To be explicit, suppose we represent vertices in the plane by (iθ , iφ ), with
iθ and iφ integers in the range 0 ≤ iθ , iφ ≤ 100, so that the tessellation takes
the form of intersecting horizontal and vertical lines. Each line will contain 101
vertices and the kinetic, potential, and overlap matrices will be 1012 × 1012 .
However, the two vertices (0, iφ ) and (100, iφ ) will map to the same vertex in
the torus (iφ 6= 0, 100), and the four vertices (0, 0), (0, 100), (100, 0), (100, 100)
(the four corners of the square in the θ-φ plane) will all map to a single vertex
in the torus. When this occurs the periodic boundary conditions are imposed
as follows. If a and b are two vertices in the flat parameterizing plane that map
to the same vertex in the torus, the row labeled b is added to the row labeled
a and row b is eliminated. Ditto for columns. In this way the 99 rows (iθ , 0)
and (iθ , 100), 0 < iθ < 100 are consolidated. Similarly for (0, iφ ) and (100, iφ ).
Finally, the contents of the three rows/columns (0, 100), (100, 0), (100, 100) of
the kinetic, potential and overlap matrices are added to the contents of the
row/column (0, 0) and removed. The resulting matrix has size (1012 − 2 × 99 −

12
3) × 1002 . Solution of the generalized eigenvalue problem (K + V − λO)c = 0
yields the values of the eigenstates ψ at the 1002 vertices distributed over the
torus. The periodic boundary conditions are enforced by this simple procedure.

4.4.1 Implementation
Creata a list of vertices that are to be identified, for example [a, b, c] with a <
b < c. There can be a complementary list of “independent” vertices: those that
are not matched to others under periodic boundary conditions. Now create an
N1 × N2 matrix S of zeroes. N1 is the number of vertices after the identification
is made (ex: 1002 ). Now start with the first list, and if [a, b, c] is the first set
of identifications in this list, overwrite a +1 in positions (1, a), (1, b), (1, c). For
the second member of this list go to the second column of S and overwrite the
zeroes with 1 (ex: (2, d), (2, e). Cointinue until the list of identified vertices is
exhausted. Now go down the singlet list and place a +1 in the next column and
appropriate row.
Then the mapping

H → S t HS = H 0 (54)
converts the initial tessellation/computation to one satisfying the appropriate
boundary conditions. When the smaller hamiltonian H 0 is diagonalized, the
column vectors [ψ(α)] of length N2 on the surface with periodic boundary con-
ditions is converted (if desirable) to a column vector of length N1 by multiplying
by S. Doing this places the same amplitude at each of the vertices that are iden-
tified. Plotting over the initial tessellation than looks more continuous.

4.5 Genus-g
Any torus of genus g > 1 can be decomposed into unions of “trinions” (particle
physics speak). The canonical decomposition contains g − 1 pairs of trinions,
and each pair consists of a “joining” trinion and a “splitting” trinion (nonlin-
ear dynamics-speak). In turn, each trinion can be decomposed into a pair of
hexagons. Thus, a torus of genus g can be tessellated by 4(g − 1) hexagons (see
Fig. ?). Tessellating a hexagon is a piece of cake. Fig. ?? shows a hexagon
tessellated with a hexagonal close packed lattice (solid state physics speak). If
there are n equally spaced vertices per
Pedge, the total number of vertices in the
2n−1
tessellation of each hexagon is T = 2 n k − (2n − 1) = 3n2 − 3n + 1.
After 4 hexagons have been assembled into a genus-2 torus, the number
of independent vertices, edges, and faces (v, e, f of dimensions 0, 1, 2) satisfy
Euler’s beautiful equality:

v − e + f = χ(g = 2) = 2 − 2g = −2 (55)
Every edge shares two faces and every face has three edges, so

13
> a:=3:b:=4:c:=5:x0[1]:=cos(a*t):x0[2]:=cos(b*t):x0[3]:=cos(c*t+Pi/2
):
> for i from 1 to 3 do x1[i]:=diff(x0[i],t):x2[i]:=diff(x1[i],t):od:
> ip00:=x0[1]*x0[1]+x0[2]*x0[2]+x0[3]*x0[3]:
> ip11:=x1[1]*x1[1]+x1[2]*x1[2]+x1[3]*x1[3]:
> ip12:=x1[1]*x2[1]+x1[2]*x2[2]+x1[3]*x2[3]:
> ip22:=x2[1]*x2[1]+x2[2]*x2[2]+x2[3]*x2[3]:
> vgeom:=(ip11*ip22-ip12*ip12)/sqrt(ip11):
> plot(vgeom,t=0..2*Pi);

2500

2000

1500

1000

500
0 1 2 3 4 5 6
t

> plot(ip11,t=0..2*Pi);
 40

35

30

25

20

15

10

0 1 2 3 4 5 6
t

> plot(ip12,t=0..2*Pi);
 80

60

40

20

0 1 2 3 4 5 6
t
–20

–40

–60

–80

> plot(ip22,t=0..2*Pi);
 800

600

400

200

0 1 2 3 4 5 6
t

> plot(ip12^2,t=0..2*Pi);
 7000

6000

5000

4000

3000

2000

1000

0 1 2 3 4 5 6
t

> plot(1/ip11,t=0..2*Pi);ke:=matrix(20,20);
 0.16

0.14

0.12

0.1

0.08

0.06

0.04

0 1 2 3 4 5 6
t

ke := array(1 .. 20, 1 .. 20, [ ])

>
>
>
> for i from 1 to 5 do for j from 1 to 5 do
ke[i,j]:=evalf(int(i*sin(i*t)*j*sin(j*t)/ip11,t=0..2*Pi));print(i,
j,ke[i,j]):od:od:
1, 1, 0.1804363349
1, 2, -0.1110423055 10 -12
1, 3, 0.01301598313
1, 4, 0.1491545251 10 -12
1, 5, -0.1941269700
2, 1, -0.1110423055 10 -12
 2, 2, 0.7390999839
2, 3, -0.1098290210 10 -12
2, 4, -0.2758938636
2, 5, 0.
3, 1, 0.01301598313
3, 2, -0.1098290210 10 -12
3, 3, 1.313546418
3, 4, 0.3990922049 10 -12
3, 5, -0.7093532830
4, 1, 0.1491545251 10 -12
4, 2, -0.2758938636
4, 3, 0.3990922049 10 -12
4, 4, 2.130337857
4, 5, 0.
5, 1, -0.1941269700
5, 2, 0.
5, 3, -0.7093532830
5, 4, 0.
5, 5, 5.531479229
> ###
> ###
> ###
> ###
> pe:=matrix(20,20):
> for i from 1 to 5 do for j from 1 to 5 do
pe[i,j]:=evalf(-1*int(cos(i*t)*cos(j*t)*vgeom/4.0,t=0..2*Pi)):prin
t(i,j,pe[i,j]):od:od:
1, 1, -1415.425823
1, 2, -0.5852048669 10 -11
1, 3, -361.8977670
1, 4, -0.5931163371 10 -11
1, 5, -109.3177836
2, 1, -0.5852048669 10 -11
2, 2, -1342.626940
2, 3, -0.5887470617 10 -11
2, 4, -182.1166668
2, 5, -0.1743810662 10 -9
3, 1, -361.8977670
3, 2, -0.5887470617 10 -11
 3, 3, -1162.845839
3, 4, -0.6013436116 10 -11
3, 5, -261.7869334
4, 1, -0.5931163371 10 -11
4, 2, -182.1166668
4, 3, -0.6013436116 10 -11
4, 4, -1242.516106
4, 5, -0.5952799506 10 -11
5, 1, -109.3177836
5, 2, -0.1743810662 10 -9
5, 3, -261.7869334
5, 4, -0.5952799506 10 -11
5, 5, -901.6795583
>