import mmf_setup;mmf_setup.nbinit()
%pylab inline --no-import-all
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.
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
Travelling Waves#
We start with a brief review of Galilean Covariance. For details see:
Boosts#
Consider a solution to the GPE in 1D:
\[
\I\hbar \dot{\psi}(x, t) =
-\frac{\hbar^2\psi''(x, t)}{2m}
+ g\abs{\psi(x, t)}^2\psi(x, t).
\]
Suppose we wish to study a solution that travels with velocity \(v\). I.e. something which should be static when considered as a function of \(x_v = x - vt\). As discussed in Galilean Covariance, we may consider several transformations of the form:
\[
\psi(x, t) = e^{\I\phi(x_v, t)}\psi_v(x_v, t)
\]
Simple boost (notation \(\psi_v(x_v, t)\)):
\[\begin{split} \hbar\phi(x_v, t) = 0, \\ \I\hbar\dot{\psi}_v(x_v, t) = -\frac{\hbar^2\psi_v''(x_v, t)}{2m} + \overbrace{\I\hbar v \psi_v'(x_v)}^{-v\op{p}\psi_v} + g\abs{\psi_v(x_v, t)}^2\psi_v(x_v, t). \end{split}\]Full Galilean Transformation (notation with a capital \(V=v\), \(\psi_V(x_V, t)\)):
\[\begin{split} \hbar\phi(x_V, t) = mVx_V + \tfrac{m}{2}V^2t, \\ \I\hbar\dot{\psi}_V(x_V, t) = -\frac{\hbar^2\psi_V''(x_V, t)}{2m} + g\abs{\psi_V(x_V, t)}^2\psi_V(x_v, t). \end{split}\]
The full transformation makes the covariance of the original equation manifest - the form of the equation stays the same. However, it does not simplify if the original equations are not Galilean covariant (e.g. if there is a modified dispersion due to a spin-orbit coupling.)
Arbitrary Dispersion#
The previous relationships can be expressed for arbitrary dispersion, which reduce to the previous equations when \(E(p) = p^2/2m\):
\[\begin{split}
\I\hbar \dot{\psi}(x, t) =
E(\op{p})\psi(x, t)
+ g\abs{\psi(x, t)}^2\psi(x, t),\\
\psi(x, t) = e^{\I\phi(x_v, t)}\psi_v(x_v, t) = e^{\I\phi(x_V, t)}\psi_V(x_V, t).
\end{split}\]
Simple boost:
\[\begin{split} \hbar\phi(x_v, t) = 0, \\ \I\hbar\dot{\psi}_v(x_v, t) = \bigl(E(\op{p})-v\op{p}\bigr)\psi_v(x_v,t) + g\abs{\psi_v(x_v, t)}^2\psi_v(x_v, t). \end{split}\]Full Galilean Transformation:
\[\begin{split} \hbar\phi(x_V, t) = mVx_V + \tfrac{m}{2}V^2t, \\ \I\hbar\dot{\psi}_V(x_V, t) = \left(E(\op{p} + mV)-V\op{p} - \frac{mV^2}{2}\right)\psi_V(x_V, t) + g\abs{\psi_V(x_V, t)}^2\psi_V(x_V, t). \end{split}\]
Madelung Transform: Quantum Hydrodynamics#
Implementing a Madelung transformation we obtain the hydrodynamic equivalent by equating real and imaginary parts (multiplying through by \(-2\sqrt{\rho}\) and \(\sqrt{\rho}/\rho\) in each case:
\[\begin{split}
\psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad u(x, t) = \tfrac{\hbar}{m} \phi'(x, t),\\
\dot{\rho} + (\rho u)' = 0, \qquad
\mu - \hbar\dot{\phi} = g\rho + \frac{m u^2}{2}
- \frac{\hbar^2}{2m}\frac{\sqrt{\rho}''}{\sqrt{\rho}}.
\end{split}\]
Taking the gradient of the last equation, multiplying through by \(\rho/m\) and using the first continuity equation to replace \(\rho\dot{u} = (\rho u)_{,t} + u (\rho u)'\) we obtain:
\[
(\rho u)_{,t} + \left(\rho u^2 + \frac{g\rho^2}{2m}\right)' = \frac{\rho}{m}\left(
\frac{\hbar^2}{2m}\frac{\sqrt{\rho}''}{\sqrt{\rho}}
\right)'
= \frac{\hbar^2}{4m}\bigl(\rho\, (\log\rho)''\bigr)'.
\]
In the first form, we can identify the roles of the force \(F = -\nabla V\) which includes both the mean-field pressure and any external potential, and the effect of the quantum pressure:
\[
(\rho u)_{,t} + (\rho u^2)' = \frac{\rho}{m}F,
\qquad
F = - \nabla\left(
(g\rho + V_{\mathrm{ext}})
- \frac{\hbar^2}{2m}\frac{\sqrt{\rho}''}{\sqrt{\rho}}
\right).
\]
The left-hand-side is simply the covariant fluid derivative of the momentum \(m \rho^{-1}\d \vect{p}/\d t\) so Newton’s law is manifest. The second form is the hydrodynamic form used in Hoefer and El.
Moving Frame#
If we implement the full Galilean transformation, then the equations of motion are invariant, but if we implement the partial transform then the equations are slightly modified:
\[
\dot{\rho}_v
+ \Bigl((u_v\overbrace{-v}^{\text{boost}})\rho_v\Bigr)' = 0, \qquad
\mu \overbrace{+ \frac{mv^2}{2}}^{\text{boost}}
- \hbar \dot{\phi}_v
= g \rho_v
+ m\frac{(u_v\overbrace{-v}^{\text{boost}})^2}{2}
- \frac{\hbar^2}{2m}\frac{\sqrt{\rho_v}''}{\sqrt{\rho_v}}.
\]
In contrast, the full Galilean transformation also shifts \(u_v \rightarrow u_v + v\) and the chemical potential \(\mu \rightarrow \mu - mv^2/2\) so that the equations revert to their original form.
Note: for pedagogical purposes, we derive the first equation explicitly. The partial Galilean transformation just changes the variable \(x \rightarrow x_v\):
\[
x_v = x - vt, \qquad
\rho_v(x_v, t) = \rho(x, t) = \rho(x_v + vt, t), \qquad
u_v(x_v, t) = u(x, t) = u(x_v+vt, t),
\]
hence, the partials are related by
\[
\dot{\rho}_v = \dot{\rho} + v\rho', \qquad
\rho' = \rho_v', \qquad
u_v' = u'.
\]
This follows from the following notation:
\[
\dot{\rho} = \left.\pdiff{\rho(x, t)}{t}\right|_{x},\qquad
\rho' = \left.\pdiff{\rho(x, t)}{x}\right|_{t}, \qquad
\dot{\rho}_v = \left.\pdiff{\rho_v(x_v, t)}{t}\right|_{x_v}, \qquad
\rho_v' = \left.\pdiff{\rho_v(x_v, t)}{x_v}\right|_{t}.
\]
Thus, the continuity equation is as advertised:
\[
\overbrace{\dot{\rho}}^{\dot{\rho}_v - v \rho_v'} + \overbrace{(\rho u)'}^{(\rho_v u_v)'} =
\dot{\rho}_v + \Bigl((u_v-v) \rho_v\Bigr)' = 0.
\]
Arbitrary Dispersion#
In the presence of arbitrary dispersion \(E(p)\), the second quantum hydrodynamic equation (the Bernoulli equation) is significantly complicated outside of homogeneous matter (see Khamehchi:2017), however, the continuity equation remains the same. Thus, in a simply boosted frame, we still have
\[
\dot{\rho}_v + \Bigl((u_v-v)\rho_v\Bigr)' = 0.
\]
Traveling Waves#
Consider a solution where the density profile moves but maintains its shape. We can then determine the phase from the continuity equation:
\[\begin{split}
\rho(x, t) = \rho_v(\overbrace{x - vt}^{x_v}), \qquad
\dot{\rho} + (\rho u)' = \Bigr(\overbrace{\bigl(u_v(x_v) - v\bigr)\rho_v(x_v)}^{-C}\Bigr)' = 0,\\
u(x, t) = u_v(x_v) = v - \frac{C}{\rho(x, t)} = v - \frac{C}{\rho_v(x_v)}.
\end{split}\]
This is simply a re-derivation of the boosted continuity equation above. The wave-function can thus be rendered time-independent with an appropriately chosen chemical potential \(\mu\):
\[\begin{split}
\psi(x, t) = \sqrt{\rho_v(x_v)}e^{\I[\phi(x_v) - \mu t/\hbar]} = \psi_v(x_v)e^{-\I\mu t/\hbar}, \qquad
\phi'(x_v) = u_v(x_v) = v - \frac{C}{\rho_v(x_v)},\\
\mu\psi_v(x_v) = \bigl(E(\op{p}) - v\op{p}\bigr)\psi_v(x_v)
+ g\abs{\psi_v(x_v)}^2\psi_v(x_v).
\end{split}\]
In the case of the Galilean covariant GPE we can go further and use the full Galilean transformation:
\[\begin{split}
\psi(x, t) = \sqrt{\rho_V(x_V)}e^{\I[\phi(x_V) - \mu t/\hbar]}e^{\I\bigl[mVx_V + \tfrac{m}{2}V^2t \bigr]/\hbar}
= \psi_{V}(x_V)e^{\I\bigl[mVx_V + \bigl(\tfrac{m}{2}V^2 - \mu\bigr)t \bigr]/\hbar},\\
\mu\psi_V(x_V) = -\frac{\hbar^2\psi_V''(x_V)}{2m}
+ g\abs{\psi_V(x_V)}^2\psi_V(x_V).
\end{split}\]
It might look like one can thus express traveling waves in a Galilean covariant theory such as the GPE as purely real stationary solutions in an appropriately moving frame, but this only describes a subset of possible solutions because the non-linear piece \(g\abs{\psi_V}^2\) couples both real and imaginary parts.
This subset of solutions can be found by direct integration starting from the initial condition that \(\psi_v(0)>0\) is the maximum of the wave with \(\psi_V'(0) = 0\). Since the wave must be convex at this point, we must have:
\[
\frac{\hbar^2\psi_V''(0)}{2m \psi_V(0)} = g\psi_V^2(0) - \mu \leq 0.
\]
We thus have a continuous family of solutions with \(0 < \psi_V(0) < \sqrt{\mu/g}\).
Numerical Example#
Here we numerically solve this equation using the scipy IVP solver, then compare this with our exact solution as implemented in the gpe.exact_solutions.TravellingWave class.
from scipy.integrate import solve_ivp
class ODE:
hbar = 1.0
g = 1.0
mu = 1.0
m = 1.0
@property
def psi_max(self):
if self.mu/self.g > 0:
return np.sqrt(self.mu/self.g)
else:
return 1.0
def f(self, x, y):
psi, dpsi = y
ddpsi = 2*self.m*(self.g*np.abs(psi)**2 - self.mu)*psi/self.hbar**2
return (dpsi, ddpsi)
def event(self, x, y):
"""Termination when dpsi=0 again"""
if x <= 0:
return -1
psi, dpsi = y
return dpsi
event.terminal = True
event.direction = -1
def solve(self, fraction=0.9, max_step=0.1):
x_max = 10.0
psi0 = fraction*self.psi_max
ivp = solve_ivp(self.f, (0, x_max), (psi0, 0),
max_step=max_step, events=[self.event])
return ivp
ode = ODE()
ivp = ode.solve(0.2)
x, psi = ivp.t, ivp.y[0]
plt.figure(figsize=(10,3))
ax1 = plt.subplot(121)
ax1.plot(x, psi)
ax1.set(xlabel="x", ylabel=r'$\psi$')
ax2 = plt.subplot(122)
ax2.plot(x, psi**2)
ax2.set(xlabel="x", ylabel=r'$n=\psi^2$')
L = ivp.t_events[0]
Here we construct the same solution from our exact solution.
# Note: there is some error in the period for longer waves.
from gpe.imports import *
from gpe.exact_solutions import TravellingWaves
n1 = (psi**2).max()
n0 = 1e-8 # Should be zero, but for numerical reasons we just make small.
s = TravellingWaves(n1=n1, n0=n0, g=ode.g, Lx=ivp.t[-1]/2.0)
s.plot()
plt.plot(ivp.t-s.Lx/2, ivp.y[0]**2)
print("psi_0/sqrt(mu/g) = {}".format(np.sqrt(n1/(s.get_mu().real/s.g))))
[I 03:38:41 numexpr.utils] NumExpr defaulting to 2 threads.
[I 03:38:41 root] Patching zope.interface.document.asReStructuredText to format code
psi_0/sqrt(mu/g) = 0.19999999812364935
Error in callback <function flush_figures at 0x714d04214040> (for post_execute), with arguments args (),kwargs {}:
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib_inline/backend_inline.py:126, in flush_figures()
123 if InlineBackend.instance().close_figures:
124 # ignore the tracking, just draw and close all figures
125 try:
--> 126 return show(True)
127 except Exception as e:
128 # safely show traceback if in IPython, else raise
129 ip = get_ipython()
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib_inline/backend_inline.py:90, in show(close, block)
88 try:
89 for figure_manager in Gcf.get_all_fig_managers():
---> 90 display(
91 figure_manager.canvas.figure,
92 metadata=_fetch_figure_metadata(figure_manager.canvas.figure),
93 )
94 finally:
95 show._to_draw = []
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/matplotlib/backend_bases.py:2186, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2182 try:
2183 # _get_renderer may change the figure dpi (as vector formats
2184 # force the figure dpi to 72), so we need to set it again here.
2185 with cbook._setattr_cm(self.figure, dpi=dpi):
-> 2186 result = print_method(
2187 filename,
2188 facecolor=facecolor,
2189 edgecolor=edgecolor,
2190 orientation=orientation,
2191 bbox_inches_restore=_bbox_inches_restore,
2192 **kwargs)
2193 finally:
2194 if bbox_inches and restore_bbox:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/backend_bases.py:2042, in FigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>(*args, **kwargs)
2038 optional_kws = { # Passed by print_figure for other renderers.
2039 "dpi", "facecolor", "edgecolor", "orientation",
2040 "bbox_inches_restore"}
2041 skip = optional_kws - {*inspect.signature(meth).parameters}
-> 2042 print_method = functools.wraps(meth)(lambda *args, **kwargs: meth(
2043 *args, **{k: v for k, v in kwargs.items() if k not in skip}))
2044 else: # Let third-parties do as they see fit.
2045 print_method = meth
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/backends/backend_agg.py:481, in FigureCanvasAgg.print_png(self, filename_or_obj, metadata, pil_kwargs)
434 def print_png(self, filename_or_obj, *, metadata=None, pil_kwargs=None):
435 """
436 Write the figure to a PNG file.
437
(...) 479 *metadata*, including the default 'Software' key.
480 """
--> 481 self._print_pil(filename_or_obj, "png", pil_kwargs, metadata)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/backends/backend_agg.py:429, in FigureCanvasAgg._print_pil(self, filename_or_obj, fmt, pil_kwargs, metadata)
424 def _print_pil(self, filename_or_obj, fmt, pil_kwargs, metadata=None):
425 """
426 Draw the canvas, then save it using `.image.imsave` (to which
427 *pil_kwargs* and *metadata* are forwarded).
428 """
--> 429 FigureCanvasAgg.draw(self)
430 mpl.image.imsave(
431 filename_or_obj, self.buffer_rgba(), format=fmt, origin="upper",
432 dpi=self.figure.dpi, metadata=metadata, pil_kwargs=pil_kwargs)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/backends/backend_agg.py:382, in FigureCanvasAgg.draw(self)
379 # Acquire a lock on the shared font cache.
380 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
381 else nullcontext()):
--> 382 self.figure.draw(self.renderer)
383 # A GUI class may be need to update a window using this draw, so
384 # don't forget to call the superclass.
385 super().draw()
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
96 renderer.stop_rasterizing()
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/figure.py:3264, in Figure.draw(self, renderer)
3261 # ValueError can occur when resizing a window.
3263 self.patch.draw(renderer)
-> 3264 mimage._draw_list_compositing_images(
3265 renderer, self, artists, self.suppressComposite)
3267 renderer.close_group('figure')
3268 finally:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3223 if artists_rasterized:
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
3230 self.stale = False
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
132 if not_composite or not has_images:
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
137 image_group = []
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/axis.py:1408, in Axis.draw(self, renderer)
1405 tlb1, tlb2 = self._get_ticklabel_bboxes(ticks_to_draw, renderer)
1407 for tick in ticks_to_draw:
-> 1408 tick.draw(renderer)
1410 # Shift label away from axes to avoid overlapping ticklabels.
1411 self._update_label_position(renderer)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/axis.py:276, in Tick.draw(self, renderer)
273 renderer.open_group(self.__name__, gid=self.get_gid())
274 for artist in [self.gridline, self.tick1line, self.tick2line,
275 self.label1, self.label2]:
--> 276 artist.draw(renderer)
277 renderer.close_group(self.__name__)
278 self.stale = False
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
68 if artist.get_agg_filter() is not None:
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
73 if artist.get_agg_filter() is not None:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/text.py:778, in Text.draw(self, renderer)
775 self.update_bbox_position_size(renderer)
776 self._bbox_patch.draw(renderer)
--> 778 gc = renderer.new_gc()
779 gc.set_foreground(self.get_color())
780 gc.set_alpha(self.get_alpha())
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/backend_bases.py:606, in RendererBase.new_gc(self)
604 def new_gc(self):
605 """Return an instance of a `.GraphicsContextBase`."""
--> 606 return GraphicsContextBase()
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/backend_bases.py:693, in GraphicsContextBase.__init__(self)
691 self._rgb = (0.0, 0.0, 0.0, 1.0)
692 self._hatch = None
--> 693 self._hatch_color = colors.to_rgba(rcParams['hatch.color'])
694 self._hatch_linewidth = rcParams['hatch.linewidth']
695 self._url = None
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/__init__.py:802, in RcParams.__getitem__(self, key)
799 from matplotlib import pyplot as plt
800 plt.switch_backend(rcsetup._auto_backend_sentinel)
--> 802 return self._get(key)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/__init__.py:701, in RcParams._get(self, key)
682 """
683 Directly write data bypassing deprecation and validation logic.
684
(...) 697 :meta public:
698 """
699 dict.__setitem__(self, key, val)
--> 701 def _get(self, key):
702 """
703 Directly read data bypassing deprecation, backend and validation
704 logic.
(...) 718 :meta public:
719 """
720 return dict.__getitem__(self, key)
KeyboardInterrupt:
Velocity of a Travelling Wave#
In the previous discussion, \(v\) is the phase velocity of the traveling wave. This is the velocity with which the peaks of the wave travel. In contrast, the group velocity determines how fast a wave packet will move, which is a different concept. This velocity, however, is frame dependent, and in general, we would prefer to characterize this with respect to some specific frames.
For waves with infinite wavelength on a non-zero background density (such as dark solitons, and some bright solitons), a preferred frame is defined by the background at infinity which is taken to be at rest.
For infinite-wavelength bright solitons with no background density, no preferred frame exists. By boosting, one can realize these solutions with any velocity.
For small amplitude oscillations (phonons), a preferred reference frame is defined by a stationary background density.
The question is how to define the preferred frame for large-amplitude, finite-wavelength modes. One possibility is to define the frame with no net current:
\[
\bar{p} = \frac{1}{L}\int_0^L\d{x}\; m \rho(x,t)u(x,t)
= \frac{m}{L}\int_0^L\d{x_v} \left(v\rho_v(x_v) - C\right)
= mv\bar{\rho} - m C = 0.
\]
With this definition, we have phase velocity:
\[
v = \frac{C}{\bar{\rho}}.
\]
Hoefer and El#
We now consider the solution of Hoefer and El. They set \(\hbar=m=g=1\) so that we have the following units (in \(d=3\) dimensions):
\[\begin{split}
[\hbar] = MD^2/T, \qquad
[g] = MD^{2+d}/T^2, \qquad
[m] = M\\
M = m = 1, \qquad
mg/h^2 = D^{d-2}
D = \left(\frac{mg}{\hbar^2}\right)^{1/(d-2)} = 1, \qquad
T = \frac{m}{\hbar}\left(\frac{mg}{\hbar^2}\right)^{2/(d-2)} = 1,\\
M = m = 1, \qquad
D = \frac{mg}{\hbar^2} = 1, \qquad
T = \frac{m^3g^2}{\hbar^5} = 1,\\
\end{split}\]
\[
\I\dot{\psi} = -\frac{1}{2}\pdiff[2]{\psi}{x} + \abs{\psi}^2\psi,\qquad
\frac{\psi_V''}{2} = \left(\abs{\psi_V}^2 - \mu + \frac{V^2}{2}\right)\psi_V.
\]
Or, using the Madelung transform:
(25)#\[\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}\]
They express their solution in terms of four constants \(r_1\leq r_2\leq r_3\leq r_4\) which we can related to some more physical parameters such as the amplitude \(a\), minimum and maximum densities \(\rho_{\min,\max}\), period \(L\), and phase velocity \(V\);
(26)#\[\begin{align}
\DeclareMathOperator{\sn}{sn}
\rho_\min &= (r_1-r_2-r_3+r_4)^2, \\
V &= \frac{r_1+r_2+r_3+r_4}{2}, \tag{2.118}\\
a = \rho_\max - \rho_\min &= (r_2-r_1)(r_4-r_3), \tag{2.117}\\
\frac{1}{l^2} &= (r_4-r_2)(r_3-r_1),\tag{2.117}\\
m = al^2 &= \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)}, \tag{2.119}\\
\rho(x, t) &= \rho_\min + a \sn^2\left(\frac{x_v}{l}, m\right), \tag{2.117}\\
L &= 2lK(m), \tag{2.121}
\end{align}\]
We disagree with their definition of \(C\) which is:
\[
C = \frac{1}{8}(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)\overbrace{(r_1-r_2-r_3+r_4)}^{\sqrt{\rho_\min}}.
\]
Instead, we find
\[
C = \frac{\sqrt{\rho_\min\rho_\max(1+\rho_\min l^2)}}{l}
=\sqrt{(-r_1-r_2+r_3+r_4)^2 - 3a}
\sqrt{(-r_1+r_2-r_3+r_4)^2 - \frac{3}{l^2}}
(r_1 - r_2 - r_3 + r_4)
\]
In our code, we would like to specify \(\rho_\min\), \(\rho_\max\), the period \(L\) and the phase velocity \(V\). To do this, we must solve for the Jacobi parameter \(m\):
\[
L = 2lK(m) = 2\sqrt{\frac{m}{a}}K(m).
\]
import gpe.utils;reload(gpe.utils)
import gpe.bec;reload(gpe.bec)
from gpe.utils import evolve
import gpe.exact_solutions;reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves
# v_p is the velocity of the constructed wave.
# v_x is the velocity of the frame.
s = TravellingWaves(Nx=64*2, Lx=10.0, n0=1.0, n1=1.1, v_p=0.5, v_x=0.5)
t_max = 10.0
for y in evolve(s, steps=500, dt_t_scale=0.5, t_max=t_max):
plt.clf()
y.plot()
The current is:
\[
\vect{j} = \rho u = V\rho - C.
\]
In the moving frame:
\[
(\rho_V u_V)' = 0, \qquad
\mu + \frac{V^2}{2} = \rho_V + \frac{u_V^2}{2}
- \frac{1}{2}\frac{\sqrt{\rho_V}''}{\sqrt{\rho_V}}.
\]
Thus, if \(\rho_V(x_v)\) is specified, then \(u_V(x_V) = C/\rho_V(x_V)\) will ensure that the conservation equation is satisfied if:
\[
\mu + \frac{V^2}{2} = \rho_V + \frac{C^2}{2\rho_V^2}
- \frac{1}{2}\frac{\sqrt{\rho_V}''}{\sqrt{\rho_V}}.
\]
The solution of Hoefer and El can be can be expressed as (Note: \(k=\sqrt{m}\) for \(\sn(z,k)\) in some CASs. We use \(\sn(z;m)\) and \(K(m)\) here):
\[\begin{split}
\psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad
u(x, t) = \phi'(x, t),\\
m = al^2, \qquad
\rho(x, t) = \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)}, \qquad
\rho_{\max} = a + \rho_{\min}, \qquad
C = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}{l},\qquad
L = 2lK(m), \qquad
Q = \frac{2\pi}{L},\\
\bar{\rho} = \frac{1}{L}\int_0^{L}\d{x}\;\rho(x, t)
= \rho_{\min} + \frac{a}{m}
= \left(1 - \frac{1}{m}\right)\rho_{\min} + \frac{1}{m}\rho_\max
= \rho_{\min} + \frac{1}{l^2},\\
V_p = \frac{C}{\bar{\rho}}
= \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}
{l\rho_{\min} + \frac{1}{l}}
\end{split}\]
Note that, averaging over a complete period \(L = 2lK(m)\) is:
\[
\bar{\rho} = \frac{1}{L}\int_0^L \rho(x)\d{x}
= \rho_\min + a\frac{1}{L}\int_0^L \sn^2(\tfrac{x}{l}, k)\d{x}
= \rho_\min + \frac{a}{m}
= \left(1-\frac{1}{m}\right)\rho_\min + \frac{1}{m}\rho_\max.
\]
\[
V_p = \frac{C}{\bar{\rho}}
= \frac{\sqrt{1+\rho l^2}}{l + \frac{1}{l\rho}}
\]
import gpe.utils;reload(gpe.utils)
import gpe.bec;reload(gpe.bec)
from gpe.utils import evolve
import gpe.exact_solutions;reload(gpe.exact_solutions)
from gpe.exact_solutions import BrightSoliton
s = BrightSoliton(Nx=64, Lx=10.0,
v=0.0, # v=10.0
v_x=10.0)
t_max = s.basis.Lxyz[0]/abs(s.v_x)
for y in evolve(s, steps=500, t_max=t_max):
plt.clf()
y.plot()
Reformulation#
Here we present a reformulation of the coefficients into more physically meaningful parameters. We start with the
\[
x_v = x - vt, \qquad
u(x_v) = v - \frac{C}{\rho(x_v)}.
\]
This form immediately satisfies the continuity equation as long as \(\rho(x_v) = \rho(x-vt)\):
\[
\rho_{,t} + (u\rho)_{,x} = -v\rho'(x_v) + (v\rho(x_v) - C)' = 0
\]
Next we have:
\[
\rho(x_v) = \rho_\min + a\sn^2\left(\frac{x_v}{l}; m\right).
\]
%pylab inline --no-import-all
from gpe.exact_solutions import sn, K
m = 0.999999999
L = 2*K(m)
x = np.linspace(-L,L,100)
plt.plot(x, sn(x, m)**2)
import numpy as np
np.random.seed(1)
r = np.random.random(5)
r[4] = -r[1]-r[2]-r[3]
C = 1./8*(-r[1]-r[2]+r[3]+r[4])*(-r[1]+r[2]-r[3]+r[4])*(r[1]-r[2]-r[3]+r[4])
V = (r[1]+r[2]+r[3]+r[4])/2
m = (r[2]-r[1])*(r[4]-r[3])/(r[4]-r[2])/(r[3]-r[1])
k = np.sqrt(m)
a = (r[4]-r[3])*(r[2]-r[1])
l = 1./np.sqrt((r[4]-r[2])*(r[3]-r[1]))
rho_min = (r[1]-r[2]-r[3]+r[4])**2
rho_max = rho_min + a
V, np.sqrt(rho_min*rho_max*(1+rho_min*l**2))/l, C
Phonon Limit#
In the limit of small amplitude \(a\) we have the following reduction:
\[ \begin{align}\begin{aligned}
a\rightarrow 0, \qquad
\bar{\rho} \rightarrow \rho, \qquad
C \rightarrow \frac{\rho}{l}\sqrt{1+\rho l^2}}
\epsilon = \rho_\max - \rho_\min = \epsil we obtain\\In the phonon limit we should have\end{aligned}\end{align} \]
\omega(q) = \pm\sqrt{q^2(\rho_0 + q^2/4)}, \qquad V_p = \frac{\omega(q)}{q} = \pm\sqrt{\rho_0 + q^2/4}, \qquad V_g = \omega’(q) = \frac{q^2+2\rho_0}{\pm\sqrt{q^2 + 4\rho_0}} $$
Thus, one can specify the following four parameters (and an overall arbitrary initial phase):
\(\rho_\min\) and \(\rho_\max\): These determine the amplitude \(a = \rho_\max - \rho_\min \geq 0\) of the wave.
\(L\): The period of the wave can be chosen implicitly by setting \(l\) which will determine the period through \(m=al^2\) and \(L = 2l K(m)\).
\(V\): The speed of the wave may also be arbitrarily chosen. This may seem strange since the speed of the waves should follow from the dispersion relationship, however, this is simply a reflection of the fact that we can look at the system in an arbitrarily moving frame.
The rest frame is defined as the frame in which there is no net momentum flow. Thus, the average current should be zero:
\[\begin{split} 0 = \bar{\vect{j}} = \frac{1}{L}\int_0^L \vect{j}(x, t)\d{x} = \frac{1}{L}\int_0^L \rho(x, t)u(x, t)\d{x} = \frac{1}{L}\int_0^L \Bigl(V\rho(x, t) - C\Bigr)\d{x}\\ V_p = \frac{CL}{\int_0^L \rho(x, t)\d{x}} = \frac{C}{\bar{\rho}} \end{split}\]where \(\bar{\rho}\) is the average density across the cell. Note: this depiction will fail for bright solitons on the vacuum since \(\bar{\rho} = 0\) but in this case, the notion of the fixed background has no meaning.
Soliton Limit#
In the long wavelength limit \(m \rightarrow 1\):
\[
K(m) \rightarrow \ln \frac{4}{\sqrt{1-m}}, \qquad
l = \frac{1}{\sqrt{a}}, \qquad
C = \rho_\max\sqrt{\rho_{\min}},\qquad
V_p = \frac{C}{\rho_\max} = \sqrt{\rho_{\min}}.
\]
Thus, the speed of the traveling waves is set by the minimum density in the long wavelength limit. This is consistent, for example, with the speed of a grey soliton which becomes stationary when it becomes dark (i.e. when the minimum density in the soliton becomes zero.)
\[\begin{split}
\hbar\omega(q) = \pm \sqrt{\frac{\hbar^2q^2}{2m}\left(\frac{\hbar^2q^2}{2m} + 2\mu\right)}\\
\omega(q) = \pm \abs{q}\sqrt{\rho_0 + \frac{q^2}{4}}\\
\end{split}\]
\[ L = 2\sqrt{\frac{m}{a}}K(m) \]
\[\begin{split}
\partial_t[e^{\I\phi}f(x-vt)] = \I\dot{\phi}e^{\I\phi}f(x-vt) -ve^{\I\phi}f'(x-vt)\\
\partial_t[\psi] = \I\partial_t\phi \psi - v\psi'(x-vt)
\end{split}\]
from gpe.imports import *
import gpe.exact_solutions; reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves
s = TravellingWaves(v_p=10.0, v_x=10.0, n0=0.1, Nx=32*4, Lx=5)
s.plot()
t_max = s.basis.Lxyz[0]/abs(s.v_p)
for y in evolve(s, steps=400, t_max=t_max):
plt.clf()
y.plot()
from gpe.imports import *
#from mmfutils.math.special import ellipkinv
from scipy.special import ellipj, ellipk
import scipy.integrate
import scipy as sp
import gpe.bec
def sn(u, m):
return ellipj(u, m)[0]
def K(m):
return ellipk(m)
def get_QV(m=0.1, a=0.1, rho_min=1.0, N=256):
rho_max = rho_min + a
#m = a*l**2
l = np.sqrt(m/a)
L = 2*l*K(m)
dx = L/N
x = np.arange(N)*dx
rho = rho_min + a*sn(x/l, m)**2
C = np.sqrt(rho_min*rho_max*(1+rho_min*l**2))/l
V = np.trapz(C/rho, x)/L
Q = 2*np.pi/L
return Q, V
rho_min = 1.0
c = np.sqrt(rho_min)
ms = 1 - 10**np.linspace(-10, -0.001, 1000)
plt.figure(figsize=(20,5))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
for a in [0.01, 0.2, 1.0]:
Q, Vp = np.array([get_QV(m, a=a, rho_min=rho_min) for m in ms]).T
l, = ax1.plot(Q**2, Vp/c, '-+', label=a)
Q_ = (Q[1:]+Q[:-1])/2
Vp_ = (Vp[1:]+Vp[:-1])/2
Vg_ = Q_*np.diff(Vp)/np.diff(Q) +Vp_
ax1.plot(Q_**2, Vg_/c, '--', c=l.get_c())
ax2.plot(Q, Q*Vp, '-+', label=a)
q_max = 2.0
q = np.sqrt(np.linspace(0, q_max**2, 100))
Vp = np.sqrt(rho_min + q**2/4)
w = Vp*q
Vg = (q**2 +2*rho_min)/np.sqrt(q**2 +4*rho_min)
l, = ax1.plot(q**2, Vp/c, label='phonon')
ax1.plot(q**2, Vg/c, '--', c=l.get_c())
ax1.set_xlim(0, q_max**2)
ax1.set_ylim(0.99, 1.5)
plt.sca(ax1)
plt.axis([0,1.5,1,1.5])
ax1.set_xlabel("$q^2$")
ax1.legend()
plt.sca(ax2)
ax2.plot(q, w, label='phonon')
ax2.set_xlabel("$q$")
ax2.set_ylabel("$\omega(q)$")
ax2.set_xlim(0, 1.25)
ax2.set_ylim(0, 2)
#plt.axis([0,0.5,0,0.5])
ax2.legend()
get_QV(m=0.9999999999999999, a=1.0)
1.15*2.0
from gpe.soc import u
L = 2.5/3.0*u.micron
g = 0.00077
m = 87*u.u
hbar = u.hbar
u_length = g/hbar**2*m
Q = 2*np.pi/(L/u_length)
rho_max = 10530.0/u.micron**3*u_length**3
rho_min = 452.0/u.micron**3*u_length**3
a = rho_max - rho_min
get_QV(m=0.999999992, a=a)
from gpe.imports import *
import scipy.optimize
from scipy.special import ellipj, ellipk
import scipy.integrate
import scipy as sp
import gpe.bec
def sn(u, m):
return ellipj(u, m)[0]
def cn(u, m):
return ellipj(u, m)[0]
def K(m):
return ellipk(m)
def get_m(a, L):
def f(m):
return 2*np.sqrt(m/a)*K(m) - L
m0 = m1 = 0
f0 = f1 = f(m0)
while f1 < 0:
m1 = (m1 + 1)/2.0
f1 = f(m1)
return sp.optimize.brentq(f, m0, m1)
class State(gpe.bec.State):
def __init__(self, Nx=256, cells_x=10,
l=1.1757305, b=1.0, a=0.1, V=0.0):
m = a*l**2
L = 2*l*K(m) * cells_x
mu = ((a+3*b)*l**2+1)/2/l**2
C = np.sqrt((a+b)*b*(1+b*l**2)/l**2)
dx = L/Nx
gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(L,),
m=1.0, g=1.0, hbar=1.0)
x = self.basis.xyz[0]
rho = b + a*sn(x/l, m)**2
u = V - C/rho
psi = np.sqrt(rho)*np.exp(1j*sp.integrate.cumtrapz(u, x, initial=0))
self.data[...] = psi
# Hpsi = np.fft.ifft(k**2*np.fft.fft(psi))/2 + (abs(psi)**2-mu)*psi
#plt.plot(x, abs(Hpsi))
def get_Vext(self):
return 0.0
s = State(Nx=256)
e = EvolverABM(s, dt=0.5*s.t_scale)
with NoInterrupt() as interrupted:
while not interrupted:
e.evolve(100)
e.y.plot()
display(plt.gcf())
plt.close('all')
clear_output(wait=True)
import scipy.integrate
import scipy as sp
l = 1.1858
l = 1.1757305
b = 1.0
a = 0.1
V = 0.0
m = a*l**2
L = 2*l*K(m)
mu = ((a+3*b)*l**2+1)/2/l**2
C = np.sqrt((a+b)*b*(1+b*l**2)/l**2)
N = 256
dx = 10*L/N
x = np.arange(N)*dx - L/2
k = 2*np.pi * np.fft.fftfreq(N, dx)
rho = b + a*sn(x/l, m)**2
u = V - C/rho
#plt.plot(x, rho)
psi = np.sqrt(rho)*np.exp(1j*sp.integrate.cumtrapz(u, x, initial=0))
Hpsi = np.fft.ifft(k**2*np.fft.fft(psi))/2 + (abs(psi)**2-mu)*psi
plt.plot(x, abs(Hpsi))
plt.plot(x, rho*u**2 + rho**2/2)
plt.plot(x[1:-1], 2.36+rho[1:-1]*np.diff(np.diff(np.log(rho)))/dx**2/4)
from gpe.imports import *
from gpe.bec import State
s =
#
Phonons#
Phonons appear in the small amplitude limit \(a\rightarrow 0\) in which case:
\[\begin{split}
\sn(z;0) = \sin(z),\qquad
K(m=0) = \frac{\pi}{2}, \qquad
L = \pi l, \qquad
Q = \frac{2}{l}\\
\rho(x, t) = \rho_{\min} + a \sin^2\frac{Q x_v}{2}
= \rho_{0} - \frac{a}{2}\cos(Q x_v), \qquad
\end{split}\]
\[\begin{split}
C = \rho_0\sqrt{l^{-2}+\rho_{0}}\\
u = V - \sqrt{\frac{Q^2}{4}+\rho_{0}}
\end{split}\]
The phase velocity is \(V = \omega(Q)/Q\)
and present the density as:
\[\begin{split}
\rho = \frac{1}{4}(r_4-r_3-r_2+r_1)^2
+ u^2\sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad
x_v = x - vt\\
v = \frac{n_1+n_2+n_3+n_4}{2}, \qquad
a = u^2 = (r_4-r_3)(r_2-r_1), \\
m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)}
= u^2l^2, \qquad
l = \frac{1}{\sqrt{(r_4-r_2)(r_3-r_1)}},\\
c^2 = .
\end{split}\]
The soliton limit is \(m=1\) (i.e. \(r_2=r_3\)) where \(\sn(z, 1) = \tanh(z)\) where the solution has the form:
\[\begin{split}
\psi = \left[\I v + u \tanh\bigl(\frac{x_v}{l}\bigr)\right]e^{-\I (u^2+v^2) t}, \qquad
u = 1/l,\\
\rho = v^2 + u^2 \tanh^2\bigl(\frac{x_v}{l}\bigr)
= v^2 + u^2 - u^2\sech^2\bigl(\frac{x_v}{l}\bigr)\\
v^2 = \frac{1}{4}(r_4-2r_2+r_1)^2\\
v^2 + u^2 = c^2 = \bar{\rho} = \frac{(r_4-r_1)^2}{4}.
\end{split}\]
\[
(r_4-r_2)(r_2-r_1) + \frac{(r_4 - 2r_2 + r_1)^2}{4} = \frac{(r_1-r_4)^2}{4}
\]
A traveling wave with velocity \(v\) has the following form:
\[\begin{split}
\psi(x, t) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt) - \I\omega t}, \\
\I\hbar \dot{\psi} = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.\\
-\I\hbar v\psi' + \hbar\omega\psi = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.
\end{split}\]
By appropriately adjusting \(\mu\) we may set \(\omega = 0\), hence we may drop the factor \(e^{-\I\omega t}\) from the wavefunction:
\[
\psi(x - vt) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt)}.
\]
The last equation above now represents a stationary solution in a moving frame with coordinate \(y = x-vt\):
\[
0 = -\frac{\hbar^2\psi''(y)}{2m} + \I\hbar v\psi'(y) + [g\rho(y)- \mu]\psi(y)
= \left[\frac{\op{p}^2}{2m} -\op{p} v + g\rho(y) - \mu\right]\psi(y)
= \left[\frac{(\op{p} - p_v)^2}{2m} - \frac{p_v^2}{2m} + g\rho(y) - \mu\right]\psi(y)
\]
where \(p_v = mv\). Finally, we may recast this in terms of the original GPE by transforming \(\psi(y)\) as follows:
\[\begin{split}
\psi(y) = e^{\I p_v x/\hbar}\tilde{\psi}(y)\\
(\op{p}-p_v)\psi(y) = e^{-\I p_v x/\hbar}\op{p}\tilde{\psi}(y)\\
\left[\frac{\op{p}^2}{2m} + g\rho(y) - \left(\mu + \frac{p_v^2}{2m}\right)\right]\tilde{\psi}(y)
\end{split}\]
If the original function \(\psi(y+L) = \psi(y)\) was periodic, then \(\tilde{\psi}(y+L) = e^{\I p_v L/\hbar}\tilde{\psi}(y)\). In other words, twisted boundary conditions should be applied with a twist \(\theta = p_vL/\hbar\).
GPE#
\[\begin{split}
\psi(x-vt) = e^{\I\phi(x-vt)}f(x-vt), \qquad f(x) = \sqrt{\rho(x)}\\
\psi' = \I\phi'\psi + e^{\I\phi}f', \qquad
\psi'' = \I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi'
+ e^{\I\phi}\I\phi'f' + e^{\I\phi}f'',\\
\I\dot{\psi} = -v\I\psi' = v\phi'\psi - v\I e^{\I\phi}f'
\end{split}\]
\[\begin{split}
2v\phi'\psi - 2v\I e^{\I\phi}f' =
\I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi' + e^{\I\phi}\I\phi'f' + e^{\I\phi}f'' + 2e^{\I\phi}f^3\\
[3\phi'+2v]f' = -\phi''f\\
[3\phi' + 2v]\phi'f = f'' + 2f^3\\
-\left(\frac{(\phi')^2}{2}\right)'
= \frac{f'f''}{f^2} + (f^2)'
= \frac{f'f''}{f^2} + \rho'\\
\end{split}\]
\[
f^2 = \rho, \qquad
2ff' = \rho', \qquad
2f'f' + 2ff'' = \rho'', \qquad
\]
Here we consider the problem of finding traveling wave solutions in a BEC. From a numerical perspective, it is highly beneficial if the problem can be stated in terms of a well-defined minimization problem. To start, we consider the available analytic solution for the conventional GPE. These solutions are presented in [El:2016]:
\[\begin{split}
\DeclareMathOperator{\sn}{sn}
\I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
\rho = \abs{\psi}^2 = \frac{(r_4-r_3-r_2+r_1)^2}{4}
+ (r_4-r_3)(r_2-r_1)\sn^2(\sqrt{(r_4-r_2)(r_3-r_1)}\xi; m)\\
C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
u = v - \frac{C}{\rho}, \qquad
\xi = x-vt - \xi_0, \qquad
v = \frac{1}{2}\sum r_i, \\
a = (r_4-r_3)(r_2-r_1), \qquad
r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
L = \frac{1}{2}\oint \frac{\d\lambda}
{\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
= \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
\end{split}\]
The physical interpretations are:
\(a\): Amplitude
\(v\): Phase velocity
\(u\):
\[\begin{split}
\DeclareMathOperator{\sn}{sn}
\I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
\rho = \abs{\psi}^2 = \alpha + a\sn^2(z; m)\\
\alpha = \frac{(r_4-r_3-r_2+r_1)^2}{4}, \\
D = \sqrt{(r_4-r_2)(r_3-r_1)}, \qquad
z = D\xi\\
C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
u = v - \frac{C}{\rho}, \qquad
\xi = x-vt - \xi_0, \qquad
v = \frac{1}{2}\sum r_i, \\
a = (r_4-r_3)(r_2-r_1), \qquad
r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
L = \frac{1}{2}\oint \frac{\d\lambda}
{\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
= \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
\end{split}\]
\[\begin{split}
\DeclareMathOperator{\cn}{cn}
\DeclareMathOperator{\dn}{dn}
\sn'(z) = \cn(z)\dn(z), \qquad
\cn'(z) = -\sn(z)\dn(z), \qquad
\dn'(z) = -k^2\sn(z)\cn(z)\\
\sn''(z) = -\sn(z)[\dn^2(z) +k^2\cn^2(z)], \qquad
\cn'(z) = -\sn(z)\dn(z), \qquad
\dn'(z) = -k^2\sn(z)\cn(z)\\
\cn^2 + \sn^2 = \dn^2 + k^2\sn^2 = 1.
\end{split}\]
A special limit is when \(m=1\):
\[
\sn(z;m=1) = \tanh(z), \qquad
\cn(z;m=1) = \dn(z;m=1) = \sech(z).
\]
\[\begin{split}
\psi = \sqrt{1+\alpha\sn^2(z)}\\
\psi' = \frac{\alpha\sn(z)\sn'(z)}{\psi}\\
\psi'' = \frac{\alpha}{\psi}\left(
\sn'(z)\sn'(z) + \sn(z)\sn''(z)
-\frac{\alpha\sn^2(z)[\sn'(z)]^2}{\psi^2}\right)\\
\end{split}\]
\[\begin{split}
\psi' = \frac{\alpha\sn(z)\cn(z)\dn(z)}{\psi}\\
\psi'' = \frac{\alpha}{\psi}\left(
\cn^2(z)\dn^2(z) - \sn^2(z)\dn^2(z) - k^2\sn^2(z)\cn^2(z)
-\frac{\alpha\sn^2(z)\cn^2(z)\dn^2(z)}{\psi^2}\right)\\
\end{split}\]
First we test these. To match units we set \(\hbar = m = g = 1\).
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec
def sn(u, m):
return ellipj(u, m)[0]
def cn(u, m):
return ellipj(u, m)[0]
def K(m):
return ellipk(m)
class State(gpe.bec.State):
def __init__(self, Nx=32, rs=[0.1, 0.2, 0.3, 0.4], xi0=0):
r1, r2, r3, r4 = rs = sorted(rs)
k = np.sqrt((r4-r2)*(r3-r1))
m = (r2-r1)*(r4-r3)/(r4-r2)/(r3-r1)
v = sum(rs)/2.0
L = 2.*K(m)/k
gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(L,),
m=1.0, g=1.0, hbar=1.0)
x = self.xyz[0]
self.T = abs(L/v)
t = 0
xi = x - v*t - xi0
a = (r4-r3)*(r2-r1)
#a = m**2*k**2
rho_0 = (r4-r3-r2+r1)**2/4.0
rho_0 = (r4-r3+r2-r1)**2/4.0
rho = rho_0 + a*sn(k*xi, m)**2
self.a_ = self.a = a
self.m_ = m
self.k_ = self.k = k
self.rho_0_ = self.rho_0 = rho_0
k_ = self.m*v/self.hbar
self[...] = np.exp(1j*k_*x) * np.sqrt(rho)
def get_Vext(self):
return 0.0
s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])
n = s.get_density()
x = s.xyz[0]
x_ = x[1:-1]
psi_ = s[1:-1]
ddpsi_ = np.diff(np.diff(s[...]))/np.diff(x)[1:]**2
plt.plot(x_, ddpsi_/psi_/2 + abs(psi_)**2)
Check the soliton limit (2.122)
s = State(Nx=128, rs=[-4.000001, 1.0, 1.000001, 2.0])
x = s.xyz[0]
s.plot()
rho_min = abs(s[...]).min()**2
rho_max = abs(s[...]).max()**2
plt.axhspan(rho_min, rho_min + s.a, fc='y', alpha=0.5)
#plt.plot(x, rho_max-s.a/np.cosh(np.sqrt(s.a)*x)**2, '+:')
rho_0 = s.rho_0
s[...] = np.sqrt(rho_0)*np.tanh(np.sqrt(rho_0)*x)
s.plot()
s.cooling_phase = 1
e = EvolverABM(s, dt=0.9*s.t_scale)
with NoInterrupt(ignore=True) as interrupted:
while e.t < e.y.T and not interrupted:
e.evolve(100)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
Solitons#
The GPE admits grey solitons with infinite period. These solitons arise from the previous solutions in the limit where \(m\rightarrow 1\) so that \(\sn\rightarrow\tanh\) and \(\cn\rightarrow\dn\rightarrow \sech\).
The soliton solution is:
\[\begin{split}
\psi(x, t) = \sqrt{\bar{\rho}}e^{-\I\mu t/\hbar}\left[
\I\frac{v}{c}
+ \frac{u}{c}\tanh\left(\frac{x-vt}{l}\right)\right], \qquad
\rho = \bar{\rho}\left[\frac{v^2}{c^2}
+ \frac{u^2}{c^2}\tanh^2\left(\frac{x-vt}{l}\right)\right]\\
u^2+v^2 = c^2, \qquad \mu=g\bar{\rho} = mc^2,\qquad l=\hbar/m u.
\end{split}\]
This can be expressed as
\[
\psi(x, t) = \sqrt{\rho}e^{-\I\mu t/\hbar + \I\phi(x)}, \qquad
\tan\phi(x) = \frac{v}{u\tanh\left(\frac{x-vt}{l}\right)}, \qquad
\phi_{\mathrm{twist}} = \phi(\infty) - \phi(-\infty) = 2\tan^{-1}\frac{v}{u}
\]
\[
\rho(0) = \bar{\rho}\frac{v^2}{c^2}
\]
or, in the same units \(m=\hbar = g = 1\) as El and Hoefer:
\[\begin{split}
\rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2
\end{split}\]
To compare with (2.117):
\[
\rho = \frac{(r_4-r_3-r_2+r_1)^2}{4} +(r_4-r_3)(r_2-r_1)\sn^2[\sqrt{(r_4-r_2)(r_3-r_1)}(x-Vt), m]
\]
we have \(m=1\), \(\sn = \tanh\) and
\[\begin{split}
\frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)} = m = 1,\\
\frac{(r_4-r_3-r_2+r_1)^2}{4} = \bar{\rho}\frac{v^2}{c^2} = v^2\\
a = (r_4-r_3)(r_2-r_1) = \bar{\rho}\left(1-\frac{v^2}{c^2}\right) = \bar{\rho} - v^2\qquad
\sqrt{(r_4-r_2)(r_3-r_1)} = u, \\
V = \frac{r_1+r_2+r_3+r_4}{2}
\end{split}\]
Everything here is reasonable except the definition of the phase velocity \(V\) which differs from the soliton velocity \(v\).
The background density is
\[
\bar{\rho} = \frac{(r_4-r_3-r_2+r_1)^2}{4} + (r_4-r_3)(r_2-r_1)
= \frac{(r_4-r_3+r_2-r_1)^2}{4}.
\]
To compare with (2.122):
, the stationary solution must have \(\bar{\rho} = a_s\) or:
\[
(r_4-r_1)^2 = 4(r_4-r_2)(r_2-r_1)
\]
If \(r_1 = -r_4-2r_2\), then
\[
4(r_4+r_2)^2 = 4(r_4-r_2)(r_4+3r_2)
\]
Issues#
I was having some issues finding solutions so I considered stationary solutions:
\[
\psi = \sqrt{\rho_0 + a\sn^2(kx,m)}.
\]
Symbolically solving for solutions to the GPE (in Maple) gives the following solutions:
(27)#\[\begin{align}
a&=k^2m^2, &&& \mu &= \rho_0 - \frac{k^4m^2}{2\rho_0}, \\
a&=0, &&& \mu &= \rho_0, & \psi &= \sqrt{\rho_0}.\tag{Constant}\\
a&=\rho_0 m^2 , & k&=\sqrt{\rho_0}, & \mu &= \rho_0- \frac{\rho_0m^2}{2}, & \psi &= \sqrt{\rho_0}\sqrt{1+m^2\sn^2(\sqrt{\rho_0}x,m)}.\tag{Constant}\\
a&=\rho_0 , & k&=\frac{\sqrt{\rho_0}}{m}, & \mu &= \rho_0 - \frac{\rho_0}{2m^2}, & \psi &= \sqrt{\rho_0}\sqrt{1+\sn^2\left(\frac{\sqrt{\rho_0}}{m}x,m\right)}.\tag{Constant}\\
\end{align}\]
(assuming \(a,k\neq 0\) and real \(k\)) gives:
\[
m^2 = \frac{a}{k^2}, \qquad
\mu = \rho_0 - \frac{ak^2}{2\rho_0}, \qquad
\rho = \rho_0[1 - \sn^2(kx,m)].
\]
restart;
X_:=JacobiSN(k*x,sqrt(m2));
psi:=sqrt(rho[0] + a*X_^2);
Rho:=psi^2:
res:=collect(numer(simplify(-diff(psi, x, x)/2/psi + Rho - mu)), X_):
mu:=expand(solve(coeff(res, X_, 0), mu));
m2:=solve(coeff(res, X_^6), m2);
solve([coeffs(res, X_)], [a,k^2]);
The dark soliton solution has \(m=1\), \(a=k^2\)
s.m
\[\begin{split}
a = (r_4-r_3)(r_2-r_1)\\
a = mk^2\\
k^2 = (r_4-r_2)(r_3-r_1)\\
\end{split}\]
s = State(rs=[-4.1, 1.0, 1.1, 2.])
x = s.xyz[0]
m = s.m_
a = s.a_
k = s.k_
rho_0 = s.rho_0_
print(m, a, m**2*k**2, rho_0, k)
k = np.sqrt(rho_0)
plt.plot(x, 1+sn(k*x,m)**2)
#s[...] = np.sqrt(rho_0*(1-sn(x,m)))
\[
m^2 = \frac{a\alpha^2}{D^2}, \qquad
(3m^2 - a) = 6\frac{a\alpha^2}{D^2}, \qquad
-2a(m^2 + 1) = 6\frac{a\alpha^2}{D^2}, \qquad
\mu = \frac{-D^2a}{2} + \alpha^2
\]
\[\begin{split}
-2a(m^2 + 1) = (3m^2 - a) = 6m^2\\
a = -3m^2\\
m^2 + 1 = 1
\end{split}\]
To do this, we first consider boosting to a moving frame with velocity \(v\) so that in this frame, the traveling wave solution is stationary and periodic. We second Here we are looking for solutions that satisfy:
\[
\psi(x + L) = e^{\I m v L / \hbar}\psi(x)
\]
Minimization#
Here we consider the hypothesis that traveling waves can be found as minimum energy solutions in a box of period \(L\) with twisted boundary conditions holding both the total particle number fixed and the value of the wavefunction at one point.
We start with the soliton solution:
\[\begin{split}
\rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2.
\end{split}\]
At the core, the density is:
\[
\bar{\rho}\frac{v^2}{c^2}
\]
\[\begin{split}
\psi(x, t) = \sqrt{n_0}\left[\frac{v}{c}\I + \sqrt{1-\frac{v^2}{c^2}}\tanh\left(\frac{x-vt}{l}\right)\right]e^{-\I\mu t/\hbar}\\
\mu = gn = mc^2\\
\psi_v(x_v) = \sqrt{n_0}e^{-\I mvx_v/\hbar}\left[\frac{v}{c}\I + \sqrt{1-\frac{v^2}{c^2}}\tanh\left(\frac{x_v}{l}\right)\right]e^{-\I(\mu + m v^2/2) t/\hbar}\\
\mu = gn = mc^2
\end{split}\]
from gpe.imports import *
import gpe.Examples.traveling_waves;reload(gpe.Examples.traveling_waves)
from gpe.Examples.traveling_waves import (
StateTravellingWave, MinimizeStateTravellingWave)
from gpe.exact_solutions import TravellingWaves
n1 = 1.0
n0 = 0.5
tw0 = TravellingWaves(n1=n1, n0=0.5)
n_avg = tw0.get_density().mean()
s = StateTravellingWave(Nx=tw0.basis.Nx, Lx=tw0.basis.Lx,
psi_0=tw0[0], ind=0,
n_avg=n_avg, v_p=tw0.v_x, twist=tw0.twist)
s.set_psi(tw0.get_psi())
#s[...] = ts0[...]
s.cooling_phase = 1
e = EvolverABM(s, dt=0.1*s.t_scale, normalize=False)
with NoInterrupt() as interrupted:
while not interrupted:
e.evolve(200)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
plt.plot(tw0.xyz[0], s[...] - tw0[...])
plt.plot(tw0.get_Vext()- s.get_Vext())
dy0 = tw0.empty()
tw0.compute_dy_dt(dy=dy0, subtract_mu=False)
dy = s.empty()
s.compute_dy_dt(dy=dy, subtract_mu=False)
dy0[...] - dy[...]
s.plot()
dy = s.empty()
s.compute_dy_dt(dy=dy)
s.cooling_phase = 1j
e = EvolverABM(s, dt=0.5*s.t_scale, normalize=True)
with NoInterrupt() as interrupted:
while not interrupted:
e.evolve(100)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
from scipy.optimize import brentq
n1 = 1.0
n_avg = 0.9
def f(n0):
print(n0)
tw = TravellingWaves(n1=n1, n0=n0)
return tw.get_density().mean() - n_avg
brentq(f, 0, n1*0.999)
from gpe.imports import *
import gpe.Examples.traveling_waves;reload(gpe.Examples.traveling_waves)
from gpe.Examples.traveling_waves import StateTravellingWave, MinimizeStateTravellingWave
vs = np.linspace(-1,1,10)
Es = []
for v in vs:
s0 = StateTravellingWave(Nx=128, L=10.0, v_p=v, twist=0.1)
m = MinimizeStateTravellingWave(s0)
s = m.minimize()
Es.append(s.get_energy())
plt.clf()
s.plot()
display(plt.gcf())
clear_output(wait=True)
plt.plot(vs, Es)
Here is a stationary dark soliton.
%%time
dark_soliton = Soliton()
s0 = StateTravellingWave(psi_0=0, twist=np.pi, v_p=0)
s0[...] = dark_soliton.get_psi(s0.xyz[0])
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
plt.plot(x, abs(dark_soliton.get_psi(x))**2, '--', label='exact')
plt.legend(loc='best')
grey_soliton = Soliton(v=0.999)
plt.plot(np.angle(grey_soliton.get_psi(s0.xyz[0])))
v_c = 0.1
np.arctan2(v_c, np.sqrt(1-v_c**2))
v_c = grey_soliton.v/grey_soliton.c
theta = np.arctan2(v_c, np.sqrt(1-v_c**2))
s0 = StateTravellingWave(psi_0=0, twist=np.pi+2*theta, v_p=grey_soliton.v)
s0[...] = grey_soliton.get_psi(s0.xyz[0])
s0.psi_0 = abs(s0[...]).min()
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
plt.plot(x, abs(grey_soliton.get_psi(x))**2, '--', label='exact')
plt.legend(loc='best')
Here is a moving grey soliton.
plt.plot(np.angle(s._twist_phase_x)/np.pi)
%%time
rho_min = 0.01
v_p = np.sqrt(s.g*rho_min/s.m)
s0 = StateTravellingWave(psi_0=np.sqrt(rho_min), twist=np.pi, v_p=v_p)
s0[...] = s[...]
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
Exact Solution#
\[\begin{split}
\rho = \rho_{\min} + a \sn^2\left(\frac{x_v}{l};al^2\right), \\
L = 2lK(al^2)\\
\end{split}\]
from gpe.imports import *
#from mmfutils.math.special import ellipkinv
from scipy.special import ellipj, ellipk
import scipy.integrate
import scipy as sp
import gpe.bec
def sn(u, m):
return ellipj(u, m)[0]
def K(m):
return ellipk(m)
def get_QV(m=0.1, a=0.1, rho_min=1.0, N=256):
rho_max = rho_min + a
#m = a*l**2
l = np.sqrt(m/a)
L = 2*l*K(m)
dx = L/N
x = np.arange(N)*dx
rho = rho_min + a*sn(x/l, m)**2
C = np.sqrt(rho_min*rho_max*(1+rho_min*l**2))/l
Vp = np.trapz(C/rho, x)/L
Q = 2*np.pi/L
return Q, Vp
rho_min = 1.0
a = 0.1
m = 0.63315
l = np.sqrt(m/a)
L = 2*l*K(m)
Q, Vp = get_QV(m=m, a=a, rho_min=rho_min)
x0 = m.pack(s0[...])
s0[...] = m.unpack(x0)
s0.plot()
s0 = State(Nx=64, L=L, mu=10.0, psi_0=np.sqrt(rho_min),
v=Vp, m=1.0, hbar=1.0, g=1.0)
s0.plot()
plt.figure()
m = MinimizeState(s0, fix_N=True)
m.check()
s = m.minimize(psi_tol=1e-12, E_tol=1e-12)
s.plot()
s = State(Nx=128*8, L=120.0, mu=1.0, psi_0=0.5, v=v)
psi_s = s.exact_psi()
v = (np.angle(psi_s[-1]) - np.angle(psi_s[0])+np.pi)*s.hbar/s.m/(s.Lxyz[0] - s.Lxyz[0]/s.Nxyz[0])
s = State(Nx=128*8, L=120.0, mu=1.0, psi_0=0.5, v=v)
s[...] = s.exact_psi()
print(v)
plt.plot(s[...]/s.twist_phase)
v
s = s1
s.cooling_phase = 1.0
e = EvolverABM(s, dt=0.5*s.t_scale)
with NoInterrupt(ignore=True) as interrupted:
while not interrupted:
e.evolve(100)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
s = State(Nx=128*4, L=20.0, mu=1.0, psi_0=0.1, v=0.1)
s[...] = 1.0
m = MinimizeState(s, fix_N=False)
s1 = m.minimize(use_scipy=True)
plt.plot(s1.xyz[0], s1.get_density() - abs(s1.exact_psi())**2)
#s1.plot()
abs(s1.get_density() - abs(s1.exact_psi())**2).max()
plt.plot(s1.xyz[0], np.angle(s1[...]))
plt.plot(s1.xyz[0], np.angle(s1.exact_psi()))
vs = np.linspace(0, 0.4, 10)
errs = []
for v in vs:
s = State(Nx=128*4, L=20.0, mu=1.0, psi_0=0.2, v=v)
m = MinimizeState(s, fix_N=False)
s1 = m.minimize(use_scipy=True)
errs.append(abs(s1.get_density() - abs(s1.exact_psi())**2).max())
New#
\[
i \dot{\psi} = \left(E(\hat{k}) + gn \right)\psi, \qquad \psi(x,t) = \psi_v(x-vt)e^{i\phi(x,t)}
= \psi_v(x-vt)e^{-i\mu t}
\]
\[
\left(E(\hat{k}) - v \cdot k + gn \right)\psi_v = 0
\]
from gpe.Examples import traveling_waves
class StateCustomDisp(traveling_waves.StateTwist_x):
def __init__(self, Nx=128, L=10.0, mu=1.0, psi_0=1.0, ind=None,
v_p=0.0, twist=0.0, m=1.0, hbar=1.0, g=1.0):
"""
Arguments
---------
"""
self.v_p = v_p
if ind is None:
ind = Nx//2
self.ind = ind
self.psi_0 = psi_0
self.p_v = p_v = m*v_p
twist_x = twist # + p_v * L / hbar
traveling_waves.StateTwist_x.__init__(self, Nxyz=(Nx,), Lxyz=(L,),
m=m, hbar=hbar, g=g,
mu=mu, twist_x=twist_x)
self[self.ind] = psi_0
def init(self):
traveling_waves.StateTwist_x.init(self)
self._kx2 = (self.get_dispersion() * 2. * self.m / self.hbar**2)**2
def get_Vext(self):
return 0
def get_dispersion(self):
k = self._kx
disp = (k * self.hbar)**2 / 2. / self.m - self.v_p * k
return disp
s.get_Vext()
s0 = StateCustomDisp(Nx=2**6, L=5, mu=1.0, psi_0=np.sqrt(0), v_p=0.1, twist=0)
m = MinimizeStateTravellingWave(s0)
s_init = m.minimize(E_tol=1e-12, psi_tol=1e-12)
def get_err(v_p, s_guess=[s_init]):
s0 = StateCustomDisp(Nx=2**6, L=5, mu=1.0, psi_0=np.sqrt(0.1), v_p=v_p, twist=0)
if s_guess[0] is not None:
s0[...] = s_guess[0][...]
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s[s.ind] = s.psi_0
s_guess[0] = s
err = abs(s.compute_dy_dt(dy=s.copy())[...]).max()
plt.clf()
s.plot()
plt.ylabel(v_p)
plt.xlabel(err)
display(plt.gcf())
clear_output(wait=True)
return err
v_ps = np.linspace(0,1,100)
errs = map(get_err, v_ps)
plt.plot(v_ps[1:], errs[1:])
v_p = 0.19
s0 = StateCustomDisp(Nx=2**6, L=10.0, mu=1.0, psi_0=np.sqrt(0.1), v_p=v_p, twist=0)
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
sp.integrate.cumtrapz?
n = s.get_density()
L = s.basis.Lx
x = s.basis.xyz[0]
rho_min = n.min()
rho_max = n.max()
def get_rho_v_p(L=5.0, rho_min=0.1, rho_max=1.2, x=None):
if x is None:
x = np.linspace(-L/2., L/2.0, 1000)
a = rho_max - rho_min
m = get_m(a=a, L=L)
l = np.sqrt(m/a)
rho = rho_min + a*sn(x/l, m)**2
C = np.sqrt(rho_min*rho_max*(1.+rho_min*l**2))/l
plt.plot(x, rho)
v_p = np.trapz(C/rho, x)/L
u = v_p - C/rho
theta = sp.integrate.cumtrapz(u, x, initial=0)
return lambda x:rho_min + a*sn(x/l, m)**2, v_p, theta
plt.plot(x, theta)
plt.plot(x, s0[...].real)
plt.plot(x, s0[...].imag)
L = 3.0
rho_max = 1.2
rho_min = 0.5
rho, v_p, theta = get_rho_v_p(L=L, rho_min=rho_min, rho_max=rho_max)
s0 = StateCustomDisp(Nx=2**7, L=L, mu=rho_max, psi_0=np.sqrt(rho_min), v_p=v_p, twist=0)
x = s0.basis.xyz[0]
rho, v_p, theta = get_rho_v_p(L=L, rho_min=rho_min, rho_max=rho_max, x=x)
s0[...] = np.sqrt(rho(s0.basis.xyz[0]))*np.exp(1j*theta)
N_tot = s0.get_N()
Nx = s.basis.Nx
#s0[...] = np.sqrt((N_tot*Nx/L - rho_min)/(Nx-1))*np.exp(1j*s.m*x*v_p/s.hbar)
#s0[s0.ind] = s0.psi_0
#print(N_tot, s0.get_N())
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
m.check()
s.plot()
#print s[s.ind]/s.psi_0 - 1
dy = s0.compute_dy_dt(s0.copy())
plt.plot(x, abs(dy[...]))
abs(dy[...]).max()
s[s.ind] = s.psi_0
s_guess[0] = s
err = abs(s.compute_dy_dt(dy=s.copy())[...]).max()
plt.clf()
s.plot()
plt.ylabel(v_p)
plt.xlabel(err)
display(plt.gcf())
clear_output(wait=True)
return err
v_ps = np.linspac
s0 = StateCustomDisp(Nx=2**6, L=10.0, mu=1.0, psi_0=np.sqrt(1.), v_p=0.1, twist=0.0)
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-6)
s.plot()
plt.figure()
dy = s.copy()
dy = s.compute_dy_dt(dy=dy)
dy.plot()
dy = s.copy()
dy = s.compute_dy_dt(dy=dy)
plt.semilogy(abs(dy[...]))
\[\begin{split}
\psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}
\psi'(x, t) = \left(
\frac{\rho'(x,t)}{2\sqrt{\rho(x, t)}}
+ \I\phi'(x,t)\sqrt{\rho(x, t)}
\right)e^{\I\phi(x, t)}\\
\psi''(x, t) = \left(
\frac{\rho''(x,t)}{2\sqrt{\rho(x, t)}}
- \frac{[\rho'(x,t)]^2}{4\sqrt{\rho(x, t)}^3}
+ \frac{2\I\phi'(x,t)\rho'(x,t)}{2\sqrt{\rho(x, t)}}
+ \I\phi''(x,t)\sqrt{\rho(x, t)}
- [\phi'(x,t)]^2\sqrt{\rho(x, t)}
\right)
e^{\I\phi(x, t)}\\
\psi''(x, t) = \left(
\frac{\rho''(x,t)}{2\sqrt{\rho(x, t)}}
- \frac{[\rho'(x,t)]^2}{4\sqrt{\rho(x, t)}^3}
+ \frac{2\I u(x,t)\rho'(x,t)}{2\sqrt{\rho(x, t)}}
+ \I u'(x,t)\sqrt{\rho(x, t)}
- [u(x,t)]^2\sqrt{\rho(x, t)}
\right)
e^{\I\phi(x, t)}\\
\end{split}\]
\[\begin{split}
, \qquad
u(x, t) = \phi'(x, t),\\
m = al^2, \qquad
\rho(x, t) = \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)}, \qquad
\rho_{\max} = a + \rho_{\min}, \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{split}\]