Part II, Chapter 5

Quantized Vortices

Topological defects in superfluids: quantized circulation, vortex dynamics, and the Kosterlitz-Thouless transition

5.1 Circulation Quantization

In a superfluid described by a macroscopic wavefunction $\Psi(\mathbf{r}) = \sqrt{n_s}\,e^{i\phi(\mathbf{r})}$, the superfluid velocity is determined by the gradient of the phase:

$$\mathbf{v}_s = \frac{\hbar}{m}\nabla\phi$$

Since the wavefunction must be single-valued, the phase $\phi$ can only change by integer multiples of $2\pi$ around any closed loop. This immediately gives the quantization of circulation:

Quantized Circulation

$$\kappa = \oint \mathbf{v}_s \cdot d\mathbf{l} = n\,\frac{h}{m}, \qquad n = 0, \pm 1, \pm 2, \ldots$$

The quantum of circulation is $h/m \approx 9.97 \times 10^{-4}\;\text{cm}^2/\text{s}$ for $^4$He.

The integer $n$ is a topological winding numberthat cannot change continuously — it is robust against smooth perturbations, making vortices genuinely topological defects.

Key Insight: The irrotationality constraint $\nabla \times \mathbf{v}_s = 0$ holds everywhere except on singular vortex lines where the phase is undefined. Vortices are the only mechanism by which a superfluid can carry angular momentum.

5.2 Vortex Line Structure

For a singly-quantized vortex ($n = 1$) aligned along the $z$-axis, the velocity field in cylindrical coordinates is purely azimuthal:

Vortex Velocity Field

$$v_s(r) = \frac{\hbar}{mr}$$

The $1/r$ divergence at $r = 0$ signals the breakdown of the superfluid description at the vortex core.

The kinetic energy density diverges as $r \to 0$, so the superfluid density must vanish at the vortex core to regularize the energy.

Vortex Core Structure from the Gross-Pitaevskii Equation

Writing $\Psi(r,\theta) = \sqrt{n_0}\,f(r)\,e^{in\theta}$, the GP equation gives the dimensionless radial profile:

$$-\frac{\xi^2}{r}\frac{d}{dr}\!\left(r\frac{df}{dr}\right) + \frac{n^2\xi^2}{r^2}f + f^3 - f = 0$$

Here $\xi = \hbar/\sqrt{2mgn_0}$ is the healing length (or coherence length).

Boundary conditions: $f(0) = 0$ and $f \to 1$ as $r \to \infty$. Near the core $f(r) \sim (r/\xi)^{|n|}$, recovering bulk density over a distance $\sim \xi$. In $^4$He, $\xi \approx 1\;\text{\AA}$; in dilute BECs, $\xi \sim 100\text{--}500\;\text{nm}$.

$n = 1$ Vortex

Energetically stable. The density vanishes linearly as $|\Psi|^2 \propto r^2$ near the core. This is the fundamental topological excitation of the superfluid.

$n \geq 2$ Vortex

Energetically unstable toward splitting into $|n|$ singly-quantized vortices, since $E \propto n^2$ while $|n|$ single vortices carry energy $\propto |n|$.

5.3 Energy of a Vortex Line

The kinetic energy of a vortex line of length $L$ in a cylindrical container of radius $R$ is obtained by integrating $\frac{1}{2}\rho_s v_s^2$ over the volume, with a lower cutoff at the vortex core radius $\sim \xi$:

Vortex Line Energy

$$E_v = \pi n_s \frac{\hbar^2}{m^2}\,L\,\ln\!\left(\frac{R}{\xi}\right)$$

For a vortex with winding number $n$: $E_v \propto n^2 \ln(R/\xi)$.

The logarithmic dependence on $R/\xi$ has key consequences:

An isolated vortex in an infinite 2D system has infinite energy; vortices exist only as bound pairs, or are stabilized by finite size or rotation.
A $n = 2$ vortex costs $4\ln(R/\xi)$ while two $n = 1$ vortices cost $2\ln(R/\xi)$ — multiply-quantized vortices are unstable.
The energy per unit length defines a line tension resisting bending (Kelvin waves).

