Bogoliubov spectrum and Bragg spectroscopy of elongated Bose-Einstein condensates

The behaviour of the momentum transferred to a trapped Bose-Einstein condensate by a two-photon Bragg pulse reflects the structure of the underlying Bogoliubov spectrum. In elongated condensates, axial phonons with different numbers of radial nodes give rise to a multibranch spectrum which can be resolved in Bragg spectroscopy, as shown by Steinhauer et al (2003 Phys. Rev. Lett. 90 060404). Here we present a detailed theoretical analysis of this process. We calculate the momentum transferred by numerically solving the time-dependent Gross-Pitaevskii (GP) equation. In the case of a cylindrical condensate, we compare the results with those obtained by linearizing the GP equation and using a quasiparticle projection method. This analysis shows how the axial-phonon branches affect the momentum transfer, in agreement with our previous interpretation of the observed data. We also discuss the applicability of this type of spectroscopy to typical available condensates, as well as the role of nonlinear effects.


Introduction
In a recent paper [1] we showed that interesting features emerge in the momentum transferred to an elongated Bose-Einstein condensate by a two-photon Bragg pulse when the duration of the pulse is long enough. On the basis of numerical simulations with the Gross-Pitaevskii (GP) equation, we interpreted those features as due to the structure of the underlying Bogoliubov spectrum: the Bragg pulse can excite axial phonons with different numbers of radial nodes, each one having its own dispersion law. Since the energy difference between these branches is of the order of the radial trapping frequency, they can be resolved only by Bragg pulses with duration comparable with the radial trapping period, as observed in [1]. In this work, we present a more detailed theoretical analysis that supports our previous interpretation and gives further information about the role played by Bogoliubov excitations in Bragg spectroscopy.
A useful insight into this problem is obtained by studying the case of an infinite condensate, unbound along z and harmonically trapped in the radial direction ρ. We show that the response of such a cylindrical condensate retains all the relevant properties needed to interpret the observed behaviour of a finite elongated condensate. In addition to the numerical solution of the full time-dependent GP equation, we explicitly determine the time evolution of the order parameter in the linear (small amplitude) limit, using the quasiparticle projection method of [2,3]. The cylindrical geometry allows us to simplify the calculations and, more important, to make the connection between the momentum transferred and the dynamic structure factor more direct than for a finite condensate. By comparing the predictions of this approach with the results of the numerical integration of the time-dependent GP equation, one can also distinguish the different effects of linear and nonlinear dynamics. This analysis of cylindrical condensates is presented in the first sections of this paper, while section 5 is devoted to the behaviour of finite elongated condensates, such as those in the experiments of [1,4,5]. Finally, in the appendix we give a short description of the algorithms used to solve the relevant equations.

Bogoliubov spectrum of a cylindrical condensate
First we investigate the excitation spectrum of a cylindrical condensate. We use the meanfield GP theory at zero temperature, assuming the condensate to be dilute and cold enough to neglect thermal and beyond mean-field contributions. The starting point is the time-dependent 54.3 GP equation for the order parameter (r, t) of a condensate of N bosonic atoms subject to an external potential V (r) [6]: where g = 4πh 2 a/m and a is the s-wave scattering length, that we assume to be positive. The order parameter is normalized according to dr | | 2 = N. What is usually named Bogoliubov spectrum is the result of an expansion of in terms of a quasiparticles basis in the form [7]- [9]: where the c j are constants and µ is the chemical potential. The quasiparticle amplitudes u and v obey the following orthogonality and symmetry relations: This expansion can be inserted into equation (1). At zero order in u and v, one gets the stationary GP equation: Among the solutions of this equation one finds the ground state order parameter, which can be chosen real and equal to the square root of the particle density [10]. The next order gives the coupled equations where H 0 = −(h 2 /2m)∇ 2 + V (r). The solutions of these equations provide the spectrum of the excited states in the linear (small amplitude) regime. Now we treat the specific case of a cylindrical condensate, which is unbound along z and harmonically trapped along the radial direction, ρ = [x 2 + y 2 ] 1/2 . Let us write the trapping potential as V (r) = V T (ρ) = (1/2)mω 2 ρ ρ 2 . The ground state order parameter is uniform along z, while the stationary GP equation for the radial part becomes where N/L is the number of bosons per unit length and the function φ 0 is chosen to be real and subject to the normalization condition 2π ∞ 0 dρ ρφ 2 0 = 1. In this geometry, the quasiparticle amplitudes u and v can also be factorized. We are actually interested in axially symmetric states (azimuthal angular momentum equal to zero), since they are the only ones excited by a Bragg process in which the momentum is imparted along the z-axis. These states, characterized by the axial wavevector, k, and the number of radial nodes, n, can be written in terms of new functions u n,k (ρ) and v n,k (ρ) in the form

