Nonlinear Schr\"odinger Equations for Bose-Einstein Condensates

TL;DR

This paper analyzes the nonlinear Schrödinger equation, specifically the Gross-Pitaevskii equation, to understand Bose-Einstein condensates' behavior under various conditions through analytical and numerical methods.

Contribution

It provides an analytical study of the NLS with $L^2$ initial data and includes numerical simulations for 2D GPE with anisotropic potentials, advancing understanding of BEC dynamics.

Findings

Analytical results on wave propagation in BECs.

Numerical simulations demonstrating wave behavior under anisotropic potentials.

Insights into the effects of electromagnetic fields on BECs.

Abstract

The Gross-Pitaevskii equation, or more generally the nonlinear Schr\"odinger equation, models the Bose-Einstein condensates in a macroscopic gaseous superfluid wave-matter state in ultra-cold temperature. We provide analytical study of the NLS with $L^2$ initial data in order to understand propagation of the defocusing and focusing waves for the BEC mechanism in the presence of electromagnetic fields. Numerical simulations are performed for the two-dimensional GPE with anisotropic quadratic potentials.

Tables

Table 1. Table 1 . Spatial discretization error analysis ‖ ψ e ​ x ​ a ​ c ​ t − ψ a ​ p ​ p ​ r ​ o ​ x ‖ L 2 subscript norm subscript 𝜓 𝑒 𝑥 𝑎 𝑐 𝑡 subscript 𝜓 𝑎 𝑝 𝑝 𝑟 𝑜 𝑥 superscript 𝐿 2 \|\psi_{exact}-\psi_{approx}\|_{L^{2}} at t = 1 𝑡 1 t=1 on R = [ − 8 , 8 ] 2 𝑅 superscript 8 8 2 R=[-8,8]^{2} (Defocusing κ = 1 𝜅 1 \kappa=1 , ε = 1 𝜀 1 \varepsilon=1 , Δ ​ t = 0.00005 Δ 𝑡 0.00005 \Delta t=0.00005 , ψ 0 = g 1 subscript 𝜓 0 subscript 𝑔 1 \psi_{0}=g_{1} )

Potential $V$	$h=\frac{1}{4}$	$h=\frac{1}{8}$	$h=\frac{1}{16}$	$h=\frac{1}{32}$
0	2.5353e-05	1.2107e-11	3.4148e-12	1.3345e-11
$\frac{x^2+y^2}{2\epsilon }$	1.8215e-05	4.0089e-10	1.7532e-10	8.5637e-11
$-\frac{x^2+y^2}{2\epsilon }$	8.7515e-04	9.0129e-06	3.3705e-06	1.5221e-06
$\frac{x^2+10y^2}{2\epsilon }$	1.7456e-01	1.3289e-03	7.5563e-10	7.9029e-11
$\frac{x^2-10y^2}{2\epsilon }$	2.7557e+00	4.8414e+00	3.5978e+00	5.0977e-02

Table 2. Table 2 . Temporal discretization error analysis ‖ ψ e ​ x ​ a ​ c ​ t − ψ a ​ p ​ p ​ r ​ o ​ x ‖ L 2 subscript norm subscript 𝜓 𝑒 𝑥 𝑎 𝑐 𝑡 subscript 𝜓 𝑎 𝑝 𝑝 𝑟 𝑜 𝑥 superscript 𝐿 2 \|\psi_{exact}-\psi_{approx}\|_{L^{2}} at t = 1 𝑡 1 t=1 on R = [ − 8 , 8 ] 2 𝑅 superscript 8 8 2 R=[-8,8]^{2}   (Defocusing κ = 1 𝜅 1 \kappa=1 , ε = 1 𝜀 1 \varepsilon=1 , Δ ​ x = Δ ​ y = 1 64 Δ 𝑥 Δ 𝑦 1 64 \Delta x=\Delta y=\frac{1}{64} , ψ 0 = g 1 subscript 𝜓 0 subscript 𝑔 1 \psi_{0}=g_{1} )

Potential $V$	$\Delta t=0.01$	$\Delta t=0.005$	$\Delta t=0.0025$	$\Delta t=0.00125$	$\Delta t=0.000625$
0	2.5615e-03	6.3592e-04	1.5871e-04	3.9662e-05	9.9139e-06
$\frac{x^2+y^2}{2\epsilon }$	1.3647e-02	3.4068e-03	8.5140e-04	2.1283e-04	5.3203e-05
$-\frac{x^2+y^2}{2\epsilon }$	4.9640e-02	1.2426e-02	3.1075e-03	7.7695e-04	1.9425e-04
$\frac{x^2+10y^2}{2\epsilon }$	2.7675e-01	6.7647e-02	1.6747e-02	4.1805e-03	1.0447e-03
$\frac{x^2-10y^2}{2\epsilon }$	1.3843e+01	4.1819e+00	1.1328e+00	3.1327e-01	5.5413e-02

Equations

$$i{u}_t=-\frac{1}{2}(\nabla -iA{)}^{2}u+Vu+\mu |u{|}^{p-1}u(t,x)\in {\mathbb{R}}^{1+n}$$

$$u(0,x)=u_0,$$

$$i\hslash \frac{\partial u}{\partial t}=\frac{\delta H}{\delta \bar{u}}$$

$$i{\partial }_{t}u=-\frac{1}{2}{\nabla }^{2}u+\mu |u{|}^{p-1}u+\tilde{V}u-\Omega \cdot Lu,$$

