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 03:38:02 root] Patching zope.interface.document.asReStructuredText to format code
[I 03:38:02 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'

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 Nx, Ny = state.basis.Nxyz
      2 x = state.x
      3 n = state.get_density()
      4 plt.plot(x, n[Nx//2, :])

NameError: name 'state' is not defined

state.evolve();

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 state.evolve();

NameError: name 'state' is not defined