Numerical implementation of the multiscale and averaging methods for quasi periodic systems

0.55 MB·English

by  Tal Kachman

NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND
AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS
TALKACHMAN,SHMUELFISHMAN,ANDAVYSOFFER
6
1 Abstract. We consider the problem of numerically solving the Schr¨odinger
0 equation with a potential that is quasi periodic in space and time. We in-
2 troduceanumericalschemebasedonanewlydevelopedmulti-timescaleand
averaging technique. We demonstrate that with this novel method we can
l
u solve efficiently and with rigorous control of the error such an equation for
J long times. A comparison with the standard split-step method shows sub-
stantialimprovementincomputationtimes,besidesthecontrollederrors. We
5
apply this method for a free particle driven by quasi-periodic potential with
2
manyfrequencies. ThenewmethodmakesitpossibletoevolvetheSchro¨dinger

equationfortimesmuchlongerthanwaspossiblesofarandtoconcludethat
]
D thereareregimeswheretheenergygrowthstopsin-spiteofthedriving.
C
.
n Numerical solution Time dependent potentials Multiscale averaging
i
l
n
[ 1. Introduction

4 A method for the study of the dynamics for the Schr¨odinger equation with time
v dependentpotentials[1]isimplementednumerically. Thepotentialisquasiperiodic
2 in both space and time. The power of this new method is demonstrated. Explo-
3
rationofthedynamicsforsuchpotentialsismotivatedbyexperimentsinoptics[2]
1
wherehypertransport, namelytransportfasterthenballisticwasfoundexperimen-
3
0 tally and numerically in some regimes. In theoretical work that followed [2, 3, 4, 5]
. a classical theory was developed for the potentials relevant for these optics exper-
1
0 iments. In particular, it was found in the framework of classical mechanics, that
5 for smooth potentials the spreading in momentum as function of time stops. The
1 calculations of the present paper are for such potentials. For short times, relevant
:
v for the existing experiments it was found that wave or quantum and classical dy-
i namics agree in general features.
X
For long times it turned out impossible to compute numerically the quantum dy-
r
a namics using the standard methods [6, 7], while with the method introduced in
[1] and implemented here, calculations for such long times are feasible as will be
demonstrated in this paper. Such calculations and comparison with the classical
results, is of fundamental importance for the issue of quantum classical correspon-
dence. The main objective of the present paper is to demonstrate the power of the
method introduced in [1] for a physically relevant example.
Potentials which are quasiperodic both in space and time, can manifest a high
degree of complexity and are subject of many studies over the last decade, mostly
in the framework of classical physics [8, 9, 10, 11, 12, 13, 14]. With different
experimental realizations of such potentials there is also a need for a numerical
approach to investigate them and their asymptotic behavior. The standard way
to numerically solve problems with time dependent potentials is based on either
1
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS2
spectral methods [15, 16] or explicit/implicit finite difference schemes [7, 17]. In
this paper we implement numerically a recently developed rigorously controlled
multi-time scale averaging technique [1]. The above mentioned method has two
distinct advantages. The first is that at each step of the averaging hierarchy there
is a well defined and completely known bound on the numerical error. The second
advantage and one which has far greater impact is the reduction in computational
time accompanied with each of the hierarchical averaging steps. This reduction
enablesustogotolongtimescales, impossiblebythemethodsweareawareof. In
this paper we describe and present an implementation for a specific problem.
The outline of the paper is as follows. In Section 2, we present the model
for which the method of multiple scale and averaging (MSA) is implemented. This
modelisofphysicalinterestandimportance. TheMSAmethodthatwasdeveloped
in Ref. [1] and details of its implementation are presented in Secs. 3 and 4. The
MSAmethodinvolvesahierarchyofcomputationswherethefirstlevelisdescribed
is Sec. 3 and the following levels are presented in Section 4. The needed level is
determinedbytherequiredprecisionandthetimeoverwhichthesystemisevolved.
In Section 5 we demonstrate that the MSA method is superior compared to the
standardsplitstepmethod,sinceforthesameprecisionitismuchfaster. Moreover
for the split step method there is only an empirical estimate on the error while for
the MSA there is a rigorous bound. Using it in Section 6, we show that with the
helpofthismethodwecouldsolvetheSchr¨odingerequationforthemodelproblem
presented in Section 2, for an extremely long time; we conclude that for this model
theenergydoesnotgrowtoinfinityin-spiteofthetimedependentdrivingpotential.
The spreading in the quantum case is wider than in the corresponding classical
system. This results from the fact that initially we observe a lot of spreading in
the quantum case, while the classical case does not show spreading.
2. The model that will be studied
In this section we introduce the model for which our MSA method developed
in [1] is implemented. In the first subsection its relation to physical systems is
explained, while in the second one it is reduced to a form for which the MSA
method can be applied(see Eqs. 2.11-2.13).
2.1. The physical model. The random potentials which are prepared in optics
[18, 2] and atom optics [19, 20, 21] experiments, are described by a sum of random
Fourier components. In experiments, potentials which are composed out of a large
number of random independent Fourier components N, are created .
More specifically here, the Schr¨odinger equation for a potential
N
1 (cid:88)
(2.1) V (x(cid:48),τ)= √ A exp(i(k x(cid:48)−ω(cid:48)τ))+C.C
n n n
N
n=1
isused. TheA areindependent,identicallydistributedcomplexrandomvariables,
m
where C.C stands for complex conjugate. The expectation values of these variables
satisfy
(2.2) <A >=<A A >=0 <A A ∗ >=σ2δ
m m n m n nm
WewillstudythespecificmodelwhereAn =Aeiφm withφm uniformlydistributed
in the interval [−π,π] and A > 0. The distribution of ω and k is specified in
n n
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS3
Section 5.
The equation of motion is the time dependent Schr¨odinger equation
(2.3) i∂ ψ(x(cid:48),τ)=Hψ(x(cid:48),τ).
τ
Where
1
(2.4) H = p(cid:48)2+V(x(cid:48),τ),
2
p(cid:48) =−i∇ .
x(cid:48)
In the following section we will introduce the reduction of the problem to a form
where the MSA technique is applicable.
2.2. Reduction of the problem. In the model with potential given by (2.1), the
particle is expected to be accelerated most effectively in the regime of velocities
where the Chirikov resonances
ω(cid:48)
(2.5) v(r) = n
n k
n
are formed [22, 23].
Inthispaperwechoosev(r)andk uniformlydistributedintheintervals(cid:2)v(min),v(max)(cid:3)
n n
and (cid:2)k(min),k(max)(cid:3) respectively. The crucial point is that the Chirikov resonant
velocities are bounded in a phase space strip. This is typically the case for smooth
potentials, and will be assumed in this paper.
Here we would like to study what is the acceleration of a particle prepared with
momentum or velocity v (we assume unit mass), so that all v(r) are far from v.
n
For the classical corresponding system we found that the acceleration is negligible
[3, 5, 4]. Here we study the corresponding quantum mechanical system. For this
purpose it is convenient to work in a frame of reference where the initial velocity
of the particle vanishes. For this we perform the Galilean transformation
(2.6) x=x(cid:48)−vτ,
where v is the velocity of the moving frame (in later stages we will relate this
quantity to the required small parameter) on the potential (2.1)
N
A (cid:88)
(2.7) V (x+vτ,τ)= √ cos(k (x+vτ)−ω(cid:48)τ +φ )
n n n
N
n=1
N
A (cid:88)
= √ cos(k x+(k v−ω(cid:48))τ +φ ).
n n n n
N
n=1
Now, re-scaling time as
(2.8) t=τv,
the potential takes the form
(2.9) V (x,t)= √A (cid:88)N cos(cid:18)k x− (ωn(cid:48) −knv)t+φ (cid:19)
N n v n
n=1
and defining
(ω(cid:48) −k v)
(2.10) ω = n n
n v
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS4
The potential in the moving frame takes the form
N
A (cid:88)
(2.11) V (x,t)= √ cos(k x−ω t+φ )
n n n
N
n=1
The time dependent Schr¨odinger equation is
1
(2.12) i∂ ψ(x,t)= H(p,x,t)ψ(x,t)
t v
Introducing the small parameter β = 1:
v
(2.13) i∂ ψ(x,t)=βH(p,x,t)ψ(x,t)
t
From this point on we will use the small parameter β to perform the averaging
steps introduced in the next section. In what follows we will also relate the small
parameter β to the time scale on which we average. The Hamiltonian will be
approximated by a finite matrix (in space and momentum) and it will be verified
that the spreading never reaches the boundaries set by this basis.
3. The averaging scheme
The multiscale averaging method is based on replacing the original Hamiltonian
byahierarchicalsetofaveragedHamiltonians. Ineachstepweperforma“peel-off”
transformation and average a part of the Hamiltonian, for a chosen time interval
√
of length T0 = √1β = v. This choice is not unique, but as shown in [1] leads to
effective error bounds. We use the fact that Eq (2.13) is of the form of (2.1) in [1].
In this section the implementation of the MSA method of [1] for equation (2.13)
will be presented. In Appendix A we summarize the main results of Ref. [1]
3.1. Zero order. The zero order average on the jth time interval has the form
(V (t)≡V (x,t))
ˆ
1 (j+1)T0
(3.1) V¯(j) = V (t)dt
0 T
0 jT0
In the case of the potential (2.1), the zero order averaging can be performed ana-
lytically
ˆ ˆ
V¯(j) = 1 (j+1)T0V (t)dt= √A (j+1)T0(cid:88)N cos(k x−ω t+φ )=
0 T0 jT0 NT0 jT0 n=1 n n n
N
A (cid:88) 1
= √ [sin(k x−ω (j+1)T +φ )
NT ω n n 0 n
0 n=1
−sin(k x−ω jT +φ )]=
n n 0 n
N (cid:18) (cid:19) (cid:18) (cid:18) (cid:19) (cid:19)
A (cid:88) 1 1 1
(3.2) = 2√ sin ω T cos k x−ω j+ T +φ
NT0 n=1ωn 2 n 0 n n 2 0 n
This defines the Hamiltonian on one time interval; accordingly we can write the
Hamiltonian of one interval as
1 ∂2
(3.3) H¯(j)(x)=− +V¯(j)(x)
0 2∂x2 0
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS5
The global Hamiltonian corresponding to (2.1) of [1], which gives the zeroth order
approximation to H(t), is generated by:
(3.4) H¯g =H¯(j)(t) for jT ≤t<(j+1)T .
0 0 0 0
Using this notation, a general time evolution can be written as the product of the
interval propagators since H¯g is piecewise constant in time,
0
(3.5) U0(t)=e−iβH¯0(jmax)(t−jmaxT0)...e−iβH¯0(1)T0e−iβH¯0(0)T0
where j −1 is the integer part of t/T The evolution in this order is
max 0
(3.6) ψ (x,t)=U (t)ψ(x,0).
0 0
The propagator satisfies
∂
(3.7) i U (t)=βH¯(g)U (t)
∂t 0 0 0
corresponding to (2.10) of [1]. This propagator is numerically implemented and
used to solve the time dependent equation of motion. The error in this order is
bounded by β21 up to times of order T0, the reasoning behind this is shown in [1],
Eq. (2.49) there. The error in diagonalization of H¯(j) is negligible (see discussion
0
in Sec. 5)
3.2. First order. The first order averaging is based upon a “peel-off” transfor-
mation of the zero order. In such a transformation the next order Hamiltonian is
constructed from the zero order in the following way: Let H (t) be defined as
1
(3.8) H (t)=U−1(t)(cid:2)H(t)−H¯ g(t)(cid:3)U (t)
1 0 0 0
Hence H (t) is the Heisenberg dynamics of the full problem with U (t) dynamics
1 0
peeled off. An important point to note is that the Laplacian term drops in H !
1
Hence H (t) is a bounded operator, for which the results of [1] directly apply. Its
1
average in the j-th interval is
ˆ
1 (j+1)
(3.9) H¯(j) = H (t)dt
1 T 1
0 jT0
The H¯ (t) dynamics is given by the propagator
1
(3.10) U1(t)=e−iβH¯1(jmax)(t−jmaxT0)...e−iβH¯1(0)T0
√
WeturnnowtocalculateH¯(j),by(2.49)of[1]itisoforder β. Beforediagonalizing
1
thisoperatorinordertocalculatethetimeevolution,thereareseveralstepsneeded
to be taken. First we write explicitly H¯ (j) using integration by parts
1
(3.11)
ˆ
H¯ (j) = 1 (j+1)T0U−1(t)(cid:2)H(t)−H¯g(t)(cid:3)U (t)dt=H¯ (j,I)+H¯ (j,II)+H¯ (j,III),
1 T 0 0 0 1 1 1
0 jT0
where
ˆ
(3.12) H¯ (j,I) = 1 (cid:20)U−1(t)(cid:18) t(cid:2)H(t(cid:48))−H¯g(t(cid:48))(cid:3)dt(cid:48)(cid:19)U (t)(cid:21) (j+1)T0
1 T0 0 0 0 0 t=jT0
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS6
(3.13)
H¯ (j,II) =− 1 ˆ (n+1)T0dt(cid:48)(cid:18) ∂ U−1(t(cid:48))(cid:19)(cid:34)ˆ t(cid:48)(cid:2)H(s)−H¯g(s)(cid:3)ds(cid:35)U (t(cid:48)),
1 T ∂t(cid:48) 0 0 0
0 nT0 0
and
(3.14)
H¯ (j,III) =− 1 ˆ (j+1)T0U−1(t(cid:48))(cid:34)ˆ t(cid:48)(cid:2)H(s)−H¯g(s)(cid:3)ds(cid:35)(cid:18)∂ U (cid:19)(t(cid:48))dt(cid:48).
1 T 0 0 ∂t 0
0 jT0 0
We note that for any integer t/T ,
0
ˆ
(3.15) (j+1)T0(cid:2)H(t(cid:48))−H¯g(t(cid:48))(cid:3)dt(cid:48) =0
0
jT0
by construction. Therefore for such t H¯ (j,I)(t) = 0 and H¯ (j,II)(t) can be simpli-
1 1
fied. Moreover the expression in H¯ (j,III) is just the hermitian conjugate of H¯ (j,II)
1 1
so we only need to analyze the second term:
H¯ (j(,I3I.)16)= 1 ˆ (j+1)T0dt(cid:48)(cid:18) ∂ U−1(t(cid:48))(cid:19)(cid:34)ˆ t(cid:48) (cid:104)H (s)−H¯ (g)(s)(cid:105)ds(cid:35)U (t(cid:48))
1 T ∂t(cid:48) 0 j 0 0
0 jT0 jT0
= 1 ˆ (j+1)T0dt(cid:48)(cid:18) ∂ U−1(t(cid:48))(cid:19)(cid:34)ˆ t(cid:48) H (s)ds−(t(cid:48)−jT )H¯ (j)(cid:35)U (t(cid:48))
T ∂t(cid:48) 0 0 0 0 0
0 jT0 jT0
Toevaluatethisoperatorwecanusearecursiveconstructionofstatesbasedonthe
zeroth order eigenstates of the averaged Hamiltonian. It is important to remember
that this is actually a matrix. For convenience of notation we will define
(3.17) U0(t)=e−iβH¯0(j)(t−jT0)......e−iβH¯0(1)T0e−iβH¯0(0)T0,
For t that is an integer multiple of T .
0
j
(cid:89)
(3.18) U (t)= W
0 j−j(cid:48)
j(cid:48)=0
where
(3.19) Wj(t)=e−iβH¯0(j)T0.
First we assume t/T is integer and then generalize for any arbitrary t. The eigen-
0
values and eigenfunctions of H¯(j) are
0
(3.20) H¯(j)|ϕk(cid:105)=Ek|ϕk(cid:105)
0 j j j
Here we use a base of size M that is assumed to be finite (in this work we take
M = 64). Since β is small the eigenfunctions are approximately eigenfunctions of
a free particle, namely
(3.21) <x|ϕk >≈eik˜x.
j
(cid:12) (cid:12)
WewillhaveM eigenvaluesandeigenvectors,thelargestvalueof(cid:12)k˜(cid:12)is50. Between
(cid:12) (cid:12)
each pair of propagators W ,W we can insert the identity resolution in the
j i−1
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS7
corresponding basis (3.20)
M
(3.22) Iˆ =(cid:88)|ϕk(cid:105)(cid:104)ϕk|
j j j
k=1
resulting in
(3.23) U (t)=W ···W W =W Iˆ ···W I W I
0 j 1 0 j j 1 1 0 0
As an example let us take just the first two terms
M M
(cid:88) (cid:88)
(3.24) W I W I = W |ϕk1(cid:105)(cid:104)ϕk1| W |ϕk0(cid:105)(cid:104)ϕk0|
1 1 0 0 1 1 1 0 0 0
k1=1 k0=1
the brakets expressions (cid:104)ϕk1|(cid:80)M W |ϕk0(cid:105) are just scalars, and since ϕkj are
1 k0=1 0 0 j
eigenfunctions of Hj (3.20)
0
M M
(cid:88) (cid:88)
(3.25) W I W I = W |ϕk1(cid:105)(cid:104)ϕk1|W |ϕk0(cid:105)(cid:104)ϕk0|
1 1 0 0 1 1 1 0 0 0
k1=1 k0=1
M M
= (cid:88) e−iβT0E1k1 |ϕk1(cid:105)(cid:104)ϕk1| (cid:88) e−iβT0E0k0 |ϕk0(cid:105)(cid:104)ϕk0|
1 1 0 0
k1=1 k0=1
using the notation
(3.26) α =(cid:104)ϕk1|ϕk0(cid:105)
k1k0 1 0
this becomes
M M
(3.27) W1Iˆ1W0Iˆ0 = (cid:88) (cid:88) αk1k0e−iβT0E1k1e−iβT0E0k0 |ϕk11(cid:105)(cid:104)ϕk00|
k0=1k1=1
In the same way for the complete sequence of propagators (3.18)
M M
(cid:88) (cid:88)
U (t)= ··· α α ...
0 kjkj−1 kj−1kj−2
k0=1 kj=1
...αk1k0e−iβ(t−jT0)Ejkj ...e−iβT0E1k1e−iβT0E0k0 |ϕkjj(cid:105)(cid:104)ϕk00|.
(3.28)
and it satisfies (3.7).
Here jT <t<(j+1)T .
0 0
The inverse takes the form
(3.29)U−1(t) = eiβH¯0(0)T0eiβH¯0(1)T0......eiβH¯0(t−jT0)
0
M M
= (cid:88) ··· (cid:88) α∗ α∗ eiβ(t−jT0)E0eiβ(j−1)E0|ϕk0(cid:105)(cid:104)ϕkj|
k0k1 k1k2 0 j
kj=1 kj−1=1
By (3.7) or direct differentiation of (3.28)
M M M M
(3.30)∂ U (t) = − (cid:88) ··· (cid:88) (cid:88) ··· (cid:88) iβEkjα α ...α
∂t 0 j kjkj−1 kj−1kj−2 k1k0
k0(cid:48)=1 kj(cid:48)=1k0=1 kj=1
e−iβ(t−jT0)Ejkj ...e−iβT0E0k0 |ϕkj(cid:105)(cid:104)ϕk0|
j 0
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS8
and in a similar way one finds
(3.31)
M M
∂∂tU0−1(t)= (cid:88) ··· (cid:88) αk∗0k1αk∗1k2...αk∗j−1kjeiβT0E0k0 ...eiβ(t−jT0)EjkjiβEjkj|ϕk00(cid:105)(cid:104)ϕkjj|
kj=1 kj−1=1
substitution of (3.30) and (3.28) into (3.16) leads to
M M M M
H¯ (j,II)(t) = − 1 (cid:88) ··· (cid:88) (cid:88) ··· (cid:88) βEkjα∗ ...α∗ α ...α
1 T0 j k0(cid:48)k1(cid:48) kj(cid:48)−1kj(cid:48) kjkj−1 k1k0
ˆ k0(cid:48)=1 kj(cid:48)=1k0=1 kj=1
(j+1)T0dt(cid:48)eiβT0E0k0(cid:48)e−iβT0E0k0(cid:48) ...eiβ(t−jT0)Ejkj(cid:48)e−iβ(t(cid:48)−jT0)Ejkj(cid:48)
jT0 (cid:34)ˆ t(cid:48) (cid:35)
|ϕk0(cid:48)(cid:105)(cid:104)ϕkj(cid:48)| H(s)ds−(t(cid:48)−jT )H¯ (j) |ϕkj(cid:105)(cid:104)ϕk0|=
0 j 0 0 j 0
jT0
M M M M
= − 1 (cid:88) ··· (cid:88) (cid:88) ··· (cid:88) βEkjα∗ ...α∗ α ...α
T0 j k0(cid:48)k1(cid:48) kj(cid:48)−1kj(cid:48) kjkj−1 k1k0
ˆ k0(cid:48)=1 kj(cid:48)=1k0=1 kj=1
(j+1)T0dt(cid:48)e−iβT0E0k0(cid:48) ...eiβ(t−jT0)Ejkj(cid:48)e−iβ(t−jT0)Ejkj(cid:48)
(3.32) (cid:42)jϕT0kj(cid:12)(cid:12)ˆ t(cid:48) H(s)ds−(t−jT )H¯ (j)(cid:12)(cid:12)ϕkj(cid:48) (cid:43)|ϕk0(cid:48)(cid:105)(cid:104)ϕk0|
j (cid:12) 0 0 (cid:12) j 0 0
jT0
Finally the full expression for the first order averaged Hamiltonian is, where (3.11)
and (3.15) were used,
(cid:16) (cid:17)†
(3.33) H¯ (j)(t)=H¯ (j,II)(t)+ H¯ (j,II)(t) .
1 1 1

H¯ (j)(t) = 2(cid:60)− 1 (cid:88)M ··· (cid:88)M (cid:88)M ··· (cid:88)M βEkjα∗ ...α∗
1  T0 j k0(cid:48)k1(cid:48) kj(cid:48)−1kj(cid:48)
k0(cid:48)=1 ˆkj(cid:48)=1k0=1 kj=1
(3.34) αkjkj−1...αk1k0 (j+1)T0dt(cid:48)e−iβT0E0k0(cid:48) ...eiβ(t−jT0)Ejkj(cid:48)e−iβ(t−jT0)Ejkj(cid:48)
(cid:42)ϕkj(cid:12)(cid:12)ˆ t(cid:48) dsH(sj)T0−(t−jT )H¯ (j)(cid:12)(cid:12)ϕkj(cid:48) (cid:43)(cid:41)|ϕk0(cid:48)(cid:105)(cid:104)ϕk0|
j (cid:12) 0 0 (cid:12) j 0 0
jT0
The integral (3.34) is performed numerically as follows. The domain of integra-
tion is derived into squares of size ∆t×∆t. The integral is approximated by a sum
ofthevaluesoftheintegrandofthemiddlepointsofthesquares,multipliedby∆t2.
The error in each term is of the order of ∆t2. This can be improved substantially.
Hereitisnotrequiredsincewecanobtaintherequiredprecisioninthissimpleway.
At first sight this expression (3.34) might seem very complicated but can be
understood quite easily; in fact what we have here is a matrix constructed from
the sum of matrix products; the terms of the matrix involve only the products of
eigenvalues and eigenvectors of H¯(j). The benefit of calculating the propagator in
0
this manner is that one only needs to diagonalize the Hamiltonian in each interval
√
only once. As a result of (3.2) the averaged potential is of order β and α are
kj,ki
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEMS9
elements of matrices that are almost diagonal as will be verified aposteriori in Sec.
5. Consequently the error in the calculation of (3.34) is of the order
(3.35) δI =∆t2
t
By the general theory see (2.15) of [1] the first order propagator takes the form
(3.5) with H¯(j) replaced by H(j) namely with
0 1
U1 =e−iβH¯1(j)(t−jT0)......e−iβH¯1(1)T0e−iβH¯1(0)T0
3.3. Normal form transformation. To improve the accuracy we implement the
normal form transformation; we will use the form given in Eq. (3.10) of [1]. In our
notation this becomes
ˆ
(cid:18) t (cid:104) (cid:105)(cid:19)
(3.36) U˜ (t)=1+iβU−1(t) dt(cid:48) H (t(cid:48))−H¯(g)(t(cid:48)) U (t).
1 1 1 1 1
0
The manner in which we will simplify the above expression will be similar to the
method used to calculate U and H given in Eq (3.10) and (20), splitting into
1 1
intervals of length T and using (3.15), replacing H by H :
0 0 1
 ˆ
(3.37) U˜1(t) = 1+iβU1−1(cid:88)j (j(cid:48)+1)T0dt(cid:48)(cid:20)H1(t(cid:48))−H¯1(j(cid:48))(cid:21)
j(cid:48)=0 j(cid:48)T0
ˆ
t (cid:104) (cid:105) (cid:19)
+ H (t(cid:48))+H¯(j)(t(cid:48)) dt(cid:48) U (t).
1 1 1
jT0
by definition of H¯(j)
1
ˆ
(j(cid:48)+1)T0dt(cid:48)(cid:20)H (t(cid:48))−H¯(j(cid:48))(cid:21)=0
1 1
j(cid:48)T0
ˆ
t (cid:104) (cid:105)
(3.38) U˜ (t) = 1+iβU−1 dt(cid:48) H (t(cid:48))−H¯(j)(t(cid:48)) U (t)=
1 1 1 1 1
jTˆ0
(cid:20)(cid:18) t (cid:19) (cid:21)
1+iβU−1 H (t(cid:48))dt(cid:48) −(t−jT )H¯(j) U (t)
1 1 0 1 1
jT0
and explicitly
(3.39)
ˆ
(cid:20)(cid:18) t (cid:19) (cid:21)
U˜(t)=1+iβeiβH¯1(0)T0...eiβH¯1(j)(t−jT0) H1(t(cid:48))dt(cid:48) −(t−jT0)H¯1(j) U1(t)
jT0
U˜(t)=e−iβH¯1(j)(t−jT0)...e−iβH¯1(0)T0.
If also the normal form transformation is performed the error is of order β23 .
4. Iterative application and error analysis
To reduce the error we introduce an iterative process. After the normal form
transformation is performed we introduce a new Hamiltonian
1 (cid:104) (cid:105)
(4.1) H˜ = √ U˜−1 H −H(g) U˜
1 β 1 1 1 1
where
(4.2) U˜ =U U U˜
1 0 1
NUMERICAL IMPLEMENTATION OF THE MULTISCALE AND AVERAGING METHODS FOR QUASI PERIODIC SYSTEM1S0
ThewavefunctionatthenewlevelsatisfiesaSchr¨odingerequationlike(3.12)of[1]
with H playing the role of A˜(t) there, and β is replaced by β3/2. Now one starts
1
from an equation like (2.13) with β replaced by β3/2 and H by H˜ . At each step
1
the effective value of β is reduced (cid:0)β →β3/2(cid:1), and so is the error. The process is
repeated until the bound on the error is satisfactory, as will be shown below.
Assume the process repeated l times. The resulting approximation for the wave
function is
(4.3) ψ(t)=U˜ ..........U˜ U˜ ψ(0)
1 l−1 l
where U˜ is the propagator at the l(cid:48) level of the hierarchy, corresponding to the
l(cid:48)
Hamiltonian H˜ .
l
Weturnnowtoestimaterigorouslytheerrorsusingthemulti-scaleandaveraging
method, assuming l-levels of the hierarchy.
With β small on a time interval of order β12 ≡ T0. Let us denote by tmax the
maximum time we want to simulate dynamics. After introducing a normal form
transformation we eliminate the cβ12 error term on (3.10) of [1], and we get for the
evolution
(4.4) ψ(t)=U(t)ψ(0)=U (t)R (t)ψ(0)
0 0
leading to
(4.5) U(t)=U (t)R (t) for 0≤t≤t
0 0 max
(cid:16) (cid:17)
withR0(t)=1+O β32tmax for0≤t≤tmax withU0(t)isgivenby(aproductof)
averaged dynamics followed by a normal form unitary transformation. Moreover
∂R
(4.6) i ∂t0 =β32H1(t)R0(t) ,
with
(4.7) (cid:107)H (t)(cid:107)=O(1).
1
If we then use the same method of averaging on R (t) in Eq. (4.6) we get R (t)=
0 0
U (t)R (t) and
1 1
(4.8) U(t)=U˜ (t)U (t)R (t)
0 1 1
where now U (t) is the averaged approximate solution for R (t), and
1 0
(cid:16) (cid:17) (cid:18)(cid:16) (cid:17)3 (cid:19)
(4.9) R1(t)=1+O β23 +O β32 2 tmax .
(cid:16) (cid:17)
Againafternormalformtransformation, theO β23 correctiondropsandwehave
(cid:18)(cid:16) (cid:17)3 (cid:19)
(4.10) U(t)=U˜0(t)U˜1(t)+O β23 2 tmax .
After l-such iterations, we get the exact solution
(4.11) U(t)=U˜0(t)U˜1(t)···U˜l(t)+O(cid:16)β(23)ltmax(cid:17).
So the convergence of the scheme is super exponentially fast, close to the Newton
type iteration.