Chapter 4

Gross-Pitaevskii Equation

Mean-Field Approximation

For a dilute gas of $N$ weakly interacting bosons at $T = 0$, nearly all particles occupy a single quantum state. We write the many-body field operator as a classical field plus fluctuations:

$$\hat{\Psi}(\mathbf{r}, t) = \Psi(\mathbf{r}, t) + \delta\hat{\Psi}(\mathbf{r}, t)$$

where $\Psi(\mathbf{r}, t) = \langle\hat{\Psi}\rangle$ is the condensate wave function (order parameter), normalized so that $\int|\Psi|^2\,d^3r = N_0 \approx N$. The condensate density is $n_0(\mathbf{r}) = |\Psi(\mathbf{r})|^2$.

For a dilute system the interaction is characterized by a single parameter — the s-wave scattering length $a_s$. The two-body potential is replaced by a contact pseudopotential $V(\mathbf{r} - \mathbf{r}') = g\,\delta(\mathbf{r} - \mathbf{r}')$ with coupling constant:

$$g = \frac{4\pi\hbar^2 a_s}{m}$$

For $^{87}$Rb: $a_s \approx 5.3\,\text{nm}$, giving repulsive interactions ($g > 0$). For $^{7}$Li at certain field values,$a_s < 0$ (attractive), enabling bright solitons.

The Gross-Pitaevskii Equation

Minimizing the energy functional with respect to $\Psi^*$ subject to the constraint $\int|\Psi|^2 = N$, or equivalently applying the Heisenberg equation of motion to the field operator in the mean-field approximation, yields the time-dependent Gross-Pitaevskii equation:

$$i\hbar\frac{\partial\Psi}{\partial t} = \left[-\frac{\hbar^2\nabla^2}{2m} + V_{\text{ext}}(\mathbf{r}) + g|\Psi|^2\right]\Psi$$

This is a nonlinear Schrödinger equation. The three terms on the right represent:

◆Kinetic energy: $-\hbar^2\nabla^2/(2m)$ — quantum pressure that resists compression
◆External potential: $V_{\text{ext}}$ — magnetic or optical trapping potential
◆Mean-field interaction: $g|\Psi|^2$ — density-dependent self-interaction

For stationary states, $\Psi(\mathbf{r},t) = \psi(\mathbf{r})\,e^{-i\mu t/\hbar}$, the time-independent GP equation reads:

$$\mu\,\psi = \left[-\frac{\hbar^2\nabla^2}{2m} + V_{\text{ext}} + g|\psi|^2\right]\psi$$

Here $\mu$ is the chemical potential (Lagrange multiplier for particle number conservation), not the single-particle energy.

Healing Length

Consider a uniform condensate with density $n$ perturbed by a hard wall at $x=0$. The condensate density must vanish at the wall and recover to $n$ over a characteristic distance — the healing length:

$$\xi = \frac{\hbar}{\sqrt{2mng}} = \frac{1}{\sqrt{8\pi n a_s}}$$

This length scale emerges from balancing kinetic energy $\sim\hbar^2/(2m\xi^2)$ against interaction energy $\sim ng$. The density profile near the wall is:

$$n(x) = n\,\tanh^2\!\left(\frac{x}{\sqrt{2}\,\xi}\right)$$

For typical BEC experiments, $\xi \sim 0.1\text{--}1\,\mu\text{m}$, much larger than the interparticle spacing but smaller than the cloud size. The healing length also sets the vortex core size in rotating condensates.

Thomas-Fermi Approximation

When the number of atoms is large, the interaction energy dominates over the kinetic energy (except near the boundary). Dropping the Laplacian from the time-independent GP equation gives the Thomas-Fermi density profile:

$$n_{\text{TF}}(\mathbf{r}) = \frac{\mu_{\text{TF}} - V_{\text{ext}}(\mathbf{r})}{g}, \quad \text{where } V_{\text{ext}} < \mu_{\text{TF}}$$

For a harmonic trap $V = m\omega^2 r^2/2$, this gives an inverted parabola with Thomas-Fermi radius $R_{\text{TF}} = \sqrt{2\mu_{\text{TF}}/(m\omega^2)}$. The chemical potential is fixed by normalization:

