E ( u ( β )) ∗ Eigenvalues λ ∗ 1 ( β ) ,λ ∗ 2 ( β )Energy

Verticillate Structures in Multi-Component Bose-Einstein Condensates

Shu-Ming Chang^a

Department of Mathematics National Tsing Hua University

Hsinchu, Taiwan

December 18, 2004

aThis is a joint work with Chang-Shou Lin (NCCU), Tai-Chia Lin (NTU), Wen-Wei Lin (NTHU) & Shih-Feng Shieh (NTHU).

Outline

• Introduction of Bose-Einstein Condensation (BEC)

• Coupled Nonlinear Schr¨odinger Equations and Coupled Gross-Pitaevskii Equations (CGPEs)

• Nonlinear Algebraic Eigenvalue Problems (NAEPs)

• Fixed Point Iteration for NAEPs

• Numerical Results

• Conclusions

Results

• Theoretical

– Propose a Jacobi-type fixed point iteration (J-FPI) and a Gauss-Seidel-type fixed point iteration (GS-FPI) to solve Multi-Component BEC.

– Prove that the GS-FPI method converges locally and linearly to a fixed point if and only if the associated

minimized energy functional problem has a strictly local minimum at the feasible fixed point.

• Numerical

– Simulate multi-component BEC.

– A New phenomenon: verticillate multipling structure.

1 Introduction of BEC

• What is BEC?

A new form of matter at the coldest temperatures in the universe...

• Theoretical prediction 1924 ...

– S. Bose: derived Planck’s black body radiation law from considering the cavity radiation as an ideal photon gas and worked out Bose statistics for photons.

– A. Einstein: generalized Bose statistics to other Bosonic particles and atoms (Bose-Einstein statistics) and predicted if the atoms were cold enough, almost all of the particles would congregate in the ground states (BEC).

A. Einstein S. Bose

• Physical experiments – Superfluid He⁴ 1938:

P. L. Kapitza, Allen and Misener: discovered the superfluidity of liquid helium.

F. London: proposed that the superfluid fraction consisting of those atoms which have “condensed” to the ground state.

P. L.Kapitza F. London

• – E. A. Cornell & C. E. Wieman (JILA, 1995):

first observed BEC of rubidium (⁸⁷Rb) atoms at 20 nK, i.e.

0.000 000 02 K.

C. E. Wieman & E. A. Cornell BEC at 400, 200, and 50 nK

– W. Ketterle (MIT, 1995):

observed BEC of sodium (²³Na) atoms.

W. Ketterle Two-Component BEC

• Experimental implementation

– The BEC named Science Magazine’s ”Molecule of the Year 1995”!

– Nobel Prize in Physics (2001), E. A. Cornell, C. E. Wieman (JILA), W. Ketterle (MIT):

for the achievement of BEC in dilute gases of alkali atoms, and for early fundamental studies of the properties of the condensates.

• Applications of BEC: atom laser, quantum computer, MEMS.

• Mathematical model: nonlinear Schr¨odinger equations, Gross-Pitaevskii equations (GPEs), coupled nonlinear

Schr¨odinger equations, coupled Gross-Pitaevskii equations (CGPEs).

2 Coupled Nonlinear Schr¨ odinger Eqs.

and CGPEs

ι~∂ψ_j(x, t)

∂t = − ~²

2m_a∇²ψ_j+ bV_jψ_j+µ_jj|ψ_j|²ψ_j+X

j6=i

µ_ij|ψ_i|²ψ_j, j = 1, . . . , m.

• Coupled Gross-Pitaevskii equations (CGPEs):

ι~∂ψ(x, t)

∂t = − ~²

2m_a ∇²ψ(x, t) + bV (x) ◦ ψ(x, t) + bB(ψ) ◦ ψ(x, t), (2.1) ψ(x, t) = (ψ₁(x, t), . . . , ψ_m(x, t))^⊤: wave ft.,

Vb (x) = ( bV₁(x), . . . , bV_m(x))^⊤: trap potential,

Bb (ψ) = ( bB₁(ψ), . . . , bB_m(ψ))^⊤: intra- and inter-species scattering

Vb_j(x) = m_a

2 ω_x,j² (x − ˆx_0,j)² + ω_y,j² (y − ˆy_0,j)² + ω_z,j² (z − ˆz_0,j)²
. µ_jℓ = ^4π~_m^2bjℓ

