---
jupytext:
  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
---

# Pancake

```{code-cell}
import mmf_setup;mmf_setup.nbinit()
import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL)
import matplotlib.pyplot as plt
import numpy as np
%load_ext autoreload
```

## Derivation

+++

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

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

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

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

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

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

+++

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

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

+++

The same result can be obtained by expanding the action:

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

For the usual energy density, this becomes

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

+++

## Gaussian axial profile

+++

A common strategy is to assume an axial Gaussian:

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

+++

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

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

+++

Ignore the term:

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

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

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

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

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

+++

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

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

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

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

+++

## Central Density

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

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

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

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

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

+++

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

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

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

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

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

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

+++

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

+++

## Energy density

+++

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

\begin{align}
    \frac{\partial \sigma^2}{\partial \phi^\dagger} &= \frac{2 \sqrt{2} a_\perp^4}{\sigma^2} \phi,\\
    \frac{\partial \sigma^4}{\partial \phi^\dagger} &= 4 \sqrt{2} a a_\perp^4 \phi,\\
    \frac{\partial \sigma^6}{\partial \phi^\dagger} &= 6 \sqrt{2} a a_\perp^4 \sigma^2 \phi.
\end{align}

+++

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

+++

## Demonstration

```{code-cell}
%autoreload
from gpe import disc

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

```{code-cell}
Nx, Ny = state.basis.Nxyz
x = state.x
n = state.get_density()
plt.plot(x, n[Nx//2, :])
```

```{code-cell}
state.evolve();
```

```{code-cell}

```