pith. sign in

arxiv: 2309.02946 · v5 · submitted 2023-09-06 · 🌀 gr-qc

Numerical stability of the Hyperbolic Formulation of the Constraint equations for mathbb{T}³ cosmological space-times

Alejandro Estrada-Llesta , Cristhian Martinez-Duarte , Leon Escobar-Diaz This is my paper

Pith reviewed 2026-05-24 06:56 UTC · model grok-4.3

classification 🌀 gr-qc
keywords numerical relativityconstraint equationscosmological initial datahyperbolic formulationslinear stability analysisT^3 topologypseudo-spectral methods
0
0 comments X

The pith

Linear analysis shows algebraic-hyperbolic constraint equations are numerically unstable near FLRW for T^3 cosmologies.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper examines an algebraic-hyperbolic version of the Einstein constraint equations for constructing initial data on T^3 cosmological spacetimes. A pseudo-spectral method of lines based on the discrete Fourier transform is applied to the system and produces growing instabilities. Linear stability analysis establishes that these instabilities cannot be removed for any spacetime close enough to a homogeneous FLRW background. The same discretization can remain stable for certain Gowdy spacetimes when the initial time is chosen appropriately, and selected subclasses of the formulation appear to support stable evolution.

Core claim

Through linear stability analysis we prove that the instabilities are unavoidable for any space-time sufficiently close to FLRW while we find that this approach can be stable for Gowdy space-times depending on the initial time choice. Additionally, we present numerical evidence that certain subclasses of the algebraic-hyperbolic formulation, when combined with a Fourier-based method of lines, are numerically stable.

What carries the argument

The algebraic-hyperbolic formulation of the Einstein constraint equations, discretized by a pseudo-spectral method of lines and subjected to linear stability analysis.

If this is right

  • The standard algebraic-hyperbolic formulation cannot be used reliably for initial data construction near FLRW with Fourier-based discretizations.
  • Stability for Gowdy spacetimes requires a suitable choice of initial time slice.
  • Certain subclasses of the formulation remain candidates for stable numerical initial data sets in inhomogeneous cosmologies.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Alternative spatial discretizations such as finite differences might evade the reported instabilities even if the continuous system is marginally stable.
  • The dependence on background geometry suggests that hyperbolic constraint formulations must be re-examined for each class of cosmological symmetry.
  • The stable subclasses identified numerically could be tested on spacetimes with different topologies or matter sources to check broader applicability.

Load-bearing premise

The spacetime is close enough to FLRW that linear perturbations around it capture the dominant numerical behavior.

What would settle it

A long-term stable numerical evolution of the full algebraic-hyperbolic system for a spacetime arbitrarily close to FLRW, using the same Fourier method of lines, would falsify the unavoidability claim.

Figures

Figures reproduced from arXiv: 2309.02946 by Alejandro Estrada-Llesta, Cristhian Martinez-Duarte, Leon Escobar-Diaz.