a , b_jℓ: scattering length. (+: repulsive, −: attractive) Z

D

|ψ_j(x, t)|²dx = N_j⁰ > 0, j = 1, . . . , m.

Assume that

(i) ω_x,1 ≤ · · · ≤ ω_x,m ≤ ω_y,1 ≤ · · · ≤ ω_y,m ≤ ω_z,1 ≤ · · · ≤ ω_z,m. (ii) b_jℓ = b_ℓj, j, ℓ = 1, . . . , m.

• Dimensionless

˜t = t

t_s, x˜ = x

x_s, ψ˜_j(˜x, ˜t) = x^3/2_s qN_j⁰

ψ_j(x, t).

t_s = 1

ω_x,1 : dimensionless “time”.

x_s =

s ~

m_aω_x,1 : dimensionless “length”.

• Dimensionless CGPEs

ι∂ψ(x, t)

∂t = −1

2∇²ψ(x, t) + V (x) ◦ ψ(x, t) + B(ψ) ◦ ψ(x, t), (2.2)

ψ(x, t) = 0, x ∈ D, (2.3)

D: A smooth bounded domain in R^d, d = 2, 3, n(ψ_j) :=

Z

D

|ψ_j(x, t)|²dx = 1, j = 1, . . . , m,

B(ψ) = (B₁(ψ), . . . , B_m(ψ))^⊤, B_j(ψ) = β_j1|ψ₁|² + · · · + β_jm|ψ_m|².

ψ(x, t) = (ψ₁(x, t), . . . , ψ_m(x, t))^⊤, V (x) = (V₁(x), . . . , V_m(x))^⊤,

V_j(x) = 1

2(γ_x,j² (x − x_0,j)² + γ_y,j² (y − y_0,j)² + γ_z,j² (z − z_0,j)²), γ_x,j = ω_x,j

ω_x,1, γ_y,j = ω_y,j

ω_x,1, γ_z,j = ω_z,j ω_x,1 , x_0,j = xˆ_0,j

b₀ , y_0,j = yˆ_0,j

b₀ , z_0,j = zˆ_0,j b₀ , β_j,ℓ = µ_jℓN_ℓ⁰

b³₀~ω_x,1 = 4π~²b_jℓN_ℓ⁰

m_ab³₀~ω_x,1 = 4πb_jℓN_ℓ⁰ b₀ .

Energy

E(ψ) =

Xm j=1

N_j⁰

N₀ E_j(ψ),

where N_j⁰ > 0 is the number of particles with P^m

j=1

N_j⁰ = N⁰ and

E_j(ψ) = Z

D

"

1

2|∇ψ_j|² + V_j(x)|ψ_j|² + 1 2

Xm k=1

β_j,k|ψ_j|²|ψ_k|²

# dx, for j = 1, . . . , m.

Let ψ(x, t) = e^−ιλ(c)t ◦ φ(x), where λ^(c) = (λ^(c)₁ , . . . , λ^(c)m )^⊤, φ(x) = (φ(x), . . . , φ_m(x))^⊤.

Substituting ψ(x, t) into CGPEs gives the NEP for (λ, φ):

λ^(c) ◦ φ(x) = −1

2∇²φ(x) + V (x) ◦ φ(x) + B(φ) ◦ φ(x), (2.4) with

Z

D

|φ_j(x)|²dx = 1, j = 1, . . . , m, where

B(φ) = (B₁(φ), . . . , B_m(φ))^⊤, B_j(φ) =

Xm k=1

β_jk|φ_k|².

Multiplying the j-th eq. in NEP (2.4) by φ_j(x), the eigenvalue λ^(c)_j and the corresp. eigenfunction φ_j for (2.4) satisfy

λ^(c)_j = Z

D

"

1

2|∇φ_j|² + V_j(x)|φ_j|² +

Xm k=1

β_jk|φ_j|²|φ_k|²

# dx

= E_j(φ) + 1 2

Z

D

Xm k=1

β_jk|φ_j|²|φ_k|²dx.

The ground state φ_g(x) of multi-comp. BEC can be found by minimizing E(φ):

Minimize

φ=(φ1,...,φm)^⊤E(φ) subject to

Z

D

