# Nonlinear Modes of a Macroscopic Quantum Oscillator

## Abstract

We consider the Bose-Einstein condensate in a parabolic trap as a macroscopic quantum oscillator and describe, analytically and numerically, its collective modes–a nonlinear generalisation of the (symmetric and antisymmetric) Hermite-Gauss eigenmodes of a harmonic quantum oscillator.

The recent observation of different types of Bose-Einstein condensation (BEC) in atomic clouds  led to the foundation of a new direction in the study of macroscopic quantum phenomena. From a general point of view, the dynamics of gases of cooled atoms confined in a magnetic trap at very low temperatures can be described by an effective equation for the condensate wave function known as the Gross-Pitaevskii (GP) equation . This is a classical nonlinear equation that takes into account the effects of the particle interaction through an effective mean field, and therefore it can be treated as a nonlinear generalization of a text-book problem of quantum mechanics, i.e. as a macroscopic quantum oscillator.

Similar models of the confined dynamics of macroscopic quantum systems appear in other fields, e.g. in the case of an electron gas confined in a quantum well , or optical modes in a photonic microcavity . In all such systems, confined single-particle states are restricted to a set of discrete energies that form a set of eigenmodes. A classical and probably most familiar example of such a system is a harmonic quantum oscillator with equally spaced energy levels .

When, instead of single-particle states, we describe quasiclassically a system of interacting bosons in a macroscopic ground state confined in an external potential, a standard application of the mean-field theory allows us to introduce a macroscopic wave function as a classical field $`\mathrm{\Psi }(𝐑,t)`$ having the meaning of the order parameter. The equation for the function $`\mathrm{\Psi }(𝐑,t)`$ looks similar to that of a single-particle oscillator, but it also includes the effect of interparticle interaction, taken into account as a mean-field nonlinear term. Then, the important questions are: Does the physical picture of eigenmodes remain valid in the nonlinear case, and what is the effect of nonlinearity on the modes ? In this paper we analyse nonlinear eigenmodes of a macroscopic quantum oscillator as a set of nonlinear stationary states that extend the well-known Hermite-Gauss eigenfunctions. We also make a link between seemingly different approximations, the well-known Thomas-Fermi approximation and the perturbation theory developed here for the case of weak nonlinearity. For both attractive and repulsive interaction, we demonstrate a close connection between the nonlinear modes and (bright and dark) multi-soliton stationary states.

We consider the macroscopic dynamics of condensed atomic clouds in a three-dimensional, strongly anisotropic, external parabolic potential created by a magnetic trap. The BEC collective dynamics can be described by the GP equation,

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V(𝐑)\mathrm{\Psi }+U_0|\mathrm{\Psi }|^2\mathrm{\Psi },$$
(1)   
where $`\mathrm{\Psi }(𝐑,t)`$ is the macroscopic wave function of a condensate, $`V(𝐑)`$ is a parabolic trapping potential, and the parameter $`U_0=4\pi \mathrm{}^2(a/m)`$ characterises the two-particle interaction proportional to the s-wave scattering length $`a`$. When $`a>0`$, the interaction between the particles in the condensate is repulsive, whereas for $`a<0`$ the interaction is attractive. In fact, the scattering length $`a`$ can be continuously detuned from positive to negative values by varying the external magnetic field near the so-called Feshbach resonances .

First of all, we derive from Eq. (1) an effective one-dimensional model, assuming the case of a highly anisotropic (cigar-shaped) trap of the axial symmetry $`V(𝐑)=\frac{1}{2}m\omega _{}^2(R_{}^2+\lambda X^2)`$, where $`R_{}=\sqrt{Y^2+Z^2}`$. This means that $`\lambda \omega _{}^2/\omega _{}^21`$, and the transverse structure of the condensate, being close to a Gaussian in shape, is mostly defined by the trapping potential .

Measuring the spatial variables in the units of the longitudinal harmonic oscillator length $`a_{ho}=(\mathrm{}/m\omega \sqrt{\lambda })^{1/2}`$, and the wavefunction amplitude, in units of $`(\mathrm{}\omega _{}/2U_0\sqrt{\lambda })^{1/2}`$, we obtain the following dimensionless equation:

$$i\frac{\mathrm{\Psi }}{t}+^2\mathrm{\Psi }[\lambda ^1(y^2+z^2)+x^2]\mathrm{\Psi }+\sigma |\mathrm{\Psi }|^2\mathrm{\Psi }=0,$$
(2)   
where time is measured in the units of $`(2/\omega _{}\sqrt{\lambda })`$, $`(x,y,z)=(X,Y,Z)/a_{ho}`$, and the sign $`\sigma =\mathrm{sgn}(a)=\pm 1`$ in front of the nonlinear term is defined by the sign of the s-wave scattering length of the two-body interaction.

