Condensate: an interactive simulator of two-dimensional quantum fluids

Introduction

Quantum fluids (and fluids in general) are best understood by watching them move, even better poking and prodding them to get a feel for how they respond. Indeed, phenomena such as vortex nucleation, the lattice of a rotating condensate, the unbinding of vortex pairs at the Berezinskii–Kosterlitz–Thouless (BKT) transition and the break-up of a dipolar gas into droplets are all dynamical and crucially can be captured in a two-dimensional model. Static figures convey this poorly, and the analysis that explains them is often beyond what students can do unaided. In a number of fields, interactive simulations have bridged a similar gap. However, typically research codes in the cold-atom community require installation, programming and offline visualisation, which puts them out of reach of new students.

This paper describes Condensate, an interactive simulator designed to let students produce and start to grasp these phenomena themselves, in real time, on a device they already carry. Its design follows a growing family of browser-based teaching tools: VisualPDE [2] solves a broad class of user-defined partial differential equations on the GPU in real time; its authors argue that interactive “play” builds intuition that pen-and-paper analysis alone does not, especially for behaviour that is analytically intractable at undergraduate level. Earlier examples include the PhET simulations for physics teaching [3], Complexity Explorables [4], the WebGL simulations of cardiac and nonlinear dynamics of Kaboudian, Cherry and Fenton [5], and, for quantum mechanics, the QuVis visualisations [6]. These tools share three features that Condensate adopts: no installation, immediate response to the user’s input, and a curated set of models with sensible parameter ranges.

For quantum fluids, existing tools fall into two groups. Research codes such as GPELab [7], XMDS2 [8], GPUE [9] are powerful and well validated, but they are tools for research rather than teaching. Interactive tools exist but are narrower. The WebGL solver [10] of Stagg simulates the single-component two-dimensional GPE in the browser, and the super_hydro explorer of Forbes and co-workers [11] offers real-time interaction with superfluids through a Python and Jupyter environment that must be installed locally. Condensate sits between these groups. It provides:

• eight quantum-fluid models in one interface, sharing one numerical method and one set of units, so that students can compare, for example, a contact-interacting gas with a dipolar or self-gravitating one;

• a WebGPU solver fast enough for a phone to run a $256^2$ grid interactively, packaged as one HTML file with no dependencies that also works offline;

• interaction as part of the physics: the user’s finger is a repulsive optical potential, and a tap imprints a quantised vortex;

• teaching material for each model, with activities whose quantitative predictions have been checked against the simulation, and references keyed to the textbook of Barenghi, Bland and Parker [1].

The intended users are lecturers and students in advanced undergraduate and postgraduate courses on superfluidity, cold atoms and nonlinear waves, and audiences for outreach. Researchers may also find it useful for building intuition about a model before writing a production code.

Using Condensate

Figure 1 shows the interface. The simulation fills a square panel. Brightness shows the density $n = |\psi|^2$ and an optional view shows the phase; for two-component models, the two species are shown in red and teal, and the phase view shows their relative phase. Touching the panel switches on a repulsive Gaussian beam at the touch point; dragging moves it, and two fingers make two beams. In “Add vortices” mode a tap instead imprints a quantised vortex of chosen sign. A “Mark vortices” button locates vortices from the phase winding around each grid cell and reports their number and signs.

The controls choose the model, the trap (harmonic with adjustable aspect ratio, box, ring, square optical lattice or none) and each model’s main parameter, together with model-specific controls such as the dipole tilt, the spin–orbit coupling strength or the polariton noise and disorder. In most cases, parameters take effect immediately, so the fluid responds while it evolves. On browsers with WebGPU the grid size can be chosen from $128^2$ to $1024^2$; it defaults to $256^2$ on phones and $512^2$ on computers.

Below the simulation, a panel denoted “Learn” describes the current model (§ 5). It updates when the model changes, and its “Set this up” buttons configure the simulation for each suggested activity.

Physical models

Common framework

All models are written in trap units, $\hbar= m = \omega= 1$, with lengths in $a_{\mathrm{ho}} = \sqrt{\hbar/m\omega}$ and times in $1/\omega$ [1 §3.8]. The fluid is quasi-two-dimensional: the wavefunction is assumed frozen in a Gaussian profile of width $l_z = 1.5$ along $z$, which sets the 2D contact coupling $g$ and the averaging of the long-range kernels below. Each model evolves a generalised, optionally damped, GPE,

