# The moment method in general Nonlinear Schrödinger Equations

## Abstract

In this paper we develop a new approximation method valid for a wide family of nonlinear wave equations of Nonlinear Schrödinger type. The result is a reduced set of ordinary differential equations for a finite set of parameters measuring global properties of the solutions, named momenta. We prove that these equations provide exact results in some relevant cases and show how to impose reasonable approximations that can be regarded as a perturbative approach and as an extension of the time dependent variational method.

The solution of nonlinear wave equations representing physically relevant phenomena is a task of the highest interest. However, it is not possible to find exact solutions except for a few simple cases where one is lucky to integrate the equations involved. In particular, in the 70’s some mathematical techniques were discovered which allowed the integration of several relevant nonlinear wave equations . So the development of rigorous approximation methods is of interest in those cases where the equations are known to be non-integrable.

One family of nonlinear wave equations with lots of practical applications is that of Nonlinear Schrödinger Equations (NSE) , which arise in plasma physics, biomolecule dynamics, fundamentals of quantum mechanics, beam physics, etc., but specially in the fields of Nonlinear Optics  and Bose–Einstein condensation . In the last two fields a great variety of these equations appear involving different spatial dimensionalities, nonlinear terms (saturation, polynomial, nonlocal, losses, etc.) and number of coupled equations.

One common approximate theoretical approach to the analysis of the dynamics involved in those problems is to assume a fixed shape for the solution with a finite set of time dependent parameters so that the dimensionality is made finite at the cost of loss of information on the solution (many degrees of freedom are lost). This method receives many denominations depending on the context: collective coordinate technique, time-dependent variational method, equivalent particle approach, energy balance equations, etc.

Although not explicitly stated most of those methods can be reduced to a more elegant formulation which is the time dependent variational technique, originally developed by Anderson  for one dimensional problems based on Ritz’s optimization procedure. This approximate technique is a good tool to study the propagation of distributions having simple shape. If the shape of the actual solution is close to the trial function, the outcome of variational method will be in good agreement with the real solutions, otherwise it may be very rough or even fail . Despite of this fact the technique has been used in many physical situations. In Nonlinear Optics it has been applied to many problems some of them being listed in Ref. . The method has been applied to many other physical problems where nonlinear wave equations (in particular NSEs) arise including random perturbations , nonlocal equations , collapse phenomena , propagation and scattering of nonlinear waves , etc. A review of the application of the technique with emphasis on problems with different scales (focused in condensed matter) is given in Ref. . This technique has been also used in the last years in the framework of Bose-Einstein condensation (BEC) applications to explain the low energy excitation spectrum of single  and double  condensates, collapse dynamics , and many other problems , etc.

In this letter we develop a completely different technique called the moment method . This method is based on the definition of several integral parameters whose evolution can be computed in closed form and has been used to obtain exact results in particular applications . We provide here a general framework for its application as well as several ways to treat it systematically as a perturbative technique.

The general NSE and moment equations.- Let us consider the $`n`$-dimensional NSE

$$i\frac{\psi }{t}=\frac{1}{2}\mathrm{}\psi +V(\stackrel{}{r})\psi +g(|\psi |^2,t)\psi i\sigma (|\psi |^2,t)\psi $$
(1)   
where $`\mathrm{}=_k^2/x_j^2`$, and $`V(\stackrel{}{r})=_k\frac{1}{2}\omega _k(t)^2(t)x_k^2`$ is a time and spatially dependent parabolic potential which has been included because it is always present in BEC problems. To aid readability, we will separate the solution in modulus and phase, $`\psi =\sqrt{\rho }e^{i\varphi }`$, and define an interaction energy density $`G(\rho )`$ as $`g=G/\rho `$. The nonlinear terms will be analytic functions of $`\rho `$ such that $`g(\rho ),G(\rho )0,\rho 0.`$

In Table I we define the so called momenta of $`\psi `$. Some of them are related to the momenta of the distribution $`\rho =|\psi |^2`$, and they all have a physical meaning (see Ref.  for their interpretation in Optics). The evolution equations for these quantities are