We assume that in Eq. (1) the nonlinear interaction is weak relative to the trapping potential force in the transverse dimensions, i.e. $`\lambda 1`$. Then, it follows from Eq. (2) that the transverse structure of the condensate is of order of $`\lambda `$, and the condensate has a cigar-like shape. Therefore, we can look for solutions of Eq. (2) in the form,

$`\mathrm{\Psi }(r,x,t)=\mathrm{\Phi }(r)\psi (x,t)e^{2i\gamma t},`$

where $`r=\sqrt{y^2+z^2}`$, and $`\mathrm{\Phi }(r)`$ is a solution of the auxiliary problem for the 2D radially symmetric quantum harmonic oscillator

$`_{}^2\mathrm{\Phi }+2\gamma \mathrm{\Phi }(r^2/\lambda )\mathrm{\Phi }=0,`$

which we take in the form of the no-node ground state, $`\mathrm{\Phi }_0(r)=C\mathrm{exp}(\gamma r^2/2)`$, where $`\gamma =1/\sqrt{\lambda }`$. To preserve all the information about the structure of the 3D condensate in an asymmetric trap describing its properties by the effective GP equation for the longitudinal profile, we impose the normalisation for $`\mathrm{\Phi }_0(r)`$ that yields $`C^2=\gamma /\pi `$.

After substituting such a factorized solution into Eq. (2), dividing by $`\mathrm{\Phi }`$ and integrating over the transverse cross-section of the cigar-shaped condensate, we finally derive the following 1D nonstationary GP equation

$$i\frac{\psi }{t}+\frac{^2\psi }{x^2}x^2\psi +\sigma |\psi |^2\psi =0.$$
(3)   
The number of the condensate particles $`N`$ is now defined as $`N=(\mathrm{}\omega /2U_0\sqrt{\lambda })Q`$, where

$$Q=_{\mathrm{}}^{\mathrm{}}|\psi (x,t)|^2𝑑x$$
(4)   
is the integral of motion for the normalised nonstationary GP equation (3).

Equation (3) includes all the terms of the same order, and it describes a longitudinal profile of the condensate state in a highly anisotropic trap. In the linear limit, i.e. when formally $`\sigma 0`$, Eq. (3) becomes the well-known equation for a harmonic quantum oscillator. Its stationary localised solutions,

$$\psi (x,t)=\varphi (x)e^{i\mathrm{\Omega }t},$$
(5)   
exist only for discrete values of $`\mathrm{\Omega }`$, such that

$`\mathrm{\Omega }_n=1+2n,n=0,1,2,\mathrm{},`$

and they are defined through the Hermite-Gauss polynomials, $`\varphi _n(x)=c_n\mathrm{e}^{x^2/2}H_n(x)`$, where $`c_n=(2^nn!\sqrt{\pi })^{1/2}`$, and

$$H_n(x)=(1)^n\mathrm{e}^{x^2/2}\frac{d^n(\mathrm{e}^{x^2/2})}{dx^n},$$
(6)   
so that $`H_0=1`$, $`H_1=2x`$, etc.

In general, the localised solutions of Eq. (3) for $`\sigma 0`$ can be found only numerically. All such solutions can be characterised by the dependence of the invariant (4) on the effective nonlinear frequency $`\mathrm{\Omega }`$. However, in some particular limits we can employ different approximate methods to find the localised solutions analytically.

First of all, to describe the effect of weak nonlinearity, we use the perturbation theory based on the expansion of the general solution of Eq. (3) in the infinite set of the eigenfunctions (6). A similar approach has been used earlier in the theory of the dispersion-managed optical solitons . To apply such a perturbation theory, we look for solutions of Eq. (3) in the form \[cf. Eq. (5)\]

$$\psi (x,t)=e^{i\mathrm{\Omega }t}\underset{n=0}{\overset{\mathrm{}}{}}B_n\varphi _n(x),$$
(7)   
where $`\varphi _n(x)`$ are the eigenfunctions of the linear equation for a harmonic oscillator that satisfy the equation

$$\frac{d^2\varphi _n}{dx^2}x^2\varphi _n+\mathrm{\Omega }_n\varphi _n=0.$$
(8)   
Inserting the expansion (7) into Eq. (3), multiplying by $`\varphi _n`$ and averaging, we obtain a system of algebraic equations for the coefficients