|φ_j(x)|²dx = 1, j = 1, . . . , m,

(2.5)

where E(φ) = Pm j=1

N_j⁰

N⁰E_j(φ) with E_j(φ) =

Z

D

"

1

2|∇φ_j|² + V_j(x)|φ_j|² + 1 2

Xm k=1

β_jk|φ_j|²|φ_k|²

# dx.

The NEP (2.4) can be regarded as the Euler-Lagrange eq. of the opt. problem (2.5).

3 Nonlinear Algebraic Eigenvalue Problems (NAEPs)

For computational purpose, we derive the discretization of NEP and the associated opt. problem. We consider D ⊆ R² a bounded domain.

The central finite difference discretizes −∇²φ_j(x) into Au_j = A[u_j1, . . . , u_jl, . . . , u_jN]^⊤, A ∈ R^{N ×N}, where u_j is an approx. of the j-th wave ft. φ_j(x).

Ab and bA^⊤ are irreducible and diag.-dominant with positive diag.

and nonpositive off-diag. elements. bA is symmetrizable to a s.p.d.

A by a D > 0, i.e.,

Ab = D⁻¹AD, A^⊤ = A ≻ 0.

1

2 1 ³₂ 2 ⁵₂ 3 ⁷₂ 4

Note that D = diag(d_l1,l2) with d²_l1,l2 = (l₁ − ¹₂)δr²δθ is equal to the area of the (l₁ + ν(l₂ − 1))-th sector.

Applying −∇² ≈ A to NEP and rewriting u_j ≡ hu_j, β_jk = ^β_h^jk2 , the discretization of NEP, referred as a NAEPs, can be formulated by

1

2Au_j + V _j ◦ u_j +

Xm k=1

β_jku _k² ◦ u_j = λ^(c)_j u_j, (3.1) u^⊤_j u_j = 1, j = 1, . . . , m, (3.2) where V _j = [V_j, . . . , V_j], j = 1, . . . , m. h is the grid size.

Let u = (u^⊤₁ , . . . , u^⊤_m)^⊤. Since the j-th kinetic energy Z

D

1

2|∇φ_j|²dx = − Z

D

φ_j(∇²φ_j)dx,

we approximate it by ¹₂u^⊤_j Au_j. Then the discretized eq. of the j-th energy E_j(φ) becomes

E_j(u) = 1

2u^⊤_j Au_j + V ^⊤_j u _j² + 1 2

Xm k=1

β_jku _k^{2 ⊤}u _j² .

Multiplying NAEPs (3.1) by u^⊤_j , the eigenvalues

λ^(c) = (λ^(c)₁ , . . . , λ^(c)m )^⊤ and the assoc. EVs u = (u^⊤₁ , . . . , u^⊤_m)^⊤ satisfy

λ^(c)_j = 1

2u^⊤_j Au_j + V ^⊤_j u _j² +

Xm k=1

β_jku _k^{2 ⊤}u _j²

= E_j(u) + 1 2

Xm k=1

β_jku _k^{2 ⊤}u _j² , j = 1, . . . , m.

Furthermore, E(u) =

Xm j=1

N_j⁰

N⁰ E_j(u) =

Xm j=1

N_j⁰

N⁰ λ^(c)_j − 1 2

Xm k=1

β_jku _k^{2 ⊤}u _j²

! .

For convenience, we let N_j⁰/N⁰ = 1/m. Then

E(u) = 1 m

Xm j=1

E_j(u) = 1 2m

Xm j=1

u^⊤_j Au_j + 1 m

Xm j=1

V ^⊤_j u _j²

+ 1 2m

Xm j=1

β_jju _j^{2 ⊤}u _j² + 1 m

X

1≤j<k≤m

β_jku _k^{2 ⊤}u _j² . The discretization of the opt. problem (2.5) becomes

Minimize

u=(u^⊤₁ ,...,u^⊤_m)^⊤E(u)

subject to u^⊤_j u_j = 1, j = 1, . . . , m.

(3.3)

Applying the optimality condition to the opt. problem (3.3), a local minimum ((λ^(L)₁ , . . . , λ^(L)_m ), (u^⊤₁ , . . . , u^⊤_m)^⊤) satisfies the Karash-Kuhn-Tucker (KKT) eq.

1 m



A + [[V _j]] + 2β_jj[[u _j² ]]