$$\tilde{V}(x)=\frac{1}{2}|A(x){|}^2+V(x)+\frac{i}{2}\mathrm{div}A(x)$$

$$i\Omega \cdot (x\wedge \nabla )=iA(x)\cdot \nabla .$$

$$B_{\tau }=B\cdot \frac{\text{r}}{r}=\left(\begin{array}{ccc}0& 2{\omega }_3& -2{\omega }_2\\ -2{\omega }_3& 0& 2{\omega }_1\\ 2{\omega }_2& -2{\omega }_1& 0\end{array}\right)\cdot \frac{\mathbf{r}}{r}=-\frac{2}{r}A,$$

$$A(x)=\frac{1}{4\pi }\int \frac{x-y}{|x-y{|}^3}\wedge B(y)dy.$$

$${\partial }^{\alpha }{A}_j(x)=O(1),\forall |\alpha |\ge 1,$$

$${\partial }^{\alpha }V(x)=O(1),\forall |\alpha |\ge 2$$

$$|{\partial }^{\alpha }B(x)|\le c_{\alpha }⟨x{⟩}^{-1-\epsilon },$$

$$\frac{2}{q}+\frac{n}{r}=\frac{n}{2}.$$

$$\Vert u{\Vert }_{S^0(I)}:=\underset{(q,r)admissible}{\sup }\Vert u{\Vert }_{L^q{L}^r(I\times {\mathbb{R}}^n)}<\infty .$$

$$\Vert f{\Vert }_{N^0(I)}=\underset{f=f_1+f_2}{\min }\left(\Vert f_1{\Vert }_{L^1{L}^2(I\times {\mathbb{R}}^n)}+\Vert f_2{\Vert }_{L^2{L}^{\frac{2n}{n+2}}(I\times {\mathbb{R}}^n)}\right),$$

$$\Vert u{\Vert }_{S^0((-T,T))}\le cT.$$

$$\Psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N)=\prod _{i=1}^N\varphi ({\mathbf{r}}_i).$$

$${\int }_{{\mathbb{R}}^3}|\varphi (\mathbf{r}){|}^{2}d\mathbf{r}=1.$$

$$H_N=\sum _{i=1}^N\left[\frac{{\mathbf{p}}_i^2}{2m}+V({\mathbf{r}}_i)\right]+U_0\sum _{i<j}\delta ({\mathbf{r}}_i-{\mathbf{r}}_j),$$

$$E_N=N{\int }_{{\mathbb{R}}^3}\left[\frac{{\hslash }^2}{2m}|\nabla \varphi (\mathbf{r}){|}^2+V(\mathbf{r})|\varphi (\mathbf{r}){|}^2+\frac{(N-1)}{2}{U}_0|\varphi (\mathbf{r}){|}^4\right]d\mathbf{r},$$

$${\int }_{{\mathbb{R}}^6}{U}_0\delta ({\mathbf{r}}_i-{\mathbf{r}}_j)⟨\varphi ({\mathbf{r}}_i)|\delta ({\mathbf{r}}_i-{\mathbf{r}}_i^{\prime })\varphi ({\mathbf{r}}_i^{\prime })⟩⟨\varphi ({\mathbf{r}}_j)|\delta ({\mathbf{r}}_j-{\mathbf{r}}_j^{\prime })\varphi ({\mathbf{r}}_j^{\prime })⟩d{\mathbf{r}}_{i}d{\mathbf{r}}_j$$

$$=$$

$$\psi (\mathbf{r})=N^{1/2}\varphi (\mathbf{r})$$

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +V(\mathbf{r})\psi +U_0|\psi {|}^2\psi .$$

$$H_N=\sum _{\iota =1}^N\left(-\frac{\hslash }{2m}{\nabla }_{A,\iota }^2+V({\mathbf{r}}_{\iota })\right)+\mu \sum _{\iota <j}^{N}g({\mathbf{r}}_{\iota }-{\mathbf{r}}_j),$$

$$i\hslash \frac{\partial }{\partial t}\psi =-\frac{{\hslash }^2}{2m}{\Delta }_A\psi +V\psi +\mu (g\ast |\psi {|}^2)\psi .$$

$$\Vert u(t){\Vert }_2=\Vert u_0{\Vert }_2$$

$$E(t):=\int (\mathcal{L}u)\bar{u}dx+\frac{2\mu }{p+1}\int |u{|}^{p+1}dx$$

$$=$$

$$i{u}_t=\mathcal{L}u+F(t),u(0)=u_0\in L^2({\mathbb{R}}^n),$$

$$u=(i{\partial }_t-\mathcal{L}{)}^{-1}F$$

$$=$$

$$U(t)f(x)=(2\pi it{)}^{-n/2}\int e^{iS(t,x,y)}a(t,x,y)f(y)dy,$$

$$\Vert U(t)f{\Vert }_{L^q{L}^r(I\times {\mathbb{R}}^n)}\le c_q\Vert f{\Vert }_2$$

$$\Vert {\int }_0^{t}U(t-s)F(s,\cdot )ds{\Vert }_{L^q{L}^r(I\times {\mathbb{R}}^n)}\le c_{q,\tilde{q}}\Vert F{\Vert }_{L^{{\tilde{q}}^{\prime }}{L}^{{\tilde{r}}^{\prime }}(I\times {\mathbb{R}}^n)},$$

Full text

1. Introduction

Consider the nonlinear Schrödinger equation (NLS)

[TABLE]