$`{\displaystyle \frac{dN}{dt}}`$ $`=`$ $`2{\displaystyle \sigma \rho },`$ (3)   
$`{\displaystyle \frac{dX_i}{dt}}`$ $`=`$ $`V_i2{\displaystyle \sigma x_i\rho },`$ (4)   
$`{\displaystyle \frac{dV_i}{dt}}`$ $`=`$ $`\omega _iX_i2{\displaystyle \sigma \frac{\varphi }{x_i}\rho },`$ (5)   
$`{\displaystyle \frac{dW_i}{dt}}`$ $`=`$ $`B_i2{\displaystyle \sigma x_i^2\rho },`$ (6)   
$`{\displaystyle \frac{dB_i}{dt}}`$ $`=`$ $`4K_i2\omega _i^2W_i2{\displaystyle DG}4{\displaystyle \sigma x_i\frac{\varphi }{x_i}\rho },`$ (7)   
$`{\displaystyle \frac{dK_i}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\omega _i^2B_i{\displaystyle DG\frac{^2\varphi }{x_i^2}}`$ (8)   
$`+`$ $`{\displaystyle \sigma \left[\sqrt{\rho }\frac{^2\sqrt{\rho }}{x_i^2}\rho \left(\frac{\varphi }{x_{\widehat{i}}}\right)^2\right]},`$ (9)   
$`{\displaystyle \frac{dJ}{dt}}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle DG\frac{^2\varphi }{x_i^2}}2{\displaystyle \sigma g\rho }+{\displaystyle \frac{G}{t}}.`$ (10)   
Here $`DG`$ is a shorthand notation for $`G(\rho )g(\rho )\rho `$. Through this paper we will concentrate in the most common case $`\sigma (|\psi |^2,t)=\sigma (t)`$ for which the equations are

$`{\displaystyle \frac{dN}{dt}}`$ $`=`$ $`2\sigma N,`$ (12)   
$`{\displaystyle \frac{dX_i}{dt}}`$ $`=`$ $`V_i2\sigma X_i,`$ (13)   
$`{\displaystyle \frac{dV_i}{dt}}`$ $`=`$ $`\omega _iX_i2\sigma V_i,`$ (14)   
$`{\displaystyle \frac{dW_i}{dt}}`$ $`=`$ $`B_i2\sigma W_i,`$ (15)   
$`{\displaystyle \frac{dB_i}{dt}}`$ $`=`$ $`4K_i2\omega _i^2W_i2{\displaystyle DG}2\sigma B_i,`$ (16)   
$`{\displaystyle \frac{dK_i}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\omega _i^2B_i{\displaystyle DG\frac{^2\varphi }{x_i^2}}2\sigma K_i,`$ (17)   
$`{\displaystyle \frac{dJ}{dt}}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle DG\frac{^2\varphi }{x_i^2}}2J+{\displaystyle \frac{G}{t}}.`$ (18)   
As stated before, some of these laws can be found in other treatments  which mostly concentrate on particular cases or treat them as basis for perturbation methods choosing one particular shape for the solution. From Eqs. (13-14) we find *exact* closed equations for zeroth and first order momenta,

$$\frac{d^2X_i}{dt^2}=\omega _i(t)X_i2\sigma (t)\frac{dX_i}{dt}2\sigma (t)X_i.$$
(19)   
In conservative systems with any type of potential the classical result of Quantum Mechanics $`dx_i/dt=V/x_i`$ is obtained as described in .

Once Eqs. (12-14) are integrated one is left with the problem of solving the remaining $`3n+1`$ equations. Typically these equations do not form a closed set but involve integral quantities which are not included in the definitions of the momenta. To close them, i.e. to equal the numbers of equations and of unknowns, one must either restrict the problem or impose some kind of approximation.

Exact closure of moment equations.- We have found only two relevant cases in which the closure is exact for the rest of the equations. Both simplified problems correspond to conservative, $`\sigma =0`$, spherically symmetric potentials, $`\omega _i=\omega (t)`$ where

$`{\displaystyle \frac{dR}{dt}}`$ $`=`$ $`B_r,`$ (21)   
$`{\displaystyle \frac{dB_r}{dt}}`$ $`=`$ $`4K2\omega ^2(t)R2D,`$ (22)   
$`{\displaystyle \frac{dK}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\omega ^2(t)B_r{\displaystyle \frac{dJ}{dt}},`$ (23)   
being $`D=[Gg\rho ]`$ and $`R=\sqrt{W}`$ the radial width. The first integrable case corresponds to $`n=2`$ and $`G=U\rho ^2`$. For it Eqs. (The moment method in general Nonlinear Schrödinger Equations) simplify to

$$\frac{d^2R}{dt^2}=\omega (t)R+\frac{M}{R^3}.$$
(24)   
Here $`M`$ a constant that depends only on the initial data and interaction strength $`U`$. This equation has been used to prove the existence of extended resonances in Ref. .