where $\gamma\ge0$ is a phenomenological damping [12], $V_{\mathrm{beam}} = V_0 e^{-\lvert\mathbf{r}-\mathbf{r}_b(t)\rvert^2/2\sigma^2}$ is the user’s beam, with $V_0 = 60$ well above the chemical potential, and $\Phi[\psi]$ collects the model-specific terms. At the default $128^2$ grid, $g = 1$ and $N = 475$, so the Thomas–Fermi chemical potential in the harmonic trap is $\mu= \sqrt{gN/\pi} \approx12.3$, the healing length is $\xi= 1/\sqrt{g n_0} \approx0.28$ and the central sound speed is $c \approx3.6$. Table 1 summarises the eight models, and Figure 2 shows a characteristic state of each.

Model	Additional physics	Controls	Phenomena
Standard	none	damping $\gamma$	critical velocity, vortex shedding, vortex dipoles, sound
Rotating	$-\Omega L_z$ (rotating frame)	$\Omega/\omega\in[0,0.85]$	vortex entry, Abrikosov lattice, Feynman’s rule
Dipolar	quasi-2D dipolar kernel with tilt $\alpha$; Lee–Huang–Yang term	$\varepsilon_{dd} \in[0,2.4]$, tilt $\alpha$, azimuth $\varphi$	roton, droplet crystals, magnetostriction, magic angle
Two-component	cross-interaction $g_{12}$; mixture Lee–Huang–Yang term	$g_{12}/g \in[-1.5,1.6]$, species share	miscibility, domains, filled cores, quantum droplets
Rabi / spin–orbit	$\frac{1}{2}(p_x-k_{so}\sigma_z)^2+\frac{1}{2}p_y^2+\frac{\Omega}{2}\sigma_x$	$\Omega$, $k_{so}$, $g_{12}/g$	stripe, plane-wave and zero-momentum phases; confined vortices
Self-gravitating	Schrödinger–Poisson, quasi-2D	$G \in[0,0.4]$	Jeans instability, solitonic cores
Polariton	driven-dissipative GPE with reservoir; noise; disorder	$P/P_{th}$, noise, disorder	threshold, outflow, spontaneous and pinned vortices
Thermal	simple-growth SPGPE at fixed $\mu$	temperature $T$	thermal vortex pairs, BKT physics

Model-specific terms

Dipolar. For dipoles tilted by $\alpha$ from the $z$ axis, with their in-plane projection along $\hat{\mathbf{e}}_{\parallel} = (\cos\phi,\sin\phi)$, averaging the dipole–dipole interaction over the axial profile gives the Fourier-space kernel [16, 17]

with $F = \sqrt{\pi}q e^{q^2}\operatorname{erfc}q$ and $q = k l_z/\sqrt{2}$. At long wavelengths this reduces to $g[1+\varepsilon_{dd}(3\cos^2\alpha- 1)]$, so the isotropic dipolar contribution vanishes at the magic angle $\alpha_m = \arccos(1/\sqrt{3}) \approx54.7^\circ$ [18]. A Lee–Huang–Yang (LHY) term $\gamma_{\mathrm{LHY}}(1+\frac{3}{2}\varepsilon_{dd}^2)|\psi|^3$ with $\gamma_{\mathrm{LHY}} = 0.02$, of the three-dimensional form [19], prevents collapse and stabilises droplets [20, 21].

Two-component. Two species with equal self-interaction $g$ and cross-interaction $g_{12}$ mix for $g_{12} < g$ and separate for $g_{12} > g$ [22]. Each species also feels $\gamma_m(|\psi_1|^2 + |\psi_2|^2)^{3/2}$ with $\gamma_m = 0.04$. For $\delta= -(g + g_{12}) > 0$ the balanced mixture has energy density $-\frac{1}{4}\delta n^2 + \frac{2}{5}\gamma_m n^{5/2}$, and zero pressure gives a self-bound liquid of density $n_0 = (5\delta/12\gamma_m)^2$, Petrov’s quantum droplet [23, 24, 25]. In strictly two dimensions the LHY term becomes logarithmic [26]; the three-dimensional form used here describes a quasi-2D gas whose axial width is large compared with the scattering lengths.

