Adiabatic Splitting, Transport, and Self-Trapping of a Bose-Einstein Condensate in a Double-Well Potential

We show that the adiabatic dynamics of a Bose-Einstein condensate (BEC) in a double well potential can be described in terms of a dark variable resulting from the combination of the population imbalance and the spatial atomic coherence between the two wells. By means of this dark variable, we extend, to the non-linear matter wave case, the recent proposal by Vitanov and Shore [Phys. Rev. A 73, 053402 (2006)] on adiabatic passage techniques to coherently control the population of two internal levels of an atom/molecule. We investigate the conditions to adiabatically split or transport a BEC as well as to prepare an adiabatic self trapping state by the optimal delayed temporal variation of the tunneling rate via either the energy bias between the two wells or the BEC non-linearity. The emergence of non-linear eigenstates and unstable stationary solutions of the system as well as their role in the breaking down of the adiabatic dynamics is investigated in detail.


I. INTRODUCTION
Bose Einstein condensates (BEC) in double well potentials have drawn a lot of attention both theoretically and experimentally for the possibilities they offer to study fundamental quantum mechanical effects at the macroscopic level as well as for potential applications like interferometry, high precision measurements or thermometry [1]. Most of the experimental realizations lead to the splitting of the condensate into two independent parts. Nevertheless, recently, also weakly linked parts of a BEC in a double well potential forming a single Josephson junction [2] have been achieved [3]. In contrast to Josephson junctions realized in superconductors and superfluids [4], in BEC the non-linear interatomic interactions play a crucial role. In the presence of the non-linearity two dynamical regimes have been predicted: (i) anharmonic Josephson oscillations [5], if the initial population imbalance of the two wells is below a critical value and (ii) macroscopic quantum self-trapping [6] i.e, inhibition of large amplitude Josephson oscillations above the threshold for the population imbalance. This threshold corresponds to the population imbalance for which the difference between the two on-site interaction energies becomes larger than the tunneling energy splitting. Both dynamical regimes have been explored experimentally in a single Josephson junction [3] and in arrays [7].
Recently, a lot of attention has been devoted to explore techniques to coherently control the non-linear dynamics of a BEC in double well potentials by modulating in time either the potential [8] or the non-linearity [9]. The two main studied regimes are the non-linear Landau-Zener [10][11][12][13][14][15][16] and the Rosen-Zener [17] regimes. In the former the tunneling rate is fixed while the energy bias is linearly varied in time, in the later the energy bias is varied in time while the tunneling rate is switching on/off following e.g., a temporal Gaussian profile. Both these regimes have been deeply investigated in double [10][11][12][15][16][17] and triple well [13,14] potentials yielding a wide variety of dynamical scenarios ranging from robust population transfer to quantum blocking, non-linear oscillations and, even, breaking down of the adiabatic dynamics.
Following a different perspective, there have been several recent proposals to coherently manipulate single atoms [18][19][20] and BECs [13,21,22] in triple-well potentials by adiabatically following a particular energy eigenstate of the system, the so-called spatial dark state. This spatial dark state only involves the two ground states of the extreme traps in a close analogy to the well known quantum optical Stimulated Raman Adiabatic Passage (STIRAP) technique [23]. Accordingly, these tools were named Three-Level Atom Optics (TLAO) techniques [18]. In this paper, following the work by Vitanov and Shore [24] on the adiabatic passage on two-level atoms, we extend the TLAO techniques to the two-level matter wave case showing that the adiabatic dynamics of a BEC, in a double well potential, can be described in terms of a dark variable resulting from the combination of the population imbalance and the spatial atomic coherence. This dark variable can be tailored varying in time the tunneling rate, the energy bias, and the non-linearity with the goal of (i) adiabatically splitting of the BEC, (ii) achieving complete BEC transfer from one well to the other; and (iii) preparing of an adiabatic self-trapping state.
The paper is organized as follows. In Section II we present the physical system consisting of a BEC in a double-well potential. We assume the two-level approximation and describe the BEC dynamics in terms of a dark variable. The conditions for the adiabatic control of tunneling by means of this dark variable are discussed in Section III from a non-linear dynamics perspective. In Section IV we present detailed numerical simulations on the adiabatic splitting, transport, and trapping of a BEC. Section V summarizes the main conclusions of the paper and briefly discusses the validity of the two-mode approximation. It also presents some possible extensions of the present work such as its formulation in second quantization.

