## 國立交通大學

### 應

### 用數學系

### 碩 士 論 文

### 橢圓與光學晶格位能井之非線性薛丁格方程

### 的數值研究

### Numerical Studies in Nonlinear Schr¨

### odinger Equations

### with Elliptic and Optical Lattices of Trap Potentials

### 研 究 生

### ：莊佩錚

### 指導教授：張書銘 博士

### 橢圓與光學晶格位能井之非線性薛丁格方程

### 的數值研究

### Numerical Studies in Nonlinear Schr¨

### odinger Equations

### with Elliptic and Optical Lattices of Trap Potentials

### 研 究 生

_{：莊佩錚}

### 指導教授：張書銘 博士

### Student: Pei-Cheng Chuang

### Advisor: Dr. Shu-Ming Chang

### 國立交通大學

### 應

### 用數學系

### 碩士論文

### A thesis

### Submitted to Department of Applied Mathematics

### college of Science

### National Chiao Tung University

### in Partial Fulfillment of the Requirements

### for the Degree of

### Master

### in

### Applied Mathematices

### July 2009

### Hsinchu, Taiwan, Republic of China

### 橢圓與光學晶格位能井之非線性薛丁格方程

### 的數值研究

### 學生：莊佩錚

### 指導教授：張書銘 博士

### 國立交通大學應用數學系（研究所）碩士班

### 摘

### 要

本論文主要運用數值方法來探討非線性薛丁格方程的解之穩定性分析，其中這邊考慮 的非線性薛丁格方程式是含有位能井（trap potential）。 首先，我們先介紹非線性薛 丁格方程式以及兩種不同的位能井。接著，簡述我們所運用的數值方法。最後，我們 呈現一些數值模擬結果。 經由我們的數值模擬，發現在一些解的形式與位能井的選擇 是可_{以觀察到解的穩定性變化。}

### Numerical Studies in Nonlinear Schr¨

### odinger

### Equations with Elliptic and Optical Lattices of Trap

### Potentials

### Student: Pei-Cheng Chuang

### Advisors: Dr. Shu-Ming Chang

Department (Institute) of Applied Mathematics National Chiao Tung University

Abstract

In this thesis, we would like to find out the stability of the solution which linearized around the special solution of the nonlinear Schr¨odinger equation. Here, we consider nonlinear Schr¨odinger equations with the potential part, optical lattices form and elliptic form. In the first part of the thesis, we introduce roughly the nonlinear Schr¨odinger equation and two kinds of trap potentials. Then, we show the mathematical analysis of the numerical method that we used. At the end, we give some numerical results of our numerical experience.

### 誌

### 謝

### 本篇論文的完成，最要感謝的是我的指導老師張書銘博士，於這兩年

### 來的用心指導，讓我得以對數值計算有更進一步的了解。除了研究上的指

### 導

_{，亦在未來人生的規劃與為人處事上給予我許多關懷指引。老師，謝謝}

### 您！於口試期間，亦承蒙黃聰明教授與吳金典博士給予我的建議及對疏漏

### 處

### 的導正，使本論文得以更加嚴謹完整，在此學生致上由衷的感謝。

### 另遠在英國的芳嫻姐姐及我的室友們：郁凌、廷芳、維彧，謝謝你們在

### 我完成論文的階段中，給予我的幫助、鼓勵與歡笑，有了你們的陪伴，讓

### 我的研究生活更加順利。

### 最後，謝謝家人對我的支持與關心，有了您們的鼓勵與幫助，讓我得

### 以專心完成本篇論文。僅將此小小成果與所有關心、幫助過我的長輩、老

### 師

### 、同學們一同分享這份喜悅，謝謝大家！

### Contents

中 中 中文文文摘摘摘要要要 i 英 英 英文文文摘摘摘要要要 ii 誌 誌 誌謝謝謝 iii Contents iv List of Tables v List of Figures vi 1 Introduction 11.1 Nonlinear Schr¨odinger Equations . . . 1

1.2 Trap Potential V . . . 2 1.2.1 Optical Lattices . . . 2 1.2.2 Elliptic Potential . . . 2 1.3 Goals . . . 3 2 Mathematical Analysis 6 2.1 Perturbation in NLS . . . 6 2.2 Stability . . . 8 3 Numerical Method 9 3.1 Matrix Form of Laplace operator . . . 10

3.1.1 Square Domain . . . 10

3.1.2 Unit disk Domain . . . 11

3.2 Algorithm . . . 12 4 Numerical Simulation 13 4.1 Setting . . . 14 4.2 ω case . . . 19 4.3 µ case . . . 25 5 Conclusions 30 Reference 31

### List of Tables

### List of Figures

1 Optical Lattices. . . 2

2 Elliptic potential and m(x) = 1 +1_{2}sin(2πx1). . . 4

3 Optical lattices potential and m(x) = 1 + 1_{2}cos(2πx1). . . 5

4 Boundary φ in ω case for λ = 0. . . 15

5 Boundary φ in ω case for λ = 1. . . 16