Rabi and spin–orbit coupling. The single-particle Hamiltonian is that of the NIST experiments [27], whose lower band has minima at $\pm k_1 = \pm k_{\mathrm{so}}\sqrt{1-(\Omega/2k_{\mathrm{so}}^2)^2}$ for $\Omega< 2k_{\mathrm{so}}^2$. Interactions select the stripe phase ($g_{12} < g$), the plane-wave phase ($g_{12} > g$) or, above $\Omega= 2k_{\mathrm{so}}^2$, the zero-momentum phase [28, 29]. The spin–orbit term only shifts each species’ kinetic energy to $\frac{1}{2}k^2 \mp k_{\mathrm{so}}k_x$, and the Raman coupling is an exact local rotation, so both remain exact in the split-step scheme. With $k_{so}=0$ the model is a coherently coupled mixture, in which a vortex in one species drags a domain wall of relative phase [31].

Self-gravitating. Newtonian gravity averaged over the axial profile gives

with the uniform background removed as in cosmological codes. This is the Schrödinger–Poisson system used for ultralight dark matter [32, 34], here with an optional contact term.

Polariton. The driven-dissipative GPE of Wouters and Carusotto [35] couples the condensate to an incoherent reservoir,

with threshold $P_{th}=\gamma_c\gamma_R/R$ and, for uniform pumping, steady density $n=(P-P_{th})/\gamma_c$ [36]. The trap buttons shape the pump. Optional truncated-Wigner noise [37] and a Gaussian-correlated disorder potential model a real microcavity.

Rotating. The contact GPE is solved in the frame rotating at $\Omega$, with $-\Omega L_z$ added to the Hamiltonian; with no trap the rotating frame is unbounded, so the trap is required.

Thermal. The simple-growth stochastic projected GPE [38, 39, 40],

with $L=-\frac{1}{2}\nabla^2+V+g|\psi|^2$, couples the low-energy modes $|\mathbf{k}|<k_c=6.5$ to a reservoir at temperature $T$ and chemical potential $\mu=12$. As in all classical-field methods, the temperature at which vortex pairs unbind depends on the cutoff.

Scaling with grid size

The grid spacing is fixed at $17/128$ in trap units, so a larger grid simulates a larger system with the same local physics rather than a finer one. On a grid $s = N/128$ times wider, the harmonic frequency becomes $1/s$, the box and ring radii and the polariton spot grow by $s$, and the atom number grows as $s^2$, which keeps the density, chemical potential and healing length fixed. The harmonic trap can be made anisotropic, $V = \frac{1}{2}(\omega_x^2x^2 + \omega_y^2y^2)$ with $\lambda= \omega_y/\omega_x \in[1, 6]$; since the two-dimensional Thomas–Fermi atom number is $N = \pi\mu^2/(g\omega_x\omega_y)$, the atom number is divided by $\lambda$, which again keeps $\mu$ and $\xi$ fixed and produces cigar-shaped clouds of fixed length. Local features keep their physical size: the box wall, the lattice period, the disorder correlation length and the ring’s cross-section. The ring therefore grows only in circumference, and its atom number scales as $s$. Rotation and gravity are scaled by $1/s$ for the same reason.

Trapping potentials

The trap buttons choose the external potential $V(\mathbf{r})$ in Eq. (1), written here for a grid $s = N/128$ times the default width, with $r = \lvert\mathbf{r}\rvert$. In every case the potential is capped at $V_{\max} = 100$. The cap only affects regions far outside the fluid, where the density is negligible, and it bounds the phase that $V$ adds in a single time step.

Harmonic. A harmonic trap, optionally anisotropic,

where the aspect slider sets $\lambda= \omega_y/\omega_x \in[1, 6]$. At $\lambda= 1$ this is the isotropic trap that defines the units, $\frac{1}{2}r^2$ on the default grid. The atom number is divided by $\lambda$ (§ 3), so that the Thomas–Fermi radii $R_i = \sqrt{2\mu/\omega_i}$ change while $\mu$ does not.

Box. A flat-bottomed disc with a smooth wall,

modelled on the uniform optical box traps of recent experiments [13]. In the box the density is $n \approx N/\pi R_{\mathrm{b}}^{2}$, so $\mu\approx5$ on the default grid; the wall is about eight times higher than $\mu$ and narrower than the local healing length, $\xi\approx0.45$, so it is effectively hard, while the tanh profile avoids the ringing that a step would produce in a spectral method.

