Non-Polynomial Schrödinger Equation (NPSEQ) in 1D (Tube)#
Here we summarize the approach and formulation for a 1D non-polynomial Schrödinger equation representing an atomic gas trapped in a tube. We start by noting that the GPE is almost separable:
\[\begin{gather*}
\I\hbar \dot{\Psi}
=
\left(
\frac{-\hbar^2}{2m}\nabla^2_{\perp} + V(x, y)
\right)\Psi
+
\left(
\frac{-\hbar^2}{2m}\nabla^2_{z} + V(z)
\right)\Psi
+
\mathcal{E}'(n)\Psi.
\end{gather*}\]
The idea is to factor the wavefunction as
\[\begin{gather*}
\Psi(\vect{x}, t) = \phi\bigl(x, y, n_1(z, t)\bigr)\psi(z, t),\\
n_1(z, t) = \abs{\psi(z, t)}^2, \qquad
n(x, y, z, t) = n_1(z, t)\Abs{\phi\bigl(x, y, n_1(z, t)\bigr)}^2,\\
\iint\d{x}\d{y}\Abs{\phi\bigl(x, y, n_1(z, t)\bigr)}^2 = 1.
\end{gather*}\]
We then make the “adiabatic” approximation of neglecting the \(z\) and \(t\) dependence in \(\phi\) – the assumption is that the dynamics along \(z\) are slow, so that the radial wavefunction \(\phi\) adjusts instantaneously. Within this approximation, we can separate the GPE:
\[\begin{gather*}
\I\hbar \phi \dot{\psi}
=
\psi
\left(
\frac{-\hbar^2}{2m}\nabla^2_{\perp} + V(x, y)
\right)\phi
+
\phi
\left(
\frac{-\hbar^2}{2m}\nabla^2_{z} + V(z)
\right)\psi
+
\mathcal{E}'(n)\phi\psi.
\end{gather*}\]
Our assumption is that the radial wavefunction \(\phi(x, y)\) at a fixed \(z\), \(t\), and density \(n_1\) adjusts instantaneously to the ground state in the transverse direction
\[\begin{gather*}
\left(
\frac{-\hbar^2}{2m}\nabla^2_{\perp}
+ V(x, y)
+ \mathcal{E}'\Bigl(n_1\abs{\phi}^2\Bigr)
\right)\phi
= \mu_{1}(n_1)\phi.
\end{gather*}\]
Note
The key outcomes of this step are the chemical potential \(\mu_{1}(n_1)\) and the solutions \(\phi_{n_1}(x, y)\) from which we will compute various moments and related properties. Here it is important that \(V(x, y, z, t) = V_{\perp}(x, y) + V_z(z, t)\) separate so that these quantities depend only on \(n_1\) and not on \(z\) or \(t\). In principle, a more complicated dependence could be done, but this would not allow us to tabulate everything with 1-dimensional splines like we will do later.
To obtain the NPSEQ, we plug this solution back into the separated GPE above, multiply by \(\phi^\dagger\) and integrate across the transverse directions.
\[\begin{gather*}
\I\hbar \dot{\psi}
=
\left(
\frac{-\hbar^2}{2m}\nabla^2_{z} + V(z) + \mu_{1}(\abs{\psi}^2)
\right)\psi
= \frac{\delta E_1[\psi]}{\delta\psi^\dagger},\\
E_1[\psi] = \int \d{x}\left(
\frac{\hbar^2}{2m}\abs{\nabla_z\psi}^2 + V(z)\abs{\psi}^2
+ \mathcal{E}_1\bigl(\abs{\psi}^2\bigr)
\right).
\end{gather*}\]
Our goal is to tablulate \(\mathcal{E}_1(n_1)\) as an effective 1D equation of state so.
Note
The importance of tabulating \(\mathcal{E}_1(n_1)\) rather than \(\mu_{1}(n_1)\) becomes clear when we have more than one component. To maintain the variational property, one requires a consistency condition on the chemical potentials to ensure that one can consistently integrate:
\[\begin{gather*}
\mu_{i} = \diff{\mathcal{E}}{n_i}, \qquad
\pdiff{\mu_{i}}{n_{j}} = \pdiff{\mu_{j}}{n_{i}} =
\frac{\partial^2\mathcal{E}}{\partial n_i \partial n_j}=
\frac{\partial^2\mathcal{E}}{\partial n_j \partial n_i}.
\end{gather*}\]
I have not explored, but I think it might be somewhat challenging to construct a consistent \(\mathcal{E}(\vect{n})\) from a numerical tabulation of \(\mu_{i}\).
Example: The Standard dr-GPE#
As a check, consider the original formulation of the dr-GPE for axially-symmetric harmonic traps:
\[\begin{gather*}
V(x, y) = \frac{m\omega_\perp^2}{2}(x^2 + y^2), \qquad
\mathcal{E}(n) = \frac{gn^2}{2}.
\end{gather*}\]
The final result with the gaussian approximation is
\[\begin{gather*}
\mu_{\perp}(n_1) =
\frac{\hbar^2}{2m\sigma^2}
+ \frac{m\omega_\perp^2\sigma^2}{2}
+ a\frac{gn_1}{4\pi\sigma^2},\qquad
\frac{m\omega_\perp^2}{2}\sigma^4 = \frac{\hbar^2}{2m} + a\frac{gn_1}{4\pi},
\end{gather*}\]
where the constant \(a=1\) for energy minimization or \(a=2\) for the chemical-potential minimization recommended by [Mateo and Delgado, 2009].
from gpe.Examples.tutorial import StateHOConvergence2
s0 = StateHOConvergence2(g=100, Lx=30.0, Nx=128, dmu=10.0)
axs = s0.plot(label="TF")
s = s0.get_initialized_state()
axs = s.plot(axs=axs, label="Minimized", plot_convergence=True)
axs[0].legend(loc='right');
[I 20:07:38 root] Patching zope.interface.document.asReStructuredText to format code
[I 20:07:38 numexpr.utils] NumExpr defaulting to 2 threads.
Our code solves the 2D GPE with
\[\begin{gather*}
( + g_{2D} \abs{\Phi}^2)\Phi = \mu \Phi\\
\phi = \frac{\Phi}{\sqrt{N}}\\
( + g_{2D} N \abs{\phi}^2)\phi = \mu \phi\\
g_{2D} N = g n_1
\end{gather*}\]
Here we show the relative error between \(\mu(n_1)\) computed numerically and the formula above from the dr-GPE:
dmus = np.linspace(0.0001, 10.0, 100)
ss = [StateHOConvergence2(g=1, Lx=30.0, Nx=128, dmu=dmu).get_initialized_state()
for dmu in dmus]
mus = np.array([s.get_mu() for s in ss])
Ns = np.array([s.get_N() for s in ss])
gtildes = np.array([s.get_gtilde() for s in ss])
g = 1.
n1s = s.g * Ns / g
s = ss[0]
a = 2
n1 = 1
hbar2_2m = s.hbar**2 / 2 / s.m
mw2_2 = s.m * s.w**2 / 2
gn1s = s.g * Ns
agn_4pi = a * gn1s / 4 / np.pi
sigma2 = np.sqrt((hbar2_2m + agn_4pi) / mw2_2)
mus_ = (hbar2_2m + agn_4pi)/sigma2 + mw2_2 * sigma2
#mu_TFs =
fig, ax = plt.subplots()
ax.plot(n1s, mus/mus_ - 1)
ax.set(xlabel=r"$\mu_1$", ylabel=r"$\mu/\mu_{dr-GPE}-1$")
ax1 = plt.twinx()
ax1.plot(n1s, gtildes)
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[3], line 2
1 dmus = np.linspace(0.0001, 10.0, 100)
----> 2 ss = [StateHOConvergence2(g=1, Lx=30.0, Nx=128, dmu=dmu).get_initialized_state()
3 for dmu in dmus]
4 mus = np.array([s.get_mu() for s in ss])
5 Ns = np.array([s.get_N() for s in ss])
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/Examples/tutorial.py:155, in StateHOConvergence1.get_initialized_state(self, fix_N, minimize_kw)
153 m = MinimizeState(s0, fix_N=fix_N)
154 m.check()
--> 155 s1 = m.minimize(**minimize_kw)
156 s = self.copy()
157 s.set_psi(s1.get_psi())
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:770, in MinimizeState.minimize(self, psi_tol, E_tol, callback, _debug, **kw)
768 return super().minimize(_debug=_debug, **kw)
769 else:
--> 770 x = super().minimize(**kw)
772 state = self.unpack(x)
774 if self.fix_N:
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:284, in Minimize.minimize(self, plot, callback, method, polish, broyden_alpha, broyden_opts, f_tol, x_tol, use_scipy, ignore_f, _test, _debug, _log, use_cache, bounds, **kw)
282 if Version(sp.__version__) > Version("1.15.0"):
283 options.pop("disp", None)
--> 284 res = self._minimize(
285 f=_f,
286 df=_df,
287 x0=_x[0],
288 method=method,
289 callback=callback_,
290 bounds=bounds,
291 options=options,
292 )
294 self.minimize_results = res
296 if not res.success:
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:313, in Minimize._minimize(self, f, df, x0, method, callback, bounds, options)
311 def _minimize(self, f, df, x0, method, callback, bounds, options):
312 """Interface to the scipy minimizer."""
--> 313 res = sp.optimize.minimize(
314 fun=f,
315 jac=df,
316 x0=x0,
317 method=method,
318 callback=callback,
319 bounds=bounds,
320 options=options,
321 )
322 res.f = f
323 res.df = df
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_minimize.py:784, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
781 res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
782 **options)
783 elif meth == 'l-bfgs-b':
--> 784 res = _minimize_lbfgsb(fun, x0, args, jac, bounds,
785 callback=callback, **options)
786 elif meth == 'tnc':
787 res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
788 **options)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_lbfgsb_py.py:420, in _minimize_lbfgsb(fun, x0, args, jac, bounds, maxcor, ftol, gtol, eps, maxfun, maxiter, callback, maxls, finite_diff_rel_step, workers, **unknown_options)
412 _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr, pgtol, wa,
413 iwa, task, lsave, isave, dsave, maxls, ln_task)
415 if task[0] == 3:
416 # The minimization routine wants f and g at the current x.
417 # Note that interruptions due to maxfun are postponed
418 # until the completion of the current minimization iteration.
419 # Overwrite f and g:
--> 420 f, g = func_and_grad(x)
421 elif task[0] == 1:
422 # new iteration
423 n_iterations += 1
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:412, in ScalarFunction.fun_and_grad(self, x)
410 if not np.array_equal(x, self.x):
411 self._update_x(x)
--> 412 self._update_fun()
413 self._update_grad()
414 return self.f, self.g
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/optimize/_differentiable_functions.py:362, in ScalarFunction._update_fun(self)
360 def _update_fun(self):
361 if not self.f_updated:
--> 362 fx = self._wrapped_fun(self.x)
363 self._nfev += 1
364 if fx < self._lowest_f:
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/_lib/_util.py:545, in _ScalarFunctionWrapper.__call__(self, x)
542 def __call__(self, x):
543 # Send a copy because the user may overwrite it.
544 # The user of this class might want `x` to remain unchanged.
--> 545 fx = self.f(np.copy(x), *self.args)
546 self.nfev += 1
548 # Make sure the function returns a true scalar
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:160, in Minimize.minimize.<locals>._f(x)
154 if (
155 not use_cache
156 or _cache[1] is None
157 or not np.allclose(x, _cache[0], atol=1e-32, rtol=_EPS)
158 ):
159 _cache[0] = x.copy()
--> 160 _cache[1:] = self.f_df(x)
161 if _log:
162 self._calls.append(tuple(_cache))
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/minimize.py:688, in MinimizeState.f_df(self, x)
685 if self.fix_N:
686 s, N = psi.normalize()
--> 688 Hpsi = psi.get_Hy(subtract_mu=self.fix_N)
690 if self.fix_N:
691 Hpsi *= s
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec.py:895, in _StateBase.get_Hy(self, subtract_mu)
893 """Return `H(y)` for convenience only."""
894 dy = self.empty()
--> 895 self.compute_dy_dt(dy=dy, subtract_mu=subtract_mu)
896 Hy = dy / self._phase
897 return Hy
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/bec.py:796, in _StateBase.compute_dy_dt(self, dy, subtract_mu)
794 Vy = y.copy()
795 Vy.apply_V(V=self.get_V_GPU())
--> 796 Hy = Ky + Vy
798 if subtract_mu:
799 if self.constraint == "N":
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:127, in StateMixin.__add__(self, y)
125 assert isinstance(y, self.__class__)
126 res = self.copy()
--> 127 res.axpy(y)
128 return res
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:661, in ArrayStateMixin.axpy(self, x, a)
659 """Perform `self += a*x` as efficiently as possible."""
660 assert self.writeable
--> 661 self._axpy(y=self.data, x=x.data, a=a)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:762, in ArrayStateMixin._get_axpy.<locals>._axpy_blas(y, x, a, _axpy)
760 def _axpy_blas(y, x, a, _axpy=_axpy):
761 assert np.isscalar(a)
--> 762 _axpy(x.ravel(), np.asarray(y).ravel(), a=a)
KeyboardInterrupt:
The maximum error is less than 1%.