Anderson localization in weakly interacting Bose gases

(1)

SECOND CYCLE, 30 CREDITS STOCKHOLM SWEDEN 2020 ,

A numerical investigation of

Anderson localization in weakly interacting Bose gases

CRYSTAL UGARTE

KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF ENGINEERING SCIENCES

(2)

(3)

A numerical investigation of Anderson localization in

weakly interacting Bose gases

CRYSTAL UGARTE

Degree Projects in Scientific Computing (30 ECTS credits)

Degree Programme in Applied and Computational Mathematics (120 credits) KTH Royal Institute of Technology year 2020

Supervisor at KTH: Patrick Henning Examiner at KTH: Michael Hanke

(4)

TRITA-SCI-GRU 2020:026 MAT-E 2020:007

Royal Institute of Technology School of Engineering Sciences KTH SCI

SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

(5)

The ground state of a quantum system is the minimizer of the total energy of that sys-

tem. The aim of this thesis is to present and numerically solve the Gross-Pitaevskii eigen-

value problem (GPE) as a physical model for the formation of ground states of dilute Bose

gases at ultra-low temperatures in a disordered potential. The first part of the report in-

troduces the quantum mechanical phenomenon that arises at ground states of the Bose

gases; the Anderson localization, and presents the nonlinear eigenvalue problem and the

finite element method (FEM) used to discretize the GPE. The numerical method used to

solve the eigenvalue problem for the smallest eigenvalue is called the inverse power iteration

method, which is presented and explained. In the second part of the report, the small-

est eigenvalue of a linear Schrödinger equation is compared with the numerically computed

smallest eigenvalue (ground state) in order to evaluate the accuracy of a linear numerical

scheme constructed as first step for numerically solving the non-linear problem. In the

next part of the report, the numerical methods are implemented to solve for the eigen-

value and eigenfunction of the (non-linear) GPE at ground state (smallest eigenvalue). The

mathematical expression for the quantum energy and smallest eigenvalue of the non-linear

system are presented in the report. The methods used to solve the GPE are the FEM

and inverse power iteration method and different instances of the Anderson localization

are produced. The study shows that the error of the smallest eigenvalue approximation

for the linear case without disorder is satisfying when using FEM and Power iteration

method. The accuracy of the approximation obtained for the linear case without disor-

der is satisfying, even for a low numbers of iterations. The methods require many more

iterations for solving the GPE with a strong disorder. On the other hand, pronounced

instances of Anderson localizations are produced in a certain scaling regime. The study

shows that the GPE indeed works well as a physical model for the Anderson localization.

(6)

(7)

Sammanfattning

En numerisk undersökning av Anderson-lokalisering i svagt interagerande Bose- gaser

Syftet med denna avhandling är att undersöka hur väl Gross-Pitaevskii egenvärdesekva- tion (GPE) passar som en fysisk modell för bildandet av stationära elektronstater i utspädda Bose-gaser vid extremt låga temperaturer. Fenomenet som skall undersökas heter Anderson lokalisering och uppstår när potentialfältets styrka och störning i systemet är tillräckligt hög.

Undersökningen görs i denna avhandling genom att numeriskt lösa GPE samt illustrera olika

utfall av Anderson lokaliseringen vid olika numeriska värden. Den första delen av rapporten

introducerar det icke-linjära matematiska uttrycket för GPE samt de numeriska metoderna

som används för att lösa problemet numerisk: finita elementmetoden (FEM) samt egenvärde-

salgoritmen som heter inversiiteration. Finita elementmetoden används för att diskretisera

variationsproblemet av GPE och ta fram en enkel algebraisk ekvation. Egenvärdesalgoritmen

tillämpas på den algebraiska ekvation för att iterativt beräkna egenfunktionen som motsvarar

det minsta egenvärdet. Det minsta egenvärdet av en fullt definierad (linjär) Schrödinger ek-

vation löses i rapportens andra del. Den linjära ekvationen löses för att ta fram en förenklad

numerisk algoritm att utgå ifrån innan den icke-linjära algoritmen tas fram. För att försäkra

sig att den linjära algoritmen stämmer bra jämförs det exakta egenvärdet för problemet med

ett numeriskt framtaget värde. Undersökningen av den linjära algoritmen visar att vi får en

bra uppskattning av egenvärdet - även vid få iterationer. Vidare konstrueras den ickelinjära

algoritmen baserat på den linjära. Ekvationen löses och undersökes. Egenfunktionen som

motsvarar minsta egenvärdet framtas och beskriver kvantsystemet i lägsta energitillstån-

det, så kallade grundtillståndet. Undersökningen av GPE visar att de numeriska metoderna

kräver många fler iterationer innan en tillräckligt bra uppskattning av egenvärdet fås. Å

andra sidan fås markanta Anderson lokaliseringar för ett skalningsområde som beskrivs av

styrkan av potentialfältet i relation till dess störning. Slutsatsen är att Gross-Pitaevskii

egenvärdesekvation passar bra som en fysisk modell för detta kvantsystem.

(8)

(9)

Acknowledgement

Firstly I would like to pay my special regards and appreciation to my supervisor Patrick Henning, who has been the biggest support and who has convincingly guided me to do the right thing throughout my thesis study. Without his persistent support and considerate encouragement to plan my own work, the goal of this project would have been hard to realize. I also appreciate that my supervisor always answered all my questions, despite being trivial at times, with a lot of details and enthusiasm - it made the study extra fun to work on, and made a big difference. I wish to show my profound gratitude to my talented cousin Monika Diaz, who has painted and designed most of the illustrations in this report. Thank you for what you create, your paintings are beautiful.

I would also like to acknowledge the support and great love of my family; my father

Dante, my mother Mirella and my sister Alma-Mia. My parents have always believed in me

and inspired me to do what I love to do, and have also been excellent academic role models

that I honor so much. I give my profound thanks to my sister for cheering on me. She has

passion and discipline for the projects she takes on, and is someone I can rely on for advice

at all times.

(10)

(11)

Anderson localization describes the general phenomenon that waves that propagate through a sufficiently disordered medium are extremely localized in few spots in space. This phe- nomenon can arise in many different types of physical systems since it applies to electromag- netic waves, acoustic waves, quantum waves and spin waves. The most famous example of this is when the electron transportation completely halts in a sufficiently disordered medium.

This thesis is the study of the Gross-Pitaevskii equation (GPE) as a physical model for the formation of stationary states of dilute Bose gases at ultra low temperatures in a sufficiently disordered potential field. Worth to mention is that there is a unique solution for the GPE.