u_j + 2 m

X

k6=j

β_jku _k² ◦ u_j = 2λ^(L)_j u_j, (3.4) where {λ^(L)_j }^m_j=1 are Lagrange multipliers. Multiplying (3.4) by

m/2 gives 1

2Au_j + V _j ◦ u_j +

Xm k=1

β_jku _k² ◦ u_j = mλ^(L)_j u_j.

We now define

A_j := A + 2[[V _j]] + 2β_jj[[u _j² ]],

λ_j := 2mλ^(L)_j , β_jk := 2β_jk, j 6= k, for j, k = 1, . . . , m. Then the NAEPs becomes

A_ju_j + X

k6=j

β_jku _k² ◦ u_j = λ_ju_j, j = 1, . . . , m and the associated opt. problem (3.3) becomes

Minimize

u=(u^⊤₁ ,...,u^⊤_m)^⊤E(u)

subject to u^⊤_j u_j = 1, j = 1, . . . , m, where

E(u) ≡

Xm  1

2u^⊤_j Au_j + V ^⊤_j u _j²



+ 1 2

X β_jku _k^{2 ⊤}u _j² .

4 Fixed Point Iteration for NAEPs

Define

M = {v ∈ R^N|v^⊤v = 1, v ≥ 0}, M= interior of M.^◦ We suppose

β_jj > 0 small, β_jk = β_kj > 0 (j 6= k), j, k = 1, . . . , m.

A in (3.1) is diagonal dominant and Ae
0, where e = (1, . . . , 1)^⊤. For V _j ≥ 0 and (u₁, . . . , u_m) ∈ ×^m

j=1 M, the matrix A¯ _j ≡ A_j +

Xm k=1

[[β_jku _k² ]],

−1

By Perron-Frobenious Theorem there is a unique positive eigenvector ¯u_j > 0 with ¯u^⊤_j u¯_j = 1 corresp. to the maximal eigenvalue µ^max_j of ¯A⁻¹_j . i.e., ¯u_j > 0 is uniquely determined by (u₁, . . . , u_m) and satisfies

A¯_ju¯_j ≡ A_j +

Xm k=1

[[β_jku _k² ]]

!

¯

u_j = λ^min_j u¯_j, where λ^min_j = 1/µ^max_j and ¯u^⊤_j u¯_j = 1, for j = 1, . . . , m.

We now define a function f : ×^m

j=1 M → ×^m

j=1 M by f(u₁, . . . , u_m) = (¯u₁, . . . , ¯u_m), where ¯u_j > 0 is well-defined, j = 1, . . . , m.

Theorem 4.1 The function f given has a fixed point in ×^m

j=1

M. In◦

other words, there is a point (u^∗₁, . . . , u^∗_m) ∈ ×^m

j=1

M and◦

λ = (λ^∗₁, . . . , λ^∗_m) which solve the NAEPs, that is, A_ju^∗_j +

Xm k=1

β_jku^∗ _k ² ◦ u^∗_j = λ^∗_ju^∗_j, j = 1, . . . , m.

We define the restricted Lagragian function of the opt. problem by L(u) = E(u) − 1

2

Xm j=1

λ_j(u^⊤_j u_j − 1),

where

E(u) ≡ 1 2

Xm j=1

u^⊤_j A_ju_j + 1 2

X

1≤j<k≤m

β_jku _k^{2 ⊤}u _j² .

Theorem 4.2 Let u^∗ = (u^∗₁, . . . , u^∗_m) be a KKT point of the opt.

problem assoc. with the Lagrangian multipliers (λ^∗₁, . . . , λ^∗_m).

Denote the Hessian of L(u) at u^∗ by ∇²L(u^∗) = 

∇²L(u^∗)_ijm i,j=1, where

∇²L(u^∗)_jj = A_j +

Xm k=1

[[β_jku^∗ _k ² ]] − λ^∗_jI_N

!

and

∇²L(u^∗)_ij = ∇²L(u^∗)_ji = 2[[β_jiu^∗_i ◦ u^∗_j]], j 6= i, Let d = (d^⊤₁ , . . . , d^⊤_m)^⊤ ∈ R^{N m}. The positivity condition

d^⊤(∇²L(u^∗))d > 0

