Dephasing and breakdown of adiabaticity in the splitting of Bose–Einstein condensates

We study the quantum dynamics of a Bose–Einstein condensate trapped in a double-well potential with a rising interwell barrier. We analytically find the characteristic timescales of the splitting process and compare our results with numerical analyses available in the literature. In the first stage of the dynamics, the relative phase of the two condensates evolves adiabatically. At a critical time tad, small-amplitude phase fluctuations around the average trajectory increase exponentially fast, signalling the breakdown of adiabaticity. After this stage, it is still not possible to neglect the tunnelling exchange between the two wells. We find highly non-trivial dependences of the dephasing time on tad and on the ramping time of the interwell barrier.


Introduction
From the early observations of a Bose-Einstein condensate (BEC), there has been a growing interest in theoretical and experimental studies of a condensate in a double-well trap. There are several goals related to these studies: (i) to understand the analogue of the Josephson effect in this type of system; (ii) to clarify the meaning of the phase in quantum mechanics; and (iii) to create interferometers working at the Heisenberg limit. The recent experimental realizations of a stable double-well trap [1,2] is renewing the interest in these topics. Even though residual sources of noise might still be limiting the coherence lifetime of the two coupled systems [1,3], the experimental progress is quite promising for future developments and possible technological applications.
In this paper we study, in a two-mode approximation, the dynamical splitting of a condensate trapped in a double-well while ramping up the interwell barrier. This problem has produced some controversy in the past [4]- [6]. In [4], Javanainen and Wilkens (JW) analysed the condensate splitting in a two-stage model: first, an initial condensate is partially split by ramping a potential barrier. The ramping rate is R I ω p , where ω p is the Josephson plasma frequency, the process remains strictly adiabatic. In the second stage, the barrier is raised to infinity at a rate R II ω p , and the system is eventually left alone for a time t. There has been a debate between JW and Legget and Sols (LS) [5] about the rate at which the two halves of the condensate lose their phase memory in this second stage. The crucial point was the estimation of the ground state phase fluctuations. In [6], Javanainen and Ivanov (JI) clarified the issue by analysing numerically a two-mode model of a continuous splitting.
In the present paper, we study the quantum dynamics of a BEC trapped in a double-well potential with a rising barrier, within an analytically solvable model. We recover simple analytical estimates of the various timescales of the problem, which agree quite well with the numerical analyses. We emphasize that, even after the loss of adiabaticity, it is still not possible to neglect the tunnelling exchange between the wells. As a consequence, the dephasing time, in which the system loses memory about the relative phase between the two condensates, has a complicated dependence on the ramping time of the interwell barrier. This has important consequences when studying the realizability of a BEC interferometer.