$$(\mathrm{\Omega }\mathrm{\Omega }_m)B_m\sigma \underset{n,l,k}{}V_{m,n,l,k}B_nB_lB_k=0,$$
(9)   
where

$`V_{m,n,l,k}={\displaystyle _{\mathrm{}}^+\mathrm{}}\varphi _m(x)\varphi _n(x)\varphi _l(x)\varphi _k(x)𝑑x.`$

Equation (9) can be also rewritten in the traditional form $`\delta (H+\mathrm{\Omega }Q)=0`$, and it allows us to develop a perturbation theory for small nonlinearities, for any Hermite-Gaussian eigenmode. For example, let us consider the ground state mode at $`n=0`$. Assuming $`B_0B_m`$ for $`m0`$ and the condition of the symmetric solution $`\varphi (x)=\varphi (x)`$, i.e. $`B_{2k+1}=0`$ for any $`k`$, we find the corrections, $`\mathrm{\Omega }\mathrm{\Omega }_0+\sigma V_{0,0,0,0}|B_0|^2`$, and

$`B_{2k}={\displaystyle \frac{\sigma V_{2k,0,0,0}}{(\mathrm{\Omega }\mathrm{\Omega }_{2k})}}|B_0|^2B_0,k0.`$

This allows us to calculate the asymptotic expansion of the invariant $`Q`$ for small nonlinearities,

$$Q=\underset{k}{}|B_{2k}|^2\sigma a_0(\mathrm{\Omega }\mathrm{\Omega }_0)[1+b_0(\mathrm{\Omega }\mathrm{\Omega }_0)^2],$$
(10)   
where the coefficients are:

$`a_0=\sqrt{2\pi },b_0={\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(2k)!}{(4k)^2(k!\mathrm{\hspace{0.17em}2}^{2k})^2}}.`$

Higher-order modes can be considered in a similar way, and the results are similar to Eq. (10) where $`a_0`$ and $`b_0`$ change to $`a_n`$ and $`b_n`$ respectively, where $`a_n`$ and $`b_n`$ depend on the mode order.

In the opposite limit, i.e. when the nonlinear term or potential are large in comparison with the kinetic term given by the second-order derivative, we can use two different approximations for describing localized modes. For $`\sigma =+1`$ and large negative $`\mathrm{\Omega }`$, localised modes are described by the stationary solutions of the nonlinear Schrödinger (NLS) equation, that appears when we neglect the trapping potential. The NLS one-soliton solution is $`\varphi _s(x)=\sqrt{2\mathrm{\Omega }}\mathrm{sech}(x\sqrt{\mathrm{\Omega }})`$, so that the dependence $`Q(\mathrm{\Omega })`$ coincides with the soliton invariant $`Q_s=4\sqrt{\mathrm{\Omega }}`$. For $`\sigma =1`$ and large positive $`\mathrm{\Omega }`$, the ground-state solution can be obtained by using the so-called Thomas-Fermi approximation, based on neglecting the kinetic term - this yields $`\varphi _{\mathrm{TF}}(x)\sqrt{\mathrm{\Omega }x^2}`$.

In general, we should solve Eq. (3) numerically. Figures 1(a) to 1(b) present examples of the numerically found ground-state solutions of Eq. (3) in the form (5) as continuous functions of the dimensionless parameter $`\mathrm{\Omega }`$, for both negative (a) and positive (b) scattering length. For $`\mathrm{\Omega }1`$, i.e. in the limit of the harmonic oscillator ground-state mode, the solution is close to Gaussian for both the cases. When $`\mathrm{\Omega }`$ deviates from 1, the solution profile is defined by the type of nonlinearity. For attraction ($`\sigma =+1`$), the profile approaches the sech-type soliton, whereas for repulsion ($`\sigma =1`$) the solution flattens, and it is better described by the Thomas-Fermi approximation, that is good except at the edge points.

