Non-Polynomial Schrödinger Equation (NPSEQ)

Non-Polynomial Schrödinger Equation (NPSEQ)#

Hide code cell content

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

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.

Introduction#

The goal here is to derive a dimensionally-reduced dynamical description for trapped quantum gasses of the GPE form

\[\begin{gather*} \I\hbar \psi(\vect{x}, t) = \Bigl[ \frac{-\hbar^2\nabla^2}{2m} + V(\vect{x}, t) + \mu_{\text{eff}}\bigl(n(\vect{x}, t)\bigr)\Bigr] \psi(\vect{x},t)\,, \qquad n(\vect{x}, t) = \abs{\psi(\vect{x},t)}^2\,, \end{gather*}\]

where \(\mathcal{E}_{\mu}(n)\) is an effective non-linear (and generally non-polynomial) interaction mocking up the transverse structure of the cloud. We will consider two cases here:

  • Tube: \(\psi(\vect{x}, t) = \psi_1(z, t)\) where \(n = n_1\) is a one-dimensional density (dimensions \(1/L\)):

    \[\begin{gather*} n_1(z, t) = \iint\d{x}\d{y}\; n_{3D}(x, y, z; t). \end{gather*}\]

    This is relevant for elongated tube geometries where the dynamics occur primarily along the \(z\) axis.

We start with a high-level formulation of the GPE as a minimum action principle from an energy density functional of the form:

\[\begin{gather*} S[\Psi(\vect{x}, t)] = \int \d{t}\d^3{\vect{x}}\left\{ \I\hbar\frac{\Psi^\dagger\dot{\Psi} - \dot{\Psi}^\dagger\Psi}{2} - \mathcal{E}[\Psi] \right\}, \\ \frac{\delta S}{\delta \Psi^*(\vect{x}, t)} = 0 = \I\hbar \dot{\Psi} - \mat{H}[n]\cdot\Psi = \I\hbar \dot{\Psi} - \int \d^3{\vect{x}'}H_{\vect{x},\vect{x}'}[n](t)\Psi(\vect{x}', t),\\ \I\hbar \dot{\Psi} = \mat{H}[n]\cdot\Psi. \end{gather*}\]

Here we assume that the non-linear portions of \(\mat{H}\) depend only on the local density \(n(\vect{x}, t) = \abs{\Psi(\vect{x}, t)}^2\). The usual GPE follows from:

\[\begin{gather*} \mathcal{E}[\Psi] = \frac{\hbar^2}{2m}\abs{\vect{\nabla}\Psi(\vect{x}, t)}^2 + V(\vect{x}, t)n(\vect{x}, t) + \frac{g n^2(\vect{x}, t)}{2}, \\ \I\hbar \dot{\Psi}(\vect{x}, t) = \left[-\frac{\hbar^2}{2m}\nabla^2 + V(\vect{x}, t) + gn(\vect{x}, t)\right]\Psi(\vect{x}, t). \end{gather*}\]

The purpose of this notebook is to consider a reformulation of the problem with a factorized wavefunction:

\[\begin{gather*} \Psi(\vect{x}, t) = \phi(\vect{x}, t)\psi(z, t) \end{gather*}\]

where assumptions will be made about the transverse wavefunction \(\phi\) that yield approximate 1D effective theories for \(\psi(z, t)\).

Summary#

Effective 1D Theory#

The procedure here is to write the wavefunction as

\[\begin{gather*} \Psi(\vect{x}, t) = \phi(\vect{x}, t) \psi(z, t), \qquad \phi(\vect{x}, t) = \frac{1}{\sqrt{\pi\sigma_x(z, t)\sigma_y(z,t)}} \exp\left(-\frac{x^2}{2\sigma_x^2} - \frac{y^2}{2\sigma_y^2}\right), \end{gather*}\]

and then minimize the action with respect to \(\sigma_{x,y}(z, t)\) assuming that:

  1. Terms proportional to \(\nabla_z\sigma_{x,y}\) can be neglected.

  2. That \(\phi\) responds instantaneously adjusting to the size that minimizes the actions in the radial direction. This is the reason that this step must be applied after removing the radial scaling. Consider for example the complete removal of the traps. The instantaneous response of the radial function would immediately send \(\sigma \rightarrow \infty\).

Effective 2D Theory (pancake code)

Analogously, we can consider

\[\begin{gather*} \Psi(\vect{x}, t) = \phi(\vect{x}, t) \psi(x, y, t), \qquad \phi(\vect{x}, t) = \frac{1}{\sqrt{\sqrt{\pi}\sigma(x, y, t)}}e^{-z^2/2\sigma^2(x, y, t)}, \end{gather*}\]

and then minimize the action with respect to \(\sigma(x, y, t)\). We will include results for this “pancake” approximation in boxes like this where needed.

We start with the simplest results for the GPE where \(m_x=m_y=m_z=m\), \(\omega_x = \omega_y = \omega_\perp\), and \(\sigma_{x} = \sigma_{y} = \sigma\), which gives the following effective theory for \(\psi(z, t)\):

\[\begin{gather*} \mathcal{E} = \frac{\hbar^2\abs{\nabla_z\psi}^2}{2m} + \Biggl( \overbrace{ \frac{\hbar^2\sigma^2_{,z}(z, t)}{2m\sigma^2} }^{\text{neglected}} + \overbrace{ \frac{\hbar^2}{2m\sigma^2} + \frac{m\omega_\perp^2\sigma^2}{2} + \frac{g\abs{\psi}^2}{4\pi\sigma^2} }^{\text{energy to minimize}} \Biggr)\psi^\dagger\psi,\\ \I\hbar \dot{\psi} = \Biggl( \frac{-\hbar^2\nabla_z^2}{2m} + \underbrace{ \frac{\hbar^2}{2m\sigma^2} + \frac{m\omega_\perp^2\sigma^2}{2} + \frac{g\abs{\psi}^2}{2\pi\sigma^2} }_{\text{chemical potential to minimize}} \Biggr)\psi,\\ \frac{m\omega_\perp^2}{2}\sigma^4 = \underbrace{ \frac{\hbar^2}{2m} + \frac{g\abs{\psi}^2}{4\pi} }_{\text{energy min.}} , \quad \text{or} \quad \frac{m\omega_\perp^2}{2}\sigma^4 = \underbrace{ \frac{\hbar^2}{2m} + \frac{g\abs{\psi}^2}{2\pi} }_{\text{chemical potential min.}}. \end{gather*}\]

Note

The final equation for \(\sigma^4\) has two forms, depending on whether one minimizes the energy or the chemical potential. Minimizing the energy is certainly the correct thing to do if one minimizes over all forms of the radial function \(\phi(r_{\perp})\), but [Mateo and Delgado, 2009] argues that minimizing over the chemical potential can give better results when one only considers restricted forms of \(\phi(r_\perp)\), such as the gaussian considered here.

Summary Formula for HO + Gaussian

  • Single component.

  • Harmonic transverse trap:

    \[\begin{gather*} V_{\perp}(\vect{r}_{\perp}) = \frac{m\omega_x^2 x^2 + m\omega_y^2 y^2}{2}. \end{gather*}\]
\[\begin{gather*} E_1 = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{Q_1(n_0)}{n_0},\\ \mu_1 = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{R_1(n_0)}{n_0},\\ n_1 = \abs{\phi}^2, \qquad n_0 = \frac{n_1}{\pi\sigma_x\sigma_y}, \qquad \pdiff{n_0}{\sigma_{x,y}^2} = -\frac{n_0}{2\sigma_{x,y}^2},\\ Q_1(n_0) = \int_0^{n_0}\d{n}\; \frac{\mathcal{E}(n)}{n}, \qquad R_1(n_0) = \int_0^{n_0}\d{n}\; \mathcal{E}'(n) = \mathcal{E}(n_0) - \mathcal{E}(0). \end{gather*}\]

Pancake Version

\[\begin{gather*} E_2 = \frac{\hbar^2}{2m}\frac{1}{2\lambda_z^2\sigma^2} + \frac{m}{2}\frac{\omega_z^2\sigma^2}{2} + \frac{Q_2(n_0)}{n_0},\\ \mu_2 = \frac{\hbar^2}{2m}\frac{1}{2\lambda_z^2\sigma^2} + \frac{m}{2}\frac{\omega_z^2\sigma^2}{2} + \frac{R_2(n_0)}{n_0},\\ n_2 = \abs{\phi_2}^2, \qquad n_0 = \frac{n_2}{\sqrt{\pi}\sigma}, \qquad \pdiff{n_0}{\sigma^2} = -\frac{n_0}{2\sigma^2},\\ Q_2(n_0) = \int_0^{n_0}\d{n}\; \frac{\mathcal{E}(n)}{\sqrt{\pi}nr}, \qquad R_2(n_0) = \int_0^{n_0}\d{n}\; \frac{\mathcal{E}'(n)}{\sqrt{\pi}r},\\ r = \sqrt{\ln n_0 - \ln n}. \end{gather*}\]

The equations of motion follow by varying either \(\mathcal{E}\) or \(\mu\) with respect to the parameters \(\sigma^2\).

Analytic Examples: Power-law EoS

We can get quite far analytically if we assume a power-law equation of state

\[\begin{gather*} \mathcal{E}(n) = \frac{\alpha n^\beta}{\beta}\,,\qquad Q_1(n_0) = \frac{\alpha n_0^{\beta}}{\beta^2}\,, \qquad R_1(n_0) = \frac{\alpha n_0^{\beta}}{\beta}\,. \end{gather*}\]

This includes the common cases of the GPE (\(\beta = 2\), \(\alpha=g\)) and the UFG (\(\beta=5/3\)).

\[\begin{gather*} E_1 = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{\alpha n_0^{\beta-1}}{\beta^2},\\ \mu_1 = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{\alpha n_0^{\beta-1}}{\beta},\\ n_1 = \abs{\phi}^2, \qquad n_0 = \frac{n_1}{\pi\sigma_x\sigma_y}, \qquad \pdiff{n_0}{\sigma_{x,y}^2} = -\frac{n_0}{2\sigma_{x,y}^2}. \end{gather*}\]

We can unify these if we introduce the parameter

\[\begin{gather*} a = \begin{cases} \frac{1}{\beta-1} & \text{(minimize energy)}, \\ \frac{1}{\beta(\beta-1)} & \text{(minimize chemical potential)}, \end{cases} \end{gather*}\]

so that the objective function \(f \in \{E_1, \mu_1\}\) is

\[\begin{gather*} f = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{\alpha n_0^{\beta-1}}{a(\beta-1)\beta^2}\,. \end{gather*}\]

Varying this, we obtain the following coupled relations for \(\sigma_{x,y}^2\)

\[\begin{gather*} -\frac{\hbar^2}{2m\sigma_{x,y}^2 \lambda_{x,y}^2} + \frac{m\omega_{x,y}^2\sigma_{x,y}^2}{2} = \frac{\alpha n_0^{\beta-1}}{a\beta^2}\,, \qquad n_0 = \frac{n_1}{\pi\sigma_x\sigma_y}\,. \end{gather*}\]

These equations are coupled through \(n_0\), which makes solving them non-trivial. Once a solution is found, the energy and effective chemical potential can be written as

\[\begin{gather*} E_1 = (1 - a)\frac{\hbar^2}{4m}\left( \frac{1}{\sigma_x^2\lambda_x^2} + \frac{1}{\sigma_y^2\lambda_y^2}\right) + (1 + a)\frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4},\\ \mu_1 = (1 - \beta a)\frac{\hbar^2}{4m}\left( \frac{1}{\sigma_x^2\lambda_x^2} + \frac{1}{\sigma_y^2\lambda_y^2}\right) + (1 + \beta a)\frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4}. \end{gather*}\]

A simple analytic solution exists for the GPE with axial symmetry. In the GPE with \(\beta = 2\) we can optimize because \(a \in \{1, 1/2\}\) and one of the two terms vanish in each case. These are implemented in gpe.tube.StateGPEdrZ.

\[\begin{gather*} \mathcal{E}(n) = \frac{gn^2}{2}\,, \qquad \alpha = g\,, \qquad \beta = 2\,, \qquad \sigma_x = \sigma_y = \sigma,\\ \lambda^{-2} = \frac{\lambda_x^{-2} + \lambda_x^{-2}}{2}, \qquad \omega^2 = \frac{\omega_x^2 + \omega_y^2}{2}\,,\\ m\omega^2\sigma^2 = \sqrt{\frac{g m \omega^2 n_1}{2\pi a} + \frac{\hbar^2\omega^2}{\lambda^2}}\,,\\ \mu_1 = (1 - 2 a)\frac{\hbar^2}{2m\sigma^2\lambda^2} + (1 + 2 a)\frac{m\omega^2\sigma^2}{2}. \end{gather*}\]

If we follow [Mateo and Delgado, 2009] and minimize the chemical potential, then \(a=1/2\) and the first term cancels giving us the effective 1D chemical potential:

\[\begin{gather*} \mu_1 = \frac{\hbar\omega}{\lambda}\sqrt{1 +\frac{g m \lambda n_1}{\hbar^2 \pi}} = \frac{\hbar\omega}{\lambda}\sqrt{1 + 4 a_s \lambda n_1}\,. \end{gather*}\]

This agrees with Eq. (4) [Mateo and Delgado, 2009] if we have no scaling \(\lambda = 1\).

Remaining Challenges:

  • In the case of pancake traps, the expressions for \(Q_2(n_0)\) and \(R_2(n_0)\) are non-analytic, even for the GPE.

  • Without axial symmetry or outside of the GPE, the self-consistency equations have a more complicated form. Within the GPE it remains polynomial, but becomes a general quartic equation.

We hope to solve these cases as described in Solving for Sigmas: Globally Convergent Newton’s Method but have not completed this general analysis yet.

Two Components#

If we have two components, then we must have a different \(\sigma\) for each component and the resulting equations become:

\[\begin{gather*} \mathcal{E} = \sum_{i=a,b}\left\{ \frac{\hbar^2\abs{\nabla_z\psi_i}^2}{2m} + \Biggl( \overbrace{ \frac{\hbar^2\sigma^2_{i,z}(z, t)}{2m\sigma_i^2} }^{\text{neglected}} + \frac{\hbar^2}{2m\sigma_i^2} + \frac{m\omega_\perp^2\sigma_i^2)}{2} \Biggr)\psi_i^\dagger\psi_i \right\} + \frac{g_{aa}\abs{\psi_a}^4}{4\pi\sigma_a^2} + \frac{g_{bb}\abs{\psi_b}^4}{4\pi\sigma_b^2} + \frac{g_{ab}\abs{\psi_a}^2\abs{\psi_b}^2}{\pi(\sigma_a^2 + \sigma_b^2)},\\ \I\hbar \dot{\psi}_a = \left( \frac{-\hbar^2\nabla_z^2}{2m} + \frac{\hbar^2}{2m\sigma_a^2} + \frac{m\omega_\perp^2\sigma_a^2}{2} + \frac{g_{aa}\abs{\psi_a}^2}{2\pi\sigma_a^2} + \frac{g_{ab}\abs{\psi_b}^2}{\pi(\sigma_a^2 + \sigma_b^2)} \right)\psi_a,\\ \frac{m\omega_\perp^2\abs{\psi_a}^2}{2} = \frac{g_{aa}\abs{\psi_a}^4}{4\pi\sigma_a^4} +\frac{g_{ab}\abs{\psi_a}^2\abs{\psi_b}^2}{\pi(\sigma_a^2 + \sigma_b^2)^2} +\frac{\hbar^2\abs{\psi_a}^2}{2m\sigma_a^4}. \end{gather*}\]

Most of the terms are just an independent sum over the components, but we must be careful with the non-linear interactions which may now appear as

\[\begin{gather*} \frac{g_{aa}}{2}(\abs{\psi_a}^2\abs{\phi_a}^2)^2 + \frac{g_{b}}{2}(\abs{\psi_b}^2\abs{\phi_b}^2)^2 + g_{ab}\abs{\psi_a}^2\abs{\phi_a}^2\abs{\psi_b}^2\abs{\phi_b}^2,\\ \frac{g_{aa}\abs{\psi_a}^2}{4\pi \sigma_a^2}\psi_a^\dagger\psi_a + \frac{g_{bb}\abs{\psi_b}^2}{4\pi \sigma_b^2}\psi_b^\dagger\psi_b + g_{ab},\\ \int \d{x_\perp^2} \frac{e^{-x_\perp^2/2\sigma_a^2 - x_\perp^2/2\sigma_b^2}}{\pi\sigma_a\sigma_b} = \frac{1}{\pi(\sigma_a^2 + \sigma_b^2)} \end{gather*}\]

To solve for \(\sigma_i\) we rewrite the equations as follows:

\[\begin{gather*} \vec{F}(x) = \begin{pmatrix} (A_a - C_ax_c) x_a^2 - B_a\\ (A_b - C_bx_c) x_b^2 - B_b\\ x_c (x_a + x_b)^2 - 1 \end{pmatrix}, \\ \mat{J} = \begin{pmatrix} 2(A_a - C_ax_c) x_a & 0 & -C_ax_a^2\\ 0 & 2(A_b - C_bx_c) x_b & -C_bx_b^2\\ 2x_c (x_a + x_b) & 2x_c (x_a + x_b) & (x_a + x_b)^2 \end{pmatrix} \\ x_a = \sigma_a^2, \qquad x_c = \frac{1}{(\sigma_a^2+\sigma_b^2)^2}, \\ A_a = \frac{m\omega_\perp^2}{2}, \qquad C_a = \frac{g_{ab}\abs{\psi_b}^2}{\pi}, \qquad B_a = \frac{g_{aa}\abs{\psi_a}^2}{4\pi} + \frac{\hbar^2}{2m}. \end{gather*}\]
\[\begin{gather*} \frac{m\omega_\perp^2}{2}\sigma_a^4 = \frac{g_{aa}\abs{\psi_a}^2}{4\pi} +\frac{g_{ab}\abs{\psi_b}^2\sigma_a^4}{\pi(\sigma_a^2 + \sigma_b^2)^2} +\frac{\hbar^2}{2m}. \end{gather*}\]

Again, details for solving this are presented in Solving for Sigmas: Globally Convergent Newton’s Method.

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 \(\psi(z,t)\) which gives:

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

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

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

so that the Hamiltonian has the form

\[\begin{gather*} 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), \\ \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). \end{gather*}\]

