---
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-1D)=
# Non-Polynomial Schrödinger Equation (NPSEQ) in 1D (Tube)

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

Here we summarize the approach and formulation for a 1D non-polynomial Schrödinger
equation representing an atomic gas trapped in a tube.  We start by noting that the GPE
is almost separable:
\begin{gather*}
  \I\hbar \dot{\Psi}
  =
  \left(
    \frac{-\hbar^2}{2m}\nabla^2_{\perp} + V(x, y)
  \right)\Psi
  +
  \left(
    \frac{-\hbar^2}{2m}\nabla^2_{z} + V(z)
  \right)\Psi
  + 
  \mathcal{E}'(n)\Psi.
\end{gather*}
The idea is to factor the wavefunction as
\begin{gather*}
  \Psi(\vect{x}, t) = \phi\bigl(x, y, n_1(z, t)\bigr)\psi(z, t),\\
  n_1(z, t) = \abs{\psi(z, t)}^2, \qquad
  n(x, y, z, t) = n_1(z, t)\Abs{\phi\bigl(x, y, n_1(z, t)\bigr)}^2,\\
  \iint\d{x}\d{y}\Abs{\phi\bigl(x, y, n_1(z, t)\bigr)}^2 = 1.
\end{gather*}
We then make the "adiabatic" approximation of neglecting the $z$ and $t$ dependence in
$\phi$ -- the assumption is that the dynamics along $z$ are slow, so that the radial
wavefunction $\phi$ adjusts instantaneously.  Within this approximation, we can
separate the GPE:
\begin{gather*}
  \I\hbar \phi \dot{\psi}
  =
  \psi
  \left(
    \frac{-\hbar^2}{2m}\nabla^2_{\perp} + V(x, y)
  \right)\phi
  +
  \phi
  \left(
    \frac{-\hbar^2}{2m}\nabla^2_{z} + V(z)
  \right)\psi
  + 
  \mathcal{E}'(n)\phi\psi.
\end{gather*}
Our assumption is that the radial wavefunction $\phi(x, y)$ at a fixed $z$, $t$, and
density $n_1$ adjusts instantaneously to the ground state in the transverse direction
\begin{gather*}
  \left(
    \frac{-\hbar^2}{2m}\nabla^2_{\perp} 
    + V(x, y)
    + \mathcal{E}'\Bigl(n_1\abs{\phi}^2\Bigr)
  \right)\phi
  = \mu_{\perp}(n_1)\phi.
\end{gather*}
:::{note}
The key outcomes of this step are the chemical potential $\mu_{\perp}(n_1)$ and the solutions
$\phi_{n_1}(x, y)$ from which we will compute various moments and related properties.
Here it is important that $V(x, y, z, t) = V_{\perp}(x, y) + V_z(z, t)$ separate so that
these quantities depend only on $n_1$ and not on $z$ or $t$.  In principle, a more
complicated dependence could be done, but this would not allow us to tabulate everything
with 1-dimensional splines like we will do later.
:::

To obtain the NPSEQ, we plug this solution back into the separated GPE above, multiply
by $\phi^\dagger$ and integrate across the transverse directions.
\begin{gather*}
  \I\hbar \dot{\psi}
  =
  \left(
    \frac{-\hbar^2}{2m}\nabla^2_{z} + V(z) + \mu_\perp(\abs{\psi}^2)
  \right)\psi.
\end{gather*}

## Example: The Standard dr-GPE

As a check, consider the original formulation of the dr-GPE for axially-symmetric
harmonic traps:
\begin{gather*}
  V(x, y) = \frac{m\omega_\perp^2}{2}(x^2 + y^2), \qquad
  \mathcal{E}(n) = \frac{gn^2}{2}.
\end{gather*}
The final result with the gaussian approximation is
\begin{gather*}
  \mu_{\perp}(n_1) = 
  \frac{\hbar^2}{2m\sigma^2}
  + \frac{m\omega_\perp^2\sigma^2}{2} 
  + a\frac{gn_1}{4\pi\sigma^2},\qquad
  \frac{m\omega_\perp^2}{2}\sigma^4 = \frac{\hbar^2}{2m} + a\frac{gn_1}{4\pi},
\end{gather*}
where the constant $a=1$ for energy minimization or $a=2$ for the chemical-potential
minimization recommended by {cite}`Mateo:2009`.

```{code-cell}
from gpe.Examples.tutorial import StateHOConvergence2
s0 = StateHOConvergence2(g=100, Lx=30.0, Nx=128, dmu=10.0)
axs = s0.plot(label="TF")
s = s0.get_initialized_state()
axs = s.plot(axs=axs, label="Minimized", plot_convergence=True)
axs[0].legend(loc='right');
```

