# Extended parametric resonances in Nonlinear Schrödinger systems

## Abstract

We study an example of exact parametric resonance in a extended system ruled by nonlinear partial differential equations of nonlinear Schrödinger type. It is also conjectured how related models not exactly solvable should behave in the same way. The results have applicability in recent experiments in Bose-Einstein condensation and to classical problems in Nonlinear Optics.

Resonances are one of the recurrent themes of Physics. After reading the current PACS sectioning scheme one finds that the word resonance appears in fifty different items which is only a naive way to measure one very important concept in Physics.

When speaking of resonances we use to refer to an “anomalously large response to a (maybe small) external perturbation”. The best known model is the undamped harmonic oscillator driven by an external force, which is textbook material. In this trivial system one easily proves existence of unbounded linearly increasing terms in the solution. A more complicated kind of resonant behavior are parametric resonances, where a relevant parameters of the system is modulated to achieve resonance. This phenomenon has been known since the Middle Ages , and its simplest examples are the parametrically forced harmonic oscillator modeled by Hill’s equation $`\ddot{x}+p(t)x=0`$  and a particular version, $`p(t)=1+ϵ\mathrm{cos}\omega t`$, known as Mathieu’s equation.

Since those are linear problems much can be said on their solutions. For instance, Mathieu’s equation can be studied by means of Floquet’s theory for periodic coefficient equations and resonance existence can be proven analytically -even in the presence of linear dissipation- for various parameter regions on the $`(ϵ,\omega )`$ plane. However resonance phenomena are not restricted to linear problems but also appear in nonlinear finite dimensional problems, e.g. in many Hamiltonian chaotic systems where torus resonances are the reason for the origin of chaos, in impact oscillators and in nuclear magnetic and spin resonances, to cite only a few examples. There is finally the framework of systems modeled by partial differential equations. Here some results are known in simple linear cases, but in general the study of resonances in extended (i.e. ruled by partial differential equations) nonlinear systems poses many open questions.

Our aim in this letter is to prove that unbounded resonances are possible in simple, experimentally realizable Hamiltonian wave problems, specifically in the framework of nonlinear Schrödinger (NLS) equations . This is a very surprising result since one is tempted to think (as it usually happens) that a resonant perturbation can only excite a small number of modes of the infinite present in an extended problem. Since the nonlinearities mix different modes the amplification of one of them is stopped by energy transfer to the others, specially in conservative systems where the conservation laws control the number and amplitude of active modes. Thus, it could be thought that once the energy is transfered to the non-resonant modes and being their growth limited by nonlinear constraints, the action of the perturbation should be controlled and unbounded growth of the relevant quantities should be inhibited. Indeed this is the common mechanism in many systems. We will show that it is not the only possibility and truly resonant behavior is possible even in simple physically relevant models.

The model.- Let us consider the following Nonlinear Schrödinger equation in $`n`$ spatial dimensions.

$$i\frac{u}{t}=\frac{1}{2}\mathrm{}u+|u|^2u+V(𝐫,t)u,$$
(1)   
where $`\mathrm{}=_{j=1}^n^2/x_j^2`$. This equation is the adimensional description of many different phenomena . One of the hottest topics where it appears is in the mean field model of nonuniform trapped Bose-Einstein condensates (BEC)  where $`V(𝐫)=\frac{1}{2}_{j=1}^n\lambda _j(t)x_j^2`$. Time dependent coefficients reflect experiments in which the trap is perturbed to obtain the response spectrum of the condensate  (For a theoretical analysis of these experiments see Ref. ). A more detailed study of the condensate response to time dependent perturbations was performed in  and the existence of strong resonances was proposed on the ground of approximate variational methods and numerical simulations of the Eq. (1). Nonetheless, the nature of the resonance could not be completely assured on the ground of approximate or numerical techniques since both have limited applicability, specially for the strong response phenomena considered here.

Let us first consider in this letter the two-dimensional case with time-dependent parabolic potential

$$i\frac{u}{t}=\frac{1}{2}\mathrm{}u+|u|^2u+\frac{1}{2}\lambda (t)\left(x^2+y^2\right)u,$$
(2)   
where $`\mathrm{}=^2/x^2+^2/y^2`$. Although Bose-Einstein condensation problems are usually three dimensional, this model can be related to pancake type traps, and to recent experiments with quasi two-dimensional condensates or coherent atomic systems . The model (2) is also known in nonlinear Optics where it is used to study the propagation of paraxial beams in fibers with a (modulated) parabolic profile of the refraction index.