A core energy of order $\pi n_0 g \xi^2 L$ is typically negligible for $R \gg \xi$.

5.4 Vortex-Antivortex Pairs & the Kosterlitz-Thouless Transition

In two-dimensional superfluids, isolated vortices have logarithmically divergent energy. However, a vortex-antivortex pair (vortex with $n = +1$ and antivortex with $n = -1$) has a finite interaction energy:

$$E_{\text{pair}}(d) = 2\pi n_s \frac{\hbar^2}{m^2}\,\ln\!\left(\frac{d}{\xi}\right)$$

where $d$ is the separation between the vortex and antivortex.

The logarithmic interaction is analogous to 2D Coulomb charges. At low temperature, pairs are tightly bound; as temperature rises, a remarkable phase transition occurs:

The Berezinskii-Kosterlitz-Thouless (BKT) Transition

The entropy of placing a vortex in a system of area $A$ is $S \sim k_B \ln(A/\xi^2)$. The free energy of a single vortex is:

$$F = E - TS = \left(\pi n_s \frac{\hbar^2}{m^2} - 2k_B T\right)\ln\!\left(\frac{R}{\xi}\right)$$

Free vortices proliferate when $F < 0$, giving the BKT transition temperature:

BKT Critical Temperature

$$k_B T_{\text{BKT}} = \frac{\pi}{2}\,n_s\,\frac{\hbar^2}{m^2}$$

Below $T_{\text{BKT}}$, bound pairs give quasi-long-range order with algebraic correlations. Above $T_{\text{BKT}}$, pairs unbind, destroying superfluidity with a universal density jump: $n_s(T_{\text{BKT}}^-) = 2mk_BT_{\text{BKT}}/(\pi\hbar^2)$.

Key Insight: The BKT transition is topological — driven by vortex proliferation rather than a local order parameter. It has been observed in superfluid helium films, ultracold gases, and superconducting thin films.

5.5 Rotating Superfluids & the Feynman Relation

Since $\nabla \times \mathbf{v}_s = 0$ away from vortex cores, a superfluid cannot achieve solid-body rotation directly. Instead, it mimics rotation by nucleating an array of quantized vortices. The Feynman relation gives the required areal vortex density:

Feynman Relation

$$n_v = \frac{2m\Omega}{h} = \frac{2\Omega}{\kappa}$$

where $n_v$ is the number of vortices per unit area and $\kappa = h/m$ is the circulation quantum.

This follows from matching the coarse-grained vorticity to $2\Omega$. At typical rotation rates ($\Omega \sim 1\;\text{rad/s}$), the inter-vortex spacing is $\ell \sim 0.2\;\text{mm}$ in superfluid helium.

Historical Note: Feynman predicted vortex arrays in 1955, confirmed by Vinen (1961) via second-sound attenuation and directly imaged by Yarmchuk, Gordon, and Packard (1979) using ion trapping.

5.6 Vortex Lattice (Abrikosov Lattice)

Multiple vortices in a rotating superfluid arrange into a triangular lattice that minimizes interaction energy — the superfluid analogue of the Abrikosov lattice in type-II superconductors. The interaction energy between two parallel vortex lines at separation $d$ is:

$$E_{\text{int}}(d) = \pi n_s \frac{\hbar^2}{m^2}\,L\,\ln\!\left(\frac{R}{d}\right), \qquad d \gg \xi$$

The logarithmic repulsion drives the formation of a triangular lattice with lattice constant $a = (2/\sqrt{3}\,n_v)^{1/2}$. Lattices with hundreds of vortices have been observed in rapidly rotating BECs. At extreme rotation, vortex cores overlap and the system enters the lowest Landau level regime with connections to fractional quantum Hall physics.

Superfluid Helium

Vortex core size $\xi \sim 1\;\text{\AA}$. Lattice spacing at typical rotation$a \sim 0.1\text{--}1\;\text{mm}$. Ratio $a/\xi \sim 10^6$: extremely dilute vortex array with well-separated cores.

Atomic BEC