II. MODEL
We consider a BEC trapped in a double well potential ( Fig. 1(a)), whose dynamics at zero temperature is described by the time-dependent Gross-Pitaevskii equation where V ( r, t) is the external trapping potential and the non-linearity is given by g = 4N πa s /m, with N the total atom number, a s the s-wave scattering length, and m the atomic mass. Assuming the two-level approximation [6,25], the wave function or classical order parameter of the BEC under study can be written as ψ( accounts for the ground state of the corresponding isolated trap. The amplitudes c L and c R obey the non-linear two-mode dynamical equations given by: with the Hamiltonian: where U R,L are the atomic self-interaction energies, Ω is the tunneling rate between the two wells, and ǫ R,L are the on-site energies. In terms of the φ L,R ( r) overlaps, the former parameters read: . Notice that the two-mode approximation assumes that the parameters of the problem (tunneling rate, non-linear interaction, and energy bias) can be varied independently. In general, however, the temporal modification of any of these parameters will result in the modification of the spatial mode functions affecting the rest of the parameters. From the experimental point of view, a double-well optical potential can be created with well separations on the range of few microns by means of two focused laser beams [26], or by the superposition of a 3D Three-level correspondence of (a) for the density matrix variables and coupling strengths. w is the population difference, u/2 (v/2) is the real (imaginary) part of the spatial coherence, and U is the non-linear self-interaction energy.
(h = 1) crossed beam dipole trap with a 1D optical lattice [3]. In the later, the intensity of the standing wave field that creates the optical lattice and its displacement with respect to the center of the dipole trap could be used to manipulate both the tunneling and the energy bias. In addition, the scattering length could be controlled playing around of a Feshbach resonance [27].
Within the density matrix formalism, assuming U R ∼ U L (≡ U ), and taking ǫ ≡ ǫ R − ǫ L as the energy bias between the two wells, the coherent dynamics of the BEC in the two-well potential can be written as follows [28]: i and σ ij = c i c * j with i, j = L, R the corresponding trap populations and spatial coherence (see Fig. 1(b)). Note that the conservation of the norm implies w 2 + u 2 + v 2 = 1.
Inspired by the work of Vitanov and Shore [24] on STI-RAP in a system of two internal atomic levels, we extend their proposal to the external degrees of freedom of matter waves. It is easy to verify that Eqs. (7) are equivalent to the discrete non linear Schrödinger equation (DNLSE) for a BEC in a triple well potential (see, for instance, Ref. [13]) under the following identifications: where C i with i = L, M, R are the probability amplitudes to find the BEC in the left, middle, or right trap, respectively. Ω LM and Ω MR denote, respectively, the tunneling interaction between the left-middle and middle-right traps. Note however that for the two level system the variables u, v and w are real quantities, while for a three level system the variables C i are, in general, complex numbers.
In the TLAO techniques within the three-level approximation being the left and right levels resonant, the transfer process is based on adiabatically follow-ing one of the three energy eigenstates of the system |D(Θ) = cos Θ |L − sin Θ |R , with the mixing angle defined as Θ(t) = tan −1 [Ω LM (t)/Ω MR (t)], and |L and |R being the ground states of the left and right wells. In the two-level system under investigation, the analogue state will be given by a combination of the population difference and the coherence d(θ) = cos θ × w − sin θ × u, with the mixing angle given by ]. This analogy opens the possibility to extend the TLAO techniques to two-level systems, by appropriately engineering the time dependence of the tunneling rate Ω(t), the energy bias ǫ(t), and/or the non-linear interaction U (t). Note that the non-linear interaction parameter can be modified in time by the temporal variation of the scattering length a s using either magnetic [30] or optical [31] Feshbach resonances or varying the trap frequency, leading to a modification of the BEC spatial profile according to Eq. (5). For a review on the manipulation of Feshbach resonances see [32].