Ring. An annular trap, harmonic in the radial direction,

the geometry of persistent-current experiments [14]. The radial frequency, and so the ring’s cross-section, is fixed, so only the circumference grows on larger grids, and the atom number grows as $s$ rather than $s^{2}$.

Lattice. A square optical lattice on a harmonic envelope,

with lattice period $a$, the recoil energy $E_{\mathrm{R}}=\hbar^{2}k_{\mathrm{L}}^{2}/2m$ of the lattice wavenumber $k_{\mathrm{L}}=\pi/a$, and a depth $V_{0}^{\mathrm{L}}\in[0,10]$ in recoil energies set by a slider [1 §3.9.3], [15]. On the default grid the period equals the oscillator length, and on every grid an integer number of periods ($17s$) fits the box, so the lattice is consistent with the periodic boundary conditions. With $E_{\mathrm{R}}=\pi^{2}/2\approx4.9$, the default depth of $3E_{\mathrm{R}}$ is comparable to $\mu$, between the nearly free and tight-binding regimes.

None. $V=0$. The periodic box then approximates an infinite uniform gas, which is the natural setting for droplets, the Jeans instability and thermal vortex pairs. The rotating model requires a trap, since the rotating frame is unbounded without one.

Polariton pump profiles. In the polariton model the same buttons also shape the pump, $P(\mathbf{r})=pP_{\mathrm{th}}f(\mathbf{r})$, where $p$ is the pump slider and

The condensate forms where $P(\mathbf{r})>P_{\mathrm{th}}$ and, in the local-density approximation, has density $n(\mathbf{r})\approx(\gamma_{\mathrm{R}}/R)[P(\mathbf{r})/P_{\mathrm{th}}-1]$, so the pumped region appears as the bright part of the fluid, and the condensate ends roughly where the pump falls below threshold. The trap potential still acts alongside the pump, and the disorder potential (§ 4) is added to $V$ when its slider is non-zero.

Numerical method

Discretisation

The wavefunction lives on an $N\times N$ periodic grid covering a square of side $L=17s$, with $s=N/128$, so the grid spacing is fixed at $\Delta x=17/128\approx0.133$ (§ 3). Spatial derivatives are spectral: the discrete Fourier transform $\mathcal{F}$ represents $\psi$ by the wavenumbers $k_{x,y}=2\pi m/L$ with $-N/2 \le m < N/2$, up to the Nyquist wavenumber $k_{\max}=\pi/\Delta x \approx23.7$. The healing length at the default density, $\xi\approx0.28$, spans about two grid spacings, and a vortex core is resolved by several points. The wavefunction is normalised to the atom number, $\sum_j |\psi_j|^2\Delta x^2=N_{\mathrm{at}}$.

Split-step integrator

Every model is advanced by Strang splitting [42, 41] of the generator into a local part, diagonal in real space, and a kinetic part, diagonal in Fourier space. Specifically, we integrate

which is Eq. (1) with time rescaled by $1+\gamma^2$, a negligible difference for $\gamma\le0.1$. With $H=T+W$, where $T=-\frac{1}{2}\nabla^2$ and $W$ collects every local term, one step of length $\Delta t$ is

with

The local potential is

where $W_{\mathrm{LHY}}$ is the relevant Lee–Huang–Yang term and $W_{\mathrm{model}}$ holds the cross-interaction, reservoir or disorder terms of the multi-component and polariton models. Each half-step recomputes $|\psi|^2$ before applying $\mathcal{L}$, so the contact nonlinearity is treated to second order in $\Delta t$; for $\gamma=0$ it is exact, because a local phase rotation leaves $|\psi|^2$ unchanged. Each factor is unitary when $\gamma=0$, so the scheme conserves the norm and has no stability limit from the kinetic or contact terms. The real-time step is $\Delta t=0.008$, except for tilted dipoles (§ 4.4).

Chemical potential and normalisation. In closed systems the wavefunction is rescaled after every step to restore the atom number, which makes Eq. (11) the usual damped GPE at fixed $N_{\mathrm{at}}$. The chemical potential is estimated as the density-weighted mean of the local potential, $\mu\approx\sum_j n_jW_j/\sum_j n_j$, which keeps the phase rotation of the condensate slow; it is updated every tenth step on the CPU and every step on the GPU, from the previous step’s reductions. Because the norm is restored anyway, $\mu$ only sets the reference energy for the damping.