54.4
The Bogoliubov equations (5)-(6) then becomē Numerical solutions of the Bogoliubov equations for cylindrical condensates have already been obtained in the context of theoretical studies of the Landau critical velocity [11], the Landau damping of low-energy collective oscillations [12], and the behaviour of solitary waves [13].
Here we solve equations (9) and (10) in order to use the functions u and v as a basis for a time-dependent calculation, as explained in the next section.
It is worth noticing that expressing all lengths and energies in units of the harmonic oscillator length, a ρ = [h/(mω ρ )] 1/2 , and energy,hω ρ , respectively, the solutions of GP and Bogoliubov equations for the cylindrical condensate scale with a single parameter, aN/L. The same parameter can be expressed in terms of the chemical potential in the Thomas-Fermi limit, µ TF . One has in fact, η ≡ µ TF /(hω ρ ) = [4aN/L] 1/2 [11] 4 . We perform calculations for different values of η, which are representative of available elongated condensates. The lowest value that we consider is η = 9.4, which simulates the condensate of Davidson and co-workers [1,5]; the largest value is η = 70, for the condensate of [14]. We also use an intermediate value, η = 26.5, which corresponds to the condensate of the calculations in [15,16].
In figures 1(a) and (b) we show the frequency ω n,k of the Bogoliubov branches as a function of k, for η = 9.4 and 26.5 respectively. The lowest branch starts linearly at low k and then becomes quadratic for large k, as expected for the textbook Bogoliubov spectrum of a uniform condensate. The transition between the two regimes occurs at k of the order of where ξ is the healing length. The slope of this branch has been the object of the investigation of [11]. The second branch starts at ω 2,0 = 2ω ρ , where it reduces to the transverse monopole (breathing) mode; its frequency is model independent as a consequence of a scaling property of the GP equation in two dimensions [17]. For k ξ −1 the dispersion law of the lowest branches can be obtained, with great accuracy, also within a hydrodynamic approach [18].
Our purpose is now to explore the effects of the discrete multibranch spectrum of figures 1(a) and (b) in the response of the condensate to a Bragg pulse.

Time-dependent quasiparticle amplitudes and momentum transferred
Let us suppose that a Bragg pulse is switched on at a time t = 0. In the experiments, this is obtained by illuminating the condensate with two laser beams with wavevectors k 1 and k 2 and frequencies ω 1 and ω 2 . The action of the stimulated light-scattering from atoms of the condensate is equivalent to a moving optical potential with wavevector q = k 1 −k 2 and frequency ω = ω 1 − ω 2 . Within the GP theory, this is accounted for by adding an extra term to the external potential, that becomes [19]  where θ(t) is the Heaviside step function and we have chosen q along z. We calculate the response of the condensate in two ways: (i) by using the quasiparticle projection method of [2,3]; (ii) by numerically solving the full time-dependent GP equation (1) for given values of η, V B , q and ω.
The latter case will be discussed in the next section. Here we briefly sketch the quasiparticle projection method and its application to infinite cylindrical condensates. The starting point is a generalization of the Bogoliubov expansion (2) that allows the population of quasiparticle states to vary with time [2]. The simplest way is to expand the 54.6 order parameter in terms of a quasiparticle basis as in (2), but taking the coefficients c j to be time dependent: The functions ψ 0 , u j and v j are the solutions of equations (4)- (6) and are assumed to be static. Let us suppose that at t = 0 the condensate is in the ground state and thus all c j (0) are zero. For t > 0 the Bragg potential starts populating the quasiparticle states. The dynamics is entirely fixed by the evolution of the coefficients c j (t), that can be obtained by inserting the expansion (12) into the GP equation (1), with the external potential given in (11). One finds [3] In order to get an expression for theċ, one can multiply by u * i and v * i the last equation and its conjugate, respectively, and integrate in dr. By using the orthogonality and symmetry relations (3) one getṡ Finally, assuming that the quasiparticle occupation is always much smaller than N , one can replace with the ground state order parameter and integrate in time. The results is [3] which is expected to be valid in the limit of small perturbations.
In the case of a cylindrical condensate, we can use the factorization (8) to rewrite the previous expression in the form The last integral simply selects out the quasiparticles with k = ±q, so that where This means that once the ground state and the Bogoliubov spectrum are known at t = 0 the order parameter at time t > 0 can be calculated in a rather simple way. An important quantity that can easily be calculated within this scheme is the total momentum imparted to the condensate, which can be defined by integrating the current density associated with the order parameter:

54.7
By inserting the expansion (12) and using the orthogonality and symmetry relations (3), one gets the simple expression Now, let us take c n,k (t) from equation (17) and use the fact that the eigenfrequencies and eigenfunctions of the Bogoliubov equations (9)-(10) are even functions of k. One gets the final result In the next section we will show typical results obtained with this expression, but some relevant features are already evident. Notice first that, for positive ω and a given q, the momentum transferred (21) is clearly resonant at the frequencies ω = ω n,q . In other words, a peak occurs in P z whenever a vertical line at k = q crosses the Bogoliubov branches of figures 1(a) and (b). The separation between these branches is roughly 1-2 ω ρ , while the width of each peak goes like 2π/t. Thus, in order to resolve the different peaks, one has to wait at least a time of the order of trapping period T ρ = 2π/ω ρ . For t in this range, the second term in the square bracket gives certainly a negligible contribution. Finally, the contribution of each axial branch to the total momentum depends on the quantities W n,q , defined in equation (18). Typically, for a given q, the lowest branches have a greater weight in the summation (21), since the corresponding quasiparticle amplitudes have a greater overlap with the ground state. This behaviour of P z reflects its direct connection with the dynamic structure factor. At zero temperature, the latter is defined as whereρ q is the density fluctuation operator, |0 is the ground state, |i are the excited states and ω i = (E i − E 0 )/h. In terms of the quasiparticle amplitudes one has [20,21] If q is taken along z, and the system is translational invariant in the same direction, one can use equations (8) and (18) to rewrite S(q, ω) as Comparing this result with equation (21) one finds [3,15] As discussed in [3], these two quantities are connected in such a simple way only for a cylindrical condensate. If the condensate is also trapped along z, the relation between S(q, ω) and P z is less direct, involving a further integration on the duration time of the Bragg pulse. However, if the axial trapping frequency is significantly smaller than the radial one, there exists a wide range of time where the connection between the momentum transferred and the dynamic structure factor is roughly the same as for an infinite cylinder. In that range, the Bogoliubov branches, which correspond to δ peaks in S(q, ω), are also observable as distinguishable peaks in P z (t), as we will see in section 5.

