Exact dynamics of two ultra-cold bosons confined in a one-dimensional double-well potential

The dynamics of two ultra-cold bosons confined in a one-dimensional double-well potential is studied. We compare the exact dynamics governed by a full two-body Hamiltonian with the dynamics obtained in a two-mode model approximation. We show that for sufficiently large interactions the two-mode model breaks down and higher single-particle states have to be taken into account to describe the dynamical properties of the system correctly.


Introduction
Double-well confinement is one of the simplest examples where quantum dynamics manifests its nonintuitive properties [1,2].Already on the single-particle level, the tunneling through a classically impenetrable barrier has many spectacular consequences like electron tunneling through p-n junctions [3,4,5,6] or the Josephson effect [7].The physics of double-well systems becomes even more interesting whenever interactions between particles are considered.
In view of recent experimental progress with ultra-cold atoms forming a Bose-Einstein condensate, double-well systems are one of the most commonly exploited schemes studied [8,9,10,11,12,13,14,15,16,17,18].Typically, in this context one assumes that weakly interacting bosons occupying different wells can be described with two independent single-particle orbitals and that the dynamics is governed by two mechanisms: contact two-body interactions acting locally and single-particle tunneling between wells.Then, in the mean-field limit, a corresponding Gross-Pitaevskii equation is introduced and numerically solved for different initial conditions [19,20,21].Generalized two-mode models, taking into account additional terms originating from long-range interactions or occupation-dependent tunnelings, are also considered in the literature and relevant corrections to the dynamics are studied [22,23,24].Although the validity of these simplified two-mode models was confirmed experimentally for weak interactions between particles, they were extended beyond the range of their applicability and adopted for strongly interacting systems, i.e. in Correspondence to: tomsow@ifpan.edu.plsituations when the local interaction energy is much larger than the single-particle tunneling energy.For example, it was shown that for initially imbalanced occupations the dynamics is heavily affected by strong interactions [25].Unfortunately, the validity of the model used was not discussed and its predictions were not compared with the exact dynamics governed by a general model.
It is quite obvious that for sufficiently strong interactions between particles any two-mode model has to break down.This comes from the observation that interactions always introduce some multi-particle correlations that exist locally at each site of the potential.Therefore, a simple assumption that all particles occupying a single site can be described with a single effective orbital cannot be valid.The situation is very similar to the problem of ultracold bosons confined in a single harmonic trap.As it was shown in the case of two strongly interacting bosons in a harmonic trap, a single-mode description is not valid [26].In the case studied here, whenever interactions are significantly larger than typical tunneling energies between wells, local correlations induced by interactions are produced much faster than correlations between wells.Just from this simple observation it is quite obvious that any two-mode approximation is inconsistent.
Inspired by this simple observation, in this article we study the dynamical properties of two bosons confined in a one-dimensional double-well potential and initially occupying a chosen site.We numerically compare the exact, many-body dynamics of the system with the dynamics governed by simplified two-mode Hamiltonians.The comparison is performed for different interaction strengths and different depths of the modeled double-well potential.
The paper is organized as follows.In Section 2 the single-particle Hamiltonian with a double-well potential is introduced and its spectral properties are described.In Section 3 we discuss the full many-body Hamiltonian describing identical and spinless ultra-cold bosons and we rewrite it in a two-site multi-orbital Bose-Hubbard form.Then, in Section 4 we focus on the problem of two bosons and we explain the simplifications of the multi-orbital Hamiltonian commonly used in literature.In Section 5 the main results are presented.We compare the evolution governed by a general Hamiltonian with those predicted by simplified models in different regimes of interactions and double-well depths.Finally, we conclude in Section 6.

Double-well model
For clarity, in this article we assume that particles of mass m are confined in a one-dimensional harmonic potential with frequency Ω.The double-well configuration is forced by an additional repulsive gaussian potential with a controlled intensity which is centered in the middle of the trap.The external potential is modeled by the following function: The dimensionless parameter λ > 0 gives the intensity of the barrier between double-well sites (See the left panel in Fig. 1).Depending on λ, the complete spectrum of the corresponding single-particle Schrödinger equation where , is found numerically via an exact diagonalization procedure.The diagonalization is done in the position representation on a dense grid.The resulting spectrum as a function of λ is presented in the right panel of Fig. 1.Obviously, for λ = 0 the standard spectrum of the harmonic oscillator is obtained.For increasing λ the single-particle spectrum changes and a characteristic two-fold quasi-degeneracy between even and odd eigenstates appears.For double-well problems it is convenient to introduce another single-particle basis {ϕ Li (x), ϕ Ri (x)} of states localized in a given well.The basis is constructed directly form the eigenstates of the Hamiltonian (2) as follows: Obviously, the states {ϕ σi (x)}, where the index σ = {L, R} enumerates wells of the potential, are not eigenstates of the Schrödinger equation (2).However, for each state ϕ σi (x) there is only one off-diagonal term coupling the state with the state localized in the other site.The corresponding matrix elements of the Hamiltonian Fig. 1.Left: The shape of the potential V (x) for different values of the parameter λ.For λ = 0 the standard harmonic oscillator shape is restored.For larger λ, two symmetric minima are separated by a potential barrier.Right: The spectrum of the single-particle Hamiltonian (2) as a function of the doublewell parameter λ.For λ = 0, the standard spectrum of the harmonic oscillator is obtained.For larger λ's two-fold degeneracy between even and odd eigenstates appear.In both figures, all the quantities are given in harmonic oscillator units, i.e., energy in hΩ and position in h/mΩ.
are called tunnelings.Moreover, just from the construction it follows that the diagonal terms E i = ϕ * σi (x)H 0 ϕ σi (x) do not depend on the site index σ.The tunnelings J i together with energies E i are related to eigenvalues E i in equation ( 2):