where $1\leq p<\infty$, $\mu\in\mathbb{R}$, $V:\mathbb{R}^{n}\to\mathbb{R}$ induces the electric field $-\nabla V$, and $A=(A_{1},\dots,A_{n}):\mathbb{R}^{n}\to\mathbb{R}^{n}$ induces the magnetic field $B=\nabla\wedge A=(\partial_{j}A_{k}-\partial_{k}A_{j})_{n\times n}$. Denote $\Delta_{A}=\nabla_{A}^{2}:=\sum_{j=1}^{n}(\frac{\partial}{\partial x_{j}}-iA_{j})^{2}$. Then $\mathcal{L}:=-\frac{1}{2}\Delta_{A}+V$ is an essentially selfadjoint Schrödinger operator with an electromagnetic potential $(A,V)$ that is gauge-invariant. The nonlinear term $F(u)=\mu|u|^{p-1}u$ has the property $\Im(\bar{u}F(u))=0$.

The physical significance for NLS in a magnetic field is well-known in nonlinear optics and Bose-Einstein condensate (BEC), where the magnetic structure is involved in scattering, superfluid, quantized vortices as well as DNLS in plasma physics [40, 29, 46]. There have been produced BEC where Bosons, Femions or other quasi-particles are trapped with atomic lasers in order to observe the macroscopic coherent wave matter in ultra-cold temperature.

The Hamiltonian $H:=\int\frac{\hbar^{2}}{2m}|\nabla_{A}u|^{2}+\frac{2\mu}{p+1}|u|^{p+1}$ generates the nonlinear system in (1):

[TABLE]

with $\hbar=m=1$ ($\hbar$ being the Planck constant and $m$ the mass of a particle), where we note that the adjoint of the covariant gradient $\nabla_{A}$ is $-\nabla_{A}$. When $p=3$, we obtain the Gross-Pitaevskii equation (GPE), which is regarded as a Ginzburg-Landau model in string theory. In general, the operator $\nabla_{A}^{2}=(\nabla-iA)^{2}$ contains components of both the (trapping) angular momentum and (attractive/repulsive) potential that can affect the dispersion of NLS. The equations $\eqref{e:u_A.V-mu.p}=\eqref{e:hu-hamilton.variation}$ play the role of Newton’s law in classical mechanics [49].

In the state of superfluid, the gaseous BEC has the vortices phenomenon which arises from (in the focusing case) the bound states of the form $u=e^{i\gamma t}Q$, $Q(x)=e^{im\theta}R_{m}(r)$ being an excited state. Another situation where it appears is when we test or manipulate the BEC by a magnetic trap with rotation. The wave function for the condensate is the solution of the following NLS

[TABLE]

where the rotation term $\Omega\cdot L=-i\Omega\cdot(x\wedge\nabla)$, $\Omega=(\omega_{1},\dots,\omega_{n})\in\mathbb{R}^{n}$ and $L$ denotes the angular momentum operator [31, 16, 4]. Comparing (1) and (4) one finds

[TABLE]

If $n=3$, then $\mathrm{div}\,A=0$ with $A=(\omega_{2}x_{3}-\omega_{3}x_{2},\omega_{3}x_{1}-\omega_{1}x_{3},\omega_{1}x_{2}-\omega_{2}x_{1})$. As a $2$-form $B=\mathrm{curl}\;A$ is constant. A simple calculation shows that the tangential component of $B$ is

[TABLE]

where ${\bf r}=(x_{1},x_{2},x_{3})$ and $r=|{\bf r}|$. This tells that $B$ is a “trapping" field whenever $\Omega\neq 0$. Heuristically $B_{\tau}\neq 0$ indicates an obstruction to the dispersion [23]. In $\mathbb{R}^{3}$, if the Coulomb gauge $\mathrm{div}\,A=0$, then one can recover $A$ as a “weighted wedge" of $x$ and $B$

[TABLE]

Let $\tilde{V}(x)=\frac{1}{2}\sum_{j}\gamma_{j}^{2}x_{j}^{2}$ and $|\gamma|=(\sum_{j}\gamma_{j}^{2})^{1/2}$. When $|\Omega|\ll|\gamma|$, the rotation action is negligible, and the potential $\tilde{V}$ is more predominant so that $A\approx 0$, $V\approx\tilde{V}=\frac{1}{2}\sum_{j}\gamma_{j}^{2}x_{j}^{2}$. In this case we anticipate trapping. When $|\Omega|\gg|\gamma|$, the rotation is much stronger than $\tilde{V}$ so that the effect of $\tilde{V}\approx 0\Rightarrow$ $V\approx-|A|^{2}/2$. This suggests that the wave function of a rotating BEC may be subject to an anisotropic repulsive potential. In this case the dispersion might hold global in time so that the (focusing) nonlinearity turns to be “short range" resulting in scattering [37, 12, 4]. Geometrically, the $x^{2}$ potential affects the wave like the trapping condition, which is stable, on a spherical portion of a manifold, while the $-x^{2}$ potential affects the wave like the scattering (non-trapping) condition on a hyperbolic portion of a manifold, which can be unstable locally in time but stable global in time.

The main analytical interest of this paper is to study the $L^{2}$ solution of (1) under the following assumptions on $A$ and $V$ throughout this section.

Assumption 1**.**

*Let $A_{j}$ and $V$ be real-valued and belong to $C^{\infty}(\mathbb{R}^{n})$. Let $V$ be bounded from below.

Assume $A=(A_{j})_{j=1}^{n}$ is sublinear and $V$ subquadratic, namely,*

[TABLE]

as $|x|\to\infty$. In addition, assume there exists some $\varepsilon>0$ such that for all $|\alpha|\geq 1$

