1. Introduction

Interference of matter waves formed from Bose-Einstein condensates (BECs) [1, 2, 3, 4, 5] has been demonstrated experimentally [1, 2]. By virtue of the nonlinear nature of the self-interaction term in the Gross-Pitaevskii(GP) equation, which describes the dynamics of such systems at zero temperature, one may expect nonlinear phenomena to occur in BEC dynamics. The equivalent of the self-focusing nonlinearity in optical Kerr media [6, 7, 8] is actually self-defocusing for the case of positive scattering length. This example of nonlinear behavior has been observed in the expansion of the condensate due to the mean field energy when the trap is released [9]. Goldstein et al. [10, 11] have proposed that phase conjugation of matter waves should be possible in analogy to this phenomenon in nonlinear optics, including the case of multiple spin-component condensates [12]. They consider the case where a “probe” BEC wavepacket interacts with two counterpropagating “pump” wavepackets to generate a fourth that is phase conjugate to the “probe”; the “probe” is weak and causes negligible depletion of the “pump”. Law et al. [13] also suggest analogies between interactions in multiple spin-component condensates and four-wave mixing. Here we examine the four-wave mixing (4WM) in a single spin-component condensate that occurs as a result of the nonlinear self-interaction term in the GP equation when three BEC wavepackets with momenta k ₁, k ₂, and k ₃ collide and interact. Nonlinear 4WM can generate a new BEC wavepacket with a new momentum k ₁+k ₂-k ₃. Our assumptions on geometry and number of atoms in the wavepackets are less restrictive than those of Goldstein et al. We suggest that experiments with such wavepackets should be feasible, for example, using Raman output coupling techniques which have been demonstrated experimentally by the NIST group [14, 15]

2. Theory of four-wave mixing

The nature of 4WM in BEC collisions is unlike 4WM for optical pulse collisions in dispersive media [8, 16, 17], since the momentum and energy constraints imposed are different in the two cases. This is because the energy-momentum dispersion relation for massive particles is quadratic in k, whereas it is linear in k for the case of light. Moreover, in dispersive optical media, the momentum of light waves is proportional to the product of the frequency of the light and the refractive index, and the refractive index depends upon frequency (and the propagation direction in non-isotropic media -hence conservation of energy does not in general guarantee conservation of momentum in optical 4WM). This complication involving the properties of the medium does not arise in the BEC case. For 4WM in BEC collisions, when the magnitude of the momenta of each of the wavepackets are identical, i.e., |k _i| = |k _j| for i ≠ j, conservation of momentum and energy. However, in general, when |k _i| ≠ |k _j| for i ≠ j, conservation of momentum does not imply conservation of energy for 4WM in BEC collisions. Clearly, creation of new BEC wavepackets in 4WM is limited to cases when momentum and energy conservation are simultaneously satisfied.

Here we will study 4WM of BECs by numerical calculations on a model with three BEC wavepackets. There are two possible choices of initial conditions: (1) a “whole collision” in which initial spatially separated wavepackets come together at the same time, or (2) a “half collision” in which the wavepackets are initially formed in the same condensate at (nearly) the same time. Although we will assume initial condition (1), similar to optical 4WM experiments with light pulses, the “half collision” version (2) will be feasible exeprimentally using the methods of [14, 15]; similar conclusions regarding 4WM will apply to such a case. We assume the initial condensates are comprised of magnetically confined atoms in the same F, M_F state, and any trapping potentials are turned off before propagation begins. The three condensate wavepackets are given initial momenta and positions as in Fig. 1 so that they collide at a given point. We have carried out calculations in 1D, 2D, and 3D.

 

Figure 1. Schematic view of the initial positions and momenta of the three BECs wavepackets. The inset shows the momentum of the additional wavepacket formed by the 4WM process.

Download Full Size | PDF

The Gross-Pitaevskii(GP) equation for a single component BEC can be written as [18],

where $T_x=\frac{-{\mathit{ħ}}^2}{2m}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right)$ is the kinetic energy operator, V(x, t) is the potential imposed on the atoms and $U_0=\frac{4\pi a_0{\mathit{ħ}}^2}{m}{N}_T$ is the atom-atom interaction strength, proportional to the s-wave scattering length, a ₀, atomic mass, m, and the total number of atoms in all the wavepackets, N_T . The initial wavefunction is comprised of three BEC wavepackets,