6 Boundary φ in µ case for λ = 0. . . 17

7 Boundary φ in µ case for λ = 1. . . 18

8 Solution of φ for ω1 = 1, ω2 = 1, = 0.001, λ = 1. . . 19

9 Solution of φ for ω1 = 1, ω2 = 1, = 0.001, λ = 0. . . 20

10 φ moves for λ = 1, = 0.0002. . . 21

11 Spectrum of ω case for λ = 1. . . 22

12 ∗ of ω case for λ = 1. . . 23

13 Value of λ in ω case. . . 24

14 Spectrum of µ case for λ = 0. . . 25

15 ∗ of µ case for λ = 0. . . 26

16 Spectrum of µ case for λ = 1. . . 27

17 ∗ of µ case for λ = 1. . . 28

### 1

### Introduction

In this thesis, we consider the nonlinear Schr¨odinger equation (NLS) with focusing pa-rameter ,

i∂tψ = −∆ψ + V (x)ψ − m(x)|ψ|2ψ,

where ψ(t, x) : R × R2 → C is the wavefunction, V (x) is the trap potential, m(x) is the mass and s-wave scattering length term, and > 0. We would like to consider two kinds of trap potentials and domains in this problem.

Firstly, we introduce nonlinear Schr¨odinger equations and the two kinds of trap po-tentials first.

### 1.1

### Nonlinear Schr¨

### odinger Equations

The Schr¨odinger equation is a basic equation in quantum theory. The study of Schr¨odinger
equations plays an important role in modern physics. In 1926, an Austrian physicist Erwin
Schr¨odinger constructed the Schr¨odinger equation for explaining the active of particles.
The famous Schr¨_{odinger equation is the form: Hψ = i~}∂ψ_{∂t}, where H is a Hamiltonian
operator, ψ is the wavefunction, ~ represents Planck’s constant over 2π, i is the
imagi-nary unit. The Schr¨odinger equation is an equation to describe the possible distribution
of atoms. Solutions to Schr¨odinger’s equation describe not only atomic and subatomic
systems, electrons and atoms, but also macroscopic systems, possibly even the whole
universe.

The nonlinear Schr¨odinger equation (NLS) is a nonlinear version of Schr¨odinger’s equation in theoretical physics. In mathematical point of view, the NLS equation is a Schr¨odinger equation with nonlinear term. The nonlinear Schr¨odinger equation is the partial differential equation

i∂tψ = − 1 2∂ 2 xψ + c|ψ| 2 ψ

for the complex field ψ. There are many kinds of nonlinear Schr¨odinger equations. One general form of such equations [12] would be

i∂tψ + ∆ψ = f (ψ, ¯ψ).

These equations (particularly the cubic NLS equation) arise as model equations from several areas of physics: nonlinear optics, quantum condensates.

Here, we review some properties of NLS equations. Let’s consider the equation i∂tψ + ∆ψ = ±|ψ|2ψ.

The sign “−” on the right hand side is focusing nonlinearity, and the sign “+” is defo-cusing. When we add the potential part to the NLS equation, then it becomes

i∂tψ + ∆ψ = −m|ψ|p−1ψ + V ψ, (1)

where V is real and time independent. The equation (1) is sometimes referred to as a Gross-Pitaevskii equation. We can look at the dynamics of the Bose-Einstein condensate from the time-dependent Gross-Pitaevskii equation.

The behavior of ψ is determined by the potential V of the NLS equation. There is a
special case of potential V = ±|x|2_{, this can be used to model a confining magnetic trap}

### 1.2

### Trap Potential V

In the previous subsection, we have introduced some property of NLS equation. Now, we introduce the trap potential. Trap potential is to trap atoms on the minima potential. The following are two forms of trap potentials in this study, which are optical lattices and elliptic form.

1.2.1 Optical Lattices

An optical lattice is essentially an artificial crystal of light. A periodic intensity pattern that is formed by the interference of two or more laser beams. The shape of optical lattices looks like an egg carton in Fig. 1(figure reference: [6]). It is called an optical lattice, since the periodic arrangement of trapping sites resembles a crystalline lattice. In Fig. 1, atoms are cooled and congregate in the minimal potential. The well depth and the periodicity are two important parts to affect the potential shape.

Figure 1: Optical Lattices.

Besides trapping cold atoms, optical lattices are also use in sorting microscopic parti-cles [11] recently.

1.2.2 Elliptic Potential

Let d × d matrix M be symmetric, positive define. Let v ∈ Rd _{and c ∈ R be arbitrary.}

The elliptic potential [2] defines as

V (u) = 1

2u

T

M u + uTv + c,

where u is a column vector in Rd. Because M is symmetric and psoitive definite, M12

exists. Therefore, we can also write V (u) = 1_{2}||M12u||2 + uTv + c. In R2, the elliptic

potential can be rewritten as

V (x, y) = ax2+ bxy + cy2_{+ dx + ey + f, a, b, c, d, e, f ∈ R}
with b2_{− 4ac < 0. The shape of this potential looks like a bowl.}

