\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2010(2010), No. 89, pp. 1--12.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2010 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2010/89\hfil Numerical computation of soliton dynamics] {Numerical computation of soliton dynamics \\ for NLS equations in a driving potential} \author[M. Caliari, M. Squassina\hfil EJDE-2010/89\hfilneg] {Marco Caliari, Marco Squassina} % in alphabetical order \address{Marco Caliari \newline Dipartimento di Informatica \\ Universit\`a degli Studi di Verona\\ C\'a Vignal 2, Strada Le Grazie 15, I-37134 Verona, Italy} \email{marco.caliari@univr.it} \address{Marco Squassina \newline Dipartimento di Informatica \\ Universit\`a degli Studi di Verona\\ C\'a Vignal 2, Strada Le Grazie 15, I-37134 Verona, Italy} \email{marco.squassina@univr.it} \thanks{Submitted March 5, 2010. Published June 25, 2010.} \subjclass[2000]{35Q40, 58E30, 81Q05, 81Q20, 37N30} \keywords{Nonlinear Schr\"odinger equations; ground states; \hfill\break\indent soliton dynamics in an external potential; numerical computation of ground states; \hfill\break\indent semi-classical limit} \begin{abstract} We provide numerical computations for the soliton dynamics of the nonlinear Schr\"odinger equation with an external potential. After computing the ground state solution $r$ of a related elliptic equation we show that, in the semi-classical regime, the center of mass of the solution with initial datum built upon $r$ is driven by the solution to $\ddot x=-\nabla V(x)$. Finally, we provide examples and analyze the numerical errors in the two dimensional case when $V$ is a harmonic potential. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{ana}{Analytical Property}[section] \newcommand{\abs}[1]{\lvert#1\rvert} \newcommand{\norm}[2][]{\|#2\|_{#1}} \section{Introduction} \subsection{Soliton dynamics behaviour} The main goal of the present paper is to provide a numerical investigation of the so-called {\em soliton dynamics} for the nonlinear Schr\"odinger equation with an external time independent smooth potential $V$. That is, we study the qualitative behaviour of \begin{equation} \label{probMF} \begin{gathered} \mathrm{i}\varepsilon\partial_t\phi_\varepsilon=-\frac{\varepsilon^2}{2}\Delta \phi_\varepsilon +V(x)\phi_\varepsilon-|\phi_\varepsilon|^{2p}\phi_\varepsilon,\quad x\in\mathbb{R}^N,\, t>0, \\ \phi_\varepsilon(x,0)=\phi_0(x),\quad x\in\mathbb{R}^N, \end{gathered} \end{equation} in the semi-classical regime. Namely, for $\varepsilon$ (which of course plays the r\^ole of Planck's constant) going to zero, we take as initial datum a (bump like) function of the form \begin{equation} \label{initialD} \phi_0(x)=r\Big(\frac{x-x_0}{\varepsilon}\Big) e^{\frac{\mathrm{i}}{\varepsilon}x\cdot\xi_0},\quad x\in\mathbb{R}^N. \end{equation} We shall assume that $N\geq 1$, $00$. Finally, $x_0$ and $\xi_0$ are given vectors in $\mathbb{R}^N$ that should be conveniently thought (in the transition from quantum to classical mechanics) as corresponding to the {\em initial position} and {\em initial velocity} of a point particle. In this framework, since~\eqref{probMF} has a conservative nature, the typical expected behaviour is that the solution travels with the shape of $r((x-x(t))/\varepsilon)$ (hence its support shrinks, as $\varepsilon$ gets small) along a suitable concentration line $x(t)$ merely depending on the potential $V$ and starting at $x_0$ with initial slope $\xi_0$. On the basis of the analytical results currently available in literature (see the discussion in Section \ref{theoryfacts}), we believe that providing some numerical study is useful to complete the overall picture of this phenomenon and furnish some practical machinery for the computation of the solutions of~\eqref{seMF} and, in turn, of~\eqref{probMF}-\eqref{initialD}. The authors are not aware of any other contribution in the literature on this issue. For the linear Schr\"odinger, some results can be found in \cite{jinyang}. \subsection{Facts from the theory} \label{theoryfacts} It is well-known that, given a positive real number $m$, the afore mentioned (ground state) solution $r$ of~\eqref{seMF} (where the value of $\lambda$ depends on $m$) can be obtained through the following variational characterization on the sphere of $L^2(\mathbb{R}^N)$ \begin{equation} \label{gsminchar} \mathcal{E}(r)=\inf\{\mathcal{E}(u): u\in H^1(\mathbb{R}^N),\,\,\|u\|_{L^2}^2=m\}, \end{equation} where $\mathcal{E}:H^1(\mathbb{R}^N)\to\mathbb{R}$ is the $C^2$ energy functional \begin{equation} \label{energyfunct} \mathcal{E}(u)=\frac{1}{2}\int_{\mathbb{R}^N}|\nabla u|^2 dx -\frac{1}{p+1}\int_{\mathbb{R}^N}|u|^{2p+2}dx. \end{equation} Furthermore, there exists a suitable choice of $m$ yielding $\lambda=1$ as eigenvalue in \eqref{seMF}. The restriction to the values of $p$ below $2/N$ is strictly related to the global well-posedness of \eqref{probMF} for any choice of initial data $\phi_0$ in $H^1$. If $p$ is larger than or equal to $2/N$, then the solution can blow-up in finite time (see e.g.\ the monograph by Cazenave~\cite{cazenave}). In particular, in the two dimensional case, $p$ will be picked in $(0,1)$. From the analytical side, it has been rigorously known since 2000 that the solution $\phi_\varepsilon(t)$ of~\eqref{probMF} remains close to the ground state $r$, in the sense stated here below, locally uniformly in time, as $\varepsilon$ is converges to zero. As we said, this dynamical behaviour is typically known as soliton dynamics (for a recent general survey on solitons and their stability, see the work of T.\ Tao~\cite{tao-solitons}). For the nonlinear equation \eqref{probMF}, rigorous results about the soliton dynamics were obtained in various papers by Bronski, Jerrard~\cite{bronski} and Keraani~\cite{Keerani1,Keerani2}. We also refer to \cite{squamagn} for a complete study of the problem with the additional presence of an external time independent magnetic vector potential $A:\mathbb{R}^N\to\mathbb{R}^N$, and to~\cite{mopelsqu} for a study of a system of two coupled nonlinear Schr\"odinger equations, a topic which is rapidly spreading in the last few years. The arguments are mainly based on the following ingredients: the energy convexity estimates proved by Weinstein~\cite{weinsteinMS,weinstein2} to get the so called modulational stability, the use of conservation laws (mass and energy) satisfied by the equation, and the associated Hamiltonian system in $\mathbb{R}^N$ built upon the guiding external potential $V$, that is the classical {\em Newton law} \begin{equation} \label{newtonsenza} \begin{gathered} \ddot x(t)=-\nabla V(x(t)), \\ x(0)=x_0, \quad \dot x(0)=\xi_0. \end{gathered} \end{equation} Under reasonable assumptions on $V$ (e.g.\ uniform boundedness of the second order partial derivatives), equation \eqref{newtonsenza} admits a unique global solution $(x(t),\xi(t))$ which satisfies the following conservation law $$ \mathcal{H}(t)=\frac{1}{2}|\xi(t)|^2+V(x(t)),\quad \mathcal{H}(t)=\mathcal{H}(0),\quad t\geq 0. $$ Let us now define a suitable scaling of the standard norm of $H^1(\mathbb{R}^N)$, $$ \|\phi\|_{\mathbb{H}_\varepsilon}^2=\varepsilon^{2-N}\|\nabla \phi\|_{L^2}^2+\varepsilon^{-N}\|\phi\|_{L^2}^2,\quad\varepsilon>0. $$ The precise statement of the soliton dynamics reads as follows. \begin{ana}[Cf.\ e.g.\ \cite{bronski,Keerani2}] \label{anal11} Let $\phi_\varepsilon(t)$ be the solution to problem \eqref{probMF} corresponding to the initial datum~\eqref{initialD}. Then there exists a family of shifts $\theta_{\varepsilon}:\mathbb{R}^+\to [0,2\pi)$ such that, as $\varepsilon$ tends to zero, $\phi_\varepsilon(x,t)$ is equal to the function \begin{equation} \label{mainconclusintro} \phi_\varepsilon^r(x,t)= r\Big(\frac{x-x(t)}{\varepsilon}\Big)e^{\frac{{\rm i}}{\varepsilon}\left[x\cdot \dot x(t)+\theta_{\varepsilon}(t)\right]}, \quad x\in\mathbb{R}^N,\,\,t>0, \end{equation} up to an error function $\omega_\varepsilon(x,t)$ such that $\|\omega_\varepsilon(t)\|_{\mathbb{H}_\varepsilon}\leq {\mathcal O}(\varepsilon)$, locally uniformly in the variable $t$. \end{ana} It is important to stress that, in the particular case of {\em standing wave solutions} of~\eqref{probMF}, namely special solutions of~\eqref{probMF} of the form $$ \phi_\varepsilon(x,t)=u_\varepsilon(x) e^{-\frac{\mathrm{i}}{\varepsilon}\theta t},\quad x\in\mathbb{R}^N,\, t\in\mathbb{R}^+,\quad(\theta\in\mathbb{R}), $$ where $u_\varepsilon$ is a real-valued function, there is an enormous literature regarding the semi-classical limit for the corresponding elliptic equation $$ -\frac{\varepsilon^2}{2}\Delta u_\varepsilon+V(x)u_\varepsilon =|u_\varepsilon|^{2p}u_\varepsilon,\quad x\in\mathbb{R}^N. $$ See the recent book~\cite{ambook} by Ambrosetti and Malchiodi and the references therein. To this regard notice that, if $\xi_0=0$ (null initial velocity) and $x_0$ is a critical point of the potential $V$, as equation~\eqref{newtonsenza} admits the trivial solution $x(t)=x_0$ and $\dot x(t)=0$ for all $t\in\mathbb{R}^+$, formula~\eqref{mainconclusintro} reduces to \begin{equation*} \phi_\varepsilon^r(x,t)=r\Big(\frac{x-x_0}{\varepsilon}\Big) e^{\frac{{\rm i}}{\varepsilon}\theta_{\varepsilon}(t)},\quad x\in\mathbb{R}^N,\,\,t>0, \end{equation*} so that the concentration of $\phi_\varepsilon(t)$ is {\em static} and takes place at $x_0$, instead occurring along a smooth concentration curve in $\mathbb{R}^N$. This is consistent with the literature for the standing wave solutions mentioned above. For other achievements about the full dynamics of~\eqref{probMF}, see also~\cite{gril1,gril2} (in the framework of orbital stability of standing waves) as well as~\cite{kaup,keener} (in the framework of non-integrable perturbation of integrable systems). Similar results were investigated in geometric optics by a different technique (WKB method), namely writing formally the solution as $u_\varepsilon=U_\varepsilon(x,t) e^{\mathrm{i}\theta(x,t)/\varepsilon}$, with $U_\varepsilon=U_0 +\varepsilon U_1+\varepsilon^2 U_2\cdots,$ where $\theta$ and $U_j$ are solutions, respectively, of a Hamilton-Jacobi type equation (the eikonal equation) and of a system of transport equations. In the presence of a constant external potential, the orbital stability issue for problem~\eqref{probMF} was investigated by Cazenave and Lions~\cite{cl}, and by Weinstein in~\cite{weinsteinMS,weinstein2}. Then, Soffer and Weinstein proved in~\cite{soffer1} the asymptotic stability of nonlinear ground states of~\eqref{probMF}. See also the following important contributions: Buslaev and Perelman~\cite{buslaev1}, Buslaev and Sulem~\cite{buslaev3}, Fr\"ohlich, Gustafson, Jonsson, Sigal, Tsai and Yau~\cite{frolich1,frohl3,aboufrosig}, Holmer and Zworski~\cite{holmer}, Soffer and Weinstein~\cite{soffer2,soffer3}, Tsai and Yau~\cite{tsai1}. Another interesting problem concerns the case where the initial datum is multibump (for simplicity two bumps), say, \begin{equation} \label{initialD2} \phi_0(x)=r_1\Big(\frac{x-x_0}{\varepsilon}\Big) e^{\frac{\mathrm{i}}{\varepsilon}x\cdot\xi_0}+ r_2\Big(\frac{x-y_0}{\varepsilon}\Big) e^{\frac{\mathrm{i}}{\varepsilon}x\cdot\eta_0},\quad x\in\mathbb{R}^N. \end{equation} where $r_i$ are solutions to the problem $$ \mathcal{E}(r_i)=\inf\{\mathcal{E}(u): u\in H^1(\mathbb{R}^N),\quad \|u\|_{L^2}^2=m_i\}, $$ for some fixed $m_i>0$, $i=1,2$ and $x_0,y_0,\xi_0,\eta_0$ are taken as initial data for \begin{equation} \label{newtonsenza22} \begin{cases} \ddot x(t)=-\nabla V(x(t)), \\ x(0)=x_0, \quad \dot x(0)=\xi_0. \end{cases} \quad \begin{cases} \ddot y(t)=-\nabla V(y(t)), \\ y(0)=y_0, \quad \dot y(0)=\eta_0. \end{cases} \end{equation} Then we state the following result. \begin{ana}[Cf.\ e.g.\ \cite{aboufrosig}] Let $\phi_\varepsilon(t)$ be the solution to~\eqref{probMF} corresponding to the initial datum~\eqref{initialD2}. Then there exist two families of shifts $\theta_{\varepsilon}^i:\mathbb{R}^+\to [0,2\pi)$ such that, as $\varepsilon$ goes to zero, $\phi_\varepsilon(x,t)$ is equal to the function \begin{equation} \label{mainconclusintro2} \phi_\varepsilon^r(x,t)= r_1\Big(\frac{x-x(t)}{\varepsilon}\Big)e^{\frac{{\rm i}}{\varepsilon}\left[x\cdot \dot x(t)+\theta_{\varepsilon}^1(t)\right]} + r_2\Big(\frac{x-y(t)}{\varepsilon}\Big)e^{\frac{{\rm i}}{\varepsilon}\left[x\cdot \dot y(t)+\theta_{\varepsilon}^2(t)\right]}, \end{equation} up to an error function $\omega_\varepsilon(x,t)$ depending both on $\varepsilon$ and on the initial relative velocity $v=|\xi_0-\eta_0|$ (the larger is $v$ the smaller is the error), locally uniformly in time. \end{ana} See figure \ref{fig:3} in the final section for a movie showing this behaviour. \section{Numerical computation of the soliton dynamics} In the numerical simulations included in the last section of the paper, we shall consider the two dimensional case. On the other hand, here we consider the general case. \subsection{Overview of the method} Our purpose is to solve the Schr\"odinger equation \begin{equation}\label{eq:schroedinger} \begin{gathered} \mathrm{i} \varepsilon \partial_t \phi_\varepsilon(x,t)= -\frac{\varepsilon^2}{2}\Delta\phi_\varepsilon(x,t)+ V(x)\phi_\varepsilon(x,t)-\abs{\phi_\varepsilon(x,t)}^{2p} \phi_\varepsilon(x,t),\quad x\in\mathbb{R}^N,\\ \phi_\varepsilon(x,0)=r_\varepsilon(x-x_0),\quad x\in\mathbb{R}^N, \end{gathered} \end{equation} where $r_\varepsilon(x)=u(x/\varepsilon)$, so that $\phi(x,t)=u(x)e^{-\mathrm{i} \lambda t}$ is the solution of \begin{equation} \label{eq:gs} \mathrm{i} \partial_t \phi(x,t)=-\frac{1}{2}\Delta \phi(x,t)- \abs{\phi(x,t)}^{2p}\phi(x,t) \end{equation} being $u$ real, positive and minimizing the energy~\eqref{energyfunct} under the constraint $\norm[L^2]{u}^2=m$. Instead of a direct minimization of the energy (see, e.g., \cite{BT03,CORT09}), here we consider the parabolic differential equation \begin{equation}\label{eq:parabolic} \begin{gathered} \partial_t r(x,t)=\frac{1}{2}\Delta r(x,t)+r^{2p+1}(x,t)+ \lambda(r(x,t))r(x,t),\quad x\in\mathbb{R}^N,\; t>0\\ r(x,0)=r_0(x),\quad \norm[L^2]{r_0}^2=m,\quad x\in\mathbb{R}^N \end{gathered} \end{equation} with vanishing boundary conditions, where the map $t\mapsto \lambda(r(\cdot,t))$ is defined by \begin{equation*} \lambda(r(x,t))=\frac{\frac{1}{2}\int_{\mathbb{R}^N}\abs{\nabla r(x,t)}^2\mathrm{d} x- \int_{\mathbb{R}^N}\abs{r(x,t)}^{2p+2}\mathrm{d} x}{\norm[L^2]{r}^2}. \end{equation*} This approach is known as \emph{continuous normalized gradient flow} (see, e.g., \cite{BD04}), i.e. the continuous version of the propagation of the Sch\"rodinger equation along imaginary time $-{\rm i} t$ and projection to the $L^2$ sphere of radius $\sqrt{m}$. In equation~\eqref{eq:parabolic}, projection is not necessary and the energy decreases: in fact, if we multiply equation~\eqref{eq:parabolic} by $r(x,t)$ and integrate over $\mathbb{R}^N$, we easily get \begin{equation*} \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t}\norm[L^2]{r(\cdot,t)}^2 =\int_{\mathbb{R}^N} r(x,t)\partial_t r(x,t)\mathrm{d} x = 0 \end{equation*} and if we multiply equation~\eqref{eq:parabolic} by $\partial_t r(x,t)$ and integrate over $\mathbb{R}^N$, we get \begin{align*} \frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t} \mathcal{E}(r(\cdot,t))&= -\int_{\mathbb{R}^N}\abs{\partial_t r(x,t)}^2\mathrm{d} x+ \lambda(r(x,t))\int_{\mathbb{R}^N} r(x,t)\partial_t r(x,t)\mathrm{d} x \\ &=-\int_{\mathbb{R}^N}\abs{\partial_t r(x,t)}^2\mathrm{d} x\le0 \end{align*} Hence, the steady-state solution $r_\infty(x)=r(x,t\to\infty)$ of \eqref{eq:parabolic} satisfies $\norm[L^2]{r_\infty}^2=m$ and has a minimal energy. In fact, notice that by the results in \cite{kwong}, for any $\lambda>0$ there exists a unique (up to translations) positive and radially symmetric solution $r=r_\lambda$ of \eqref{seMF}. In turn, given $\lambda_1,\lambda_2>0$, if $r_1,r_2:\mathbb{R}^N\to\mathbb{R}$ denote, respectively, the positive radial solutions of the equations \begin{equation*} -\frac{1}{2}\Delta r_1+\lambda_1 r_1=r_1^{2p+1},\quad -\frac{1}{2}\Delta r_2+\lambda_2 r_2=r_2^{2p+1}, \end{equation*} then it is readily verified that \begin{equation*} r_2(x)=\mu r_1(\gamma x),\quad \gamma=\Big(\frac{\lambda_2}{\lambda_1}\Big)^{1/2},\quad \mu=\Big(\frac{\lambda_2}{\lambda_1}\Big)^{1/(2p)}, \end{equation*} which tells us that that, up to a scaling, the solution corresponding to different values of $\lambda$ is unique. Notice now that, due to the choice of the bump like initial datum (Gaussian like, see~\eqref{initialparabolic}) in the iterations to compute $r_\infty$ (see the discussion below), it turns out that $\lambda_\infty$, defined as $\lambda(r_\infty)$, is negative and $r_\infty$ is positive, radially symmetric (see figure~\ref{fig:1}) and solves \begin{equation*} -\frac{1}{2}\Delta r_\infty+\hat\lambda_\infty r_\infty=r_\infty^{2p+1}, \end{equation*} where $\hat\lambda_\infty=-\lambda_\infty>0$. If $r_m$ denotes the ground state solution (with the corresponding positive eigenvalue denoted by $\lambda_m$), then we have \begin{equation} \label{realvscomp} r_m(x)=\mu r_\infty(\gamma x),\quad \gamma=\Big(\frac{\lambda_m}{\lambda_\infty}\Big)^{1/2},\quad \mu=\Big(\frac{\lambda_m}{\lambda_\infty}\Big)^{1/(2p)}. \end{equation} \begin{figure}[htb] \includegraphics[scale=0.6]{fig1} % gs2 \caption{The positive, radially symmetric and radially decreasing ground state solution $r_\infty$ of~\eqref{gsminchar} with $m=1$ and $p = 0.2$. Of course, in the computation of $r_\infty$, there is a spurious imaginary part of maximum value around $10^{-16}$, since the complex FFT algorithm is involved. The corresponding value of $\lambda_\infty$ is $\lambda_\infty = -0.37921$.} \label{fig:1} \end{figure} On the other hand, by construction, we have $$ m=\|r_m\|_{L^2}^2=\int_{\mathbb{R}^N} r_m^2(x)dx =\mu^2\gamma^{-N}\int_{\mathbb{R}^N} r_\infty^2(x)dx=m\mu^2\gamma^{-N}, $$ namely $\mu^2\gamma^{-N}=1$. Finally, by the definition of $\gamma$ and $\mu$ in~\eqref{realvscomp}, we get $\lambda_m=\hat \lambda_\infty$ and $\gamma=\mu=1$, yielding from~\eqref{realvscomp} the desired conclusion; that is, $$ r_\infty=r_m. $$ Moreover, $r_\infty(x) e^{-\mathrm{i} \lambda(r_\infty(x))t}$ is a solution of \eqref{eq:gs}. We will take $r_\varepsilon(x-x_0)=r_\infty((x-x_0)/\varepsilon)$ as our candidate initial condition for the time-dependent nonlinear Schr\"odinger equation~\eqref{eq:schroedinger}. From a numerical point of view, it is convenient to compute directly $r_\infty(x/\varepsilon)$ instead of $r_\infty(x)$ and to apply the change of variable $\Phi(X,t)=\sqrt[4]{\varepsilon^N}\phi_\varepsilon(x,t)$, $\sqrt{\varepsilon}X=x$, to the nonlinear Schr\"odinger equation~\eqref{eq:schroedinger}, and hence to equation~\eqref{eq:parabolic}. Altogether, we need to solve \begin{equation}\label{eq:parabolictosolve} \begin{gathered} \partial_t R(X,t)=\frac{\varepsilon}{2}\Delta R(X,t)+ \varepsilon^{-Np/2}R(X,t)^{2p+1}+ \Lambda(R(X,t))R(X,t),\quad X\in\mathbb{R}^N\\ R(X,0)=R_0(X),\quad \norm[L^2]{R_0}^2=m\varepsilon^N,\quad X\in\mathbb{R}^N \end{gathered} \end{equation} with \begin{equation*} \Lambda(R(X,t))= \frac{\frac{\varepsilon}{2}\int_{\mathbb{R}^N}\abs{\nabla R(X,t)}^2\mathrm{d} X- \varepsilon^{-Np/2}\int_{\mathbb{R}^N}\abs{R(X,t)}^{2p+2}\mathrm{d} X}{\norm[L^2]{R}^2} \end{equation*} where $R(X,t)=r(x/\varepsilon,t)\sqrt[4]{\varepsilon^N}$. Since it is not possible to numerically integrate the equation up to an infinite time, we will consider $R(X,\bar{t})$ the steady-state as soon as $\mathcal{E}(R(X,\bar{t}))$ is stabilized within a prescribed tolerance. The initial condition $R_0(X)$ can be arbitrarily chosen (in the class of bump like functions), but an initial solution with small energy will shorten the ``steady-state'' time $\bar{t}$. Among the family of the Gaussian functions parameterized by $\sigma$ \begin{equation} \label{initialparabolic} R^{\sigma}(X)=\sqrt{m}\sigma^{N/2}e^{-\abs{\sigma X/\sqrt{\varepsilon}}^2/2} \sqrt[4]{\big(\frac{\varepsilon}{\pi}\big)^N} \end{equation} with $\norm[L^2]{R^\sigma}^2=m\varepsilon^N$ it is possible to choose the one with minimal energy. In fact \begin{equation*} \mathcal{E}(R^\sigma)=\sigma^2\frac{\varepsilon}{2}\int_{\mathbb{R}^N}\abs{\nabla R^1(X)}^2\mathrm{d} X -\sigma^{Np}\frac{\varepsilon^{-Np/2}}{p+1}\int_{\mathbb{R}^N} \abs{R^1(X)}^{2p+2}\mathrm{d} X. \end{equation*} If we define \begin{equation*} A=\frac{\varepsilon}{2}\int_{\mathbb{R}^N} \abs{\nabla R^1(X)}^2\mathrm{d} X,\quad B=\frac{\varepsilon^{-Np/2}}{p+1}\int_{\mathbb{R}^N} \abs{R^1(X)}^{2p+2}\mathrm{d} X. \end{equation*} the minimum for $\mathcal{E}(R^\sigma)$ is attained for \begin{equation*} \sigma=\Big(\frac{BNp}{2A}\Big)^{\frac{1}{2-Np}}. \end{equation*} The quantities $A$ and $B$ can be analytically computed and give \begin{gather*} A=\begin{cases} \frac{m\varepsilon}{4}& N=1\\ \frac{m\varepsilon^2}{2}& N=2\\ \frac{3m\varepsilon^3}{4}& N=3, \end{cases}\\ B=\frac{m^{p+1}\varepsilon^N}{\pi^{Np/2}(p+1)^{1+N/2}} \end{gather*} \subsection{Numerical discretization} With the normalization introduced above, the nonlinear Schr\"odinger equation to solve is \begin{equation} \label{eq:schroedingertosolve} \begin{gathered} \mathrm{i}\partial_t\Phi(X,t)= -\frac{1}{2}\Delta\Phi(X,t)+ \frac{V(\sqrt{\varepsilon}X)}{\varepsilon}\Phi(X,t) -\frac{\abs{\Phi(X,t)}^{2p}}{\sqrt{\varepsilon^{2+Np}}}\Phi(X,t),\quad X\in\mathbb{R}^N\\ \Phi(X,0)=R(X-X_0,\bar{t}),\quad X\in\mathbb{R}^N \end{gathered} \end{equation} A well-established numerical method for the cubic Schr\"odinger equation (focusing or defocusing case) is the Strang splitting \cite{BJM03,BJaM03,CNT09}. It is based on a split of the full equation into two parts, in which the first is spectrally discretized in space and then exactly solved in time and the second has an analytical solution. We used the Strang splitting method as well. The first part is \begin{subequations} \begin{equation}\label{eq:schro1} \mathrm{i}\partial_t\Phi_1(X,t)= -\frac{1}{2}\Delta\Phi_1(X,t). \end{equation} Thus, the Fourier coefficients of $\Phi_1(X,t)$ restricted to a sufficiently large space domain satisfy a linear and diagonal system of ODEs, which can be exactly solved. The second part is \begin{equation}\label{eq:schro2} \mathrm{i}\partial_t\Phi_2(X,t)= \frac{V(\sqrt{\varepsilon}X)}{\varepsilon}\Phi_2(X,t) -\frac{\abs{\Phi_2(X,t)}^{2p}}{\sqrt{\varepsilon^{2+Np}}}\Phi_2(X,t). \end{equation} \end{subequations} It is easy to show that the quantity $\abs{\Phi_2(X,t)}^{2p}$ is constant in time for this equation. Hence it has an analytical solution. Given the approximated solution $\Phi_n(X)\approx\Phi(X,t_n)$ of equation~\eqref{eq:schroedingertosolve}, a single time step of the Strang splitting Fourier spectral method can be summarized as: \begin{enumerate} \item take $\Phi_n(X)$ as initial solution at time $t_n$ for~\eqref{eq:schro1} and solve for a time step $k/2$, obtaining $\Phi_1(X,t_n+k/2)$; \item take $\Phi_1(X,t_n+k/2)$ as initial solution at time $t_n$ for \eqref{eq:schro2} and solve for a time step $k$, obtaining $\Phi_2(X,t_n+k)$;\label{realspace} \item take $\Phi_2(X,t_n+k)$ as initial solution for~\eqref{eq:schro1} and solve for a time step $k/2$, obtaining $\Phi_{n+1}(X)$. \end{enumerate} The result $\Phi_{n+1}(X)$ is an approximation of $\Phi(X,t_n+k)$. Since the solutions of the first part and the second part are trivial to compute in the spectral space and in the real space, respectively, it is necessary to transform the solution from spectral space to real and from real space to spectral before and after step~\eqref{realspace} above, respectively. All the transformations can be carried out by the FFT algorithm. The method turns out to be spectrally accurate in space and of the second order in time. Therefore, we used the Fourier spectral decomposition for the solution of equation~\eqref{eq:parabolictosolve}, too. Together with the Galerkin method, it yields a nonlinear system of ODEs \begin{equation}\label{eq:ODEs} \begin{gathered} \hat R'(t)=\frac{\varepsilon}{2}D \hat R(t)+f(\hat R(t)),\quad t>0,\\ \hat R(0)=\hat R_0, \end{gathered} \end{equation} where $\hat R$ is the vector of Fourier coefficients, $D$ the diagonal matrix of the eigenvalues of the Laplace operator and $f$ the truncated Fourier expansion of the whole nonlinear part of equation~\eqref{eq:parabolictosolve}. For the solution of equation~\eqref{eq:ODEs} we used an exponential Runge--Kutta method of order two (see, e.g., \cite{HO05}), with the embedded exponential Euler method. Given the approximation $\hat R_n\approx \hat R(t_n)$, a single time step of the method is \begin{enumerate} \item set $A_{n+1}=k_{n+1}\frac{\varepsilon}{2}D$ and $R_{n1}=R_n$; \item compute $\hat R_{n2}=\exp(A_{n+1})\hat R_{n1}+ k_{n+1}\varphi_1(A_{n+1})f(\hat R_{n1})$ (exponential Euler method); \item compute $\hat R_{n+1}=\hat R_{n2}+k_{n+1}\varphi_2(A_{n+1})(-f(\hat R_{n1})+ f(\hat R_{n2}))$\label{error} \end{enumerate} where $\varphi_1(z)$ and $\varphi_2(z)$ are the analytic functions \begin{align*} \varphi_1(z)&=\frac{e^z-1}{z},\ z\ne 0, & \varphi_2(z)&=\frac{e^z-1-z}{z^2},\ z \ne 0,\\ \varphi_1(0)&=1, & \varphi_2(0)&=\frac{1}{2}. \end{align*} The result is an approximation of $\hat R(t_n+k_{n+1})$. Exponential integrators are explicit and do not suffer of time step restrictions. However, they require the computation of matrix functions. In our case, the matrices involved $A_{n+1}$ are diagonal and the computation of the matrix functions $\exp(A_{n+1})$, $\varphi_1(A_{n+1})$ and $\varphi_2(A_{n+1})$ is trivial. In order to compute the terms $f(\hat R_{n1})$ and $f(\hat R_{n2})$, it is necessary to recover the functions in the real space corresponding to the Fourier spectral coefficients $\hat R_{n1}$ and $\hat R_{n2}$, respectively, then to compute the nonlinear part of equation~\eqref{eq:parabolictosolve} and finally to compute its Fourier transform. All the transformations can be carried out by the FFT algorithm. The term $\hat R_{n+1}-\hat R_{n2}$ in step~(\ref{error}) above can be used as an error estimate for $R(t_{n+1})-R_{n+1}$ and then it is possible to derive a variable time step integrator. This is particularly useful for our aim of computing the steady-state of the equation: in fact, we expect that the as soon as the solution approaches the steady-state it is possible to enlarge the time step, thus reducing the computational cost. The method turns out to be spectrally accurate in space and of the second order in time. \begin{figure}[htb] \centering \href{run:fig2a.avi}{\includegraphics[scale=0.72]{fig2a}} % movie1bh \href{run:fig2b.avi}{\includegraphics[scale=0.72]{fig2b}} % movie1bhirr \caption{In both simulation movies we set $\varepsilon=0.01$, $p = 0.02$, $m = 1$, $(x_0,y_0)=(-3.0,-3.0)$, $v_0 = (0,0)$. In the left movie, we chose $\omega_1=1.4$ and $\omega_2=1$ (rational ratio). In the right movie, we chose $\omega_1=\sqrt{2}$ and $\omega_2=1$ (irrational ratio). Notice that, although the ratios $\omega_2/\omega_1$ are very close in the two examples, the soliton dynamics is ergodic in the right movie. Of course the figures refer to the (squared modulus of the) solution at the time $t=0$ and contain the concentration paths (admitting an analytic expression) that the soliton travels.} \label{fig:2} \end{figure} \begin{figure}[htb] \centering \href{run:fig3.avi}{\includegraphics[scale=0.8]{fig3}} %movie2bh \caption{In the simulation movie, we set $\varepsilon=0.01$, $p = 0.02$, $m = 1$, $(x_0^1,y_0^1)=(-3,-3)$, $(x_0^2,y_0^2)=(1,1)$, $v_0^1=(2,0)$, $v_0^2=(0,0)$, $\omega_1=1.1$ and $\omega_2=1$.} \label{fig:3} \end{figure} \begin{figure}[ht] \centering \includegraphics[scale=0.45]{fig4} % herr \caption{For the error analysis, we set $p=0.02$, $(x_0,y_0)=(-0.5,-0.5)$, $m=1$, $\omega=(2,1)$ and a final time $t=\pi$. With the change of variable we used, $\sqrt{\varepsilon}X=x$ and $U(X)=\sqrt{\varepsilon}u(x)$, we have $\lVert u\rVert_{L^2}^2=\lVert U\rVert_{L^2}^2$ and $\lVert \nabla_x u\rVert_{L^2}^2=\frac{1}{\varepsilon}\lVert \nabla_X U\rVert_{L^2}^2$. Hence, the numerical error is computed through formula (written for the 2D case) $\| u\|_{\mathbb{H}_\varepsilon}= \sqrt{\varepsilon^{-1}\| \nabla_X U\|_{L^2}^2+ \varepsilon^{-2}\| U\|_{L^2}^2}$. As predicted by the Analytical Property~\ref{anal11}, the error in $\| \cdot\|_{\mathbb{H}_\varepsilon}$ is below the order ${\mathcal O}(\varepsilon)$.} \label{fig:4} \end{figure} \section{Two dimensional examples and error analysis} In this section, in order to provide some examples, we consider the two dimensional setting and focus on the physically relevant case of harmonic potential \begin{equation*} V(x,y)=\omega_1^2 x^2+\omega_2^2 y^2,\quad \omega_1,\omega_2>0, \end{equation*} well-established in the theory of Bose-Einstein condensates. In the two movies starting with figure~\ref{fig:2} we show the dynamics of the solitary wave along two Lissajous curves, periodic in the left side and ergodic for the right side. In the movie starting from figure~\ref{fig:3} we report the soliton dynamics in the case of an initial datum exhibiting a double bump behaviour (with a sufficiently large distance between the centers) up to the collision time. It is important to stress that in these figures the paths have an analytic expression and are plotted before the dynamics starts. The movies will then show that the centers of mass of the solitons follow adherently these curves up to the final computation time. An analysis of the error (in the single bump case) arising when the modulus of the solution $|\phi_\varepsilon(x,t)|$ is replaced by the modulus of the expression in the representation formula~\eqref{mainconclusintro}, namely $r((x-x(t))/\varepsilon)$, is indicated in figure~\ref{fig:4}. As predicted by the analytical property~\ref{anal11}, the error in the $\| \cdot\|_{\mathbb{H}_\varepsilon}$ is below the order ${\mathcal O}(\varepsilon)$. \subsection*{Acknowledgements} The first author was supported by the GNCS ``Programma Giovani Ricercatori''. The second author was supported by the Italian PRIN Project 2007: Metodi Variazionali e Topologici nello Studio di Fenomeni non Lineari. \clearpage \begin{thebibliography}{00} \bibitem{aboufrosig} Abou Salem W.K., Froehlich J., Sigal I.M., Colliding solitons for the nonlinear Schr\"odinger equation, {\em Comm. Math. Phys.} {\bf 291} (2009), 151--176. \bibitem{ambook} Ambrosetti A., Malchiodi A., Perturbation methods and semilinear elliptic problems on $\mathbb{R}^n$, Progress in Mathematics {\bf 240}, Birkh\"auser Verlag, Basel, 2006, xii+183 pp. \bibitem{BD04} Bao W., Du Q., Computing the ground state solution of Bose--Einstein condensates by a normalized gradient flow, {\em SIAM J.\ Sci.\ Comput.} {\bf 25} (2004), 1674--1697. \bibitem{BJaM03} Bao W., Jaksch D., Markowich P.A., Numerical solution of the Gross--Pitaevskii Equation for Bose--Einstein condensation, {\em J.\ Comput.\ Phys.} {\bf 187} (2003), 318--342. \bibitem{BJM03} Bao W., Jin S., Markowich P.A., Numerical study of time-splitting spectral discretizations of nonlinear Schr\"odinger equations in the semi-classical regimes, {\em SIAM J.\ Sci.\ Comput.} {\bf 25} (2003), 27--64. \bibitem{BT03} Bao W., Tang W.J., Ground state solution of Bose--Einstein condensate by directly minimizing the energy functional, {\em J.\ Comput.\ Phys.} {\bf 187} (2003), 230--254. \bibitem{berlions1} Beresticki H., Lions P.L., Nonlinear scalar fields equation I. Existence of a ground state, {\em Arch.\ Ration.\ Mech.\ Anal.} {\bf 82} (1983), 313--346. \bibitem{bronski} Bronski J., Jerrard R., Soliton dynamics in a potential, {\em Math.\ Res.\ Letters} {\bf 7} (2000), 329--342. \bibitem{buslaev1} Buslaev V.S., Perelman G.S., On the stability of solitary waves for nonlinear Schr\"odinger equations. Nonlinear evolution equations, 75-98, Amer.\ Math.\ Soc.\ Transl.\ Ser.\ 2, {\bf 164}, AMS, 1995. \bibitem{buslaev3} Buslaev V.S., Sulem C., On asymptotic stability of solitary waves for nonlinear Schr\"odinger equations, {\em Ann.\ Inst.\ H.\ Poincar\'e Anal.\ Non Lin\'eaire} {\bf 20} (2003), 419--475. \bibitem{CNT09} Caliari M., Neuhauser Ch., Thalhammer M., High-order time-splitting Hermite and Fourier spectral methods for the Gross–-Pitaevskii equation, {\em J. Comput. Phys.} {\bf 228} (2009), 822--832. \bibitem{CORT09} Caliari M., Ostermann A., Rainer S., Thalhammer M., A minimisation approach for computing the ground state of Gross--Pitaevskii systems, {\em J. Comput. Phys.} {\bf 228} (2009), 349--360. \bibitem{c} Carles R., WKB analysis for nonlinear Schr\"odinger equations with potential, {\em Commun.\ Math.\ Phys.} {\bf 269} (2007), 195--221. \bibitem{cazenave} Cazenave T., An introduction to nonlin\-ear Schr\"o\-din\-ger equation, Text.\ Metod.\ Mat., {\bf 26} Univ.\ Fed.\ Rio de Janeiro, 1993. \bibitem{cl} Cazenave T., Lions P.L., Orbital stability of standing waves for some nonlinear Schr\"odinger equations, {\em Comm.\ Math.\ Phys.} {\bf 85} (1982), 549--561. \bibitem{fw} Floer A., Weinstein A., Nonspreading wave packets for the cubic Schr\"odinger equation with a bounded potential, {\em J.\ Funct.\ Anal.} {\bf 69} (1986), 397--408. \bibitem{frolich1} Fr\"ohlich J., Gustafson S., Jonsson B.L.G., Sigal I.M., Dynamics of solitary waves external potentials, {\em Comm.\ Math.\ Phys.} {\bf 250} (2004), 613--642. \bibitem{frohl3} Fr\"ohlich J.; Tsai T.-P., Yau H.-T., On the point-particle (Newtonian) limit of the non-linear Hartree equation, {\em Comm.\ Math.\ Phys.} {\bf 225} (2002), 223--274. \bibitem{gril1} Grillakis M., Shatah, J., Strauss, W., Stability theory of solitary waves in the presence of symmetry. I, {\em J.\ Funct.\ Anal.} {\bf 74} (1987), 160--197. \bibitem{gril2} Grillakis M., Shatah, J., Strauss, W., Stability theory of solitary waves in the presence of symmetry. II, {\em J.\ Funct.\ Anal.} {\bf 94} (1990), 308--348. \bibitem{HO05} Hochbruck M., Ostermann A., Explicit exponential Runge-Kutta methods for semilinear parabolic problems, {\em SIAM J. Numer. Anal.} \textbf{43} (2005), no.~2--4, 323--339. \bibitem{holmer} Holmer J., Zworski M., Soliton interaction with slowly varying potentials, {\em Int.\ Math.\ Res.\ Not.} (2008), no. 10, 36 pp. \bibitem{frolich2} Jonsson B.L.G., Fr\"ohlich J., Gustafson S., Sigal I.M., Long time motion of NLS solitary waves in a confining potential, {\em Annals Henri Poincare} {\bf 7} (2006), 621--660. \bibitem{kaup} Kaup D.J., Newell A.C., Solitons as particles and oscillators and in slowly changing media: a singular perturbation theory, {\em Proc.\ Roy.\ Soc.\ London A.} {\bf 361} (1978), 413--446. \bibitem{keener} Keener J.P., McLaughlin D.W., Solitons under perturbation, {\em Phys.\ Rev.\ A} {\bf 16} (1977), 777--790. \bibitem{Keerani1} Keraani S., Semiclassical limit of a class of Schr\"odinger equation with potential. {\em Comm.\ Partial Differential Equations} {\bf 27} (2002), 693--704. \bibitem{Keerani2} Keraani S., Semiclassical limit for nonlinear Schr\"odinger equation with potential. II {\em Asymptotic\ Anal.} {\bf 47} (2006), 171--186. \bibitem{kwong} Kwong M.K., Uniqueness of positive solutions of $\Delta u-u+u\sp p=0$ in $\mathbb{R}^n$, {\em Arch.\ Rational Mech.\ Anal.} {\bf 105} (1989), 243--266. \bibitem{jinyang} Jin S., Yang X., Computation of the semiclassical limit of the Schr\"odinger equation with phase shift by a level set method, {\em J.\ Sci.\ Comput.} {\bf 35} (2008), 144--169. \bibitem{pllcc} Lions P.L., The concentration-compactness principle in the calculus of variations. The locally compact case. Part II, {\em Annales Inst.\ H.\ Poincar\'e Anal.\ Nonlin.} {\bf 1} (1984), 223--283. \bibitem{mopelsqu} Montefusco E., Pellacci B., Squassina M., Soliton dynamics for CNLS systems with potentials, {\em Asymptotic\ Anal.} {\bf 66} (2010), 61--86. \bibitem{soffer1} Soffer A., Weinstein M.I., Multichannel nonlinear scattering for nonintegrable equations, {\em Comm.\ Math.\ Phys.} {\bf 133} (1990), 119--146. \bibitem{soffer2} Soffer A., Weinstein M.I., Multichannel nonlinear scattering for nonintegrable equations. II. The case of anisotropic potentials and data, {\em J.\ Differential Equations} {\bf 98} (1992), 376--390. \bibitem{soffer3} Soffer A., Weinstein M.I., Selection of the ground state for nonlinear Schr\"odinger equations, {\em Rev.\ Math.\ Phys.} {\bf 16} (2004), 977--1071. \bibitem{squamagn} Squassina M., Soliton dynamics for the nonlinear Schr\"odinger equation with magnetic field, {\em Manuscripta Math.} {\bf 130} (2009), 461--494. \bibitem{sulemsulem} Sulem C., Sulem P.L., The nonlinear Schr\"odinger equation. Self-focusing and wave collapse. Applied Mathematical Sciences, {\bf 139}. Springer-Verlag, New York, 1999, +350pp. \bibitem{tao-solitons} Tao T., Why are solitons stable? {\em Bull.\ Amer.\ Math.\ Soc.} {\bf 46} (2009) 1--33. \bibitem{tsai1} Tsai T.-P., Yau H.-T., Asymptotic dynamics of nonlinear Schr\"odinger equations: resonance-dominated and dispersion-dominated solutions, {\em Comm.\ Pure Appl.\ Math.} {\bf 55} (2002), 153--216. \bibitem{weinsteinMS} Weinstein M., Modulation stability of ground state of nonlinear Schr\"odinger equations, {\em SIAM J.\ Math.\ Anal.} {\bf 16} (1985), 472--491. \bibitem{weinstein2} Weinstein M., Lyapunov stability of ground states of nonlinear dispersive evolution equations, {\em Comm.\ Pure Appl.\ Math.} {\bf 39} (1986), 51--67. \end{thebibliography} \end{document}