Vortex core size $\xi \sim 0.2\;\mu\text{m}$. Lattice spacing$a \sim 5\text{--}20\;\mu\text{m}$. Ratio $a/\xi \sim 25\text{--}100$: vortex cores can be directly imaged via absorption imaging.

5.7 Vortex Dynamics

Magnus Force

A vortex moving with velocity $\mathbf{v}_L$ through a superfluid with background velocity $\mathbf{v}_s$ experiences a Magnus forceper unit length:

$$\mathbf{f}_M = \rho_s \boldsymbol{\kappa} \times (\mathbf{v}_L - \mathbf{v}_s)$$

where $\boldsymbol{\kappa} = \kappa\,\hat{z}$ is the quantized circulation vector.

Without dissipation, $\mathbf{v}_L = \mathbf{v}_s$: the vortex moves with the local superfluid flow. Same-sign vortex pairs orbit their center; a vortex-antivortex pair translates perpendicular to the line joining them.

Mutual Friction

In real superfluids at finite temperature, the vortex core interacts with the normal fluid component, giving rise to mutual friction. The equation of motion for a vortex line becomes:

$$\rho_s \boldsymbol{\kappa} \times (\mathbf{v}_L - \mathbf{v}_s) = -\gamma_0\,\hat{\boldsymbol{\kappa}} \times [\hat{\boldsymbol{\kappa}} \times (\mathbf{v}_L - \mathbf{v}_n)] - \gamma_0'\,\hat{\boldsymbol{\kappa}} \times (\mathbf{v}_L - \mathbf{v}_n)$$

Here $\gamma_0, \gamma_0'$ are dissipative and reactive friction coefficients and$\mathbf{v}_n$ is the normal fluid velocity. Mutual friction governs the decay of superfluid turbulence and normal-superfluid coupling in rotating helium.

Kelvin Waves

Small helical displacements of a vortex line propagate as Kelvin waves with dispersion $\omega(k) = (\kappa k^2 / 4\pi)[\ln(1/k\xi) + c]$. These waves drive the energy cascade in quantum turbulence, transferring energy to small scales where it is radiated as phonons.

Detailed Derivation: Energy of a Quantized Vortex Line

Step 1: Kinetic Energy Density

For a singly-quantized vortex ($n = 1$) aligned along the $z$-axis, the velocity field is $v_s(r) = \hbar/(mr)$. The kinetic energy per unit length in a cylindrical shell of radius $r$ and thickness $dr$ is:

$$dE = \frac{1}{2}\rho_s\,v_s^2\cdot 2\pi r\,dr = \frac{1}{2}\rho_s\left(\frac{\hbar}{mr}\right)^2 2\pi r\,dr = \frac{\pi\rho_s\hbar^2}{m^2}\,\frac{dr}{r}$$

Step 2: Integration Over the Volume

Integrating from the vortex core radius $\xi$ (where the superfluid density vanishes) to the container radius $R$:

$$\frac{E_v}{L} = \frac{\pi\rho_s\hbar^2}{m^2}\int_{\xi}^{R}\frac{dr}{r} = \frac{\pi\rho_s\hbar^2}{m^2}\,\ln\!\left(\frac{R}{\xi}\right)$$

For a vortex with winding number $n$, $v_s = n\hbar/(mr)$, so$E_v/L = \pi n^2\rho_s\hbar^2\ln(R/\xi)/m^2$. The energy scales as $n^2$, making multiply-quantized vortices unstable.

Step 3: Instability of Multiply-Quantized Vortices

Compare the energy of a single $n = 2$ vortex with two $n = 1$ vortices separated by distance $d$. The $n = 2$ vortex has energy (per unit length):

$$E_{n=2} = 4\,\frac{\pi\rho_s\hbar^2}{m^2}\ln\!\left(\frac{R}{\xi}\right)$$

Two $n = 1$ vortices contribute self-energies plus an interaction energy:

$$E_{2\times(n=1)} = 2\,\frac{\pi\rho_s\hbar^2}{m^2}\ln\!\left(\frac{R}{\xi}\right) + \frac{\pi\rho_s\hbar^2}{m^2}\ln\!\left(\frac{R}{d}\right)$$

