Pancake

Pancake#

import mmf_setup;mmf_setup.nbinit()
import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL)
import matplotlib.pyplot as plt
import numpy as np
%load_ext autoreload

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.

Derivation#

We start by deriving an effective 1D theory expressing \(\Psi(\vect{x}, t) = \phi(\vect{x},t)\psi(z, t)\) and varying the action wrt \(\phi(\vect{x},t)\) which gives:

\[\begin{split} \Psi(\vect{x}, t) = \phi(\vect{x}, t)\psi(z, t), \qquad \int\abs{\psi(z, t)}^2 \d{z} = 1, \\ \frac{\delta S}{\delta\phi^\dagger(\vect{x}, t)} = 0 = \int\d{z}\left\{ \I\hbar \abs{\psi}^2(z,t)\dot{\phi}(\vect{x}, t) + \I\hbar \psi^\dagger(z,t)\dot{\psi}(z,t)\phi(\vect{x}, t) - \psi^\dagger(z, t)\int \d^3{\vect{x}'}H_{\vect{x},\vect{x}'}[n](t)\phi(\vect{x}', t)\psi(z', t) \right\},\\ \I\hbar\dot{\phi}(\vect{x}, t) = E_\psi(\vect{x}, t) \phi(\vect{x}, t) + \int\d{z}\left\{ \psi^\dagger(z, t)\int \d^3{\vect{x}'}H_{\vect{x},\vect{x}'}[n](t)\phi(\vect{x}', t)\psi(z', t) \right\}, \qquad E_\psi(\vect{x}, t) = -\I\hbar\int\d{z}\left\{ \psi^\dagger(z,t)\dot{\psi}(z,t)\right\} \end{split}\]

To be more specific, we now consider a general energy density of the form:

\[ \mathcal{E}[\Psi] = \frac{\hbar^2}{2m}\abs{\nabla\Psi}^2 + f\bigl(n(\vect{x}, t)\bigr) \]