The thesis will aim to reproduce this solution numerically based on the GPE using the Galerkin finite-element method (Galerkin FEM) and an iterative eigenvalue algorithm called inverse power iteration. The solution represents the phenomenon Anderson localization and the aim is thus to evaluate the GPE as a physical model for the formation of Anderson localization. The following physical questions will be answered in this report; How does the Anderson localization depend on the disorder and strength of the potential field applied to the system? What happens to the Anderson localization phenomena when we change the number and the type of particles in it (change the chemical element)? For that reason, there is also an emphasis on quantum mechanical theory in the introduction and the results section of the report to answer the questions posed. Before presenting the mathematical equation and the numerical methods that will be used, the two first chapters of the report will provide a physical background.

Figure 1: Anderson localization. Drawing by Monika Diaz, February 2020

In this section, an introduction to the Anderson localization phenomenon from the basis

(14)

of electric properties of disordered solids is presented. The reason is that the Anderson local- ization phenomenon also manifests in electrons in sufficiently disordered crystals (conductors) with high concentrations of defects. This quantum system was the basis of PW Anderson’s seminal work on the phenomenon that was published in 1958. (See Ad Lagendijk, Bart van Tiggelen, and Diederik S. Wiersma article in Physics today 2019, page 24) The Anderson localization phenomena is illustrated in Fig 1.

In PW Anderson’s ground-breaking research he stated that an increase in the degree of lattice disorder decreases the mean free path of the electron. It even ceases completely when the order of disorder is sufficiently high, changing the behavior of the metal completely. As nearly all electrons become localized to few small regions, the originally conducting metal transform into an insulator. This of course intrigued physicists for decades and has even caused many disputes (See Ad Lagendijk, Bart van Tiggelen, and Diederik S. Wiersma article in Physics today 2019, page 25). The finding was unfeasible to explain before Anderson’s research work and most importantly the discovery of quantum wave propagation. The free electron model (or Drude model) in Classical mechanics assumes that freely moving electrons move randomly across the metal lattice by electrons scattering from immobile positive ions in the metal lattices. It also assumes that electric conductivity is proportional to the mean free path (the average length that an electron travels before it collides with an ion). This neglects the impact of electron interactions. An illustration of the free electron model is shown in Fig 2,

Figure 2: Electron travelling freely across a conductor and scattered by immobile positive ions along the way. Drawing by Monika Diaz, February 2020.

This simple model fails to describe the Anderson localization that arises in disordered

regular metal lattices in a (sufficiently) disordered potential. The Drude model namely

claims that the conductivity is proportional to the mean free path. This is not the case since

electric conductivity halts in sufficiently disordered metal lattices. This generates an absence

of diffusion in the metal all together. The reason that the Drude model is not sufficient is

that it is based on the simplest form of kinetic theory and treats the particles as identical

(15)

solid spheres. For complex systems such as the Anderson localization phenomenon in weakly interacting Bose-gases (at ultra-low temperatures), particle interactions cannot be ignored.

The Gross-Pitaevskii equation (GPE) is the model we will use in this study to describe the behavior of the quantum states of such Bose-gases. This model not only considers the kinetic theory of quantum particles - but also the potential energy and the interactions energy of the quantum particles. To avoid confusion between the report’s Bose-gas quantum state and the example in P.W. Anderson’s work, remark that the GPE cannot be used to described the behavior of electrons. The purpose of the study of the GPE is to evaluate it as a physical model for the quantum states of Bose-gases in particular.

Motivations

It is given that the Anderson localization can be observed in a system of Bose gases in a (sufficiently) disordered potential as the temperature of the gas is lowered close to absolute zero. In ground state, the gas forms a peculiar state of matter called the Bose-Einstein con- densate (BEC) or superfluid state. This state of matter is unlike all other state of matters.

The individual particles in this state of matter overlap and behave as if they where one single super atom that occupies the whole domain the particles are in. This makes the localization phenomenon that arises difficult to explain. The behavior of the system cannot be described with the theories we are most familiar with, such as classical mechanics. Despite this, the behavior can be reproduced using the mathematical model GPE and numerical methods to solve it. This will be demonstrated in this thesis by investigating on the impact of different parameters in the model. The impact of the parameters such as the disorder of the poten- tial field such align with the theory. How does each parameter influence the phenomenon?

Another interesting and strongly connected property of Bose-Einstein condensates is super-

fluidity which describes fluids that flow with zero friction and zero viscosity. The simplest

form of superfluidity can be observed by the hydrogen-4 in ultra-low temperatures (or abso-

lute zero). In the superfluid state the hydrogen-4 behaves in a very unusual manner. It does

not dissipate any energy and can also be described by the Gross-Pitaevskii equation. The

applications of phenomenon are very broad. For example, superfluids may serve as agents

to trap and drastically reduce the speed of light (See Lise Bix’ article in ScienceNordic 2014,

October). A development of far-reaching physics with strong connections to BEC (and ultra-

slow light) is the pulsed-atom lasers with spatially controlled output modes used in atom

holography (Christopher Slowe, Naomi S. Ginsberg, Trygve Ristroph, Anne Gooodsell and

Lene Vestergaard Hau in a scientific paper from Harward Physics department on ’Ultra-slow

lights and Bose-Einstein Condensates’ 2005). Research of BEC and the GPE model is thus

important for both understanding superfluidity and pulsed-atom lasers too.

(16)

1.1 Problem statement

In this thesis we consider the Gross-Pitaevskii eigenvalue problem (GPE) as a model for Bose gases in ground states in a disordered potential at ultra-low temperatures. The quan- tum state represents a Bose-Einstein condensate and can be described by a standing wave function strongly localized at few spots in space. We will reproduce the phenomenon An- derson localization and the shape of the potential numerically using Galerkin FEM. The investigation aims to answer how the Anderson localization behaves depending on

• the disorder of the potential field

• the strength of the potential field

• the number and the type of particles.

As the physical interpretation of the numerical investigation will be used to validate the results of the report, the physical background section will include quantum mechanical con- cepts. It will refreshes the reader on the interpretation of the wave function as well as the existence of integer- and half spins, and the BEC (Sec 2). Furthermore, since numerical schemes will be constructed to model the mathematical problem in this study, a Finite Ele- ment Method background section will be presented next. It will cover the numerical basics needed in Sec 3. An introduction to the inverse power iteration method is provided there as well. The methodology (the approach) for solving the problem step-wise is covered in Sec 5.

