Collective tunnelling of strongly interacting cold atoms in a double well potential
D. Sokolovski X. Gutiérrez de la Cal M. Pons






⟨
  

 
⟩
ℰ
ϵ
τ

Ĥ









^1Departmento de Química-Física, Universidad del País Vasco, UPV/EHU, Leioa, Spain
^2 IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, 48013, Bilbao, Spain

^3Departmento de Física Aplicada I, Universidad del País Vasco, UPV/EHU, Bilbao, Spain

It is known that under resonance conditions, a group of strongly interacting bosonic atoms, trapped in a double-well potential, 
mimics a single particle, performing Rabi oscillations between the wells.
By implication, all atoms need to tunnel at roughly the same time, even though 
the Bose-Hubbard Hamiltonian accounts only for one-atom-at-a-time transfers.
We analyse the mechanism of this collective behaviour,
evaluate the Rabi frequencies in the process, and discuss the limitation of this simple picture.
In particular, it is shown that the small rapid oscillations superimposed 
on the slow Rabi cycle, result from splitting the transferred cluster at the sudden onset 
of tunnelling, and disappear if tunnelling is turned on gradually.






03.65.-w, 03.65.Yz, 03.75.Nt
Collective tunnelling of strongly interacting cold atoms in a double well potential
    D.  Sokolovski^1,2, X. Gutieŕrez de la Cal^2, and  M. Pons^3
    December 30, 2023
===================================================================================






§ INTRODUCTION

For a quantum particle, which can occupy several quantum states with roughly the same energy, 
even a small perturbation is capable of causing transitions between the states.
If only a pair of states satisfies this resonance condition, these transitions would, in general, be in the form Rabi oscillations <cit.>, 
moving the particle periodically from one state to the other. 
A more interesting case is the one where the resonance occurs between many-body states of interacting particles, 
which correspond to different spatial configurations of the system.
A direct matrix element may or may not connect a pair of such states, say, |Ψ_1 and |Ψ_2. In the latter case the system
would need to reach the final state via a pathway, passing through the states of the system, whose energies can lie far from the resonance. The passage must be, in some sense, rapid, since a measurement will almost 
always find the system either in |Ψ_1, or in |Ψ_2, and only rarely detect it elsewhere.



Examples of collective behaviour can be found, for example, in cold atom physics.
In <cit.> the authors observed coherent many-body Rabi oscillations in 
electronics transitions of interacting rubidium atoms. More relevant to our analysis 
is the direct observation of 
 
correlated tunnelling of pairs of strongly 
interacting cold atoms,  reported in <cit.>. 


Recently developed laser techniques are capable of trapping
such atoms in quasi one-dimensional traps <cit.>, and the dynamics of interacting atomic systems has been extensively studied both by solving the Schrödinger equation numerically (see, e.g. <cit.>, <cit.>, <cit.>), and by using a simplified Bose-Hubbard model,

as was done, for example, in <cit.>, <cit.>, <cit.>. 
In <cit.>, <cit.> and <cit.> experimental and theoretical analysis of  larger systems are presented, and the contributions from higher bands are included. In these systems, effects such as breathing are cradle modes can be observed as a consequence of this high-bands contributions, and tunnelling between wells as an effect of the lower modes. 
Similar dynamics may also be observed, for example, using mixtures of atoms, where different inter-species interaction regimes were studied
 
 <cit.>, 
 or spin chains <cit.>. Other theoretical studies, with a large number of particles,  address the dynamics of a BEC in a double-well using the Gross-Pitaevskii equations including many-body interactions studying the self-trapping effect induced <cit.> or different regimes, from coherent oscillations to their suppression when the number of bosons is high <cit.>. Important experimental results have also been obtained for BEC in double-well potentials, such as the first realization of a single Josephson junction <cit.> or the more recent <cit.> where the dynamical control of correlated tunnelling processes of strongly interacting particles is presented.