[TABLE]

where $B=(b_{jk})_{n\times n}$, $b_{jk}=\partial_{j}A_{k}-\partial_{k}A_{j}$.

Define the $\mathcal{L}$-Sobolev space $\mathscr{H}^{s,r}:=\{u:\nabla^{s}u\in L^{r},\langle x\rangle^{s}u\in L^{r}\}$, where $\langle x\rangle=(1+|x|^{2})^{1/2}$. When $r=2$, we will also use the abbreviation $\mathscr{H}^{1}=\mathscr{H}^{1,2}$. For $u_{0}$ in $\mathscr{H}^{1}$, local wellposedness of (1) was proven for $1\leq p<1+4/(n-2)$ e.g., in [9, 43, 41] based on the fundamental solution constructed in [54]. The $\mathscr{H}^{s}$ subcritical result was considered in [58] for $1\leq p<1+4/(n-2s)$. When $s=1$, the following are known: Let $u_{0}\in\mathscr{H}^{1}$, $r=p+1$ and $q=\frac{4p+4}{n(p-1)}$.

1. (1)

Let $1\leq p<1+4/(n-2)$. Then in the defocusing case $\mu>0$, (1) has an $\mathscr{H}^{1}$-bounded global solution in $C(\mathbb{R},\mathscr{H}^{1})\cap L^{q}_{loc}(\mathbb{R},\mathscr{H}^{1,r})$. In the focusing case $\mu<0$, if $1\leq p<1+4/d$, then (1) has an $\mathscr{H}^{1}$-bounded global solution in $C(\mathbb{R},\mathscr{H}^{1})\cap L^{q}_{loc}(\mathbb{R},\mathscr{H}^{1,r})$. 2. (2)

Let $p=1+4/(n-2)$, $n\geq 3$. If $\|u_{0}\|_{\mathscr{H}^{1}}<\varepsilon$ for some $\varepsilon=\varepsilon(n,|\mu|)$ sufficiently small, then (1) has a unique local solution in $C((-T,T),\mathscr{H}^{1})\cap L^{q}((-T,T),\mathscr{H}^{1,r})$ for some $T>0$.

In two and three dimensions similar results on the $\mathscr{H}^{1}$ subcritical problem for (4) have been obtained in [30, 31, 4]. The main theorem (Theorem 1) we state below is the global wellposedness of (1) for $L^{2}$ initial data by virtue of the maximal Strichartz norm. This strengthens Theorem 3.3 in [58].

Definition 1**.**

We call $(q,r)=(q,r,n)$ an admissible pair if $q,r\in[2,\infty]$ satisfy $(q,r,n)\neq(2,\infty,2)$ and

[TABLE]

Definition 2**.**

Let $I\subset\mathbb{R}$ be an interval. The Strichartz space $S^{0}(I):=S^{0}(I\times\mathbb{R}^{n})$ is a Banach space consisting of functions in $\cap_{(q,r)\;admissible}L^{q}L^{r}(I\times\mathbb{R}^{n})$ satisfying

[TABLE]

Define $N^{0}(I)$ to be the linear span of $\cup_{(q,r)\;admissible}L^{q^{\prime}}L^{r^{\prime}}(I\times\mathbb{R}^{n})$, where $q^{\prime}=q/(q-1)$ is the Hölder conjugate of $q$. If $n\geq 3$, the admissible pairs include the endpoint pair $(2,\frac{2n}{n-2})$, which allows us to identify $S^{0}(I)$ with $L^{\infty}L^{2}\cap L^{2}L^{\frac{2n}{n-2}}(I\times\mathbb{R}^{n})$ through interpolation. In this case $N^{0}(I)=L^{1}L^{2}+L^{2}L^{\frac{2n}{n+2}}(I\times\mathbb{R}^{n})$ and $N^{0}(I)$ is endowed with the norm

[TABLE]

where the infimum is taken over all possible $f_{1}\in L^{1}L^{2}(I\times\mathbb{R}^{n})$ and $f_{2}\in L^{2}L^{\frac{2n}{n+2}}(I\times\mathbb{R}^{n})$ such that $f=f_{1}+f_{2}$ [36].

Theorem 1**.**

Let $A$ and $V$ satisfy the conditions in Assumption 1. Suppose $u_{0}\in L^{2}(\mathbb{R}^{n})$.

1. (1)

If $1\leq p<1+4/n$, then equation (1) has a unique solution $u$ in $C(\mathbb{R},L^{2}(\mathbb{R}^{n}))\cap S^{0}_{loc}(\mathbb{R}\times\mathbb{R}^{n})$. Furthermore, for any $R>0$ there exists $T_{R}>0$ such that the flow $u_{0}\mapsto u$ is Lipschitz continuous from $\mathcal{B}_{R}$ into $S^{0}((-T_{R},T_{R}))$. 2. (2)

If $p=1+4/n$, then there exists an $\varepsilon>0$ such that $\|u_{0}\|_{2}<\varepsilon$ implies that equation (1) has a unique solution $u$ in $C(\mathbb{R},L^{2}(\mathbb{R}^{n}))\cap S^{0}_{loc}(\mathbb{R}\times\mathbb{R}^{n})$. The flow $u_{0}\mapsto u$ is Lipschitz continuous from $\mathscr{B}_{\epsilon/2}$ into $S^{0}((-T_{0},T_{0}))$.

In both cases, it holds that for all $T>0$,

[TABLE]

