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.
Travelling Waves#
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\).
K(0.9)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 1
----> 1 K(0.9)
NameError: name 'K' is not defined
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec
def sn(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.0 * 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])
[I 20:02:16 numexpr.utils] NumExpr defaulting to 2 threads.
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[3], line 1
----> 1 from gpe.imports import *
2 from scipy.special import ellipj, ellipk
3 import gpe.bec
4
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/imports.py:28
26 from mmfutils.plot import imcontourf # noqa: E402
27 from gpe.minimize import MinimizeState # noqa: E402
---> 28 from gpe.utils import evolve_to, evolve, evolves # noqa: E402
29 from gpe.plot_utils import MPLGrid # noqa: E402
30 from pytimeode.evolvers import EvolverSplit, EvolverABM # noqa: E402
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/utils.py:35
32 from persist.objects import Archivable
33 from persist.archive import Archive
---> 35 from pytimeode.evolvers import EvolverABM
36 from pytimeode.mixins import ArrayStateMixin
37 from pytimeode.interfaces import implementer, IStateForABMEvolvers
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/__init__.py:2
1 from . import interfaces
----> 2 from . import mixins
3 from . import evolvers
5 __all__ = ["interfaces", "mixins", "evolvers"]
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:39
30 from .interfaces import (
31 IState,
32 IStateApply,
(...) 36 implementer,
37 )
38 from . import interfaces
---> 39 from .utils import expr
41 __all__ = [
42 "StateMixin",
43 "StatesMixin",
(...) 47 "ArraysStateWithBraketMixin",
48 ]
51 @implementer(IState)
52 class StateMixin:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/utils/expr.py:18
15 import numpy as np
17 try:
---> 18 import sympy
19 except ImportError:
20 sympy = None
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/__init__.py:77
70 from .logic import (to_cnf, to_dnf, to_nnf, And, Or, Not, Xor, Nand, Nor,
71 Implies, Equivalent, ITE, POSform, SOPform, simplify_logic, bool_map,
72 true, false, satisfiable)
74 from .assumptions import (AppliedPredicate, Predicate, AssumptionsContext,
75 assuming, Q, ask, register_handler, remove_handler, refine)
---> 77 from .polys import (Poly, PurePoly, poly_from_expr, parallel_poly_from_expr,
78 degree, total_degree, degree_list, LC, LM, LT, pdiv, prem, pquo,
79 pexquo, div, rem, quo, exquo, half_gcdex, gcdex, invert,
80 subresultants, resultant, discriminant, cofactors, gcd_list, gcd,
81 lcm_list, lcm, terms_gcd, trunc, monic, content, primitive, compose,
82 decompose, sturm, gff_list, gff, sqf_norm, sqf_part, sqf_list, sqf,
83 factor_list, factor, intervals, refine_root, count_roots, all_roots,
84 real_roots, nroots, ground_roots, nth_power_roots_poly, cancel,
85 reduced, groebner, is_zero_dimensional, GroebnerBasis, poly,
86 symmetrize, horner, interpolate, rational_interpolate, viete, together,
87 BasePolynomialError, ExactQuotientFailed, PolynomialDivisionFailed,
88 OperationNotSupported, HeuristicGCDFailed, HomomorphismFailed,
89 IsomorphismFailed, ExtraneousFactors, EvaluationFailed,
90 RefinementFailed, CoercionFailed, NotInvertible, NotReversible,
91 NotAlgebraic, DomainError, PolynomialError, UnificationFailed,
92 GeneratorsError, GeneratorsNeeded, ComputationFailed,
93 UnivariatePolynomialError, MultivariatePolynomialError,
94 PolificationFailed, OptionError, FlagError, minpoly,
95 minimal_polynomial, primitive_element, field_isomorphism,
96 to_number_field, isolate, round_two, prime_decomp, prime_valuation,
97 galois_group, itermonomials, Monomial, lex, grlex,
98 grevlex, ilex, igrlex, igrevlex, CRootOf, rootof, RootOf,
99 ComplexRootOf, RootSum, roots, Domain, FiniteField, IntegerRing,
100 RationalField, RealField, ComplexField, PythonFiniteField,
101 GMPYFiniteField, PythonIntegerRing, GMPYIntegerRing, PythonRational,
102 GMPYRationalField, AlgebraicField, PolynomialRing, FractionField,
103 ExpressionDomain, FF_python, FF_gmpy, ZZ_python, ZZ_gmpy, QQ_python,
104 QQ_gmpy, GF, FF, ZZ, QQ, ZZ_I, QQ_I, RR, CC, EX, EXRAW,
105 construct_domain, swinnerton_dyer_poly, cyclotomic_poly,
106 symmetric_poly, random_poly, interpolating_poly, jacobi_poly,
107 chebyshevt_poly, chebyshevu_poly, hermite_poly, hermite_prob_poly,
108 legendre_poly, laguerre_poly, apart, apart_list, assemble_partfrac_list,
109 Options, ring, xring, vring, sring, field, xfield, vfield, sfield)
111 from .series import (Order, O, limit, Limit, gruntz, series, approximants,
112 residue, EmptySequence, SeqPer, SeqFormula, sequence, SeqAdd, SeqMul,
113 fourier_series, fps, difference_delta, limit_seq)
115 from .functions import (factorial, factorial2, rf, ff, binomial,
116 RisingFactorial, FallingFactorial, subfactorial, carmichael,
117 fibonacci, lucas, motzkin, tribonacci, harmonic, bernoulli, bell, euler,
(...) 138 Znm, elliptic_k, elliptic_f, elliptic_e, elliptic_pi, beta, mathieus,
139 mathieuc, mathieusprime, mathieucprime, riemann_xi, betainc, betainc_regularized)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/__init__.py:124
118 from .orthopolys import (jacobi_poly, chebyshevt_poly, chebyshevu_poly,
119 hermite_poly, hermite_prob_poly, legendre_poly, laguerre_poly)
121 from .appellseqs import (bernoulli_poly, bernoulli_c_poly, genocchi_poly,
122 euler_poly, andre_poly)
--> 124 from .partfrac import apart, apart_list, assemble_partfrac_list
126 from .polyoptions import Options
128 from .rings import ring, xring, vring, sring
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/partfrac.py:13
9 from sympy.polys.polytools import parallel_poly_from_expr
10 from sympy.utilities import numbered_symbols, take, xthreaded, public
---> 13 @xthreaded
14 @public
15 def apart(f, x=None, full=False, **options):
16 """
17 Compute partial fraction decomposition of a rational function.
18
(...) 67 apart_list, assemble_partfrac_list
68 """
69 allowed_flags(options, [])
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/utilities/decorator.py:82, in xthreaded(func)
65 def xthreaded(func):
66 """Apply ``func`` to sub--elements of an object, excluding :class:`~.Add`.
67
68 This decorator is intended to make it uniformly possible to apply a
(...) 80
81 """
---> 82 return threaded_factory(func, False)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/utilities/decorator.py:19, in threaded_factory(func, use_add)
17 """A factory for ``threaded`` decorators. """
18 from sympy.core import sympify
---> 19 from sympy.matrices import MatrixBase
20 from sympy.utilities.iterables import iterable
22 @wraps(func)
23 def threaded_func(expr, *args, **kwargs):
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/matrices/__init__.py:8
6 from .exceptions import ShapeError, NonSquareMatrixError
7 from .kind import MatrixKind
----> 8 from .dense import (
9 GramSchmidt, casoratian, diag, eye, hessian, jordan_cell,
10 list2numpy, matrix2numpy, matrix_multiply_elementwise, ones,
11 randMatrix, rot_axis1, rot_axis2, rot_axis3, rot_ccw_axis1,
12 rot_ccw_axis2, rot_ccw_axis3, rot_givens,
13 symarray, wronskian, zeros)
14 from .dense import MutableDenseMatrix
15 from .matrixbase import DeferredVector, MatrixBase
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/matrices/dense.py:15
13 from .exceptions import ShapeError
14 from .decompositions import _cholesky, _LDLdecomposition
---> 15 from .matrixbase import MatrixBase
16 from .repmatrix import MutableRepMatrix, RepMatrix
17 from .solvers import _lower_triangular_solve, _upper_triangular_solve
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/matrices/matrixbase.py:40
38 from sympy.core.decorators import call_highest_priority
39 from sympy.core.logic import fuzzy_and, FuzzyBool
---> 40 from sympy.tensor.array import NDimArray
41 from sympy.utilities.iterables import NotIterable
43 from .utilities import _get_intermediate_simp_bool
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/tensor/__init__.py:7
5 from .index_methods import get_contraction_structure, get_indices
6 from .functions import shape
----> 7 from .array import (MutableDenseNDimArray, ImmutableDenseNDimArray,
8 MutableSparseNDimArray, ImmutableSparseNDimArray, NDimArray, tensorproduct,
9 tensorcontraction, tensordiagonal, derive_by_array, permutedims, Array,
10 DenseNDimArray, SparseNDimArray,)
12 __all__ = [
13 'IndexedBase', 'Idx', 'Indexed',
14
(...) 22 'Array', 'DenseNDimArray', 'SparseNDimArray',
23 ]
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/tensor/array/__init__.py:251
1 r"""
2 N-dim array module for SymPy.
3
(...) 248
249 """
--> 251 from .dense_ndim_array import MutableDenseNDimArray, ImmutableDenseNDimArray, DenseNDimArray
252 from .sparse_ndim_array import MutableSparseNDimArray, ImmutableSparseNDimArray, SparseNDimArray
253 from .ndim_array import NDimArray, ArrayKind
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/tensor/array/dense_ndim_array.py:8
6 from sympy.core.singleton import S
7 from sympy.core.sympify import _sympify
----> 8 from sympy.tensor.array.mutable_ndim_array import MutableNDimArray
9 from sympy.tensor.array.ndim_array import NDimArray, ImmutableNDimArray, ArrayKind
10 from sympy.utilities.iterables import flatten
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/tensor/array/mutable_ndim_array.py:1
----> 1 from sympy.tensor.array.ndim_array import NDimArray
4 class MutableNDimArray(NDimArray):
6 def as_immutable(self):
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/tensor/array/ndim_array.py:591
586 raise ValueError('Dimension of index greater than rank of array')
588 return index
--> 591 class ImmutableNDimArray(NDimArray, Basic):
592 _op_priority = 11.0
594 def __hash__(self):
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/core/basic.py:218, in Basic.__init_subclass__(cls)
212 def __init_subclass__(cls):
213 # Initialize the default_assumptions FactKB and also any assumptions
214 # property methods. This method will only be called for subclasses of
215 # Basic but not for Basic itself so we call
216 # _prepare_class_assumptions(Basic) below the class definition.
217 super().__init_subclass__()
--> 218 _prepare_class_assumptions(cls)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/core/assumptions.py:624, in _prepare_class_assumptions(cls)
622 for k in _assume_defined:
623 attrname = as_property(k)
--> 624 v = cls.__dict__.get(attrname, '')
625 if isinstance(v, (bool, int, type(None))):
626 if v is not None:
KeyboardInterrupt:
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.2 * s.t_scale)
s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])
with NoInterrupt(ignore=True) as interrupted:
while e.t < e.y.T and not interrupted:
e.evolve(1000)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
%debug
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:
(28)#\[\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.0])
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 following hypothesis:
The 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}
\]
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec, gpe.minimize
reload(gpe.bec)
class State(gpe.bec.StateTwist):
def __init__(
self,
Nx=32,
L=10.0,
mu=1.0,
psi_0=1.0,
ind=None,
v=0.0,
twist=None,
m=1.0,
hbar=1.0,
g=1.0,
):
if ind is None:
ind = Nx // 2
self.ind = ind
self.psi_0 = psi_0
self.p_v = p_v = m * v
if twist is None:
twist = np.pi + p_v * L / hbar
gpe.bec.StateTwist.__init__(
self, Nxyz=(Nx,), Lxyz=(L,), m=m, hbar=hbar, g=g, mu=mu, twist=(twist,)
)
self[self.ind] = psi_0
def exact_psi(self):
"""Return the analytic soliton"""
v = abs(self.psi_0)
u = np.sqrt(self.mu - v ** 2)
return 1j * v + u * np.tanh(u * self.xyz[0])
def get_Vext(self):
return -(self.mu + self.p_v ** 2 / 2.0 / self.m)
def _compute_dy_dt(self, dy, subtract_mu=False):
dy = gpe.bec.StateTwist.compute_dy_dt(self, dy=dy, subtract_mu=subtract_mu)
# dy[self.ind] = 0
return dy
class MinimizeState(gpe.minimize.MinimizeState):
def __init__(self, state, **kw):
self.psi_0 = state.psi_0
self.ind = state.Nxyz[0] // 2
gpe.minimize.MinimizeState.__init__(self, state, **kw)
def pack(self, psi):
psi[self.ind] = self.psi_0
fact = 1 if self.real else 2
x = gpe.minimize.MinimizeState.pack(self, psi)
return np.concatenate([x[: self.ind * fact], x[(self.ind + 1) * fact :]])
def unpack(self, x, state=None):
fact = 1 if self.real else 2
x = np.concatenate([x[: self.ind * fact], [0, 0], x[self.ind * fact :]])
state = gpe.minimize.MinimizeState.unpack(self, x, state=state)
assert np.allclose(state[self.ind], 0)
state[self.ind] = self.psi_0
# plt.clf()
# state.plot()
# plt.twinx()
# plt.plot(state.xyz[0], np.angle(state[...]/state.twist_phase), 'r--')
# display(plt.gcf())
# clear_output(wait=True)
return state
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())
plt.plot(vs, errs, "-+")
s = gpe.bec.State(Nxyz=(64,), Lxyz=(10.0,))
s[...] = 1.0
m = gpe.minimize.MinimizeState(
fix_N=False,
)
s1 = m.minimize(use_scipy=True)
s1.plot()
Es = []
m = MinimizeState(s)
m.check()
s1 = m.minimize(use_scipy=True, fix_N=False, callback=callback)
# plt.plot(s1.xyz[0], np.log10(abs(abs(s1[...])**2-s.mu)))
s1.plot()
plt.plot(s1.xyz[0], abs(s1.exact_psi()) ** 2)
plt.plot(Es)
s1.get_density().max(), 1 + (s1.k_B[0] ** 2) / 2
-s1.mu, s1.get_Vext()
Es = []
twists = np.linspace(0, 2 * np.pi, 10)
for twist in twists:
s = State(Nx=128, L=50.0, mu=1.0, psi_0=0.05, twist=twist)
m = MinimizeState(s)
s1 = m.minimize(use_scipy=True, fix_N=False)
Es.append(s1.get_energy())
print twist, Es[-1]
s1.plot()
plt.plot(twists, Es)
m = MinimizeState(s)
s1 = m.minimize(use_scipy=True)
plt.clf()
s1.plot()
s1.get_energy()