Two component bosonic Josephson junctions in elongated traps

We study binary mixtures of Bose–Einstein condensates confined in a two-well potential within the framework of the Gross–Pitaevskii equation. We re-examine both the single component and the binary mixture cases for such a potential. We investigate the most usual dimensional reductions used to solve the Gross–Pitaevskii equations, including the one proposed by Reatto and collaborators. To this end, we compare numerical simulations of the 1D reductions with the full 3D numerical solutions of the Gross–Pitaevskii equation. Our analysis considers an experimentally feasible binary mixture of an F = 1 spinor condensate where two of its Zeeman manifolds are populated.


Introduction
The first evidence of phase coherence of a BEC was obtained in early experiments [1] where clean interference patterns appeared in the overlapping region of two expanding condensates. It has been only recently that clear evidence of an external bosonic Josephson junction in a weakly linked scalar (single component) BEC has been experimentally reported, first by the group of M. Oberthaler in Heidelberg [2], followed by the group of J. Steinhauer [3]. Internal Josephson dynamics has also been experimentally achieved [4].
In the experiment of [2], two condensates are confined in a double-well potential with an initial population imbalance between both sides which triggers the Josephson oscillations. The tunnelling of particles leads to a coupled dynamical evolution of the two conjugate variables, the phase difference between the two weakly linked condensates and their population imbalance. In spite of the system being dilute, the atom-atom interaction plays a crucial role in the Josephson dynamics, leading to new regimes beyond the standard Josephson effect, e.g. macroscopic quantum self-trapping (MQST).
The Gross-Pitaevskii (GP) theory provides an appropriate framework for investigating Josephson dynamics in weakly interacting systems provided: (a) the number of atoms is large enough so that quantum fluctuations can be neglected, and (b) the initial many-body state is of mean-field type. Josephson oscillations in scalar Bose-Einstein condensates were initially studied in [5,6], and since then they have been studied using different techniques, for a review see [7]. For N9100, the exact Josephson dynamics has been recently shown to depart from the mean-field Gross-Pitaevskii, as seen in the evolution of the population imbalance or with the appearance of fragmentation [8]. For small N the transition from Josephson to self-trapped dynamics has also been shown to involve the appearance of strongly correlated quantum states [9]. For a larger number of atoms, N01000, the full three-dimensional time-dependent Gross-Pitaevskii equation (GP3D) provides a reasonable agreement with the experimental data of [2], also see [7,10], where N $ 1150 and the system is allowed to evolve for less than a Rabi time.
However, since 3D dynamics need in general rather involved calculations, one can benefit from the fact that the barrier is created along one direction and the tunnelling of particles is mainly one dimensional (1D) in order to investigate the Josephson dynamics by means of effective 1D GP-like equations. Among these reduced GP equations, the non-polynomial nonlinear Schro¨dinger equation (NPSE) proposed by Prof. Reatto and collaborators in [11] has provided the best agreement with the experimental results in scalar condensates, whereas another effective 1D Gross-Pitaevskii equation (GP1D) fails to describe the dynamics for a large number of trapped atoms under the same trapping conditions as in the Heidelberg experiment [7].
Josephson oscillations in binary mixtures confined in double-well potentials have been addressed in a number of recent articles. The case of two-component BECs with density-density interactions has been studied within two-mode approaches in [12][13][14][15][16][17][18][19]. In [17] Mazzarella et al. go one step further and also consider GP1D simulations. Spin-dependent interactions have been addressed in [20][21][22]. In [20,23] the interest of studying Josephson dynamics in binary mixtures has been emphasized as it can give access to information on the different scattering lengths present in the system. The paper is organized as follows. The coupled Gross-Pitaevskii equations for a binary mixture are presented in Section 2, as well as the different onedimensional reductions of the GP3D equations for the mixture: GP1D and NPSE. In Section 3 we review results for the single component case. And in Section 4, we discuss the dynamics of a bosonic binary mixture in a double-well potential with the same parameters as in the experiment [2]. Finally, the conclusions are presented in Section 5.

