Extended Gross–Pitaevskii Framework

• Extended Gross–Pitaevskii framework is a set of generalized equations that relax standard assumptions to model finite interaction ranges, nonlocal effects, and quantum fluctuations.
• It incorporates modifications like derivative corrections, LHY terms, and pairing-field dynamics to capture beyond-mean-field behaviors in Bose systems.
• The framework underpins rigorous derivations and numerical methods, enabling precise modeling of condensate dynamics in systems ranging from dipolar droplets to self-gravitating superfluids.

The expression extended Gross–Pitaevskii framework denotes a family of non-equivalent generalizations of the Gross–Pitaevskii equation constructed to relax one or more of the standard GP assumptions: a dilute Bose gas at $T\approx 0$, dominance of $s$-wave scattering, zero-range contact interactions modeled by a $\delta$-pseudopotential, and a mean-field Hartree description in which all atoms occupy the same orbital $\psi$. Across the literature, the same label is used for derivative corrections that encode finite interaction range, Lee–Huang–Yang-augmented equations for droplets, rigorous condensate-plus-Bogoliubov descriptions in the Gross–Pitaevskii regime, matrix-valued quantum equations on the continuous matrix product state manifold, nonlocal self-interaction-subtracted theories for excitonic solids, and Gross–Pitaevskii–Poisson systems for self-gravitating superfluids (Veksler et al., 2014, Haegeman et al., 2015, Caraci et al., 2023, Conti et al., 27 Jul 2025, Shukla et al., 12 Dec 2025, Huang et al., 6 Apr 2026).

1. Baseline Gross–Pitaevskii structure and the scope of extension

In its canonical form, the time-dependent Gross–Pitaevskii equation is

$$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$

with energy functional

$$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$

The order parameter may be normalized either to the particle number $N$ or to $1$, with the factor $N$ absorbed into the nonlinear coefficient. This baseline is the reference point from which all “extended” GP constructions depart (Veksler et al., 2014).

The extensions in the literature are not a single hierarchy. Some retain a scalar order parameter but modify the interaction term; some keep the mean-field PDE but change its equation of state; some replace the scalar wavefunction by matrix-valued fields; and some treat the GP equation only as the condensate sector of a larger condensate-plus-fluctuation dynamics. A compact classification is therefore useful.

Extension class	Defining modification	Representative sources
Finite-range or effective-range GP	Adds $g_2\nabla^2(|\psi|^2)\psi$ or equivalent short-range nonlocal correction	(Veksler et al., 2014)
LHY-augmented eGPE	Adds a local beyond-mean-field term proportional to $s$0 or $s$1	(Stringari, 2018, Huang et al., 6 Apr 2026)
Rigorous GP-plus-Bogoliubov dynamics	Keeps GP evolution for the condensate and adds microscopic correlation dressing and quadratic fluctuation dynamics	(Caraci et al., 2023, Brooks, 2023)
Multi-component hard-core GP derivation	Derives coupled GP equations with uniform-in-time convergence bounds	(Chong et al., 30 Jan 2025)
Pairing-field droplet theory	Interprets the eGPE as the equation of motion for a pairing field rather than the condensate	(Hu et al., 2020)
Quantum/cMPS GP	Replaces the scalar field by matrix-valued cMPS variables and derives a non-commutative QGPE	(Haegeman et al., 2015)
Nonlocal correlated and self-interaction-subtracted GP	Uses a nonlocal kernel and explicitly excludes intra-cell self-interactions	(Conti et al., 27 Jul 2025)
Self-gravitating GPPE	Couples GP dynamics to Poisson gravity, and in neutron-star applications to Maxwell fields and multi-component order parameters	(Shukla et al., 12 Dec 2025)

This plurality is the central feature of the subject. A statement about an “extended GP equation” is therefore only meaningful together with the mechanism of extension.

2. Finite-range, nonlocal, and long-range interaction extensions

