🔢 mathematics

A fully decoupled and structure-preserving relaxation Crank--Nicolson finite element method for Gross--Pitaevskii--Poisson model

This paper proposes a fully decoupled, structure-preserving relaxation Crank--Nicolson finite element method for the Gross--Pitaevskii--Poisson system that achieves second-order temporal accuracy, optimal spatial convergence, and strict preservation of mass and modified discrete energy invariants.

Original authors: Dongqian Li, Huini Liu, Yin Yang, Peimeng Yin

Published 2026-03-24

📖 4 min read🧠 Deep dive

CC0 1.0

Original authors: Dongqian Li, Huini Liu, Yin Yang, Peimeng Yin

Original paper dedicated to the public domain under CC0 1.0 (http://creativecommons.org/publicdomain/zero/1.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to simulate a complex dance between two groups of tiny, invisible dancers (quantum particles) in a crowded room. These dancers don't just move on their own; they influence each other's steps, and they also create a "mood" in the room (an electric field) that changes how everyone else moves.

This is the Gross-Pitaevskii-Poisson (GPP) model. It's a set of equations used by physicists to understand ultracold plasmas—super-cold clouds of charged particles that behave like a single quantum wave.

The problem? These equations are incredibly messy. They are non-linear (small changes cause huge, unpredictable reactions) and coupled (everything depends on everything else). Solving them on a computer is like trying to predict the path of a hurricane while also calculating the wind speed of every single air molecule at the same time.

Here is what the authors of this paper did, explained simply:

1. The Problem: A Tangled Knot

Traditionally, to solve these equations, computers had to use "iterative" methods. Imagine trying to untangle a knot by pulling on one string, then checking if it's loose, then pulling another, then checking again. You have to repeat this loop thousands of times for every single step of time.

The Result: It's slow, computationally expensive, and sometimes the computer gets stuck in a loop or crashes.
The Old Way: To simulate two types of particles, the computer might have to solve seven different linear systems (groups of equations) for every single time step.

2. The Solution: The "Relaxation" Trick

The authors invented a new way to untangle the knot. They used a technique called Relaxation combined with a classic method called Crank-Nicolson.

Think of it like this:
Instead of trying to solve the whole tangled mess at once, they introduced two "Shadow Variables" (called $Z_+$ and $Z_-$ ).

The Metaphor: Imagine the dancers are wearing heavy, confusing costumes. Instead of calculating the weight of the costume inside the dance step, they take a snapshot of the costume's weight before the step, treat it as a fixed object, and then let the dancer move.
The Magic: By doing this, the complex, tangled equations become linear (straight lines) and fully decoupled (unlinked).

3. The Result: A Simple, Fast Dance

Because of this trick, the computer no longer needs to loop or guess.

The New Way: For every time step, the computer only needs to solve two simple, independent systems (one for each type of dancer) and one equation for the room's "mood" (the electric field).
The Benefit: It's like going from solving seven puzzles simultaneously to solving three easy puzzles one after another. It is much faster and much more stable.

4. Keeping the Physics Real (Structure-Preserving)

In physics simulations, a common problem is that computers make tiny errors that add up over time. Eventually, the simulation might show energy appearing out of nowhere or particles vanishing, which breaks the laws of physics.

The authors designed their method to be "Structure-Preserving."

The Metaphor: Imagine a bank account. Even if you make a tiny rounding error on a transaction, a good accounting system ensures that the total money in the bank never magically increases or decreases.
The Guarantee: Their method guarantees that the total mass (number of particles) and the total energy of the system remain exactly the same throughout the simulation, just like in the real world.

5. The Proof: It Works!

The authors didn't just guess; they did the math to prove:

Accuracy: The method is highly accurate (second-order in time, meaning if you cut the time step in half, the error drops by four times).
Efficiency: It is significantly faster than previous methods.
Reliability: They ran computer experiments (simulations) that showed the method perfectly preserved energy and mass, and the results matched theoretical predictions.

Summary

This paper presents a smart, fast, and reliable recipe for simulating complex quantum plasmas.

Old Recipe: Slow, tangled, requires constant checking and re-checking.
New Recipe: Uses "shadow variables" to untangle the math, making it fast and linear, while strictly obeying the laws of conservation (energy and mass).

It's a significant step forward for scientists who want to model ultracold matter without waiting days for a computer to finish the calculation or worrying that the simulation will drift away from reality.

1. Problem Statement

The paper addresses the numerical simulation of the coupled Gross–Pitaevskii–Poisson (GPP) system, which models the dynamics of ultracold plasmas consisting of two interacting quantum species (e.g., cations and anions). The system is governed by:

Two coupled nonlinear Schrödinger equations (Gross–Pitaevskii equations) describing the wave functions $\psi_+$ and $\psi_-$ .
A Poisson equation describing the self-consistent electrostatic potential $\phi$ generated by the charge densities $|\psi_+|^2$ and $|\psi_-|^2$ .

Key Challenges:

Strong Nonlinearity: The system contains cubic self-interaction terms ( $|\psi|^2\psi$ ), cross-interaction terms ( $|\psi_{\mp}|^2\psi_{\pm}$ ), and coupling to the Poisson potential.
Coupling: The wave functions and the potential are tightly coupled, making standard explicit or implicit schemes computationally expensive or unstable.
Structure Preservation: Physical models of this type conserve mass (for each species) and total energy. Standard numerical methods often fail to preserve these invariants at the discrete level, leading to long-term instability or unphysical results.
Efficiency: Existing methods often require solving large, coupled nonlinear algebraic systems or multiple linear systems per time step (e.g., splitting methods requiring up to 7 linear systems).

2. Methodology

The authors propose a Fully Decoupled and Structure-Preserving Relaxation Crank–Nicolson Finite Element Method (FEM).

A. Reformulation via Auxiliary Variables

To handle the nonlinearities, the authors introduce two auxiliary variables:
$Z_+ = |\psi_+|^2, \quad Z_- = |\psi_-|^2$
This reformulates the original nonlinear interaction terms (e.g., $|\psi_+|^2\psi_+$ ) into linear forms involving the auxiliary variables (e.g., $Z_+\psi_+$ ). This transforms the system into an equivalent form where the nonlinearities are "relaxed."

B. The Numerical Scheme

The method combines:

Time Discretization: A Relaxation Crank–Nicolson scheme. It uses a midpoint approximation for the time derivative and averages the auxiliary variables and wave functions.
Space Discretization: Standard Finite Element Method (FEM) using piecewise polynomial spaces ( $P_k$ ).

Key Features of the Scheme:

Fully Decoupled: At each time step, the system splits into independent linear problems.
- Two independent linear systems for $\psi_+$ and $\psi_-$ (solved in parallel).
- One Poisson equation for $\phi$ .
- Two linear updates for the auxiliary variables $Z_+$ and $Z_-$ .
No Iteration: Unlike traditional implicit methods, this scheme does not require iterative solvers (like Newton-Raphson) to handle nonlinearities because the relaxation technique linearizes the system.
Structure Preservation: The scheme is designed to exactly preserve:
- Discrete Mass: $M(\psi_{\pm}) = \int |\psi_{\pm}|^2 dx$ is conserved for each species.
- Modified Discrete Energy: A specific discrete energy functional (incorporating the auxiliary variables) is conserved exactly.

3. Key Contributions

Novel Algorithm: Development of a relaxation-based Crank–Nicolson FEM specifically for the coupled GPP system with multiple species. It reduces the computational complexity significantly compared to existing splitting or iterative methods.
Computational Efficiency: The method requires solving only three linear systems per time step (two for wave functions, one for the potential), whereas comparable splitting methods may require up to seven. It avoids the cost of nonlinear iterations entirely.
Rigorous Error Analysis:
- Proved the well-posedness and uniqueness of the fully discrete scheme.
- Established optimal error estimates:
  - Second-order accuracy in time ( $O(\tau^2)$ ).
  - Optimal $(k+1)$ -th order accuracy in space ( $O(h^{k+1})$ ) in the $L^2$ -norm for $P_k$ elements.
- The proof utilizes mathematical induction, projection-based error splitting, and careful handling of the interaction between multiple species and the Poisson coupling.
Conservation Laws: Rigorous proof that the discrete scheme preserves mass and a modified energy invariant, ensuring long-term stability and physical fidelity.

4. Results

The authors validated the method through extensive numerical experiments:

Convergence Rates:
- Spatial Convergence: Experiments with $P_1$ and $P_2$ elements confirmed the theoretical convergence rates of $O(h^2)$ and $O(h^3)$ respectively.
- Temporal Convergence: Experiments confirmed second-order convergence in time ( $O(\tau^2)$ ).
Conservation Properties: Numerical simulations showed that the total mass and energy errors remained at the level of machine precision (round-off error) over long time intervals, confirming the structure-preserving nature of the method.
Physical Phenomena:
- Successfully captured density wave patterns and instability phenomena (modulational instability) inherent in the GPP system.
- Demonstrated the ability to simulate complex dynamics, such as the splitting of solitons and the formation of dipole/quadrupole states under varying interaction parameters.

5. Significance

This work provides a robust, efficient, and theoretically sound numerical framework for simulating ultracold plasmas and multi-component Bose-Einstein condensates.

Practical Impact: By decoupling the system and eliminating the need for nonlinear iterations, the method significantly reduces computational cost, making it feasible to simulate larger domains and longer time scales.
Theoretical Impact: It extends the relaxation technique to coupled systems with Poisson coupling, a non-trivial extension that preserves the physical invariants crucial for long-time dynamics.
Generalizability: The methodology is shown to naturally extend to the coupled Gross–Pitaevskii equations (without the Poisson term), offering a unified approach for various quantum fluid models.

In summary, the paper presents a breakthrough in numerical methods for coupled quantum systems, balancing high-order accuracy, computational efficiency, and strict adherence to physical conservation laws.