Semiclassical analysis of Bose-Hubbard dynamics

In this work the two site Bose-Hubbard model is studied analytically in the limit of weak coupling u and large number of particles N . The semiclassical approximation where \frac{1}{N} plays the role of Planck's constant was used and perturbation theory to order u^{2} was applied. In particular, the difference in the occupation between the two sites, where initially all particles are at one site was calculated analytically. Excellent agreement with the exact numerical solution was found. This quantity exhibits collapses and revivals that superimpose rapid oscillations. The occupation difference was calculated also for the case where initially both sites are occupied provided that the difference in occupation is sufficiently large. It provides an analytical description of results that were so far found only numerically. Similar behavior and analysis are expected for a large variety of physical situations in optics, atom optics and quantum dynamics of electrons in Rydberg atoms.


I. INTRODUCTION
The physics of Bose-Einstein condensates (BECs) is extensively studied in the recent years [1][2][3]. For weakly coupled atoms in a large variety of systems the Gross-Pitaevskii equation (GPE) [1,3] describes well the static properties. For the dynamics, the situation is more complicated and it is instructive to study simple paradigmatic systems. For the double well potential, the GPE does not reproduce the correct dynamics. In particular, the collapse of the amplitude of Rabi oscillations as a function of time is not reproduced, in contradiction with the numerical solutions for the many-body system [4,5]. The double well is a paradigmatic model system that was extensively studied experimentally [6][7][8][9].
Much interest was in the Bosonic Josephson effect [6,7]. This is a clear manifestation of macroscopic quantum coherence. It encourages theoretical exploration of this and related systems [4,[10][11][12][13]. If the inter-particle interaction is sufficiently weak so that the coupling between the lowest levels of the double well and higher levels can be ignored, the system may be described by the two site Bose Hubbard (BH) model where bosons can occupy only two sites [14]. This model has been studied numerically and analytically [4,10,[15][16][17][18][19][20].
Fascinating phenomena that were explored are collapse and revival of the difference in the occupation of the sites [4]; and the semiclassically-related statistics of the fluctuations that are associated with the occupation [10,21,22]. The latter are important for the fringe visibility in interference experiments [9]. The understanding of the revivals is crucial for understanding of the coherence in the system.
Collapses and revivals were observed in experiments where a BEC was confined to a lattice. The interference pattern of the matter wave field originating in different lattice sites showed collapses and revivals as a function of time [23,24]. These were found also experimentally for other condensates [25,26] and Rydberg atoms [27,28]. Coherence was explored for models of dynamics of atoms on optical lattices in [29,30]. Collapses and revivals were found in numerical calculations following heuristic arguments and in direct numerical studies of the double well problem [5,15,16]. For the two site Bose-Hubbard model, collapses were found in exact numerical calculations [15,16]. Collapses and revivals were found theoretically for interacting bosons in a harmonic well and the relevant times were estimated [31][32][33][34][35]. These were found also for wave packets in harmonic wells with small nonlinearities [36]. In quantum optics these were found in the Jaynes-Cummings model [37] analytically and numerically [38]. It is the closest to the one found in the present paper for interacting bosons. In optics this phenomenon is well understood and is known as the Talbot effect [39,40] (see also [41,42]). A related phenomenon is the "Quantum carpet" [42,43].
Collapses and revivals can be found in many situations. A generic picture is outlined in [2,44].
A semiclassical picture for the two site Bose-Hubbard model was developed and studied in some detail [4,10,17]. The dimensionless parameter that controls the corresponding classical behavior is where U is the inter-particle interaction, J is the strength of the hopping between the sites, while N is the number of particles. In this picture, 1 N plays the role of Planck's constant and the thermodynamic limit N → ∞ plays the role of the classical limit. The various regimes are [10]: 1. Rabi regime u < 1 2. Josephson regime 1 < u < N 2 3. Fock regime N 2 < u.
The Josephson regime is the most extensively studied one [10,17,18]. It exhibits an interesting phase space, with dynamics related to the experimental observations [6,7]. We confine ourselves to the Rabi regime where the classical behavior is very simple. It enables us to study analytically quantum collapses and revivals that are crucial in the understanding of coherence.
The present work will follow the formulation presented in detail in the work of D. Cohen and coworkers [10]. In the limit N ≫ 1 we find an analytical formula for the difference between the occupation of the two sites. It is rare to find such results for interacting systems and to the best of our knowledge it is the first time such an analytic expression is found in the present context, namely for interacting bosons.
The outline of the paper is as follows: In section II the model is defined and the semiclassical picture is presented. In section III a transformation to angle action variables is performed and used to find the energies within the WKB approximation to the order u 2 , and in section IV it is calculated in standard quantum perturbation theory. In section V the difference in occupation between the two sites is calculated for an initial condition where all atoms are on one site, while in section VI the initial condition where both sites are occupied is used. The results are discussed in Sec. VII.

II. THE TWO STATES HUBBARD MODEL
The two states Hubbard model we study is defined by the Hamiltonian The sites are denoted by L (Left) and R (Right). The creation and annihilation operators on the sites are a † L , a † R and a L , a R . The number operators for the two sites are n L = a † L a L and n R = a † R a R . The commutation relations are a L , a † L = 1, a R , a † R = 1, and the units are such that = 1 N . It is assumed that the on site energies on the two sites are identical. The total number of particles n L + n R = N is conserved. The first term in (2) represents the hopping between the two sites while the second one is the energy of the interparticle interaction that in the present work is assumed to be small. The Hamiltonian (2) can be written in the form By using the angular momentum operators (see for example [10]) the Hamiltonian (3) can be written up to a constant as (shown in App. A) Namely, where C N = 1 2 N 2 U − NU. The operators (4) satisfy the commutation relations of angular momentum operators as can be easily verified. These are the standard commutation relations of the angular momentum operators. It is convenient to measure the energy in units of 2JN, and work with the Hamiltonian where u ≡ U N J (see (2)). Since (4) are angular momentum operators, and the eigenvalues of For large N the semi-classical limit is justified. In the classical limit, the equation of motion can be obtained by replacing N i [f, g] −→ {f, g} where {f, g} are the Poisson's brackets. These are the Hamilton equations obtained from (8). As the total number of particles N is conserved, the total angular momentum S 2 is conserved as well. Therefore, the vector lies on the Bloch sphere of radius 1 2 and it is possible to write where 0 < ϕ < 2π is an angle circling the S x axis. Now, the Hamiltonian (8) takes the form

III. THE SEMI-CLASSICAL CALCULATION OF THE SPECTRUM
In the absence of inter-particle interactions (u = 0), the bosons undergo Rabi oscillations and the phase space trajectories circle around the S x axis with frequency 2J. We would like to study the dynamics in the Rabi regime u ≪ 1 (weak inter-particle interactions) by using semi-classical methods. Our aim is to find the spectrum of (10) by using WKB quantization for the action variable (a similar approach was adopted by [10] for the Josephson regime 1 < u < N 2 ). S x and ϕ are canonically conjugate variables. Their variation is given by Hamilton's equations generated by H of (8). It was verified that these are identical to the equations satisfied by the components of S. We turn now to calculate the action variable via [45] For this we use the relation between S x and ϕ given by where only the + solution is consistent with (10) for u = 0. In the first order in u, the action can be calculated from and by (11), Now, one can write the Hamiltonian in terms of I as where only the + solution satisfies (14) for u = 0. Therefore, to the first order in u, The action variable is quantized [45] so that where n = − N 2 , ..., N 2 are integers. Note thatφ = − ∂H ∂Sx ≈ −1 for small u and thereforė ϕ never vanishes. Consequently the Maslov index vanishes. Hence, the spectrum of the Hamiltonian (8) is In order to compare the energies (18) to the exact spectrum of the BH model (3) (which can be obtained by diagonalizing the Hamiltonian matrix), we should multiply it by 2JN and add the constants which were omitted in (5) (see (93)) , namely, For u < 1, This spectrum is a good approximation to the exact BH spectrum (see Fig. 1).
In second order in u, one finds: which leads to an action variable of the form In order to find the corrections to the spectrum (18), we substitute H = E (1) n + u 2 · δ n in (21) and keep terms up to the second order in u, resulting in and leading to Numerical calculations ( Fig. 1) verify that the spectrum E (BH2) n is indeed closer to the BH spectrum than E (BH1) n . Although the second order correction is extremely small compared to the first order, it turns out to be of great importance for the dynamics and in particular for the shape of the revival peaks as will be shown in Sec. V. Note that the agreement is very good even for u that is not much smaller than 1. There are predictions based on low order perturbation theory that hold even when the perturbations are not very small [46][47][48][49].
For the present work, it is particularly instructive to note Eq. (A.2) and (A.3) of Ref. [49].
The result of (24) is actually not the one of quantum perturbation theory but the expansion to the order u 2 of the leading semiclassical result. This expansion is convergent in general for u < 1 2 , and for small n used in the paper, it is sufficient that u < 1 as can be seen from (12). In the next section, standard quantum perturbation theory is used and we note that it requires uN < 1. A natural question is what are the corrections to the leading order in the Semiclassical expansion presented here. In App. D it is shown that the correction is of order 1 N 2 . Hence it is of the form 1 N 2 f (H). Then, this term should be added to the RHS of (21) leading to an additional correction to the energy levels. In App. D. we estimate this correction for representative values of the parameters and find it to be extremely small.

IV. THE PERTUBATIVE CALCULATION OF THE SPECTRUM
It is possible to calculate the spectrum of (8) by using standard quantum perturbation theory for small u. The perturbation series is likely to converge if uN < 1 since the energy differences are of order 1 N , see (33). In the first order in u, The matrix element k |S 2 z | n can be calculated easily by using the relation where S ± = (S z ∓ iS y ) given by (4) are ladder operators satisfying Hence, Assuming N ≫ 1, we expand n |S 2 z | n to the second order in 1 N and get resulting in which is equivalent to the semiclassical correction calculated in (18), if 1 N is ignored compared to 1.
The energies to the second order in u are The energy differences are and the matrix elements k |S 2 z | n does not vanish only for k = n ± 2. Therefore, and, which is equivalent to the semiclassical correction calculated in (23), when 1 N is ignored compared to 1.

V. DYNAMICS
In this section, our aim is to derive an analytic expression for the expectation value of S z (t) that is the difference in occupation of the two sites where the initial condition is that all the bosons occupy the state L and S z = 1 2 (it is the north pole of the phase space Bloch sphere). In the framework of the BH model (3), it is possible to calculate S z (t) numerically [4]. The resulting S z (t) is a series of collapses and revivals, superimposed on rapid oscillations. We would like to utilize the spectrum (24) in order to study analytically the dynamics in the Rabi regime u < 1.
In absence of inter-particle interactions (u = 0), the operator S x commutes with the Hamiltonian (8). Hence, for u < 1, the eigenstates of (8) can be approximated by the eigenstates of S x (corrections of higher order will be discussed later), namely by where The reason for (36) is that In what follows, we calculate the evolution of the operator for the initial condition It is useful to expand |ψ (t = 0) in the basis of (36), where For N ≫ 1 and n ≪ N 2 , the binomial coefficients can be approximated by a Gaussian, We note that the normalized difference between the occupation of the two sites is In the basis {|n }, S + is a raising operator, therefore, we used (37), (38) and (39).The expectation value of S + at time t for the initial condition (43) is First, note that According to (24), Since u < 1, we can neglect the term 3 2N 2 u 2 n which is much smaller than 2 N un. For large N and n ≪ N 2 , (48) can be written in the form and We note that e −iφt is a rapidly oscillating function of t with a period that is approximately π J . We turn now to explore the envelope of S. Since n is an integer, in first order in u, the envelope of the sum (48) is a periodic function of t with period (revival time) of Actually T R is the inverse of the coefficient of the linear term in n in the RHS of (50), The estimate (54) assumes u N ≪ 1. The terms proportional to u 2 in (50) are ignored for the same reason, taking into account that in what follows only terms where n ≪ N are important. Around the m-th revival, we write t = m · T R + τ with We approximate the sum by an integral What enables to approximate the sum over n by an integral is the fact that in the vicinity of a revival J N uτ + 3 N 2 u 2 nτ is small (while J N unt is typically large). The integral was calculated in App. B (where we should take β = γ = 1 andm = m for the calculations of the present section) using ψ in the order u 0 and in App. C the corrections of the order u and u 2 were added. In App. E it is verified that the semiclassical wave function gives the same result.
The result is where D R , D D and φ s are given by (98), (100), (106), and R and φ ′ s will be calculated in what follows. and namely, the revivals start to mix at time For times t < m max T R , in the leading order in u, S can be approximated by Therefore, 1 + 9 16 u 2 m 2 π 2 1/4 exp and (see (46)) where and In the expression for φ we neglected 1 N 2 compared to 1 in (52). The other phase variable is neglecting 1 N compared to 1, one finds The evolution of the expectation of the normalized difference in occupation of the two sites is the main result of the present work. In Fig. 2 it is compared to exact results found by numerical diagonalization of the Hamiltonian (3), for u = 1 2 , N = 100 and J = 1. In Fig. 2 as well as in Figs. 3 the expressions (69) and (71) for the phases were used. We checked that if (52) and (70) are used instead, the results cannot be distinguished in the plots. We note remarkable agreement of the envelope with the exact numerical result. The rapid oscillations, exhibit good agreement for short times (Fig. 2(c)) but it deteriorates for longer times (Fig. 2(d)).
In Fig. 3 the evolution of the difference in occupation between the two sites is presented for u = 1 20 , N = 50 and J = 1. We note also the remarkable agreement between the analytical and numerical results found for the envelope. The prediction for the rapid oscillations agrees with the exact results for longer times and more revivals than in Fig. 2.  For short times (m = 0), the dynamics is described by and Both the expectations of S z and S y oscillate rapidly with the Rabi frequency π J , and at the short time scale have a Gaussian envelope which is in the leading order in u and 1 N . Namely, it decays on the time scale Note that correction term J 1 8 u 2 + u N to the phase in (69) improves the agreement with the exact numerical results compared to Rabi's phase 2Jt (see Fig. 2c).
For m > m max the revival peaks overlap and the picture presented in Figs. 2a and 3a is blurred as demonstrated in Fig. 2b.

VI. INITIAL CONDITIONS WHERE BOTH SITES ARE OCCUPIED
It is interesting to study the dynamics of a double well where the initial condition is different occupation of the two wells. Such situation is encountered, for example, if a condensate is suddenly separated into two unequal parts, as was done in [7]. The initial condition is of the form Expansion of (76) as a sum N n=0 c n |n (with |n given by (36)) yields The coefficients c n of substantial magnitude are distributed around so that and where β = 1 cos 2 (2α) .
The expectation value ∆ (t) = S z (t) is calculated in a similar way to what was done in the previous section. The differences are:  (47)) .
3. The β in the exponent of (80) affects the result of the integral S m of (56), see App. B.
4. For n ≈ n max , it is possible that 3 2N 2 u 2 n 2 in (50) is not negligible compared to 2 N un. Consequently, the revival time T R will be modified as described in what follows. We substitute in (50) n = n max + ∆n and write E Therefore, for the initial condition (76), the expectation value ∆ (t) takes the form (as can be seen by modifying (66)), and φ is given by (69).
We turn now to estimate the conditions for the validity of the approximation (83). The width of the Gaussian (80) is N β , therefore it is required that therefore by (78), The result (88) is demonstrated in Fig. 4.
In addition, the perturbation theory in u adds to S m a term of the form unmax N S m (see App. C (120), where the first order correction is calculated). Therefore, the higher orders can be neglected only if un max ≪ N, namely (see (78)), Furthermore, the spectrum (24) is more accurate for small values of |n| (see Fig. 1) where H in (12) is small. If one wants to describe the dynamics for ε 2 1, higher orders in the expansion of (12) might be needed.

VII. SUMMARY AND DISCUSSION
In the present work the dynamics of the two site Bose Hubbard model defined by (2) and (3) were analyzed. We analyzed it for weak coupling u (1) and for large number of particles N. The calculation was preformed to order u 2 and to the leading order in 1 N , using a semiclassical method where 1 N plays the role of the Planck's constant. It is important to note that this is not the standard quantum perturbation theory that requires uN < O (1) but here it is requires only that u ≤ O (1).
In particular, the normalized difference in the occupation of the sites ∆ (t) = Comparison between the approximate result and the exact numerical calculation demonstrates that the result obtained indeed requires the terms in order u 2 and 1 N . The classical approximation (10) reproduces correctly the rapid oscillations for short times. Such a behavior is found also for the GPE in double well [4,11]. Quantization is essential for the collapses and revivals. The collapse and revival times are predicted correctly by the first order in the interaction u, however for the width of the peaks the order u 2 is required, since the width depends on m via the combination m 2 u 2 .
Collapses and revivals were found in various situations [4, 23, 31-34, 36, 38-41, 43]. To the best of our knowledge, (66) is the only complete analytic description of this situation for a specific model of interacting bosons. It is reminiscent of the dynamics of the Jaynes-Cummings model [38].
We studied also the case where initially both sites are populated and found out an approximation that is good if the initial difference in occupation is sufficiently large. For finer details see Eqs. (88), (89) and Fig. 4. In this case, collapses and revivals are found as well where also here the collapse time T c and the width of the reviving peaks are proportional to √ N and the revival time is proportional to N. However, the revival time depends on the initial condition, as is seen from (82).
The generalization to other situations is left for further work.

Acknowledgments
This work resulted of a discussion with D. Cohen on ref. [10]. We thank him for motivating this direction of research and many critical discussions and communications. We The first term of H ′ is identical to the first term of H BH . The second term is We find that Therefore, That reduces to (5) up to a constant.

Appendix B
In this appendix, we calculate the integral (56) and the corresponding integral required in Sec. VI, which are of the form where In Sec. V we consider the case β = 1, γ = 1 while in section VI, β = 1 cos 2 (2α) and γ = 1 − 3 4 u sin (2α) −1 . There, α determines the initial conditions, see (76). In order to write (94) explicitly, we preform some manipulations where for each order of τ , only the dominant order in u is taken into account.
After multiplying the numerator and the denominator by the complex conjugate of the denominator, To the leading order in u, D R can be written as Now we turn to calculate π A which appears in (94). According to (100), where Therefore, and where To the order u 2 , and φ s = u Jτ N + 3 8 The matrix element k |S 2 z | n can be calculated easily by using (29) up to the second order in 1 N . The result is where The energy difference is (33) and therefore For small n and u relevant for the present work, it agrees with the semiclassical result (160).
We would like to expand the wavefunction in basis |n ′ . For this purpose, we define the expansion coefficients According to (112), while the coefficients c n are given by (79) (β and n max are defined by (81) and (78), in Sec.
V, β = 1 and n max = 0, resulting in (45)). Therefore, In the leading order in 1 and in the first order in u, The resulting correction to S m is a of the form where A, B are presented explicitly in App. B, Eq. (96) and x = n − n max . The integral can be solved by usingˆ∞ Therefore, Since A and B are of the same order of magnitude (see (59) and (60) m is negligible. However, in other cases (discussed in Sec. VI) it might be important and then namely, Therefore, Expending the exponent to the second order in 1 N yields The second order correction to S m is of the form where A,B are defined in (96). The integral can be calculated by using (119) and Therefore, and it is a small correction.

Appendix D
In this appendix we calculate higher orders of the WKB expansion and show that its contribution to the spectrum is not important. In the WKB expansion [45], one makes the ansatz ψ ∝ e iS(ϕ)/ where S (ϕ) is the series and Here, = 1 N . In Sec. III, we used the Bohr-Sommerfeld quantization, namely, we demanded S (ϕ) = S (ϕ + 2π) + 2πn in order to find the spectrum. In the present work, = 1 N is understood. There, we replaced S by S 0 which is justified only in the leading order in .
Finally, it turned out that the spectrum contains terms of higher orders of (21)-(24) and therefore, the effects of S 1 and S 2 should be taken into account as well. Fortunately, the contributions of S 1 and S 2 are negligible as described in what follows. S x is periodic in ϕ so that S 1 does not contribute to the spectrum.
In order to find the contribution of S 2 , we first calculate the derivatives of S 1 : and Hence, In what follows, all calculations are performed to the order 2 . We substitute S x of (20) and find and Therefore, to the second order in u, − H 2 sin 2 (2ϕ) + 2u cos (2ϕ) + 8u 2 H 3 sin 2 ϕ cos 2 ϕ − sin 4 ϕ . (146) Assuming u ≪ H we find 1 − H 2 sin 2 (2ϕ) + 8 cos (2ϕ) sin 2 ϕ (147) +2u cos (2ϕ) + 8u 2 H 3 sin 2 ϕ cos 2 ϕ − sin 4 ϕ . Therefore,δ This should be added to the right hand side of (21), resulting in a contribution of to the spectrum (23). The approximation leading to this term is not valid for small n (see (18) where H is of order u). To find an estimate for the correction in this regime we repeat the calculation for I = 0. If n = 0, S x ≈ − u 8 + 1 4 u sin 2 ϕ (151) and the derivatives are ∂Sx ∂ϕ = 1 4 u sin (2ϕ) , ∂ 2 Sx ∂ϕ 2 = 1 2 u cos (2ϕ). Therefore, to the second order in u, This expression is antisymmetric with respect to 2ϕ = π 2 + α → 2ϕ = π 2 − α. Therefore the integral for S 2 vanishes. The above estimates are only for part of the spectrum. Therefore, we turn to a numerical estimate.
In Fig. 5, we present the numerically calculated deviations in the spectrum originating of S 2 and show that it is small for the parameters of Figs. 2-3. The calculation of the spectrum presented in Fig. 5 was carried out by iterations as described in what follows: 1. For each n, S x (ϕ) was calculated according to (12) where H is replaced by the spectrum E (2) n of (23).
2. S 2 (ϕ) was found by substitution of S x in (143) and integration over ϕ.
3. The term 2 2π (S 2 (ϕ + 2π) − S 2 (ϕ)) was added to the RHS of (21), which we solved numerically to obtain a corrected spectrum E (2) n . 4. We repeated steps 1-3 where H in S x is replaced by E (2) n until conversion. 5. We Multiplied the resulting spectrum by 2JN and added the constant C N to be able to compare with the exact BH spectrum.

Appendix E
In this appendix, we calculate the eigenstates in the semiclassical approximation and show that it can be approximated by the eigenstates of S x as was done in Sec. V and VI. According to (138) and (13), in the first order in u, The eigenstates of S x (obtained by substituting (155) with u = 0 in (154)) are e inϕ . These are denoted by |n of (36). The overlap between |k ′ and |n is n|k ′ = 1 2πˆ2 . In order to solve the integral, we expand to series of Bessel functions: and obtain n| (n + 2l) ′ = J l C 2 (158) for positive integer l. Since C 2 is small, the Bessel functions can be approximated by , so that the overlap is substantial only for small values of l and This result reduces to (112) for small n and contribute only small corrections to the dynamics, as was shown in App. C.