Mean field approach: Gross-Pitaevskii equations
We consider a binary mixture of weakly interacting bosons at zero temperature, confined by a symmetric double-well potential, VðrÞ. For dilute systems with large enough number of particles, the Gross-Pitaevskii equation provides a suitable framework to study the dynamics. In the mean field approximation, each condensate is described by the corresponding wave function C i ðr; tÞ, with i ¼ a, b denoting each of the two components of the binary mixture. In most situations, the system will behave as if there were four weakly linked Bose-Einstein condensates, two per each component of the binary mixture per each side of the potential barrier. The mean field description will reflect this feature by the homogeneous quantum phase of C i ðr; tÞ at each side of the potential barrier.
The dynamical evolution of the two wave functions can be obtained by solving the two coupled GP equations: i h @C a ðr; tÞ @t ¼ For each component, the condensate wave function C i ðr; tÞ is normalized to 1, m i is the atomic mass, and g ii ¼ 4p h 2 a i =m i is the effective atomic interaction between atoms of the same species, with a i the corresponding s-wave scattering length. The coupling between both components is governed by the interspecies interaction g ab g ba , which depends on the specific nature of the binary mixture. The total number of atoms in the mixture is N ¼ N a þ N b . We will consider binary mixtures made of F ¼ 1 87 Rb atoms populating the m ¼ AE1 Zeeman sublevels. This implementation simplifies the dynamics as the inter-and intra-species couplings are similar in magnitude. Of course this choice limits the phenomena which can be observed. On the other hand its simplicity allows us to discuss the different dimensional reductions of the 3D equations taken in the literature.
We consider the same setup and the same trap parameters as in the experiments of the Heidelberg group [2]. There, a condensate of 87 Rb with 1150 atoms is confined to a fairly small region of $ 5 mm through the potential, In the full GP3D simulations we define the number of atoms in the left well as: The number of atoms in the right well is computed as N R ðtÞ ¼ N À N L ðtÞ. From these values, the population imbalance reads zðtÞ ¼ ðN L ðtÞ À N R ðtÞÞ=N. Analogous definitions are used in the GP1D and NPSE equations.
The phase at each side of the potential barrier is obtained in the following way. The phase at each point at a certain time, ðx, y, z; tÞ, is: Cðx, y, z; tÞ ¼ ½ðx, y, z; tÞ 1=2 exp½iðx, y, z; tÞ, where the local density is ðx, y, z; tÞ ¼ jCðx, y, z; tÞj 2 . Averaged densities, integrating over the z component, are defined as ðx, y; tÞ ¼ R 1 À1 dz ðx, y, z; tÞ. To visualize the phase coherence along some of the planes we define an average phase, integrating the z component, as ðx, y; tÞ ¼ ½1=ðx, y; tÞ R 1 À1 dz Â ðx, y, z; tÞ ðx, y, z; tÞ. The phase on the left, L ðtÞ, is defined as dzðx,y,z; tÞðx,y,z; tÞ: The phase on the right, R ðtÞ, is defined accordingly. The phase difference between each side of the barrier is computed as ðtÞ ¼ R ðtÞ À L ðtÞ.
Assuming that most of the dynamics takes place in the direction which contains the barrier, the x direction in our case, one can approximate the wave function of the system by Cðx, y, z; tÞ $ C 1D ðx; tÞ ' g:s: ð yÞ ' g:s: ðzÞ, where ' g:s: are the corresponding ground state wave functions for the trapping potential in the y or z direction in the absence of interactions. In this way it can be shown [24] that in the case of the one-component system, C 1D ðx; tÞ fulfills a Gross-Pitaevskii-like 1D equation, where the corresponding 1D coupling constant is obtained rescaling the 3D one, The extension to binary mixtures (and also to spinor condensates [25]) may be written down readily, where the rescaled couplings are g ij;1D ¼ g ij =ð2pa 2 ? Þ.

Non-polynomial Schro¨dinger equation (NPSE)
A more sophisticated reduction that includes to some extent the transverse motion of the elongated BEC in the corresponding potential is the so-called nonpolynomial Schro¨dinger equation, proposed for a scalar BEC by Prof. Reatto and collaborators [11].
The NPSE recovers the previously discussed 1D reduction in the weakly interacting limit, but it has been shown to provide the best agreement with the experimental results on Josephson oscillations between two coupled BECs [22]. The NPSE for the onecomponent case reads i h @Cðx; tÞ @t The generalization of the NPSE for two components in a binary mixture of BECs has been addressed in [26].
The system of equations, which become rather involved, can be further simplified in the case of equal interactions, both intra-and inter-species: i h @C j ðx; tÞ @t , and, as before, g 1D ¼ g=ð2pa 2 ? Þ, g g aa ¼ g bb ¼ g ab ¼ g ba , and R dxjC j ðxÞj 2 ¼ 1.

