EXACT SOLUTION OF THE TWO-MODE MODEL OF MULTICOMPONENT BOSE-EINSTEIN CONDENSATES

We find the explicit solution of the so-called two-mode model for multicomponent Bose-Einstein condensates. We prove that all the solutions are constants or periodic functions and give explicit formulae for them.

1. Introduction. One of the most remarkable achievements of contemporary physics has been the realization of Bose-Einstein condensation with ultracold atomic gases [1]. In the regime of very low temperatures, the gas dynamics is modelled by the so called Gross-Pitaevskii equation, which reads where m is the mass of the atoms, is the Planck constant, U is a parameter related to the probability of collisions between individual atoms and ∆ = D j=1 ∂ 2 ∂x 2 j is the Laplace operator in D dimensions. The quantity |Ψ(t, r)| 2 has the physical meaning of the local density of the gas and the total number of atoms, which is given by N := Ψ 2 L 2 = d D r|Ψ(t, r)| 2 , is conserved under time evolution. The potential V (r) = 1 2 mω 2 r 2 , where r 2 = D j=1 r 2 j and ω ∈ R + , is an external potential generated by magnetic fields. The presence of this potential is responsible for the physical confinement of the atoms which would otherwise spread out due to the combined effect of dispersion (term with Laplacian) and the nonlinear interaction (which is repulsive).
These type of systems are naturally three-dimensional, however, in many situations in which the atoms are tightly confined along one of the dimensions, it has been shown that a two-dimensional description of the problem is possible. The only change is that the scattering parameter U and number of particles N are then understood as effective parameters proportional to the real ones [2].
Although our analysis will be related to dynamical phenomena in Bose-Einstein condensates equation (1) is an equation of nonlinear Schrödinger type and thus our results to be described later will have a much larger range of applicability because of the well-known universality of this model equation [3].
In this paper we concentrate on the analysis of multicomponent systems. When two different type of atoms of the same mass (usually two different hyperfine levels of the same atomic specie), which occupy the same spatial region are simultaneously cooled down below the critical temperature, each type of atom is described by a single function Ψ 1 (r, t), Ψ 2 (r, t) so that the system is described by the equations and the number of particles of the jth type is given by N j = d D r |Ψ j (t, r)| 2 . These equations are a generalization of Eq. (1) to the case of two atomic species. Here U jk are constants controlling the nonlinear behavior proportional to the mutual collision amplitudes and U 21 = U 12 .
It is convenient to work with a set of adimensional units. To this end we scale coordinates and time by r → r 0 r,t → t/ω with r 0 = /mω, (hereafter r and t stand for dimensionless coordinates and time) and introduce new functions ψ j by with u jk = U jk r D 0 / ω. These type of physical systems were first generated experimentally in 1997 [4]. One of the most remarkable achievements using these systems was the first experimental generation of vortices in Bose-Einstein condensates reported in Ref. [6].
Let us now present the concept of solitary wave or stationary solution of Eq. (1) which in adimensional units can be written as By definition solitary waves are solutions of the form ψ(t, r) = φ µ (r)e −iµt where the pair µ ∈ R and φ µ (r) satisfy the elliptic problem with boundary conditions lim r→∞ φ µ (r) = 0. These solutions are stationary points of the functional of fixed L 2 -norm. Of the full family of solutions φ µ (r) of Eq. (5) the one which corresponds to the global minimum of the functional (6) is called the ground state and we will denote it by φ 1 (r).
In the linear case (U = 0) this system has a numerable set of solutions. Say, in the D = 2 case this set can be chosen as φ kl (r) = H k (x)H l (y)e −(x 2 +y 2 )/2 , where µ kl = k + l + 1, k, l ∈ N and H k,l are the Hermite polynomials. For the nonlinear repulsive case (U > 0) it is has been proven in limited situations [7] and conjectured for the general case that the situation is the same, i.e., for each value of the norm N there is a numerable set of solutions {µ n , φ n }.
It is known that the ground state of Eq. (5) is a unique positive radial solution [8] which is finite when r → 0. General properties of the solutions of Eq. (5) have been studied in a more general context in [9,10].
When the nonlinear parameter is large enough [17] it has been checked numerically that the first excited state, φ 2 (r), is a solution of the so-called vortex type. The simplest vortex solutions of Eq. (5) in the two-dimensional case is of the form φ 2 (r) = R 2 (r) e iθ , where µ 2 ∈ R, R 2 (r) : R + → R and (r, θ) are the polar radius and angle. The solutions µ 2 , R 2 (r) are to be determined from the equations It is possible to prove that this equation has regular solutions satisfying the condition that R 2 (r) ∼ r when r → 0 and R 2 (r) ≤ re −r 2 , r → ∞. The name vortex comes from the fact that if the Bose-Einstein condensate is a superfluid then the gradient of the nonhomogeneous phase term is physically associated to a closed circulation of the quantum fluid (see e.g. [11] for details). If Θ = arg ψ the term vortex refers to solutions for which m = 1 2π ∇Θd is an integer nonzero value for any closed curve γ containing a point -the so called vortex core -(x 0 , y 0 ) where ψ(x 0 , y 0 ) = 0.
There has been a great interest on multicomponent systems in which one of the atomic species has initially a vortex, i.e., the corresponding ψ j (r, 0) is of vortex type and ψ(r, 0) with = j is a nodeless function. This interest was triggered by the experimental results reported in Ref. [6] where one of the possible configurations with j = 1, 2 was found to be stable and the other one unstable. In recent works [12,13] numerical simulations of Eqs. (3) have been used to show that the dynamics of the unstable configuration can be understood within the framework of mean field theories for the double-condensate system. Later [14], it was shown that a simple model was able to describe accurately the dynamics of the multiple-condensate system in a particular two dimensional situation and an specific parameter region close to the experimental parameters of Ref. [6].
In this paper we concentrate on the analysis of the simple model reported in Ref. [14] proving that the system is integrable and that all the solutions are either periodic or constants (equilibria). This proves that the tendency of the system to exhibit periodic exchange of vortices between both atomic species shown in Refs. [12,13,14] is not the consequence of a fortunate choice of parameters. We also study the period of the solutions as a function of the parameters of the problem and compare quantitatively the analytical solutions with numerical simulations.