The difference is $E_{n=2} - E_{2\times(n=1)} = (\pi\rho_s\hbar^2/m^2)[2\ln(R/\xi) - \ln(R/d)] > 0$for all $d > \xi$, confirming the instability. Physically, vortex splitting is driven by the $n^2$ scaling of the self-energy.

Detailed Derivation: The BKT Transition Temperature

Step 1: Free Energy of an Isolated Vortex

The energy of a single vortex in a 2D system of area $A = \pi R^2$ is:

$$E_v = \pi n_s \frac{\hbar^2}{m^2}\,\ln\!\left(\frac{R}{\xi}\right)$$

where $n_s$ is the 2D superfluid density (mass per unit area). The vortex core can be placed anywhere in the area $A$, with $A/(\pi\xi^2)$ distinguishable positions, giving an entropy:

$$S_v = k_B\ln\!\left(\frac{A}{\pi\xi^2}\right) = 2k_B\ln\!\left(\frac{R}{\xi}\right)$$

Step 2: Free Energy Balance

The free energy of an isolated vortex is:

$$F_v = E_v - TS_v = \left(\pi n_s\frac{\hbar^2}{m^2} - 2k_BT\right)\ln\!\left(\frac{R}{\xi}\right)$$

In the thermodynamic limit ($R \to \infty$), free vortices are thermodynamically favourable when $F_v < 0$, i.e.:

$$\pi n_s\frac{\hbar^2}{m^2} < 2k_BT$$

Step 3: The Universal Jump and Critical Temperature

The transition occurs when the free energy changes sign, giving:

$$k_BT_{\text{BKT}} = \frac{\pi}{2}\,n_s\,\frac{\hbar^2}{m^2}$$

Equivalently, using the 2D superfluid density $\rho_s^{(2D)} = n_s$ (mass per area) and the circulation quantum $\kappa = h/m$:

$$\frac{\rho_s^{(2D)}(T_{\text{BKT}}^-)}{T_{\text{BKT}}} = \frac{2m^2 k_B}{\pi\hbar^2}$$

This is the Nelson-Kosterlitz universal jump (1977): the superfluid density drops discontinuously to zero at $T_{\text{BKT}}$ from a universal value determined solely by the particle mass and temperature. This remarkable prediction has been confirmed experimentally in thin helium films and in quasi-2D ultracold gases.

Derivation: Feynman Relation for Vortex Density

Step 1: Coarse-Grained Vorticity

A superfluid rotating at angular velocity $\Omega$ mimics solid-body rotation by distributing quantized vortices uniformly. Each vortex carries circulation $\kappa = h/m$. Coarse-graining over a region containing many vortices, the average vorticity is:

$$\langle\nabla\times\mathbf{v}_s\rangle = n_v\,\kappa\,\hat{z}$$

where $n_v$ is the areal vortex density (number per unit area perpendicular to the rotation axis).

Step 2: Matching to Solid-Body Rotation

For solid-body rotation, $\mathbf{v} = \boldsymbol{\Omega}\times\mathbf{r}$, so$\nabla\times\mathbf{v} = 2\boldsymbol{\Omega}$. Matching the coarse-grained vorticity:

$$n_v\kappa = 2\Omega \quad \Longrightarrow \quad n_v = \frac{2\Omega}{\kappa} = \frac{2m\Omega}{h}$$

The inter-vortex spacing in a triangular lattice is:

$$\ell = \left(\frac{1}{n_v}\right)^{1/2} = \sqrt{\frac{\kappa}{2\Omega}} = \sqrt{\frac{h}{2m\Omega}}$$

For He-4 at $\Omega = 1$ rad/s: $\ell = \sqrt{9.97\times 10^{-8}/(2\times 1)} \approx 0.22$ mm. For $^{87}$Rb at $\Omega = 2\pi\times 50$ Hz: $\ell \approx 5\;\mu$m, comparable to the healing length, making vortex cores resolvable by direct imaging.

Historical Context

The concept of quantized circulation in superfluids was introduced by Lars Onsager in a brief remark at a 1949 conference and independently by Richard Feynman in his 1955 paper on the application of quantum mechanics to liquid helium. Feynman provided the detailed physical picture of quantized vortex lines, predicted the formation of vortex arrays in rotating superfluids, and estimated the critical angular velocity for vortex nucleation.

