import mmf_setup
mmf_setup.nbinit()
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
Two-Component BECs (BEC2)#
Here and in the file bec2.py we present code that solves the GPE
for a two-component system. The energy density and equations of motion are:
\[\begin{gather*}
\mathcal{E}[\Psi] = \Psi^\dagger\begin{pmatrix}
\frac{-\hbar^2\nabla^2}{2m} + V_{a} & V_{ab}\\
V_{ab}^\dagger & \frac{-\hbar^2\nabla^2}{2m} + V_{b}\\
\end{pmatrix}\Psi
+
\mathcal{E}(n_a, n_b),
\qquad
\Psi = \begin{pmatrix}
\psi_a\\
\psi_b
\end{pmatrix}\\
-\I\hbar \pdiff{}{t}\Psi =
\begin{pmatrix}
\frac{-\hbar^2\nabla^2}{2m} + V_{a}
+ \mathcal{E}_{,n_a}(n_a, n_b) & V_{ab}\\
V_{ab}^\dagger & \frac{-\hbar^2\nabla^2}{2m} + V_{b}
+ \mathcal{E}_{,n_b}(n_a, n_b)\\
\end{pmatrix}\Psi.
\end{gather*}\]
For the standard GPE, the equation of state is:
\[\begin{gather*}
\mathcal{E}(n_a, n_b) = \frac{g_{aa}}{2} n_a^2 + \frac{g_{aa}}{2} n_a^2 + g_{ab} n_an_b.
\end{gather*}\]
This structure is supported in an augmented state.data array of shape (2, N_x, N_y, N_z) where the first index refers to indices a or b. Note now that functions like integrate(), get_N(), and braket() now return arrays of shape (2,) instead of scalars.
Twisted Boundary Conditions#
Instead of a periodic box, one can work with twisted boundary conditions whereby
\[\begin{gather*}
\psi(x+L) = e^{\I\phi}\psi(x).
\end{gather*}\]
To implement such boundary conditions, we use the periodic function
\[\begin{gather*}
\tilde{\psi}(x) = e^{-\I\phi x/L}\psi(x)
\end{gather*}\]
when computing the FFT. To compute derivatives, we must thus first apply this shift to the wavefunction to get a periodic function, then take the FFT, then use momenta that have been shifted by \(-\phi/L\) etc.
%pylab inline --no-import-all
from IPython.display import display, clear_output
from gpe import bec2
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
[I 03:36:21 root] Patching zope.interface.document.asReStructuredText to format code
[I 03:36:21 numexpr.utils] NumExpr defaulting to 2 threads.
reload(bec2)
from gpe.bec2 import State, u
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 reload(bec2)
2 from gpe.bec2 import State, u
NameError: name 'reload' is not defined
s = State(Nxyz=(128,), Lxyz=(40 * u.micron,), N=1e5)
s.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 s = State(Nxyz=(128,), Lxyz=(40 * u.micron,), N=1e5)
2 s.plot()
NameError: name 'State' is not defined
Testing#
To test the code, we will set \(g=0\) and use the exact solution for the Harmonic Oscillator:
\[\begin{gather*}
\psi(x) \propto e^{-(x/a)^2/2}, \qquad
a^2 = \frac{\hbar}{m\omega}
\end{gather*}\]
def get_err(N=128, L=24 * u.micron):
s = State(Nxyz=(N,), Lxyz=(L,), N=1e5)
s.gs = 0, 0, 0
a = np.sqrt(u.hbar / u.m / s.ws[0])
x = s.xyz[0]
psi_0 = np.exp(-((x / a) ** 2) / 2.0)
s[...] = psi_0[None, :]
s.normalize()
dy = s.empty()
s.compute_dy_dt(dy=dy, subtract_mu=True)
return abs(dy[...]).max()
Ns = 2 ** np.arange(2, 8)
errs = map(get_err, Ns)
plt.semilogy(Ns, errs, "-+")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 def get_err(N=128, L=24 * u.micron):
2 s = State(Nxyz=(N,), Lxyz=(L,), N=1e5)
3 s.gs = 0, 0, 0
4 a = np.sqrt(u.hbar / u.m / s.ws[0])
NameError: name 'u' is not defined
Why are \(L=23\)microns and \(N=2^6\) optimal?
s = State(Nxyz=(46,), Lxyz=(23 * u.micron,))
a = np.sqrt(u.hbar / u.m / s.ws[0])
L, N = s.Lxyz[0], s.Nxyz[0]
k_max = np.pi * (N - 2) / L # For Khalid...
print k_max, np.max(s.kxyz[0])
print np.exp(-((L / 2 / a) ** 2) / 2) # Wavefunction drops by factor of macheps
Cell In[7], line 5
print k_max, np.max(s.kxyz[0])
^
SyntaxError: Missing parentheses in call to 'print'. Did you mean print(...)?
psi_0 = s.xyz[0] * np.exp(-((s.xyz[0] / a) ** 2) / 2)
plt.semilogy(
np.fft.fftshift(s.kxyz[0]).ravel(),
np.fft.fftshift(abs(np.fft.fft(psi_0))).ravel(),
"-+",
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 psi_0 = s.xyz[0] * np.exp(-((s.xyz[0] / a) ** 2) / 2)
2 plt.semilogy(
3 np.fft.fftshift(s.kxyz[0]).ravel(),
4 np.fft.fftshift(abs(np.fft.fft(psi_0))).ravel(),
NameError: name 's' is not defined
So we see that for the ground state \(k\) needs to go up to \(6\).
Exact Solution with Interactions#
Consider now a state with \(\psi_a = e^{-(x/a)^2/2}\) and \(\psi_b = e^{-(x/b)^2/2}\).
\[\begin{gather*}
\psi_a''(x)
\end{gather*}\]
%pylab inline --no-import-all
from IPython.display import display, clear_output
from gpe import bec2
reload(bec2)
from gpe.bec2 import State, u
s = State(Nxyz=(64,), Lxyz=(23 * u.micron,), N=1e5)
a = np.sqrt(u.hbar / u.m / s.ws[0])
x = s.xyz[0]
psi_0 = np.exp(-((x / a) ** 2) / 2)
class State1(State):
def __init__(self, *v, **kw):
State.__init__(self, *v, **kw)
a = np.sqrt(u.hbar / u.m / self.ws[0])
x = self.xyz[0]
k = 1.0 / 2.0 / a ** 2
psi_0 = 4.0 * np.exp(-((x / a) ** 2) / 2)
n_0 = abs(psi_0) ** 2
g_aa, g_bb, g_ab = self.gs
self._V_ext = (
u.hbar ** 2 / 2.0 / u.m * (4 * (k * x) ** 2 - 2 * k) - g_aa * n_0,
u.hbar ** 2 / 2.0 / u.m * (4 * (k * x) ** 2 - 2 * k) - g_bb * n_0,
0 * n_0,
)
self.data[...] = psi_0
self.get_Vext = lambda: self._V_ext
self.pre_evolve_hook()
s = State1(Nxyz=(64,), Lxyz=(23 * u.micron,))
s.plot()
plt.plot(x, s.get_Vext()[0])
dy = s.empty()
s.compute_dy_dt(dy=dy, subtract_mu=False)
abs(dy[...]).max()
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 5
1 get_ipython().run_line_magic('pylab', 'inline --no-import-all')
2 from IPython.display import display, clear_output
3 from gpe import bec2
4
----> 5 reload(bec2)
6 from gpe.bec2 import State, u
7
8 s = State(Nxyz=(64,), Lxyz=(23 * u.micron,), N=1e5)
NameError: name 'reload' is not defined
from mmfutils.contexts import NoInterrupt
from pytimeode.evolvers import EvolverSplit, EvolverABM
from IPython.display import display, clear_output
s = State1(Nxyz=(64 * 4,), Lxyz=(23 * u.micron,))
assert np.allclose(s._N, s.get_N())
s[...] = 1.0
s.normalize()
s.cooling_phase = 1j
E_max = u.hbar ** 2 * np.abs(s.kxyz).max() ** 2 / 2.0 / u.m
# e = EvolverSplit(s, dt=0.01*u.hbar/E_max, normalize=True)
e = EvolverABM(s, dt=0.1 * u.hbar / E_max, normalize=True)
with NoInterrupt(ignore=True) as interrupted:
while e.y.t < 4 * u.ms and not interrupted:
e.evolve(100)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 5
1 from mmfutils.contexts import NoInterrupt
2 from pytimeode.evolvers import EvolverSplit, EvolverABM
3 from IPython.display import display, clear_output
4
----> 5 s = State1(Nxyz=(64 * 4,), Lxyz=(23 * u.micron,))
6 assert np.allclose(s._N, s.get_N())
7
8 s[...] = 1.0
NameError: name 'State1' is not defined
from mmfutils.contexts import NoInterrupt
from pytimeode.evolvers import EvolverSplit, EvolverABM
from IPython.display import display, clear_output
s = State1(Nxyz=(64 * 4,), Lxyz=(23 * u.micron,))
s *= np.sign(s.xyz[0] - 0.5)
s.cooling_phase = 1 + 0.01j
E_max = u.hbar ** 2 * np.abs(s.kxyz).max() ** 2 / 2.0 / u.m
# e = EvolverSplit(s, dt=0.01*u.hbar/E_max, normalize=True)
e = EvolverABM(s, dt=0.5 * u.hbar / E_max, normalize=True)
with NoInterrupt(ignore=True) as interrupted:
while e.y.t < 40 * u.ms and not interrupted:
e.evolve(100)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 5
1 from mmfutils.contexts import NoInterrupt
2 from pytimeode.evolvers import EvolverSplit, EvolverABM
3 from IPython.display import display, clear_output
4
----> 5 s = State1(Nxyz=(64 * 4,), Lxyz=(23 * u.micron,))
6 s *= np.sign(s.xyz[0] - 0.5)
7 s.cooling_phase = 1 + 0.01j
8
NameError: name 'State1' is not defined
2D Scissor Modes (incomplete)#
We base the physical parameters on a system like [Marago:2001] so we can explore damping of the normal modes of the system. They have trapping frequencies of \(\omega_y=\omega_z = 128\)Hz and \(\omega_x = \sqrt{8}\omega_y\) and \(N = 2\times 10^4\) particles. In the TF approximation at \(T=0\) this corresponds to \(\mu = \)
[Marago:2001]: http://dx.doi.org/10.1103/PhysRevLett.86.3938 (Onofrio Marag`o, Gerald Hechenblaikner, Eleanor Hodby, and Christopher Foot, “Temperature Dependence of Damping and Frequency Shifts of the Scissors Mode of a Trapped Bose-Einstein Condensate”, Phys. Rev. Lett. 86, 3938–3941 (2001) )
%pylab inline --no-import-all
from IPython.display import display, clear_output
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
from gpe import bec
reload(bec)
from pytimeode.evolvers import EvolverABM
from mmfutils.contexts import NoInterrupt
from gpe.bec import State, u
s = State()
s.cooling_phase = 1j
s.t = -100 * u.ms
e = EvolverABM(s, dt=0.001)
with NoInterrupt(ignore=True) as interrupted:
while not interrupted:
e.evolve(500)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 3
1 from gpe import bec
2
----> 3 reload(bec)
4 from pytimeode.evolvers import EvolverABM
5 from mmfutils.contexts import NoInterrupt
6 from gpe.bec import State, u
NameError: name 'reload' is not defined
s = e.get_y()
n = s.get_density()
x, y, z = s.xyz
plt.plot(x.ravel(), n.sum(axis=-1).sum(axis=-1))
s.mu, s._mu
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 s = e.get_y()
2 n = s.get_density()
3 x, y, z = s.xyz
4 plt.plot(x.ravel(), n.sum(axis=-1).sum(axis=-1))
NameError: name 'e' is not defined
Bright Solitons#
There is nice integrable 1D model that is useful for testing code in the case of attractive interactions (the so-called focusing NLSEQ):
\[\begin{gather*}
\I\hbar \partial_t \Psi = - \frac{\hbar^2\nabla^2}{2m}\Psi
+ \begin{pmatrix}
-gn_{+} & \Omega(t)\\
\Omega(t) & -gn_{+}
\end{pmatrix}\cdot\Psi
\end{gather*}\]
where \(n_{+} = n_a + n_b\). This corresponds to a coupled set of GPEs with equal attractive interactions \(g_{aa} = g_{bb} = g_{ab} = -g\).
\[\begin{gather*}
\DeclareMathOperator{\sech}{sech}
\Psi(x, t) = \frac{1}{2}\begin{pmatrix}
\cos\theta e^{-\I\Gamma(t)} + \sin\theta e^{\I\Gamma(t)}\\
\cos\theta e^{-\I\Gamma(t)} - \sin\theta e^{\I\Gamma(t)}
\end{pmatrix}e^{\I g b^2 t/4} b\sech\frac{b \sqrt{gm} x}{\hbar}, \qquad
-\frac{\hbar^2 \psi''(x)}{2m} - g\abs{\psi}^2\psi = -\frac{b^2 g}{2}\psi, \\
\Gamma(t) = \int_{0}^{t}\Omega(t')\d{t'}
\end{gather*}\]
\[\begin{gather*}
\Psi(x, t) = \left(\frac{\cos\theta e^{-\I\Gamma(t)}}{2}
\begin{pmatrix}1\\1\end{pmatrix}
+
\frac{\sin\theta e^{\I\Gamma(t)}}{2}
\begin{pmatrix}1\\ -1\end{pmatrix} \right)
e^{\I g b^2 t/4} b\sech(b x\sqrt{gm/\hbar^2}), \qquad
-\frac{\hbar^2 \psi''(x)}{2m} - g\abs{\psi}^2\psi = -\frac{b^2 g}{2}\psi, \\
\Gamma(t) = \int_{0}^{t}\Omega(t')\d{t'}
\end{gather*}\]
We have the following dimensions:
\[\begin{gather*}
[\hbar] = \frac{M D^2}{T}, \qquad
[2m] = M, \qquad
[gn] = [E] = \frac{M D^2}{T^2}, \qquad
[g] = [VE] = \frac{M D^5}{T^2}, \qquad
[\Omega] = [E] = \frac{M D^2}{T^2}\\
\left[\frac{2mg}{\hbar^2}\right] = D, \qquad
\left[\frac{2m\Omega}{\hbar^2}\right] = \frac{1}{D^2}, \qquad
\left[\frac{t\hbar}{2m}\right] = D^2.
\end{gather*}\]
In the paper, units are scaled so that \(2m = \hbar = 1\). This means that \(g\) provides the length scale for the problem.
from numpy import cos, sin, exp
theta, Gamma = np.random.random(2)
(cos(theta) * exp(-1j * Gamma) + sin(theta) * np.exp(1j * Gamma)) / 2, (
cos(theta + Gamma) - 1j * sin(theta - Gamma)
) / 4
(np.complex128(0.5510033158108356-0.020159308102452346j),
np.complex128(0.24700595507245932-0.01593005327747076j))
%pylab inline --no-import-all
from IPython.display import display, clear_output
from pytimeode.evolvers import EvolverABM
from mmfutils.contexts import NoInterrupt
from gpe import bec2
reload(bec2)
from gpe.bec2 import State2, u
class StateBS(State2):
t0 = u.micron ** 2 * 2 * u.m / u.hbar
def get_Vext(self):
Va = Vb = 0 * self.xyz[0]
return (Va, Vb, self.Omega)
# Paper dimensions
g_ = 2.0
Omega_ = 2.0
theta_ = np.pi / 4.0
b_ = 2.0
s = StateBS(Nxyz=(128,), Lxyz=(4 * u.micron,), N=1e5)
s = StateBS(
Nxyz=(64 * 4,),
Lxyz=(16 * u.micron,),
N=1e5,
gs=(-g_ * u.hbar ** 2 / 2.0 / u.m,) * 3,
Omega=Omega_ * u.hbar ** 2 / 2.0 / u.m,
cooling_phase=1.0j,
)
# Fig 9
g = -g_ * u.hbar ** 2 / 2.0 / u.m
s = StateBS(
Nxyz=(64,),
Lxyz=(16 * u.micron,),
N=1e5,
gs=(g * 1.01, g * 1.01, g * 0.99),
Omega=Omega_ * u.hbar ** 2 / 2.0 / u.m,
cooling_phase=1.0,
)
x = s.xyz[0]
s[0, ...] = (
(np.cos(theta_) + np.sin(theta_)) / 2.0 * b_ / np.cosh(b_ * np.sqrt(g_) * x / 2.0)
)
s[1, ...] = (
(np.cos(theta_) - np.sin(theta_)) / 2.0 * b_ / np.cosh(b_ * np.sqrt(g_) * x / 2.0)
)
# s[...] += np.random.random(s.shape)*0.1
s.plot()
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 7
3 from pytimeode.evolvers import EvolverABM
4 from mmfutils.contexts import NoInterrupt
5 from gpe import bec2
6
----> 7 reload(bec2)
8 from gpe.bec2 import State2, u
9
10
NameError: name 'reload' is not defined
e = EvolverABM(s, dt=0.5 / s.E_max / u.hbar)
with NoInterrupt(ignore=True) as interrupted:
while not interrupted:
e.evolve(10)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 e = EvolverABM(s, dt=0.5 / s.E_max / u.hbar)
2 with NoInterrupt(ignore=True) as interrupted:
3 while not interrupted:
4 e.evolve(10)
NameError: name 's' is not defined
import sympy
sympy.init_session()
var("a, theta, b, g, x")
psi = a / cosh(x * b * sqrt(g) / 2)
n = psi ** 2
r1, r2 = (-psi.diff(x, x) / 2).trigsimp().simplify(), (g * n * psi).simplify()
r1, r2
IPython console for SymPy 1.14.0 (Python 3.14.4-64-bit) (ground types: python)
These commands were executed:
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()
Documentation can be found at https://docs.sympy.org/1.14.0/
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[20], line 6
2 var("a, theta, b, g, x")
3 psi = a / cosh(x * b * sqrt(g) / 2)
4 n = psi ** 2
5 r1, r2 = (-psi.diff(x, x) / 2).trigsimp().simplify(), (g * n * psi).simplify()
----> 6 r1, r2
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/decorator/__init__.py:247, in decorate.<locals>.fun(*args, **kw)
245 if not kwsyntax:
246 args, kw = fix(args, kw, sig)
--> 247 return caller(func, *(extras + args), **kw)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/interactive/printing.py:179, in _init_ipython_printing.<locals>._print_latex_png(o)
177 s = '$\\displaystyle %s$' % s
178 try:
--> 179 return _preview_wrapper(s)
180 except RuntimeError as e:
181 debug('preview failed with:', repr(e),
182 ' Falling back to matplotlib backend')
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/interactive/printing.py:90, in _init_ipython_printing.<locals>._preview_wrapper(o)
88 exprbuffer = BytesIO()
89 try:
---> 90 preview(o, output='png', viewer='BytesIO', euler=euler,
91 outputbuffer=exprbuffer, extra_preamble=extra_preamble,
92 dvioptions=dvioptions, fontsize=fontsize)
93 except Exception as e:
94 # IPython swallows exceptions
95 debug("png printing:", "_preview_wrapper exception raised:",
96 repr(e))
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/preview.py:311, in preview(expr, output, viewer, euler, packages, filename, outputbuffer, preamble, dvioptions, outputTexFile, extra_preamble, fontsize, **latex_settings)
308 raise RuntimeError("latex program is not installed")
310 try:
--> 311 _check_output_no_window(
312 ['latex', '-halt-on-error', '-interaction=nonstopmode',
313 'texput.tex'],
314 cwd=workdir,
315 stderr=STDOUT)
316 except CalledProcessError as e:
317 raise RuntimeError(
318 "'latex' exited abnormally with the following output:\n%s" %
319 e.output)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/preview.py:26, in _check_output_no_window(*args, **kwargs)
24 else:
25 creation_flag = 0 # Default value
---> 26 return check_output(*args, creationflags=creation_flag, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/subprocess.py:473, in check_output(timeout, *popenargs, **kwargs)
470 empty = b''
471 kwargs['input'] = empty
--> 473 return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
474 **kwargs).stdout
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/subprocess.py:557, in run(input, capture_output, timeout, check, *popenargs, **kwargs)
555 with Popen(*popenargs, **kwargs) as process:
556 try:
--> 557 stdout, stderr = process.communicate(input, timeout=timeout)
558 except TimeoutExpired as exc:
559 process.kill()
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/subprocess.py:1208, in Popen.communicate(self, input, timeout)
1206 self._stdin_write(input)
1207 elif self.stdout:
-> 1208 stdout = self.stdout.read()
1209 self.stdout.close()
1210 elif self.stderr:
KeyboardInterrupt:
sympy.expand_trig