FACt: FORTRAN toolbox for calculating fluctuations in atomic condensates
FFACt: FORTRAN toolbox for calculating fluctuations inatomic condensates
Arko Roy a,b , Sukla Pal a, ∗ , S. Gautam c , D. Angom a , P. Muruganandam d,e a Physical Research Laboratory, Navarangpura, Ahmedabad 380 009, Gujarat, India b Max-Planck-Institut f¨ur Physik komplexer Systeme, N¨othnitzer Straße 38, D-01187 Dresden, Germany. c Department of Physics, Indian Institute of Technology Ropar, Rupnagar, Punjab 14001, India d Department of Physics, Bharathidasan University, Tiruchirappalli 620 024,India e Department of Medical Physics, Bharathidasan University, Tiruchirappalli 620 024,India
Abstract
We develop a FORTRAN code to compute fluctuations in atomic condensates (FACt)by solving the Bogoliubov-de Gennes (BdG) equations for two component Bose-Einsteincondensate (TBEC) in quasi two dimensions. The BdG equations are recast as ma-trix equations and solved self consistently. The code is suitable for handling quantumfluctuations as well as thermal fluctuations at temperatures below the critical point ofBose-Einstein condensation. The code is versatile, and the ground state density pro-file and low energy excitation modes obtained from the code can be easily adapted tocompute different properties of TBECs — ground state energy, overlap integral, quasiparticle amplitudes of BdG spectrum, dispersion relation and structure factor and otherrelated experimental observables.
Keywords:
Gross-Pitaevskii equation; Hartree-Fock-Bogoliubov theory;Bogoliubov-de Gennes equations; Quasiparticle spectra; Goldstone mode;Kohn/Slosh mode; miscibility-immiscibility transition;
PROGRAM SUMMARY
Program Title:
FACt
Journal Reference:Catalogue identifier:Licensing provisions: none
Programming language:
FORTRAN 90
Computer:
Intel Xeon,
Operating system:
General
RAM: at least 1.5Gbytes per core.
Number of processors used: Supplementary material: none
Classification: ∗ Corresponding author.
E-mail address: [email protected]
Preprint submitted to Computer Physics Communications June 5, 2018 a r X i v : . [ c ond - m a t . qu a n t - g a s ] J un xternal routines/libraries: ARPACK
Subprograms used:Journal reference of previous version: * Nature of problem:
Compute the ground state density profile, ground state energyand chemical potential for individual species, evaluate the quasiparticle mode ener-gies and corresponding amplitudes which can capture the transformation of the modesagainst the change of the parameters (intraspecies interaction, interspecies interaction,anisotropy parameter etc.) using Hartree-Fock Bogoliubov theory with the Popov ap-proximation Calculate the overlap integral, dispersion relation and structure factor.
Solution method:
In the first step, the pair of coupled Gross-Pitaevskii (GP) equationsare solved using split step Crank-Nicolson method to compute the condensate density.To solve BdG equations, as a basic input the first N b harmonic oscillator eigenstates arechosen as a basis to generate the BdG matrix with dimension of N b + 1)( N b + 1) × N b + 1)( N b + 1) . Since the matrix size rapidly increases with N b , A RPACK routinesare used to diagonalise the BdG matrix efficiently. To compute the fluctuation and non-condensate density, a set of the low energy quasiparticle amplitudes above a thresholdvalue of the Bose factor are considered. The equations are then solved iteratively tillthe condensate, and non-condensate densities converge to predefined accuracies. Toaccelerate the convergence we use the method of successive under-relaxation (SUR).
Restrictions:
For a large system size, if the harmonic oscillator basis size is also taken to be large,the dimension of the BDG matrix becomes huge. It may take several days to computethe low energy modes at finite temperature and this package may be computationallyexpensive.
Additional comments:
After successful computation of this package, one should obtain the equilibriumdensity profiles for TBEC, low energy Bogoliubov modes and the corresponding quasi-particle amplitudes. In addition, one can calculate the dispersion relation, structure fac-tor, overlap integral, correlation function, etc. using this package with minimal modi-fications. In the theory section of the manuscript, we have provided the expressions tocompute the above quantities numerically.
Running time: ∼ minutes for the sample case. For self consistent calculation with 15 iterations,it could take approximately 2 days for the parameters specified in the manuscript.
1. Introduction
The self-consistent Hartree-Fock-Bogoliubov theory with the Popov (HFB-Popov)approximation is an effective model to examine the fluctuations of equilibrium state so-lutions of trapped BEC at zero temperature as well as finite temperatures. The theory isin particular well suited to examine the evolution of the low-lying modes as a functionof the interaction parameters, temperature or trapping parameters. It has been used ex-tensively in single-species BEC to study finite temperature effects and mode energies[1–4], and the results are in good agreement with experimental results [5] at low tem-peratures. The detailed and systematic information about the quasiparticle spectrum,2oth of single and multispecies condensate, are described by the HFB formalism. Intwo-species BECs (TBECs), where the phenomenon of phase-separation is important[6, 7], the HFB-Popov approximation has been used in the miscible [8] and immiscibledomain [9–11] to compute the low-lying modes.In the present work we report the development of a FORTRAN code which im-plements the HFB-Popov theory to compute the low energy elementary excitations ofthe TBECs. At T = 0 K, where only the quantum fluctuations are present in the sys-tem, the code captures the essence of quantum fluctuations. These are important inthe stabilization of quantum droplets in binary BEC mixtures [12–15]. In our recentworks [16, 17] we have investigated the elementary excitations in radially symmetricand and anisotropic TBECs using the present version of FACt. However, the mainstrength of HFB-Popov approximation is in encapsulating properties of trapped BECat finite temperatures, which is more realistic and experimentally relevant. It must beemphasized that our code provides high precision and converged results for T (cid:28) T c and computes the low energy excitation modes for TBECs in quasi two dimension.It must also be mentioned that the HFB-Popov has been used to study quantum andthermal fluctuations in optical lattices [18, 19].An important feature of our implementation, which optimizes the computationalrequirements, is the absence of any constraints on the symmetry. That is, we imple-ment the code in Cartesian coordinates. The basic and important advantage of thisapproach is that, our code is very general and applicable to the anisotropic cases wherethe frequency of the trap in x and y directions are different.
2. Finite temperature theory for two component BEC
In the dilute limit, when the interparticle interactions are weak, the nonlinear Sch¨odingerequation (NLSE), also known as the Gross-Pitaevskii equation (GPE) provides a gooddescription of BECs. To incorporate the statics and dynamical properties of TBECs,this equation can be generalized to a pair of coupled GPEs (CGPEs). This, however, isa description valid at zero temperature T = 0 and they form the basis of our compu-tational scheme. Neglecting the quantum fluctuations, the condensed state of TBEC at T = 0 can be described by the macroscopic wave function φ ( x, y, t ) ( φ ( x, y, t ) ) withenergy functional E [ φ ] ( E [ φ ] ) for the first (second) species. The energy functionalof the total system is E = E + E + E = (cid:90) (cid:90) dxdy (cid:20) (cid:88) i =1 (cid:18) (cid:126) m i |∇ φ i | + V i ( x, y ) | φ i | + 12 U ii | φ i | (cid:19) + U | φ | | φ | (cid:21) . (1)where E is the contribution from the interspecies interaction, m i is the mass of thebosonic atom of species i , and V i ( x, y ) is the external harmonic trapping potential. Theinteraction strengths are given by U ij = 2 π (cid:126) a ij /m ij , where m − ij = m − i + m − j is3he reduced mass for an atom i and an atom j . Using these definitions and the mean-field theory, the static and dynamical properties of TBEC, albeit at T = 0 , can beexamined through the time-independent CGPE − (cid:126) m i ∇ + V i ( x, y ) + (cid:88) j =1 U ij | φ j | φ i = µ i φ i , (2)which are obtained by variational minimization of the energy functional E = E − (cid:80) i µ i N i with φ ∗ i as the parameter of variation. The Eq. (2) forms the starting point ofour analysis of TBECs at finite temperatures ( T (cid:54) = 0 ). At equilibrium, depending uponthe relative strengths of intra- ( U ii ) and inter-species ( U ) interactions, the TBECsmay either be in miscible or immiscible phase. The latter is also referred to as phase-separated and we use these two terms interchangeably. The emergence of these phasesrenders the physics of TBEC drastically different from single-species BEC. And, thenatural question is the role of fluctuations, both quantum and thermal, on these phases.For this, the first step is to solve Eqns. (2), and then use the HFB-Popov approximationto calculate the thermal cloud densities.For T (cid:54) = 0 , along with the two coherent condensate clouds, there exist the incoher-ent non-condensate clouds of both the species. This introduces additional interparticleinteractions, the intra- and inter-species interactions between the condensate and non-condensate clouds. The presence of larger number of interaction terms complicates thegoverning equations, and poses difficulty to theoretically model it. In the present work,we have assumed that the thermal clouds of both the species are static, and consider T less than the lower critical temperature among the two. To obtain the Hartree Fock Bogoliubov equation we consider the grand-canonicalHamiltonian for TBECs in a quasi-2D trap, ˆ H = (cid:88) i =1 , (cid:90) (cid:90) dxdy ˆΨ † i ( x, y, t ) (cid:20) − (cid:126) m i ( ∂ ∂x + ∂ ∂y ) + V i ( x, y ) − µ i + U ii † i ( x, y, t ) ˆΨ i ( x, y, t ) (cid:21) ˆΨ i ( x, y, t )+ U (cid:90) (cid:90) dxdy ˆΨ † ( x, y, t ) ˆΨ † ( x, y, t ) ˆΨ ( x, y, t ) ˆΨ ( x, y, t ) , (3)where i = 1 , is the species index, ˆΨ i ’s are the Bose field operators of the two species,and µ i ’s are the chemical potentials. The intra- and interspecies interactions strengthsare U ii = 2 a ii √ πλ and U = 2 a √ πλ (1 + m /m ) , respectively, where λ =( ω z /ω ⊥ ) is the anisotropy parameter. Here, a ii , a represent the s -wave scatteringlengths of intra and inter species interactions respectively. The requirement of havinga quasi-2D geometry is satisfied through the following inequalities: λ (cid:29) , (cid:126) ω z (cid:29) µ i [20, 21] and (cid:126) ω z (cid:29) k B T (at finite temperature T ) [10, 22]. Under these constraintconditions, the motion of the trapped atoms will be confined strongly along z direction4nd the atoms will remain frozen in the ground state providing a quasi-2D confinement.The Heisenberg equation of motion for the Bose field operators ˆΨ i in two-componentnotation is i (cid:126) ∂∂t (cid:18) ˆΨ ˆΨ (cid:19) = (cid:18) ˆ h + U ˆΨ † ˆΨ U ˆΨ † ˆΨ U ˆΨ † ˆΨ ˆ h + U ˆΨ † ˆΨ (cid:19)(cid:18) ˆΨ ˆΨ (cid:19) , (4)where ˆ h i = ( − (cid:126) / m i )( ∂ /∂x + ∂ /∂y ) + V i ( x, y ) − µ i . Using Bogoliubov ap-proximation, the field operators can be written as ˆΨ i ( x, y, t ) = φ i ( x, y ) + ˜ ψ i ( x, y, t ) ,where φ i ( x, y ) is a c -field and represents the condensate, and ˜ ψ i ( x, y, t ) is the fluctu-ation operator corresponding to the i th species. We can write the total field operatoras (cid:18) ˆΨ ˆΨ (cid:19) = (cid:18) φ φ (cid:19) + (cid:18) ˜ ψ ˜ ψ (cid:19) , ⇒ ˆΨ = Φ + ˜Ψ , (5)where Φ and ˜Ψ are the condensate and fluctuation operator in two-component nota-tions. Using the expression of ˆΨ i , we can separate the Hamiltonian into terms of dif-ferent orders in fluctuation operators i.e., ˆ H = (cid:80) i =1 , (cid:80) n =0 ˆ H in , where (cid:54) n (cid:54) denotes the order of the fluctuation operators. The explicit forms of ˆ H in are ˆ H = (cid:90) (cid:90) dxdy φ ∗ (cid:18) ˆ h − µ + U | φ | + U | φ | (cid:19) φ , ˆ H = (cid:90) (cid:90) dxdy φ ∗ (cid:18) ˆ h − µ + U | φ | + U | φ | (cid:19) φ , ˆ H = (cid:90) (cid:90) dxdy (cid:104) φ ∗ (cid:16) ˆ h − µ + U | φ | + U | φ | (cid:17) ˜ ψ + ˜ ψ † (cid:16) ˆ h − µ + U | φ | + U | φ | (cid:17) φ (cid:105) , ˆ H = (cid:90) (cid:90) dxdy (cid:104) φ ∗ (cid:16) ˆ h − µ + U | φ | + U | φ | (cid:17) ˜ ψ + ˜ ψ † (cid:16) ˆ h − µ + U | φ | + U | φ | (cid:17) φ (cid:105) , ˆ H = (cid:90) (cid:90) dxdy (cid:20) ˜ ψ † (cid:16) ˆ h − µ + 2 U | φ | + U | φ | (cid:17) ˜ ψ + U (cid:16) φ ∗ ˜ ψ ˜ ψ + φ ˜ ψ † ˜ ψ † (cid:17) + U (cid:16) φ ∗ φ ∗ ˜ ψ ˜ ψ + φ ∗ φ ˜ ψ † ˜ ψ + φ φ ∗ ˜ ψ † ˜ ψ + φ φ ˜ ψ † ˜ ψ † (cid:17) (cid:21) , ˆ H = (cid:90) (cid:90) dxdy (cid:20) ˜ ψ † (cid:16) ˆ h − µ + 2 U | φ | + U | φ | (cid:17) ˜ ψ + U (cid:16) φ ∗ ˜ ψ ˜ ψ + φ ˜ ψ † ˜ ψ † (cid:17) + U (cid:16) φ ∗ φ ∗ ˜ ψ ˜ ψ + φ ∗ φ ˜ ψ † ˜ ψ + φ φ ∗ ˜ ψ † ˜ ψ + φ φ ˜ ψ † ˜ ψ † (cid:17) (cid:21) , ˆ H = (cid:90) (cid:90) dxdy (cid:20) U (cid:16) φ ∗ ˜ ψ † ˜ ψ ˜ ψ + φ ˜ ψ † ˜ ψ † ˜ ψ (cid:17) + U (cid:0) φ ∗ ˜ ψ † ˜ ψ ˜ ψ + φ ∗ ˜ ψ † ˜ ψ ˜ ψ + φ ˜ ψ † ˜ ψ † ˜ ψ + φ ˜ ψ † ˜ ψ † ˜ ψ (cid:1)(cid:21) , H = (cid:90) (cid:90) dxdy (cid:20) U (cid:16) φ ∗ ˜ ψ † ˜ ψ ˜ ψ + φ ˜ ψ † ˜ ψ † ˜ ψ (cid:17) + U (cid:0) φ ∗ ˜ ψ † ˜ ψ ˜ ψ + φ ∗ ˜ ψ † ˜ ψ ˜ ψ + φ ˜ ψ † ˜ ψ † ˜ ψ + φ ˜ ψ † ˜ ψ † ˜ ψ (cid:1)(cid:21) , ˆ H = (cid:90) (cid:90) dxdy (cid:20) U ψ † ˜ ψ † ˜ ψ ˜ ψ + U ψ † ˜ ψ † ˜ ψ ˜ ψ (cid:21) , ˆ H = (cid:90) (cid:90) dxdy (cid:20) U ψ † ˜ ψ † ˜ ψ ˜ ψ + U ψ † ˜ ψ † ˜ ψ ˜ ψ (cid:21) . (6)Using the definition of field operator from Eq. (5) and putting it in Eq. (4), theHeisenberg equation of motion for the first species ( i = 1 ) is i (cid:126) ∂ ( φ + ˜ ψ ) ∂t = (cid:20) − (cid:126) m ∇ φ − (cid:126) m ∇ ˜ ψ + V φ + V ˜ ψ + U ˆΨ † ˆΨ ˆΨ + U ˆΨ † ˆΨ ˆΨ − µ φ − µ ˜ ψ (cid:105) . (7)The interaction terms in the equation can be written in terms of c -number and fluctua-tion operators as ˆΨ † ˆΨ ˆΨ = | φ | φ + 2 | φ | ˜ ψ + 2 φ ˜ ψ † ˜ ψ + φ ∗ ˜ ψ ˜ ψ + φ ˜ ψ † + ˜ ψ † ˜ ψ ˜ ψ , (8a) ˆΨ † ˆΨ ˆΨ = | φ | φ + | φ | ˜ ψ + φ ∗ ˜ ψ φ + φ ∗ ˜ ψ ˜ ψ + ˜ ψ † φ φ + ˜ ψ † φ ˜ ψ + ˜ ψ † ˜ ψ φ + ˜ ψ † ˜ ψ ˜ ψ . (8b)Since all the atomic fluctuations (quantum and thermal) associated in this theory arewhite noise (cid:104) ˜ ψ i (cid:105) = (cid:104) ˜ ψ i † (cid:105) = 0 . Hence the expectation value of the product of operatorsare (cid:104) ˆΨ † ˆΨ ˆΨ (cid:105) = | φ | φ + φ ∗ (cid:104) ˜ ψ ˜ ψ (cid:105) + 2 φ (cid:104) ˜ ψ † ˜ ψ (cid:105) + (cid:104) ˜ ψ † ˜ ψ ˜ ψ (cid:105) , (9a) (cid:104) ˆΨ † ˆΨ ˆΨ (cid:105) = | φ | φ + φ ∗ (cid:104) ˜ ψ ˜ ψ (cid:105) + φ (cid:104) ˜ ψ † ˜ ψ (cid:105) + φ (cid:104) ˜ ψ † ˜ ψ (cid:105) + (cid:104) ˜ ψ † ˜ ψ ˜ ψ (cid:105) . (9b)Considering that the fluctuations of the two species are uncorrelated (cid:104) ˜ ψ ˜ ψ (cid:105) = (cid:104) ˜ ψ † ˜ ψ (cid:105) =0 , the equation of motion of the condensate of the first species is obtained by taking theaverage of Eq. (7) as i (cid:126) ∂φ ∂t = (cid:20) − (cid:126) m ∇ + V − µ (cid:21) φ + U [ n c + 2˜ n ] φ + U ˜ m φ ∗ + U [ n c + ˜ n ] φ + (cid:104) ˜ ψ † ˜ ψ ˜ ψ (cid:105) + (cid:104) ˜ ψ † ˜ ψ ˜ ψ (cid:105) . (10)Similarly, the equation of motion for the condensate of the second species is i (cid:126) ∂φ ∂t = (cid:20) − (cid:126) m ∇ + V − µ (cid:21) φ + U [ n c + 2˜ n ] φ + U ˜ m φ ∗ + U [ n c + ˜ n ] φ + (cid:104) ˜ ψ † ˜ ψ ˜ ψ (cid:105) + (cid:104) ˜ ψ † ˜ ψ ˜ ψ (cid:105) , (11)6here we have introduced the local densities: n ic ≡ | φ i | , ˜ n i ≡ (cid:104) ˜ ψ † i ˜ ψ i (cid:105) , ˜ m i ≡ (cid:104) ˜ ψ i ˜ ψ i (cid:105) as the condensate, non-condensate, and anomalous densities, respectively. The equa-tion of motion for the non-condensate density of the first species is i (cid:126) ∂ ˜ ψ ∂t = i (cid:126) ∂∂t ( ˆ ψ − φ ) . (12)Using Eq. (7) and Eq. (10) and applying mean-field approximation, ˜ ψ † i ˜ ψ j (cid:39) (cid:104) ˜ ψ † i ˜ ψ j (cid:105) , ˜ ψ i ˜ ψ j (cid:39) (cid:104) ˜ ψ i ˜ ψ j (cid:105) , ˜ ψ † ˜ ψ ˜ ψ (cid:39) (cid:104) ˜ ψ † ˜ ψ (cid:105) ˜ ψ + (cid:104) ˜ ψ ˜ ψ (cid:105) ˜ ψ † , ˜ ψ † ˜ ψ ˜ ψ (cid:39) (cid:104) ˜ ψ † ˜ ψ (cid:105) ˜ ψ , theequation of motion of the fluctuation operator for the first species is i (cid:126) ∂ ˜ ψ ∂t = (cid:18) − (cid:126) m ∇ + V + 2 U ( n c + ˜ n ) − µ + U | φ | + U ˜ n (cid:19) ˜ ψ + U (cid:0) φ + ˜ m (cid:1) ˜ ψ † + U φ φ ∗ ˜ ψ + U φ φ ˜ ψ † . (13)where for the same species i = j , the fluctuation operators are (cid:104) ˜ ψ † i ˜ ψ i (cid:105) = ˜ n i , and (cid:104) ˜ ψ i ˜ ψ i (cid:105) = ˜ m i . However, as mentioned earlier (cid:104) ˜ ψ † i ˜ ψ j (cid:105) = (cid:104) ˜ ψ i ˜ ψ j (cid:105) = 0 .Similarly, the equation of motion of the fluctuation operator of the second speciesis, i (cid:126) ∂ ˜ ψ ∂t = (cid:18) − (cid:126) m ∇ + V + 2 U ( n c + ˜ n ) − µ + U | φ | + U ˜ n (cid:19) ˜ ψ + U (cid:0) φ + ˜ m (cid:1) ˜ ψ † + U φ ∗ φ ˜ ψ + U φ φ ˜ ψ † . (14)For compact notation, we have used the definitions n i = n ic + ˜ n i , and m i = φ i + ˜ m i .The next step is to diagonalise the Hamiltonian matrix and obtain the quasiparticleamplitude functions u s and v s. Incorporating the Bogoliubov transformation, the fluc-tuation operators have the following form ˜ ψ i = (cid:88) j (cid:104) u ij ˆ α j e − iE j t/ (cid:126) − v ∗ ij ˆ α † j e iE j t/ (cid:126) (cid:105) , (15a) ˜ ψ † i = (cid:88) j (cid:104) u ∗ ij ˆ α † j e iE j t/ (cid:126) − v ij ˆ α j e − iE j t/ (cid:126) (cid:105) . (15b)Here, j is the index representing the sequence of quasiparticle excitations. We takethe operators α and α † as common to both the species which is consistent in describingthe coupled multispecies dynamics. Furthermore, this reproduces the standard coupledBdG equations at T = 0 and in the limit a → , the quasiparticle spectra separatesinto two distinct sets: one set for each of the condensates. On substituting Eq. (15) inEqns. (13) and (14) we obtain the BdG equations for TBEC. And, in scaled units theBdG equations are ˆ L u j − U φ v j + U φ ( φ ∗ u j − φ v j ) = E j u j , (16a) ˆ L v j + U φ ∗ u j − U φ ∗ ( φ v j − φ ∗ u j ) = E j v j , (16b) ˆ L u j − U φ v j + U φ ( φ ∗ u j − φ v j ) = E j u j , (16c) ˆ L v j + U φ ∗ u j − U φ ∗ ( φ v j − φ ∗ u j ) = E j v j , (16d)7here ˆ L = (cid:0) ˆ h + 2 U n + U n ) , ˆ L = (cid:0) ˆ h + 2 U n + U n (cid:1) and ˆ L i = − ˆ L i and the quasiparticle amplitudes are normalized as (cid:90) (cid:90) dxdy (cid:88) i ( | u ij ( x, y ) | − | v ij ( x, y ) | = 1 . (17)Under time-independent HFB-Popov approximation for a TBEC, φ i s are the static so-lutions of the coupled generalized GP equations ˆ h φ + U [ n c + 2˜ n ] φ + U n φ = 0 , (18a) ˆ h φ + U [ n c + 2˜ n ] φ + U n φ = 0 . (18b)To solve Eq. (16) we define u ij and v ij ’s as linear combination of N b harmonic oscil-lator eigenstates, u j ( x, y ) = N b (cid:88) κ,l =0 p jκl ϕ κj ( x ) ϕ lj ( y ) , v j ( x, y ) = N b (cid:88) κ,l =0 q jκl ϕ κj ( x ) ϕ lj ( y ) ,u j ( x, y ) = N b (cid:88) κ,l =0 r jκl ϕ κj ( x ) ϕ lj ( y ) , v j ( x, y ) = N b (cid:88) κ,l =0 s jκl ϕ κj ( x ) ϕ lj ( y ) , (19)where ϕ kj s and ϕ lj s are the j th harmonic oscillator eigenstates and p jκl , q jκl , r jκl and s jκl are the coefficients of linear combination. Using this expansion Eq. (16) is reducedto a matrix eigenvalue equation and solved using standard matrix diagonalization algo-rithms. The matrix has a dimension of N b + 1)( N b + 1) × N b + 1)( N b + 1) andis non-Hermitian, non-symmetric and may have complex eigenvalues. Considering theorthogonality of harmonic oscillator basis, the matrix becomes sparse. Due to the N b scaling of the BdG matrix, the matrix size rapidly increases with the basis size, and itis essential to use algorithms capable of large matrix diagonalization. For this reason,we use A RPACK [23]. The eigenvalue spectrum obtained from the diagonalization ofthe matrix has an equal number of positive and negative eigenvalues E j ’s. Using thequasiparticle amplitudes obtained, the number density ˜ n i of the non-condensate atomsis ˜ n i = (cid:88) j { [ | u ij | + | v ij | ] N ( E j ) + | v ij | } , (20)where (cid:104) ˆ α † j ˆ α j (cid:105) = ( e βE j − − ≡ N ( E j ) is the Bose factor of the quasiparticle statewith real and positive energy E j . The coupled Eqns. (16) and (18) are solved itera-tively till the solutions converge to desired accuracy. We use this theory to investigatethe evolution of Goldstone modes and mode energies as a function of the interactionstrengths and temperature. Although, HFB-Popov does have the advantage vis-a-viscalculation of the modes, it is nontrivial to get converged solutions. A measure of phase separation is the overlap integral,
Λ = [ (cid:82)(cid:82) n ( x, y ) n ( x, y ) dxdy ] [ (cid:82)(cid:82) n ( x, y ) dxdy ][ (cid:82)(cid:82) n ( x, y ) dxdy ] . (21)8he TBEC is in the miscible phase when Λ ≈ and signifies complete overlap betweenthe two species when Λ has unit value. The TBEC is completely phase separated when Λ = 0 [24]. The other important measure is the response of the TBEC when subjectedto external perturbations, and one which defines this is the dispersion relation. Todetermine the dispersion relation we compute the root mean square of the wave number k rms of each quasiparticle mode [25, 26] k rms j = (cid:26) (cid:80) i (cid:82) d k k [ | u ij ( k ) | + | v ij ( k ) | ] (cid:80) i (cid:82) d k [ | u ij ( k ) | + | v ij ( k ) | ] (cid:27) / . (22)It is to be noted here that k rms j are defined in terms of the quasiparticle modes corre-sponding to each of the constituent species defined in the k or momentum space throughthe index i = 1 , . It is then essential to compute u ij ( k ) and v ij ( k ) , the Fourier trans-form of the Bogoliubov quasiparticle amplitudes u ij ( x, y ) and v ij ( x, y ) , respectively.Once we have k rms j for all the modes we obtain a discrete dispersion curve. It is to bementioned that to obtain k rms j , we consider 2D Fourier transform and k = ( k x , k y ) andthe integration in Eq. 22 is carried over in 2D Fourier space. The dynamical correlation function or the dynamic structure factor (DSF) charac-terizes the dynamic properties of a quantum many body system and it is a quantity ofconsiderable experimental interest. Unlike other quantum systems where DSF providesinformations ranging from low (characterized by spectrum of collective excitations) tohigh momentum transfer (characterized by momentum distribution), for BECs of di-lute Bose gases DSF is of importance in exploring the domain of high momenta, wherethe response of the system is not affected by its collective features [27]. Rather it isdetermined by the momentum distribution of condensate atoms. In experiments DSF ismeasured by the inelastic light scattering [28] and Bragg spectroscopy [29]. Followingrefs. [27, 30, 31], the dynamic structure factor in terms of j th quasi particle amplitudes u ji ( x, y ) and v ji ( x, y ) for a TBEC is S d ( q x , q y , E ) = (cid:88) j,i (cid:12)(cid:12)(cid:12) (cid:90) (cid:90) dxdy [ u ∗ ji ( x, y )+ v ∗ ji ( x, y )] e i ( xq x + yq y ) / (cid:126) ψ i ( x, y ) (cid:12)(cid:12)(cid:12) δ ( E − (cid:15) j ) , (23)where i corresponds to the species index and for TBEC system i = 1, 2. φ i ( x, y ) is thecondensate order parameter for i th species.Another important measure of the TBEC which is related to the coherence of thesystem is the first-order or the off-diagonal correlation function g (1) i ( x, y, x (cid:48) y (cid:48) ) = (cid:104) ˆΨ † i ( x, y ) ˆΨ i ( x (cid:48) , y (cid:48) ) (cid:105)(cid:104) ˆΨ † i ( x, y ) ˆΨ i ( x, y ) (cid:105)(cid:104) ˆΨ † i ( x (cid:48) , y (cid:48) ) ˆΨ i ( x (cid:48) , y (cid:48) ) (cid:105) , (24)which is also measure of the phase fluctuations. It can also be expressed in terms ofoff-diagonal condensate and noncondensate densities as g (1) i ( x, y, x (cid:48) y (cid:48) ) = n ci ( x, y ; x (cid:48) , y (cid:48) ) + ˜ n i ( x, y ; x (cid:48) , y (cid:48) ) (cid:112) n i ( x, y ) n i ( x (cid:48) , y (cid:48) ) , (25)9here n ci ( x, y ; x (cid:48) , y (cid:48) ) = φ ∗ i ( x, y ) φ i ( x (cid:48) , y (cid:48) ) , (26) ˜ n i ( x, y ; x (cid:48) , y (cid:48) ) = (cid:88) j { [ u ∗ ij ( x, y ) u ij ( x (cid:48) , y (cid:48) ) + v ∗ ij ( x, y ) v ij ( x (cid:48) , y (cid:48) )] N ( E j )+ v ∗ ij ( x, y ) v ij ( x (cid:48) , y (cid:48) ) } (27)At T = 0 , when the entire system is coherent and characterized by the presence of acondensate only, then g (1) i = 1 within the extent of the condensate, whether it is inthe miscible or in the immiscible regime. So, one cannot distinguish between the twophases from the nature of the correlation functions of the individual species. However,at T (cid:54) = 0 , a clear signature of a miscible-immiscible transition of the density profiles isreflected in the form of correlation functions.
3. Details of implementation
As a first step to compute the BdG matrix and derive the BdG equations, we solvethe pair of coupled GP Eqs. (18) using split-step Crank-Nicolson method [32, 33]adapted for binary condensates. The method when implemented with imaginary timepropagation is appropriate to obtain the stationary ground state wave function of theTBEC. To represent the quasiparticle amplitudes u s and v s as a linear combinationof N b direct product states ϕ ( x ) ⊗ ϕ ( y ) as defined in Eq.(19), ϕ ( x ) and ϕ ( y ) areconsidered to be the harmonic oscillator eigenstates. To generate ϕ ( x ) and ϕ ( y ) , westart with the ground ϕ ( x ) and first excited state ϕ ( x ) , and higher excited states aregenerated using the following recurrence relations H n +1 ( x ) = 2 xH n ( x ) − nH n − ( x ) (28) ϕ n ( x ) = (cid:112) /nxϕ n − ( x ) − (cid:114) n − n ϕ n − ( x ) (29)where H n ( x ) is the n th order Hermite polynomial. The computation of basis functionis implemented in the subroutine basis.f90 and stored on a grid. The BdG matrix from the set of BdG Eqs.(16) can be written as E (cid:18) pqrs (cid:19) = (cid:18) BdG BdG BdG BdG (cid:19) (cid:18) pqrs (cid:19) , (30)10here the submatrices in the above matrix equation are defined as BdG = A · · · A N b ... . . . ... A N b · · · A N b N b B · · · B N b ... . . . ... B N b · · · B N b N b E · · · E N b ... . . . ... E N b · · · E N b N b F · · · F N b ... . . . ... F N b · · · F N b N b , (31) BdG = C · · · C N b ... . . . ... C N b · · · C N b N b D · · · D N b ... . . . ... D N b · · · D N b N b G · · · G N b ... . . . ... G N b · · · G N b N b H · · · H N b ... . . . ... H N b · · · H N b N b , (32) BdG = I · · · I N b ... . . . ... I N b · · · I N b N b J · · · J N b ... . . . ... J N b · · · J N b N b M · · · M N b ... . . . ... M N b · · · M N b N b N · · · N N b ... . . . ... N N b · · · N N b N b , (33) BdG = K · · · K N b ... . . . ... K N b · · · K N b N b L · · · L N b ... . . . ... L N b · · · L N b N b O · · · O N b ... . . . ... O N b · · · O N b N b P · · · P N b ... . . . ... P N b · · · P N b N b , (34) pq = p ... p N b N b q ... q N b N b , (35) rs = r ... r N b N b s ... s N b N b . (36)11he BdG matrix is non-Hermitian and non-symmetric with a dimension of N b +1) × N b + 1) , so it can have both real and complex eigenvalues depending on the physicalparameters of the system under study.The eigenvalue spectrum obtained from the diagonalization of the matrix has anequal number of positive and negative eigenvalues E j ’s. From the structure of thematrix elements, we can identify 16 blocks ( A , B , C , D , ..., P ) in the BdG matrix inEq. (30) and in subroutine hfb2d2s.f90 , we compute the matrix elements for theseblocks. In subroutine hfb2d2s.f90 , the blocks A , B , C , D , ..., P correspond to block1, 2 3, 4,..., 16. The elements of each block have the following general expressions A pq = (cid:90) (cid:90) ϕ p ( x, y )[ h + 2 U ( n c + ˜ n ) + U ( n c + ˜ n )] ϕ q ( x, y ) dxdy, B pq = (cid:90) (cid:90) ϕ p ( x, y )[ − U φ ] ϕ q ( x, y ) dxdy, C pq = (cid:90) (cid:90) ϕ p ( x, y )[ U φ φ ∗ ] ϕ q ( x, y ) dxdy, D pq = (cid:90) (cid:90) ϕ p ( x, y )[ − U φ φ ] ϕ q ( x, y ) dxdy, E pq = (cid:90) (cid:90) ϕ p ( x, y )[ U φ ∗ ] ϕ q ( x, y ) dxdy, F pq = − (cid:90) (cid:90) ϕ p ( x, y )[ h + 2 U ( n c + ˜ n ) + U ( n c + ˜ n )] ϕ q ( x, y ) dxdy, G pq = (cid:90) (cid:90) ϕ p ( x, y )[ U φ ∗ φ ∗ ] ϕ q ( x, y ) dxdy, H pq = (cid:90) (cid:90) ϕ p ( x, y )[ − U φ ∗ φ ] ϕ q ( x, y ) dxdy, I pq = − (cid:90) (cid:90) ϕ p ( x, y )[ − U φ ∗ φ ] ϕ q ( x, y ) dxdy, J pq = (cid:90) (cid:90) ϕ p ( x, y )[ − U φ φ ] ϕ q ( x, y ) dxdy, K pq = (cid:90) (cid:90) ϕ p ( x, y )[ h + 2 U ( n c + ˜ n ) + U ( n c + ˜ n )] ϕ q ( x, y ) dxdy, L pq = (cid:90) (cid:90) ϕ p ( x, y )[ − U φ ] ϕ q ( x, y ) dxdy, M pq = (cid:90) (cid:90) ϕ p ( x, y )[ U φ ∗ φ ∗ ] ϕ q ( x, y ) dxdy, N pq = − (cid:90) (cid:90) ϕ p ( x, y )[ U φ φ ∗ ] ϕ q ( x, y ) dxdy, O pq = (cid:90) (cid:90) ϕ p ( x, y )[ − U φ ∗ ] ϕ q ( x, y ) dxdy, P pq = − (cid:90) (cid:90) ϕ p ( x, y )[ h + 2 U ( n c + ˜ n ) + U ( n c + ˜ n )] ϕ q ( x, y ) dxdy (37)The BdG matrix is sparse as the harmonic oscillator basis are orthonormal. So, we use12parse matrix representation to store the matrix, and diagonalized using A RPACK [23]in the subroutine hfbpopov.f . Depending on the parameters, from the diagonaliza-tion we compute the lowest D eigenvalues and corresponding V eigenvectors. u and v From the eigenvectors of the BdG matrix, we compute the quasiparticle amplitudes u and v in the subroutine hfb2d2s.f90 . Considering the array of eigenvectors V ,from Eq.(19) the quasiparticle amplitudes are computed as u j ( x, y ) = N b (cid:88) p,q =0 N b − (cid:88) jj =0 v jj ϕ pj ( x ) ϕ qj ( y ) , (38) v j ( x, y ) = N b (cid:88) p,q =0 2 N b − (cid:88) jj = N b v jj ϕ pj ( x ) ϕ qj ( y ) , (39) u j ( x, y ) = N b (cid:88) p,q =0 3 N b − (cid:88) jj =2 N b v jj ϕ pj ( x ) ϕ qj ( y ) , (40) v j ( x, y ) = N b (cid:88) p,q =0 4 N b − (cid:88) jj =3 N b v jj ϕ pj ( x ) ϕ qj ( y ) . (41)The non degenerate u s and v s are orthonormal. However, to make the degenerate u sand v s orthonormal, we use the Gram Schmidt orthogonalization scheme. Once the eigenvalues ( E j ) of the BdG matrix are obtained after diagonalization,the Bose factor of the j th state in Eq.(42) is N ( E j ) = 1 e βE j − , (42)and the corresponding thermal or non-condensate components are computed using thedefinition of ˜ n i in Eq.(20). As mentioned earlier, for the degenerate states to render the u s and v s orthonormal we use the Gram Schmidt orthogonalization. Among the low-energy collective modes, a few are zero energy, and these are the the Nambu-Goldstone(NG) modes. For TBEC, there exists two NG modes for each of the condensate speciesdue to the breaking of U (1) global gauge symmetry when BEC is formed. These NGmodes do not contribute to ˜ n i , and must be skipped while computing ˜ n i . This is imple-mented through the parameter SKIP = 4 in the main subroutine. In the subroutine hfb2d2s.f90 , we compute the quasi particle amplitudes corresponding to these NGmodes separately.The solutions are iterated until n ic and ˜ n i converge to a predefined accuracy pa-rameter. For T (cid:54) = 0 , the convergence is either very slow due to the thermal fluctuationsor tend to diverge. To accelerate the convergence and ameliorate divergence, we use13he method of Successive under relaxation (SUR)[34], and choose the underrelaxationparameter S = 0 . . The new solution at the k th iteration is then φ new k ( x, y ) = Sφ k ( x, y ) + (1 − S ) φ k − ( x, y ) , (43)where k is the iteration index. To compute ˜ n i we consider the modes with N ( E j ) larger than a threshold value, say − . For parameters relevant to experiments, this isachieved by considering the first 250 or less number of modes.To show the structure of the code, we show a flowchart which describes the howdifferent modules of the code are related..14nitialize Φ i Generate ϕ s SolvecoupledGGPEsEvaluateBdG matrix u s & v s in terms of ϕ sDiagonalizeBdG matrix Use A RPACK
Update Φ i , n ic & ˜ n i with SUR Compute n ic & ˜ n i Converged n ic & ˜ n i ?Stopno yes
4. Description of FACt
This package requires a single input data file input.dat . It consists of ten lines,and description of the input parameters are provided in the contents of the sample file input.dat given below for
Cs - Rb TBEC in miscible regime shown below.15
Where, the parameters are related to various physically significant parameters and theseare as follows:
G011, G022 : s -wave scattering lengths of intraspecies interaction for species 1 and species 2 respectively, G012, G021 : s -wave scattering lengths of interspecies interaction between species 1 and 2, M1, M2 : Mass of species 1 and species 2 respectively,
NUR : Frequency along radial direction, AL : Anisotropy parameter in quasi-2D confinement. ( AL = ω y /ω x ), LAMBDA : Anisotropy parameter to create quasi-2D confinement. (
LAMBDA = ω z /ω x ), TN01, TN02 : Total number of atoms of species 1 and 2 respectively,
SUNDER : Under relaxation parameter to ensure convergence,
NBX, NBY : Number of harmonic oscillator basis taken into account to construct BDG matrix,
NEV, NCV : Number of eigenvalues and eigen vectors ARPACK will print in output file,
TEMPK : Temperature of the system in Kelvin,
SKIP : Number of Goldstone modes,
NUR : Number of HFB Popov self consistent iteration that will ensure convergence,where, the scattering lengths are in the units of Bohr radius ( a ) and the masses are inthe units of amu (atomic mass unit) The above sample input file corresponds to thecase of radially symmetric (AL = 1) Cs - Rb TBEC at zero temperature. Toexamine the effect of anisotropy in the trapping parameter one can consider
AL < 1 (corresponding to ω y (cid:28) ω x , the TBEC is elongated along y axis) or AL > 1 (corre-sponding to ω x (cid:28) ω y , the TBEC is elongated along x axis). In our recent work [17],we have considered the effect of anisotropy in Rb - Rb TBEC at zero temperaturefor
AL > 1 . To make the system quasi-2D a large value of anisotropy parameter alongaxial direction
LAMBDA = 12.5 is chosen so that the condition µ (cid:28) (cid:126) ω z is satisfied.With this condition the atoms are strongly confined along axial ( z ) direction and theyare frozen in the ground state. The size of the harmonic oscillator basis ϕ i chosen toexpand u s and v s is determined by NBX=NBY=55 . This optimal basis size is chosento produce very low ( ∼ O (10 − ) residuals while diagonalising the BdG matrix usingA RPACK . Initially, the total the number of atoms in each species are chosen to be 2000each (
TN01 = TN02 = 2000 ). The under relaxation parameter
SUNDER is keptfixed at 0.1, and number of NG modes to skip is set to 4 (
SKIP = 4 ), this avoidsdivergence associated with the NG modes. The parameter
ITMAX is the maximumnumber of iterations to check the self consistency through HFB-Popov iterations of theBdG equations. 16n addition to the parameters entered from the input.dat , there are other pa-rameters and variables which are defined through modules in the main subroutine hfb main.f90 . The modules
COMM DATA , GPE DATA , and
CN DATA are fromthe original GPE solver code [32, 33]. Solving the HFB-Popov equations require ad-ditional data and variables. For this we introduce two modules
HFB 2D DATA and
ARPK DATA . The former consist of arrays and constants pertaining to the BdG ma-trix and HFB-Popov approximation. These include arrays to store harmonic oscillatorstates ϕ , kinetic energy and potential energy contribution to BdG matrix, etc. The lattermodule has arrays and constants pertaining to A RPACK . Following input files are considered to show a testrun which takes ≈
10 min tocomplete.
On successful completion of computation, the package generates the eigenval-ues and eigen vectors of the BdG matrix. The eigenvalues are stored in data file eigenvalue.out and their corresponding quasiparticle amplitudes are stored in file uv***.dat . Where, *** can take any value between 001 to 200. The details relatedto the computation are given in the output file hfb2d2s.out . Also, the eigen valuesand number of atoms at each HFB Popov iterations are written in hfb2d2s.out . Tocheck for convergence in HFB-Popov iterations, one needs to follow the contents ofoutput file converge.out . The contents of the hfb2d2s.out file for
Cs- Rbat temperature is written below where
Norm1 and
Norm2 check the normaliza-tion,
Psi1ˆ2(0) and
Psi2ˆ2(0) state the density at the center of the confiningpotential for species 1 and 2 respectively. ------------------------------------------------------------------------------Trapping potential, mass and temperature of the quasi-2D TBEC------------------------------------------------------------------------------ALPHA = 1.000, LAMBDA = 12.500NUR = 8.000M1 = 133.000, M2 = 87.000
011 = 280.000, G012 = 100.000G021 = 100.000, G022 = 100.000BETA = Infinity------------------------------------------------------------------------------Derived constants, basis size and spatio-temporal grid information------------------------------------------------------------------------------Oscillator Length = 0.308263D-05MRATIO(MASS1/MASS2) = 1.529No. of basis X = 20No. of basis Y = 20No of spatial points NX = 200No of spatial points NY = 200Spatial step size DX = 0.050000Spatial step size DY = 0.050000Tempotral step size DT = 0.001000Total number of atomsTN01 = 200.00, TN02 = 200.00Number of iterationsNPAS = 5000 NRUN = 1000------------------------------------------------------------------------------iter Norm1 Chem1 Ener
N01 and
N02 correspond to the number of condensate atoms forspecies 1 and 2 respectively. Though the eigenvalues are printed in hfb2d2s.out ,for other detailed computations like the mode evolution as a function of anisotropyand interaction parameters, the energy eigenvalues are also stored in the output file eigenvalue.out . Such data is useful in studies like our previous works [16, 17],where we have shown the mode evolution as a function of various parameters usingthese package. It is to be mentioned that, the energy eigen values, chemical potentialsand total energy of the system, calculated in this package are in units of (cid:126) ω x .
5. Numerical results
In this section, we describe the results from our code in different parameter regimeat zero temperature as well as in finite temperature. At zero temperature, the self-consistent HFB-Popov iterations do not produce significant changes in density pro-files. Since, HFB Popov iterations are computationally expensive and takes time, theresults of zero temperature calculations are provided after single HFB-Popov iterations(
ITMAX = 1). Whereas for finite temperature we consider
ITMAX = 15 which providesrequired convergence.In TBEC, the unique and easily observable effect is phase separation, where thedensity peaks of the component BECs are separate. Alternatively, we can say the mis-cible TBEC phase separates, and enters into immiscible configurations. Numerically,this is quantifiable from the overlap integral Λ as well as the quasi particle amplitudes.In two dimensional (as well as in quasi 2D) systems, the phase separation of TBEC canoccur in two ways. First, the density peaks of the BECs get shifted either along x -axisor along y -axis in x - y plane. This type of phase separation is referred to as side-by-side phase separation. And second possibility arises when one species occupies the core re-gion while the second species surrounds the first one like an annular ring. This type ofphase separated density profile is termed as shell structured density profile. In earlierkind of phase separation, the symmetry of the confining potential is broken where as itis preserved in the later case. In this section we describe the zero temperature condensate density profiles n ic andthe Bogoliubov quasi particle amplitudes u and v in miscible and immiscible regions.In Fig.1, we show the density of condensate atoms n ic ( x, . This figure is obtainedby plotting column 1, 3 and 5 of file den00x.dat for three different inter speciesinteraction strengths. If otherwise mentioned, in all the figures the species 1 and 2correspond to Cs and Rb, respectively. For Fig.1(a) and Fig.1(c) we consider20otal 2000 of atoms where as in Fig.1(b) we consider total 5000 atoms. To obtainequilibrium ground states and avoid metastable states for side by side phase separatedTBEC, it is essential to start the iterations with the initial guess wave functions havingspatially separated peaks. This is implemented in the subroutine initialize.f90 by setting
SHIFT1 = 5.0D0. This also ensures rapid convergence. For other densityconfigurations,
SHIFT1 = 0.0D0 is considered and implies complete overlap of theinitial guess wave functions. -4 0 400.030.06 n c (a) a CsRb =100 Cs Rb -4 0 4 (b) a CsRb =200 -4 0 4 (c) a CsRb =220 x Figure 1: Equilibrium ground state of
Cs- Rb TBEC at zero temperature for three different values ofinterspecies interaction strength (a) a CsRb = 100 a : TBEC is in miscible domain (b) a CsRb = 200 a :TBEC is in shell structures domain and (c) a CsRb = 220 a : TBEC is side by side phase separated. n c ismeasured in units of a − and the spatial coordinates x is measured in units of a osc . From Fig.1(b) it is clear that the TBEC shell structured for the chosen set of pa-rameters, where
Cs BEC is at the core and with the Rb BEC surrounding it. InFig.1(c),
Cs and Rb BECs occupy right and left sides, respectively. Here, the po-sitions of the BECs are not unique, and can interchange depending on the shift in initialguess wave functions. Below we provide content of the input file to corresponding toFig.1(a). input file corresponding to Fig.1(a):
The formation of BEC is associated with the spontaneous symmetry breaking (SSB)of U (1) global gauge. Due to this SSB, in trapped quasi-2D TBEC, the low-energyBdG spectrum has two Goldstone modes for each of the condensate species. In otherwords, the excitation spectrum of the BEC is gapless, and the two lowest energy modes21ith finite energies are the dipole modes. The dipole modes which oscillate out-of-phase with each other are called slosh modes. The in-phase slosh modes with center-of-mass motion are called the Kohn modes and have frequency identical to the naturalfrequency of the harmonic confining potential. Thus the frequency of the Kohn mode isindependent of the type of interactions and interaction strength as well. For this reasonthe getting Kohn mode energy close to 1 serves as an important consistency check ofour FACt package.The Bogoliubov quasi particle amplitudes corresponding to low energy modes areshown in Fig.2, 3 and 4 for miscible, side-by-side and shell structured TBEC respec-tively. -404 (a)E =0.557 u (x,y) (b) v (x,y) (c) u (x,y) (d) v (x,y) -404 (e)E =0.890 (f) (g) (h) -4 0 4-404 (i)E =1.0 -4 0 4 (j) -4 0 4 (k) -4 0 4 (l) -0.100.1-0.1300.06-0.1300.13 x y Figure 2: Quasiparticle amplitudes corresponding to miscible
Cs- Rb TBEC at zero temperature. (a)-(b) show slosh modes for species 1 and (c) - (d) corresponds to those of species 2.(e)-(f) show quadrupolemodes for species 1 and (g) - (h) are those for species 2. (i)-(j) describe the Kohn mode corresponding tospecies 1 and (k) -(l) are those due to species 2. u s and v s are in units of a − and spatial coordinates x and y are in units of a osc . The quasiparticle amplitudes of the selected low-energy modes in the miscible do-main obtained with a CsRb = 100 a are shown in Fig. 2. The images in Fig.2 (a)-(d)correspond to the slosh mode of the system. To obtain the quasiparticle amplitudes, weplot column 3, 4, 5 and 6 of file uv005.dat . In Fig. 2(e)-(h), the quasiparticle ampli-tudes from the file uv010.dat are shown, and these correspond to quadrupole modeof the system. And, the Kohn modes, from the data in the file uv013.dat , are shownin Fig.2(i)-(l). Here, the numerical value in file name uv013.dat indicates thatit is the 13th excited state. For each of the quasiparticle amplitudes the correspondingenergies, taken from the output file eigenvalue.dat , are given in the bottom leftcorner. 22 (a) E =0.000769 u (x,y) (b) v (x,y) (c) u (x,y) (d) v (x,y) -404 (e)E =0.6658 (f) (g) (h) -4 0 4-404 (i)E =0.999 -4 0 4 (j) -4 0 4 (k) -4 0 4 (l) -0.2500.25-0.0500.05-0.2300.13 x y Figure 3: Quasiparticle amplitudes corresponding to side by side phase separated
Cs- Rb TBEC at zerotemperature. (a)-(d) show quasiparticle amplitudes corresponding to NG mode for each of the species.(e)-(h)show those for interface mode for each species. (i)-(l) describe those corresponding to the Kohn mode foreach of the species. Subscript indexes 1 and 2 refer to species 1 and 2 respectively. u s and v s are in units of a − and spatial coordinates x and y are in units of a osc . For the case of side-by-side immiscible phase, with a CsRb = 220 a , the quasipar-ticle amplitudes of low-lying modes are shown in Fig. 3. The images in Fig.3 (a)-(d)correspond to the NG modes of the system which in general resemble n ic , and are basedon the data in the output file uv005.dat . Due to the rotational symmetry breakingassociated with the miscible to side-by-side immiscible phase transition, each specieshas two additional NG modes. The Fig.3(e)-(h) show the quasiparticle amplitudes from uv009.dat , and these correspond to interface mode of the system. In the immiscibledomain the interface modes, as the name suggests, are localized at the interface of thetwo species. The Kohn modes of the system are shown in Fig. 3(i)-(l) and correspondto the data in uv018.dat . 23 (a)E =0.0002 u (x,y) (b) v (x,y) (c) u (x,y) (d) v (x,y) -404 (e)E =0.267 (f) (g) (h) -4 0 4-404 (i)E =0.999 -4 0 4 (j) -4 0 4 (k) -4 0 4 (l) -0.01800.018-0.0500.05-0.2300.23 x y Figure 4: Quasiparticle amplitudes corresponding to shell structured
Cs- Rb TBEC at zero temperature.(a)-(d) show quasiparticle amplitudes corresponding to NG mode for each of the species.(e)-(h) show thosefor interface mode for each species. (i)-(l) describe those corresponding to the Kohn mode for each of thespecies. Like Fig.2 and Fig.3, subscript indexes 1 and 2 refer to species 1 and 2 respectively. u s and v s arein units of a − and spatial coordinates x and y are in units of a osc . For shell structured TBEC, the quasiparticle amplitudes corresponding to NG modes,quadrupole modes and Kohn modes are shown in Fig.4(a)-(d), (e)-(h) and (i)-(l) respec-tively.
For finite temperature computations, solving the HFB-Popov equations require it-erations and we consider
ITMAX = 15 for all the finite temperature computationsreported in this work. The density profiles of n ic corresponding to each HFB-Popoviterations are stored in the file den00x.dat where x runs from 0 to ITMAX. When T (cid:54) = 0 , at each iteration, the number of condensate atoms decreases, whereas thenumber of thermal (non condensate) atoms increases. Fig. 5 shows the equilibriumprofiles of n ic and ˜ n ic for three different temperatures in miscible domain. The plotsin Fig. 5(a) correspond to n ic at T = 0nK , and hence in Fig.5(d) ˜ n ic are negligiblysmall. The plots in Fig. 5(b) and (c) correspond to n ic at T = 5nK and T = 10nK ,respectively. To obtain the plots in the top row, we plotted column 1, column 3 andcolumn 5 file of den00x.dat with column 3 and column 5 multiplied by number ofcondensate atoms N and N (taken from hfb2d2s.out ), respectively. Although,the changes in n ic are not dramatic, there is a large change in ˜ n ic as shown in Fig.5(e)-(f). From Fig. 5, there is a notable feature of ˜ n ic : it has a minimum where n ic hasmaximum value. 24 n c (a) T =0nK Cs Rb (b)
T =5nK (c)
T =10nK -4 0 402.04.0 ˜ n (d) Cs Rb -4 0 4 (e) -4 0 4 (f) x Figure 5: Equilibrium ground state density of
Cs- Rb TBEC in miscible domain for three differentvalues of temperature (a) T = 0nK (b) T = 5nK and (c) T = 10nK . interspecies interaction strength isfixed at a CsRb = 100 a . n c and ˜ n are measured in units of a − and the spatial coordinates x is measuredin units of a osc . For the side by side configuration the density profiles at finite temperature areshown in Fig. 6. Like in the miscible domain, here as well, we observe growth in ˜ n ic with the increase of temperature and thereby lowering the number of condensateatoms. It is to be noted that at the interface of two species, where the n ic are low, ˜ n ic have maximum value. 25 n c (a) T =2nK Cs Rb (b)
T =5nK (c)
T =10nK -4 0 402.04.0 ˜ n (d) Cs Rb -4 0 4 (e) -4 0 4 (f) x Figure 6: Equilibrium ground state density of
Cs- Rb TBEC in immiscible (side by side) domain forthree different values of temperature (a) T = 2nK (b) T = 5nK and (c) T = 10nK . Interspeciesinteraction strength is fixed at a CsRb = 220 a . n c and ˜ n are measured in units of a − and the spatialcoordinates x is measured in units of a osc . Acknowledgments
The example results shown in the paper are based on the computations using theHPC cluster Vikram-100 at Physical Research Laboratory, Ahmedabad.
References [1] A. Griffin, Phys. Rev. B , 9341 (1996).[2] D. A. W. Hutchinson, E. Zaremba, and A. Griffin, Phys. Rev. Lett. , 1842(1997).[3] R. J. Dodd, M. Edwards, C. W. Clark, and K. Burnett, Phys. Rev. A , R32(1998).[4] C. Gies, B. P. van Zyl, S. A. Morgan, and D. A. W. Hutchinson, Phys. Rev. A ,023616 (2004).[5] D. S. Jin, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell, Phys.Rev. Lett. , 764 (1997).[6] S. Gautam and D. Angom, J. Phys. B , 025302 (2011).[7] S. Gautam and D. Angom, Phys. Rev. A , 053616 (2010).[8] M. O. C. Pires and E. J. V. de Passos, Phys. Rev. A , 033606 (2008).269] A. Roy, S. Gautam, and D. Angom, Phys. Rev. A , 013617 (2014).[10] A. Roy and D. Angom, Phys. Rev. A , 011601 (2015).[11] A. Roy and D. Angom, New Journal of Physics , 083007 (2016).[12] C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarru-ell, Science , 301 (2018).[13] D. S. Petrov, Phys. Rev. Lett. , 155302 (2015).[14] D. S. Petrov and G. E. Astrakharchik, Phys. Rev. Lett. , 100401 (2016).[15] G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi,M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, ArXiv e-prints (2017),1710.10890 .[16] S. Pal, A. Roy, and D. Angom, J. Phys B. , 195301 (2017).[17] S. Pal, A. Roy, and D. Angom, J. Phys. B , 085302 (2018).[18] K. Suthar, A. Roy, and D. Angom, Phys. Rev. A , 043615 (2015).[19] K. Suthar and D. Angom, Phys. Rev. A , 063608 (2016).[20] D. S. Petrov, M. Holzmann, and G. V. Shlyapnikov, Phys. Rev. Lett. , 2551(2000).[21] L. Salasnich, A. Parola, and L. Reatto, Phys. Rev. A , 043614 (2002).[22] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Phys. Rev. Lett. , 3745(2000).[23] R. Lehoucq, D. Sorensen, and C. Yang, ARPACK Users’ Guide (Society forIndustrial and Applied Mathematics (Philadelphia), 1998).[24] P. Jain and M. Boninsegni, Phys. Rev. A , 023602 (2011).[25] C. Ticknor, Phys. Rev. A , 053601 (2014).[26] R. M. Wilson, S. Ronen, and J. L. Bohn, Phys. Rev. Lett. , 094501 (2010).[27] F. Zambelli, L. Pitaevskii, D. M. Stamper-Kurn, and S. Stringari, Phys. Rev. A , 063608 (2000).[28] D. M. Stamper-Kurn, A. P. Chikkatur, A. G¨orlitz, S. Inouye, S. Gupta, D. E.Pritchard, and W. Ketterle, Phys. Rev. Lett. , 2876 (1999).[29] J. Stenger, S. Inouye, A. P. Chikkatur, D. M. Stamper-Kurn, D. E. Pritchard, andW. Ketterle, Phys. Rev. Lett. , 4569 (1999).[30] W.-C. Wu and A. Griffin, Phys. Rev. A , 4204 (1996).2731] C. Menotti, M. Kr¨amer, L. Pitaevskii, and S. Stringari, Phys. Rev. A , 053609(2003).[32] P. Muruganandam and S. Adhikari, Comput. Phys. Commun. , 1888 (2009).[33] D. Vudragovi´c, I. Vidanovi´c, A. Balaˇz, P. Muruganandam, and S. K. Adhikari,Comput. Phys. Commun. , 2021 (2012).[34] T. Simula, S. Virtanen, and M. Salomaa, Comput. Phys. Commun. , 396(2001).[35] R. Shankar, Principles of Quantum Mechanics (Springer, 2008).[36] S. Wilson, Journal of Molecular Structure: THEOCHEM , 279 (2001).[37] C.-H. Zhang and H. A. Fertig, Phys. Rev. A , 013601 (2007).[38] P. ¨Ohberg, Phys. Rev. A , 013601 (1999).[39] N. P. Proukakis and B. Jackson, J. Phys. B , 203002 (2008).[40] B. D. Esry, C. H. Greene, J. P. Burke, Jr., and J. L. Bohn, Phys. Rev. Lett. ,3594 (1997).[41] P. ¨Ohberg and S. Stenholm, Phys. Rev. A , 1272 (1998).[42] B. D. Esry and C. H. Greene, Phys. Rev. A , 1457 (1999).[43] Y. N. M. de Escobar, P. G. Mickelson, M. Yan, B. J. DeSalvo, S. B. Nagel, andT. C. Killian, Phys. Rev. Lett. , 200402 (2009).[44] S. Stellmer, M. K. Tey, B. Huang, R. Grimm, and F. Schreck, Phys. Rev. Lett. , 200401 (2009).[45] S. Kraft, F. Vogt, O. Appel, F. Riehle, and U. Sterr, Phys. Rev. Lett.103