# Barrier resonances in Bose–Einstein condensation

## I Introduction

The recent experimental realization of Bose–Einstein condensation (BEC) in ultra-cold atomic gases,  has triggered the theoretical exploration of the properties of Bose gases. Specifically there has been a great interest in the development of applications which make use of the properties of this new state of matter. Perhaps, the recent development of the so–called atom laser  is the best example of the interest of these applications.

The current model used to describe a system with a fixed mean number $`N`$ of weakly interacting bosons, trapped in a parabolic potential $`V(\stackrel{}{r})`$ is the following Nonlinear Schrödinger Equation (NLSE) (which in this context is called the Gross–Pitaevskii equation (GPE))

$$i\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2}{2m}^2\psi +V(\stackrel{}{r},t)\psi +U_0|\psi |^2\psi ,$$
(1)   
which is valid when the particle density and temperature of the condensate are small enough. Here $`U_0=4\pi \mathrm{}^2a/m`$ characterizes the interaction and is defined in terms of the ground state scattering length $`a`$. The normalization for $`\psi `$ is $`N=|\psi |^2d^3\stackrel{}{r},`$ and the trapping potential is given by

$$V(\stackrel{}{r},t)=\frac{1}{2}m\nu ^2\left(\lambda _x^2(t)x^2+\lambda _y^2(t)y^2+\lambda _z^2(t)z^2\right),$$
(2)   
The $`\lambda _\eta ,(\eta =x,y,z)`$ are, as usual, functions that describe the anisotropies of the trap . In real experiments with stationary systems they are constants and the geometry of the trap imposes usually the condition $`\lambda _x=\lambda _y=1`$. $`\lambda _z=\nu _z/\nu `$ is the quotient between the frequency along the $`z`$-direction $`\nu _z`$ and the radial one $`\nu _r\nu `$.

At this point we want to emphasize here that our analysis does not restrict to the stationary case. Instead it focuses on situations where $`\lambda _\eta `$ are periodic sinusoidal functions of the time. To be more precise we will adopt the notation

$$\lambda _i(t)=\lambda _{i,0}(1+ϵ_icos(\omega _it))$$
(3)   
with $`i=x,y,z`$.

In its original derivation Eq. (1) was supposed to be strictly valid in the $`T=0`$ and low density limits. Recent theoretical work extends the applicability of the GPE to the high density limit . On the other hand linearized stability analysis based on perturbative expansions on $`1/N`$ seems to point that the validity of the equation is restricted to the cases where no exponential separation of nearby orbits appear as it happens for example in chaotic pulsations of the atom cloud .

A lot of work has concentrated in the analysis of the resonance structure when a periodic time dependent perturbation is applied to the magnetic field, mainly because of the availability of experimental data . The theoretical advances include both semi-classical analysis of Eq. (1) in the hydrodynamic limit , variational methods  and numerical simulations . Recent research based on Eq. (1) also includes the study of the different problems using tools from nonlinear science  and the analysis of its multicomponent extensions .

The variational approach provides us with very simple equations for the evolution of the widths of the Gaussian atom cloud together with a simple picture of the movement of the condensate . However, those equations were initially derived for the static field case and only applied to the analysis of the normal modes and frequencies of the condensate motion in the weak perturbation regime. Similar equations are found using various scaling arguments  or moment theory .

Evolution of the condensate in a time dependent trap has been addressed in many different ways. There are papers where the GPE is solved numerically, either by watching the time evolution of a suitable initial condition  or by using a linear expansion in normal modes . And there are other works  where the authors derive a set of ordinary differential equations using scaling arguments plus the Thomas–Fermi approximation. Among all the papers that treat time dependent potentials, the one which is closer to what we present here is . In that work the authors study the excitations of condensed and non–condensed clouds under a periodic parametric perturbation and find that the condensate depletes and that chaotic evolution comes out. We will comment more on this at Sec. VII.

It is our intention in this paper to explore the behavior of the condensate under periodic perturbations of the trap strength using analytical and numerical tools. Analytical techniques include exact results and a variational analysis extended to situations with a time dependent potentials. Comparison of the results of both methods will allow us to state rigorously and derive a simple model for the resonant behavior and to qualitatively predict the evolution of the system under other conditions (less symmetry, dissipation regimes, non–condensed corrections, etc). Other conclusions as well as experimental implications will be found in the course of the analysis.

Our plan is as follows: In Sec. II we obtain a set of exact equations for the evolution of the center of mass of the condensate. The result is found to be three uncoupled Mathieu equations. We show that a complete perturbative analysis is possible and obtain the resonance structure for the center of mass. In Section III we introduce the time dependent variational model for the time dependent GP equation. In Sec. IV we solve the GPE numerically in the radially symmetric case. We show the resonance structure for the condensate in all regimes – small and large amplitude oscillations –. Next we turn to the variational equations for the widths and show that numerical simulations agree qualitatively with the preceding picture. Finally, we use those ordinary differential equations to predict the locations of the resonances in a simple model that connects the linear and the nonlinear cases using tools from Sec. II. This is the main result of our paper. In Sec. V we consider the possibility of extending the analysis to the multidimensional case and study the effect of the coupling between the different variational parameters. The resulting predictions are found to agree with a posterior set of 3D simulations of the GPE equation. In Sec. VI we comment on the effect of losses. We show that these resonances are of a persistent nature and show how chaotic evolution may arise from losses effects. Finally in Sec. VII we discuss the experimental implications of our results and summarize our conclusions.

Note that, unless otherwise stated, all magnitudes are in adimensional units. To adimensionalize these figures we have employed the change of variables that is introduced in Sec. III

## II Exact analytical results

### A Variational form of the GP equation

It can be proved that every solution of Eq. (1) is a stationary point of an action corresponding, up to a divergence, to the Lagrangian density

