{'fraction_non_alphanumeric': 0.04345505310412789, 'fraction_numerical': 0.016273404138837183, 'mean_word_length': 4.21074026529739, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weakly linked binary mixtures of F=1 87-Rb Bose-Einstein condensates

We present a study of binary mixtures of Bose-Einstein condensates confined in a double-well potential within the framework of the mean field Gross-Pitaevskii equation. We reexamine both the single component and the binary mixture cases for such a potential, and we investigate in which situations a simpler two-mode approach leads to an accurate description of their dynamics. We also estimate the validity of the most usual dimensionality reductions used to solve the Gross-Pitaevskii equations. To this end, we compare both the semi-analytical two-mode approaches and the numerical simulations of the 1D reductions with the full 3D numerical solutions of the Gross-Pitaevskii equation. Our analysis provides a guide to clarify the validity of several simplified models that describe mean field non-linear dynamics, using an experimentally feasible binary mixture of an F=1 spinor condensate with two of its Zeeman manifolds populated, m=+1 and m=-1.

The phase coherence of a Bose-Einstein Condensate (BEC) is an important and characteristic property of ultracold bosonic gases that leads to fascinating macroscopic phenomena such as interference effects or Josephson-type oscillations. Two condensates trapped in a double-well potential exhibit interference fringes when the barrier is released and the two expanding condensates, with a well-defined quantum phase, overlap. Instead, if the barrier is not switched-off and is large enough to ensure a weak link between both condensates in each side of the trap, the quantum phase difference will drive Josephson-like effects, which consist on fast oscillating tunneling, much faster than the single particle tunneling, of atoms through the potential barrier [1,2]. The first evidence of the phase coherence of a BEC was obtained in early interference experiments [3] where clean interference patterns appeared in the overlapping region of two expanding condensates. It has been only recently that a clear evidence of a bosonic Josephson junction in a weakly linked scalar BEC ‡ has been experimentally reported by the group of M. Oberthaler in Heidelberg [4]. In this experiment, two condensates are confined in a double-well potential with an initial population imbalance between both sides which triggers the Josephson oscillations. The tunneling 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 very dilute, the inter-species 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) mean field theory provides a natural framework for investigating Josephson dynamics in weakly interacting systems at very low temperature. Josephson oscillations in scalar Bose-Einstein condensates have been theoretically studied by using different techniques [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. As expected, the full three dimensional time-dependent Gross-Pitaevskii equation (GP3D ) provides an excellent agreement with the experimental data [9,10,20]. 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 tunneling of particles is mainly one dimensional (1D) to investigate the Josephson dynamics by means of effective 1D GP-like equations. Among these reduced GP equations, the non-polynomial nonlinear Schrödinger equation (NPSE ) proposed in Ref. [21] has provided the best agreement with the experimental results in scalar condensates [20], whereas another effective 1D Gross-Pitaevskii equation ( GP1D ) fail to describe the dynamics for large number of trapped atoms in the same trapping conditions as in the Heidelberg's experiment [10,20].
Interactions are important to understand the different regimes of the tunneling dynamics. Therefore, multi-component BECs in double-well potentials offer an interesting extension to study phenomena related to phase coherence. In particular the Josephson dynamics will become richer due to the interplay between intra-and inter-species interactions.
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 Refs. [22,23,24,25,26,27,28,29]. Refs. [26,29] go one step further and also consider GP1D simulations. Spin-dependent interactions have been addressed in Refs. [30,31]. Josephson dynamics in spinor condensates confined in double-wells, characterized by an exchange of population between different Zeeman components, has also been investigated in Refs. [32,33]. In Refs. [26,30] the interest of studying Josephson dynamics in binary mixtures has been emphasized as it can give access to information of the different scattering lengths present in the system.
Recently, the equations of the tunneling dynamics in a binary mixture within the two-mode approximation to the GP equations have been derived in Ref. [25]. However, the authors have not compared their two-mode analysis to direct numerical resolutions of the GP equation, and have also not provided microscopic values to the parameters of the two-mode equations. Their main result is the description of a symmetry breaking pattern occurring when the inter-and intra-species interactions differ substantially. In Ref. [26] a comparison of the standard two-mode approach and the coupled GP1D equations for the mixture has been presented for one specific double-well potential which allows an analytical treatment.
For single component BECs, it has been already studied the range of validity of the different approximations to the Josephson dynamics by comparing with GP3D calculations and with the experimental results. However, no comparison with the full GP3D dynamics has been yet performed for a binary mixture in a double-well potential.
The aim of this paper is to investigate systematically the tunneling dynamics of a binary mixture of BECs trapped in a double-well potential, as well as the validity of the different mean-field approximations. We consider a mixture of two components obtained by populating two Zeeman states of an F = 1 87 Rb condensate confined in the same double-well potential as in the experiments [4]. This system corresponds to a natural extension of the experimental work of Ref. [4], where only one of the Zeeman components was populated.
We provide a general overview of the different techniques used to investigate Josephson dynamics within the two-mode model (standard and improved two-mode) and within the Gross-Pitaevskii framework (one-dimensional reductions of the GP equation, GP1D and NPSE ). To this end, we solve the full 3D time-dependent GP equation for the mixture as a reference to assess and analyze the validity of the previous approximations.
The paper is organized as follows. The general framework of the coupled Gross-Pitaevskii equations for a binary mixture is presented in Sec. 2. In Sec. 3, we derive analytic two-mode models both for single and two component systems. First we recall the standard two-mode model (S2M ). Then we derive the equations of the improved two-mode model (I2M ) for a binary mixture, generalizing the work for a single component BEC performed in Ref. [9]. We discuss the stability of the dynamical equations and look for the stationary points for a binary mixture. In Sec. 4, we analyze the different one-dimensional reductions of the GP3D equations for the mixture: GP1D and NPSE . In Sec. 5, we revisit the dynamics of a single component condensate in a double-well potential with the same parameters as in the experiment [4]. The tunneling dynamics in two-component systems is accurately discussed in Sec. 6. We obtain the dynamics by solving the coupled GP3D equations for the mixture and show that for certain conditions there exists a good agreement between I2M and GP3D , as well as for NPSE and GP3D . The range of validity of the two-mode models is explored, paying special attention to situations which fall beyond the two-mode approximation. Finally, we discuss cases that present characteristic features arising from the mixture, with no analog in the tunneling dynamics of a single component BEC. Conclusions are given in Sec. 7.

Mean field approach: Gross-Pitaevskii equations
We consider a binary mixture of weakly interacting atoms at zero temperature, confined by the same double-well potential, V (r). For dilute systems with sufficiently large 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 Ψ i (r; t), with i = a, b denoting each of the two components of the binary mixture. To avoid any misunderstanding let us remind the reader that we are describing two different kind of atoms, a and b, which evolve on a double-well external potential. 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 Ψ i (r; t) at each side of the potential barrier, as will be discussed in great detail in the following sections.
The dynamical evolution of the two wave functions can be obtained by solving the two coupled GP equations: For each component, the condensate wave function Ψ i (r; t) is normalized to 1, m i is the atomic mass, and g ii = 4π 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 inter-species 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 . There are many experimental possibilities to study the dynamics of binary mixtures of BECs. We will restrict our study to one of them, which is experimentally feasible. We will consider binary mixtures made of F = 1 87 Rb atoms populating the m = ±1 Zeeman sublevels [34]. This implementation greatly simplifies the dynamics as the interand intra-species couplings are very similar in magnitude. Of course this choice limits the phenomena which can be observed, e.g. the interesting symmetry breaking pattern discussed in Ref. [25], which relies on the inter-species coupling being larger than the intra-species one, will not take place, see Sec. 6.5.
On the other hand its simplicity allows to discuss in detail the different approaches taken in the literature, e.g. two-mode models of the GP equations, one dimensional reductions, etc. As occurred in the scalar case, the dynamical features contained in Eqs. (1) can to a large extent be described by a simplified two-mode model for each component. In the next section we follow Refs. [5], [9] and [25] and derive two-mode expressions for the scalar and binary case. The usual assumption of neglecting the overlaps involving the right and left modes gives rise to the so called standard two-mode (S2M ) equations, while retaining them one also gets a closed system of equations, the improved two-mode (I2M ). Both the S2M and I2M are also derived for the binary mixture case.

Two-mode approaches
The two-mode approximation allows to study the dynamics of weakly linked Bose-Einstein condensates, without solving the full GP3D nor reducing the dimensionality of the GP equation [5,9]. Depending on the specific double-well potential, e.g. on the energy gap between the first two levels of the single particle Hamiltonian and the next two, it can provide an excellent description of the full GP solution. The relevant physical quantity is the ratio between the energy gap between the ground state and first excited state of the double-well potential, δE 0;1 = E 1 − E 0 § and the energy difference of the ground state and the second excited state, δE 0;2 = E 2 − E 0 . The smaller the ratio δE 0;1 /δE 0;2 the more accurate the two-mode approach. The two mode description characterizes the dynamics of the scalar condensate in a double-well potential with only two variables: the relative population and the phase difference between the left and right side of the potential barrier. § Note that this is zero if the barrier is infinitely high.

Standard two-mode model for the single component case
The GP equation for the scalar case corresponds to a particular limit of the GP equations for the binary mixture, Eqs. (1), We will make use of the following notation, H 0 = −( 2 /2m)∇ 2 + V (r), and H g [Ψ] = gN|Ψ(r; t)| 2 . Let us recall the two-mode approximation for a single component condensate in a double-well potential. We consider N interacting atoms with atomic mass m, and coupling constant g, trapped in a symmetric double-well potential V (r). When both sides of the potential barrier are weakly linked, the total wave function can be approximately written as a superposition of two time-independent spatial wave functions Φ L(R) (r) mostly localized at the left (right) side of the trap: The left and right modes, can be expressed as linear combinations of the ground (+) and the first excited (−) states of the double-well potential including the interaction term. They satisfy, (H 0 + H g [Φ ± ])Φ ± = µ ± Φ ± , and the left/right modes can be written as [9]: We observe that in a symmetric double-well, Φ ± have a well defined parity: Φ ± (r) = ±Φ ± (−r) , and therefore Φ i Φ j = δ ij with i, j = +, −. Since they are stationary solutions of the GP equation, Φ ± are real functions, and so are the left and right modes Φ L(R) . The integrated density in the z-direction, ρ(x, y) = dz|Ψ(x, y, z)| 2 , associated to the ground (Φ + ), and first excited (Φ − ) states are depicted in Fig. 1 together with the densities associated to the left and right modes. The plots correspond to the experimental set up of Ref. [4]. From the phase coherence properties of a BEC, one can assume that the wave function in each side of the trap has a well defined quantum phase φ j (t), which is independent of the position but changes during the time evolution. We can write, where N L(R) (t) corresponds to the number of atoms on the left (right) side of the trap, and the total number of atoms is As a first step, we consider the so-called standard two-mode approximation (S2M ), which neglects a certain set of overlapping integrals involving mixed products of Φ L and Φ R . This approximation yields essentially the correct qualitative results in the scalar condensate although it may lead to incorrect quantitative predictions depending on the specific barrier properties [4,9] .
Inserting the two-mode ansatz (3) in the GP equation for a single component condensate (2) and neglecting terms involving mixed products of Φ L and Φ R of order larger than one, yields into a system of equations for the two localized modes which can be written in terms of two dynamical variables: the population imbalance z(t) = [N L (t) − N R (t)]/N and the phase difference δφ(t) = φ R (t) − φ L (t) between each side of the barrier: where ω R = 2K/ is the Rabi frequency and For a symmetric double-well, E 0 L = E 0 R and U L = U R ≡ U, therefore ∆E = 0. Moreover, the Rabi frequency only appears as a scale in the problem and thus can be absorbed in the time by rescaling t → ω R t. Then, together with the definition Λ ≡ NU/( ω R ), we obtain,ż cos δφ(t) .
Note that, Λ > 0 and Λ < 0 correspond to repulsive and attractive atom-atom interactions, respectively. There are different regimes depending on the initial values of the population imbalance and phase difference, z(0) and δφ(0), Sec. 3.3. From the energy functional of the GP equation (2): and using the two-mode ansatz (3), we can define the conserved energy per particle of the system as: where C is a rescaling constant. If we consider again a symmetric double-well we have: Note that the equations of motion (8) can be written in the Hamiltonian form: being z and δφ canonical conjugate variables.

Improved two-mode model for the single component case
Ananikian and Bergeman [9] noticed that for a symmetric double-well there was no need to neglect any of the overlapping integrals to obtain a closed set of equations relating z and δφ. Thus, remaining in the two-mode approximation but retaining all the overlaps it is straightforward to write down the following set of equations (cf. Eqs. (22) in Ref. [9]), called the "improved two-mode" (I2M ) equations, Defining, γ ij = g dr Φ 2 i (r)Φ 2 j (r) , i, j = +, −, we have, A = N (10γ +− − γ ++ − γ −− ) /4, B = 2K + N(γ −− − γ ++ )/2, and C = gN dr Φ 2 L (r)Φ 2 R (r). As discussed in detail in Ref. [9], the physics arising from the I2M is similar to the one present in the S2M . The I2M , however, is in much better agreement with the GP3D for a broader set of double-well potentials, as we will see in Sec. 5. In particular, for the double-well considered in the experimental setup of the Heidelberg group [4]  the S2M (with the corresponding microscopic parameters, K and U, computed from the GP3D with the experimental 3D potential and experimental coupling, g) does not correctly predict the physics of the experiment, mainly due to the importance of the neglected overlaps, and not to a dynamics far from a two-mode one.

Stability analysis
In this section we use the S2M to analyze the stability of the single component system, and focus on the case of repulsive interactions Λ > 0. Using the Hamiltonian (11) and the equations of motion (12), the stationary points (z 0 , δφ 0 ) can be found by solving the equations: To asses the stability of these points, we need to study the Hessian matrix of the system, which for the possible values of the phase difference, δφ 0 = 0 or π, is always diagonal and its eigenvalues are ∂ 2 z H| z 0 ,δφ 0 and ∂ 2 δφ H| z 0 ,δφ 0 . Depending on the sign of these eigenvalues the stationary points will be maxima, saddle points or minima. The stationary points and their stability are summarized in Table 1.
The evolution of the system can be represented on a z − δφ plane, where the system follows trajectories with constant energy, H, see curves in Fig. 2. Note that oscillations around a stationary point, closed curves, occur only if the central point is either a maximum or a minimum of the energy, but not a saddle point. As we will see in the following sections, these orbits will give rise to the Josephson oscillations and to the Figure 3. Different regimes for a set of initial conditions, imbalance z(0) in the y-axis and phase difference δφ(0) in the x-axis. The upper panels correspond to repulsive interactions while the lower ones to attractive interactions. The values of |Λ| are 0.5, 1.5, and 5. for the left, middle and right panels respectively. Grey regions correspond to Josephson oscillations, blue regions to π-modes (upper panels) and zero-modes (lower panels), and red regions to running phase modes.

Symmetry between attractive and repulsive interactions
The stability analysis has been presented only for repulsive interactions, but from the system (8) we can see that if we change the interactions, Λ −→ −Λ, we recover the same system of equations if δφ −→ π − δφ: d dt which means that the dynamics of the system and the different regimes are the same for both interactions, with a phase-shift of π. This can be seen in Fig. 3, that shows (z 0 , δφ 0 ) stationary minimum saddle maximum (0, 0) Table 1. Stationary points of the system for repulsive interactions, Λ > 0, and their stability. the behavior of the system for a given set of initial conditions. The upper panels are for repulsive interactions Λ > 0 and the lower ones for attractive interactions Λ < 0. The grey regions correspond to Josephson oscillations, the blue regions to zero-and π-modes, and the red regions to running phase modes, as it will be seen in the following sections.

Josephson dynamics
This regime is characterized by fast oscillating tunneling of population across the potential barrier. Plotted in a z − δφ map, the system evolves following closed trajectories around a minimum or a maximum (z 0 = 0, δφ 0 ) configuration, with a zero time-average of the population imbalance, < z > t = 0. The stability analysis shows that for Λ > −1, which corresponds to repulsive or slightly attractive interactions, the stationary point (z 0 = 0, δφ 0 = 0) is a minimum permitting Josephson oscillations around it. Analogously, for Λ < 1, either attractive or slightly repulsive interactions, the stationary point (z 0 = 0, δφ 0 = π) becomes a maximum, and therefore also allows for closed orbits around it. For |Λ| > 1, there are Josephson oscillations around only one point: (z 0 = 0, δφ 0 = 0), or (z 0 = 0, δφ 0 = π). However, in the region of weak interaction, |Λ| < 1, the oscillations around both points are allowed.
It is also interesting to study the behavior of the system for small oscillations around these two stationary points of zero imbalance, smallest orbits in Fig. 4 (a). In this limit, the system (8) can be linearized giving the dynamical equation:z(t) = −z(t)(1 + Λ cos δφ 0 ) with cos δφ 0 = ±1. The population imbalance performs sinusoidal oscillations with a frequency ω J = ω R 1 + Λ cos δφ 0 , independent of the initial conditions. Note that this frequency only exists when these points are either maxima or minima. The phase difference oscillates with the same frequency but with a phase-shift of π/2 with respect to the imbalance. If the initial population imbalance increases, the dynamics of the system changes substantially to non-sinusoidal oscillations, and the frequency becomes dependent on the initial conditions.

Macroscopic quantum self trapping
In the case of repulsive interactions, we have seen that for Λ > 1, the stationary point (z 0 = 0, δφ 0 = π) becomes a saddle point and there appear two maxima, A similar behavior is found for attractive interactions. These stationary points allow for oscillations around them with < z > t = 0. In fact, in this regime, the imbalance has the same sign during the evolution, and therefore one of the wells is always overpopulated.
This regime is called macroscopic quantum self trapping, as the tunneling is strongly suppressed and the particles remain mostly trapped in one of the wells. This is a phenomenon arising from the atom-atom interaction, which appears as a non-linearity in the Gross-Pitaevskii equation.
The critical condition for the existence of the MQST regime can be found by imposing that the system remains on one side of the trap [5]. For a given set of initial conditions, z(0) = 0 and δφ(0), the system will remain trapped if, where the limits of the interaction parameter are due to the fact that only when |Λ| > 1 the (z 0 = 0, δφ 0 ) stationary points exist. In this regime however, there are two different kind of MQST depending on whether the phase difference evolves bounded, giving the so-called zero-and π-modes, or whether it evolves unbounded, increasing (or decreasing) always in time, giving rise to the running phase modes.
For values of the interaction parameter of 1 < |Λ| < 2 the only MQST regime that one can have is the zero-mode for attractive interactions and the π-mode for repulsive interactions (which are plotted in blue dotted lines in panel (b) of Fig. 4). In these regimes the phase difference evolves bounded around δφ = 0 and δφ = π, respectively.
On the other hand, for values of |Λ| > 2 one can have both classes of MQST. In general however, for a given set of initial conditions, the system will evolve following a running phase mode (dashed-red lines of panel (c) of Fig. 4), because the values of z 0 = ± 1 − 1/Λ 2 , that allow closed orbits, are very close to 1 (see the small π-modes of panel (c) in blue dotted lines).
In panel (c), one can see that the broadest closed orbit around (z 0 = 0, δφ 0 = π), for Λ > 2, is the one that goes through (z = ±1, δφ = 0). Notice that an orbit that crosses the δφ = 0 axis in any other point, z = ±1, would correspond to a running phase mode. The case of attractive interactions can be understood by taking into account the phase-shift of π in δφ. The latter can be used to find the condition to have bounded or running phase difference modes. For a given set of initial conditions (z(0), δφ(0)) fulfilling the self-trapping condition (16), the system will evolve in a bounded phase mode only if: Moreover, in a zero-or a π-mode MQST, we can study small oscillations around the corresponding minima or maxima, z(t) = z 0 + δz and δφ(t) = δφ 0 +δφ(t), so the linearized system (8) becomes: which gives a sinusoidal behavior with a frequency:

Standard two-mode model for the binary mixture
Let us recall the two-mode approximation for weakly linked binary mixtures [22,23,24,25]. The total wave function of each component is written as a superposition of two time-independent spatial wave functions localized in each well: with Φ iα |Φ jβ = δ ij δ αβ , i, j = a, b and α, β = L, R. For a given component, the condensates in each side of the trap are weakly linked. Then, as in the scalar case, one can assume that the wave function in each side of the trap has a well defined quantum phase φ j,α (t), which is independent of the position but which changes during the time evolution. Thus, N j,L(R) (t) corresponds to the population of the j-component on the left (right) side of the trap, with N j = N j,L (t) + N j,R (t). Inserting the two-mode ansatz (20) in the coupled GP equations for the mixture (1), retaining up to first order crossed terms yields the following system of coupled equations: cos δφ a (t) where, with i, j = a, b and α = L, R. Let us consider a mixture with the same atomic mass for both components M ≡ m a = m b , which are trapped in the same symmetric double-well potential. Then, the localized modes are the same for both components but depend on the site: Defining for each component the population imbalance and phase difference between both sides of the barrier, the above equations can be rewritten as: where ω R = 2K/ is the Rabi frequency, the same for both species. It is useful to define, These equations correspond to two coupled nonrigid pendulums. The stability of these systems of equations have been analyzed recently in Ref. [24].

Improved two-mode model for the binary mixture
As was noted for the scalar case in Ref. [9], it is not mandatory to neglect any of the overlaps to obtain a closed set of equations relating the population imbalances and phase differences for a symmetric double-well potential. The complete set of two-mode equations were called the improved two-mode (I2M ) equations. In principle, if the experimental setup is appropriately chosen the left and right modes may be quite well localized at each side of the trap. In this case, the S2M equations are expected to provide quantitative agreement with the experimental data. When the two-modes are not so well localized, then it becomes necessary to consider the I2M to have quantitative agreement. In [9] the authors considered explicitly the set up of the Heidelberg group and showed that the I2M is necessary in the single component case to provide a quantitative understanding of the experimental data.
Following similar steps as in the previous section and assuming the double-well potential to be symmetric as in the experiment, then the wave functions for the ground state and first excited state, Φ j± (r), have a well defined parity. The symmetry properties and the ortho-normalization conditions are capital to derive the coupled equations within the I2M model: provides an exact description of the dynamics in the symmetric doublewell potential, with no approximations beyond the assumption of a two-mode ansatz of the total wave function Ψ j (r; t), Eq. (20).
The resulting system of equations relating the population imbalance and phase difference for each component within the I2M approximation reads : where we have defined Our system of equations differs slightly with the previously derived ones, cf. appendix of Ref. [25]. We believe their system has some minor errors, which do not affect their discussion which is based on the S2M equations.
µ j + and µ j − are the chemical potentials of the ground and first excited state of the j component, that can be calculated from the time-independent GP equation for Φ j± , respectively. Analogously one can define 2K ba by exchanging the subindex a and b in the previous expression.
Notice that we have kept the full 3D dependence of the wave functions Φ j± (r), instead of averaging the transverse spatial dependence as in Refs. [9,25]. Thus, the coupling parameters g ij in Eqs. (29) are the 3D ones and are not renormalized.
The equations for the I2M are essentially similar to the S2M . The main difference is that the tunneling term, K ab (t), is time dependent and contains effects due to the interactions. As expected, if the localization of the modes is increased, i.e. by increasing the barrier height, K ab (t) approaches the constant value, 2K ab → µ a − − µ a + , which equals 2K of Eq. (23). The coupled equations obtained in the I2M model reduce to well-known dynamical equations in two limiting cases: i) Setting to zero the overlapping integrals that involve mixed products of left and right modes of order larger than 1, the I2M equations reduce to the S2M model for the mixture, Eqs. (25).
ii) Assuming a noninteracting mixture, the inter-species interaction is g ab = 0, and the I2M equations for the mixture reduce to a two non-coupled system of equations, that are the dynamical equations of the I2M for a single component, Sec. 3.2.
As discussed in the introduction we are interested in the particular case of a binary mixture made of atoms populating two different hyperfine states. Then, both components have the same mass M, and are trapped in the same symmetric double-well potential. We initially restrict to the case in which the inter-species interaction is also almost equal to the intra-species one, g ≡ g aa = g bb ∼ g ab . This is the situation for F = 1 m = ±1 of 87 Rb. This case allows straightforward comparisons between the results of the I2M and the ones obtained by solving the NPSE or GP1D for a mixture explained in Sec. 4.
The ground and first excited states in a symmetric double-well potential are the same for both components. Moreover, since g = g ab the overlap integrals (29) reduce to: and the chemical potentials µ a α = µ b α ≡ µ α with α = +, −. This yields the following relations: K ab = K ba and ∆ a = ∆ b . The I2M system reduces to: In this case both components obey the same system of coupled differential equations. Then, if the initial conditions are the same for both, z a (0) = z b (0) and δφ a (0) = δφ b (0), they will evolve with the same imbalance and phase, and no mixture effects will be observed.