The validation section, Sec 7, will present the numerical parameter values that will be used to reproduce different instances of the Anderson localization phenomena. As will be described in section 5, we will first solve a linear case of the Schrödinger equation before solving the non-linear GPE. The discussion and result sections are thus divided into two different sec- tions. The results and discussions for the linear case are presented in Sec 8 and the results for the nonlinear case (the GPE) are presented in Sec 9. The discussions for the nonlinear case (the GPE) are presented in Sec 10 separately. Section 11 will cover the conclusions of the thesis. Before presenting all the work done and the literature study needed for the report, we will shortly describe the limitations and scope of the report and provide an outline and useful notations next.

1.2 Limitations and Scope

The limitations of this report are set to the time-independent equations of the model prob-

lems. The time dependent equations exist but are not in the scope of this thesis. The finite

element analysis (FEA) that will be done in the report can in general terms be described by

the following steps,

(17)

1. Pre-processing

a) Establishing the coordinates, and finite set of nodes and elements.

b) Discretizing the problem and identifying the primary unknowns.

2. Solver

Constructing the numerical scheme of the discretized problem and solve it 3. Post-processing

Illustrating the results and make sense of the data obtained.

4. Validation by analytical and experimental methods.

(18)

2 Physical Background

The discovery of quantum mechanics (QM) and its wave-particle duality is as mentioned before essential for explaining all particle behavior. When applying classical mechanics to explain particle phenomena the results clearly conflict with experiment. The common example provided in quantum mechanics literature is the two slits experiment where an electron beam passages through one slit and forms a continuous pattern on a screen placed behind the slit. A second passage of electron beam next pass through two slits, one above the other, that according to classical mechanics should lead to a superposition of their independent path patterns. Instead, a diffraction pattern caused by wave interference is observed in such experiments. Landau, L.D. (1965, page 3) As a result numerous applications in microscopic mechanical systems in modern world successfully use the theory of quantum mechanics and the wave-particle duality. In computer and smartphones for example we find both the quantum devices: the laser and packs of transistors at nano-meter scale that make the computer chip. Samuel J. Ling, Jeff Sanny, and Bill Moebs. (2019, Chapter 7.2) The phenomenon Anderson localization is thus explained with quantum mechanics as opposed to classical mechanics.

2.1 Mathematical notations

Symbol Expression or values Name

∆ or ∇

² ^∂_x²

1

+ .... +

_x^∂²

n

Laplacian

∇ (

_∂x^∂

1

, ...,

_∂x^∂

n

)

^>

Gradient

x Space variable in one dimension

p Momentum variable in one dimension

~ r [x

₁

, x

₂

, ..., x

_n

] Space vector (multiple dimension)

~

p [p

₁

, p

₂

, ..., p

_n

] Momentum vector (multiple dimensions)

δ~ r Infinitesimal distance in space (multiple dimensions)

i i

²

= −1 Imaginary number

h 6.626068 m

²

· kg/s Planck’s constant Ψ(~ r) Ψ(x

1

, x

2

, ..., x

n

) Wave function

Table 1: List of Mathematical Notations

(19)

2.2 The Uncertainty principle

The uncertainty principle is a key concept in QM developed by Heisenberg. In rough terms, the principle states the exact position and the exact momentum of any quantum object is impossible to measure simultaneously. The more we know about the latter, the less we know of the former, and vice versa. An illustration that can explain it is shown in Fig 3.

Figure 3: Several plane waves on the left hand side and a wave packet on the right hand side.

Acrylic painting by Monika Diaz, February 2020)

As quantum objects and particles behave as waves, they are expected to be found any- where in their problem space according to the uniform probability density that describes them. We are uncertain of the position but completely certain of the particles momentum because we know that a particle has a definite value of wavelength. By linearly superposing different plane waves, we form a localized particle confined to the region ∆x (illustrated in Fig 3). The smaller the confined area ∆x, the more combinations of momenta we have. This is summarised by the relationship suggested by Heisenberg’s principle,

∆x∆p ≥ ~

2 , (1)

where x here stands for displacement, p for the momentum, ~ =

_2π^h

is the reduced Planck’s constant, and ∆ the range of uncertainty for the given quantities. The minimum value of product is given by

∆x∆p = ~

2 . (2)

If we for example have an uncertainty in speed ∆u of an electron, we can compute the

uncertainty in momentum (∆p = m∆u where m is mass of the electron) and use equation

(2) to compute the uncertainty in space. Samuel J. Ling, Jeff Sanny, and Bill Moebs.

(20)

(2019, Chapter 7.3). A mathematical explanation for the uncertainty principle is that the exact momentum and the exact position are conjugates of each other, which implies that we cannot find them in exact terms simultaneously. This is justified by the wave mechanics interpretation,

Ψ(~ r) = 1

√ 2π~

Z

∞

−∞

φ(~ p)e

^i~^p~^~^r

, (3)

where i is the imaginary number and ~ r is the displacement vector. Here φ(~ p) represents the mode amplitude and is also called the wave function in momentum space. Ψ(~ r) represents the wave function of the quantum particle described by the integral over all possible modes and ~ =

_2π^h

the Planck’s constant. Landau L.D. and Lifshitz E.M. (1965, page 3) As the φ(~ p) is the Fourier transform of Ψ(~ r) we have that ~ p and ~ r are conjugates variables and impossible to determine simultaneously.

The Wave function

The statistical interpretation of the wave function Ψ(~ r, t) and the interpretation of its square

|Ψ|

²

, called Born interpretation, will be presented here. The probability density function of a quantum particle describes the probability of finding the particle in a narrow interval [~ r, ~ r + δ~ r] at time t. The probability is uniform and also time-independent.

P ([~ r, ~ r + ~ δr], t) = |Ψ(~ r)|

²

δ~ r. (4)

If the wave function varies rapidly, the square of the wave function has to be integrated over the narrow interval to consider the variations. Furthermore, if we are interested in finding the probability over the complete enclosed domain Ω, the expression also includes the integral over the square of the wave function,

P ([~ r, ~ r + δ~ r], t) = Z

Ω

|Ψ(~r)|

²

d~ r. (5)

Samuel J. Ling, Jeff Sanny, and Bill Moebs. (2019, Chapter 7.1, Eqn 7.1.1)

The variations of the wave function across the enclosed domain Ω have to be considered, which will be the case for the Bose gases in this thesis. Moreover, the total probability of finding a quantum particle within the enclosed domain Ω at an arbitrary time equals 1, meaning it is impossible to find the particle outside of Ω. This constraint is described by the expression,

Z

Ω

|Ψ(~r)|

²

d~ r = 1, (6)

(21)

and will thus be added as a conditional constraint to the model problem of the thesis.

A probability distribution for a wave is shown in Fig 4 (b) and examples of wave functions and their squares are given in Fig 5.

Figure 4: The probability density distribution |Ψ

_n

