Eindhoven University of Technology

BACHELOR

Energy levels in one-dimensional hydrogen atoms and Rydberg crystals

van der Weerden, T.H.P.

Award date:
2014

Link to publication

Disclaimer
This document contains a student thesis (bachelor's or master's), as authored by a student at Eindhoven University of Technology. Student
theses are made available in the TU/e repository upon obtaining the required degree. The grade received is not published on the document
as presented in the repository. The required complexity or quality of research of student theses may vary by program, and the required
minimum study period may vary in duration.

General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners
and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.

• Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
• You may not further distribute the material or use it for any profit-making activity or commercial gain
 Energy levels in one-dimensional
hydrogen atoms and Rydberg crystals

T.H.P. van der Weerden

August 2014

CQT 2014-06

Bachelor Thesis

Supervisor Physics: dr. ir. Servaas

Kokkelmans
Supervisor Mathematics: dr. Georg Prokert

Applied Physics: Coherence and Quantum

Technology
Applied Mathematics: Centre for Analysis,
Scientific computing and Applications

Eindhoven University of Technology

 Abstract
We theoretically study the energy levels and wave functions of the electrons in Rydberg
crystals, a collection of Rydberg atoms on a lattice. A Rydberg atom can be represented
as an electron with a very high principal quantum number around a shielded nucleus. We
can simplify the system to a toy model by using the Born-Oppenheimer approximation
and assuming a shielded nucleus generates a Coulomb potential. This model is the
simplest in one dimension with a single Rydberg atom, in which case the model represents
a one-dimensional hydrogen atom. This happens to be a controversial subject and after
considering different views, it is found that the singularity in the Coulomb potential acts,
in one dimension, as an impenetrable barrier, which means there are two independent
regions for the electron. Other than that, the problem is very similar to three-dimensional
Hydrogen. For two Rydberg atoms the singularities cause the model to cease representing
the physical system. Regularizing, and thus removing the singularity, fixes this problem.
Electron-electron repulsion has straightforward effects on the system.
Contents

1 Introduction 1

2 Physical model 3
2.1 The Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 One-dimensional hydrogen . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2.1 General solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2.2 Different views . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.3 Anomalous states . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Numerical methods 13
3.1 Single well . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.1.1 Energy spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.2 Wave functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Double well . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2.1 Merge method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.2 Partition method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4 Regularization 27
4.1 Single electron . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.2 Two electrons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2.1 Perturbational . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.2.2 Numerical . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

5 Conclusion 35

Bibliography 37
1. Introduction

Classical computers, like the one this thesis is written on, are based on binary digits
(bits). A single bit can at any given time be in only one of two possible states, most
commonly represented as 0 and 1. Quantum bits (qubits), on the other hand, can, by
the laws of quantum mechanics, be in a superposition of states, which roughly means
that a qubit can simultaneously be 0 and 1. Subsequently, three qubits each represent
two possible values at once (both 0 and 1) so that they hold as much information as
eight (two times two times two) classical bits. In general N qubits correspond to 2N
bits, an exponential increase in information storage. For example, all the data on a hard
disk of 1 TB (nearly 1013 bits) can be stored in just 43 qubits and furthermore, 300
qubits represent more bits than the amount of atoms in the universe. Because of this
promising outlook for computing tasks, there has been a lot of research in this field of
quantum computing.
It was Richard Feynman, who proposed a basic model of the quantum computer in
1982 [1], later described by David Deutsch in 1985, when he came up with the quantum
Turing machine [2]. For the same reason that quantum computers, using the quantum
phenomena superposition and entanglement, can handle exponentially more states, they
cannot be simulated by a classical computer. An example of a computation that can
be done faster on a quantum computer is Shor’s algorithm, an integer factorization al-
gorithm that can run in polynomial time on quantum computers [3], in contrast to the
general number field sieve, which is classically the fastest integer factorization algorithm
and runs in exponential time [4]. These long factorization times are the basis for cryptog-
raphy, the field concerning, among other things, the protection of money transfers and
secure authentication of passwords. The present cryptography will be rendered useless
with the introduction of a working quantum computer, so that more secure techniques
are needed, such as those studied in quantum cryptography.
Instead of classical supercomputers, quantum simulators can be used to simulate
the behavior of quantum computers [5]. Quantum simulators are real many particle
systems that exhibit quantum behavior and that can be controlled by experimentalists
(quantum computers can be seen as a more practical quantum simulator). Examples of
possible quantum simulators are cold atoms in optical lattices, trapped ions, photons
and quantum dots, each possessing a controllable quantum mechanical property in a
collection of particles. This thesis is concerned with Rydberg crystals, which falls under

1
CHAPTER 1. INTRODUCTION

the category of cold atoms in optical lattices.

A Rydberg crystal consists of ultracold Rydberg atoms, placed in a lattice [6]. Ul-
tracold means that the atoms are in the temperature regime of micro- and nano-Kelvins
so that they are approximately standing still and are therefore easily controllable, which
makes ultracold atoms a good candidate for a quantum simulator. Rydberg atoms are
atoms that have a single electron with a very high principal quantum number [7], which
gives rise to properties such as a high atom radius and a very high dipole moment for
very strong long-range interactions. Another interesting aspect of Rydberg atoms is
that they exhibit a blockade effect, which means that a Rydberg atoms shifts the energy
levels of nearby atoms so that they cannot be excited to a Rydberg state with the same
laser wavelength [8].
When several (more than two) Rydberg atoms are brought together to form a crystal
due to the blockade effect for example, it is possible to excite an electron in such a way
that its wave function is still confined inside the lattice, but the electron is free to move
around all the nuclei. This state is called a quantum plasma and the corresponding
energies of the electrons will be derived from Schrödinger’s equation. From the energy
levels the effects of changing the principal number and the interatomic distance on the
interaction behaviour between the Rydberg atoms in a crystal can be deduced, which is
the goal of this project.
The idea is to make a one-dimensional model of a Rydberg crystal and see how
far Schrödinger’s equation can describe the behaviour of the electrons and their energy
levels. The simplest case of this model is a crystal with a single atom, which is nothing
else than the one-dimensional hydrogen atom, a subject on which many different views
have been given the last decades [10–13]. However, there is still no definite conclusion
on which everyone agrees, so it is crucial to explore this problem before advancing to
the Rydberg crystals.

2
2. Physical model

In this section we will establish a simplified physical model based on known properties of
Rydberg atoms. The idea is to make crude assumptions for the terms in the Hamiltonian
and formulate a toy model based on Schrödinger’s equation. The single atom case
happens to be identical to the one-dimensional hydrogen atom, which will be investigated
analytically. This is mainly done by discussing different views from various articles
[9–13].

2.1 The Hamiltonian

A Rydberg atom can be represented as a single electron, charge −q and mass m, in a