Spectral filter. The kinetic factor includes the filter $f(k)=\exp[-1.5(k/k_{\max})^{12}]$. It is above $0.99$ for $k<0.66\,k_{\max}$, falls to $0.22$ at $k_{\max}$, and removes the corners of the Fourier square, where $|\mathbf{k}|$ reaches $\sqrt{2}\,k_{\max}$. Applied once per step, it acts as a weak hyperviscosity that controls the build-up of aliased power from the nonlinear terms, which generate wavenumbers beyond those the grid can represent.

Nonlocal interactions

The dipolar and gravitational potentials are convolutions of the density with a kernel, and are evaluated with the convolution theorem,

at a cost of two extra two-dimensional FFTs per step. The kernels $\widetilde{U}(\mathbf{k})$ of § 3 are tabulated on the grid when a parameter changes. The function $\operatorname{erfcx} q = e^{q^2}\operatorname{erfc} q$ in both kernels is evaluated with a rational approximation for $q < 2$ and a continued fraction for $q \ge2$. For the dipolar kernel the $k = 0$ mode, $g\varepsilon_{dd}(3\cos^2\alpha- 1)$, is retained (the mean-field shift of a uniform dipolar gas); for self-gravity the $k = 0$ mode is set to zero, which removes the interaction with the uniform background, as is standard in cosmological codes.

The transform is periodic, so each cloud also interacts with its periodic images. The dipolar interaction falls off as $1/r^3$ and the clouds fill only part of the box, so the image interaction is small; for gravity the leading error is a tidal term of relative order $(R/L)^3$ for a cloud of radius $R$. Whilst zero padding would remove these errors (at four times the computational cost), we judged it unnecessary for a teaching tool.

$\Phi$ is computed once per step, from the density at the start of the step, and held fixed through both half-steps. This saves two FFTs per step but makes the nonlocal coupling first order in time.

Time step for tilted dipoles

Tilting the dipoles makes the kernel strongly repulsive at short wavelengths along their in-plane direction, where $\widetilde{U} \to g + 2g\varepsilon_{dd}$, whereas at $\alpha= 0$ it tends to $g - g\varepsilon_{dd}$. Because $\Phi$ is treated explicitly, the scheme is then only conditionally stable, so the dipolar step is reduced adaptively to

which leaves $\alpha= 0$ unchanged and reaches 0.0027 at $\varepsilon_{dd} = 2.4$, $\alpha= 90^\circ$. The same step is stable for every azimuth $\varphi$, including the grid diagonal, where the largest wavenumbers are $\sqrt{2}$ times those along an axis. Strongly tilted, strongly dipolar runs therefore advance up to three times more slowly.

Rotation

In the rotating frame the Hamiltonian gains $-\Omega L_z = \Omega y p_x - \Omega x p_y$. Combined with the kinetic energy, this term is diagonal in neither real nor Fourier space, so the kinetic step is replaced by an alternating-direction split [43],

in which $A$ is diagonal in $(k_x,y)$ and $B$ in $(x,k_y)$. The rotating kinetic step is then

nested inside the outer Strang split. Each factor is applied with one-dimensional transforms: $e^{-(i+\gamma)Ah}$ by a forward FFT along $x$ for every row of fixed $y$, multiplication by $\exp[-(i+\gamma)(\frac{1}{2}k_x^2+\Omega yk_x)h]$ and an inverse FFT along $x$; $B$ likewise along $y$ for every column of fixed $x$. The three factors cost six sets of one-dimensional transforms, the equivalent of one and a half two-dimensional FFT pairs, and each is unitary when $\gamma= 0$, so the scheme remains norm-preserving and unconditionally stable. The spectral filter is applied in each direction separately. Because $[A,B] \ne0$ the inner split has its own second-order error, of the same order as the outer one.

Two-component and coupled models

In the two-component model each species has its own kinetic step and local potential, the latter including $g_{12}|\psi_{3-j}|^2$ and the mixture LHY term $\gamma_m(n_1+n_2)^{3/2}$; both local potentials are evaluated from the densities before either species is updated. Each species is normalised separately.