(~ r)|

²

for a quantum particle at state for: a) ground state n=1; b) the first excited state, n=2; and , c) the nineteenth state, n=20. Acrylic painting by Monika Diaz, February 2020.

Figure 5: Wave functions and their squares. Acrylic painting by Monika Diaz, February 2020.

(22)

Copenhagen interpretation

An interesting interpretation is the Copenhagen interpretation, which aims to explain where a particle is located. Consider a box with an enclosed problem space Ω and a quantum particle in it. The first explanation to the question "where is the particle located?" is that when the observer is not looking, the particle exist everywhere in the enclosed space. When the observer is looking, the second explanation says that the particle jumps to a position state (~ r, ~ r +δ~ r) with a probability |Ψ(~ r, t)|

²

in a process called wave-function collapse or state reduction. The explanation of the Copenhagen interpretation thus depends on whether or not the observer is looking. The assumption that a particle can only be at one location and only have one momentum and energy at a time (according to classical mechanics) is thus abandoned in quantum physics. It is assumed that the quantum particle has an "indefinite position" meaning that it, for example, can be in one location when the observer is looking and in another one when the observer is not looking. Landau L.D. and Lifshitz E.M. (1965, page 3)

2.3 Identical Particles

Distinguishable particles are able to be labeled and followed forever. Identical particles however are not able to be labeled because they are indistinguishable in behavior and ap- pearance. Identical particles have equal mass, charge, and spin. If you consider four dis- tinguishable particles and a set of four identical particles for instance. We can observe that the first set of particles can be ordered in 24 unique ways (4!). For the second set it is not possible to tell. It is said that the second set has a multiplicity of 1 and the first set has a multiplicity of 24. Beatty, Adam. (2013) The concept of identical particles is an important concept for statistical mechanics as it is based on probabilistic arguments.

Additionally, there are two different types of identical particles that exist according to quan- tum field theory. They are called fermions and bosons. Fermions only have half-integer spin and obey Fermi-Dirac statistics while bosons have integer spin and obey Bose-Einstein statistics. In the report, the second type of identical particles are described by the GPE (Bose gases). Since bosons have an integer spin they have the capability of sharing quantum states, which fermions do not (Pauli exclusion principle).

2.4 Bose-Einstein Condensate

With this in mind, the state of matter Bose Einstein Condensate (BEC) is introduced next.

BEC is simply put a group of identical particles that behave as if it were one big super

(23)

particle resulting in an extreme wave-like behavior. Whitfield, John (2013) The individual particles are hard to distinguish in this state of matter as their wave-packets overlap making them oscillate in unison. This state of matter is unlike solid, gas and liquid difficult to imagine. However a good comparison is to remark what happens as the temperature of atoms increases and decreases. As we increase temperature, the atoms form gases. As we decrease temperature the atoms undergo a transformation and become liquids. Lowering the temperature of liquids, the atoms undergo another transformation to become solids.

As we decrease the temperature to extremely low temperatures next, the vibrations of the atoms are slowed down extremely and reach lower quantum states while the wave-packets (which characterize the atoms) also increase in size. As bosons, quantum particles that are capable of sharing states, reach temperatures close to absolute zero, (-273.14 C

^◦

= 0K), they eventually share the lowest energy state, called the ground state. At the point of reaching the ground state the bosons wave-packets have become so big that they overlap and are no longer distinguishable. The wave-packets (the particles themselves) have all collapsed into the lowest possible energy state (the ground state) and hence all occupy the same quantum state simultaneously. In this state, the atoms have transformed into a peculiar state of matter, the BEC, where they form a giant wave vibrating in unison. It is in this state of matter that the Anderson localization may arise, which we will further investigate in the report.

Landscape analogy

A helpful analogy of the Anderson phenomenon in BEC is comparing the particle localiza- tions with a landscape full of pointy mountains hidden under sea. Say we are figuring out the landscape under a sea by tuning the sea level, or observing how the sea-level is altered with time as the average temperature alters. As we alter the sea-level up and down, eventually we understand the shape of underlying landscape. Similarly, the localizations of the dilute Bose gases at ultra low temperatures are revealed as the energy level of the system is altered.

As the Bose gases reach ground state energy level, the strength of the particle interactions vanishes and the system behaves like one big super-atom, as explained in the concept of BEC (Sec 2.4). To understand the behavior of the dilute Bose gases as we alter the energy level, it is thus needed to analyze the numerical approximation of the mathematical model GPE for different parameters.

Quantum tunneling

Another analogy worth considering is that of physical barriers and the level of the potential.

If we make this analogy and think in classical terms, how is it then possible that the Bose

gases localize at extremely localized spots instead of diffusing given the random nature of

the potential? The explanation for this is given by the phenomenon quantum tunneling. It

describes particles penetrating potential energy barriers with heights greater than the energy

(24)

of the particles themselves although they should not ’physically’ be able to according to classical mechanics. This question is justified by the Born interpretation and the Copenhagen interpretation discussed in Sec 2.2. It is not impossible for the quantum particle to be located on ’the other side’ of a potential barrier in a sufficiently random potential, it is only less probable that the particle will tunnel/ relocate to the other side. The taller and wider the barriers are, the less probable the relocation is. An illustration of quantum tunneling is found in Fig 6.

Figure 6: Quantum tunneling representation. Acrylic painting by Monika Diaz, February 2020.

Although we can understand the physical concepts behind the Anderson localization, it is still difficult to predict were the localizations will be located exactly. Understanding the P.W.

Anderson localization requires more than looking into its simple components. The approach

for solving complex systems is suggested by Anderson in his paper "More is Different" (1972)

where he argues that the issue with solving complex systems lies in that of the reductionist

hypothesis: the belief that all systems can be reduced to simple fundamental laws. The

reductionist hypothesis does not ’by any means imply that of the constructionist one: the

ability to start from those (fundamental equations) and construct the universe’ Anderson

says. Numerical investigations of the Anderson localization and other complex systems are

argued to be needed in order to reveal the true behavior that otherwise is impossible to

predict theoretically or by ’induction’.

(25)

3 Finite Element Background

The Finite element method (FEM) is a numerical method used to solve continuous differ- ential equations with boundary conditions by approximating the variational formulation of the differential equation with a linear system of algebraic equations based on nodal values.

Imagine that we are working to solve a problem defined over a two-dimensional problem space Ω that we partition using a finite set of nodes across the entire problem space. The number of nodes then equals the number of degrees of freedom that we solve for using the finite element equation obtained using the method, i.e. the number of algebraic equations.