Many-body Hamiltonian
In the following we consider two ultra-cold bosons confined in the potential V (x) and interacting via short-range forces modeled with a point-like potential gδ(x−x ′ ), where g is an effective interaction constant related to the s-wave scattering length.Note that, in contrast to higher dimensions, in the one-dimensional case the Dirac δ function is a well defined self-adjoint operator and any regularization is not necessary [27].Although we consider only two particles, it is still convenient to work in the formalism of the second quantization.Therefore, we introduce the field operator Ψ (x) annihilating a particle at position x.The field operator fulfills natural bosonic commutation relations Ψ (x), Ψ † (x ′ ) = δ(x − x ′ ) and Ψ (x), Ψ (x ′ ) = 0.In this language the many-body Hamiltonian of the system can be written in the form The standard route to analyze the Hamiltonian ( 5) is to decompose the field operator Ψ (x) in a chosen singleparticle basis.Here we decompose the bosonic field in the basis of left-right single-particle states of the double-well potential where an operator âσi annihilates a boson at site σ in level i, i.e. a boson in a single-particle state described by the wave function ϕ σi (x).These operators fulfill standard bosonic commutation relations âσi , â † σ ′ j = δ σσ ′ δ ij and [â σi , âσ ′ j ] = 0.After substitution of the decomposition (6) into (5) one rewrites the Hamiltonian in a Hubbard-like form where for simplicity we have introduced in the interaction part of the Hamiltonian a super-index A = (σ, i) indicating two quantum numbers: site index σ and excitation index i respectively.Values of the interaction terms U ABCD can be calculated straightforwardly from the shapes of single-particle states The Hamiltonian ( 7) is completely equivalent to the initial Hamiltonian (5).To make the notation simpler, in what follows we also introduce operators of the total number of bosons in a given energy level ni and the total number of bosons in a given site Nσ : Note that the operators Nσ are defined very differently to a common definition of a sum over occupations i â † σi âσi .In our definition, the spreading of the basis wave functions ϕ σi (x) to the neighboring well is automatically taken into account.Moreover, the definition used here treats singleparticle superpositions appropriately, i.e. a left/right probability is calculated directly from the full single-particle density and not as a sum of partial probabilities.