holds, for all d with u^∗_j^⊤d_j = 0, j = 1, . . . , m, if and only if u^∗ is a

Jacobi Iteration (JI)

Define f : ×^m

j=1 M → ×^m

j=1 M by

f(u₁, . . . , u_m) = (¯u₁, . . . , ¯u_m), where ¯u_j > 0 is well-defined, j = 1, . . . , m.

Theorem 4.3 Let (λ^∗, u^∗) = ((λ^∗₁, . . . , λ^∗_m), (u^∗₁, . . . , u^∗_m)) be a

fixed point of NAEPs. Suppose β_jj > 0 suff. small, j = 1, . . . , m. If the JI converges to (λ^∗, u^∗) locally and linearly with an initial in

×m j=1

M, then u◦ ^∗ = (u^∗₁, . . . , u^∗_m) is a strictly local min. of the opt.

problem.

Gauss-Seidel Iteration (GSI)

Define g : ×^m

j=1 M → ×^m

j=1 M by

g(u₁, . . . , u_m) = (¯u₁, . . . , ¯u_m), where

¯

u₁ = g₁(u₁, . . . , u_m) = f₁(u₁, u₂, . . . , u_m),

¯

u₂ = g₂(u₁, . . . , u_m) = f₂(¯u₁, u₂, u₃, . . . , u_m),

... ...

¯

u_m = g_m(u₁, . . . , u_m) = f_m(¯u₁, ¯u₂, . . . , ¯u_m−1, u_m),

in which {f_j}^m_j=1 are given in JI. The ft. g defines a Gauss-Seidel

Theorem 4.4 Let (λ^∗, u^∗) = ((λ^∗₁, . . . , λ^∗_m), (u^∗₁, . . . , u^∗_m)) be a fixed point of the NAEPs. Suppose the matrix Z^⊤∇²L(u^∗)Z is nonsingular. Suppose β_jj > 0 suff. small, j = 1, . . . , m. The GSI converges to (λ^∗, u^∗) locally and linearly with an initial in ×^m

j=1

M iff◦

u^∗ = (u^∗₁, . . . , u^∗_m) is a strictly local min. of the opt. problem.

J^∗_s = Ω^∗12J^∗_fΩ^∗−12 = −2 [ ∗ ],









0 β12Ω

∗− 12 1 Z∗T

1 [[u∗ 1 ◦ u∗

2 ]]Z∗ 2Ω

∗− 12

2 · · · · · · β1mΩ

∗− 12 1 Z∗T

1 [[u∗ 1 ◦ u∗

m ]]Z∗ mΩ

∗− 12 m

0 β2mΩ

∗− 12 2 Z∗T

2 [[u∗ 2 ◦ u∗

m ]]Z∗ mΩ

∗− 12 m ..

.

.. .

Symm.

..

. βm−1,mΩ

∗− 12 m−1Z∗T

m−1[[u∗ 1 ◦ u∗

2 ]]Z∗ mΩ

∗− 12 m 0









where

Ω^∗12 = diag{Ω^∗₁¹², . . . , Ω^∗m¹²}, Ω^∗_j = diagn

1

ζ_j2^∗ −λ^∗_j , . . . , _ζ^{∗ 1}

jN−λ^∗_j

o.

5 Numerical Results

Two-component BEC

(a) green: β^∗ = 1000, λ^∗₁ = λ^∗₂ = 7.07, E(u^∗) = 7.02

(b) red: β^∗ = 1000, λ^∗ = 10.34, λ^∗ = 14.54, E(u^∗) = 12.43

0 4.1 10.8 20 40 2

4 6 8 10 12 14 16

10⁶10⁷10⁸10⁹10¹⁰

β₁ β₂

(a) m = 2

β

Eigenvaluesλ

∗ 1(β),λ

∗ 2(β)

0 4.1 10.8 20 40

2 4 6 8 10 12 14

10⁶10⁷ 10⁸ 10⁹ 10¹⁰

m = 2

β₁ β₂

(b)

β EnergyE(u∗ (β))

Figure 5.1: (a): Eigenvalue curves, (b): energy curves, vs β.

Three-component BEC

(a) green: β^∗ = 1000, λ^∗₁ = λ^∗₂ = λ^∗₃ = 9.57, E(u^∗) = 9.52