### 1.3

### Goals

In this thesis, we consider nonlinear Schr¨odinger equations with focusing parameter ,

i∂tψ = −∆ψ + V (x)ψ − m(x)|ψ|2ψ, (2)

where ψ(t, x) : R × R2 _{→ C is the wavefunction, V (x) is the trap potential, m(x) is}

the mass and s-wave scattering length term, and > 0. Such equation occurs in many physics, including nonlinear optics, quantum physics, and water waves.

In this thesis, we consider two kinds of different trap potentials and domains as follows: Case 1. (ω case)

V (x) = (ω1x1)2+ (ω2x2)2,

where ω1, ω2 ≥ 0 and ω1, ω2 = 1, 2, ..., 5, and m(x) = 1 + 1_{2} sin(2πx1) with unit disk

domain Ω = {(x1, x2)|x21+ x22 < 1}.

Case 2. (µ case)

V (x) = 1 + sin2(µ1x1) + sin2(µ2x2),

where µ1, µ2 = 1, 2, ..., 5 and m(x) = 1 + 1_{2}cos(2πx1) with square domain Ω =

{(x1, x2)| − 1 < x1 < 1, −1 < x2 < 1}.

Here, we have two kinds of trap potentials and function m(x) in Fig. 3 and Fig. 2 respec-tively. We focus on different situation in ω case and µ case. We concentrate on the well depth of the potential by varying the parameters ω1 and ω2 for the elliptic trap potential.

And then we control the periodic for the optical lattices trap potential by changing the two parameters, µ1 and µ2. For keeping the shapes of trap potentials, we use different

types of domains for ω case and µ case.

The NLS equation has a special solution: ψ(t, x) = eiλt_{φ(x). The aim of this thesis}

is to find out the spectra of linearized operator which arises when the equation (2) is linearized around the special solution ψ. The goal is to study the stability of the solution by the spectrum of its linear operator.

There are many literature on NLS equation. The necessary condition for orbital sta-bility and instasta-bility of single-spike bound state can be obtained in [9, 10].

Next, in Section 2 we will present the mathematical analysis about the NLS equation. And then in Section 3 we will discuss the numerical method. At last, some numerical results will be shown in Section 4.

### 2

### Mathematical Analysis

### 2.1

### Perturbation in NLS

Recall the nonlinear Schr¨odinger equation:

i∂tψ = −∆ψ + V (x)ψ − m(x)|ψ|2ψ. (3)

We consider the special solutions of NLS equation (3): ψ(t, x) = eiλt_{φ(x), where λ =}

0, 1, and general case. φ(x) is a real-valued function independent of time. Substitution ψ(t, x) = eiλtφ(x) into (3) and the φ satisfy the nonlinear elliptic equation:

−∆φ + V φ − m|φ|2_{φ = −λφ,}

that is

(−∆ + V − m|φ|2)φ = −λφ (4)

or

(−∆ + (V + λ) − m|φ|2)φ = 0. (5)

Equation (4) and (5) will be used in different solution form and algorithm.

To study the stability of the special solution form, we consider solutions of the NLS equation of the form

ψ(t, x) = eiλt(φ(x) + h(t, x)). (6)

The perturbation h(t, x) ∈ R satisfies an equation:

∂th = Lh + (nonlinear terms),

where

Lh = −i{(−∆ + (V + λ) − 3mφ2_{)h}.} _{(7)}

Remark. To explain how to obtain (7):

In equation (3), replace ψ(x) by eiλt(ψ(x) + h(t, x)) then the left hand side of (3): i∂tψ = i∂t(eiλt(φ + h))

= i(iλeiλt(φ + h) + eiλt∂th)

= eiλt(−λ(φ + h) + i∂th).

And the right hand side of (3):

−∆ψ + V (x)ψ − m(x)|ψ|2_{ψ = e}iλt_{(−∆φ − ∆h + V (x)ψ + hφ − m|φ + h|}2_{(φ + h))}

= eiλt{−∆φ + V φ − ∆h + hφ

−m(|φ|2_{φ + φ}2_{(¯}_{h + h) + φh¯}_{h + φ}2_{h + φ(h + ¯}_{h)h + h}2_{¯}_{h)}}

= eiλt(−∆φ + V φ − m|φ|2φ)

+eiλt(−∆h + V h − m|φ|2(h + ¯h) − m|φ|2h) + eiλtO(h2).
Notice that |φ + h|2 _{= (φ + h)(φ + h) = |φ|}2_{+ φ(¯}_{h + h) + h¯}_{h. Then}

eiλt(−λ(φ + h) + i∂th) = eiλt(−∆φ + V φ − m|φ|2φ)

So

i∂th = (−∆φ + (V + λ)φ − m|φ|2φ)

+(−∆h + (V + λ)h − m|φ|2(h + ¯h) − m|φ|2h) + O(h2) = −∆h + (V + λ)h − m|φ|2(h + ¯h) − m|φ|2h + O(h2). If we consider h(t, x) as a complex perturbation, ie. h(t, x) ∈ C, then Lh = L(Re h + i Im h)