The quantum-phase model
In this section, we review the quantum-phase model (QPM) of a condensate trapped in a symmetric two-well trap at zero temperature. The second quantization Hamiltonian of a system of bosons interacting with a δ-pseudo-potential is given bŷ whereˆ is the bosonic field operator, V(z, t) is the time-dependent external double-well potential and g = 4πh 2 a/m is the strength of the interparticle interaction with a being the s-wave scattering length. The two-mode ansatz readŝ where ψ 1,2 (z, t) can be constructed as a sum and difference of the first symmetric and antisymmetric Gross-Pitaevskii dynamical wavefunctions in the double-well trap. The operator a † 1,2 (â 1,2 ) creates (destroys) a particle in the modes 1,2 respectively. In the following, we will decouple the 'internal dynamics' (distribution of the particles between the two modes) and the 'external dynamics' (evolution of spatial wave function of the modes) [7]. The important timescale of the external dynamics is given by the trapping oscillation period τ z = 2π/ω z , where ω z is the trap frequency. If the barrier is raised on a timescale t τ z , then the final wave function ψ 1,2 (z, t) will be very close to the ground state of the 1,2 well, respectively. If t τ z , on the other hand, the splitting process excites the system. Menotti et al in [7] indicate the revival time of coherence, τ r [8] as the adiabaticity timescale for the internal dynamics. In the current experiments, τ z τ r , with τ r longer than the life time of the condensate. Therefore, raising the potential barrier at a rate τ z t τ r , we can decouple the internal and external dynamics. The two-mode model has been extensively discussed in the literature [4]- [10], [12,17]. In general, it is expected to be fairly realistic when the ground state and the first excited state are close to each other in energy and well separated from higher energy modes. It has limited validity in the case of a low potential barrier, when it is not allowed to neglect higher excitation modes and strong atom-atom interactions. However, it becomes increasingly accurate by raising the potential barrier. In this case, in fact, the two lower lying modes become closer in energy and separated from the higher ones.
In the two-mode approximation, the Hamiltonian of the system iŝ The operatorN =n 1 +n 2 =â † 1â 1 +â † 2â 2 is the total number of particles and it commutes withĤ. The quantity E j is the 'Josephson coupling energy', and E c is the 'one-site energy': It is convenient to study equation (3) in the Bargmann phase-states representation [9]. We write a general state in the Hilbert space of the two-mode system as where φ is the relative phase between the two modes, and are non-normalized vectors of the overcomplete Bargmann base, written for the relative number of particles n. In this representation, the action of any operator on |ψ can be represented in terms of differential operators acting on the associated phase amplitude (φ, t). The main consequence of the overcompleteness is the non-standard inner product between Bargmann vectors (7) φ | θ ≈ cos N ( φ−θ 2 ) which affects the inner product between states (6). In the limit E j NE c , the dynamical equation for the 2π-periodic phase amplitude In this paper, we study the dynamical evolution of the wave function (φ, t) when the initial condensate is split with a symmetric double-well potential. This experiment has been recently realized [1], with the two wells ramped apart linearly in time. The distance between the centre of the wells evolves according to d(t) = d 0 + d fin t/ t R , where d 0 and d fin are the initial and final distances, respectively, and t R is the total ramping time. With this setup, it is possible to find both in WKB approximation [11,17] and by a numerical 1D simulation, that the Josephson coupling energy evolves in time with the exponential law E j (t) = E j (0) e −t/τ , where the effective ramping time τ = t Rh / 2m(V 0 − µ)d 2 fin depends on the particle mass m, the initial height of the potential barrier V 0 and the chemical potential µ. In the two-mode approximation, the one-site energy E c remains approximately constant during the dynamics. We notice that E j scales exponentially with the interwell distance only when the condensates are well separated.
During the initial splitting, when the chemical potential is close to the interwell barrier, the dynamics remains adiabatic. The adiabaticity will break down at a large separation of the two condensates, in the deep tunnelling regime, which will be the focus of the next sections.

Numerical solution
We first numerically solved equation (8). Figure 1 presents a plot of the normalized | (φ, t)| 2 at different times. Superimposed on the phase amplitude (blue line), the red line presents the cosine potential in arbitrary units. At t = 0 (figure 1(A)) the phase amplitude (φ, t) is in its ground state, which, for sufficiently high values of E j , is well approximated by a Gaussian. At t > 0, the height of the potential barrier decreases exponentially (here we choose τ = 5 ms) and the phase amplitude spreads (figure 1(B)) until it touches the borders at ±π. The wiggles in figure 1(C) arise from the interference between the two overlapping tails in the region around φ ∼ π, which eventually spread over the full region ( figure 1(D)). The period of oscillations of the interference pattern depends on the velocity of spreading of the phase amplitude: the more adiabatic the expansion, the smaller the number of oscillations.

Variational approach
We study equation (8) with a time-dependent variational approach [12]. The time evolution of the variational parameters is characterized by the minimization of an action with the effective Lagrangian with q i being the time-dependent parameters of the phase amplitude (φ, q i (t)). This provides the familiar Lagrange equations d dt We choose the time-dependent variational phase amplitude with the condition that the width σ φ (t) 2π during the dynamics. At the beginning of the dynamics, the phase amplitude is rather narrow and, feeling only the quadratic part of the cos(φ) potential, it can be approximated by the Gaussian (11). During the first stage of the dynamics, the phase amplitude will follow the instantaneous ground state of the system. A breakdown of adiabaticity will occur at time t ad , which will depend on the ramping time τ, and the initial conditions E c and E j (0). The Gaussian ansatz will fail when the wavefunction touches the borders φ = ±π, namely when σ φ (t) 1. With the variational ansatz (11), the equation of motion for the width of the phase amplitude becomes 3 :