$$\mu_{\text{TF}} = \frac{\hbar\omega}{2}\left(\frac{15Na_s}{a_{\text{ho}}}\right)^{2/5}$$

The Thomas-Fermi approximation breaks down within a distance $\sim\xi$ of the cloud edge, where the kinetic energy smooths the sharp density cutoff into the exponential tail characteristic of quantum tunneling.

Bogoliubov Transformation & Excitation Spectrum

To study excitations above the condensate, we write $\hat{\Psi} = (\sqrt{n_0} + \delta\hat{\Psi})\,e^{-i\mu t/\hbar}$ and linearize the Heisenberg equation of motion. The fluctuation operator is expanded as:

$$\delta\hat{\Psi} = \sum_k \left(u_k\,\hat{b}_k\,e^{i\mathbf{k}\cdot\mathbf{r} - i\omega_k t} - v_k\,\hat{b}_k^\dagger\,e^{-i\mathbf{k}\cdot\mathbf{r} + i\omega_k t}\right)$$

where $\hat{b}_k^\dagger$ creates a Bogoliubov quasiparticle. The amplitudes satisfy$u_k^2 - v_k^2 = 1$ and yield the celebrated Bogoliubov dispersion relation:

$$\hbar\omega_k = \sqrt{\varepsilon_k(\varepsilon_k + 2ng)}$$

where $\varepsilon_k = \hbar^2k^2/(2m)$ is the free-particle energy. Two important limits emerge:

Phonon Regime ($k\xi \ll 1$)

$$\omega_k \approx c\,k$$

Linear dispersion — sound waves with velocity $c = \sqrt{ng/m}$. Quasiparticles are collective density oscillations.

Free-Particle Regime ($k\xi \gg 1$)

$$\omega_k \approx \varepsilon_k/\hbar + ng/\hbar$$

Quadratic dispersion shifted by the mean-field energy $ng$. Quasiparticles are dressed single particles.

The crossover between these regimes occurs at $k \sim 1/\xi$, i.e., at the inverse healing length. The gapless, linear spectrum at low $k$ is a Goldstone mode associated with broken $U(1)$ gauge symmetry.

Sound Velocity & Landau Criterion

The speed of Bogoliubov sound is:

$$c = \sqrt{\frac{ng}{m}} = \sqrt{\frac{4\pi\hbar^2 n a_s}{m^2}}$$

This is also related to the compressibility: $mc^2 = n\,\partial\mu/\partial n$. The Landau critical velocity for superfluidity is:

$$v_c = \min_k\frac{\omega_k}{k} = c$$

Since the Bogoliubov spectrum is linear at low $k$ with slope $c$, the minimum of $\omega_k/k$ equals $c$. An object moving through the condensate at velocity $v < c$ cannot create excitations — it experiences no drag. This is the microscopic basis for superfluidity in weakly interacting Bose gases.

Soliton Solutions

The 1D GP equation supports exact soliton solutions that propagate without changing shape due to a balance between dispersion and nonlinearity.

Dark Solitons ($g > 0$)

$$\Psi(x,t) = \sqrt{n}\left[i\frac{v}{c} + \sqrt{1-\frac{v^2}{c^2}}\tanh\!\left(\frac{x - vt}{\xi_v}\right)\right]e^{-i\mu t/\hbar}$$

A density notch moving at velocity $v < c$ with width $\xi_v = \xi/\sqrt{1-v^2/c^2}$. At $v = 0$ the soliton is stationary with a complete density dip — a black soliton.

Bright Solitons ($g < 0$)

$$\Psi(x,t) = \frac{\sqrt{n_0}}{\cosh\!\left(\frac{x - vt}{\xi_b}\right)}e^{i(mvx - Et)/\hbar}$$

A localized density peak that exists without any external trap, held together by attractive interactions. Width $\xi_b = \hbar^2/(m|g|n_0)$.

Dark solitons have been observed experimentally in $^{87}$Rb condensates by phase-imprinting techniques, while bright solitons have been created in $^{7}$Li near a Feshbach resonance where $a_s$ can be tuned to negative values.

Time-Dependent GP & Condensate Dynamics