The system studied
To study the dynamical properties of the system we assume that initially two bosons occupy the lowest state of a chosen (left) site of the double-well potential It is worth noticing that this kind of state can be prepared experimentally for a few particles.For example, a similar state for two distinguishable particles was obtained in recent experiments [28].Therefore, the problem analyzed here has not only a theoretical meaning but also may have some importance in few-body problems considered recently in the field of ultra-cold atomic systems [29,30].
It is also worth noting that analogous systems were studied recently with engineered optical waveguide lattices [31,32].
The initial state |ini⟩ is not an eigenstate of the manybody Hamiltonian (7) and its time evolution is not trivial.Our aim is to compare the exact dynamics of the state governed by the full many-body Hamiltonian (7) with the dynamics predicted by a simplified model in a two-mode approximation, i.e. when the decomposition ( 6) is cut down to the two lowest single-particle states, i = 0.In this approximation the Hamiltonian (7) has the form where for simplicity we introduce standard notations for relevant parameters: and J = J 0 defined in (4).Note that in definitions (12) we directly exploit the fact that single-particle wave functions are real functions of positions.The additional interactioninduced tunneling controlled by the amplitude T depends directly on the overlap between wave functions localized in neighboring wells.Therefore, it becomes important for shallow barriers.
It is also worth noticing that in the model studied, pair-tunneling and inter-site density-density interaction terms are controlled by the same parameter V .This is a direct consequence of the short-range interactions that are assumed.Therefore, a typical approximation made in this context (see for example [25]) of neglecting pair-tunneling without neglecting density-density interactions seems to be inconsistent.
It is known that in the case of an extended Hubbard model describing bosons interacting via long-range interactions, additional tunnelings, i.e. pair-tunneling and densityinduced tunneling, may play a crucial role and may lead to the appearance of exotic quantum many-body phases in the system [33,34,35].
Typically, on this level of simplification, further approximations are performed and terms controlled by am-plitudes T and V are neglected.Then, the Hamiltonian is reduced to the standard two-site Bose-Hubbard form [10,36] In the case of two ultra-cold bosons the simplified models (11) and ( 13) have a very simple matrix representation.Indeed, in this case the Hilbert space is spanned by three two-body Fock states Therefore, the Hamiltonian (11) can be rewritten in a simple 3 × 3 matrix form The matrix form of the reduced Hamiltonian ( 13) can be derived straightforwardly from ( 14) by setting T = V = 0.It is worth noticing that from a theoretical point of view, the latter model ( 13) is controlled effectively by only one dimensionless parameter U/J.However, from experimental side the two parameters U and J are controlled independently, i.e. the single-particle tunneling J is controlled only by an external potential parameter λ; the on-site interaction U is controlled by λ and a coupling strength g which can be tuned independently.That means that in this approximation the same value of the dimensionless parameter U/J can be obtained for different parameters λ by changing the coupling g.Since U/J is the only parameter in the model, in all these cases the evolution of the system is the same provided that time is measured in the units of h/J.Of course for shallow barriers, when the interaction energy is comparable with the energy gap to higher orbitals, terms neglected in the original Hamiltonian start to influence the dynamics and an evolution for the same ratio U/J starts to depend on λ.This would be the sign that the simplified model breaks down.

