Mechanical response functions of finite temperature Bose-Einstein Condensates Stephen Choi1, Vladimir Chernyak2, and Shaul Mukamel1,3 1 Department of Chemistry, University of Rochester, Box 270216, Rochester, New York 14627-0216 2 Corning Incorporated, Process Engineering and Modeling, Corning, New York 14831 3 Department of Physics and Astronomy, University of Rochester, Box 270216, Rochester, New York 14627-0216 UsingtheLiouvillespaceframeworkdevelopedinnonlinearopticswecalculatethelinearresponse functionsandsusceptibilitiesofBose-Einsteincondensates(BEC)subjecttoanarbitrarymechanical 3 force. Distinct signatures of the dynamics of finite temperature BEC are obtained by solving the 0 0 Hartree-Fock-Bogoliubov theory. Numerical simulations of the position dependent linear response 2 functions of one dimensional trapped BEC in thetime and thefrequency domains are presented. n a J 8 ] t f I. INTRODUCTION o s . Since the first experimental observation of atomic Bose-Einstein Condensates (BEC), extensive theoretical and t a experimentalefforthasfocusedonunderstandingtheirproperties[1]. BECpresentsmanypossibleavenuesofresearch, m as it may be viewed in a variety of ways such as mesoscopic “superatom”, gaseous “superfluid”, or as a source of - coherent atoms used in an Atom Laser. In particular, the close analogy between the nonlinear interaction of BEC d matter-waves and photon waves in nonlinear optics gave rise to fascinating atom optics applications [2]. n Optical spectroscopy has been a major tool for studying the properties of matter since the early days of quatum o c physics. And, with the advent of the laser, nonlinear optical spectroscopy has become an important technique for [ studying the properties of matter that are not accessible with incoherent light sources. The primary theoretical tool usedinnonlinearopticalspectroscopytoanalyzethestructureanddynamicprocessesinmanybodyquantumsystems 1 is the linear and higher order optical response functions [3]. In this paper we extend the systematic formalism of v 8 optical response functions to probe the response of trapped, atomic BEC to an arbitrary mechanical force coupled 0 to the atomic density. We calculate the first order (linear) suseptibilities for the condensate, non-condensate density, 1 and non-condensate correlationin the time and the frequency domains. 1 A major difference that arises when the formalism of nonlinear spectroscopy is applied to the mesoscopic BEC is 0 that it is generally not possible to apply the dipole approximation commonly made in the atom-light interaction. It 3 should also be noted that there have been a number of calculations of the linear response of BEC in the past [4]; 0 the response functions calculated here provide a more fundamental Green’s function description of the response to / t a force with arbitrary spatial and temporal profiles. The time domain response functions or their frequency domain a m counterparts, the susceptibilities, provide unique signatures of the dynamics of the system under consideration. The dynamics of zero temperature atomic BEC is commonly described using the time-dependent and time- - d independent Gross-Pitaevskii Equation (GPE). Numerical solutions of the GPE for various properties of zero tem- n perature BEC as well as for atom optics and four wave mixing applications have been reported [5–8]. BEC at finite o temperatures, on the other hand, require sophisticated theories that go beyond the GPE, and various approachesin- c cluding the time-dependent Bogoliubov-deGennes equations [13],the Hartree-Fock-Bogoliubov(HFB) theory [9–12], : v Quantum Kinetic Theory [14–19], and Stochastic methods [21–25] have been employed. i X In order to go beyond the GPE, two or more particle correlations must be taken into account. One of the original contributions was made by Bogoliubov with his introduction of the Bogoliubov transformation [26,27] which shows r a how the condensed state of an interacting homogeneous gas differs from that of the noninteracting gas. That result was extended to inhomogeneous gases by de Gennes [27]. Recently, a time-dependent version of the Bogoliubov-de Gennes equations was employed by Castin and Dum [13] to describe the dynamics of BEC in time dependent traps. Their approach was based on an expansion of the evolution equations for the atomic field operator which is valid if the number of noncondensate particles is small. Their result has been used for analyzing the stability and depletion of a strongly driven BEC [13]. Thetime-independentHFBtheory(TIHFB)hasbeenusedtocalculatesomeoftheimportantequilibriumproperties such as the quasiparticle excitation frequencies and the equilibrium condensate and non-condensate density profiles [11]. The dynamical properties of finite temperature BEC predicted by the TDHFB [12] have not been studied extensively, because the full solution to these equations is computationally prohibitive. The response functions 1 computedhereofferasystematicperturbativeapproachforexploringtheTDHFBdynamics,atanaffordablenumerical cost. The quantum kinetic theory of dilute interacting Bose gas was derived long ago by Kadanoff and Baym [14] in terms of nonequilibrium real-time Green’s functions that parametrize the condensating gas. Their equations were latergeneralizedbyHohenbergandMartintoincludecondensates[15]. Amorecontemporaryversionofthequantum kinetic theory applicable to the experimentally produced BEC has been developed in a series of papers of Gardiner and Zoller [16]. They describe a system composed of interacting a condensate and a noncondensate vapor, where the vaporisdescribedbyaquantumkineticmasterequation,equivalenttoaquantumBoltzmannequationoftheUehling- Uhlenbeckform[20]. Themasterequationwhichdescribesthetransferofenergyandparticlesbetweenthevaporand condensate has been used to describe the formation of BEC. A similar formalism was also put forward by Walser et al. [18] who derived, using a method reminiscent of the classical Bogoliubov-Born-Green-Kirkwood-Yvon(BBYGK) technique, a generalized kinetic theory for the coarse-grained Markovian many-body density operator. Their result incorporates second order collisional processes, and describes irreversible evolution of a condensed bosonic gas of atoms towards thermal equilibrium. Important experimental observations such as the damping of elementary excitations has prompted attempts to develop a theory that works in the hydrodynamic regime. These have led to various versions of Quantum Kinetic theorysuchasthe twofluidhydrodynamicdescriptionoffinite temperatureBEC[17]developedbyZaremba,Nikuni, and Griffin (ZNG). ZNG is a semiclassical Hartree Fock description where one neglects the off-diagonal distribution functions. It combines a non-Hermitian generalized Gross-PitaevskiiEquation for the condensate wave function and the Boltzmann kinetic equation for the noncondensate phase density. This theory was recently shown to properly describe damping [19]. Itwill be helpful to clarifyhowthe variouskinetic theoriesarerelated. The Green’s functionapproachofKadanoff and Baym [14] are equivalent to the equations of Walser et al. [18]. From the equations of Walser et al., the ZNG hydrodynamic equations may be obtained by neglecting the anomalous fluctuations. This greatly simplifies the equations and facilitates the numerical solution. The TDHFB equations are obtained by dropping the second order collisional terms from the kinetic equations of Walser et al. while keeping the anomalous fluctuations. The theory of GardinerandZoller[16]is basedonthe quantumBoltzmannmasterequationsreminiscentofthe quantumstochastic methods used in Quantum Optics; their quantum kinetic master equations also contain the higher order collisional processes described by Walser et al. and ZNG. A number of approximations are common to these kinetic theories: the BornandMarkovapproximation,ergodicassumption,andGaussianinitialreferencedistribution. The validity of theseassumptionswillultimatelybeconfirmedbyexperiments;sofartheexperimentshavesupportedthepredictions of these equations. The stochastic method is another approach used for the description of finite temperature BEC beyond GPE. A description of a Bose gas in thermal equilibrium has been developed using the quantum Monte Carlo techniques, based on Feynman path integral formulation of quantum mechanics [21,22]. A stochastic scheme corresponding to the dynamical evolution of density operator in the positive P representation has been applied to BEC [23,24]. An alternative approach that describes the dynamics of the gas based on a stochastic evolution of Hartree states and avoids some of the instability problems of earlier works was proposed recently [25]. The stochastic methods make the simulation of the exact dynamics of N boson system numerically feasible; this, at present, requires assuming an initial state such as Hartree Fock state which is not very realistic. The TDHFB theory is a self-consistent theory of BEC in the collisionless regime that progresseslogically from the Gross-Pitaevskii Equation by taking into account higher order correlations of noncondensate operators. Although TDHFB does not take into account higher order correlations that are done in the various quantum kinetic theories, the TDHFB equations are valid at temperatures near zero, even down to the zero temperature limit, and are far simpler than the kinetic equations which can only be solved using approximations such as ZNG. Another attractive featureofTDHFB fromapurelypragmaticpointofviewis thatthe Fermionicversionofthe theoryhasalreadybeen well-developed in Nuclear Physics [9]. We therefore work at the TDHFB level in this paper and our approach draws upon the analogywith the time-dependent Hartree-Fock(TDHF) formalismdeveloped for nonlinear opticalresponse of many electron systems [28]. The paper is organized as follows: in Section II we show how to systematically solve the TDHFB equations for externallydrivenfinite temperatureBECbyorderbyorderexpansionofthe dynamicalvariablesintheexternalfield. In Section III we define the n’th order response function in time and frequency domains and calculate the linear (n = 1) response function. Numerical results for the linear response functions and susceptibilities of a 2000 atom condensate in a one dimensional harmonic trap, and its variation with position, time and frequency, are presented in Section IV at zero and finite temperatures. Our main findings are finally summarized in Section V. 2 II. THE TDHFB EQUATIONS A. Equations of motion Ourtheorystartswithatime-dependentmany-bodysecondquantizedHamiltoniandescribingasystemofexternally driven, trapped, structureless bosons with pairwise interactions. Introducing the boson operators aˆ† and aˆ that i i respectively create and annihilate a particle from a basis state i with wave functions φ (r), the Hamiltonian is: i Hˆ =Hˆ +Hˆ′(t) (1) 0 where 1 Hˆ = (H −µ)aˆ†aˆ + V aˆ†aˆ†aˆ aˆ . (2) 0 ij i j 2 ijkl i j k l Xij iXjkm The matrix elements of the single particle Hamiltonian H are given by ij ¯h2 H = d3rφ∗(r) − ∇2+V (r) φ (r), (3) ij Z i (cid:20) 2m trap (cid:21) j where V (r) is the magnetic potential that confines the atoms and µ is the chemical potential. The basis state trap φ (r) is arbitrary; a convenient basis for trapped BEC is the eigenstates of the trap since H is then diagonal. The i ij symmetrized two particle interaction matrix elements in Eq. (2) are 1 V = hij|V|kli+hji|V|kli , (4) ijkl 2 h i where hij|V|kli= d3rd3r′φ∗(r)φ∗(r′)V(r−r′)φ (r′)φ (r), (5) Z i j k l withV(r−r′)beingageneralinteratomicpotential. H′(t)describesthecouplingofanexternalfieldwiththesystem: H′(t)=η E (t)aˆ†aˆ . (6) ij i j Xij η is a bookkeepingexpansionparameter(to be setto 1 atthe endofthe calculation). The matrixelements E (t)are ij given by E = d3rφ∗(r)V (r,t)φ (r), (7) ij i f j Z where V (r,t) denotes a time- and position- dependent external potential that exerts an arbitrary mechanical force f on the system. More general forms of H′(t) could include, for example, terms of the form F (t)aˆ aˆ in addition to the term ij ij i j given in Eq. (6) which couples to atomic density. In this work we shall focuPs on the most experimentally relevant perturbation; for instance, the time-dependent modulation of the trap spring constant which mimics the mechanical force that was recently applied experimentally [29,30]. The dynamics of the system is calculated by deriving equations of motion for the condensate mean field, z ≡haˆ i, i i the non-condensate density ρ ≡ haˆ†aˆ i− haˆ†ihaˆ i, and the non-condensate correlations κ ≡ haˆ aˆ i− haˆ ihaˆ i. ij i j i j ij i j i j Closed equations are derived starting with the Heisenberg equations of motion for the operators, aˆ , aˆ†aˆ , and aˆ aˆ , i i j j j andassumingacoherentmanybodystate. Theresultingmanybodyhierarchyisthentruncatedusingthegeneralized Wick’s theorem for ensemble averages[31–33]: hA i=6 0, (8) i hA A i=hA ihA i+hhA A ii (9) 1 2 1 2 1 2 hA A A i=hA ihA ihA i+hA ihhA A ii+hA ihhA A ii+hA ihhA A ii, (10) 1 2 3 1 2 3 1 2 3 2 1 3 3 1 2 3 and similarly for products involving higher number of operators. A denote boson creation or annihilation operators, i aˆ† or aˆ , and are assumed to be normally ordered. We follow the convention that normal ordered operators are time i i ordered. ThedoubleangularbracketshhA A iidenoteirreducibletwopointcorrelations. Usingthisnotationwehave i j ρ ≡hhaˆ†aˆ ii and κ ≡hhaˆ aˆ ii. ij i j ij i j This procedure yields the TDHFB equations of motion [9,10,12] dz i¯h =[H +ηE(t)]z+H z∗, (11) z z∗ dt dρ i¯h =[h,ρ]−(κ∆∗−∆κ∗)+η[E(t),ρ], (12) dt dκ i¯h =(hκ+κh∗)+(ρ∆+∆ρ∗)+∆+η[E(t),κ] , (13) + dt where [...] denotes the anticommutator. Here, H , H , h, and ∆ are n×n matrices with n being the basis set + z z∗ size used: [H ] =H −µ+ V [z∗z +2ρ ], (14) z i,j ij iklj k l lk Xkl [H ] = V κ , (15) z∗ i,j ijkl kl Xkl h =H −µ+2 V [z∗z +ρ ], (16) ij ij iklj k l lk Xkl ∆ = V [z z +κ ] . (17) ij ijkl k l kl Xkl h is knownas the Hartree-FockHamiltonianand∆ as the pairing field[10]. µis the chemicalpotential introducedin Eq. (2). Eqs. (12) and (13) may also be recast in a more compact matrix form [10]: dγ˜R i¯h =γ˜[Hγ˜,γ˜R]+ηE, (18) dt where H, R, E and γ˜ are 2n×2n matrices defined as follows: h−µ ∆ ρ κ [E(t),ρ] [E(t),κ] H = R= E = + , (19) (cid:18) ∆∗ h∗−µ(cid:19) (cid:18)κ∗ ρ∗+11(cid:19) (cid:18)−[E(t),κ]∗ −[E(t),ρ]∗ (cid:19) + 11 0 γ˜ = , γ˜2 =γ˜, (20) (cid:18) 0 −11(cid:19) B. Solutions of the TDHFB Direct numerical solution of the TDHFB equations [Eqs (11 - 17)] for arbitrary field strength is complicated by their highly nonlinear character. A practical way to proceed is expanding all dynamical variables in powers of η: z (t)=z(0)+ηz(1)(t)+η2z(2)(t)+··· (21) i i i i ρ (t)=ρ(0)+ηρ(1)(t)+η2ρ(2)(t)+··· (22) ij ij ij ij κ (t)=κ(0)+ηκ(1)(t)+η2κ(2)(t)+··· (23) ij ij ij ij and upon substituting these expansions in Eqs. (11 - 17), collecting all terms to a given power in η. This results in a hierarchy of equations such that the equation of motion for the n’th order solution α(n), where α denotes the variablesz,ρorκ,isexpressedasafunctionof0,...(n−1)thordersolutions,α(0),α(1),...,α(n−1). Whilethezero’th 4 order equation is nonlinear, each of the higher order equations are linear; the n’th order solutions, α(n) are therefore obtained by solving the sequence of linear equations. This is analogous to the TDHF [28] or time dependent density functional algorithms used for many electron systems. Findingthezero’thordersolution{z(0),ρ(0),κ(0)}shouldbethefirststepinsolvingtheTDHFB.Thisdescribesthe system at equilibrium in the absence of the external driving field, and requires the solution of the TIHFB equations: H(0)z(0)+H(0)z(0)∗ =0 (24) z z∗ [h(0),ρ(0)]−(κ(0)∆(0)∗−∆(0)κ(0)∗)=0 (25) (h(0)κ(0)+κ(0)h(0)∗)+(ρ(0)∆(0)+∆(0)ρ(0)∗)+∆(0) =0 (26) Here, H(0) , H(0), h(0), and ∆(0) are n×n matrices defined in Eq. (14-17) with the variables {z,ρ,κ} replaced by z z∗ {z(0),ρ(0),κ(0)}. It follows from Eq. (18) that Eqs. (25) and (26) can be written in the compact form: γ˜[H(0)γ˜,γ˜R(0)]=0. (27) Self-consistent numerical methods for solving the TIHFB are given in the Appendix [10]. So far we have worked in the trap basis since this enables us to maintain a general form for the interatomic interactions V and makes the numerical solution feasible. However, it may be more interesting to recast the ijkl solution in real space α(n)(r,t) (α=z, ρ, or κ) by transforming the solution to TDHFB in the trap basis α(n)(t): z(n)(r,t)= z(n)(t)φ (r), (28) j j Xj ρ(n)(r,t)= ρ(n)(t)φ∗(r)φ (r), (29) ij i j Xij κ(n)(r,t)= κ(n)(t)φ (r)φ (r). (30) ij i j Xij In general, real space non-condensate density and non-condensate correlations are nonlocal functions of two spatial points ρ(r′,r) and κ(r′,r). We only computed these quantities for r = r′ in this paper since these are the most physically accesible. Measuring these quantities with r 6= r′ involves observing atomic correlations which is much more difficult than photon correlations. C. Liouville space representation We next introduce the Liouville space notation [3,34] that will be used in the folllowing sections. One well-known example for which the Liouville space formalism is used is in the optical Bloch equations, where the 2×2 density matrix is recast as a 4 component vector, and the Liouvillian is written accordingly as a 4×4 matrix superoperator. We rearrange the TIHFB equations Eqs. (24-26) by writing ρ and κ as vectors in Liouville space, and introduce ij ij the following set of n2×n2 matrices (superoperators): H(−) =h(0)δ −h(0)δ , (31) ij,mn im jn nj im (+) (0) (0) H =h δ +h δ +V , (32) ij,mn im jn nj im ijmn D =∆(0)δ , (33) ij,mn im jn D∆ =−∆∗(0)δ , (34) ij,mn nj im where the n×n matrices h and ∆ were defined in Eqs. (16-17). We further define the n2×1 matrices Λκ = V z z . (35) ij ijkl k l Xkl Using this notation, the TIHFB assume the form: 5 H(0) H(0) 0 0 0 0 ~z(0) 0 z z∗ H(0)∗ H(0)∗ 0 0 0 0  ~z(0)∗   0  z∗ z  0000 0000 −HDD0(−∗∆)∗ HDD0(∆+∗) H−(DD0−∆)∗ HD(D0∆+∗)∗ ~κρ~~κρ~((((0000))))∗∗ +ΛΛ~~00κκ∗ =0. (36)      Hereafter, we shall denote the zero’th order solution in Liouville space notation as ψ~(0)(t) ≡ [~z(0),~z(0)∗,ρ~(0),~κ(0),ρ~(0)∗,~κ(0)∗]T and refer to the 2n(2n+1) by 2n(2n+1) matrix multiplying ψ~(0) as the zero’th order Liouville operator, L . 0 Substituting the expansion for z , ρ , and κ [Eqs. (21 - 23)] into Eqs. (11-13), we obtain to first order in η: i ij ij dψ~(1)(t) i¯h =Lψ~(1)(t)+ζ(t). (37) dt Here ψ~(1)(t) = [~z(1),~z(1)∗,ρ~(1),~κ(1),ρ~(1)∗,~κ(1)∗]T i.e. a 2n(2n+1)×1 vector with the variables in first order as its components. Also, L≡L +L (38) 0 1 where L is the total Liouvillian, L was introduced in Eq. (36), and L and ζ(t), obtained by equating all first order 0 1 terms in η are given in Appendix B. L is given by the sum of the original TIHFB matrix L plus a perturbation L . 0 1 This perturbation induces a shift in the excitation frequencies, as will be shown below. Adopting Liouville space notation, the position-dependent nth order solution ψ~(n)(r,t) can be defined using the relations Eq. (28 - 30) and introducing a 2n(2n+1)×2n(2n+1) square matrix Υ˜(r): ψ~(n)(r,t)≡Υ˜(r)ψ~(n)(t) (39) where Υ˜(r)=diag φ˜(r),φ˜∗(r),Φ (r),Φ (r),Φ∗(r),Φ∗(r) . (40) ρ κ ρ κ h i Here“diag[···]”denotesthatΥ˜(r)isablockdiagonalsquarematrixmadeofn×nblocksφ˜(r),φ˜∗(r)andn2×n2blocks Φ (r),Φ (r),Φ∗(r),Φ∗(r). φ˜(r)isadiagonalmatrixwiththei-thdiagonalelementgivenbythebasisstatesφ (r),and ρ κ ρ κ i Φ (r) and Φ (r) are also diagonal matrices whose ij’th diagonal element are given by [Φ (r)] = φ∗(r)φ (r), and ρ κ ρ ij,ij i j [Φ (r)] = φ (r)φ (r) respectively. The real space variables z(n)(r,t), ρ(n)(r,t), and κ(n)(r,t) are finally obtained κ ij,ij i j by summing over the appropriate elements of the vector ψ~(n)(r,t): n 2n+n2 2n+2n2 z(n)(r,t)= ψ~(n)(r,t), ρ(n)(r,t)= ψ~(n)(r,t), κ(n)(r,t)= ψ~(n)(r,t). (41) i i i Xi=1 i=X2n+1 i=2nX+n2+1 III. THE RESPONSE FUNCTIONS The n’th order response function K(n)(t,t ,...t ,r,r ,...r ) is defined by the relation: αρ 1 n 1 n α(n)(r,t)= K(n)(t,t ,...t ,r,r ,...r ),V (r ,t )···V (r ,t )dt ···dt dr ···d3r , (42) Z Z αρ 1 n 1 n f 1 1 f n n 1 n 1 n whereα(n)(r,t) withα=z,ρ,orκ areposition-dependentn’thordersolutions. TheρsubscriptofK(1) indicatethat αρ this is the response to an external field that couple to the atomic density. The n’th order susceptibility is defined as the Fourier transform of the response function to the frequency domain: ∞ ∞ χ(n)(Ω,Ω ,...,Ω ,r,r ,...r )= dt dt ···dt K(n)(t,t ,...,t ,r,r ,...r ) αρ 1 n 1 n 1 n αρ 1 n 1 n Z Z 0 0 ×exp(iΩt+iΩ t +···+iΩ t ). (43) 1 1 n n 6 A. The time-domain linear response function The solution of the matrix equation Eq. (37) is 1 t i ψ~(1)(t)= exp − L(t−t′) ζ(t′)dt′. (44) i¯hZ (cid:20) ¯h (cid:21) 0 The corresponding position-dependent solution can then be written: 1 t ψ~(1)(r,t)= dr′ Υ˜(r)U(t−t′)Φ˜(r′)ψ~(0)V (r′,t′)dt′ (45) f i¯hZ Z 0 1 t = dr′ K~(1)(t−t′,r,r′)V (r′,t′)dt′, (46) i¯hZ Z ψ~ f 0 where Υ˜(r) is as given in Eq. (40) and we have defined the 2n(2n+1)×2n(2n+1) matrices i U(t−t′)=θ(t−t′)exp − L(t−t′) , (47) (cid:20) ¯h (cid:21) Φ˜(r)=diag Φ(r),Φ∗(r),Φ(−)(r),Φ(+)(r),Φ(−)∗(r),Φ(+)∗(r) , (48) h i with θ(t−t′) being the Heaviside function and as discussed above,the notation “diag[···]” is used to denote Φ˜(r′) as a 2n(2n+1)×2n(2n+1) block diagonal square matrix with the blocks consisting of n×n square matrices having the ith row and jth column given by [Φ(r)] =φ∗(r)φ (r), (49) ij i j and n2×n2 square matrices Φ(±)(r) =φ∗(r)φ (r)δ ±φ∗(r)φ (r)δ . (50) i m jn n j im h iij,mn The function K~(1)(t−t ,r,r ) of Eq. (46), ψ~ 1 1 K~(1)(t−t ,r,r )≡Υ˜(r)U(t−t )Φ˜(r )ψ~(0), (51) ψ~ 1 1 1 1 may be viewed as the linear response function for the position dependent vector ψ~(1)(r,t) i.e. for all the variables z(1), ρ(1) and κ(1) in the trap basis. The real space response functions for the condensate z, non-condensate density ρ,andthe non-condensatecorrelationκ arethereforegivenbysumming overappropriateindices inK~(1)(t−t ,r,r ), ψ~ 1 1 using the relations Eq. (28-30): n K(1)(t−t ,r,r )= K~(1)(t−t ,r,r ), (52) zρ 1 1 ψ~i 1 1 Xi=1 2n+n2 K(1)(t−t ,r,r )= K~(1)(t−t ,r,r ), (53) ρρ 1 1 ψ~i 1 1 i=X2n+1 2n+2n2 K(1)(t−t ,r,r )= K~(1)(t−t ,r,r ). (54) κρ 1 1 ψ~i 1 1 i=2nX+n2+1 To analyze the physical significance of the response functions it will be useful to expand them in the basis of the eigenvectors ξ~ of matrix L ν Lξ~ =ω ξ~ , ν =1,2,...2n(2n+1). (55) ν ν ν We define the Green’s function 7 i G (t−t′)=θ(t−t′)exp − ω (t−t′) , (56) ν ν (cid:20) ¯h (cid:21) and further introduce µ , η (r), and δ (r) as the expansion coefficients of the following vectors in the basis of ν ν ν eigenvectors ξ~ ν 2n(2n+1) 2n(2n+1) 2n(2n+1) ψ~(0) = µ ξ~ ; Φ˜(r)ξ~ = η (r)ξ~ ; Υ˜(r)ξ~ = δ (r)ξ~ . (57) ν ν ν ν ν ν ν ν Xν=1 νX=1 νX=1 We then have K~(1)(t−t ,r,r )= K(1)(t,t′,r,r )ξ~ (58) ψ~ 1 1 νρ 1 ν Xν where Kν(1ρ)(t,t′,r,r1)= δν(r)ην′(r1)µν′′Gν′(t−t′). (59) νX′,ν′′ Here,G (t−t′)isdefinedinEq. (56),µ ,η (r), andδ (r)aredefinedasthe expansioncoefficientsinEq. (57),while ν ν ν ν the matrices Φ˜(r) and Υ˜(r) are given by Eqs. (40) and (48). Eqs. (58-59) express the linear response function Eq. (51) as an expansion in quasiparticle modes. B. The frequency domain response function The linear susceptibility is defined as the Fourier Transform of the response function to the frequency domain: ∞ ∞ χ~(1)(Ω,Ω ,r,r )= dt dt K~(1)(t−t ,r,r )exp(iΩt+iΩ t ). (60) 1 1 1 1 1 1 1 Z Z 0 0 It is possible to change variables t−t → τ and define χ~(1)(Ω,r,r ) since K~(1)(t−t ,r,r ) only depends on the 1 1 1 1 time difference, t−t . Below, we shall keep the notation K~(1)(Ω,Ω ,r,r ) i.e. a function of both frequency of the 1 1 1 signal, Ω, and frequency of perturbation, Ω , rather than K~(1)(Ω,r,r ), in line with Bloembergen’s notation [3]. 1 1 In order to evaluate the Fourier Transform, we note that it follows from causality 1 ∞ 1 U(t)=− dω exp(−iωt) (61) 2πiZ ω−L+iǫ −∞ ∞ = dωU(ω)exp(−iωt). (62) Z −∞ Substituting Eq. (62) into the expression for K~(1)(t−t ,r,r ) Eq. (51), and taking the Fourier Transform, we ψ~ 1 1 obtain a vector: 1 ∞ ∞ χ~(1)(Ω,Ω ,r,r )= dω dtdt Υ˜(r)U(ω)exp(−iωt+iωt )exp(iΩt+iΩ t )Φ˜(r )ψ~(0) (63) ψ~ 1 1 2πiZ Z 1 1 1 1 1 −∞ 0 1 ∞ = dωΥ˜(r)U(ω)Φ˜(r )ψ~(0)δ(Ω−ω)δ(Ω +ω). (64) 1 1 2πiZ −∞ This implies that Ω=−Ω , and we have the susceptibility in vector form 1 χ~(1)(−Ω ;Ω ,r,r )=Υ˜(r)U(Ω )Φ˜(r )ψ~(0). (65) ψ~ 1 1 1 1 1 Thez,ρandκsusceptibilitiesareobtainedbysummingovertheappropriateelementsinthevectorχ~(1)(−Ω;Ω,r,r ): ψ~ 1 8 n χ(1)(−Ω;Ω,r,r )= χ~(1)(−Ω;Ω,r,r ) (66) zρ 1 ψ~i 1 Xi=1 2n+n2 χ(1)(−Ω;Ω,r,r )= χ~(1)(−Ω;Ω,r,r ) (67) ρρ 1 ψ~i 1 i=X2n+1 2n+2n2 χ(1)(−Ω;Ω,r,r )= χ~(1)(−Ω;Ω,r,r ). (68) κρ 1 ψ~i 1 i=2nX+n2+1 Transforming the trap basis to the basis of eigenstates ξ~ as before, one has: ν χ~(1)(−Ω;Ω,r,r )= K(1)(−Ω;Ωr,r )ξ~ , (69) ψ~ 1 νρ , 1 ν Xν where K(1)(−Ω;Ω,r,r )= δν(r)ην′(r1)µν′′, (70) νρ 1 νX′ν′′ Ω−ων′ +iǫ and δ (r), η (r ), and µ are as defined in Eq.(57). ν ν 1 ν As was done in Section III A for the time domain, we have recast the linear susceptibility in two forms: matrix form of Eqs. (65) and an expansion in modes of Eq. (69-70). IV. NUMERICAL SIMULATIONS FOR CONTACT POTENTIAL A. Zero’th order solution (TIHFB) and the frequency shifts So far, all our results hold for a general pairwise interatomic interaction potential. In the following numerical calculations, we approximate the interatomic potential V(r−r′) in Eq. (5) by a contact potential, as is standard in BEC applications: 4π¯h2a V(r−r′)→U δ(r−r′), U = , (71) 0 0 m whereaisthes-wavescatteringlengthandmistheatomicmass. Thisapproximationmaybejustifiedsincethewave functions at ultracoldtemperatures havevery long wavelengthscomparedto the rangeof interatomic potential. This implies that details of the interatomic potential become unimportant and the potential may be approximated by a contact potential. The tetradic matrices V are then simply given by: ijkl 4π¯h2a V = φ∗(r)φ∗(r)φ (r)φ (r)dr. (72) ijkl m Z i j k l We assume a 2000 atom one dimensional condensate in a harmonic trap. The parameters used for our numerical calculationsare: U = 4πh¯2a =0.01,andtemperatures0h¯ω /k and10h¯ω /k whereω isthetrapfrequency 0 m trap B trap B trap and k is the Boltzmann constant. We used 256 grid points for position and the basis set of n = 5 states. We keep B the trap units throughout. We have followed the prescription of Griffin for solving the TIHFB for ψ~(0) in terms of the Bogoliubov-de Gennes EquationswhichareobtainedfromEq. (19)bytransformingtorealspaceandusingthecontactinteratomicpotential [11]. In the Appendix we show the equivalence between the Bogoliubov-de Gennes Equations and the matrix H(0), and summarize the numerical procedure. The eigenvalues of the non-Hermitian matrix L required for computing the response functions were calculated using the Arnoldi algorithm [35]. TheeigenvaluesoftheLiouvillianaretrasitionfrequenciesratherthanstateenergies. Thefrequenciescomeinpairs ofpositiveandnegativefrequencies;thisindicatesthattheLiouvillianmaybemappedontoaharmonicoscillatorspace for whichthere arealwayspositive andnegativefrequencysolutions. The eigenvaluesofL≡L +L areshifted with 0 1 respecttothecorrespondingeigenvaluesofL (theTIHFBequations). Welistsomeoftherepresentativeeigenvalues, 0 namely the lowest few positive eigenvalues of L and L in Table I for both zero and nonzero temperatures. Similar 0 9 frequency shifts were noted by Giorgini [4]. Physically it is easy to understand how the TDHFB frequencies may be shifted from the TIHFB: a dynamical system contains the effect of the interacting condensate and non-condensate atomswhich,bydefinition,isnotpresentintheequilibriumsystem. Thisisanalogoustoopticalexcitationsoffermions where the time-dependent Hartree-Fock theory shows excitonic shifts which are lacking by the time-independent Hartree-Fock states [28]. B. Time domain response The linear response functions were calculated by substituting the numerical solution ψ~(0) evaluated at zero and finite temperatures into Eq. (51) together with Eqs. (52-54). Eq. (51) is a matrix multiplication of 2n(2n+1)×1 vector ψ~(0) with 2n(2n+1)×2n(2n+1) matrices Υ˜, U(t), and Φ˜ which are defined in Eqs. (48) and (40). Φ˜ and Υ˜ are constructed in terms of the harmonic oscillator basis states which are calculated numerically from the recursive formulathatinvolvestheGaussianfunctionmultiplyingtheHermitepolynomials[36]. ThematrixU(t)wascalculated using a MATLAB function that uses the Pad´e approximation for matrix exponentiation [37]. We first present the dependence of the linear response function on r and r′ at fixed times t − t′. This gives a snapshot of the position-dependent correlations across the condensate. Such dependence is important since the experimentally produced condensates are mesoscopic in size; in contrast, the dipole approximationusually applies in optical spectroscopy and consequently the spatial dependence of the response is not observable. In Figs 1 we display the real space response functions at times t − t′ = 0/ω ,7.2/ω ,15.7/ω at zero trap trap trap temperature. Fig. 2showsthecorrespondingfinitetemperatureresults. Physically,t′ andr′ arerespectivelythetime andpostionatwhich the externalperturbationis appliedand t andr arethe correspondingcoordinatesat whichthe measurementis made. We found that for longer times t−t′ >5π/ω , the plot maintains generally the same shape trap asthe thirdcolumnoffigures;it maythereforebe possible to experimentallyobservethe correlationsatlongertimes. At zero temperature, the correlationattains this stable shape faster than at finte temperature. To model a uniform perturbation applied acrossthe condensate, we have integratedthe zero temperature response function over r′ and plotted its absolute value i.e. | K(1)(r,r′,t−t )dr′| vs. time t−t′ in Fig. 3. The response 1 functionsKzρ,Kρρ,Kκρ growexponentiallywithinarRelativelyshorttimespanofaround5π/ωtrap,sothatthedetails ofthe structureforearliertimesupto∼2π/ω arenotvisibleinthe plot. Thisrapidgrowthisshownmoreclearly trap in the right hand column which depicts the response functions integrated over both r and r′, i.e. K(1)(t−t ). The 1 figure shows the real and imaginary parts as well as the absolute values of K(1)(t−t ). The response function grows 1 rapidly by around an order of magnitude over the time scale ∼ 3π/ω . The response function keeps growing over trap time since TDHFB does nothaveany dissipativeterm; this shouldbe areasonablemodelfor BECin the collisionless regime. It also shows that, even in the absence of a dissipative term, it takes some time after the initial impulse at t′ = 0 before the effect of the force is reflected appreciably in the response functions. From the plots we note that a mechanical force applied on the condensate can be seen to “generate” noncondensate atoms and anomalous correlations even at zero temperature. ThecalculationsofFig. 3arerepeatedatfinitetemperatureinFig. 4. Themaindifferenceisthesmallermagnitude of K , K and K . The fact that the BEC is less responsive at finite temperature may be attributed to the fact zρ ρρ κρ that the condensate to non-condensate interaction is greater at finite temperatures where additional collisions shield the effect of the applied perturbation. C. Frequency domain response The susceptibilities were computed using Eq. (65) in conjunction with Eqs. (66-68) Eq. (65) is a matrix multi- plication of 2n(2n+1)×1 vector ψ~(0) with 2n(2n+1)×2n(2n+1) matrices Υ˜, U(ω), and Φ˜. The matrix U(ω) is calculated as follows: 1 ξ ζ† U(ω)= = ν ν (73) ω−L+iǫ ω−ω +iǫ Xν ν where ξ is the right eigenvalue of L with eigenvalues ω (Lξ = ω ξ ) and ζ are the left eigenvectors such that ν ν ν ν ν ν νξνζν† =11. The eigenvalues ων of L were calculated using the Arnoldi algorithm [35]. PTo clearly display the resonance structure, we present in Fig. 5 the absolute value of zero and finite temperature linearresponse integratedoverr andr inthe frequency domain,| χ(1)(−Ω,Ω,r,r′)drdr′|, α=z,ρ,κona logarith- 1 αρ mic scale. The eigenvalues of L and hence the resonant frequencieRs come in positive and negative pairs; we present 10