III. ADIABATIC CONTROL OF TUNNELING
In this section we study different scenarios for the coherent control of the external degrees of freedom of a BEC by adiabatically following the dark variable d(θ). In particular, we show how to adiabatically split, transport, and inhibit tunneling of a BEC via the matter wave analogues of both STIRAP and double-STIRAP techniques in two-level systems.
Let us assume that the BEC is initially prepared in the left trap with Ω(t = −∞) = 0. If so, note then that θ = 0 meaning w = −1 and u = v = 0 will be the initial state. Then, to coherently split the condensate, θ should vary adiabatically from θ = 0 to θ = π/2 corresponding to u = 1 and w = v = 0 (see Fig. 1(b)) by appropriately changing Ω, ǫ, and U in time. To transfer the BEC from the left to the right trap, the mixing angle should be slowly increased up to θ = π. For the adiabatic selftrapping state case, the evolution will consist in varying θ from 0 to e.g., π/2 and then back to 0. In all cases, in order to guarantee that the BEC follows the dark variable d(θ) during the whole process two conditions must be fulfilled: and U (t) should be smoothly varied in time to adiabatically follow the dark variable d(θ) which, in turn, means that the mixing angle θ(t) = tan −1 2 [Ω(t)/(ǫ(t) + U (t)w)] must be slowly changed from θ(t = −∞) = 0 to its expected final value. Gaussian profiles for the temporal variations of the control parameters are assumed: with either the energy bias ǫ(t), or the BEC non-linearity U (t), preceding the tunneling interaction Ω(t), i.e., the counter intuitive sequences ∆t ǫ > 0 or ∆t U > 0 will be assumed. n ǫ,U = 0, ±1 is a switch that takes the value 0 for the two-level matter wave analogue of the STIRAP sequence and +1 (−1) for the symmetric (antisymmetric) double-STIRAP sequences. Adiabaticity means that, at any time,|θ(t)| should be much smaller than the energy separation between the selected eigenstate and the one energetically closest. For a weak non-linear interaction, |U (t)/Ω(t)| ≪ 1, this adiabaticity condition readṡ is not a parameter of the system but a dynamical variable since it contains the population difference w in its definition.