William Vinen (1961) provided the first experimental evidence for quantized circulation by measuring the circulation around a fine wire suspended in rotating He-II using its vibration frequency. He confirmed that the circulation was quantized in units of$h/m = 9.97\times 10^{-4}$ cm$^2$/s. The direct imaging of vortex arrays was achieved by Yarmchuk, Gordon, and Packard (1979) using trapped ions as markers for vortex positions in rotating He-II.

The Berezinskii-Kosterlitz-Thouless theory was developed in stages: Vadim Berezinskii (1971) in the Soviet Union and J. Michael Kosterlitz and David Thouless (1973) in the UK independently recognized that vortex-antivortex unbinding could drive a phase transition in two dimensions without conventional symmetry breaking. Kosterlitz and Thouless shared half of the 2016 Nobel Prize in Physics (with Duncan Haldane) for “theoretical discoveries of topological phase transitions.” The BKT transition has been confirmed in thin helium films, 2D superconductors, and ultracold atomic gases, making it one of the most broadly validated predictions in theoretical physics.

Applications of Quantized Vortices

Type-II superconductor technology: The magnetic flux vortices (Abrikosov vortices) in type-II superconductors are the direct electromagnetic analogue of superfluid vortices. Understanding vortex pinning, creep, and dynamics is essential for the design of high-current superconducting magnets, cables, and energy storage systems.
Quantum turbulence: Tangles of quantized vortex lines in superfluid helium provide a minimal model of turbulence with quantized degrees of freedom. The Kolmogorov $k^{-5/3}$ energy spectrum has been observed in quantum turbulence, suggesting universal features shared with classical turbulence. This research informs our understanding of turbulence in neutron stars and early-universe cosmology.
Rapidly rotating BECs and quantum Hall physics: At extreme rotation rates where the vortex density approaches the particle density, the system enters the lowest Landau level regime. This provides a route to creating bosonic fractional quantum Hall states and studying topological order with neutral atoms.
Neutron star interiors: The superfluid neutrons in neutron star cores are expected to form vast arrays of quantized vortex lines. Vortex pinning to the nuclear lattice and subsequent unpinning events (vortex avalanches) are believed to cause pulsar glitches — sudden spin-up events observed in rotating neutron stars.
Topological quantum computation: Non-Abelian vortices (such as those in the A-phase of superfluid $^3$He or in $p + ip$superconductors) support Majorana zero modes that could serve as topologically protected qubits, resistant to local decoherence. Braiding operations on these vortices implement fault-tolerant quantum gates.

Interactive Simulations

Vortex Velocity Field, Density Profile & Vortex Lattice

Python

script.py124 lines

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))
fig.patch.set_facecolor('#0a0a0a')

# --- Panel 1: Velocity field of a single vortex ---
ax = axes[0]
ax.set_facecolor('#0a0a0a')

N = 200
x = np.linspace(-5, 5, N)
y = np.linspace(-5, 5, N)
X, Y = np.meshgrid(x, y)
R = np.sqrt(X**2 + Y**2)
R[R < 0.1] = 0.1  # avoid singularity

# velocity components (hbar/m = 1 units)
vx = -Y / R**2
vy = X / R**2
speed = np.sqrt(vx**2 + vy**2)

# Plot speed as background color
im = ax.pcolormesh(X, Y, speed, cmap='inferno', shading='auto',
                   norm=LogNorm(vmin=0.05, vmax=5))
# Streamlines
seed_points = np.array([[np.cos(t)*r, np.sin(t)*r]
                        for r in [0.5, 1.0, 1.5, 2.5, 4.0]
                        for t in [0]])
ax.streamplot(X, Y, vx, vy, color='cyan', linewidth=0.7,
              density=1.5, arrowsize=1.2, start_points=seed_points)