The tunnelling frequencies for a symmetric trap were first evaluated in <cit.>, where the authors relied on the time independent perturbation theory of <cit.>, in order to obtain energy splitting between the resonant states.
In their follow up paper <cit.>, the authors of <cit.> studied time evolution of the average difference of the wells' populations,
δ n(t), for various initial conditions, and  analysed the frequency spectrum of quantum fluctuations, superimposed on the Rabi cycle. 
However, this is not the whole story, and certain aspects of the collective tunnelling phenomenon require a further discussion.
In particular, in the Rabi oscillations one would  find all the transferred atoms in the same well, at all times. 
This suggests that the atoms must tunnel together, almost instantaneously, or at least during a time much shorter than the Rabi period  <cit.>.
Neither the analysis of <cit.>, <cit.>, nor the form of the Bose-Hubbard Hamiltonian, which contains only single atom transfer terms, give an immediate clue as to how this may be possible. 
A study of the  mechanism of this rapid collective transfer, and identification of the relevant time scales is the first of our aims.
Furthermore, the picture in which a cluster of transferred atoms behaves as a single particle, performing Rabi oscillations, is only approximate,
and deserves further attention. The mean populations difference <cit.>, <cit.>,  δ n(t) = ∑_n=0^N[ p(n,t)-(N-n)(1-p(n,t)], where p(n,t) are the probabilities for finding n out of N bosons in the right well at a time t, is a rather crude averaged quantity, and may not be best suited for such an analysis.
Quantum fluctuations, evident in the δ n(t) <cit.>, come from the directly measurable individual probabilities p_n(t).
These, in turn, are absolute squares of the sums of the probability amplitudes, corresponding to elementary processes, such as transfer of a single atom  from one well to the other. Identification of interfering scenarios, responsible for the additional oscillatory patterns, specific to many body Rabi oscillations, is the second main subject of this paper. 
Fortunately, with the tunnelling matrix element small, and the Rabi period large, the required analysis can be carried out already in the first non-vanishing order of the time dependent perturbation theory. With its help, we will show that the largest contribution to the probability of the resonance transfer of n atoms will come from a process, in which all n bosons "jump"  together, 
almost instantly if compared to the Rabi period. 
We will also demonstrate that the additional oscillations  result from the processes, in which one or more atoms are split from the tunnelling cluster.  This breakup of the cluster can be related to a sudden onset of tunnelling at the start of the experiment. Our prediction that the oscillations will be quenched if tunnelling is turned on slowly, (compared to the characteristic "jump time"), can be subjected to experimental verification.
The rest of the paper is organised as follows. In Section II we formulate many-body resonance conditions  
for atoms in an asymmetric trap. In Sec. III we use time-dependent perturbation theory to analyse 
a rapid transfer of a group of atoms along an indirect pathway, connecting the resonance states, 
and obtain expressions for the Rabi frequencies.  In Sec. IV we relate additional oscillatory patterns, seen in the Rabi probabilities,
to sudden switching of the tunnelling. In Sec. V we briefly consider a detuned regime, and Section VI contains our conclusions.
A further detailed discussion of the interference mechanism of the transfer can be found in Appendix A. Our  findings are tested in the Appendix B on some of the exactly solvable cases.



§ TRAPPED ATOMS IN THE BOSE-HUBBARD APPROXIMATION


We consider  N strongly interacting identical bosons, contained  in an asymmetric double-well
potential, as shown in Fig.1. An experimental realisation of such a system can be achieved, e.g.,  using ^87Rb atoms (see, for example <cit.>).
The energy of the system is modelled by the Bose-Hubbard Hamiltonian  (we use ħ=1)


    Ĥ_(t)= U(c_L^+)^2(c_L)^2/2+U(c_R^+)^2(c_R)^2/2+β c_R^+c_R
    +Ω (c_L^+c_R+c_R^+c_L),

where c_L^+(c_L) and c_R^+(c_R) create (annihilate)
a boson in the left (right) state in Fig. 1.

In the first two terms in Eq.(<ref>), we take  U to describe short range interactions between the bosons in the same well,  which is neglected for atoms placed on different sides of the barrier,

β is the difference 
between one-particle energies on the right and on the left, and Ω is the tunnelling amplitude,
which allows for a transfer of a particle between the wells.
Note that both U and β can, in principle, be positive or negative, although in what follows 
U,β >0 will be considered.
It is convenient to describe the system 
by the number of bosons populating the right well, 0≤ n ≤ N.
If no tunnelling is possible, Ω=0, the eigenstates of Ĥ_


    |n = (c_L^+)^N-n(c_R^+)^n|0/[n!(N-n)!]^1/2

correspond to the energies 

    (n)=U(n-ν)^2+_0(N,β),  n=0,1,...,N,

where

    ν≡ (N-β/U)/2, and_0≡ U[N(N-1)/2-ν^2].


The dependence of (n) on n is quadratic and, since (n)=(2ν-n), 
whenever the minimum of the parabola, at n=ν is an integer, 
or half integer, 

    (N-β/U)=K,  K=1,2,...

there are several pairs of doubly degenerate states. 
Fig. 2 shows the energy levels for an asymmetric potential with N=8 and β/U=1, so the minimum energy corresponds to ν=3.5. We have degeneracy between several pairs of states, for instance, between |1 and |6, that share the same energy, (1)=(6).
 In general, for ν≤ N/2 the number of degenerate pairs
is given by  that of the integers in the interval [0,ν), and 
for ν > N/2, by the number of integers inside (ν/2,N].
 We note that if ν happens to be an integer, there is  an unpaired  non-degenerate ground state
|n=ν, and ( ν)= _0. 
Next we switch the tunnelling on, in such a manner that the tunnelling matrix element in Eq.(<ref>) will remain small, 
compared to the interatomic interaction
Ω /U <<1.
The perturbation will lift the degeneracy
between the levels 
|n and |2ν-n, and introduce pairs of new eigenstates 

    |n_±=[|n± |2ν-n]/√(2)

with the energies

    (n_±)≡(n)±ω_n = (2ν -n)±ω_n,

with the  splitting 2ω_n small, compared to other energy
differences shown in Fig. 2.



 If as a result, the atoms, initially prepared in the
state |Ψ(0)=|n with n atoms in the right well, will perform collective Rabi oscillations
between the states |n and |2ν-n, we expect their state at a time t to be given by

    |Ψ(t)≈ [cos(ω_nt)|n+αsin(ω_nt)|2ν-n]exp(-i(n) t),     |α|=1.

For a symmetric setup, the frequencies ω_n were obtained in <cit.> within time-independent 
perturbation theory, using a rather complicated procedure <cit.> for diagonalising a 
tridiagonal matrix with degenerate eigenvalues. 
The authors of <cit.> correctly note that the shortest way to reach the state |N-n from |n
is by "moving one boson per step". There is, however, room for a further clarification.
In the approximation (<ref>), the atoms will never be observed
 in any other state
|m, m n,2ν-n  (or rather the probability of such an observation
will be negligibly small).
For example, in the case shown in Fig. 2, for |Ψ(0)=|1, counting atoms
in the right well will  almost always give a result either 1, or 6. This suggests,
and we will show below that this suggestion is correct, that the 5 atoms 
making the difference, will need to tunnel, in some sense,  all at the same time.
Next we  demonstrate this by using time-dependent perturbation theory,
and obtain, while we are at it, collective tunnelling frequencies for an asymmetric trap (β 0).  







§ PERTURBATION THEORY AND THE COLLECTIVE FREQUENCIES

Writing the wave function of the system as


    |Ψ(t) = ∑_n=0^N b_n(t) |n,



and introducing dimensionless time and energy,

    τ≡ Ut,  E_n≡(n) /U ,


we obtain the equations for the coefficients  b_n(τ) 

    i∂_τb_n=∑_m=0^N(E_nδ_nm+γ_n-1δ_n,m+1+γ_nδ_n,m-1)b_m


where


    γ_n≡ U^-1Ω [(n+1)(N-n)]^1/2.

A probability to find n atoms in the right well is, therefore, given by

    p(n)=|b_n(t)|^2.

Equations (<ref>) can be solved analytically only for N≤ 3 in general, and for N≤ 7 if the potential is symmetric, as discussed in Appendix B. In order to evaluate the Rabi frequencies for an arbitrary N,
we note that since no matrix element connects the states |n and |2ν-n directly,
and the latter can only be reached from the former via 2(ν-n) intermediate 
steps, shown in Fig. 2. 



For simplicity, we first consider the symmetric potential, β=0, 
assume, for the moment,


that the γs depend on time, let the system start with n bosons in the right well,
and employ the time-dependent perturbation theory.
To the leading approximation, the transition amplitude between the states |n and |N-n

 is given by (M≡ N-2n, n≤ N/2)

    A_n,N-n(τ)≡ N-n|exp[-i∫_0^τ H(τ')dτ']|n =
    (-i)^M∫_0^τ d_M...∫_0^_3d_2 ∫_0^_2d_1×exp[-iE_N-n(-_M)]γ_N-n-1(_M)×exp[-iE_N-n-1(_M-_M-1)]...
    ×γ_n+1(_2)exp[-iE_n+1(_2-_1)]×γ_n(_1)exp(-iE_n_1)+O(Ω^M+1/U^M+1).


This expression has the standard interpretation <cit.>:
the first of the N-2n atoms jumps into the right well at t_1,
the second at t_2, and so on. The transition amplitude is then found
by summing over all t_j, leaving the precise moments of jumps
indeterminate.  Perturbative treatment will limit us to 
times much shorter that the Rabi period, yet as we will see below, 
it captures the essential features of the collective transfer mechanism.
Next we demonstrate that for Ω/U << 1,  and Ω(t) slowly varying 
compared to the inverse of the separations between the system's levels in Fig.2, the sum is dominated 
by the process in which all N-2n atoms jump roughly at the same time.
Returning to the original unscaled time variable, t, 
we have

    A_n,N-n(T) =
    (-i)^Mexp[-i(N-n)T]×∫_0^Tdt_M W_N-n-1(t_M)exp[iΔ_Mt_M]...
    ×∫_0^t_3dt_2W_n+1(t_2)exp[iΔ_2 t_2]×∫_0^t_2dt_1W_n(t_1)exp(iΔ_1 t_1)+O(Ω^M+1/U^M+1)

 where Δ_m is the energy, separating two adjacent states, 

    Δ_m≡(n+m)-(n+m-1)
    =U[1-2(n+m)-N],

and W_n(t)≡γ_nU= Ω(t) [(n+1)(N-n)]^1/2 is the scaled 
tunnelling amplitude.



Thus, we have to evaluate M oscillatory integrals. 
Oscillations of exp(iΔ_n t) become more rapid as U increases,
so that, for a large U/Ω,  the main contributions to the integrals will come from the
endpoints (for details see Appendix A). For example, we may write

    ∫_0^t_2dt_1W_n(t_1)exp(iΔ_1 t_1)= 
    [W_n(t_2)exp(iΔ_1 t_2)-W_n(0)]/iΔ_1+O(Ω/U).

Here the first term corresponds to the first particle jumping at the 
same time as the second, t_1=t_2. The second term clearly corresponds
to the first jump occurring at t_1=0, i.e., immediately after the tunnelling
is switched on.
Jumps at 0 < t_1<t_2  are suppressed,  for large U/Ω, due to destructive 
interference.
Continuing  in the same vein, we obtain a total of 2^M terms ranging 
from all particles jumping at t=0 to all particles jumping at the same time.
In the case of an exact resonance,
∑_m=1^M Δ_m=0,
 the contribution from all fully coordinated jumps between t=0 and t=T is readily seen to be

    i(-1)^M∫_0^Tdt
    ∏_j=0^M-1W_n+j(t)×
    [Δ_1(Δ_1+Δ_2)....
    (Δ_1+Δ_2+...+ Δ_M-1)]^-1.

Note that since the amplitude results from the interference between all t's, 
the exact moment in which the collective transfer takes place remains 
indeterminate, much like the number of the slit chosen by an electron 
in Young's double-slit experiment.
Returning to the case of constant Ω, we, therefore, obtain

    A_n,N-n(T) =i(-1)^MT
    Ω^M/U^M-1exp[-i(N-n)T]×(n+1)...(N-n)/(E_n+1-E_n)
    
    ....(E_N-n-1-E_n) +R,

where the remainder R
 contains the terms corresponding to, at least, one atom jumping 
 immediately after tunnelling is turned on.
 In the next Section we will demonstrate that  R will vanish if 
 tunnelling is switched on sufficiently slowly. 
 
For a sufficiently small T, Eq.(<ref>) predicts

    A_n,N-n(T) =αω_n T.

Comparing Eq.(<ref>) with Eq.(<ref>), and evaluating the products, yields







    ω_n=Ω^N-2n/U^N-2n-1(N-n)!/n![(N-2n-1)!]^2,
      α=(-1)^ni,

which agrees with the result obtained by a different method in <cit.>.
A calculation  for an asymmetric well can be done in exactly the same way, and here we will 
only quote the final result. As discussed in Sec. II, the resonances between
many-body states occur provided ν=1/2,1, ...,N/2. The Rabi 
frequency for a process in which the atoms start in a state |n, and 2(ν-n),
[0≤ n≤ν-1 if ν is an integer, and 0≤ n≤ν-1/2, if ν is odd]
where the atoms are transferred simultaneously, is given by [ν≡ (N-β/U)/2]

    ω_n=Ω^2(ν-n)/U^2(ν-n-1/2)[2(ν-n)-1]!^2
    √((N-n)!(2ν-n)!/n!(N+n-2ν)!).

For an asymmetric well, the leading probabilities p(n), to find n atoms in the right well, 
after switching the tunnelling on at t=0,
are shown 
in Fig.3, together with the Rabi oscillations at the frequency (<ref>).

Superimposed upon the Rabi oscillations, there are much faster oscillations (see inset), which will be discussed in the next Section.








§ SUDDEN VS. ADIABATIC SWITCHING ON OF THE TUNNELLING

These small rapid  oscillations are much more pronounced in a system containing just N=2 atoms, as shown in Fig. 4a.

We have already used perturbation theory to obtain slow Rabi frequencies ω_n, and next we will
try to use it in order to explore the origin of these additional oscillatory patterns. 
To the leading order in Ω, the amplitude to find one atom in each well is given 
by Eq.(<ref>),  and we have (from (<ref>), (0)=U, and (1)=0)
 
    A_0,1(T)  ≈
    -i √(2)Ω∫_0^T dt 
    exp[-i(1)(T-t)]exp[-i(0)t] 
    = 2√(2)Ω/Uexp(-iUT/2)sin(iUT/2),

so that the corresponding probability rapidly oscillates,
 
    p(1)  ≈8 Ω^2/U^2sin^2(iUT/2).

As discussed in Appendix A, by turning the tunnelling on gradually, we can avoid splitting the 2-atom cluster at 
t=0, and quench the oscillations. This can be checked directly, e.g., by choosing 

    Ω(t)=Ω[1-exp(-Γ t)],

where 1/Γ is the characteristic switching time, over which the tunnelling matrix
element is brought to its stationary value. Evaluation of the integral in Eq.(<ref>) gives

    A_0,1(T) ≈
    
    √(2)Ω{exp(-iU T)/U-exp[-i(U-iΓ) T]/(U-iΓ)
    -1/U
    +1/(U-iΓ)}.

For Γ /U<<1, and Γ T>>1, this reduces to

    A_0,1(T) ≈√(2)Ω/U  {exp(-iU T), and p(1) ≈ 2Ω^2/U^2.

Oscillations of p(1) disappear, and its largest value is reduced by a factor of four.
The rapid oscillations also disappear from  the amplitude to find 2 atoms in the right well.
This, as described in Appendix A, can be written as a sum of three terms,

    A_0,2(T)  
    ≈ A^(I)_0,2(T) +A^(II)_0,2(T)+A^(III)_0,2(T)=
    i2Ω^2/Uexp(-iUT)T
    -2Ω^2/U^2
    +2Ω^2/U^2exp(-iUT),
 
related to the three scenarios shown in Fig.5. The first term corresponds to both atoms being transferred
together at some unspecified time between t=0 and t=T. The second term comes from the process 
in which one of the atoms jumps immediately, at t=0, while the second one is transferred later.
The third term is the amplitude for both atoms jumping together at t=0
(see Appendix A). With tunnelling switched on gradually from zero we, therefore, have 

    A_0,1(T) ≈ A^(I)_0,2(T) ,and p(1) ≈ 4Ω^4T^2/U^2.

What has been said so far should apply to times much shorter than the large Rabi period,
T_Rabi=2π/ω_0, 1/U<< T<< T_Rabi, and we still need to check that
this analysis remains correct for much longer times, shown in Fig. 4. We note first that the exact 
amplitude [cf. Eq.(<ref>)] has a form, similar to A_0,2 in Eq.(<ref>)
 
    A_0,1(T)  = 
     -2√(2)Ω/κ Uexp(-iUT/2)sin(iκ UT/2),

where  κ =√(1+16Ω^2/U^2). The amplitude remains small, ∼Ω, 
at all times, since the system never stays in the one-atom-in-each-well configuration for long.
We can, therefore, expect  that, in the slow switching-on regime, the probability p(1) would reduce to 
 
    p(1) = 2Ω^2/√(U^2+16Ω^2),

and the oscillations will be eliminated also from p(0) and p(2), 
which depend on p(1), through the condition p(0)+p(1)+p(2)=1.
The ultimate proof consists in solving numerically Eqs.(<ref>) with an Ω(t) given in Eq.(<ref>).
The results are shown in Fig.4b. The Rabi oscillations, delayed until  Ω reaches its final magnitude, 
lack the rapid oscillations seen in Fig. 4a, while p(1) does indeed have a constant value given by (<ref>).
This approach can be extended, with similar results, to more than two atoms.
Before giving our conclusions, we briefly discuss detuned Rabi oscillations.



§ DETUNED RABI OSCILLATIONS
 
At the exact resonance, (n)=(2ν -n)
Rabi oscillations between the corresponding states 
|n and |2ν-n, ν=(N-β/U)/2,
can be described by an effective Hamiltonian,   acting in a reduced two-dimensional Hilbert 
space, 
^eff=ω_nσ_x, with σ_x denoting the Pauli matrix.  
In a slightly asymmetric well, with β replaced by β +δβ,

the two levels are detuned by Δ, |Δ|=|2δβ (n-ν)|

and the effective Hamiltonian acquires an additional term, proportional to σ_z, 
^eff=ω_nσ_x+Δσ_z/2. Here, as in the general case, the 
probabilities are given by the formulae, known for detuned Rabi oscillations. 
In particular, starting with n atoms in the right well we should have

    P(2ν-n)=1-P(n)≈ω_n^2/ω_n^2+δβ^2 (ν-n)^2sin^2[√(ω_n^2+δβ^2 (ν-n)^2t)],
    P(m)<< 1, for m n, 2ν-n

where ν and ω_n are given by Eqs.(<ref>) and (<ref>), respectively.
Detuned oscillations for N=5, n=0, in a slightly perturbed symmetric potential




are shown in Fig. 6.


§ CONCLUSIONS AND DISCUSSION

In summary, we show that a system of strongly interacting bosonic atoms in a not necessarily symmetric double-well potential is capable 
of performing collective Rabi oscillations of a frequency ω_n given by Eq.(<ref>).
This will happen provided the resonance condition (<ref>) is satisfied for the many body states containing n and n+M atoms in one of the wells. 
The frequencies of the Rabi oscillations are easily predicted with the help of time-dependent perturbation theory, 
which recovers the result obtained in <cit.> for the symmetric case.
The physical picture is the one in which 
all M atoms are transferred together, as a single cluster, from one well to another, during a short period of time δ t ∼ U^-1,
while the precise moment of the transfer remains indeterminate, in accordance with the uncertainty principle.
This is achieved through destructive interference between other scenarios.
Compared to  the observation time, typically of the order of the large Rabi period,
the rapid transfer appears to be almost instantaneous. 
Accordingly, an observation  made in either well will almost certainly find there either n, or n+M atoms, with all other counts occurring only rarely. While it is difficult to probe directly the duration of the collective transfer, a less direct proof
of this transfer mechanism is available.
 While the behaviour of the group of M atoms is broadly similar to that of a 
single particle,  the composite nature of the group can be probed  if  tunnelling is switched on suddenly.
If so, the immediate transfer of individual atoms, or smaller groups of atoms, will result in a 
rapid oscillatory pattern, superimposed on the much slower Rabi cycle.
An analysis of the transition amplitude in Appendix A, as well as numerical results shown in Fig.4, 
show that these oscillations disappear, if tunnelling is turned on slowly, compared to δ t,
and no splitting of the transferred group of atoms occurs at the start of the evolution.
For example, the one-atom-in-each well probability p(1), shown in Fig. 4, ceases to oscillate, and assumes a
constant value in the adiabatic limit. This justifies the approximation of the integral in Eq.(<ref>), and of all other 
oscillatory integrals in Appendix A, by a sum of end point contributions, and leaves the transition amplitude dominated
by the apparently simultaneous transfer of all atoms.
We note that the spectrum of these rapid oscillations typically contains combinations of all internal frequencies of the 
system, and does not lend itself to a simple analysis, except in the case M=N=2, discussed in Sec. IV. 
 Finally, we expect the present theory, which treats collective tunnelling beyond the single-cluster picture, 
to find applications in 
experimental studies of such subtle aspects of the phenomenon, as quantum fluctuations in the observed averages, or the system's response to 
 temporal variations of its parameters.


Acknowledgements



Financial support  of
MINECO and the European Regional Development Fund FEDER, through the grant
FIS2015-67161-P (MINECO/FEDER,UE), the Basque Government through Grant No IT986-16.
is also acknowledged  (DS and M.P), and of FIS2016-80681-P (MP).



§ APPENDIX A. OSCILLATORY INTEGRALS AND THE TRANSITION AMPLITUDE

To illustrate the use of oscillatory integrals, discussed in Sec. III, we first consider a
 transition between states |1 and |2, with the energies E_1 and E_2, 
caused by a time-dependent perturbation Ŵ(t)=W(t)[|2 1|+h.c.]. 
If the system is in |1 at t=0, the amplitude for finding it in |2 at t=T, 
to the leading order in Ŵ(t), is given by

    A=-i∫_0^T exp[-iE_2(T-t)]W(t)exp(-iE_1t)dt≡ 
     -iTexp(-iE_2T)∫_0^1 W(τ) exp[iT Δτ]W(t)]dτ,

where τ≡ t Δ. 


 If the observation time 
is large, TΔ >> 1, and W(t) varies slowly, compared to δ t∼ 1/Δ =1/|E_2-E_1|, 

the main contributions to the integral
comes from the incomplete oscillations at the endpoints, 
i.e., from the regions  around 0 and T, of a width ∼δ t.

    A≈ -Δ^-1[W(T) exp(-E_1T)W(0)exp(-iE_2T)].


Now the amplitude A
has lost all the information about the overall 
behaviour of W(t), and depends only on its initial value at t=0, 
and its current value at t=T. 
The short time it takes for the system to "jump" is, clearly ∼δ t.
Moreover, if W(t) is gradually switched
on from zero, the likelihood of finding the system in |2, will 
depend only on the perturbation's current value W(T). 


 §.§ Even number of atoms

Next in complexity is the case in  which two atoms are transferred via a third state.
For simplicity, we consider a symmetric well, β=0, and start transferring an
even number of atoms M=N=2K, so that the minimum of the parabola in Fig. 2 is at an integer value of n.
 The transfer takes two steps, and the corresponding second-order amplitude
is given by 

    A≈ (-i)^2exp(-iE_3T) ∫_0^T dt_2 ∫_0^t_2dt_1×exp(iΔ_2 t_2)W_2(t_2)
    exp(iΔ_1 t_1)W_1(t_1),

where, as before, Δ_i=E_i+1-E_i, the resonance condition reads E_1=E_3, 
and the "jump time" is δ t ∼ 1/(E_2-E_1). The diagram in Fig.7 helps to visualise
the structure of the double integral in Eq.(<ref>). Every integral, containing an exponential exp(iΦ t), 
splits into contributions from its endpoints. An exception is the second integral in the left branch, 
where the phase Φ=Δ_1+Δ_2=0 vanishes, and the integration has to be carried out 
over the whole interval 0≤ t_2≤ T. This process, in which both atoms are transferred 
together at an unspecified moment t_2, contributes to the amplitude in Eq.(<ref>) an amount

    A^(I)≈ iΔ_1^-1exp(-iE_3T) ∫_0^T dt_2 W_2(t_2)
    W_1(t_1).

There is also a process, in which the first atom jumps at t=0, and is followed by the second one at t=T
Its  contribution is

    A^(II)≈ -Δ_1^-2exp(-iE_2T) W_2(T)
    W_1(0).

(Note that the phase is E_2t, since the system continues in the state E_2, 
until the second jump at t=T).
Finally, it is possible for both atoms to jump at t=0, and we have

    A^(III)≈Δ_1^-2exp(-iE_3T) W_2(0)
    W_1(0),

so that A=A^(I)+A^(II)+A^(III). We note that if W_1(t) is switched on 
from zero, slowly compared to the jump time  δ t, the last two terms vanish, 
and the amplitude is fully determined by collective transfer of two atoms, 
A=A^(I) at an unspecified time t_2.
The scheme in Fig. 7 is easily extended to collective transfer of more than two atoms, M>2
provided N remains even. In general, such an amplitude is seen to be dominated by nearly 
instantaneous transfer of  M-atom cluster,


if  W_i(t) varies slowly, compared to δ t=1/min{Δ_m}.
 With the observation times typically of order of the large Rabi period  T_Rabi∼ W^-M,
 the pictures should hold in most practical cases.








 §.§ Odd number of atoms


Although Eq.(<ref>) is valid for any N, 
the situation for M=2K+1, K=1,2,..., is slightly different, since now the 
two lowest states in Fig. 2, through which the system has to pass, are also in resonance. Thus, we have an 
integral

    A= (-i)^3exp(-iE_4T) ∫_0^Tdt_3 ∫_0^t_3dt_2 ∫_0^t_2dt_1×exp(iΔ_3 t_3)W_3(t_3)exp(iΔ_2 t_2)W_2(t_2)exp(iΔ_1 t_1)W_1(t_1),

where Δ_2=0, and Δ_1+Δ_3=0.
Four possible processes, which contribute to the amplitude (<ref>) are shown in Fig.9.
The left branch on the diagram in Fig.7 corresponds, as before, to all atoms jumping at the 
same, yet unspecified, time, for which we have, 

    A^(I)=-iexp(-iE_4T)Δ_1^-2×∫_0^T W_1(t_3)W_2(t_3)W_3(t_3)dt_3.

The process in which the first two atoms jump at t=0, and the third at t=T, 
contributes 

    A^(II)=iexp(-iE_3T)Δ_1^-2Δ_3^-1×
     W_1(0)W_2(0)W_3(T),

while all three atoms jumping at t=0 add

    A^(III)=-iexp(-iE_4T)Δ_1^-2Δ_3^-1×
     W_1(0)W_2(0)W_3(0).

Finally, it is possible for the first atom to jump at t=0, 
for the second one to jump at some 0≤ t_2≤ T, 
and for the third one to complete the transition at t=T.
The amplitude for this is

    A^(IV)=-iexp(-iE_4T)Δ_1^-1Δ_3^-1×
     W_1(0)W_3(T)exp(iΔ_3T)∫_0^TW_2(t_2)dt_2.

The amplitude for the fifth process in Fig.9 vanishes, since it involves
∫_0^0W(t_2)dt_2=0, and the full amplitude is the sum of four terms
A=A^(I)+A^(II)+A^(III)+A^(IV).
Typically, all W_i are of the same order of magnitude, W_1∼ W_2∼ W_3∼ W, 
and in the strong interaction limit, Wδ t << 1, the second and the third terms, 
which contain an extra factor of Δ_3^-1, can be omitted. 
For all W_i chosen to be constant, both A^(I) and A^(IV) grow linearly with time, 

    A^(I)∼ T,  A^(IV)∼ T exp(iΔ_3T),

as shown in the inset of Fig.10.
However, A^(IV) rapidly oscillates, compared to the large Rabi period, and can be omitted, 
when the Rabi frequency is calculated from the rate of change of A. 
Although not immediately seen from Eq.(<ref>), this becomes clear from the graphs in
the main panel of Fig. 10. The oscillations do not increase proportional to time, and the linear
growth of the amplitude is determined, as one would expect, only by collective transfer of all 
three atoms. As before, switching tunnelling on slowly from zero, helps avoiding the splitting of the 
tunnelling cluster, and eliminates A^(II), A^(III), and A^(IV), all proportional to W_1(0).
The diagram in Fig.9 can be extended to the transfer of more than three atoms, by adding integrations and branches, as appropriate. In general, the amplitude would contain the leading term, describing nearly simultaneous 
transfer of all M atoms,  as well as the additional sub-amplitudes, if the cluster is split at t=0  by sudden onset 
of tunnelling.







§ APPENDIX B. EXACTLY SOLVABLE CASES
 
To write down a formal analytic solution to Eqs.(<ref>),
 we need to diagonalise the Hamiltonian matrix
 H_mn≡ m|Ĥ|n, m,n=0,1,...,N. 
 Since algebraic equations have analytical solutions up to the 
 fourth order, this can be done for no more than N=3 atoms,
 in the case of a general asymmetric well.
 For a symmetric well, β =0, a solution can be found for up to N=6 
 atoms. Additional symmetry of H_mn,
H_mn=H_N-m,N-n
requires that its eigenvectors a^n≡
(a_0^n,a_1^n,...,a_N^n)
have a definite parity,

    a^n_m=α a^n_N-m,  α = ± 1.

This, in turn, allows us to reduce the size of the matrix 
we need to diagonalise in order to obtain the eigenvalues of H_mn,
and an analytical solution can, in principle, be obtained for N≤ 7.
Below we analyse the N=2 and N=3 cases, assuming β = 0.



 §.§ The two-particle case (N=2)

For N=2, we write A_n,N-n(t) =exp(-iUt)b_n(t), 

where b_n(t) satisfy

    i∂_t b_0=√(2)Ω b_1
    
    i∂_tb_1=√(2)Ω b_0-Ub_1+√(2)Ω b_2
    
    i∂_t b_2=√(2)Ω b_1

with the initial condition
b_i(t=0)=δ_i0,    i=0,1,2.
This is equivalent to a second order equation for b_1

    ∂_ttb_1-iU∂_t b_1+4Ω^2b_1=0

with the initial condition
b_1(0)=0,   ∂_t b_1(0)=-√(2)iΩ.
The corresponding solution is



    b_1(t)=-√(2)Ω/ϵ_1-ϵ_2
    [exp(iϵ_1t)-exp(iϵ_2t)]

where.
ϵ_1,2=U/2±√(U^2/4+4Ω^2).
For two other coefficients we then have

    b_2(t)=-√(2)iΩ∫_0^t b_1(s)ds = 
    2Ω^2/(ϵ_1-ϵ_2)×
    [exp(iϵ_1t)/ϵ_1-exp(iϵ_2t)/ϵ_2]
    +2Ω^2/ϵ_1ϵ_2

and

    b_0(t)=1+b_2(t).

Returning to the strong interaction case, Ω/U <<1, we note then that

    ϵ_1-ϵ_2≈ϵ_1≈ U,
    ϵ_2 ≈- 4Ω^2/U,ϵ_1ϵ_2≈ -4Ω^2.

With this we have

    b_0=exp(-iω_0t)cos(ω_0 t), 
     b_2=-iexp(-iω_0t)sin(ω_0 t),ω_0=2Ω^2/U,

which agrees with Eq.(<ref>).
For the condensate prepared initially with one boson in each well
we must change the initial condition for b_1(0) to
b_1(0)=1,   ∂_τ b_1(0)=iU.
Hence we have

    b_1(t)=U-ϵ_2/ϵ_1-ϵ_2exp(iϵ_1t)
    +ϵ_1-U/ϵ_1-ϵ_2exp(iϵ_2t)

and
b_0(t)=b_2(t)=-i√(2)Ω∫_0^tb_1(s)ds.
It is readily seen that in the strong interaction limit
Ω/U<<1, when the energy separation from the two states 
corresponding to two bosons in the same  well grows as U, the 
system stays in its ground state |1. Indeed, rapid oscillations
of  b_1(t)≈exp(-iUt), guarantee that both 
b_0(t) and b_2(t) remain of order of  Ω/U).



 §.§ The three-particle case (N=3)

For three bosons, N=3, in a symmetric well, β=0, 
to evaluate  b_n=exp(2iUt) A_n,N-n, we need to solve 

    i∂_t b_0=√(3)Ω b_1
    
    i∂_tb_1=√(3)Ω b_0-2Ub_1+2Ω b_2
    
    i∂_tb_2=2Ω b_1-2Ub_2+√(3)Ω b_3
    
    i∂_t b_3=√(3)Ω b_2.

This is equivalent to 

    i∂_t b_0=√(3)Ω b_1
    ∂_ttb^±+2iU^±∂_tb^±+3Ω^2b^±=0
    
    
    
    
    i∂_t b_3=√(3)Ω b_2.
 
where

    b^±≡ b_1± b_2   and   U^±≡±Ω -U
 
which are solved with appropriate initial conditions discussed below.
Initially, no particles in the right well (n=0): b_n(0)=δ_n0,   n=0,1,..3.
Introducing (subscripts 1 and 2 correspond to the 
+ and - signs of the square root) 

    ^±_1,2=[U^± ]±√((U^±)^2+3Ω^2)

we have

    b_0=
    
    -3/2Ω^2{exp(i_1^+t)/_1^+(_1^+-_2^+)
    -exp(i_2^+t)/_2^+(_1^+-_2^+)
    +exp(i_1^-t)/_1^-(_1^–_2^-)
    -exp(i_2^-t)/_2^-(_1^–_2^-)}
    b_1=-√(3)/2Ω{exp(i_1^+t)-exp(i_2^+t)/_1^+-
    _2^++
    exp(i_1^-t)-exp(i_2^-t)/_1^–_2^-}
    b_2=-√(3)/2Ω{exp(i_1^+t)-exp(i_2^+t)/_1^+-
    _2^+-
    exp(i_1^-t)-exp(i_2^-t)/_1^–_2^-}
    b_3=-3/2Ω^2{exp(i_1^+t)/_1^+(_1^+-_2^+)
    -exp(i_2^+t)/_2^+(_1^+-_2^+)
    -exp(i_1^-t)/_1^-(_1^–_2^-)
    +exp(i_2^-t)/_2^-(_1^–_2^-)}
 
In the strong interaction limit
Ω/U<<1 we have

    _1^±-_2^±≈^±_1≈ 2U, 
    e^±_2≈ -3Ω^2/2U∓ 3Ω^3/2U^2

so that

    b_0(t) ≈ -exp(-3iΩ^2t/2U) cos(3Ω^3t/2U^2),   
    b_3(t) ≈ iexp(-3iΩ^2t/2U) sin(3Ω^3t/2U^2),   
    b_1(t)≈ b_2(t) =
    o(Ω/U)

in accordance with ω_0 ≈ 3Ω^3/2U^2,  predicted by Eq.(<ref>).

Initially, one particle in the right well (n=1):

    b_n(0)=δ_n1,   n=0,1,..3.




Solving Eqs.(<ref>) with this initial condition we have

    b_0=-√(3)Ω/2×
    {-_2^++2U^+/_1^+(_1^+-
    _2^+)exp(i_1^+t)+ 
    _1^++2U^+/_2^+(_1^+-
    _2^+)exp(i_2^+t)
    -_2^-+2U^-/_1^-(_1^–_2^-)exp(i_1^-t)+
    _1^-+2U^-/_2^-(_1^–_2^-)exp(i_2^-t)},
    b_1=-1/2
    {-_2^++2U^+/_1^+-
    _2^+exp(i_1^+t)+
    _1^++2U^+/_1^+-
    _2^+exp(i_2^+t)
    -_2^-+2U^-/_1^–_2^-exp(i_1^-t)+
    _1^-+2U^-/_1^–_2^-exp(i_2^-t)}


    b_2=-1/2
    {-_2^++2U^+/_1^+-
    _2^+exp(i_1^+t)+
    _1^++2U^+/_1^+-
    _2^+exp(i_2^+t)
    +_2^-+2U^-/_1^–_2^-exp(i_1^-t)-
    _1^-+2U^-/_1^–_2^-exp(i_2^-t)}
    b_3=-√(3)Ω/2×
    {-_2^++2U^+/_1^+(_1^+-
    _2^+)exp(i_1^+t)+
    _1^++2U^+/_2^+(_1^+-
    _2^+)exp(i_2^+t)
    +_2^-+2U^-/_1^-(_1^–_2^-)exp(i_1^-t)-
    _1^-+2U^-/_2^-(_1^–_2^-)exp(i_2^-t)}.


(Note that the expressions for b_0 and b_3 do not contain the 
time-independent terms obtained in integrating b_1 and b_2.
Since the hamiltonian matrix has four non-zero
eigenvalues, b_0 and b_3 cannot contain an additional
zero frequency.)
In the strong interaction limit
Ω/U<<1 we have



    b_1(t) ≈ -exp(2iUt)cos(2Ω t), 
    b_2(t) ≈ -iexp(2iUt)sin(2Ω t), 
     b_0(t)≈ b_3(t) =O(Ω/U),

in accordance with ω_1 ≈ 2Ω, obtained from Eq.(<ref>).



 10

 Rabi  I. I. Rabi, S. Millman, P. Kusch, J. R. Zacharias, Phys. Rev. 1939, 55, 526.
   
 Dud  Y. O. Dudin, L. Li, F. Bariani, A. Kuzmich, Nature Physics 2012, 8, 790.
  
Bloch S. Foelling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, A. Widera, T. Mueller, I. Bloch,
Nature 2007, 448, 1029.
 

 
 Raiz M. G. Raizen, S. Wan, C. Zhang, Q. Niu. Phys. Rev. A 2009, 80, 030302(R).
 
 MB1 
K. Sakmann, A. I. Streltsov, O. E. Alon, L. S. Cederbaum, Phys. Rev. Lett. 2009, 103, 220601.
 
 MB2 
R. l. Beinke, S. Klaiman, L. S. Cederbaum, A. I. Streltsov,1,  O. E. Alon,  Phys. Rev. A  2015, 92, 043627.
 
 MB3 
J. Erdmann, S.I. Mistadikis, P. Schmelcher,  Phys. Rev. A 2018,  98, 053614.
 
 BH1 R. Ma, M. E. Tai, P. M. Preiss, W. S. Bakr, J. Simon,  M. Greiner, Phys. Rev. Lett. 2011, 107, 095301.
  
 BH2 F. Meinert,  M. J. Mark, E. Kirilov, K. Lauber, P. Weinmann, M. Groebner, A. J. Daley, H. C. Naegerl, Science 2014, 344, 1259
 Ref2_5 S. I. Mistakidis, L. Cao, P. Schmelcher, Mol. and Opt. Phys. 2014, 47, 225303. 
 
 Ref2_6 S. I. Mistakidis, L. Cao, P. Schmelcher, Phys. Rev, A 2015 91, 033611. 
 
Ref2_7 J. Neuhaus-Steinmetz, S. I. Mistakidis, P. Schmelcher, Phys. Rev. A 2017, 95, 053610. 

 Ref1_1 L. Cao, I. Brouzos, B. Chatterjee, P. Schmelcher, New. J. Phys. 2012, 14, 093011.
  
 Ref1_2 R. E. Barfknecht, A. Foerster, N. T. Zinner, New. J. Phys. 2018, 20, 063014.
  
 Ref2_2 G.J. Milburn, J. Corney, E.M. Wright, D. F. Walls, Phys. Rev. A 1997, 55, 4318.
   
Ref2_3 A. Smerzi, S. Fantoni, S. Giovanazzi, S. R. Shenoy, Phys. Rev. Lett., 1997,79, 4950. 
   
Ref2_1 M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, M. K. Oberthaler, Phys. Rev. Lett.  2005, 95, 010402.
  
Ref2_4 Y. Chen, S. Nascimbène, M. Aidelsburger, M. Atala, S. Trotzky, I. Bloch, Phys. Rev. Lett., 2011, 107, 210405.
   




BH0 G. Kalosaka, A. R. Bishop,  V. M. Kenkre,  J. Phys. B., 2003, 36, 3233.

Bern L. Bernstein, J. C. Eilbeck, A. C. Scott, Nonlinearity, 1990, 3, 293.

BH00 G. Kalosaka, A. R. Bishop,  V. M. Kenkre,  Phys. Rev. A, 2003, 68, 023602.

FOOTall 
Imagine there are two rooms, one of which is full of people. The people move from one room to the other, and each time one checks, all people are found in the same room. This  implies that whenever people decide to move, they do it not one by one, but all together. The same must be the case with the bosons.

Feyn R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York 1965. 










In particular, for an odd number of particles, N=2K-1,
we need to diagonalise two K× K matrices obtained by 
taking the K× K left upper corner sub-matrix of H_mn and
replacing H_KK with  H_KK± H_K,K+1 for α equal to 
1 and -1, respectively.

For an even number of particles, N=2K, and α=1 we need to 
diagonalise the K+1× K+1 matrix obtained by taking the K+1× K+1
 left upper corner sub-matrix of H_mn  and replacing 
H_K+1,K with H_K+1,K+H_K+1,K+2.
For the negative parity, α=-1, we simply need to diagonalise
the K× K left upper corner sub-matrix of H_mn.
Bearing in mind that algebraic equations only have analytical solutions
up to the fourth order, we may conclude that an analytical solutions
can, in principle, be obtained for N=2,3,4,5,6.

To conclude the Section we compare predictions of Eq.(<ref>) with some known results.

a). The  2-boson case, the frequency of two-boson transfer between the states
|0 and |2 (N=2,n=0, E_0=0,E_1=1:

    ω_0= 2Ω^2/U


b). The  3-boson case, the frequency of three-boson transfer between the states
|0 and |3 (N=3,n=0, E_0=0,E_1=2, E_2=2:

    ω_0= 3Ω^3/2U^2


c). The  3-boson case, the frequency of one-boson transfer between the states
|1 and |2 (N=3,n=1, E_0=0,E_1=2, E_2=2:

    ω_1= 2Ω


    b_0=1/2i
    {-_2^++2U^+/_1^+(_1^+-
    _2^+)exp(i_1^+t)+ 
    _1^++2U^+/_2^+(_1^+-
    _2^+)exp(i_2^+t)
    -_2^-+2U^-/_1^-(_1^–_2^-)exp(i_1^-t)+
    _1^-+2U^-/_2^-(_1^–_2^-)exp(i_2^-t)},
    b_1=1/2
    {-_2^++2U^+/_1^+-
    _2^+exp(i_1^+t)+
    _1^++2U^+/_1^+-
    _2^+exp(i_2^+t)
    -_2^-+2U^-/_1^–_2^-exp(i_1^-t)+
    _1^-+2U^-/_1^–_2^-exp(i_2^-t)}
    b_2=1/2
    {-_2^++2U^+/_1^+-
    _2^+exp(i_1^+t)+
    _1^++2U^+/_1^+-
    _2^+exp(i_2^+t)
    +_2^-+2U^-/_1^–_2^-exp(i_1^-t)-
    _1^-+2U^-/_1^–_2^-exp(i_2^-t)}
    b_3=1/2i
    {-_2^++2U^+/_1^+(_1^+-
    _2^+)exp(i_1^+t)+
    _1^++2U^+/_2^+(_1^+-
    _2^+)exp(i_2^+t)
    +_2^-+2U^-/_1^-(_1^–_2^-)exp(i_1^-t)-
    _1^-+2U^-/_2^-(_1^–_2^-)exp(i_2^-t)},



we consider a system of N strongly interacting bosonic atoms, trapped in a double-well potential.
Its energy spectrum in the absence of tunnelling is shown in Fig.2.
As usual (see, e.g., <cit.>), the strong interaction regime is realised if the strength of the attractive, or repulsive
interaction, is large compared to the one-particle 
tunnelling matrix element. 
With the resonance condition (<ref>) satisfied for the many body states containing n and n+δ n atoms,
tunnelling the system will perform Rabi oscillations with a frequency ω_n in Eq.(<ref>).
In consequence, an observation made on the right well will almost certainly find there either n, or n+δ n atoms,
with all other counts occurring only rarely. The picture is the one in which all of δ n atoms are transferred together, 
during a short period of time ∼ U^-1, which is also roughly the range of times contributing to the oscillatory integrals (<ref>).
The precise moment of the transfer remains indeterminate, since transfers at all times equally contribute to the amplitude in Eq.(<ref>). 

Next we extend this approach to the M-th order of the perturbation theory
in Eq.(<ref>), where the short "jump" time is now δ t=1/min{Δ_m}, 
m=1,2,..,M, T>>δ t, and W(t) varies slowly, compared to δ t
W(t+δ t)≈ W(t). 
Now there are many endpoint contributions,
and instead of analysing the general case, we consider here transfer of only 
three atoms, M=3 
(The simplest case M=2 is studied in Sec. IV.)
We consider, to the leading approximation, a transition
|4 |3 |2 |1, between the satiates with the energies E_i,
We have an integral 

    A= (-i)^3exp(-iE_4T) ∫_0^Tdt_3 ∫_0^t_3dt_2 ∫_0^t_2dt_1×exp(iΔ_3 t_3)W_3(t_3)exp(iΔ_2 t_2)W_2(t_2)exp(iΔ_1 t_1)W_1(t_1).
,
where W_i(t) are matrix elements of a time dependent perturbation W, 
and Δ_i=E_i+1-E_i.
 Following the outmost left-hand branch of the diagram in Fig.7, 
and imposing the resonance condition  Δ_1+Δ_2+ Δ_3 =E_4-E_1=0, 
we obtain the results in Eq.(<ref>) [cf. Eq.(<ref>)]

    A_1≈ -exp(-iE_3T)∫_0^TW_1(t)W_2(t)W_3(t) /Δ_1(Δ_1+Δ_2)(Δ_1+Δ_2+Δ_3),

corresponding to all three atoms making a transition within a small vicinity ∼δ t of
an unspecified moment t. Following the outmost right-hand branch gives an amplitudes of a process
in which all three atoms jump at t=0. In Fig.7 there are two other branches, corresponding to the same 
process and, adding their contributions, we have
 
    A_2≈ -exp(-iE_3T)W_1(0)W_2(0)W_3(0)[Δ_1^-1Δ_2^-1Δ_3^-1+
    -Δ_1^-1Δ_2^-1(Δ_2+Δ_3)^-1
    +Δ_1^-1(Δ_1+Δ_2)^-1Δ_3^-1].

Then there is  a process, in which two atoms jump at t=0, and the remaining one 
does so at t=T, with an amplitude
 
    A_3≈exp(-iE_2T)W_1(0)W_2(0)W_3(T)×
    [(Δ_1Δ_2Δ_3)^-1-
    Δ_1^-1(Δ_1+Δ_2)^-1Δ_3)^-1]

