Numerical stability of the Hyperbolic Formulation of the Constraint equations for $\mathbb{T}^3$ cosmological space-times

This paper investigates the numerical stability of an algebraic-hyperbolic formulation of the Einstein constraint equations for $\mathbb{T}^3$ cosmological space-times using a pseudo-spectral method, revealing that while the scheme suffers from unavoidable instabilities near FLRW space-times, it can be stable for Gowdy space-times and specific subclasses, offering a potential new approach for computing inhomogeneous initial data.

The Gist

The Big Picture: Building a Universe from Scratch

Imagine you are an architect trying to build a house. Before you can start the construction crew (the "evolution" of the universe), you need to lay down the perfect foundation. In the world of Einstein's General Relativity, this foundation is called Initial Data.

The problem is that the blueprints for this foundation (the "Constraint Equations") are incredibly complicated. They are like a massive, tangled knot of rules that say, "If the floor is curved this way, the walls must be curved that way, or the whole house collapses."

For decades, architects have used a specific tool called the Conformal Method to untangle this knot. It works great for building isolated houses (like binary black holes floating in empty space). But when you try to build a universe (which is a closed loop, like a video game world where you walk off the right edge and appear on the left), that old tool gets stuck. It's like trying to use a straight ruler to measure the curve of a donut.

The New Tool: The "Hyperbolic" Approach

The authors of this paper tried a different tool: the Algebraic-Hyperbolic Formulation (AHF).

Instead of trying to solve the whole knot at once (like the old method), this new approach treats the problem like a river flowing downstream. You start at one bank (a specific slice of space) and "evolve" the data step-by-step toward the other bank. It's much more flexible for closed universes because it doesn't need to know the shape of the whole river at once; it just needs to know how to flow from the next step.

They decided to test this new tool using a computer method called Pseudo-Spectral Fourier, which is like using a super-precise musical tuner to analyze the shape of the universe.

The Problem: The "Shaky Ladder"

When they ran their computer simulations, they hit a wall.

• The Good News: When they tested the tool on a specific type of universe called Gowdy (which is full of gravitational waves, like a bumpy road), the tool worked perfectly. The foundation was solid.
• The Bad News: When they tried to use it on a Perturbed FLRW universe (a universe that is mostly smooth and uniform, like a calm ocean with tiny ripples), the computer simulation exploded. The numbers went wild, the errors grew instantly, and the "foundation" collapsed.

It was as if they tried to climb a ladder, but every time they reached for the next rung, the ladder vibrated so violently they fell off.

The Investigation: Why Did the Ladder Shake?

The authors decided to act like detectives. They performed a Linear Stability Analysis, which is essentially asking, "If I nudge the system slightly, does it settle back down, or does it go into a panic?"

They discovered the culprit was a mismatch between the tool and the terrain:

1. The Tool: They were using a "Fourier" method. Think of this as a method that assumes everything is a smooth, repeating wave (like a sine wave). It's incredibly fast and accurate for smooth things.
2. The Terrain: The "Perturbed FLRW" universe they were trying to model has a specific mathematical property. When they looked at the "spectrum" (the list of frequencies) of their equations, they found that the numbers were sitting in a "danger zone."

The Analogy: Imagine trying to push a swing. If you push it at the right time (the stable zone), it goes higher smoothly. But if you push it at the wrong time (the unstable zone), it fights against you.
The authors proved that for the smooth, uniform universe (FLRW), the math of their "river flow" method combined with the "Fourier tuner" meant they were always pushing the swing at the wrong time. No matter how careful they were, the math guaranteed the simulation would become unstable. It wasn't a bug in their code; it was a fundamental law of the system they were using.

The Solution: Tweaking the Rules

Since they couldn't fix the math of the universe itself, they decided to tweak the rules of their construction game. They realized that the instability was caused by a specific part of the foundation (a field called $Y_i$, which represents how the layers of space are sliding past each other).

They proposed two new ways to build the foundation:

1. The "Fix the Current" Method: They forced the sliding part to be zero and calculated what the "energy current" must be to make that happen. It's like saying, "Okay, the floor can't slide, so the water pipes must be arranged in this specific way to compensate." This worked, but it meant they couldn't choose the water flow freely anymore.
2. The "Minimal Surface" Method: They imposed a rule that the layers of space must be "minimal" (like a soap bubble that has the least amount of surface area). This allowed them to solve the equations in a different way (using a "parabolic relaxation" method, which is like slowly heating a metal rod until it settles into the right shape). This worked beautifully and allowed them to keep the energy flows free.

The Conclusion

What did they learn?

• The new "River Flow" method (AHF) is a promising way to build initial data for universes, especially those with complex shapes.
• However, if you try to use it on a perfectly smooth, expanding universe with a specific computer method (Fourier), it will fail. The math is just too shaky.
• But! If you add a few extra rules (like forcing the layers to be minimal or fixing the sliding), the method becomes stable and works great.

The Takeaway:
The authors didn't just find a broken tool; they figured out why it broke and invented two new "safety harnesses" to make it work. This opens the door for scientists to simulate complex, inhomogeneous universes (universes with clumps of galaxies and voids) without getting stuck in the old, rigid methods. They turned a shaky ladder into a sturdy bridge, provided you follow the new safety guidelines.

Technical Summary

1. Problem Statement

The construction of initial data sets for the Einstein constraint equations in cosmological settings (specifically with compact $T^3$ topology) is a significant challenge.

• Limitations of Standard Methods: The traditional Lichnerowicz–York conformal method relies on elliptic partial differential equations (PDEs) and asymptotic flatness. In compact cosmological topologies, the lack of spatial infinity makes boundary conditions ill-defined, and the conformal decomposition becomes under-constrained.
• The Alternative: The Algebraic-Hyperbolic Formulation (AHF), introduced by Rácz, recasts the constraints into a first-order hyperbolic system. This allows for localized evolution and is naturally suited to periodic boundary conditions via spectral methods.
• The Core Issue: While AHF is theoretically promising, its numerical implementation using Fourier-based Method of Lines (MoL) exhibits pathological instabilities, particularly for perturbed Friedmann-Lemaître-Robertson-Walker (PFLRW) universes. The paper aims to diagnose the source of these instabilities and explore modifications to achieve stable numerical solutions.

2. Methodology

The authors employ a rigorous combination of numerical experimentation and linear stability analysis:

• Numerical Implementation:

  • Discretization: A pseudo-spectral method using the Discrete Fourier Transform (DFT) for angular derivatives ($x^1, x^2$) and a 4th-order Runge-Kutta (RK4) scheme for radial evolution ($r$).
  • Domain: $T^3$ topology ($S^1 \times T^2$) with periodic boundary conditions in all coordinates.
  • Test Cases: The scheme is tested against two benchmark spacetimes:
    1. $T^3$-Gowdy Space-times: Vacuum solutions with gravitational waves.
    2. Perturbed FLRW (PFLRW): Models representing inhomogeneous cosmologies with small scalar perturbations.
  • Error Metrics: Constraint violations (Hamiltonian and Momentum) are evaluated using 3D Fourier differentiation to measure the deviation from the exact constraints.

• Stability Analysis:

  • Linearization: The nonlinear AHF system is linearized around analytical solutions (FLRW and Gowdy).
  • Spectral Analysis: The authors analyze the eigenvalues of the semi-discretization matrix $L$ (arising from Fourier differentiation).
  • Stability Criteria: They apply the Lax Equivalence Theorem and Eigenvalue Stability criteria. Specifically, they check if the eigenvalues of $\Delta r L$ fall within the stability region of the RK4 time-stepping scheme.
  • Pseudospectrum Analysis: For non-normal matrices (common in hyperbolic systems), they utilize $\epsilon$-pseudospectrum analysis to assess transient growth and stability robustness.

3. Key Contributions and Results

A. Diagnosis of Instability in PFLRW

