Many-body interband tunneling as a witness for complex dynamics in the
  Bose-Hubbard model
Many-body interband tunneling as a witness for complex dynamics in the
Bose-Hubbard model
Andrea Tomadin,1 Riccardo Mannella,1 and Sandro Wimberger1,2
Dipartimento di Fisica, Unversità degli Studi di Pisa, Largo Pontecorvo 3, 56127 Pisa, Italy
CNISM, Dipartimento di Fisica del Politecnico, C. Duca degli Abruzzi 24, 10129 Torino, Italy
(Dated: November 4, 2018)
A perturbative model is studied for the tunneling of many-particle states from the ground band
to the first excited energy band, mimicking Landau-Zener decay for ultracold, spinless atoms in
quasi-one dimensional optical lattices subjected to a tunable tilting force. The distributions of the
computed tunneling rates provide an independent and experimentally accessible signature of the
regular-chaotic transition in the strongly correlated many-body dynamics of the ground band.
PACS numbers: 03.65.Xp,32.80.Pj,05.45.Mt,71.35.Lk
The experimental advances in atom and quantum op-
tics allow the experimentalist to directly study a plethora
of minimal models which have been developed to de-
scribe usually much more complex phenomena occurring
in solid states [1, 2, 3]. Bose-Einstein condensates loaded
into optical lattices, which perfectly realize spatially peri-
odic potentials, are used, e.g., to implement the Wannier-
Stark problem [4, 5, 6] as a paradigm of quantum trans-
port where atoms move in a tilted lattice. Up to now all
experiments on the Wannier-Stark system with ultracold
atoms have been performed in a regime where atom-atom
interactions are either negligible [4] or reduce to an effec-
tive mean-field description [5, 7]. State-of-the-art setups
are, however, capable to achieve small filling factors of
the order of one atom per lattice site [2]. Moreover, the
atom-atom interactions can be tuned by the transversal
confinement and by Feshbach resonances [3, 8], resulting
in strong interaction-induced correlations.
The regime of strong correlations in the Wannier-Stark
system was addressed in [9, 10], revealing the sensitive
dependence of the system’s dynamics on the Stark force
F . The single-band Bose-Hubbard model of [9, 10] is
defined by the following Hamiltonian with the creation
l,1, annihilation âl,1, and number operators n̂ l,1 for the
first band of a lattice l = 1 . . . L:
Fl n̂ l,1−
â l+1,1
† â l,1 + h.c.
n̂ l,1 ( n̂ l,1 − 1) .
A transition from a regular dynamical (dominated by F )
to a quantum chaotic regime (with comparable values of
J1, U1 , F ) was found [9, 10]. The transition was quan-
titatively studied using the distribution of the spacings
between next nearest eigenenergies of the Hamiltonian
(1). This analysis [9, 10] verifies that the normalized level
spacings s ≡ ∆E/∆E obey a Poisson (P(s) = exp(−s) )
and a Wigner-Dyson (WD: P(s) = sπ/2 exp(−πs2/4) )
distribution in the regular and chaotic case, respectively
[11]. P(s) and the cumulative distribution functions
(CDF: C(s) ≡
ds′ P(s′)) are shown for typical cases in
Fig. 1, where we scanned F to emphasize the crossover
0 1 2 3 4 5
1.8 2 2.2 2.4 2.6 2.8 3
(2π/F)
0 1 2 3 4 5
0 1 2 3 4 5
0 1 2 3 4 5
0 5 10 15 20 25
FIG. 1: (a,b) CDF (stairs) and P(s) (stairs in insets) for N =
5 atoms, L = 8, lattice depth V = 10 recoil energies (fixing
J1 = 0.038), U1 = 0.032, F ≃ 0.063 (a) and 0.021 (b), with
WD (solid) and Poisson distributions (dashed). (c) χ2 test
with values close to zero for good WD statistics. The dashed
line marks the transition to quantum chaos as F is tuned. (d)
variance of the number of levels in intervals of length dE (with
normalized mean spacing), for the cases of (a) (squares) and
(b) (circles), with the random matrix predictions for Poisson
(dashed) and WD (solid) [11].
between the regular and the chaotic regime. Statistical
tests are also shown which confirm the analysis of [9, 10]
in a more systematical manner [12].
As shown in [9], the strong correlations in the quantum
chaotic regime induce a fast and irreversible decay of the
Bloch oscillations, which otherwise would persist in the
ideal, non-interacting case. Therefore, the crossover be-
tween the two regimes discussed above could be measured
in experiments by observing just the mean momentum as
a function of time. Here we introduce a new, robust and
hence also experimentally accessible prediction for this
crossover. In the presence of strong interactions parame-
terized by U1 , the single-band model should be extended
to allow for interband transitions [13], as recently realized
http://arxiv.org/abs/0704.0233v1
at F = 0 in experiments with fermionic interacting atoms
[3]. Instead of using a numerically hardly tractable com-
plete many-bands model, we introduce a perturbative de-
cay of the many-particles modes in the ground band to
a second energy band. Our novel approach to study the
Landau-Zener-like tunneling between the first and the
second band [1, 5, 7, 14, 15] leads to predictions for the
expected decay rates and their statistical distributions.
As we will show, the latter are drastically affected by the
dynamics in the ground band, and they therefore provide
a measurable witness for the regular-chaotic transition.
We first derive the individual decay rates of the dom-
inating interband coupling channels. These decay rates
will serve to effectively open the single-band model (1)
for mimicking losses arising from the interband coupling.
Our analysis starts from the following “unperturbed”
Hamiltonian for the first two bands:
ε1 n̂ l,1 + ε2 n̂ l,2 − J22 ( â l+1,2
† â l,2 + h.c.)
+Fl ( n̂ l,1 + n̂ l,2) +
n̂ l,1 ( n̂ l,1 − 1)
. (2)
For a moment, we neglect the hopping in the lower band,
where the single-particle Wannier functions [14] are more
localized than in the upper band. In the latter we neglect
the interactions, since initially only a few particles pop-
ulate the excited levels. A closer analysis of the full two-
bands system [12] shows that there are two dominating
mechanisms that promote particles to the second band.
The first one is a single-particle dipole coupling arising
from the force term:
H1 = F · D
â l,2
â l,1 + â l,1
â l,2
, (3)
where D depends only on the lattice depth V (measured
in recoil energies according to the definition in [7]). The
second one is a many-body effect, describing two particles
of the first band entering the second band together:
â l,2
â l,2
â l,1 â l,1 + (1 ↔ 2)
. (4)
The cross-band interaction is characterized by the pa-
rameter U× ≡ ãs
dxχ21χ
2 ≃ 0.5U1 (for V = 3 . . . 10)
[12], for U1 ≡ ãs
dxχ41, with renormalized scattering
length ãs [8, 12] and the Wannier functions χ1,2 local-
ized in each well for the first or second band. To justify
the following perturbative approach, it is crucial to real-
ize that the terms (3) and (4) must be small compared
with the band gap ∆ ≡ ε2 − ε1 (not necessarily small
with respect to the single band terms in (1)), and indeed
FD,U×, U1 ≪ ∆ for the parameters considered here.
For the first perturbation, the decay channel of a given
unperturbed Fock state labelled |b〉 (with a total number
of atoms N and nh atoms in an arbitrary well h) is
|b;N〉 ⊗ |vac〉 → |b′;N − 1〉 ⊗ |w〉 , n′h = nh − 1. (5)
Here, |w〉 =
m=−∞ Jm−w(|J2|/F ) âm,2
†|vac〉 is the
single-particle eigenstate for the Wannier-Stark problem,
localized around the site w in the second band, with the
Bessel function of the first kind Jm(x) [14].
The expectation value of (3) for |b;N〉 of the first band,
equal to the first-order δE(b), is zero because the opera-
tor does not conserve the number of particles within the
bands. The decay width at first-order is given by the ma-
trix element of the perturbation between the initial and
final state according to Fermi’s Golden Rule, and only
the first term in (3) gives a nonzero contribution [12]:
〈k|〈b′|
l=1 â l,2
† â l,1 |b〉|vac〉 =
l=1 Jl−w(|J2|/F )
·δ(n′
, nl − 1)
m 6=l δ(n
m, nm). (6)
The δ(·, ·) functions act as a selection rule for the Fock
states that are coupled by the perturbation. The tun-
neling mechanism does not include any income of en-
ergy from an external source, so the initial and final
energies E0(b) = 〈vac|〈b|H0|b〉|vac〉 and E0(b′, w) =
〈w|〈b′|H0|b′〉|w〉, respectively, must be equal as required
by the Golden Rule. The condition on the energy conser-
vation is, however, relaxed to account for the uncertainty
∆E(b) of the unperturbed energy levels of the initial and
final states in the lower band arising from the hopping in
this band initially neglected in (2). A detailed derivation
is given in [12], and here we only state the result:
∆E(b) = 2π (J1 /2)
∆E(b → b′) =
2π (J1 /2)
∆l=±1
n2l δ(nl+∆l + 1, nl). (7)
The level density ρ(E, b) around the unperturbed en-
ergy E0(b) of a Fock state |b〉 is then approximated by
a rectangular profile, of width ∆E(b) and unit area:
ρ(E, b) = χ {|E − E0(b)| ≤ ∆E(b)/2} / ∆E(b). The re-
laxed energy conservation rule selects from (5) the set K
of permitted decay channels (h,w) parameterized by the
two indices h,w such that:
′, w)− E0(b) = ∆− F (h− w)− U1 (nh − 1)
−∆E(b)+∆E(b
∆E(b)+∆E(b′)
. (8)
Hence the energy ∆ required to promote a particle to the
second band is supplied by the decrease of the interaction
(∝ U1 ) and by the work of the force (∝ F ) exerted on
the promoted particle.
The total width Γ1(b) for the decay via the allowed
channels K, is proportional to the square of the matrix
element and to the level density ρ(E, b):
Γ1(b) = 2π(FD)
(h,w)∈K
Jh−w( |J2|F ) ·
∆E(b)∆E(b′)
. (9)
Jm(x) significantly contributes only for |m| <∼ |x|. If
U1 ,∆E(b) ≪ ∆, the energy conservation is roughly
given by |∆| ≃ F (h − w). Requiring that the Bessel
function in (9) is substantially larger than zero, we ob-
tain the inequality |∆| ≤ |J2|. The last condition does
not depend on F , since a twofold effect is at work: a
stronger force produces a larger energy gain when a par-
ticle moves along the lattice, but the extension |J2/F | of
the single-particle state shrinks. Therefore, increasing F
results in an increased energy matching and a strongly
reduced “geometrical” matching. For 3 < V < 26, we
have |∆|− |J2| > 1.0 [12], such that the energy matching
cannot be realized by just tuning the lattice depth. The
decay can, however, be activated by an increase of the in-
teractions, which can be experimentally achieved by act-
ing on the transversal confining potential of a quasi-one
dimensional lattice, or by a Feshbach resonance [8]. In
the calculations presented below, we augmented U1 used
in [9, 10] by a factor of order 10, and as noted in the in-
troduction, a similar increase of the interaction strength
was used in the experiment to promote fermions to higher
bands [3], in close analogy to the here described field- and
interaction-induced interband coupling of bosons.
The second term (4) is treated in a similar way, with
the difference that two particles are promoted to the sec-
ond band, and the position of the second single-particle
state |w′〉 is an additional degree of freedom for the tran-
sition. The decay channels are:
|b,N〉⊗ |vac〉 → |b′, N − 2〉 ⊗ |w,w′〉 ; n′h = nh − 2. (10)
The energy matching selects a set K of decay channels,
parameterized by the three site indices h,w,w′:
(h,w,w′) ∈ K s.t. E0(b′, w, w′)− E0(b) =
= 2∆− F (2h− w − w′)− U1 (2nh − 3)
∆E(b) + ∆E(b′)
∆E(b) + ∆E(b′)
. (11)
The computation of the matrix element yields [12]:
Γ2(b) = 2π
(h,w,w′)∈K
Jh−w( |J2|F ) ·
Jh−w′( |J2|F )
· 4nh (nh − 1) 1∆E(b)∆E(b′)
. (12)
With respect to (9), the additional degree of freedom
w′ results in a summation over all possible values of
w−w′. This follows from the possibility to conserve the
energy even if a particle is pushed far, if the other parti-
cle is pushed almost equally far in the opposite direction.
Since the decay widths in (12) depend on the product of
two (rapidly decaying) Bessel functions – again a “geo-
metrical” matching condition – we apply the truncation
|w−w′| ≤ |J2/F |, to reduce the formula to a finite form.
We can now compute the total width ΓF(b) = Γ1(b) +
Γ2(b) defined by the two analyzed coupling processes for
0 1 2 3 4 5
-6.5 -6 -5.5
0 1 2 3 4 5
-5 -4.5 -4 -3.5
Log10Γ
0 1 2 3 4 5
-4 -3.5 -3 -2.5
-6.5 -6
-4 -3.5
FIG. 2: (a,c,e) CDF from Re {Ej} (stairs), together with
WD (solid) and Poisson predictions (dashed). (b,d,f) distri-
butions of the logarithm of the rates. In (a,b), (c,d), (e,f)
F ≃ 0.17, 0.31, 0.47, respectively, with (N,L) = (7, 6), V =
3, U1 = 0.2 (fixing U× ≃ 0.1). In the regular regime (f),
a log-normal distribution (dotted) well fits the data, with a
scaling P(Γ) ∝ Γ−x for the largest Γ (dashed line in the inset
of (f) with x = 1). In the chaotic case, a global power-law
behavior with x ≈ 2 is found (dashed line in the inset of (b)).
each basis state |b〉 of the single-band problem given in
(1). The ΓF(b) are inserted as complex potentials in the
diagonal of the single-band Hamiltonian matrix. After
a gauge transform that recovers the translational invari-
ance of the problem (see [10, 12] for details), the latter
matrix is used to compute the evolution operator over one
Bloch period TB, which is finally diagonalized to obtain
its eigenphases exp (−iEj TB). Along with the statistics
of the level spacings defined by Re {Ej}, the Figs. 2 and 3
analyze the statistical distributions of the tunneling rates
Γj = −2Im{Ej} for some paradigmatic cases. All rates
are much smaller than unity, which a posteriori is fully
consistent with our perturbative approach.
To observe what happens at the regular-chaotic tran-
sition (c.f. Fig. 1), we scan F in Fig. 2, and as F in-
creases, the average decay increases by orders of mag-
nitude, while the distributions broaden. The large in-
crease of the rates is due to an improved energy match-
ing, when F supplies the necessary energy to promote
particles to the second band. For the parameters of
Fig. 2, the single-particle Landau-Zener formula [14] gives
ΓLZ = F/(2π) exp
−π2∆2/(8F )
∼ 10−23, 10−12, 10−8
for (b,d,f). This huge variation, typical of semiclassical
formulae, implies that there are possibly parameters for
which our results are comparable to the single-particle
prediction, but, in general, the many-particle effects can-
not be neglected. Moreover, mean-field treatments of
the Landau-Zener tunneling at best predict a shift of Γ
[7, 15], but cannot account for their distributions.
In the chaotic regime, the Fock states are strongly
mixed by the dynamics [9, 10, 12] and a fast decaying
-12 -11 -10 -9 -8
-5.5 -5 -4.5 -4
10a) b)
c) d)
FIG. 3: (a,c) rate distributions in the chaotic regime with F ≃
0.17, U1 = 0.2 (U× ≃ 0.1), together with the corresponding
unscaled P(Γ) in (b,d). In (a,b) (N,L) = (7, 6), V = 4, and
in (c,d) (N,L) = (9, 8), V = 3. Power-laws P(Γ) ∝ Γ−x are
found with x ≈ 2 (dashed lines in (b,d)).
Fock state can act as a privileged decay channel for many
eigenstates. Many states then share similar rates, lead-
ing to thinner distributions. Therefore, the thinner dis-
tribution of Fig. 2 (b) is a direct signature of the chaotic
dynamics evidenced in (a), as compared with the regular
case in (e,f). In Fig. 2 (f), we found a good agreement
with the expected log-normal distribution of decay rates
[16] (or of the similarly behaving conductance [17]) in
the regular regime. There the system shows nearly per-
fect Bloch oscillations [9], and the motion of the atoms
is localized along the lattice [14]. We can even detect a
qualitative crossover to a power-law P(Γ) ∝ Γ−1 in the
right tail of the distribution, as predicted from localiza-
tion theory [16, 18, 19]. The distributions in Figs. 2 (b)
and 3 follow the expected power-law for open quantum
chaotic systems in the diffusive regime [18]. The expo-
nents x are, however, nonuniversal and depend on the
opening of the system. In our case, the decay channels
are defined by the interband coupling, which in a sense
attaches “leads” to all lattice sites within the sample.
Going along with the regular-to-chaotic transition in the
lower band of our model (from Fig. 2 (f) to (b), or to
Fig. 3) the Γ distributions transform from a log-normal
to a power-law with x ≈ 2, in close analogy to the tran-
sition from Anderson-localized to diffusive dynamics in
open disordered systems [18, 20].
In summary, our perturbative opening of the single-
band Wannier-Stark system allows one to study Landau-
Zener-like interband tunneling within a many-body de-
scription of the dynamics of ultracold atoms. The statis-
tical characterization of the tunneling rates (mean values
and form of the distributions) provides clear and robust
signatures of the regular-to-chaotic transition for future
experiments. A more detailed analysis of the interband
coupling in a full-blown model, in which at least two
bands are completely included, calls for huge computa-
tional resources to access the complete quantum spec-
tra. Nonetheless, our results are a first step in the di-
rection of studies for which “horizontal” and “vertical”
quantum transport along the lattice are simultaneously
present and influence each other in a complex manner.
We thank the Centro di Calcolo, Dipartimento di
Fisica, Università di Pisa, for providing CPU, and the
Humboldt Foundation, MIUR-PRIN, and EU-OLAQUI
for support.
[1] O. Morsch et al., Rev. Mod. Phys. 78, 179 (2006).
[2] M. Greiner et al., Nature 415, 39 (2002); T. Stöferle et
al., Phys. Rev. Lett. 92, 130403 (2004); S. Fölling et al.,
ibid. 97, 060403 (2006).
[3] M. Köhl et al., Phys. Rev. Lett. 94, 080403 (2005).
[4] M. BenDahan et al., Phys. Rev. Lett. 76, 4508 (1996);
S.R. Wilkinson et al., ibid. 4512; B.P. Anderson et al.,
Science 282, 1686 (1998).
[5] O. Morsch et al., Phys. Rev. Lett. 87, 140402 (2001).
[6] G. Roati et al., Phys. Rev. Lett. 92, 230402 (2004).
[7] S. Wimberger et al., Phys. Rev. A 72, 063610 (2005).
[8] T. Bergeman et al., Phys. Rev. Lett. 91, 163201 (2003).
[9] A. Buchleitner et al., Phys. Rev. Lett. 91, 253002 (2003).
[10] A.R. Kolovsky et al., Phys. Rev. E 68, 056213 (2003).
[11] M.L. Mehta, Random Matrices and the Statistical Theory
of Energy Levels (Academic Press, New York, 1991).
[12] A. Tomadin, Master’s thesis, Università di Pisa (2006).
[13] V.W. Scarola et al., Phys. Rev. Lett. 95, 033003 (2005).
[14] M. Glück et al., Phys. Rep. 366, 103 (2002).
[15] B. Wu et al., Phys. Rev. A 61, 023402 (2000); O. Zobay
et al., ibid. 61, 033603 (2000); D. Witthaut et al., ibid.
75, 013617 (2007); S. Wimberger et al., J. Phys. B 39,
729 (2006).
[16] M. Terraneo et al., Eur. Phys. J. B 18, 303 (2000).
[17] C.W.J. Beenakker, Rev. Mod. Phys. 69, 731 (1997).
[18] T. Kottos, J. Phys. A 38, 10761 (2005).
[19] G. Casati et al., Phys. Rev. Lett. 82, 524 (1999); S. Wim-
berger et al., ibid. 89, 263601 (2002).
[20] S. Wimberger et al., J. Phys. A 34, 7181 (2001).
