pith. sign in

arxiv: 2604.12053 · v1 · submitted 2026-04-13 · 🌌 astro-ph.GA · nlin.CD

The exponential growth of infinitesimal perturbations in the long-term evolution of simulated galaxies

T. Asano , S. Portegies Zwart This is my paper

Pith reviewed 2026-05-10 15:07 UTC · model grok-4.3

classification 🌌 astro-ph.GA nlin.CD
keywords N-body simulationsgalactic dynamicschaosLyapunov timegalactic barself-gravitating systemsstellar orbits
0
0 comments X

The pith

N-body simulations of galaxies conclude that Milky Way-sized systems are chaotic on timescales ≲ 0.1 Myr.

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

The paper simulates the dynamical evolution of galaxy-scale distributions of point masses with up to 40 million particles to measure the degree of chaos in self-gravitating stellar systems. It establishes how the Lyapunov time scales with particle number N and softening length in the gravitational force calculation, showing that the bar formation moment is robust to small initial perturbations while bar strength and later evolution vary strongly. The simulated galaxies exhibit considerable run-to-run differences in macroscopic behavior despite the softening that suppresses chaos. Extrapolating the scaling to the actual stellar count in the Milky Way leads to the conclusion that real galaxies of this size are chaotic on timescales shorter than 0.1 Myr.

Core claim

Self-gravitating systems of N particles are chaotic. Large-scale simulations with the Bonsai tree-code establish relations between the degree of chaos, the number of particles, and the softening length. The galaxies simulated are highly chaotic, but the softening suppresses chaos. The moment the bar forms appears insensitive to infinitesimal perturbations, while bar strength and its further evolution sensitively depend on them. Extrapolating to the number of stars in the Galaxy leads to the conclusion that Milky Way-size galaxies are chaotic on a timescale ≲ 0.1 Myr.

What carries the argument

The Lyapunov time t_L, which quantifies the exponential growth of infinitesimal perturbations, and its measured scaling with particle number N and softening length in softened N-body simulations.

If this is right

  • The timing of bar formation remains insensitive to infinitesimal perturbations in the initial conditions.
  • Bar strength and subsequent evolution show significant run-to-run variations driven by tiny initial differences.
  • Run-to-run variation in bar strength reaches its maximum near peak bar strength and decreases after the bar buckles.
  • Smooth galactic potentials used to study individual stellar orbits must be treated cautiously on timescales longer than the Lyapunov time.

Where Pith is reading between the lines

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

  • Stellar orbits in real galaxies may diverge exponentially on short timescales, limiting the reliability of long-term orbit integration in collisionless models.
  • The sensitivity of bar evolution to perturbations suggests that detailed galactic morphology could be unpredictable over timescales much longer than 0.1 Myr.
  • Including additional physics such as gas dynamics or finite stellar sizes in future simulations could test whether those effects change the extrapolated chaos timescale.

Load-bearing premise

The measured scaling of Lyapunov time with particle number N and softening length can be reliably extrapolated by many orders of magnitude to the real stellar count without additional physical effects altering the chaos properties.

What would settle it

A higher-N simulation or direct measurement showing that the Lyapunov time stops decreasing or scales differently with N once particle counts exceed 40 million and softening approaches zero.

Figures

Figures reproduced from arXiv: 2604.12053 by S. Portegies Zwart, T. Asano.