= −i{(−∆ + (V + λ) − mφ2_{) Re h + i(−∆ + (V + λ) − mφ}2_{) Im h − 2mφ}2 _{Re h}}

= (−∆ + (V + λ) − mφ2) Im h + i{−(−∆ + (V + λ) − 3mφ2) Re h}.

And we rewrite L as a matrix acting on Re h

Im h , L = 0 L− −L+ 0 , (8) where L+ = −∆ + (V + λ) − 3mφ2, L−= −∆ + (V + λ) − mφ2. (9)

L+ and L− are self-adjoint.

In the following Lemma, we will explain how we could rewrite the linerization term as
(8). Are eigenvalues and eigenvectors of (8) and the original form the same or not?
Lemma 1. Let L =
0 L−
−L+ 0
and Lh = −∆h + (V + λ)h − m|φ|2(h + ¯h) −
m|φ|2_{h, where L}

+ and L− are defined by (9). If ρ is an eigenvalue of L and [uT, vT]T

be the eigenvector, then ρ is also an eigenvalue of L and u + iv is an eigenvector of L corresponding to ρ.

Proof. Since ρ is an eigenvalue of L, that is,
L
u
w
= ρ
u
w
. (10)
So
0 L−
−L+ 0
u
w
= ρ
u
w
.
Then
−(−∆ + (V + λ) − 3mφ2_{)u = ρw} _{(11)}
and
(−∆ + (V + λ) − mφ2)w = ρu. (12)
(12)+i×(11), therefore

ρ(u + iw) = −i(−∆ + (V + λ) − 3mq2)u + (−∆ + (V + λ) − mq2)w

= −i{(−∆ + (V + λ) − 3mq2)u + i(−∆ + (V + λ) − mq2)w}