and also a process in which one atom jumps at t=0, and two at t=T, 
 
    A_4≈ -exp(-iE_1T)W_1(0)W_2(T)W_3(T)×Δ_1^-1Δ_2^-1(Δ_2+Δ_3)^-1.

Note that if W(t) is switched on in such a way that W_i(0)=0, 
the processes in which one of the atoms breaks away from the other to disappear, 
A_3=A_4=0, and only collective behaviour of all three is possible.
For a constant W, W_i W_i(t), only A_1 in Eq.(<ref>) grows ∼ T, 
while the rest of the amplitudes are responsible for a small oscillatory structure, superimposed 
on the linear growth, and visible also in the transition probability P(4 1)=|A|^2.
The scheme shown in Fig. 7 is easily extended to any number of atoms transferred,
split into groups making transitions at t=0 and t=T.



The frequencies of the Rabi oscillations are easily predicted with the help of time dependent perturbation theory. 
In general, there is no direct matrix element connecting the resonance states and, if the atoms repel each other, 
the system has to go down first, and then climb the ladder of intermediate states, as shown in Fig. 2. 
(If the interaction is attractive, parabola in Fig.1 is turned upside down, and the system has to climb first).
In either case, the transitions first occur in the higher order of perturbation theory, and the Rabi periods grow with the number of atoms
transferred, as seen from Eq.(<ref>).
At the exact resonance, for just two atoms in the strong coupling regime, N=2, 

the time evolution of the states, containing initially  0 and 2 
atoms in the right well, can be written as [cf. (<ref>)] 

    |
    |ψ_0(t) = exp(-iω_0t)[cos(ω_0 t)|0 -isin(ω_0 t)|2] 
    |ψ_n(t) = exp(-iω_0t)[cos(ω_0 t)|2 -isin(ω_0 t)|0].

This can be used to demonstrate that for an arbitrary state 
|Ψ = a_0|0 + a_2(t)|2, the coefficients a_0(t) and 
a_2(t) evolve according to a Schroedinger equation with an effective 2× 2 
Hamiltonian ^eff_BEC=ω_0[I+σ_x], where σ_x is the Pauli matrix.