From the mathematical point of view, Eq. (2) is a parametrically forced NLS equation that has not been studied rigorously before. Related systems are the externally driven damped NLS in one dimension , the AC-driven damped sine-Gordon system  and one type of non-self-adjoint parametrically driven NLS . However those works concentrate on one dimensional equations and mostly on soliton stability problems.

Moment method.- We will first study the radially symmetric version of Eq. (2), searching solutions of the form $`\psi (r,\theta ,t)=u(r,t)e^{im\theta }`$, which includes both the typical radially symmetric problem corresponding to $`m=0`$, and vortex line solutions, with $`m0`$. The simplified equation for $`u`$ is

$$i\frac{u}{t}=\frac{1}{2r}\frac{}{r}\left(r\frac{u}{r}\right)+\left(\frac{m^2}{2r}+|u|^2+\frac{\lambda (t)}{2}r^2\right)u,$$
(3)   
Eq. (3) is non-integrable and has no exact solutions even in the constant $`\lambda `$ case. Its solutions are stationary points of the action $`S=_{t_0}^{t_1}(t)`$, where

$`(t)`$ $`=`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \left(u\frac{u^{}}{t}u^{}\frac{u}{t}\right)d^2x}+{\displaystyle \frac{\lambda (t)}{2}r^2|u|^2d^2x}`$ (4)   
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\left|\frac{u^{}}{r}\right|^2+\frac{m^2}{r^2}|u|^2+|u|^4\right)d^2x}.`$ (5)   
Let us define the following integral quantities 

$`I_1(t)`$ $`=`$ $`{\displaystyle |u|^2d^2x},`$ (7)   
$`I_2(t)`$ $`=`$ $`{\displaystyle |u|^2r^2d^2x},`$ (8)   
$`I_3(t)`$ $`=`$ $`i{\displaystyle \left[u\frac{u^{}}{r}u^{}\frac{u}{r}\right]rd^2x},`$ (9)   
$`I_4(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\left|u\right|^2+\frac{m^2}{r^2}|u|^2+|u|^4\right)d^2x},`$ (10)   
where $`d^2x=2\pi rdr`$ because of the symmetry. These magnitudes are related physically to the norm (intensity or number of particles), width, radial momentum and energy of the wave packet. In the Optical interpretation of Eq. (3) these quantities are known as moments and are used in (usually approximate) calculations related to beam parameters evolution . It is remarkable that the $`I_j`$ satisfy simple, and most important, closed evolution laws

$`{\displaystyle \frac{dI_1}{dt}}`$ $`=`$ $`0,`$ (12)   
$`{\displaystyle \frac{dI_2}{dt}}`$ $`=`$ $`I_3,`$ (13)   
$`{\displaystyle \frac{dI_3}{dt}}`$ $`=`$ $`2\lambda (t)I_2+4I_4,`$ (14)   
$`{\displaystyle \frac{dI_4}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\lambda (t)I_3.`$ (15)   
The first equation comes from the phase invariance of Eq. (3) under global phase transformations and corresponds to the $`L^2`$-norm conservation (in BEC it is interpreted as the particle number conservation in the mean field model). The other equations can also be obtained in connection with the invariance of the action (4) under symmetry transformations.

Singular Hill’s equation.- Eqs. (Extended parametric resonances in Nonlinear Schrödinger systems) form a linear non-autonomous system for the unknowns, $`I_j,j=1,\mathrm{},4`$ that has several positive invariants under time evolution, of which the most important is

$$Q=2I_4I_2I_3^2/4>0.$$
(16)   
With the help of this quantity, the system (Extended parametric resonances in Nonlinear Schrödinger systems) can be reduced to a single equation for the most relevant parameter $`I_2(t)`$, which is

$$\frac{d^2I_2}{dt^2}\frac{1}{2I_2}\left(\frac{dI_2}{dt}\right)^2+2\lambda (t)I_2=\frac{Q}{I_2}.$$
(17)   
If we were able to solve Eq. (17) then the use of Eqs. (Extended parametric resonances in Nonlinear Schrödinger systems) would allow us to track the evolution of the other ones. We can do so by defining $`X(t)=\sqrt{I_2}`$, whose physical meaning is the wavepacket width, and substituting it into (17). This procedure gives us