Bogoliubov branches and multipeak spectrum
In figures 2(a) and (b) we show the momentum transferred to two cylindrical condensates with η = 9.4 and 26.5 by a Bragg pulse with wavevector q = 2.3 a −1 ρ (a) and q = 4 a −1 ρ (b). Recalling that ξ −1 = √ η, the values of q are chosen in order to have roughly the same value of qξ in the two cases. As solid curves we plot the quantity P z (t)/(NV 2 Bh q) as a function of ω and for three different durations: t = 2, 6 and 10 ω −1 ρ . According to equation (21), this quantity is independent of V B .
For a pulse shorter than the trapping period, the response is distributed over a single broad 54.9  (18), for the two condensates in figure 2. In the linear response regime, these functions provide the height of the δ-peaks in the dynamic structure factor (24), as well as of the peaks in the momentum transferred (21). The values of q used in figure 2 are shown as vertical dotted lines.
peak. Short pulses were used in the first experiments of [4,5]. In this regime, accurate predictions can be found by using the local density approximation (LDA) [15,21], i.e. assuming the system to be locally uniform. If the excitation spectrum is that of a uniform gas, but with a dispersion law fixed by the local density n(r), the dynamic structure factor S(q, ω) turns out to be non-zero only within the frequency range, ω r < ω < ω r [1 + 2µ/(hω r )] 1/2 , wherehω r =h 2 q 2 /(2m) is the free recoil energy [21]. This range corresponds to the shaded area in figures 1(a) and (b). The LDA approximation for S(q, ω) can be used to estimate the momentum transferred P z (t) as in [15]. The corresponding curves are shown as dashed curves in figures 2(a) and (b). The agreement between dashed and the solid curves for lowest curve at t = 2 ω −1 ρ is reasonable; for shorter times we reproduce the results of [15], and the agreement with the LDA becomes better and better.
For t of the order of T ρ = 2π ω −1 ρ new structures begin to appear, strongly deviating from LDA. In particular, one finds a multipeak spectrum that reflects the existence of Bogoliubov branches, i.e. the same as plotted in figure 1. The peak at the lowest frequency is the collective axial mode with no nodes in the radial direction (n = 0). The n = 1 and 2 branches are also visible for the condensate with η = 9.4 in (a), and even more for η = 26.5 in (b). The corresponding frequencies, ω n,q , are indicated by triangles on the horizontal axis of figure 2.
From equation (21) one sees that the height of each peak, for t T ρ , becomes proportional to t 2 |W n,q | 2 . In figure 3 we plot |W n,q | 2 as a function of q for the same condensates as figure 2. This figure tells us how many branches give significant contributions to P z and S(q, ω) at a given q. One notices that the contribution of each branch is sizeable when the corresponding ω n,q is, roughly speaking, within the shaded area in figure 1. The first branch is the dominant one at small q. Its weight |W 0,q | 2 also vanishes at q → 0, in agreement with the limiting behaviour of long-wavelength phonons, for which u 0,q −v 0,q ∝ q −1/2 , so that the integral (18) vanishes. One also notices that, for a given q, the number of branches that contribute to the summation (21) increases with η. This effect is rather dramatic for the condensates used by Ketterle and co-workers [4], where η 10 (see figure 4). In this case, even a small extra broadening in the experimental detection of P z might significantly reduce the visibility of the Bogoliubov branches.
The momentum transferred can also be calculated by direct integration of the time-dependent GP equation (1), using the definition (19) for P z . In the linear response regime the results must coincide with those of the quasiparticle projection method so far discussed. An example is shown in figure 5. The solid curve corresponds to the momentum calculated with equation (21) at t = 6 ω −1 ρ , as in figure 2(a). The empty squares are the results of GP simulations for a weak Bragg pulse (V B = 10 −3h ω ρ ). As expected, the two results are in good agreement, except for small differences that are compatible with the accuracy of our GP simulations. We also checked that, for such small values of V B , the quantity P z (t) scales with V 2 B as expected for the linear response.
A major advantage of using both approaches is that one can now increase V B beyond the linear regime and look at the differences between the linear response, given by equation (21), and the numerical GP calculation, which is valid also in the nonlinear regime [3,22]. A typical result is shown in figure 5, where we compare the small V B limit (solid curves and empty squares) with the results of GP simulations for V B = 0.5hω ρ (solid circles). With such a pulse intensity, the fraction of condensate atoms which are excited is roughly 25% at resonance. This value is slightly larger than that of typical experiments, where a fraction of the order of 10-20% is required for the optical detection of the excited atoms. The figure shows that, even for such highly excited condensates, the main features of the multipeak response associated with the underlying Bogoliubov spectrum are still visible. The most significant effects of nonlinearity, originating from the mean-field interaction in the GP equation, are that (i) P z does not scale with V 2 B , (ii) the peak frequencies are slightly shifted downwards, and (iii) additional structures appear superimposed on the original shape of the peaks.
This behaviour is a consequence of typical nonlinear processes like mode mixing and higher harmonic generation. Similar effects have already been predicted for spherical [3] and onedimensional [22] condensates.

Bragg scattering from a prolate elongated condensate
Here we show that the main features found for the response of a cylindrical condensate to a Bragg pulse also remain observable in a finite elongated condensate.
Let us take the trapping potential in the form where λ = ω z /ω ρ . If λ < 1 the condensate at equilibrium is a prolate ellipsoid. The ground state at t = 0 can be found as the stationary solution of equation (1). Then, the time-dependent GP equation can be solved at t > 0 to simulate the Bragg process. In figures 6(a) and (b) we show typical results obtained for P z in the linear response (small V B ) regime. We simulate the Bragg scattering from two different condensates, whose chemical potential, in units ofhω ρ , is equal to 9.4 and 26.5, as for the cylindrical condensates of figures 1 and 2. These two values correspond to the condensate used in the experiments of [1], made of N = 10 5 atoms of 87 Rb in a trap with ω ⊥ = 2π (220 Hz) and λ = 0.114, and to the one used in the calculations of [15,16], with λ = 0.125, respectively. With this choice, the condensates in figures 2(a), (b) and 6(a), (b) have also the same density profile in the radial direction for z = 0. As one can see, the momentum transferred to the trapped condensates of figure 6 has the same multipeak structure as the corresponding cylindrical condensates of figure 2.
The response of an ellipsoidal (axially confined) condensate is expected to differ from that of a cylindrical (axially unbound) condensate for two main reasons. First, their Bogoliubov spectra are different: the excitations of a confined condensate are also discretized along z and the quasiparticle amplitudes u and v have a non-trivial dependence on ρ and z. Second, the 54.12