Numerical results for a single component condensate
First, let us consider the case of a single component condensate. In Figure 1 we present the time evolution of the population imbalance for Josephson, and selftrapped regimes. We compare the full GP3D (solid red) with the two 1D reductions, GP1D (dotted black) and NPSE (dashed blue). The first panel, (a), contains simulations performed with zero initial phase difference, i.e. Josephson oscillations and self-trapping cases. For the Josephson cases, zð0Þ ¼ 0:1, 0:35, the imbalance oscillates with a frequency which is mostly independent of the initial imbalance (for small imbalances). With zð0Þ ¼ 0:1 the oscillations are almost sinusoidal, while as we increase the initial imbalance their shape becomes more involved but remaining periodic. In the self-trapped case, zð0Þ ¼ 0:6, the atoms remain mostly on their initial side of the trap and there are short and small periodic oscillations as predicted by two-mode models [5]. At longer times, the imbalance is seen to decrease smoothly, implying a departure from the predicted two-mode dynamics.
The two 1D reductions give in most situations qualitatively similar results as the GP3D, but not quantitatively in all cases. The NPSE is seen to reproduce very well the GP3D in all the runs up to times near $40 ms. Above those times, the period of oscillation predicted by the NPSE is slightly shorter than the GP3D one. The GP1D on the contrary only captures the amplitude of oscillation in the Josephson dynamics, failing in all cases to give the same period as the GP3D or the NPSE. Moreover, the GP1D departs notably from the GP3D for the self-trapped case.
Panel (b) is computed very close to the critical value for the appearance of self-trapping in the full GP3D, zð0Þ ¼ 0:39 for ð0Þ ¼ 0. The GP1D and NPSE predict a critical initial imbalance close to the GP3D value.
Panel (c) contains two self-trapped cases obtained with an initial ð0Þ ¼ p and zð0Þ ¼ 0:2, and 0.4. Notice that for ð0Þ ¼ p the critical imbalance is smaller than for ð0Þ ¼ 0. The discussion is similar to the Josephson case, i.e. the NPSE captures most of the dynamical features of the GP3D while the GP1D only provides a qualitative understanding of the problem.
To further explore the quality of the 1D reductions, we present in Figure 2 the density profiles in the x direction after integrating the y and z components, ðx; tÞ ¼ R 1 À1 dy R 1 À1 dz jCðx, y, z; tÞj 2 at t ¼ 50 ms. The agreement between the NPSE and the GP3D is very good in most situations, except for the critical case, as expected. In all cases the density profiles show a clear bi-modal structure. The GP1D does not predict the correct density profiles and, as seen in the self-trapped case,  contribution of higher modes. The critical initial imbalance starting with no phase difference that we find numerically by means of the GP3D is the same as found in [10], z c ¼ 0:39, and differs from the one reported in [2], z c ¼ 0:5. The agreement of the NPSE with GP3D results justifies the use of NPSE in [2] to analyse their experiment.
4. Numerical solutions of the 3D Gross-Pitaevskii equations: binary mixture As discussed in Section 2, one feasible way of experimentally preparing binary mixtures of BECs is to consider a number of atoms populating the m ¼ AE1 Zeeman components of an 87 Rb F ¼ 1 spinor. The experimental observation of Josephson tunnelling phenomena by the Heidelberg group seems to be possibly extended to trap both Zeeman components. In this case the two components of the mixture have the same mass, M m a ¼ m b , and equal intra-species interactions, g aa ¼ g bb g. With respect to the interspecies interaction we will consider the case of 87 Rb which implies g ab $ g.
The mean field GP3D system of equations governing the dynamics of the three components of an F ¼ 1 spinor BEC can be written as [27] i h with H s ¼ À h 2 =ð2MÞ J 2 þ V þ c 0 n being the spinindependent part of the Hamiltonian. The density of the mth component is given by n m ðrÞ ¼ j m ðrÞj 2 , while nðrÞ ¼ P m j m ðrÞj 2 is the total density normalized to the total number of atoms N. The couplings are c 0 ¼ 4p h 2 ða 0 þ 2a 2 Þ=ð3MÞ and c 2 ¼ 4p h 2 ða 2 À a 0 Þ=ð3MÞ, where a 0 and a 2 are the scattering lengths describing binary elastic collisions in the channels of total spin 0 and 2, respectively. Their values for 87 Rb are a 0 ¼ 101:8a B and a 2 ¼ 100:4a B [28]. Since the spindependent coupling, c 2 , is much smaller than the spinindependent one, c 0 , and the total number of atoms that we will consider is relatively small, N ¼ 1150, the population transfer between the different components can be neglected [20]. Therefore, in our calculation the number of atoms in each sublevel remains constant in time allowing us to treat the system as a binary mixture of components a and b. Comparing the system of Equations (2) and (8) the value of the couplings can be read off, g aa ¼ g bb ¼ c 0 þ c 2 and g ab ¼ g ba ¼ c 0 À c 2 .
Once the total number of atoms is fixed, we investigate the Josephson-like dynamics for different numbers of atoms populating each component, N a ¼ f a N and N b ¼ f b N, and for different initial conditions z a (0), z b (0), a ð0Þ and b ð0Þ.

