PHYSICAL REVIEW A 82, 053601 (2010) Dimensional reduction of a binary Bose-Einstein condensate in mixed dimensions Luis E. Young-S.,1,* L. Salasnich,2,3,† and S. K. Adhikari1,‡ 1Instituto de Fı́sica Teórica, UNESP–Universidade Estadual Paulista, 01.140-070 São Paulo, São Paulo, Brazil 2INO-CNR, Research Unit at the Dipartimento di Fisica “Galileo Galilei”, Università di Padova, Via Marzolo 8, I-35131 Padova, Italy 3CAMTP, University of Maribor, Krekova 2, SLO-2000 Maribor, Slovenia (Received 30 August 2010; published 1 November 2010) We present effective reduced equations for the study of a binary Bose-Einstein condensate (BEC), where the confining potentials of the two BEC components have distinct asymmetry so that the components belong to different space dimensions as in a recent experiment [G. Lamporesi et al., Phys. Rev. Lett. 104, 153202 (2010).]. Starting from a binary three-dimensional (3D) Gross-Pitaevskii equation (GPE) and using a Lagrangian variational approach we derive a binary effective nonlinear Schrödinger equation with components in different reduced dimensions, for example, the first component in one dimension and the second in two dimensions as appropriate to represent a cigar-shaped BEC coupled to a disk-shaped BEC. We demonstrate that the effective reduced binary equation, which depends on the geometry of the system, is quite reliable when compared with the binary 3D GPE and can be efficiently used to perform numerical simulation and analytical calculation for the investigation of static and dynamic properties of a binary BEC in mixed dimensions. DOI: 10.1103/PhysRevA.82.053601 PACS number(s): 03.75.Mn, 03.75.Hh, 03.75.Kk I. INTRODUCTION It is now well established that, at ultralow temperature, dilute Bose-Einstein condensates (BECs) can be accurately described by the three-dimensional (3D) Gross-Pitaevskii equation (GPE) [1,2]. Some years ago it was found that, starting from the 3D GPE and using a Gaussian variational ap- proach [3,4], it is possible to derive effective one-dimensional (1D) or two-dimensional (2D) wave equations which describe the dynamics of quasi-1D or quasi-2D BECs as appropriate for cigar or disk shapes [5]. For instance, in that derivation of the effective 1D equation the trapping potential is harmonic and isotropic in the transverse direction and generic in the axial one. The effective 1D equation, which is a time-dependent nonpolynomial Schrödinger equation (NPSE) [5], has been used to model a cigar-shaped condensate by many experimen- tal and theoretical groups [6–10]. A similar reduction scheme also exists for a disk-shaped BEC leading to a 2D NPSE. (The lowest-order approximation of the NPSE leads to a GP-type nonlinear Schrödinger equation (NLSE) equation with cubic nonlinearity. Different variants of this reduction scheme have been suggested, which offer certain advantages in some cases [7,8].) Moreover, by using the 1D NPSE, analytical and numerical solutions have been found for solitons and vortices [11]. Recently, a discrete version of the 1D NPSE (1D DNPSE) has been obtained [12]. The 1D DNPSE has been used to model a BEC confined in a combination of a cigar-shaped trap and a deep optical lattice acting in the axial direction [12]. More recently the 2D NPSE, obtained for a disk-shaped BEC subject to strong confinement along the axial direction, has been used to predict the density profiles and dynamical stability of repulsive and attractive BECs with zero and finite topological charge in various planar trapping configurations, including the *lyoung@ift.unesp.br. †luca.salasnich@cnr.it. ‡adhikari44@yahoo.com; [http://www.ift.unesp.br/users/adhikari/]. axisymmetric harmonic confinement and 1D periodic potential [13]. The reduction scheme has also successfully been applied [8,14] to the study of a density-functional formulation [15] for a Fermi superfluid at unitarity as well as to the study of Anderson localization in a BEC [16]. In a very recent experiment Lamporesi et al. [17] investi- gated the mixed-dimensional scattering occurring when the collisional atoms of two species of a binary BEC live in different space dimensions. This achievement opens the way to the experimental and theoretical studies of a binary BEC with the components belonging to mixed dimensions. In the present article we consider the dimensional reduction of a binary 3D GPE appropriate for the study of a binary BEC in mixed dimensions. Starting from a binary 3D GPE we derive a binary nonpolynomial Schrödinger equation (BNPSE) in reduced dimensions by using a Lagrangian variational approach. (The BNPSEs represent a set of coupled nonlinear differential equations with nonpolynomial non-power-law nonlinearities.) We also present a simplified version of the reduction scheme first, where the reduced equations in the form of a binary NLSE (BNLSE) have simpler power-law nonlinearities. The BNPSEs depend on the geometry of the system and also on the interspecies interaction. In fact, the main difference with respect to the case of a single BEC is the coupling between the two BECs, which gives rise, quite generally, to n integral term in the reduced wave equations. (In the past, a binary 1D NPSE has been derived and used [18] to study a binary cigar-shaped BEC, rather than a binary BEC in mixed dimensions as considered here.) We use the BNLSE and BNPSE to study numerically the static and dynamic properties of a 1D-2D binary BEC where the first quasi-1D component lies along the z axis and the quasi- 2D component lies in the (x,y) plane. For small nonlinearities both the BNLSE and BNPSE yielded results for density in good agreement with the full 3D result. At larger values of nonlinearities the BNPSE provided better approximation to the full 3D result. We also studied small oscillations of the 1D-2D binary BEC initiated by a sudden change of the 1050-2947/2010/82(5)/053601(9) 053601-1 ©2010 The American Physical Society http://dx.doi.org/10.1103/PhysRevLett.104.153202 http://dx.doi.org/10.1103/PhysRevLett.104.153202 http://dx.doi.org/10.1103/PhysRevA.82.053601 LUIS E. YOUNG-S., L. SALASNICH, AND S. K. ADHIKARI PHYSICAL REVIEW A 82, 053601 (2010) radial angular frequency of the quasi-2D component. The BNPSE provided a better description of this oscillation over the BNLSE, when compared with the full 3D model, especially for large nonlinearities and large times. The reduced effective BNLSEs with simple power-law nonlinearity, but with different dimensionalities for the com- ponents, are presented in Sec. II. These equations are derived by considering the ansatz for the wave function that sets the transverse spatial dependence as the Gaussian ground states in transverse harmonic traps. The transverse spatial dependence is finally integrated out and the reduced equations are derived in the usual fashion using a Gaussian variational approach. In Sec. III the BNPSEs with different dimensionalities for the components are derived by taking the ansatz for the wave function that sets the transverse spatial dependence as generic Gaussian states which are not the ground states in transverse harmonic traps. The reduced effective BNPSEs are then derived by a Lagrangian variational approximation which determines the widths of the Gaussian functions in a self- consistent fashion. In Sec. IV we present the numerical results, where we compare the densities and chemical potentials of the 3D binary GPE with those of the BNLSE and the BNPSE. The binary BEC has a disk-shaped component along the z direction coupled to a cigar-shaped component in the (x,y) plane. For smaller nonlinearities of the 3D binary GPE, we establish good agreement between the results for density and chemical potential of the effective BNLSE and the full 3D binary GPE. For larger nonlinearities a better result was obtained with the BNPSE. We also studied small oscillations of the binary 1D-2D BEC initiated by a sudden change of the radial angular frequency of the quasi-2D disk-shaped BEC using the BNPSE, the BNLSE, and the full 3D equations. Again the results obtained with the BNPSE were in better agreement with those of the BNLSE. Finally, a summary of our investigation is presented in Sec V. II. THEORETICAL FORMULATION FOR THE BNLSE We shall now explicitly demonstrate the reduction pro- cedure in three different types of binary mixtures: 1D-3D, 1D-2D, and 2D-3D. We shall consider a binary BEC where the trap anisotropy of the two components are distinct. When the component 2 has a fully asymmetric trap, the component 1, due to its specific trap symmetry, will be assumed to have a quasi-1D or quasi-2D configuration. The quasi-1D shape emerges when the traps in the two transverse directions are much stronger than that in the linear direction. The quasi-2D shape emerges when the trap in the direction perpendicular to the 2D plane is much stronger than those in the 2D plane. In this fashion one could have effective 1D-3D and 2D-3D mixed configurations for the traps acting on the binary BEC described by effective 1D-3D and 2D-3D BNLSEs. There remains one distinct case where one of the components has a quasi-1D shape and the other component has a quasi-2D shape due to specific trap symmetry of the two components. In that case the binary mixture is described by an effective 1D-2D BNLSE. We discuss these possibilities in the following and illustrate the corresponding reduction schemes. The action of a system of two BECs is given by [19] A = ∫ (L1 + L2 − g12|ψ1|2|ψ2|2) d3rdt, (1) where the Lagrangian density is Lj = { i h̄ 2 (ψ∗ j ∂tψj − ψj∂tψ ∗ j ) − h̄2 2mj |∇rψj |2 −Uj (r)|ψj |2 − gj 2 |ψj |4 } , (2) with Uj (r) being the external potential acting on the j th BEC component (j = 1,2), which is described by the macroscopic wave function ψj (r,t). Here ∂t denotes the time derivative, gj = 4πh̄2aj/mj , g12 = 2πh̄2a12/mR , and ∫ |ψj (r,t)|2 = Nj , where aj is the intraspecies scattering length for species j and a12 is the interspecies scattering length, mj is the mass of species j , and mR ≡ m1m2/(m1 + m2) is the reduced mass. By extremizing the action A with respect to ψ∗ 1 (r,t) and ψ∗ 2 (r,t) we obtain the following binary GPE [19]: ih̄∂tψ1 = [ − h̄2 2m1 ∇2 r + U1(r) + g1|ψ1|2 + g12|ψ2|2 ] ψ1, (3) ih̄∂tψ2 = [ − h̄2 2m2 ∇2 r + U2(r) + g2|ψ2|2 + g12|ψ1|2 ] ψ2. (4) These two coupled partial differential equations are both in 3 + 1 (space-time) dimensions. In the general case, these can be solved numerically in a routine fashion [20,21]. However, in many situations of interest, when one or both of the BEC components have extreme distinct cigar and disk shapes and when we are interested in density and dynamics in axial and radial directions, respectively, the following reduction procedures are of particular use. A. 1D-3D binary BEC Let us suppose that the confining potential of the first BEC is U1(r) = 1 2m1 ( ω2 1xx 2 + ω2 1yy 2 ) + V1(z), (5) where V1(z) is a generic potential in the z direction, and ω1x and ω1y are the angular frequencies of the harmonic trap along the transverse x and y directions. The transverse harmonic potentials are assumed to be so strong that the first BEC is quasi-1D, that is, in an elongated cigar-shaped configuration. We can then use the variational wave function [5] ψ1(r,t) = e−x2/(2a2 1x ) (π1/2a1x)1/2 e−y2/(2a2 1y ) (π1/2a1y)1/2 f1(z,t), (6) for the first BEC, where a1x = √ h̄/(m1ω1x) and a1y =√ h̄/(m1ω1y) are the characteristic harmonic lengths of the os- cillators in the x and y directions, respectively. Consequently, the action of the system becomes A = ∫ L1dzdt + ∫ L2d 3rdt + ∫ L12d 3rdt, (7) 053601-2 DIMENSIONAL REDUCTION OF A BINARY BOSE- . . . PHYSICAL REVIEW A 82, 053601 (2010) where L1 = i h̄ 2 (f ∗ 1 ∂tf1 − f1∂tf ∗ 1 ) − h̄2 2m1 |∂zf1|2 −V1(z)|f1|2 − E1|f1|2 − g1 4πa1xa1y |f1|4, (8) L12 = − g12 πa1xa1y e−(x2/a2 1x+y2/a2 1y )|f1|2|ψ2|2, (9) E1 = h̄2 4m1a 2 1x + m1ω 2 1xa 2 1x 4 + h̄2 4m1a 2 1y + m1ω 2 1ya 2 1y 4 . (10) We now extremize the action A with respect to f ∗ 1 (z,t) and ψ∗ 2 (r,t) to obtain the effective BNLSE for f1(z,t) and ψ2(r,t): ih̄∂tf1 = [ − h̄2 2m1 ∂2 z + V1(z) + E1 + g1 2πa1xa1y |f1|2 + g12 πa1xa1y ∫ e−(x2/a2 1x+y2/a2 1y )|ψ2|2dxdy ] f1, (11) ih̄∂tψ2 = [ − h̄2 2m2 ∇2 r + U2(r) + g2|ψ2|2 + g12 πa1xa1y e−(x2/a2 1x+y2/a2 1y )|f1|2 ] ψ2. (12) Equations (11) and (12) are simplified if we impose cylindrical symmetry in the harmonic potential acting on the first BEC, that is, ω1x = ω1y = ω1ρ , and also in the generic potential acting on the second BEC, that is, U2(r) = U2(ρ,z) with ρ = (x2 + y2)1/2. In this way we can choose ψ2(r,t) = ψ2(ρ,z,t) and find the following effective BNLSE for f1(z,t) and ψ2(ρ,z,t): ih̄∂tf1 = [ − h̄2 2m1 ∂2 z + V1(z) + E1 + g1 2πa2 1ρ |f1|2 + 2g12 a2 1ρ ∫ ∞ 0 e−ρ2/a2 1ρ |ψ2|2ρdρ ] f1, (13) ih̄∂tψ2 = [ − h̄2 2m2 ∇2 ρz + U2(ρ,z) + g2|ψ2|2 + g12 πa2 1ρ e−ρ2/a2 1ρ |f1|2 ] ψ2, (14) where a1ρ ≡ a1x ≡ a1y = √ h̄/(m1ω1ρ). In the coupled set of Eqs. (13) and (14), the first is integro-differential in 1 + 1 dimensions and the second is differential in 3 + 1 dimensions (but in 2 + 1 variables). B. 1D-2D binary BEC 1. Parallel BECs Let us suppose that the confining potential of the first BEC is still given by Eq. (5), while the confining potential of the second BEC is U2(r) = 1 2m2ω 2 2xx 2 + V2(y,z), (15) where V2(y,z) is a generic potential acting in the (y,z) plane. In addition we impose that the harmonic potential of angular frequency ω2x along the x axis is so strong that the second BEC is quasi-2D, that is, a very flat two-dimensional configuration in the (y,z) plane. Consequently, we have binary 1D-2D BEC where the first BEC is cigar shaped along the z axis and the second disk shaped in the (y,z) plane. We can then use the variational wave function of Eq. (6) for the first BEC, while for the second BEC we adopt the following variational wave function [5] ψ2(r,t) = e−x2/(2a2 2x ) (π1/2a2x)1/2 φ2(y,z,t), (16) where a2x = √ h̄/(m2ω2x) is the characteristic harmonic length of the oscillator in the x direction acting on the second BEC. In this way the action of the system becomes A = ∫ L1dzdt + ∫ L̃2dydzdt + ∫ L̃12dydzdt, (17) where L1 is still given by Eq. (8), while L̃2 = i h̄ 2 (φ∗ 2∂tφ2 − φ2∂tφ ∗ 2 ) − h̄2 2m2 |∇yzφ2|2 − Ẽ2|φ2|2 −V2(y,z)|φ2|2 − g2 2(2π )1/2a2x |φ2|4, (18) Ẽ2 = h̄2 4m2a 2 2x + m2ω 2 2xa 2 2x 4 , (19) L̃12 = − g12 πa1y √ a2 1x + a2 2x |f1|2|φ2|2e−y2/a2 1y . (20) We now extremize the action A and find the following effective BNLSE for f1(z,t) and φ2(y,z,t): ih̄∂tf1 = [ − h̄2 2m1 ∂2 z + V1(z) + E1 + g1 2πa1xa1y |f1|2 + g12 πa1y √ a2 1x + a2 2x ∫ ∞ −∞ e−y2/a2 1y |φ2|2dy ] f1, (21) ih̄∂tφ2 = [ − h̄2 2m2 ∇2 yz + V2(y,z) + Ẽ2 + g2√ 2πa2x |φ2|2 + g12e −y2/a2 1y πa1y √ a2 1x + a2 2x |f1|2 ] φ2. (22) These are two coupled equations, where the first is integro- differential in 1 + 1 dimensions and the second is differential in 2 + 1 dimensions. Equations (21) and (22) describe the binary system where the axis of the cigar-shaped BEC along the z direction is in the (y,z) plane of the two-dimensional BEC. 2. Perpendicular BECs Another interesting possibility is a binary system where the cigar-shaped BEC is perpendicular to the plane of the two- dimensional BEC. This kind of configuration can be obtained by considering the following trapping potential acting on the second BEC: U2(r) = W2(x,y) + 1 2m2ω 2 2zz 2, (23) 053601-3 LUIS E. YOUNG-S., L. SALASNICH, AND S. K. ADHIKARI PHYSICAL REVIEW A 82, 053601 (2010) where W2(x,y) is a generic potential in the (x,y) plane, while the harmonic potential of angular frequency ω2z along the z axis is so strong that the second BEC is quasi-2D in the (x,y) plane. For the first BEC we still use the potential of Eq. (5) and the variational wave function of Eq. (6), but for the second BEC we adopt the variational wave function ψ2(r,t) = h2(x,y,t) e−z2/(2a2 2z) (π1/2a2z)1/2 , (24) where a2z = √ h̄/(m2ω2z) is the characteristic harmonic length of the potential acting on the second BEC along the z direction. From Eqs. (11) and (12), using Eq. (24), it is straightforward to derive the effective BNLSE for f1(z,t) and h2(x,y,t) [5]: ih̄∂tf1 = [ − h̄2 2m1 ∂2 z +V1(z) + E1 + g1 2πa1xa1y |f1|2 + g12e −z2/a2 2z π3/2a1xa1ya2z ∫ ∞ −∞ e−(x2/a2 1x+y2/a2 1y )|h2|2dxdy ] f1, (25) ih̄∂th2 = [ − h̄2 2m2 ∇2 xy + W2(x,y) + g2 (2π )1/2a2z |h2|2 +E2 + g12 e−(x2/a2 1x+y2/a2 1y ) π3/2a1xa1ya2z ∫ ∞ −∞ e−z2/a2 2z |f1|2dz ] h2, (26) E2 = h̄2 4m2a 2 2z + m2ω 2 2za 2 2z 4 . (27) Equations (25) and (26) are simplified if we impose cylindrical symmetry in the harmonic potential acting on the first BEC, that is, ω1x = ω1y = ω1ρ , and also in the generic potential acting on the second BEC, that is, W2(x,y) = W2(ρ) with ρ = (x2 + y2)1/2. In this way we can choose h2(x,y,t) = h2(ρ,t) and find the following effective BNLSE for f1(z,t) and h2(ρ,t): ih̄∂tf1 = [ − h̄2 2m1 ∂2 z + V1(z) + E1 + g1 2πa2 1ρ |f1|2 + E1 + 2g12e −z2/a2 2z π1/2a2 1ρa2z ∫ ∞ 0 e−ρ2/a2 1ρ |h2|2ρdρ ] f1, (28) ih̄∂th2 = [ − h̄2 2m2 ∇2 ρ + W2(ρ) + E2 + g2 (2π )1/2a2z |h2|2 + g12 e−ρ2/a2 1ρ π3/2a2 1ρa2z ∫ ∞ −∞ e−z2/a2 2z |f1|2dz ] h2. (29) These are two coupled equations, where the first equation is integro-differential in 1 + 1 dimensions and the second is integro-differential in 2 + 1 dimensions (but in 1 + 1 variables). C. 3D-2D binary BEC Let us suppose now that the first BEC is fully 3D and it is under the action of the generic potential U1(r), while the second BEC is trapped by the potential (23). Again we assume that the harmonic potential of angular frequency ω2z along the z axis is so strong that the second BEC is quasi-2D in the (x,y) plane. The first BEC is described by the generic wave function ψ1(r,t) while the second BEC is modeled by the variational wave function of Eq. (24). Following the same procedure above, the wave functions ψ1(r,t) and h2(x,y,t) are found to satisfy the following effective BNLSE ih̄∂tψ1 = [ − h̄2 2m1 ∇2 r + U1(r) + g1|ψ1|2 + g12 π1/2a2z e−z2/a2 2z |h2|2 ] ψ1, (30) ih̄∂th2 = [ − h̄2 2m2 ∇2 xy + W2(x,y) + g2 (2π )1/2a2z |h2|2 +E2 + g12 π1/2a2z ∫ ∞ −∞ e−z2/a2 2z |ψ1|2dz ] h2. (31) As in the previous section, these equations are simplified if the confining potentials have cylindrical symmetry: U1(r) = U1(ρ,z) and W2(x,y) = W2(ρ). In this way we can choose ψ1(r,t) = ψ1(ρ,z,t) and h2(x,y,t) = h2(ρ,t) and the effective BNLSEs become ih̄∂tψ1 = [ − h̄2 2m1 (∇2 ρ + ∂2 z ) + U1(ρ,z) + g1|ψ1|2 + g12 π1/2az e−z2/a2 2z |h2|2 ] ψ1, (32) ih̄∂th2 = [ − h̄2 2m2 ∇2 ρ + W2(ρ) + E2 + g2 (2π )1/2a2z |h2|2 + g12 π1/2a2z ∫ ∞ −∞ e−z2/a2 2z |ψ1|2dz ] h2. (33) These are two coupled equations, where the first is differential in 3 + 1 dimensions (but in 2 + 1 variables) and the second is integro-differential in 2 + 1 dimensions (but in 1 + 1 variables). III. DIMENSIONAL REDUCTION BEYOND THE EFFECTIVE BNLSE: THE BNPSE In Sec. II we considered dimensional reduction of the GPE for binary BEC in the form of effective BNLSEs in reduced dimensions. Such equations are simpler than the 3D binary GPE and they could be good approximations to the original equations. It is possible to improve on the effective BNLSE presented in Sec. II through a binary nonpolynomial Schrödinger equation (BNPSE). This improvement is obtained by considering a more flexible variational wave function in place of Eqs. (6), (16), and (24), where in the transverse direc- tion the wave function is taken as the ground state in the res- pective harmonic potential(s) with Gaussian form. In the following we maintain the Gaussian form but with a flexible width to obtain the required BNPSEs in two cases: (i) the axially symmetric 1D-3D case described by Eqs. (13) and (14) and (ii) the axially symmetric perpendicular 1D-2D case described by Eqs. (28) and (29). The other cases presented in Sec. II can be treated in a similar and straightforward fashion and the respective BNPSE obtained. 053601-4 DIMENSIONAL REDUCTION OF A BINARY BOSE- . . . PHYSICAL REVIEW A 82, 053601 (2010) A. 1D-3D with cylindrical symmetry Here we would like to introduce the next-order correction to Eqs. (13) and (14). In this case the confining potential (5) of the first BEC can be written as U1(r) = 1 2m1ω 2 1ρρ 2 + V1(z). (34) We can then use the variational wave function ψ1(r,t) = e−ρ2/[2σ 2(z,t)] π1/2σ (z,t) f1(z,t) (35) in place of Eq. (6), for the first BEC, where σ ≡ σ (z,t) is its width, which can be quite different from the characteristic harmonic length aρ = √ h̄/(m1ω1ρ) of the transverse con- finement. We assume, consistent with the philosophy of the reduction scheme, that in Eq. (35) the z and t derivatives act only on f1. Then the action of the system becomes A = ∫ L̂1dzdt + ∫ L̂2d 3rdt + ∫ L̂12d 3rdt, (36) where L̂1 = i h̄ 2 (f ∗ 1 ∂tf1 − f1∂tf ∗ 1 ) − h̄2 2m1 |∂zf1|2 − V1(z)|f1|2 − ( h̄2 2m1σ 2 + m1ω 2 1ρσ 2 2 ) |f1|2 − g1 4πσ 2 |f1|4, (37) L̂2 = i h̄ 2 (ψ∗ 2 ∂tψ2 − ψ2∂tψ ∗ 2 ) − h̄2 2m2 |∇ψ2|2 −U2(ρ,z)|ψ2|2 − g2 2 |ψ2|4, (38) L̂12 = − g12 πσ 2 e−ρ2/σ 2 |f1|2|ψ2|2, (39) where we have neglected, in the spirit of the reduction scheme, the spatial derivative of σ . We now extremize the action A with respect to f ∗ 1 (z,t), ψ∗ 2 (ρ,z,t), and σ (z,t), finding the equations for f1(z,t), ψ2(ρ,z,t), and σ (z,t): ih̄∂tf1 = [ − h̄2 2m1 ∂2 z + V1(z) + g1 2πσ 2 |f1|2 + ( h̄2 2m1σ 2 + m1ω 2 1ρσ 2 2 ) + 2g12 σ 2 ∫ ∞ 0 e−ρ2/σ 2 |ψ2(ρ,z,t)|2ρdρ ] f1, (40) ih̄∂tψ2 = [ − h̄2 2m2 ∇2 ρz + U2(ρ,z) + g2|ψ2|2 + g12 πσ 2 e−ρ2/σ 2 |f1(z,t)|2 ] ψ2, (41) m1ω 2 1ρσ 4 = h̄2 m1 + g1 2π |f1|2 + 4g12 ∫ ∞ 0 dρρe−ρ2/σ 2 |ψ2|2 ( 1−ρ2 σ 2 ) . (42) Equations (40) and (41) are the BNPSEs, where the first is integro-differential in 1 + 1 dimensions, and the second is differential in 3 + 1 dimensions (but in 2 + 1 variables). Equation (42) is integro-algebraic and determines the width σ as a function of |f1|2 and |ψ2|2. If we take σ = a1ρ , the harmonic length in the transverse ρ direction for the first component, in place of Eq. (42), we get the previously obtained effective 1D-3D Eqs. (13) and (14). If we set the coupling term g12 = 0 in Eq. (42) we obtain the model previously obtained in Ref. [5] for a cigar-shaped BEC. B. Perpendicular 1D-2D with cylindrical symmetry Here we would like to introduce the next-order correction to Eqs. (28) and (29). In this case the confining potential is given by Eq. (34) and the variational wave function of the first BEC will be taken as Eq. (35). The potential of the second BEC will be taken as U2(r) = W2(ρ) + 1 2m2ω 2 2zz 2, (43) and its variational wave function as ψ2(r,t) = e−z2/(2η2(ρ,t)) π1/4η1/2(ρ,t) h2(ρ,t), (44) where the width η ≡ η(ρ,t) can be quite different from the characteristic harmonic length along the z direction and the ρ and t derivatives are supposed to act only on h2. The action of the system can now be written as A = ∫ L̂1dzdt + 2π ∫ L̄2ρdρdt + ∫ L̄12d 3rdt, (45) where L̂1 is given by Eq. (37) and L̄2 = i h̄ 2 (h∗ 2∂th2 − h2∂th ∗ 2) − h̄2 2m2 |∇ρh2|2 − W2(ρ)|h2|2 − 1 2 ( h̄2 2m2η2 + m2ω 2 2zη 2 2 ) |h2|2 − g2 2η √ 2π |h2|4, (46) L̄12 = − g12 π3/4ησ 2 e−ρ2/σ 2 e−z2/η2 |f1|2|h2|2. (47) Again we have neglected the space derivatives of the widths σ and η. We now extremize the action A with respect to f ∗ 1 , h∗ 2, σ , and η to get the desired equations: ih̄∂tf1 = [ − h̄2 2m1 ∂2 z + V1(z) + g1 2πσ 2 |f1|2 + ( h̄2 2m1σ 2 + m1ω 2 1ρσ 2 2 ) + 2g12e −z2/η2 π1/2σ 2η ∫ ∞ 0 e−ρ2/σ 2 |h2(ρ,t)|2ρdρ ] f1, (48) ih̄∂th2 = [ − h̄2 2m2 ∇2 ρ + W2(ρ) + g2 (2π )1/2η |h2|2 + 1 2 ( h̄2 2m2η2 + m2ω 2 2zη 2 2 ) + g12 e−ρ2/σ 2 π3/2σ 2η ∫ ∞ −∞ e−z2/η2 |f1(z,t)|2dz ] h2, (49) 053601-5 LUIS E. YOUNG-S., L. SALASNICH, AND S. K. ADHIKARI PHYSICAL REVIEW A 82, 053601 (2010) m1ω 2 1ρσ 4 = h̄2 m1 + g1 2π |f1|2 + 4g12π 1/4 ∫ ∞ 0 dρρ 1 η × e−ρ2/σ 2 e−z2/η2 |h2|2 ( 1 − ρ2 σ 2 ) . (50) m2ω 2 2zη 4 = h̄2 m2 + g2η 2 √ 2π |h2|2 + g12η π3/4 ∫ ∞ −∞ dz 1 σ 2 × e−ρ2/σ 2 e−z2/η2 |f1|2 ( 1 − 2 z2 η2 ) . (51) Equations (48) and (49) constitute the BNPSE, where the first one is integro-differential in 1 + 1 dimensions and the second is integro-differential in 2 + 1 dimensions (but in 1 + 1 variables). Equations (50) and (51) are integro-algebraic and determine the widths σ and η as functions of |f1|2 and |h2|2. If we take σ = a1ρ and η = a2z, the harmonic lengths in the transverse directions for the first and the second components, respectively, in place of Eqs. (50) and (51), we get back Eqs. (28) and (29). If we set the coupling term g12 = 0 in Eqs. (50) and (51), Eqs. (48) and (49) reduce to previously studied models in Ref. [5] for cigar- and disk-shaped BECs, respectively. IV. NUMERICAL RESULTS We present results for the coupled 1D-2D case when the cigar-shaped BEC is perpendicular to the disk-shaped BEC. The potentials V1(z) and W2(ρ) in Eqs. (28) and (29) are taken as V1(z) = 1 2m1ω 2 1zz 2; W2(ρ) = 1 2m2ω 2 2ρρ 2. (52) The cigar-shaped BEC, labeled 1, lies along the z axis and has trapping frequencies satisfying ω1ρ/ω1z = 10. The disk- shaped BEC, labeled 2, lies in the (x,y) plane and has trapping frequencies satisfying ω2z/ω2ρ = 10. To keep the algebra under control, we work in dimension- less units and set m1 = m2 = h̄ = 1 and take N1 = N2 = N . For the cigar-shaped BEC we take ω1ρ = 10 and ω1z = 1; for the disk-shaped BEC we take ω2z = 10 and ω2ρ = 1. For the cigar-shaped BEC 1, the harmonic oscillator lengths are a1z ≡ √ h̄/(mω1z) = 1 and a1ρ ≡ √ h̄/(mω1ρ) = 1/ √ 10. For the disk-shaped BEC 2, the harmonic oscillator lengths are a2z ≡ √ h̄/(mω2z) = 1/ √ 10 and a2ρ ≡ √ h̄/(mω2ρ) = 1. The linear and radial densities are calculated from the solution of Eqs. (3) and (4) via |f1(z)|2 = 2π ∫ ∞ 0 ρdρ|ψ1(r)|2, (53) |h2(ρ)|2 = ∫ ∞ −∞ dz|ψ1(r)|2, (54) to be compared with the solutions of Eqs. (28) and (29). To solve the coupled GPEs in three and reduced dimen- sions we employ the split-step Crank-Nicolson discretization technique [20] with real and imaginary time propagation. We essentially use the Fortran programs provided in Ref. [20]. The convergence of the numerical result was tested by varying the 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 0 0.04 0.08 0.12 0.16 0.2 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 16π (3D) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 0 0.2 0.4 0.6 0.8 1 -3 -2 -1 0 1 2 3 4 lin ea r de ns ity |f 1( z) |2 z Ng1 = 0 Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 (a) Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 16π (3D) FIG. 1. (Color online) Plot of (a) linear and (b) radial densities of a binary perpendicular cigar-disk mixture as obtained from the 3D GPEs (3) and (4) (denoted 3D in the figure) and the reduced 1D-2D BNLSEs (28) and (29) (denoted NL in the figure) in dimensionless units for Ng1 = 0 and Ng2 = 40π for different Ng12. Parameters used: m1 = m2 = h̄ = 1, ω1ρ = ω2z = 10,ω1z = ω2ρ = 1, a1z = a2ρ = 1, and a1ρ = a2z = 1/ √ 10. Both densities are normalized to unity. space (typically ∼0.05) and time (typically ∼0.001) steps and the total number of spatial and temporal discretization points. A. Stationary states In Figs. 1(a) and 1(b) we plot the linear and radial densities |f1(z)|2 and h2(ρ)|2, respectively, obtained from the 3D equations, Eqs. (3) and (4), as well as the BNLSEs (28) and (29) for Ng1 = 0 and Ng2 = 40π and for different Ng12. In this case in the absence of coupling (g12 = 0) the two linear densities are identical and the two radial densities are very close to each other. The BNLSEs (28) and (29) provide a faithful representation of the densities even in the presence of coupling as can be seen from Fig. 1. In Figs. 2(a) and 2(b) we plot the linear and radial densities |f1(z)|2 and h2(ρ)|2, respectively, as in Figs. 1(a) and 1(b), for Ng1 = 4π and Ng2 = 40π . In this case even for g12 = 0, there is some difference between the results of the binary 3D GPE and the BNLSE. This difference is consistent with the previous studies [5,8]. The agreement between the densities obtained from the binary 3D GPE and the BNLSE continues to be good. The small discrepancy between the results of the binary 3D GPE and the BNLSE continues as a nonzero g12 is introduced. This discrepancy does not increase as g12 is increased to moderate values, as we see from Fig. 2. In Figs. 3(a) and 3(b) we plot the linear and radial densities for Ng1 = 40π and Ng2 = 400π . In this case there is a larger difference between the results of 3D and reduced 1D-2D models. What is interesting is that the agreement between the results of the two models does not deteriorate with the increase of g12 up to 128π . 053601-6 DIMENSIONAL REDUCTION OF A BINARY BOSE- . . . PHYSICAL REVIEW A 82, 053601 (2010) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 0 0.08 0.16 0.24 0.32 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 16π (3D) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 0 0.02 0.04 0.06 0.08 0.1 0 1 2 3 4 5 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 4π Ng2 = 40π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -4π (NL) -4π (3D) 0 (NL) 0 (3D) 4π (NL) 4π (3D) 16π (NL) 16π (3D) FIG. 2. (Color online) The same as in Fig. 1 for Ng1 = 4π and Ng2 = 40π . From the results displayed in Figs. 1–3 we find that the BNLSEs (28) and (29) provide an excellent account of the actual state of affairs, when compared with the full 3D binary GPEs (3) and (4), especially for small nonlinearities Ng1 and Ng2. We next attempt an improvement over the BNLSEs (28) and (29) of Sec. II by applying BNPSEs of Sec. III. In other words, we attempt a solution of Eqs. (48)–(51), in place of Eqs. (28) and (29), to reduce the discrepancy noted in Figs. 3(a) and 3(b) between the results of the 3D equation and the BNLSE. The results so obtained are plotted in Figs. 4(a) and 4(b) for Ng1 = 40π and Ng2 = 400π , which show a significant improvement over the results in Figs. 3(a) and 3(b). The large discrepancy between the results of the 3D binary GPE and the BNLSE for 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 16π (3D) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 16π (3D) 128π (NL) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 16π (3D) 128π (NL) 128π (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 16π (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 16π (3D) 128π (NL) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NL) -8π (3D) 0 (NL) 0 (3D) 16π (NL) 16π (3D) 128π (NL) 128π (3D) FIG. 3. (Color online) The same as in Fig. 1 for Ng1 = 4π and Ng2 = 400π . 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 16π (3D) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 16π (3D) 128π (NP) 0 0.04 0.08 0.12 0.16 -8 -6 -4 -2 0 2 4 6 8 10 lin ea r de ns ity |f 1( z) |2 z (a) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 16π (3D) 128π (NP) 128π (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 16π (3D) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 16π (3D) 128π (NP) 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 ra di al d en si ty |h 2( ρ) |2 ρ (b) Ng1 = 40π Ng2 = 400π ω1ρ/ω1z = 10 ω2z/ω2ρ = 10 Ng12 = -8π (NP) -8π (3D) 0 (NP) 0 (3D) 16π (NP) 16π (3D) 128π (NP) 128π (3D) FIG. 4. (Color online) The same as in Fig. 3 obtained using the reduced BNPSEs (48)–(51) (NP). the linear density of component 1 has virtually disappeared in Fig. 4(a) and that for the radial density of component 2 is significantly reduced in Fig. 4(b). This shows that although the results from the effective BNLSE of Sec. II are quite good for small nonlinearities, significant improvement can be obtained from the BNPSE of Sec. III, especially for larger values of nonlinearities. Finally, we consider the chemical potential µ1 and µ2 of the two components calculated using the 3D binary GPEs 0 10 20 30 40 0 100 200 300 400 500 600 ch em p ot g2N/π Ng12 = 128π Ng1 = 40π (a) µ1 (NL) 0 10 20 30 40 0 100 200 300 400 500 600 ch em p ot g2N/π Ng12 = 128π Ng1 = 40π (a) µ1 (NL) µ1 (NP) 0 10 20 30 40 0 100 200 300 400 500 600 ch em p ot g2N/π Ng12 = 128π Ng1 = 40π (a) µ1 (NL) µ1 (NP) µ1 (3D) 0 10 20 30 40 0 100 200 300 400 500 600 ch em p ot g2N/π Ng12 = 128π Ng1 = 40π (a) µ1 (NL) µ1 (NP) µ1 (3D) µ2 (NL) 0 10 20 30 40 0 100 200 300 400 500 600 ch em p ot g2N/π Ng12 = 128π Ng1 = 40π (a) µ1 (NL) µ1 (NP) µ1 (3D) µ2 (NL) µ2 (NP) 0 10 20 30 40 0 100 200 300 400 500 600 ch em p ot g2N/π Ng12 = 128π Ng1 = 40π (a) µ1 (NL) µ1 (NP) µ1 (3D) µ2 (NL) µ2 (NP) µ2 (3D) 0 10 20 30 40 0 20 40 60 80 ch em p ot g1N/π Ng12 = 128π Ng2 = 400π (b) µ1 (NL) 0 10 20 30 40 0 20 40 60 80 ch em p ot g1N/π Ng12 = 128π Ng2 = 400π (b) µ1 (NL) µ1 (NP) 0 10 20 30 40 0 20 40 60 80 ch em p ot g1N/π Ng12 = 128π Ng2 = 400π (b) µ1 (NL) µ1 (NP) µ1 (3D) 0 10 20 30 40 0 20 40 60 80 ch em p ot g1N/π Ng12 = 128π Ng2 = 400π (b) µ1 (NL) µ1 (NP) µ1 (3D) µ2 (NL) 0 10 20 30 40 0 20 40 60 80 ch em p ot g1N/π Ng12 = 128π Ng2 = 400π (b) µ1 (NL) µ1 (NP) µ1 (3D) µ2 (NL) µ2 (NP) 0 10 20 30 40 0 20 40 60 80 ch em p ot g1N/π Ng12 = 128π Ng2 = 400π (b) µ1 (NL) µ1 (NP) µ1 (3D) µ2 (NL) µ2 (NP) µ2 (3D) FIG. 5. (Color online) (a) Chemical potential of the two com- ponents µ1 and µ2 calculated using the 3D GPEs (3) and (4) (3D), the reduced BNPSEs (48)–(51) (NP), and the reduced BNLSEs (28) and (29) (NL) versus g2N/π for Ng1 = 40π and Ng12 = 128π . (b) Chemical potential of the two components µ1 and µ2 calculated using the 3D GPEs (3) and (4) (3D), the reduced BNPSEs (48)–(51) (NP), and the reduced BNLSEs (28) and (29) (NL) versus g1N/π for Ng2 = 400π and Ng12 = 128π . 053601-7 LUIS E. YOUNG-S., L. SALASNICH, AND S. K. ADHIKARI PHYSICAL REVIEW A 82, 053601 (2010) (3) and (4), the BNPSEs (48)–(51) as well as the BNLSEs (28) and (29). The equations for the chemical potentials are obtained by replacing the time derivative ih̄∂t by µ1 and µ2, respectively, in the dynamical equations for components 1 and 2, respectively. The chemical potentials are plotted in Fig. 5(a) versus g1N/π for Ng2 = 400π and Ng12 = 128π . In Fig. 5(b) we plot the chemical potentials versus g2N/π for Ng1 = 40π and Ng12 = 128π . The chemical potentials obtained from the BNLSEs (28) and (29) are good but those calculated from the BNPSEs (48)–(51) are better approximations to the results obtained from the 3D binary GPEs (3) and (4) in all cases. B. Dynamics Now we investigate if the reduced equation presented here can satisfactorily describe the dynamics of the perpendicular 1D-2D system. For this purpose first we consider the coupled system presented in Figs. 2(a) and 2(b) with Ng1 = 4π,Ng2 = 40π, and Ng12 = 16π . After the stationary condensate is created during time evolution of the numerical routine using real-time propagation, we reduce the radial angular frequency of the disk-shaped quasi-2D BEC suddenly by 5% (the new angular frequencies being ω1ρ = ω2z = 10,ω1z = 1,ω2ρ = 0.95) and study the subsequent dynamical evolution of the system. The relevant dynamics is best illustrated by plotting the evolution of the rms axial size of the cigar-shaped quasi-1D BEC and the rms radial size of the disk-shaped quasi-2D BEC 1.5 2 2.5 0 5 10 15 20 25 rm s si ze ( z or ρ ) time t Ng1 = 4π Ng2 = 40π Ng12 = 16π (a) Ng1 = 4π Ng2 = 40π Ng12 = 16π (a) Ng1 = 4π Ng2 = 40π Ng12 = 16π (a) Ng1 = 4π Ng2 = 40π Ng12 = 16π (a) ρ2(NL) ρ2(NP) ρ2(3D) z1(NL) z1(NP) z1(3D) 2.5 3 3.5 4 0 5 10 15 20 25 rm s si ze ( z or ρ ) time t Ng1 = 40π Ng2 = 400π Ng12 = 128π (b) ρ2(NL) ρ2(NP) ρ2(3D) z1(NL) z1(NP) z1(3D) FIG. 6. (Color online) Evolution of axial and radial rms sizes of the cigar- and disk-shaped components of the binary BEC after a sudden reduction of the radial angular frequency of the disk-shaped component by 5%, for the nonlinearities (a) Ng1 = 4π , Ng2 = 40π , and Ng12 = 16π and (b) Ng1 = 40π , Ng2 = 400π , and Ng12 = 128π . The cigar-shaped BEC is perpendicular to the disk-shaped BEC. in Figs. 6(a) and 6(b) for two sets of nonlinearities exhibited in Figs. 2 and 3. Here we show the results of the 3D GPEs (3) and (4) (3D), the reduced BNPSEs (48)–(51) (NP), and the reduced BNLSEs (28) and (29) (NL). For both sets of nonlinearities illustrated in Fig. 6(a) Ng1 = 4π , Ng2 = 40π , and Ng12 = 16π and Fig. 6(b) Ng1 = 40π , Ng2 = 400π , and Ng12 = 128π , the NP reduction produced far better results compared to the NL reduction. In both cases the angular frequency of the quasi-2D disk-shaped BEC was modified and both the NL and NP schemes as well as the full 3D calculation exhibited sinusoidal oscillation for the axial size of the disk-shaped BEC right from t = 0. It took some 5 units of time for similar sinusoidal oscillation to initiate in the unperturbed cigar-shaped component due to the nonzero coupling term g12. The oscillations in rms sizes in the NP scheme were always in phase with those of the full 3D calculation. But for larger nonlinearity and at large times the oscillations in rms sizes in the NL scheme are not always in phase with those of the full 3D calculation. The differences between the rms sizes of the NL and NP schemes compared to the full 3D calculation were consequences of the same differences, existing at t = 0, in the corresponding stationary results in the absence of any perturbation in the angular frequency. V. SUMMARY In this article we derived and studied simple BNPSEs for the description of a binary BEC where the two components subject to distinct trap symmetries belong to two different spatial dimensions. We considered three possibilities: the first component is in 3D and the second in (a) 1D (cigar shape) or (b) 2D (disk shape), and (c) the first component is in 1D and the second is in 2D. The 1D and 2D configurations were achieved in an axially symmetric setting with a strong harmonic trap in the transverse radial and axial directions, respectively. The BNPSEs were obtained by a Lagrangian variational scheme after approximating the spatial wave-function component in the transverse direction by a generic Gaussian function. A simpler set of reduced equations, the BNLSEs with power-law nonlinearity, were obtained with a simple ansatz for the 3D wave function where the wave-function component in the transverse direction is taken to be the ground state in the corresponding harmonic trap. We tested the accuracy of the aforementioned reduction scheme by numerically solving the 3D binary GPE in the coupled 1D-2D case where the cigar-shaped component is in the z direction and the disk-shaped component lies in the (x,y) plane. We calculated the density profiles and the chemical potentials of the two components. For smaller values of inter- and intraspecies nonlinearities, the BNLSEs (28) and (29) with power-law nonlinearity produced results for density and chem- ical potential in good agreement with the binary 3D GPEs. For larger nonlinearities the BNPSEs (48)–(51) produced signif- icant improvement over the BNLSEs (28) and (29). We also studied small oscillation of the binary 1D-2D BEC, initiated by a sudden change of the radial angular frequency of the 2D com- ponent. For larger nonlinearities and large times, the BNPSEs (48)–(51) produced better results for dynamics compared with the BNLSEs (28) and (29). To conclude, the BNLSEs (28) 053601-8 DIMENSIONAL REDUCTION OF A BINARY BOSE- . . . PHYSICAL REVIEW A 82, 053601 (2010) and (29) as well as the BNPSEs (48)–(51) are shown to be very useful in investigating a binary 1D-2D BEC, where the cigar-shaped BEC is perpendicular to the disk-shaped one. Similar reduced equations are derived for other symmetries, such as the 1D-3D and 2D-3D cases as well as the 1D-2D case where the cigar-shaped BEC is in the plane of the disk-shaped one. ACKNOWLEDGMENTS CNPq and FAPESP (Brazil) provided partial support. [1] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Clarendon Press, Oxford, 2003). [2] A. J. Leggett, Quantum Liquids: Bose Condensation and Cooper Pairing in Condensed-Matter Systems (Oxford University Press, Oxford, 2006). [3] V. M. Pérez-Garcı́a, H. Michinel, and H. Herrero, Phys. Rev. A 57, 3837 (1998). [4] A. D. Jackson, G. M. Kavoulakis, and C. J. Pethick, Phys. Rev. A 58, 2417 (1998). [5] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 65, 043614 (2002); L. Salasnich, Laser Phys. 12, 198 (2002). [6] P. Massignan and M. Modugno, Phys. Rev. A 67, 023614 (2003); M. Modugno, C. Tozzo, and F. Dalfovo, ibid. 70, 043625 (2004); C. Tozzo, M. Kramer, and F. Dalfovo, ibid. 72, 023613 (2005); M. Modugno, ibid. 73, 013606 (2006). [7] A. Munoz Mateo and V. Delgado, Phys. Rev. A 75, 063610 (2007); 77, 013617 (2008); Ann. Phys. 324, 709 (2009). [8] C. A. G. Buitrago and S. K. Adhikari, J. Phys. B 42, 215306 (2009). [9] G. Theocharis, P. G. Kevrekidis, M. K. Oberthaler, and D. J. Frantzeskakis, Phys. Rev. A 76, 045601 (2007); A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, M. K. Oberthaler, D. J. Frantzeskakis, G. Theocharis, and P. G. Kevrekidis, Phys. Rev. Lett. 101, 130401 (2008). [10] M. L. Chiofalo and M. P. Tosi, Phys. Lett. A 268, 406 (2000); A. M. Kamchatnov and V. S. Shchesnovich, Phys. Rev. A 70, 023604 (2004); W. Zhang and L. You, ibid. 71, 025603 (2005); F. K. Abdullaev and R. Galimzyanov, J. Phys. B 36, 1099 (2003). [11] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A 66, 043603 (2002); L. Salasnich, ibid. 70, 053617 (2004); L. Salasnich, A. Parola, and L. Reatto, ibid. 70, 013606 (2004); A. Parola, L. Salasnich, R. Rota, and L. Reatto, ibid. 72, 063612 (2005). [12] A. Maluckov, L. Hadzievski, B. A. Malomed, and L. Salasnich, Phys. Rev. A 78, 013616 (2008); L. Salasnich, J. Phys. A: Math. Theor. 42, 335205 (2009); G. Gligoric, A. Maluckov, L. Hadzievski, L. Salasnich, and B. A. Malomed, Chaos 19, 043105 (2009). [13] L. Salasnich and B. A. Malomed, Phys. Rev. A 79, 053620 (2009). [14] S. K. Adhikari and L. Salasnich, New J. Phys. 11, 023011 (2009). [15] Y. E. Kim and A. L. Zubarev, Phys. Rev. A 70, 033612 (2004); N. Manini and L. Salasnich, ibid. 71, 033625 (2005); G. Diana, N. Manini, and L. Salasnich, ibid. 73, 065601 (2006); W. Wen and G. X. Huang, ibid. 79, 023605 (2009); S. K. Adhikari and L. Salasnich, ibid. 78, 043616 (2008); S. K. Adhikari, Laser Phys. Lett. 6, 901 (2009); Phys. Rev. A 77, 045602 (2008); J. Phys. B 43, 085304 (2010); L. Salasnich and F. Toigo, Phys. Rev. A 78, 053626 (2008); L. Salasnich, F. Ancilotto, and F. Toigo, Laser Phys. Lett. 7, 78 (2010). [16] S. K. Adhikari and L. Salasnich, Phys. Rev. A 80, 023606 (2009); S. K. Adhikari, ibid. 81, 043636 (2010); Y. Cheng and S. K. Adhikari, ibid. 81, 023620 (2010); 82, 013631 (2010); Laser Phys. Lett. 7, 824 (2010). [17] G. Lamporesi, J. Catani, G. Barontini, Y. Nishida, M. Inguscio, and F. Minardi, Phys. Rev. Lett. 104, 153202 (2010). [18] L. Salasnich and B. A. Malomed, Phys. Rev. A 74, 053610 (2006); L. Salasnich, B. A. Malomed, and F. Toigo, ibid. 81, 045603 (2010). [19] G. P. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonzlez, eds., Emergent Nonlinear Phenomena in Bose-Einstein Conden- sates (Springer, Berlin, 2008). [20] P. Muruganandam and S. K. Adhikari, Comput. Phys. Commun. 180, 1888 (2009). [21] P. Muruganandam and S. K. Adhikari, J. Phys. B 36, 2501 (2003); S. K. Adhikari and P. Muruganandam, ibid. 35, 2831 (2002). 053601-9 http://dx.doi.org/10.1103/PhysRevA.57.3837 http://dx.doi.org/10.1103/PhysRevA.57.3837 http://dx.doi.org/10.1103/PhysRevA.58.2417 http://dx.doi.org/10.1103/PhysRevA.58.2417 http://dx.doi.org/10.1103/PhysRevA.65.043614 http://dx.doi.org/10.1103/PhysRevA.65.043614 http://dx.doi.org/10.1103/PhysRevA.67.023614 http://dx.doi.org/10.1103/PhysRevA.70.043625 http://dx.doi.org/10.1103/PhysRevA.72.023613 http://dx.doi.org/10.1103/PhysRevA.73.013606 http://dx.doi.org/10.1103/PhysRevA.75.063610 http://dx.doi.org/10.1103/PhysRevA.75.063610 http://dx.doi.org/10.1103/PhysRevA.77.013617 http://dx.doi.org/10.1016/j.aop.2008.10.002 http://dx.doi.org/10.1016/j.aop.2008.10.002 http://dx.doi.org/10.1088/0953-4075/42/21/215306 http://dx.doi.org/10.1088/0953-4075/42/21/215306 http://dx.doi.org/10.1103/PhysRevA.76.045601 http://dx.doi.org/10.1103/PhysRevLett.101.130401 http://dx.doi.org/10.1103/PhysRevLett.101.130401 http://dx.doi.org/10.1016/S0375-9601(00)00197-3 http://dx.doi.org/10.1103/PhysRevA.70.023604 http://dx.doi.org/10.1103/PhysRevA.70.023604 http://dx.doi.org/10.1103/PhysRevA.71.025603 http://dx.doi.org/10.1088/0953-4075/36/6/302 http://dx.doi.org/10.1088/0953-4075/36/6/302 http://dx.doi.org/10.1103/PhysRevA.66.043603 http://dx.doi.org/10.1103/PhysRevA.66.043603 http://dx.doi.org/10.1103/PhysRevA.70.053617 http://dx.doi.org/10.1103/PhysRevA.70.013606 http://dx.doi.org/10.1103/PhysRevA.72.063612 http://dx.doi.org/10.1103/PhysRevA.72.063612 http://dx.doi.org/10.1103/PhysRevA.78.013616 http://dx.doi.org/10.1088/1751-8113/42/33/335205 http://dx.doi.org/10.1088/1751-8113/42/33/335205 http://dx.doi.org/10.1063/1.3248269 http://dx.doi.org/10.1063/1.3248269 http://dx.doi.org/10.1103/PhysRevA.79.053620 http://dx.doi.org/10.1103/PhysRevA.79.053620 http://dx.doi.org/10.1088/1367-2630/11/2/023011 http://dx.doi.org/10.1103/PhysRevA.70.033612 http://dx.doi.org/10.1103/PhysRevA.70.033612 http://dx.doi.org/10.1103/PhysRevA.71.033625 http://dx.doi.org/10.1103/PhysRevA.71.033625 http://dx.doi.org/10.1103/PhysRevA.73.065601 http://dx.doi.org/10.1103/PhysRevA.73.065601 http://dx.doi.org/10.1103/PhysRevA.79.023605 http://dx.doi.org/10.1103/PhysRevA.79.023605 http://dx.doi.org/10.1103/PhysRevA.78.043616 http://dx.doi.org/10.1103/PhysRevA.78.043616 http://dx.doi.org/10.1002/lapl.200910090 http://dx.doi.org/10.1103/PhysRevA.77.045602 http://dx.doi.org/10.1088/0953-4075/43/8/085304 http://dx.doi.org/10.1103/PhysRevA.78.053626 http://dx.doi.org/10.1002/lapl.200910107 http://dx.doi.org/10.1002/lapl.200910107 http://dx.doi.org/10.1103/PhysRevA.80.023606 http://dx.doi.org/10.1103/PhysRevA.81.043636 http://dx.doi.org/10.1103/PhysRevA.81.023620 http://dx.doi.org/10.1103/PhysRevA.82.013631 http://dx.doi.org/10.1002/lapl.201010063 http://dx.doi.org/10.1002/lapl.201010063 http://dx.doi.org/10.1103/PhysRevLett.104.153202 http://dx.doi.org/10.1103/PhysRevA.74.053610 http://dx.doi.org/10.1103/PhysRevA.74.053610 http://dx.doi.org/10.1103/PhysRevA.81.045603 http://dx.doi.org/10.1103/PhysRevA.81.045603 http://dx.doi.org/10.1016/j.cpc.2009.04.015 http://dx.doi.org/10.1016/j.cpc.2009.04.015 http://dx.doi.org/10.1088/0953-4075/36/12/310 http://dx.doi.org/10.1088/0953-4075/36/12/310 http://dx.doi.org/10.1088/0953-4075/35/12/317 http://dx.doi.org/10.1088/0953-4075/35/12/317