In the above, $\mathscr{B}_{R}:=\{u:\|u\|_{2}\leq R\}$, $T_{0}=T_{0}(A,V)$, and $\varepsilon$ and $c$ are constants depending on $n$, $\mu$ and $\|u_{0}\|_{2}$ only.

In the focusing ($\mu<0$), $L^{2}$ critical or supercritical but energy subcritical regime $1+4/n\leq p<1+4/(n-2)$, there can occur finite time blowup solutions for (1), see e.g., [9, 11, 49]. Such situation is more complicated, where the occurrence of wave collapse is equivalent to the existence of soliton, which depends on the interaction between linear and nonlinear energies, the expectation of momentum as well as the profile of the initial data. For the rotating problem (4), wave collapse can occur for either cases where $|\Omega|\ll|\gamma|$ or $|\Omega|\gg|\gamma|$. In [4] blowup conditions are given in terms of $(\Omega\cdot L)V$. More recently, Garcia [25] obtained a general blowup criteria for (1) based on spectral properties of $A$ and $V$.

It is desirable to observe numerical results that can experimentally verify the theory. In Section 4 we apply the Strang splitting scheme to find numerical solutions for the GPE (1) in 2D (a cubic NLS) where we take $A=0$ and $V(x_{1},x_{2})=\frac{1}{2}\sum_{j=1}^{2}\delta_{j}\gamma_{j}^{2}x_{j}^{2}$, $\delta_{j}\in\{\pm 1\}$. Our algorithm and implementations are based on time-splitting Fourier-spectral methods developed in [6] and GPELab [27]. Such scheme is stable and has higher accuracy under appropriate conditions on $V$ and initial data, see [39, 38]. Numerical schemes typically use spectral or pseudo-spectral method to approximate the solution by discretizing spacial dimensions and then advancing a time step, while physicists have used e.g., Crank-Nicholson method via Lagrangian variational techniques [47, 18].

The organization of the remaining of the paper is as follows. In Section 2, the time dependent Gross-Pitaevskii equation, in particular the electromagnetic GPE, is introduced and formally derived as mean field approximation for the $N$-particle state of the BEC. In Section 3, we prove Theorem 1, mainly in the $L^{2}$-critical case, concerning global wellposedness of (1). In Section 4, we present numeral simulations to illustrate the focusing and defocusing nonlinear effects on the wave function of BEC subject to various anisotropic harmonic potentials.

2. Formal derivation of the Gross-Pitaevskii equation

In the early stage of quantum mechanics there arose questions concerning fundamental aspects of decoherence and measurement theory as well as understanding the correlation between classical and quantum scattering models. In 1924, Satyendra Nath Bose published a paper describing the statistical nature of light [8]. Using Bose’s paper, Albert Einstein predicted a phase transition in a gas of noninteracting atoms could occur due to these quantum statistical effects [20, 21]. This phase transition period, Bose-Einsten Condensation, would allow for a macroscopic number of non-interacting bosons to simultaneously occupy the same quantum state of lowest energy.

It wasn’t until 1938, with the discovery of superfluidity in liquid helium, that F. London conjectured that this superfluidity may be one of the first manifestations of BEC. The real breakthrough came in 1995, when the BEC were produced from a vapor of rubidium, and of sodium atoms [2, 17].

The Gross-Pitaevskii equation (1), $p=3$ describes the macroscopic wave functions $u$ of the condensate in the presence of the magnetic and electric potentials $A$ and $V$. The nonlinear term results from the mean field interaction between atoms. The constant $\mu$ accounts for the attractive ($\mu<0$) or repulsive ($\mu>0$) interaction, whose sign depends on the chemical elements.

Nowadays BEC can be simulated in the computer and the lab. The rotating BEC, for instance, involves the decoherence $\leftrightarrow$ coherence phase. The angular momentum operator breaks up the beams, hence split the spectral lines when performing the experiment on silver atoms in normal state. It can help create quasi-particles so to manipulate or observe not only the macroscopic atoms, but also individual particle. There are potential applications in higher degree precision for measurement, navigation, computing and communications.

2.1. A formal derivation

We follow a mean-field approach to derive the time-dependent GPE for the $N$-body system of bosons. At ultra low temperatures, all bosons exist in identical single-particle state $\phi({\bf r})$, ${\bf r}\in\mathbb{R}^{3}$ and so we can write the wave function of the $N$-particle system as

[TABLE]

The single-particle wave function $\phi({\bf r})$ obeys the typical normalization condition

[TABLE]

Due to the fact that we are dealing with dilute gases, the distance between any two particles in positions ${\bf r}$ and ${\bf r}^{\prime}$ is such that the only interaction term is $U_{0}\delta({\bf r}-{\bf r}^{\prime})$, where $\delta$ is the usual Dirac function and $U_{0}=\frac{4\pi\hbar^{2}a}{m}$ is the strength of effective contact interaction ($a$ being the scattering length). Thus the Hamiltonian reads

[TABLE]

where ${\bf p}=-i\hbar\nabla$ stands for the momentum and $V({\bf r})$ the external potential. Meanwhile the $N$-state (7) has energy

[TABLE]

where the nonlinear energy term is attributed to the inherent self-interaction and interaction between a pair of bosons on the same state

[TABLE]

This is equivalent to an expression in terms of the expectation of the collision contact.

Introduce the wave function for the condensed state

[TABLE]

so that $N=\int|\psi|^{2}d{\bf r}$. By a variation argument for $E_{N}$, similar to (3) we formulate the GPE as $N\to\infty$

[TABLE]

Derivation of magnetic NLS.