ax.set_title('Single Vortex Velocity Field', color='cyan', fontsize=13, fontweight='bold')
ax.set_xlabel('x / \u03be', color='gray')
ax.set_ylabel('y / \u03be', color='gray')
ax.tick_params(colors='gray')
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.set_aspect('equal')
cb = plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
cb.set_label('|v| (\u0127/m\u03be)', color='gray')
cb.ax.tick_params(colors='gray')

# --- Panel 2: Density profile through vortex core ---
ax2 = axes[1]
ax2.set_facecolor('#0a0a0a')

r_vals = np.linspace(0, 8, 500)
# Approximate GP solution for n=1: f(r) ~ r/sqrt(r^2 + 2*xi^2)
xi = 1.0
f1 = r_vals / np.sqrt(r_vals**2 + 2*xi**2)
# n=2 vortex
f2 = r_vals**2 / np.sqrt(r_vals**4 + 8*xi**4)

ax2.plot(r_vals, f1**2, color='cyan', linewidth=2.5, label='n = 1 vortex')
ax2.plot(r_vals, f2**2, color='#ff6b6b', linewidth=2.5, linestyle='--', label='n = 2 vortex')
ax2.axhline(y=1, color='gray', linestyle=':', alpha=0.5, label='Bulk density')
ax2.axvline(x=xi, color='cyan', linestyle=':', alpha=0.4)
ax2.annotate('\u03be', xy=(xi, 0.02), color='cyan', fontsize=14, ha='center')
ax2.set_xlabel('r / \u03be', color='gray', fontsize=12)
ax2.set_ylabel('|\u03a8|\u00b2 / n\u2080', color='gray', fontsize=12)
ax2.set_title('Vortex Core Density Profile', color='cyan', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10, facecolor='#1a1a1a', edgecolor='cyan', labelcolor='lightgray')
ax2.tick_params(colors='gray')
ax2.set_xlim(0, 8)
ax2.set_ylim(0, 1.15)
ax2.grid(True, alpha=0.15)

# --- Panel 3: Vortex lattice in rotating condensate ---
ax3 = axes[2]
ax3.set_facecolor('#0a0a0a')

# Generate triangular lattice positions
Nv = 7  # layers
vortex_x = []
vortex_y = []
a_lat = 2.0  # lattice spacing
for i in range(-Nv, Nv+1):
    for j in range(-Nv, Nv+1):
        px = a_lat * (i + 0.5*j)
        py = a_lat * (np.sqrt(3)/2 * j)
        if px**2 + py**2 < (Nv*a_lat*0.7)**2:
            vortex_x.append(px)
            vortex_y.append(py)

vortex_x = np.array(vortex_x)
vortex_y = np.array(vortex_y)

# Compute condensate density with vortex depletion
Ng = 300
xg = np.linspace(-12, 12, Ng)
yg = np.linspace(-12, 12, Ng)
Xg, Yg = np.meshgrid(xg, yg)

# Thomas-Fermi envelope
R_TF = 11.0
density = np.maximum(1 - (Xg**2 + Yg**2)/R_TF**2, 0)

# Deplete density at each vortex core
xi_v = 0.4
for vx_i, vy_i in zip(vortex_x, vortex_y):
    r2 = (Xg - vx_i)**2 + (Yg - vy_i)**2
    density *= r2 / (r2 + xi_v**2)