2.
The two-mode model in two dimensions. We will derive here the twomode model following a different formalism than that of Ref. [14]. The idea of the method is to use the reduced Lagrangian approach. Equations (3) can be obtained from the variational problem where (hereafter we restrict ourselves to the two-dimensional case D = 2) with In terms of some basis, associated with the linear problem the Lagrangian can be rewritten as follows: In what follows we will use a limited expansion, i.e. we keep in (11) only two modes: where φ 1 is the mode resembling the ground state and φ 2 corresponds to a state with a vortex. Is this approximation accurate enough to describe the dynamics contained in Eqs. (3)? It is evident that in a general situation the answer is negative but in this paper our interest is in the dynamics of special configurations having initially a vortex in any of the components plus ground state in the other one. Thus, our choice for is the simplest choice which is able to approximate the initial data of interest for our physical situation. Moreover, numerical stability analysis of the full partial differential equations ruling the phenomena [12] shows that when vortex configurations are unstable the destabilization is caused by the growth of gaussian (ground-state) modes. This numerical evidence makes us expect that our choice will be able to describe at least the stationary states and the first stages of the unstable cases. Finally, the fact that the dynamics obtained from the dynamics of the full PDE system has some sort of recurrent dynamics leads us to the hope that the phenomena could more accurately described by this crude approximation than one could expect in principle. In fact, despite its simplicity and crudeness this approximation has been shown to retain many qualitative features of the dynamics of the system [14]. Further comparisons with the numerical results obtained from simulations of the full problem (3) to be presented later in this paper will also support our statement.
In this paper we make a more detailed derivation and mathematical analysis of this model together with a deeper physical analysis of its range of validity and implications.
Inserting (15) into Eq. (14) and noting that φ 1 is axially symmetric while φ 2 ∝ exp(iθ) one can obtain that of all 16 integrals in the nonlinear part of the Lagrangian only 6 are nonzero, namely the ones which correspond to the following values of the indices: (j, k, j , k ) = (1, 1, 1, 1), (1, 2, 1, 2), (1, 2, 2, 1), (2, 1, 1, 2), (2, 1, 2, 1), (2, 2, 2, 2). Thus the truncated Lagrangian becomes Namely this Lagrangian we consider as describing the two-mode model in two dimensions. As discussed before, this model allows a very complete description of the dynamics of the system and will be the starting point for our discussion. In particular we will try to understand why most of the simulations of the two-mode model and the full PDE system (3) lead frequently to periodic solutions. Specifically our main result, apart from obtaining the explicit form of the solutions, will be the proof that all solutions of this system are periodic functions of time and the analysis of the dependence of this period on the relevant parameters of the problem.
3. Analysis of the two-mode model. (17) with respect tō a nj , n, j = 1, 2 we get the evolution equations for the coefficients a nj :

