Exact Solutions

Exact Solutions#

Here we present the various exact solutions described in gpe/exact_solutions.py

Harmonic Oscillator#

Here we construct an exact solution for a trapped BEC in a harmonic oscillator:

\[\begin{gather*} \psi_0(x) = \sqrt{n_0}e^{-x^2/2\sigma^2}, \qquad V(x) = \frac{m\omega^2x^2}{2} - \frac{\hbar^2}{2m\sigma^2} - gn_0e^{-x^2/\sigma^2}, \qquad \sigma = \sqrt{\frac{\hbar}{m\omega}}. \end{gather*}\]

This potential has the non-linear piece subtracted so that the constructed solution is an exact solution with zero chemical potential (i.e. a stationary state) and energy:

\[\begin{gather*} E = \int_{-\infty}^{\infty}\d{x}\;\frac{-gn_0^2}{2}e^{-2x^2/\sigma^2} = -\sqrt{\frac{\pi}{2}}\frac{gn_0^2\sigma}{2} \end{gather*}\]

%pylab inline --no-import-all
from gpe.imports import *
from gpe.minimize import MinimizeState
from gpe.exact_solutions import HarmonicOscillator

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
[I 03:36:35 numexpr.utils] NumExpr defaulting to 2 threads.
[I 03:36:36 root] Patching zope.interface.document.asReStructuredText to format code

s = HarmonicOscillator(n0=1.0)
s.plot()
assert abs(s.compute_dy_dt(s.copy())[...]).max() < 1e-14
assert np.allclose(s.get_energy(), -s.g * s.n0**2 * s.sigma * np.sqrt(np.pi / 2) / 2)

../_images/3b9da4d1db1ce59d3bc8af8c5a586b782011158e19f3a011cb6e917da654ed50.png

We can also do this in higher dimensions:

\[\begin{gather*} \psi_0(\vect{x}) = \sqrt{n_0}\exp\left(-\sum_{i}\frac{x_i^2}{2\sigma_i^2}\right), \qquad \frac{-\hbar^2}{2m}\nabla^2\psi_0 = \frac{\hbar^2}{2m}\sum_{i} \left(\frac{1}{\sigma_i^2} - \frac{x_i^2}{\sigma_i^4}\right) = \sum_{i} \left(\frac{\hbar^2}{2m\sigma_i^2} - \frac{m\omega_i^2x_i^2}{2}\right) \\ \qquad V(x) = \sum_i\left(\frac{m\omega_i^2x_i^2}{2} - \frac{\hbar^2}{2m\sigma_i^2}\right) - g\abs{\psi_0}^2, \qquad \sigma_{i}^2 = \frac{\hbar}{m\omega_i}. \end{gather*}\]

This potential has the non-linear piece subtracted so that the constructed solution is an exact solution with zero chemical potential (i.e. a stationary state) and energy:

\[\begin{gather*} E = \int_{-\infty}^{\infty}\d{x}\;\frac{-gn_0^2}{2}e^{-2x^2/\sigma^2} = -\sqrt{\frac{\pi}{2}}\frac{gn_0^2\sigma}{2} \end{gather*}\]

Travelling Waves#

The standard GPE admits a family of periodic traveling waves with an analytic form:

\[\begin{gather*} \psi(x, t) = e^{-\I \mu t /\hbar}\psi_v(x_v), \qquad x_v = x - vt \end{gather*}\]

where the wave is moving with speed \(v\). We base our solution on the presentation of Hoefer and El El:2016.

The GPE can be expressed in terms of \(\psi(x, t)\) or in terms of \(\psi_v(x_v) = e^{-\I\mu t/\hbar} \psi_v(x- vt)\):