One of the most direct generalizations of the GP equation replaces the zero-range contact approximation by an effective-range correction. A simple model uses a repulsive core plus an attractive thin shell, which in the zero-thickness limit generates a sum of $s$2-functions and captures the essential structure of van der Waals interactions while remaining analytically tractable. In the resulting stationary equation, the leading finite-range correction is a derivative term independent of the microscopic details of the interaction potential,

$s$3

with

$s$4

In an elongated trap, the same physics can be recast in a simple 1D three-$s$5 model,

$s$6

so that effective-range information is encoded in a single length scale $s$7 (Veksler et al., 2014).

In the Thomas–Fermi regime of a 1D harmonic trap, this extension yields closed-form leading-order corrections. The density acquires a constant offset, the chemical potential shifts by

$s$8

and the per-particle energy gains

$s$9

The relative correction is

$\delta$0

A distinctive feature is that $\delta$1, whereas the standard GP Thomas–Fermi energy scales as $\delta$2 in the 1D setting analyzed there. This scaling separation is the proposed experimental handle for finite-range physics (Veksler et al., 2014).

The same derivative viewpoint fits into a broader nonlocal perspective. For short-range nonlocal kernels, the leading effect can be written as $\delta$3, whereas for genuinely long-range kernels one keeps an explicit convolution. A different long-range extension appears in self-gravitating superfluids, where the Gross–Pitaevskii–Poisson equation takes the form

$\delta$4

This framework is used for boson and axion stars, dark-matter haloes, and minimal neutron-star models; for axions it is extended further to a cubic–quintic GPPE, and in neutron-star applications it is coupled to a proton superconductor, the Maxwell equations for the vector potential $\delta$5, and Poisson gravity (Shukla et al., 12 Dec 2025).

These interaction-based extensions remain closest to the original GP philosophy: a macroscopic order parameter is preserved, but the contact Hartree term is replaced by a more structured interaction functional.

3. Quantum-fluctuation corrections and droplet eGPEs

A second major branch of the extended GP literature incorporates beyond-mean-field quantum fluctuations through a local Lee–Huang–Yang correction. In uniform matter, linearized time-dependent GP theory reproduces the Bogoliubov dispersion,

$\delta$6

the static structure factor $\delta$7, the momentum distribution

$\delta$8

and the quantum depletion fraction

$\delta$9

Within the same framework, the chemical potential acquires the LHY form

$\psi$0

which corresponds to an additional local term proportional to $\psi$1 in the extended GP equation (Stringari, 2018).

In dimensionless notation, a widely used $\psi$2-dimensional eGPE with LHY correction is

$\psi$3

with energy

$\psi$4

For $\psi$5 and free space, the ground-state problem exhibits dimension-dependent existence properties: if $\psi$6, $\psi$7 and there is no minimizer; if $\psi$8, minimizers exist for all $\psi$9 in $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$0, while in $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$1 there is a threshold $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$2 below which no minimizer exists and above which minimizers exist. Under a confining potential $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$3 with $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$4 as $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$5, a minimizer exists for every $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$6 and $$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$7 (Huang et al., 6 Apr 2026).

This equation supports distinct free-space regimes that the paper organizes as no-ground-state, soliton-like, and droplet-like. In the droplet regime, a simple flat-top approximation neglects kinetic energy and yields

$$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$8

which captures the high-density core of self-bound droplets (Huang et al., 6 Apr 2026).

Dipolar quantum droplets provide a stringent test of this eGPE logic. The eGPE used for dysprosium droplets adds dipolar mean-field interactions and a local LHY chemical potential

$$i\hbar\,\partial_t\psi(\mathbf{r},t) = \left[ -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g\,|\psi(\mathbf{r},t)|^2 \right]\psi(\mathbf{r},t), \qquad g=\frac{4\pi\hbar^2 a_s}{m},$$9

together with a three-body loss term in dynamical simulations. Although the central density of a saturated Dy droplet corresponds to $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$0, the measured critical atom numbers are systematically lower than the eGPE prediction, while quantum Monte Carlo calculations agree well with experiment. This directly questions the validity of the current theoretical state-of-the-art description of dilute dipolar droplets within the eGPE framework and indicates that correlations are significant (Böttcher et al., 2019).