Adiabatic variational solution
We can calculate the adiabatic solution from equation (12) by imposing the adiabatic condition σ φ = 0. In the limit E j E c , the dynamics is governed by the particle exchange through the two wells, obtaining where LW(x) is the Lambert-W-function [13]. In figure 2, we compare the adiabatic solution (blue line) with the numerical one (red line).

Linear approximation
We now solve equation (12) seeking a solution of the form σ φ (t) = σ φ, ad (t) + ε(t). We expect that there will be a time at which the solution ε(t) becomes a significant correction to 3 The differential equations for the parameters σ φ and δ of the variational wave function (11) are [12]: from which we obtain equation (12).  (25)) and the dephasing time t D = 200.6 ms (equation (40)). Here the parameters are τ = 5 ms, E c = 10 −3h kHz and E j = 10 3h kHz. σ φ, ad (t): this will give the criterion for the breakdown of adiabaticity. Replacing σ φ (t) = σ φ, ad (t) + ε(t) in equation (12), where σ φ, ad (t) is the adiabatic solution (13), we obtain a secondorder differential equation for ε(t) where we neglect quadratic terms in ε(t). We are interested in studying this equation for times 0 t t ad assuming that σ 2 ad (t) 1 and, equivalently, E j (t) E c . In this limit, after replacing equation (13) in equation (14), we obtain This is the equation of a harmonic oscillator with a time-dependent (exponentially decreasing) frequency. It is possible to find the solution of this equation in terms of Bessel functions: where J 0 is the zero-order Bessel function of the first kind and Y 0 is the zero-order Bessel function of the second kind [14], while C and D are constants which depend on the initial conditions ε(0), andε(0). Choosing the initial conditions as ε(0) = 0 andε(0) = 0, we have  (23) and (24), respectively.
where A = 4E c E j (0)/h 2 . Figure 3 presents a plot of the functions J 0 (2 √ Aτ e −t/2τ ) and Y 0 (2 √ Aτ e −t/2τ ). Initially, the solutions oscillate with a π/2 phase difference: The amplitude of the oscillations increase exponentially as e −t/4τ . For large times (t 4τ), the functions approach their asymptotic behaviour: J 0 tends to a constant value independently of A and τ, and Y 0 diverges linearly in time The times t J and t Y in figure 3 indicate the last oscillation of the Bessel functions Y 0 and J 0 . They are defined by the equations We interpret this divergence at t > 4τ as the signature of breakdown of adiabaticity. From equations (24) and (23), it is natural to define the time of breakdown of adiabaticity t ad by the relation where c is an arbitrary number ∼1. Since A = 4E c E j (0)/h 2 , E J (t) = E J (0) e −t/τ and introducing the Josephson oscillation frequency ω j (t) = E c E j (t)/h, we can rewrite equation (25) as The explicit expression for t ad is Equation (26) coincides with the definition of breakdown of adiabaticity suggested with a different heuristic argument by Javanainen in [6] with α = 4/c ≈ 2π. This value, found in [6] with numerical calculations, is in reasonable agreement with our analytical calculations. Equation (26) has a simple physical meaning. Initially the system is in its ground state, and we assume E j (0) E c ; the phase dispersion is σ 2 φ (0) = 1/2 E c /E j (0) 1, and the system feels only the quadratic part of the cos(φ) potential. By ramping the two wells, the quantity E j (t) decreases with a timescale τ. As long as ω j (t) 1/τ, the system adjusts itself in such a way that it is always in the ground state. After the time t ad defined by the equation (26), the frequency ω j (t) becomes of the order of 1/τ, and it is impossible for the system to adjust in the ground state following the decreasing of the tunnelling rate. The evolution of the system is no longer adiabatic.