The nodal values that are solved for can in structural problems for example represent the displacements in space. They can for example also represent the alternation in temperature across a plate. In this report the nodal values will represent point values of the wave func- tion that describes the quantum state of the ultracold Bose gase that is trapped in a random potential.

3.1 The development and advantage of FEM

The story of the development of the numerical method Finite element method (FEM) can be traced back to the early 1940s in the work of the structural engineer A.Hrennikoff and the mathematician R.Courant. A.Hrennikoff (1941, page 169) R.Courant (1943, page 1) It is a success story on its own as it had physicists and mathematicians across the world collaborating together through open source in order to solve structural problems of aerospace engineering during the computer age. Strang, Gilbert. (2013, page 1132) Long before this time Galerkin contributed with the Galerkin method which was used obtain high accuracy with minimal computational effort using cleverly chosen basis functions, as was first presented in his paper on elasticity published in 1915. Eriksson, K. (1996) It is thus difficult to attribute the discovery of FEM to one single person as many different types and aspects of FEM were developed during many centuries. Characteristic for all FE methods however is that they are based on variational principles rather than differential equations. Strang, Gilbert. (2013, page 1132)

What is so powerful about the Galerkin Finite Element Method is that it produces the same algebraic structure for a large set of physical problems. Approximating the solutions and the equation using the same interpolation functions produces a simple and easily solved algebraic system of equations, where the coefficients of the basis functions are the unknowns.

Since the discretization process of the FEM is based on the variational formulation of a dif-

ferential equation, the process also weakens the differentiability- and satisfaction of boundary

requirements of the original physical problem. Using FEM as a method for solving a contin-

(26)

uous differential equation thus means finding an approximation of the continuous solution in the discrete space rather than in the continuous space.

3.2 Procedure

The FEM will be used to solve the model problem of the report that will be presented in Sec 4, and before we present the model problem we will both present the numerical setting and go through some technical background i.e. ’Finite Element background’ needed.

To begin with, an example of a continuous differential equation with boundary conditions (in 1D) is (D):

(D)







−

αu

⁰

(x)

⁰

= f (x) x ∈ (0, 1) u(0) = u(1) = 0.

To solve the differential equation the first step is to produce the variational formulation which will be the basis of the Finite element method. Using Galerkin method, which means multiplying the problem (D) with an arbitrary test function v and integrating over a problem domain, we obtain the variational formulation or weak problem (W),

(W)



 

  R

1

0

αu

⁰

v

⁰

dx = R

1

0

f vdx ∀v ∈ V, u(0) = u(1) = 0,

α(x) ≥ α

₀

> 0 ∀x

were the solution of (W) belongs to the test space V which coincides with the Sobolev space H

₀¹

, defined by

V = H

₀¹

(0, 1) = {v ∈ C

⁰

(0, 1) : kvk

_L²_(0,1)

+ kv

⁰

k

_L²_(0,1)

< ∞ v(0) = v(1) = 0}. (7) We note that since the weak formulation is a variational integral formulation of problem (D), the differentiability requirements of the original problem are reduced. The solution to the weak formulation (W) is thus a generalized solution typically called the weak solution, whereas the solution to the problem (D) – if it exits – is typically called a strong solution.

We refer to the weak solution as the exact solution. In the expressions for (D) u represents

the strong solution while in (W), u represents the weak solution in V. The test function v

belongs to the same test space V and is arbitrary. Another remark is that a strong solution

is always a weak solution, but a weak solution is not always a strong solution. A weak

solution is only a strong solution if the it fulfills additional requirements. We say that both

concepts are equivalent if the setting allows for enough regularity. In addition it can be

shown with the Lax-Milgram lemma that for every constant α 6= 0, the weak solution always

(27)

exists and is unique. Strong solutions do not always exist, because the requirements (two- times differentiable) might be too strong in some cases. The next step of the Galerkin finite element discretization is assuming that the solution can be approximated with a combina- tion of a finite set of basis functions that can be global polynomials, piece-wise polynomials, trigonometric polynomials or other functions. In this report the FEM equation is based on piece-wise polynomial basis functions (presented in Sec 3.4). Before further discussing how we can develop the Galerkin FEM formulation for (D) we will introduce some additional FEM concepts.

3.3 Function spaces

In this section we will present the function spaces that are needed for the Galerkin FEM.

Expressions for the L

²

norm and the H

¹

-norm for a given problem domain Ω are,

||v||

_L²_(Ω)

= Z

Ω

|v|

²

dx

¹₂

,

||v||

_H¹_(Ω)

= Z

Ω

|v|

²

+ |∇v|

²

dx

¹₂

.

These will be used in the report to define spaces of the solutions of the weak problem and the FEM solution. The H

²

-norm equals,

kvk

_H²_(Ω)

= X

²

i,j=0,i+j≤2

||∂

_x_i_x_j

v||

²_L2(Ω)

¹₂

. (8)

In one-dimension this reduces to

||u||

_H²_(Ω)

= q

||u||

²_L2(Ω)

+ ||u

⁰

||

²_L2(Ω)

+ ||u

⁰⁰

||

²_L2(Ω)

. (9)

3.4 Basis functions

In the same way that every vector in a vector space can be described as a linear combination

of the basis vectors of that space, every continuous function of a function space can be

described as a linear combination of the basis functions of that space. The basis functions

are linearly independent and can be used in interpolation, and are used in Galerkin FEM to

describe the FEM solution in a structured way by making it possible to simplify the FEM

(28)

equation. Common basis functions used for solving many structural problems are called the Lagrange basis functions. The Lagrange basis functions of order 1 (also called hat-functions), are given by

φ

_j

(x) =



 

 

0 x

_j+1

< x, x ≤ x

_j−1

(x−xj−1)

xj−x_j−1

x

_j−1

≤ x ≤ x

_j

1 −