The normalization condition for \(\phi\), i.e. \(\braket{1}_\phi = 1\), is independent of \(z\) allowing 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{gather*} \I\hbar\dot{\psi}(z,t) = \left( \frac{-\hbar^2\nabla_z^2}{2m} + \Braket{ - \I\hbar\partial_t - \frac{\hbar^2\nabla^2_\vect{x}}{2m} + f'\bigl(n(\vect{x}, t)\bigr)}_{\phi} \right)\psi(z, t),\\ \Braket{\op{A}}_\phi = \int \phi^\dagger(\vect{x}, t)\op{A}\phi(\vect{x}, t)\;\d{x}\d{y}, \end{gather*}\]

where the average \(\braket{}_{\phi}\) is over \(\phi(\vect{x}, t)\). For the NPSEQ, we will approximate this, but a complete extremization of the action defines \(\phi(\vect{x}, t)\) as the solution to the coupled radial equation:

\[\begin{gather*} \I\hbar\dot{\phi}(\vect{x}, t) = \left( \frac{-\hbar^2\nabla^2_\vect{x}}{2m} + f'\bigl(n(\vect{x}, t)\bigr)\right) \phi(\vect{x}, t). \end{gather*}\]

The same result can be obtained by expanding the action:

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

For the usual energy density, this becomes

\[\begin{gather*} S = \int\d{t} \d{z}\left\{ \I\hbar \psi^\dagger\dot{\psi} + \psi^\dagger\left(\frac{\hbar^2\nabla_z^2}{2m}\right)\psi + \overbrace{\psi^\dagger\psi}^{n_1(z,t)}\Braket{\I\hbar\partial_t + \frac{\hbar^2\nabla_\vect{x}^2}{2m} -\frac{f(n)}{n_1} }_{\phi} \right\} \end{gather*}\]

Pancake version

This is a straightforward generalization

\[\begin{multline*} S = \int\d{t}\d{x}\d{y}\Biggl\{ \I\hbar \psi_2^\dagger\dot{\psi}_2 + \psi_2^\dagger\left(\frac{\hbar^2(\nabla_{x}^2+\nabla_{y}^2)}{2m}\right)\psi_2 +\\ + \underbrace{\psi_2^\dagger\psi_2}_{n_2(x,y,t)} \Braket{\I\hbar\partial_t + \frac{\hbar^2\nabla_\vect{x}^2}{2m} -\frac{f(n)}{n_2} }_{\phi_2} \Biggr\} \end{multline*}\]

In what follows, we will use the notation \(\phi_2\) and \(\psi_2(x, y)\) for the pancake version, but keep \(\phi \equiv \phi_1\) and \(\psi \equiv \psi_1\) in the main text.

This is exact, but to proceed further we must make some assumptions about the form and behavior of \(\phi(\vect{x}, t)\) and \(f(n)\). Before proceeding, we define two auxiliary quantities that will play a key role in our minimization procedure:

\[\begin{align*} \mathcal{E}_1 &= \Braket{-\frac{\hbar^2\nabla_\vect{x}^2}{2m} + \frac{f(n)}{n_1}}_{\phi}, & \mu_{1} &= \Braket{- \frac{\hbar^2\nabla^2_\vect{x}}{2m} + f'(n)}_{\phi},\\ \end{align*}\]

Pancake version

\[\begin{align*} \mathcal{E}_2 &= \Braket{-\frac{\hbar^2\nabla_\vect{x}^2}{2m} + \frac{f(n)}{n_2}}_{\phi_2}, & \mu_{2} &= \Braket{- \frac{\hbar^2\nabla^2_\vect{x}}{2m} + f'(n)}_{\phi_2},\\ \end{align*}\]

Adiabatic Approximation#

The adiabatic approximation is that the variations of \(\psi(z, t)\) are slow so that the radial equation adjusts quickly to become stationary states of the radial equation. Recall that we normalize our states so that

\[\begin{gather*} n(\vect{z};t) = \abs{\phi}^2(\vect{x},t)\overbrace{\abs{\psi}^2(z, t)}^{n_1(z,t)},\\ n_1(z,t) = \int\d^2{\vect{r}_{\perp}}\; n(\vect{x};t) = \overbrace{\int\d^2{\vect{r}_\perp}\, \abs{\phi}^2(\vect{x},t)}^{1}\; \abs{\psi}^2(z,t) = \abs{\psi}^2(z,t). \end{gather*}\]

Here \(n_1(z,t)\) is the integrated density along the cloud. The adiabatic approximation is that

\[\begin{gather*} \phi(\vect{x}, t) = \phi\bigl(\vect{r}_\perp; n_1(z,t)\bigr)e^{\I\mu_1 t/\hbar}. \end{gather*}\]

I.e., up to a phase, all of the time dependence and \(z\) dependence come through the integrated density \(n_1(z, t)\). As part of the adiabatic approximation, we neglect spatial (along \(z\)) and temporal variations of \(n_1\) in the equation for \(\phi(\vect{r}_\perp; n_1)\), which “instantaneously” adjusts to satisfy the radial equation for the given axial density. This becomes a problem in some cases, such as the common case of suddenly turning off the external trap to let a cloud expand. In this case, the radial function will suddenly expand to fill all space, which is a very poor approximation. The Dynamically Rescaled GPE (dr-GPE) fixes this by transforming to an expanding frame where there is an effective trapping potential (from the expanding frame), even when the trap is suddenly turned off.

Variational Approach#

Note that the equation for \(\phi(\vect{r}_\perp; n_1)\) follows from minimizing the following energy functional

\[\begin{gather*} \mathcal{E}_1[\phi] = \frac{E_1[\phi]}{n_1} = \int\d^2{\vect{r}_\perp}\left( \frac{\hbar^2\abs{\vect{\nabla}_\perp\phi}^2}{2m} + V_\perp(\vect{r}_\perp)\abs{\phi}^2 + \frac{\mathcal{E}\bigl(n_1 \abs{\phi}^2\bigr)}{n_1} \right). \end{gather*}\]

In general, we cannot solve this radial equation analytically, so we typically make an approximation on the form of \(\phi\). The variational approach is to stay within a restricted set of functional forms, but now to minimize with respect to a few parameters, such as the widths of a gaussian \(\sigma_{x,y}\). [Mateo and Delgado, 2009] suggests two alternatives, minimizing the quantities we defined above:

  • Minimizing the energy per unit length:

    \[\begin{gather*} \mathcal{E}_1[\phi] %= \frac{E_1[\phi]}{n_1} = \overbrace{\int \d^2\vect{r}_{\perp}\;\left( \frac{\hbar^2\abs{\nabla_\perp \phi}^2}{2m} + V_\perp(\vect{r}_\perp)\abs{\phi}^2 + \frac{\mathcal{E}\bigl(n_1\abs{\phi}^2\bigr)}{n_1} \right)}^{\Braket{-\frac{\hbar^2}{2m}\nabla_\vect{x}^2 + f(n)/n_1}_{\phi}}. \end{gather*}\]

  • Minimizing the chemical potential:

    \[\begin{gather*} \mu_1[\phi] = \pdiff{\mathcal{E}_1[\phi]}{n_1} = \overbrace{\int \d^2\vect{r}_{\perp}\; \phi^*\left( \frac{-\hbar^2\nabla_\perp^2}{2m} + V_\perp(\vect{r}_\perp) + \mathcal{E}'\bigl(n_1\abs{\phi}^2\bigr) \right)\phi}^{\Braket{-\frac{\hbar^2}{2m}\nabla^2_\vect{x} + f'(n)}_{\phi}}. \end{gather*}\]

They argue that, while minimization of \(\mathcal{E}_1[\phi]\) always yields the correct result when the minimum is taken over the whole space of functions \(\phi\), the second approach of minimizing \(\mu_1[\phi]\) can yield better approximation in the case of restricted forms of \(\phi\), and that these two approaches are equivalent when the non-linear interactions vanish.

Pancake version

Everything is similar for reduction along only \(z\) to produce effective 2D versions for pancake”-shaped clouds.

  • Minimizing the energy per unit area:

    \[\begin{gather*} \mathcal{E}_2[\phi_2] = \int_{-\infty}^{\infty}\d{z}\;\left( \frac{\hbar^2\abs{\nabla_z \phi_2}^2}{2m} + V_z(z)\abs{\phi_2}^2 + \frac{\mathcal{E}\bigl(n_2\abs{\phi_2}^2\bigr)}{n_2} \right). \end{gather*}\]

  • Minimizing the chemical potential:

    \[\begin{gather*} \mu_2[\phi_2] = \pdiff{\mathcal{E}_2[\phi_2]}{n_2} = \int_{-\infty}^{\infty}\d{z}\; \phi_2^*\left( \frac{-\hbar^2\nabla_z^2}{2m} + V_z(z) + \mathcal{E}'\bigl(n_2\abs{\phi_2}^2\bigr) \right)\phi_2. \end{gather*}\]

Gaussian Approximation#

As mentioned above, a common approach is to use a gaussian:

\[\begin{gather*} n = n_1\abs{\phi_1(\vect{r}_\perp;n_1)}^2 = \overbrace{\frac{n_1}{\pi\sigma_x\sigma_y}}^{n_0} \exp\left(-\frac{x^2}{\sigma_x^2} - \frac{y^2}{\sigma_y^2}\right) \end{gather*}\]

where \(\sigma_{x,y}(n_1)\) are functions of \(n_1\). Here we also introduce the central density \(n_0\) which will play an important role below.

Pancake formula

For the effective 2D theory, the transverse function is a single gaussian:

\[\begin{gather*} n = n_2\abs{\phi_2(z; n_2)}^2 = \overbrace{\frac{n_2}{\sqrt{\pi}\sigma}}^{n_0} \exp\left(-\frac{z^2}{\sigma}\right) \end{gather*}\]

where \(\sigma(n_2)\) is a function of \(n_2\).

Depending on whether we choose to minimize the energy or the chemical potential, we must be able to compute \(R(n_0)\) or \(Q(n_0)\):

\[\begin{gather*} n = n_1\abs{\phi}^2 = n_0e^{-r^2} = n_0e^{-x^2/\sigma_x^2 -y^2/\sigma_y^2},\qquad n_0 = \frac{n_1}{\pi\sigma_x\sigma_y},\\ n_1\mathcal{E}_1 = \cdots + \int \d^2\vect{r}_{\perp}\; \mathcal{E}\bigl(n_1\abs{\phi}^2\bigr) = \cdots + \pi \sigma_x\sigma_y Q_1(n_0),\\ n_1\mu_1 = \cdots + \int \d^2\vect{r}_{\perp}\; n_1\abs{\phi}^2 \mathcal{E}'\bigl(n_1\abs{\phi}^2\bigr) = \cdots + \pi \sigma_x\sigma_y R_1(n_0),\\ %Q_1'(n) = \frac{\mathcal{E}(n)}{n}, \qquad %R_1'(n) = \mathcal{E}'(n),\\ Q_1(n_0) = \int_0^{n_0}\d{n}\; \frac{\mathcal{E}(n)}{n}, \qquad R_1(n_0) = \int_0^{n_0}\d{n}\; \mathcal{E}'(n) = \mathcal{E}(n_0) - \mathcal{E}(0). \end{gather*}\]

Depending on the form of the equation of state, \(Q_1(n_0)\) may also have an analytic expression.

Pancake Formulae

The analogous formula for pancake clouds (2D) are not so simple due to the extra factor of \(r\) in the denominator of the integrands:

\[\begin{gather*} n = n_2\abs{\phi}^2 = n_0e^{-r^2} = n_0e^{-z^2/\sigma^2},\qquad n_0 = \frac{n_2}{\sqrt{\pi}\sigma},\\ n_2\mathcal{E}_2 = \cdots + \int_{-\infty}^{\infty}\d{z}\; \mathcal{E}\bigl(n_2\abs{\phi}^2\bigr) = \cdots + \sqrt{\pi}\sigma Q_2(n_0),\\ n_2\mu_2 = \cdots + \int_{-\infty}^{\infty}\d{z}\; n_2\abs{\phi}^2 \mathcal{E}'\bigl(n_2\abs{\phi}^2\bigr) = \cdots + \sqrt{\pi}\sigma R_2(n_0),\\ %Q_2'(n) = \frac{\mathcal{E}(n)}{rn\sqrt{\pi}}, \qquad %R_2'(n) = \frac{\mathcal{E}'(n)}{r\sqrt{\pi}},\\ Q_2(n_0) = \int_0^{n_0}\d{n}\; \frac{\mathcal{E}(n)}{\sqrt{\pi}nr}, \qquad R_2(n_0) = \int_0^{n_0}\d{n}\; \frac{\mathcal{E}'(n)}{\sqrt{\pi}r},\\ r = \sqrt{\ln n_0 - \ln n}. \end{gather*}\]

Integration Details

For gaussian radial functions, we can compute these integrals with a change of variables and rescaling:

\[\begin{gather*} n \equiv n_1\abs{\phi}^2 = n_0 e^{-r^2}, \qquad r^2 = \frac{x^2}{\sigma_x^2} + \frac{y^2}{\sigma_y^2}, \qquad n_0 = \frac{n_1}{\pi\sigma_x\sigma_y},\\ \d{n} = -2nr\d{r},\\ \int \d^2\vect{r}_{\perp}\; \mathcal{E}\bigl(n_1\abs{\phi}^2\bigr) = \pi \sigma_x\sigma_y \int_{0}^{\infty} 2r\d{r}\; \mathcal{E}(n) = \pi \sigma_x\sigma_y \int_{0}^{n_0} \d{n}\; \frac{\mathcal{E}(n)}{n} = \pi \sigma_x\sigma_y Q_1(n_0),\\ \int \d^2\vect{r}_{\perp}\; n\mathcal{E}'(n) = \pi \sigma_x\sigma_y \int_{0}^{n_0} \d{n}\; \mathcal{E}'(n) = \pi \sigma_x\sigma_y \Bigl(\mathcal{E}(n_0) - \mathcal{E}(0)\Bigr) = \pi \sigma_x\sigma_y R_1(n_0). \end{gather*}\]

In , we will need the following version of these:

\[\begin{gather*} n \equiv n_2\abs{\phi}^2 = n_0 e^{-r^2}, \qquad r^2 = \frac{z^2}{\sigma^2}, \qquad n_0 = \frac{n_2}{\sqrt{\pi}\sigma},\\ r = \sqrt{-\ln n/n_0}, \qquad \sigma\d{r} = \d{z}, \qquad \d{n} = -2nr\d{r},\\ \int \d{z}\; \mathcal{E}\bigl(n_2\abs{\phi}^2\bigr) = 2 \sigma \int_{0}^{\infty} \d{r}\; \mathcal{E}(n) = \sigma \int_{0}^{n_0} \d{n}\; \frac{\mathcal{E}(n)}{n\sqrt{-\ln n/n_0}} = \sqrt{\pi}\sigma Q_2(n_0),\\ \int \d{z}\; n\mathcal{E}'(n) = \sigma \int_{0}^{n_0} \d{n}\; \frac{\mathcal{E}'(n)}{\sqrt{-\ln n/n_0}} = \sqrt{\pi}\sigma Q_2(n_0). \end{gather*}\]

These are not so nice.

At the end of the day we need the following integrals

\[\begin{gather*} n = n_0e^{-r^2}, \qquad r = \sqrt{\ln n_0 - \ln n},\\ Q_{1}(n_0) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{n}, \qquad R_{1}(n_0) = \int_0^{n_0}\d{n}\;\mathcal{E}'(n), \\ \sqrt{\pi}Q_{2}(n_0) = \int_0^{\infty}2\d{r}\;\mathcal{E}(n) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{nr} = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{n\sqrt{\ln n_0 - \ln n}}, \\ \sqrt{\pi}R_{2}(n_0) = \int_0^{\infty}2\d{r}\;n\mathcal{E}'(n) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}'(n)}{r} = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}'(n)}{\sqrt{\ln n_0 - \ln n}}. \end{gather*}\]

These can be computed for power-law EoS

\[\begin{gather*} \mathcal{E}(n) = \alpha \frac{n^\beta}{\beta} = \frac{\alpha n_0^\beta e^{-\beta r^2}}{\beta}, \qquad n\mathcal{E}'(n) = \alpha n^\beta = \alpha n_0^\beta e^{-\beta r^2},\\ Q_{1}(n_0) = \frac{\alpha n_0^{\beta}}{\beta^2}, \qquad R_{1}(n_0) = \frac{\alpha n_0^{\beta}}{\beta} = \beta Q_1(n_0), \\ \sqrt{\pi}Q_{2}(n_0) = 2\int_0^{\infty}\d{r}\;\mathcal{E}(n), \qquad \sqrt{\pi}R_{2}(n_0) = 2\int_0^{\infty}\d{r}\;n\mathcal{E}'(n) \end{gather*}\]

These follow from

\[\begin{gather*} \sqrt{\pi}Q_{2}(n_0) = 2\int_0^{\infty}\d{r}\;\mathcal{E}(n) = 2\frac{\alpha n_0^\beta}{\beta} \overbrace{\int_0^{\infty}\d{r}\;e^{-\beta r^2}}^{\sqrt{\pi}/2\sqrt{\beta}} = \frac{\sqrt{\pi} \alpha n_0^\beta}{\beta^{3/2}},\\ \sqrt{\pi}R_{2}(n_0) = 2\int_0^{\infty}\d{r}\;n\mathcal{E}'(n) = 2 \alpha n_0^{\beta}\int_0^{\infty}\d{r}\;e^{-\beta r^2} = \frac{\sqrt{\pi}\alpha n_0^{\beta}}{\sqrt{\beta}} = \beta \sqrt{\pi}Q_{2}(n_0). \end{gather*}\]

Hide code cell content

 
# Check these relations

import numpy as np
from scipy.integrate import quad

alpha = 1.23
beta = 2.34
n0 = 1.11
s = 2.11

def E(n, d=0):
    if d == 0:
        return alpha * n**beta / beta
    elif d == 1:
        return alpha * n**(beta - 1)


def n1(r):
    return n0*np.exp(-r**2)

def Q1_r(r):
    n = n1(r)
    return 2*np.pi * r * E(n,d=0) / np.pi

def Q1_n(n):
    return E(n,d=0)/n

def R1_r(r):
    n = n1(r)
    return 2*np.pi * r * n * E(n,d=1) / np.pi

def R1_n(n):
    return E(n,d=1)

def n2(z):
    return n0*np.exp(-z**2/s**2)

def Q2_z(z):
    n = n2(z)
    return E(n,d=0) / s / np.sqrt(np.pi)

def Q2_n(n):
    r = np.sqrt(np.log(n0) - np.log(n))
    return E(n,d=0)/n/r/np.sqrt(np.pi)

def R2_z(z):
    n = n2(z)
    return n*E(n,d=1) / s / np.sqrt(np.pi)

def R2_n(n):
    r = np.sqrt(np.log(n0) - np.log(n))
    return E(n,d=1)/r/np.sqrt(np.pi)

Q1 =  alpha*n0**beta / beta**2
R1 = beta * Q1
Q2 = alpha*n0**beta / beta**(3/2)
R2 = beta * Q2

ans_z = quad(Q1_r, 0, np.inf)[0]
ans_n = quad(Q1_n, 0, n0)[0]
assert np.allclose([ans_z, ans_n], Q1)

ans_z = quad(R1_r, 0, np.inf)[0]
ans_n = quad(R1_n, 0, n0)[0]
assert np.allclose([ans_z, ans_n], R1)

ans_z = quad(Q2_z, -np.inf, np.inf)[0]
ans_n = quad(Q2_n, 0, n0)[0]
assert np.allclose([ans_z, ans_n], Q2)

ans_z = quad(R2_z, -np.inf, np.inf)[0]
ans_n = quad(R2_n, 0, n0)[0]
assert np.allclose([ans_z, ans_n], R2)

# Now check with general EoS
def E(n, d=0):
    if d == 0:
        return n**2*np.exp(n)
    elif d == 1:
        return (2*n + n**2)*np.exp(n)

ans_z = quad(Q1_r, 0, np.inf)[0]
ans_n = quad(Q1_n, 0, n0)[0]
assert np.allclose(ans_z, ans_n)

ans_z = quad(R1_r, 0, np.inf)[0]
ans_n = quad(R1_n, 0, n0)[0]
assert np.allclose(ans_z, ans_n)

ans_z = quad(Q2_z, -np.inf, np.inf)[0]
ans_n = quad(Q2_n, 0, n0)[0]
assert np.allclose(ans_z, ans_n)

ans_z = quad(R2_z, -np.inf, np.inf)[0]
ans_n = quad(R2_n, 0, n0)[0]
assert np.allclose(ans_z, ans_n)

In summary, including the definition of \(n_1\) and \(n_2\) in terms of \(n_0\), the non-linear pieces have a nice form:

\[\begin{gather*} \mathcal{E}_1 = \cdots + \frac{Q_1(n_0)}{n_0}, \qquad \mu_1 = \cdots + \frac{R_1(n_0)}{n_0},\\ \mathcal{E}_2 = \cdots + \frac{Q_2(n_0)}{n_0}, \qquad \mu_2 = \cdots + \frac{R_1(n_0)}{n_0}. \end{gather*}\]

Harmonic Potential#

The remaining portion of these equations can be computed if we assume that the transverse potential is harmonic:

\[\begin{gather*} V_{\perp}(\vect{r}_\perp) = \frac{m\omega_x^2 x^2 + m\omega_y^2y^2}{2},\qquad \text{or for pancakes}\qquad V_{z}(z) = \frac{m\omega_z^2 z^2}{2}. \end{gather*}\]

With these we can compute all of the remaining expectation values.

Expectation Values

Explicitly, here are the various expectation values:

\[\begin{gather*} \Braket{\I\hbar \partial_t}_{\phi} = -\mu_1, \qquad \Braket{\nabla^2_{x,y}}_{\phi} = \frac{-1}{2\sigma_{x,y}^2}, \qquad \Braket{x^2}_{\phi} = \frac{\sigma_{x,y}^2}{2}\\ \Braket{\nabla^2_{z}}_{\phi} = -\frac{(\nabla_z\sigma_{x})^2}{2\sigma_{x}^2} -\frac{(\nabla_z\sigma_{y})^2}{2\sigma_{y}^2} \qquad \text{(neglected in our approx.)},\\ \Braket{\abs{\phi}^2} = \frac{1}{2\pi\sigma_x\sigma_y}, \qquad \Braket{\abs{\phi}^\nu} = \frac{2}{(2+\nu)(\sqrt{\pi\sigma_x\sigma_y})^{\nu}}. \end{gather*}\]

The corresponding pancake expectation values are:

\[\begin{gather*} \Braket{\I\hbar \partial_t}_{\phi} = -\mu_2, \qquad \Braket{\nabla^2_{z}}_{\phi} = \frac{-1}{2\sigma^2}, \qquad \Braket{z^2}_{\phi} = \frac{\sigma^2}{2}\\ \Braket{\nabla^2_{x,y}}_{\phi} = -\frac{(\nabla_{x,y}\sigma)^2}{2\sigma^2} \text{(neglected in our approx.)},\\ \Braket{\abs{\phi}^2} = \frac{1}{\sqrt{2\pi}\sigma}, \qquad \Braket{\abs{\phi}^\nu} = \frac{\sqrt{2}}{\sqrt{2+\nu}(\sqrt{\pi}\sigma)^{\nu/2}}. \end{gather*}\]

Hide code cell source

 
import sympy;sympy.init_printing()
from sympy import exp, sqrt, pi, oo, S
x, y, z = sympy.var('x, y, z', real=True)
sigma_x, sigma_y = sympy.var('sigma_x, sigma_y', positive=True)
dsigma_x, dsigma_y = sympy.var("sigma'_x, sigma'_y", positive=True)
ddsigma_x, ddsigma_y = sympy.var("sigma''_x, sigma''_y", positive=True)
sx, sy = sympy.Function('sigma_x'), sympy.Function('sigma_y')
phi = exp(-x**2/2/sx(z)**2-y**2/2/sy(z)**2)/sqrt(pi*sx(z)*sy(z))
subs = {sx(z): sigma_x, sy(z): sigma_y, 
        sx(z).diff(z): dsigma_x, sy(z).diff(z): dsigma_y,
        sx(z).diff(z,z): ddsigma_x, sy(z).diff(z,z): ddsigma_y}

def braket(f):
    return (phi*f(phi)).subs(subs).integrate((x, -oo, oo), (y, -oo, oo)).simplify()

assert braket(lambda phi:phi) == 1
assert braket(lambda phi:x**2*phi) == sigma_x**2/2
assert braket(lambda phi:y**2*phi) == sigma_y**2/2
assert braket(lambda phi:phi.diff(x, x)) == -1/sigma_x**2/2
assert braket(lambda phi:phi.diff(y, y)) == -1/sigma_y**2/2
assert braket(lambda phi:phi.diff(z, z)) == (-dsigma_x**2/sigma_x**2/2 - dsigma_y**2/sigma_y**2/2)

x, y, z = sympy.var('x, y, z', real=True)
sigma = sympy.var('sigma', positive=True)
dsigma_x, dsigma_y = sympy.var("sigma'_x, sigma'_y", positive=True)
ddsigma_x, ddsigma_y = sympy.var("sigma''_x, sigma''_y", positive=True)
s = sympy.Function('sigma')
phi1 = exp(-z**2/2/s(x,y)**2)/sqrt(sqrt(pi)*s(x,y))
subs = {s(x,y): sigma, 
        s(x,y).diff(x): dsigma_x, s(x,y).diff(y): dsigma_y}

def braket(f):
    return (phi1*f(phi1)).subs(subs).integrate((z, -oo, oo)).simplify()

assert braket(lambda phi1:phi1) == 1
assert braket(lambda phi1:z**2*phi1) == sigma**2/2
assert braket(lambda phi1:phi1.diff(x, x)) == -dsigma_x**2/2/sigma**2
assert braket(lambda phi1:phi1.diff(y, y)) == -dsigma_y**2/2/sigma**2
assert braket(lambda phi1:phi1.diff(z, z)) == -1/sigma**2/2
assert braket(lambda phi1:phi1**2*phi1) == sqrt(2/pi)/sigma/2
for nu in map(S, [1, 2, 3, 4, 5]):
    assert braket(lambda phi1:phi1**nu*phi1) == sqrt(2*nu + 4)/(nu+2)/(pi*sigma**2)**(nu/4)

Putting these together, we have

\[\begin{gather*} E_1 = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{Q_1(n_0)}{n_0},\\ \mu_1 = \frac{\hbar^2}{2m}\left( \frac{1}{2\sigma_x^2 \lambda_x^2} + \frac{1}{2\sigma_y^2 \lambda_y^2} \right) + \frac{m}{2}\frac{\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2}{2} + \frac{R_1(n_0)}{n_0},\\ n_1 = \abs{\phi}^2, \qquad n_0 = \frac{n_1}{\pi\sigma_x\sigma_y}, \qquad \pdiff{n_0}{\sigma_{x,y}^2} = -\frac{n_0}{2\sigma_{x,y}^2}. \end{gather*}\]

Pancake Version

\[\begin{gather*} E_2 = \frac{\hbar^2}{2m}\frac{1}{2\lambda_z^2\sigma^2} + \frac{m}{2}\frac{\omega_z^2\sigma^2}{2} + \frac{Q_2(n_0)}{n_0},\\ \mu_2 = \frac{\hbar^2}{2m}\frac{1}{2\lambda_z^2\sigma^2} + \frac{m}{2}\frac{\omega_z^2\sigma^2}{2} + \frac{R_2(n_0)}{n_0},\\ n_2 = \abs{\phi_2}^2, \qquad n_0 = \frac{n_2}{\sqrt{\pi}\sigma}, \qquad \pdiff{n_0}{\sigma^2} = -\frac{n_0}{2\sigma^2},\\ \end{gather*}\]

Note: in the kinetic energy, we have allowed for different scaling factors \(\lambda_x\) and \(\lambda_y\) in each direction. This will help with the book-keeping needed for the Dynamically Rescaled GPE (dr-GPE) later. Additional factors and the effective potential induced by the co-scaling frame (fictitious forces) can be absorbed into the effective trapping frequencies \(\omega_{x,y}\).

Varying these with respect to \(\sigma_x\) and \(\sigma_y\), we obtain the following equations:

\[\begin{gather*} \overbrace{\frac{m\omega_{x,y}^2}{2}}^{B_{x,y}}\sigma_{x,y}^2 - \overbrace{\frac{\hbar^2}{2m\lambda_{x,y}^2}}^{A_{x,y}}\frac{1}{\sigma_{x,y}^2} = Q_1'(n_0) - \frac{Q_1(n_0)}{n_0} \tag{energy}\\ \frac{m\omega_{x,y}^2\sigma_{x,y}^2}{2} - \frac{\hbar^2}{2m\lambda_{x,y}^2}\frac{1}{\sigma_{x,y}^2} = R'_1(n_0) - \frac{R_1(n_0)}{n_0}. \tag{chemical potential} \end{gather*}\]

Pancake Version

\[\begin{gather*} \overbrace{\frac{m\omega_{z}^2}{2}}^{B}\sigma^2 - \overbrace{\frac{\hbar^2}{2m\lambda_{z}^2}}^{A}\frac{1}{\sigma^2} = Q_2'(n_0) - \frac{Q_2(n_0)}{n_0} \tag{energy}\\ \frac{m\omega_{z}^2\sigma^2}{2} - \frac{\hbar^2}{2m\lambda_{z}^2}\frac{1}{\sigma^2} = R'_2(n_0) - \frac{R_2(n_0)}{n_0}. \tag{chemical potential} \end{gather*}\]

Solving these is a bit tricky, and we present details in Solving for Sigmas: Globally Convergent Newton’s Method.

Power-law Equation of State

Things simplify somewhat for power-law EoS:

\[\begin{gather*} \mathcal{E}(n) = \alpha \frac{n^{\beta}}{\beta}, \quad R_{i}(n_0) = \beta Q_{i}(n_0) \propto n_0^{\beta}. \end{gather*}\]

Explicitly:

\[\begin{gather*} R_{1}(n_0) = \beta Q_{1}(n_0) = \alpha\frac{n_0^{\beta}}{\beta}, \qquad n_0 = \frac{n_1}{\pi\sigma_{x}\sigma_{y}},\\ R_{2}(n_0) = \beta Q_{2}(n_0) = \alpha\frac{n_0^{\beta}}{\sqrt{\beta}}, \qquad n_0 = \frac{n_2}{\sqrt{\pi}\sigma}. \end{gather*}\]

Old Stuff: Needs cleaning

Power-law EoS (Old)

For example, consider a pow-law equation of state, and a separable external potential:

\[\begin{gather*} \mathcal{E}(n)= \alpha \frac{n^\beta}{\beta}, \qquad Q(n) = \int_{0}^{n}\d{n}\frac{\mathcal{E}(n)}{n} = \frac{\alpha n^{\beta}}{\beta^2}. \end{gather*}\]

This includes the form of the usual GPE when \(\alpha = g\) and \(\beta = 2\). We then have

\[\begin{gather*} n_1\mathcal{E}_1 = \cdots + \pi \sigma_x\sigma_y \frac{\alpha n_0^\beta}{\beta^2} = \cdots + \frac{\alpha n_1^\beta} {\beta^2(\pi \sigma_x\sigma_y)^{\beta - 1}},\\ n_1\mathcal{\mu}_\perp = \cdots + \pi\sigma_x\sigma_y \alpha \frac{n_0^{\beta}}{\beta} = \cdots + \frac{\alpha n_1^{\beta}} {\beta (\pi \sigma_x\sigma_y)^{\beta - 1}} \end{gather*}\]

This allows us to write the equations cleanly separating the radial from axial components:

\[\begin{gather*} \mu_1 = \Braket{ \frac{-\hbar^2\nabla^2_\perp}{2m} + V_\perp(\vect{r}_\perp) + \alpha \abs{\phi}^{2(\beta-1)}n_1^{\beta-1} }_{\phi} \end{gather*}\]
\[\begin{gather*} n_1(z,t) = \int\d^2{\vect{r}_{\perp}}\; n(\vect{x};t) = \overbrace{\int\d^2{\vect{r}_\perp}\, \abs{\phi}^2(\vect{x},t)}^{1}\; \abs{\psi}^2(z,t) = \abs{\psi}^2(z,t),\\ \I\hbar\dot{\psi}(z,t) = \left( \frac{-\hbar^2\nabla_z^2}{2m} + V_z(z) + \mu_1[n_1(z, t)] \right)\psi(z, t),\\ \mu_1[n_1] \phi(\vect{r}_\perp; n_1) = \left( \frac{-\hbar^2\nabla^2_\perp}{2m} + V_\perp(\vect{r}_\perp) + \alpha \abs{\phi}^{2(\beta-1)}n_1^{\beta-1}\right) \phi(\vect{r}_\perp; n_1),\\ \mu_1[n_1] = \Braket{ \frac{-\hbar^2\nabla^2_\perp}{2m} + V_\perp(\vect{r}_\perp) + \alpha \abs{\phi}^{2(\beta-1)}n_1^{\beta-1} }_\phi. \end{gather*}\]

Variational Approach#

To solve for \(\phi(r_\perp, n_1)\) exactly is still challenging, so one typically restricts \(\phi(\vect{r}_{\perp}, n_1)\) to some restricted subspace of functions, minimizing an appropriate quantity. [Mateo and Delgado, 2009] suggests two alternatives:

  • Minimizing the energy per unit length:

    \[\begin{gather*} \mathcal{E}_1[\phi] = \frac{E_1[\phi]}{n_1} = \int \d^2\vect{r}_{\perp}\;\left( \frac{\hbar^2\abs{\nabla_\perp \phi}^2}{2m} + V_\perp(\vect{r}_\perp)\abs{\phi}^2 + \frac{\mathcal{E}(n_1\abs{\phi}^2)}{n_1} \right)\\ = \int \d^2\vect{r}_{\perp}\;\left( \frac{\hbar^2\abs{\nabla_\perp \phi}^2}{2m} + V_\perp(\vect{r}_\perp)\abs{\phi}^2 + \frac{Q_1(n_0)}{n_1} \right). \end{gather*}\]

  • Minimizing the chemical potential:

    \[\begin{gather*} \mu_1[\phi] = \pdiff{E[\phi]}{n_1} = \int \d^2\vect{r}_{\perp}\; \phi^*\left( \frac{-\hbar^2\nabla_\perp^2}{2m} + V_\perp(\vect{r}_\perp) + \mathcal{E}'(n_1\abs{\phi}^2) \right)\phi. \end{gather*}\]

Power-law forms (old)

  • Minimizing the energy per unit length:

    \[\begin{gather*} \mathcal{E}_1[\phi] = \frac{E_1[\phi]}{n_1} = \int \d^2\vect{r}_{\perp}\;\left( \frac{\hbar^2\abs{\nabla_\perp \phi}^2}{2m} + V_\perp(\vect{r}_\perp)\abs{\phi}^2 + \alpha\frac{\abs{\phi}^{2\beta}}{\beta}n_1^{\beta-1} \right). \end{gather*}\]

  • Minimizing the chemical potential:

    \[\begin{gather*} \mu_1[\phi] = \pdiff{E[\phi]}{n_1} = \int \d^2\vect{r}_{\perp}\; \phi^*\left( \frac{-\hbar^2\nabla_\perp^2}{2m} + V_\perp(\vect{r}_\perp) + \alpha\abs{\phi}^{2(\beta-1)}n_1^{\beta-1} \right)\phi. \end{gather*}\]

They argue that, while minimization of \(\mathcal{E}_1[\phi]\) always yields the correct result when the minimum is taken over the whole space of functions \(\phi\), the second approach of minimizing \(\mu_1[\phi]\) can yield better approximation in the case of restricted forms of \(\phi\), and that these two approaches are equivalent when \(\alpha \rightarrow 0\).

HO Potential and Gaussian Wavefunction#

Here we derive the results for harmonic oscillator potentials with the Gaussian approximation:

\[\begin{gather*} V_{\perp}(\vect{r}_\perp) = \frac{m\omega_x^2x^2 + m\omega_y^2y^2}{2}, \qquad \phi\bigl(x, y; \sigma(z, t)\bigr) = \frac{e^{\I\mu t/\hbar}e^{-x^2/2\sigma_x^2 - y^2/2\sigma_y^2}} {\sqrt{\pi\sigma_x(z, t)\sigma_y(z, t)}}. \end{gather*}\]

The expectations values are:

\[\begin{gather*} \Braket{ \frac{-\hbar^2}{2m}\frac{\nabla_x^2\phi/\lambda_x^{2} + \nabla_y^2\phi/\lambda_y^{2}}{\phi} }_{\phi} = \frac{\hbar^2}{2m}\left( \frac{\sigma_x^2 - x^2}{\lambda_x^2\sigma_x^4} + \frac{\sigma_y^2 - y^2}{\lambda_y^2\sigma_y^4} \right)\\ \Braket{\frac{\hbar^2}{2m\lambda_x^2}\nabla_{x}^2}_{\phi} = -\frac{\hbar^2}{4\lambda_x^2\sigma_x^2}, \qquad \Braket{\frac{\hbar^2}{2m\lambda_y^2}\nabla_{y}^2}_{\phi} = -\frac{\hbar^2}{4\lambda_y^2\sigma_y^2},\\ \Braket{\abs{\phi}^2} = \frac{1}{2\pi\sigma_x\sigma_y}, \qquad \Braket{\abs{\phi}^\nu} = \frac{2}{(2+\nu)(\sqrt{\pi\sigma_x\sigma_y})^{\nu}}, \\ \Braket{x^2} = \frac{\sigma_x^2}{2}, \qquad \Braket{y^2} = \frac{\sigma_y^2}{2}, \\ \Braket{\nabla_z^2}_\phi = -\frac{(\sigma'_x)^2}{2\sigma_x^2} -\frac{(\sigma'_y)^2}{2\sigma_y^2}, \qquad \Braket{\I\hbar\partial_t} = -\mu. \end{gather*}\]

Note: in the kinetic energy, we have allowed for different scaling factors in each direction. This will help with the book-keeping needed for the Dynamically Rescaled GPE (dr-GPE) later.

\[\begin{gather*} E_1[\phi] = \frac{\hbar^2}{4m}\left( \frac{1}{\lambda_x^2\sigma_x^2} + \frac{1}{\lambda_y^2\sigma_y^2} \right) + \frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4} + \overbrace{Q(n_0)/n_0} ^{\frac{\alpha}{\beta^2}\frac{n_1^{\beta-1}}{(\pi\sigma_x\sigma_y)^{\beta - 1}}},\\ \mu_1[\phi] = \frac{\hbar^2}{4m}\left( \frac{1}{\lambda_x^2\sigma_x^2} + \frac{1}{\lambda_y^2\sigma_y^2} \right) + \frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4} + \underbrace{\mathcal{E}(n_0)/n_0}_{\frac{\alpha}{\beta}\frac{n_1^{\beta-1}}{(\pi\sigma_x\sigma_y)^{\beta-1}}},\\ n_0 = \frac{n_1}{\pi \sigma_x\sigma_y}, \qquad \pdiff{n_0}{\sigma_{x,y}^2} = \frac{-n_0}{2\sigma^2_{x,y}}. \end{gather*}\]

Varying these with respect to \(\sigma_x\) and \(\sigma_y\), we obtain the following equations:

\[\begin{gather*} \overbrace{\frac{m\omega_{x,y}^2}{2}}^{B_{x,y}}\sigma_{x,y}^2 - \overbrace{\frac{\hbar^2}{2m\lambda_{x,y}^2}}^{A_{x,y}}\frac{1}{\sigma_{x,y}^2} = Q'(n_0) - \frac{Q(n_0)}{n_0}, \tag{energy} \\ \frac{m\omega_{x,y}^2\sigma_{x,y}^2}{2} - \frac{\hbar^2}{2m\lambda_{x,y}^2}\frac{1}{\sigma_{x,y}^2} = \mathcal{E}'(n_0) - \frac{\mathcal{E}(n_0)}{n_0}. \tag{chemical potential} \end{gather*}\]

For the power-law equation of state, we have

\[\begin{gather*} Q'(n_0) - \frac{Q(n_0)}{n_0} = \frac{\alpha(\beta-1)}{\beta^2}n_0^{\beta-1} = \overbrace{\frac{\alpha(\beta - 1)n_1^{\beta-1}}{\beta^2\pi^{\beta - 1}}}^{C} \frac{1}{(\sigma_x\sigma_y)^{\beta - 1}}, \\ \mathcal{E}'(n_0) - \frac{\mathcal{E}(n_0)}{n_0} = \frac{\alpha(\beta-1)}{\beta}n_0^{\beta-1} = \underbrace{\frac{\alpha\beta (\beta-1)n_1^{\beta-1}}{\beta^2\pi^{\beta - 1}}}_{\beta C} \frac{1}{(\sigma_x\sigma_y)^{\beta - 1}}. \end{gather*}\]

These are a little tricky to solve and we defer the discussion of their solution to Solving for Sigmas: Globally Convergent Newton’s Method. Once we solve for \(\sigma_x\) and \(\sigma_y\), we can express the energy and chemical potential as follows, eliminating the interaction term:

\[\begin{gather*} E_1[\phi] = (1-a)\frac{\hbar^2}{4m}\left( \frac{1}{\lambda_{x}^2\sigma_x^2} + \frac{1}{\lambda_{y}^2\sigma_y^2}\right) + (1+a)\frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4},\\ \mu_1[\phi] = (1-\beta a)\frac{\hbar^2}{4m}\left( \frac{1}{\lambda_{x}^2\sigma_x^2} + \frac{1}{\lambda_{y}^2\sigma_y^2}\right) + (1+\beta a)\frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4}. \end{gather*}\]

The value of \(a\) depends on whether we minimize the energy, or chemical potential:

\[\begin{gather*} \frac{\alpha}{\beta^2}\frac{n_1^{\beta-1}}{(\pi\sigma_x\sigma_y)^{\beta - 1}} = \frac{1}{2}\sum_{x,y} (-a) \left( \frac{\hbar^2}{2m}\frac{1}{\lambda_{x,y}^2\sigma_{x,y}^2} - \frac{m\omega_{x,y}^2\sigma_{x,y}^2}{2} \right),\\ a = \begin{cases} \frac{1}{\beta-1} & \text{(energy)},\\ \frac{1}{\beta(\beta - 1)} & \text{(chemical potential)}. \end{cases} \end{gather*}\]

GPE with \(\beta = 2\), Minimize Energy, \(a=1\):

\[\begin{gather*} E_1[\phi] = \frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{2},\\ \mu_1[\phi] = -\frac{\hbar^2}{4m}\left( \frac{1}{\lambda_{x}^2\sigma_x^2} + \frac{1}{\lambda_{y}^2\sigma_y^2} \right) + 3\frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4}. \end{gather*}\]

GPE with \(\beta = 2\), Minimize Chemical Potential, \(a=1/2\):

\[\begin{gather*} E_1[\phi] = \frac{1}{2} \frac{\hbar^2}{4m}\left( \frac{1}{\lambda_{x}^2\sigma_x^2} + \frac{1}{\lambda_{y}^2\sigma_y^2} \right) + \frac{3}{2} \frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{4},\\ \mu_1[\phi] = \frac{m(\omega_x^2\sigma_x^2 + \omega_y^2\sigma_y^2)}{2}. \end{gather*}\]

Old Results

Here are the same formulae in the symmetric limit \(\lambda_{x}=\lambda_{y}=1\), \(\omega_x = \omega_y = \omega_\perp\), and \(\sigma_x = \sigma_y = \sigma\):

\[\begin{gather*} E_1[\phi] = \frac{\hbar^2}{2m\sigma^2} + \frac{m\omega_\perp^2\sigma^2}{2} + \frac{\alpha}{\beta^2} \frac{n_1^{\beta-1}}{(\pi\sigma^2)^{\beta-1}},\\ \mu_1[\phi] = \frac{\hbar^2}{2m\sigma^2} + \frac{m\omega_\perp^2\sigma^2}{2} + \frac{\alpha}{\beta} \frac{n_1^{\beta-1}}{(\pi\sigma^2)^{\beta-1}},\\ \frac{\hbar^2}{2m\sigma^4} = \frac{m\omega_{\perp}^2}{2} + \frac{\alpha(1-\beta)}{\beta^2}\frac{n_1^{\beta-1}}{(\pi\sigma^2)^{\beta - 1}\sigma^2}, \tag{Energy} \\ \frac{\hbar^2}{2m\sigma^4} = \frac{m\omega_{\perp}^2}{2} + \frac{\alpha(1-\beta)}{\beta}\frac{n_1^{\beta-1}}{(\pi\sigma^2)^{\beta - 1}\sigma^2}, \tag{Chemical Potential} \end{gather*}\]

Explicit form for GPE

For the GPE with \(\alpha = g\), and \(\beta = 2\) we have:

\[\begin{gather*} \mu_1[n_1] = \frac{\frac{\hbar^2}{2m} + \frac{gn_1}{2\pi}}{\sigma^2} + \frac{m\omega_\perp^2\sigma^2}{2}, \qquad g = \frac{4\pi a \hbar^2}{m}, \\ \sigma^2 = \sqrt{\frac{\hbar^2}{m^2\omega_\perp^2} + \frac{gn_1}{2\pi m\omega_\perp^2}},\\ \mu_1[n_1] = \hbar\omega_\perp\left( \frac{1 + 3an_1}{\sqrt{1 + 2a n_1}} \right), \tag{Energy}\\ \sigma^2 = \sqrt{\frac{\hbar^2}{m^2\omega_\perp^2} + \frac{gn_1}{\pi m\omega_\perp^2}},\\ \mu_1[n_1] = m\omega_\perp^2\sigma^2 = \hbar \omega_\perp\sqrt{1 + 4a n_1}, \tag{Chemical Potential} \end{gather*}\]

HO Potential and TF Wavefunction (Incomplete)

An alternative strategy is to use the TF approximation radially:

\[\begin{gather*} \phi\bigl(x, y; \sigma(z, t)\bigr) = \sqrt{\frac{2\mu_1 - m\omega_\perp^2 (x^2+y^2)}{2gn_1}}, \qquad \mu_1[n_1] = \hbar\omega_\perp\sqrt{4 a n_1}. \end{gather*}\]

This is missing the zero-point energy of the radial motion. We do not use this.

HO Potential and Exact Wavefunction (Incomplete)

The exact solution requires solving radial equation. Here we present a general solution for this (numerically) for the GPE after first performing some scaling transformations. We start with a symmetric trap:

\[\begin{gather*} \left(\frac{-\hbar^2}{2m}\frac{1}{r}\diff{}{r}\left(r\diff{}{r}\right) + m\frac{\omega_\perp^2r^2}{2} + gn_1\abs{\phi}^2\right)\phi = \mu_1\phi,\\ u = \sqrt{r}\phi, \qquad \tilde{r} = \sqrt{\frac{m\omega_\perp}{\hbar}} r, \qquad \int 2\pi r\d{r}\; \abs{\psi}^2 = \int 2\pi \d{r}\; \abs{u}^2 = 1,\\ \left( - \nabla^2_{\tilde{r}} - \frac{1}{4\tilde{r}^2} + \tilde{r}^2 + \overbrace{2gn_1\sqrt{\frac{m}{\omega_\perp\hbar^3}}}^{\alpha}\frac{\abs{u}^2}{\tilde{r}} \right)u(\tilde{r}) = \frac{2\mu_1}{\hbar\omega_\perp} u(\tilde{r}),\\ \end{gather*}\]
\[\begin{gather*} \left(\frac{-\hbar^2(\nabla_x^2 + \nabla_y^2)}{2m} + m\frac{\omega_x^2x^2 + \omega_y^2y^2}{2} + gn_1\abs{\phi}^2\right)\phi = \mu_1\phi,\qquad \int\d{x}\d{y}\; \abs{\psi}^2 = 1,\\ \tilde{x} = \sqrt{\frac{m\omega_x}{\hbar}} x, \qquad \tilde{y} = \sqrt{\frac{m\omega_y}{\hbar}} y, \qquad \end{gather*}\]
\[\begin{gather*} \left(\frac{-\hbar(\omega_x\nabla_\tilde{x}^2 + \omega_y\nabla_\tilde{y}^2)}{2} + \hbar\frac{\omega_x\tilde{x}^2 + \omega_y\tilde{y}^2}{2} + gn_1\abs{\phi}^2\right)\phi = \mu_1\phi,\qquad \int\d{x}\d{y}\; \abs{\psi}^2 = 1,\\ \end{gather*}\]

Central Density#

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

\[\begin{gather*} n(0,z) = \frac{n_{1D}(z)}{\Lambda \pi \sigma_{x}\sigma_{y}}. \end{gather*}\]

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:

Let

\[\begin{gather*} \tilde{R}^2 = \frac{2\mu}{m} - \omega_z^2z^2 = \frac{2\alpha}{m} n_0^{\beta-1}. \end{gather*}\]

Then

\[\begin{gather*} \alpha n^{\beta-1}(\vect{x}) = \mu - \frac{m}{2}\sum_i\omega_i^2x_i^2, \qquad \alpha n^{\beta-1}_0 = \alpha n^{\beta-1}(0, z) = \mu - \frac{m}{2}\omega_z^2z^2,\\ n_{1D}(z) = \int\d{x}\d{y}\; n(\vect{x}) = \int\d{x}\d{y}\; \left( \frac{\mu - \frac{m}{2}\sum_i\omega_i^2x_i^2}{\alpha} \right)^{1/(\beta-1)}\\ = \frac{2\pi}{\omega_x\omega_y}\left(\frac{m}{2\alpha}\right)^{1/(\beta-1)} \int_0^{\tilde{R}}\tilde{r}\d{\tilde{r}}\; (\tilde{R}^2 - \tilde{r}^2)^{1/(\beta-1)}\\ = \frac{2\pi}{\omega_x\omega_y}\left(\frac{m}{2\alpha}\right)^{1/(\beta-1)} \frac{\tilde{R}^{2/(\beta-1)+2}}{2/(\beta-1)+2},\\ = \frac{2\pi}{\omega_x\omega_y}\left(\frac{m}{2\alpha}\right)^{1/(\beta-1)} \frac{\left(\frac{2\alpha}{m} n_0^{\beta-1}\right)^{1/(\beta-1)+1}}{2/(\beta-1)+2} = \frac{2\pi\alpha(\beta-1)}{m\beta\omega_x\omega_y}n_0^{\beta}. \end{gather*}\]

Hence:

\[\begin{gather*} n(0,z) = \frac{1}{\Lambda} \left(\frac{m\beta\omega_x\omega_y}{2\pi\alpha(\beta-1)}n_{1D}(z)\right)^{1/\beta}. \end{gather*}\]

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

GPE with \(\beta=2\) and \(\alpha = g\):

\[\begin{gather*} n(\vect{x}) = \frac{\mu - \frac{m}{2}\sum_i\omega_i^2x_i^2}{g}, \qquad n_0 = n(0, z) = \frac{\mu - \frac{m}{2}\omega_z^2z^2}{g},\\ n_{1D}(z) = \int\d{x}\d{y} n(\vect{x}) = \frac{m}{2g}\frac{2\pi}{\omega_x \omega_y} \int_0^{\tilde{R}}\tilde{r}\d{\tilde{r}}\; (\tilde{R}^2 - \tilde{r}^2) = \frac{\pi gn_0^2}{m \omega_x\omega_y}, \qquad \tilde{R}^2 = \frac{2\mu}{m} - \omega_z^2z^2 = \frac{2gn_0}{m},\\ n(0,z) = \frac{1}{\Lambda}\sqrt{\frac{m \omega_x\omega_y n_{1D}(z)}{\pi g}}. \end{gather*}\]

Thomas-Fermi Profiles#

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

\[\begin{gather*} n(z) = \abs{\psi}^2 = \frac{\mu - V(z)}{g} = \frac{V(R_{TF}) - V(z)}{g}. \end{gather*}\]

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 our effective 1D model, this must be modified somewhat. We start with the definition of the Thomas Fermi where \(n_{1D}(R) = 0\). As we shall see shortly, at this point \(\sigma_{x,y}^2(R) = \hbar/m\omega_{x,y}\) and the chemical potential must be augmented with the transverse energy (zero-point energy \(\hbar\omega/2\) in each direction):

\[\begin{gather*} \mu = V(R_{TF}) + \frac{\hbar\omega_x + \hbar\omega_y}{2}. \end{gather*}\]

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

Recall that the Thomas-Fermi approximation has the form

\[\begin{align*} \mu - V(z) &= \frac{A}{2}\left(\frac{1}{\sigma_x^2} + \frac{1}{\sigma_y^2}\right) + \frac{B_x\sigma_x^2 + B_y\sigma_y^2}{2} + \frac{a\beta C}{(\sigma_x\sigma_y)^{\beta-1}},\\ B_{x,y}\sigma_{x,y}^2 - \frac{A}{\sigma_{x,y}^2} & = \frac{C}{(\sigma_x\sigma_y)^{\beta - 1}}, \end{align*}\]

where

\[\begin{gather*} A = \frac{\hbar^2}{2m},\qquad C = \frac{\alpha}{a\beta^2} \left(\frac{n_1}{\pi}\right)^{\beta - 1},\qquad B_{x,y} = \frac{m\omega_{x,y}^2}{2},\\ a = \begin{cases} \frac{1}{\beta-1} & \text{(energy)},\\ \frac{1}{\beta(\beta - 1)} & \text{(chemical potential)}. \end{cases} \end{gather*}\]

While this can be solved in principle, the expressions are cumbersome, so for now, we consider only the case \(\omega_x = \omega_y \equiv \omega_\perp\) so that \(B_x = B_y = B\) and \(\sigma_x = \sigma_y = \sigma\). Then we can solve:

\[\begin{gather*} \mu - V(z) = \frac{A}{\sigma^2} + B\sigma^2 + \frac{a\beta C}{(\sigma^2)^{\beta-1}}, \qquad B\sigma^2 - \frac{A}{\sigma^2} = \frac{C}{(\sigma^2)^{\beta - 1}},\\ \frac{a\beta C}{(\sigma^2)^{\beta-1}} = \mu - V(z) - \frac{A}{\sigma^2} - B\sigma^2 = a\beta \left(B\sigma^2 - \frac{A}{\sigma^2}\right),\\ \mu - V(z) = (1 + a\beta) B\sigma^2 + (1 - a\beta)\frac{A}{\sigma^2}. \end{gather*}\]

This gives a quadratic equation for \(\sigma^2\) which in turn can be used to find \(C\) and \(n_1\):

\[\begin{gather*} B\sigma^2 = \frac{\Bigl(\mu - V(z)\Bigr) + \sqrt{\Bigl(\mu - V(z)\Bigr)^2 - 4AB(1 - a^2\beta^2)} }{2(1 + a\beta)},\\ C = \frac{\alpha}{a\beta^2} \left(\frac{n_1}{\pi}\right)^{\beta - 1} = (\sigma^2)^{\beta - 1}\left(B\sigma^2 - \frac{A}{\sigma^2}\right). \end{gather*}\]

Expanding the constants, we have

\[\begin{gather*} \frac{m\omega_\perp\sigma^2}{\hbar} = \frac{\frac{\mu - V(z)}{\hbar\omega_\perp} + \sqrt{\left(\frac{\mu - V(z)}{\hbar\omega_\perp}\right)^2 - (1 - a^2\beta^2)} }{1 + a\beta},\\ n_{TF}(z) = \pi\sigma^2 \left( \frac{a\beta^2 \hbar \omega_\perp}{2\alpha} \left[\frac{m\omega_\perp\sigma^2}{\hbar} - \frac{\hbar}{m\omega_\perp\sigma^2}\right]\right)^{1/(\beta-1)}. \end{gather*}\]

GPE (\(\alpha=g\), \(\beta = 2\)), Minimize Energy, \(a=1\):

\[\begin{gather*} \frac{m\omega_\perp\sigma^2}{\hbar} = \frac{1}{3}\left(\frac{[\mu - V(z)]}{\hbar\omega_\perp} + \sqrt{\frac{[\mu - V(z)]^2}{\hbar^2\omega_\perp^2} + 3} \right) ,\\ \frac{4\pi C}{g} = n_1 = \frac{2\pi \hbar^2}{mg}\left(\left(\frac{m\omega_\perp\sigma^2}{\hbar}\right)^2 - 1\right). \end{gather*}\]

GPE (\(\alpha=g\), \(\beta = 2\)), Minimize Chemical Potential, \(a=1/2\):

\[\begin{gather*} \frac{m\omega_\perp \sigma^2}{\hbar} = \frac{\mu - V(z)}{\hbar\omega_\perp},\\ \frac{2\pi C}{g} = n_1 = \frac{\pi\hbar^2}{mg}\left(\left(\frac{m\omega_\perp\sigma^2}{\hbar}\right)^2 - 1\right) = \frac{[\mu - V(z)]^2 - \hbar^2\omega_\perp^2}{gm\omega_\perp^2/\pi}. \end{gather*}\]

2D Version (Pancakes)#

We can repeat the same analysis for clouds that are tightly confined in a single dimension (pancakes) to obtain an effective 2D NPSEQ:

\[\begin{gather*} \Psi(\vect{x},t) = \phi(\vect{x},t)\psi(x, y, t), \qquad \int \abs{\psi(\vect{x},t)}^2\d{z} = 1\\ \phi(\vect{x}, t) = \frac{1}{\sqrt{\pi \sigma}}e^{-z^2/2\sigma},\\ n(\vect{x},t) = \abs{\phi}^2(\vect{x},t)\overbrace{\abs{\psi}^2(x,y,t)}^{n_2(z,y,t)},\\ \mathcal{E}_2[\phi] = \int \d{z}\left( \frac{\hbar^2 \abs{\nabla_z \phi}^2}{2m} + V_z(z)\abs{\phi}^2 + \frac{\mathcal{E}(n_2\abs{\phi}^2)}{n_2} \right),\\ \mathcal{\mu}_2[\phi] = \int \d{z}\; \phi^*\left( \frac{-\hbar^2 \nabla_z^2}{2m} + V_z(z) + \mathcal{E}'(n_2\abs{\phi}^2) \right)\phi,\\ n_1\mathcal{E}_1 = \cdots + 2\sigma Q_2(n_0),\\ n_1\mu_1 = \cdots + 2\sigma R_2(n_0),\\ Q_2(n_0) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{n\sqrt{\ln n_0 - \ln n}}, \qquad R_2(n_0) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{\sqrt{\ln n_0 - \ln n}},\\ n_0 = \frac{n_2}{\sqrt{\pi}\sigma}, \qquad Q_2'(n) = \frac{\mathcal{E}(n)}{n\sqrt{\ln n_0 - \ln n}}, \qquad R_2'(n) = \frac{\mathcal{E}(n)}{\sqrt{\ln n_0 - \ln n}}. \end{gather*}\]

For power-law EoS:

\[\begin{gather*} \mathcal{E}(n) = \alpha \frac{n^{\beta}}{\beta}\\ Q_2(n_0) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{n\sqrt{\ln n_0 - \ln n}}, \qquad R_2(n_0) = \int_0^{n_0}\d{n}\;\frac{\mathcal{E}(n)}{\sqrt{\ln n_0 - \ln n}},\\ \end{gather*}\]

Demonstration#

{code-cell}
from importlib import reload
import gpe.tube;reload(gpe.tube)
from gpe.tube import StateGPEdrZ
from gpe.minimize import MinimizeState


class State(StateGPEdrZ):
    _mu = 0.0

    def get_ws_perp(self, t):
        return self.ws[1:]

    def _get_Vext_(self):
        mu = self._mu if self._mu else 0.0
        return super()._get_Vext_() - mu


s = State(
    Lxyz=(10.0,),
    Nxyz=(64,),
    ws=(1.1, 1.2, 1.3),
    g=10.0,
)
x = s.xyz[0]
x_TF = 2.0
V_TF = s.get_V_TF(x_TF=x_TF)
s._mu = 0.0
s._mu = s.get_mu_from_V_TF(V_TF=V_TF)
n_TF = s.get_n_TF(V_TF=V_TF)
s = MinimizeState(s, fix_N=False).minimize()
# plt.plot(x, n_TF)
plt.plot(x, s.get_density() - n_TF);

{code-cell}
m, g, w = s.m, s.g, s.w0_perp
w2 = w**2
assert s.hbar == 1
sigma2 = s.get_sigma2()
n = s.get_density()
# plt.plot(np.sqrt(1 + m*g*n/2/np.pi)/m/w_perp - sigma2)
V = s.get_Vext()

def get_sigma2_TF(mu):
    mu_eff = mu - s.get_Vext()
    return mu_eff / 3 / m / w2 + 1.0 / 3.0 / m / w * np.sqrt(mu_eff**2 / w2 + 3)

def get_n_TF(mu):
    s2 = get_sigma2_TF(mu)
    return 2 * np.pi * np.maximum((m * w * s2) ** 2 - 1, 0) / m / g

mu = ((3 * m * w2 * sigma2**2 / 2 - 1.0 / 2.0 / m) / sigma2 + V).min()
# plt.plot(get_sigma2_TF(mu))
# plt.plot(sigma2)
plt.plot((get_n_TF(mu) - n) / n.max());

{code-cell}
from gpe import utils, bec, tube, minimize

reload(bec)
reload(tube)


class State(bec.State):
    # Goal:
    # healing length = 1.0
    # trap length = 2.0
    # central density = 1.0
    # g = 1.0
    # cloud radius R = 5.0
    # Box size = 20.0
    # mw^2R^2/2 = mu = gn
    def __init__(self, dim=3, g=1.0, trap_a=0.0, trap_T=1.0):
        m = 1.0
        n0 = 1.0
        mu = g * n0
        Rxyz = np.array([10.0, 2.0, 2.0])
        ws = np.sqrt(2 * mu / m) / Rxyz
        dx = 0.25
        Lxyz = 3.0 * Rxyz
        Nxyz = list(map(utils.get_good_N, Lxyz / dx))
        self.trap_a = trap_a
        self.trap_T = trap_T
        args = dict(
            g=g,
            m=m,
            ws=ws,
            Lxyz=Lxyz[:dim],
            Nxyz=Nxyz[:dim],
            mu=mu,
        )
        super().__init__(**args)
        self.mu = mu

    def get_ws_perp(self, t=None):
        if t is None:
            t = self.t
        if t <= 0:
            return self.ws[1:]
        else:
            w_trap = 2 * np.pi / self.trap_T
            ws = 1 * self.ws[1:]
            ws *= 1 + self.trap_a * np.sin(w_trap * t)
            return ws

    def get_ws(self, t=None):
        ws = 1 * self.ws
        ws[1:] = self.get_ws_perp(t=t)
        return ws

    def _get_Vext_(self):
        return (
            0.5 * self.m * sum((_w * _x) ** 2 for _w, _x in zip(self.get_ws(), self.xyz))
            - self.mu
        )


class State1(tube.StateGPEdrZ):
    def __init__(self, state):
        args = dict(
            Nxyz=state.basis.Nxyz[:1],
            Lxyz=state.basis.Lxyz[:1],
            ws=state.ws,
            g=state.g,
            m=state.m,
            mu=state.mu,
        )
        self.trap_a = state.trap_a
        self.trap_T = state.trap_T
        super().__init__(**args)

    def get_ws_perp(self, t=None):
        if t is None:
            t = self.t
        if t <= 0:
            return self.ws[1:]
        else:
            w_trap = 2 * np.pi / self.trap_T
            ws = 1 * self.ws[1:]
            ws *= 1 + self.trap_a * np.sin(w_trap * t)
            return ws

    def get_ws(self, t=None):
        ws = 1 * self.ws
        ws[1:] = self.get_ws_perp(t=t)
        return ws

    def _get_Vext_(self):
        return (
            0.5 * self.m * sum((_w * _x) ** 2 for _w, _x in zip(self.get_ws(), self.xyz))
            - self.mu
        )


def get_states(g=1.0):
    """Returns three states `(s, s_1D, s_tube)` for comparison:

    Returns
    -------
    s : Full 3D state.
    s_1D : 1D state representing a translationally invariant system in y and z.
       This should match the central density of s in the TF limit.
    s_tube : 1D state but in the dr-GPE-Z approximation where the cross-section
       of the cloud is modelled as a Gaussian.
    """
    s = minimize.MinimizeState(State(g=g), fix_N=False).minimize()
    Nx, Ny, Nz = s.basis.Nxyz
    n_central = s.get_density()[:, Ny // 2, Nz // 2]

    # w_perp = np.prod(s.ws[1:])**(1./len(s.ws[1:]))
    # a_perp = np.sqrt(s.hbar/s.m/w_perp)
    # zeta_1_2 = -1.46035450880958681288949915251529801246722933101258
    # scattering_length = s.g/s.m/4/np.pi/s.hbar**2
    # g = (a_perp * s.hbar * w_perp * 2 * scattering_length/a_perp
    #              / (1 - scattering_length/a_perp * abs(zeta_1_2)))
    # print(g)
    s_1D = State(dim=1, g=g)
    s_1D[...] = np.sqrt(n_central)
    s_1D.mu = s.mu
    s_1D = minimize.MinimizeState(s_1D, fix_N=False).minimize()

    s_tube = State1(state=s)
    s_tube[...] *= np.sqrt(s.get_N() / s_tube.get_N())
    s_tube = minimize.MinimizeState(s_tube, fix_N=True).minimize()
    return s, s_1D, s_tube

{code-cell}
s_g1 = get_states(g=1.0)
s_g2 = get_states(g=2.0)
s_g10 = get_states(g=10.0)
s_g1[0].plot()
plt.figure()
s_g10[0].plot()

{code-cell}
plt.figure(figsize=(10, 7))


def plot_central_density(s, s_1D, s_tube):
    Nx, Ny, Nz = s.basis.Nxyz
    n_central = s.get_density()[:, Ny // 2, Nz // 2]
    plt.plot(s_1D.xyz[0], s_1D.get_density(), label="1D")
    plt.plot(s.xyz[0].ravel(), s.get_density().max(axis=-1).max(axis=-1), label="3D")
    plt.plot(s_tube.xyz[0], s_tube.get_central_density(), label="dr-GPE")
    plt.plot(s_tube.xyz[0], s_tube.get_central_density(TF=True), label="dr-GPE-TF")
    plt.ylabel("Central Density")
    plt.legend(loc="best")


def plot_1D_density(s, s_1D, s_tube):
    Nx, Ny, Nz = s.basis.Nxyz
    plt.plot(s_1D.xyz[0], s_1D.get_density(), label="1D")
    dxyz = s.basis.Lxyz / s.basis.Nxyz
    plt.plot(
        s.xyz[0].ravel(),
        s.get_density().sum(axis=-1).sum(axis=-1) * np.prod(dxyz[1:]),
        label="3D",
    )
    plt.plot(s_tube.xyz[0], s_tube.get_density(), label="dr-GPE")
    plt.ylabel("1D integrated density")
    plt.legend(loc="best")


for n, (s, s_1D, s_tube) in enumerate([s_g1, s_g2, s_g10]):
    plt.subplot(231 + n)
    plt.title("g={}".format(s.g))
    plot_central_density(s=s, s_1D=s_1D, s_tube=s_tube)

    plt.subplot(234 + n)
    plot_1D_density(s=s, s_1D=s_1D, s_tube=s_tube)

Here we see that in the TF limit (large \(g\), small healing length, right) that the 1D model and the 3D model have very similar behavior, but the tube deviates. This is because the cross-section in the trap is a Gaussian which is only valid for tightly confining traps. For tightly confining traps, the 1D model does not work well. We also compare the extraction of the central density assuming a TF profile dr-GPR-TF which does work well in the right plot. Which version you use depends on the trap parameters.

In the bottom, we compare the integrated line density integrating across the trap in the \(x\) and \(y\) directions. This is well-modeled by the gr-GPE-Z model as expected because this is how the Gaussian is fit. The 1D model does not represent this well at all (also as expected).

Here we add some dynamics in the form of a sinusoidally varying \(\omega(t)\).

{code-cell}
from IPython.display import display, clear_output
from pytimeode.evolvers import EvolverABM
from mmfutils.contexts import NoInterrupt

ss = s_g2
s = ss[0]
dt = 0.2 * s.t_scale

es = []
for _s in ss:
    s = _s.copy()
    x = s.xyz[0]
    s.trap_a = 0.1
    # s.trap_T = 10.0
    # s[...] *= (np.tanh(0.1*(x - 2))**2 + 10.0)/11.0
    s[...] *= np.tanh(0.1 * (x - 2))
    es.append(EvolverABM(s, dt=dt))

with NoInterrupt() as interrupted:
    while not interrupted:
        for e in es:
            e.evolve(200)
        plt.clf()
        plt.subplot(211)
        plot_1D_density(*[_e.y for _e in es])
        plt.subplot(212)
        plot_central_density(*[_e.y for _e in es])

        plt.suptitle("t={:0.2f}, g={}".format(e.t, e.y.g))
        display(plt.gcf())
        clear_output(wait=True)

This works well until about 29ms. I am not sure what the source of the errors are at this point - possibly due to the Z approximation, or maybe due to neglecting the derivatives of \(\sigma_{,z}\)?

TF Approximation#

Now we consider the TF approximation and chemical potentials. The previous states were chosen so that the total particle number is the same. Here we compare how this relates to the chemical potential \(\mu\) and \(V_{TF} = V(x_{TF})\) at the boundary of the cloud:

{code-cell}
plt.figure(figsize=(10, 3))
for n, (s, s_1D, s_tube) in enumerate([s_g1, s_g2, s_g10]):
    plt.subplot(131 + n)
    plt.title("g={}".format(s.g))
    for _s in (s, s_tube):
        _s0 = _s.copy()
        _s0.set_initial_state(x_TF=9)
        (l,) = plt.plot(_s.xyz[0].ravel(), _s.get_density_x())
        plt.plot(_s0.xyz[0].ravel(), _s0.get_density_x(), c=l.get_c(), ls=":")
        # plot_central_density(s=s, s_1D=s_1D, s_tube=s_tube)

        # ## References

        # * [Massignan:2003]: One-dimensional model for the dynamics and expansion of elongated Bose-Einstein condensates.  This is where they introduce the so-called dynamically rescaled GPE (dr-GPE).
        # * [Mateo:2009]: Clearest presentation of the results from [Mateo:2006], [Mateo:2007], and [Mateo:2008].  Here they consider two forms of the NPSEQ.
        #
        # [Massignan:2003]: http://dx.doi.org/10.1103/PhysRevA.67.023614 'Pietro Massignan and Michele Modugno, "One-dimensional model for the dynamics and expansion of elongated Bose-Einstein condensates", Phys. Rev. A 67, 023614 (2003)'
        # [Mateo:2006]: http://dx.doi.org/10.1103/physreva.74.065602 'A. Muñoz Mateo and V. Delgado, "Extension of the Thomas-Fermi approximation for trapped Bose-Einstein condensates with an arbitrary number of atoms", pra 74(6),  (2006)'
        # [Mateo:2007]: https://doi.org/10.1103/PhysRevA.75.063610 'A. Mu\~noz Mateo and V. Delgado, "Ground-state properties of trapped Bose-Einstein condensates: Extension of the Thomas-Fermi approximation", pra 75, 063610 (2007)'
        # [Mateo:2008]: https://doi.org/10.1103/PhysRevA.77.013617 'A. Mu\~noz Mateo and V. Delgado, "Effective mean-field equations for cigar-shaped and disk-shaped Bose-Einstein condensates", pra 77, 013617 (2008) '
        # [Mateo:2009]: https://doi.org/10.1016/j.aop.2008.10.002 'A. Muñoz Mateo and V. Delgado, "Effective one-dimensional dynamics of elongated Bose--Einstein condensates", Annals of Physics 324(3), 709 - 724 (2009) '

        # ## Separable Ansatz

        # The other portion of the analysis is to factorize the wavefunction as
        #
        # \begin{gather*}
        #   \Psi(x, y, z; t) = \phi(x, y, z, t)\psi(z, t), \qquad
        #   \int \abs{\phi}^2(x, y, z, t)\;\d{x}\d{y} = 1.
        # \end{gather*}
        #
        # The normalization condition gives the interpretation of $\abs{\psi}^2(z,t)$ as the density of particles per unit length along $z$.
        #
        # Noting that the GPE derives from a principle of minimal action:
        #
        # \begin{gather*}
        #   S = \int \Psi^\dagger \left(
        #   \I\hbar\partial_t
        #   - \frac{\hbar^2\overleftarrow\nabla\overrightarrow\nabla}{2m}
        #   - V_\perp(x, y) - V_z(z) - \frac{g}{2}\abs{\Psi}^2\right)\d{x}\d{y}\d{z},
        # \end{gather*}
        #
        # we obtain the following equation for $\psi(z, t)$:
        #
        # \begin{gather*}
        #   \newcommand{\Braket}[1]{\left\langle{#1}\right\rangle}
        #   \I\hbar\dot{\psi} = \left(
        #   - \frac{\hbar^2\nabla_z^2}{2m}
        # + V_ {"incorrectly_encoded_metadata": "{1D}(z) + g_{1D} \\abs{\\psi}^2\\right)\\psi\\\\"}
        #   \braket{A} = \frac{\int \phi^\dagger A\phi\; \d{x}\d{y}}
        #                     {\int \abs{\phi}^2\;\d{x}\d{y}}, \qquad
        #   V_{1D}(z) = V_z(z) + \Braket{
        #     -\I\partial_t
        #     - \frac{\hbar^2\nabla^2}{2m}
        # + V_\perp(x, y)}, \qquad
        #   g_{1D} = \braket{\abs{\phi}^2}.
        # \end{gather*}
        #
        # The normalization condition for $\phi$ ensures that the linear term $C\nabla_z$ with $C\propto\braket{\nabla_Z} = 0$ vanishes.

        # Likewise, minimizing wrt $\phi$ gives:
        #
        # \begin{gather*}
        #   \I\hbar \dot{\phi} = \left(
        #     -\frac{\hbar^2\nabla^2}{2m} + V_\perp(x, y) + V_z(z) - \mu + g\abs{\Psi}^2
        # + \frac {"incorrectly_encoded_metadata": "{\\hbar^2\\abs{\\nabla_z\\psi}^2}{2m\\abs{\\psi}^2}"}
        #     \right)\phi,\\
        #   V_{2D}(x, y, z) = V_\perp(x, y) + V_z(z).
        # \end{gather*}

        # So far this is exact.  To proceed, one assumes that the last term above is small, and that $\phi$ instantly adjusts so that the time-dependence can be dropped.

        # Proceeding further, we factor the wavefunction as:
        #
        # \begin{gather*}
        #   \Phi(x, y, Z;t) = \phi(x, y, t; Z)\psi(Z, t)
        # \end{gather*}
        #
        # and make the "slowly varying approximation" that the dependence of $\phi$ on $Z$ is slowly varying so that we may ignore derivatives.  This allows $\phi$ to commute with the term $E(\op{p}_Z)$ so that the equations separate:
        #
        # \begin{gather*}
        #   \I\hbar\frac{\dot{\phi}}{\phi} =
        #     \frac{\frac{-\hbar^2(\nabla_X^2+\nabla_Y^2)}{2m}\phi}{\phi} + \epsilon(\abs\Psi^2)
        # + \frac {"incorrectly_encoded_metadata": "{E(p_Z)\\psi}{\\psi} + V(Z, t) - \\I\\hbar \\frac{\\dot{\\psi}}{\\psi}."}
        # \end{gather*}
        #
        # If we make the ansatz
        #
        # \begin{gather*}
        #   \phi\bigl(x, y; \sigma(Z, t)\bigr) = \frac{1}{\sqrt{\pi}\sigma(Z, t)}e^{-(x^2+y^2)/2\sigma^2}, \\
        #   \frac{\frac{-\hbar^2(\nabla_X^2+\nabla_Y^2)}{2m}\phi}{\phi} = \frac{\hbar^2}{2m}\frac{2\sigma^2-x^2-y^2}{\sigma^4}
        # \end{gather*}

        # ## Effective $g$ in Axial Confinement

        # [Olshanii:1998] gives the relationship:
        #
        # \begin{gather*}
        #   a_{1D} = \frac{-a_\perp^2}{2a}\left(1 - C\frac{a}{a_{\perp}}\right), \qquad
        #   a_\perp = \sqrt{\frac{\hbar}{m\omega_\perp}}.\\
        #   g_{1D} = \frac{\hbar \omega_\perp 2 a}{1-C\frac{a}{a_\perp}}
        #          = \frac{2\hbar^2}{m}\frac{a}{a_\perp^2}\frac{1}{1-C\frac{a}{a_\perp}}
        #          = \frac{-a_\perp^2\hbar\omega_\perp}{a_{1D}}
        #          = \frac{-\hbar^2}{ma_{1D}},\\
        #          C = -\zeta(1/2) = 1.46035450880958681288949915251529801246722933101258.
        # \end{gather*}
        #
        #
        # but see [Zhang:2014] for a correction with SOC.
        #
        # [Olshanii:1998]: http://dx.doi.org/10.1103/PhysRevLett.81.938 'M. Olshanii, "Atomic Scattering in the Presence of an External Confinement and a Gas of Impenetrable Bosons", Phys. Rev. Lett. 81, 938--941 (1998)'
        #
        # [Zhang:2014]: http://dx.doi.org/10.1038/srep04992 'Yi-Cai Zhang, Shu-Wei Song, and Wu-Ming Liu, "The confinement induced resonance in spin-orbit coupled cold atoms with Raman coupling", Sci. Rep. 4, 4992 (2014)'
        #

        # In the weak-confinement limit $a \ll a_\perp$ this becomes:
        #
        # \begin{gather*}
        #   g_{1D} = 2 a \hbar \omega_\perp = \frac{2 \hbar^2}{m}\frac{a}{a_\perp^2}.
        # \end{gather*}

        zeta_1_2 = -1.46035450880958681288949915251529801246722933101258
        gs_1d = (
            a_perp
            * u.hbar
            * w_perp
            * 2
            * a_s
            / a_perp
            / (1 + a_s / a_perp * abs(zeta_1_2))
        )