Regimes for binary mixtures
We proceed now to analyze the stability of the system of equations (26), cf. see the appendix of Ref. [22]. As in the single component case, and in order to get analytical results that allow for a physical insight, we perform the study in the framework of the S2M approximation. First we note that an stationary point, defined by the equations: and the following system of equations, Therefore there are four different cases: noting that in all of them there is an obvious stationary point, z 0 a = z 0 b = 0. These stationary points will be referred to as "trivial stationary points". We need to find the conditions for non-trivial solutions in each case.
The stability of the system is analyzed by considering small variations around the stationary points for each of the four situations. Defining the displacements η i , the following system of equations for the η's can be derived from Eqs.  Table 3. Same as Table 2 but retaining up to the first order in β, whereΛ = Λ (1 + β). where, In Table 2 we give the explicit values of the eigenfrequencies of Ω for the trivial stationary points, z 0 i = 0. These are obtained for the case under consideration, where Λ > 0 and Λ > 0.
Approximate simpler expressions for the same eigenfrequencies can be derived for the case whenΛ ∼ Λ. DefiningΛ = Λ(1 + β) and retaining up to terms of order β, one obtains the frequencies listed in Table 3.
3.6.1. Stationary points with (δφ 0 a = 0, δφ 0 b = 0) In this case, the condition for the existence of non-trivial solutions to the equations (34) depends on the slope at the origin of the two curves (34) [22]. The condition guarantees the existence of two additional solutions besides the trivial one. If we restrict ourselves to the case ofΛ ∼ Λ > 0, we have that (38) cannot be fulfilled and therefore the only stationary point is the trivial one, z 0 a = z 0 b = 0.
In this case, it is straightforward to substitute the stationary point, z 0 a = z 0 b = 0 and δφ 0 a = δφ 0 b = 0 into Eq. (37) to get, which has two eigenvalues, listed in Table 2.
In the very polarized case, f a ∼ 1, f b ∼ 0, the population imbalance of the most populated component decouples from the less populated one and oscillates with the Josephson frequency w J = ω 1 . The less populated component is driven by the other component and follows its dynamics, thus giving rise to "anti-Josephson" oscillations. The smaller frequency oscillation seen in the population imbalance of the less populated component is ω 2 , which in this case is very similar to ω R [30].
Also interesting is the non-polarized case, f a = f b = 1/2, then (assumingΛ ∼ Λ, which is the case for 87 Rb), and defining ∆z( Therefore, ∆z behaves as the single component case, oscillating with the usual Josephson frequency, w J = ω R √ 1 + Λ while δz oscillates with the Rabi frequency, as would a single component case in the absence of atom-atom interactions. This mode can be further enhanced by imposing that z a (0) = −z b (0) thus forcing both imbalances to oscillate with the same frequency.
We have proposed in Ref. [30] to use these two configurations to extract the frequencies governing the dynamics of the system in order to obtain the microscopic atom-atom interaction. The idea was to profit from the fact that the difference between the inter-and intra-species interaction is small for the case of 87 Rb,Λ = Λ(1 + β), so we can use the expressions listed in Table 3, Note that in the anti-Josephson case the oscillation with larger period is ω 2 2 = ω 2 R (1 + O(βf b )) and the shorter is ω 2 1 = ω 2 R (1 + Λ + O(βf b )), with β << 1 and f b << 1, allowing to extract both the Rabi and Josephson frequencies with good precision. The second configuration only has one frequency which is ω 2 1 = ω 2 R (1 + Λ(1 + β/2)) which allows to isolate the value of β.
3.6.2. Stationary points with (δφ 0 a = π, δφ 0 b = 0) In this case, the condition for the existence of three stationary points is, For the case considered here,Λ ∼ Λ, and, in most applications, Λ > 1. Therefore, an appropriate choice of f a can ensure the existence of three stable points. The stability of the trivial solution is checked by studying, whose eigenvalues are listed in Table 2. The stability of the other two solutions is easy to study with the same tools. Simple analytic expressions are only attainable for the caseΛ = Λ. Then we have, whose eigenvalues are, 3.6.3. Stationary points with (δφ 0 a = π, δφ 0 b = π) The condition for the existence of three stationary points is in this case [22], The eigenvalues corresponding to small oscillations around the trivial point are listed in Table 2. Its dynamical stability depends on the specific values of f i ,Λ, Λ and ω R . For the caseΛ = Λ, it is stable provided that ω R > Λ. The eigenfrequencies for the non-trivial solution are the same as for the case (δφ a = 0, δφ b = π). For the simplest case,Λ = Λ, they are,

Effective 1D mean field approaches
In the experimental realization [4] the condensate is confined by an asymmetric harmonic trap, characterized by ω x , ω y , and ω z , with a barrier on the x direction. Thus, in a first approximation one can assume that the dynamics takes place mostly along the x axis and derive descriptions of the system where the other two dimensions have been integrated out reducing the GP3D equation to an effective 1D equation. There are different procedures to derive effective one dimensional GP-like equations starting from the three dimensional one. Their generalization to binary mixtures, with two coupled GP equations, or spinor BEC, with three or more coupled GP equations, is presented below together with the single component case.

One dimensional Gross-Pitaevskii-like equations (GP1D )
Assuming that most of the dynamics occurs in the direction which contains the barrier, the x direction in our case, one can approximate the wave function of the system by where ϕ g.s. are the corresponding ground state wave functions for the trapping potential in the y or z direction in absence of interactions (in the case of harmonic traps they are Gaussian). In this way it can be shown [35] that Ψ 1D (x; t) fulfills a Gross-Pitaevskii-like 1D equation, where the corresponding 1D coupling constant is obtained rescaling the 3D one, g 1D = g/(2πa 2 ⊥ ), with a ⊥ the transverse oscillator length, a ⊥ = /mω ⊥ , with ω ⊥ = √ ω z ω y .
The extension to binary mixtures (and also to spinor condensates [36]) may be written down readily, where, the rescaled couplings are g ij;1D = g ij /(2πa 2 ⊥ ).

Non-polynomial Schrö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 non-polynomial Schrödinger equation, proposed for a scalar BEC in Ref. [21]. 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 [20]. The NPSE for the scalar case reads, The generalization of the NPSE for two components in a binary mixture of BECs has been addressed in Ref. [37]. The system of equations, which become rather involved, can be greatly simplified in the case when all the interactions, both intra-and inter-species, are equal: , and, as before, g 1D = g/(2πa 2 ⊥ ), g ≡ g aa = g bb = g ab = g ba , and dx|Ψ j (x)| 2 = 1 .

Numerical solutions of the 3D Gross-Pitaevskii equation: single component
Before analyzing the binary mixtures in the next section, we will present here numerical results for the single component to illustrate the main differences between the various two-mode models and 1D reductions.
As discussed in the introduction, we consider the same setup and the same trap parameters as in the experiments of the Heidelberg group [4]. There, a condensate of 87 Rb with 1150 atoms is confined to a fairly small region of ∼ 5 µm through the potential, with ω x = 2π × 78 Hz, ω y = 2π × 66 Hz, ω z = 2π × 90 Hz, q 0 = 5.2µm, and V 0 = 413 h Hz. In Fig. 5 we show the potential in the x direction together with the first four energy levels of the single particle Hamiltonian and the corresponding modes. The energy levels of the single particle Hamiltonian show a clear separation between the two first eigenvalues, ground and first excited state, which are almost degenerate, and the next two.
The atom-atom interaction strength is in this case, g = 4π 2 a/M. The scattering length for 87 Rb is a = 100.87a B , therefore g/ =0.04878 KHz µm 3 . Noting that the number of atoms is known up to 10% in the experiment, the relevant product, gN/ is in the range [51.22, 60.98] KHz µm 3 . Ref. [9] uses a value of 58.8 KHz µm 3 to simulate the experimental setup. This large value of gN corresponds to a situation similar to panel (c) of Fig. 4, where the possible dynamical situations we can have are: Josephson oscillations, i.e. closed orbits around the stationary point (z 0 , δφ 0 ) = (0, 0), and self trapping regimes, usually funning phase modes.
In the experiments, the system is prepared in a slightly uneven double-well potential which produces an initial population imbalance between both sides of the barrier. At t = 0 the asymmetry is removed and the BEC is left to evolve in a symmetric doublewell potential. In our numerical simulations the initial states with either δφ(0) = 0 or π are constructed in a different way than in the experiment. We build initial states which are by construction two-mode-like. First, we obtain numerically the ground and first excited states of the condensate in the double-well potential by solving the time independent GP equation (both for the 1D reductions and the 3D case), then use those to build the left and right modes, Eq. (4), and finally construct initial states of any given initial imbalance, z 0 : Ψ z 0 (r; t = 0) = αφ L (r) + e lıπ βφ R (r), with α 2 + β 2 = 1 and α 2 − β 2 = z 0 . The ground and first excited states are obtained by a standard imaginary time evolution of the equation from an initial state with the proper parity. The density profiles of the ground, first excited and left and right modes computed numerically are plotted in Fig. 1. As can be seen, the left/right modes are indeed well localized at each side of the barrier.
From these ground and first excited states we compute all the parameters entering in the S2M and I2M descriptions presented in Secs. 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 difference between both sides of the potential barrier is computed in the following way. The phase at each point at a certain time, φ(x, y, z; t), is: where the local density, ρ(x, y, z; t) = |Ψ(x, y, z; t)| 2 . Averaged densities are defined as, i.e. integrating over the z component, To visualize the phase coherence along some of the planes we define, e.g. integrating the z component, φ(x, y; t) = 1 ρ(x, y; t) ∞ −∞ dz ρ(x, y, z; t) φ(x, y, z; t) .
(55) ¶ These values compare reasonably well with the ones provided in page 33 of Albiez PhD thesis [20], there they are given in units of ω x : A/ω x = 2.43572, B/ω x = 0.0751497, C/ω x = 0.0048, and K/ω x = 0.0163.  Figure 6. The two smaller plots above depict in solid-black line the GP3D time evolution of z (left) and δφ (right), computed as explained in the text compared to the I2M predictions in dashed-red. Then we show 3D pictures complemented with contour plots, left, of ρ(x, y; t), ρ(x, z; t) and ρ(y, z; t) at three different times, 0.5 ms (upper), 30 ms (middle) and 75 ms (lower), respectively. On the right of each plot we present a contour plot of the averaged quantum phase φ(x, y; t), φ(x, z; t) and φ(y, z; t) at the same times. They correspond to the first run presented in Fig. 8, z(0) = 0.1 and δφ(0) = 0.
The phase on the left, φ L (t), is defined as, The phase on the right is defined accordingly. The way to implement the above averages over the phase has been done in the following way, φ(x, y; t) = arctan We plot ρ(x, y; t) and contour plots. On the right panels we present contour plots of the averaged quantum phase, φ(x, y; t). The phase coherence of the condensates at each side of the barrier is clearly seen.

GP3D results
In Figs. 6 and 7 we present full GP3D simulations for a Josephson regime and a running phase mode self-trapped case, respectively. These two figures clearly show two relevant aspects of the problem. First, it is clear that during the full time evolution, which covers up to t = 80 ms in the figure, the system remains mostly localized on the two minima of the potential. Therefore, the density has a two-peaked structure over the considered time period. Secondly, the atoms in each of the two wells remain to a large extent in a coherent phase during all times. This can be seen from the uniform color, constant phase, at each side of the barrier in the right panels of the figures. These two characteristics of the time evolution of the 3D Gross-Pitaevskii equation support the use of two-mode approximations. The modulation of the density profiles on the transverse direction is seen to be small, with a mostly constant quantum phase in the region populated by the atoms. This indicates that the transverse dynamics can be integrated out to a large extent, as is done in the 1D reductions discussed in Sec. 4.
The Josephson dynamics, Fig. 6, is clearly seen in the small upper panels depicting z(t) and δφ(t). They both oscillate with the same period but with a phase-shift of π/2.
A self-trapped case is shown in Fig. 7. The atoms remain trapped mostly on the left side of the trap (they start with an imbalance of z(0) = 0.7) and remain trapped in this potential-well during the considered time evolution. The coherence of the phase at each side of the potential barrier can also be appreciated in the figure, although here we should note that the right side of the barrier, being less populated, is concentrated on a smaller (x, y) domain.

Comparison between the different models
The GP3D cases described above indicate that within the configuration considered here the two commonly employed two-mode models and 1D equations are expected to be reasonable. In this section we present comparisons between the different approaches described in the previous sections: 1D reductions (NPSE , GP1D ) and two-mode models, S2M and I2M .

GP3D vs 1D reductions: GP1D and NPSE
In Fig. 8 we present the time evolution of the population imbalance for the different dynamical conditions described in Sec. 3.3, i.e. Josephson, and self-trapping. We compare the full GP3D (solid red) with the two previously described 1D reductions, GP1D (dotted black) and NPSE (dashed blue).
First, we note that the dynamics emerging from the GP3D is indeed similar to what was predicted by analyzing the S2M equations in Sec. 3 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 the initial side of the trap and there are short and small periodic oscillations as predicted by the two-mode models. At longer times, the imbalance is seen to decrease smoothly, implying a departure from the predicted two-mode dynamics [38]. The two 1D reductions give qualitatively similar results in most situations to 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 cases, failing in all cases to give the same period as the GP3D or the NPSE . Moreover, the GP1D departs notably from two-mode for the selftrapped case. It does predict self trapping, but more than two modes contribute to the time evolution. These results justify the use of the NPSE in Ref. [4] to analyze their experiment.
To further explore the quality of the 1D reductions, we present in Fig. 9 the density profiles in the x direction after integrating the y and z components, ρ(x; t) = ∞ −∞ dy ∞ −∞ dz |Ψ(x, y, z; t)| 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 , as could be inferred from the previous results, does not predict the correct density profiles and, as seen in the self-trapped case, (z(0) = 0.6, δφ(0) = 0), do show the 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 Ref. [9], z c = 0.39, and differs from the one reported in Ref. [4], z c = 0.5.

GP3D vs two-mode approximations, S2M and I2M
As explained above, the use of two-mode models is suggested by the GP3D results, see Figs. 6 and 7. What is, a priori, not clear, is whether the extra assumption used in deriving the S2M (which are the most commonly employed equations) will work for each specific double-well potential. As discussed in Sec. 3.2, the conditions of the Heidelberg experiment are such that the S2M predictions are not good. This does not mean that the dynamics is not two-mode but that the overlaps involving high powers of the two localized modes are not negligible as assumed in deriving the S2M equations.
In Fig. 10 we compare GP3D (solid red), the S2M (dotted black) and the I2M (dashed blue) results using the parameters calculated microscopically from the ground and first excited state of the GP3D . Both two-mode schemes predict the same phenomenology and thus qualitatively capture the dynamics of the system. At the quantitative level, however, the I2M is clearly better. In the run with z(0) = 0.1 and δφ(0) = 0. (panel (a)), both the S2M and I2M predict a similar behavior with the correct amplitude and oscillation period close to the GP3D one. As the imbalance is increased, e.g. ( [4] considers z(0) = 0.28), the S2M fails to describe the correct period and predicts smaller amplitudes. This is analyzed in full detail in Ref. [9]. The critical initial imbalances determined by both two-mode approaches are smaller than the GP3D one. For the latter they predict a self-trapped case, see panel (b). Finally, for the self-trapped cases with δφ(0) = π (panel (c)) the I2M give similar oscillation amplitudes with shorter periods than the GP3D . The S2M fails both in reproducing the amplitudes and the periods.

Numerical solutions of the 3D Gross-Pitaevskii equations: binary mixture
As discussed in Sec. 2, one feasible way of experimentally prepare binary mixtures of BECs is to consider a number of atoms populating the m = ±1 Zeeman components of an 87 Rb F = 1 spinor. The experimental observation of Josephson tunneling phenomena by the Heidelberg group seems to be possibly extended to trap both Zeeman components [34]. 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 inter-species 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 [39], with H s = − 2 /(2M) ∇ 2 + V + c 0 n being the spin-independent part of the Hamiltonian. The density of the m-th component is given by n m (r) = |ψ m (r)| 2 , while n(r) = m |ψ m (r)| 2 is the total density normalized to the total number of atoms N. The couplings are c 0 = 4π 2 (a 0 + 2a 2 )/(3M) and c 2 = 4π 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 [40]. Since the spin-dependent coupling, c 2 , is much smaller than the spin-independent 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 [30]. Therefore, in our calculation the number of atoms in each sublevel remains constant in time allowing to treat the system as a real binary mixture of components a and b. Comparing the system of Eqs. (1) and (58) 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 want to investigate the Josephson-like  Table 3 for the specific conditions considered in the numerical simulations as a function of the fraction of atoms in the a component, f a . The notation is as follows, ω i;αβ , with i = 1, 2 and α, β = 0, π. (right) The conditions for the existence of the non-trivial equilibrium points given in Eqs. (38,41,45), upper panel, as function of f a for the conditions described in the text. The lower panel contains the explicit equilibrium points z 0 a , z 0 b as a function of f a obtained by solving equations (34). Note that each equilibrium point has a trivial partner which is obtained by flipping the sign of z 0 a , z 0 b .
dynamics for different number 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). The values of Λ = NU/ ω R andΛ = NŨ / ω R are Λ = 74.278 andΛ = 74.968. With Λ/Λ = 0.99. These are obtained from the microscopic 3D parameters computed in the scalar case, with the same total number of particles, Sec. 5. This is reasonable for the case we are considering where g aa = g bb ∼ g ab , which implies that the ground state wave functions for the GP equations of the mixture do not depend on f a and f b for a fixed total number of particles. This would certainly not be the case if g aa = g bb = g ab , in such case one would need to recompute the ground state wave functions for a and b for each value of f a .
Following the discussion in Sec. 3.6, where the predictions of the S2M were discussed in detail, the system has the trivial equilibrium points, listed in Table 3 with β = 0.009. In Fig. 11 we show the values of the two eigenfrequencies for each of the trivial equilibrium points listed in Table 3 for the specific conditions described above. The figure shows a number of important features about the stability of the trivial equilibrium points. First, the (z 0 a , δφ 0 a , z 0 b , δφ 0 b ) = (0, 0, 0, 0) is always stable regardless of the total polarization of the system (measured by f b −f a ). Second, the (z 0 a , δφ 0 a , z 0 b , δφ 0 b ) = (0, π, 0, π) mode is always unstable, as seen by the negative value taken by the square of the frequencies. Third, the (z 0 a , δφ 0 a , z 0 b , δφ 0 b ) = (0, 0, 0, π) mode should be stable for f a 0.43, correspondingly the (z 0 a , δφ 0 a , z 0 b , δφ 0 b ) = (0, π, 0, 0) is stable for f b 0.43 and therefore there is a range of polarizations, given by 0.43 f a 0.57 where the only trivial mode which is stable is the (z 0 a , δφ 0 a , z 0 b , δφ 0 b ) = (0, 0, 0, 0). The non-trivial equilibrium points in this case can be obtained by analyzing the conditions given in Sec. 3.6. For (δφ 0 a , δφ 0 b ) = (0, 0) there are no equilibrium points apart from the trivial one, due to Λ ∼Λ. In the other three cases there are nontrivial equilibrium points depending on the specific values of f a . In Fig. 11(right) we analyze their existence. First, we note that there are non-trivial points corresponding to (δφ 0 a , δφ 0 b ) = (0, π) provided f a 0.37, correspondingly there are also equilibrium points for (δφ 0 a , δφ 0 b ) = (0, π) if f b 0.37. There is also a non-trivial equilibrium point for (δφ 0 a , δφ 0 b ) = (π, π) regardless of f a . As can be seen in the figure, all these nontrivial equilibrium points correspond to fairly imbalanced conditions and can in most cases be understood in simple terms from the analysis of the scalar case. For instance, the equilibrium point for (δφ 0 a , δφ 0 b ) = (π, π) corresponds to z 0 a ∼ z 0 b ∼ 1 (or −1), which can be understood as having both components locked in a π-mode. Similarly, the equilibrium points in the (0, π) or (π, 0) cases exist whenever the most abundant component is populated enough to drive the dynamics close to being π locked.