large orbit around a shielded nucleus, which, very roughly, is an ion with charge +q for
states where the electron is never close to the nucleus. An example of a precise model
for the spherically symmetric potential at a radial distance r exhibited by the shielded
nucleus is given by Marinescu et al. [14],
q 2 −Zl (r) αc h 6
i
VM ar (r) = − 4 1 − e−(r/rc ) , r > 0, (2.1)
4πε0 r 2r
where ε0 is the permittivity of vacuum and Zl (r) is a model nuclear charge, depending
on the nuclear charge number Z, given by
Zl (r) = 1 + (Z − 1)e−a1 r − r(a3 + a4 r)e−a2 r , r > 0, (2.2)
with parameters a1 , . . . , a4 and αc that are determined by experiment and are dependent
on the orbital quantum number l. As we just want to establish a toy model, we can
instead use a simpler approximation for (2.1), namely the Coulomb potential,
q2 1
VC (r) = , r > 0, (2.3)
4πε0 r
which is especially true for states in which the electron does not penetrate the nucleus.
This can be generalized for N Rydberg atoms on a one-dimensional lattice to
N
q2 X 1
V (x) = − , x ∈ R\{xi | i = 1, . . . N }, (2.4)
4πε0 |x − xi |
i=1

3
CHAPTER 2. PHYSICAL MODEL

Figure 2.1: Example energy levels of a single electron, charge −q, in the Coulomb
potential of five ions, charge +q, at xi , an unspecified distance 2a apart. Source: [15]

where the xi are the positions of the nuclei.

In the sub-Kelvin regime and on the time scale experiments are performed, the posi-
tions of the atoms can be regarded as fixed. This helps in making a simple model in two
ways. First, we can apply the Born-Oppenheimer approximation and neglect the kinetic
energies of the ions. Second, it makes the potentials of the ions time-independent (the
xi in (2.4) are time-independent). A simple case of the problem is when all the atoms
are placed on an equidistant lattice and then N − 1 atoms are ionized and a single atom
is brought into a Rydberg state. The Hamiltonian Ĥ for the single electron case then
becomes
N
~2 d q2 X 1
Ĥ = − − , x ∈ R\{xi | i = 1, . . . N }, (2.5)
2m dx2 4πε0 |x − xi |
i=1

with ~ the reduced Planck constant and xi − xi−1 =: 2a the constant interionic distance
for i = 2, . . . , N . This Hamiltonian has a discrete set of eigenenergies dependent on the
principal quantum number n, shown in Figure 2.1 for N = 5 and an arbitrary a.
We can distinguish two kinds eigenstates. For low energies the corresponding wave
function of the electron is bound to one (or by symmetry two) ion(s), while for higher
energies the wave function can span over the entire range of the lattice. The latter case
is what is called a quantum plasma, a state in which the electron is free to move around
all ions, but not far outside the lattice.

4
 2.2. ONE-DIMENSIONAL HYDROGEN

2.2 One-dimensional hydrogen

We start by analytically solving the simplest one-dimensional case, a single electron

in the proximity of a single ion, charge +q, at x = 0. According to the toy model
Schrödinger’s equation for the wave function ψ(x) is given by

~ d2 ψ(x) q2 1
− − ψ(x) = Eψ(x), x ∈ R\{0}. (2.6)
2m dx2 4πε0 |x|

This is the exact same equation as for the one-dimensional hydrogen atom, a contro-
versial subject without any decisive arguments [9–13, 18]. It is therefore necessary to
be precautious with assumptions. Surely, we are interested in bound (E < 0) solutions
ψ that are square-integrable functions, so ψ is an element of the Hilbert space L2 (R).
Scattering states (E > 0) are less interesting from a Rydberg crystal point of view. Ac-
tually we are looking for solutions that are in a smaller space, a weighted L2 (R) if you
will, because we would also like to demand that the expectation value of the potential
energy,
Z b
hV i = ψ ∗ (x)V (x)ψ(x) dx, a, b ∈ R, b > a, (2.7)
a

exists and is finite. By (2.6) and the fact that ψ ∈ L2 (R) it follows that the expectation
value of the kinetic energy,
b
~2 d2
Z
hT i = − ψ ∗ (x) ψ(x) dx, a, b ∈ R, b > a, (2.8)
2m a dx2

exists and is finite as well. If hV i is infinite, it would mean that hT i is also infinite, but
in such a way that their sum, hT i + hV i = hHi = E, is finite. Because the integrand
q2
in hV i, 4πε 0
|ψ(x)|2 /|x|, is strictly negative, it suffices to find a single interval [a, b] on
which hV i is infinite to show that a solution does not meet our requirements. This finite
potential energy restriction is also used by Andrews [9]. For the solution of (2.6) we
start by closely following Loudon [10], a well-respected approach.

2.2.1 General solution

By introducing a coordinate transformation x → y and a positive constant λ,

√ r
−2mE q2 −m
y=2 x and λ = > 0, (2.9)
~ 4πε0 ~ 2E
we obtain a new function φ,
√
2 −2mE
φ(y) = ψ( x), (2.10)
~

5
CHAPTER 2. PHYSICAL MODEL

that has to satisfy similar restrictions as ψ(x) and the equation

d2 φ(y)
 
λ 1
+ − φ(y) = 0, y ∈ R\{0}. (2.11)
dy 2 |y| 4

First we examine this equation on only R+ = {y ∈ R| y > 0}, which is a special case of
Whittaker’s equation on R+ with µ = 1/2 [16],

d2 φ(y) λ 1 1/4 − µ2
 
+ − + φ(y) = 0, y ∈ R+ . (2.12)
dy 2 y 4 y2

For µ = 1/2 and y > 0 the solutions of this equation can be expressed in terms of
Kummer’s confluent hypergeometric functions U and M ,

φU,λ (y) = y e−y/2 U (1 − λ, 2; y) and φM,λ (y) = y e−y/2 M (1 − λ, 2; y). (2.13)