Phase coherence and localization
The numerical solutions of the GP3D for the single component system showed two features. First, the atoms remained mostly localized in the two minima of the potential well and secondly, each group of atoms had to a large extent the same quantum phase. This, clearly supported the picture of having two BECs, one at each side of the barrier, with a well-defined phase at each side during the dynamical evolution.
As in the scalar case, our exact GP3D numerical solutions of the dynamics of the binary mixture in several initial conditions of population imbalances and phase differences show the same two distinctive features, see Figure 3. First, the density of atoms for each component is always bi-modal, with the two atom bunches centred around the minima of the potential well. Secondly, the phase of the wave function is mostly constant for each species at each side of the potential trap. Thus, we find that the GP3D does predict the dynamics to be mostly bi-modal also for the binary mixture case.
At the end of the section we will consider some deviations from the bi-modal behaviour that are found in very specific conditions.

Rabi and anti-Josephson modes
The two predictions of the two-mode model described in [20,22], namely the 'anti-Josephson' behaviour, and the enhancement of the Rabi mode, are confirmed by the NPSE and GP1D simulations as can be seen in Figure 4. In Figure 4 (left panels) we consider a very polarized case with f a ¼ 0:8. As expected from twomode analysis the dynamics of the most populated component should to a large extent decouple from the less populated one and perform fast Josephson oscillations with a frequency close to the corresponding one for the scalar case, ! J ¼ ! R ½1 þ NU=ð h! R Þ 1=2 (where ! R is the Rabi frequency, that is, the frequency of oscillation of a single atom in the double-well potential). The GP3D simulation is seen to confirm the above predictions. The less abundant component is strongly driven by the most populated one and shows an anti-Josephson behaviour as described in [20].
Another prediction is related to the behaviour of z a þ z b and z a À z b in the non-polarized case, f a ¼ f b .
As explained in [22] in this case the difference, z a À z b , should enhance the long mode which oscillates with the Rabi frequency of the system, while the sum z a þ z b should mostly oscillate with the Josephson frequency, ! J . In the right part of Figure 4 we present the extreme case when z a ð0Þ ¼ Àz b ð0Þ computed with GP3D and NPSE. In this case, both population imbalances and phase differences oscillate mostly with the Rabi frequency of the system, keeping during the time evolution z a þ z b $ 0.
Both 1D reductions produce qualitatively similar physics. The only important difference is that the frequency of the Josephson oscillations is higher in the GP1D, as occurred already for the single component [22].
Interestingly, they predict different Josephson oscillations while the Rabi frequencies are similar. The long oscillation corresponding to the Rabi mode is seen to agree well with the corresponding long oscillation seen in the right panels of Figure 4. The Josephson-like oscillations of binary mixtures of spinor F ¼ 1 87 Rb BECs around the ðz 0 a , 0 a , z 0 b , 0 b Þ ¼ ð0, 0, 0, 0Þ are therefore essentially controlled by two frequencies, ! R and ! J .

