Transcript Numerical Simulations for Bose

Mathematical Analysis and Numerical
Simulation for Bose-Einstein Condensation
Weizhu Bao
Department of Mathematics
& Center of Computational Science and Engineering
National University of Singapore
Email: [email protected]
URL: http://www.math.nus.edu.sg/~bao
Collaborators
External
–
–
–
–
–
–
–
–
P.A. Markowich, Institute of Mathematics, University of Vienna, Austria
D. Jaksch, Department of Physics, Oxford University, UK
Q. Du, Department of Mathematics, Penn State University, USA
J. Shen, Department of Mathematics, Purdue University, USA
L. Pareschi, Department of Mathematics, University of Ferarra, Italy
I-Liang Chern, Department of Mathematics, National Taiwan University, Taiwan
C. Schmeiser & R.M. Weishaeupl, University of Vienna, Austria
W. Tang & L. Fu, Beijing Institute of Appl. Phys. & Comput. Math., China
Internal
– Yanzhi Zhang, Hanquan Wang, Fong Ying Lim, Ming Huang Chai
– Yunyi Ge, Fangfang Sun, etc.
Outline
Part I: Predication & Mathematical modeling
–
–
–
–
Theoretical predication
Physical experiments and results
Applications
Gross-Pitaevskii equation
Part II: Analysis & Computation for Ground states
–
–
–
–
Existence & uniqueness
Energy asymptotics & asymptotic approximation
Numerical methods
Numerical results
Outline
Part III: Analysis & Computation for Dynamics in BEC
– Dynamical laws
– Numerical methods
– Vortex stability & interaction
Part IV: Rotating BEC & multi-component BEC
–
–
–
–
–
BEC in a rotational frame
Two-component BEC
Spinor BEC
BEC at finite temperature
Conclusions & Future challenges
Part I
Predication
&
Mathematical modeling
Theoretical predication
Particles be divided into two big classes
–
Bosons: photons, phonons, etc
•
•
•
–
Integer spin
Like in same state & many can occupy one obit
Sociable & gregarious
Fermions: electrons, neutrons, protons etc
•
•
Half-integer spin & each occupies a single obit
Loners due to Pauli exclusion principle
Theoretical predication
For atoms, e.g. bosons
–
Get colder:
•
–
Very cold:
•
–
Behave more like waves & less like particles
Overlap with their neighbors
Extremely cold:
•
•
•
Most atoms behavior in the same way, i.e gregarious
quantum mechanical ground state,
`super-atom’ & new matter of wave & fifth state
Theoretical predication
S.N. Bose: Z. Phys. 26 (1924)
–
–
–
Study black body radiation: object very hot
Two photons be counted up as either identical or different
Bose statistics or Bose-Einstein statistics
A. Einstein: Sitz. Ber. Kgl. Preuss. Adad. Wiss. 22 (1924)
–
–
Apply the rules to atoms in cold temperatures
Obtain Bose-Einstein distribution in a gas
f ( i ) 
1
e  i  / kBT  1
Experimental results
JILA (95’, Rb, 5,000):
Science 269 (1995)
–Anderson et al.,
Science, 269 (1995),
198: JILA Group; Rb
–Davis et al., Phys.
Rev. Lett., 75 (1995),
3969: MIT Group; Rb
–Bradly et al., Phys.
Rev. Lett., 75 (1995),
1687, Rice Group; Li
Experimental results
Experimental implementation
–
–
–
–
–
JILA (95’): First experimental realization of BEC in a gas
NIST (98’): Improved experiments
MIT, ENS, Rice,
ETH, Oxford,
Peking U., …
2001 Nobel prize in physics:
– C. Wiemann: U. Colorado
– E. Cornell: NIST
– W. Ketterle: MIT
ETH (02’, Rb, 300,000)
Experimental difficulties
Low temperatures  absolutely zero (nK)
Low density in a gas
Experimental techniques
Laser cooling
Magnetic trapping
Evaporative Cooling
($100k—300k)
Possible applications
Quantized vortex for studying superfluidity
Square Vortex
lattices in spinor
BECs
Giant
vortices
Vortex
lattice
dynamics
Test quantum mechanics theory
Bright atom laser: multi-component
Quantum computing
Atom tunneling in optical lattice trapping, …..
Mathematical modeling
N-body problem
– (3N+1)-dim linear Schroedinger equation
Mean field theory: T  Tc  O(nK)
– Gross-Pitaevskii equation (GPE):
– (3+1)-dim nonlinear Schroedinger equation (NLSE)
Quantum kinetic theory
– High temperature: QBME (3+3+1)-dim
– Around critical temperature: QBME+GPE
– Below critical temperature: GPE
Gross-Pitaevskii equation (GPE)
Physical assumptions
– At zero temperature
– N atoms at the same hyperfine species (Hartree ansatz)


Ψ (x, x
N
1
2

,, x

, t)   ψ( x , t)
N
N
i
i 1
– The density of the trapped gas is small
| as |
a
 1,
3
 |as|
 1 .
– Interatomic interaction is two-body elastic and in Fermi form
Second Quantization model
The second quantized Hamiltonian:
– A gas of bosons are condensed into the same single-particle state
– Interacting by binary collisions
– Contained by an external trapping potential
1
H (  )    ( x ') H 0 ( x ')dx '   † ( x ')  † ( x ) Vint ( x ', x ) ( x ') ( x )dx ' dx
2
†
 ( x ) :  ( x , t ) : Bose field operator
P2
H0 
 Vext ( x ') : single particle Hamiltonian
2m
P  i   ( px , p y , pz )T : moentum operator
 † ( x ) &  ( x ) : creation & annihilation of a particle at position x
Second quantization model
– Crucial Bose commutation rules:
 ( x '),  ( x )   ( x ' x ),  ( x '), ( x )   † ( x '), † ( x )   0

 



†
– Atomic interactions are low-energy two-body s-wave
collisions, i.e. essentially elastic & hard-sphere collisions
Vint ( x ', x)  U0  ( x  x ') with U0  4 2as / m
– The second quantized Hamiltonian
U0
†
†
H (  )    ( x ') H 0 ( x ')dx ' 

(
x
')

( x ') ( x ') ( x ')dx '

2
†
Second quantization model
The Heisenberg equation for motion:

i
 ( x )    ( x ), H ( )  :  H 0  U 0  † ( x ) ( x )   ( x )
t
For a single-particle state with macroscopic occupation
( x )  N  ( x, t )   ( x, t )
 ( x , t ) : marcoscopic wave function (=  ( x ) / N : expectation value of  ( x ))
 ( x , t ) : fluctation operator satisfies  ( x , t )  0
– Plugging, taking only the leading order term
– neglecting the fluctuation terms (i.e., thermal and quantum depletion of the
condensate)
– Valid only when the condensate is weakly-interacting & low tempertures
Gross-Pitaevskii equation
The Schrodinger equation (Gross, Nuovo. Cimento., 61; Pitaevskii, JETP,61 )
 ψ(x,t)
 H ( )
i

,
x  ( x, y, z )
t
  *( x, t )
– The Hamiltonian:

2
H ( )    *( x , t ) [
2m
2
 V ( x , t )] ( x , t ) dx
1
  *( x, t )  *( x , t )( x  x)  ( x, t )  ( x , t ) dx dx
2
– The interaction potential is taken as in Fermi form

( x)  ( N 1)
U
0

 ( x),
U 0  4 
2
a
s
m
.
Gross-Pitaevskii equation
The 3d Gross-Pitaevskii equation (

x  ( x, y, z ) )
2




 2 

2

i   ( x, t )  
  ( x, t )  V ( x )  ( x, t )  N U 0 |  ( x, t ) |  ( x, t )
t
2m
– V is a harmonic trap potential

V (x) 
m
(
2

2
x
x2  y y 2  z z 2 )
2
2
– Normalization condition
 2 
 3 |  ( x, t ) | dx  1.
R
Gross-Pitaevskii equation
Scaling (w.l.o.g.
x  y  z
)
– Dimensionless variables
x
t   t,
x ,
 ( x , t )  a  ( x , t ),
a
a
– Dimensionless Gross-Pitaevskii equation
3/ 2
x
0
0

0
m x



 2 

1 2 
i  ( x, t )     ( x, t )  V ( x ) ( x, t )   |  ( x, t ) |  ( x, t )
t
2
– With
2
2
 1 2
2
V ( x )  ( x   y   z 2 ),
y
z
2

y



y
x
,
 

z
z
x
,

4 N as
a0
Gross-Pitaevskii equation
Typical parameters (
  1.051034 [Js]
)
– Rb Used in JILA
87
 
m  1.44 1025 [kg],
as  5.1 [nm],
–
23
s

a0 
m x
y
 10 2 [1 / s],

z
 8x
 0.3407105 [m],  
4 N a s
a
 0.01881N
0
NaUsed in MIT
 
m  3.8 1026 [kg],
a
x
 2.75 [nm],
x
a
0


m z
y
 360 2 [1 / s],

z
 3.5  2 [1 / s]
 1.120910 [m],  
5
4 N as
a
0
 0.003083N
Gross-Pitaevskii equation
Reduction to 2d (disk-shaped condensation)
– Experimental setup x  y , z  x   y  1,  z  1
– Assumption: No excitations along z-axis due to large energy
 ( x, y, z , t )   12 ( x, y, t ) 3 ( z ) with
3 ( z )  (  |g ( x, y, z ) |2 dxdy)1/ 2
R
1
 
 ho ( z )   z 
 
1/ 4
e  z z
2
/2

2d Gross-Pitaevskii equation ( x  ( x, y),   12 )
x 2   y2 y 2

1
i  ( x , t )    
   2 | |2  ,
t
2
2


 2    ( z )dz   ho4 ( z )dz  
4
3


z
:  2a
2
Numerical Verification
Numerical Results
Bao, Y. Ge, P. Markowich & R. Weishaupl, 06’
Gross-Pitaevskii equation

General form of GPE ( x  R d )



 2 

1 2 
i  ( x, t )     ( x, t )  Vd ( x ) ( x, t )   d |  ( x, t ) |  ( x, t )
t
2
with

 y z
4
  2  ( y, z ) dydz  
,
23
R
2



4

z
 d      3 ( z ) dz  
,
2



,



Normalization condition

Rd
1 2
2 x ,
2
  1 2
Vd ( x )   ( x   y 2 ),
y
2
 1 ( x 2   2 y 2   2 z 2 ),
y
z
 2
 2 
|  ( x , t ) | dx  1.
d 1
d 2
d 3
Gross-Pitaevskii equation
Two kinds of interaction
– Repulsive (defocusing) interaction
 0
a
s

d  0
– Attractive (focusing) interaction
a  0   
Four typical interaction regimes:
s
d
0
– Linear regime: one atom in the condensation
d  0
– Weakly interacting condensation
|  d |  1
Gross-Pitaevskii equation
– Strongly repulsive interacting condensation
d
 1
– Strongly attractive interaction in 1D
1  0 & | 1 |  1
Other potentials
–
–
–
–
Box potential
Double-well potential
Optical lattice potential
On a ring or torus
Gross-Pitaevskii equation
Conserved quantities
– Normalization of the wave function
N ( (t )) 
– Energy



| ( x , t ) | 2 dx  N ( (0))  1
Rd
d
1
2
2
4
[
|


(
x
,
t
)
|

V
(
x
)|

(
x
,
t
)
|

|

(
x
,
t
)
|
] dx
d
Rd
2
2
 E  ( (0))
E  ( (t ))  
Chemical potential


( (t ))  
Rd
 2

 2
 4 
1
[ |  ( x , t ) | Vd ( x ) |  ( x , t ) |   d | ( x , t ) | ] dx
2
Semiclassical scaling
When
, re-scaling x  x 
d  1
1/ 2
    d / 4   1/ d2/( d 2)




  
2 2  
i   ( x, t )     ( x, t )  Vd ( x )  ( x, t )  |  ( x, t ) |2   ( x, t )
t
2
With
E  (  )  
Rd
[
2


1
|   |2 Vd ( x ) |   |2  |  |4 ] dx  O(1)
2
2
Leading asymptotics (Bao & Y. Zhang, Math. Mod. Meth. Appl. Sci., 05’)
E ( )   1E (  )  O  1   O  d2/( d 2) 
 ( )   1  (  )  O  1   O  d2/( d 2) 
Comparison of two scaling
Quanties
ts
xs
Thomas-Fermi scaling
1/  x
a0 
/ m x
Semiclassical scaling
1/ x
a0 1/ 2  a0  d1/( d  2)  a0 (
s
a03/ 2
a03/ 2 d / 4
Es
mxs2
 x : 2
ts
 x 1
Energy E
O   d2 /( d  2) 
O (1)
Chemical potential 
O   d2 /( d  2) 
O(1)
length of wave function
O   d1/( d  2)   O( 2  )
O(1)
height of wave function O   d d / 2( d  2)   O(  /  d )
O(1)
4 Nas 1/( d  2)
)
a0
Quantum Hydrodynamics
Set


   e
iS  / 


J    v  ,   1 /  d2 /(d  2 )

, v   S  ,
Geometrical Optics: (Transport + Hamilton-Jacobi)
 t      (   S  )  0,
1
2
 2

 t S  S  Vd ( x )   
2
2

1

 
Quantum Hydrodynamics (QHD): (Euler +3rd dispersion)
t      (   v  )  0

t ( J )    (
J  J




)  P(  )   Vd 
P(  )   2 / 2
2
4
(    ln   )
Part II
Analysis & Computation
for
Ground states
Stationary states
Stationary solutions of GPE


i  t
 ( x, t )  e  ( x )
Nonlinear eigenvalue problem with a constraint

 
 2 
1 2 
  ( x )     ( x )  Vd ( x ) ( x )   d |  ( x ) |  ( x ),
2
 2 
 d |  (x) | dx  1

x Rd
R
Relation between eigenvalue and eigenfunction
    ( )  E ( ) 
d
2

Rd
 4 
|  (x) | dx
Stationary states
Equivalent statements:
– Critical points of E ( ) over the unit sphere S   |  L  1, E ( )  
– Eigenfunctions of the nonlinear eigenvalue problem
– Steady states of the normalized gradient flow:(Bao & Q. Du, SIAM J. Sci. Compu., 03’)
 ( )
1 2
2
t ( x, t )  [   V ( x )   |  | 
] ,
2
||  ||2
 ( x, 0)  0 ( x )
with
|| 0 ( x ) || 1.
2
Minimizer/saddle points over the unit sphere :
– For linear case d  0
(Bao & Y. Zhang, Math. Mod. Meth. Appl. Sci., 05’)
• Global minimizer vs saddle points
– For nonlinear case d  0
• Global minimizer, local minimizer vs saddle points
Ground state
Ground state:
E (g )  min E ( ),  g    (g )  E (g ) 
|| || 1
d
2

R
4
|

(
x
)
|
dx
g
d
Existence and uniqueness of positive solution :
d  0
– Lieb et. al., Phys. Rev. A, 00’
Uniqueness up to a unit factor
g  g ei
0
with any constant 0
Boundary layer width & matched asymptotic expansion
– Bao, F. Lim & Y. Zhang, Bull. Institute of Math., Acad. Scinica , 07’
Excited & central vortex states
Excited states: 1 , 2 , 3 , 
i 

(
x
,
y
,
t
)

Central vortex states:
e
m
t
1 d  dm (r )   m2 r 2 
2
m m (r )  
r
   2   m (r )   2 | m | m ,
2r dr 
dr   2r
2
2

 m (r ) r dr  1,
2
i m t
m ( x, y )  e
m ( r ) e
0r
m (0)  0
0
Central vortex line states in 3D:
Open question: (Bao & W. Tang, JCP, 03’; Bao, F. Lim & Y. Zhang, TTSP, 06’)
g ,
1 ,
2 ,

E ( g )  E (1 )  E ( 2 )  
  ( g )    (1 )    ( 2 )   ???????
im 
Approximate ground states
Three interacting regimes
– No interaction, i.e. linear case
– Weakly interacting regime
– Strongly repulsive interacting regime
Three different potential
– Box potential
– Harmonic oscillator potential
– BEC on a ring or torus
Energies revisited
Total energy:
d
1
2
2
[
|


(
x
)
|

V
(
x
)|

(
x
)
|

| ( x ) |4 ] dx : E kin( )  E pot( )  E int( )
d
d
R
2
2
1
2
(

)

|


(
x
)
|
dx
Kinetic energy: E kin
d
2 R
E  ( )  
–
– Potential energy: E ( )   V ( x)|  ( x) |
– Interaction energy: E ( )    | ( x ) |
2
pot
R
d
int
Chemical potential
1
  ( )   [ 2 |  ( x ) |
2
Rd
2
dx
4
dx
d
d
R
d
Vd ( x )|  ( x ) |2   d | ( x ) |4 ] dx
=E ( )  Eint ( )  Ekin ( )  Epot ( )  2 Eint ( )
Box Potential in 1D
 0, 0  x  1,
V ( x)  
, otherwise.
The potential:
The nonlinear eigenvalue problem
1
  ( x)    ( x)   |  ( x) |2  ( x),
2
 (0)   (1)  0
0  x  1,
1
2
|

(
x
)
|
dx  1

with
Case I: no interaction, i.e.   0
0
– A complete set of orthonormal eigenfunctions
l ( x)  2 sin(l x),
1
2
l  l 2 2 ,
l  1, 2,3,
Box Potential in 1D
– Ground state & its energy:
g ( x)   ( x)  2 sin( x),
0
g
Eg : E0 ( ) 
0
g
2
2
– j-th-excited state & its energy
 j ( x)   ( x)  2 sin(( j  1) x),
0
j
  g : 0 (g0 )
( j  1)2  2
E j : E0 ( ) 
  j : 0 ( j0 )
2
0
j
Case II: weakly interacting regime, i.e. |  | o(1)
– Ground state & its energy:
3
2
0
g ( x)   ( x)  2 sin( x), Eg : E (g )  E ( ) 

,  g :  (g )   (g ) 
 3
2
2
2
0
g
0
g
2
– j-th-excited state & its energy
( j  1) 2  2 3
 j ( x)   ( x)  2 sin(( j  1) x), E j : E ( j )  E ( ) 

,
2
2
( j  1) 2  2
0
 j :  ( j )    ( j ) 
 3
2
0
j
0
j
Box Potential in 1D
Case III: Strongly interacting regime, i.e.

1
– Thomas-Fermi approximation, i.e. drop the diffusion term
 gTF gTF ( x)   | gTF ( x) |2 gTF ( x),
0  x  1,
TF

g
 gTF ( x) 

1

TF
2
|

(
x
)
|
dx  1
 g
0
g (x)  gTF ( x)  1,
E g  E TF
g 

2
,
g   gTF   ,
• Boundary condition is NOT satisfied, i.e.
• Boundary layer near the boundary
gTF (0)  gTF (1)  1  0
Box Potential in 1D
– Matched asymptotic approximation
• Consider near x=0, rescale
• We get
x
1
g
1
( X )   ( X )   3 ( X ), 0  X  ;
2
• The inner solution
( X )  tanh( X ), 0  X  
X,
 ( x) 
g
 ( x)

(0)  0,
lim ( X )  1
X 
g
 g ( x ) 
tanh(  g x), 0  x  o(1)

• Matched asymptotic approximation for ground state
MA

g
g ( x)  gMA ( x) 

 tanh(  MA x)  tanh(  MA (1  x))  tanh(  MA )  , 0  x  1
g
g
g


1
1   | gMA ( x) |2 dx   g   gMA    2   1  2   gTF  2   1  2, 
0
1.
Box Potential in 1D
• Approximate energy
4
 2
MA
  1  2, Eint, g  Eint,


  1,
g
2 3
2 3
2
MA
 Ekin, g 
 1  2
3
Eg  EgMA 
Ekin, g


• Asymptotic ratios:
Eg
1
lim
 ,
  
2
g
lim
 
Eint, g
Eg
 1,
lim
 
Ekin, g
Eg
 0,
• Width of the boundary layer:
O(1/  )
Numerical observations:
g  gMA

L
 O(e3
 /2
Eg  EgMA  O(1/  ),
),
g  gMA
2
L
 O(e3
 /2
MA
Ekin, g  Ekin,
g  O(1/  ),
),
g  gMA  O(e3
 /2
MA
Eint, g  Eint,
g  O(1/  )
)
Box Potential in 1D
• Matched asymptotic approximation for excited states
MA

j
 j ( x)   jMA ( x) 

[ j / 2]
[( j 1) / 2]
[

tanh( 
l 0
 tanh(

l 0
MA
g
MA
g
2l
(x 
)) 
j 1
2l  1
(
 x))  C j tanh(  gMA )]
j 1
• Approximate chemical potential & energy
 j   MA
   2( j  1)   ( j  1)2  2( j  1)2 ,
j

4
 ( j  1)   ( j  1) 2 ,
2 3
 2
MA
Eint, j  Eint,
 ( j  1)   ( j  1) 2 ,
j 
2 3
2
MA
Ekin, j  Ekin,
( j  1)   ( j  1) 2  2( j  1) 2
j 
3
E j  E MA

j
Fifth excited states
Energy & Chemical potential
Box potential in 1D
• Boundary layers & interior layers with width
O(1/  )
• Observations: energy & chemical potential are in the same order
E (g )  E (1 )  E (2 ) 
  (g )   (1 )   (2 ) 
• Asymptotic ratios:
Ej
1
lim
 ,
  
2
j
lim
 
Ej
Eg
 1,
lim
 
Eint, j
Ej
 1,
j
 1,
  
g
lim
• Extension to high dimensions
lim
 
lim
 
Ej
g
Ekin, j
Ej
 0,
 0,
Harmonic Oscillator Potential in 1D
The potential:
The nonlinear eigenvalue problem
x2
V ( x) 
2
1
  ( x)    ( x)  V ( x) ( x)   |  ( x) |2  ( x), with
2
Case I: no interaction, i.e.   0

2
|

(
x
)
|
dx  1


– A complete set of orthonormal eigenfunctions
l ( x)  (2l l !) 1/ 2
1
 1/ 4
e x
2
/2
H l ( x),
l 
l 1
,
2
l x
l x d e
H l ( x)  (1) e
: Hermite polynomials with
l
dx
H 0 ( x)  1, H1 ( x)  2 x, H 2 ( x)  4 x 2  2,
2
2
l  0,1, 2,3,
Harmonic Oscillator Potential in 1D
– Ground state & its energy:
g ( x)  g0 ( x) 
1
 1/ 4
e x
2
/2
Eg : E0 (g0 ) 
,
– j-th-excited state & its energy
 j ( x)   j0 ( x)  (2 j j !) 1/ 2
1
 1/ 4
e x
2
/2
H j ( x),
1
  g : 0 (g0 )
2
E j : E0 ( 0j ) 
( j  1)
  j : 0 ( j0 )
2
Case II: weakly interacting regime, i.e. |  | o(1)
– Ground state & its energy:
g ( x)  g0 ( x) 
1
 1/ 4
e x
2
/2
, Eg : E (g )  E (g0 ) 
1 
1
 C0 ,  g :   (g )    (g0 )    C0
2 2
2
– j-th-excited state & its energy
 j ( x)   j0 ( x), E j : E ( j )  E ( j0 ) 
( j  1)
 j :  ( j )   ( ) 
 Cj
2
0
j
( j  1) 
 Cj,
2
2

with C j =  | j0 ( x) |4 dx
-
Harmonic Oscillator Potential in 1D
Case III: Strongly interacting regime, i.e.

1
– Thomas-Fermi approximation, i.e. drop the diffusion term

TF
g
 (  TF  x 2 / 2) /  , | x | 2  TF
g
g
 ( x )  V ( x )  ( x )   |  ( x ) |  ( x)   ( x)  
0,
otherwise


2(2  gTF )3/ 2
1 3
TF
2
1   | g ( x) | dx 

 g   gTF  ( ) 2 / 3
3
2 2
-
TF
g
TF
g
TF
g
2
TF
g
TF
g
1/3
O
(

)
– Characteristic length:
TF
x


2

g
– It is NOT differentiable at
– The energy is infinite by direct definition:
E (gTF )  , Ekin (gTF )  
Harmonic Oscillator Potential in 1D
– A new way to define the energy
1 3 2 / 3
TF
TF
Eint, g  Eint,

E
(

)

( ) ,
g
int
g
5 2
1 3 2 / 3
TF
TF
Epot, g  Epot,

E
(

)

( )
g
pot
g
10 2
3 3 2 / 3
TF
TF
Eg   g  Eint, g   gTF  Eint,

(
)
:

E
g
g ,
10 2
TF
TF
Ekin, g  EgTF  Eint,
g  Epot, g  0
– Asymptotic ratios
Eg
3
lim
 ,
  
5
g
lim
 
Eint, g
Eg
2
 ,
3
lim
 
Epot, g
Eg
1
 ,
3
lim
 
Ekin, g
Eg
 0,
Numerical observations:
g  gTF
 O(
ln 
),
2/5

ln 
Eg  EgMA  O( 2 / 3 ),

L
g  gMA
L2
 O(
MA
Ekin, g  Ekin,
g  O(
ln 

),
2/5
ln 

),
2/3
 g   gMA  O(
ln 

MA
Eint, g  Eint,
g  O(
2/3
)
ln 

2/3
)
Harmonic Oscillator Potential in 1D
– Thomas-Fermi approximation for first excited state
1TF 1TF ( x)  V ( x) 1TF ( x)   | 1TF ( x) |2 1TF ( x)
sign( x) (  TF  x 2 / 2) /  , 0 | x | 2  TF
1
1
 ( x)  
0,
otherwise

TF
1

2(21TF )3/ 2
1   |  ( x) | dx 
3
-
TF
1
2
• Jump at x=0!
• Interior layer at x=0

1 3 2 / 3
)
2 2
1  1TF  (
Harmonic Oscillator Potential in 1D
– Matched asymptotic approximation

MA
1

1MA
sign( x) ( 1MA  x 2 / 2) /  , 0 | x | 21MA
MA
( x) 
[tanh( 1 x)  sign( x)]  

0,
otherwise


– Width of interior layer:
O(1/ 1MA )  O(1/  1/ 3 )  1MA  O(  2 / 3 )
– Ordering:
E (g )  E (1 )   (g )   (1 )
Harmonic Oscillator Potential
Extension to high dimensions
Identity of energies for stationary states in d-dim.
2 Ekin  2 Epot  d Eint  0
– Scaling transformation
 ( x)  (1  )d / 2 0 ((1  ) x)
with  0 ( x): a stationary state
– Energy variation vanishes at first order in 
E ( ( x ))  (1  )2 Ekin ( 0 )  (1  )2 Epot ( 0 )  (1  ) d Eint ( 0 )
d
E ( ( x )) | 0  0
d
BEC on a ring
The potential: V ( x)  0 on an interal
The nonlinear eigenvalue problem
1
  ( )    ( )   |  ( ) |2  ( ),
2
 (  2 )   ( ),
0    2 ,
2
2
|

(

)
|
d  1

with
For linear case, i.e.   0
0
– A complete
set of orthonormal
eigenfunctions
1
1
1
0 ( ) 
2
, 2l ( ) 
l2
0  0, 2l  2l 1  ,
2

cos(l ), 2l 1 ( ) 
l  1, 2,3,

sin(l );
BEC on a ring
– Ground state & its energy:
g ( x)  g0 ( x) 
1
,
2
Eg : E0 (g0 )  0   g : 0 (g0 )
– j-th-excited state & its energy
 j ( x )   ( x) 
0
j
Some properties
1

cos(l ),
j2
E j : E0 ( ) 
  j : 0 ( j0 )
2
|  | o(1)
– Ground state & its energy
g ( x)  g0 ( x) 
1
,
2
0
j
Eg : E (g0 ) 


,  g :  (g0 ) 
8
4
– With a shift:
 ( ) is a solution   (  0 ) is also a solution
– Interior layer can be happened at any point in excited states
Numerical methods for ground states
Runge-Kutta method: (M. Edwards and K. Burnett, Phys. Rev. A, 95’)
Analytical expansion: (R. Dodd, J. Res. Natl. Inst. Stan., 96’)
Explicit imaginary time method: (S. Succi, M.P. Tosi et. al., PRE, 00’)
Minimizing E ( ) by FEM: (Bao & W. Tang, JCP, 02’)
Normalized gradient flow: (Bao & Q. Du, SIAM Sci. Comput., 03’)
– Backward-Euler + finite difference (BEFD)
– Time-splitting spectral method (TSSP)
Gauss-Seidel iteration method: (W.W. Lin et al., JCP, 05’)
Spectral method + stabilization: (Bao, I. Chern & F. Lim, JCP, 06’)
Imaginary time method
Idea: Steepest decent method + Projection


1  E ( ) 1 2
t ( x , t )  
    V ( x )   |  |2  , tn  t  tn 1
2 
2

 

 ( x , t n1)

ˆ

 ( x , tn1 ) 
,
n  0,1,2,
 
||  ( x , t n 1) ||



 ( x ,0)  0 (x) with || 0 ( x ) || 1.
0
1
1
2
E (ˆ1 )  E (0 )
E (ˆ )  E ( )
1
1
E (1 )  E (0 ) ??
g
Physical institutive in linear case

– Solution of GPE:  ( x, t )   a e  ( x)
– Imaginary time dynamics:   i t
j 0

 ( x , )   ( x , t )   a j e
j 0
 j 
 j ( x)
0  1 

i  j t
j
j
 ( x , )
a0 e
 0 

with  ( x,0)   0 ( x )   a j  j ( x )
j 0
 
 0 ( x ) : grond state
Mathematical justification
For gradient flow
For linear case:
(Bao & Q. Du, SIAM Sci. Comput., 03’)
(Bao & Q. Du, SIAM Sci. Comput., 03’)
E0 (  (., tn1 ) )  E0 (  (., tn ) )    E0 (  (., 0) )
For nonlinear case: ???
Mathematical justification
Normalized gradient flow
Idea: (Bao & Q. Du, SIAM Sci. Comput., 03’)
– The projection step is equivalent to solve an ODE
t ( x, t )   (t , tn ) ( x, t ), tn  t  tn1 with  (t , k )  
1
ln  ( x, tn1 )
2tn
2
&  ( x, tn )   ( x, tn1 )
– Gradient flow with discontinuous coefficients:
1
2
t ( x , t )   2  V ( x )    |  |2    (t , tn )  , t  0,
– Letting time step go to 0
 ( (., t ))
1 2
2
t ( x, t )     V ( x )    |  |  
 , t  0,
2
2
||  (., t ) ||
 ( x,0)  0 ( x )
with
|| 0 ( x ) || 1.
– Mass conservation & Energy diminishing
||  (.,t ) |||| 0 || 1,
d
E ( (.,t ))  0,
dt
t 0
Fully discretization
Consider in 1D:
1
2
t ( x, t )  xx  V ( x)   |  |2  , x    (a, b), tn  t  tn 1 ,  (a, t )   (b, t )  0

 ( x, tn 1 ) 
 ( x, t n 1)

||  ( x, t n 1) ||
,
 ( x, 0)  0 (x)
with
|| 0 ( x) || 1.
Different Numerical Discretizations
– Physics literatures: Crank-Nicolson FD or Forward Euler FD
– BEFD: Energy diminishing & monotone
(Bao & Q. Du, SIAM Sci. Comput., 03’)
– TSSP: Spectral accurate with splitting error (Bao & Q. Du, SIAM Sci. Comput., 03’)
– BESP: Spectral accuracy in space & stable (Bao, I. Chern & F. Lim, JCP, 06’)
– Crank-Ncolson FD for normalized gradient flow
Backward Euler Finite Difference
Mesh and time steps:
x j  a  j h,
j  0,1,
h  x 
ba
;
M
k  t  0;
, M ; tn  n k , k=0,1,2, ;
 jn  (x j ,t n )
BEFD discretization
2nd order in space; unconditional stable; at each step,
only a linear system with sparse matrix to be solved!
Backward Euler Spectral method
Discretization
Spectral order in space; efficient & accurate
Ground states
Numerical results (Bao&W. Tang, JCP, 03’; Bao, F. Lim & Y. Zhang, TTSP, 06’)
– In 1d
• Box potential: V ( x)  0 0  x  1;  otherwise
– Ground state; excited states: first fifth
• Harmonic oscillator potential:
V(x)  x 2 / 2
– ground & first excited & Energy and chemical potential
• Double well potential :
V ( x)  (4  x 2 )2 / 2
– Ground & asymptotic excited state
• Optical lattice potential:
V ( x)  x 2 / 2  12sin 2 (4x)
– Ground & asymptotic excited state
with potential
next
back
back
back
back
back
1
0
E ( g ) E (1 ) E ( 2 ) E (3 )   ( g )   (1 )   ( 2 )   (3 )
0.5000 1.5000 2.5000 3.500 0.5000 1.5000 2.5000 3.500
3.1371 1.0441 1.9414 2.8865 3.8505 1.5266 2.3578 3.2590 4.1919
31.371 3.9810 4.7438 5.5573 6.4043 6.5527 7.2802 8.0432 8.8349
156.855 11.464 12.191 12.944 13.719 19.070 19.784 20.512 21.252
313.71 18.171 18.891 19.629 20.383 30.259 30.971 31.691 32.419
Observations
E ( g )  E (1 )  E ( 2 )  ,   ( g )    (1 )    ( 2 )  ,
for any fixed
lim
1  
back
E ( j )
E ( g )
1
lim
1  
  ( j )
1
  ( g )
1
back
back
back
back
back
Ground states
Numerical results (Bao&W. Tang, JCP, 03’; Bao, F. Lim & Y. Zhang, BIM, 07’)
– In 2d
• Harmonic oscillator potentials:
– ground
• Optical lattice potential:
– Ground & asymptotic excited states
– In 3D
• Optical lattice potential:
ground
asymptotic excited states
next
back
back
back
back
Part III
Analysis & Computation
for
Dynamics in BEC
Dynamics of BEC
Time-dependent Gross-Pitaevskii equation



 2 

1 2 
i  ( x , t )     ( x , t )  Vd ( x ) ( x , t )   d | ( x , t ) |  ( x , t )
t
2


 ( x ,0)   0 ( x )
Dynamical laws
–
–
–
–
–
Time reversible & time transverse invariant
Mass & energy conservation
Angular momentum expectation
Condensate width
Dynamics of a stationary state with its center shifted
Angular momentum expectation
Definition:


Lz (t ) :  * Lz dx  i  * ( y x  x y ) dx , t  0
Rd
Lemma
Dynamical laws
d Lz (t )
dt
Rd
(Bao & Y. Zhang, Math. Mod. Meth. Appl. Sci, 05’)
 2 
 (   )  xy | ( x, t ) | dx, t  0
2
x
2
y
Rd
For any initial data, with symmetric trap, i.e.  x   y , we have
Lz (t )  Lz (0),
Numerical test
E ,0 ( )  E ,0 ( 0 ),
next
t  0.
Angular momentum
expectation
Energy
back
Dynamics of condensate width
Definition:




 r (t )   ( x 2  y 2 ) |  ( x , t |2 dx ,   (t )    2 |  ( x , t |2 dx
Rd
Rd
Dynamic laws (Bao & Y. Zhang, Math. Mod. Meth. Appl. Sci, 05’)
– When
d  2 & x   y
for any initial data:
d 2 r (t )
2

4
E
(

)

4

 r (t ),

0
x
2
dt
im
with initial data  0 ( x, y)  f (r) e
1
 x (t )   y (t )   r (t ),
t0
Numerical Test
2
– For any other cases:
– When
next
t 0
d  2 & x   y
d 2 (t )
2

4
E
(

)

4

 (t )  f (t ),

0

2
dt
t 0
Symmetric trap
Anisotropic trap
back
Dynamics of Stationary state with a shift

 
Choose initial data as:  0 ( x)  s ( x  x0 )
The analytical solutions is: (Garcia-Ripoll el al., Phys. Rev. E, 01’)
 ( x, t )  ei t s ( x  x(t )) eiw( x ,t ) ,
– In 2D:
x(t )   x2 x(t )  0,
s
w( x, t )  0, x(0)  x0
y (t )   y2 y (t )  0,
x(0)  x0 ,
y (0)  y0 , x(0)  0,
y (0)  0
– In 3D, another ODE is added
z(t )   z2 z(t )  0,
z(0)  z0 , z(0)  0
Solution of the center of mass
Center of mass: Bao & Y. Zhang, Appl. Numer. Math., 2006
x (t ) :  x |  ( x , t ) | dx   x |  ( x  x (t )) | dx  x (t )
2
2
s
d
d
In a non-rotating BEC:
x(t )  x0 cos( x t ), y(t )  y0 cos( y t ), t  0
– Trajectory of the center
– Pattern Classification:
•
•
•
•
Motion of the solution
next
Each component of the center is a periodic function
In a symmetric trap, the trajectory is a straight segment
2 p
If  y /  x is a rational #, the center moves periodically with period
If  y /  x is an irrational #, the center moves chaotically, envelope is a rectangle
back
back
Numerical methods for dynamics
Lattice Boltzmann Method (Succi, Phys. Rev. E, 96’; Int. J. Mod. Phys., 98’)
Explicit FDM (Edwards & Burnett et al., Phys. Rev. Lett., 96’)
Particle-inspired scheme (Succi et al., Comput. Phys. Comm., 00’)
Leap-frog FDM (Succi & Tosi et al., Phys. Rev. E, 00’)
Crank-Nicolson FDM (Adhikari, Phys. Rev. E 00’)
Time-splitting spectral method (Bao, Jaksch&Markowich, JCP, 03’)
Runge-Kutta spectral method (Adhikari et al., J. Phys. B, 03’)
Symplectic FDM (M. Qin et al., Comput. Phys. Comm., 04’)
Time-splitting spectral method (TSSP)
Time-splitting:
Step 1:
Step 2:
1
i  t ( x, t )   2  ,
2
i  t ( x , t )  Vd ( x ) ( x , t )   d | ( x , t ) |2  ( x , t )
 | ( x , t ) || ( x , tn ) |
 ( x , tn 1 )  e
 i (Vd ( x )   d | ( x ,tn )|2 ) t
 ( x , tn )
For non-rotating BEC
– Trigonometric functions (Bao, D. Jaksck & P. Markowich, J. Comput. Phys., 03’)
– Laguerre-Hermite functions (Bao & J. Shen, SIAM Sci. Comp., 05’)
Time-splitting spectral method
Properties of TSSP
–
–
–
–
–
–
Explicit, time reversible & unconditionally stable
Easy to extend to 2d & 3d from 1d; efficient due to FFT
Conserves the normalization
Spectral order of accuracy in space
2nd, 4th or higher order accuracy in time
Time transverse invariant


Vd ( x)  Vd ( x)  


|  ( x, t ) |2
unchanged
– ‘Optimal’ resolution in semicalssical regime
h  O   ,
k  O   ,
  1/ d 2 /(2d )
Dynamics of Ground states
1d dynamics: 1  100 at t  0, x  4x
2d dynamics of BEC (Bao, D. Jaksch & P. Markowich, J. Comput. Phys., 03’)
– Defocusing: 2  20, at t  0 x  2x ,  y  2 y
– Focusing (blowup): At t  0  2  40  50
3d collapse and explosion of BEC (Bao, Jaksch & Markowich,J. Phys B, 04’)
– Experiment setup leads to three body recombination loss



1 2
i  ( x, t )      V ( x )   | |2   i 0  2 | |4 
t
2
– Numerical results:
• Number of atoms , central density & Movie
next
back
back
back
Collapse and Explosion of BEC
back
Number of atoms in condensate
back
Central density
back
back
Central quantized vortices
Central vortex states in 2D:
i m t
 ( x, y , t )  e
with
i m t
m ( x, y )  e
1 d  dm (r )   m 2 r 2 
2
 m m ( r )  
r
   2   m (r )   2 | m | m ,
2r dr 
dr   2r
2
2
m ( r ) e
0r

 m (r ) r dr  1
2
0
Vortex Dynamics
– Dynamical stability
– Interaction
• Pattern I
• Pattern II



V ( x )  V ( x )  W ( x, t )
N
 0 ( x, y)   n ( x  x j , y  y j ) / ||  ||
j 1
N
j
 0 ( x, y)  n ( x  x j , y  y j ) / ||  ||
j 1
j
im 
Central Vortex states
Central Vortex states
Vortex stability & interaction
Dynamical stability (Bao & Y. Zhang, Math. Mod. Meth. Appl. Sci., 05’)
– m=1: stable velocity
– m=2: unstable velocity
Interaction (Bao & Y. Zhang, Math. Mod. Meth. Appl. Sci., 05’)
– N=2: Pair: velocity trajectory phase phase2
Anti-pair: phase trajectory angular trajectory2
– N=3: velocity trajectory
– Pattern II: Linear
nonlinear
Interaction laws:
– On-going with Prof. L. Fu & Miss Y. Zhang
next
back
back
back
back
back
Linear case
back
Noninear case: BEC
back
Linear case
back
back
Linear case
back
Linear case
back
back
back
back
Some Open Questions
Dynamical laws for vortex interaction
dx j (t )
dt
 ?????
With a quintic damping, mass goes to zero
N (t ) 
 | ( x, t ) |
2
d
t 
dx

0  0 ????
Convergence & error estimate of the TSSP?
Energy diminishing of the gradient flow in
nonlinear case & error estimate ?
Part IV
Rotating BEC
&
multi-component BEC
Rotating BEC

The Schrodinger equation ( x  ( x, y, z ) )

 ψ(x, t)
 H ( )
i

t
*
– The Hamiltonian:






H ( )   * ( x , t ) [
 V ( x )   Lz ] ( x , t ) dx
2m


 


 
1



  * ( x , t )  * ( x , t ) ( x  x )  ( x , t )  ( x , t ) dx dx
2
2
2
– The interaction potential is taken as in Fermi form

( x)  ( N 1)
U
0

 ( x),
U 0  4 
2
a
s
m
.
Rotating BEC
The 3D Gross-Pitaevskii equation ( x  ( x, y, z) )
2



2

i   ( x , t )  [
  V ( x )   Lz  N U 0 |  |2 ]
t
2m
– Angular momentum rotation
   
Lz : xpy  ypx  i( x y  y x )  i , L  x  P, P  i
– V is a harmonic trap potential

V (x) 
m
(
2

2
x
x2  y y 2  z z 2 )
2
2
– Normalization condition

R3
 2 
|  ( x , t ) | dx  1
Rotating BEC
General form of GPE (

x Rd
)




1 2
2
i  ( x, t )  [   Vd ( x )   Lz   d | | ] ( x, t )
t
2
with
Lz : i( x y  y x )  i 


d  



z
(
z
)
dz


,

 3
2
,
4
2
1 2

2
(
x

y
),
d 2

 
y
2
Vd ( x )  
2
2
1
 ( x 2   y 2   z 2 ), d  3
y
z
2
Normalization condition
 2 
 d |  ( x, t ) | dx  1.
R
Rotating BEC
Conserved quantities
– Normalization of the wave function
N ( (t )) 



| ( x , t ) | 2 dx  N ( (0))  1
Rd
– Energy


d
1
2
2
4
[
|


|

V
(
x
)
|

|



*
L


|

|
]
d
x
d
z
Rd
2
2
 E  , ( (0))
E  , ( (t ))  
Chemical potential

 ,
( (t ))  
Rd


1
[ |  |2 Vd ( x ) |  |2  * Lz   d | |4 ] dx
2
Semiclassical scaling
When
With
d  1
, re-scaling


x   1/ 2 x      d / 4   1/ d2 /(d 2)

  
2 2
i   ( x, t )  [   Vd ( x )    Lz  |  |2 ] 
t
2
E , ( )  

Rd
[
2


1
|   |2 Vd ( x ) |   |2   (  ) * Lz    |  |4 ] dx
2
2
 O(1)
Leading asymptotics
  


E , ( )   1E , (  )  O  1  O d2 /(d 2) ,  , ( )  O d2 /(d 2)

Quantum Hydrodynamics
Set


   e
iS  / 

, v   S  ,


J    v  ,   1 /  d2 /(d  2 )
Geometrical Optics: (Transport + Hamilton-Jacobi)
 t      (   S  )  Lˆ z    0,
Lˆ z : ( x y  y x )  
2

1

 2


 t S  S  Vd ( x )    Lˆ z S 
2
2

1

 
Quantum Hydrodynamics (QHD): (Euler +3rd dispersion)

 t      (   v  )  Lˆ z    0
 


  2
J J


t (J )    (
)  P(  )   Vd  Lˆ z J  A J  (    ln   )


4

 0 1


P(  )   2 / 2 ,
J   v ,
A  
 1 0 
Stationary states
Stationary solutions of GPE


i  t
 ( x, t )  e  ( x )
Nonlinear eigenvalue problem with a constraint


 2

1 2
  ( x )  [   Vd ( x )   Lz   d |  ( x ) | ]  ( x ),
2
 2 
 d |  (x) | dx  1
R
Relation between eigenvalue and eigenfunction
    , ( )  E , ( ) 
d
2

Rd
 4 
|  (x) | dx
Stationary states
Equivalent statements:
– Critical points of E , ( ) over the unit sphere
– Eigenfunctions of the nonlinear eigenvalue problem
– Steady states of the normalized gradient flow:
 ( )


1
t ( x, t )  [  2  V ( x )   Lz   |  |2   , 2 ]  ,
2
||  ||



 ( x,0)  0 ( x )
with
|| 0 ( x ) || 1.
Minimizer/saddle points over the unit sphere :
– For linear case
d  0 &   1
• Global minimizer vs saddle points
– For nonlinear case d  0 &   1
• Global minimizer, local minimizer (?) vs saddle points
Ground state
Ground state:
E , ( g )  min E , ( ),
|| || 1
 g    , ( g )  E , ( g ) 
d
2

R
 4 
| g ( x ) | dx
d
Existence: |  | 1 & d  0
– Seiringer (CMP, 02’)
Uniqueness of positive solution:   0 & d  0
– Lieb et al. (PRA, 00’)
Energy bifurcation:
0  c  1
– Aftalion & Du (PRA, 01’); B., Markowich & Wang 04’
Numerical results
–
Ground states:
in 2D
in 3D
isosurface
Quantized vortex generation in 2D
–
•
–
surface
contour
Vortex lattice
•
–
Symmetric trapping
anisotropic trapping
Giant vortex generation in 2D
•
–
surface
contour
Giant vortex
•
In 2D
In 3D
next
back
back
back
back
back
back
back
back
back
back
back
Numerical & Asymptotical results
Critical angular frequency: symmetric state vs quantized vortex state
0  c  1, c  E ,0 (1 )  E ,0 (0 )
Asymptotics of the energy:
Eg,  Eg,1  O(1 ),   1; E , (s )  O( 2 /(2d ) ),   
Ratios between energies of different states
lim
 d 
Es ,
Eg,
 s ,
 1, lim g  1,   0;
  
 ,
d
lim
 d 
Es ,
Eg,
 s ,
 const., lim g  const, 0 |  | 1
  
 ,
d
Rank according to energy and chemical potential
– Stationary states are ranked according to their energy, then their chemical
potential are in the same order.
Next
back
back
back
Dynamical laws of rotating BEC
Time-dependent Gross-Pitaevskii equation



1 2
i  ( x, t )  [   Vd ( x )   Lz   d | |2 ]
t
2


 ( x,0)   0 ( x ),
Lz : i( x y  y x )  i
Dynamical laws
–
–
–
–
–
Time reversible & time transverse invariant
Conservation laws
Angular momentum expectation
Condensate width
Dynamics of a stationary state with its center shifted
Conservation laws
Conserved quantities
– Normalization of the wave function
N ( (t )) 



| ( x , t ) | 2 dx  N ( (0))  1
Rd
– Energy


d
1
2
2
4
[
|


|

V
(
x
)
|

|



*
L


|

|
]
d
x
d
z
Rd
2
2
 E  , ( (0))
E  , ( (t ))  
Chemical potential

 ,
( (t ))  
Rd


1
[ |  |2 Vd ( x ) |  |2  * Lz   d | |4 ] dx
2
Angular momentum expectation


Lz (t ) :  * Lz dx  i  * ( y x  x y ) dx , t  0
Definition:
Lemma The dynamics of
Rd
d Lz (t )
dt
Rd
Lz (t )
satisfies
 2 
 (   )  xy | ( x, t ) | dx, t  0
2
x
2
y
Rd
For any initial data, with symmetric trap, i.e.  x   y , we have
Lz (t )  Lz (0),
Numerical test
E ,0 ( )  E ,0 ( 0 ),
next
Bao, Du & Zhang, SIAM J. Appl. Math., 66 (2006), 758
t  0.
Angular momentum
expectation
Energy
back
Dynamics of condensate width
Definition:
Bao, Du & Zhang, SIAM J. Appl. Math., 66 (2006), 758




 r (t )   ( x 2  y 2 ) |  ( x , t |2 dx ,   (t )    2 |  ( x , t |2 dx
Rd
Rd
Dynamic laws
– When
d  2 & x   y
for any initial data:
d 2 r (t )
2

4
E
(

)

4

 r (t ),

,

0
x
2
dt
im
with initial data  0 ( x, y)  f (r) e
1
 x (t )   y (t )   r (t ),
t0
Numerical Test
2
– For any other cases:
– When
next
t 0
d  2 & x   y
d 2  (t )
2

4
E
(

)

4

 (t )  f (t ),

,

0

2
dt
t 0
Symmetric trap
Anisotropic trap
back
Dynamics of Stationary state with a shift

 
Choose initial data as:  0 ( x)  s ( x  x0 )
The analytical solutions is: Bao, Du & Zhang, SIAM J. Appl.Math., 2006
 ( x, t )  ei t s ( x  x(t )) eiw( x ,t ) ,
w( x, t )  0, x(0)  x0
– In 2D:
x(t )  2y (t )  ( x2   2 ) x(t )  0,
s
y (t )  2x(t )  ( y2   2 ) y (t )  0,
x(0)  x0 ,
y (0)  y0 , x(0)  y0 ,
y (0)  x0
– In 3D, another ODE is added
z(t )   z2 z(t )  0,
z(0)  z0 , z(0)  0
Solution of the center of mass
Center of mass: Bao & Zhang, Appl. Numer. Math., 2006
x (t ) :  x |  ( x , t ) | dx   x |  ( x  x (t )) | dx  x (t )
2
2
s
d
d
In a non-rotating BEC:   0
x(t )  x0 cos( x t ), y(t )  y0 cos( y t ), t  0
– Pattern Classification:
•
•
•
•
Each component of the center is a periodic function
In a symmetric trap, the trajectory is a straight segment
2 p
If  y /  x is a rational #, the center moves periodically with period
If  y /  x is an irrational #, the center moves chaotically, envelope is a rectangle
Solution of the center of mass
In a rotating BEC with a symmetric trap:
x0
|  | y0
cos(a t )  cos(bt ) 
sin(a t )  sin(bt ) ,
2
2
y
|  | x0
y (t )  0  cos(a t )  cos(bt ) 
  sin(a t )  sin(bt ) ,
2
2
x(t ) 
a  x ||
b  x||
| x (t ) |: x 2 (t )  y 2 (t )  x02  y02 | cos( x t ) |
–
–
–
–
Trajectory of the center
Distance between the center and trapping center
Motion of the solution: 0.5
1 2 4
Pattern Classification:
next
1/5,
back
4/5,
1
3/2, 6,
Pi
back
back
back
back
back
Pattern Classification
Pattern Classification: Bao & Zhang, Appl. Numer. Math., 2006
– The distance between the center and trap center is periodic function
  q/ p
– When
is a rational #
• The center moves periodically
• The graph of the trajectory is unchanged under a rotation
– When   q / p is an irrational #,
• The center moves chaotically
• The envelope of the trajectory is a circle
back
– The solution of GPE agrees very well with those from the ODE system
Solution of the center of mass
In a rotating BEC with an anisotropic trap
– When
|  |  x or  y
results
• The trajectory is a spiral coil to infinity
• The trajectory is an ellipse
– Otherwise
result1
result2
• The center moves chaotically & graph is a bounded set
• The center moves along a straight line to infinity
next
back
back
back
Total density with dissipation
Time-dependent Gross-Pitaevskii equation

1 2
i



(
x
,
t
)

[

  Vd ( x )  W ( x, t )   Lz  d |  |2 ]
 
t
 ( x, 0)   0 ( x ),
2
Lz : i( x y  y x )  i
Lemma The dynamics of total density satisfies
d t d
1  2
N ( )(t ) 
| ( x, t ) |2 d x  
 , ( )  0, t  0
d
2
– The total density decreases when   0&|  | min{ x ,  y }
density function
energy
next
back
back
Numerical Methods
Time-splitting pseudo-spectral method (TSSP)
Step 1 :
Step 2 :

1
i  t ( x , t )    2    Lz  ,
Lz : i( x y  y x )  i
2







i  t ( x , t )  Vd ( x ) ( x , t )   d | ( x , t ) |2  ( x , t )  | ( x , t ) || ( x , t n ) |
– Use polar coordinates (B., Q. Du & Y. Zhang, SIAP 06’)
– Time-splitting + ADI technique (B. & H. Wang, JCP, 06’)
– Generalized Laguerre-Hermite functions (B., J. Shen & H. Wang, 06’)
Numerical methods for rotating BEC
Numerical Method one: (Bao, Q. Du & Y. Zhang, SIAM, Appl. Math. 06’)
– Ideas
• Time-splitting
• Use polar coordinates: angular momentum becomes constant coefficient
• Fourier spectral method in transverse direction + FD or FE in radial direction
• Crank-Nicolson in time
– Features
•
•
•
•
•
Time reversible
Time transverse invariant
Mass Conservation in discretized level
Implicit in 1D & efficient to solve
Accurate & unconditionally stable
Numerical methods for rotating BEC
Numerical Method two: (Bao & H. Wang, J. Comput. Phys. 06’)
– Ideas
• Time-splitting
• ADI technique: Equation in each direction become constant coefficient
• Fourier spectral method
– Features
•
•
•
•
•
Time reversible
Time transverse invariant
Mass Conservation in discretized level
Explicit & unconditionally stable
Spectrally accurate in space
Dynamics of ground state
Choose initial data as:  100,   0.8,


 0 ( x)  g ( x) : ground state
 y   z 1
Change the frequency in the external potential:
Case 1: symmetric:  x :1  2 &  y :1  2
surface
contour
– Case 2: non-symmetric:  x :1 1.8 &  y :1  2.2
surface
contour
– Case 3: dynamics of a vortex lattice with 45 vortices:

image contour   1000,   0.9, V ( x, t ) : anisotropic
next
–
back
back
back
back
back
back
Interaction of two vortices in linear
0
  1/ 3
Interaction of two vortices in linear
  1/ 2
 1
Interaction of two vortices in linear
  1/ 
 
Interaction of vortices in nonlinear
  1/ 
 
Interaction of vortices in nonlinear
0
  1/ 2
Interaction of vortices in nonlinear
 1
4
Interaction of vortices in nonlinear
  1/ 
 
Some Open Questions
Dynamical laws for vortex interaction
dx j (t )
dt
 ?????
With a quintic damping, mass goes to zero
N (t )   | ( x, t ) | d x

0  0 ????
Semiclassical limit when initial data has
vortices???
Vortex line interaction laws, topological change?
What is a giant vortex?
2
d
t 
Two-component BEC
The 3D coupled Gross-Pitaevskii equations
2
i

 1 ( x , t )  [
 2  V ( x )   Lz  U11 |  1 |2 U12 |  2 |2 ]  1    2
t
2m
i

 2 ( x , t )  [
 2  V ( x )   Lz  U 21 |  1 |2 U 22 |  2 |2 ]  2    1
t
2m
2
Normalization
conditions
2
N (t )    | j ( x, t ) |2 dx  N10  N20 : N
j 1
with
3
3
Intro- & inter-atom Interactions
4  2 a jl
U jl 
m
N 0j   |  j ( x,0) |2 dx,
with
a12  a21
Two-component BEC
Nondimensionalization

1
 1 ( x , t )  [  2  V ( x )   Lz  11 |  1 |2  12 |  2 |2 ] 1   2
t
2

1
i
 2 ( x , t )  [  2  V ( x )   Lz   21 |  1 |2   22 |  2 |2 ] 2   1
t
2
i
Normalization conditions
– There is external driven field   0
N (t )   |  ( x , t ) | dx   |  ( x , t ) | dx  1
2
2
1
3
2
3
– No external driven field   0
N10
3 |1 ( x, t ) | dx  N ,
2
N20
3 | 2 ( x, t) | dx  N
2
Two-component BEC
Energy
E ()  
Rd
2 
1
jl
2
2
[ ( |  j | V ( x )| j |  j * Lz j  
| j |2 | l |2 )  2 Re( 1* 2 )]dx
j=1 2
l 1 2
2
Reduction to one-component:
N20
N2 (t )   | 2 ( x, t ) | dx 
: 
N
3
2
i
  0, N10
N20 , N10  O( N )
N10
1, N1 (t )   |  1 ( x, t ) | dx 
: 1    1
N
3
2

1
 ( x, t )  [ 2  V ( x )   Lz   |  |2 ] ( x, t ),
t
2
 ( x, t )  N / N  1 ( x, t ) &  =N 11 / N
0
1
0
1
| E ()  Es ( ) |
 O( )
Es ( )
Two-Component BEC
Semiclassical scaling

2 2
i   1 ( x , t )  [   V ( x )   Lz  11 |  1 |2 12 |  2 |2 ] 1   2
t
2

2 2
i   2 ( x , t )  [   V ( x )   Lz   21 |  1 |2  22 |  2 |2 ] 2   1
t
2
Semiclassical limit
– No external field:   0
• WKB expansion, two-fluid model
– With external field:   0
• WKB expansion doesn’t work, Winger transform
Ground state
No external field:
min
||1 ||  ,||2 ||  
 0
E (1 , 2 ) with     1
Nonlinear eigenvalue problem
1
2
1
2 2 ( x )  [  2  V ( x )   Lz   21 | 1 |2   22 | 2 |2 ] 2
2
1 1 ( x )  [  2  V ( x )   Lz  11 | 1 |2  12 | 2 |2 ] 1
Existence & uniqueness of positive solution
Numerical methods can be extended
Ground states
crater
Ground state
With external field:
0
min 2 E (1 , 2 )
||1 || ||2 || 1
2
Nonlinear eigenvalue problem
1
2
1
 2 ( x )  [  2  V ( x )   Lz   21 | 1 |2   22 | 2 |2 ] 2  1
2
 1 ( x )  [   2  V ( x )   Lz  11 | 1 |2  12 | 2 |2 ] 1  2
Existence & uniqueness of positive solution ???
Numerical methods can be extended????
Dynamics
Dynamical laws:
–
–
–
–
Conservation of Angular momentum expectation
Dynamics of condensate width
Dynamics of a stationary state with a shift
Dynamics of mass of each component, they are
periodic function when 11  12  22
– Vortex can be interchanged!
Numerical methods
– Time-splitting spectral method
Dynamics
Dynamics
Spinor BEC
Spinor F=1 BEC
2
i

 1  [
 2  V ( x )   Lz  g n  ] 1  g s ( 1   0   1 ) 1  g s *1 02
t
2m
2
i

 0  [
 2  V ( x )   Lz  g n  ] 0  g s ( 1   1 ) 1  2 g s 1 1 0*
t
2m
2
i

 1  [
 2  V ( x )   Lz  g n  ] 1  g s (  1   0  1 ) 1  g s 1* 02
t
2m
With
4 2 a0  2a2
4 2 a2  a0
  1  0  1 ,  j | j | , gn 
, gs 
m
3
m
3
a0 , a2 : s-wave scattering length with the total spin 0 and 2 channels
2
Spinor BEC
Total mass conservation
N (t )    | ( x, t ) | dx  N  N  N
1
2
j
j 1
0
1
0
0
0
1
: N
with
3
N 0j   |  j ( x,0) |2 dx,
3
Total magnetization conservation
M (t ) 
2
2
0
0
|

(
x
,
t
)
|
dx

|

(
x
,
t
)
|
dx

N

N
1
1 : M
 1
 1
3
3
Energy conservation
E ( )  
1
Rd
[ (
j=-1
2
2m
|  j |2 V ( x )|  j |2  j * Lz j ) 
gn 2
 
2
gs 2
( 1   21  2 1 0  2  1 0  2 1 1 )  g s ( *1 02 1*  1 ( 0* ) 2 1 )]dx
2
Spinor BEC
Dimension reduction
Ground state
– Existence & uniqueness of positive solution??
– Numerical methods ???
Dynamics
– Dynamical laws
– Numerical methods: TSSP
Semiclassical limit & hydrodynamics equation??
BEC at Finite Temperature
Condensate coexists with noncondensed thermal cloud
Coupled equations of motion
for condensate and thermal
cloud
Mean-field theory in
collisionless regime
ZGN theory in collision
dominated regime
Mean-field Theory
Evolution of quantum field operator
ˆ  2 2
ˆ

ˆ †
ˆ
ˆ
i
 
  Vext    g
t  2m

g  4 2as / m

ˆ ( x, t ) is the annihilation field operator
where 
ˆ † ( x, t ) is the creation field operator
and 
Mean-field description
 
 
~ 
 ( x , t )   ( x , t ) aˆ 0   ( x , t )
 

( x, t )  ( x, t )
~ 
( x, t )  0
Condensate wavefunction

   2 2
~†~ ~
~*  g 
i
 
  Vext  gnc  2 gnT    gm

t  2m

Mean-field Theory
Generalized GPE for condensate wavefunction

   2 2
~†~ ~
~*  g 
i
 
  Vext  gnc  2 gnT    gm

t  2m


2
nc ( x , t )  
condensat edensit y

~ ~
nT ( x , t )   † 
non - condensat edensit y
~~
~ ( x, t )  
m

off - diagonal non - condensat edensit y
~ ~~
 † 
hree
t
- field correlat ion funct ion
Temperature-dependent fluctuation field for non-condensate
~
~
   2 2
~
~ ~~
i
 
  Vext  2 gn   gm †  g  † 
t  2m

~ ~
~ ~



 †   †
n( x, t )  nc ( x, t )  nT ( x, t )
~~
~~


  
2 
~
m( x, t )   ( x, t )  m( x, t )
~ ~~
~ ~ ~
~~ ~
 †   2  †      †
Hartree-Fock Bogoliubov Theory
Ignore the three-field correlation function
~ ~~
 † 

   2 2
~ *
i
 
  Vext  g nc  2nT    gm
t  2m

~
2
~
  
~
i
 
 2  Vext  2 gn   gm †
t  2m

Bogoliubov transformation



~ 
 ( x , t )   u j ( x , t )ˆ j  v*j ( x , t )ˆ j†
j




~ 
 † ( x , t )   u *j ( x , t )ˆ j†  v j ( x , t )ˆ j
j

whereˆ (ˆ j ) creates (annihilates) a Bogoliubov quasiparticle of
energy εj
The quasiparticles are non-interacting
†
j
Hartree-Fock Bogoliubov Theory
Bogoliubov equations for non-condensate
u j
 2 2

i
 
  Vext  2 gnu j  gm vj
t  2m

v j   2 2

 i
 
  Vext  2 gn v j  gm*u j
t  2m

where

2
nc ( x , t )  
2
2

nT ( x , t )   u j N j  v j ( N j  1)

j

*
m( x , t )   2   u j v j (1  2 N j )
j
N j  ˆ †j ˆ j 
1
exp( j / kT )  1

Time-independent Hartree-Fock Bogoliubov Theory
Stationary states


 ( x , t )   ( x ) e i t / 

 i  t / 
u j ( x , t )  u j ( x ) e j e i t / 

 i  t / 
v j ( x , t )  v j ( x ) e j ei t / 
Time-independent generalized GPE and Bogoliubov equations
 2 2

~ *  





V

g
n

2
n


g
m
ext
c
T
 2m



 2 2




V



2
gn
ext
 2m
u j  gm vj   j u j


 2 2

*



V



2
gn
v

gm
u j   j v j
ext
 2m
 j


HFB-Popov Approximation
HFB produces an energy gap in the excitation spectrum
~
Solution: leave out m
Generalized GPE and Bogoliubov equations within Popov
approximation (gapless spectrum)
 2 2
 *











i
Vext g nc 2nT  

t  2m

u j
 2 2

i
 
  Vext  2 gnu j  g 2 v j
t  2m

v j   2 2

2
 i
 
  Vext  2 gn v j  g  * u j
t  2m

 
Hartree-Fock Approximation
Approximate Bogoliubov excitations with single-particle excitations,
i.e. let v j  0

   2 2
i
 
  Vext  g nc  2nT   *
t  2m

 j   2 2

i
 
  Vext  2 gn j
t  2m


2
nc ( x , t )  
2

nT ( x , t )    j N j
j
Restricted to finite temperature close to Tc, where the non-condensed
particles have higher energies
ZGN Theory
Mean-field theory deals with BEC in collisionless region (low
density thermal cloud):
l >> 
l is the collisonal mean-free-path of excited particles
 is the wavelength of excitations
In collision-dominated region l <<  (higher density thermal
cloud), the problem becomes hydrodynamic in
nature
ZGN theory (E. Zaremba, A. Griffin, T. Nikuni, 1999)
describes finite-T BEC with interparticle collisions in
the semi-classical limit
kBT >> ħ0 0 : trap frequency)
kBT >> gn
ZGN Theory
~ ) but include the three-field
Apply Popov approximation (ignore m
~†~ ~
correlation function 

GPE for condensate wavefunction

   2 2
~ ~~
i
 
  Vext  g nc  2nT    g  † 
t  2m


   2 2
i
i
 
  Vext  g nc  2nT  
12 [ f ] 
t  2m
2nc

Quantum Boltzmann equation for phase-space distribution function of
non-condensate

 
 
f ( x, p, t )
 p 


 C12 [ f ]  C22 [ f ]
 t  m   x  Vext  2 gn  p  f ( x, p, t ) 
t
collision

source term 12 [ f ]  


2g
dp
~ ~~
Im *  †   
C12 [ f ]
3

(2)
ZGN Theory
Thermal cloud density


 
dp
nT ( x , t )  
f
(
x
, p, t )
3
2 
Collision between condensate and non-condensate
-- transfer atoms from/to the condensate
2 g 2 nc



  

C12 [ f ] 
d
p
d
p
d
p

(
m
v

p

p

p
1
2
3
s
1
2
3)
(2 ) 2  4 
 
 
 
  ( c   p1   p 2   p 3 ) ( p  p1 )   ( p  p2 )   ( p  p3 )
 (1  f1 ) f 2 f 3  f1 (1  f 2 )(1  f 3 )
Collision between non-condensate particles
 
    
2g 2
C22[ f ] 
d
p
d
p
d
p
1
2
3 ( p  p1  p2  p3 )
5 7 
(2 ) 
  ( p   p1   p 2   p 3 ) (1  f )(1  f1 ) f 2 f 3  f f1 (1  f 2 )(1  f 3 )


ZGN Theory
Energy of condensate atoms


1
2

 c ( x , t )   c ( x , t )  mv s2 ( x , t )
Local chemical potential

 2  nc
c ( x, t )  
 Vext  gnc  2 gnT
2m
nc
2
Superfluid velocity
 


v s ( x , t )   ( x , t )
m



( x, t )  nc ( x, t ) expi ( x, t )
Energy of non-condensate atoms
– Hartree-Fock energy
2

p
 p ( x, t ) 
 Vext  2 gn
2m
Limited to high temperature (close to Tc)
For lower temperature, the spectrum of excited atoms should be
described by Bogoliubov approximation
Open questions
Mathematical theory
– Quantum Boltzmann Master equation (QBE)
– GPE with damping term
– Coupling QBE +GPE
Numerical methods
–
–
–
–
For QBE: P. Markowich & L. Pareschi (Numer. Math., 05’)
For QBE+GPE
Comparison with experiments
Rotational frame
Conclusions
–
–
–
–
–
–
–
–
Review of BEC
Experiment progress
Mathematical modeling
Efficient methods for computing ground & excited states
Efficient methods for dynamics of GPE
Comparison with experimental results
Vortex dynamics
Quantized vortex stability & interaction
Future Challenges
–
–
–
–
–
–
–
–
–
Multi-component BEC for bright laser
Applications of BEC in science and engineering
Precise measurement
Fermions condensation, BEC in solids & waveguide
Dynamics in optical lattice, atom tunneling
Superfluidity & dissipation, quantized vortex lattice
Coupling GPE & QBE for BEC at finite temperature
Mathematical theory for BEC
Interdisciplinary research: experiment,physics, mathematics, computation, ….
References
[1] M.H. Anderson, J.R. Ensher, M.R. Matthews, C.E. Wieman and E.A. Cornell, Science
269 (1995) 198-201.
[2] W. Bao, J. Shi and P.A. Markowich, J. Comput. Phys. , Vol. 175, pp. 487-524, 2002.
[3] W. Bao and W.J. Tang, J. Comput. Phys., Vol. 187, No. 1, pp. 230 - 254, 2003.
[4] W. Bao, D. Jaksch and P.A. Markowich, J. Comput. Phys., Vol. 187, No. 1, pp. 318 342, 2003.
[5] W. Bao, S. Jin and P.A. Markowich, SIAM J. Sci. Comput., Vol. 25, No. 1. pp. 27-64,
2003.
[6] W. Bao and D. Jaksch, SIAM J. Numer. Anal., Vol. 41, No. 4. pp. 1406-1426, 2003.
[7] W. Bao, D. Jaksch and P.A. Markowich, J. Phys. B: At. Mol. Opt. Phys., Vol. 37, No. 2,
pp. 329-343, 2004.
[8] W. Bao, Multiscale Modeling and Simulation: a SIAM Interdisciplinary Journal, Vol. 2,
No. 2. pp. 210-236, 2004.
[9] W. Bao and Q. Du, SIAM J. Sci. Comput. , Vol. 25, No. 5. pp. 1674-1697, 2004.
References
[9] W. Bao and Q. Du, SIAM J. Sci. Comput. , Vol. 25, No. 5. pp. 1674-1697, 2004.
[10] W. Bao, H.Q. Wang and P.A. Markowich, Comm. Math. Sci. , Vol. 3, No. 1, pp. 57-88,
2005.
[11] W. Bao, P.A. Markowich, C. Schmeiser and R. M. Weishaupl, Math. Mod. Meth. Appl.
Sci. , Vol. 15, No. 5, pp. 767-782, 2005.
[12] W. Bao and J. Shen, SIAM J. Sci. Comput. , Vol. 26 , No. 6, pp. 2010-2028, 2005.
[13] W. Bao and Y.Z. Zhang, Math. Mod. Meth. Appl. Sci. , Vol. 15 , No. 12, pp. 1863-1896,
2005.
[14] W. Bao, Qiang Du and Yanzhi Zhang, SIAM J. Appl. Math., Vol. 66 , No. 3, pp. 758786, 2006.
[15] W. Bao and H. Wang, J. Comput. Phys., Vol. 217, No. 2, pp. 612-626, 2006.
[16] W. Bao, I-L. Chern and F. Y. Lim, J. Comput. Phys., Vol. 219, No. 2, pp. 836-854, 2006
[17] W. Bao and Y. Zhang, Appl. Numer. Math., Vol. 57, No. 5-7, pp. 697-709, 2007.
References
[18] W. Bao, F. Y. Lim and Y. Zhang, Bulletin of the Institute of Mathematics, Academia
Sinica, Vol. 2, No. 2, pp. 495-532, 2007.
[19] W. Bao, H.L. Li and Y. Zhang, Physica D: Nonlinear Phenomena, Vol. 234, pp. 49-69,
2007.
[20] W. Bao, Y. Ge, D. Jaksch, P. A. Markowich and R. M. Weishaeupl, Comput. Phys.
Comm., Vol. 177, No. 11, pp. 832-850, 2007.
[21] W. Bao and H. Wang, A mass and magnetication conservative and energy diminishing
numerical method for computing ground state of spin-1 Bose-Einstein condensates,
SIAM J. Numer. Anal., Vol. 45, No. 5, pp. 2177-2200, 2007.
[22]. A. Klein, D. Jaksch, Y. Zhang and W. Bao, Dynamics of vortices in weakly interacting
Bose-Einstein condensates, Phys. Rev. A, Vol. 76, article 043602, 2007.
[23]. W. Bao and M.-H. Chai, A uniformly convergent numerical method for singularly
perturbed nonlinear eigenvalue problems, Commun. Comput. Phys., to appear.
[24]. W. Bao and F. Y. Lim, Computing Ground States of Spin-1 Bose-Einstein
Condensates by the Normalized Gradient Flow, arXiv: 0711.0568; SIAM J. Sci.
Comput., to appear.
References
[25] Bradly et al., Phys. Rev. Lett., 75 (1995), 1687.
[26] Davis et al., Phys. Rev. Lett., 75 (1995), 3969.
[27] A.L. Fetter and A. A. Svidzinsky, Vortices in a trapped dilute Bose-Einstein condensate
(topical review), J. Phys.: Condens. Matter 13 (2001), 135-194.
[28] A.J. Leggett, Bose-Einstein condensation in the alkali gases: some fundamental
concepts, Rev. Modern Phys., 73 (2001), 307-356.
[29] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation, Oxford University Press,
2003.
[30]E.H. Lieb, R. Seiringer, J.P. Solovej and J. Yngvason, The Mathematics of the Bose
Gas and its Condensation, Birkhauser, 2000.
[31] A. Aftalion, Vortices in Bose-Einstein Condensates, Birkhauser, 2006.
[32] F. Dalfovo, S. Giorgini, L.P. Pitaevskii and S. Stringari, Theory of Bose-Einstein
condensation in trapped gases, Rev. Modern. Phys., 71 (1999), 463-512.
[33] C.J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases, Cambridge
University Press, 2002.