_x^(x−x^j⁾

j+1−x_j

x

_j

≤ x ≤ x

_j+1

,

for nodal points j = 1, ..., n. These simple basis functions of order one are commonly used for solving physical problems because of the low regularity of such problems. Many physical problems do not need to have many derivatives to be fully defined and thus using polynomials of higher polynomial degree as basis functions would not improve the accuracy of problems with low regularity. In this report, the basis functions will thus be of order one. It is interesting to mention however that for problems with high regularity, increasing the order of the polynomial would also increase the accuracy of the approximation.

Characteristic for all basis functions, despite their polynomial order is that they take the value one at their nodal point j (x

_j

), and value zero at its neighboring nodal points (and outside of this neighborhood it is extended by zero). We also say, in both the onedimensional case and twodimensional case, that two basis functions are almost orthogonal to each other over a problem domain Ω if they are far enough apart. This is because only the basis functions φ

_i

and φ

_j

sharing support on a common triangle (in 2D) or string element (in 1D) give non-zero value in expressions such as

Z

Ω

φ

_i

(x)φ

_j

(x) dx or Z

Ω

αφ

⁰_i

(x)φ

⁰_j

(x) dx.

In one dimension the basis function takes the shape of a triangle with value one in nodal

point j and value zero at its neighbouring points and beyond. In two dimensions the shape

is slightly more complicated. An example of a piecewise linear Lagrange basis function in

two dimensions is provided in Fig 7.

(29)

Figure 7: The support of a piecewise linear basis function φ

_j

in 2D. The function takes the value one in nodal point j and value zero at the neighbouring nodal points. Drawing by Monika Diaz, February 2020.

3.5 Trial function and trial space

To proceed with developing the FEM equation for (D) we next introduce a discrete trial space which is a finite dimensional subspace of V. We refer to the set of functions V, in which we seek for the solution u, as the continuous solution space and the set of functions V

_h

, where we seek the discrete solution u

_h

, as the (discrete) trial space. The trial functions will be the basis functions of the discrete space V

_h

(on the interval Ω = [0, 1]), defined by

V

_h

= {v ∈ C

⁰

(Ω) : v|[κ] ∈ P

¹

(κ) for all κ ∈ T

_h

, v ≡ 0 on ∂Ω, } (10) where T

_h

= {κ} is a partition of the problem domain Ω. To obtain the Galerkin FEM formulation we assume that the solution approximation is a combination of piecewise linear basis functions,

u

_h

= X

j

ζ

_j

φ

_j

(x). (11)

The next step is substituting the expression v with a discrete basis function (first order

Lagrange function) φ

_j

(x). This can be done because the test function v can be chosen arbi-

trarily and because V

_h

is a subspace of V. Adjusting the variational formulation with the new

(30)

requirements (piecewise linearity and substitute the test function with the basis function) we finally obtain the algebraic Galerkin FEM formulation:

Find u

_h

∈ V

_h

(i.e. find ζ

₀

, ..., ζ

_n

∈ R) such that

(W

^GF⁰

)



 

  R

1

0

αu

⁰_h

φ

⁰_i

dx = R

1 0

f φ

_i

dx u

_h

(0) = u

_h

(1) = 0, α(x) > 0 ∀x or respectively

(W

^GF

)



 

  P

n

j=1

ζ

_j

R

1

0

αφ

⁰_i

φ

⁰_j

dx = R

1 0

f φ

_i

dx u

_h

(0) = u

_h

(1) = 0,

α(x) > 0 ∀x

(12)

for all 1 ≤ i ≤ n. Since φ

_i

forms a basis of V

_h

, this implies that Z

1

0

αu

⁰_h

v

⁰_h

dx = Z

1

0

f v

_h

dx for all v

_h

∈ V

_h

.

Furthermore, since basis functions are almost orthogonal, the integral on the left hand side of equation (12) represents a matrix with entries

a

_ij

= X

κ∈Th

a

^κ_ij

= X

κ∈Th

Z

κ

αφ

⁰_i

φ

⁰_j

dx, (13)

with value zero for most of the entries. The integral on the right hand side represents a vector with entries

b

i

= X

κ∈Th

Z

κ

b

^κ_i

= X

κ∈Th

Z

κ

f φ

i

(x)dx. (14)

3.6 Discretization

An essential first step in solving the Galerkin FEM equation is dividing the problem space

Ω into a finite number of elements that we call a triangulation T

_h

= {κ}. A triangulation

is in other words a subdivision of the problem domain, where none of the elements overlap

each other. The domain of the Galerkin FEM problem W

^GF

presented in the previous

section is for example an interval from zero to one. Its triangulation, i.e. subdivisions are

(31)

subintervalls that are connected to each other by a node. In a two-dimensional case, a uniform triangulation is illustrated in Fig 8. The corners of these triangles are the nodes x

_j

and there is a basis function φ

j

for each node.

Figure 8: 2D Example of a triangulation with equidistant intervals in both axes produced in MATLAB)

3.7 The unknowns of the FEM equation

Our finite element solution is of the form

u

_h

= X

j

ζ

_j

φ

_j

(x), (15)

where the coefficients in front of the basis functions, i.e. ζ

₀

, ..., ζ

_n

, will be the unknowns of

an algebraic system of equations produced by the Galerkin FEM, see (W

^GF

).

(32)

3.8 Discretization notations

In the following, we fix the notation.

Symbol Name

n + 1 Number of degrees of freedom

i Discretization index / index of a node x

_i

or basis function φ

_i

T

_h

The non-overlapping partition of Ω into triangles (or subintervals).

κ An arbitrary triangle of T

h

.

h Maximum diameter of a triangle in T

_h

N

_n

Node at degree of freedom ’n’

E Energy level of a particle

L

²

(Ω) The space of square-integrable functions on Ω C

¹

(Ω) Continuously differentiable functions

C

^k

(Ω) The function and its derivative up to order k exist and are continuous

H

₀¹

(Ω) Sobolev space of weakly differentiable functions which are zero on the boundary

∂

_r

Partial derivative with respect to variable r

∇ Gradient operator

Ω Problem domain

∂Ω Boundary of the problem domain

Table 2: List of FEM and mathematical notations

3.9 One-dimensional example

Say we are interested in applying the Galerkin FEM to the following differential equation in one dimension that we name (D),