Another ample family of systems for which moment equations are closed and exact is those with a time-independent interaction strength, $`G/t=0`$ (the usual case), and a divergenceless velocity distribution (given by the phase gradient) $`\text{div}\left(\varphi \right)=\mathrm{}\varphi =0.`$ This condition imposes no restriction on the density distribution $`\rho `$ and is automatically satisfied by the well known vortex-line solutions, which in $`n`$ spatial dimensions read

$`\psi `$ $`=`$ $`\rho _B(x_1,\mathrm{},x_n;t)e^{i\varphi _B},`$ (26)   
$`\varphi _B`$ $`=`$ $`\stackrel{}{\alpha }\stackrel{}{X}+\mathrm{arctan}{\displaystyle \frac{x_kX_k}{x_lX_l}}.`$ (27)   
Here $`\stackrel{}{X}`$ are free, time dependent parameters, and $`\rho _0`$ is arbitrary. With those simple conditions we get an infinite number of constants of evolution named “supermomenta”, $`Q(F)F(\rho )`$, built up from differentiable functions of the modulus, $`F(\rho )`$, that satisfy the regularity conditions. In these cases one can prove that $`D`$ is a constant and that Eqs. (The moment method in general Nonlinear Schrödinger Equations) become equivalent to Eq. (24).

Uniform divergence approximation.- Intuition dictates that the zero-Laplacian phase condition is related to (a) the configuration of the cloud, i.e. $`\rho `$, does not change, (b) the soliton or the wavepacket is either stationary or at most suffers *displacements and rotations* and (c) all of the “supermomenta” depend solely on the norm, $`Q(F)=f(N)`$. General solutions have nonzero divergence of the velocity field, thus the zeroth order approximation, $`\mathrm{}\varphi =0`$, fails to describe the system. The next possibility is a first order approximation in which the Laplacian of the phase is uniform. As we will see, this extension now allows for *changes in the shape* of the cloud and introduces three new independent variables in the supermomenta $`Q(F)`$. Mathematically, the first order approximation the phase is

$$\varphi =\varphi _B(\stackrel{}{x},t)+\underset{j}{}\beta _j(t)x_j^2,$$
(28)   
where $`\varphi _B(\stackrel{}{x},t)`$ is any function satisfying $`\mathrm{}\varphi _B=0`$. This approximation will be called uniform divergence approximation on the phase in what follows. A limited version of this approximation that uses a linear function (which has zero divergence) in place of $`\varphi _B`$ was first applied to radially symmetric problems in Ref.  and to the study of resonances in general 3D forced NSE problems in Refs. . In that case, were $`\varphi (\stackrel{}{x},t)=\varphi _0(t)+\stackrel{}{\alpha }\stackrel{}{x}+_j\beta _j(t)x_j^2`$ it is evident that the approximation consists on a Taylor expansion of the phase or even better a polynomial fitting with time dependent parameters. Here we provide a general framework for the application of the technique to arbitrary NSE problems. In a certain sense this is a generalization of the usual time dependent variational method but now no assumption on the shape of the amplitude of the wave is needed and the phase has a free is approximated by a least-squares type fitting with time-dependent parameters and only very general restrictions on the phase.

Let us first take the case with $`\varphi _B=\stackrel{}{\alpha }\stackrel{}{x}`$. We will assume for simplicity that the nonlinearity is or can be approximated by a polynomial

$$G(\rho )=\underset{k}{}\alpha _k\rho ^k.$$
(29)   
The dissipative terms can be removed from the equations by rescaling the solution with $`\gamma =_0^t\sigma (t^{})𝑑t^{},`$ so that $`\stackrel{~}{\rho }=e^\gamma \rho `$. We will denote the momenta obtained using $`\stackrel{~}{\rho }`$ by a tilde, as in $`\stackrel{~}{W}_i`$. It is important to stress that Eqs. (The moment method in general Nonlinear Schrödinger Equations) do not close immediately and it is necessary to study the monomial supermomenta

$$\stackrel{~}{Q}^{(m)}=\stackrel{~}{\rho }^md^nx.$$
(30)   
Their evolution laws are

$$\frac{d\stackrel{~}{Q}^{(m)}}{dt}=(m1)\underset{i}{}\beta _i\stackrel{~}{Q}^{(m)}.$$
(31)   
Since the parameters in the phase can be expressed as

$$\beta _i=\frac{d}{dt}\mathrm{log}\sqrt{\stackrel{~}{W}_i},$$
(32)   
one obtains a closed form for each of the supermomenta