\[\begin{gather*} \I\hbar\dot{\psi} = -\frac{\hbar^2}{2m}\pdiff[2]{\psi}{x} + g\abs{\psi}^2\psi,\qquad 0 = -\frac{\hbar^2\psi_v''}{2m} + v \I \hbar \psi_v' + (g\abs{\psi_v}^2 - \mu)\psi_v = \left(\frac{\op{p}^2}{2m} - v\op{p} + g\abs{\psi_v}^2 - \mu\right)\psi_v. \end{gather*}\]

Note: this is not the typical Galilean transformation for quantum mechanical systems which includes an extra factor of \(e^{\I (mvx +mv^2t/2)\hbar}\). This additional phase allows one to shift the momentum, completing the square of the kinetic energy, and removing the \(v\op{p}\) linear term at the expense of shifting the chemical potential. The present form, however, is more general, and works with arbitrary dispersion, so we maintain it.

Hoefer and El El:2016 set \(\hbar=m=g=1\) which defines the following dimensions:

\[\begin{gather*} m_\mathrm{unit} = m, \qquad l_\mathrm{unit} = \frac{mg}{\hbar^2}, \qquad t_\mathrm{unit} = \frac{m^3g^2}{\hbar^5}. \end{gather*}\]

Using the Madelung transform:

(6)#\[\begin{gather} \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad u(x, t) = \phi'(x, t),\\ \dot{\rho} + (\rho u)' = 0, \qquad (\rho u)_{,t} + \left(\rho u^2 + \frac{\rho^2}{2}\right)' = \frac{1}{4}\bigl(\rho\, (\log\rho)''\bigr)'. \tag{2.111} \end{gather}\]

The solution can be expressed as:

\[\begin{gather*} \DeclareMathOperator{\sn}{sn} \psi(x, t) = e^{-\I\mu t/\hbar}\psi_v(x_v), \qquad \psi_v(x_v) = \sqrt{\rho_v(x_v)}e^{\I\phi_v(x_v)}, \qquad u(x, t) = \phi_v'(x, t),\\ m = al^2, \qquad \rho_v(x_v) = \rho_{\min} + a \sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad x_v = x - vt, \qquad u(x, t) = v - \frac{C}{\rho(x, t)}, \\ a = \rho_{\max} - \rho_{\min}, \qquad \mu = \rho_{\min} - \frac{v^2}{2} + \frac{\rho_\max + a/m}{2}\qquad C = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}{l},\qquad L = 2lK(m), \qquad Q = \frac{2\pi}{L}. \end{gather*}\]

(Note: \(k=\sqrt{m}\) for \(\sn(z,k)\) in some CASs. We use \(\sn(z;m)\) and \(K(m)\) here)

The full solution (with proper coefficients) is thus:

\[\begin{gather*} \psi(x, t) = e^{-\I\mu t/\hbar}\psi_v(x_v), \qquad \psi_v(x_v) = \sqrt{\rho_v(x_v)}e^{\I\phi_v(x_v)}, \qquad u(x_v) = \frac{\hbar}{m}\phi_v'(x_v) = v - \frac{C}{\rho_v(x_v)},\\ \rho_v(x_v) = \rho_{\min} + a\sn^2\bigl(\frac{x_v}{l};\tfrac{mg}{\hbar^2}al^2\bigr), \qquad x_v = x - vt, \qquad a = \rho_{\max} - \rho_{\min}, \\ \mu = g\rho_\min - \frac{mv^2}{2} + g\frac{\rho_\max +a/m}{2}\qquad C = \frac{\sqrt{\rho_{\min}\rho_{\max}(\hbar^2+mgl^2\rho_{\min})}}{ml}. \end{gather*}\]

\[\begin{gather*} m = al^2, \qquad l^{-1} = \frac{2K(m)}{L} \end{gather*}\]

%pylab inline --no-import-all
from gpe.imports import *
import gpe.exact_solutions

reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves

# Stationary wave
args = dict(n0=0.9, n1=1.0, Lx=10.0, Nx=16)
s0 = TravellingWaves(v_p=0, **args)

# Travelling wave in the lab frame
s1 = TravellingWaves(v_p=0.6047197603, v_x=0, **args)

# Travelling wave in a co-moving frame.
s2 = TravellingWaves(v_p=0.6047197603, **args)

for ss in evolves([s0, s1, s2], t_max=50.0):
    plt.clf()
    for s in ss:
        s.plot()

# s.twist_x
# s.mu, s.get_mu().real, s._twist

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[4], line 16
     12 # Travelling wave in the lab frame
     13 s1 = TravellingWaves(v_p=0.6047197603, v_x=0, **args)
     14 
     15 # Travelling wave in a co-moving frame.
---> 16 s2 = TravellingWaves(v_p=0.6047197603, **args)
     17 
     18 for ss in evolves([s0, s1, s2], t_max=50.0):
     19     plt.clf()

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/exact_solutions.py:191, in TravellingWaves.__init__(self, Nx, Lx, n0, n1, m, g, hbar, v_p, v_x, twist, **kw)
    189 if twist is None:
    190     x = np.linspace(-Lx / 2.0, Lx / 2.0, 10000)
--> 191     _psi, twist = self.psi_exact(x=x, Lx=Lx)
    192 self._twist = twist
    194 if v_x is None:

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/exact_solutions.py:230, in TravellingWaves.psi_exact(self, x, Lx)
    228 theta = [0]
    229 for _n in range(1, len(x)):
--> 230     theta.append(theta[-1] + scipy.integrate.quad(dtheta, x[_n - 1], x[_n])[0])
    231 theta = np.array(theta)
    232 theta -= theta[len(x) // 2]

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/integrate/_quadpack_py.py:479, in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst, complex_func)
    476     return retval
    478 if weight is None:
--> 479     retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
    480                    points)
    481 else:
    482     if points is not None:

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/integrate/_quadpack_py.py:626, in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    624 if points is None:
    625     if infbounds == 0:
--> 626         return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    627     else:
    628         return _quadpack._qagie(func, bound, infbounds, args, full_output, 
    629                                 epsabs, epsrel, limit)

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/exact_solutions.py:225, in TravellingWaves.psi_exact.<locals>.dtheta(x)
    222 _a = self.v_p / self.v_unit
    223 C = self._C * self.n_unit
--> 225 def dtheta(x):
    226     return _a + C / self.n_exact(x)
    228 theta = [0]

KeyboardInterrupt: 

As a simple test, we consider the phonon limit \(a = \epsilon \ll 1\). In this limit, we have

\[\begin{gather*} K(m) \approx \frac{\pi}{2}\left(1 + \frac{m}{4}\right) + \order(m^2), \qquad m = \frac{L^2}{\pi^2}\epsilon, \qquad \sn(\theta, m) = \sin(\theta) + \order(m), \qquad l = \frac{L}{\pi}\left(1 - \frac{L^2}{4\pi^2}\epsilon\right) + \order(\epsilon^2)\\ \end{gather*}\]

\[\begin{gather*} \frac{L = \approx l\pi\left(1 + \frac{m}{4}\right) \end{gather*}\]

consider a \(v=0\) solution with \(\rho_\min = \rho_\max = 1\):

\[\begin{gather*} \rho = 1, \qquad a = m = 0, \qquad \frac{L}{l} = \pi, \qquad u(x_v) = \frac{C}{\rho}, \qquad \theta(x_v) = \frac{C}{\rho}x_v \end{gather*}\]

from gpe.exact_solutions import K

s._C

s = TravellingWaves(n0=1.0, n1=1.0, Lx=np.pi, Nx=64, v_p=0.0)
s.plot()
s.twist_x
s._mu, s.get_mu().real, s._twist

e = EvolverABM(s, dt=0.2 * s.t_scale)
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(1000)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        plt.close("all")
        clear_output(wait=True)

x = s.xyz[0]
c = np.sqrt(s.g * s.n1 / s.m)
v = 0.5 * c
u = np.sqrt(c**2 - v**2)
n0 = s.n1 * (v**2 / c**2)
l = s.hbar / s.m / u
psi_soliton = v / c * 1j + u / c * np.tanh(x / l)
twist = np.angle(psi_soliton)[-1] - np.angle(psi_soliton)[0]
s1 = TravellingWaves(n0=n0, n1=s.n1, Lx=s.basis.Lx, Nx=s.basis.Nx, v_p=0, twist=twist)
s1[...] = psi_soliton
s1.plot()

e = EvolverABM(s, dt=0.2 * s.t_scale)
e1 = EvolverABM(s1, dt=0.2 * s1.t_scale)
pe = None
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(1000)
        e1.evolve(1000)
        plt.clf()
        e.y.plot()
        e1.y.plot()
        display(plt.gcf())
        plt.close("all")
        clear_output(wait=True)

Bright Soliton#

Here we demonstrate the analytic bright-soliton for a GPE with attractive interactions. This object moves at a specified speed \(v\), and we consider the solution in three frames: comoving (the soliton is stationary), lab frame (the soliton makes a single oscillation in time \(T=L/v\)) and a frame moving backwards at the same speed (the soliton crosses the box twice in time \(T\)).

This forms the basis of gpe/tests/test_bec.py::test_StateTwist_x_v_x and verifies that the frame velocity argument of StateTwist_x works. Because the density vanishes at the boundaries, this does not test the twist.

\[\begin{gather*} \I\hbar\dot{\psi}(x, t) = \frac{-\hbar^2\psi''(x, t)}{2m} + g\abs{\psi}^2\psi,\\ \psi(x, t) = \frac{\sqrt{n_0}}{\cosh\bigl(\eta (x-vt)\bigr)} \exp\left\{\frac{\left(\mu + \tfrac{mv^2}{2}\right)t - vx}{\I\hbar}\right\}, \qquad \mu = \frac{-\hbar^2\eta^2}{2m}, \\ gn_0 = 2\mu = \frac{-\hbar^2\eta^2}{m}, \qquad \sigma = \frac{1}{\eta}. \end{gather*}\]

import sympy
from sympy import Eq, var, I, exp, cosh, sqrt, Abs

x, t, g = var(["x", "t", "g"], real=True)
eta, n0, mu, m, hbar = var(["eta", "n_0", "mu", "m", "hbar"], positive=True)
mu = -(hbar**2) * eta**2 / 2 / m
g = 2 * mu / n0
psi = exp(mu * t / I / hbar) * sqrt(n_0) / cosh(eta * x)
n = Abs(psi) ** 2
display(Eq(sympy.S("psi"), psi))
Hpsi = -(hbar**2) * psi.diff(x, x) / 2 / m + g * n * psi
display(Eq(sympy.S("H*psi"), Hpsi))
(Hpsi - I * hbar * psi.diff(t)).simplify()

%pylab inline --no-import-all
from gpe.imports import *
import gpe.exact_solutions

reload(gpe.exact_solutions)
from gpe.exact_solutions import BrightSoliton

Nx = 32
Lx = 10.0
v = 2.0

args = dict(Nx=Nx, Lx=Lx, v=v, sigma=1)
s0 = BrightSoliton(v_x=0, **args)
s1 = BrightSoliton(v_x=-v, **args)
s2 = BrightSoliton(v_x=v, **args)

for ss in evolves([s0, s1, s2], t_max=Lx / v):
    plt.clf()
    for s in ss:
        s.plot()

1D GPE (NLSEQ)#

Here we consider the exact solutions to the 1D GPE via the inverse scattering method (ISM). These follow from a so-called Lax representation, which can be expressed as follows. Consider two linear operators \(\op{L}(\lambda)\) and

\[\begin{gather*} \op{M}(\lambda) = \I\diff{}{t} + V(x, t, \lambda) \end{gather*}\]

that commute:

\[\begin{gather*} [\op{L}(\lambda), \op{M}(\lambda)] = 0. \end{gather*}\]

\[\begin{gather*} \mat{L} = \I\mat{1}\pdiff{}{x} + \begin{pmatrix} 0 & u^*\\ u & 0 \end{pmatrix},\\ \mat{M} = -\mat{1}\pdiff[2]{}{x} + \begin{pmatrix} \abs{u}^2 & \I u_{,x}^*\\ -\I u_{,x} & -\abs{u}^2 \end{pmatrix}. \end{gather*}\]

differential equation

\[\begin{gather*} \I \dot{\mat{L}} = [\mat{L}, \mat{A}] \end{gather*}\]

We start by expressing the problem as

\[\begin{gather*} \I \dot{u} + u'' + 2 \abs{u}^2 u = 0, \qquad u(x, t=0) = u_0(x). \end{gather*}\]

Define

\[\begin{gather*} \mat{q}_0(x) = \begin{pmatrix} 0 & u_0(x)\\ u_0^*(x) & 0 \end{pmatrix}. \end{gather*}\]