( −u

⁰⁰

(x) = f (x) Ω = [0, 1]

u(0) = u(1) = 0.

The weak formulation of the equation is obtained by multiplying the test function v ∈ V = H

₀¹

(Ω) = {v ∈ C

⁰

(Ω) : kvk

_L²_(Ω)

+ kv

⁰

k

_L²_(Ω)

< ∞}, and integrating over the problem space. The weak formulation is

Find u ∈ V such that, Z

1

0

u

⁰⁰

(x)v(x)dx = Z

1

0

f (x)v(x)dx ∀v ∈ V(0, 1),

(33)

where v(0)=v(1)=0. Integrating by parts and using the boundary condition yields, Z

1

0

u

⁰

(x)v

⁰

(x)dx = Z

1

0

f (x)v(x)dx ∀v ∈ V(0, 1).

The next step is introducing the finite dimensional subspace of V(Ω) = H

₀¹

(0, 1):

V

_h

= {v ∈ C

⁰

(Ω) : v|

_[x_i−1_,x_i_]

∈ P

¹

, v(0) = v(1) = 0, }, i = 1, 2, ..., n

where x

_i−1

, x

_i

are the corners of the element κ

_i

= [x

_i−1

, x

_i

], which is just an interval in 1D.

The variational formulation with piecewise linear functions is then Find u

_h

∈ V

_h

such that,

Z

1 0

u

⁰_h

(x)v

⁰

(x)dx = Z

1

0

f (x)v(x)dx ∀v ∈ V

_h

(0, 1). (16) The final steps for producing the algebraic system of equations is inserting the two sub- stitutions

• u

_h

= P

n

j=0

ζ

_j

φ

_j

(x) and

• v = φ

_i

(x)

into the discrete variational formulation (16). The Lagrange polynomial of order one will be the basis function in this problem and reads

φ

_i

(x) =



 

 

0 x

_i+1

< x, x ≤ x

_i−1

(x−xi−1)

xi−xi−1

x

_i−1

≤ x ≤ x

_i

1 −

_x^(x−xⁱ⁾

i+1−x_i

x

i

≤ x ≤ x

i+1

.

These two substitutions are the final steps needed to form the Galerkin FEM formulation of the original problem which we denote (W

^GF

). The linear system of equations is,





 R

1

0

φ

⁰₀

(x)φ

⁰₀

(x)dx · · · R

1

0

φ

⁰_n

(x)φ

⁰₀

(x)dx

.. . .. .

R

1

0

φ

⁰₀

(x)φ

⁰_n

(x)dx · · · R

1

0

φ

⁰_n

(x)φ

⁰_n

(x)dx











 ζ

₀

.. . ζ

_n





 =





 R

1

0

f (x)φ

₀

(x)dx .. .

R

1

0

f (x)φ

n

(x)dx





 .

Since each of the basis functions φ

_j

equal one at the nodal point j, the coefficient we solve

for in the Galerkin FEM equation approximate the exact solution at that node. This is an

important remark which can be compared with the analytical (exact solution, if it exists) of

the problem at those nodes to obtain an error.

(34)

The coefficients (and exact nodal values) are represented by the vector ~ ζ. The left hand side integral represents a matrix A with entries

a

_ij

= X

κ∈Th

Z

κ

φ

⁰_i

φ

⁰_j

dx (17)

and the right hand integral represents the vector ~b with entries b

_i

= X

κ∈T_h

Z

κ

f φ

_i

dx. (18)

where both i, j = 0, ..., n. The manipulations done on the original differential equation up to this point yields an algebraic system of equations,

A~ ζ = ~b. (19)

The only unknowns are the coefficients of the FEM expression (if the right hand side is given).

The load function can for example represent the exterior loads applied to a mechanical system and the coefficients of the basis functions can represent the nodal displacements of a rigid body.

3.10 h- and q- method

The h- and q- methods are two methods used to increase the accuracy of a numerical ap- proximation of a Galerkin FEM solution. The h- method is the classical way of increasing the accuracy and means that we refine the triangulation, which also means increasing the number of nodes and thus the computational time. The other method is called the q- method and means increasing the degree of the polynomial of the basis function to increase accuracy.

It suggests that the error in the H

¹

-norm is bounded by the product of a constant value C

_q

, the exponent h

^q

and the H

^q+1

- norm of the solution. It holds an estimate such as

kek

_H¹

0(Ω)

= k∇u − ∇u

h

k

_L²_(Ω)

≤ C

q

h

^q

kuk

_H^q+1_(Ω)

,

where u represents the exact solution and u

_h

the approximate solution produced by FEM.

The expression suggests that the higher the degree q of the used polynomial, the higher is the expected convergence rate in h, where we obtain the rate h

^q

. This however also requires that the solution is (q+1)-times differentiable. Since we deal with a physical problem with low regularity in this report, we do not observe a second order convergence rate in the H

¹

-norm as the q- method claims (for q = 2). It will be sufficient to consider finite elements with a polynomial order one, q = 1. The estimate for the error then becomes

k∇u − ∇u

_h

k

_L²_(Ω)

≤ Chkuk

_H²_(Ω)

,

(35)

meaning that we can expect a linear order convergence in h. The conclusions is thus that any higher order basis functions only add unnecessary computations and do not allow us to see the second order convergence the method claims we would observe. As a result, we will use piecewise linear basis functions in the report.

An illustration of the degrees of freedom in a triangular Lagrange elements of different orders is provided in Fig 9.

Figure 9: P

^q

Lagrange elements on a triangle. Drawing by Monika Diaz, February 2020

(36)

4 The model problem

A good approach to solve any complex problem is to look into easier similar problems. In this report, the Gross Pitaevskii equation (GPE) is an eigenvalue problem where the linear part of the equation is the linear and time independent Schrödinger equation, which is a type of elliptic eigenvalue problem. The approach for solving the GPE will thus be to first solve the linear Schrödinger equation. The approach to evaluate the numerical scheme of the Schrödinger equation before working on the numerical scheme of the GPE will be to compute the difference between the approximate eigenvalue and the exact eigenvalue for a fully defined example problem. The accuracy is evaluated by a tolerance criteria. We will also look into the shape of the solution produced in order to verify that the numerical scheme is doing the right thing. The Schrödinger equation is a type of diffusion equation. (Boris Baeumer, Mark M. Meerschaert and Mark Naber. 2010, page 1) We should therefore expect that the solution of a quantum system that is only influenced by constant potential field has a symmetrical wave-like shape.

4.1 Schrödinger Equation

As mentioned in Sec 2, the square of the wave function can be interpreted as the probability distribution of a particle’s position over problem domain Ω. Because of this wave-particle duality, the words ’particle’ and ’wave’ are used interchangeably in this section. The principle of uncertainty (Sec 2.2) can explain why particle mechanics does not confine with the theories of classical mechanics and to investigate mechanical systems at particle level we are forced to look for a model beyond the classical mechanical equations that are insufficient for these types of problems. The theories of interest is of course described by or based on Quantum Mechanics (QM). Even though we cannot compute the exact position and exact momentum of a particle simultaneously, we can compute the exact energy and its corresponding wave function of the quantum object. Consider a particle in a potential field. Schrödinger was a physicist who helped develop the theory behind the mathematical equation which states that the energy and wave function of particles (in a potential V

⁰

) behave according to the relationship,

i∂

_t

u(~ r, t) =

− ~

²

2m ∆ + V

⁰

(~ r)

u(~ r, t). (20)

This is called the Schrödinger equation. Because this PDE generates a solution u(~ r, t) that

obeys Huygens principle of diffraction, i.e. has an oscillating and smooth shape, the solution

describes a standing wave that propagates in time t and space ~ r. If we assume that the

(37)

wave function is a product of an exponential function in time and a time independent wave function in space,

u(~ r, t) = e

^−iEt

Ψ(~ r), (21)

and substitute this assumption into the Schrödinger equation, we obtain the ’time-independent’

form of the Schrödinger equation EΨ(~ r) =

− ~

²

2m ∆ + V

⁰

(~ r)

· Ψ(~r), (22)

where the notation E represents the quantum energy, V

⁰

(~ r) represents a time independent potential, ~ =

_2π^h

a scaled value of the Planck’s constant. The notation m represents the particle mass and ~ r the position vector. The time independent wave function denoted Ψ represents the shape of the standing wave at quantum energy level E and space ~ r and does not depend on time t. Moreover, the equations yields predictions about the future state of a quantum particle given its current state, similar to Newton’s second law of motion for mass objects. We can for instance calculate the momentum of a quantum particle using only the energy and wave function that we obtain with the Schrödinger equation just as we can calculate the momentum of a mass object using Newton’s second law of motion. The terms within the parenthesis in Eq (22) have a central meaning to the equation and is called the Hamiltonian,

H = − ˆ ~

²

2m ∇

²

+ V

⁰

(~ r). (23)

It represents the total energy enclosed in the conservative quantum system (no dissipation of energy). Furthermore, the Schrödinger is intriguing not only by all the information we can obtain about a quantum system but also by how the equation is capable of simplify a complex nature into a linear expression. This is not trivial at all. There are limitations however. The Schrödinger equation was applied to a hydrogen atom when it was developed and it is thus an excellent model for making accurate predictions about quantum particle(s) that are not affected by their interaction to each other. As soon as there is interaction however, the approximations are not correct. In this report we are interested in considering the influence of the weakly interacting Bose gases on the Anderson localization phenomenon.

This model will not be sufficient to model the quantum system we will investigate but gives us a basic idea of the type of how quantum systems are described.

Before presenting a mathematical equation that considers the interaction and that will be used in this report, a last remark worth considering regarding both Eq (20) and Eq (22), is that Eq (22) represents the eigenvalue problem of (20) where the wave function Ψ(~ r) corresponds to the eigenfunction with corresponding eigenvalue λ

⁰

,

λ

⁰

= E. (24)

(38)

The eigenvalue coincides with the quantum energy. When the eigenfunction satisfies a pre- defined boundary conditions on top of the eigenvalue problem, it is also called system’s mode function.

4.2 GPE or the Total Energy Schrödinger Equation

The Gross-Pitaevskii equation (GPE) was applied to identical bosons fully condensed in ground state when developed (Bose-Einstein Condensate described in Sec 2.4), and considers the interactions between bosonic particles. This mathematical equation is sufficient and will be used as a model when reproducing the Anderson localization phenomenon numerically in this report. The Gross-Pitaevskii equation is,

µΨ(~ r) = − ~

²

2m ∆Ψ(~ r) + V

⁰⁰

(~ r)Ψ(~ r) + N 2π~

²

a

_s

m |Ψ(~r)|

²

Ψ(~ r), (25) where:

µ: eigenvalue of the GPE

~: reduced Planck’s constant [J · s], V

⁰⁰

(~ r) =

^m₂

V ( ~ ¯ s

_k

r): the potential [J]

a

_s

: scattering length [m] ∼ 10

⁻⁹

,

N : number of bosons (dimensionless) ∼ 10

³

− 10

⁷

. We obtain the simple expression,

µΨ(~ r) =

− ~

²

2m ∆ + V

⁰⁰

(~ r) + δV

⁰⁰

(~ r)

Ψ(~ r), (26)

by substituting the coefficient of the non-linear term with, δV

⁰⁰

(~ r) = 2N π~

²

a

_s

m |Ψ(~r)|

²

. (27)

A way of comparing the time independent Schrödinger equation (22) with the Gross-Pitaevskii equation (GPE) (25) is by using a conversation of energy perspective. If the Schrödinger equation considers the kinetic- and potential energy of a quantum system with many quan- tum particles, then the GPE considers not only the kinetic- and potential energy, but also by the interaction energy between bosonic particles (represented by the nonlinear term). Note that if the interaction energy is insignificant the GPE is equal to the Schrödinger equation.

Therefore, the GPE is also called the ’total energy Schrödinger equation’. The mathemat-

ical appearance is to that extend very similar, but the solutions produced by the different

models are the contrary. They have little in common. While the Schrödinger equation is

(39)

fundamental for describing simple quantum mechanical systems, the GPE describes a differ- ent regime of matter and complex quantum systems, that historically was considered purely theoretical. As mentioned in Sec 2.4 weakly interacting Bose gases at ultra-low temperatures transform into the state of matter Bose-Einstein Condensate (BEC) due to the integer spin property. This state of matter describes Bose atoms forming a giant wave/particle with peculiar physical properties. As pointed out in the physical background section, this is be- cause the atoms share the same ground state energy level in BEC, where there is no energy left to move relative to any other particles. In other words, this means that the solution of the GPE describes one wave composed by different interacting quantum particles (that are interacting to the extent that they act as if they were one giant particle). To observe the Anderson localization phenomenon, the Bose gases need to be in the BEC state of matter.

The disorder and strength of the potential field also need to be sufficiently strong.

How does the interaction energy between the Bose gases affect the Anderson localization phenomena? By comparing the GPE (25) with the Schrödinger equation (22), one thing we can directly point out is that the non-linear interaction term of the dilute Bose gases intervenes with the potential field, with a magnitude δV

⁰⁰

(~ r). This is clearly observed by rewriting Eqn (26) such that,

µΨ(~ r) =

− ~

²

2m ∇

²

Ψ(~ r) +

V

⁰⁰

(~ r) + δV

⁰⁰

(~ r) Ψ(~ r).

The notations related to the expression are:

~: reduced Planck’s constant [J · s], V

⁰⁰

(~ r) =

^m₂

V ¯

⁰⁰

( ~ s

_k

r): the potential [J]

V ¯

⁰⁰

(~ r): the trapping potential.

ω

k

: trapping frequencies [rad s

⁻¹

] k = x, y, z N : number of bosons (dimensionless) ∼ 10

³

− 10

⁷

; m: mass of the particles [kg] ∼ 10

⁻²⁵

,

a

s

: scattering length [m] ∼ 10

⁻⁹

, Ψ(~ r): wave function [m

^d²

]

∆: Laplace operator (∇

²

).

A harmonic trapping potential can for example look like: ω

_x²

x

²

+ ω

_y²

y

²

+ ω

_z²

z

²

Anderson localization in weakly interacting Bose gases

SECOND CYCLE, 30 CREDITS STOCKHOLM SWEDEN 2020 ,

A numerical investigation of