A plausible implication is that “extended GP” in the LHY sense is quantitatively reliable only within a limited correlation window. The formal addition of a local $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$1 energy density is not automatically sufficient when nonlocality, finite-range effects, or strong dipolar correlations become important.

4. Microscopic derivations, rigorous many-body theory, and coupled GP limits

A different meaning of extended GP emerges when the GP equation is derived as only one sector of a more complete many-body description. In the three-dimensional Gross–Pitaevskii regime, one rigorous framework lets the condensate evolve according to a time-dependent GP equation, factors out the microscopic two-body correlation structure through the zero-energy scattering solution, and describes the orthogonal excitations by a quadratic Bogoliubov dynamics. In this setting the many-body dynamics is approximated in norm with rate $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$2, the Bogoliubov dynamics converges to an $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$3-independent limit with rate $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$4, and a central limit theorem for bounded observables is obtained with rate $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$5 (Caraci et al., 2023).

A complementary rigorous development diagonalizes dilute Bose gases algebraically in the GP regime and beyond. In the GP case it recovers the leading-order energy $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$6, yields the discrete Lee–Huang–Yang correction, identifies the excitation dispersion

$$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$7

and provides an upper bound with a logarithmic remainder of the conjectured next-order scale. This same analysis motivates an extended GP functional with an LHY term and an associated Euler–Lagrange equation valid in the dilute regime $$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$8 (Brooks, 2023).

For mixtures, the extension can be made fully quantitative. In the hard-core GP limit of a two-component Bose gas, the effective equations are

$$E[\psi] = \int d^3r\, \left[ \frac{\hbar^2}{2m}|\nabla\psi|^2 + V(\mathbf{r})|\psi|^2 + \frac{g}{2}|\psi|^4 \right].$$9

with a modified $N$0-dependent system incorporating Neumann scattering correctors used in the derivation. The main theorem proves a uniform-in-time trace-norm bound

$N$1

and the same strategy extends to the one-component case and to any finite number of components $N$2 (Chong et al., 30 Jan 2025).

Binary quantum droplets introduce a further microscopic refinement. In a homonuclear Bose–Bose mixture with attractive inter-species interaction, the condensate wavefunction obeys

$N$3

where $N$4 is a bosonic pairing field generated by a Hubbard–Stratonovich decoupling of the inter-species channel. The paper’s central claim is that the “extended Gross–Pitaevskii equation” used earlier for droplets should be understood as the equation of motion for the pairing field, rather than the condensate. In uniform matter this coupled theory yields a gapless branch

$N$5

and a gapped branch

$N$6

so that low-energy collective excitations are fluctuations of the pairing field, while the atomic Bogoliubov equations describe single-particle-like excitations (Hu et al., 2020).

This many-body branch of the subject therefore extends GP in a stricter sense than a modified nonlinear coefficient. It retains the condensate PDE but embeds it within a fluctuation theory, a scattering renormalization, or a coupled condensate–pairing system.

5. Strong-correlation and non-mean-field reformulations

In one dimension, the failure of the ordinary GPE is not merely quantitative. A non-commutative quantum Gross–Pitaevskii equation is obtained by applying the time-dependent variational principle to the variational manifold of continuous matrix product states. The dynamical variables are matrix-valued fields $N$7, $N$8, and $N$9, together with left and right density matrices $1$0 and $1$1. The gauge-covariant QGPE reads

$1$2

supplemented by an equation for $1$3 and by gauge-density equations. When the bond dimension is $1$4, all commutators vanish and the QGPE reduces to the ordinary GPE; for $1$5, it retains entanglement and correlations. Linearization yields quantum Bogoliubov–de Gennes equations whose steady-state response shows a pronounced feature near $1$6, consistent with the Tonks–Girardeau limit and type II excitations (Haegeman et al., 2015).

A still more radical extension appears in the theory of an excitonic incompressible Bose solid in double-layer semiconductor heterostructures. There the effective interaction kernel is nonlocal and correlated,

$1$7