$``$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}\left(\psi {\displaystyle \frac{\psi ^{}}{t}}\psi ^{}{\displaystyle \frac{\psi }{t}}\right)`$ (4)   
$`+`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}|\psi |^2+V(r)|\psi |^2+U_0|\psi |^4,`$ (5)   
where the asterisk denotes complex conjugation. That is, instead of working with the NLSE we can treat the action,

$$S=d^3r𝑑t=_{t_i}^{t_f}L(t)𝑑t,$$
(6)   
and study its invariance properties and extrema, which are in turn solutions of Eq. (1).

For instance, from the invariance of Eq. (4) under global phase transformations, one can assure the conservation of the norm of the wave function

$$N=|\psi |^2d^3r,$$
(7)   
which in this context is interpreted as the number of particles in the Bose condensed state. We can also define another quantity,

$$E=\left\{\frac{\mathrm{}^2}{2m}|\psi |^2+V(r,t)\psi +\frac{U_0|\psi |^4}{2}\right\}d^3r,$$
(8)   
that can be thought of as the energy, and whose time evolution is simply

$$\frac{dE}{dt}=\frac{dV}{dt}|\psi |^2d^3r.$$
(9)   
Thus, when the potential is not time dependent, the energy is another conserved quantity. And when the potential has the form (2), the evolution of energy can be easily connected to that of the mean square radii of the cloud

$$\frac{dE}{dt}=\frac{1}{2}m\nu ^2\underset{\eta =x,y,z}{}\frac{d\lambda _\eta }{dt}<\eta ^2>.$$
(10)   
Eqs. (7) and (10) are also useful to test the stability of the numerical scheme we will use to simulate Eq. (1).

### B Newton’s equations for the center of mass in a general GPE

Let us consider the following function

$$\psi (\stackrel{}{r})=\varphi (\stackrel{}{r}\stackrel{}{r}_0),$$
(11)   
where $`\varphi `$ is a solution of Eq. 1. Substituting it in Eq. (4) and calculating the averaged Lagrangian, we obtain

$$L=d^3r=L_{free}[\varphi ]+L_{cm}[\varphi ,\stackrel{}{r}_0].$$
(12)   
The Lagrangian has been splitted in two parts: one, the “free” contribution $`L_{free}`$, which depends only on $`\varphi `$

$`L_{free}[\varphi ]={\displaystyle \frac{i\mathrm{}}{2}\left(\varphi _t\varphi ^{}\varphi ^{}_t\varphi \right)d^3r}`$ (13)   
$`+{\displaystyle \frac{\mathrm{}^2}{2m}|\varphi |^2d^3r}+{\displaystyle \frac{U_0}{2}|\varphi |^4d^3r},`$ (14)   
and another one, $`L_{cm}`$, that includes both the potential and the displacement $`\stackrel{}{r}_0`$

$$L_{cm}[\varphi ,\stackrel{}{r}_0]=\left\{V(\stackrel{}{r}+\stackrel{}{r}_0)|\varphi (\stackrel{}{r})|^2i\mathrm{}\varphi ^{}\varphi \dot{\stackrel{}{r}}_0\right\}d^3r.$$
(15)   
If we impose that the action be stationary for some $`\stackrel{}{r}_0(t)`$, i.e. if we use Lagrange’s equations (37), we get

$$\frac{d}{dt}<i\mathrm{}>=<\frac{}{r_i}V(\stackrel{}{r}+\stackrel{}{r}_0)>.$$
(16)   
Here the brackets denote, as usual, the mean value of an operator over the unperturbed wave function, $`\varphi `$, $`<A>=\varphi ^{}(\stackrel{}{r})A\varphi (\stackrel{}{r})d^3r.`$

As a final step let us use that $`\varphi `$ is a solution of the GPE. Then Eq. (16) must be satisfied at least for $`\stackrel{}{r}_0=0`$

$$\frac{d}{dt}<i\mathrm{}>=<\frac{}{r_i}V(\stackrel{}{r})>.$$
(17)   
It is now easy to prove the relation between the mean value of the moment operator, $`i\mathrm{}`$, and the speed of the center of mass. We start from

$$\frac{d}{dt}<\stackrel{}{r}>=\stackrel{}{r}(\varphi ^{}\frac{}{t}\varphi +\varphi \frac{}{t}\varphi ^{})d^3r.$$
(18)   
Next we replace the time derivatives with spatial ones using Eq. (1) and its complex conjugate

$$\frac{d}{dt}<\stackrel{}{r}>=\frac{1}{i\mathrm{}}\stackrel{}{r}(\varphi ^{}\frac{\mathrm{}^2}{2m}^2\varphi \varphi \frac{\mathrm{}^2}{2m}^2\varphi ^{})d^3r,$$
(19)   
and finally integrate this expression to obtain

$$\frac{d}{dt}<\stackrel{}{r}>=<i\mathrm{}>.$$
(20)   
Eqs. (17) and (20) are the quantum equivalent of Newton’s Second Law and are exact for functions $`\varphi `$ that satisfy the GPE, in fact these equations coincide with the ones appearing in the linear Schrödinger equation.

### C The condensate in a harmonic trap. Mathieu equations for the center of mass

In our setup $`V(\stackrel{}{r})`$ is a harmonic potential with trap strengths of the form of Eq. (3). Thus, Eqs. (17) become a set of three decoupled ODE

$`{\displaystyle \frac{d^2}{dt^2}}<x>={\displaystyle \frac{1}{2}}m\omega ^2\lambda _x^2(t)<x>,`$ (22)   
$`{\displaystyle \frac{d^2}{dt^2}}<y>={\displaystyle \frac{1}{2}}m\omega ^2\lambda _x^2(t)<y>,`$ (23)   
$`{\displaystyle \frac{d^2}{dt^2}}<z>={\displaystyle \frac{1}{2}}m\omega ^2\lambda _x^2(t)<z>.`$ (24)   
After a change of scale, all of the preceding equations are equivalent to a model one that we will write as