The time-dependent GP equation conserves the total energy, particle number, and momentum. Writing $\Psi = \sqrt{n}\,e^{i\phi}$, the GP equation decomposes into hydrodynamic form:

$$\frac{\partial n}{\partial t} + \nabla\cdot(n\mathbf{v}_s) = 0$$

$$m\frac{\partial \mathbf{v}_s}{\partial t} = -\nabla\left(\frac{1}{2}mv_s^2 + V_{\text{ext}} + gn - \frac{\hbar^2}{2m\sqrt{n}}\nabla^2\sqrt{n}\right)$$

where $\mathbf{v}_s = (\hbar/m)\nabla\phi$ is the superfluid velocity. The continuity equation and Euler-like equation are identical to irrotational, inviscid fluid dynamics except for the quantum pressure term $\propto \nabla^2\sqrt{n}$, which is important only on scales $\lesssim \xi$.

Key dynamical phenomena include: collective oscillations of trapped condensates (breathing and dipole modes), free expansion after trap release, interference of overlapping condensates, and the nucleation of vortices in stirred or rotating traps.

Detailed Derivation: GP Equation from Energy Functional

Step 1: The Energy Functional

The total energy of the condensate in the mean-field approximation is:

$$E[\Psi] = \int d^3r \left[\frac{\hbar^2}{2m}|\nabla\Psi|^2 + V_{\text{ext}}|\Psi|^2 + \frac{g}{2}|\Psi|^4\right]$$

The three terms are kinetic energy, potential energy, and mean-field interaction energy. The factor of $1/2$ in the interaction term avoids double-counting.

Step 2: Constrained Variation

We minimize $E[\Psi]$ subject to the constraint $\int|\Psi|^2\,d^3r = N$ by introducing the Lagrange multiplier $\mu$ and varying $\Psi^*$ independently:

$$\frac{\delta}{\delta\Psi^*}\left[E - \mu\int|\Psi|^2\,d^3r\right] = 0$$

$$-\frac{\hbar^2}{2m}\nabla^2\Psi + V_{\text{ext}}\Psi + g|\Psi|^2\Psi = \mu\Psi$$

This is the time-independent Gross-Pitaevskii equation. The Lagrange multiplier $\mu$is the chemical potential: $\mu = \partial E/\partial N$.

Step 3: Time-Dependent GP via the Action Principle

The time-dependent GP equation follows from the action:

$$S[\Psi,\Psi^*] = \int dt\,d^3r\left[\frac{i\hbar}{2}\left(\Psi^*\frac{\partial\Psi}{\partial t} - \Psi\frac{\partial\Psi^*}{\partial t}\right) - \mathcal{E}\right]$$

where $\mathcal{E}$ is the energy density. The Euler-Lagrange equation$\delta S/\delta\Psi^* = 0$ yields:

$$i\hbar\frac{\partial\Psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2\Psi + V_{\text{ext}}\Psi + g|\Psi|^2\Psi$$

Detailed Derivation: Bogoliubov Excitation Spectrum

Step 1: Linearization About the Uniform Condensate

Write $\Psi(\mathbf{r},t) = [\sqrt{n_0} + \delta\psi(\mathbf{r},t)]\,e^{-i\mu t/\hbar}$where $\mu = gn_0$ for a uniform condensate. Substituting into the time-dependent GP equation and keeping terms linear in $\delta\psi$:

$$i\hbar\frac{\partial\delta\psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2\delta\psi + gn_0(\delta\psi + \delta\psi^*)$$

The coupling between $\delta\psi$ and $\delta\psi^*$ is a crucial feature: it arises from the broken $U(1)$ symmetry of the condensate.

Step 2: Bogoliubov Transformation

Expand in plane waves: $\delta\psi = \sum_{\mathbf{k}} [u_k\,e^{i(\mathbf{k}\cdot\mathbf{r} - \omega_k t)} + v_k^*\,e^{-i(\mathbf{k}\cdot\mathbf{r} - \omega_k t)}]$. Substituting and equating coefficients of $e^{i\mathbf{k}\cdot\mathbf{r}}$ and $e^{-i\mathbf{k}\cdot\mathbf{r}}$:

$$\hbar\omega_k\,u_k = \varepsilon_k\,u_k + gn_0(u_k + v_k)$$

$$-\hbar\omega_k\,v_k = \varepsilon_k\,v_k + gn_0(u_k + v_k)$$

where $\varepsilon_k = \hbar^2k^2/(2m)$ is the free-particle energy.

Step 3: Solving the Eigenvalue Problem

Writing this as a matrix eigenvalue problem:

$$\begin{pmatrix} \varepsilon_k + gn_0 & gn_0 \\ -gn_0 & -\varepsilon_k - gn_0 \end{pmatrix}\begin{pmatrix} u_k \\ v_k \end{pmatrix} = \hbar\omega_k \begin{pmatrix} u_k \\ v_k \end{pmatrix}$$

The eigenvalues are found from $\det(M - \hbar\omega_k I) = 0$:

$$(\hbar\omega_k)^2 = (\varepsilon_k + gn_0)^2 - (gn_0)^2 = \varepsilon_k^2 + 2gn_0\varepsilon_k = \varepsilon_k(\varepsilon_k + 2gn_0)$$

This gives the Bogoliubov dispersion:

$$\boxed{\hbar\omega_k = \sqrt{\varepsilon_k(\varepsilon_k + 2gn_0)}}$$

The amplitudes are $u_k^2 = \frac{1}{2}\left(\frac{\varepsilon_k + gn_0}{\hbar\omega_k} + 1\right)$and $v_k^2 = \frac{1}{2}\left(\frac{\varepsilon_k + gn_0}{\hbar\omega_k} - 1\right)$, satisfying $u_k^2 - v_k^2 = 1$. The quantum depletion at $T = 0$ is$(N - N_0)/N = (8/3)\sqrt{na_s^3/\pi}$ (Lee-Huang-Yang result).

Derivation: Healing Length and Density Profile Near a Wall

Step 1: Dimensionless GP Equation in 1D

Consider a semi-infinite condensate bounded by a hard wall at $x = 0$. The time-independent GP equation in 1D (no trapping potential, $\mu = gn$ for bulk density $n$) is:

$$-\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2} + g|\psi|^2\psi = gn\,\psi$$

Write $\psi(x) = \sqrt{n}\,f(\tilde{x})$ with $\tilde{x} = x/(\sqrt{2}\,\xi)$where $\xi = \hbar/\sqrt{2mgn}$. After substitution:

$$-\frac{d^2 f}{d\tilde{x}^2} + f^3 - f = 0$$

Step 2: Solving the ODE

Multiply by $df/d\tilde{x}$ and integrate: $\frac{1}{2}(df/d\tilde{x})^2 = \frac{1}{4}(f^2 - 1)^2$. Taking the square root with $f(0) = 0$ and $f(\infty) = 1$:

$$\frac{df}{d\tilde{x}} = \frac{1}{\sqrt{2}}(1 - f^2)$$

This is separable: $\int df/(1 - f^2) = \tilde{x}/\sqrt{2}$, giving$\text{arctanh}(f) = \tilde{x}/\sqrt{2}$, so $f = \tanh(\tilde{x}/\sqrt{2})$. Returning to physical units:

$$n(x) = n\,\tanh^2\!\left(\frac{x}{\sqrt{2}\,\xi}\right)$$

The density recovers from zero to $n$ over a characteristic scale $\sim\xi$, confirming $\xi$ as the healing length. The same profile describes the core structure of a dark soliton at rest.

Historical Context

The equation was derived independently by Eugene P. Gross(1961) at Brandeis University and Lev P. Pitaevskii (1961) at the Institute for Physical Problems in Moscow. Both were motivated by the need for a self-consistent mean-field description of a weakly interacting Bose gas that could describe spatially inhomogeneous phenomena such as vortices and boundary effects.

The theoretical foundations were laid earlier by Nikolay Bogoliubov(1947), who developed the transformation that bears his name to treat the excitation spectrum of a weakly interacting Bose gas. Bogoliubov's work showed that interactions convert the quadratic free-particle spectrum into a linear (phonon) spectrum at low momenta, providing the microscopic basis for superfluidity in the Landau sense.

The GP equation gained enormous practical importance after the experimental realization of BEC in 1995. It successfully describes the ground-state density profiles, collective oscillation frequencies, vortex nucleation, soliton dynamics, and interference patterns of dilute atomic condensates. Extensions such as the stochastic GP equation (including thermal fluctuations) and the multi-component GP equations (for spinor condensates) remain active areas of theoretical development.

Applications of the Gross-Pitaevskii Equation

Trapped condensate modeling: The GP equation quantitatively predicts the density profile, chemical potential, and Thomas-Fermi radius of BECs in magnetic and optical traps. It is the standard workhorse for experimental design and data interpretation in ultracold atom laboratories worldwide.
Nonlinear atom optics and solitons: The GP equation describes dark and bright solitons in BECs, analogous to optical solitons in nonlinear fibers. Soliton collisions, soliton-vortex interactions, and soliton trains have been experimentally observed and are used to probe nonlinear dynamics in quantum fluids.
Vortex lattice formation: Numerical solutions of the rotating GP equation predict the critical rotation frequency for vortex nucleation, the triangular vortex lattice structure, and the transition to the quantum Hall regime at extreme rotation rates.
Spinor condensates and topological excitations: Multi-component GP equations describe spinor BECs (e.g., spin-1 $^{87}$Rb or $^{23}$Na), which support rich topological structures including half-quantum vortices, skyrmions, monopoles, and knot solitons.
Analogue gravity: The GP equation for a flowing condensate provides a laboratory realization of analogue spacetime metrics. Phonons propagating in a condensate with a supersonic flow region experience an effective acoustic horizon, enabling studies of analogue Hawking radiation.

Python Simulation: GP Ground State

Imaginary-time evolution of the 1D GP equation in a harmonic trap using the split-step Fourier method. The converged ground state density is compared with the Thomas-Fermi approximation.

GP Ground State via Imaginary-Time Evolution

Python

Split-step Fourier solution of the 1D Gross-Pitaevskii equation with Thomas-Fermi comparison

script.py126 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────────────
# Gross-Pitaevskii Equation: 1D Harmonic Trap
# Imaginary-time evolution for the ground state
# Comparison with Thomas-Fermi approximation
# ──────────────────────────────────────────────────────

# --- Dimensionless units: lengths in a_ho, energies in hbar*omega ---
# a_ho = sqrt(hbar / (m * omega))
# GPE: -0.5 d^2/dx^2 psi + 0.5 x^2 psi + g_1d |psi|^2 psi = mu psi

N_particles = 500       # number of atoms
a_s_over_aho = 0.01     # scattering length / oscillator length
g_1d = 2.0 * N_particles * a_s_over_aho  # effective 1D coupling

# Spatial grid
Nx = 512
x_max = 12.0
x = np.linspace(-x_max, x_max, Nx)
dx = x[1] - x[0]

# Harmonic trap
V_trap = 0.5 * x**2

# Kinetic energy operator in momentum space
dk = 2.0 * np.pi / (Nx * dx)
k = np.fft.fftfreq(Nx, d=dx) * 2.0 * np.pi
T_k = 0.5 * k**2

# Initial guess: Gaussian
psi = np.exp(-x**2 / 4.0).astype(complex)
norm = np.trapezoid(np.abs(psi)**2, x)
psi /= np.sqrt(norm)

# Imaginary-time evolution parameters
dt_imag = 0.002
n_steps = 8000

# Split-step method in imaginary time
for step in range(n_steps):
    # Half-step potential + interaction
    psi *= np.exp(-0.5 * dt_imag * (V_trap + g_1d * np.abs(psi)**2))

# Full-step kinetic (in k-space)
    psi_k = np.fft.fft(psi)
    psi_k *= np.exp(-dt_imag * T_k)
    psi = np.fft.ifft(psi_k)

# Half-step potential + interaction
    psi *= np.exp(-0.5 * dt_imag * (V_trap + g_1d * np.abs(psi)**2))

# Renormalize
    norm = np.trapezoid(np.abs(psi)**2, x)
    psi /= np.sqrt(norm)

# Compute chemical potential: mu = <psi|H|psi>
density = np.abs(psi)**2
kinetic = np.trapezoid(0.5 * np.abs(np.gradient(psi, dx))**2, x).real
potential = np.trapezoid(V_trap * density, x)
interaction = np.trapezoid(g_1d * density**2, x)
mu = kinetic + potential + interaction
energy = kinetic + potential + 0.5 * np.trapezoid(g_1d * density**2, x)

# Thomas-Fermi approximation: n_TF(x) = (mu_TF - V(x)) / g_1d
mu_TF = 0.5 * (1.5 * g_1d)**( 2.0/3.0)
R_TF = np.sqrt(2.0 * mu_TF)
n_TF = np.where(mu_TF > V_trap, (mu_TF - V_trap) / g_1d, 0.0)

# Healing length at center
n_center = density[Nx // 2]
xi = 1.0 / np.sqrt(2.0 * g_1d * n_center) if n_center > 0 else np.inf

# ---- Plotting ----
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#1e3a5f')

# Left panel: density profile
axes[0].fill_between(x, density, alpha=0.3, color='#06b6d4')
axes[0].plot(x, density, color='#22d3ee', linewidth=2.5, label='GP ground state')
axes[0].plot(x, n_TF, color='#f59e0b', linewidth=2, linestyle='--', label='Thomas-Fermi')
axes[0].axhline(0, color='#334155', linewidth=0.5)
axes[0].axvline(-R_TF, color='#f59e0b', linewidth=0.8, linestyle=':', alpha=0.6)
axes[0].axvline(R_TF, color='#f59e0b', linewidth=0.8, linestyle=':', alpha=0.6)
axes[0].set_xlabel('x / a_ho', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('Density |psi(x)|^2', color='#94a3b8', fontsize=12)
axes[0].set_title('Condensate Density in Harmonic Trap', color='white', fontsize=14)
axes[0].legend(fontsize=10, facecolor='#0f172a', edgecolor='#1e3a5f', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].set_xlim(-8, 8)

# Right panel: effective potential and energy levels
V_eff = V_trap + g_1d * density
axes[1].plot(x, V_trap, color='#64748b', linewidth=1.5, linestyle='--', label='V_trap = x^2/2')
axes[1].plot(x, V_eff, color='#22d3ee', linewidth=2, label='V_eff = V + g|psi|^2')
axes[1].axhline(mu, color='#f43f5e', linewidth=1.5, linestyle='-.', label=f'mu = {mu:.3f}')
axes[1].axhline(mu_TF, color='#f59e0b', linewidth=1, linestyle=':', label=f'mu_TF = {mu_TF:.3f}')
axes[1].fill_between(x, V_trap, V_eff, alpha=0.15, color='#06b6d4')
axes[1].set_xlabel('x / a_ho', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('Energy / (hbar omega)', color='#94a3b8', fontsize=12)
axes[1].set_title('Effective Potential & Chemical Potential', color='white', fontsize=14)
axes[1].legend(fontsize=9, facecolor='#0f172a', edgecolor='#1e3a5f', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')
axes[1].set_xlim(-8, 8)
axes[1].set_ylim(0, mu * 2.5)

plt.tight_layout()
plt.savefig('gp_ground_state.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

print(f"1D coupling constant g_1d = {g_1d:.4f}")
print(f"Chemical potential (GP): mu = {mu:.4f} hbar*omega")
print(f"Chemical potential (TF): mu_TF = {mu_TF:.4f} hbar*omega")
print(f"Thomas-Fermi radius: R_TF = {R_TF:.4f} a_ho")
print(f"Healing length at center: xi = {xi:.4f} a_ho")
print(f"Total energy E = {energy.real:.4f} hbar*omega")
print(f"Kinetic / Total = {(kinetic/energy).real:.4f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Simulation: Bogoliubov Spectrum

Computation of the Bogoliubov excitation spectrum $\omega(k)$ showing the crossover from the phonon regime ($\omega \propto k$) at low momentum to the free-particle regime ($\omega \propto k^2$) at high momentum. Bogoliubov amplitudes $u_k, v_k$ are also computed.

Bogoliubov Excitation Spectrum

Fortran

Phonon-to-free-particle crossover in the quasiparticle dispersion relation

program.f9091 lines

program bogoliubov_spectrum
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: Nk = 200

real(dp), parameter :: pi = 4.0_dp * atan(1.0_dp)

! ── Physical parameters (dimensionless units) ──
  ! Energy in units of ng, momentum in units of 1/xi
  ! where xi = hbar / sqrt(2*m*n*g) is the healing length
  ! and ng is the interaction energy per particle

real(dp) :: k, dk, k_max
  real(dp) :: eps_k, E_bog, E_free, c_sound
  real(dp) :: k_arr(Nk), E_arr(Nk), Efree_arr(Nk), Ephonon_arr(Nk)
  real(dp) :: u_k, v_k, depletion
  integer  :: i

! Sound velocity in these units: c = 1 (in units of sqrt(ng/m))
  c_sound = 1.0_dp

k_max = 6.0_dp
  dk = k_max / real(Nk - 1, dp)

write(*,'(A)') '=== Bogoliubov Excitation Spectrum ==='
  write(*,'(A)') 'Units: energy in ng, momentum in 1/xi'
  write(*,'(A)') ''
  write(*,'(A)') 'Bogoliubov dispersion: E(k) = sqrt(eps_k * (eps_k + 2*ng))'
  write(*,'(A)') 'where eps_k = hbar^2 k^2 / (2m)'
  write(*,'(A)') ''

depletion = 0.0_dp

write(*,'(A6, A14, A14, A14, A14, A14)') &
       'k/k_xi', 'eps_k', 'E_Bog', 'E_phonon', 'u_k', 'v_k'
  write(*,'(A)') repeat('-', 84)

do i = 1, Nk
    k = real(i - 1, dp) * dk
    if (i == 1) k = 0.01_dp  ! avoid k=0

k_arr(i) = k

! Free particle energy (in units of ng):
    ! eps_k = (k*xi)^2 * ng = k^2  (since k is in units of 1/xi, and eps in units of ng)
    ! Actually: eps_k = hbar^2 k^2/(2m). In units where ng=1 and 1/xi=1:
    ! xi = hbar/sqrt(2mng), so hbar^2/(2m) = xi^2 * ng = 1
    ! Hence eps_k = k^2
    eps_k = k**2

! Bogoliubov dispersion
    ! E = sqrt(eps_k * (eps_k + 2))  [since ng = 1 in our units]
    E_bog = sqrt(eps_k * (eps_k + 2.0_dp))

! Free particle limit
    E_free = eps_k

E_arr(i)      = E_bog
    Efree_arr(i)  = E_free + 1.0_dp  ! shift: free particle is eps_k + ng
    Ephonon_arr(i) = c_sound * k      ! phonon: E = c*k, c = sqrt(ng/m) = 1 in our units

! Bogoliubov amplitudes
    u_k = sqrt(0.5_dp * (eps_k / E_bog + 1.0_dp))
    v_k = -sqrt(0.5_dp * (eps_k / E_bog - 1.0_dp))

! Quantum depletion integrand: v_k^2
    if (i > 1) then
      depletion = depletion + v_k**2 * dk
    end if

if (mod(i-1, 20) == 0 .or. i == Nk) then
      write(*,'(F6.2, F14.6, F14.6, F14.6, F14.6, F14.6)') &
           k, eps_k, E_bog, c_sound * k, u_k, v_k
    end if
  end do

write(*,'(A)') ''
  write(*,'(A)') '=== Key Results ==='
  write(*,'(A,F10.6)') 'Sound velocity c (units sqrt(ng/m)): ', c_sound
  write(*,'(A,F10.6)') 'Quantum depletion (1D integral):     ', depletion
  write(*,'(A)') ''
  write(*,'(A)') 'Low-k regime (phonon): E ~ c*k   (linear, sound waves)'
  write(*,'(A)') 'High-k regime (free):  E ~ k^2   (quadratic, free particles)'
  write(*,'(A)') 'Crossover at k ~ 1/xi where eps_k ~ ng'
  write(*,'(A)') ''
  write(*,'(A)') 'The Bogoliubov spectrum is gapless (Goldstone mode)'
  write(*,'(A)') 'ensuring superfluidity via the Landau criterion:'
  write(*,'(A)') 'v_c = min(E_k/k) = c > 0'

end program bogoliubov_spectrum

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Chapter 3: Two-Fluid Model Chapter 5: Quantized Vortices →