so that the Hamiltonian has the form $\( H_{\vect{x}, \vect{x}'} = \delta^3(\vect{x}-\vect{x}')\left( \frac{-\hbar^2}{2m}\nabla_{\vect{x}'}^2 + f'\bigl(n(\vect{x}, t)\bigr) \right), \qquad \int \d^3{\vect{x}'}H_{\vect{x},\vect{x}'}[n](t)\phi(\vect{x}', t)\psi(z', t) = \frac{-\hbar^2}{2m}\left( \psi(z,t)\nabla^2\phi(\vect{x}, t) +\phi(\vect{x}, t)\nabla_z^2\psi(z,t) \right) + f'\bigl(n(\vect{x}, t)\bigr)\Psi(\vect{x}, t). \)$

The normalization condition for \(\phi\) being independent of \(z\) allows us to neglect the cross-term \(2\nabla_z\phi(\vect{x}, t)\nabla_z\psi(z,t)\) (it can be removed from the Lagrangian by a total derivative which vanishes under physical boundary conditions).

The final effective 1D equation of motion can be expressed as:

\[\begin{split} \I\hbar\dot{\phi}(\vect{x},t) = \left( \frac{-\hbar^2\nabla_\vect{x}^2}{2m} + \Braket{ - \I\hbar\partial_t - \frac{\hbar^2\nabla^2_z}{2m} + f'\bigl(n(\vect{x}, t)\bigr)}_{\psi} \right)\phi(\vect{x}, t), \\ \Braket{\op{A}}_\psi = \int \psi^\dagger(z, t)\op{A}\psi(z, t)\;\d{z}. \end{split}\]

The same result can be obtained by expanding the action:

\[ S = \int \d{x}\d{y}\left\{ \I\hbar \phi^\dagger\dot{\phi} +\Braket{\I\hbar\partial_t}_{\psi}\phi^\dagger\phi + \int\d{z}\; \mathcal{E}(\Psi) \right\} \]

For the usual energy density, this becomes

\[ S = \int \d{x}\d{y}\left\{ \I\hbar \phi^\dagger\dot{\phi} +\phi^\dagger\Braket{\I\hbar\partial_t + \frac{\hbar^2\nabla_z^2}{2m} -f'(n) }_{\psi}\phi \right\} \]

Gaussian axial profile#

A common strategy is to assume an axial Gaussian:

\[\begin{split} \psi\bigl(z; \sigma(\vect{x}, t)\bigr) = \frac{e^{\I\mu t/\hbar}}{\sqrt{\pi}\sigma(\vect{x}, t)}e^{-z^2/2\sigma^2}, \\ \frac{\frac{-\hbar^2 \nabla_z^2}{2m}\psi}{\psi} = \frac{\hbar^2}{2m}\frac{\sigma^2-z^2}{\sigma^4},\\ \Braket{\frac{\hbar^2}{2m}\nabla_z^2}_{\psi} = -\frac{\hbar^2}{4m\sigma^2}, \qquad \Braket{\abs{\psi}^2} = \frac{1}{\sqrt{2}}\frac{1}{(\sqrt{\pi}\sigma)^2}, \qquad \Braket{z^2} = \frac{\sigma^2}{2}, \qquad \Braket{\nabla_\vect{x}^2}_\psi = -\frac{\sigma^2_{,\vect{x}}(\vect{x}, t)}{\sigma^2}, \qquad \Braket{\I\hbar\partial_t} = -\mu. \end{split}\]

For \(f(n) = m\omega_z^2 z^2n/2 + gn^2/2\), we thus have

\[ S = \int \d{x}\d{y}\left\{ \I\hbar \phi^\dagger\dot{\phi} +\phi^\dagger\left(-\mu -\frac{\hbar^2\sigma^2_{,\vect{x}}(\vect{x}, t)}{2m\sigma^2} - \frac{\hbar^2}{4m\sigma^2} - \frac{m\omega_z^2\sigma^2}{4} - \frac{g\abs{\phi}^2}{2\sqrt{2}\pi\sigma^2} + \frac{\hbar^2\nabla_\vect{x}^2}{2m} \right)\phi \right\} \]

Ignore the term:

\[ \frac{\hbar^2\sigma_{,\vect{x}}^2}{2m\sigma^2} \]

assuming that the radial variations of \(\sigma(\vect{x}, t)\) are small compared with the axial energies (something that can be tested during a simulation). Varying this with respect to \(\sigma^2\) and multiplying by \(\sigma^4\), we obtain the equation:

\[ \frac{m\omega_z^2}{4}\sigma^4 = \frac{g\abs{\phi}^2}{2\pi\sqrt{2}}+\frac{\hbar^2}{4m} +\frac{\hbar^2(\nabla_\vect{x}\sigma^2)^2}{4m\sigma^4} + \frac{\hbar^2\sigma^2}{4m\phi^\dagger\phi}\nabla_\vect{x}\frac{\phi^\dagger\phi\nabla_\vect{x}\sigma^2}{\sigma^4}, \]

which, after neglecting the \(\sigma_{,\vect{x}}\) terms gives the usual NPSEQ condition:

\[ \frac{m\omega_z^2}{4}\sigma^4 = \frac{g\abs{\phi}^2}{2\pi\sqrt{2}}+\frac{\hbar^2}{4m}. \]

We can use a few more definitions to make this result simpler.

Recall the definition of the scattering length and the harmonic oscillator length as $\( a = \frac{m g}{4 \pi \hbar^2},\\ a_\perp = \sqrt{\frac{\hbar}{m \omega_\perp}} \)$ respectively.

With some algebra, we get: $\( \sigma^4 = a_\perp^4 \left( 1 + 4\sqrt{2} a \left| \phi \right|^2 \right), \)$

and arguing that only the positive root is relevant; we arrive at the simplified NPSEQ condition: $\( \sigma^2 = a_\perp^2 \sqrt{\left( 1 + 4\sqrt{2} a \left| \phi \right|^2 \right)}. \)$

Central Density#

If the Gaussian is a good approximation, then the central density is related to the 1D integrated density by:

\[ n(x, y, 0) = \frac{n_{2D}(x, y)}{\lambda_\perp \sqrt{\pi} \sigma}. \]

As demonstrated below, if the radial direction is large, the Gaussian ansatz may be poor even though the resulting dynamics might be well captured by this model. In this case, a different approach is required to get the central density based on the TF approximation:

\[\begin{split} n(\vect{x}) = \frac{\mu - \frac{m}{2}\sum_i\omega_i^2x_i^2}{g}, \qquad n_0 = n(x, y, 0) = \frac{\mu - \frac{m}{2}\omega_x^2x^2 - \frac{m}{2}\omega_y^2y^2}{g},\\ n_{2D}(x, y) = \int\d{z} n(\vect{x}) = \frac{m}{2g}\frac{1}{\omega_z} \int_{-L_z/2}^{L_z/2}\d{\tilde{z}}\; (\tilde{R}^2 - \tilde{z}^2) = \frac{2 L_z gn_0^2}{m \omega_z}, \qquad \tilde{R}^2 = \frac{2\mu}{m} - \omega_x^2x^2 - \omega_y^2y^2 = \frac{2gn_0}{m},\\ n(x, y, 0) = \frac{1}{\lambda_\perp}\sqrt{\frac{m \omega_z n_{2D}(x, y)}{2 L_z g}}. \end{split}\]

This latter option is provided in by get_central_density(TF=True).

Thomas Fermi Profiles#

To initialize states we typically use a Thomas Fermi approximation. For the usual 2D case, this leads to a profile:

\[ n(x, y) = \abs{\phi}^2 = \frac{\mu - V(x, y)}{g} = \frac{V(R_{TF}) - V(x, y)}{g}. \]

where \(R_{TF}\) is the Thomas Fermi “radius” (for a harmonic potential). (Of course, only positive densities should be used, if the expression is negative, the density should be set to zero.)

For the pancake effective 2D model, this must be modified somewhat. We start with the definition of the Thomas Fermi where \(n_{2D}(R) = 0\), hence \(\sigma^2(R) = \hbar/m\omega_\perp\). At this point, the chemical potential must still be augmented with the transverse energy:

\[ \mu = V(R_{TF}) + \hbar\omega_\perp. \]

(Note: we have set \(\lambda_\perp = 1\) here since we only consider initializing a state prior to expansion.)

\[ \sigma^2\left(\mu - V(x, y)\right) = + \frac{\hbar^2}{4m} + \frac{m\omega_\perp^2\sigma^4}{4} + \frac{g\abs{\Phi}^2}{2\sqrt{2}\pi}, \qquad \sigma^2 = a_\perp^2 \sqrt{\left( 1 + 4\sqrt{2} a \left| \phi \right|^2 \right)}. \]
\[\begin{split} \frac{3m\omega_\perp^2\sigma^4}{2} - \sigma^2[\mu - V(x, y)] - \frac{\hbar^2}{2m} = 0, \qquad n_{1D} = \abs{\Phi}^2 = 2\pi m \frac{(\omega_\perp\sigma^2)^2 - \frac{\hbar^2}{m^2}}{g}\\ \omega_\perp \sigma^2 = \frac{\hbar}{3m}\left( \frac{[\mu - V(x, y)]}{\hbar\omega_\perp} + \sqrt{\frac{[\mu - V(x, y)]^2}{\hbar^2\omega_\perp^2} + 3} \right) \end{split}\]

Energy density#

Since we have an expression for the variance in the NPSEQ formalism, $\( \sigma^2 = a_\perp^2 \sqrt{\left( 1 + 4\sqrt{2} a \left| \phi \right|^2 \right)} \)$ we can take it’s first 3 derivatives to see if we can express the energy as a total derivative of some quantity, which we can then interpret as the energy density.

(22)#\[\begin{align} \frac{\partial \sigma^2}{\partial \phi^\dagger} &= \frac{2 \sqrt{2} a_\perp^4}{\sigma^2} \phi,\\ \frac{\partial \sigma^4}{\partial \phi^\dagger} &= 4 \sqrt{2} a a_\perp^4 \phi,\\ \frac{\partial \sigma^6}{\partial \phi^\dagger} &= 6 \sqrt{2} a a_\perp^4 \sigma^2 \phi. \end{align}\]

The energy is given as:

(23)#\[\begin{align} \mathcal{E} = \frac{\hbar^2 |\vect{\nabla} \phi|^2}{2 m} + \left\{ V(x, y, t) + \frac{1}{2} m \left( \omega_x^2 x^2 + \omega_y^2 y^2 \right) \right\} |\phi|^2 + \frac{m \omega_\perp^2 \sigma^6}{12 \sqrt{2} a a_\perp}. \end{align}\]

Demonstration#

%autoreload
from gpe import disc

expt = disc.Experiment()
#state = expt.get_initial_state()
state = expt.get_state()
state.plot()
#plt.imshow(state.get_density())

[I 20:01:26 root] Patching zope.interface.document.asReStructuredText to format code
[I 20:01:26 numexpr.utils] NumExpr defaulting to 2 threads.

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 2
      1 get_ipython().run_line_magic('autoreload', '')
----> 2 from gpe import disc
      3 
      4 expt = disc.Experiment()
      5 #state = expt.get_initial_state()

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/disc.py:27
     22 # from mmfutils.math.bases import PeriodicBasis
     23 # from mmfutils.optimize import bracket_monotonic
     25 from pytimeode.evolvers import EvolverABM
---> 27 import pytest
     29 piston = pytest.importorskip("piston")
     30 from piston.mixins import EvolveMixin

ModuleNotFoundError: No module named 'pytest'

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/stack_data/core.py:391, in Line.render(self, markers, strip_leading_indent, pygmented, escape_html)
    389 if pygmented and self.frame_info.scope:
    390     assert_(not markers, ValueError("Cannot use pygmented with markers"))
--> 391     start_line, lines = self.frame_info._pygmented_scope_lines
    392     result = lines[self.lineno - start_line]
    393     if strip_leading_indent:

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/stack_data/utils.py:145, in cached_property.cached_property_wrapper(self, obj, _cls)
    142 if obj is None:
    143     return self
--> 145 value = obj.__dict__[self.func.__name__] = self.func(obj)
    146 return value

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/stack_data/core.py:824, in FrameInfo._pygmented_scope_lines(self)
    821     ranges = []
    823 code = atext.get_text(scope)
--> 824 lines = _pygmented_with_ranges(formatter, code, ranges)
    826 start_line = self.source.line_range(scope)[0]
    828 return start_line, lines

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/stack_data/utils.py:166, in _pygmented_with_ranges(formatter, code, ranges)
    164 lexer = MyLexer(stripnl=False)
    165 try:
--> 166     highlighted = pygments.highlight(code, lexer, formatter)
    167 except Exception:
    168     # When pygments fails, prefer code without highlighting over crashing
    169     highlighted = code

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pygments/__init__.py:82, in highlight(code, lexer, formatter, outfile)
     77 def highlight(code, lexer, formatter, outfile=None):
     78     """
     79     This is the most high-level highlighting function. It combines `lex` and
     80     `format` in one function.
     81     """
---> 82     return format(lex(code, lexer), formatter, outfile)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pygments/__init__.py:64, in format(tokens, formatter, outfile)
     62 if not outfile:
     63     realoutfile = getattr(formatter, 'encoding', None) and BytesIO() or StringIO()
---> 64     formatter.format(tokens, realoutfile)
     65     return realoutfile.getvalue()
     66 else:

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pygments/formatters/terminal256.py:250, in Terminal256Formatter.format(self, tokensource, outfile)
    249 def format(self, tokensource, outfile):
--> 250     return Formatter.format(self, tokensource, outfile)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pygments/formatter.py:124, in Formatter.format(self, tokensource, outfile)
    121 if self.encoding:
    122     # wrap the outfile in a StreamWriter
    123     outfile = codecs.lookup(self.encoding)[3](outfile)
--> 124 return self.format_unencoded(tokensource, outfile)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pygments/formatters/terminal256.py:261, in Terminal256Formatter.format_unencoded(self, tokensource, outfile)
    258 while ttype and not_found:
    259     try:
    260         # outfile.write( "<" + str(ttype) + ">" )
--> 261         on, off = self.style_string[str(ttype)]
    263         # Like TerminalFormatter, add "reset colors" escape sequence
    264         # on newline.
    265         spl = value.split('\n')

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pygments/token.py:44, in _TokenType.__repr__(self)
     43 def __repr__(self):
---> 44     return 'Token' + (self and '.' or '') + '.'.join(self)

KeyboardInterrupt: 

Nx, Ny = state.basis.Nxyz
x = state.x
n = state.get_density()
plt.plot(x, n[Nx//2, :])

state.evolve();