GP3D calculations: phase coherence and localization
The numerical solutions of the GP3D presented in Sec. 5 for the single component case 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 BEC, one at each side of the barrier, with a well defined phase at each side during the dynamical evolution. Essentially those are the premises used to derive the two-mode models, both for single component and for binary mixtures, as we did in Sec. 3.
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 two distinctive features, see Fig. 12. First, the density of atoms for each component is always bi-modal, with the two atom bunches centered 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 behavior that are found in very specific conditions, e.g. for very large population imbalances and also when analyzing a case with g ab = g aa = g bb .
6.2. Small oscillations around z 0 a,b and δφ 0 a,b = 0 The two predictions of the S2M described in Sec. 3.6 are confirmed by the NPSE and GP1D simulations as can be seen in Figs. 13 and 14. In Fig. 13 (left panels) we consider a very polarized case, f a = 0.8. As expected from the two-mode 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 + Λ. The GP3D simulation is seen to confirm the above and follow closely the predictions of the I2M . The less abundant component is strongly driven by the most populated one and shows an anti-Josephson behavior as described in Ref. [30].
Another prediction is related to the behavior of z a + z b and z a − z b in the nonpolarized case, f a = f b . As explained in Sec. 3.6, 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. In the right part of Fig. 13 we present the extreme case when z a (0) = −z b (0) computed with GP3D , NPSE and I2M . 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. As seen in Fig. 14 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, see Sec. 5.
Interestingly, they predict different Josephson oscillations while the Rabi frequencies are similar. In panel (c) of Fig 14 the long oscillation corresponding to the Rabi mode is seen to agree well with the corresponding long oscillation seen in the right panels of Fig. 13. 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 .
As a general statement, in the conditions of the Heidelberg experiment, as occurred for the scalar case, the I2M produces more reliable results than the S2M model, which are not shown in the figures. Notice that the parameters that we use for the I2M are extracted from the GP3D calculation as given in Sec. 5. Other representative cases with (δφ a (0), δφ b (0)) = (0, 0) but with larger initial imbalances, z i (0) ∼ 0.5 are shown in Fig. 15 initial population imbalances. In this case, both frequencies ω J and ω R show up in the evolution. The I2M provides a satisfactory description of the dynamics.
6.3. Small oscillations around z 0 a,b , δφ 0 a = 0 and δφ 0 b = π As explained above, for these conditions there can exist up to three stationary points depending on the specific value of f a considered. The trivial equilibrium point exists provided f a 0.43, see Fig. 11. This prediction of the two-mode models is observed both in the GP3D and NPSE as it can be seen in Fig. 16. In the figure, we consider a simulation with z a (0) = 0.1, z b (0) = −0.15, and f a = 0.2 < 0.43 (left panels). The population imbalance (upper panel) of both components oscillates in the usual Josephson regime. At the same time, the phase difference oscillates with its characteristic phaseshisft of π/2 with respect to the imbalance (lower panel). The phase of the a component oscillates around δφ a = 0 while δφ b does oscillate around δφ b = π.
A completely different picture emerges when the fraction of atoms in both components is exchanged, f a = 0.8 > 0.43 (right panels), with most of the atoms populating the a component. In this case, the oscillation amplitude is large, both components remain trapped on their original sides and the phase difference becomes unbounded. This should be considered as a genuine effect of the binary mixture as each component follows a running phase mode at each side of the potential barrier. 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 80 ms in the run without equilibrium point (right panel).
The I2M gives a good qualitative picture of both cases but fails to provide predictions as accurate as the NPSE , as happened in the scalar case, see for instance Figs. 8 and 10. In particular the predicted periods of oscillation are much longer than the actual ones.
An example of simulations around non-trivial equilibrium points is presented in Fig. 17. As explained previously, these involve very large and opposite initial population imbalances for both components. In Fig. 17 we consider a case with initial conditions very close to the predicted equilibrium point using the standard two-mode, and described in Fig.11, 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. These simulations will be discussed again in Sec. 6.5 as they exhibit effects which clearly go beyond a two-mode approximation. 6.4. Small oscillations around z 0 a,b and δφ 0 a,b = π The trivial equilibrium point is not stable in the considered conditions as seen in Fig. 11. The non-trivial one, however, is only attainable if extremely imbalanced configurations for both components are considered. This case would correspond essentially to having both components in a π mode state, which in our conditions only exists for z ∼ 1 as can be seen in the blue spots in panel (c) of Fig. 4. In Fig. 18 we present two simulations with different initial conditions. First, we consider a simulation with z a (0) = 0.4 and z b (0) = −0.2, with f a = 0.9. The behavior is understood in simple terms, the most populated component remains self-trapped while the other component is forced by the other one. The phase evolves unbounded. The figure again contains GP3D and NPSE simulations.
The second simulation (right panels) is closer to a non-trivial equilibrium point, we consider z a (0) = 0.9 and z b (0) = 0.85 with f a = 0.9. In this case, both components remain self trapped, the phase difference is unbounded, but we do not get the expected behavior of two π modes because the initial imbalances are not close enough to z 0 ∼ 1.