where Ψ(x - x _i) is the solution to the GP equation with a locally harmonic potential centered around x = x _i, i = 1,2,3; the normalization constant N chosen so that the norm of Ψ is unity. We assume the three inital positions x _i are spatially separated so the initial wavepackets are non-overlapping (one could also consider the “half-collision” case where the three x _i are the same and the wavepackets are generated in situ from the same initial condensate). Although the initial wavepackets can have arbitrary phases multiplying the amplitudes Ψ(x - x _i) [3, 5], for simplicity we take the initial relative phase between wavepackets to be zero, since 4WM will occur for an arbitrary set of initial relative phases.

The nonlinear term in the GP equation will have terms with the factor exp[i(k _i + k _j - k _l) · x] where i, j and l can be 1, 2 or 3 respectively. These terms can generate a wavepacket with a central momentum that is not in the initial wavefunction Ψ(x, t = 0). For example, if k ₂ = -k ₁ (see Fig. (2)), then it is possible to produce a wavepacket with central momentum k ₄ = k ₁ - k ₂ + k ₃ = 2k ₁ + k ₃.

We set up a numerical calculation where the initial state, Eq. (2), evolves according to the GP equation (1). The time evolution is carried out using a split-operator Fourier transform method [19, 20]. We have verfied numerically that energy and momentum are conserved during the course of our calculation. Thus, dE(t)/dt = 0, where $E\left(t\right)=〈\Psi \left(t\right)\mid \left(T_x+\frac{1}{2}{U}_0{\mid \Psi \mid }^2\right)\mid \Psi \left(t\right)〉$ is the energy per particle in the BEC, and d P(t)/dt = 0, where P(t) = -iħ⟨Ψ(t)|∇|Φ(t)⟩ is the momentum per particle.

In order to estimate the importance of the various terms in the GP equation, it can be written in terms of characteristic time scales t_DF,t_NL in the following manner [19, 20]:

Here the diffraction time t_DF = 2${\mathit{\text{mw}}}_0^2$/ħ, and the nonlinear interaction time t_NL = (U ₀|Ψ_m|²/ħ)^-l, where ${\mid {\Psi }_m\mid }^2={\left(\frac{4}{3}\pi w_0^3\right)}^{-1}$ approximates the mean value of |Φ(x)|², and w ₀ stands for initial halfwidth of the colliding wavepackets. The smaller the characteristic time, the more important the corresponding term in the GP equation. We also define the collision duration time t_col = (2w ₀)/v, where v = k ₁/m is the initial group velocity of a wavepacket. The ratio t_col /t_NL gives an indication of the strength of the nonlinearity during the collision. The larger the ratio of t_col /t_NL , the stronger the effects of the nonlinearity during the overlap of the wavepackets. These characteristic times stand in the ratios $t_{\mathit{DF}}:t_{\mathit{col}}:t_{\mathit{NL}}=1:\frac{\lambda }{2\pi w_0}:\frac{w_0}{6a_0{N}_T}$, where λ is the De Broglie wavelength associated with the wavepacket velocity v. Experimental condensates with t_col /t_NL ≫ 1 can be readily achieved. Thus, the nonlinear term will have time to act while the BEC wavepackets remain physically overlapped during a collision.

3. Numerical simulations

We solve Eq. (3) numerically in reduced form by chosing the units of length x ₀ and time t ₀ so that (t ₀/t_DF )(w ₀/x ₀)² = 1/2; once x ₀ is chosen, t ₀ is given by t ₀ = ${\mathit{\text{mx}}}_0^2$/ħ. Here we choose the unit of length x ₀ to be x ₀ = 10μm, so that t ₀ = 36.2 ms for ²³Na atoms. Consequently, the unit of energy E ₀ = m(x ₀/t ₀)² = ħ/t ₀ is E ₀ = 2.91 × 10^-33 J = h(4.39 Hz), and the unit of momentum p ₀ = m(x ₀/t ₀) = ħ/x ₀ is p ₀ = 1-05 × 10^-29 Kg m/s. For comparison purposes, the recoil energy and momentum for a 589 nm photon (the Na resonance transition) are 5690E ₀ and 107p ₀ respectively.