im3 = ax3.pcolormesh(Xg, Yg, density, cmap='viridis', shading='auto')
ax3.scatter(vortex_x, vortex_y, c='red', s=8, zorder=5, marker='o')
ax3.set_title('Vortex Lattice in Rotating BEC', color='cyan', fontsize=13, fontweight='bold')
ax3.set_xlabel('x / \u03be', color='gray')
ax3.set_ylabel('y / \u03be', color='gray')
ax3.tick_params(colors='gray')
ax3.set_xlim(-12, 12)
ax3.set_ylim(-12, 12)
ax3.set_aspect('equal')
cb3 = plt.colorbar(im3, ax=ax3, fraction=0.046, pad=0.04)
cb3.set_label('|\u03a8|\u00b2 / n\u2080', color='gray')
cb3.ax.tick_params(colors='gray')

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
print("Vortex velocity field, density profile, and lattice structure plotted.")
print(f"Number of vortices in lattice: {len(vortex_x)}")
print(f"Lattice constant: {a_lat:.1f} xi")
print(f"Healing length: {xi:.1f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Vortex Energy vs System Size and Healing Length

Fortran

program.f9087 lines

program vortex_energy
  implicit none

! Compute vortex line energy E_v = pi*n_s*(hbar/m)^2*L*n^2*ln(R/xi)
  ! Units: pi*n_s*(hbar/m)^2*L = 1

integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: N_R = 20
  integer, parameter :: N_xi = 5

real(dp) :: R_min, R_max, R, dR_log
  real(dp) :: xi_values(N_xi), energy
  real(dp), parameter :: pi = 3.14159265358979323846_dp
  integer :: i, j, k, n_wind

xi_values = (/0.1_dp, 0.5_dp, 1.0_dp, 2.0_dp, 5.0_dp/)
  R_min = 2.0_dp
  R_max = 1000.0_dp
  dR_log = (log(R_max) - log(R_min)) / real(N_R - 1, dp)

! Table 1: Energy vs R for fixed xi = 1.0, winding numbers n=1,2,3
  write(*,'(A)') "============================================="
  write(*,'(A)') " VORTEX LINE ENERGY vs SYSTEM SIZE"
  write(*,'(A)') " (xi=1.0, E in units of pi*n_s*(hbar/m)^2*L)"
  write(*,'(A)') "============================================="
  write(*,'(A12, 3A14)') "R/xi", "n=1", "n=2", "n=3"
  write(*,'(A)') "---------------------------------------------"

do i = 1, N_R
    R = exp(log(R_min) + real(i-1, dp) * dR_log)
    write(*,'(F12.2)', advance='no') R
    do k = 1, 3
      energy = real(k**2, dp) * log(R)
      write(*,'(F14.6)', advance='no') energy
    end do
    write(*,*)
  end do

! Table 2: Energy vs healing length for R = 100
  write(*,*)
  write(*,'(A)') "============================================="
  write(*,'(A)') " VORTEX ENERGY vs HEALING LENGTH (R=100, n=1)"
  write(*,'(A)') "============================================="
  write(*,'(A12, A14, A14)') "xi", "R/xi", "E_v"
  write(*,'(A)') "---------------------------------------------"

R = 100.0_dp
  do j = 1, N_xi
    energy = log(R / xi_values(j))
    write(*,'(F12.4, F14.2, F14.6)') xi_values(j), R/xi_values(j), energy
  end do

! Table 3: Vortex-antivortex pair energy
  write(*,*)
  write(*,'(A)') "============================================="
  write(*,'(A)') " VORTEX-ANTIVORTEX PAIR ENERGY"
  write(*,'(A)') " E_pair = 2*ln(d/xi)"
  write(*,'(A)') "============================================="
  write(*,'(A12, A14)') "d/xi", "E_pair"
  write(*,'(A)') "---------------------------------------------"

do i = 1, 15
    R = 1.0_dp + real(i, dp) * 3.0_dp
    write(*,'(F12.2, F14.6)') R, 2.0_dp * log(R)
  end do

! Table 4: Feynman relation
  write(*,*)
  write(*,'(A)') "============================================="
  write(*,'(A)') " FEYNMAN RELATION: n_v = 2*Omega/kappa"
  write(*,'(A)') "============================================="
  write(*,'(A12, A14, A14)') "Omega", "n_v*R^2", "l/R"
  write(*,'(A)') "---------------------------------------------"

do i = 1, 10
    R = 0.2_dp * real(i, dp)
    energy = 2.0_dp * R / (2.0_dp * pi)
    write(*,'(F12.4, F14.6, F14.6)') R, energy, &
      1.0_dp / sqrt(max(energy, 1.0e-10_dp))
  end do

write(*,*)
  write(*,'(A)') "Key result: E_vortex ~ n^2 * ln(R/xi)"
  write(*,'(A)') "Multiply-charged vortices (n>1) unstable"
  write(*,'(A)') "against decay into |n| single vortices."
end program vortex_energy

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Chapter 4: Gross-Pitaevskii Part III: Cooper Pairs →