$$\ddot{x}+(1+ϵ\mathrm{cos}(\omega t))x=0.$$
(25)   
This equation is known as Mathieu’s equation. It is a well know problem which appears frequently in the study of parametrically forced oscillators and where one can obtain a lot of information by analytical means .

First, Floquet’s theory for linear ODE with periodic time dependent coefficients  shows that Eq. (25) has an infinite set of instability regions in the parameter space. The limits of these zones can be found and have the shape of wedges that start on the points $`(\omega _{min},ϵ_{min})=(2,0),(1,0),(2/3,0),\mathrm{},`$ and widen as $`ϵ`$ is increased up from zero. Inside this regions at least one branch oscillates with an exponentially increasing amplitude.

Either with an asymptotic method, or by making use of the singular perturbation theory, we can also locate those resonances and study the evolution of the system around them. For a perturbation frequency close enough to the first resonance, that is for $`|\omega 2|=o(1)`$, an asymptotic method  yields up to first order

$`\sigma =\pm \sqrt{{\displaystyle \frac{ϵ^2}{4\omega ^2}}\delta ^2},`$ (27)   
$`rce^{\sigma t}\mathrm{cos}(\omega t/2+\theta _0).`$ (28)   
Here we see that for some values of $`\delta `$ and $`ϵ`$ the exponent $`\sigma `$ is a positive real number and the amplitude of the oscillations grows unlimitedly. Also, the strength of the resonance is maximum for a value of

$$\delta _{max}=1+\sqrt{1\frac{ϵ^2}{4}}ϵ^2+𝒪(ϵ^4).$$
(29)   
A second order Taylor expansion in Eq. (27) lays the following limits

$$|\omega 2|\frac{ϵ}{2}+\frac{ϵ^2}{32}.$$
(30)   
The treatment of other resonances is more difficult as they are caused by higher order terms –at least of second order in the $`\omega =1`$ case–. In practice this means that they have a smaller region of influence and that they are not so strong. One has to choose large values of $`ϵ`$, initial conditions $`(x,\dot{x})`$ not too close to the equilibrium point and an excellent numerical integration method, in order to find real instabilities. However, if one concentrates not just in looking for exponential divergences, but on efficient pumping, it will be easily observed that on top of these subharmonics there are peaks of the energy gain speed (See Fig 4b).

Finally, we wish to point out that these resonances are of a peculiarly persistent nature: as we will precise in section VI, they do resist even the presence of dissipation. This has a serious and inmediate consequence which is that feeding the condensate in a resonant regime can result in a large amplitude oscillation of the center of mass. On the other hand, as we will outline later in Sec. VI, a measure of this effect can give us some insight in the losses effects, as well as of any additional terms that could be added to Eq. (1) so as to better modeling the condensate.

## III Variational equations for the time–dependent GP equation with a periodic parametric perturbation

Although the center of mass of the wave packet satisfies very simple equations it does not happen the same to other parameters such as the width. Only in the two dimensional case it is possible to apply moment techniques and find analytically its time evolution  but in the fully three dimensional problem no exact results have been derived yet using the moment technique.

To simplify the problem, we restrict the shape of the function $`\psi `$ to a convenient family of trial functions and study the time evolution of the parameters that define that family. A natural choice, which corresponds to the exact solution in the linear limit ($`U_0=0`$) and provided quite good results in our previous works  is a three dimensional Gaussian-like function with sixteen free parameters

$$\psi (x,y,z,t)=A\underset{\eta =x,y,z}{}\mathrm{exp}\left\{\frac{[\eta \eta _0]^2}{2w_\eta ^2}+i\eta \alpha _\eta +i\eta ^2\beta _\eta \right\}.$$
(31)   
For a matter of convenience and ease of interpretation we will make here a change of parameters, from $`A`$ and $`A^{}`$ to $`N`$ (the norm of the wave function) and $`\varphi `$ (its global phase):

$$A=\frac{N}{\pi ^{3/2}w_xw_yw_z}e^{i\varphi }.$$
(32)   
The rest of the parameters are $`w_\eta `$ (width), $`\alpha _\eta `$ (slope), $`\beta _\eta `$ (square root of the curvature radius) and $`\eta _0`$ (center of the cloud).

This trial function must now be placed in Eq. (6) to obtain an averaged Lagrangian per particle

$`{\displaystyle \frac{L}{N}}={\displaystyle \frac{1}{N}}{\displaystyle _{\mathrm{}}^+\mathrm{}}d^3r=`$ (33)   
$`\mathrm{}\dot{\varphi }+{\displaystyle \underset{\eta }{}}\left\{{\displaystyle \frac{w_\eta ^2}{2}}+\eta _0^2\right\}\left\{\mathrm{}\dot{\beta }_\eta +{\displaystyle \frac{2\mathrm{}^2}{m}}\beta _\eta ^2+{\displaystyle \frac{1}{2}}m\nu ^2\lambda _\eta ^2(t)\right\}`$ (34)   
$`+{\displaystyle \underset{\eta }{}}\left\{\mathrm{}\dot{\alpha }_\eta +{\displaystyle \frac{\mathrm{}^2}{m}}2\alpha _\eta \beta _\eta \right\}`$ (35)   
$`+{\displaystyle \frac{\mathrm{}^2}{m}}{\displaystyle \underset{\eta }{}}\left\{{\displaystyle \frac{1}{2w_\eta ^2}}+\alpha _\eta ^2\right\}+{\displaystyle \frac{U_0}{4\sqrt{2}}}{\displaystyle \frac{N}{\pi ^{3/2}w_xw_yw_z}}.`$ (36)   
The evolution of the parameters is ruled by the corresponding set of Lagrange equations