In a similar way we can formally derive (1) for $p=3$. Let $A\in L_{loc}^{2}(\mathbb{R}^{3},\mathbb{R}^{3})$, $V:\mathbb{R}^{3}\to\mathbb{R}$. Assume an $N$-particle weakly interacting condensate of non-relativistic bosons without spin in the mean field. The Hamiltonian in the electromagnetic frame has the form on $\mathbb{R}^{3N}$

[TABLE]

where $\nabla_{A}=\nabla-iA$ is the covariant gradient on $\mathbb{R}^{3}$, $V$ represents the external potential, $\mu g$ the inherent potential for a two-body bosons, that is, the interaction between two particles is given by $\mu g({\bf r}-{\bf r}^{\prime})$. Using the fact that the expectation at $(t,{\bf r})$ of the interaction from the $\iota$-th particle is $\mu\int_{\mathbb{R}^{3}}g({\bf r}-{\bf r}_{\iota})|\psi(t,{\bf r}_{\iota})|^{2}d{\bf r}_{\iota}$ we arrive at the GPE that decries the wave function of the condensate

[TABLE]

In the case $g=\delta$ where only local contact interaction from collision is accounted for while other interactions are neglected in a dilute gas, the equation becomes the standard magnetic cubic NLS. ∎

Remark.  The derivation above relies on the fact that the $N$ particles of a dilute gas are condensed in the same state for which the wave function minimized the energy. The note [28] contains derivation and discussions of the magnetic GPE in the physical setting. For rigorous derivation of the mean field limit of the $N$-particle coherent state as $N\to\infty$ as well as $t\to\infty$ involving ground state trapping and scattering (dispersion) we refer to [24, 22].

GPE with harmonic potential and angular momentum

In (9), $|\psi(t,x)|^{2}$ denotes the probability density of the condensate at $(t,x)$. The coefficient $\mu$ measures the strength of interaction and depends on a quantity called the $s$-scattering length. It has positive sign (defocusing) for $\tensor*[^{87}]{Rb}{}$, $\tensor*[^{23}]{Na}{}$, $\tensor*[^{1}]{H}{}$ atoms, but negative sign (focusing) for $\tensor*[^{7}]{Li}{}$, $\tensor*[^{85}]{Rb}{}$, $\tensor*[^{133}]{Cs}{}$ [53, 11]. The typical example $V=\frac{1}{2}\sum_{j}\gamma_{j}^{2}x_{j}^{2}$ represents an external trapping potential imposed by a system of laser beams, where $\gamma_{1},\gamma_{2},\gamma_{3}$ are the magnitudes of the frequencies of the oscillator in three directions. It works as an anisotropic trap that allows one to observe the behavior of macroscopic waves traveling along a waveguide with varying width or excitations when a BEC is released from a trap.

With the addition of a rotation term we arrive at the GPE in (4). This equation is viewed as a conservation for the angular momentum on a quantum level that involves Newton’s law and Lorentz force where the magnetic field is divergence free. The momentum operator $L_{\Omega}:=i\Omega\cdot(x\wedge\nabla)$ with non-vanishing angular velocity $\Omega$ gives rise to vortex lattices in a condensate that supports it in turn, e.g., one can obtain the vortex lattices of a BEC by setting the $Na$ condensate in rotation using laser beams [3, 42].

The study of BEC as a rotating superfluid leads to the quantization of circulation and quantized vortices. Physically this makes it impossible for a superfluid to rotate as a rigid body: In order to rotate, it must swirl [3]. The existence of quantized vortices with such particular pattern has been verified by experiments and numerics, see e.g., [40, 7] and [1, 16]. They can be observed in a condensate with either optical or magnetic traps.

3. The $L^{2}$ solution using maximal Strichartz norm

This section is devoted to the proof of Theorem 1. We let $A$ and $V$ satisfy Assumption 1. A priori, note that equation (1) has the conservation of mass and energy on its lifespan

[TABLE]

Let $u$ and $F$ be $L^{2}\cap L^{r}(\mathbb{R}^{n})$-valued functions in $t\in I$, $I$ an interval in $\mathbb{R}$. If $u$ solves

[TABLE]

then the solution can be expressed in an integral form according to Duhamel principle

[TABLE]

From [54] we know there exists $T_{0}$ such that for $0<|t|<T_{0}$ the propagator $U(t):=e^{-it\mathcal{L}}$ is given as

[TABLE]

where $S(t,x,y)$ is a real solution of the Hamilton-Jacobi equation, both $S$ and $a$ are $C^{1}$ in $(t,x,y)$ and $C^{\infty}$ in $(x,y)$, with $|\partial_{x}^{\alpha}\partial_{y}^{\beta}a(t,x,y)|\leq c_{\alpha\beta}$ for all $\alpha,\beta$. Write $I:=I_{T_{0}}=[-T_{0},T_{0}]$ and $L^{q}L^{r}(I\times\mathbb{R}^{n})=L_{t}^{q}(I,L_{x}^{r}(\mathbb{R}^{n}))$.

Lemma 1** (Strichartz estimates [9, 58]).**

If $A$ and $V$ satisfy Assumption 1, then we have for $I=[-T_{0},T_{0}]$, there exist constants $c_{q},c_{q,\tilde{q}}$ such that

• [TABLE]
• [TABLE]

where $(q,r),(\tilde{q},\tilde{r})$ are any admissible pairs, and $q^{\prime}$ is the Hölder conjugate of $q$.

The Strichartz estimates yield the following lemma, consult [51, Chapter 3].