Reformulation of the model. Variating Lagrangian
where while the equations for phases are +C 12 u nn ρ 2 n + C 12 u 12 These equations must be complemented with the appropriate initial conditions for ρ kl and Φ. ρ nj (t 0 ) =ρ nj θ nj (t 0 ) =θ nj n, j = 1, 2 (22) The form of the two-mode equations [Eqs. (19), (21)] presented here is slightly different from that of Ref. [14], but is more suitable for our purposes in this paper.
and the initial data areΓ ≡ Γ(0) = nj γ njρ 2 nj . The following relations are then evident: withΓ = nk γ nkρ 2 nk . These equations imply that all the relevant densities may be obtained from a single quantity, Γ(t), and our goal now will be to find an equation ruling its dynamics. Inserting representation (33) where and the constants Γ nk , n, k = 1, 2 are defined as Γ 11 =Γ + γ * ρ equations, allows to obtain the expressions for ρ jk (t). In particular Eq. (38) leads to .
Let us rewrite the polynomial (39) in terms of its roots as follows where P * = 4C 2 12 u 4 12 /γ 2 * − 1/4 and P 1 , P 2 , Q 1 , Q 2 are the roots of P 4 . There are at least two real roots P 1 , P 2 of P 4 satisfying thatΓ is either a root or is located between two of them, i.e. P 1 <Γ < P 2 .
To prove the last affirmation let us first compute P 4 (Γ) = 4C 2 12 γ 2 * ρ11ρ12ρ21ρ22 sin 2Φ > 0. Note that if the initial data are such that ifΦ = 0 thenΓ is a real root of P 4 (Γ).
To understand what is the behavior when Φ = 0 we evaluate This means that there is at least one real root located in each of the intervals As a corollary we get that P 4 (Γ) > 0 for P 1 < Γ < P 2 . Concerning the other two roots Q 1,2 (Q 2 > Q 1 if real), nothing may be said in general. However, for the case P * > 0 it is clear that P 4 (±∞) > 0 and then Q 1,2 are real numbers. When P * < 0 it happens that Q 1 , Q 2 could be (and, in fact, they are) complex conjugate roots for certain parameter ranges. In any case, the boundedness of Γ ∈ [P 1 , P 2 ] ensures the periodicity of the solutions.
Let us now proceed to find an explicit form for the solutions. First, to present polynomial P 4 (Γ) in the canonical form we make the transform: where Equivalently we may also write with the quantities ∆ and P 0 being given by It is clear that P 1 < P 0 < P 2 and thus P(P 0 ) > 0. In terms of Y the polynomial P 4 (Γ) can be written as where Thus, the function Y (t) satisfies the equation with Equation (52) can be solved as follows: where the constantt should be determined from the initial conditions by solving the equation This allows us to get a explicit solution for our problem. To do so let us first use Eq. (54) to find Thus, inserting this explicit expression for Γ(t) and Eqs. (33) one may get explicit expressions for ρ nj (t), n, j = 1, 2. The function sn is periodic, and then the period of our solutions is where the quantity K(k) = π/2 0 dθ/ 1 − k 2 sin 2 θ is the complete elliptic integral of first kind [16].

