{'fraction_non_alphanumeric': 0.041112109087939906, 'fraction_numerical': 0.012125418469829346, 'mean_word_length': 4.337764218784049, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical results for Josephson dynamics of ultracold Bosons

We study the dynamics of ultracold Bosons in a double-well potential within the two-mode Bose-Hubbard model by means of semiclassical methods. By applying a WKB quantization we find analytical results for the energy spectrum, which are in excellent agreement with numerical exact results. They are valid in the energy range of plasma oscillations, both in the Rabi and the Josephson regime. Adopting the reflection principle and the Poisson summation formula we derive an analytical expression for the dynamics of the population imbalance depending on the few relevant parameters of the system only. This allows us to discuss its characteristic dynamics, especially the oscillation frequency, and the collapse- and revival time, as a function of the model parameters, leading to a deeper understanding of Josephson physics. We find that our fomulae match previous experimental observations.


I. INTRODUCTION
Fundamental issues of non-equilibrium physics of interacting many-body quantum systems and of phase coherence and phase stability, in particular, have a long history. A simple yet relevant model, the two-site Bose-Hubbard Hamiltonian, features phase and fluctuation decay, and also revivals and thus, over the years, many thorough investigations of its quantum dynamics have appeared. Most remarkably, recent experiments involving ultracold Bose gases trapped in an effectively one-dimensional double-minimum potential represent an almost ideal realizations of this fundamental model [1,2], with the fascinating possibility to vary relevant model parameters over a wide range.
A full many-body calculation of the dynamics of an interacting, trapped ultracold Bose gas is only possible for a very small number of particles, even for weakly interacting Bosons.
Most often a mean field approximation in form of the Gross-Pitaevskii equation is applied, which provides good results for low temperatures and for a large number of particles N, if only for a limited time and a limited set of observables. These limits are intensively studied.
Once the field operators are replaced by a c-number field, some truly quantum phenomena, e.g. wave function revivals, cannot be described. The double-well potential provides an ideal playground to analyze these issues. Thus, a purely classical field approach quickly comes to its limits, and the question arises whether semiclassical methods can improve the theoretical treatment of such bosonic systems, allowing us in the future to study more challenging problems whose many-body Schrödinger equation can no longer be solved fully numerically.
A number of articles deal with the discussion of the consequences of the mean-field approximation and many-body quantum corrections [3,4] and the many-body quantum and classical dynamics in phase space [5]. Furthermore, semiclassical methods were applied to the double-well system. In [6,7,13,[23][24][25] a WKB quantization is adopted to analyze the energy spectrum and the wave functions in certain parameter regions.
Despite this fair amount of investigations, it is remarkable to realize that -leaving some fairly straightforward cases aside -no analytical expressions for the relevant dynamical quantities appear to be known. Thus, the purpose of this article is to find a generally applicable analytical description of the population imbalance dynamics of an ultracold Bose gas in a double-well potential by applying semiclassical methods. Since the full quantum dynamics can be determined numerically up to many thousand particles, we are able to compare to exact results. Clearly, the interesting case of very large N → ∞ can no longer be investigated numerically, yet our analytical approach is suited to study this very limit in detail.
At low temperatures a Bose-Einstein condensate in a double-well potential can be described by a two-mode approximation. The corresponding second quantized many particle two-site Bose-Hubbard Hamiltonian is written aŝ with the creation and annihilation operators for a boson in the ith well denoted byâ † i , a i with [a i , a † j ] = δ ij . Thus, the particle number operator of the ith site isn i =â † iâ i . U is a measure for the on-site two-body interaction strength, T is a tunneling amplitude, which in the experiments can be controlled by varying the barrier hight. The tilt parameter δ leads to an asymmetry in the one-particle site energies of the two wells and is used to initiate the dynamics. Note that in the standard notation adapted in Josephson physics we have E J = NT and E C = 4U [11].
It has been shown that the Bose-Hubbard Hamiltonian describes the dynamics of the bosons in the double-well potential properly [8], provided that the interaction energy U is small compared to the level spacing of the trap potential, such that only the two lowest lying modes have to be taken into account. Transverse modes should also be suppressed. It should be mentioned that there are finer descriptions of the two-mode limit that also take into account tunnel coupling energies depending explicitly on the nonlinear two-body interaction term [9]. In this work, however, we restrict ourselves to the standard Bose-Hubbard Hamiltonian (1).
First, there are three qualitatively quite different regimes [10,11] with respect to crucial features of the energy spectra. They are best explained by introducing the which thus separates the Rabi-(Λ < 1) from the so-called Josephson regime for which 1 < Λ ≪ N 2 and the Fock regime with Λ ≫ N 2 .
The Rabi regime is the non-interacting limit Λ ≪ 1, when the system consists of N independent particles leading to an almost harmonic oscillator energy spectrum and thus, after an initial tilt, to plasma oscillations with the known plasma frequency ω p = 2T √ 1 + Λ ≈ 2T [2,15].
In the Fock regime all eigenenergies are grouped in doublets with a quasi-degenerate symmetric and antisymmetric state. Thus, the dynamics of the mean population imbalance follows an extremely slow evolution in time which is called self trapping.
The Josepshon regime combines the two characteristics of the spectrum just discussed.
We distinguish the self trapping regime E > 2NT from the plasma oscillating regime, where E < 2NT holds. In the former, the energy eigenstates appear as doublets again leading to self trapping. In the latter the energy eigenstates correspond to an (an-harmonic) oscillator spectrum and the population imbalance oscillates around zero.
Thus, in the Josephson regime the dynamics will depend on the energy of the initial state. For low energies -the subject of this work -the dynamics undergoes plasma oscillations, for higher energies we see self-trapping, which is beyond the scope of this paper.
In this article we have in mind an experiment as in reference [1], so the double-well system is initially prepared in the ground state ψ 0 of a tilted potential, i.e. δ = 0 in (1).
Then, at t = 0 it is quickly switched to a symmetric potential, i.e. δ = 0. Starting from an initial population imbalance unequal to zero the system is left to evolve in time.
In our paper we first discuss the spectrum using the semiclassical WKB-or Bohr-Sommerfeld quantization. We find a way to systematically obtain an approximate, useful expression for energies in the plasma oscillating regime. In order to describe imbalance dynamics, we need to explore overlap matrix elements in the following section, which we do with the help of the reflection principle. We then apply the Poisson summation formula, which has a long history in semiclassical approaches to quantum dynamics. As a result, we find a useful expression for the time evolution of the imbalance, containing parameters that can be obtained analytically on the basis of the classical Hamiltonian. We then compare exact calculations with our new formula and find remarkable agreement over the whole relevant range of Λ, covering the known Rabi-but also the plasma oscillating Josephson region. In particular, the oscillation frequency, the collapse and revival times are reproduced astonishingly well. We finally discuss the corresponding analytical expressions. It should be noted that the experimentally observed oscillation frequency in [1] of about 40ms follows directly from our formula.

II. SEMICLASSICAL DESCRIPTION
We will follow mainly Braun [12] and his discrete WKB method, as already applied to the double-well problem by Korsch et al. [13]. The two-mode Bose-Hubbard Hamiltonian can be written in the Schwinger spin representation by transforming to angular momentum whereN is the total particle number operator. For fixed N a change from basis |n, N − n to the angular momentum states |l, j is useful, with l = N 2 and j = n 1 −n 2 2 . With w j = 2Ul 2 −2Ul+2Uj 2 +2δj and p j = −T l(l + 1) − j(j − 1), the eigenvalues of the Hamiltonian are determined by an equation of the form as discussed in [12]. By introducing the "coordinate" operator φ = i ∂ ∂j (note [14]), equation (4) can be written as a Schrödinger equation for the function c j with eigenvalue E and Hamilton operatorĤ = w(j) + p(j)e −iφ + p(j + 1)e iφ . In the classical limit the operators turn to canonically conjugate coordinate φ and momentum j (population imbalance), where φ turns out to be the phase difference between the two wells. Since p(j) is a slowly varying function of j in the classical limit (N → ∞) one can replace both p j and p j+1 by p j+ 1 2 and one finds the Hamilton function which can also be found from the mean-field Gross-Pitaevskii functional in the two-mode limit [15,16]. The classical dynamics of the population imbalance and the relative phase (for δ = 0) is then determined by Hamilton's equations of motion: The rich dynamics in this "classical picture" have been studied by several groups [8,17,18], focusing on the differences between the classical and the quantum description of the dynamics [19][20][21]. Clearly a purely classical description cannot picture the collapses and the revivals of the population imbalance, but it is able to shed light on the transition from the tunneling to the self-trapping regime. Recently, the phase space region near the classical bifurcation was also investigated experimentally with ultracold Bosons [22].

A. Semiclassical energy spectrum: Bohr-Sommerfeld quantization
An analytical approach to the energy spectrum relies on the WKB method following Braun [12] and others [13,[23][24][25]. In [13], only the noninteracting case is investigated analytically. In [23,24] the authors concentrate on energies close to the extremal points, and in [25] the case of an attractive gas for the single value of Λ = 1 is studied. We here concentrate on the plasma oscillating regime and aim for solutions over the whole range of For the Hamilton function (6) it is convenient to introduce two potential-energy curves such that the classically allowed energies lie in the region confined by the two potential curves V + and V − . The minimum energy is chosen to be V − (j = 0) = 0. The potential curves display the transition from the Rabi-to the Josephson regime very nicely, as shown in Fig. 1. The energy eigenvalues change from a (non-harmonic) oscillator like spectrum for Λ < 1 to a spectrum with doublets for Λ > 1 due to tunneling, which can be seen from the potential curves. For Λ > 1, V + attains a local minimum which leads to doublets in the spectrum for energies E > V + (0). The deeper the minimum, the bigger this so-called Fock-fraction of the spectrum. Since we are interested in plasma oscillations, the Fock-fraction will not be investigated here, but a semiclassical analysis along similar linesif only more involved -is possible, see for instance [26].
In the WKB approximation the eigenenergies E n are obtained from the quantization where n is the quantum number, and φ(j) is determined by the Hamilton function (6) at fixed energy E (recall that zero energy E = 0 corresponds to S(E = 0) = 0). The integration limit j + is the (positive) classical turning point as obtained from which leads to a quadratic equation in j 2 with solutions Recall that ω p = 2T √ 1 + Λ is the plasma frequency. For the plasma oscillating regime the relevant turning point is j + . Note that j + → 0 for E → 0, while (j 2 − ) approaches the negative constant (j 2 − ) → −(ω p /(2U)) 2 as E → 0. The integral in eqn. (10) can be solved numerically and the results agree very well with the exact quantum results even for quite small numbers of particles as has already been noticed in [13]. It is impossible to solve the action integral analytically without approximation. As we aim at the plasma oscillation regime, we expand in powers of E. First, however, we take the derivative with respect to energy and rescale to find .
with κ 2 = (j + ) 2 /|j − | 2 . Since κ 2 → 0 for E → 0, and 0 < λ < 1, an expansion of (1+κ 2 λ 2 ) −1/2 in powers of κ 2 λ 2 leads to a series in powers of E. The corresponding integrals are known analytically. Finally, a systematic expansion of κ 2n and 1/|j − | in E leads to which is one of the important results of this paper. Apparently, the formal expansion in E is an expansion in the dimensionless parameter The expression on the right hand side clearly shows that our results are expected to be valid in the plasma oscillating regime E < V + (0), irrespectively of the value of Λ. From a simple integration together with the Bohr-Sommerfeld-quantization condition (10) we find The time evolution of the population imbalanceĵ = (n 1 −n 2 )/2 is then with the matrix The dynamics of the population imbalance thus depends on the energy spectrum through the differences E n − E m , and on the matrix A nm , which contains the initial condition and matrix elements φ m |ĵ|φ n . Figure 3 shows the matrix A nm for increasing Λ, obtained from a numerically exact calculation. Due to parity with respect toĵ, A nm is zero for an even number n − m, as can be seen in Fig. 3: 0 = A nn = A nn±2 = A nn±4 + . . .. Among the non-zero matrix elements, there is a strong hierarchy, in particular for small Λ, which will be important later. It is worth noting that in the limit U → 0 (and therefore Λ → 0, for fixed N) the dynamics is well described by a harmonic oscillator. In that case it is easy to prove (the φ n (j) are Hermite polynomials) that only the A nn±1 are in fact different from zero.
Along the diagonals, the matrix elements A nn±k (with k = 1, 3, 5, . . .) have a Gaussianlike n-dependence. This is due to the n-dependence of the overlap c n = φ n |ψ 0 , which will be discussed in the next section. By contrast, the n-dependence of the matrix elements φ n±k |ĵ|φ n is weak. Thus, it is save to assume the form with n-independent parameters d k ≈ φn ±k |ĵ|φn (with the most relevantn), for which, following (20), we expect

IV. SEMICLASSICAL DYNAMICS OF THE POPULATION IMBALANCE
For a semiclassical evaluation of j(t) according to equations (18) and (19) we need semiclassical expressions for E n − E n±k and the overlap coefficients c n . While the spectrum was discussed in section II, we start here with the latter.

A. Reflection principle
The problem to find overlap integrals of an initial wavepacket ψ 0 (j) localized near j ≈ j 0 with eigenstates |φ E of the Hamiltonian with potential V (j) = V − (j) is often encountered in molecular photo-dissociation [27]. The semiclassical solution (reflection principle) states that with some constant c. It is important to note that here the eigenstates are understood to be energy normalized, i.e. φ E |φ E ′ = δ(E − E ′ ), since in typical applications these are scattering states. For the coefficients c n we therefore find c n = φ n |ψ 0 = dE dn φ E |ψ 0 . The normalization condition 1 = n |c n | 2 ≈ dn|c n | 2 yields c = 1/ V ′ (j 0 ), and we get In our calculations, following the experiments, the initial wave function is prepared as the ground state of the tilted trap potential (achieved through the term δ(n 1 −n 2 ) = 2δĵ in the Bose-Hubbard Hamiltonian (1)). In a harmonic approximation near the potential minimum of the tilted potential we find the Gaussian density with j 0 uniquely determined by the tilting strength δ and Clearly, the shape of the initial wave function determines the shape of the c n as a function of n. On closer inspection of equation (24), however, we observe that for the initial state (25), due to the nonlinear relation between E n and n, the coefficients c n are gaussian in E n but not in n.

B. Population imbalance
Having all the ingredients at hand we can now aim at a semiclassical expression for the dynamics of the population imbalance j(t) which we choose to write as with k = 1, 3, 5, . . . taking into account the diagonal structure of A nm as discussed in the last section. Replacing A nn±k by expression (21) and using the Poisson summation formula we find with This rather complicated expression is readily simplified by changing the integration variable from n to E. Further, as only very small k (k = 1, 3) are relevant (see equ. (22)), it is safe to replace E n − E n−k ≈ dE dn k = 2πk S ′ (E) and neglect the k-dependence in the reflection principle, i.e. c * n±k ≈ c * n . Finally, we replace 2πn = S(E) − π according to the semiclassical quantization rule (10). With τ = kt we find This expression, together with equ. (28) is one of the main results of our paper. As we will see, even with further simplifications, the formula captures all essential details of the dynamics, allows for a thorough understanding of decay and revival dynamics, and, most importantly, is the starting point for analytical expressions.
Due to the localization of the initial state ψ 0 (j), the energy integration in (30) is confined to a relatively small interval near E ≈ V (j 0 ), which we assume to be in the plasma oscillating regime (E < V + (0)). Therefore, for the evaluation of the overall phase mS(E) − 2πkt/S ′ (E) we can rely on our semiclassical series expansions (14) and (16).
With a Gaussian initial state as in (25) and expanding the overall phase up to second order around E ≈ V (j 0 ) allows us to take the Gaussian integral and leads us to the analytical result with τ = kt. In the following we want to discuss the structure of this central result. The most important features are the plasma oscillations (ω p ), their collapse (T collapse ) and their revivals (T rev ).
The phaseφ =φ(τ, m) can be ignored for a qualitative discussion -it is a complicated expression and can be found in the appendix. Importantly,φ varies slowly with time and thus needs only be taken into account when quantitative agreement with exact calculations over extremely long time scales is sought.
The parameter A = A(τ, m) = τ · Σ τ − m · Σ m (expressions for the constants Σ τ and Σ m can be found in the appendix) describes an additional slow broadening and decay of the signal. As for the phaseφ, the inclusion of A leads to quantitative agreement with exact calculations as shown later, but need not be discussed further here.
Thus we concentrate on the important plasma oscillations (ω p ), their collapse (T collapse ) and their revivals (T rev ).
The analytical formula for the generalized plasma frequency for arbitrary Λ is For the revival time we find and for the collapse time with being the width of the wavepacket in energy in units of the mean excited energy. Again, a more elaborate discussion of these results will be done in section VI.

C. Simple Rabi limit
In the well studied Rabi limit, i.e. when Λ ≪ 1, our results simplify. In particular, ω p → ω p , T rev → π U , and T collapse → (2g∆V 0 ω p ) −1 . Moreover, only the main off-diagonal contribution k = 1 of the matrix A nn±k needs to be taken into account. Thus, in the Rabi limit, the dynamics of the population imbalance is governed by the simple expression a result that with an appropriate identification of the parameters can also be found in the literature [28].

A. Plasma oscillation frequency
The plasma oscillation frequency was found to bẽ with . For very small Λ, the correctω p approaches the standard plasma frequency ω p , since the constants c 1 and c 2 tend to one in this limit, and the parameter g approaches zero: with V (j 0 ) ≈ ω 2 p j 2 0 /(2NT ) (harmonic approximation of the potential), it is worth writing the latter parameter in the which shows that g tends to zero linearly in Λ for fixed initial imbalance (j 0 /N). However, for increasing Λ the correction terms inω p become more and more relevant, especially for large j 0 , as can be seen from (38). Figure 9 shows a comparison of the plasma frequency [2], this discrepancy becomes by all means relevant. Sure enough, with the parameters given in [2], our formula leads to 2π/ω p = 39ms -which is the experimentally observed value. By contrast, without our corrections one would find 2π/ω p = 30ms.

B. Collapse time
According to (31), the collapse time is given by The expression in front of the brackets (which can be identified with the collapse time in the Rabi regime, i.e. for small Λ) can be approximated as N/(Λω p (2j 0 /N))σ. Thus, assuming Λ and (j 0 /N) are kept fixed so that σ is proportional to √ N, the collapse time is proportional to √ N. This √ N-behavior has been stated before in [25,29]. Our semiclassical formula shows, however, that this statement is only correct for the special case of fixed Λ and j 0 /N, or in the Rabi limit (Λ ≪ 1). In all the other cases, the collapse time depends in a nontrivial way on N through Λ and j 0 /N. Figures 4, 5, 6, 7, and in detail Fig. 8 show that our semiclassical expression for the collapse time is remarkably reliable.

C. Revival time
Following (31), the revival time is Most interestingly, in the Rabi limit it becomes independent of the number of particles and in fact independent of any other system parameter except the interaction strength U.
The revival time was already discussed in [21] where it was found to be equal to 4π, with an interaction strength of 1 4 (considering the different definition of parameters), which we confirm here, in the Rabi limit. Furthermore, it is stated in [21,25] that the revival time grows linearly with the number of particles N. This is obviously true for those investigations with UN = const only, as can be seen from our expression (40). However, note that even only slightly away from the Rabi limit, when Λ approaches or becomes greater than one, the constants c 1 and c 2 and the parameter g become relevant. This can be seen from Fig.   10. Thus, for Λ > 1 no simple scaling law for the revival time exists. The figure shows that for increasing Λ the exact revival times differ strongly from the revival time π/U predicted by the Rabi limit formula. On the other hand, it can be seen that the improved semiclassical expression (40) reproduces the exact revival times very nicely even for Λ > 1. For increasing values of Λ the self-trapping fraction of the phase space is increasing as well, such that for large initial excitations, e.g. j 0 = 20 for N = 100, the semiclassical analysis ceases to give reliable results.

VII. CONCLUSION
We applied semiclassical methods to the well-known two-mode Bose-Hubbard model, in order to investigate in detail BEC tunneling in a double-well trap. Within the plasma oscillation regime we found analytical expressions for the energy spectrum and the initial state agreeing nicely with numerically exact results. Employing the reflection principle and the Poisson summation formula led us to an analytical expression for the time evolution of the population imbalance of the Bose gas in the double well. This allows us to discuss the dependence of characteristic quantities of the dynamics, like plasma oscillation frequency, collapse and revival times, on the relevant system parameters. Remarkably enough, despite a wealth of publications on the two-mode model, such detailed understanding has not been achieved before. Finally, our generalized formula for the plasma oscillation frequency agrees perfectly well with experimental findings. Challenging as it may be, we hope that our predictions for collapse and revival times will be confirmed experimentally, too.
Semiclassical methods are well suited to study the non-equilibrium dynamics of a Bosonic interacting many-body quantum system. For systems with more degrees of freedom, an explicitly time dependent approach might prove useful.
parameters. The phase of the oscillation reads where the dominantly m-dependent part is defined separately as