$$\ddot{X}+\lambda (t)X=\frac{Q}{X^3}.$$
(18)   
The resulting equation is a singular (nonlinear) Hill equation. In general it cannot be solved explicitly (see Ref.  for other examples of these equations), however in this case we are fortunate that the general solution of Eq. (18) can be obtained  and is given by

$$X(t)=\sqrt{u^2(t)+\frac{Q}{W^2}v^2(t)},$$
(19)   
where $`u(t)`$ and $`v(t)`$ are the two linearly independent solutions of the equation

$$\ddot{x}+\lambda (t)x=0,$$
(20)   
which satisfy $`u(t_0)=X(t_0),\dot{u}(t_0)=X^{}(t_0),v(t_0)=0,v^{}(t_0)0`$, and $`W`$ is the Wronskian $`W=u\dot{v}\dot{u}v=\text{const.}0.`$

The conclusion is that the evolution of the width of the wave packet, which is given by (19) is closely related to the solutions of the Hill equation (20). This is a very well studied problem  which is explicitly solvable only for particular choices of $`\lambda (t)`$, but whose solutions are well characterized and many of its properties are known.

A practical application of Eq. (19) is that one may design $`\lambda (t)`$ starting from the desired properties for the wave packet width evolution. This can be done in BEC applications by controlling the trapping potential and in Optics by the precise control of the $`z`$dependent refraction index of the waveguide.

If it is supposed that $`\lambda (t)`$ depends on a parameter $`\lambda (t)=1+\stackrel{~}{\lambda }(t)`$ with $`\stackrel{~}{\lambda }(t)`$ a periodic function with zero mean value and peak value $`ϵ`$ (not necessarily small), then there is a complete theory describing intervals of $`ϵ`$ for which all solutions of Eq. (20) are bounded (stability intervals) and intervals for which all solutions are unbounded (instability intervals). Both types of intervals are ordered in a natural way. Let us finally concentrate in the physically relevant case $`\lambda (t)=1+ϵ\mathrm{cos}\omega t`$, and where $`ϵ`$ needs not to be a perturbative parameter. In that case the parameter regions where exact resonances occur can be studied by several means. First, for any fixed $`ϵ`$, there exist infinite ordered sequences $`\{\omega _n\},\{\omega _n^{}\}`$ tending to zero such that Eq. (20) is resonant if $`\omega `$ belongs to $`(\omega _n,\omega _n^{})`$ for some $`n`$. Second, when $`\omega `$ is fixed, resonances appear for $`ϵ`$ big enough. Further, a stability diagram can be drawn in the $`ϵ\omega `$ plane as shown in Fig. 1. The boundaries of these regions are the so-called characteristic curves, which can not be explicitly derived but whose existence can be proven analytically: if $`D(ϵ,\omega )`$ is the discriminant of the equation, characteristic curves are obtained by solving the equations $`D(ϵ,\omega )=2`$ and $`D(ϵ,\omega )=2`$. In particular, instability regions start on the frequencies $`\omega =2,1,1/2,\mathrm{},2/n^2,\mathrm{}`$  \[Fig. 1\].

The resonant behavior depends only on the mutual relation of the parameters but not on the initial data. Since a resonance in Eq. (20) automatically implies unbounded behavior in the solution \[Eq. (19)\] we obtain the result that for resonant parameters all solutions have divergent width provided they start from finite values of the moments (which is the generic case). With respect to stability, it is deduced from Massera’s Theorem that if $`(ϵ,\omega )`$ belongs to a stability region of Eq. (20), there exists a periodic solution of Eq. (18) which is Lyapunov stable.

We emphasize that our analysis of the cilindrically symmetric problem performed up to now is completely rigorous. Now we will turn to approximations to discuss other related problems.

Non-symmetric problems.- Throughout this paper we have made use of a special symmetry for the confining potential and for the solution in order to make our theory exact. When these constrains are removed and we turn to a general two or three-dimensional model (1), a corresponding set of coupled Hill’s equations may still be obtained by some kind of approximation, which can be scaling laws  or a variational ansatz .