Comparison of the two-mode model with the results from the PDEs.
The explicit forms of the solutions given by (56) together with (33) and the formula for the period (57) are the main results of this paper. Although we tried to justify a priori the particular choice that we have used to model the initial problem, a more detailed comparison is in order.
In general, we must say that for the parameter ranges which we have explored in more detail the agreement of the two-mode model with the simulations of the partial differential equations (3) is good. For instance, in Fig. 1 it is shown the period of the vortex transfer as a function of the initial conditions ρ jk (0) for a specific choice of the parameters with a very good result for the two-mode model. The specific values of ρ jk (t) in particular situations are not in such an excellent quantitative agreement but at least the qualitative behavior is captured by the simple model as seen in Fig. 2. In Fig. 2(a) (upper solid curve) we observe that the two modes used capture most of the norm of the solution, which is consistent with the very good approximation to the period provided by them. However, there are some inaccuracies which also are present in other simulations which we have performed: (i) the instability tends to set in faster in the two-mode model than in the full system and (ii) the fraction of particles transferred is smaller in the PDE model. In any case, even though these details differ, the inter-peak distance which gives the period is very similar between both models, as was expected from the results shown in Fig. 1.
Of course, when the ρ jk (0) or u jk constants are very large, corresponding to a strongly nonlinear situation then the model, which is based on an expansion on the basis of modes of the linear problem, is expected to fail. Another situation where we have found the model to fail is for medium to large perturbation values for which the periods are overestimated by a factor as large as 4.

Analysis of the two-mode equations.
A systematic analysis of the dependence of the period of the solutions of the two mode model is a formidable problem due to the large number of parameters involved. In what follows we study several specific examples.
dependence of the period on the phase Φ ( Fig. 3(b)). Since in normal experimental situations, the perturbations are not controlled but generated by imperfect initial data, noise, etc., this sensitivity could limit the applicability of the model (3) to predict the specific value of the transfer period. On the other side, by controlling the distance of the initial configuration to the stationary (unstable) one, it is possible to achieve a wide variety of periods for given physical parameters.
We have also studied, as shown in Fig. 4, the dependence of the period as a function of the values of physical constants u ij for different values of the initial data ρ ij (0). The most remarkable feature is the strong increase of the period when u 22 is much smaller than u 11 . This phenomenon is due to the small coupling between the two populations, which leads to slower transfers. 4.3. Other applications. One interesting application of the formulae from the two-mode model is to give the period for small coupling values. Due to the presence of terms of the type e iµt for µ of the order of one in the solutions ψ 1 , ψ 2 the time step in the numerical simulation cannot be chosen to be too large. On the other hand, for small coupling values the period of the transfer becomes very large (as can be seen for instance in Fig. 1). Thus, the solution has two time scales one small and other very large, which makes it difficult to extract the period from numerical simulations of Eq. (3). In this situation Eqs. (56) and (57) can be a helpful tool to provide results. Another interesting applications of the two-mode formulae is the analysis of the vortex transfer mechanisms in three-dimensional systems. Again, in such problems the predictive power of Eqs. (3) is limited due to computational limitations. On the other hand the two mode model can be used by just changing the coefficients E j , C jk by their three-dimensional counterparts.

5.
Conclusions. In this paper we have found the exact solution of the two-mode model developed for the explanation of the vortex transfer mechanisms described in Refs. [6,12,14]. The most remarkable result is that all the solutions are periodic functions and thus the periodic transfer mechanisms of the vorticity are natural within the range of validity of the model, which was previously found to be at least that of pancake traps and some regimes of fully three-dimensional traps. It has been also shown the relevance of the different parameters on the frequency of the oscillations. In general, this dependence is nontrivial and given by the explicit formulae here found