---
execution:
  timeout: 30
jupytext:
  cell_metadata_json: true
  encoding: '# -*- coding: utf-8 -*-'
  formats: md:myst,ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

(sec:NPSEQ)=
# Non-Polynomial Schrödinger Equation (NPSEQ)

```{code-cell}
:init_cell: true
:tags: [hide-cell]

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

## 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)$.

:::{margin}
A major use-case is to deal with expansion, which we do with a dynamic rescaling of the
coordinates.  This induces an effective radial potential (the equivalent of the
centrifugal force introduced by a rotating frame) and the remaining theory can be
reduced to an effective 1D theory as we discuss here. 
:::
## 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$.
   
:::{admonition} 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
{cite}`Mateo: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.
:::

:::{margin}
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
{ref}`sec:drGPE` 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}$.

In this expressions $n_0$ is the central density.
:::
:::::{admonition} 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*}
:::{admonition} 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$.

:::{admonition} 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 {class}`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 {cite}`Mateo: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) {cite}`Mateo: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 {ref}`sec:sigmas` 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 {ref}`sec:sigmas`.

## 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).

:::{margin}
Since $\phi$ is normalized, the expectation $\braket{}_{\phi}$ here is really just the
integral over the transverse coordinates $\int\d{x}\d{y}$, but using the solution $\phi$
that minimizes the action.
:::
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*}
:::{admonition} 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*}
:::{admonition} 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 {ref}`sec:drGPE` 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.

:::{margin}
Formally we should include here all of the potential terms:
\begin{gather*}
  f(n) = V(\vect{x}) n + \mathcal{E}(n).
\end{gather*}
We assume, however, that the potential is separable:
\begin{gather*}
  V(\vect{x}) = V_z(z) +V_{\perp}(x, y),
\end{gather*}
and then we drop the $V_z(z)$ term as it will just give a constant shift.

Note also that we only include the transverse derivatives here in the kinetic term
$\abs{\vect{\nabla}_{\perp}\phi}^2$.  Formally, we should include all derivatives, but,
as per the adiabatic approximation, we neglect variations $\nabla_z\phi$. 
:::
### 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}$. {cite}`Mateo: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.

:::{admonition} 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.
:::{admonition} 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$.
:::

:::{margin}
Is there a name/symbol for $Q(n)$ whose derivative is the energy per particle?
\begin{gather*}
  Q(n) = \int_0^{n}\d{n}\frac{\mathcal{E}(n)}{n}, \\
  Q'(n) = \frac{\mathcal{E}(n)}{n} = \frac{E}{N} = \epsilon(n).
\end{gather*}
:::
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.

:::::{admonition} 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*}
:::::

:::::{admonition} Integration Details
:class: dropdown

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 {rec}`sec:NPGPE2`, 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*}
:::::

```{code-cell}
:tags: [hide-cell]

# 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)
```

:::{margin}
We have included the factor of $\sqrt{\pi}$ in the definitions of $Q_2(n_0)$ and
$R_2(n_0)$ to keep this symmetric form.
:::
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.
:::::{admonition} Expectation Values
:class: dropdown

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*}
:::::

```{code-cell}
:tags: [hide-input]

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*}
:::{admonition} 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*}
:::
:::{margin}
In previous versions of the code where we assumed symmetry $\lambda_x = \lambda_y$, we
absorbed these scaling factors into a redefined mass $m$.
:::
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
{ref}`sec:drGPE` 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}$.

:::{margin}
The coefficients $A_{x,y}$, and $B_{x,y}$ will be used in {ref}`sec:sigmas`.
:::
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*}
:::{admonition} 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 {ref}`sec:sigmas`.

:::{admonition} 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




















:::::{admonition} 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. {cite}`Mateo: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*}

:::{admonition} 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 {ref}`sec:drGPE` 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*}

:::{margin}
The coefficients $A_{x,y}$, $B_{x,y}$, and $C$ will be used in {ref}`sec:sigmas`.  Note that,
in terms of the parameter $a$ defined below:
\begin{gather*}
  C = \frac{\alpha}{a\beta^2}\left(\frac{n_1}{\pi}\right)^{\beta-1}.
\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
{ref}`sec:sigmas`.  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*}

:::{margin}
We choose the sign of $a$ here so that the value is typically positive.
:::
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*}

:::{admonition} 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*}
:::
  
:::{admonition} 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*}
:::


:::::{admonition} Old Results
:class: dropdown

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*}

```{code-cell}
from mmfutils.math.bases import CylindricalBasis
b = CylindricalBasis(Nxr=(1, 32), Lxr=(1, 10.0))
x, r = b.xyz
V = r**2
K = b._get_K()[0]
H = -K + np.diag(V)
```

\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

:::{margin}
The factor of $1/\Lambda = 1/\lambda_{x}\lambda_{y}\lambda_{z}$ is for the dynamic rescaling as
discussed in {ref}`sec:drGPE`.
:::
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:
:::{margin}
Note that
\begin{gather*}
  \int_0^{R}\d{r}\; r(R^2-r^2)^a = \frac{R^{2a+2}}{2a+2}.
\end{gather*}
:::
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)`.

:::{admonition} GPE with $\beta=2$ and $\alpha = g$:
:class: dropdown

\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.)*

:::{margin}
Using the results below, at $R_{TF}$, the density vanishes, so $C=0$ and we have
\begin{gather*}
  \sigma_{x,y}^2 = \sqrt{\frac{A}{B_{x,y}}} = \frac{\hbar}{m\omega_{x,y}}.
\end{gather*}
Inserting this we find
\begin{gather*}
  \mu - V(z)= \frac{\hbar\omega_{x} + \hbar\omega_{y}}{2}.
\end{gather*}
:::
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*}

:::::{admonition} 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*}

:::::
:::::{admonition} 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))
        )
```