$$\stackrel{~}{Q}^{(m)}=C^{(m)}\left(\sqrt{\stackrel{~}{W}_1\mathrm{}\stackrel{~}{W}_d}\right)^{m+1},$$
(33)   
where the constants $`C^{(m)}`$ must be determined from initial data. We can use this property to estimate all the integrals in which $`G(\rho )`$ appears, as a function of powers of the mean square widths. Using this, and defining as before the natural widths, $`\stackrel{~}{R}_i=\sqrt{\stackrel{~}{W}_i}`$, in the general nonsymmetric case we arrive to

$$\frac{d^2\stackrel{~}{R}_i}{dt}=\omega _i\stackrel{~}{R}_i+\frac{M_i}{\stackrel{~}{R}_i^3}+\frac{𝒢(\stackrel{~}{R}_1\mathrm{}\stackrel{~}{R}_d,t)}{\stackrel{~}{R}_1}$$
(34)   
where $`M_i`$ is a constant to be determined from the initial data, and $`𝒢(R_1\mathrm{}R_d,t)`$ is a function that may be calculated from Eqs. (29) and (33). Finally, to interpret these results one must remember that the actual width of the cloud is actually $`R_i(t)=\stackrel{~}{R}_i(t)e^{2\gamma }`$.

Eq. (34) means that given an initial data it is possible to compute the evolution of the width (and all the momenta of the initial datum) provided the parabolic phase is a good description for the solution. The method is very powerful as it accounts for the evolution of any solution that can be described in terms of a finite number of momenta (since higher order momenta are functions of lower order ones). Though much more general than the collective coordinate method the one presented here is simpler to apply and extend since it only involves computing the integrals in Eqs. (The moment method in general Nonlinear Schrödinger Equations).

We must also remark that the phase is not restricted to be a polynomial in $`\stackrel{}{r}`$ as in Ref. . Instead one still obtains information when the solution has many other forms. For example, if phase is of type (The moment method in general Nonlinear Schrödinger Equations), then one can combine the two widths from the plane on which the vorticity is present, $`W_k`$ and $`W_l`$, into a radial one, $`RW_k+W_l`$ and the evolution of $`R`$ can be studied together with the widths from the remaining spatial directions.

Independent moment approximation.- We have seen that the uniform divergence condition leads to a closed expression for every $`Q^{(m)}`$ in terms of the widths $`W_i`$. Although this approach is more powerful than the usual time dependent ansatz it is possible to improve the moment method to higher precision. The basic idea is to assume that only a finite set of momenta are independent, the number of those independent momenta being related to the accuracy of the solution and finding expressions for higher order momenta as functions of lower order ones. For instance if only the $`\{N,W_i\}`$ momenta are truly independent (as it occurs with the uniform divergence approximation) it is possible to find expressions for the rest of the magnitudes in terms of these momenta. By scaling the wave with respect to one of the coordinates

$$\psi (x_1,\mathrm{},x_i,\mathrm{})\frac{1}{\sqrt{1+ϵ}}\psi (x_1,\mathrm{},\frac{x_i}{1+ϵ},\mathrm{}),$$
(35)   
and relating the first order changes in $`Q^{(m)}`$ and $`W_k`$ one arrives to

$$(m1)Q^{(m)}=\underset{k}{}2\frac{Q^{(m)}}{W_k}W_k.$$
(36)   
This equation has a solution of the form reflected in Eq. (33) and also similar expressions may be derived for the rest of the unknown integrals in Eq. (3). Thus, in the first order case discussed here the method can provide a closed set of equations for any type of nonlinear term. In principle the procedure could be extended to higher precision approximation of the evolution of initial data. Work on this point is in progress and will be reported in detail elsewhere.

Finally we would like to point out that Eqs. (The moment method in general Nonlinear Schrödinger Equations) can be used straightforwardly by replacing $`\psi (\stackrel{}{r})`$ with an appropriate ansatz. In the conservative case, $`\sigma =0`$, this procedure is equivalent to Ritz’s optimization procedure, with the advantages that one needs not build a huge Lagrangian integral, and that our central equations (3) can also be applied to dissipative systems that lack a Lagrangian density at all. It must also be remarked that unless the ansatz has a phase different from (28), one will always arrive to Eq. (34).

In conclusion we have presented the moment method in a general framework and discussed under which conditions it leads to closed equations for a finite set of parameters. The uniform divergence ansatz has been introduced as a way to improve usual collective coordinate methods and obtain a closed set of equations for any NSE with polynomial nonlinearity and linear dissipation (with or without a parabolic time dependent spatial potential). The second method takes the zeroth, first and second order momenta to be the only independent ones and builds analytical approximations for other magnitudes. It is an extension of the uniform divergence approximation that seems to work for non polynomic self-interaction energies.

This work has been partially supported by CICYT under grant PB96-0534.