Lemma 2**.**

Let $u$ be a solution of (12). Then for any admissible pairs $(q,r)$, $(\tilde{q},\tilde{r})$ as in Definition 1 we have

[TABLE]

Moreover,

[TABLE]

Now we begin to prove part (2) in Theorem 1.

Proof of (2) in Theorem 1.

(I) Let $p=1+4/n$. According to (18), we have

[TABLE]

Since $N^{0}(I)\supset\cup_{(q,r)\;admissible}L^{q}L^{r}(I\times\mathbb{R}^{n})$ and $q=r={(2n+4)}/{n}$ are admissible, it follows that

[TABLE]

Hence we obtain

[TABLE]

(II) Let $\|u_{0}\|_{2}\leq\varepsilon:=\eta\gamma=(2c_{n})^{-1-n/4}\min(1,(5|\mu|)^{-n/4})$, where we choose $\eta=(2c_{n})^{-1}$ and $\gamma=\min((2c_{n})^{-n/4},(10c_{n}|\mu|)^{-n/4})$. In view of (13) we need to prove that the mapping

[TABLE]

is a contraction on the closed set $E_{\gamma}=\{u\in S^{0}(I):\|u\|_{S^{0}(I)}\leq\gamma\}$.

(a) In doing so first we show $\Phi$: $E_{\gamma}\to E_{\gamma}$. According to (19) we have, for $u\in E_{\gamma}$

[TABLE]

(b) Then we show that $\Phi$ is contraction on $E_{\gamma}$. Note the following inequality: For all $p>1$

[TABLE]

Hölder inequality gives, with $p=1+4/n$,

[TABLE]

The same type of inequality holds with $|u|^{4/n}(u-v)$ replaced by $|v|^{4/n}(u-v)$.

Hence, applying Lemma 2 we obtain, for $p=1+4/n$ and ${\tilde{q}}={\tilde{r}}=(2n+4)/n$

[TABLE]

It follows that for $u,v\in E_{\gamma}$

[TABLE]

by the choice of $\gamma$ above, where $C_{n}=2c_{n}p|\mu|$ . Therefore we have proved that $\Phi$ has a fixed point in the set $E_{\gamma}$. We conclude that if $\|u_{0}\|_{2}\leq\varepsilon$, there exists a unique solution $u$ in $L^{\infty}([-T_{0},T_{0}],L^{2}(\mathbb{R}^{n}))\cap S^{0}([-T_{0},T_{0}],\mathbb{R}^{n})$. The global in time existence follows from the conservation of the $L^{2}$ norm (10) by observing that $\varepsilon$ only depends on $n$ and $\mu$.

(III) The Lipschitz continuity is based on iteration of the contraction $\Phi$, see e.g., Proposition 1.38 and Proposition 3.17 in [51]. ∎

Proof of (1) in Theorem 1.

Let $p<1+4/n$. The proof for the subcritical case follow the same line as for the critical case but use the following: Choose $(q,r)=(\tilde{q},\tilde{r})=(\frac{4p+4}{n(p-1)},p+1)$ to arrive at

[TABLE]

if choosing $T=T(\|u_{0}\|_{2})>0$ sufficiently small. Here we notice that $\alpha=\frac{4-n(p-1)}{4}>0\iff p<1+4/n$. ∎

Remarks:  For $A=V=0$ the analogous result was proven in [52, 15, 13] using $L_{t}^{q}L_{x}^{r}$ norm. The case where $A=0$ and $V$ is subquadratic or quadratic was treated in [44, 11, 12]. The proof presented here is a modification of the standard argument, see [51].

When $1+4/n\leq p<1+4/(n-2)$, Carles [12, Theorem 1.4] shows that if $A=0$ and $V=-\frac{1}{2}|x|^{2}$ (more generally, $V$ has a stronger repulsive component), then global in time existence and scattering hold in $\mathscr{H}^{1}$. Carles’ proof relies on global in time Strichartz estimate where the repulsive component of $V$ produces exponential decay for $U(t)$ that balances the confining force from its attractive component to control the nonlinear effects. In the energy critical case $p=1+4/(n-2)$, Killip, Visan and Zhang proved the GWP and scattering for radial initial data in $\mathscr{H}^{1}$ [55, 34].

4. Numerical simulations for GPE

The Strang splitting method [50] deals with hyperbolic model problems with second order accuracy for finite difference schemes, which initially appeared in [35]. For a general nonlinear system one can write

[TABLE]

to obtain the following two equations for which $u=v+w$ and

[TABLE]

In the NLS case this method also has second order stability [39]. Since the splitting scheme can preserve the structure of the PDE, it also preserves the same conservation quantities (10) and (11) for the numerical solution as well as the analytic solution.

In this section we apply the Strang splitting algorithm to find numerical solution of the GPE in two dimensions. By truncation we consider the following equation defined on the rectangle $R:=[a,b]\times[c,d]$ with periodic boundary conditions

[TABLE]

The algorithm is implemented based on time-splitting trigonometric spectral approximations with fine mesh grids and time steps. The solutions are computed mainly using GPELab [27] adapted to various cases where $V(x,y)={(\pm\gamma_{1}^{2}x^{2}\pm\gamma^{2}_{2}y^{2})}/2$, $\kappa>0$ or $\kappa<0$. The initial data is taken as either a gaussian in $C^{\infty}(\mathbb{R}^{2})$ or a hat function in $H^{1}(\mathbb{R}^{2})$.