In the coupled model the spin–orbit term is diagonal in Fourier space and only shifts the kinetic energy of the two species to $T_{\pm}=\frac{1}{2}k^{2}\mp k_{\mathrm{so}}k_{x}$, which the kinetic factor includes exactly. The Raman coupling $\frac{1}{2}\Omega\sigma_{x}$ is a local $2 \times2$ operator whose exponential is known in closed form,

which reduces to a rotation, $\cos(\frac{1}{2}\Omega h)1-i\sin(\frac{1}{2}\Omega h)\sigma_x$, when $\gamma=0$. It is applied after the local factor in the first half-step and before it in the second, so the full step remains second order; the coupling moves atoms between species, only the total atom number is restored, and the two species share one chemical potential.

Ground states

Initial states are found by imaginary-time evolution [44]: the same split-step factors with $i+\gamma$ replaced by 1 and $\mu=0$, and the norm restored after every step, a normalised gradient flow. The imaginary step is $\Delta\tau=0.004$ ($0.01$ for the coupled model), for 500 steps (900 rotating, 600 coupled). The initial guess is the Thomas–Fermi profile $n=\max(0,\mu-V)/g_{\mathrm{eff}}$, with $\mu$ found by bisection and an effective coupling that includes the long-wavelength interaction: $g[1+\varepsilon_{\mathrm{dd}}(3\cos^{2}\alpha-1)]$ for dipoles and $\frac{1}{2}g(1+g_{12}/g)$ for balanced mixtures. A random perturbation of up to $\pm3\%$ breaks symmetries. Three models need a more specific seed. The rotating model imprints the Feynman number of vortices, $\Omega R^{2}$ for a cloud of radius $R$, at random positions, so that cooling relaxes them into a lattice instead of waiting for them to enter through the surface. The coupled model starts from lower-band spinors at $\pm k_{1}$ with unequal weights, so that cooling can select a stripe or plane-wave phase. A self-bound droplet starts as a flat-topped disc at the predicted density $n_{0}$. After cooling, a small amount of noise is added so that real-time instabilities can develop.

Open and stochastic models

Polariton. The gain and loss enter the local factor as the real exponential $\exp[\frac{1}{2}(Rn_{R}-\gamma_{c})h]$, alongside the potential. The reservoir obeys a linear equation once $|\psi|^{2}$ is frozen over a step, and is advanced exactly,

an exponential integrator that is stable for any step. The kinetic factor carries a weak damping, $\gamma=0.01$, which models energy relaxation. The wavefunction is not renormalised: its amplitude is set by the balance of gain and loss. The condensate grows from a smooth random seed, built from modes with $|\mathbf{k}|<3$, so that it forms in real time.

Stochastic terms. The thermal model and the optional polariton noise add complex Gaussian white noise once per step, in Fourier space after the kinetic factor, to the modes with $|\mathbf{k}|<k_{c}=6.5$. The noise has real-space variance $2D\Delta t/\Delta x^{2}$ per grid point, with $D=\gamma T$ for the thermal model, matching $\langle dW^{*}dW\rangle=2\gamma T\,dt$, and $D=\sigma\gamma_{c}/2$ for polaritons. This is an Euler–Maruyama treatment of the noise within the split step. In the thermal model the projector $\mathcal{P}$ replaces the spectral filter: the kinetic factor is multiplied by the indicator of $|\mathbf{k}|<k_{c}$, so the noise and the nonlinear term never populate modes above the cutoff. The damping $\gamma=0.05$ appears in both the local and kinetic factors, through $i+\gamma$, and $\mu$ is fixed at 12, so the atom number is set by the reservoir rather than imposed. Both noisy models display a running average of the density, $\bar n\leftarrow\bar n+(\Delta t/\tau)(|\psi|^{2}-\bar n)$ with $\tau=0.25$, which plays the role of a camera’s exposure time.

Disorder. The polariton disorder potential is a Gaussian random field with correlation

It is drawn in Fourier space with amplitude $e^{-k^{2}\ell^{2}/8}$ and random phase, transformed back, and rescaled to zero mean and unit variance; the slider sets its rms amplitude.

Interaction