• Observation: The standard AHF + Fourier MoL approach fails to converge for PFLRW metrics. Constraint violations grow rapidly (anti-convergence) regardless of grid resolution.
• Theoretical Proof (Theorem 4.1): The authors prove that for any spacetime sufficiently close to FLRW, the instability is unavoidable under the specific combination of:
  1. Using Fourier-based MoL.
  2. Applying it to the standard AHF PDE system.
  3. Restricting the problem to PFLRW-like metrics.
• Mechanism: The linearized system for FLRW yields a semi-discretization matrix $L$ whose eigenvalues lie on the real axis (specifically, some have positive real parts). Since the stability region of any consistent Runge-Kutta method excludes a segment of the positive real axis near the origin, the numerical scheme is inherently unstable. This is directly linked to the hyperbolicity condition ($XZ < 0$) of the AHF, which is violated for FLRW metrics.

B. Conditional Stability in Gowdy Space-times

• Dependence on Parameters: Unlike FLRW, the stability of the Gowdy system depends on the time parameter $t$ and the radial coordinate.
• Hyperbolicity Link: The eigenvalues of the discretization matrix depend on the sign of the product $XZ$.
  • If $XZ > 0$ (violating hyperbolicity), eigenvalues are real and positive, leading to instability.
  • If $XZ < 0$ (satisfying hyperbolicity), eigenvalues are purely imaginary, falling within the stability region of implicit schemes (like Crank-Nicolson) and potentially stable for explicit schemes depending on step size.
• Result: The scheme is stable for Gowdy spacetimes only when the evolution stays within the domain where the hyperbolicity condition holds.

C. Proposed Modifications for Stable Initial Data

Recognizing the limitations of the standard formulation, the authors propose two modifications to the AHF system to enforce stability ($Y_i = 0$) and construct new initial data sets:

1. Type 1: Algebraic Determination of Currents

   • Strategy: Impose the condition $Y_i = 0$ (vanishing tangential extrinsic curvature). This converts the momentum constraint equation into an algebraic relation that determines the tangential current $J^{(||)}_i$ directly from the geometry.
   • Outcome: This reduces the system to one differential equation and three algebraic equations. Numerical tests show constraint violations of $\sim 10^{-13}$ for PFLRW, proving viability, though it restricts the freedom of the energy-momentum sources.

2. Type 2: Parabolic Relaxation

   • Strategy: Impose $Y_i = 0$, $\beta^i = 0$, and minimal surface foliation ($H^i_i = 0$). This transforms the problem of finding the initial field $X$ into a second-order elliptic PDE for an arbitrary function $F$.
   • Implementation: To avoid solving the elliptic PDE directly, they introduce a "parabolic relaxation" method (solving a heat-equation-like evolution in a fictitious time $t$) to converge to the steady-state solution.
   • Outcome: This approach allows for fully free energy sources ($\rho, J$) while maintaining stability. Constraint violations are found to be $\sim 10^{-12}$.

4. Significance and Implications

• Theoretical Insight: The paper provides a definitive explanation for why hyperbolic formulations of constraints fail in standard cosmological perturbation theory when paired with spectral methods. It highlights that the spectral properties of Fourier differentiation matrices are incompatible with the eigenvalue structure of the linearized AHF for FLRW backgrounds.
• Practical Utility: The proposed modifications (Type 1 and Type 2) offer a viable path forward for generating initial data in compact cosmologies without resorting to elliptic solvers.
• Future Directions: The work suggests that while the standard AHF is unstable for generic inhomogeneous cosmologies, constrained versions (fixing specific geometric degrees of freedom) can yield stable, periodic initial data sets. This opens new avenues for studying fully relativistic inhomogeneous universes and backreaction effects in numerical relativity.

In summary, the paper demonstrates that while the direct application of the Algebraic-Hyperbolic Formulation to perturbed FLRW universes is numerically unstable due to spectral incompatibility, strategic geometric constraints can restore stability, enabling the construction of valid initial data for compact cosmological models.