We summarize the numerical results in Figures 1 to 7 and then provide error analysis in Tables 1 and 2 with progressively finer and finer mesh sizes and time steps. These errors are relatively very small and yield quite high accuracy. Relevant error estimates can be found in [39] and equation (26) in [38]. Corresponding to two type of nonlinear regimes, we will select $\kappa=1$ and $\kappa=-1.9718$ for repulsive and attractive interactions, respectively. Let $a=c=-8,b=d=8$. Let $h=\Delta x=\Delta y=(b-a)/M$ be the meshgrid size and $\Delta t=T/N$ the time step.

A. Defocusing case: $\kappa=1>0$.  Set the initial data $\psi_{0}(x,y)=g_{\sigma}(x,y)$ with $\sigma=1$, where

[TABLE]

The following show the figures for the numerical solution $\psi_{approx}$ of (21) on $R=[-8,8]^{2}$ at different times in the presence of anisotropic quadratic potentials. The numerical results are in consistence with the theory that attractive $V$ confines the waves to the ground state while the repulsive $V$ enhances the dispersion or scattering.

B. Focusing case: $\kappa=1.9718<0$ with the same gaussian initial data (22).  In the mass-critical case $p=1+4/n=3$, the focusing NLS may have finite blowup solution. The physics dictates that a positive harmonic potential is attractive and confines the cooled bosonic atoms. On the other hand, an inverted (negative) harmonic potential is repulsive and supports the dispersion which offsets the impact of the focusing effect.

Now we observe from Figure 4 that if $V$ changes from $(x^{2}+y^{2})/{2\varepsilon}$ to $-(x^{2}+y^{2})/{2\varepsilon}$, then the blowup time has a slight delay at approximately $t=0.35$.

C. Initial data equal to the “hat function” $\psi_{0}=h\in H^{1}(\mathbb{R}^{2})$.  In Figures 5(a) and 5(b) we compare the density function $|\psi|^{2}$ for two different initial data, one is given by the gaussian (22) and the other is given by the “hat function”

[TABLE]

The solutions show that if the magnitude of the frequency is large enough then the inverted harmonic potential $V$ counteracts the focusing effect which leads to global in time existence. On a quite long time interval, they both reveal self-similarity (“multifractal-like”) for the density function although $h\in H^{1}$ has a larger magnitude with low regularity. However with $\psi_{0}$ being the hat function, $|\psi(t,x,y)|^{2}$ is more irregular in temporal and much more singular in spatial variables, see Figure 6.

The numerical results agree with Theorem 1 and [58, Theorem 3.2]. Note that on local time interval the amplitude of $\psi$ is higher than in the free case $V=0$. The lack of the long time decay or scattering may be due to the fact that equation (21) has a “truncation” version.

Tables 1 and 2 show the error analysis between the approximate solution and the exact solution. In both cases the results have reached good accuracy as well as efficiency. However when we test on the case where $\kappa=-1.9718$ ($\varepsilon=0.3$, $\Delta t=0.00005$), the spatial and the temporal error analysis seem to indicate quite big difference between the use of the gaussian and the use of the hat function. Numerical result shows that in the focusing case, if $\psi_{0}=h$, then the approximation solution $\psi_{approx}$ along with the error $\|\psi_{exact}-\psi_{approx}\|_{L^{2}}$ becomes larger in short time and the blowup comes sooner with more singularities. Note that the relative error is not small since $\|\psi(t)\|_{2}=1024/3$. This might suggests that for numerical purpose one needs to use smoother initial data in order to maintain the prescribed accuracy, see the discussions in [39, 38].

5. Conclusion

In the study of the NLS for Bose-Einstein Condensation, the analytic and numerical tools and results we have applied, discovered and reviewed provide good understanding of the modeling equations. On the numerical aspect, the Strang splitting method has been shown to be very accurate in many cases [6]. In the literature this type of splitting schemes apply to a wide range of nonlinear problems including KdV, Maxwell-Dirac and Zakharov systems, Boltzmann equation and Landau damping [5, 38, 10, 19].

Recent theory informs that when the quadratic potential has only positive frequency, the wave function exists locally in time and stable, and when $V$ has large negative frequency, then it can counteract the nonlinear effect. The outcome of the simulations agree with the physics of the BEC under trapping conditions. In the case where the anisotropic quadratic potential has sufficiently higher negative coefficients we observe a dissipative pattern over time, similar to that of the defocusing nonlinearity [12, 32]. The focusing nonlinearity causes an attractive effect on the condensate that can cause it to “blowup”. These are true when $\psi_{0}$ is a gaussian. In Figure 5(a) after short time the linear $V$ starts to take over and there shows scattering like in the linear periodic case. However, when the initial data has low regularity, we observe singularities over very short time. Nevertheless, over a quite long time interval the solutions exhibit multi-fractal structure similar to the linear case [33]. Thus it may be worthwhile to look into the post-blowup behavior of the solutions.

The general understanding is that the BEC mechanism decouples into two states: The ground state from focusing effect and the dispersion from the repulsive interaction. Considering the recent work on BEC with rotation or more generally, the NLS with magnetic fields [37, 4, 58], where some questions are quite open, it would be of interest to continue to study such model under more critical conditions. This investigation, on the analytic and numerical levels, would help explain how the excited sates are formed and how dispersion or scattering can be achieved by appropriately manipulating BEC with potentials, the nonlinearities, and actions of symmetries.

Acknowledgment This work is funded in part by the COSM pilot interdisciplinary project at GSU. The authors would like to thank Dr. M. Edwards for constructive discussions.