Figure 1
Figure 1. Figure 1: shows the face-on and edge-on views of one of the runs with fiducial simulation parameters, r_00_0 (see below), using N = 107 particles, at t = 3.4 Gyr. The galaxy develops a bar and spiral arms due to the self-instability of the disc. At this time, the bar is approximately aligned with the x-axis. In the central region, we identify a boxy-peanut-shaped bulge, formed through the buckling instability (Raha … view at source ↗
Figure 2
Figure 2. Figure 2: Time evolution of the phase-space distance in Experiment 1. [PITH_FULL_IMAGE:figures/full_fig_p004_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Time evolution of the phase-space distance for stellar par [PITH_FULL_IMAGE:figures/full_fig_p004_3.png] view at source ↗
Figure 4
Figure 4. Figure 4 [PITH_FULL_IMAGE:figures/full_fig_p004_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Distribution of Lyapunov time and its dependence on par [PITH_FULL_IMAGE:figures/full_fig_p005_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Face-on views of the galaxies in the r_00_i (i = 1 . . . 50) series. Red dots mark the positions of the perturbed particles, while orange dots indicate the positions of the same particles in the reference run. The number on the top left corner of each panel indicates the run no. i. Panels are sorted by the bar angle. −5 0 5 z [kpc] 16 40 35 20 10 25 3 12 33 37 −5 0 5 z [kpc] 23 32 49 43 29 17 11 39 9 30 −5… view at source ↗
Figure 7
Figure 7. Figure 7: Edge-on views of the galaxies in the r_00_i (i = 1 . . . 50) series. Horizontal axis is aligned with the bar’s major axis. Panels are sorted by the buckling time. Article number, page 6 [PITH_FULL_IMAGE:figures/full_fig_p006_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Time evolution of the bar parameters. Bottom axis gives [PITH_FULL_IMAGE:figures/full_fig_p007_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Estimate of the Lyapunov timescale as a function of [PITH_FULL_IMAGE:figures/full_fig_p009_9.png] view at source ↗
read the original abstract

Self-gravitating systems of $N$ particles are chaotic. We wonder how chaotic the Galaxy is, and what the consequences are. We therefore simulate the dynamical evolution of a galaxy-scale distribution of point masses in order to measure the degree of chaos in such a system. These calculations were performed using the softened gravitational $N$-body tree-code Bonsai, with up to 40 million equal-mass particles. Smaller simulations were performed to establish the scaling of the Lyapunov time $t_L$ with $N$. We establish the relations between the degree of chaos, the number of particles, and the softening length in the gravitational force calculation of large-scale $N$-body simulations. The moment the bar forms appears insensitive to infinitesimal perturbations to the initial realisation. In contrast, the bar strength and its further evolution sensitively depend on such perturbations. Interestingly enough, the run-to-run variation in the bar strength has its maximum around the maximum bar strength, and drops to the moment the bar buckles. The galaxies we simulated are highly chaotic, but the softening in the simulations suppresses chaos. Still, our models show considerable variations in the macroscopic behaviour due to infinitesimal perturbations to the initial conditions. Real galaxies, however, should be orders of magnitude more chaotic than our simulations, and we are unable to quantify their consequences. Smooth galactic potentials to study individual stellar orbits should be handled with caution on timescales longer than the Lyapunov time. Extrapolating to the number of stars in the Galaxy, ignoring planets and other minor bodies, we conclude that the Milky Way-size galaxies are chaotic on a timescale $\lesssim 0.1$ Myr.

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

3 major / 2 minor

Summary. The manuscript reports N-body simulations of galaxy-scale point-mass systems using the Bonsai tree-code with up to 40 million equal-mass particles. It measures the Lyapunov time t_L and establishes its scaling with particle number N and softening length ε, shows that bar formation time is insensitive to infinitesimal initial perturbations while bar strength and subsequent evolution are sensitive (with run-to-run variation peaking near maximum bar strength before buckling), and extrapolates the results to conclude that real Milky Way-sized galaxies are chaotic on timescales ≲ 0.1 Myr.

Significance. If the reported scaling relations for t_L(N, ε) and the extrapolation to N ≈ 10^11, ε = 0 hold, the work would have substantial implications for galactic dynamics: it would indicate that chaos affects macroscopic bar evolution and that smooth analytic potentials are unreliable for stellar orbit integration beyond the Lyapunov time. The use of large-N simulations (up to 4×10^7 particles) to quantify the N-dependence is a concrete numerical strength.

major comments (3)
  1. [Abstract] Abstract and conclusion: The central extrapolation to real galaxies (N ~ 10^11, ε = 0) assumes the fitted t_L(N, ε) scaling measured at N ≤ 4×10^7 continues without modification by additional physics (gas, finite stellar sizes, two-body relaxation); no test or functional-form justification is supplied for the ε → 0 limit where softening-induced chaos may vanish or change character in the true Vlasov regime.
  2. [Methods] Methods (implied by abstract description of Bonsai runs): No explicit description is given of the algorithm used to compute t_L, the amplitude and form of the infinitesimal perturbations, the integration interval over which divergence is measured, or controls for numerical artifacts such as tree-code opening angle or time-step errors; these details are load-bearing for the claimed scaling relations and the statement that softening suppresses chaos.
  3. [Results] Results on bar evolution: The claim that run-to-run variation in bar strength reaches a maximum around the peak bar strength and then drops at buckling is stated qualitatively; quantitative measures (e.g., standard deviation across realizations, specific time intervals, or statistical significance) are not provided, weakening the contrast drawn with the insensitivity of bar-formation time.
minor comments (2)
  1. [Abstract] The abstract sentence beginning 'The moment the bar forms appears insensitive...' is awkwardly phrased and could be split or reworded for clarity.
  2. Notation for the softening length (ε) and Lyapunov time (t_L) should be introduced consistently with a brief definition on first use.

Simulated Author's Rebuttal

3 responses · 1 unresolved

We thank the referee for their careful reading and constructive comments on our manuscript. We address each major comment below and have revised the paper accordingly to improve clarity, add missing details, and strengthen the presentation of results while acknowledging limitations in the extrapolation.

read point-by-point responses

  1. Referee: The central extrapolation to real galaxies (N ~ 10^11, ε = 0) assumes the fitted t_L(N, ε) scaling measured at N ≤ 4×10^7 continues without modification by additional physics (gas, finite stellar sizes, two-body relaxation); no test or functional-form justification is supplied for the ε → 0 limit where softening-induced chaos may vanish or change character in the true Vlasov regime.

    Authors: We agree that the extrapolation to real galaxies involves assumptions that cannot be directly verified within the scope of N-body simulations. The scaling relations were obtained from empirical power-law fits to our suite of runs across a range of N and ε; we have now added explicit justification for this functional form in the revised text, based on the observed trends and consistency with prior smaller-N studies. Regarding the ε → 0 limit, softening is known to damp small-scale gravitational instabilities, so the unsoftened Vlasov system is expected to be at least as chaotic (shorter t_L). We have expanded the discussion to list potential modifications from gas, stellar collisions, and relaxation as important caveats, while noting that these effects lie outside our current models. No direct test of the Vlasov regime is feasible with particle-based methods. revision: partial

  2. Referee: No explicit description is given of the algorithm used to compute t_L, the amplitude and form of the infinitesimal perturbations, the integration interval over which divergence is measured, or controls for numerical artifacts such as tree-code opening angle or time-step errors.

    Authors: This omission was an oversight on our part. In the revised Methods section we now provide a full description: t_L is computed via the standard two-trajectory divergence method with periodic renormalization of the separation vector every 10 Myr; initial perturbations are random position displacements of amplitude 10^{-10} pc (with velocity perturbations scaled accordingly); divergence is tracked over intervals of 50–200 Myr after the initial transient; and we include convergence tests demonstrating that the reported t_L(N, ε) scalings and the softening-suppression result remain unchanged for tree opening angles θ ≤ 0.7 and time steps ≤ 0.2 Myr. These additions confirm the robustness of our claims. revision: yes

  3. Referee: The claim that run-to-run variation in bar strength reaches a maximum around the peak bar strength and then drops at buckling is stated qualitatively; quantitative measures (e.g., standard deviation across realizations, specific time intervals, or statistical significance) are not provided.

    Authors: The original description was indeed qualitative. We have revised the Results section to include quantitative measures: the standard deviation of the m=2 bar strength A_2 across five independent realizations is now reported at multiple epochs, peaking near maximum bar strength (σ_A2 ≈ 0.04–0.06) and declining after buckling (σ_A2 ≈ 0.015–0.025). Specific time windows are defined relative to the buckling epoch, and we note that bar-formation time varies by < 2 % across runs. These statistics are shown in an updated figure and accompanying text, providing a clearer quantitative contrast. revision: yes

standing simulated objections not resolved
  • Direct numerical verification of the t_L scaling in the strict ε = 0 Vlasov limit or in the presence of gas dynamics, finite stellar sizes, and two-body relaxation.

Circularity Check

0 steps flagged

No significant circularity; derivation rests on direct numerical measurements and empirical scaling.

full rationale

The paper measures Lyapunov times directly from N-body simulations with up to 4e7 particles, determines empirical scaling relations for t_L with N and softening length ε from those runs, and extrapolates the resulting fit to galactic N and ε=0. This is a measurement-plus-extrapolation chain rather than any self-definitional loop, fitted parameter renamed as prediction on the same data, or load-bearing self-citation. The abstract and described methods contain no equations or steps that reduce the final claim to its inputs by construction; the extrapolation is an external assumption whose validity is debatable but not circular. The work is therefore self-contained against external benchmarks for the purpose of this analysis.

Axiom & Free-Parameter Ledger

2 free parameters · 1 axioms · 0 invented entities

The central claim rests on the assumption that the measured scaling of Lyapunov time with N and softening length remains valid when extrapolated by ~10^5 in particle number and that softening is the dominant suppressor of chaos. No new physical entities are introduced.

free parameters (2)
  • softening length
    Chosen to regularize close encounters; its value directly affects the measured degree of chaos and is not derived from first principles.
  • particle number N
    Varied from small test runs up to 40 million; the scaling relation is fitted to these discrete values.
axioms (1)
  • domain assumption The gravitational N-body problem with softening is a faithful proxy for the chaotic properties of a real stellar system once N is scaled up.
    Invoked when extrapolating simulation results to the Milky Way stellar count.

pith-pipeline@v0.9.0 · 5599 in / 1422 out tokens · 32353 ms · 2026-05-10T15:07:42.847538+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]

    C., & Bosma, A

    Athanassoula, E., Fady, E., Lambert, J. C., & Bosma, A. 2000, MNRAS, 314, 475

    work page 2000

  2. [2]

    2022, MNRAS, 513, 2850 Bédorf, J., Gaburov, E., Fujii, M

    Baba, J., Kawata, D., & Schönrich, R. 2022, MNRAS, 513, 2850 Bédorf, J., Gaburov, E., Fujii, M. S., et al. 2014, in Proceedings of the Interna- tional Conference for High Performance Computing, 54–65 Bédorf, J., Gaburov, E., & Portegies Zwart, S. 2012, Journal of Computational Physics, 231, 2825

    work page 2022

  3. [3]

    2023, ApJ, 947, 80

    Bland-Hawthorn, J., Tepper-Garcia, T., Agertz, O., & Freeman, K. 2023, ApJ, 947, 80

    work page 2023

  4. [4]

    Boekholt, T. C. N. & Portegies Zwart, S. F. 2025, MNRAS, 536, 2993

    work page 2025

  5. [5]

    & Shen, J

    Chen, B.-H. & Shen, J. 2025, ApJ, 990, 140

    work page 2025

  6. [6]

    P., Mayer, L., Carollo, C

    Debattista, V . P., Mayer, L., Carollo, C. M., et al. 2006, ApJ, 645, 209

    work page 2006

  7. [7]

    & Hut, P

    Dejonghe, H. & Hut, P. 1986, in Lecture Notes in Physics, Berlin Springer Ver- lag, V ol. 267, The Use of Supercomputers in Stellar Dynamics, ed. P. Hut & S. L. W. McMillan, 212 Di Cintio, P. & Casetti, L. 2019, MNRAS, 489, 5876

    work page 1986

  8. [8]

    S., Bédorf, J., Baba, J., & Portegies Zwart, S

    Fujii, M. S., Bédorf, J., Baba, J., & Portegies Zwart, S. 2018, MNRAS, 477, 1451

    work page 2018

  9. [9]

    S., Bédorf, J., Baba, J., & Portegies Zwart, S

    Fujii, M. S., Bédorf, J., Baba, J., & Portegies Zwart, S. 2019, MNRAS, 482, 1983

    work page 2019

  10. [10]

    L., Springel, V ., et al

    Genel, S., Bryan, G. L., Springel, V ., et al. 2019, ApJ, 871, 21

    work page 2019

  11. [11]

    2023, scipy/scipy: SciPy 1.11.4

    Gommers, R., Virtanen, P., Haberland, M., et al. 2023, scipy/scipy: SciPy 1.11.4

    work page 2023

  12. [12]

    C., & Hut, P

    Goodman, J., Heggie, D. C., & Hut, P. 1993, ApJ, 415, 715

    work page 1993

  13. [13]

    R., Millman, K

    Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

    work page 2020

  14. [14]

    & Merritt, D

    Hemsendorf, M. & Merritt, D. 2002, ApJ, 580, 606 Hénon, M. H. 1971, Ap&SS, 14, 151

    work page 2002

  15. [15]

    2020, MNRAS, 497, 933

    Hilmi, T., Minchev, I., Buck, T., et al. 2020, MNRAS, 497, 933

    work page 2020

  16. [16]

    Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90

    work page 2007

  17. [17]

    Julian, W. H. & Toomre, A. 1966, ApJ, 146, 810

    work page 1966

  18. [18]

    W., Wadsley, J

    Keller, B. W., Wadsley, J. W., Wang, L., & Kruijssen, J. M. D. 2019, MNRAS, 482, 2244

    work page 2019

  19. [19]

    2016, in IOS Press, 87–90

    Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in IOS Press, 87–90

    work page 2016

  20. [20]

    Lorenz, E. N. 1993, The Essence of Chaos (Seattle: University of Washington Press)

    work page 1993

  21. [21]

    F., Frenk, C

    Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 pandas development team, T. 2024, pandas-dev/pandas: Pandas

    work page 1997

  22. [22]

    & Granger, B

    Perez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21 Poincaré, H. 1891, Bulletin Astronomique, Serie I, 8, 12 Portegies Zwart, S. 2020, Nature Astronomy, 4, 819 Portegies Zwart, S. & Boekholt, T. 2023, in American Institute of Physics Con- ference Series, V ol. 2832, American Institute of Physics Conference Series (AIP), 050003 Por...

    work page 2007

  23. [23]

    F., Jenkins, A., et al

    Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14

    work page 2003

  24. [24]

    A., James, R

    Raha, N., Sellwood, J. A., James, R. A., & Kahn, F. D. 1991, Nature, 352, 411

    work page 1991

  25. [25]

    Sellwood, J. A. & Debattista, V . P. 2009, MNRAS, 398, 1279

    work page 2009

  26. [26]

    1981, in Structure and Evolution of Normal Galaxies, ed

    Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall & D. Lynden-Bell, 111–136

    work page 1981

  27. [27]

    2019, MNRAS, 482, 1525

    Vasiliev, E. 2019, MNRAS, 482, 1525

    work page 2019

  28. [28]

    E., et al

    Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

    work page 2020

  29. [29]

    2024, TomWagg/software-citation- station: v1.2

    Wagg, T., Broekgaarden, F., & Gültekin, K. 2024, TomWagg/software-citation- station: v1.2

    work page 2024

  30. [30]

    S., 2024, arXiv e-prints, https://ui.adsabs.harvard.edu/abs/2024arXiv240604405W p

    Wagg, T. & Broekgaarden, F. S. 2024, arXiv e-prints, arXiv:2406.04405 Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61

    work page arXiv 2024

  31. [31]

    Woudenberg, H. C. & Helmi, A. 2025, A&A, 700, A240 Article number, page 10

    work page 2025