Condition 2: Avoiding adverse bifurcation points.
For large enough values of the non-linearity, the interaction between the atoms of the BEC results in a non-linear temporal coupling producing additional non-linear stationary states yielding loop structures and a rich variety of level crossing scenarios [10,13,33,34]. In the matter wave STIRAP case for a triple-well potential [13], it has been shown that even in the adiabatic limit given by|θ(t)| → 0, the appearance of non-linear stationary states breaks down, in some cases, the adiabatic evolution. To address this issue in the double-well potential, we start first looking for those critical U values giving rise to non-linear stationary states by solving the eigenvalue equation: with H given in Eq. (3) and µ being the chemical potential. After some algebra one obtains the following fourth-order eigenvalue equation for µ: For ǫ = 0 and U 2 /(2Ω) 2 < 1 (> 1) this quartic equation gives two (four) real roots [10]. Additional information on the BEC dynamics can be obtained looking for the stationary solutions of the density matrix equations (7) and analyzing their stability (see Appendix A). The stationary solutions can be written as {u ss , v ss = 0, w ss } with u ss /w ss = tan θ ss and (u ss ) 2 + (w ss ) 2 = 1. The energy of these stationary solutions, up to four, is given by Eq. (12). Note that, as the dynamics is conservative, the eigenvalues of the linear stability matrix must satisfy Σ 3 i=1 λ i = 0 for each stationary solution. In fact, a Linear Stability Analysis (LSA) around any of these solutions yields one eigenvalue equal to zero together with either two pure imaginary eigenvalues corresponding to a stable fixed point (or center) or two real eigenvalues accounting for an unstable saddle point. In the limit of a large non-linear interaction, i.e., for U 2 /(2Ω) 2 > 1, the system presents four stationary solutions, one of them being an unstable saddle point while the other three are stable elliptical ones.
If the number of stationary solutions remains constant during the whole adiabatic dynamics, the system will successfully follow the selected energy eigenstate. However, depending on the interplay between the non-linear interaction and the tunneling rate, during the dynamics the number of stationary solutions will change through a bifurcation point from four to two or vice versa. If the selected stationary solution is not involved in the bifurcation, the system will again adiabatically follow the corresponding eigenstate. On the opposite case, whether the adiabatic dynamics will be affected will depend on the particular bifurcation scenario: Case 1. From 4 to 2 stationary solutions through a backward pitchfork bifurcation point in which the stationary solution to be followed merges two more solutions (one unstable) and yields one stable solution. In this case, the system will evolve adiabatically. In what follows, this case will be termed PB42 case, even if the pitchfork bifurcation is perturbed and becomes imperfect. Case 2. From 2 to 4 solutions through a forward pitchfork bifurcation point in which the stationary solution to be followed becomes unstable and yields two stable solutions corresponding, in the limit of vanishing tunneling, to w = ±1. In this case, the system after reaching the bifurcation point will split into a combination of the two stationary solutions. In this case, named PB24 in what follows, the adiabatic dynamics will, in general, break down. Case 3. From 4 to 2 solutions through a saddle-node bifurcation in which the stationary solution to be followed is annihilated by an unstable solution. The adiabatic dynamics breaks down even in the adiabatic limit. After reaching the bifurcation the system will split into a combination of two non-degenerate energy eigenstates and, therefore, will oscillate at the corresponding energy splitting. Case SNB42 in what follows.

IV. NUMERICAL SIMULATIONS
This section is devoted to the numerical simulations of the BEC adiabatic splitting, transport, and self-trapping by means of the temporal variation of the energy bias, the non-linearity, and the tunneling rate. We will show, for the most relevant cases, the eigenvalues and stationary states predicted by Eqs. (12) and (7) as well as their linear stability. The main goal of these numerical simulations will consist in illustrating the different dynamical scenarios to adiabatically control the BEC while looking for those parameter values that prevent the system to reach the two unwanted bifurcation cases PB24 and SNB42 previously described.  (7). In (b) and (c), the dashed curve accounts for the unstable solution, while the short arrows indicate the adiabatic solution (in red) to be followed by the system. PB42 and PB24 denote the two pitchfork bifurcations that the system suffers at Ω0t = 777 and Ω0t = 1750, respectively. Note that the pitchfork bifurcation PB24 yields one unstable solution (dashed curve) plus two energetically degenerated stable solutions (solid curve) corresponding, in the limit of vanishing tunneling, to w = ±1.