$$\frac{d}{dt}\left(\frac{L}{\dot{q_j}}\right)=\frac{L}{q_j},$$
(37)   
which give us equations for the conservation of the norm,

$$\frac{dN}{dt}=0;$$
(38)   
the movement of the center of mass,

$$\ddot{\eta }_0+m\nu ^2\lambda _\eta (t)\eta _0=0;$$
(39)   
the evolution of slope and curvature,

$`\beta _\eta `$ $`=`$ $`{\displaystyle \frac{m\dot{w}_\eta }{2\mathrm{}w_\eta }},`$ (41)   
$`\alpha _\eta `$ $`=`$ $`{\displaystyle \frac{m}{\mathrm{}}}\dot{\eta }_02\beta _\eta \eta _0;`$ (42)   
and finally the evolution of the widths,

$`\ddot{w_x}+\nu ^2\lambda _x^2(t)w_x`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{m^2}}{\displaystyle \frac{1}{w_x^3}}+{\displaystyle \frac{U_0}{2\sqrt{2}m}}{\displaystyle \frac{N}{\pi ^{3/2}w_x^2w_yw_z}},`$ (44)   
$`\ddot{w_y}+\nu ^2\lambda _y^2(t)w_y`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{m^2}}{\displaystyle \frac{1}{w_y^3}}+{\displaystyle \frac{U_0}{2\sqrt{2}m}}{\displaystyle \frac{N}{\pi ^{3/2}w_xw_y^2w_z}},`$ (45)   
$`\ddot{w_z}+\nu ^2\lambda _z^2(t)w_z`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{m^2}}{\displaystyle \frac{1}{w_z^3}}+{\displaystyle \frac{U_0}{2\sqrt{2}m}}{\displaystyle \frac{N}{\pi ^{3/2}w_xw_yw_z^2}}.`$ (46)   
As one can see, the introduction of a time dependent potential does not affect the form of the equations, which remain the same as those of .

Let us introduce the constants $`P=\sqrt{2/\pi }Na/a_0`$ (strength of the atom-atom interaction) and $`a_0=\sqrt{\mathrm{}/(m\nu )}`$ (harmonic potential length scale), as well as a set of rescaled variables for time, $`\tau =\nu t`$, and the widths, $`w_\eta =a_0v_\eta ,(\eta =x,y,z)`$. This leads us to

$`\ddot{v}_x+\lambda _x^2(t)v_x`$ $`=`$ $`{\displaystyle \frac{1}{v_x^3}}+{\displaystyle \frac{P}{v_x^2v_yv_z}},`$ (48)   
$`\ddot{v}_y+\lambda _y^2(t)v_y`$ $`=`$ $`{\displaystyle \frac{1}{v_y^3}}+{\displaystyle \frac{P}{v_xv_y^2v_z}},`$ (49)   
$`\ddot{v}_z+\lambda _z^2(t)v_z`$ $`=`$ $`{\displaystyle \frac{1}{v_z^3}}+{\displaystyle \frac{P}{v_xv_yv_z^2}}.`$ (50)   
This set of rescaled variables and units is the one that we will use throughout the paper, unless otherwise stated.

## IV Analysis of the radially symmetric case

### A The equations

In this section we will analyse the case where the trap has the same strength on all directions, that is

$$V(x,y,z)=\frac{1}{2}m\nu ^2\lambda ^2(t)(x^2+y^2+z^2).$$
(51)   
and the solutions are supposed to be radially symmetric. This high degree of symmetry simplifies the equations considerably. First, the variational model for the widths of the condensate (44) reduces to a single ODE for the radial width $`v(t)`$

$$\ddot{v}=\lambda ^2(t)v+\frac{1}{v^3}+\frac{P}{v^4}.$$
(52)   
And secondly the following change of function

$$\psi (r,t)=A\frac{u(r)}{r},$$
(53)   
with the constrains

$`u(r)0,r0`$ (55)   
$`{\displaystyle _0^{\mathrm{}}}u(r)𝑑r=1`$ (56)   
$`|A|^2=4\pi N,`$ (57)   
transforms Eq. (1) into the one-dimensional PDE

$$i\mathrm{}_tu=\frac{\mathrm{}^2}{2m}_r^2u+\left\{\frac{1}{2}m\nu ^2\lambda ^2(t)r^2+4\pi \frac{U_0N}{4\pi }\frac{|u|^2}{r^2}\right\}u.$$
(58)   
This is the equation that we have actually solved numerically.

### B Numerical study of the equations

The numerical solution of Eqs. (52) and (58) is not a trivial job. We can see without much effort that both equations are stiff  when the width of the cloud becomes very small, due to the presence of strong singular potentials. Thus, we need numerical methods that account for the nonlinearities and are stable enough to be trusted when close to a resonance.

To solve Eq. (52) we have used an adaptive step size Runge–Kutta–Fehlberg method, a Dormand–Prince pair , the ODE Suite of MATLAB , and finally Vazquez’s conservative scheme  –a finite differences scheme that conserves a discretized version of the energy and is unconditionally stable–. All of them gave the same accurate results about the frequencies and amplitude of the oscillations, the regions of divergence, etc.

To solve Eq. (58) we have utilized a modification of a second order accurate finite difference scheme developed in . This new scheme is time reversible, conserves the norm and has a discrete analogue for Eq. (9) which provides enhanced stability . On the other hand, even using the best methods, we face another important difficulty which is the finite size of either the spatial grid –in finite differences schemes– or the momentum space –in spectral or pseudospectral methods–. This size effect becomes particularly important in the case of parametrical perturbations and imposes a severe limit on the time for which simulations may be trusted.

Using all this computational machinery we have achieved several important results. First we have checked our programs with low amplitude oscillations. In the PDE we imposed a gaussian initial condition of width $`v0`$ and used this same value as initial condition for Eq. (52). By this procedure we obtained the linearized excitation frequencies of the condensate, concluding that there is a significant coincidence of both models as was shown in . In the large amplitude simulations we found a surprising result which is that the variational model still follows the evolution of the PDE with a 90% of accuracy even in situations where the condensate remains no longer gaussian but gains an important contribution from the first and second modes. Another confirmed fact is that the width of the condensate develops fast tough bounces against the “origin” (Fig. 1a-b), with its excitation frequency deviating from the linearized predictions of  and approaching those of the harmonic trap.

We have also simulated the system with a periodic time dependent perturbation of the form

$$\lambda ^2(t)=1+ϵ\mathrm{cos}(\omega t).$$
(59)   
Using the variational equation Eq. 52 we have scanned the parameter space $`(\omega ,ϵ)`$. We have found one region where the radial width diverges exponentially (Fig. 2a) and causes most numerical methods to fail after a finite time, and two more where the width grows almost exponentially up to a point where the simulation cannot account for this growth. All of these zones have the form of wedges, with a peak at $`(\omega _{min},ϵ_{min})`$, and a growing width as $`ϵ`$ is increased. The most important resonance stands on $`\omega _{min}=2.04`$ (See Fig. 2b). The other ones are weaker and rest on $`\omega _{min}=1.02,0.68`$. We have studied a wide range of setups and found that these frequencies change no more than a $`0.5\%`$ depending on the initial conditions and the nonlinearity. On the other hand, the lowest perturbation amplitude for which the resonance exists $`ϵ_{min}`$ does exhibit a strong dependence on the initial conditions. Indeed, around the weaker resonances, instability is never reached for points close to the minimum of the potential. However, as we already pointed in Sec. II this is probably a numerical effect.

In order to study the location of the resonances we have also made some plots of the efficiency of the energy absorption process against the perturbation frequency for different values of the perturbation amplitude for the variational system and the partial differential equation (Figs. 4a-b, 3a-b). The way we have measured “efficiency” is by letting the system evolve for a fixed time and then computing the maximum value of the mean square radius of the cloud. For the sake of simplicity and to approach the experimental setups, we have used an equilibrium state of the static GPE as initial condition.

In a rather complete inspection of the parameter space using the efficiency plots, we have found that, though the resonance regions cannot be precisely delimited because of the dependence on the initial data, they do exist and behave much like the variational model predicted. As we see in Figs. 4 and 3 there are two important features in the response of the condensate. The main one is the width of the resonances and the dependence of that width on the strength of the perturbation. The second important feature is the change of the frequency for which the perturbation is most efficient. This peak is centered on the frequency of the linearized model only for very small perturbation amplitudes, and switches to the trap natural frequencies very quickly as the amplitude is made stronger – of about $`10\%`$ or so–. For even stronger amplitudes the optimal frequency decreases slowly.

Another consequence of this work is that response of the cloud is stronger in the PDE than in the variational simplification. For instance, it is possible to find (See Fig. 2a) perturbation amplitudes that do not cause an significant growth in Eq. (52) but make the cloud width increase quite linearly in the PDE. We have also studied the case where an exponential growth of the width is present in both equations –(52) and (58)– and we have seen that the growth is qualitatively similar though there is a tendency in the exact model to exhibit slightly larger amplitude motion. These discrepancies are originated by high modes which are not present in the variational treatment. In fact, even though an important part of the cloud remains close to the origin, it is observed that for long times a long tail appears which is hard to appreciate but is responsible for the growth of the mean square width (and the presence of higher order modes).

### C Analysis of the Mathieu equation with a singular potential

In this subsection we are to develop a simplified model that explains why do resonances appear in the perturbed GPE equation. This model makes use of the variational equations for the cloud width but does not care for the actual shape of the variational ansatz. The reason for it is that both the approximate model and the exact one agree for short times, the response of the latter being always stronger.

Let us limit ourselves to Eq. (52). In the previous subsection we said that this system is very stiff and that the origin acts as an elastic wall. In view of this, it is intuitively appealing to replace the singular but differentiable potential $`1/v^3+P/v^4`$ with a discontinuous bounce condition on the origin, i.e. an impact oscillator. Numerical simulations confirm that this approximation is good for large amplitude oscillations. Also, the connection between soft singular potentials and impact oscillators has been widely studied and these kind of models have found ample application in real–life situations in the fields of Physics and Engineering (See  for many references).

Since the variational model is a lossless one, our boundary condition must be elastic. We replace Eq. (52) with the following one

$`\ddot{v}+\lambda ^2(t)v`$ $`=`$ $`0,`$ (60)   
$`\underset{tt_c^{}}{lim}(v,\dot{v})=(0^+,V_c)`$ $``$ $`\underset{tt_c^+}{lim}(v,\dot{v})=(0^+,V_c),`$ (61)   
where $`t_c`$ denotes any isolated instance when the system bounces against the $`v=0`$ singularity.

Let us show that this equation is in turn equivalent to an elastic oscillator without barrier conditions. We introduce the change of variables

$$v=|u|,$$
(62)   
where $`u`$ is an unrestricted real number and satisfies the following one–dimensional harmonic oscillator equation

$$\ddot{u}+\lambda ^2(t)u=0,$$
(63)   
It is now easy to prove that every solution of Eq. (63) provides a solution of Eq. (60). And vice versa, from every solution of Eq. (60) it is possible to construct a solution of Eq. (63), unique up to a sign.

So, what do we have now? We have proved that for large amplitude motion Eq. (52) behaves as Mathieu’s equation. This implies that in the regime of medium to large amplitude oscillations, or in a situation of large amplitude perturbations, Eq. (52) will have instability regions that are more or less centered on Mathieu’s frequencies $`\omega =2,1,1/3,\mathrm{}`$. This prediction is indeed confirmed by the numerical simulations: the main resonance is capable of causing an exponential divergence, while the other ones are harder to track down but do appear in plots of pumping “efficiency”.

Another, important but minor result of this equivalence is that these resonance regions must become wider and move on to smaller frequencies as the perturbation amplitude is increased. This result is also obtained in the numerical simulations. Indeed, from the graphs we can estimate how fast the optimal frequency decreases as $`ϵ`$ grows, and we will see that the order of magnitude corresponds to that of Eq. (29).

Summing up, what we have seen here is that for low amplitude oscillations the condensate moves in an effective potential that is parabolic. Thus, the frequency that excites the condensate most efficiently in a parametric way is the one that results from the linearization of Eq. (52). On the other hand, when we start to consider large amplitude oscillations, we find that the harmonic trap gains importance over the details of the well. It is in this situation that the instability arises, and it happens for perturbations that oscillate according to multiples of the frequency of the trap.

Even more, higher modes are little influenced by the nonlinearity, just because they are more spread and the value of $`|\psi (r,t)|^2`$ is smaller. As a consequence, the energies of these modes come even closer to those of the linear harmonic oscillator, resulting in the fact that the response of the GPE is stronger than what the variational model, limited to a gaussian shape, predicts.

## V Analysis of the nonsymmetric case

After finishing the study of the radially symmetric problem, it seems a natural step to proceed with the nonsymmetric one. However this step is not simple for several reasons, the main one being that the simulation of full three dimensional GPE is computationally very expensive work. Thus, it would be unwise to directly attack the full problem without gaining some insight on what is to be expected from the simulations by cheaper means.

In our case the cheapest tool is the variational model. We have seen that it describes rather well the behaviour of a radially symmetric condensate. And, as we pointed out before, there are exact analytical studies  where the moment method reduce exactly the two dimensional GPE to a set of ODE which are similar to the ones we have.

### A Predictions of the variational model

When we remove the radial symmetry in the variational equations, we are left with two to three coordinates, and the perturbation can bear many different forms. However, trying to follow the experimental setups , we should once more take a sinusoidal time dependence for every $`\lambda _\eta (t)`$ coefficient, as in Eq. (3). This choice accounts both for the $`m=0`$ mode $`(ϵ_x=ϵ_y,ϵ_z=0)`$ and the $`m=2`$ perturbations $`(ϵ_x=ϵ_y,ϵ_z=0)`$ from the JILA experiment . In the latter case the potential is a parabolic one, with fixed frequencies on a rotating frame. However, for our trial function (Eq. (1) it behaves just as a potential of the form of Eq. (2)).

Substituting our effective perturbation frequencies into Eq. (48) we get a set of three coupled Mathieu equations with a potential that is singular on the $`v_x=0`$, $`v_y=0`$ and $`v_z=0`$ planes. The singularities are at least as strong as $`1/v^3`$, and the numerical simulations again confirm that they act as elastic walls, so we now proceed with a change of variables formally equivalent to that of Eq. (62):

$`w_\eta `$ $`=`$ $`|u_\eta |`$ (65)   
$`\ddot{u}_\eta `$ $`=`$ $`\lambda _\eta ^2(t)u_\eta `$ (66)   
for $`\eta =x,y,z`$.

Now the situation is a bit more complex. The first new feature is the existence of several sets of instability regions. Due to having three a priori different constants $`\lambda _{0\eta }`$, the three oscillators in Eq. (V A) are not equivalent and we may get three sets of resonances in the $`(ϵ_\eta ,\omega )`$ space, each one containing the instability regions that start on the frequencies

$$\frac{\omega _{\eta ,min}}{\lambda _{0\eta }}=2,1,2/3,\mathrm{}$$
(67)   
Numerical simulations of the variational equations for the $`m=0`$ type excitation confirm this prediction with a relative accuracy around $`0.5\%`$ in the frequencies. The results show again that, opposite to the pure Mathieu equation, these “wedges” rest on a nonzero value of the perturbation amplitude $`ϵ_{min,\eta }`$. In Fig. 5 we see the evolution of the condensate width with parameters close to the main resonance region.

Another new feature is the possibility of coupling between the widths of the condensate. This coupling is seen both in the $`m=0`$ and $`m=2`$ excitations setups. In the first one, the perturbed width feeds the unperturbed one. Figure 3 demonstrates that efficiency is not very high. In the $`m=2`$ case two widths are associated to the same trap frequencies and the perturbations are of equal magnitude and opposite sense. The explanation is that both widths block each other, eliminating the resonance and leaving just a bounded “movement” of small amplitude.

As a side note, we must say that the variational model itself predicts the lack of resonances in the $`m=2`$ setup. We already pointed out that the JILA $`m=2`$ perturbation corresponds to a situation where the potential of the trap is rotated while maintaining its shape. It can be easily proved that if we choose any family of trial functions with enough degrees of freedom for a general rotation, the variational solution will always stick to the potential and rotate at the same speed. This also implies that the trial function (31) is not suitable for describing the condensate when this kind of perturbation is applied.

### B Numerical study of the three dimensional GPE

To perform the simulation of the full GPE we have used a fourier pseudospectral method, using typically a grid of $`108^3`$ collocation points and integrating in time with a second order, symmetrized split step method .

Our numerical study is based on a $`O((\mathrm{\Delta }t)^2)`$ scheme which behaves extremely well for long time runs. However, as in the radial PDE, no matter how accurate the scheme is, there’s a limit in the time during which simulations can be trusted and this limit is imposed by the growth rate of the condensate width and the size of the grid. This limit is specially important for the $`m=0`$ perturbation, where the condensate develops a large tail in the unperturbed direction, and this tail breaks the simulation after a certain time. The more asymmetric the trap is, the sooner this effect comes out. In our simulations the grid was a box of $`108^3`$ equally spaced points and whose sides measured from 20 to 40 length units, depending on the symmetry and intensity of the trap. This allowed us to track the condensate for about 12 periods in truly resonant setups, and for much more in nonresonant ones.

We have applied the algorithm to many different problems. First we checked our programs against stationary and radially symmetric problems. Secondly we introduced time dependent traps and reproduced the calculations of Section IV. In both cases we got the expected results.

The third set of experiments consisted in a resonant time dependent radially symmetric trap applied to several slightly asymmetric gaussian wave packets. The initial asymmetry was slight enought to treat it as a weak perturbation and we saw that it departured little from the symmetric case –i.e., no modes with higher energy or angular moments break the exponential growth.

The fourth and probably most important group of simulations consisted in the study of the $`m=0`$ perturbation from the JILA  experiment. Here we confirmed the predictions of the previous subsection, that is we obtained at least one resonance region where the variational model predicts.

We also observed that the response of the condensate as modeled by the GPE is stronger and exhibits a slightly more intense growth rate than what the variational model predicts. This and other similar results from the radial case (Sect. IV) favor the theory of a cooperative staircase effect, where the higher modes contribute to the energy absorption process without interfering in the evolution of the width. Figs. 8 and 7 show just two examples of the kind of evolution of have seen for this perturbation model on the resonant regions.

The reason for this cooperative effect seems to lay in the energy level structure of the GPE equation. We have studied the evolution of the correlation a condensate wave against its initial data. In all cases the initial data was a displaced and deformed gaussian cloud, while the environment corresponded to a stationary trap. In the linear case it is easy to show that the spectrum of the correlation must reveal a subset of the eigenvalues of the hamiltonian, and indeed that is what we got. In a nonlinear context it is not clear what the frequencies of the correlation mean –they may be or may be not eigenvalues of the GPE–, but at least we know that they must rule the energy absorption process somehow. What numerical experiments show is that for an extensive family of initial conditions these generalized “spectra” can be approximated by the formula $`E_n=\omega n+E_0`$, where $`E_0`$ depends on the nonlinearity and the level spacing is a regular, harmonic one.

This picture of equally spaced levels is intuitively appealing to explain the existence of such strong resonances. On one side we have the fact that, in any other system, a continous parametric perturbation would become bounded as the higher state become populated. This is probably what one would first think when facing this system. On the other side we see that due to being equally spaced these higher energy states are themselves sensitive to the same resonant frequency and offer no resistence to the particle promotion process. In the end, this is more or less what happens in the variational equations, considered from the point of view of Classical Mechanics.

However, it is not relevant for the existence of resonances whether they cause a sustained growth or not. What is confirmed without any doubt is that the perturbations that are most efficient in the variational model are also the most efficient in the full GPE.

## VI Analysis of the effect of losses

We now want to show the effect of a dissipative term in the variational equations. This term will be introduced in a phenomenological way so as to model the damping of the oscillations of the condensate in regimes where the number of noncondensed particles is small. We will choose a viscous damping term that models well the behaviour of the condensate in the experiments . This choice introduces a significant loss of energy in the oscillations while preserving the number of condensed particles. Using this term we will see that the resonance regions of the pure Mathieu persist.

Let us see what Eq. (25) looks like once we add damping:

$$\ddot{x}+(1+ϵ\mathrm{cos}(\omega t))x+\gamma \dot{x}=0.$$
(68)   
A simple change of variables $`x(t)=\rho (t)e^{\gamma t}`$ makes this new term disappear, transforming it back to a pure Mathieu equation

$$\ddot{\rho }+(1\gamma ^2+ϵ\mathrm{cos}(\omega t))\rho =0.$$
(69)   
With the introduction of the damping we are shifting the resonances to lower values, given by

$$\frac{\omega }{\nu (\gamma )}=2,1,2/3\mathrm{}$$
(70)   
where $`\nu ^2(\gamma )=1\gamma ^2`$ is the new effective frequency for the trap. We can approximately solve Eq. (69) around the first resonance, obtaining

$$x(t)ce^{(\sigma \gamma )t}\mathrm{cos}(\frac{\omega t}{2}+\theta _0),$$
(71)   
where $`\lambda `$ is given by Eq. (27).

This shows that the resonance regions in the parameter space are constrained to values of $`(\omega ,ϵ)`$ for which the strength of the resonance, $`\lambda `$, is larger that the strength of the dissipative term. The new regions have a larger, nonzero, value of $`ϵ_{min}`$, and are typically thinner, but do not disappear unless $`\gamma `$ is very large.

This effect is reproduced in the variational model when we introduce similar viscous damping terms. For instance, taking the data from the JILA experiment , we can estimate a condensate lifetime of about 110 ms and a value for $`\gamma `$ of $`0.15`$ in natural units of the condensate. Such damping makes the $`ϵ_{min}`$ value raise from $`0.09`$ to $`0.18`$ for the $`P=9.2`$ case. Thus, the instability should not be appreciated unless the perturbation amplitude exceeds the 20$`\%`$.

An interesting effect of damping is that the evolution of a continuously perturbed condensate outside the instability regions becomes more ordered than in the undamped mode since the motion is constrained to limit cycle synchronized to the frequency of the parametric perturbation, and with a size that depends only on the perturbation parameters, $`(\omega ,ϵ)`$. In Fig. 6 we show the different limit cycles that appear under peridic perturbations. The largest one is always the one with its frequency on top of the peak of the resonance as shown by Fig. 3 and the size of the limit cycle decreases as the frequency is detuned from this value.

Finally we wish to point out that the appearance of a limit cycle opens the door to a wide family of phenomena, from chaotic motion to bifurcation theory . This limit cycle would exist under a great variety of dissipative terms, and is not exclusive of linear damping. Also, the dependence of the limit cycle on the damping constant can be useful from the experimental viewpoint to separate the condensed and noncondensed clouds as will be pointed out in the last section.

## VII Conclusion and discussion.

In this work we have analysed the resonant dynamics of the parametrically forced time dependent GPE using exact analytical techniques, approximate time dependent variational techniques and numerical simulations. All the results point to the existence of various resonant behaviors associated to the same parameter regions concerning the motion of the center of mass and the width oscillations of the wave packet. The role played by the nonlinearity is to provide a strong repulsive term at the origin which acts as a barrier, which depending on the dimensionality and the symmetry of the external forcing could be stronger than the repulsive term related to the linear dispersion (kinetic energy term).

We have developed a simplified version of the variational equations for the spherically symmetric condensate under periodic sinusoidal change of the trapping potential. This model is based on an impact oscillator, i.e. a harmonic oscillator with an elastic barrier condition which allows to find explicit solutions. When applied to the parametrically perturbed condensate, our model predicts that medium to large oscillations approach the harmonic trap frequencies, not the ones resulting from the linearization of the variational equations. The model also predicts that the response of the condensate to an external perturbation is ruled by the harmonic trap frequencies. In the case of a sinusoidal parametric perturbation it accurately predicts the existence of a family of resonances on multiples of the trap natural frequencies.

In the end, what we get from this work is a precise picture of the resonances for a wide family of equations that include the Eq. (1) and the harmonic oscillator. In this description for the radially symmetric case we have one base frequency that is essentially the same for both equations (linear and nonlinear one) and which corresponds to the energy separation between the ground state and the first excited state of the harmonic oscillator up to a high precision. The invariance of this base frequency has been checked for a nonlinearity constant going from $`P=9.2`$ –the JILA  experiment– to 20 times this value. It differs from the predictions of , where formulas regarding the $`P0`$ and $`P\mathrm{}`$ are derived. Also there is a whole set of subharmonics of this frequency, all of which are capable of exciting the cloud quite efficiently. At least three of them have been found, both with significant responses. These subharmonics are not predicted in .

We have also found that for this kind of parametric drive the resonances are wide. The width grows with the strength of the interaction and decreases with the effect of dissipation. Both facts can be checked in the experiments by forcing the system for a longer time than what it is currently done.

Both predictions are exact in the linear limit, $`U_0=0`$, and have been confirmed with simulations coming both from the exact variational model and the GPE. The discretization of the GPE has also allowed us to scan the parameters space, studying the efficiency of the perturbation process. In this study we have only found peaks centered on Mathieu’s frequencies.

For the general case with more than one degree of freedom (axial and nonsymmetric cases) we obtain a set of two to three decoupled pure Mathieu equations. We have shown that, due to having more than one frequency, the predicted Mathieu resonances do exist in a larger number. On the other hand, we have also seen that some of these resonances may disappear due to the locking of “equivalent” variables, an effect that our decoupled equations do not account for.

This resonance scheme for the nonsymmetric case has been confirmed with accurate simulations of the full GPE for the $`m=0`$ perturbation. An analysis of the correlation of a state against its initial data shows that both the linear and the nonlinear problems exhibit a spectral structure which is likely to present such behavior.

Finally, damping has been shown to limit the effect of the parametric perturbation. Once more we have proved that only frequencies close to the Mathieu resonance regions do excite the condensate as a whole in an efficient way, causing the appearance of a stable limit cycle. All other frequencies are inefficient in the sense that the system stays extremely close to the equilibrium configuration, which acts as a focus.

A main conclusion of this work is that for this set of resonances to exist one only needs a singularity that prevents from collapse. The variational method showed that the kinetic terms in the evolution equations guarantee a $`1/x^3`$ singularity as far as we impose a repulsive interaction between the atoms in the cloud. This is the reason why we say that we have a family of systems that behave much the same. An immediate result of this is that the response of the noncondensed atoms under the parametrical perturbation will be qualitatively similar to that of the condensed ones, with the only difference that the former are subject to a more intense dissipation. But as we already saw in Sec. VI this dissipation can be enough to distinguish both kind of fluids: while the condensed part might suffer an exponential growth, the uncondensed part might develop low amplitude bounded oscillations.

We have also demonstrated that these resonances show up in the movement of the center of mass as well, causing any initial displacement of the center of mass to be exponentially amplified while the perturbation works. Opposite to our models for the widths, this is an exact prediction based solely on the GPE, and it shows that the parametrical perturbation may also have a disturbing effect in the experiments. On the other hand, a measure of this effect can give us information about the intensity of dissipation and collision effects.

All of the preceding statements are based solely on the GPE. In a few words, they include the existence of resonance regions both for the widths and the center of mass, the shape and the location of those regions, and its intensity as possible measure of damping. The failure of any prediction should be interpreted as a failure of the GPE to describe the condensate. Thus we provide with simple experimental checks to perform a quantitative study of the regimes for which the GPE properly describes the Bose-Einstein condensates in time dependent traps.

Finally, we must mention that throughout this work we have concentrated on regular motion regions in the parameter space. These regions can be “safely” reached in the experiments. There are many other cases where chaos appears in the variational equations and complex behavior is seen in the numerical simulations of Eq. (1). Although the study of those disordered regions could be interesting from the nonlinear science point of view they seem not to be of interest for Bose–Einstein condensation since the exponential separation of nearby orbits which is characteristic of chaotic behavior has been shown to induce instabilities and take the system out of the regime where it can be described using the mean field GP equation .

## ACKNOWLEDGMENTS

This work has been supported in part by the Spanish Ministry of Education and Culture under grants PB95-0389, PB96-0534 and AP97-08930807.