In 1D, only head-on collisions of the condensates are possible. With the energy-momentum dispersion relation, E = ħ ² k ² /(2m), the constraints imposed by conservation of energy and momentum during the collision do not permit additional wavepackets to be created in 1D. Let us consider the following initial conditions: x ₁ = -x ₂, and x ₃ = 0, and k ₂ = -k ₁, and k ₃ = 0. Two wavepackets move symmetrically towards the central wavepacket whose center is at x ₃ = 0. Since the nonlinear term in the propagation equation is of third order in Ψ, and Ψ is a superposition of condensates with momenta k ₁, -k ₁ and 0, the nonlinear term could become a source of wavepackets propagating with momentum 0, ±k ₁, ±2k ₁ and ±3k ₁; in addition the collision could transfer population between the condensate wavepackets. We carried out numerical experiments with different values of nonlinearity U ₀, up to the value of the τ_col /τ_NL = 10. As expected, the only effect observed in 1D simulations was a slight delay of the maximum of the moving wavepacket peaks. No transfer of population, or additional peak creation was present, i.e., no wavepackets of momentum ±2k ₁ and ±3k ₁ were created. Even if we had taken k ₃ ≠ 0, or |k ₁| ≠ |k ₂|, no new wavepacket would appear in a 1D calculation.

 

Figure 2. Probability distribution |Φ(x,y,t)|² versus x and y in length units of x ₀ = 10 μm. Panels (a), (b), and (c) are for respective times t = -1, 0, and 1 t ₀ before, during, and after the collision, where t ₀ = 36.2 ms. The initial wavepacket started at t = -2t ₀, and expanded about 8-fold by the time of panel (a). Panels (d), (e), and (f) show the corresponding momentum distribution |Φ(k_x, k_y, t)|² versus k_x and k_y in $x_0^{-1}$ units.

Download Full Size | PDF

Next we consider the two dimensional case with initial configurations such that x ₁ = -x ₂ = (20,0), k ₁ = -k ₂ = (10,0), and x ₃ = (20α,-40α), k ₃ = (-10α, 20α), with several different values of the parameter α; α = 0.5, 0.7, 1, 1.25 and 1.5. We take N_T = 1.4 × 10^{6 23}Na atoms equally partitioned between wavepackets with initial w ₀ = 10 μm, for a typical tight trap with a mean trap frequency of 200 Hz. For this case t_DF = 72 ms, t_col = 7.2 ms, and t_NL = 0.031 ms. The wavepacket momenta k _i are significantly larger than the internal momentum. Fig. 2a shows the initial configuration of the colliding wavepackets. The initial conditions were selected such that the three wavepackets collide at t = 0 at the origin of the reference frame. Fig. 2 shows results of our 2D calculations for the case of α = 1.0. Fig. 2a–c shows the probability distribution |ΨΦ(x,y,t)|² before, during, and after the collision, and Fig. 2 d–f shows its Fourier transform, the momentum distribution |ΨΦ(k_x, k_y,t)|². Fig. 2b illustrates the wavepacket interference during the collision.

The striking feature in Figs. 2c and 2f is an additional 4WM wavepacket created in the collision with momentum k ₄ = k ₁ - k ₂ + k ₃. The three other wavepackets seen in Fig. 2c are the ones that pass through the collision region without changing their central momenta. Redistribution of population between different wavepackets must satisfy conservation of energy and momentum. Before or after the collision, when the wavepackets are separated in space, the population of the i-th one is given by the integral N_i = N_T ∫_Vi d x⟨Ψ_i|Ψ_i⟩ where the integration region V_i is selected to include the region around the i-th wavepacket. If we use N_i and N′_i to denote inital and final populations, conservation of energy and momentum show that N′ ₄ = ΔN ₃ = ΔN ₁ = -ΔN ₂, where ΔN_i = N′_i - N_i . These relationships are satisfied in our numerical simulations.

4. Interpretation