Dephasing time
We define the dephasing time t D as the time needed for the phase dispersion to become of order 1, which roughly implies that the phase amplitude (φ, t) reaches the borders φ = ±π. It might be useful to recall that in a single experiment, a well-defined phase will actually be measured. However, in different experiments, repeated in identical conditions, the measured phases would differ, with a mean-squared fluctuation σ φ . The main result of this section will be the analytical estimate of t D . We first notice that (comparing the red and blue lines in the figure 2) we can approximate well the numerically exact phase amplitude with a Gaussian even for t > t ad . We can therefore expect that the Gaussian variational ansatz can give a good estimate of t D .
Thus, we approximate the QPMĤ φ with the Hamiltonian of a harmonic oscillator with a time-dependent frequency: The quantum evolution of the system governed by equation (28) can be calculated exactly in the Wigner phase space. We first need to calculate the solution of the classical equation of motion where ω 2 j (0) = E c E j (0)/h 2 is the Josephson plasma frequency. The solution of equation (29) can be expressed in terms of Bessel functions of the first and the second kind. We have (similar to equations (15) and (16)): where C and D depend on the initial conditions φ 0 ≡ φ(0) andφ 0 ≡φ(0). We introduce the quantity χ ≡ 2ω j (0)τ. Depending on χ we can distinguish between three regimes: (i) if χ 1, the dynamics of the system will be approximately free; (ii) if χ ≈ 1, the adiabaticity will be broken almost immediately; and (iii) if χ 1, the initial dynamics of the system will be adiabatic, while the breakdown of adiabaticity will occur at time (27). We can write where the quantities J n (χ) and Y n (χ) are Bessel functions of the first and second kind (of degree n), respectively, and is the Wigner transform of the phase probability amplitude (φ, t). For an harmonic oscillator, the (normalized) Wigner transform [15] becomes where σφ 0 = E c /(2hσ φ 0 ) = σ φ 0 ω j (0) is the width of the distribution P(φ 0 ,φ 0 , t) in the variablė φ 0 . Notice that, in equation (33), the integral is performed between ±∞ and not ±π. This is consistent with the condition σ φ 0 , σφ 0 2π. Since the Wigner transform (33) is symmetric in φ 0 andφ 0 , we have that φ(t, φ 0 ,φ 0 ) = 0, and where A direct study of the quantities (36)-(38) in the variable χ shows that, for χ 1, CD is negligible with respect to C 2 and D 2 . Notice that in the limit τ = 0, we recover the free evolution dynamics as given by equation (42). As shown in figure 3, the harmonic oscillator approximation remains valid even for t > t ad . The dephasing of the phase amplitude occurs at t D > τ. We can therefore use the asymptotic expansion t τ for the Bessel functions J O (2ω j (0)τe −t/2τ ) → 1 and Y O (2ω j (0)τe −t/2τ ) → 2 π ln(2ω j (0)τ) − t πτ . We can rewrite equation (35) as The dephasing time t D is defined by the relation σ φ (t D ) = 1 and it can be obtained from equation (39) that Figure 4 presents the times t D and t ad as a function of τ. After a time t = t free , we approximate the evolution of σ φ (t) by the free evolution given by equation (42). The red line represents the numerical calculation corresponding to t free = ∞ which physically corresponds to an unlimited ramping of the two potential wells. The other lines refer to different choices of finite t: t free = t ad = 34.5 ms (blue line); t free = 100 ms (green line); and t free = t D = 200 ms (yellow line).

Free evolution
The free evolution can be studied, in our variational model, by neglecting the potential energy term in equation (12), which becomes The solution of this equation is given by where σ 2 φ (t free ) is the Gaussian width at t = t free when the free evolution takes place. For t t free + 2σ 2 φ (t free )h/E c , we have that σ φ (t) is almost constant, while for t t free + 2σ 2 φ (t free )h/E c , the width σ φ (t) is a linear function of time. In figure 5 we present the evolution of σ φ (t) for different values of t free . In this figure, τ = 5 ms, E c = 10 −3h kHz and E c = 10 3h kHz.
In the model of Legget and Sols [5,16], the free evolution takes place with the breakdown of adiabaticity. This model assumes that t free = t ad . We can check this model by looking at figure (5): the red line represents the numerical solution of equation (8), for τ = 5 ms; the blue line represents the free expansion occurring at the breakdown of adiabaticity t ad as defined by equation (25). We can see that the model of Legget and Sols overestimates the dephasing time in the limit of slow ramping while in the limit of fast ramping, when t ad → 0, it gives a good prediction of t D .
In [1] the two potential wells are separated with a ramping time t free . At the end of the ramping, the condensates are held in the trap for a time t hold . In order to study the phase dispersion during this holding process, we can use the variational approximation (12) with a constant parameter . If we have t free = βτ, where the β 1, we can neglect the potential energy term in equation (12), and we obtain that the phase width evolves freely during the holding time. Figure 5 presents the behaviour of σ φ for different ramping times: t free = 100 ms (green line) and t free = 200 ms (yellow line).