From the behavior of these confluent hypergeometric functions we can find the solu-
tions that satisfy our restrictions.
• For small y we have
(
1/y λ∈
/N
U (1 − λ, 2, y) ∼ and M (1 − λ, 2, y) ∼ 1, (2.14)
1 λ∈N

so that at y = 0 φU,λ (y) is zero for positive integer λ and finite for other λ, while
φM,λ (y) is always zero at y = 0.
• Furthermore they behave for large y as follows,

ey y −1−λ (−y)λ−1
 
λ−1
U (1 − λ, 2, y) ∼ y and M (1 − λ, 2, y) ∼ + , (2.15)
Γ(1 − λ) Γ(1 + λ)

where Γ(s) is the gamma function. This means that φU,λ is always square-integrable,
but φM,λ only if λ is a positive integer (a distinction Loudon did not make and
stated that φM,λ (y) diverged for all λ). The gamma function brings for those λ the
exponential down so that φM,λ (y) behaves as a polynomial multiplied by e−y/2 ,
which is square-integrable.
• Finally, we can investigate the expectation of the potential energy hV i. As stated
before, it suffices to find a single interval [a, b] to disprove a solution. As we
know that we are dealing with smooth and square-integrable solutions, only the
behaviour at small y is interesting. This leads to considering the interval [ε, b] for
some b > ε and taking the limit ε → 0. As we know what U (1 − λ, 2, y) and
M (1 − λ, 2, y) do for small y, we know that
(
 2 1/y λ ∈ /N  2
V (x)φU,λ (y) ∼ and V (x)φM,λ (y) ∼ 1. (2.16)
1 λ∈N

6
 2.2. ONE-DIMENSIONAL HYDROGEN

This means that as ε → 0, hV i → ∞ on [ε, b] for φU,λ (y) unless λ is a positive

integer. Loudon did not take into account the finiteness of the potential energy
and instead required finiteness of the derivatives of the wave functions. While
this leads to the same conclusion and may in fact be equivalent (a finite potential
energy means a finite kinetic energy, which involves the momentum operator and
thus the derivative), he did not give any arguments about why it is necessary for
the derivative to be finite.
We can thus conclude from square-integrability and the finiteness of the potential
energy that only for λ =: n ∈ N we are able to satisfy the restrictions on the wave
function.
It happens to be that φU,n (y) and φM,n (y) are linearly dependent for all integer n,
so we need to find a second linear independent solution φen (y) as we are dealing with a
second order differential equation. Because we already know one solution, φU,n (y), we
can try reduction of order. Because the differential equation contains no term involving
φ0 (y), this comes down to solving the integral
Z y
1 1
φn (y) =
e ds, s > 0, (2.17)
φU,n (y) φU,n (s)2

which is analytically difficult, if not impossible. In the literature there is no mention of

the linear dependence of φU,n (y) and φM,n (y) and thus also no correct second solution.
One can, however, use Mathematica to find φen (y), but only for fixed integer n (so
not for all integer n simultaneously). It happens to be a combination of powers of y,
exponentials of y and the exponential integral Ei(y) and it looks like that this solution
is divergent for every n. A quick check of finding the first few terms of the Taylor
d
expansion of dx φen (y)/φU,n (y) and comparing these with those of 1/φU,n (y)2 confirms
that this also would be the solution found by reduction of order. Because we want our
wave function to be square-integrable, we can discard this solution and find

φn (y) = φU,n (y) = y e−y/2 U (1 − n, 2; y), y > 0, (2.18)

for our solution of (2.11) for y > 0. One could have also taken M (1 − n, 2; y instead of
U (1 − n, 2; y) as they are equivalent for integer n.
So far we actually have only solved the radial equation of the three-dimensional
hydrogen atom with the azimuthal quantum number l = 0. Recall, see for example [17],
that the wave function can be separated in a radial part R(r) and a angular part Ylm (θ, φ)
and that by the transformation u(r) = rR(r) the radial equation becomes

~2 d2 u q2 1 ~2 l(l + 1)
 
− + − + u = Eu, r > 0, (2.19)
2m dr2 4πε0 r 2m r2

which is exactly (2.6) when l = 0. The difference with the one-dimensional problem
is that the three-dimensional problem is only defined for r > 0, while we also need a

7
CHAPTER 2. PHYSICAL MODEL

solution for x < 0. According to (2.11), the space of solutions for fixed λ is invariant
under the transformation y 7→ −y, i.e.

h y 7→ φ(y) i = h y 7→ φ(−y) i. (2.20)

Assuming this space is one-dimensional, we have,

∃α ∈ C : ∀y ∈ R, φ(−y) = αφ(y) ⇒
(2.21)
φ(y) = φ(−(−y)) = αφ(−y) = α2 φ(y) ⇒ α = ±1.

This means that φ(y) has an undetermined symmetry around y = 0. The tricky part is
finding the correct α.

2.2.2 Different views

Loudon [10] sees the problem as a limit of a regularized potential, e.g.

q2 1
− , (2.22)
4πε0 |x| + ε

for very small ε > 0, which is nonsingular. For this potential there exist both even and
odd wave functions, for which he concludes that in the limit ε → 0 these will remain to
exist. Andrews [9,11] shows that this singularity (in fact any non-integrable singularity)
in the potential acts as an impenetrable barrier so that the wave function must vanish
there, think of an infinite square well of no width. He also shows that the probability
current or particle flux,

~2
 
∗ d d ∗
j= ψ ψ−ψ ψ , (2.23)
2mi dx dx

is zero everywhere, because we only have one linear independent solution, that is real.
From this and the fact that the equation is symmetrical around x = 0, he concludes that
there is no physical connection between the regions x > 0 and x < 0 and that if one
nevertheless regards them to be a single system, then it is doubly degenerate. As there
is no point in connecting the wave functions, both α = +1 and α = −1 are correct, in
agreement with Loudon.
There has also been argued that only the odd wave functions are correct. While
Abramovici and Avishai [12] just state that by taking α = −1 the solutions and their
derivatives are continuous (supposedly natural), Xianxi et al. [13] come with an argument
based on Schrödinger’s equation. They say that integrating the equation over an interval
[−ε, ε] with ε > 0 gives

2m ε
Z
0 0
ψ (ε) − ψ (−ε) = 2 (V (x) − E) ψ(x) dx, (2.24)
~ −ε

8
 2.2. ONE-DIMENSIONAL HYDROGEN

Figure 2.2: Both the even and odd versions of the first two unnormalized wave functions
in one-dimensional hydrogen, n = 1 is solid line and n = 2 is dashed.

so that in the limit ε → 0,

Z 0+
0 + 0 − 2m
ψ (0 ) − ψ (0 ) = 2 V (x)ψ(x) dx. (2.25)
~ 0−

Because ψ(x) is known from (2.18) and (2.21) (ignoring the linear transformations for
now) it can be seen that the odd solutions (α = −1) do and the even solutions (α = 1)
do not satisfy (2.25), meaning only the odd solutions are correct. However, this is only
valid if one considers (2.6) on the entire real axis and not, like what is done here, only for
positive x. As Andrews states, there is no physical connection between the two regions
and thus no point in connecting the two domains. By following the approach of Xianxi
et al., one poses unnecessary extra restrictions on the problem.

From here on both even and odd wave functions are regarded as correct, unless we
are able to show a counterexample. The energies are, nevertheless, independent of α
and are determined by the fact that λ = n ∈ N, so that

mq 4 Ry
En = − 2 2 2 2
=− 2, (2.26)
32π ε0 ~ n n

with Ry the Rydberg constant of energy equal to 13.61 eV. This is exactly the spectrum of
the three-dimensional hydrogen atom, because, as stated before, our problem is actually
just the radial part of three-dimensional hydrogen, equation (2.19). The wave function

9
CHAPTER 2. PHYSICAL MODEL

Figure 2.3: The odd versions of the unnormalized wave functions for n = 2, 3, 4, 5 inside
the Coulomb potential V (x) and around their respective energies En .

can be obtained by applying the transformations (2.9) backwards on (2.18),

(
cn kn x e−kn x U (1 − n, 2; 2kn x) x≥0
ψn (x) = n = 1, 2, . . . , (2.27)
±cn kn x ekn x U (1 − n, 2; −2kn x) x < 0,
√
with wave number kn = −2mEn /~ and normalisation constant cn which follows from
the normalisation of the probability distribution,
Z ∞
|ψ(x)|2 dx = 1, (2.28)
−∞

so that,
√
2
cn = . (2.29)
n!
To retrieve the even solution, one has to take the plus sign in (2.27) and subsequently the
minus sign for the odd solution. In Figure 2.2 both the even and odd (unnormalized) wave
functions for n = 1 (solid line) and n = 2 (dashed line) are shown. The (unnormalized)
wave functions ψn for n = 2 to n = 5 are shown around the energies En in Figure 2.3
inside the Coulomb potential V (x) in eV. They extend slightly outside the potential (as
expected) and they also look like the wave functions of three-dimensional hydrogen with
ψn (x) having n − 1 nodes. ψ1 is not shown as by En ∼ n−2 has an energy four times as
large as ψ2 and therefore lies too far beneath this figure.

10
 2.2. ONE-DIMENSIONAL HYDROGEN

2.2.3 Anomalous states

Abramovici and Avishai [12] have also studied the one-dimensional Coulomb problem
and claim to have found another set of bound states, which they call anomalous states.
The first part is similar to what we have done here, but then they introduce a new
interlacing energy spectrum E en = En+1/2 = − Ry 2 . All they look for are square-
(n+1/2)
integrable and continuous wave functions without any extra condition at x = 0. But
in that case every λ should be allowed, not only integer and half-integer, as ψλ (x) =
cλ kλ |x|e−kλ |x|/2 U (1 − λ, 2, kλ |x|) is continuous and square-integrable for every λ. All
states with λ ∈ / N have infinite expectation values for the potential and kinetic energies
and it even happens to be that these anomalous states are not orthogonal so that the
Hamiltonian is not Hermitian (although it still has real eigenvalues). According to
Núñez et al. [18] this happens because the domain on which the Hamiltonian is defined
is not specified and extra care has to be taken for the singularity, concluding that these
anomalous states cannot be accepted as solutions to the problem, which is in agreement
which Andrews and what we have said here.

11
CHAPTER 2. PHYSICAL MODEL

12
3. Numerical methods

The single ion single electron problem lies at the end of the analytical approach. While
the single electron in the vicinity of two ions case still involves a linear differential
equation, its solutions cannot be expressed in known functions. It is therefore necessary
to turn to numerical methods. We start by developing two slightly different methods
for the single well case and try to replicate the analytical results, after which we extend
this to the double well case. It turns out that these two methods, of which one uses the
approach of Andrews, produce drastically different results.

3.1 Single well

As stated, we will start by replicating the analytical results of a single electron in a 1D

single Coulomb well potential. The Hamiltonian can by (2.6) be written as

~2 d2 q2 1
Ĥ = − − . (3.1)
2m dx2 4πε0 |x|

We want to find solutions to the Schrödinger equation, Ĥψ = Eψ, that are in L2 (R)
and have a finite expectation value for the potential energy, the same restrictions as
for the analytical approach. The finite difference method will be used to discretize our
Hamiltonian Ĥ to a matrix H and our wave function ψ(x) to a vector ψe and then find
the eigenvalues and eigenvectors of the eigenvalue problem

H ψe = E ψ.
e (3.2)

We can approximate the second derivative of ψ(x) for small ∆x by

d2 ψ ψ(x − ∆x) − 2ψ(x) + ψ(x + ∆x)

= + O(∆x2 ) (3.3)
dx2 ∆x2
and dropping the ∆x2 and higher order terms. Sampling ψ(x) on the grid x1 , x2 , . . . , xN
with equal spacing xi − xi−1 = ∆x, i = 2, 3, . . . , N , we define ψe as

ψe = (ψe1 ψe2 . . . ψeN )T = (ψ(x1 ) ψ(x2 ) . . . ψ(xN ))T . (3.4)

13
CHAPTER 3. NUMERICAL METHODS

The finite difference equation becomes

~2  e ei + ψei+1 − q
 2 1 e
− 2
ψi−1 − 2 ψ ψi = E ψei , i = 1, 2, . . . , N, (3.5)
2m∆x 4πε0 |xi |

with ψ0 and ψN +1 determined by the boundary conditions. These equations can be

represented as the matrix-eigenvalue problem in (3.2).
Two slightly different methods will be used to solve (3.5). The first method defines
an N -point grid with x1 = ∆x and xN = L − ∆x, where ∆x = L/(N + 1) and L is
large enough so we can use the boundary condition ψ(L) = 0, which follows from the
fact a smooth and square-integrable function is needed. For the boundary condition at
x = 0 we can use according to Andrews and our prior knowledge that ψ(0) = 0. This
method, which we will call the partition method, thus uses ψ0 = ψN +1 = 0. The second
method, the merge method, does not set the wave function to zero at the singularity and
instead uses a grid with x1 = −L + ∆x and xN = L − ∆x so that ∆x = 2L/(N + 1).
Again, L is large enough so that we can use the boundary conditions ψ0 = ψ(−L) = 0
and ψN +1 = ψ(L) = 0, the same as in the first method. If one takes an odd N , a grid
point is taken directly in the singularity of potential, which leads to a singular matrix.
It is therefore necessary to make N even. For both methods a sparse matrix H can
be constructed using (3.5). We can then solve (3.2) for E and ψe using the MATLAB
function eigs( ).

3.1.1 Energy spectra

First we consider the spectra of both methods, using N = 106 grid points and L = 10 nm.
The numerical results for the eigenvalues and their analytical counterparts are shown in
Table 3.1. This is exactly what we should expect, the spectra of both methods converge
to the analytical spectrum and the errors are relatively small for a moderate amount of
grid points. The partition method is of course more precise as its ∆x is twice as small
(the matrix in the merge method contains the same information twice). The merge
method, however, also gives another interlacing spectrum that goes as approximately
−Ry/(n + 0.119)2 unlike the analytical spectrum that goes as −Ry/n2 . The explanation
of this spectrum will follow after we examine the wave functions.

n 1 2 3 4 5
Analytical -13.6056927 -3.4014232 -1.5117436 -0.8503557 -0.5442277
Rel. err. partition 8.9 · 10−9 2.2 · 10−9 9.5 · 10−10 8.4 · 10−10 2.6 · 10−10
Rel. err. merge 1.1 · 10−7 6.2 · 10−8 4.4 · 10−8 3.4 · 10−8 2.7 · 10−8

Table 3.1: First five analytical eigenvalues E of the Hamiltonian in eV and the relative
errors of the numerically calculated eigenvalues (both methods with N = 106 grid points
and boundaries at L = 10 nm).

14
 3.1. SINGLE WELL

Figure 3.1: Comparison of the partition and merge method with the analytical results
(numerical) (x)−ψ (analytic) (x)
on the basis of the relative change η = ψ |ψ (analytic) (x)|
as a function of x. For
small x it can be seen that the partition method is exactly zero (non-divergent relative
change) unlike the merge method. At large x numerical errors take over as the wave
function approaches zero again. In between the relative changes are small.

3.1.2 Wave functions

The partition method solves for a wave function defined for x > 0, which we can continue
for x < 0 to form either an even or odd function. For the merge method, however, we
only obtain odd wave functions for the spectrum we obtained analytically. As expected
from the accuracy of the eigenvalues, these numerically calculated wave functions are
also accurate, see Figure 3.1 for relative changes given by
ψ (numerical) (x) − ψ (analytic) (x)
η(x) = . (3.6)
|ψ (analytic) (x)|
The partition and merge method have relative changes of the order 10−8 and 10−7
respectively, unless the wave function is (nearly) zero. Unlike the merge method, the
partition method wave function is exactly zero at x = 0, so that the relative change
does not diverge there. The other spectrum in the merge method corresponds to even
wave functions and that can be explained as follows. Without manually setting the
wave function to zero, the discretization of the singular Coulomb potential leads to a

15
CHAPTER 3. NUMERICAL METHODS

Figure 3.2: The differences between the even analytical wave functions and the even
numerical wave functions found by the merge method. The analytical even and odd
wave functions are degenerate, while those wave functions from the merge method have
different spectra due to a lack of incorporated singularity, which also causes the wave
function to differ (that is nonzero at x = 0).

regularized non-singular Coulomb potential. The information of the singularity is not

contained inside the discretization matrix and can only be achieved if the number of
grid points is infinite. Because the regularization is governed by the density of the grid,
which is relatively high in areas where the gradient of the potential is low, we expect
these even wave functions to closely resemble the analytical wave function except maybe
around x = 0. In Figure 3.2 we compare the even wave functions of the merge method
and the analytical wave functions and see that they differ quite substantially. By not
being zero at x = 0 (this might be difficult to see, but their values are around 0.1 in the
figure) the even numerical wave functions gain relatively much energy from the potential,
causing the defect of 0.119 in the eigenvalues. It is important to understand that even for
∆x = 0.02 pm this numerical regularization has a tremendous effect in the behaviour of
the even wave functions. The degeneracy between even and odd wave functions is purely
analytical, that is when there is a singularity, which cannot be numerically replicated in
this fashion.

16
 3.2. DOUBLE WELL

3.2 Double well

We are now ready to solve Schrödinger’s equation for a single electron in the vicinity
of two ions, that is with the Hamiltonian of (2.5) with N = 2. This problem is still
symmetric around x = 0, but the singularities are now positioned at x = ±a instead of
x = 0. First we will try to predict what happens when we bring two ions from infinity
together, because the wave function of a single electron in the potential of two infinitely
separated ions is known: The wave function consists of two single ion wave function (as
deduced for one-dimensional hydrogen) that are centered around x = ±a. Classically,
if the electron is in the vicinity of one ion, it does not feel the potential of the other
ion and vice versa. Note that this state is eight (four regions, so 24 divided by two for
a switch of sign) times degenerate according to Andrews and two (either even or odd
around x = 0) times degenerate if only considering locally odd wave functions.
When a is finite, there exists a smallest n0 such that ∀n > n0 the 1D hydrogen wave
function still oscillates for x > a, meaning that there is overlap in the wave functions
and we cannot just add them up to find the wave function for the two well problem
(these are the quantum plasma solutions). For n << n0 , however, this is still a good
approximation as the wave functions are still bunched around their ion and do not really
feel the other ion. For smaller a n0 decreases, as the ions influence each other more. It
is important to understand that by decreasing a the symmetries and general behavior
of the spectrum are retained. This means that every state gradually changes its energy
and wave function; no abrupt changes are expected when bringing two ions from infinity
together. This enables us to ‘follow’ the eigenstates as we bring down a so that we expect
to see some kind of band structure that accounts for the degeneracies we see at infinity.
For this approach it physically does not matter if we take the eight or two degeneracy
approach, because it is actually |ψ| that influences the (kinetic and potential) energy,
which is always two times degenerate.
We can now use both methods to calculate the eigenvalues and eigenfunctions and see
if our expectations are correct and what Andrews’ argument means for this locally asym-
metric potential. Again, a matrix for the discretized Hamiltonian can be constructed for
the two well problem that MATLAB can use to find its eigenvalues and eigenfunctions.
We start with the merge method to postpone the subtleties concerning the partition
method.

17
CHAPTER 3. NUMERICAL METHODS

Figure 3.3: First two eigenstates for the two well problem with separation a = 2 nm,
using the merge method. They resemble even and odd combinations of the odd solutions
in the single well case as in Figure 2.2 (albeit with slightly lower energies).

Figure 3.4: The thirteenth and fourteenth eigenstate for the two well problem with
separation a = 2 nm, using the merge method. There would be overlap between the
single well wave functions so that approximating the solutions for the double well case by
taking a linear combination of the single well solutions is no longer a valid approximation.

18
 3.2. DOUBLE WELL

3.2.1 Merge method

The merge method is as simple as choose a desired a, define a grid and it gives the entire
eigenstate (wave function on the entire domain and its eigenenergy). As the precision
of both methods for the single well was good, this is no longer a point to worry about.
Let a first be large, say 2 nm, so that the ground states of the single well problem
(approximately 0.5 nm broad) do not overlap. The expected band structure can be
traced to finding pairs (because we have two ions) of wave functions, one even and one
odd around x = 0. The ground state and first excited state of the double well problem
can be seen in Figure 3.3. As expected, they are the even and odd combinations of the
single well solutions centered around x = ±a. This is nearly a degenerate state, the
even state has a slightly lower energy, that has an energy just below the ground state
of the single well problem (−13.61 eV). The thirteenth and fourteenth excited state can
be seen in Figure 3.4 and are clearly not the sum of the single well wave functions.
There are now big differences between the even and odd wave function, that for large a
were degenerate. For one their energies are now relatively far apart and furthermore the
nodes of the wave functions shift more apart the further they are from the singularity.
At the singularity both are still approximately zero. The same thing also happens to the
ground state when a is smaller (shown in Figure 3.5 for a = 0.35 nm). In this figure both
the single well solutions as the effect of lifted degeneracy are visible. For the interlacing
spectrum with even single well solutions the same results are found (shown in Figure
3.6 for a = 0.35 nm). Because the quantum plasma states in Rydberg crystals, which
have interionic distances a in the order of µm, have a very high quantum number n, it is
useful to study solutions around the hydrogenic distances in the order of nm, for which
there exist quantum plasma states with very low n, as the behaviour is similar.

19
CHAPTER 3. NUMERICAL METHODS

Figure 3.5: First two eigenstates for the two well problem with separation a = 0.35
nm, using the merge method. Just as in Figure 3.4, but unlike Figure 3.3, there is a
slight overlap between the single well ground state wave functions. The nuclei are closer
together now so that the overlapping of wave functions also happens for states with lower
energies.

Figure 3.6: First two even eigenstates for the two well problem with separation a = 0.35
nm, using the merge method. They are both similar to those in Figure 3.2 (nonzero
at the singularity) and those in Figure 3.5 (slight perturbations from single well wave
functions)

20
 3.2. DOUBLE WELL

Figure 3.7: The behaviour of a part of the energy spectrum in the double well problem
using the merge method as a function of the separation a. The legend includes the local
symmetries (L) around the singularities x = ±a and the global symmetries (G) around
the centre, x = 0, for the different states. Both can independently be even (e) or odd
(o). For analysis of the figure, refer to the text.

The behaviour of the energy levels as a function of a is very interesting and can
be seen in figure 3.7. It is important to note that even though the Hamiltonian is
spin-independent, spin effects must be considered. Let the two electron wave function
(2)
ψnn0 (x1 , x2 ), where xi is the position of electron i with i ∈ {1, 2}, be an eigenstate of
the Hamiltonian

H (2) = H1 + H2 =
~2 ∂ 2 q2 ~2 ∂ 2 q2
   
1 1 1 1
− − + − − + ,
2m ∂x21 4πε0 |x1 − a| |x1 + a| 2m ∂x22 4πε0 |x2 − a| |x2 + a|
(3.7)

that lacks electron-electron repulsion. As this Hamiltonian can be split into the single
electron Hamiltonians H1 and H2 so that Hi only depends on xi , the wave function
(2)
separates into ψnn0 (x1 , x2 ) = ψn (x1 )ψn0 (x2 ) and thus is an element of L2 (R) ⊗ L2 (R) ∼
=
2 2 (2)
L (R ) and has energy Enn0 = En +En0 . Actually we need to take the linear combination

(2) 1
ψnn0 (x1 , x2 ) = √ (ψn (x1 )ψn0 (x2 ) − ψn0 (x1 )ψn (x2 )) , (3.8)
2
as electrons are fermions [17]. Because the Hamiltonian is independent of spin, the

21
CHAPTER 3. NUMERICAL METHODS

Figure 3.8: A visualization of the greater influence of changing the well separation a on
the inner region of the double well Coulomb potential V (x) than on the outer region.

spinor also separates, resulting in a total (spatial and spin) wave function
(2)
ψnn0 (x1 , x2 )χ(s1 , s2 ). (3.9)

Returning to the energies in Figure 3.7, we can see the degenerate energies for large
a that for very large a will converge to the single well energies. Then as we lower a
we see that the degeneracy is lifted, which happens for larger a (ψ9 , red) when the
energy is high and thus the wave functions are wider compared to lower energies (ψ5 ,
dark blue). Before the minimum is reached, the energies of globally even wave functions
decrease more rapidly than those of the globally odd wave functions. This effect can be
contributed to the spin of the electron. The total wave function, including spin, always
has to be odd according to Pauli, so if the wave function is globally odd, then the spin
is in an even triplet state, for which the exclusion principle generates an extra repulsion.
This repulsion then rises the energy levels of globally odd wave functions. On the other
hand, globally even wave functions have their electrons in a singlet state, meaning they
are not subject to the exclusion principle. There is also a difference between locally even
and odd wave functions. From Figures 3.5 and 3.6 we can deduce that the locally even
wave functions are more concentrated in the outer regions, |x| > a, while locally odd wave
functions are concentrated more in the inner region, |x| < a. The effects of decreasing a
has the most effect on V (x) in the inner region (see Figure 3.8) so that the energies of
the locally odd wave functions change more rapidly. Finally, after a local maximum all
energy levels drop again, which is hard to explain in our model. The model lacks any
kind of nucleus-nucleus repulsion (Coulomb nor strong nuclear) as we are working in the
Born-Oppenheimer approximation and thus does not represent the physical system we
want to study anymore for such small a. We can this ignore this behavior, although it
is interesting to see that the globally even and odd wave functions become degenerate
again for very small a while their wave functions retain their symmetry and are very
different.

22
 3.2. DOUBLE WELL

Figure 3.9: First four eigenstates for the two well problem with separation a = 0.35 nm
in the region −a < x < a, using the partition method. The wave functions are exactly
zero at x = ±a.

3.2.2 Partition method

Andrews’ argument states for the partition method that the wave function is zero at x =
±a and that there is no physical connection between the any of the three regions x < −a,
−a < x < a and x > a. We will first look at the eigenstates in the inner region. The first
four wave functions and their eigenenergies are given in Figure 3.9 for a = 0.35 nm. The
wave functions are in turn even and odd, zero at the boundaries and similar to the wave
functions in the inner region of the merge method. It is interesting to note that there
exist only a finite amount of bound states in this inner region, #ψn := |{n ∈ N|En < 0}|,
and its dependence on a is shown in Figure 3.10. From the numerical method we find
that for a < 0.027 nm even the ground state is an unbound state and thus by fitting the
data in the figure with #ψn = A(a − 0.027)B we find that B ≈ 1/2. The reason for a
finite #ψn is that in order for a smooth and normalized function to have more nodes on
a bounded interval a more rapidly changing derivative is needed. This corresponds to a
bigger second derivative and thus kinetic energy T , which is always positive and thus can
bring the total energy E = T + V to values above zero. On a (semi-)infinite interval the
wave function always has more space to ‘expand’. The reason that the merge method
has an infinite amount of solutions is that they, on closer inspection, are not zero at
the singularity. Although one can recognize the inner partition solutions in the merge
solutions, they are altered and even more so for higher energies.

23
CHAPTER 3. NUMERICAL METHODS

Figure 3.10: The amount of bound states #ψn := |{n ∈ N|En < 0}| in the inner
region, using the partition method, as a function of a. The data has been fitted with
#ψn = A(a − 0.027)B , where A = 7.8 and B = 0.52 ≈ 1/2.

The solutions on the outer regions are by itself less interesting. In Figure 3.11 it
can be seen that for a = 0.35 nm and a = 2 nm the wave functions are very similar
to the single well case. And for the energies we can once again say that they approach
the single well energies as a increases. As expected, the outer regions have an infinite
amount of eigenstates and thus energies, but they are different from the energies in the
inner region. This means that whereas we could see the single well case as a single
system because of symmetry (thus both sides had equal energies), we now lack any kind
of connection between the inner and outer regions. By requiring a finite potential energy,
and thus a vanishing wave function at the singularity, not only the wave functions are
not connected, but also their energies in this locally asymmetric potential. This is in
favour of what Andrews has stated. One could, nevertheless, smoothly attach the wave
functions in the inner and outer regions for comparison with the merge method, but
then this state has no specified energy on the full domain. The relative change between
the ground states at a = 0.35 nm,
(part) (merge)
ψ1 (x) − ψ1 (x)
η(x) = (merge)
, (3.10)
|ψ1 (x)|
is shown in Figure 3.12, from which we can see that the merge and partition method
differ by a maximum of seven percent, except when the wave function approaches zero
(at the singularities and for large |x|). This is a significant difference in the behaviour
of the electron and apparently causes the splitting of regions, while both satisfy the
numerical Schrödinger equation. For the merge method the node near the singularity is
just a node and has no other physical consequence than that there is a low probability
of finding the particle near that point. For the partition method, on the other hand,
the node in the singularity represent an impenetrable barrier that separates two regions
with completely different energies.

24
 3.2. DOUBLE WELL

Figure 3.11: First three eigenstates for the two well problem with separations a = 0.35
nm and a = 2 nm in the region x > a, using the partition method. The wave functions
are exactly zero at x = a and their shape changes very little for these a (they look like
those in Figures 2.2 and 2.3).

Figure 3.12: Comparison of the partition and merge method for the ground state at
(part) (merge)
ψ1 (x)−ψ1 (x)
a = 0.35 nm on the basis of the relative change η(x) = (merge) as a function
|ψ1 (x)|
of x. The relative change is at most seven percent, unless the wave function is nearly
zero.

25
CHAPTER 3. NUMERICAL METHODS

26
4. Regularization

So far we have used two different methods to establish a 1D toy model for Rydberg
crystals, which gave slightly different results for a single ion, but vastly different results
for two ions. It is the one-dimensionality of this toy model that makes the potential
singularities act as impenetrable barriers, just as Andrews predicted. But we are actually
dealing with a physically three-dimensional problem that allows the electron to move
around the singularity due to angular components. The partition method, which sets
the wave function to zero at the singularity, is therefore unreasonable as it disconnects
regions that should not be disconnected. The merge method is more applicable to the
physical problem, but is mathematically crude, as it actually involves regularization of
the potential that is controlled by the numerical grid. The next section will establish a
numerical method for an externally controlled regularization of the Coulomb potential,
so that we stay close to the original problem.

4.1 Single electron

We introduce the regularized potential for small regularization parameter ε as

q2 1
Vε (x) = − √ , x ∈ R. (4.1)
4πε0 x2 + ε2

This is a nonsingular, continuous, and smooth potential resembling the Coulomb poten-
1 1
tial. This potential is preferred over potentials like |x|+ε and |x+ε| because of smoothness.
For two different ε of 0.03 and 0.08 nm the potential with two of these wells is shown in
Figure 4.1 as a function of x. We can see that ε has a tremendous effect on the depth
of the wells, so it is necessary to use a very small ε. But the resolution of our grid sets
a limit, because if ε is not much larger than the distance between grid points ∆x, the
effects of regularization are negligible. We will study the behaviour around ε = 10−4
nm, which corresponds to approximately 20∆x for the used grid.
First we investigate what happens to the eigenstates when ε is varied. For an ε of
10−5 ,10−4 , and 10−3 nm and an interionic distance a of 0.3 nm the ground state wave
function, ψ1 , is shown in Figure 4.2. They are similar to the result we found in Figure
3.5 and vary little for these changes in ε. The energies decrease as the regularization

27
CHAPTER 4. REGULARIZATION

Figure 4.1: Effect of varying the regularization parameter ε on the regularized Coulomb
potential. A smaller ε greatly increases the depth of the wells.

parameter gets smaller, because the potential well becomes deeper. The locally even
states are more interesting and are shown for the same ε and a in Figure 4.3. The wave
functions now vary more noticeably, especially in the centre. The peaks of the wells
decrease as ε decreases, meaning that in the limit ε → 0 the locally even wave function
probably go to zero at the well, just as Loudon and Andrews predicted. The energies
now also change drastically with again a decrease in energy for smaller ε. They will most
likely converge to the values of the locally odd states so that the predicted degeneracy
is achieved.
For the energies we also expect the results to match the results in the merge method
of the last section. For a regularization parameter ε = 10−4 nm the dependence of the
energies of the first four eigenstates on the interionic distance a is shown in Figure 4.4.
There only little differences, such as the fact there is less case of degeneracy at low a,
when compared to the results from before (although different eigenstates are shown this
time, the behaviour is the same). The reasoning of those results is still valid.

28
 4.1. SINGLE ELECTRON

Figure 4.2: The little effects of varying the regularization parameter ε on the ground
state of the double well problem with an interionic distance a = 0.3 nm.

Figure 4.3: The large effects of varying the regularization parameter ε on the first locally
even state of the double well problem with an interionic distance a = 0.3 nm.

29
CHAPTER 4. REGULARIZATION

Figure 4.4: The behaviour of a part of the energy spectrum in the regularized double
well problem as a function of the interionic distance a. The figure and its explanation
are similar to Figure 3.7

30
 4.2. TWO ELECTRONS

4.2 Two electrons

A more realistic system of Rydberg crystals is the double well double electron problem,
meaning we are dealing with two complete Rydberg atoms instead of one being ionized.
The solution to the problem without electron-electron repulsion can be constructed from
the wave function in (3.8). Incorporating the repulsion can be done in several ways, for
example perturbative and numerical, which are discussed in this section.

4.2.1 Perturbational

We can use perturbative methods to add electron-electron repulsion in our model for
states in which the wave functions do not overlap as then the effects of this repulsion is
small and can been seen as a perturbation. The Hamiltonian without repulsion is the
two electron Hamiltonian H (2) as in (3.7) and the perturbation is given by

q2 1
H0 = , (4.2)
4πε0 |x1 − x2 |
(2)
which is strictly positive. These work on the wave function ψnn0 (x1 , x2 ) as in (3.8). The
first and second order corrections for the eigenstates of H (2) + H 0 are given by [17]
(2) (2)
(2) (2)
X |hψmm0 |H 0 |ψnn0 i|2
1
Enn0 = hψnn0 |H 0 |ψnn0 i > 0 and 2
Enn 0 = (2) (2)
< 0. (4.3)
m6=n,m0 6=n0 Enn0 − Emm0

Although everything is known to calculate these corrections, they involve an integral

over R2 with an integrand that is singular on (nearly) the whole x1 = x2 line. Without
regularizing one can use for example Mathematica’s function NIntegrate, which cannot
calculate this. With regularization the results are heavily dependent on the regulariza-
tion parameter and seem to diverge as the parameter goes to zero, meaning we can never
be certain about the answer.

4.2.2 Numerical

Another approach is to incorporate the repulsion term in our numerical Hamiltonian

matrix. This means that we have to discretize our wave function on a two dimensional
grid, which also causes the the matrix to be less sparse. Whereas for the one electron
case the distance between grid points ∆x was in the order of 10−5 nm, ∆x is in the order
of 10−2 nm for two electrons, a severe drop in resolution. With a regularized repulsion
potential as in (4.1) it is possible to find the eigenstates in the usual manner, but the
energies are dependent of the amount of grid points so that they are meaningless except
for the order they are in. For an interionic distance a = 0.3 nm and a regularization
parameter ε = 10−4 nm the ground and first excited state are shown in Figures 4.5 and

31
CHAPTER 4. REGULARIZATION

Figure 4.5: The ground state for a two electron system with an interionic distance a = 0.3
nm and regularization parameter ε = 10−4 nm. Because of repulsion the wave function
is nearly zero (green) for x1 ≈ x2 . Furthermore, the wave function is symmetric under
the change of electrons.

4.6. We see immediately that the wave function is zero (green) when x1 ≈ x2 , meaning
the repulsion causes the electrons to avoid each other. If one electron is near the well at
−a then the other electron is at the well at +a and vice versa. Because the electrons are
identical we have the symmetric and antisymmetric wave function, of which the former
has a lower energy due to it having a singlet spin state as explained before.
More insight can be gained from higher excited states with the same parameters, for
example the ninth excited state that is shown in Figure 4.7. This state shows that the
wave function is zero if both |x1 | and |x2 | are larger than |a|. If they have the same sign
this is of course caused by the electron-electron repulsion and if they have different signs
the system can gain energy by bringing one electron between the wells as this results
in more energy from the ions than it loses from repulsion. Similarly, if one electron is
between the two wells, then the other electron has a greater probability of being found
outside of the wells, which is also slightly larger for the well that is furthest from the
former electron (seen from the triangle like shapes just outside the square with vertices
at (±a, ±a)). Lastly, if one electron is, again, far outside either well, then the other
electron has a relatively unperturbed (in this case odd) ground state single well wave
function around the other well.

32
 4.2. TWO ELECTRONS

Figure 4.6: The first excited state for a two electron system with an interionic distance
a = 0.3 nm and regularization parameter ε = 10−4 nm. Because of repulsion the wave
function is nearly zero for x1 ≈ x2 . Furthermore, the wave function is antisymmetric
under the change of electrons.

Figure 4.7: The ninth excited state for a two electron system with an interionic distance
a = 0.3 nm and regularization parameter ε = 10−4 nm. For analysis of the figure, refer
to the text.

33
CHAPTER 4. REGULARIZATION

34
5. Conclusion

The main idea was to establish a one-dimensional toy model for energy levels Rydberg
crystals based on Schrödinger’s equation. We use the Born-Oppenheimer approximation
and represented Rydberg atoms as an ion (shielded nucleus) plus an electron, so that
the model is nothing more than an electron in a lattice of point charges. The single
atom case is exactly the same as one-dimensional hydrogen, a controversial subject. On
a half-plane, so on either side of the singularity in the Coulomb potential, the situation
is similar to the radial part of three-dimensional hydrogen, as, among other things,
the energies go as n−2 . The problem, however, lies at connecting the two sides. The
strongest arguments are given by Loudon [10] and Andrews [9, 11], who say that there
exist both even and odd eigenstates around the singularity, that have the same energy.
Andrews actually states that the problem on the entire axis is not a single system and
the two sides must be viewed independently. The singularity acts as an impenetrable
barrier that does not allow an electron to transfer from one side to the other.
In order to study the system for more ions numerical methods were used, with which
we started to replicate the analytical results with two different methods, the partition
and merge method. In the partition method we manually set the wave function to zero
at the singularity, so that we can form the even and odd solutions ourselves. The merge
method, on the other hand, only has boundary conditions at infinity and uses a potential
that is actually not singular. The odd solutions are the same as in the partition method,
but the even solutions are not. They are nonzero at the singularity and have different
energies. It is expected that for an infinitely small grid the even wave functions will
converge to the analytical solutions.
The situation gets complicated for two atoms (of which one is ionized, so still one
electron), for which the two method produce vastly different results. The partition
method still splits the regions by setting the wave function to zero at the singularity,
but now the problem is not symmetrical around the singularity, resulting in different
energies in different regions. For a wide separation the merge method gives all globally
(around the centre) even and odd combinations of the locally (around the singularity)
either even or odd wave functions of the single atom case. When the atoms are brought
together symmetries are retained, but wave functions start to overlap. If not extremely
close together, the behaviour of the energies can be explained by symmetries and spin
effects.

35
CHAPTER 5. CONCLUSION

The partition method does not adequately represent the physical system, as in the
three-dimensional situation there is no barrier. Furthermore, the merge method is math-
ematically crude, as it is in fact a grid-controlled regularization of the Coulomb potential.
By manually controlling the regularization we can see that the odd solutions are more
or less independent of the regularization parameter (if it is small enough), while the
change is even solutions (both wave functions and energies) is large. The convergence
to the degenerate situation is slow but visible. The study of two electron systems with
electron-electron repulsion has not given sufficient results yet due to computational prob-
lems involving grid sizes and integrals involving singularities. This is the main point to
continue the research of this toy model for Rydberg crystals. Another point of continu-
ation lies at increasing the amount wells, which is relatively simple, and the amount of
electrons, which is more difficult as it also increases the dimensionality of the problem.

36
 BIBLIOGRAPHY

Bibliography

[1] R.P. Feynman, Simulating physics with computers, International Journal of Theo-
retical Physics, 21, 467-468, 1982.
[2] D. Deutsch, Quantum theory, the Church-Turing principle and the universal quan-
tum computer, Proceedings of the Royal Society A, 400, 97-117, 1985.
[3] P.W. Shor, Polynomial-Time Algorithms for Prime Factorization and Discrete Log-
arithms on a Quantum Computer, SIAM Journal of Computing, 26, 1484-1509,
1997.
[4] C. Pomerance, A Tale of Two Sieves, Notices of the AMS, 43, 1473-1485, 1996.
[5] T.H. Johnson, S.R. Clark, D. Jaksch, What is a quantum simulator?, arXiv,
1405.2831, 2014.
[6] R. van Bijnen et al., Adiabatic formation of Rydberg crystals with chirped laser
pulses, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 2011.
[7] T.F. Gallagher, Rydberg Atoms, Cambridge University Press, 1994.
[8] D. Comparat, P. Pillet, Dipole blockade in a cold Rydberg atomic sample, Journal
of the Optical Society of America B, 27, A208, 2010.
[9] M. Andrews, Singular potentials in one dimension, American Journal of Physics,
44, 1064, 1976.
[10] R. Loudon, One-Dimensional Hydrogen Atom, American Journal of Physics, 27,
649, 1959.
[11] M. Andrews, The one-dimensional hydrogen atom, American Journal of Physics,
56, 776, 1988.
[12] G. Abramovici, Y. Avishai, The one-dimensional Coulomb problem, Journal of
Physics A: Mathematical and Theoretical, 42, 285302, 2009.
[13] D. Xianxi et al., Orthogonality criteria for singular states and the nonexistence of
stationary states with even parity for the one-dimensional hydrogen atom, Phys.
Rev. A, 55, 2617, 1997.

37
BIBLIOGRAPHY

[14] M. Marinescu, H. R. Sadeghpour, and A. Dalgarno, Dispersion coefficients for

alkali-metal dimers, Physical Review A, 49, 982, 1994.
[15] PhET Interactive Simulations, Quantum Bound States,
http://phet.colorado.edu/en/simulation/bound-states.
[16] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas,
Graphs, and Mathematical Tables, New York: Dover, 504, 1965.
[17] D.J. Griffiths, Introduction to Quantum Mechanics, Second Edition, Pearson Pren-
tice Hall, 2005.
[18] H.N. Núñez-Yépez, A.L. Salas-Brito, D.A. Solı́s, Comment on ‘The one-dimensional
Coulomb problem’, Journal of Physics A: Mathematical and Theoretical, 46, 208003,
2013.

38