(b) red: β^∗ = 1000, λ^∗₁ = λ^∗₃ = 18.36, λ^∗₂ = 20.85, E(u^∗) = 19.09

(c) blue: β^∗ = 1000, λ^∗₁ = 20.84, λ^∗₂ = 24.84, λ^∗₃ = 32.14, E(u^∗) = 25.85

0 3.9 10.8 20 29.4 40 0

5 10 15 20 25 30 35

10⁶10⁷ 10⁸ 10⁹ 10¹⁰

β₁ β₂

β₃ (a)

β

Eigenvalueλ

∗ j(β),j=1,2,3

0 3.9 10.8 20 29.4 40

0 5 10 15 20 25 30

10⁶10⁷ 10⁸ 10⁹ 10¹⁰

m = 3

β₁ β₂

β₃ (b)

β EnergyE(u∗ (β))

Figure 5.2: (a): Eigenvalue curves, (b): energy curves, vs β.

Verticillate Structures

• How to distribute in multi-component BEC when the scattering length is sufficiently large?

• All positive bound state solutions may repel each other and form finitely segregated nodal domains when scattering length approaches to infinity.

• Verticillate: [Botany] leaf, arranged in verticils.

We observe that verticillate or multiple verticillate structure (i) (n₁, . . . , n_γ) depends on m and

Pγ i=1

n_i = m (β ≫ 1), (Single, Double, Triple, Quadruple verticillate, ...) (ii) 1 ≤ n₁ ≤ 5.

Figure 5.3: Single, Double, Triple, Quadruple verticillate:

2 3 4 5 6 7 5

10 15 20 25

8 9 10 11 12 13 14

20 25 30 35 40 45 50

45 50 55 60 65 70 75 80

70 80 90 100 110 120 (2)

(3)

(4)

(5)

(1,5)

(2,7) (1,7)

(1,6)

(3,9) (2,9)

(2,8)

(4,10) (4,9)

(1,6,10)

(1,5,12) (1,5,13)

(1,7,12) (1,7,13)

(2,8,12)

(2,8,14)

(3,8,14) (3,9,14)

(4,9,14) (4,10,14) (1,3)

(1,4)

(6)

(7)

(8) (1,8)

(9)

(1,20)

(1,10) (12)

(10)

(1,19) (1,18)

(1,17) (1,16)

(1,15) (1,14)

(1,13)

(1,21) (1,22)

(1,23)

(1,24) (1,25)

(1,12)

(1,27) (1,11)

(1,26) (1,9)

(a)

(d) (c)

(b)

Energy

29 30 31 32 33 70

80 90 100 110 120 130

(1,28)

(1,29)

(4,10,15)

(5,10,15)

(5,10,16)

(1,5,11,15)

(1,5,11,16)

m-component BEC (e)

Energy

(i) Single verticillate (2) occur- ring at m = 2,

(ii) Double verticillate (1,5) oc- curring at m = 6,

(iii) Triple verticillate (1,6,10) occurring at m = 17,

(iv) Quadruple verticillate (1,5,11,15) occurring at m = 32.

(a) (b)

Figure 5.4: m = 5: (a) Ground state solutions, (b) bound state

(a) (b)

Figure 5.5: m = 6: (a) Ground state solutions, (b) bound state solutions.

(a) m = 14 (b) m = 17

(a) m = 29 (b) m = 32

Figure 5.7: Ground state solutions, (a) m = 29, (b) m = 32.

0 2.75 5.37 10 0

5 10 15 20 25 30 35

Λ

Energy

(1)

(1,8) (1,7)

(2,7)

(1,8) (2,7)

(9)

Λ₁(9)

Λ₂(9)

10⁵ 10^5.5 10⁶

6 Conclusion

• Theoretical

– The JI and GSI are proposed from the viewpoint of eigenvalue approach.

– The necessary and sufficient conditions of convergence of the GSI method are proven that the energy functional has a strictly local minimum at the fixed point.

• Numerical

– GSI method converges much faster than JI, globally and linearly between 10 to 20 steps.

– New phenomenon: verticillate multipling.

• Future works

7 Numerical Algorithms

Gauss-Seidel Type Iteration (GSI(m))