The beam is the potential $V_{\mathrm{beam}} = V_{0}\sum_{b}a_{b}(t)\exp[-|\mathbf{r}-\mathbf{r}_{b}(t)|^{2}/2\sigma^{2}]$, which enters $W$. Each beam centre $\mathbf{r}_{b}$ follows the touch point at a speed of at most 4 trap units, which is faster than the central sound speed, and its amplitude $a_{b}$ ramps between 0 and 1 over 0.2 time units, so that switching a beam on or off does not deliver an impulsive kick. A tap in vortex mode multiplies the wavefunction by $r/\sqrt{r^{2}+2\xi^{2}}\,e^{\pm i\theta}$, centred on the tap, with $\xi$ from the local density, which imprints a vortex with approximately the right core profile. Vortices are located by summing the phase differences around each grid plaquette and keeping windings of $\pm2\pi$ where the smoothed density exceeds a fifth of its maximum, which rejects phase defects in the low-density edge.

Implementation

GPU engine. The GPU engine implements the whole loop in WGSL compute shaders, in single precision with complex numbers stored as pairs of 32-bit floats. We note that the built-in WGSL sine and cosine are only guaranteed to about $5 \times10^{-4}$. Hence, the GPU engine uses its own hand-coded sine and cosine, with argument reduction by $\pi/2$ in two parts followed by minimax polynomials, accurate to about $10^{-7}$. This ensured good agreement with the CPU approach.

CPU engine. The CPU engine performs the same operations in double precision at $128^{2}$, with an in-place radix-2 FFT applied along rows and, after a transpose, along columns. It caches the kinetic and rotation factors and recomputes them only when a parameter changes. It runs wherever WebGPU is unavailable and served as the reference for testing the GPU engine.

Architecture. The application is written in plain JavaScript and assembled from three sources into a single HTML file: the interface, the CPU engine and the GPU engine. Both engines expose the same interface, so the controls drive either. The Learn panels are written in LaTeX notation and converted to MathML at build time, so equations render natively in current browsers without external libraries. The page needs no network access once loaded.

Teaching with Condensate

The Learn panels. Each model has a panel, aimed at graduate students, with four parts: the model, giving the equations the application solves; the physics, giving the key arguments and results; try this, three or four activities; and read more. A collapsible Foundations section, common to all models, introduces the GPE, the hydrodynamic (Madelung) interpretation of density and phase, the Thomas–Fermi profile, the healing length and speed of sound with their values in the app, quantised vortices, and the numerical method.

Activities. Each activity pairs a prediction that students can make with pen and paper with an experiment in the app. Examples include estimating the lattice spacing from Feynman’s rule and checking it with the vortex markers; predicting the stripe period $\pi/k_{1} \approx1.8$ of the spin–orbit-coupled gas and the coupling $\Omega= 2k_{\mathrm{so}}^{2}$ at which the stripes vanish; predicting the miscibility threshold; computing the flat-top density of a quantum droplet; and deriving the Jeans threshold from the dispersion relation. Every quantitative prediction quoted in the panels was checked against the simulation. Some activities are deliberately open-ended. A “Set this up” button beside each activity configures the model, trap and parameters on the $128^2$ grid, so that the quoted numbers apply.

Links to the literature. For convenience, the panels cite sections of the standard reference Quantum Fluids, Solitons, and Vortices [1], whose chapters on the GPE, waves, superfluidity, vortices, two-component and dipolar condensates, and droplets and supersolids map directly onto the models. Where the book does not cover a model, the panels cite review and primary literature instead: for spin–orbit coupling [27, 30], ultralight dark matter [32, 33], polaritons [36] and BKT physics [45, 46, 47, 48].

Limitations and Outlook. Condensate is a teaching tool, and its design makes compromises because of this. The models are fixed rather than user-defined and two-dimensional and the GPU engine runs in single precision. We would be interested in hearing your comments or criticisms, additions to the current models or further extensions that would prove useful to the community.

Availability

Condensate is freely available from https://abag.github.io/condensate/ and runs in most modern browsers; the GPU engine needs WebGPU, available in Chrome and Edge 113 or later, Safari 26 or later on macOS, iOS and iPadOS, and Firefox on Windows.

Acknowledgements

The software the paper describes was developed with the assistance of Claude (Anthropic), a large language model, under the author’s direction. We specified the physics, reviewed the code and derivations, directed validation tests and take full responsibility for the software and its description. We acknowledge fruitful conversations with Nick Parker and Ryan Doran, which improved the depth and functionality of the software.