Repeated interference experiments
It has been shown in [18]- [20] that the phase of a condensate is established by measurement. For example two BECs, initially in number 'Fock' states, will interfere and have a definite relative phase in each single experiment. A different phase, however, will be measured in different experiments, so that, by averaging over the ensemble, no interference is observed.
Here we consider a similar case. We start with a single condensate which is split by the confining trap. Even if the relative phase of the two condensates has large quantum fluctuations, in a single experimental run, the two condensates will show a well-defined phase. In many repeated experiments, prepared with identical initial conditions, we will measure different values of the phase, with uncertainty σ φ .
We consider the interference between two independent (not overlapping) condensates in the Thomas-Fermi (TF) approximation. Initially, the two condensates are in equilibrium, the condensate 1 being centred at x = −d/2 and the condensate 2 at x = d/2. We assume that at t = 0, the order parameter is described by the linear combination where 1 ( r) and 2 ( r) are the equilibrium wave functions (order parameters) of the two condensates, respectively, which we assume to be well separated in space: At t = 0, the trap is switched off. We neglect the interaction between the two condensates but we account for the atom-atom interaction in each single condensate, which is important during the free expansion [21]. The total density ρ = | | 2 of the overlapping condensates exhibits modulation of the form where 1,2 ( r, t) = ρ 1,2 ( r, t) e S 1,2 ( r,t) . By recalling the quadratic behaviour of the phase in the TF approximation, we obtain asymptotically [21,22] In a single experiment, the interference pattern is characterized by straight line fringes which are orthogonal to the x-axis (radial axis of the two cigar-shaped parallel condensates). We obtain fringes perpendicular to the x-axis with a spacing of ht/md between two consecutive fringes. To experimentally test the quantum phase dynamics, it would be necessary to average over several identical interferometric realizations. The relative phase will be chosen randomly with a Gaussian distribution of width σ. The ensemble-averaged density is therefore ρ( r, t) = +π −π dφ ρ( r, t; φ)| (φ, t)| 2 = 1 2 [ρ 1 ( r, t) + ρ 2 ( r, t) + ρ int ( r, t)], with ρ 1,2 ( r, t) ≡ |ψ 1,2 ( r, t)| the densities of each released condensate, and ρ int ( r, t) = 2 ρ 1 ( r, t)ρ 2 ( r, t) +π −π dφ cos md ht x + φ exp −(φ 2 /2σ 2 ) √ 2πσ 2 (47) = 2 ρ 1 ( r, t)ρ 2 ( r, t) cos md ht x e −σ 2 /2 .
The ensemble-averaged density given by equation (47) is plotted in figure 6 for different values of the phase uncertainty σ. For σ 1 the interference pattern is washed out.

Conclusions
We have studied the dynamical splitting of a BEC in a two-mode approximation by solving the QPM within a variational ansatz. We have considered the case of a linear ramping of the potential wells, which corresponds to an exponential decrease of the Josephson coupling. In the limit of a very fast ramping, we recover the Legget and Sols results ( [5,16]), in which the phase spreading evolves almost freely. For slow ramping, the adiabaticity is broken when the inverse of the Josephson oscillation frequency becomes of the order of the variation rate of the Josephson coupling. This happens when the system can no longer adjust to follow the separation of the wells. We analytically found the timescale t ad characterizing the breakdown of adiabaticity, equation (27), and compared it with the previous numerical estimate. Even for times longer than t ad , the tunnelling between the wells cannot be neglected and the phase amplitude dynamics cannot be solved within a free-evolution model, making the calculation of the dephasing time highly non-trivial. We have developed a harmonic oscillator model with time-dependent frequency leading to an analytical estimate of the dephasing time, equation (40), which is in very good agreement with the numerical simulations. The results of the paper are summarized in figure 2, where we compared the numerical solution of the QPM with both the variational adiabatic solution and the harmonic oscillator approximation. Finally, we have shown how the quantum phase uncertainty manifests itself experimentally. Given an ensemble of systems prepared with the same initial conditions, we have shown how the visibility of interference fringes is degraded when the phase dispersion is ∼1. Our analysis has been mainly motivated by the need to refine and partially reconsider the previous understanding of the splitting of a BEC, in the light of the recent experiments [1,2].