(i) Given A_j = A + 2[[V _j]] + 2β_jj[[u⁽⁰⁾ _j ² ]], β_jj ≪ 0, β_jk = β_kj ≥ 0 (j 6= k), j, k = 1, . . . , m and u⁽⁰⁾_j > 0 with ku⁽⁰⁾_j k₂ = 1,

j = 1, . . . , m. Let n = 0;

(ii) Repeat n: until convergence, (ii) For j = 1, . . . , m,

Use e.g., the Shift-Invert Arnoldi algorithm or the

Jacobi-Davidson algorithm to solve the minimal positive EW.

λ⁽ⁿ⁺¹⁾_j of A⁽ⁿ⁺¹⁾_j and the assoc. EV u⁽ⁿ⁺¹⁾_j with ku⁽ⁿ⁺¹⁾_j k₂ = 1, where

A⁽ⁿ⁺¹⁾_j := A_j + X

k<j

[[β_jku⁽ⁿ⁺¹⁾_j ]] + X

k≥j

[[β_jku⁽ⁿ⁾_j ]],

Comment: we denote u⁽ⁿ⁺¹⁾_j =f_j(u⁽ⁿ⁺¹⁾₁ , . . . , u⁽ⁿ⁺¹⁾_j−1 ,u⁽ⁿ⁾_j+1, . . . , u⁽ⁿ⁾_m );

(iv) Compute the residual,

res⁽ⁿ⁺¹⁾_j = A⁽ⁿ⁺¹⁾_j u⁽ⁿ⁺¹⁾_j − λ⁽ⁿ⁺¹⁾_j u⁽ⁿ⁺¹⁾_j , j = 1, . . . , m, (7.1) (v) If kres⁽ⁿ⁺¹⁾_j k₂ < Tol, j = 1, . . . , m, then stop, else n ← n + 1, go

to Repeat.

Variant GSI(2)≡ V1-GSI(2)

(i) Given A_j = A + 2[[V _j]] + 2α_j[[u⁽⁰⁾ _j ² ]], u⁽⁰⁾_j > 0 with ku⁽⁰⁾_j k₂ = 1, j = 1, 2, α_j ≪ 1, β > 0; Let n = 0;

(ii) Repeat n: until convergence,

(iii) Compute u⁽ⁿ⁺¹⁾₁ = f₁(u⁽ⁿ⁾₂ ), u⁽ⁿ⁺¹⁾₁ ← nl(ave(u⁽ⁿ⁺¹⁾₁ )), Compute u⁽ⁿ⁺¹⁾₂ = f₂(u⁽ⁿ⁺¹⁾₁ ),

(iv) Compute the residuals as in (7.1),

(v) If converges, then stop; else n ← n + 1, go to Repeat (ii).

Variant GSI(3)

(i) Given A_j := A + 2[[V _j]] + 2α_j[[u⁽⁰⁾ _j ² ]], u⁽⁰⁾_j > 0 with ku⁽⁰⁾_j k₂ = 1, j = 1, 2, 3, α_j ≪ 1, β > 0; Let n = 0;

(ii) Repeat n: until convergence, V1-GSI(3)

(iii) Compute u⁽ⁿ⁺¹⁾₁ = f₁(u⁽ⁿ⁾₂ , u⁽ⁿ⁾₃ ), u⁽ⁿ⁺¹⁾₁ ← nl(ave(u⁽ⁿ⁺¹⁾₁ )), Compute u⁽ⁿ⁺¹⁾₂ = f₂(u⁽ⁿ⁺¹⁾₁ , u⁽ⁿ⁾₃ ),

u⁽ⁿ⁺¹⁾₃ = f₃(u⁽ⁿ⁺¹⁾₁ , u⁽ⁿ⁺¹⁾₂ ),

(iv) Compute the residuals as in (7.1),

(v) If converges, then stop, else n ← n + 1, go to Repeat (ii);

V2-GSI(3)

(iii) Compute u⁽ⁿ⁺¹⁾₁ = f₁(u⁽ⁿ⁾₂ , u⁽ⁿ⁾₃ ), u⁽ⁿ⁺¹⁾₁ ← nl(ave(u⁽ⁿ⁺¹⁾₁ )), Compute u⁽ⁿ⁺¹⁾₂ = f₂(u⁽ⁿ⁺¹⁾₁ , u⁽ⁿ⁾₃ ), u⁽ⁿ⁺¹⁾₂ ← nl(ave(u⁽ⁿ⁺¹⁾₂ )), Compute u⁽ⁿ⁺¹⁾₃ = f₃(u⁽ⁿ⁺¹⁾₁ , u⁽ⁿ⁺¹⁾₂ ),

