Exact Dynamics of Multicomponent Bose-Einstein Condensates in Optical Lattices in One, Two and Three Dimensions
aa r X i v : . [ c ond - m a t . o t h e r] N ov Exact Dynamics of Multicomponent Bose-Einstein Condensates in Optical Lattices inOne, Two and Three Dimensions
R. Mark Bradley
Department of Physics, Colorado State University, Fort Collins, CO 80523, USA
James E. Bernard and L. D. Carr
Department of Physics, Colorado School of Mines, Golden, CO 80401, USA (Dated: November 6, 2018)Numerous exact solutions to the nonlinear mean-field equations of motion are constructed formulticomponent Bose-Einstein condensates on one, two, and three dimensional optical lattices. Wefind both stationary and nonstationary solutions, which are given in closed form. Among thesesolutions are a vortex-anti-vortex array on the square optical lattice and modes in which two ormore components slosh back and forth between neighboring potential wells. We obtain a variety ofsolutions for multicomponent condensates on the simple cubic lattice, including a solution in whichone condensate is at rest and the other flows in a complex three-dimensional array of intersectingvortex lines. A number of physically important solutions are stable for a range of parameter values,as we show by direct numerical integration of the equations of motion.
PACS numbers: 03.75.Lm, 03.75.Kk, 03.75.Mn
I. INTRODUCTION
Two areas at the forefront of research in Bose-Einsteincondensates (BECs) over the last few years have beenoptical lattices and the hyperfine degree of freedom [1,2]. Optical lattices allow one to explore, in both themean-field and quantum regimes, the effects of periodicpotentials on bosons. These studies complement the vastbody of knowledge concerning the behavior of fermionsin periodic potentials coming from solid state physics.Optical lattices are free of defects and disorder and theatomic potential has a simple closed form. In contrast,an electron in a solid is subject to a complex, imperfectlyknown potential that is usually marred by defects.The hyperfine states of the atoms making up a BECallow one to construct exotic spin structures based onthe occupation of different hyperfine spin states of theform | F, m F i . BECs of this kind, which are called spinoror multicomponent , have a vector order parameter. Re-search on multicomponent BECs has been instrumentalin reaching important milestones in BEC research, in-cluding the creation of a quantum vortex and the sub-sequent demonstration that a BEC made from a weaklyinteracting alkali gas is superfluid [3, 4]. Recently, exper-imentalists have placed multicomponent BECs in opticallattices [5, 6].In this article, we construct exact solutions to themean-field equations of motion for multicomponentBECs in one, two and three dimensional optical lattices.Band theory was invented as a tool to analyze station-ary solutions of the linear Schr¨odinger equation with aperiodic potential, a problem that does not have a so-lution in closed analytic form for realistic potentials. Incontrast, for the nonlinear Schr¨odinger equation (NLS)with a sinusoidal optical potential, exact, closed-formstationary solutions have been discovered [7–14]. Herewe generalize and extend previous work to an overar- ching and rigorous treatment of certain classes of exactsolutions describing the dynamics of condensates with s components in D dimensions. Our formalism permits usto construct exact, nonstationary solutions of the vec-tor NLS. This article brings together our new work andpast treatments into a single, general, rigorous analyti-cal framework. At the same time, we elucidate the mostexperimentally relevant and aesthetically pleasing solu-tions. We also perform a full nonlinear stability analysiswhere computationally tractable. A number of surprisingstability regimes present themselves, as we will demon-strate.The basic idea that leads to our exact solutions is touse the nonlinearity to cancel the spatial variation of thepotential, leading to an effective free particle problem.(Clearly, this is not possible for the linear Schr¨odingerequation.) This idea is due to Bronski, Carr, Deconinck,and Kutz, who originally applied it to a one-component,or scalar, BEC in a one-dimensional periodic potential [7–9]. Deconinck, Frigyik and Kutz extended this work onone-component BECs to higher dimensions [10, 11], andto multicomponent condensates in one dimension [13],but not to multicomponent BECs in two and three di-mensions. Although the higher-dimensional Jacobi el-liptic periodic potentials considered by Deconinck et al. do not include the square, rectangular and simple cubicoptical potentials readily available in the lab, Hai andcoworkers have shown that the cancellation technique canbe used to construct solutions for a scalar condensate ona square optical lattice [12].Our extension of the cancellation technique to conden-sates with an arbitrary number of components s in a sin-uoidal optical potential of arbitrary dimension D leads toan enormous number of new solutions, many of great ex-perimental import. A crucial and challenging step in thecancellation technique is to find a solution ansatz that,for a given potential, allows the cancellation to occur. Atthe same time, the solution ansatz must satisfy the freeparticle Schr¨odinger equation. A key aspect of our workis the introduction of a novel, very general solution ansatzthat allows the solution of a wide range of problems.Mean-field theory, which for a single-component con-densate takes the form of the scalar NLS or Gross-Pitaevskii equation [15], has been quite successful in de-scribing experiments on multicomponent BECs in opticallattices [5, 16]. However, there are important approxima-tions underlying its use. First, the tunneling or hoppingenergy t h must be much larger than the on-site interac-tion energy U so that the system is far from the Mottinsulating regime [17, 18]. This means that the poten-tial barriers are not so high that the sites lose mutualphase coherence, and that a full many body quantumFock state treatment is not needed [20]. Second, three-body and other loss processes are neglected. Third, quan-tum fluctuations are ignored, as is always the case whenthe NLS is applied to BECs [15, 20]. Fourth, any possibleresonances induced by the lattice or dimensional confine-ment are neglected [21]. Fifth, when treating D = 1and D = 2, mean-field theory requires that the confiningpotential that reduces the effective dimensionality havea length scale smaller than or of the order of the heal-ing length but larger than the scattering length [22–24],i.e., the underlying scattering process must remain threedimensional.Three important conditions must be met if our can-cellation technique is to apply. First, we cannot includethe low frequency harmonic trap often, but not always,present in experiments. Such a trap is used to keep atomsfrom spilling off the edge of a finite lattice. We requirethe potential to be sinusoidal, although we do allow thelattice constant to be different in each direction, leadingto, e.g., a rectangular lattice in two dimensions. Sec-ond, for condensates with three or more components, themean-field theory normally includes coherent couplingsbetween different components of the vector order param-eter; for two components, such couplings are preventedby angular momentum selection rules [25, 26]. We requireincoherent couplings only, which is appropriate for s = 2.For s > F . Finally, we can-not treat a mixture of scalar BEC’s of different masses,despite the fact that such a system can be described by avector mean-field theory with incoherent couplings only.Although certain conditions must be satisfied for it tobe applicable, our cancellation technique yields a panoplyof exact solutions of great physical significance. For ex-ample, for a two-component condensate on a rectangularoptical lattice, we find temporally periodic solutions inwhich the optical lattice is divided into two sublattices,and the condensate components oscillate back and forthbetween these sublattices. For the square optical lattice,we find a vortex-anti-vortex array for a scalar conden-sate, while for two-component condensates we obtain ex- otic solutions in which the optical lattice is divided into atotal of four sublattices, and the condensate componentsmove cyclically between these sublattices. As the dimen-sion D and the number of components s are increased,the number of solutions our technique generates growsrapidly. The number of solution types is so vast in threedimensions (3D) that, for the sake of brevity, we limitour discussion to stationary solutions with a high degreeof symmetry and to two examples of non-stationary so-lutions.The article is organized as follows. In Sec. II, we intro-duce the mean-field equations of motion and the opticalpotentials we will study. In Sec. III, a complete and rigor-ous treatment of exact dynamical solutions for arbitrarydimensions D and number of components s is presented.In Sec. IV, we treat select cases in detail, giving exam-ples of how to apply the results of Sec. III; of particularinterest are the vortex-anti-vortex array we find in 2Dand the array of intersecting vortex lines we find in 3D.In Sec. V, we present detailed stability studies of im-portant solution classes, including the vortex-anti-vortexarray, and an explicit connection to experimental units.Finally, in Sec. VI, we conclude. II. THE MEAN-FIELD EQUATIONS OFMOTION
Consider an s -component BEC in D dimensions withincoherent couplings between components. The conden-sate is subject to an optical potential formed by D lin-early polarized retroflected light waves. Let k l be thewave vector of the l th light wave, where l ∈ { , . . . , D } will be called the directional index. We will restrict ourattention to optical lattices with k l · k l ′ = 0 for l = l ′ .In particular, we will study BECs in one-dimensional,square, rectangular, and cubic optical lattices. For areview of optical potentials for neutral atoms, see, forinstance, Ref. [28].The mean-field equations of motion are i ¯ h ∂ψ j ∂t = − ¯ h m ∇ + s X j ′ =1 g jj ′ | ψ j ′ | + V j ψ j , (1)where ψ j ( r , t ) is the j th component of the conden-sate order parameter, j ∈ { , . . . , s } , and the posi-tion r is a vector in D dimensions. The j th compo-nent of the order parameter may be written ψ j ( r , t ) = p n j ( r , t ) exp[ iS j ( r , t )], where n j ( r , t ) is the number den-sity of the j th component at position r and v j ( r , t ) =(¯ h/m ) ∇ S j ( r , t ) is its velocity at that point [15]. Atomsin the j th component are subject to the optical potential V j . The coefficients g jj ′ of the nonlinear terms describethe binary interaction of an atom in component j andan atom in component j ′ : explicitly, g jj ′ = 4 π ¯ h a jj ′ /m ,where a jj ′ is the s -wave scattering length and m is theatomic mass, which, as stated in Sec. I, is assumed tobe independent of the component index j . Note that thenonlinear coefficient g jj ′ is renormalized by transverseconfinement as briefly alluded to in Sec. I for D = 1 and D = 2; see the references given there for more details.The n × n real, symmetric matrix M ≡ { g jj ′ } will be re-ferred to as the interaction matrix . We will assume thatall of the diagonal elements of M are nonzero, as is thecase in experiments.The atoms of the j th component are subject to theoptical potential V j ( r ) = − p j D X l =1 e l cos ( k l · r ) , (2)where p j is the atomic polarizability of an atom in the j th hyperfine state and e l is the electric field amplitude of the l th standing light wave. The wave vector k l determinesthe lattice constant in the l th direction. For convenience,we let V jl ≡ p j e l . Equation (2) then becomes V j ( r ) = − D X l =1 V jl cos ( k l · r ) . (3)We assume that all of the p j ’s and e l ’s are nonzero, sothat none of the V jl ’s vanish; if V jl were zero, the atomsof the j th component would be subject to an optical po-tential independent of the l th spatial coordinate. Notethat the atomic polarizability p j generally depends onthe component, or hyperfine state, j .To illustrate how one arrives at the potential given byEq. (2), consider the case D = 2. The total electric field E = E + E , where E l = e l cos( k l · r + χ l ) cos( c | k l | t + φ l ) . (4)The vector e l has constant, real components and is or-thogonal to k l , and c is the speed of light in vacuum. Byshifting the location of the origin if necessary, one canarrange for both of the spatial phases χ l to vanish. Thetemporally averaged intensity is then I = 12 e cos ( k · r ) + 12 e cos ( k · r ) + δ k ,k × e · e cos( φ − φ ) cos( k · r ) cos( k · r ) . (5)The j th component of the BEC is subject to the externaloptical potential V j = − p j I . The optical potentials V j take the form (2) with D = 2 if and only if the termin I coming from the interference of the two light wavesvanishes. If k and k differ, the interference term is zeroand the optical potential has rectangular symmetry. If k = k , on the other hand, the interference term van-ishes if either e · e = 0 or cos( φ − φ ) = 0. The opticallattice is then simply a square lattice. Finally, if k = k and 0 < | e · e cos( φ − φ ) | < e e , the structure of theoptical lattice is more complex [29].The interference terms in the intensity can also bemade to vanish for D = 3: for example, we can choose the vectors e j to be orthogonal to one another. If k = k = k , the optical lattice has simple cubic sym-metry.For simplicity, in this paper we will confine ourselves tooptical potentials in which the interference terms vanish.The potential is then given by Eq. (2). It is worth noting,though, that our solution techniques can be generalizedto optical potentials with nonzero interference terms. III. METHODS OF CONSTRUCTING EXACTSOLUTIONS TO THE MEAN-FIELDEQUATIONS OF MOTION
The time evolution of the order parameter is describedby the mean-field equations of motion (1) with the opticalpotentials (3). We will seek solutions to this problem inwhich each of the effective potentials U j ( r , t ) ≡ V j ( r ) + s X j ′ =1 g jj ′ | ψ j ′ ( r , t ) | (6)is constant. Solutions of this kind will be referred to as potential-canceling (PC) solutions and our solution tech-nique will be called the cancellation method because foreach j ∈ { , , . . . , s } , the spatial variation of the op-tical potential V j ( r ) is canceled by the variation of theterm P sj ′ =1 g jj ′ | ψ j ′ ( r , t ) | , rendering the effective poten-tial U j constant. PC solutions reduce the coupled non-linear mean-field equations of motion (1) to uncoupledlinear Schr¨odinger equations with constant potential: i ¯ h ∂ψ j ∂t = (cid:18) − ¯ h m ∇ + U j (cid:19) ψ j , (7)for j ∈ { , , . . . s } .For the U j ’s to be constant, the ψ j ’s must be linearcombinations of terms that vary sinusoidally with posi-tion. To be precise, we seek solutions of the form ψ j ( r , t ) = e − i Ω j t D X l =0 A jl cos( k l · r ) e − iω l t , (8)where ¯ hω l = ¯ h k l / m is the energy of a free particle ofmass m with wave number k l in the absence of nonlin-earity. Note that the possible values of the directionalindex l have been extended to include l = 0 in orderto simplify the notation: k ≡ l = 0 termgives rise to a constant offset in the order parameter. Thecoefficients A jl in our solution ansatz (8) are in generalcomplex, while the frequencies Ω j are real. The A jl ’s areconstrained by the requirement that the effective poten-tials U , U , . . . , U s are constant. These constraints willbe discussed in detail below.Substituting Eq. (8) into Eq. (1), we see that V j = − s X j ′ =1 g jj ′ | ψ j ′ | + ¯ h Ω j , (9)for j ∈ { , , . . . s } . This means that the effective poten-tial U j takes on the constant value ¯ h Ω j . Inserting Eq. (8)into Eq. (9) and comparing the resulting expression for V j with Eq. (3), we obtain D X l =1 V jl cos ( k l · r ) = D X l =0 D X l ′ =0 s X j ′ =1 g jj ′ A j ′ l A ∗ j ′ l ′ × cos( k l · r ) cos( k l ′ · r ) e − i ( ω l − ω l ′ ) t − ¯ h Ω j , (10)for j ∈ { , , . . . , s } . Equation (10) yields a set of alge-braic equations that the coefficients A jl must satisfy. Thefirst set of equations ensure that the cross terms on theright hand side of Eq. (10) vanish. Specifically, for eachpair of integers ( l, l ′ ) with 0 ≤ l < l ′ ≤ D , one obtains aset of conditions. If ω l = ω l ′ , we must have s X j ′ =1 g jj ′ A j ′ l A ∗ j ′ l ′ = 0 (11)for j ∈ { , . . . , s } . On the other hand, if ω l = ω l ′ , it issufficient to impose the weaker conditions ℜ s X j ′ =1 g jj ′ A j ′ l A ∗ j ′ l ′ = 0 . (12)Equating the coefficients of the terms that are propor-tional to cos ( k l · r ) on either side of Eq. (10), we seethat V jl = s X j ′ =1 g jj ′ | A j ′ l | (13)for j ∈ { , , . . . , s } and l ∈ { , , . . . , D } . (Note that thisequation does not apply for l = 0.) Finally, the constantterm on the right hand side of Eq. (10) must vanish, andso ¯ h Ω j = s X j ′ =1 g jj ′ | A j ′ | (14)for j ∈ { , , . . . , s } .The task of finding solutions to the coupled nonlinearpartial differential equations (1) has now been reducedto solving a system of algebraic equations: Eq. (13) mustbe solved for the coefficients A jl subject to the condi-tions (11) or (12) for each pair ( l, l ′ ) with l < l ′ . Asolution to these equations does not necessarily exist. Ifa solution does exist, Eq. (14) yields the frequencies Ω j .The solution, if it exists, is not uniquely specified bythe system of algebraic equations. To see this, let A l ≡ ( A l , . . . , A sl ) T (15)for l ∈ { , , , . . . , D } . Equation (13) determines thenorm of the vector A l for l ∈ { , , . . . , D } but does notconstrain the magnitude of A . In addition, Eqs. (11) and (12) with l = 0 place constraints on the direction of A but not its norm. The quantity | A | is therefore afree parameter.The length of the vector A is determined if the spa-tial average of the total density is given, as we will nowestablish. At time t , the number density of the j th com-ponent at position r is n j ( r , t ) ≡ | ψ j ( r , t ) | . Let h f i denote the spatial average of an arbitrary function f ( r ).Using Eq. (8), we find that the spatial average of thetotal number density h n i ≡ s X j =1 h n j i (16)is given by h n i = | A | + 12 D X l =1 | A l | . (17)Equation (17) has a solution for | A | if and only if2 h n i ≥ D X l =1 | A l | ; (18)if this condition is met, | A | is uniquely specified.Equation (11) states that the vector( A l A ∗ l ′ , . . . , A nl A ∗ nl ′ ) T is in the kernel of M , whileif Eq. (12) applies, the real part of this vector must bein the kernel of M . For this reason, most (but not all)of the solutions we obtain will be for cases in which theatomic interactions are such that det M = 0.We will now consider three particularly interesting andphysically important special cases in which the algebraicconditions that must be solved to yield a solution sim-plify dramatically. These special cases will be referredto as Special Cases A, B and C. We will also provide anexample that shows that the formalism just developedyields solutions to the mean-field equations of motioneven when the interaction matrix is nonsingular. Thisexample appears in Subsection III D. A. Factorizable Equations of Motion
A particularly simple special case is obtained when therank of M is unity. This is true to an excellent approxi-mation for the two-component condensates first producedby the JILA group that consist of two different hyperfinespin states of Rb: g , g and g are known to the1% level, and are in the proportion 1.03 : 1 : 0.97 [30].As a result, det M/ Tr M is zero to within experimentalerror.Because M has rank 1 and is a symmetric matrix, thereare nonzero, dimensionless, real numbers λ j and a σ = ± g jj ′ = σgλ j λ j ′ (19)for all j and j ′ . The quantity g is a positive constant withdimensions of energy times volume, which is inserted inEq. (19) to render the λ j ’s dimensionless. The magnitudeof g is arbitrary but fixed and, if desired, may be takento be the typical magnitude of the interaction coefficients g jj ′ .Let L ≡ ( λ , λ , . . . , λ s ) T . Equation (19) may then bewritten M = σg LL T , (20)showing that when the rank of the interaction matrix is1, M can be factored.Let Λ be the s × s matrix with elementsΛ jj ′ ≡ λ j δ j,j ′ . (21)For each pair of directional indices ( l, l ′ ) with 0 ≤ l 1. Constructing Exact Non-stationary Solutions Using aTransformation If we have a PC solution ψ to the equation of mo-tion (35), under certain circumstances we can constructnew solutions by transforming ψ . Let P be an invertible s × s matrix, and suppose that ψ is given by Eq. (8) andsatisfies Eq. (35). We set ζ ( r , t ) = T ( P ) ψ ( r , t ) , (36)where T ( P ) ≡ exp( − i ΩΛ t ) P exp( i ΩΛ t ) . (37)From Eq. (8), it follows that if P † Λ P = Λ , (38)then ζ ( r , t ) is also a solution to the equation of mo-tion (35). An important aspect of this transformationis that even if ψ ( r , t ) is stationary, ζ ( r , t ) can turn out tobe nonstationary. Thus, by transforming a single solu-tion ψ , we obtain a set of stationary and nonstationarysolutions. We will call each such set a P-set .For a given Λ, let S (Λ) be the set of invertible s × s matrices P that satisfy Eq. (38). It is straightforwardto show that the set of matrices T ( P ) with P ∈ S (Λ)forms a group. We will call this group by G (Λ). Laterin the paper we will study two examples in which G (Λ)is a continuous group, i.e., it has an uncountably infinitenumber of elements. As a result, the P -set is uncountablyinfinite in these examples. B. Factorizable Equations of Motion with EqualAtomic Polarizabilities A particularly important factorizable problem has λ j = 1 for all j , so that Λ is the identity matrix I . In thiscase, which we will call Special Case B , all of the inter-action strengths g jj ′ have the value σg , and the atomicpolarizabilities p j are all equal.The equation of motion (35) with Λ = I and V = 0was first studied by Manakov [31] and is now known asthe Manakov equation. We will extend this terminologyby also calling Eq. (35) with Λ = I and nonzero potential V the Manakov equation.The Manakov Case, i.e., Special Case B, is of consid-erable physical interest. Provided that the atoms are nottoo close to resonance, the p j ’s are to a good approxi-mation equal [32]. The interaction strengths are nearlyequal in two-component Rb condensates [30, 33]. As aresult, the dynamics of these condensates are reasonablywell described by the Manakov equation with s = 2.Three-component Na condensates with hyperfinespin F = 1 were first studied by the MIT group [34, 35].For F = 1 spinor condensates, the interaction strengths g jj ′ are identical and there are no incoherent couplingsif l and l are equal, where l F is the s -wave scatteringlength for two colliding atoms with total hyperfine spin F [25, 26]. Since the difference l − l is small comparedto l for Na [35, 36], it is a reasonable approximationuse the Manakov equation with s = 3 to model the three-component condensates produced by the MIT group, atleast for the initial stage of the time evolution.For the Manakov case, Eq. (28) reduces to | A l | = σp g e l , (39)where l ∈ { , , . . . , D } . This shows that σp must bepositive for there to be a solution to Eq. (39). Thus, forthe remainder of the paper, when we discuss Special CaseB, we will assume that σp > 0. For convenience, let a l ≡ s | p | g e l (40)for l = 1 , , . . . , D . Equation (39) is then simply | A l | = a l . The optical potentials V j are all equal to V for SpecialCase B. Equation (31) shows that V = − σgn + ¯ h Ω , (41)where n ≡ | ψ | is the total condensate number density.The total density is independent of time and varies si-nusoidally with position. The maxima of n are locatedat the potential minima for σ = +1. In contrast, for σ = − 1, the maxima of n are located at the maximaof the potential. This leads to an obvious instability, aspointed out by Bronski et al. [9] for the single-componentcase. Accordingly, for the remainder of the paper, we willlimit our attention to the case σ = +1 whenever we studya Case B problem [37]. Since we have already assumedthat σp > 0, this means that for Case B the atomic po-larizability p will be taken to be positive throughout theremainder of the paper.The condition (38) is particularly simple for the Man-akov equation: P can be any unitary matrix. Since Λ isthe identity matrix, T ( P ) = P and ζ is a unitary trans-formation of ψ . Unitary transformations of solutions tothe Manakov equation with an external potential havebeen studied elsewhere in one spatial dimension [14, 38].The transformation given by Eqs. (36)–(38) generalizesthat work to problems in which the λ j ’s are not all iden-tical, as well as to higher spatial dimensions. C. Factorizable Equations of Motion forTwo-Component Condensates with p = − p Consider a two-component BEC with factorizableequations of motion. Recall that λ and λ have unitmodulus and are real. As a result, there are four pos-sibilities: (i) λ = λ = 1, (ii) λ = λ = − 1, (iii) λ = − λ = 1 and (iv) − λ = λ = 1. Case (i) has al-ready been discussed: it is the Manakov case with s = 2.The equations of motion for case (ii) are unchanged ifwe reverse the signs of λ , λ and p , and so case (ii) isidentical to case (i). In precisely the same way, case (iv)is equivalent to case (iii). In this section, we will studycase (iii).The case in which a two-component BEC with factor-izable equations of motion has λ = − λ = 1 will bereferred to as Special Case C or the factorizable-with-opposite-polarizabilities (FOP) case. In this case, theatomic polarizabilities p and p have opposite signs andthe interaction matrix M = σg (cid:18) − − (cid:19) . (42)If σ is positive, atoms in the same condensate compo-nent repel each other and atoms in different condensatecomponents attract. The situation is reversed if σ is neg-ative. Finally, note that by switching the labels of thetwo components if necessary, we can arrange for σp to benegative. We will always take σp to be negative when wediscuss Special Case C.For Special Case C, the invertible matrix P (∆) ≡ (cid:18) cosh ∆ − sinh ∆ − sinh ∆ cosh ∆ (cid:19) , (43)satisfies the condition (38) for arbitrary real ∆. Equa-tions (36) and (37) with P = P (∆) and real ∆ thereforedefines a P -set. D. Two-Component Condensates with aNon-singular Interaction Matrix In the special cases discussed so far, the interactionmatrix M was singular. The cancellation method, how-ever, does yield solutions even if det M is nonzero. Toillustrate this point, let us consider a two-component con-densate in one dimension with det M = 0.It follows from Eq. (11) that A A ∗ = A A ∗ = 0 . (44) A and A cannot both vanish because the potentialcoefficients V j are nonzero. If both A and A arenonzero, A = A = 0, Ω = Ω = 0, and ψ j = A j cos( k x ) e − iω t (45)for j = 1 , 2. A solution of this form exists only if theequations V j = X j ′ =1 g jj ′ | A j ′ | (46)have a solution for A and A . The equations (46) havea solution if and only if χ ≡ ( g V − g V ) / det M > χ ≡ ( g V − g V ) / det M > . (48)If the inequalities (47) and (48) hold, A j = √ χ j yieldsa solution.We next turn to the case in which only one of the A j ’s is nonzero. The equations (46) have a solution with A = 0 if p g = p g > . (49)On the other hand, Eqs. (46) have a solution with A =0 if p g = p g > . (50)If the condition (50) is satisfied, we can switch the label-ing of the two condensate components, yielding Eq. (49).It is therefore sufficient to consider the case in which the condition (49) holds. In this case, A = p V /g and A = A = 0. It follows that ψ = A e − i Ω t and ψ = A e − i ( ω +Ω ) t cos( k x ). Equation (14) becomes¯ h Ω j = g j | A | , and, without loss of generality, we maytake A to be real. We conclude that the two-componentorder parameter is given by ψ = a e − ig a t/ ¯ h (51)and ψ = s V g e − i ( ω + g a / ¯ h ) t cos k x , (52)where a ≡ A is an arbitrary real constant. IV. APPLICATION OF ANALYTICALTECHNIQUES TO SELECT CASES We will now construct solutions to the mean-field equa-tions of motion (1) using the exact analytical methodsdeveloped in the preceding section. We will start withthe simplest cases as an introduction to the applicationof our solution methods, and to make connections withthe relatively simple solutions to be found in the litera-ture. We will then move on to progressively more richand complex problems with higher dimensions and/ormore condensate components than have previously beenconsidered. A. Solutions on a One Dimensional Optical Lattice It is convenient to orient the x axis along k , so that k = k ˆ x . Since ω = ω , Eq. (11) applies with l = 0and l ′ = 1. 1. One-Component Condensates We will begin with the simplest case, D = s = 1. For s = 1, Eq. (11) reduces to g A A ∗ = 0. It followsthat A and/or A must vanish. Equation (13) reducesto V = g | A | . Since V has been assumed to benonzero, A cannot vanish, and hence A = 0. Equa-tion (14) then shows that Ω = 0. There is a solution ofthe form of (8) only if g and V have the same sign. Ifthis is the case, we have the solution ψ = s V g cos( k x ) e − iω t (53)previously found by Bronski et al. [8, 9]. Note that thereis a PC solution only if the spatially averaged condensatedensity h n i happens to be V / (2 g ). Although this is avery restrictive condition, this simple case is neverthelessuseful in the development of nonlinear band theory [39]. 2. Two-Component Condensates We only found a single stationary solution for a one-component condensate in 1D. For two-component con-densates of the Manakov and FOP types, we find a muchlarger parameter space of solutions, including nonstation-ary solutions.We briefly touched on solutions for two-componentcondensates in one dimension in Section III D. In ob-taining the solution given by Eqs. (51) and (52), we didnot use our assumption that det M = 0. Moreover, thecondition (49) holds for both case B and case C. As aresult, the solution is also valid for cases B and C. Inboth cases, the solution takes the form ψ = ψ (1) ∗ , where ψ (1) ∗ ≡ e − iσga Λ t/ ¯ h (cid:18) a a cos( k x ) e − iω t (cid:19) . (54)The solutions for s = 2 constructed to this point arestationary, and have previously been obtained by Decon-inck et al. [13]. Let us now consider the Manakov andFOP cases B and C. We will demonstrate that for thesetwo cases there are nonstationary solutions in the same P -set as Eq. (54), where a P -set is the set of solutionsconnected by a matrix transformation of the order pa-rameter as described in Sec. III A 1; in the Manakov case,the transformation is just a unitary transformation.For Case B, Eq. (22) becomes A · A ∗ = 0. By changingthe phase of ψ j if necessary, we can arrange for A j tobe real for j = 1 , 2. We can arrange for A to be real bychanging the zero of time if needed. It then follows that A is real as well, and so the vectors A and A have realcomponents. Recalling that | A | = a , we obtain A = a (cos θ, sin θ ) T and A = a ( − sin θ, cos θ ) T , where a and θ are arbitrary real constants. The correspondingorder parameter is ψ ( r , t ) = e − iga t/ ¯ h [ A + A cos( k x ) e − iω t ] (55)by Eqs. (8) and (30). Recasting this solution, we have ψ ( r , t ) = P ( θ ) ψ (1) ∗ ( r , t ), where P ( θ ) ≡ (cid:18) cos θ − sin θ sin θ cos θ (cid:19) (56)is a unitary matrix. The solution ψ ( r , t ), which has beenpreviously described by Bradley et al. [14], is therefore aunitary transformation of the stationary solution (54). Ifsin(2 θ ) is nonzero, it is a nonstationary solution becausethe condensate component densities n = | ψ | = a cos θ + a sin θ cos ( k x ) − a a sin(2 θ ) cos( k x ) cos( ω t ) (57)and n = | ψ | = a sin θ + a cos θ cos ( k x )+ a a sin(2 θ ) cos( k x ) cos( ω t ) (58) oscillate in time with period T = 2 π/ω . We concludethat by performing unitary transformations of the singlesolution ψ (1) ∗ , we generate an uncountably infinite P -setof stationary and nonstationary solutions.What is the physical meaning of the solution (55)? Theexternal potentials V and V coincide and are equal to V = − pe cos ( k x ), and so the potential minima occurat the points x = qπ/k , where q is any integer. We di-vide the lattice of potential minima into two sublattices:sublattice 1 with even q , and sublattice 2 with odd q .The total condensate density n = n + n does not varyin time, and its maxima occur at the points x = qπ/k ,where q is any integer. Suppose for the sake of speci-ficity that a is positive and that π/ < θ < π . At time t = 0, the maxima of n are on sublattice 1, while at time t = T / 2, the maxima of n are on sublattice 2. At time t = T , the maxima of n are again on sublattice 1. Themaxima of n also oscillate between sublattices 1 and 2,but the oscillations of n lag those of n by half a period,ensuring that n is time-independent.The spatial average of the total number density h n i is a + a . If h n i is given, there is a solution of theform (55) if and only if h n i ≥ a . In constrast to thesingle-component case, we obtain a solution not just fora single value of h n i , but for a whole rangle of h n i values.This concludes our discussion of the relation betweenthe solutions obtained using the general formalism ofSec. III and the existing literature for one dimension.Let us now turn to the novel FOP case, Special Case C.We can again arrange for the vectors A and A to bereal. Equations (22) and (28) have the solution A = a (cosh ∆ , − sinh ∆) T and A = a ( − sinh ∆ , cosh ∆) T valid for arbitrary real a and ∆. Since ¯ h Ω = σga , thecondensate component order parameters are ψ = e − iσga t/ ¯ h [ a cosh ∆ − a sinh ∆ cos( k x ) e − ω t (cid:3) (59)and ψ = e iσga t/ ¯ h [ − a sinh ∆+ a sinh ∆ cos( k x ) e − ω t (cid:3) . (60)In contrast to the solution (55) we constructed for theManakov case, the component densities n and n oscil-late in phase between the two sublattices, and the totaldensity n oscillates in time as well. Since the spatial aver-age of the total number density h n i = cosh(2∆)( a + a ),we obtain solutions provided that h n i ≥ a .The vector order parameter may be written ψ = exp( − i ΩΛ t ) P (∆) exp( i ΩΛ t ) ψ (1) ∗ , (61)where P (∆) is defined by Eq. (43). Once again, the non-stationary solution can be constructed by transformingthe stationary solution (54) and we have an uncountablyinfinite P -set. 3. Three-Component Condensates We will not carry out the analysis for all possible casesfor three components. However, we will touch on two newfeatures that arise when we go from two components tothree.For two-component condensates governed by the Man-akov equations of motion, we showed that the cancella-tion method only yields solutions in which the oscillationsof the two components between the sublattices are 180 ◦ out of phase. Three-component condensates have an ad-ditional degree of freedom, and this leads to solutionswith a wide range of relative phases.We will restrict our attention to Case B and to so-lutions with A = A = A > A j . By changing the zero of time if needed,we can arrange for A to be real and positive. Set A j = | A j | e iφ j for j = 1 , , φ = 0.The condition (22) gives A + A + A = 0. Recallthat | A | = A + A + A = a . Clearly, there are so-lutions in which both A and A are real. In solutionsof this type, two of the components oscillate in phasewith one another and the other component is 180 ◦ out ofphase.In addition to these solutions, there is a solution forany φ and φ satisfying the conditions0 < φ < π < φ < π (62)and 0 < φ − φ < π . (63)The condensate component densities are | ψ j | = A j + | A j | cos ( k x )+2 A j | A j | cos( k x ) cos( ω t − φ j ) . (64)Thus, a wide variety of phase relationships among thecondensate component densities are possible for s = 3.If s is increased still further, the range of possible phaserelationships grows rapidly.For three-component condensates, it is possible for theinteraction matrix M to be singular and to have rankgreater than one. Even though the equations of motionare not factorizable, the cancellation method developedat the outset of Section III can be applied to yield so-lutions for problems of this type. We will illustrate thiswith a one-dimensional example.Suppose the eigenvalues of M are Υ = 0, Υ = 0 andΥ = 0, and let the µ i ’s be the associated real eigenvec-tors. We also set V ≡ ( V , V , V ) T (65)and w ll ′ ≡ ( A l A ∗ l ′ , A l A ∗ l ′ , A l A ∗ l ′ ) T (66)for l, l ′ = 0 , 1. The conditions (11) and (13) are then M w = 0 (67) and M w = V . (68)Equation (68) has a solution only if V is a linear com-bination of µ and µ . Suppose this is indeed the case,so that V = ˜ v µ + ˜ v µ . The equation M ξ = V has the solutions ξ = (˜ v / Υ ) µ + (˜ v / Υ ) µ + ˜ ξ µ , (69)where ˜ ξ is an arbitrary real constant. We can set w = ( | A | , | A | , | A | ) T = ξ = ( ξ , ξ , ξ ) T (70)if ξ j ≥ j = 1 , , 3. Suppose that the ξ j ’s are in factall positive. Then A j = p ξ j exp( iφ j ), where φ j is real.Equation (67) shows that w = ( A A ∗ , A A ∗ , A A ∗ ) T = C µ = C ( µ , µ , µ ) T , (71)where C is an arbitrary nonzero complex constant, andso A j = Cµ j exp( iφ j ) / p ξ j for j = 1 , , 3. Since A j and A j are both proportional to exp( iφ j ), the order pa-rameter ψ j is proportional to exp( iφ j ) as well, and wemay set φ j = 0 for j = 1 , , j ’s can be readily obtained from Eq. (14). We con-clude that the order parameter for the j th condensatecomponent is ψ j = ξ − / j e − i Ω j t (cid:2) ξ j + Cµ j cos( k x ) e − iω t (cid:3) (72)for j = 1 , , 3. By changing the zero of time if necessary,we can arrange for C to be real and positive. The densityof the j th condensate component is then n j = ξ j + C µ j ξ j cos ( k x )+2 Cµ j cos( k x ) cos( ω t ) , (73)which shows that the oscillations of each pair of conden-sate components are either in phase or 180 ◦ out of phase.The phase and amplitude of the oscillations of the com-ponents between sublattices 1 and 2 are determined bythe nature of µ , the eigenvector of the interaction ma-trix M with eigenvalue zero. B. Solutions on a Square Optical Lattice For D = 2, the optical lattice is formed by two standinglight waves with orthogonal wave vectors. We take the x axis to lie along k and the y axis to lie along k . For k = k , the optical lattice has rectangular symmetry,while for k = k ≡ k , we obtain a square optical lattice.In this section, we will study solutions on the square op-tical lattice. Solutions on the rectangular optical latticewill be discussed in Section IV C.0 1. One-Component Condensates The A jl ’s must satisfy Eq. (11) for l = 0 and l ′ = 1 , l = 1 and l ′ = 2. For s = 1, these conditions reduce to A A ∗ = A A ∗ = 0and ℜ ( A A ∗ ) = 0, respectively. Equation (13) has asolution only if p and g have the same sign. We assumethat this is the case. Since | A l | = V l /g > l = 1 , 2, both A and Ω must vanish. We can arrangefor A to be real and positive. A is then imaginary.The condensate wave function is thus ψ = 12 r p g [ e cos( kx ) ± ie cos( ky )] e − iωt , (74)where ω ≡ ¯ hk / µ . This is a stationary solution on thesquare optical lattice with potential V ≡ V = − p [ e cos ( kx ) + e cos ( ky )] (75)and time-independent condensate density n = − V /g .A PC solution therefore exists only if the spatiallyaveraged condensate density h n i is precisely p ( e + e ) / (8 g ).Although the solution (74) has previously been con-structed by Hai et al. [12], its physical interpretationhas not yet been discussed. The solution is a vortex-antivortex lattice [see Fig. 1, Parts (a) and (b)], and sothe condensate is flowing even though its density is nottime-dependent. Each square of side λ/ ψ ( r , t ) is a solution, then so is thetime-reversed state ψ ∗ ( r , − t ). Using our method of solu-tion, we found both the solution with the upper sign inEq. (74) and its time-reversed version, the solution withthe lower sign. 2. Two-Component Condensates We have now finished establishing contacts betweenthe literature and the solutions obtained using the for-malism of Section III. To the best of our knowledge, thesolutions found from this point on are new.The analysis for two-component condensates on asquare optical lattice runs parallel to that given for twocomponents on a one-dimensional optical lattice, and soonly the final results will be given. If det M = 0, there isa solution of the form ψ j = ( | A j | cos kx ± i | A j | cos ky ) e − iωt (76) (f)(a) (b)(c) (d)(e) FIG. 1: (Color online) Solutions on a square optical lattice :(a) Gray scale plot of the square optical lattice potential V ( x, y ). Regions of low (high) potential are shown in black(white). (b) The current density for the one-component so-lution (74) with the upper sign. The direction (size) of thearrows indicates the direction (magnitude) of the current flow.This solution is a vortex-anti-vortex array. (c)–(f) Gray scaleplots of the density of the first component n ( x, y ) for thetwo-component solution to the Manakov case with a = a , ρ = 1 and θ = 3 π/ 4. The plots are for times t = T / t = 3 T / t = 5 T / t = 7 T / n are shown inblack (white). Each plot shows the region with − λ ≤ x ≤ λ and − λ ≤ y ≤ λ . for j = 1 , 2, provided that Eq. (13) has a solution for the A jl ’s. Each of the two condensate components moves ina vortex-antivortex lattice in this solution. If the condi-tion (49) holds, on the other hand, we have a solutionwith ψ = a e − ig a t/ ¯ h (77)1and ψ = 12 r p g ( e cos k x ± ie cos ky ) e − i ( ω + g a / ¯ h ) t , (78)where a = A is an arbitrary real constant. For bothCase B and Case C, this is a valid solution and ψ = ψ (2) ∗ ,where ψ (2) ∗ ≡ e − iσga Λ t/ ¯ h × (cid:18) a ( a cos kx ± ia cos ky ) e − iωt (cid:19) . (79)There are nonstationary solutions if the problem is fac-torizable. For the Manakov Case B, we have an uncount-ably infinite P -set: the solution ψ = P ( θ ) ψ (2) ∗ is validfor arbitrary real θ . The densities of the two componentsof the condensate are time-dependent in this solution ifsin(2 θ ) is nonzero. For example, n ( x, y, t ) = a cos θ + a sin θ [cos ( kx ) + ρ cos ( ky )] − a a sin(2 θ )[cos( kx ) cos( ωt )+ ρ cos( ky ) sin( ωt )] , (80)where ρ ≡ ± e /e . The solution with ρ = − e /e issimply the time-reversed version of the solution with ρ = e /e , and so we may restrict our attention to the case ρ > 0. The external potentials V and V coincide and areequal to V = − p ( e cos kx + e cos ky ). The potentialminima occur at the points ( x, y ) = λ ( q , q ), where q and q are integers and λ is the optical wavelength. Wedivide the lattice of potential minima into four squaresublattices with lattice spacing λ , as shown in Fig. 2.Although n and n are time-dependent, the total con-densate density n does not vary in time, and its maximaoccur at the potential minima. The time evolution of thedensity of the first component, as described by Eq. (80),is illustrated in Fig. 1(c)–(f). Let T = 2 π/ω be the periodand suppose for the sake of specificity that a is positiveand that π/ < θ < π . At time t = T / 8, the maximaof n are on sublattice 1 [Fig. 1 (c)]. One quarter periodlater, the maxima of n are on sublattice 2 [Fig. 1 (d)].They are on sublattice 3 at time t = 5 T / t = 7 T / n return to sublattice 1 at time t = 9 T / n also oscillate among the sublattices,but the oscillations of n lag those of n by half a period.For the FOP Case C, we obtain a P -set of solutions bytransforming ψ (2) ∗ : explicitly, ψ = e − i ΩΛ t P (∆) e i ΩΛ t ψ (2) ∗ is a solution for arbitrary real ∆. The external potentialsare V = − V = V . Suppose for the sake of specificitythat p is positive. The minima of V and the maxima of V then occur at the lattice of points ( x, y ) = λ ( q , q ),where q and q are integers and λ is the optical wave-length. We divide this lattice into the same four squaresublattices as we did for Case B. As time passes, the FIG. 2: (Color online) Division of the lattice of potentialminima into four square sublattices. The lattice of potentialminima are the vertices of the grid. The vertices in sublattices1, 2, 3 and 4 are indicated by circles, squares, diamonds andstars, respectively. The origin is at the center of the figure. maxima of n move periodically among the four sublat-tices, just as they do for Case B. In Case C, however, theoscillations of component 2 are in phase with those ofcomponent 1, and the total condensate density n variesin time as a result. This is analogous to what we foundfor Case C in one dimension.For the solutions just discussed, the spatial average ofthe total number density h n i = a + ( a + a ) for theManakov case, while h n i = cosh(2∆)[ a + ( a + a )]for the FOP case. In both cases, we obtain solutionsprovided that h n i ≥ ( a + a ). C. Solutions on Rectangular Optical Lattices We now turn to the case in which D = 2 and k = k ,i.e., to the rectangular optical lattice. No solutions of theform of Eq. (8) exist for a single-component condensateon a rectangular optical lattice. For a two-componentcondensate, PC solutions are obtained only for the Man-akov Case B, and so we will confine our attention to thatcase. Choosing l = 0 and l ′ = 1 , A is nonzero, A and A must be parallel.Equation (22) with l = 1 and l ′ = 2 then shows thateither A or A must vanish and, hence, Eq. (28) can-not be satisfied for both l = 1 and l = 2. It followsthat A = 0. We can take the components of A tobe real without loss of generality. Since | A | = a , wemay set A = a (cos θ, sin θ ) T . We can arrange for thecomponents of A to be real through a change in thezero of time. Because A · A = 0 and | A | = a , itfollows that A = ± a ( − sin θ, cos θ ) T . Equation (30)shows that Ω = 0, and hence we have the solution givenby ψ = a cos θ cos( k x ) e − iω t ∓ a sin θ cos( k y ) e − iω t (81)2and ψ = a sin θ cos( k x ) e − iω t ± a cos θ cos( k y ) e − iω t , (82)where the angle θ is arbitrary. Equations (81)and (82) define a P -set of solutions since ψ = P ( θ )( a cos( k x ) e − iω t , ± a cos( k y ) e − iω t ) T .Equations (81) and (82) give a nonstationary solu-tion with temporal period T = 2 π/ | ω − ω | for 0 <θ < π/ 2. The time evolution of this solution canbe understood as follows. The optical potential V = − p [ e cos ( k x )+ e cos ( k y )] has minima at the points( x, y ) = π ( q /k , q /k ), where q and q are integers. Di-vide the lattice of potential minima into two sublattices:sublattice A with even q + q and sublattice B with odd q + q . For the solution given by Eqs. (81) and (82) withthe lower signs and 0 < θ < π/ 2, the maxima of n areinitially on sublattice A , as illustrated in Fig. 3(b) for a = a and θ = π/ 4. Half a period later, maxima of n are on sublattice B [see Fig. 3(c)]. The maxima of n are on sublattice A once again at time t = T . Thesecond component oscillates between the two sublatticesin the same way, but its oscillations lag those of the firstcomponent by half a period.In the limit that k and k coincide, the period of os-cillation T tends to infinity and we obtain a stationarysolution on the square optical lattice. In this solution, themaxima of n reside on one sublattice and the maximaof n are on the other. D. Solutions on the Simple Cubic Optical Lattice For condensates with two or more components on threedimensional optical lattices, the set of solutions of theform (8) is prohibitively large. Therefore, we will makea number of simplifying assumptions and will limit our-selves to giving examples of solutions. The stationarysolutions we will discuss all have a high degree of sym-metry.For D = 3, the optical lattice is formed by three stand-ing waves with orthogonal wave vectors. The x l axis willbe taken to lie along k l for l = 1 , , 3. We will confine ourattention to the case in which the three standing waveshave the same wavelength λ , so that the lattice of po-tential minima is a simple cubic (SC) lattice with latticespacing λ/ 2. We will further simplify the problem byrestricting our attention to the Manakov Case B and byassuming that the e l ’s coincide. To simplify the notation,set k ≡ k = k = k , ω ≡ ¯ hk / m , e ≡ e = e = e , f l ( r ) ≡ cos( k l · r ) for l = 1 , , 3, and a ≡ p p/ge/ ψ = exp( − ig | A | t/ ¯ h ) × (cid:2) A + ( A f + A f + A f ) e − iωt (cid:3) (83)is a solution to the mean-field equations of motion if | A l | = a and A ∗ · A l = 0 (84) (d)(a) (b)(c) FIG. 3: (Color online) Two-component condensate on a rect-angular optical lattice with λ = 2 λ : (a) Gray scale plot ofthe 2D optical potential V ( x, y ). Regions of low (high) po-tential are shown in black (white). (b)–(c) Density of the firstcomponent n ( x, y, t ) at times t = 0 and T / 2, respectively; re-gions of high (low) n are shown in black (white). (d) Currentdensity of the first condensate component at time t = T / A to sublattice B . Thedirection (size) of the arrows indicate the direction (magni-tude) of the current flow. Each of the four plots shows theregion with − λ ≤ x ≤ λ and − λ ≤ y ≤ λ . for l = 1 , , ℜ ( A ∗ · A ) = ℜ ( A ∗ · A ) = ℜ ( A ∗ · A ) = 0 . (85)There is no solution to Eqs. (84) and (85) for s = 1, andso we will only consider condensates with two or morecomponents.For brevity, the lattice of potential minima will be re-ferred to as “the lattice.” The lattice can be divided intoeight simple cubic sublattices with lattice spacing λ . Thevector f ≡ ( f , f , f ) (86)takes on a different value on each of these sublattices.The lattice can also be divided into four body-centeredcubic (BCC) sublattices. Each of these BCC sublatticesis the union of two simple cubic sublattices with f ’s thatsum to zero. In Table I, we assign labels to each of theeight simple cubic sublattices and to each of the four BCCsublattices. These sublattices are illustrated in Fig. 4 andwill play an important role in our examples.All of the stationary solutions we will discuss have atleast one of the two symmetries we will now define. If,3 TABLE I: Labeling of the sublattices of the simple cubic lat-tice f SC sublattice label BCC sublattice label(1,1,1) 0+ 0(-1,-1,-1) 0 − − − − z y x λ / 0+ 0 − − − 3+ 1+2+ 3 − FIG. 4: (Color online) Sublattices of the simple cubic lattice :The vertices of the grid are sites of the lattice of potentialminima. The gray sites belong to the SC sublattice 0+, whilethe white sites belong to the SC sublattice 0 − . The gray andwhite sites together make up the BCC sublattice 0. Sites ofall eight SC sublattices are labelled at the corners of the cubein which x , y and z all lie between 0 and λ/ 2. The sites ofSC sublattices 1+, 1 − , 2+, 2 − , 3+ and 3 − are colored red,cyan, green, magenta, blue and yellow, respectively. after a certain lattice translation, the density n j ( r ) is un-changed by a rotation of 90 ◦ about the x , y and z axesfor j = 1 , , . . . , s , then we say that a solution has four-fold rotational symmetry . Note that the lattice transla-tion could depend on the condensate component index j and could be the null translation. If a solution has four-fold rotational symmetry, the x , y and z directions areequivalent, and so this symmetry is a type of discretizedisotropy. On the other hand, if for each pair ( j, j ′ ) thereis a sequence of lattice translations or a series of rota-tions of 90 ◦ about the x , y or z axes that maps n j ( r )onto n j ′ ( r ), then we say that the solution has compo-nent symmetry . Intuitively speaking, the s componentsof the condensate all play the same role in a solution with component symmetry. 1. Two-Component Condensates For s = 2, all solutions of the form (83) are stationarybecause Eqs. (84) and (85) do not have a solution withnonzero A . Both examples of solutions we give willtherefore be stationary solutions.One possible solution with four-fold rotational symme-try is given by ψ = 1 √ a ( f + f + f ) e − iωt ≡ ψ ( rot )1 (87)and ψ = r a ( f + e πi/ f + e πi/ f ) e − iωt ≡ ψ ( rot )2 . (88)The maxima of the total density are located at the po-tential minima. A straightforward analysis reveals that n has its maxima on the BCC sublattice 0 and that themaxima of n reside on BCC sublattices 1, 2 and 3. Thissolution does not possess component symmetry.Component 1 is at rest since the phase of ψ is inde-pendent of position. In contrast, the flow of component2 is fascinating: it flows in a three-dimensional vortexlattice of great beauty, as we will now demonstrate.Consider the cube C in which x , y and z range between0 and λ/ 2. Each of the eight corners of the cube belongto a different simple cubic sublattice (see Fig. 4). Thesecond component of the condensate flows along six ofthe twelve edges of the cube. Specifically, component 2flows from the site 1+ to the site 3 − , and then to sites2+, 1 − , 3+ and 2 − before returning to site 1+. Thereis no mass current along the remaining six edges of thecube.The cyclic flow of component 2 suggests that there isa vortex line within the cube, and this is fact the case:a vortex line has its core along the cube diagonal thatjoins the 0+ site to the 0 − site. To establish this, we willbegin by considering the behavior of ψ close to the line x = y = z . Letˆ e ′ = r (cid:18) ˆ x − 12 ˆ y − 12 ˆ z (cid:19) , (89)ˆ e ′ = r 12 (ˆ y − ˆ z ) , (90)and ˆ e ′ = r 13 (ˆ x + ˆ y + ˆ z ) . (91)The vectors ˆ e ′ , ˆ e ′ and ˆ e ′ form an orthonormal triadwith ˆ e ′ = ˆ e ′ × ˆ e ′ . We introduce the new coordinates x ′ i = ˆ e ′ i · r , where i = 1 , , 3. On the line x = y = z ,4both x ′ and x ′ vanish and x ′ is arbitrary. Rewriting ψ in terms of the new coordinates and expanding for small x ′ and x ′ , we obtain ψ ∼ = − ka sin (cid:18) kx ′ √ (cid:19) ( x ′ + ix ′ ) e − iωt . (92)Equation (92) shows that there is a vortex line with itscore along the line x = y = z . The direction of the vectorˆ e ′ and the direction of the current flow around the vortexcore are related by the right hand rule; for brevity, wewill say that the vortex line is oriented along the vectorˆ e ′ .We have just shown that there is a vortex line withinthe cube C with its core along the cube diagonal, andthat the vortex line is oriented along the vector ˆ e ′ . Todetermine the nature of the flow throughout space, firstnote that ψ is invariant under the transformation x →− x , i.e., it is invariant under reflection about the y − z plane. Naturally, ψ is also invariant under the reflections y → − y and z → − z . These invariances give the flowwithin the cubical region in which x , y and z range from − λ/ λ/ 2. Since ψ is invariant under the latticetranslation r → r + ( q ˆ x + q ˆ y + q ˆ z ) λ for all integers q , q and q , the nature of the flow in the whole of spacecan now be inferred.The following picture emerges from this analysis.There is a vortex core along every line in the BCC sub-lattice 0 that joins an infinite chain of nearest-neighborsites. At these sites, the density of the second compo-nent is at a maximum, its current density is zero andthe potential is at a minimum. Each vortex core alsopasses thorough a chain of neighboring potential max-ima which alternate with the potential minima. At thepotential maxima, the density of the second componentof the condensate n is zero. In this way, the energy ofthe vortex array is minimized. Every vortex line is ori-ented along one of the following four vectors: ˆ x + ˆ y + ˆ z ,ˆ x − ˆ y − ˆ z , − ˆ x + ˆ y − ˆ z or − ˆ x − ˆ y + ˆ z .The solution with ψ = a √ f + e iπ/ f + e − iπ/ f ) e − iωt (93)and ψ = a √ f − e − iπ/ f − e iπ/ f ) e − iωt (94)has component symmetry but not four-fold rotationalsymmetry: The x direction is not equivalent to the y and z directions, and n and n differ only by a translationthrough the distance λ/ x axis. The maximaof n are on BCC sublattice 0, while the maxima of n are on BCC sublattice 1.A natural question to ask is whether, for a given s ,there is a solution with both four-fold rotational and component symmetries. The answer to this question is“no” for both s = 2 and 3, as we show in Appendix B. 2. Three-Component Condensates For s = 3, the stationary solution given by ψ = ψ ( rot )1 , ψ = √ φ ψ ( rot )2 and ψ = √ − φ ψ ( rot )2 with 0 < φ < ψ j = 23 a (cid:18) f + f + f − f j (cid:19) e − iωt , (95)for j = 1 , , 3. The maxima of n j are on the j th sublatticeand each component of the condensate is at rest. This so-lution does not have four-fold rotational symmetry sincethe x j direction is special for condensate component j .However, it does have component symmetry. To see this,consider an arbitrary pair of indices ( l, l ′ ) with l = l ′ .Let l ′′ be the integer belonging to the set { , , } thatdiffers from both l and l ′ . A 90 ◦ rotation about the x l ′′ axis interchanges f l and f l ′ , and so maps n l onto n l ′ .A nonstationary solution that illustrates just how com-plex the solutions for s = 3 can be is given by ψ = a √ f + f + if ) e − i ( ω +Ω) t , (96) ψ = a h √ − f + f + if ) e − iωt i e − i Ω t , (97)and ψ = a h √ − ( − f + f + if ) e − iωt i e − i Ω t , (98)where Ω = ga / ¯ h . The density of the first conden-sate component is time-independent and its maxima areon BCC sublattices 0 and 3. The density maxima ofcomponent 2 are on simple cubic sublattice 1+ at time t = T π tan − (1 / ≡ τ , on simple cubic sublattice 2+ at t = T / − τ , on simple cubic sublattice 1 − at t = T / τ ,and are on simple cubic sublattice 2 − at time t = T − τ .At time T + τ , the maxima of n have returned to simplecubic sublattice 1+. The motion of condensate compo-nent 3 is identical to that of the second component, ex-cept that the oscillations of n lag those of n by half aperiod. 3. Four-Component Condensates The solution space for four-component condensates isvery large. To see this, consider an arbitrary set of or-thonormal vectors with real components in four dimen-sions, { ˆ ǫ , ˆ ǫ , ˆ ǫ , ˆ ǫ } . Setting A = a ˆ ǫ and A l = a ˆ ǫ l for l = 1 , , 3, we obtain a solution to Eqs. (84) and (85)for arbitrary nonnegative real numbers a . One such so-lution is given by ψ = 12 (cid:2) − a + a ( f + f + f ) e − iωt (cid:3) e − iga t/ ¯ h (99)5and ψ j = 12 (cid:2) a + a ( f + f + f − f j − ) e − iωt (cid:3) × e − iga t/ ¯ h (100)for j = 2, 3, and 4.For two- and three-component condensates, no station-ary solution of the form (83) has both four-fold rotationaland component symmetries. Such a solution does existfor four-component condensates, however. The solutionis given by Eqs. (99) and (100) with a = 0. In this solu-tion, the maxima of n j are on sublattice j − n is unchanged by a rotation of 90 ◦ about the x , y and z axes. After a translation through λ ˆ x j − , the density n j becomes n and so is unchanged by a rotation of 90 ◦ about the x , y and z axes for j = 2, 3, 4.As we have seen, n j can be mapped onto n by a prim-itive lattice translation for j = 2, 3 and 4. On the otherhand, n is mapped onto n j ′ by a translation through λ ˆ x j ′ − for j ′ = 2, 3, and 4. It follows that, for eachpair ( j, j ′ ), there is a sequence of at most two primitivelattice translations that carries n j onto n j ′ , and so thesolution has component symmetry.For a > 0, Eqs. (99) and (100) describe a nonstation-ary solution. In this solution, the density maxima of thefirst component of the condensate are on the simple cubicsublattice 0 − at time t = 0 and are on the simple cubicsublattice 0+ at t = T / 2. For j = 2, 3 and 4, the maximaof n j are on the simple cubic sublattice j + at time t = 0and are on the simple cubic sublattice j − at t = T / t = T . V. NONLINEAR STABILITYA. Dimensionless mean-field equations Our numerical investigations of the stability of selectedsolutions to the mean-field equations (1) are performedwith a dimensionless form of the equations, using dimen-sionless position, time, and potential-energy variables, ξ = r /x , τ = t/t , and ˜ V j = V j /E , (101)defined in terms of units x , t , and E . For opticalpotentials of the form (3), we also define dimensionlesspotential-strength coefficients˜ V jl = V jl /E . (102)The length and time units, x and t , are related to theenergy unit E by x = ¯ h/ p mE and t = ¯ h/E . (103) The elements g jj ′ of the interaction matrix must alsobe put into dimensionless form, but first we note thatthey may be renormalized, depending on the dimen-sionality of the optical lattice. If there is narrow har-monic transverse confinement for one-dimensional andtwo-dimensional optical lattices, the appropriate formsin all dimensions are [40] g (3) jj ′ = 4 π ¯ h m a jj ′ , g (2) jj ′ = (cid:18) π ¯ h ω z m (cid:19) / a jj ′ , and g (1) jj ′ = 2¯ hω ⊥ a jj ′ , (104)where a jj ′ is the low-energy s -wave scattering length forspecies j and j ′ , the superscript on g jj ′ is D , the dimen-sionality of the optical lattice, and the confining poten-tials are characterized by the angular frequencies ω z intwo dimensions and ω ⊥ in one dimension.We may choose the elements ˜ g jj ′ of the dimensionlessform of the interaction matrix to be typically of orderone, so that the elements of the scattering-length matrixand the dimensional interaction matrix decompose as a jj ′ = a ˜ g jj ′ and g ( D ) jj ′ = g ( D ) ˜ g jj ′ , (105)where a and g ( D ) are scalar, dimensional factors. Thelatter is the same as the g appearing in Eq. (19), butwith the possible need for renormalization in lower di-mensions explicitly indicated by the superscript. Thedimensionless form of g ( D ) is˜ g ( D ) = g ( D ) E x D . (106)Its values, which follow from Eqs. (103)–(106), are˜ g (3) = 4 π ¯ h p mE a , ˜ g (2) = (cid:18) πmω z ¯ h (cid:19) / a , and ˜ g (1) = 2 r mE ω ⊥ a . (107)We absorb the square root of this dimensionless scalefactor into the order parameter, which then takes thedimensionless form˜ ψ j = p ˜ g ( D ) x D/ ψ j . (108)Thus, for a PC solution, the coefficients appearing in ˜ ψ j are related to those in Eq. (8) by˜ A jl = p ˜ g ( D ) x D/ A jl . (109)The normalization of the dimensionless order parameteris given by s X j =1 Z (cid:12)(cid:12) ˜ ψ j ( ξ , τ ) (cid:12)(cid:12) d ξ = N ˜ g ( D ) ∀ τ , (110)6wherein we see that ˜ g ( D ) plays a role equivalent to thenumber of particles, N . For a PC solution, the dimen-sionless mean number density is then h ˜ n i ≡ x D h n i = 1˜ g ( D ) s X j =1 (cid:18) | ˜ A j | + 12 D X l =1 | ˜ A jl | (cid:19) , (111)where h n i is given in Eq. (17).Finally, the mean-field equations (1) take the dimen-sionless form i ∂ ˜ ψ j ∂τ = (cid:20) − ∇ ξ + (cid:18) s X j ′ =1 ˜ g jj ′ | ˜ ψ j ′ | (cid:19) + ˜ V j (cid:21) ˜ ψ j . (112) B. Details of the calculations Our numerical stability tests use Eq. (112), propagat-ing a specified initial condition ˜ ψ j ( ξ , 0) forward in timevia a fifth-order Runge-Kutta algorithm with adaptivestep-size control [41], the spatial derivatives being calcu-lated in wave-vector space via a pseudo-spectral method.We perturb the solution by adding some white noise tothe initial condition before beginning the time propaga-tion. This is accomplished for each component of the or-der parameter by adding to the real and imaginary partsof each of its Fourier components a random number froma uniform distribution in the range ± . × − times themodulus of the largest Fourier component.We work within a spatial cell comprising four periodsof the optical lattice (two optical wavelengths) in eachof the D dimensions, applying periodic boundary condi-tions to that cell. The spatial grids contain n grid = 128points for one-dimensional cases and n grid = 32 points ineach dimension for two-dimensional cases. Thereby weare able to test the stability of solutions against pertur-bations having wavelengths ranging from 2 λ/n grid to 2 λ ,where λ is the optical wavelength.As a measure of the instability of a component of asolution at a particular instant of time, we use the vari-ance of its Fourier power spectrum relative to that at τ = 0 [39, 43], σ j ( τ ) = vuuuuut X κ (cid:2) ˜ f j ( κ , τ ) − ˜ f j ( κ , (cid:3) X κ (cid:2) ˜ f j ( κ , (cid:3) . (113)Here ˜ f j ( κ , τ ) ≡ (cid:12)(cid:12) ˜ φ j ( κ , τ ) (cid:12)(cid:12) , ˜ φ j ( κ , τ ) is the Fourier trans-form of ˜ ψ j ( ξ , τ ), and κ i ≡ π/ξ i .A solution is deemed to have reached the onset of in-stability when each of the σ j has exceeded 0 . vs. solutionparameters that the times represent literal lifetimes thatwould be observed in any particular experiment.However, as will become clear below, there is a fairlywell-defined boundary between unstable solutions andstable solutions, beyond which the instability-onset timesincrease extremely rapidly. The locations of those bound-aries are largely insensitive to details of the calculations.We therefore expect the parameter boundaries delimitingnumerically stable solutions to be experimentally mean-ingful, in the sense that within the stable regions ob-served lifetimes should be at least of order one second. C. Results of the calculations To make our results more concrete, we have chosentypical values for the laser wavelengths, λ = 800 nm inall directions, the trap frequencies, ω z = 2 π × 100 Hzand ω ⊥ = 2 π × 200 Hz, the atomic mass, m = 87 u, andthe s -wave scattering length, a = 55 ˚A. We will refer tothese below as the “system parameters.”The results are displayed primarily in recoil units, set-ting E = E R = ¯ h k L m = 2 π ¯ h mλ ≈ . µ K × k B (114)and t = t R = ¯ h/E R ≈ . × − s , (115)where k L is the laser wave number, and k B is Boltz-mann’s constant. The numerical values are obtained fromour chosen system parameters above.We present below selected numerical stability analysesfor one to three condensate components in one and twodimensions. Three-dimensional cases are not included,for they are too computationally demanding at this time.It is straightforward to transform the system-specificvalues shown in the figures below, the potential strengths,the instability-onset times, and the numbers of particlesper well, to values appropriate for alternative choices ofthe system parameters. We will elaborate on this pointin Sec. V D, following the presentation of the results.7 Unstable S t a b l e l og ( t i / t R ) Particles per well V /E R l og [ t i ( m s ) ] s = 1 40 100 200 300 400 500 600 7003 3 D = 2 245 0 FIG. 5: (Color online) Instability-onset times t i for a one-component PC solution on a square optical lattice. Thedimensionless potential-strength parameter V /E R and thecorresponding particle density are shown on the horizontalaxes. The solution is stable at both low and high potentialstrengths. 1. One component on a square lattice We choose the dimensionless interaction parameterand the dimensionless potential coefficients ˜ V jl to be˜ g = 1 . , ˜ V = V /E R , and ˜ V = V /E R , (116)where the potential-strength parameter V can be var-ied. Then the dimensionless solution corresponding toEq. (74) with the upper sign has coefficients˜ A = 0 , ˜ A = q ˜ V , and ˜ A = i q ˜ V . (117)Because the spatially constant term is required to vanish,the particle density, Eq. (17), is uniquely determined bythe potential-strength parameter.The instability-onset times for this solution are shownas a function of V in Fig. 5. Two time scales are in-cluded, showing both recoil times and milliseconds, witha maximum propagation time of several seconds. Thetop scale shows the particle density in particles per wellcorresponding to the potential-strength parameter shownon the bottom scale.While the solution becomes unstable in just a shorttime over much of the range of potential strengths shown,it is stable for sufficiently weak potentials. Much moresurprisingly, it is also stable for sufficiently strong po-tentials. To test whether the solution becomes unstableagain for potentials stronger than that at the boundarynear 3 . E R , we performed propagations to t ≈ t R ≈ 68 ms of solutions having V /E R as high as 48, finding norecurrence of instability. This extends well into the Mottinsulating regime, beyond the point of physical relevanceof the mean-field equations, as discussed in Sec. I. s = 2 V /E R P a r t i c l e s p e r w e ll D = 1 FIG. 6: (Color online) Regions of stability and instability fora PC solution having two components on a one-dimensionaloptical lattice. The points denote calculated boundaries be-tween stable and unstable behavior, the lines connecting themserving as guides to the eye. The inaccessible region has noPC solutions of the form (119). 2. Two components on a one-dimensional lattice Here we study the stability for the Manakov case, inwhich the dimensionless interaction matrix has rank oneand has all elements equal to one. The potential coeffi-cients we choose are˜ V = ˜ V = V /E R , (118)and the coefficients of the dimensionless solution corre-sponding to Eq. (54) are˜ A = α , ˜ A = 0 , ˜ A = 0 , and ˜ A = q ˜ V , (119)where α is a free parameter. The components of thissolution are then mixed using P ( θ ) of Eq. (56) with θ = π/ α to control the particle density independently of thepotential-strength parameter, giving a two-dimensionaldomain in which to investigate the stability of the so-lution. As is clear from Fig. 5, the boundaries of thestable regions are approximated well by the positions atwhich the instability-onset times have exceeded aboutthree hundred times t R . Consequently, we have used thatas a threshold to define those boundaries for the presentcase in scans over the potential-strength parameter atseveral fixed values of the particle density. The resolu-tion of the potential grid was 0 . E R , easily adequate forthe graphical delimitation of the stable regions. The re-sults are shown as a map of stable and unstable regionsin Fig. 6.8The boundary of the region marked “inaccessible” cor-responds to the vanishing of the coefficient ˜ A of the con-stant term in the solution. Within that region it is impos-sible to construct a PC solution of the chosen form, sinceEq. (34) cannot be satisfied. The region marked “stable”is that portion of the parameter space where the criterionfor stability is satisfied, and the region marked “unstable”corresponds to parameters for which the instability-onsettime falls below the threshold. As in the two-dimensionalcase shown in Fig. 5, the solution becomes stable in theweak-potential limit. However, in striking contrast to thetwo-dimensional case, there is no evidence of a second re-gion of stability at high potential strengths.In order to verify that apparent absence of stabil-ity for deep potentials, we performed calculations alongthe boundary of the inaccessible region with V /E R ashigh as 48, well beyond the highest shown in Fig. 6,finding only monotonically decreasing instability-onsettimes [44]. This confirms the observation in Fig. 6 thatthere is no additional stable region at high potentialstrength. 3. Two components on a square lattice As we did in one dimension, here we study the Man-akov case, with all elements of the dimensionless interac-tion matrix equal to one. Now there are four potentialcoefficients, which we choose to be equal:˜ V = ˜ V = ˜ V = ˜ V = V /E R . (120)The coefficients of the solution corresponding to Eq. (79)are ˜ A = α , ˜ A = 0 , ˜ A = ˜ A = 0 , ˜ A = q ˜ V , and˜ A = i q ˜ V , (121)where α allows us to set the particle density. Once againwe mix the components using P ( θ ) with θ = π/ α fixed atone and varying potential strengths are shown in Fig. 7,where the conventions are similar to those used for thesingle-component case in Fig. 5. As in that case, re-gions of stability occur at both low and high potentialstrengths, but now a rather striking additional region ofstability appears at intermediate potential strengths, justbelow 10 recoil energies.To explore this behavior in more detail, we map outregions of stability in the plane of particle density andpotential strength in Fig. 8 using the same strategy ap-plied in the one-dimensional case in Fig. 6. Three areasof stability are clearly evident, though the central one issomewhat narrower than the others. The fine black linerunning parallel to the boundary of the inaccessible re-gion is the track in the parameter-space plane followed l og [ t i ( m s ) ] Unstable UnstableParticles per well V /E R l og ( t i / t R ) S t a b l e S t a b l e S t a b l e s = 2 345 5 10 15 D = 2 FIG. 7: (Color online) Instability-onset times for a two-component PC solution on a square optical lattice. Thedimensionless potential-strength parameter V /E R and corre-sponding particle density are shown on the horizontal axes.The solution is stable in three regions, having low, high, andintermediate potential strengths. by the graph shown in Fig. 7. We extended the search forrenewed instability along this line to V /E R = 48, limit-ing the propagation time to t ≈ t R ≈ 68 ms, findingno evidence of further instability beyond the crossoverinto the stable region near V /E R = 16. Stable P a r t i c l e s p e r w e ll U n s t a b l e Inaccessible U n s t a b l e S t a b l e S t a b l e V /E R s = 2 D = 2 FIG. 8: (Color online) Regions of stability and instability forthe PC solution Eq. (121), having two components on a squareoptical lattice. The points denote calculated boundaries be-tween stable and unstable behavior, the lines connecting themserving as guides to the eye. The inaccessible region has noPC solutions of the form Eq. (121). The fine line parallel tothe boundary of the inaccessible region is the track followedby the graph of instability-onset times in Fig. 7. V /E R Unstable s = 3 Stable Inaccessible P a r t i c l e s p e r w e ll D = 1 FIG. 9: (Color online) Regions of stability and stability fora PC solution having three components on a one-dimensionaloptical lattice. The points denote calculated boundaries be-tween stable and unstable behavior, the lines connecting themserving as guides to the eye. The inaccessible region has noPC solutions of the form Eq. (124). 4. Three components on a one-dimensional lattice We have also tested the stability of a class of PC so-lutions having three components. The dimensionless in-teraction matrix then has nine elements, all ones in theManakov case, and rank equal to one. The potential co-efficients are all chosen to be the same:˜ V = ˜ V = ˜ V = V /E R . (122)The coefficients of the components of the dimensionlesssolution are all of equal magnitude, those of the spatiallyconstant term being real:˜ A = ˜ A = ˜ A = α , (123)with those of the space-dependent part chosen to havethe phases of the cube roots of unity:˜ A = s ˜ V , ˜ A = s ˜ V e i π/ , and˜ A = s ˜ V e i π/ . (124)The map of stable regions in the space of the free pa-rameters of this solution is shown in Fig. 9, which is vi-sually indistinguishable from its two-component analog,Fig. 6, having a region of stability where the potential issufficiently weak. In fact, the coordinates of all the pointsplotted on the graph are identical, within the resolutionof the scans. D. Alternative choices of system parameters Those aspects of the results presented above that aredependent on the system parameters chosen in Sec. V Care readily transformed to values corresponding to al-ternative choices of those parameters. Obviously, thepotential-strength parameter V is trivially obtained bymultiplying the dimensionless value V /E R by the recoilenergy corresponding to any desired set of system param-eters, and the instability-onset time t i is easily convertedby multiplying it by the ratio t ′ R /t R of the time units t ′ R corresponding to the alternative parameters and t R corresponding to our chosen system parameters.The average number of particles per well is just themean density times the well volume, N well = (cid:18) λ (cid:19) D h n i . (125)From the dimensionless density given in Eq. (111), wesee that this can be expressed in terms of dimensionlessparameters as N well = (cid:18) λ x R (cid:19) D g ( D ) s X j =1 (cid:18) | ˜ A j | + 12 D X l =1 | ˜ A jl | (cid:19) , (126)with the length unit x set to the recoil length x R =¯ h/ √ mE R .For all of our stability figures, the dimensionlesspotential-strength parameter V /E R is a free parameter,and it determines the values of the dimensionless coeffi-cients ˜ A jl having l > 0. For Figures 6, 8, and 9, there isone additional free parameter, α , and it determines thevalues of one or more of the coefficients ˜ A j . Thus, forany given abscissa in Figure 5 or 7, or abscissa and ordi-nate in Figure 6, 8, or 9, the sum in Eq. (126) is fixed,and the value of N well corresponding to an alternativechoice of system parameters can be obtained from thatshown on the graph by simply rescaling the prefactors: N ′ well = N well (cid:18) x R λ (cid:19) D ˜ g ( D ) (cid:18) λ ′ x ′ R (cid:19) D g ( D ) ′ , (127)where the unprimed quantities correspond to our choiceof system parameters, and the primed quantities to somealternative choice. VI. CONCLUSIONS In this paper, we made a comprehensive study ofpotential-canceling (PC) solutions for an s componentBose-Einstein condensate in in a D -dimensional opticallattice. Studies of specific cases with small s and D , es-pecially in one spatial dimension, have appeared in theliterature. Our work brings these previous studies to-gether, generalizes to arbitrary s and D , and provides0intriguing new solution types and novel physical inter-pretations.Currently, there is a great deal of interest in thea Berezinskii-Kosterlitz-Thouless phase in Bose-Einsteincondensates at intermediate temperatures in 2D [45]. Insuch a phase, vortex-anti-vortex pairs become bound to-gether, in contrast to the free vortex proliferation whichoccurs at high temperatures. However, this phase is re-stricted to a truly 2D system, which is difficult to achieveexperimentally. We have shown that an optical lat-tice stabilizes vortex-anti-vortex pairs in the quasi-two-dimensional case, and that the lattice causes the array tobe tightly packed.Not only have we presented multicomponent gener-alizations of 2D vortex-anti-vortex arrays, but we havealso generalized them to three dimensions. In 3D, weconstructed a solution in which one condensate compo-nent forms a lovely and complex three-dimensional ar-ray of intersecting vortex lines. As a part of our studyof PC solutions in 3D, we gave a thorough treatment ofthe most highly symmetric solutions for condensates withtwo, three and four components. Our formalism can alsobe used to gain insight into complex systems now exper-imentally available, such as five-component condensatesin 3D optical lattices.We studied the stability of PC solutions numericallyin 1D and 2D for one-, two and three-component con-densates. We found three main results: (1) potential-canceling solutions tend to become stable as the poten-tial strength is reduced; (2) there is a remarkable differ-ence between the one-dimensional and two-dimensionalsolutions, in that the latter are also stable for deep po-tentials; and (3) for two-components in a square opticallattice, there is a fascinating third region of stability forintermediate-strength potentials. We found no evidenceof stabilization at high potential strength in one dimen-sion.Finally, we mention that the possibility of experimen-tally realizing vortex-anti-vortex arrays in a 2D lattice isprovided for in the recent experiments of Sebby-Strabley et al. [46]. In those experiments, two polarizations ofthe lasers used to create the optical lattice potential aremanipulated to create lattices that can be dynamicallycontrolled on a site-by-site basis. In this way, one canimagine creating an array of small “propellers” to stirup vortex-anti-vortex pairs. The parameter ranges inwhich such a procedure would lead to stable structureswere determined in our numerical studies. To manipu-late multicomponent condensates, one can imagine moreadvanced versions of such an experiment, in which thefact that different hyperfine components “feel” differentlattice strengths for a given optical wavelength λ can beused to one’s advantage.We thank B. Deconinck, J. N. Kutz, and J. N. Robertsfor useful discussions. LDC’s work was supported by theNational Science Foundation under Grant PHY-0547845as part of the NSF CAREER program. APPENDIX A: PROOF THAT THE λ j ’S CAN BERESCALED TO HAVE UNIT MODULUS For Special Case A, the equations of motion are fac-torizable and are given by Eq. (35). The elements of theinteraction matrix are g jj ′ = σgλ j λ j ′ and the optical po-tentials are V j = λ j V . Consider an associated “normal-ized” problem that is also factorizable. In this normal-ized problem, the elements of the interaction matrix are˜ g jj ′ = σg ˜ λ j ˜ λ j ′ and the optical potentials are ˜ V j = ˜ λ j V ,where ˜ λ j ≡ λ j / | λ j | has unit modulus. Suppose we havea solution ˜ ψ j = e − i ˜Ω j t s X l =0 ˜ A jl cos( k l · r ) e − iω l t (A1)to the normalized problem. We can then construct a cor-responding solution to the original, unnormalized prob-lem as follows. We let A jl = ˜ A jl / p | λ j | and Ω j = | λ j | ˜Ω j ,and define ψ through Eq. (8). ψ is then a solution toEq. (35). Moreover, | ψ j | = | ˜ ψ j | / | λ j | . We see thatfor each solution ˜ ψ of the normalized problem, there isa corresponding solution ψ to the original, unnormalizedproblem, and that the density of the j th condensate com-ponent simply differs by the constant factor | λ j | − in thetwo problems. As a result, we may assume without lossof generality that the λ j ’s all have unit modulus.The length of ˜ A is a free parameter at this point.However, if the average total density h n i is specified inthe original, unnormalized problem, then Eq. (17) gives h n i = s X j =1 λ j | ˜ A j | + 12 D X l =1 | ˜ A jl | ! . (A2)If Eq. (A2) has a solution, it fixes the value of | ˜ A | . Weconclude that | ˜ A | is determined if h n i is given. APPENDIX B: PROOF THAT THERE ARE NOSOLUTIONS WITH BOTH FOUR-FOLDROTATIONAL AND COMPONENTSYMMETRIES FOR TWO- ANDTHREE-COMPONENT CONDENSATES IN 3D It was stated in Sec. IV D that there are no solutionswith both four-fold rotational and component symme-tries in 3D if the condensate has two or three compo-nents. Our proof is as follows. Let s be 2 or 3. Thecondensate order parameters are ψ j = X l =1 A jl f l ! e − iωt , (B1)where j ranges from 1 to s and f l ≡ cos( k l · r ). Thedensities are n j = X l =1 | A jl | f l + 2 X ≤ l 1. Similarly,the condition ℜ ( A ∗ · A ) = 0 implies that β = π + σ β ,where σ = ± 1. If σ = σ , the condition ℜ ( A ∗ · A ) = 0becomes cos( α − β ) = 0. Referring to Eq. (B3), weobserve that if n is to have four-fold rotational symme-try, cos α and cos β must vanish as well. This is notpossible. If σ = − σ , on the other hand, the condition ℜ ( A ∗ · A ) = 0 becomes cos α cos β = 0. Equation (B3)shows that if n is to have four-fold rotational symmetry,it is required that cos α = cos β = cos( α − β ) = 0,which is an impossibility. We conclude that there is nosolution with both four-fold rotational and componentsymmetries for two components. [1] A. J. Leggett, Rev. Mod. Phys. , 307 (2001).[2] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A.Sen, and U. Sen, Advances in Physics , 243 (2007).[3] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S.Hall, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. , 2498 (1999).[4] J. E. Williams and M. J. Holland, Nature , 568(1999).[5] J. M. Higbie, L. E. Sadler, S. Inouye, A. P. Chikkatur,S. R. Leslie, K. L. Moore, V. Savalli, and D. M. Stamper-Kurn, Phys. Rev. Lett. , 050401 (2005).[6] A. Widera, F. Gerbier, S. Folling, T. Gericke, O. Mandel,and I. Bloch, Phys. Rev. Lett. , 190405 (2005).[7] J. C. Bronski, L. D. Carr, B. Deconinck, and J. N. Kutz,Phys. Rev. Lett. , 1402 (2001).[8] J. C. Bronski, L. D. Carr, B. Deconinck, J. N. Kutz, andK. Promislow, Phys. Rev. E , 036612 (2001).[9] J. C. Bronski, L. D. Carr, R. Carretero-Gonz´alez, B. De-coninck, J. N. Kutz, and K. Promislow, Phys. Rev. E ,056615 (2001).[10] B. Deconinck, B. A. Frigyik, and J. N. Kutz, Phys. Lett.A , 177 (2001).[11] B. Deconinck, B. A. Frigyik, and J. N. Kutz, J. NonlinearSci. , 169 (2002).[12] W. Hai, C. Lee, X. Fang, and K. Gao, Physica A ,445 (2004).[13] B. Deconinck, J. N. Kutz, M. S. Patterson, and B. W.Warner, J. Phys. A: Math. Gen. , 5431 (2003).[14] R. M. Bradley, B. Deconinck, and J. N. Kutz, J. Phys.A: Math. Gen. , 1901 (2005).[15] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari,Rev. Mod. Phys. , 463 (1999).[16] J. Mur-Petit, M. Guilleumas, A. Polls, A. Sanpera, M.Lewenstein, K. Bongs, and K. Sengstock, Phys. Rev. A , 013629 (2006).[17] See Ref. [19] for a detailed derivation of t h and U fromfirst principles quantum field theory.[18] Note that t h is normally called t in the condensed matterliterature and is sometimes denoted J in the case of ul-tracold quantum gases. The J notation leads to a “ J − J ”model instead of a t − J model, and so we avoid it. Wewill reserve t for time as is standard in dynamics, anduse t h for the hopping energy.[19] A. M. Rey, Ph.D. thesis, University of Maryland, 2004.[20] R. V. Mishmash and L. D. Carr, Phys. Rev. Lett. , underreview; e-print http://arxiv.org/abs/0710.0045 (2007).[21] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[22] L. D. Carr, M. A. Leung, and W. P. Reinhardt, J. Phys.B: At. Mol. Opt. Phys. , 3983 (2000).[23] D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov,Phys. Rev. Lett. , 2551 (2000).[24] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven,Phys. Rev. Lett. , 3745 (2000).[25] T.-L. Ho, Phys. Rev. Lett. , 742 (1998).[26] T. Ohmi and K. Machida, J. Phys. Soc. Japan , 1822(1998).[27] C. K. Law, H. Pu, and N. P. Bigelow, Phys. Rev. Lett. , 5257 (1998).[28] R. Grimm, M. Weidemuller, and Y. B. Ovchinnikov, Adv.Atom. Mol. Opt. Phys. , 95 (2000).[29] A. Hemmerlich, D. Schropp, T. Esslinger, and T. W.H¨ansch, Europhys. Lett. , 391 (1992).[30] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman,and E. A. Cornell, Phys. Rev. Lett. , 1539 (1998).[31] S. V. Manakov, Sov. Phys. JETP , 693 (1974).[32] J. N. Roberts, private communication.[33] C. J. Myatt, E. A. Burt, R. W. Ghrist, E. A. Cornell,and C. E. Wieman, Phys. Rev. Lett. , 586 (1997). [34] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur,S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle,Phys. Rev. Lett. , 2027 (1998).[35] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H.-J. Mies-ner, A. P. Chikkatur, and W. Ketterle, Nature , 345(1998).[36] J. P. Burke, Jr., C. H. Greene, and J. L. Bohn, Phys.Rev. Lett. , 3355 (1998).[37] Note, however, that it is a simple matter to write downthe corresponding σ = − , 063605 (2004).[39] B. T. Seaman, L. D. Carr, and M. J. Holland, Phys. Rev.A , 033622 (2005).[40] L. D. Carr, M. J. Holland, and B. A. Malomed, J. Phys.B: At. Mol. Opt. , 3217 (2005).[41] We use the driver routine odeint and the stepper routine rkqs of [42].[42] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of ScientificComputing (Cambridge Univ. Press, Cambridge, U.K.,1993).[43] B. T. Seaman, L. D. Carr, and M. J. Holland, Phys. Rev.A , 033602 (2005).[44] For very weak potentials, with V /E R below about 0 . V /E R = 0 . 16. Evidently the struc-ture of the stability map is more complicated in the farlower-left corner of Fig. 6 than in the rest of the plane.[45] P. Kr¨uger, Z. Hadzibabic, and J. Dalibard, Phys. Rev.Lett. , 040402 (2007).[46] J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V.Porto, Phys. Rev. A73