A. Adiabatic Splitting of a BEC
We start first by fixing the non-linear interaction control parameter U to a constant value U g while, counter intuitively, time-varying the energy bias and the tunneling rate (see the inset of Fig. 2(b))) with the goal of achieving equal population in the two wells, i.e., w(t out ) = 0 starting with the BEC being in the left trap, i.e., w(t in ) = −1, or, equivalently, θ(t out ) = π/2 from θ(t in ) = 0. Fig. 2(a) shows w(t out ), after the numerical integration of Eqs. (7) with w(t in ) = −1 and the rest of parameters values given in the figure caption. Clearly, there is a large set of parameters where the adiabatic splitting process takes place with a high fidelity (see the plateau in Fig. 2(a)). Fig. 2(b) shows the plateau w(t out ) = 0 for U g /Ω 0 = (−0.12, 1.5) and ǫ 0 = 1.5Ω 0 . To give more insight into the dynamics of the system, we have divided the nonlinear interaction range studied in five regions, from a to e, as shown at the top of Fig. 2(b). In the following we will discuss in detail the dynamics for each of these regions. Note that in all cases, the dynamics starts and ends with Ω(t in ) = Ω(t out ) = 0 which, for U = 0, implies that for t = t in,out the system presents four fixed points or stationary solutions. For intermediate times and large enough tunneling rates, there are only two stationary solutions. In Fig. 3 we show the temporal dynamics corresponding to the parameter values of case (i) in region (c) of Fig. 2(b) with U g = 0.5Ω 0 (see Fig. 3(a)). The adiabatic energy eigenvalues and stationary solutions w ss are shown in Figs. 3(b) and (c), respectively, with the dashed curve accounting for the unstable ones. Short arrows in Figs. 3(b) and 3(c) indicate the stationary solution to be adiabatically followed. The dynamical variable w(t), after integration of Eqs. (7), is plotted in Fig. 3(d). For this set of parameter values the adiabaticity condition is fulfilled and the selected eigenstate is involved in a single bifurcation, at Ω 0 t = 777, corresponding to the PB42 case previously described. After the bifurcation point, the system follows the eigenstate that, at the end of the process, yields w = 0 and µ = 0.25Ω 0 . This non-linear dynamical scenario holds for the whole plateau, region c in Fig. 2(b), where the adiabatic splitting succeeds.
For the parameter range U g /Ω 0 = (1.5, 2.07) corresponding to region d in Fig. 2(b), the final state of the BEC strongly depends on the parameter values. For a slight modification of U g , the final state alternates between w ∼ +1 and w ∼ −1 accounting for the BEC being located at the right or left trap, respectively. To understand the origin of this behavior, we plot in Fig. 4 the detailed dynamics for U g = 1.7Ω 0 corresponding to case (ii) in Fig. 2(b). At Ω 0 t = 1374, the solution that the system is adiabatically following is annihilated with an unstable one through a saddle-node bifurcation, i.e., case SNB42 previously described. At this point, the system splits into two non-degenerate components (being the largest one the energetically closest stable solution) and starts to oscillate at the frequency corresponding to the energy separation of the two solutions. At Ω 0 t = 1588, the system reaches a bifurcation point of the PB24 type. This bifurcation yields two energetically degenerated stable solutions corresponding to the two final states w = ±1. Thus, at Ω 0 t = 1588, the largest component of the system chooses one out of the two new stable solutions. We have numerically verified that the selected solution strongly depends on the previous oscillatory dynamics.
For |U g | > 2.07Ω 0 corresponding to regions a and e of Fig. 2(b), there are four stationary solutions during the whole dynamics that do not present any bifurcation. In this case, the non-linear coupling forces the adiabatic evolution from θ(t in ) = 0 to θ(t out ) = 0 instead of evolving to the desired θ(t out ) = π/2. Then, the initial (u, v, w) = (0, 0, −1) and the final state of the system coincide.
Finally, for −2.07Ω 0 < U g < −0.12Ω 0 , region b in Fig. 2, the energy eigenstate to be followed reaches, during the ramping down of the tunneling interaction Ω(t), a pitchfork bifurcation of the PB24 type. After this bifurcation point, the two stable stationary solutions are energetically degenerated and, therefore, the system splits into a non-oscillatory combination of them that for U g = −0.12Ω 0 (U g = −2.07Ω 0 ) yields w(t out ) = 0 (w(t out ) = −1).
Up to now, we have discussed the possibility to adiabatically split a BEC by fixing U and time varying ǫ(t) and Ω(t). Note, however, that it is also possible to adiabatically split the BEC by appropriately time varying the non-linear interaction U (t) and the tunneling rate Ω(t) with a fixed energy bias. Thus, for U g = −0.4Ω 0 , U 0 = 1.5Ω 0 , ǫ = 0 and the rest of the parameters given in the figure caption, Fig. 5 shows w(t out ) as a function of the time delay between the modulation of the tunneling rate and the non-linear interaction. We again obtain a plateau of parameter values, Ω 0 ∆t U ≃ (250, 600), where the adiabatic splitting, corresponding to w(t out ) ∼ 0, becomes successful. A detailed analysis of the temporal dynamics reveals that, as in Fig. 2, the plateau corresponds to the case where the selected solution reaches during its dynamics a single bifurcation point of the PB42 type.