The dynamics
To analyze the limitations of the simplifications of the Hamiltonian, we study the dynamics governed by Hamiltonians ( 7), (11), and (13) of the system initially prepared in the state |ini⟩.In our numerical approach we first diagonalize the matrix forms of Hamiltonians in the Fock basis built from single-particle states (3) and we obtain two-body eigenstates and their eigenenergies.In the case of the full multi-orbital Hamiltonian (7) we cut the decomposition (6) at a sufficiently large i max (in practice i max = 30).Then we decompose the initial state |ini⟩ into the eigenstates and we find the time-dependent state of the system |ψ(t)⟩.In principle, an exact form of the state |ψ(t)⟩ depends on the approximation assumed.
The simplest quantity that is accessible in experiments and can be compared for different theoretical models is an expectation value of the site occupation number ⟨ψ(t)| Nσ |ψ(t)⟩.Since the total number of particles is conserved, in practice it is sufficient to calculate the imbalance of populations between potential wells In the following, we compare the dynamics of the system in two different regimes of parameters.To make the comparison clear we compare two different barriers: λ = 3 and λ = 5.In both cases we tune the interaction coupling g in such a way that the ratio U/J represents three different regimes: soft interactions (U/J = 0.1), intermediate interactions (U/J = 1), and strong interactions (U/J = 12).With this strategy, a comparison of different cases is quite natural since simplified Hamiltonian (13) depends only on the ratio U/J (with J fixing the energy scale) and therefore it gives exactly the same dynamics for both λ's (understood as an evolution of coefficients in the decomposition of the many-body state in a Fock state basis).Whereas, in the other cases, there is also an explicit dependence on λ through the other parameters of the Hamiltonian, T and V .
In Fig. 2 we show the evolution of imbalance I(t) in the case of a very deep potential barrier, λ = 5.The solid red line represents the imbalance of particles predicted by an exact many-body Hamiltonian (7).The dotted blue and thin black lines correspond to the simplified two-mode models (11) and (13), respectively.In this case, independently of the strength of interactions between particles, both simplified models recover the exact dynamics appropriately.A small phase shift in oscillations is visible only in the very strong interaction regime.However, the general behavior of the imbalance is reproduced.
The situation changes significantly for a shallow barrier.In Fig. 3 we compare the analogous behavior for λ = 3.It is seen that only in the case of small interactions is the dynamics well captured by simplified models.However, a comparison of the dotted blue and solid red lines suggests that a full two-mode model (11) works quite well even for strong interactions and that it can predict the evolution of the imbalance of occupations in a satisfactory way for long times.Only the small oscillations around the general behavior are not reproduced correctly by this model.For example, for U/J = 12, where the simplified model ( 13) is completely wrong, the model ( 11) is sufficient to describe the oscillation of occupation between different wells with almost adequate frequency.The small discrepancy in the resulting frequency leads only to a slow growth of the phase shift between curves.
This illusory conviction that a complete two-mode Hamiltonian (11) is sufficient to describe the dynamical properties of the system in the strong interaction limit has to be revisited when, instead of densities, inter-particle correlations are considered.For example, let us consider one of the simplest correlations -the probability that bosons occupy different wells of the potential.In the case of two bosons, this probability is related to the density-density Fig. 2. Imbalance of particles I(t) = nL(t) − nR(t) as a function of normalized time in the case of a deep potential barrier (λ = 5).The solid red curve corresponds to the dynamics governed by a full multi-orbital Hamiltonian (7).The dotted blue and thin black curves show the simplified dynamics described by Hamiltonians ( 11) and ( 13), respectively.In all cases we assume that initially the system is prepared in the state |ini⟩.
Note that the exact dynamics is well approximated by the reduced Hamiltonians independently of the strength of the interactions being considered.
U/J=0.Fig. 3. Imbalance of particles I(t) = nL(t) − nR(t) as a function of normalized time in the case of a shallow potential barrier (λ = 3).The solid red curve corresponds to the dynamics governed by a full multi-orbital Hamiltonian (7).The dotted blue and thin black curves show the simplified dynamics described by Hamiltonians (11) and (13), respectively.In all cases we assume that initially the system is prepared in the state |ini⟩.In contrast to the cases shown in Fig. 2, the exact dynamics can not be approximated by the most simplified Hamiltonian (13).However, comparison of the solid red and dotted blue lines may suggest that the dynamics of occupations are well captured by the two-mode Hamiltonian (11).As explained in the main text and in Fig. 4, the two-mode Hamiltonian is not sufficient to describe correlations between particles. correlation: where nσi = â † σi âσi .The time evolution of the probability P(t) may shed some light on the differences between evolutions governed by different Hamiltonians in the strong interaction limit.This comes from the fact that for the simplest case (13), single-particle tunnelings is strongly suppressed by the conservation of energy.The energy difference between the initial state (10) and a state in which bosons occupy different wells is equal to U , and is much larger than the energy gain from the tunneling, J. Therefore, the dynamics is governed effectively by the second-order process in the tunneling, i.e. bosons tunnel between wells mainly in pairs [37].This is visible in the evolution of the probability ( 16) (black curve in Fig. 4) -it is close to 0 at all moments.Situation may change for extended models (11) and (7) where other processes are taken into account.For example, in the case of the extended two-mode model (11) the density-induced tunnelings give additional contributions to single-particle tunneling and they effectively support the breaking of a bosonic pair.On the other hand, a pair-tunneling term supports an ordinary second-order tunneling of a composed pair.After all, as it is seen from numerical results, for strong enough interactions the onsite energy U always dominates over other terms and the probability ( 16) is close to 0 also for the extended twomode model (11) (blue line in Fig. 4).
The full many-orbital model ( 7) opens additional channels for single-particle tunneling.For strong interactions, couplings to higher orbitals, in which single-particle tunnelings J i are large, become relevant.Therefore, the breaking of a bosonic pair is amplified.In consequence, particles can tunnel through the barrier almost independently.This fact is clearly visible in the evolution of the probability P(t) (red line in Fig. 4).The probability, initially being equal to 0, rapidly grows to 1/2 and oscillates around 1/4 during the entire evolution.Consequently, the probability that two bosons may be found in different wells is quite large and cannot be neglected.Moreover, as was explained above, this fact is not reproduced correctly by the simplified models.It is quite natural that for a higher number of particles the situation can be even more complicated.Therefore, one should treat all dynamical results obtained with simplified models with increased care.

Conclusions
In this article, we have studied the dynamical properties of two ultra-cold bosons confined in a one-dimensional double-well potential initially occupying the lowest state of a chosen site.We compare the exact dynamics governed by a full two-body Hamiltonian with two simplified twomode models.In particular, we compared the evolution of particle density and spatial correlations between particles.We show that for a shallow barrier and strong enough  11) and ( 13) (blue and black curves, respectively) the probability is close to 0. In the case of the complete many-orbital Hamiltonian (7) the probability is clearly nonzero, i.e. it is not so rare to find bosons in opposite wells of the potential.
interactions the simplified models break down and the correct multi-orbital description cannot be substituted with a two-mode model even if all appropriate interaction terms are taken into account.The fundamental difference between the exact and two-mode descriptions emerges when inter-particle correlations are considered.For example, the evolution of the probability that both bosons are found in opposite wells of the potential crucially depends on couplings to higher orbitals of an external potential.This fact sheds some light on recent theoretical results and opens some perspectives for further experimental explorations.

Fig. 4 .
Fig.4.Time evolution of the probability that bosons occupy different wells of the potential in the case of a shallow potential barrier (λ = 3) and strong interactions (U/J = 12).For simplified Hamiltonians(11) and (13) (blue and black curves, respectively) the probability is close to 0. In the case of the complete many-orbital Hamiltonian(7) the probability is clearly nonzero, i.e. it is not so rare to find bosons in opposite wells of the potential.