Effects beyond two-mode
Most of the dynamics described in the previous sections can to a large extent be understood within the two-mode models developed in Sec. 3. 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 double-well 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 seen in Fig. 17. There, as clearly seen in the density profiles along the x direction, the two-mode approximation is clearly not valid. The simplest way of seeing this is by noting the zero in the density of one of the components at x ∼ 2µ m. 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, (b), of effects beyond two-mode involve excitations of the transverse components. These effects are present in any binary mixture calculation whenever the intra-and inter-species interactions are not equal. To enhance this effect, and also to explore the interesting symmetry breaking phenomena described in Ref. [25], we consider a case with g aa = g bb , but with g ab = g ba = 2.3g aa . Therefore, now the inter-species interaction strength is larger than the intra-species one. The two-mode prediction for this case, S2M , which was analyzed in Ref. [25] shows a large symmetry breaking pattern during the time evolution of the system. In Fig. 19 we consider a full GP3D simulation of a representative example with z a (0) = −0.2, z b (0) = 0.1, δφ a (0) = δφ b (0) = 0, and f a = 0.7.
The qualitative prediction of the I2M also shows the symmetry breaking, and the two components do separate from each other and mostly concentrate on one of the wells as time evolves. But, as it can be seen in the 3D depictions of ρ(x, y; t) at three different times, the evolution of the system departs almost from the beginning from the twomode. At t = 1 ms we have the density distributions of each component corresponding to a small initial imbalance. Then at t = 11 ms, we can already see that the most populated component is expelling the other one from the minima of the potential. This fact can be appreciated as a four peaked distribution, ρ b (x, y; t). After that, each of the components start to accumulate on their original sides following qualitatively the prediction of the I2M and thus presenting the symmetry breaking pattern discussed in Ref. [25]. The two-mode approximation is in this case broken for a short period of time, when the first modes along the transverse directions are excited due to the large inter-species interaction. As explained in the text in this case, g aa = g bb and g ab = 2.3g aa . The upper two plots depict z a (t) (left) and z b (t) (right). Then each row contains from left to right: 3D depictions complemented by contour plots of ρ a (x, y; t), a contour plot of the averaged phase φ a (x, y; t), 3D depictions complemented by contour plots of ρ b (x, y; t), and a contour plot of the averaged phase φ b (x, y; t). Each row corresponds to a different time, 1 ms (upper), 11 ms (middle) and 51 ms (lower), respectively. In all cases, solid black lines are computed with GP3D and dashed red ones with I2M .