A new approach involves only imposing the wave-function to have a quadratic dependence in the complex phase $`\psi =|\psi |\text{exp}\left(i_{jk}\alpha _{jk}x_jx_k\right)`$, $`j,k=1,\mathrm{},3`$, which is a less restrictive ansatz than those used previously. By introducing this trial function into (4), building Lagrange’s equations for the parameters, and using some transformations, one obtains

$$\ddot{X_i}+\lambda _i(t)X_i=\frac{1}{X_i^3}+\frac{Q}{X_iX_1X_2X_3}.$$
(21)   
Here $`X_i,i=1,2,3`$ are the three root mean square radius of the solution for each spatial direction, i. e. the extensions of the previous $`X(t)`$ to the three dimensional problem. In this case the $`\lambda _i(t)`$ need not to be equal and in fact experimental problems in BEC involve different $`\lambda `$ values because of asymmetries of the traps.

Eqs. (21) are not integrable and form a six-dimensional nonautonomous dynamical system for which few things can be said analytically. Nevertheless, the numerical study of the approximate model (21) exhibits an extended family of resonances which is more or less the Cartesian product of those of (19), with minor displacements due to the coupling. There exist also parameter regions of chaotic and periodic solutions. Numerical simulations of Eq. (1) confirm the predictions of the simple model (21).

Nonlinear spectrum and resonances.- Although resonant behavior has been proven in the radially symmetric two dimensional case, one could try to make a qualitative explanation for the reason that this behavior is also present in non-symmetric problems. Let us look for stationary solutions $`u(𝐫,t)=\varphi (𝐫)e^{i\mu t}`$ for the stationary trap, $`V=V(𝐫)`$, and write the nonlinear eigenvalue problem

$$\mu \varphi =\frac{1}{2}\mathrm{}\varphi +|\varphi |^2\varphi +V(𝐫)\varphi .$$
(22)   
Let us assume that we may use these solutions to expand, in a possibly non-orthogonal and time dependent way, any solution of the non-stationary problem (1). By doing so one finds that the energy absorption process in the time dependent potential is ruled by the separation between the eigenvalues, $`\mu _i\mu _j`$, of any two different modes of (22). If these differences are distributed in a random way, the perturbation is not efficient and does not lead to resonances, but to chaos. On the other hand, if the differences may be approximated by multiples of a fixed set of frequencies, then an appropriate parametric excitation will induce a sustained process of energy gain (and width growth) such as the one we have observed.

We have studied the spectrum $`\{\mu \}`$ for the case of an axially symmetric potential $`V(𝐫)=\lambda _r(x^2+y^2)+\lambda _zz^2`$ in three dimensions using a variant of a pseudospectral scheme using the harmonic oscillator basis as described in . Our results show that the spectrum exhibits an ordered structure, with different directions of uniformity that may be excited by the parametric perturbation. We have checked the computations in one, two and three spatial dimensions without assumptions of symmetry on the solution with similar results despite the fact that the spectrum becomes more complex as dimensionality increases. Details of these calculations will be published elsewhere.

Losses effects.- When losses are included in Eq. (3), e.g. by the addition of a new term of the type $`i\sigma u`$, it is not possible to obtain a set of closed equations for the momenta and thus one cannot solve exactly the problem. However, as discussed in the preceding comments one still may further restrict the analysis to the parabolic phase approximation to obtain the modified equation $`\ddot{X}+\lambda (t)X=Q(t)/X^3`$, where $`Q(t)`$ is a decreasing function satisfying $`Q1`$ as $`t\mathrm{}`$. At least up to the range of validity of the approximation one obtains parametric resonances since the behavior of the singular term does not affect essentially the resonance relation. Of course, at the same time one has a decrease in the norm of the solution since now $`dI_1/dt=2\sigma I_1`$.

In conclusion we have studied a parametrically perturbed nonlinear extended system. By using the moment technique for the cilindrically symmetric problem we obtain a singular Hill equation which is reducible to a linear Hill equation and thus existence of resonances can be shown analytically. For the physically relevant case of a periodic perturbation it is shown that there exist strong extended resonances for the relevant parameters of the solution even when the solution is constrained by conservation laws. It is conjectured using the parabolic ansatz, analysis of nonlinear spectrum and numerical simulations that this behavior is also present in non-radially symmetric problems or when losses are present.

V.M.P-G. has been partially supported by CICYT grants PB96-0534 and PB95-0839. P. T. is supported by CICYT grant PB95-1203.