In Fig. 2 we present the dependence of the invariant $`Q`$ on the parameter $`\mathrm{\Omega }`$, for both the types of the ground-state solution, corresponding to two different signs of the scattering length. The dashed curve for the zero-order mode (marked as 0th-mode) shows the dependence $`Q_\mathrm{s}(\mathrm{\Omega })`$ for the soliton solution of the NLS equation without a trapping potential. In the asymptotic region of negative $`\mathrm{\Omega }`$, i.e. say for $`\mathrm{\Omega }<2`$, the dependence $`Q(\mathrm{\Omega })`$ for the BEC condensate in a trap approaches the curve $`Q_\mathrm{s}(\mathrm{\Omega })`$. This means that for such a narrow localised state the effect of a parabolic potential is negligible, and the condensate ground state becomes localised mostly due to an attractive interparticle interaction. In contrast, for large positive $`\mathrm{\Omega }`$ the effect of a trapping potential is crucial, and the solution of the Thomas-Fermi approximation, $`\varphi _{\mathrm{TF}}(x)=\sqrt{\mathrm{\Omega }x^2}`$, defines the common asymptotics for all the modes, $`Q_{\mathrm{TF}}\frac{4}{3}\mathrm{\Omega }^{3/2}`$ (dashed-dotted curves in Fig. 2).

As has been mentioned above, in the linear limit ($`\sigma 0`$), Eq. (3) possesses a discrete set of localised modes described by the Hermite-Gauss polynomials. We have demonstrated that all such modes can be readily calculated by the perturbation theory in the weakly nonlinear approximation, and therefore they should all exist for the nonlinear problem as well, describing an analytical continuation of the Hermite-Gauss linear modes to a set of nonlinear stationary states. In application to the BEC theory, these non-ground-state solutions were first discussed by Yukalov et al .

Figures 1(c) to 1(f) show examples of the first- and second-order modes for both negative and positive scattering length, respectively. In the limit $`\mathrm{\Omega }\mathrm{\Omega }_n`$, all those modes transform into the corresponding eigenfunctions of a linear harmonic oscillator.

It is clear that nonlinearity has a different effect for the negative and positive scattering length. For the negative scattering length (attraction), the higher-order modes transform into multi-soliton states consisting of a sequence of solitary waves with alternating phases \[see Figs. 1(c) and 1(e)\]. This is further confirmed by the analysis of the invariant $`Q`$ vs. $`\mathrm{\Omega }`$, where all the branches of the higher-order modes approach asymptotically the soliton dependencies $`Q_n(n+1)Q_\mathrm{s}`$, where $`n`$ is the order of the mode ($`n=0,1,\mathrm{}`$). From the physical point of view, in the case of attractive interaction the higher-order stationary modes exist due to a balance between repulsion of out-of-phase bright NLS solitons and attraction imposed by the trapping potential. The analysis of the global stability of such higher-order multihump multi-soliton modes is still an open problem, however the recent results indicate that, at least in some nonlinear models, multihump soliton states can be stable .

For the positive scattering length ($`\sigma =1`$), the higher-order modes transform into a sequence of dark solitons (or kinks) , so that the first-order mode corresponds to a single dark soliton, the second-order mode, to a pair of dark solitons, etc. \[see Figs. 1(d) and 1(f)\]. Again, these stationary solutions satisfy a force balance condition \- repulsion between dark solitons is exactly compensated by an attractive force of the trapping potential.

The modal structure of the condensate macroscopic states described above and summarised in Fig. 2 for both positive and negative values of the scattering length allows us to draw an analogy between BEC in a trap and the guided-wave optics where the condensate dynamics in time corresponds to the stationary mode propagation along an optical waveguide, with the parameter $`\mathrm{\Omega }`$ as the propagation constant. As is well known from different problems of guided-wave optics, in the presence of interaction the guided modes become coupled and the coupling can lead to both power exchange and nonlinear phase shifting between the modes. In application to the BEC theory, the mode coupling resembles a kind of internal Josephson effect. These issues are beyond the scope of this paper and will be analysed elsewhere .

The theory of nonlinear stationary modes of a macroscopic quantum oscillator developed above for the 1D analog of BEC can be easily extended to both 2D and 3D cases. Moreover, the coupled-mode theory for a single condensate is closely connected to the dynamics of strongly coupled two-component BECs  where excitation of an antisymmetric (or, in our notation, the first order) collective mode, in the form of collapses and revivals, has been recently observed experimentally .

At last, we would like to mention that the basic concepts and results presented above can find their applications in other fields. For example, the effect of nonlinearity can lead to a mode coupling in the so-called “photonic atom”, a micrometer-sized piece of semiconductor that traps photons inside , or two such ‘photonic atoms’ coupled together, a ‘photonic molecule’ . In such photonic microcavity structures, the macroscopic nature of the states may lead to different nonlinear effects including the mode mixing and power exchange.

In conclusion, we have analysed nonlinear stationary modes of a macroscopic quantum oscillator considering, as an example, the cigar-shaped Bose-Einstein condensate in a parabolic trap.