A symplectic Python integration of the Friedmann equations recovers the age of the Universe to within 0.1 per cent of Planck 2018

A new paper published today in European Journal of Physics reconstructs the age of the Universe from first principles using a symplectic Python implementation of the Friedmann equations, agreeing with the Planck 2018 value of 13.8 Gyr to within 0.1 per cent. The team, led by Homer Dávila Gutiérrez, FRAS, is affiliated with the Master’s Programme in Astrophysics and Astronomy at Universidad Internacional de La Rioja (UNIR). The work compares seven distinct dark-energy regimes in a single integration framework, from concordance ΛCDM to phantom energy and the Big Rip. ✨

European Journal of Physics, the journal published by IOP Publishing on behalf of the Institute of Physics (IOP) and the European Physical Society (EPS), released today, 21 May 2026, the article «Numerical implementation of flat FLRW models of cosmic expansion with Planck 2018 cosmological parameters», authored by Homer Dávila Gutiérrez, FRAS, Julio Carlos Bertua Marasca, O. A. Martínez Osorio, J. I. Mosquera Hadatty and R. C. Videla Ricci. The team is affiliated with the Master’s Programme in Astrophysics and Astronomy at Universidad Internacional de La Rioja (UNIR), with the first author holding an additional affiliation at SKYCR, San José, Costa Rica. The full reference is Eur. J. Phys. 47 (2026) 035605, DOI 10.1088/1361-6404/ae6377, open access under Creative Commons Attribution 4.0.

A 0.1 per cent agreement with Planck 2018, obtained from first principles

The headline numerical result of the paper is straightforward and quantitative. Integrating the second Friedmann equation backwards in time from the present epoch down to a scale factor a ≃ 0.01, with the ΛCDM concordance assumption w = −1 and the Planck 2018 cosmological parameters (Ω_m ≃ 0.315, Ω_r ≃ 9.24 × 10⁻⁵, Ω_Λ ≃ 0.685, H₀ ≃ 67.4 km s⁻¹ Mpc⁻¹), the code recovers an age of the Universe of approximately 13.78 Gyr. The Planck Collaboration 2020 value is 13.8 Gyr. The agreement falls within 0.1 per cent.

This result is non-trivial precisely because no free parameters are tuned a posteriori. The concordance arises directly from the physics encoded in the differential system and from the numerical scheme adopted to integrate it.

The methodological choice: a symplectic integrator over cosmological timescales

The paper does not use an explicit Euler scheme, which would accumulate systematic drift over the integration window required to reach a ≃ 0.01. Instead, the authors implement a symplectic Störmer-Verlet integrator (semi-implicit Euler), in which the velocity v = ȧ is updated first using the current acceleration, and the scale factor a is then advanced using the new velocity.

Symplectic schemes are the canonical tool for Hamiltonian systems integrated over long timescales, because they preserve the geometric structure of the underlying phase space. The choice is essential here: the integration spans timescales of order 10¹⁰ years with a time step Δt = 10⁷ yr, and the energy-density terms scale as steep negative powers of the scale factor (radiation as a⁻⁴, matter as a⁻³). Without the symplectic structure, error accumulation would corrupt the late-time behaviour.

The paper substantiates this methodological choice with an explicit convergence study (Figure 3), comparing integrations with Δt = 10⁷, 10⁶ and 10⁵ yr. The maximum relative error of the coarsest run with respect to the finest stays below 0.25 per cent over the entire backward integration interval, decreasing rapidly towards the present epoch. The morphology of a(t) remains visually indistinguishable across the three resolutions.

Seven dark-energy regimes in a unified framework

Beyond the concordance ΛCDM case, the paper performs a systematic scan of the dark-energy equation-of-state parameter w. Seven physically distinct regimes are integrated within the same framework and compared in a single figure (Figure 2) over a temporal window of approximately 30 Gyr.

The phantom regime (w = −1.5) corresponds to a hypothetical dark-energy component whose pressure is more negative than that of a cosmological constant, causing its energy density to grow with the expansion. The acceleration equation enters a runaway regime, and the scale factor diverges in finite cosmic time, the Big Rip scenario originally formalised by Caldwell, Kamionkowski and Weinberg.

The concordance regime (w = −1) reproduces the canonical ΛCDM behaviour, with late-time exponential expansion of the form a(t) ∝ exp(Ht).

The topological-defect regimes (w = −2/3, domain walls; w = −1/3, cosmic strings) produce intermediate acceleration histories between matter dominance and the cosmological constant. The value w = −1/3 sits exactly at the boundary between accelerated and decelerated expansion.

The matter and radiation regimes (w = 0 and w = 1/3 respectively) recover the analytic power-law solutions a(t) ∝ t^(2/3) and a(t) ∝ t^(1/2) within the multi-component model.

A stiff-fluid case (w = 0.6) is included as an illustrative example within the physically admissible range 1/3 < w ≤ 1, bounded above by the causality requirement that the adiabatic sound speed not exceed the speed of light.

Position relative to CLASS, CAMB and the standard tools of computational cosmology

The paper does not aim to compete with established Boltzmann codes such as CLASS or CAMB. Those tools solve the full perturbation system with sub-per-cent accuracy and remain the standard for precision cosmology. Their architecture, however, is highly opaque to anyone outside their development communities: their internal logic is not easily inspectable, and their codebases are not designed for pedagogical re-implementation.

The contribution of this paper is complementary. It produces a self-contained, end-to-end transparent implementation of the background expansion problem, in which every line of code maps unambiguously to a physical or numerical decision. The trade-off is explicit: precision is exchanged for inspectability, and the result is a code base that can serve as a foundation for further extensions, including spatial curvature (Ω_k ≠ 0), dynamical dark energy with time-varying w(t), adaptive time-stepping schemes to resolve the radiation-dominated era at a < 10⁻⁴, and higher-order integrators such as fourth-order Runge-Kutta.

Open access and reproducibility

The paper is published open access under Creative Commons Attribution 4.0, and the full Python source code is deposited on Zenodo with a persistent DOI. The script depends only on the standard scientific Python stack (math, numpy, matplotlib) and runs in a Jupyter Notebook environment. Indexing in NASA/ADS, Scopus and Web of Science is immediate upon publication.

Access to the paper and the code

🔗 Paper on IOPscience: https://iopscience.iop.org/article/10.1088/1361-6404/ae6377 🔗 Code on Zenodo: https://doi.org/10.5281/zenodo.19228187