A physical interpretation of our results can be made in analogy with 4WM in nonlinear optics. We can regard the collision of the condensates as producing a grating formed due to the nonlinear term in Eq. (1). If Bragg conditions are satisfied (i.e., if energy and momentum are conserved in the formation of new peaks), new 4WM induced wavepackets may be formed in addition to the wavepackets present initially. Consider a nonlinear term in Eq. (1) and assume that Φ(x, t) consists of three wavepackets: Ψ(x, t) = ${\mathrm{\sum }}_{i=1}^3$ Ψ_i(x, t), as in Eq. (2). Hence, the nonlinear term |Ψ(x,t)|²Ψ(x,t) in Eq. (1) becomes a sum of nine contributions. Terms homogeneous in the index i describe self phase modulation (self-focusing or actually, self-defocusing for the case of positive scattering length), in analogy with nonlinear optics, and these terms can not be a source of new wavepackets. Moreover, crossed phase modulation terms of the form |Ψ_i(x,t)|²Ψ_j(x,t) for i,j = 1,2,3 and j ≠ i also can not contribute to the formation of new peaks. Only mixed terms, containing different indices may be a source of new wavepackets, but only when the Bragg conditions (i.e., conservation of momentum and energy) are satisfied. Specifically, these conditions are: (a) k ₄ = k _i - k _j + k _l, and (b) conservation of energy, which for our initial configuration (see Fig. 2) is equivilant to |k ₄| = |k ₃| where i = 3 and j = 1, l = 2. Referring again to the analogy with nonlinear optics, one can say that wavepackets with indices j and l create a grating during the collision whose grating vector is K = -k _j + k _l;, and the wavepacket i = 3 scatters off of this grating. This is confirmed by our numerical simulations. Only when α = 1 does the central K vector satisfy the Bragg conditions stated above. However, since the wavepackets contain momentum components spread around the central K vector, formation of the additional wavepacket may still occur for α ≠ 1 due to internal momentum compensating the momentum mismatch of the central K vectors. As |α - 1| increases, the number of atoms in the additional wavepacket decreases rapidly.

 

Figure 3. Blowup of the probability distribution |Ψ(x, y, t)|² of panel (c) in Fig. 2.

Download Full Size | PDF

It is important to note that the new 4WM wavepacket with momentum k ₄ is not a phase conjugate of the wavepacket with momentum k ₃, since the wavepackets with momenta k ₁ and k ₂ do not form a static grating, but rather, form a dynamic grating which changes in time due to the evolution and interaction of the wavepackets throught the course of the collision. The new k ₄ wavepacket is cigar-shaped, whereas the initial k ₃ wavepacket is spherical. The final k ₁ and k ₃ wavepackets are mirror images of one another (about a plane containing the centers of wavepackets 2 and 4). Their final shape is distorted relative to their initial shape; a bite has been removed from wavepackets 1 and 3 by the 4WM process which (1) created the new wavepacket with momentum k ₄ and (2) added Bose atoms to the wavepacket with momentum k ₂. Fig. 3 is an enlarged view of the final wavepackets in coordinate space with the details of their shapes shown more clearly. If we were to view the final wavepackets in the reference frame in which the new wavepacket is stationary, the trailing edges of wavepackets 1 and 3 would be reduced and the trailing edge of wavepacket 2 would be enhanced. The structure of the wavepackets clearly shows that the nature of the matter waves obtained after the collision is sensitive to the details of the collision dynamics and the properties of the initial wavepackets. A static grating picture is not sufficient to explain the results.

Our 3D calculations show similar results to the 2D ones. The quantity

which indicates the z-averaged distribution, is similar to the 2D distribution and to the probability distribution cut at the collision plane z = 0. These different distributions show only a few percent difference in the ratio of number of atoms in the four final wavepackets.

The formulation developed here assumes that the BECs are at zero temperature and can therefore be treated in the mean-field approximation as given by the GP equation. At finite temperatures, the collision must include the above-the-mean-field part of the wavefunction (order parameter), and a density matrix treatment is required. Moreover, even at zero temperature, there may be effects due to non-vanishing expectation values of the above-the-mean-field part of the wavefunction involving the creation of excitations due to the collision [21]. Moreover, in our treatment we have assumed that only one value of M_F is present for atoms in the condensates (e.g., F = 1, M_F = -1) with the z-axis perpendicular to the scattering plane. Another means of carrying out 4WM experiments is to use far off-resonance light traps to confine separate BECs [22], and impart momentum boosts to the BECs. In this case, the BECs would be multi-component in nature with all values of M_F being present in the condensates. A multi-component GP equation could be used to describe such experiments, as for example, the phase-conjugation experiment proposed by [12].