B. Adiabatic Transport
To adiabatically transfer the BEC from the left to the right trap, the mixing angle must be slowly varied from θ(t in ) = 0 up to θ(t out ) = π. Taking n ǫ = −1 and the rest of parameters as in Fig. 2, Fig. 6 shows the population difference w(t out ) at the end of the antisymmetric double STIRAP process (shown in the inset) as a function of the non-linear interaction U g . Clearly, there is a wide parameter setting, the plateau U g /Ω 0 = (−2, 1.4), where the BEC transfer process takes places with a high fidelity. In the parameter region U g /Ω 0 = (−2, 0), the system starts at w(t in ) = −1 following an energy eigenstate of the system that does not involve any bifurcation point. For U g /Ω 0 = (0, 1.4) the system crosses first a PB42 bifurcation to later on reaching a PB24 bifurcation point that, for these parameters, yields eventually w(t out ) = 1. In regions U g /Ω 0 = (−4, −2) and U g /Ω 0 = (1.4, 4) the  non-linearity yields θ(t out ) = 0 forcing the initial and final state to be the same.

C. Adiabatic Self Trapping
It is also possible to completely inhibit the BEC transport between the two wells by means of the matter wave analogue to the double STIRAP process (see inset of Fig. 7). Starting with the BEC located at the left trap, w(t in ) = −1, we show in Fig. 7 the population difference at the time when the tunneling rate is maximum, w(t Ω ), and at the end of the double-STIRAP sequence, w(t out ), as a function of the energy bias ǫ g . The adiabatic self-trapping process succeeds for a wide domain of parameter values, except for the region ǫ g /Ω 0 = (0, 0.8). In this last region, the solution to be followed reaches during its dynamics a saddle-point bifurcation while, for the rest of values of the energy bias, the selected solution does not involve any bifurcation point.

V. CONCLUSIONS
By means of a dark variable that results from the combination of the population imbalance and the spatial atomic coherence, we have investigated in detail the adiabatic dynamics of a BEC in a double well potential in the framework of the two-level approximation. We have shown that it is possible to robustly split, transport or trap a BEC by appropriate temporal variation of either the energy bias or the non-linear interaction together with the tunneling rate. All these proposals have been studied from a non-linear dynamics perspective deriving the stationary solutions of the system, evaluating their linear stability, and discussing the bifurcation scenarios. For the eventual implementation of the techniques here discussed, we want to highlight the recent work by M. Bauer et al. [27] reporting an accurate control on magnetic Feshbach resonances by means of laser light.
In order to check the validity of the two-mode approximation, we have numerically integrated the onedimensional Gross-Pitaevskii (GP) equation for the nonlinear adiabatic splitting of a BEC in a double-well potential, obtaining the characteristic 'plateau' of the STI-RAP protocol shown in Fig. 2(b). However, a detailed numerical investigation of the GP equation still remains to be performed to validate that the matter wave nonlinear STIRAP techniques here derived in the two-mode approximation could be used for matter wave interferometry or coherent transport of a BEC. In performing this analysis, the global landscape provided by the results of the two-mode approximation should be a useful roadmap.
Although it is out of the scope of the present paper, it would be very interesting also to extend the present nonlinear matter wave STIRAP techniques to the secondquantization formalism [35]. Within this formalism, we could elucidate whether the adiabatic splitting of the BEC results in a macroscopic superposition where the entire condensate localizes in one of the wells, i.e. in a NOON state [36], or in a perfect 50% spatial splitting of all the atoms. Note that in both cases the final population difference reads w(t out )=0.