Small oscillations around z 0
a,b , d/ 0 a^0 and d/ 0 b^n ; effects beyond two-mode As explained in detail in [22], the dynamics for these conditions depend on the specific value of f b considered. The trivial equilibrium point at z i ¼ 0 exists  provided f b 90:43 [22]. This prediction of the twomode models is observed both in the GP3D and NPSE as can be seen in Figure 5. In the figure, we consider a simulation with z a ð0Þ ¼ 0:1, z b ð0Þ ¼ À0: 15 Again, the comparison between the NPSE and GP3D is very satisfactory. The NPSE captures almost completely the dynamics up to times of 100 ms. In all cases, the NPSE reproduces correctly both the phase difference and population imbalance. The only sizeable discrepancies occur for times 080 ms in the run without equilibrium point (right panel).
An example of simulations around non-trivial equilibrium points is presented in Figure 6. These points involve very large and opposite initial population imbalances for both components. In Figure 6 we consider a case with initial conditions very close to the predicted equilibrium point using the standard twomode, z a ð0Þ ¼ À0:78, and z b ð0Þ ¼ 0:99, with f a ¼ 0:1. Also in the same figure we consider a similar run but with f a ¼ 0:9. In both cases the NPSE and GP3D predict a very similar dynamics.
Most of the dynamics described up to now can to a large extent be understood within two-mode models [22]. There are, however, a number of situations where the two-mode fails. Some are a direct consequence of having two components evolving in the same doublewell potential, others are due to having initial configurations, mostly with large initial imbalances, producing situations where the atom-atom interaction energy per atom is comparable to the gap between the first excited state and the second/third excited states.
We can distinguish two different cases: (a) involving excitations along the coordinate which contains the barrier, (b) involving excitations of the transversal coordinates.
An example of (a) is shown in Figure 6. There, as clearly seen in the density profiles along the x direction, the two-mode approximation is no longer valid. The simplest way of seeing this is by noting the zero in the density of one of the components at x $ 2 mm. This effect beyond two-mode is well taken care of by the NPSE which reproduces the density profile quite well during most of the time evolution considered in the simulation. Thus, the excitations of higher modes along the direction which has not been integrated out in the 1D reduction do not pose a great difficulty to the 1D reductions.
The second type of effect beyond two-mode, (b), involves excitations of the transverse components. These effects are present in any binary mixture calculation whenever the intra-and inter-species interactions are not equal, see for instance the last section of [22].

Conclusions
The rich dynamical regimes which take place in binary mixtures, like double self-trapped modes, Josephson oscillations, or zero-and -bound phase modes, have been studied by performing full GP3D simulations covering all the relevant initial conditions. The 3D numerical solutions of the Gross-Pitaevskii equations have been used to critically discuss the validity of the most common 1D reductions of the GP equations, GP1D and NPSE.
To fix the conditions of the dynamics, we have focused on one particular setup which corresponds to a natural extension of the experiments reported in [2]: the case of a binary mixture made by populating two of the Zeeman states of an F ¼ 1 87 Rb condensate. As discussed in the present paper, this setup already allows us to observe and characterize a large variety of phenomena which are genuine for the binary mixture, e.g. anti-Josephson oscillations in highly polarized cases, long Rabi-like oscillating modes, zero-and -locked modes, etc.
Two commonly employed dimensional reductions of the GP3D, the GP1D and NPSE, have been shown to differ substantially among each other, with the NPSE being clearly in much better agreement with the original 3D dynamics for a broader set of conditions. In general, the GP1D describes essentially the correct physics but quantitatively far from the GP3D predictions. Also, for self-trapped cases already in the single component case, it departs from the two-mode behaviour earlier than the GP3D or the NPSE. The agreement between the NPSE and the full 3D dynamics is astonishingly good both for single component and the considered binary mixtures, where the intraand inter-species are very similar and the NPSE equations are particularly easy to handle. This agreement is not only seen for fully integrated magnitudes, for instance population imbalances, but also for the density profiles predicted along the direction hosting the barrier.