(iv) Compute the residuals as in (7.1),

(v) If converges, then stop; else n ← n + 1, go to Repeat (ii).

The energy curves E(u^∗(β)) in Figure 7.10(b) and 7.11(b) are computed by

E(u^∗(β)) = 1 2

Xm j=1

u^∗_j^⊤A_ju^∗_j + β 2

X

1≤j<k≤m

u^∗ _k ^{2 ⊤}u^∗_j ² .

Table 7.1: (g): ground states, (b): bound states.

m = 2 m = 3 green curves (g) GSI(2) GSI(3) red curves (b) V1-GSI(2) V1-GSI(3) blue curves (b) — V2-GSI(3)

(a) green: β^∗ = 1000, λ^∗₁ = λ^∗₂ = 7.07, E(u^∗) = 7.02

(b) red: β^∗ = 1000, λ^∗ = 10.34, λ^∗ = 14.54, E(u^∗) = 12.43

Table 7.2: Two-component BEC.

θ = π, m = 2 green red

(0, β₁) λ^∗₁ = λ^∗₂, u^∗₁ = u^∗₂ —

(β₁, β₂) λ^∗₁ = λ^∗₂, λ^∗₁ = λ^∗₂, u^∗₁ = u^∗₂ (β₂, β_∞) u^∗₂ = R_θ(u^∗₁) λ^∗₁ 6= λ^∗₂,

u^∗_j = ave(u^∗_j), j = 1, 2

0 4.1 10.8 20 40 2

4 6 8 10 12 14 16

10⁶10⁷10⁸10⁹10¹⁰

β₁ β₂

(a) m = 2

β

Eigenvaluesλ

∗ 1(β),λ

∗ 2(β)

0 4.1 10.8 20 40

2 4 6 8 10 12 14

10⁶10⁷ 10⁸ 10⁹ 10¹⁰

m = 2

β₁ β₂

(b)

β EnergyE(u∗ (β))

Figure 7.10: (a): Eigenvalue curves, (b): energy curves, vs β.

(a) green: β^∗ = 1000, λ^∗₁ = λ^∗₂ = λ^∗₃ = 9.57, E(u^∗) = 9.52

(b) red: β^∗ = 1000, λ^∗₁ = λ^∗₃ = 18.36, λ^∗₂ = 20.85, E(u^∗) = 19.09

(c) blue: β^∗ = 1000, λ^∗₁ = 20.84, λ^∗₂ = 24.84, λ^∗₃ = 32.14, E(u^∗) = 25.85

Table 7.3: Three-component BEC.

θ = ^2π₃ , m = 3 green red blue

(0, β1) λ^∗₁ = λ^∗₂ = λ^∗₃,

u^∗₁ = u^∗₂ = u^∗₃ — —

(β1, β2) —

(β2, β3)

λ^∗₁ = λ^∗₂ = λ^∗₃, u^∗₂ = Rθ(u^∗₁), u^∗₃ = Rθ(u^∗₂)

λ^∗₁ = λ^∗₂ 6= λ^∗₃, u^∗₁ = u^∗₂,

{u^∗_j = ave(u^∗_j)}³_j=1

(β3, β^∞)

λ^∗₁ 6= λ^∗₂ = λ^∗₃, u^∗₁ = Rπ(u^∗₁), u^∗₃ = Rπ(u^∗₂)

λ^∗₁ 6= λ^∗₂ 6= λ^∗₃, {u^∗_j = ave(u^∗_j)}³_j=1

0 3.9 10.8 20 29.4 40 0

5 10 15 20 25 30 35

10⁶10⁷ 10⁸ 10⁹ 10¹⁰

β₁ β₂

β₃ (a)

β

Eigenvalueλ

∗ j(β),j=1,2,3

0 3.9 10.8 20 29.4 40

0 5 10 15 20 25 30

10⁶10⁷ 10⁸ 10⁹ 10¹⁰

m = 3

β₁ β₂

β₃ (b)

β EnergyE(u∗ (β))

Figure 7.11: (a): Eigenvalue curves, (b): energy curves, vs β.