Figure 1
Figure 1. Figure 1: Hamiltonian constraint evaluation depending on the number of nodes considered in a [PITH_FULL_IMAGE:figures/full_fig_p011_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: First (left) and second (right) components of the Momentum constraint evaluation [PITH_FULL_IMAGE:figures/full_fig_p012_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Convergence test for Gowdy metric. (left) Convergence test for N = 32. (right) Constraint violation of the solution reconstructed from the evolution process as a function of the number of nodes and number of evolution radial steps Nr,evol = F actor N. this filtering strategy and explore it in Appendix (A). For this and PFLRW case we show results for the (1/2)−filter. However, this result cannot be fully un… view at source ↗
Figure 4
Figure 4. Figure 4: Convergence test for GowdyRX metric of eq. ( [PITH_FULL_IMAGE:figures/full_fig_p015_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Convergence test for PFLRW metric of eq. ( [PITH_FULL_IMAGE:figures/full_fig_p015_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Eigenvalues of ∆r times the Matrix L, for L = 0.5, superimposed on the stability region of (from left to right) Crank Nicholson and Implicit Euler. The eigenvalues were scaled by a factor θ = 102 for visualization purposes. As illustrated in figure 6, half of the scaled eigenvalues of L lie outside the stability region of all the integration schemes considered in this work. Consequently, since the eigenval… view at source ↗
Figure 7
Figure 7. Figure 7: show the distribution of the eigenvalues of the angular discretization operator when the coefficients in equation 37 are frozen using the two different methods mentioned above, for ∆r = 10−2N −1 . In both cases, the numerically calculated eigenvalues fall outside the stability region, regardless of the chosen radial step, just as in the unperturbed case. This result indicates that, similar to the FLRW syst… view at source ↗
Figure 8
Figure 8. Figure 8: Behaviour of the numerical solution of the linearized AHF equations using CN for [PITH_FULL_IMAGE:figures/full_fig_p020_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Graph of the inequality XZ in the domain [0, 1] × [0, 1]. The blue region indicates where the hyperbolicity condition is satisfied (XZ < 0). The horizontal lines illustrate the three possibilities for the parameter t: the white line at t = 0.2 represents a value of t for which the hyperbolicity condition is always satisfied in all the radial domain [0, 1]. Whereas the green line at t = 0.4 represents a val… view at source ↗
Figure 10
Figure 10. Figure 10: On the left, Eigenvalues of ∆r times the Matrix L, with coefficients frozen by taking their average in the mesh, superimposed on the stability region of Crank Nicholson. On the right, the eigenvalues are superimposed on the stability region of the Implicit Euler method. The eigenvalues were scaled by a factor θ = 102 for visualization purposes. (a) (b) [PITH_FULL_IMAGE:figures/full_fig_p022_10.png] view at source ↗
Figure 11
Figure 11. Figure 11: Behaviour of the numerical solution of the AHF equations with CN for Gowdy space [PITH_FULL_IMAGE:figures/full_fig_p022_11.png] view at source ↗
Figure 12
Figure 12. Figure 12: On the left, eigenvalues of ∆r times the matrix L, for L = 0.5 and t = 0.2, superim￾posed on the stability region of Crank-Nicholson and on the right superimposed on the stability region of Implicit Euler. The eigenvalues were scaled by a factor θ = 102 for visualization pur￾poses [PITH_FULL_IMAGE:figures/full_fig_p023_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: On the right, eigenvalues of random perturbations [PITH_FULL_IMAGE:figures/full_fig_p024_13.png] view at source ↗
Figure 14
Figure 14. Figure 14: Distance of ϵ−pseudospectrum to the Crank-Nicholson stability region as a function of ϵ for different grid sizes. Filled circles indicate measured distances, and the dashed line indicates the linear regression line [PITH_FULL_IMAGE:figures/full_fig_p024_14.png] view at source ↗
Figure 15
Figure 15. Figure 15: Behavior of the numerical solution of the AHF equations with CN for Gowdy spacetime [PITH_FULL_IMAGE:figures/full_fig_p025_15.png] view at source ↗
Figure 16
Figure 16. Figure 16: Solution of X at some values of r = rk = −L+k∆r for different initial conditions. The first row show stages of the solutions for the original initial condition, i.e , without modification. The other three rows show the solution for the modifications of the initial conditions presented in eqs. (50) to (52). The color bar represents the values of X in the scale of X − Xbackground where Xbackground is the un… view at source ↗
Figure 17
Figure 17. Figure 17: Plots achieved in a grid of 32 nodes, RK4 as the integrator and a [PITH_FULL_IMAGE:figures/full_fig_p030_17.png] view at source ↗
Figure 18
Figure 18. Figure 18: Plots were achieved in a grid of 32 nodes, RK4 as the integrator and a [PITH_FULL_IMAGE:figures/full_fig_p030_18.png] view at source ↗
Figure 19
Figure 19. Figure 19: (a) Values of the parameters p and α that satisfies that the maximum of the derivative of σ(k) is at k = q(N/2). For each q there is not only a point in the parameter space that satisfies this condition, but a curve in the plane p − ln α. (b) Results for q = 2/3 compared with the (2/3)−filter. Increase the value of p results in a higher value of the derivative of σ(k) producing a faster screening of the s… view at source ↗
Figure 20
Figure 20. Figure 20: Result of the filtering parameter exploration for the PFLRW metric in square grids [PITH_FULL_IMAGE:figures/full_fig_p036_20.png] view at source ↗
Figure 21
Figure 21. Figure 21: Result of the filter parameter exploration applied to the Gowdy gauge wave metric [PITH_FULL_IMAGE:figures/full_fig_p037_21.png] view at source ↗
read the original abstract

In this work, we study of the algebraic-hyperbolic formulation of the Einstein constraint equations for numerically constructing initial data sets for inhomogeneous cosmological space-times with $\mathbb{T}^3$ topology. We implement a pseudo-spectral method of lines based on the discrete Fourier transform and find that the scheme exhibits pathological instabilities. Through linear stability analysis, we prove that the instabilities are unavoidable for any space-time sufficiently close to FLRW while we find that this approach can be stable for Gowdy space-times depending on the initial time choice. Additionally, we present numerical evidence that certain subclasses of the algebraic-hyperbolic formulation, when combined with a Fourier-based method of lines, are numerically stable, thus offering a potential new path for computing initial data sets for inhomogeneous cosmological space-times.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The manuscript studies the algebraic-hyperbolic formulation of the Einstein constraint equations on T^3 topologies for cosmological initial data. A pseudo-spectral method of lines based on the discrete Fourier transform is implemented and shown to exhibit pathological instabilities. Linear stability analysis is invoked to prove that instabilities are unavoidable for any spacetime sufficiently close to FLRW, while stability is possible for Gowdy spacetimes depending on the initial time slice. Numerical evidence is presented that certain subclasses of the formulation remain stable under the same discretization.

Significance. If the linear analysis rigorously establishes persistence of instabilities under small deviations from FLRW, the result would be significant for numerical relativity: it would identify a structural limitation of the algebraic-hyperbolic system near the most common cosmological backgrounds and point to stable subclasses as a practical alternative. The combination of analytic proof and targeted numerical tests is a positive feature when the mode-coupling question is resolved.

major comments (2)
  1. [linear stability analysis (abstract and implementation paragraph)] Abstract and linear stability analysis section: the central claim that instabilities are 'unavoidable for any space-time sufficiently close to FLRW' is load-bearing. The analysis appears to be performed on the exact homogeneous FLRW background (where DFT modes remain eigenmodes). Small spatial inhomogeneities introduce off-diagonal coupling in the mode equations; without a first-order perturbation of the background coefficients to verify that unstable eigenvalues persist, the extension to nearby inhomogeneous spacetimes is not established.
  2. [numerical results section] Numerical evidence paragraph: the reported stability for certain subclasses is presented as offering a 'potential new path,' yet no quantitative error norms, convergence rates, or comparison against a known exact solution (e.g., perturbed FLRW) are referenced. This weakens the practical utility claim.
minor comments (2)
  1. [abstract] The abstract states 'we prove' but does not indicate whether the linear analysis includes a systematic expansion in spatial inhomogeneity amplitude; this should be clarified in the main text.
  2. [formulation section] Notation for the algebraic-hyperbolic system variables should be introduced once with explicit reference to the original formulation paper.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for their thorough review and valuable comments on our manuscript. We address each major comment below and indicate the revisions we plan to make.

read point-by-point responses

  1. Referee: Abstract and linear stability analysis section: the central claim that instabilities are 'unavoidable for any space-time sufficiently close to FLRW' is load-bearing. The analysis appears to be performed on the exact homogeneous FLRW background (where DFT modes remain eigenmodes). Small spatial inhomogeneities introduce off-diagonal coupling in the mode equations; without a first-order perturbation of the background coefficients to verify that unstable eigenvalues persist, the extension to nearby inhomogeneous spacetimes is not established.

    Authors: We agree with the referee that our linear stability analysis is conducted on the homogeneous FLRW background, where the discrete Fourier modes are exact eigenmodes. The extension to spacetimes sufficiently close to FLRW relies on the expectation that the unstable eigenvalues will persist under small perturbations of the coefficients. However, as noted, a first-order perturbative analysis of the coupled mode equations would provide a more rigorous justification. We will revise the manuscript to explicitly acknowledge this limitation of the current analysis and to discuss the conditions under which the instabilities are expected to persist, including a sketch of how a perturbative treatment could be carried out. revision: yes

  2. Referee: Numerical evidence paragraph: the reported stability for certain subclasses is presented as offering a 'potential new path,' yet no quantitative error norms, convergence rates, or comparison against a known exact solution (e.g., perturbed FLRW) are referenced. This weakens the practical utility claim.

    Authors: We acknowledge that the numerical evidence for the stability of certain subclasses could be strengthened with quantitative measures. In the revised version, we will include L^2 error norms, demonstrate convergence rates with respect to the number of Fourier modes, and provide comparisons to exact solutions for the Gowdy spacetimes and other stable cases where possible. This will better support the claim of numerical stability and the potential utility of these subclasses. revision: yes

Circularity Check

0 steps flagged

No circularity: linear analysis is independent of numerics

full rationale

The paper derives its main result (instabilities unavoidable near FLRW) via linear stability analysis of the algebraic-hyperbolic system and separately reports numerical evidence for stability in Gowdy cases. No quoted step shows a parameter fitted to data then renamed as prediction, a self-definitional loop, or a load-bearing self-citation chain that reduces the claim to its own inputs. The derivation remains self-contained as a direct mathematical analysis of the PDE system.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The work rests on the validity of the algebraic-hyperbolic reformulation of the Einstein constraints and on the applicability of linear stability analysis to the discretized system; no free parameters or invented entities are indicated in the abstract.

axioms (2)
  • domain assumption The Einstein constraint equations admit an algebraic-hyperbolic formulation suitable for numerical evolution.
    Central object of study; invoked throughout the abstract.
  • domain assumption Linear stability analysis of the semi-discrete system correctly captures the behavior of the full nonlinear numerical scheme near FLRW.
    Used to prove the main instability claim.

pith-pipeline@v0.9.0 · 5671 in / 1285 out tokens · 26214 ms · 2026-05-24T06:56:04.949731+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

31 extracted references · 31 canonical work pages

  1. [1]

    Alcubierre

    M. Alcubierre. Introduction to 3+ 1 numerical relativity , volume 140. OUP Oxford, 2008

    work page 2008

  2. [2]

    Alcubierre, G

    M. Alcubierre, G. Allen, C. Bona, D. Fiske, T. Goodale, F. S. Guzm´ an, I. Hawke, S. H. Hawley, S. Husa, M. Koppitz, C. Lechner, D. Pollney, D. Rideout, M. Salgado, E. Schnetter, E. Seidel, H. aki Shinkai, D. Shoemaker, B. Szil´ agyi, R. Takahashi, and J. Winicour. Towards 31 standard testbeds for numerical relativity. Classical and Quantum Gravity , 21(2...

    work page 2003

  3. [3]

    T. W. Baumgarte and S. L. Shapiro. Numerical relativity: solving Einstein’s equations on the computer. Cambridge University Press, 2010

    work page 2010

  4. [4]

    Beyer, L

    F. Beyer, L. Escobar, and J. Frauendiener. Asymptotics of solutions of a hyperbolic formu- lation of the constraint equations. Classical and Quantum Gravity , 34(20):205014, 2017

    work page 2017

  5. [5]

    Beyer, L

    F. Beyer, L. Escobar, J. Frauendiener, and J. Ritchie. Numerical construction of initial data sets of binary black hole type using a parabolic-hyperbolic formulation of the vacuum constraint equations. Classical and Quantum Gravity , 36(17):175005, 2019

    work page 2019

  6. [6]

    Beyer, J

    F. Beyer, J. Frauendiener, and J. Ritchie. Asymptotically flat vacuum initial data sets from a modified parabolic-hyperbolic formulation of the einstein vacuum constraint equations. Physical Review D, 101(8):084013, 2020

    work page 2020

  7. [7]

    Canuto, M

    C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral methods: evolution to complex geometries and applications to fluid dynamics . Springer Science & Business Media, 2007

    work page 2007

  8. [8]

    Csuk´ as and I

    K. Csuk´ as and I. R´ acz. Numerical investigations of the asymptotics of solutions to the evolutionary form of the constraints. Classical and Quantum Gravity , 37(15):155006, 2020

    work page 2020

  9. [9]

    Csuk´ as and I

    K. Csuk´ as and I. R´ acz. Is it possible to construct asymptotically flat initial data using the evolutionary forms of the constraints? Physical Review D, 107(8):084013, 2023

    work page 2023

  10. [10]

    M. P. Do Carmo and J. Flaherty Francis. Riemannian geometry, volume 6. Springer, 1992

    work page 1992

  11. [11]

    Dodelson and F

    S. Dodelson and F. Schmidt. Modern cosmology. Academic press, 2020

    work page 2020

  12. [12]

    R. Durrer. The Cosmic Microwave Background . Cambridge University Press, 2008

    work page 2008

  13. [13]

    Garfinkle and L

    D. Garfinkle and L. Mead. Cosmological initial data for numerical relativity. Phys. Rev. D , 102:044022, Aug 2020

    work page 2020

  14. [14]

    D. J. Higham and L. N. Trefethen. Stiffness of odes. BIT Numerical Mathematics , 33:285– 303, 1993

    work page 1993

  15. [15]

    D. A. Kopriva. Implementing spectral methods for partial differential equations: Algorithms for scientists and engineers . Springer Science & Business Media, 2009

    work page 2009

  16. [16]

    Ma and E

    C. Ma and E. Bertschinger. Cosmological perturbation theory in the synchronous and con- formal newtonian gauges. Astrophysical Journal, 455(1), 12 1995

    work page 1995

  17. [17]

    H. J. Macpherson, P. D. Lasky, and D. J. Price. Inhomogeneous cosmology with numerical relativity. Phys. Rev. D , 95:064028, Mar 2017

    work page 2017

  18. [18]

    H. J. Macpherson, D. J. Price, and P. D. Lasky. Einstein’s universe: Cosmological structure formation in numerical relativity. Phys. Rev. D , 99:063522, Mar 2019

    work page 2019

  19. [19]

    F. K. Morton and B. MJ. Lax-stability vs. eigenvalue stability of spectral methods. 32

    work page

  20. [20]

    P. J. E. Peebles. The large-scale structure of the universe , volume 98. Princeton university press, 2020

    work page 2020

  21. [21]

    I. R´ acz. Cauchy problem as a two-surface based ‘geometrodynamics’.Classical and Quantum Gravity, 32(1):015006, 2014

    work page 2014

  22. [22]

    I. R´ acz. Is the bianchi identity always hyperbolic? Classical and Quantum Gravity , 31(15):155004, 2014

    work page 2014

  23. [23]

    I. R´ acz. Constraints as evolutionary systems.Classical and Quantum Gravity, 33(1):015014, 2015

    work page 2015

  24. [24]

    S. C. Reddy and L. N. Trefethen. Stability of the method of lines. Numerische Mathematik, 62(1):235–267, 1992

    work page 1992

  25. [25]

    Ringstr¨ om

    H. Ringstr¨ om. Cosmic censorship for gowdy spacetimes. Living Reviews in Relativity , 13:1– 59, 2010

    work page 2010

  26. [26]

    Stoer, R

    J. Stoer, R. Bulirsch, R. Bartels, W. Gautschi, and C. Witzgall. Introduction to numerical analysis, volume 1993. Springer, 1980

    work page 1993

  27. [27]

    M. Taylor. Partial differential equations II: Qualitative studies of linear equations , volume

    work page

  28. [28]

    Springer Science & Business Media, 2013

    work page 2013

  29. [29]

    L. N. Trefethen. Spectral methods in MATLAB . SIAM, https://doi.org/10.1137/1.9780898719598, 2000

    work page doi:10.1137/1.9780898719598 2000

  30. [30]

    L. N. Trefethen and M. Embree. Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators. Princeton University Press, 2005

    work page 2005

  31. [31]

    R. M. Wald. General Relativity. Chicago Univ. Pr., Chicago, USA, 1984. 33 A Aliasing error and filtering strategy Aliasing error emerges from the non-linearities of the PDE and is due to the finiteness of the grid (see [7, 15]). Since our equations (eqs. (6) to (8)) are non-linear, determine the relevance of this error is important to obtain accurate nume...

    work page 1984