Our code solves the 2D GPE with
\begin{gather*}
  ( + g_{2D} \abs{\Phi}^2)\Phi = \mu \Phi\\
  \phi = \frac{\Phi}{\sqrt{N}}\\
  ( + g_{2D} N \abs{\phi}^2)\phi = \mu \phi\\
  g_{2D} N  = g n_1
\end{gather*}

Here we show the relative error between $\mu(n_1)$ computed numerically and  the formula
above from the dr-GPE: 
```{code-cell}
dmus = np.linspace(0.0001, 10.0, 100)
ss = [StateHOConvergence2(g=1, Lx=30.0, Nx=128, dmu=dmu).get_initialized_state()
      for dmu in dmus]
mus = np.array([s.get_mu() for s in ss])
Ns = np.array([s.get_N() for s in ss])
gtildes = np.array([s.get_gtilde() for s in ss])

g = 1.
n1s = s.g * Ns / g

s = ss[0]
a = 2
n1 = 1
hbar2_2m = s.hbar**2 / 2 / s.m
mw2_2 = s.m * s.w**2 / 2
gn1s = s.g * Ns
agn_4pi = a * gn1s / 4 / np.pi
sigma2 = np.sqrt((hbar2_2m + agn_4pi) / mw2_2)
mus_ = (hbar2_2m + agn_4pi)/sigma2 + mw2_2 * sigma2
fig, ax = plt.subplots()
ax.plot(n1s, mus/mus_ - 1)
ax.set(xlabel=r"$\mu_1$", ylabel=r"$\mu/\mu_{dr-GPE}-1$")
ax1 = plt.twinx()
ax1.plot(n1s, gtildes)
```

<!--


## Introduction

Here we specialize to the tube case, deriving 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 only the
tube case here $\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.  The full GPE takes the following 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
  + \overbrace{V(\vect{x}, t)n(\vect{x}, t) + \frac{g n^2(\vect{x}, t)}{2}}^{f(n)}, \\
  \I\hbar \dot{\Psi}(\vect{x}, t) = \Biggl[
    -\frac{\hbar^2}{2m}\nabla^2 + \underbrace{V(\vect{x}, t) + gn(\vect{x}, t)}_{f'(n)}
  \Biggr]\Psi(\vect{x}, t).
\end{gather*}

## Derivation

Consider the following ansatz:
\begin{gather*}
  \Psi(\vect{x}, t) = \phi\bigl(x, y;n_1(z, t)\bigr)\psi(z, t), \qquad
  \iint \abs{\phi}^2\d{x}\d{y} = 1.
\end{gather*}
Here we assume that the transverse wavefunction responds instantly as the density
$n_1(z, t) = \abs{\psi(z, t)}^2$ changes.  Varying the action with respect to $\phi$ and
$\psi$ we find:
:::{margin}
Here we suppress the arguments $\phi\bigl(x, y; n_1(z, t)\bigr) \rightarrow \phi$ and
$\psi(z, t) \rightarrow \psi$ to simplify the notation.
:::
\begin{gather*}
  \frac{\delta S}{\delta\psi^\dagger(z, t)} = 0
  = \int\d{x}\d{y}\Biggl\{\I\hbar \abs{\phi}^2\dot{\psi}
    + \I\hbar \phi^\dagger\dot{\phi}\psi
    - \phi^\dagger
    \int \d^3{\vect{x}'}H_{\vect{x},\vect{x}'}[n](t)\phi\psi
  \Biggr\},\\
  \I\hbar\dot{\psi}
  = E_\phi(z, t) \psi + \int\d{x}\d{y}\left\{
      \phi^\dagger\int \d^3{\vect{x}'}H_{\vect{x},\vect{x}'}[n](t)\phi\psi
  \right\}, \\
  E_\phi(z, t) = -\I\hbar\int\d{x}\d{y}\left\{\phi^\dagger\dot{\phi}\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\psi
  =
  \frac{-\hbar^2}{2m}\left(
    \psi\nabla^2\phi
    +\phi\nabla_z^2\psi
  \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).
\begin{gather*}
  \iint\d{x}\d{y}\;
    \phi^\dagger(x, y, z)
    \nabla_{z}\phi(x, y, z)\nabla_{z}\psi(z)\\
  \nabla_{z}\iint\d{x}\d{y}\; \phi^\dagger(x, y, z)\phi(x, y, z)
  = \iint\d{x}\d{y}\; \Bigl(
    \nabla_{z}\phi^\dagger(x, y, z)\phi(x, y, z) 
    + 
    \nabla_{z}\phi^\dagger(x, y, z)\phi(x, y, z)
  \Bigr)
\end{gather*}

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