= −i{(−∆(u + iw) + (V + λ)(u + iw) − mq2(u + iw) − 2mq2u}.

### 2.2

### Stability

Our aim of this thesis is to study the stability of these solutions. We show the stability by solving eigenvalue problem. So in this subsection, we recall the definition of some stability, and the relation between eigenvalues and stability [4, 5].

At first, we give some notations and definitions. Consider an ODE dx

dt = ˙x = f (x); x ∈ R

n_{,} _{(13)}

where f : Rn _{→ R}n _{and x = [x}

1, x2, ..., xn]T. If f (x∗) = 0 for all t, the point x∗ is called

an equilibrium point.

Definition 1. x∗ is a Lyapunov stable equilibrium if for every neighborhood U of x∗ there is a neighborhood V ⊆ U of x∗ such that every solution x(t) with x(0) = x0 ∈ V is defined

and remains in U for all t ≥ 0.

Definition 2. If V can be chosen above so that, in addition to the properties for stability, we have limt→∞x(t) = x∗ then we say that x∗ is asymptotically stable.

An equilibrium is called neutrally stable if it is Lyapunov stable but not asymptotically stable.

In studying the stability of x∗, we consider x∗ plus a small perturbation h(t), ie,
x(t) = x∗+ h(t), where |h(t)| << 1. Subtitute x(t) into (13) and expand f (x) by Taylor
series: ˙x∗+ ˙h = f (x∗+ h) = f (x∗) + Df (x∗)h + O(|h2_{|). The notation Df (x}∗_{) is the n × n}

Jacobian matrix of partial derivative of a vector-valued function f .

The eigenvalues and eigenvectors of the matrix Df (x∗) determine the general solution. In studying stability we want to know whether the solution grows, stays constant, or decay to 0 as t → ∞. It can be answered by evaluating the eigenvalues.

If λ is a real eigenvalue with eigenvector v, then there is a solution to the linearization

equation of the form: h(t) = cveλt. If λ = a ± ib is a complex conjugate pair with

eigenvectors v = u ± iw (where u, w are real), then h1(t) = eat(u cos bt − w sin bt) and

h2(t) = eat(u cos bt + w sin bt) are two linearly independent solutions. In both cases, the

real part of λ almost determines stability. Any solution of the linearized equation can be written as a linear combination of terms of these forms. We can obtain the following conclusions:

1. If all eigenvalues of Df (x∗) have negative real parts, then |h(t)| → 0 as t → ∞ for all solutions.

2. If there exists one eigenvalue of Df (x∗) has a positive real part, then there is a solution h(t) with |h(t)| → +∞ as t → ∞.

3. If some pair of complex-conjugate eigenvalues have zero real parts with distinct imaginary parts, then the corresponding solutions for |h(t)| oscillate as t → ∞ and neither decay nor grow as t → ∞.

Moreover, if x∗ is an equilibrium of ˙x = f (x) and all the eigenvalues of the matrix Df (x∗) have strictly negative real parts, then x∗is asymptotically stable. If at least one eigenvalue has strictly positive real part, then x∗ is unstable.

We had discussed the solution of the ODE is an equilibrium. The statement of stability may be extend to non-constant orbits of ODE.

Here, let Ot(x) = x(t) and the initial value x(0) = x0. Then the set O(x) = {Ot(x) :

0 ≤ t} is called orbit.

Definition 3. Let two orbits O(x) and O(ˆx). If there is a reparameterization of time ˆt(t) such that |Ot(ˆx) − Oˆt(t)(ˆx)| < for all t ≥ 0,then we say two orbits O(x) and O(ˆx) are

-close.

Definition 4. An orbit O(x) is orbitally stable if for any > 0, there is a neighborhood V of x so that, for all ˆx ∈ V , O(x) and O(ˆx) are -close. If additionally V may be chosen so that, for all ˆx ∈ V , there exists a constant τ (ˆx) so that |Ot(ˆx) − Ot−τ (ˆx)(ˆx)| < as

t → ∞, then Ot(x) is asymptotically stable.

The linearization skill in general orbit is similar to the the pervious linearization method. Also discuss the eigenvalues of the linearization operator.

### 3

### Numerical Method

In this section, we discretized the Laplace operator by finite-difference method first, and we will show it on two kinds of domain respectively. Then we present the numerical method for solving φ(x) from (5) first. And finally, we will describe a numerical method to compute the spectra of L from (7).

Now, we recall our main question: to study the stability of the special solution form.
For our solution form ψ(t, x) = eiλt_{φ(x), we will solving the time-inedpent term φ(x).}

By (5), φ = φ(x1, x2) satisfies − ∂ 2 ∂x2 1 φ + ∂ 2 ∂x2 2 φ + (V (x1, x2) + λ)φ − m(x1, x2)φ3 = 0. (14)

The natural boundary condition is lim

|x|→∞φ(x1, x2) = 0. (15)

Consider the same equation (5) in the unit disk domain Ω = {(x1, x2) : x21+ x22 < 1},

applying the polar coordinate transformation,

x1 = r cos θ, x2 = r sin θ,

where r = px2

1+ x22, θ = tan
−1_{(x}

2/x1). We can rewrite (14), φ = φ(r, θ) in the polar

coordinate system satisfies − 1 r ∂ ∂r(r ∂φ ∂r) + 1 r2 ∂2φ ∂θ2 + (V (r, θ) + λ)φ − m(r, θ)φ3 = 0, (16)

with 0 < r < 1, 0 ≤ θ < 2π. And the boundary condition is lim

### 3.1

### Matrix Form of Laplace operator

Now, we are going to describe the matrix form of the Laplace operator. In this study, we have two kinds of domains, unit disk domain and square domain. We will show them respectively.

3.1.1 Square Domain

The Laplace operator is a seconed order differential operator. Here, we consider the
function φ in R2. Then
∆φ = ∂
2_{φ}
∂x2 +
∂2_{φ}
∂y2.

Here we use Finite-Difference method for the boundary value problem [1].

First step is to partition x into −1 = x0 < x1 < x2 < · · · < xN −1 < xN = 1, and

partition y into −1 = y0 < y1 < y2 < · · · < yN −1 < yN = 1, where xi = −1 + i∆x,

yj = −1 + j∆y and h = ∆x = ∆y = 2/N . We use the Taylor series in the variable x, y

about xi, yj to generate the centered-difference formula

∂2φ ∂x2 =

fi−1,j − 2fi,j + fi+1,j

(∆x)2 + O((∆x)

2_{),}

∂2_{φ}

∂y2 =

fi−1,j − 2fi,j + fi+1,j

(∆y)2 + O((∆y)

2

)

for each i = 1, 2, ..., N − 1 and j = 1, 2, ..., N − 1, where fi,j = φ(xi, yj).

In difference-equation form, this result in the finite-difference method, with local
trun-cation error of order O((∆x)2_{+ (∆y)}2_{). Therefore}

−∆φ ' − fi−1,j − 2fi,j + fi+1,j

(∆x)2 +

fi,j−1− 2fi,j+ fi,j+1

(∆y)2

= −fi−1,j + fi+1,j− 4fi,j+ fi,j−1+ fi,j+1

h2 .

We can rewrite it as a matrix representation:

1 h2 ˆ A −I 0 −I Aˆ −I . .. ... ... −I Aˆ −I 0 −I Aˆ f1 f2 .. . fN−2 fN− 1 ≡ AF, where ˆ A = 4 −1 0 −1 4 −1 . .. ... ... −1 4 −1 0 −1 4 (N −1)×(N −1) , fj = [f1,j, f2,j, ..., fN −1,j]T,j = 1, 2, ..., N − 1, and I is an (N − 1)2 × (N − 1)2 identity matrix.

3.1.2 Unit disk Domain

Then we consider the equation (5) over a 2-dimensional disk. We use a discretization scheme [8] for equation (16). The Laplace operator is

∆φ(r, θ) = ∆rφ + ∆θφ = 1 r ∂ ∂r r∂φ ∂r + 1 r2 ∂2φ ∂θ2.

We discretize the Laplace operator in two parts. The first step is to discrstize the operator ∆rφ on r ∈ [0, 1], and then to discretize the operator ∆θφ on θ ∈ [0, 2π].

Firstly, we consider the operator −∆rφ on r ∈ [0, 1],

−∆rφ = − 1 r∂r r∂φ ∂r . To partition r into 0 = r0 < r1 2 < r1 < · · · < rN − 1 2 < rN = 1, rj = j∆r and ∆r = 1 N, and

to partition θ into 0 = θ0 < θ1 < θ2 < · · · < θM −1 < θM = 2π, θj = j∆θ and ∆θ = 2π_{M},

and denote fi,j = φ(ri, θj). For each θj,

−∆rφ '

"

− 1

r_{i−}1
2

ri(fi+1_{2},j− fi−1_{2},j) − ri−1(fi−1_{2},j − fi−3_{2},j)

(∆r)2
#
.
That is
1
(∆r)2
r1
r1
2
−r1
r1
2
0
−r1
r3
2
1
r3
2
(r1 + r2) −_{r}r2_{3}
2
. .. . .. . ..
−rN −2
r_{N − 3}
2
1
r_{N − 3}
2
(rN −2+ rN −1) −_{r}rN −1
N − 3_{2}
0 −rN −1
r_{N − 1}
2
1
r_{N − 1}
2
(rN −1+ rN)
f1
2,j
f3
2,j
..
.
f_{N −}3
2,j
f_{N −}1
2,j
≡ bAf ,
where
b
A = 1
(∆r)2
a1 c1 0
b1 a2 c2
. .. ... . ..
bN −2 aN −1 cN −1
0 bN −1 aN
and fj =
f1
2,j
f3
2,j
..
.
f_{N −}1
2,j
,
ai =
1
r_{i−}1
2
(ri−1+ ri) =
i − 1 + i
i −1_{2} = 2, i = 1, . . . , N
and
bi = −
ri
r_{i+}1
2
= −1 + 1
2i + 1, ci = −
ri
r_{i−}1
2
= −1 − 1
2i − 1, i = 1, . . . , N − 1.

Now we consider the operator
−∆θφ = −
1
r2
∂2_{φ}
∂θ2

on θ ∈ [0, 2π]. Then for each r_{i−}1

2, i = 1, 2, ..., N
−∆θφ '
"
− 1
r2
i−1_{2}
f (r_{i−}1
2, θj+1) − 2f (ri−
1
2, θj) + f (ri−
1
2, θj−1)
(∆θ)2
#
.
Let
B = diag 2
r2
1
2
(∆θ)2,
2
r2
3
2
(∆θ)2,
2
r2
5
2
(∆θ)2, . . . ,
2
r2
N −1_{2}(∆θ)
2
!
and
C = diag −1
r2
1
2
(∆θ)2,
−1
r2
3
2
(∆θ)2,
−1
r2
5
2
(∆θ)2, . . . ,
−1
r2
N −1
2
(∆θ)2
!
,
fj = [f1
2,j, f
3
2,j, . . . , fN −
1
2,j]
>
for j = 1, 2, ..., M.
Then we can rewrite the Laplace operator as

b A + B C C C A + Bb C . .. . .. . .. C A + Bb C C C A + Bb f1 f2 .. . fM−1 fM ≡ AF.

Then the matrix A represents the Laplace operator.

### 3.2

### Algorithm

Now, we describe a numerical method to compute the spectrum of the linear operator L defined by (7) for > 0. There are two steps in our numerical method: First, we will solve the nonlinear problem (14) for φ. Second, compute the spectrum of the discretized linearized operator around φ.

Step I. Compute φ and it is energy minimizing among all solution of (14)-(17).

Step II. Compute the spectra of the discretized linearized operator L. Since h is no
longer real, so the linearized operator is a little different from L in the equation (7).
Here, we denote some notation of the following. For A ∈ RM ×M_{, q, m ∈ R}M_{, m ∗ q}

denote the Hadamard product of m and q, and q
p _{means the p-time Hadamard product}

of q. And diag(q) represents the diagonal matrix of q. In Step I. First, we represent −∆q as

where q is an approximation of φ. As for the form of matrix A, it was shown in previous two subsections. Then the discretization of the equation (14) (and (16)) become

Aq + (V + λ) ∗ q − m ∗ q 3

= 0,

where V and m are approximation of the function V (x1, x2) and m(x1, x2), respectively.

We solve it by using an iteration method [7]:

[A + diag(V + λ)]˜qnew = m ∗ q 3

old, (18)

where ˜qnew and qold are unknown and known vector. The iteration step is shown in

Algorithm 1 .

Algorithm 1 Iterative algorithm for solving φ(x)

Step 1 Choose an initial guess of qold> 0, and ||qold||2 = 1.

Step 2 Solve (18), then obtain ˜qnew.

Step 3 Let αnew = _{||˜}_{q}_{new}1 _{||}

2, qnew = αnewq˜new.

Step 4

if (converge) then

Output the solution (αnew)

1

2q_{new}. Stop.

else

Let qold= qnew.

Go to Step 2. end if

When λ in a general case, λ is no longer a constant. λ varies from , so the numerical method has a little difference form. From (4), the discretization form becomes

[A + diag(V − m ∗ q 2

)]q = −λq. (19)

It turns into an eigenvalue problem. In Algorithm 2, it shows the iteration step. In Step II. Now we discretize L of (8) into an eigenvalue problem.

L
u
w
= ρ
u
w
, (20)
where
L =
0 A + diag(V + λ) − diag(m ∗ q2_{)}

−A − diag(V + λ) + diag(3m ∗ q2_{)} _{0}

. q and λ are obtain from Step I. We use ARPACK in MATLAB version R2007a to solve the linear algebraic eigenvalue problem and obtain eigenvalues ρ of L near origin.

### 4

### Numerical Simulation

For each potential case, we summary the numerical results for three solution forms.

Algorithm 2 Iterative algorithm for solving φ(x) in general solution Step 1 Choose an initial guess of qold> 0, and ||qold||2 = 1.

Step 2 Solve

(A + diag(V − m ∗ q 2

old))˜qnew = −λ˜qnew,

where −λ is the smallest eigenvalue of A + diag(V − m ∗ q 2

old). Then obtain ˜qnew.

Step 3 Let ˆq = (A + diag(V − m ∗ ˜q 2

new))˜qnew, αnew = |−λ|

||ˆq||2, qnew = αnewq˜new.

Step 4

if (converge) then

Output the solution qnew and λ. Stop.

else

Let qold= qnew.

Go to Step 2. end if

### 4.1

### Setting

In ω case, the discretized matrix of Laplace operator has size N T by N T with N = r/∆r
and T = 2π/∆θ. Here we use zero boundary condition and r = 1, N = 32, T = 64. And
in µ case with square domain, the discretized natrix of Laplace operator has size N2 _{by}

N2 with N = 2/h and here N = 64. And so the matrix size of the operator L are 2N T

by 2N T and 2N2 _{by 2N}2 _{respectively.}

In our numerical method, we use finite-difference method to solve NLS equation. The
finite-difference method has truncation error O(h2), where h is the grid size. In our
experiment, O(h2_{) ≈ c · 10}−3_{.}

The region of that we consider about is larger than 10−5, that is because the value is smaller than 10−5 we treat it as zero. We should control the region of to satisfy the boundary condition, which means the region of would not be large. In the region that we considered, the boundary mean values of all cases are less than 10−5. In Fig. 4, we plot the mean value of the boundary of φ for from 0 to 0.003. We can see the mean of boundary will be less than 10−5 as smaller than 0.0024. So in that case, we only try ∈ [10−5, 0.0024] in our numerical experiment.

Besides the boundary condition, we also take care of the shape of φ, here we narrow down the focus on the solitary solution for this part only. As goes larger, sometimes the shape would change. The following statement like ∗ are all in the region of .

In the following table(Table 1), we write down the region of for each case. In each potential case, for our converient, we choose the intersection of region for different parameters.

In this thesis, we focus on the parameter changing, and to study when the solution

of the NLS would be unstable or stable. The notation ∗ denote as the bifurcation of

stable and unstable.

We find that a pair of purely imaginary eigenvalues will collide at the origin and split into a pair of real eigenvalues for < ∗. That is, when less than ∗, the spectrum of the linearized operator are all pure imaginary. As larger than ∗, the spectrum of the linearized operator has at least one pair of eigenvalue with nonzero real part. That means

Table 1: The region of .

case µ case ω

λ = 0 [10−5, 0.016] [10−5, 0.0024] λ = 1 [10−5, 0.025] [10−5, 0.005] general λ [10−5, 0.001] [10−5, 0.001]

Figure 8: Solution of φ for ω1 = 1, ω2 = 1, = 0.001, λ = 1.

the solution is stable while < ∗, unstable when > ∗.

Since the numerical result of eigenvalues are approximartion of the exact mathematical values. There exist some error of the result. We determined that the eigenvalue has nonzero real part when the absolute value of the real part of eigenvalus is larger than 10−3.

Since the regions of are small, and the spectrum of lineraization operation changes so fast. We move in a small step for each case. δ = 10−5. To prevent losing some pairs of pure complex eigenvalue from turning to real eigenvalues, we also move the target of the algorithm for finding eigenvalues.

This method is not an efficient method but an easy way to catch the changing of spectra.

### 4.2

### ω case

In ω case, the maxima of φ is not at center of the disk; instead, it is on the right of the center. As go larger and the potential change, the maxima of φ seem to moving to the center of domain. In figure 10, we fix = 0.0002 and changing the well depth of potential V (x). We find that the location of maxima φ changes. It moves to the center of the domain.

We know trap potential V (x) and m(x) would effect the solution φ. The NLS equation in our model is an focusing case, that is the atoms would stay near the maxima of m(x), ie, x1 = 1/4.

1. λ = 0: Testing in the region that we mentioned before, and find that the eigen-values of the linearized operation L are all pure imaginary.

2. λ = 1: We find there is an ∗ in this solution case. The ∗ represent the bifurcation of pure imaginary to has a pair of eigenvalue with nonzero real part. In Fig. 11

shows the spectrum of L . And in Fig. 12, it shows the ∗ of the parameters

Figure 9: Solution of φ for ω1 = 1, ω2 = 1, = 0.001, λ = 0.

3. In elliptic potential case, we can find ∗in eit_{φ(X) this solution form in our numerical}

method, but we can not find the ∗ in other solution form. At first we guess ∗ in ei0tφ(x) form of solution can be found in smaller . But we do not find that in the region of as we mentioned before.

4. λ is general case: In this solution case, we can not find the ∗. When is quite small, the value of λ is already negative. That means the trap potential is not positive at all, V (x) < 0 on the center of the disk. In Figure 13 it shows the value of λ as changes.

Figure 14: Spectrum of µ case for λ = 0.

### 4.3

### µ case

1. λ = 0, we can find ∗ in this solution case. For µ1, µ2 = 1, 2, ..., 5, we find that the

eigenvalues are pure imaginary when < ∗, and turn into a pair real value. Figure 14 plot the eigenvalues around original. And in Figure 15 show the ∗ for every µ1, µ2.

2. λ = 1, we also find ∗. Figure 16, 17 show the eigenvalues and values of ∗ for each µ1, µ2.

3. Consider the same µ1 and µ2, these two solution forms can be seen as changing

the well depth of potential (from (5)). The ∗ are different. So we know that the potential well depth may change the stability of solution.

In optical lattices potential case, changing µ1or µ2for different periodic of potential.

We can see that the value of ∗ increasing when µ1, µ2 go larger. To comparison

eit_{φ(x) and e}i0t_{φ(x) these two forms of solutions, we can see that the second kind}

of ∗ is less than the other ones.

4. λ is general case, there is no ∗ in this solution case. When is quit small( less than 1.2 × 10−4), the value of λ is positive(In Figure 18). This result is consistent with the condition of Lin in [9](λ is positive for small enough ). But λ become negative when gets larger, and the value of λ gets smaller as gets larger. When < 3.3 × 10−4, λ is smaller than −1. And the minimum of trap potential is 1, that means the trap potential is no longer positive value.

### 5

### Conclusions

In this thesis, we consider the NLS equation with the special solution eiλt_{φ(x).} _{We}

mainly show the spectrum of the linerization operator L. By numerical computation, we can conclude some results: when < ∗ the eigenvalues are pure imaginary, and the eigenvalues turn to real as > ∗. That is, there is a bifurcation of stable and unstable, when λ = 0, 1 in optical lattices potential case and λ = 1 in elliptic potential case.

Here, we use Matlab code to find the spectrum of the linearization operator. It is a convenient way to solve the eigenvalue problem. Since the method of finding eigenvalue in Matlab is to search some eigenvalues which are near the target. And the distribution of eigenvalues changes fast when near ∗. For future work, we are interested in improving the method of finding the spectrum of the linearization operator.

### References

[1] R. Burden and J. D. Faires. Numerical Analysis. Brooks Cole, 7 edition, 2000. [2] N. Cesa-Bianchi and G. Lugosi. Prediction, Learning, and Games. Cambridge

Uni-versity Press, 2006.

[3] S. M. Chang, S. Gustafson, K. Nakanishi, and T. P. Tsai. Spectra of linearized oper-ators for NLS solitary waves. SIAM Journal on Mathematical Analysis, 39(4):1070– 1111, 2007.

[4] M. W. Hirsch, S. Smale, and R. L. Devaney. Differential equations, dynamical sys-tems, and an introduction to chaos. Academic Press, 2004.

[5] P. Holmes and E. T. Shea-Brown. Stability. Scholarpedia, 1(10):1838, 2006.

[6] http://en.wikipedia.org/wiki/Optical lattice. Optical lattice- wikipedia, the free en-cyclopedia.

[7] T. M. Hwang and W. Wang. Analyzing and visualizing a discretized semilinear elliptic problem with neumann boundary conditions. Numerical Methods for Partial Differential Equations, 18(3):261–279, 2002.

[8] M. C. Lai. A note on finite difference discretizations for poisson equation on a disk. Numerical Methods for Partial Differential Equations, 17:199–203, 2001.

[9] T. C. Lin and J. Wei. Orbital stability of bound states of semi-classical nonlinear

schr¨odinger equations with critical nonlinearity. SIAM Journal on Mathematical

Analysis, 40(1):365–381, 2008.

[10] T. C. Lin, Juncheng Wei, and Wei Yao. Orbital stability of bound states of nonlinear schr¨odinger equations with linear and nonlinear lattices.

[11] M. P. MacDonald, G. C. Spalding, and K. Dholakia. Microfluidic sorting in an optical lattice. Nature, 426:421–424, 2003.

[12] C. Sulem and P.-L. Sulem. The nonlinear Schrodinger equation: Self-Focusing and Wave Collapse. Springer, 1999.