with the hard-core radius $1$8 extracted from QMC structure factors of 2D dipoles. The generalized stationary GP equation is

$1$9

where

$N$0

The second term removes the particle’s own-cell contribution and explicitly excludes self-interactions. This prescription stabilizes an incompressible supersolid with exactly one boson per lattice site and a finite homogeneous background $N$1. Numerically, at fixed low density $N$2, the solutions are uniform superfluid-like at $N$3, become supersolid at $N$4, and the energy eigenvalues cross around $N$5; the stationary problem is solved by imaginary-time propagation on a $N$6 grid with convergence at relative energy change $N$7 (Conti et al., 27 Jul 2025).

This framework also clarifies a common misconception. Cluster supersolids in dipolar gases are generally compressible and arise from a competition between attractive and repulsive forces, whereas the exciton supersolid uses purely repulsive interactions together with explicit self-interaction exclusion to enforce single occupancy and incompressibility (Conti et al., 27 Jul 2025).

6. Solvable constructions, perturbative schemes, computation, and conceptual ambiguities

Some extended GP frameworks are analytic rather than physical extensions. One construction defines a family of exactly solvable stationary GP-type equations with variable coefficients and arbitrary nonlinear powers by the transform

$N$8

together with explicit coefficient maps

$N$9

The resulting “family” transfers exact solvability from one member to all others and covers the nonlinear Schrödinger equation, quintic GP, cubic–quintic GP, generalized polynomial nonlinearities, and a 3D spherically symmetric radial GP-type equation (Liu et al., 2020).

Another analytic route is perturbative. For the time-independent GP equation in $g_2\nabla^2(|\psi|^2)\psi$0 with general external potential, a Green’s-function/eigenbasis expansion yields first-order corrections to the wavefunction, density, chemical potential, internal energy, and mean-square radius. The true small parameter is $g_2\nabla^2(|\psi|^2)\psi$1, and the validity criterion is

$g_2\nabla^2(|\psi|^2)\psi$2

Applied to harmonic traps in $g_2\nabla^2(|\psi|^2)\psi$3, the method gives explicit first-order formulas and predicts breakdown once the interaction-driven admixture of excited states is no longer small; the paper notes, for example, breakdown around $g_2\nabla^2(|\psi|^2)\psi$4 in the 2D and 3D examples discussed there (Abulseoud et al., 2023).

The numerical literature extends the framework in a different direction. For steady symmetric and central vortex states of nonlocal GP equations coupled to microwaves, two accelerated solvers are proposed: the accelerated stabilized normalized gradient flow (ASGF) method and the perturbed, projected nonlinear conjugate gradient (PPNCG) method. The governing equations are discretized by a Legendre–Galerkin spectral method in polar coordinates, and the PPNCG method performs the best in the reported benchmarks, consistently reducing iteration counts and CPU times relative to ASGF and earlier normalized-gradient methods (Wang et al., 2022).

Two recurring conceptual cautions follow from this body of work. First, the label extended Gross–Pitaevskii equation is not unique: it may denote a gradient-corrected mean-field PDE, an LHY-augmented droplet equation, a pairing-field equation, a condensate-plus-Bogoliubov dynamical scheme, or a matrix-valued quantum evolution. Second, several papers explicitly challenge over-identification of the extension with a single scalar condensate equation. In binary Bose-mixture droplets, the earlier eGPE is reinterpreted as the equation of motion for the pairing field (Hu et al., 2020). In dilute dipolar droplets, disagreement between experiment, eGPE, and quantum Monte Carlo shows that the current theoretical state-of-the-art description within the eGPE framework can miss significant correlations (Böttcher et al., 2019).

Taken together, these developments suggest that the extended Gross–Pitaevskii framework is best understood not as one equation but as a research program: a set of controlled modifications of the GP paradigm, each designed to encode a specific physical ingredient—finite range, quantum fluctuations, hard-core correlations, pairing, multicomponent coupling, nonlocal kernels, self-gravity, or exact solvability—while preserving as much as possible of the variational and dynamical structure that made the original GP equation useful.