Conclusions
We have presented a thorough investigation of the mean-field dynamics of a binary mixture of Bose-Einstein condensates trapped in a double-well potential. 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 scrutinized by performing full GP3D simulations covering all the relevant initial conditions. The 3D numerical solutions of the Gross-Pitaevskii equations have been used as a benchmark to critically discuss the validity of the most common 1D reductions of the GP equations, GP1D and NPSE , and the often employed simple two-mode reductions, S2M and I2M .
The full 3D solutions of the binary mixture have shown to have a large amount of phase coherence and localization at each side of the potential barrier for both components, predicting a dynamics which is mostly bi-modal. This feature permits to speak of Bose-Einstein condensates at each side of the barrier, where the atoms mostly share a common phase, and to support the use of two-mode approximations, which analytical solutions allow to gain physical insight into the problem.
To fix the conditions of the dynamics, we have focused in one particular setup that corresponds to a natural extension of the experiments reported in Ref. [4]: 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 to observe and characterize a large variety of phenomena which are genuine of the binary mixture, e.g. anti-Josephson oscillations in highly polarized cases, long Rabi-like oscillating modes, zero and π locked modes, etc.
For the sake of completeness and to better frame the physics of the binary mixture we have provided a detailed description of the single component dynamics, with explicit expressions for all the commonly employed approximations to the 3D mean field Gross-Pitaveskii equation. The natural extension of the latter to the binary mixture, i.e. S2M and I2M equations and 1D reductions, have been consistently derived providing a self-contained reference, easy to read, with all the relevant formulae used in the article.
The standard two-mode model, with its microscopic parameters computed with the GP3D , has been used to reexamine the existence and stability of the different regimes that can occur in both single component and binary mixture condensates, describing the Josephson oscillations and the macroscopic quantum self-trapping, including running phase modes and zero-and π−modes.
The comparisons between the two-mode models and the numerical solutions of the GP3D show an excellent agreement for conditions close to the stable stationary regimes predicted by the two-mode models. As we depart from those stable points, the S2M fails to provide a quantitative agreement with the results obtained with the GP3D equations. The range of validity of the I2M is much larger, fully capturing the dynamics of single and binary mixtures for a larger set of initial conditions.
The two most 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 in 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 behavior 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 intra-and inter-species are very similar and the NPSE equations are particularly easy to handle. This agreement is not only seen on fully integrated magnitudes, for instance population imbalances, but also on the density profiles predicted along the direction hosting the barrier.
We have also considered two situations where the two-mode approximation fails. This is naturally due to the excitation of higher modes. Two different cases have been described, first the excitation of modes in the direction of the barrier and secondly, excitation of modes in the transverse direction. The NPSE has been shown to capture perfectly the excitations along the barrier direction, reproducing the integrated density profiles obtained with the GP3D . The second case has been studied in a simulation performed with different intra-and inter-species, which can be achieved in principle experimentally through Feshbach resonance modulation of the scattering lengths. In this case, the dynamics of the less populated component in each side of the trap departs notably from the two-mode with clear excitations of transverse modes, seen already in the density profiles along a transverse direction.
The present article is intended both to motivate the experimental effort to study binary mixtures of BECs, where we have shown that a large variety of phenomena related to phase coherence and localization can be observed, and to serve as a tool in the analysis of such experiments providing a concise and self-contained derivation of the most commonly used models.