Efficient High-order Mass-conserving and Energy-balancing Schemes for Schrödinger-Poisson Equations

This paper proposes and validates efficient, high-order numerical schemes that combine implicit-explicit Runge-Kutta time-stepping with Fourier collocation and relaxation-based post-processing to rigorously conserve mass and energy (or satisfy energy balance equations) for Schrödinger-Poisson systems, including those with time-varying coefficients relevant to cosmological simulations.

The Gist

Imagine you are trying to simulate the behavior of a vast, invisible cloud of "fuzzy" particles that make up the universe's dark matter. These particles don't just sit still; they ripple, crash into each other, and swirl around like a cosmic fluid. To predict how they move, scientists use a complex mathematical recipe called the Schrödinger-Poisson (SP) equations.

Think of these equations as the "laws of physics" for this digital universe. They tell us two very important things must always happen:

1. Mass Conservation: You can't create or destroy matter out of thin air. The total amount of "stuff" in your simulation must stay exactly the same.
2. Energy Balance: In a stable universe, energy stays constant. But in our real, expanding universe (like the one we live in), energy doesn't stay perfectly still; it shifts and balances in a specific, predictable way as space stretches.

The Problem: The Leaky Bucket

For decades, when scientists tried to run these simulations on computers, they faced a frustrating problem. The computer algorithms they used were like a leaky bucket.

Every time the computer took a tiny step forward in time to calculate the next moment, it would accidentally "spill" a little bit of mass or energy. At first, the spill is tiny and unnoticeable. But over millions of steps (simulating billions of years), the bucket empties. The simulation might show the universe losing all its matter or gaining infinite energy, making the results scientifically useless.

To fix this, previous methods tried to build a "perfect bucket" by making the math incredibly complex and slow. These methods were like trying to carry water in a bucket made of solid gold—very sturdy, but so heavy and slow that you could barely move.

The Solution: The "Relaxation" Trick

The authors of this paper, Manvendra Rajvanshi and David Ketcheson, proposed a clever, lightweight solution. They didn't try to build a perfect bucket from the start. Instead, they used a technique called Relaxation.

Here is the analogy:
Imagine you are walking a tightrope (the correct path of physics). You take a step forward using a standard walking method (a fast, efficient computer algorithm). Because you're human (or a computer), you might wobble slightly off the rope.

• The Old Way: Try to walk perfectly straight from the start, which is slow and exhausting.
• The New "Relaxation" Way: Take a quick, fast step (even if you wobble a bit). Then, immediately after the step, you pause, look at where you are, and gently nudge yourself back onto the tightrope.

This "nudge" is the Relaxation. It's a quick calculation that checks: "Did we lose any mass? Did we gain too much energy?" If the answer is yes, the computer makes a tiny adjustment to fix it instantly.

Two Types of Nudges

The paper tests two specific ways to do this nudge:

1. Multiple Relaxation (MR): This is like having two safety nets. You check the mass and the energy separately and adjust both. It's very accurate, but sometimes the math gets tricky, and the computer might get stuck trying to find the perfect adjustment.
2. Projection Relaxation (PR): This is the "smart" nudge. It projects the solution directly onto the "safe zone" where mass and energy are correct. It's like using a laser guide to snap you back to the line. The authors found this method is faster, more reliable, and rarely gets stuck.

Why This Matters for the Universe

The paper shows that this method works beautifully, even in the most difficult scenarios:

• Stable Systems: It keeps the "bucket" perfectly full, even after simulating complex interactions.
• Expanding Universe: It correctly handles the "energy balance" of an expanding universe (where energy isn't constant but follows a specific rule), ensuring the simulation doesn't drift off course.
• 3D Cosmology: They tested this on a massive 3D simulation of the universe. Without the fix, the simulation would have drifted and become inaccurate. With the "Relaxation" nudge, the simulation stayed true to physics, producing realistic "clumps" of matter (galaxies and dark matter halos) that look exactly like what we see in the real sky.

The Bottom Line

This paper introduces a "plug-and-play" upgrade for astrophysicists. You don't have to throw away your favorite, fast simulation tools. You just add this "Relaxation" step at the end of every calculation.

It's like giving a race car a self-correcting steering system. The car can still go fast (using efficient algorithms), but now it stays perfectly on the track (conserving mass and energy) without crashing. This allows scientists to run longer, more accurate simulations of the universe's history, helping us understand how galaxies formed and how dark matter shapes our cosmos.

Technical Summary

1. Problem Statement

The Schrödinger-Poisson (SP) equations are fundamental for modeling nonlinear phenomena in fields ranging from semiconductor physics to cosmology (specifically fuzzy dark matter and self-gravitating Bose-Einstein condensates). The system is given by:
$$\psi_t - i p(t) \nabla^2 \psi + i q(t) V \psi = 0$$
$$\nabla^2 V = |\psi|^2 - \mu$$

Key Physical Challenges:

• Conservation Laws: The system possesses a conserved total mass ($\mu = \|\psi\|^2$). The total energy ($E$) is conserved if coefficients $p$ and $q$ are constant. However, in cosmological applications (expanding universe), $p(t)$ and $q(t)$ vary, and energy satisfies a balance equation (the Layzer-Irvine equation) rather than a strict conservation law.
• Numerical Limitations: Existing numerical schemes often struggle to preserve these invariants.
  • Operator-splitting methods typically conserve mass but fail to conserve energy.
  • Fully implicit schemes (e.g., Crank-Nicolson) can conserve both but are limited to second-order accuracy in time and require solving expensive nonlinear systems.
  • High-order methods often sacrifice conservation properties.

The authors aim to develop high-order in time numerical schemes that strictly conserve mass and satisfy the energy balance equation (or conservation law) without the prohibitive cost of solving large nonlinear systems at every step.

2. Methodology

The proposed methodology combines Implicit-Explicit (ImEx) Runge-Kutta (RK) time integration with Relaxation techniques and Fourier collocation in space.

A. Spatial Discretization

• Pseudospectral Fourier Collocation: The spatial derivatives are approximated using Fourier spectral methods.
• Skew-Hermitian Property: The discrete derivative operator $D$ is constructed to be skew-Hermitian ($D^\dagger = -D$). The authors prove that this property ensures the semi-discrete system naturally satisfies discrete analogs of the mass conservation and energy balance laws, provided the time integration is handled correctly.

B. Time Integration: ImEx RK

• The stiff term (linear diffusion, $\nabla^2 \psi$) is treated implicitly, while the non-stiff term (nonlinear potential, $V\psi$) is treated explicitly.
• This allows the use of high-order schemes (e.g., 3rd and 4th order ARK methods) while only requiring the solution of linear systems (via Fast Fourier Transforms, FFTs) rather than nonlinear systems.

C. Relaxation Techniques

To enforce conservation on top of the ImEx RK step, the authors employ two relaxation strategies that perturb the solution back onto the conservative manifold after each time step:

1. Multiple Relaxation (MR):

   • Uses two embedded ImEx methods with different weights to generate update vectors.
   • Solves a system of two algebraic equations to find relaxation parameters ($\gamma$) that enforce both mass and energy conservation simultaneously.
   • Drawback: Requires solving a 2x2 system; convergence can fail for large time steps.

2. Projection Relaxation (PR):

   • First, the solution is orthogonally projected onto the mass-conserving manifold.
   • Second, a scalar relaxation parameter is found to enforce energy conservation (or the energy balance equation) while staying on the mass manifold.
   • Advantage: Requires solving only one scalar nonlinear equation, making it more robust and computationally cheaper than MR.

D. Handling Time-Varying Coefficients (Cosmology)

For cases where $p(t)$ and $q(t)$ vary (expanding universe), the method enforces the energy balance equation rather than strict energy conservation.

• The discrete energy change over a time step is calculated via an integral of the balance law terms ($\frac{dp}{dt}E_k - \frac{1}{2}\frac{dq}{dt}E_v$).
• The relaxation parameter is adjusted so that the numerical solution's energy change exactly matches this calculated integral.

3. Key Contributions

1. High-Order Conservation: The paper demonstrates that relaxation techniques can be successfully combined with high-order ImEx RK schemes to achieve mass and energy conservation (or balance) at 3rd and 4th order accuracy in time.
2. Generalizability: The approach is modular; it can be applied as a post-processing step to any time-stepping scheme, provided the spatial discretization preserves the skew-Hermitian property.
3. Efficiency: By using ImEx schemes, the method avoids solving nonlinear systems. The implicit step only requires linear solves via FFTs, which are computationally efficient.
4. Cosmological Applicability: The framework explicitly handles time-varying coefficients, enabling accurate long-term simulations of expanding universes where energy is not conserved but follows a specific balance law.

4. Results

The authors validated the method through three numerical experiments:

• 2D Self-Gravitating Gaussians (Constant Coefficients):

  • Conservation: Relaxation methods (MR and PR) maintained mass and energy to machine precision ($10^{-14}$), whereas the baseline method showed significant drift.
  • Accuracy: PR maintained the 3rd-order convergence rate of the baseline scheme while reducing absolute errors. MR showed slightly degraded density accuracy.
  • Cost: Relaxation increased computational time by a factor of ~2–3 compared to the baseline, but this is negligible compared to fully implicit nonlinear solvers. PR was faster and more robust than MR.

• 2D Sine-Wave Collapse (Time-Varying Coefficients):

  • Balance Law: The PR method successfully enforced the energy balance equation, keeping the "energy violation" near machine epsilon, while the baseline method violated the balance law significantly.
  • Convergence: The 4th-order ImEx scheme maintained its 4th-order convergence rate when combined with relaxation.
  • Solution Quality: Visually, relaxation methods produced sharper interference patterns and cores that matched the reference solution (computed with tiny time steps) better than the non-relaxed baseline.

• 3D Cosmological Simulation:

  • Scale: Simulated a 3D universe ($256^3$ grid) from redshift $z=127$ to $z=0$.
  • Impact: While visual differences in density plots were subtle, the power spectra revealed significant differences. The PR method with a coarse time step ($\Delta t = 5 \times 10^{-5}$) produced results within 1–3 orders of magnitude of the highly refined reference solution, whereas the baseline method with the same step size deviated significantly.
  • Conclusion: Conservation properties are critical for capturing the correct statistical properties (halo formation, power spectrum) in long-term cosmological simulations.

5. Significance

This work bridges the gap between high-order accuracy and structure preservation in SP simulations.

• For Cosmology: It provides a robust tool for simulating Fuzzy Dark Matter and other scalar field models where long-term energy balance is crucial for the formation of solitonic cores and halo structures.
• For Numerical Mathematics: It validates the "Relaxation" paradigm as a versatile, low-cost add-on that can retrofit existing high-order time integrators with exact conservation properties, avoiding the complexity of developing new implicit conservative schemes from scratch.
• Recommendation: The authors recommend Projection Relaxation (PR) for practical applications, particularly in 3D, due to its superior robustness (solving a scalar equation vs. a system) and lower computational overhead compared to Multiple Relaxation.