54.13
response for long pulses is affected by the presence of the axial confinement, which can produce a reflection of quasiparticles at the boundaries and a centre-of-mass motion of the condensate in the trap.
The first effect, however, is small if q L −1 , where L is the axial length of the condensate. In this case, the Bragg pulse excites quasiparticles that behave like travelling waves along z, having a wavelength much shorter than the condensate size. One can safely introduce a continuous wavevector k to identify such excitations as in the multibranch spectrum of the cylindrical condensate. The axial size of the two condensates in figures 6(a) and (b) is approximately L = 75 and 120 a ρ respectively, so that the Bogoliubov modes can be safely represented by continuous branches down to about k = 0.5 a −1 ρ . Moreover, for k L −1 , the frequencies ω n,k for the elongated ellipsoidal condensate are expected to be close to the ones of the corresponding cylindrical condensate, with a similar multibranch structure. In order to check that the peaks of P z are still located at the frequencies of the Bogoliubov branches, we make a further simulation for the condensate in figure 6(a). We excite the condensate by acting with a Bragg pulse for a given duration t and then we let it to freely evolve in the trap. We perform a Fourier analysis of the induced density variations in the axial direction. As expected, we find that, for each q and ω, the density oscillates as a superposition of waves of different frequencies ω n,k . The latter are extracted from the Fourier spectrum and reported in figure 6(a) as solid circles. As one can see, the momentum is resonantly transferred at those frequencies.
The second effect is also small if the duration time of the Bragg pulse is significantly less than the axial trapping period T z . If λ is of the order of 10 −1 , or less, there exists a sufficiently wide interval of time between T ρ and T z where the multipeak structures of P z can be resolved without being significantly affected by the reflection of quasiparticles at the boundaries and by the centre-of-mass motion.
Finally, the experimental observation of the excited atoms requires rather high Bragg intensities. As already shown in the previous section, nonlinear processes are not expected to hinder the visibility of the multibranch Bogoliubov spectrum. Conversely, the combined use of both the quasiparticle projection method and the time-dependent GP equation might allow one to extract from the observed P z information about nonlinear effects in the condensate dynamics.
In figure 7 the results of our GP simulations are compared with experimental measurements of P z . Similar examples have already been reported in [1], where the major sources of noise and broadening of the observed data were also discussed. The main message of [1] was that the measured response to long Bragg pulses exhibits clear evidence of the underlying Bogoliubov spectrum. In the present paper we have provided a deeper theoretical analysis, which supports this view.

Stationary GP equation
The ground state configuration at t = 0 is obtained by solving the stationary GP equation (4) with the external potential given in (26). The case λ = 0 corresponds to the cylindrical condensate of sections 3 and 4: the order parameter is uniform along z and the problem is reduced to the one-dimensional equation (7). The solution of this equation is equivalent to the minimization of a GP energy functional [6]. We find the minimum by mapping the function φ 0 on a grid of points and propagating it in imaginary time with an explicit first-order algorithm, starting from a trial configuration, as described in [23]. The case λ = 0 corresponds to an ellipsoidal condensate, as in section 5. We use the same minimization algorithm. Due to the axial symmetry of V T , the problem is two-dimensional. The order parameter is mapped on a N ρ × N z grid (typically, 64 × 1024 points).

Time-dependent GP equation
Let us define an effective Hamiltonian H by rewriting equation (1) in the form ih∂ t = H . We solve this equation by propagating the order parameter in real time, starting from the ground state configuration. At each time step, t, the evolution is determined by the implicit algorithm where H and n are calculated at the nth iteration. For λ = 0, the propagation is split into axial and radial parts. The former is obtained by using a fast Fourier transform algorithm to treat the kinetic energy term, while the latter is performed with a Crank-Nicholson differencing method, as in [15]. For the cylindrical case (λ = 0), we first write the order parameter as (ρ, z, t) = ν ψ ν (ρ, t) exp(iνqz), where q is the wavevector of the Bragg potential, and then we propagate in time the radial functions φ ν , for ν = 0, ±1, ±2, . . .. Only a few values of ν give significant contributions to the sum, their number depending on the intensity of the pulse.

Bogoliubov equations
In order to solve the eigenvalue problem (9)-(10) for a cylindrical condensate, we project the unknown radial functions u(ρ) and v(ρ) on a basis set of Bessel functions. The problem is thus reduced to a matrix diagonalization. We use Bessel functions with up to 100 nodes in the computational box. The radius of the box is typically two or three times the size of the condensate.