## On the Generality and Persistence of Cosmological Stasis

James Halverson and Sneh Pandya

NSF AI Institute for Artificial Intelligence and Fundamental Interactions (IAIFI) Department of Physics, Northeastern University, Boston, MA 02115, USA

Hierarchical decays of N matter species to radiation may balance against Hubble expansion to yield stasis, a new phase of cosmological evolution with constant matter and radiation abundances. We analyze stasis with various machine learning techniques on the full 2 N -dimensional space of decay rates and abundances, which serve as inputs to the system of Boltzmann equations that governs the dynamics. We construct a differentiable Boltzmann solver to maximize the number of stasis e -folds N . High-stasis configurations obtained by gradient ascent motivate log-uniform distributions on rates and abundances to accompany power-law distributions of previous works. We demonstrate that random configurations drawn from these families of distributions regularly exhibit many e -folds of stasis. We additionally use them as priors in a Bayesian analysis conditioned on stasis, using stochastic variational inference with normalizing flows to model the posterior. All three numerical analyses demonstrate the generality of stasis and point to a new model in which the rates and abundances are exponential in the species index. We show that the exponential model solves the exact stasis equations, is an attractor, and satisfies N ∝ N , exhibiting inflation-level e -folding with a relatively low number of species. This is contrasted with the N ∝ log( N ) scaling of power-law models. Finally, we discuss implications for the emergent string conjecture and string axiverse.

{j.halverson,pandya.sne}@northeastern.edu , /github

## I. INTRODUCTION

One of the hallmark characteristics of an expanding universe is the time-evolution of components that contribute to the overall energy density of the universe. It has recently been introduced [1] that this assumption is not always true, and in fact the universe could experience extended phases of cosmological 'stasis' in which the cosmological abundances of matter, radiation, and/or vacuum energy can remain steady over extended e -folds of cosmological evolution, facilitated by various physical mechanisms driving energy pumps which oppose the effects of cosmological expansion [2]. These phenomena can arise naturally from a variety of beyond the Standard Model (BSM) physical theories, for instance in those that predict towers of unstable states that inevitably decay.

rates Γ ℓ of the N particle species that matter dominate prior to stasis. Whereas [1] laid out a general theory of strict stasis and derived many analytic results in an eightdimensional power-law model, including attractor behavior, we focus on understanding stasis in a model-agnostic manner on the full 2 N -dimensional space of rates and abundances. Such an analysis seems to require numerics, for which we employ a variety of machine learning tools (to optimize stasis e -folds, not cosmological viability; see Section VI), but lets us probe the generality of stasis and understand aspects of its duration. The numerical methods uncover a new exponential model of stasis, which leads to extended periods of stasis.

There are a number of different flavors of stasis. In the original formulation [1], stasis is achieved when a tower of ϕ ℓ matter states dominate the energy density but then hierarchically decays into radiation (hereby referred to as M → γ stasis), whereas [2] further introduced mechanisms of vacuum energy to matter (Λ → M ) and vacuum energy to radiation (Λ → γ ) stasis in a similar context, as well as studies of the dynamics of 'triple stasis' with simultaneous Λ → M → γ stasis. Such alternative cosmological histories can have observational consequences [3, 4] depending crucially on the flavor, timescale, and duration of stasis. Though in many cases these provide interesting possibilities for the evolution of our universe, stasis may also constrain ultraviolet-complete theories such as string theory if it arises after nucleosynthesis. We will comment on stasis in the context of Kaluza-Klein towers, string towers, and the string axiverse.

In this paper, we present an analysis of stasis to complement [1], focusing on M → γ stasis. This flavor of stasis depends crucially on the abundances Ω (0) ℓ and decay

This paper is organized as follows. First, we review essentials of M → γ stasis from [1]. In Section II we develop a differentiable Boltzmann solver that facilitates a number of numerical analyses involving stasis. Performing gradient ascent to maximize stasis, we see the emergence of exponential models that motivate log-uniform statistics. In Section III we study random stasis where rates and abundances are drawn from both these loguniform distributions and also power-law distributions. Stasis occurs quite generally, with longer duration in the log-uniform case. In Section IV we use both types of distributions as priors for stasis-conditioned posteriors. The posteriors in this Bayesian analysis are modeled using a neural network known as a normalizing flow, which are optimized using stochastic variational inference. Posterior samples lead to more robust stasis, referring to longer a duration of stasis epochs, and again prefer an exponential model. In Section V we study the exponential model directly, demonstrating that it leads to parametricallyinN longer periods of stasis than power-law models, and discuss potential interfaces with string theory.

## A. Matter-Radiation Stasis

For a tower of states ϕ ℓ where ℓ ∈ [0 , 1 2 , , ..., N ], let ρ ℓ denote a corresponding energy density and Ω ℓ a corresponding abundance. For general types of stasis, ϕ ℓ can either be a tower of massive scalar fields contributing to the total matter abundance Ω M of the universe, or a tower of vacuum energy fields contributing to Ω Λ . Below, we will focus on the M → γ stasis formalism. Recall that for any energy density ρ i , Ω is related via i

$$\Omega _ { i } \equiv \frac { 8 \pi G } { 3 H ^ { 2 } } \rho _ { i } \, \quad \ \ \ ( 1 ) \quad t e.$$

where H is the Hubble parameter and G is Newton's gravitational constant.

Differentiating, the time evolution of Ω i is then

$$\frac { d \Omega _ { i } } { d t } = \frac { 8 \pi G } { 3 } \left ( \frac { 1 } { H ^ { 2 } } \frac { d \rho _ { i } } { d t } - 2 \frac { \rho _ { i } } { H ^ { 3 } } \frac { d H } { d t } \right ) \,. \quad \ ( 2 ) \quad \ \text{duri} \\ \text{abru}.$$

This set of Boltzmann equations is dependent on the time-evolution of individual ρ i as well as H . Using the Friedmann equation for dH/dt in a FriedmannRobertson-Walker (FRW) universe, we obtain

$$\frac { d H } { d t } = - H ^ { 2 } - \frac { 4 \pi G } { 3 } \left ( \sum _ { \ell } \rho _ { \ell } + 3 \sum _ { \ell } p _ { \ell } \right ) \, \quad ( 3 ) \quad \text{Equ} \\ \text{stasi} \quad \text{stasi}$$

where in practice this is expressed in terms of the equation-of-state parameter w ≡ p /ρ ℓ ℓ for a component. Simplifying further, we arrive at

$$\frac { d H } { d t } = - \frac { 1 } { 2 } H ^ { 2 } ( 4 - \Omega _ { M } ) \, \quad \ \ ( 4 ) \quad \text{ for } \\ \quad \text{ the }$$

having invoked w γ = 1 3 and / w M = 0 in equation 3. Integrating both sides, we arrive at

$$H ( t ) = \frac { 2 } { 4 - \langle \Omega _ { M } \rangle } \left ( \frac { 1 } { t - t ^ { ( 0 ) } } \right ) \,, \quad \quad ( 5 )$$

where we have used the approximation that H (0) ( t -t (0) ) ≫ 1 and ⟨ Ω M ⟩ is the time-averaged matter abundance defined as

$$\langle \Omega _ { M } \rangle = \frac { 1 } { t - t ^ { ( 0 ) } } \int _ { t ^ { ( 0 ) } } ^ { t } d t ^ { \prime } \Omega _ { M } ( t ^ { \prime } ) \. \quad \quad ( 6 )$$

During stasis, d ⟨ Ω M ⟩ /dt = 0, in which ⟨ Ω M ⟩ = Ω M , the asymptotic stasis abundance.

In M → γ stasis, the decays of individual matter species are what source the production of radiation in this universe. Stasis epochs necessarily require timedependent energy densities, whose equations of motion in this case are given by

$$\frac { d \rho _ { \ell } } { d t } = - 3 H \rho _ { \ell } - \Gamma _ { \ell } \rho _ { \ell } \. \quad \quad \quad ( 7 )$$

Returning to equation 2 and substituting in the result of equation 4 and 7, we arrive at the set of N ODEs which govern the time evolution of individual Ω ℓ

$$\frac { d \Omega _ { \ell } } { d t } = H \Omega _ { \ell } \left ( 1 - \Omega _ { M } \right ) - \Gamma _ { \ell } \Omega _ { \ell } \,. \quad \quad \left ( 8 \right )$$

Equation 8, in combination with the ODE for the Hubble parameter directly gives the dynamics for our system. As our Universe only contains matter and radiation, Ω M + Ω γ = 1 at all times and the dynamics for radiation are easily obtained by recognizing that d Ω M /dt = -d Ω γ /dt .

For M → γ stasis, the universe is beginning in a matter dominated state, Ω M ( t (0) ) = 1. With time-evolution, the individual matter species gradually redshift or begin decaying into radiation; the effects of both must balance to have d Ω M /dt = 0 during stasis. This will happen for all times when the decays of ϕ ℓ are exactly counterbalanced by Hubble expansion. Further, all individual Ω 's ℓ during stasis must cooperate to produce an asymptotic abundance Ω M . These combined form the two necessary and sufficient conditions for an extended period of stasis:

$$\sum _ { \ell } \Gamma _ { \ell } \Omega _ { \ell } = H ( \Omega _ { M } - \Omega _ { M } ^ { 2 } ) \quad \quad ( 9 )$$

$$\sum _ { \ell } \Omega _ { \ell } ( t ) = \overline { \Omega } _ { M } \,.$$

Equation 9 as written is actually a condition for eternal stasis, which of course cannot be physical for some finite tower of states. However, we can illuminate how the stasis epoch ends by operating under the assumption of eternal stasis. Let us assume we're in a period of stasis where Ω M has achieved its Ω M stasis value. With these conditions, we can study the explicit time dependence for many of the quantities of interest. Beginning with the Hubble parameter in equation 3, the solution is

$$H ( t ) = \left ( \frac { 2 } { 4 - \overline { \Omega } _ { M } } \right ) \frac { 1 } { t } \, \quad \quad \ \ ( 1 1 )$$

which further implies that the scale factor grows as

$$a ( t ) = a _ { * } \left ( \frac { t } { t _ { * } } \right ) ^ { 2 / ( 4 - \overline { \Omega } _ { M } ) } \quad \text{(1)}$$

for a fiducial time t ∗ . It further follows from equation 7 that

$$\rho _ { \ell } ( t ) = \rho _ { \ell } ^ { * } \left ( \frac { t } { t _ { * } } \right ) ^ { - 6 / ( 4 - \overline { \Omega } _ { M } ) } e ^ { - \Gamma _ { \ell } ( t - t _ { * } ) } \,, \quad ( 1 3 )$$

which in turn implies that

$$\Omega _ { \ell } ( t ) = \Omega _ { \ell } ^ { * } \left ( \frac { t } { t _ { * } } \right ) ^ { 2 - 6 / ( 4 - \overline { \Omega } _ { M } ) } e ^ { - \Gamma _ { \ell } ( t - t _ { * } ) } \,. \quad ( 1 4 )$$

The result of equation 11 when inserted into equation 9 while assuming a period of stasis gives

$$\sum _ { \ell } \Gamma _ { \ell } \Omega _ { \ell } = \frac { 2 \overline { \Omega } _ { M } ( 1 - \overline { \Omega } _ { M } ) } { 4 - \overline { \Omega } _ { M } } \frac { 1 } { t } \, \quad \quad ( 1 5 )$$

exhibiting a power-law dependence for t , which cannot be true for all t . Thus, this will not yield an eternal stasis epoch, but a stasis epoch which is terminated when all species decays have concluded.

The models studied in [1, 2] consider a spectrum of decay widths { Γ ℓ } and abundances { Ω ℓ } motivated by a variety of BSM models which follow a power-law scaling

$$\Gamma _ { \ell } = \Gamma _ { 0 } \left ( \frac { m _ { \ell } } { m _ { 0 } } \right ) ^ { \gamma }, \ \Omega _ { \ell } ^ { ( 0 ) } = \Omega _ { 0 } ^ { ( 0 ) } \left ( \frac { m _ { \ell } } { m _ { 0 } } \right ) ^ { \alpha } \quad ( 1 6 ) \quad \text{$\frac{R\text{urt} } { \text{dial} 1$} \\ \text{$\frac{L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{\tt{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\tt{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$S\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$S\tt{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$F\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$F\tt{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$C\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$C\tt$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text	$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text`$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text
$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text<$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text'$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text}$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text.$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text1$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\text{$L\$$$

where the mass spectrum takes the form

$$m _ { \ell } = m _ { 0 } + ( \Delta m ) \ell ^ { \delta } \, \quad \ \ ( 1 7 ) \, \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

and Ω (0) 0 is a normalization factor enforcing Ω M ( t (0) ) = 1. The parameters α γ , , and δ are further restricted to the following range:

$$- \frac { 1 } { \delta } < \alpha \leq \frac { \gamma } { 2 } - \frac { 1 } { \delta } \. \quad \quad \ \ ( 1 8 ) \quad \ g r _ { i } ^ { \ v i c } \\ \text{sol}$$

This therefore defines a 8 D model parameterized by

$$\{ \alpha, \gamma, \delta, m _ { 0 }, \Delta m, \Gamma _ { 0 }, \Omega _ { 0 } ^ { ( 0 ) }, t ^ { ( 0 ) } \} \quad \ \ ( 1 9 ) \quad \text{ise} \, \cdot$$

which is crucially a subset of the full 2 N -dimensional input parameter space that we seek to study.

It is important to note that missing from this list of parameters if the initial value of the Hubble constant, H (0) . This model of stasis has exhibited global attractor properties that were extensively studied in [1]. The initial timescale for ϕ ℓ decays is dictated by the ratio Γ N -1 /H (0) . When this ratio is small, the decays begin slowly after the starting time t (0) , and the socalled 'edge effects' in [1] are mild. Conversely, when Γ N -1 /H (0) ≫ 1, particle decays begin almost immediately and there are severe edge effects. These edge effects are indeed necessary for a stasis state to both begin and end, as the condition in equation 15 cannot be true when the decay process is just beginning or has concluded. However, due to the global attractor nature of stasis it is possible to achieve the same configuration of stasis, with the exception that the stasis state is approached from below rather than above.

## II. MAXIMIZING STASIS WITH DIFFERENTIABLE SIMULATIONS

The duration of stasis may be determined by solving the N +1 Boltzmann equations on a (2 N +1)-dimensional parameter space of decay rates, abundances, and the initial value of the Hubble parameter. This problem is difficult due to its high dimensional nature, and in general we would also like to be able to differentiate the numerical solution to the Boltzmann equations to understand how the duration of stasis responds to variations in the rates and abundances. A numerical Boltzmann solver may be thought of as a type of simulator, and we seek to differentiate through the entire simulation process.

Differentiable simulators are part of a growing trend in ML applications within the sciences, motivated by the emergence of more powerful and robust simulations. They have accelerated scientific analyses from molecular dynamics [5] and biology to physics and cosmology [6]. Further, with the advent of neural-network based probabilistic modeling, differentiable simulations are essential for implementing techniques such as stochastic variational inference [7] and Hamiltonian Monte Carlo [8]. Differentiable simulations thus serve as powerful scientific tools for which to do both purely data-driven and probabilistic modeling studies.

For N species, our task is to solve a set of N +1 coupled ODEs as given in equations 4 and 8 and compute the number of e -folds of a potential stasis epoch. We utilize diffrax [9], a jax -based [10] library that provides numerical differential equation solvers that preserve gradients and also allow backpropagation through our solutions to the Boltzmann equations. jax is a differentiable high-performance numerical computing library which utilizes just-in-time ( JIT ) compilation and vectorized computations, earning popularity in the sciences for its efficiency and modularity. The set of ODEs is stiff, meaning certain numerical techniques will be unstable without a sufficiently small resolution for the step-size. We use the Kvaerno5 [11] solver which is a 5th order explicit singly diagonal implicit Runge-Kutta (ESDIRK) method suited for stiff ODEs, with an absolute and relative tolerance atol = rtol = 10 -8 in the step size. We terminate the solver at t = t max when Ω M ( ) t = 10 -4 , indicating that the universe has passed into a radiationdominated epoch. Backpropagation through a differential equation can also be very memory expensive. We use RecursiveCheckpointAdjoint [12, 13] which utilizes a binomial checkpointing scheme to preserve memory usage during backpropagation.

We seek to utilize the gradients in our solver to guide the parameter space towards stasis configurations. It is therefore not enough to solve the Boltzmann equations and be able to backpropagate through them, we also define a differentiable algorithm to compute the stasis duration and asymptotic abundance from a given Ω M curve. To this extent, we must first define a numerical notion of numerical stasis that can be applied to numerical Ω M curves.

## A. ϵ -Stasis

Stasis as defined up to this point is a strict condition on the time-evolution of a cosmological component. In [1], the model introduced is constructed so that this can be achieved exactly. However, from a numerical perspective, configurations that yield small deviations from exact stasis can nevertheless yield significant alterations to a cosmology. A more general definition of stasis is an

epoch in which Ω M is 'flat' enough for 'long' enough. This requires a parametric notion of stasis, and we therefore introduce an ϵ -tolerance on stasis according to two different definitions. First, we develop a notion of stasis that allows us to create a differentiable (to enable the differentiable simulator) and accurate stasis finder algorithm to isolate epochs of stasis. Then we present a more intuitive notion of stasis that is utilized in all the presentations of our results.

We begin with the notion of stasis that admits a differentiable stasis finder algorithm. Recall that stasis is a phenomenon induced by the cooperative behavior of individual Ω . ℓ A period of stasis is therefore computed by analyzing the total abundance Ω M ( ) t = ∑ ℓ Ω ( ). ℓ t For a given Ω M ( ) t curve, we must compute the asymptote around which stasis occurs, and further isolate the duration of a stasis epoch. To do this in a differentiable way, we introduce an exponential weighting to compute a 'flatness score' that yields the stasis duration in t

$$\Lambda \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \dots \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \delta \, \dagma \, \dagma \, \dagma \, \dagma \, \dagma \, \dagma \, \dagma \, \dagma \, \dagma \, \dagma \, \d$$

where Θ is the Heaviside function, enforcing that only mixed-component cosmologies are considered, and the index i is that of the i th time-step in the solution Ω M ( ). t A tolerance σ = 0 02 and . δ = 0 09, which roughly corre-. sponds to a window-tolerance of ± 0 1 about the asymp-. totic stasis abundance, was found to work best for optimization.

Upon termination of the solver, we use the solution for H t ( ) to calculate the total number of e -folds of the simulation via

$$\mathcal { N } _ { m a x } = \int _ { t ^ { ( 0 ) } } ^ { t _ { \max } } H ( t ) \, d t \, \quad \ \ ( 2 1 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where it is assumed that a t ( (0) ) = 1. This integral is evaluated numerically using the composite trapezoidal rule (i.e. jax.numpy.trapezoid ). It is then straightforward to compute the number of e -folds of stasis N via normalization with respect to the total number of e -folds,

$$\mathcal { N } = \mathcal { N } _ { m a x } \cdot \frac { t _ { s } } { t _ { m a x } } \. \quad \quad ( 2 2 ) \quad \text{pa} \\ \text{un}$$

This normalization step allows us to isolate the contribution of stasis epoch relative to the entire simulation duration.

Intuitively, the stasis finder measures flatness by considering both how similar consecutive in Ω M,i values are and how close these values are to Ω M . It does so by creating a score that rewards sequences where values are close to each other (indicating flatness) and close to the target value (indicating relevance), adjusted by the parameters ϵ and δ to fine-tune sensitivity and scaling.

This differentiable calculation, however, requires that Ω M is known. This is retrieved from a separate, differentiable abundance-finder algorithm. The algorithm first marks Ω M values that are not matter dominated or radiation dominated (i.e. 0 01 . &lt; Ω M &lt; 0 99). . It then determines links between consecutive valid values that are within a specified tolerance ϵ . It proceeds to count these links for each time step and identifies the range of indices with the maximum number of links, indicating a period of stasis. Finally, it filters the valid Ω M values within this range and computes Ω M as the median of these values, representing the typical abundance during the stasis period. This algorithm, of course, also yields a prediction for N , but we find by examining solutions that it is inaccurate and does not function as well as equation 20.

A pitfall of the differentiable stasis-finder in equation 20 is that the sum of flatness scores considers all parts of Ω M ( ) that are near the asymptotic abundance, as opt posed to those only within the stasis period. As such, the differentiable finder can be biased in e -folds; however, we will see that it still properly serves the purpose of a guiding the parameter space towards stasis configurations.

For these reasons, a more accurate, non-differentiable stasis finder which uses a sliding-window algorithm to isolate the longest stasis period is used in all presented results. This sliding window stasis finder is configured for a 10%-tolerance, e.g. the number of e-folds of stasis is determined by considering a window of ± 0 1 around . Ω M as valid. All the results presented in this paper use this notion of stasis.

## B. Maximizing Stasis

We have amassed the necessary ingredients to use the differentiable simulation to optimize for stasis. Using the simulation S ( θ = Γ { ℓ , Ω (0) ℓ , H (0) } ), which returns N and Ω M after solving the Boltzmann equations and its gradients ∇ Γ ℓ and ∇ Ω (0) ℓ , we can employ gradient ascent to optimize on stasis e -folds.

We begin with a single vector of samples of Ω (0) ℓ and Γ ℓ for a given N , initialized according to draws from any appropriate distribution. The Boltzmann equation for Ω ℓ is inherently dimensionless; the only dimensionful parameters are H (0) and Γ . We will take all Γ ℓ ℓ to be in units of Planck mass M p and no greater than this scale, e.g., max(Γ ) ℓ ≤ 1 M p . Similarly, we will restrict decay rates to be no smaller than those corresponding to the current age of the universe, min(Γ ) ℓ ≥ 10 -62 M p .

We consider both power-law and standard uniform distributions for initializing Γ ℓ and Ω (0) ℓ . The power-law distribution is representative of the original model of stasis in the limit that ∆ m/m 0 ≫ 1, while the standard uniform respects minimal physical constraints (i.e. positivedefiniteness and max(Γ ) ℓ ≤ 1 M p ) on the parameters. Recall that Ω (0) ℓ are normalized upon entering the simu-

lation such that Ω M ( t (0) ) = 1. To generate samples from a power-law distribution with a relative scaling ℓ β , we draw from a Pareto distribution with a shape parameter α p = 1 /β . These samples are then inverted to represent draws from a power-law distribution. Henceforth, we denote power-law distribution samples as X ∼ ℓ β , with the understanding that they are generated as the inverse of samples X ∼ Pareto( α p = 1 /β ).

After sampling but before entering the simulation, we sort the Γ , which is without loss of generality since the ℓ ℓ subscript amounts to a species definition. However, we also often sort the Ω (0) ℓ , so that increasingly large rates correspond to increasingly large abundances. This introduces a non-trivial physics-motivated correlation between the parameters, and we henceforth refer to such samples as 'sort-correlated.' When we speak of identically and independently draw (i.i.d.) parameters, we mean prior to sort-correlation, unless otherwise stated.

Sort-correlation may be performed with any off-theshelf sorting algorithm, including but not limited to a simple jax.numpy.sort . As Γ ℓ and Ω (0) ℓ are updated with gradient ascent, we would like to continue to enforce that their spectra are sort-correlated. Doing so crucially requires that they are sorted smoothly and differentiably, to not interrupt the flow of gradients during optimization. For this reason, we implement a custom differentiable bitonic sorting algorithm inspired by [14]. This algorithm works by recursively dividing an array into smaller sub-arrays, sorting them, and then merging them using compare-and-swap operations, the latter of which is modified to be differentiable. This sorting technique is also used when doing SVI experiments.

One necessary constraint to recognize is that simply optimizing on stasis e -folds will encourage matterdominated cosmologies, as species' time spent decaying detracts from time spent redshifting, which is what contributes to the overall stasis duration. As such, a stasis optimization condition f is designed to enforce a mixed component cosmology

$$f ( \theta ) = & \, \mathcal { N } ( \theta ) - \alpha \left [ ( \overline { \Omega } _ { M } ( \theta ) - l ) ^ { 2 } + ( \overline { \Omega } _ { M } ( \theta ) - u ) ^ { 2 } \right ] \ \ ( 2 3 ) \quad \ \ \ 9 \colon \\ & \times \left ( 1 - \Theta ( \overline { \Omega } _ { M } ( \theta ) - l ) \Theta ( u - \overline { \Omega } _ { M } ( \theta ) ) \right ) \, \quad \ \ \ 9 \colon$$

where we recognize that N and Ω M are outputs of S ( θ ). Above, l = 0 2 and . u = 0 8, . defining bounds on the allowed matter abundance, and α specifies the regularization strength depending on the experiment, ensuring the constraint has a meaningful effect on the optimization landscape. The Heaviside function indicates that the penalty is only applied when Ω M is outside the allowed abundance window. It is now clear to see why the ϵ -stasis finder algorithm needs to be differentiable in both e -folds and matter abundances, as the gradients of Ω M are necessary to study optimization of stasis.

When performing gradient ascent on stasis, the sequen- tial gradient updates for Γ ℓ and Ω (0) ℓ are computed as

$$\Omega _ { \ell, i + 1 } ^ { ( 0 ) } = \Omega _ { \ell, i } ^ { ( 0 ) } + \eta ( t ) \nabla _ { \Omega _ { \ell, i } ^ { ( 0 ) } } f ( \theta ) \text{ \quad \ \ } ( 2 4 )$$

$$\Gamma _ { \ell, i + 1 } = \Gamma _ { \ell, i } + \eta ( t ) \nabla _ { \Gamma _ { \ell, i } } f ( \theta ) \text{ \quad \ \ } ( 2 5 )$$

for a given step i , where η t ( ) is a (potentially) timedependent learning rate. In subsequent experiments, a decay-factor γ is applied at epoch t ′ such that the learning rate has the functional form

$$\eta ( t ) = \begin{cases} \eta _ { 0 }, & \text{for $t<t^{\prime}$} \\ \gamma \cdot \eta _ { 0 }, & \text{for $t\geq t^{\prime}$} \end{cases}$$

Optimization is subject to an early-stopping criterion if f ( θ ) does not improve over a specified number of epochs ξ or if a NaN is encountered during optimization.

A detailed algorithm of optimizing stasis with gradient ascent is shown in Algorithm 1. In short, the algorithm 1) generates initial samples for rates and abundances; 2) sort-correlates them; 3) solves the Boltzmann equations; 4) uses the differentiable stasis-finder to compute the number of stasis e-folds; 5) compute gradients through the solution; 6) updates the parameters according to gradient ascent on stasis. This yields the parameters for the next iteration of the pipeline.

## Algorithm 1: Gradient Ascent on Stasis

Require: θ = { Γ ℓ , Ω (0) ℓ } simulation parameters, N stasis e -folds, α penalty coefficient, l lower bound on matter abundance, u upper bound on matter abundance, η t ( ) learning rate, ξ early-stopping threshold

- 1: Initialize: Γ , Ω ℓ (0) ℓ via sampling and sort-correlate their spectra
- 2: while not converged do
- 3: N , Ω M ←S ( θ )
- 4: f ( θ ) ← N ( θ ) -α [ (Ω M ( θ ) -l ) 2 +(Ω M ( θ ) -u ) 2 ] × ( 1 -Θ(Ω M ( θ ) -l )Θ( u -Ω M ( θ )) )
- 5: ∇ Ω (0) ℓ f ( θ , ) ∇ Γ ℓ f ( θ ) ▷ Compute gradients
- 6: Ω (0) ℓ,i +1 ← Ω (0) ℓ,i + ( ) η t ∇ Ω (0) ℓ,i f ( θ ) ▷ Update abundances
- 7: Γ ℓ,i +1 ← Γ ℓ,i + ( ) η t ∇ Γ ℓ,i f ( θ ) ▷ Update decay rates
- 8: Ω (0) 1 ≤ Ω (0) 2 ≤ . . . ≤ Ω (0) N ▷ Differentiably sort Ω (0) ℓ

9:

Γ

1

≤

Γ

2

≤

. . .

≤

Γ

N

▷

Differentiably sort Γ

ℓ

- 10: Γ ℓ ← clip(Γ ℓ , 0 1) ,

▷

Ensure Γ

ℓ

physical

- 11: if f ( θ ) does not improve over ξ epochs or NaN is encountered then
- 12: break

▷ Early-stopping criterion

- 13: end if
- 14: end while

We study the outcome of gradient ascent optimization for initializations from a power law distribution corresponding to Γ ℓ ∼ ℓ 3 and Ω (0) ℓ ∼ ℓ 1 and a uniform initialization where Γ ℓ , Ω (0) ℓ ∼ Uniform(0 1). , We conduct the experiments for N = 50 species and with Γ N -1 /H (0) = 0 1. . We optimize both for a total of 50000 epochs with an early stopping threshold of ξ = 2000 epochs. An initial learning rate of η = 0 01 is used. .

FIG. 1. Gradient ascent trajectories for uniform initialization with N = 50 species and Γ N -1 /H (0) = 0 1. . Parameters were optimized subject to the condition in equation 23 for 50,000 epochs of optimization with α = 10. Intermediate trajectories for this initialization are relatively noisy, as seen by the red lines.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-a1709ab5-0b02-4218-9a16-076fa120df5b.jpg)

Example gradient ascent trajectories for the uniform initialization is shown in Figure 1, where the optimization has settled on Ω M = 0 2 and has achieved 27 . e -folds of stasis. The change from initialized to optimized parameters is seen in going from the black-dashed line at initialization, through increasingly dark red intermediate trajectories as gradient ascent progresses, converging to the solid black line. From the trajectories, it can be seen that the optimization was relatively noisy. Indeed, from a numerical perspective what is essential for a robust epoch of stasis is large hierarchies in Γ ℓ and smaller hierarchies in Ω (0) ℓ , as is manifest for the power-law model of stasis and not characteristic of draws from a uniform distribution. Conceptually, this means that large changes are required to achieve stasis.

When initializing parameters as a power-law, the benefit of immediate hierarchies in the species spectra is apparent in Figure 2. The optimization is more stable, as shown by the well-behaved intermediate values in red, and results in a longer epoch of stasis lasting 34 e -folds. It is then interesting to wonder: does a power-law initialization still result in a power-law model after following gradients, or does it change qualitatively? Similarly, does the uniform initialization change qualitatively under gradient ascent?

We can further deploy the differentiable simulation to study the model dependence of the optimized Γ ℓ and Ω (0) ℓ . We optimize 50 random initializations of θ = { Γ ℓ , Ω (0) ℓ } with the differentiable simulation initialized according to a uniform distribution, Γ ℓ , Ω (0) ℓ ∼ U (0 , 1),

FIG. 2. Gradient ascent trajectories for power-law initialization with N = 50 and Γ N -1 /H (0) = 0 1. . The benefit of decompressed Γ ℓ and Ω (0) ℓ spectra associated to power-law initialization is seen by the less noisy (red) background trajectories in the right panel and the more robust epoch of stasis when compared to the uniform initialization.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-7be1855f-8d5e-4727-af21-86ceaf799db4.jpg)

and power-law distributions, Γ ℓ ∼ ℓ 3 and Ω (0) ℓ ∼ ℓ 1 . We conduct our experiments for N = 50 and Γ N -1 /H (0) = 0 1, the results of which are shown in Figure 3. .

We can quantify the model dependence of Γ ℓ and Ω (0) ℓ to a power-law model with linear fit in log-log space when looking at the functional dependence of Γ ℓ and Ω (0) ℓ with ℓ . Similarly, a linear fit in semi-log space is indicative of an exponential dependence with ℓ . We use the coefficient of determination ( R 2 score) for comparisons of model dependence across distributions and across scales, as it enjoys a scale invariance while still encoding information of fit residuals. The scale invariance is crucial, as the relative scale of optimized Γ ℓ and Ω (0) ℓ are a priori different than from their respective initializations. Traditional metrics such as mean squared error are sensitive to this, and can yield misleading results.

In Figure 3 we see that at initialization there is a clear preference towards a power-law model over exponential, as both Ω (0) ℓ and Γ ℓ generally exhibit R 2 &gt; 0 9, . while exponential fits feature R &gt; 2 0 6. . This is expected, since neither of the initialization distributions were exponential. After optimization, there is a clear drift towards perfect exponential fit ( R 2 = 1) across all experiments, indicating a clear preference towards an exponential model of stasis when optimizing on stasis e -folds, even when being initialized with a power-law distribution that is reflective of the existing model of stasis. It is also seen in Figure 3 that the optimized values are much more abundant in stasis e -folds than their respective initializations.

FIG. 3. Experiments showing a preference for the exponential model upon optimizing stasis with gradients for uniform and power-law initializations for N = 50 and Γ N -1 /H (0) = 0 1. . Gradient ascent was subject to the constraint 0 2 . &lt; Ω M &lt; 0 8. . Optimization was done for 50 random initializations for 50,000 epochs with early-stopping. An initial learning rate η 0 = 0 01 . was used with a γ = 0 1 multiplicative decay at epoch . t ′ = 10 000. , A clear bias towards an exponential model is shown for optimal Γ ℓ and Ω (0) ℓ , even when initialized with a power-law distribution similar to the original model of stasis. It is also seen that the drift shifts across the 1:1 line dividing exponential and power-law confidence equality. In some instances, the initialized parameters are shown to already be a good exponential fit due to the compression of the relative abundance and decay spectra. Even under such conditions, the drift towards a more exponential model and away from a power-law model is evident.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-09a7d329-b0d1-4dd7-a703-d19f35a3b783.jpg)

This result suggests a new model of stasis for which

$$\Omega _ { \ell } ^ { ( 0 ) } \in e ^ { \alpha \ell } \quad \Gamma _ { \ell } \in e ^ { \gamma \ell }, \quad \quad ( 2 7 ) \quad \Re$$

As discussed, we may sort the Γ ℓ without loss of generality, but in this section we keep Ω (0) ℓ unsorted.

an exponential model which is qualitatively different than the power-law model introduced in [1]. It also opens the door for additional physical mechanisms that can result in particle spectra that can induce a stasis state. We will comment on those physics models in Section V.

We emphasize that this result has very little model bias: no trained neural networks or strong prior beliefs that restrict the effective parameter dimension were used in arriving at this result. The only assumption is a prior on the full parameter space from which the initial parameter are drawn, and then we simply differentiably simulate the stasis phenomenon conditioned on increasing stasis. On statistical grounds, these results suggest a new distribution that these parameters follow; that is, a log-uniform distribution, a simple distribution for which samples are uniformly distributed in orders of magnitude. We will further study this model of stasis, and more generally the stasis parameter space from a distributional point of view, using neural networks in a Bayesian inference setting. To that end, these results suggest a new prior distribution.

## C. Stasis and Unsorted Abundances

Wehave so far exclusively operated under the setting of sort-correlated Γ ℓ and Ω (0) ℓ . It is nonetheless interesting to study the emergence of a stasis epoch without respecting strict physical correlations between the parameters.

We proceed by initializing the decay rates as Γ ℓ ∼ Log-Uniform(10 -62 , 10 ) and sorting the spectrum, and 0 Ω (0) ℓ ∼ Log-Uniform(10 -2 , 10 ) 0 for N = 50 species, evolving the parameters according to the gradient update rules given in equations 24 for 50,000 epochs with early-stopping.

The result of this experiment is shown in Figure 4. This analysis has indeed shown that robust stasis epochs are possible without strictly correlated decay rates and abundances , even though the species contributing to stasis have abundances increasing with decay rates. This realization would have been difficult to make with analytical modeling!

We see in Figure 4 that the randomly distributed abundances upon initialization become roughly monotonic up to ℓ ∼ 30, which are exactly the species in the left panel that balance to produce the stasis epoch lasting 43 e -folds. With the spectrum viewed on a semi-log plot in the right panel of Figure 4, it is clear to see that the contributing species follow an exponential model. The species ℓ &gt; 30 are those that decay at early times, setting the Ω M value. By just following gradients, it was deduced that these species did not in fact need to respect strict correlations with Γ ℓ to produce stasis.

## III. RANDOM STASIS IN PHYSICS-MOTIVATED DISTRIBUTIONS

The results of simply following gradients in the differentiable simulation motivates an exponential model,

FIG. 4. Example stasis epoch and gradient ascent trajectories for unsorted abundances with N = 50 and Γ N -1 /H (0) = 0 1. The . initializations chosen were Γ ℓ ∼ Log-U(10 -62 , 10 ) and Ω 0 (0) ℓ ∼ Log-U(10 -2 , 10 ). 0 Γ ℓ samples were sorted before optimization via gradient ascent, which is without loss of generality as ℓ is a species definition. Ω (0) ℓ are left unsorted. We see that the gradients learn to correlate lower ℓ species which contribute to the overall stasis duration, resulting in a period of stasis lasting ∼ 43 e -folds at an abundance of Ω M = 0 80. We see in the right panel that these contributing species approximately follow an . exponential, characterized by a linear dependence on ℓ on a semi-log plot. For high ℓ species, which decay at early times to set the stasis abundance Ω M , the algorithm does not learn to sort them. Indeed, ∼ 40% of abundances do not monotonically increase with ℓ .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-e1f9b9a4-1b49-4cf1-82af-2c35b6349a41.jpg)

as opposed to the original power-law model of [1]. In a statistical language, this motivates the analysis of log-uniform priors in addition to the power-law priors. Power-law models were argued to be physically relevant in [1], and we will argue in Section V that exponential models are also physically motivated.

In this section, we study stasis in random draws from these physically-motivated priors, and argue that they are also statistically well-motivated in Bayesian statistics. First we will provide a statistical perspective on the priors, and then study stasis.

## A. Scale-Invariant Priors and Decay Rates

A prior on a vector of random variables θ is a density that represents some knowledge of the system or beliefs about the way it behaves, before further evidence is introduced. There exists a particular class of priors known as 'uninformative priors,' and a subclass of such priors that are scale-agnostic.

One of these scale-invariant priors is known as the Jeffreys prior. Mathematically, it is defined to exhibit the scaling

$$p ( \theta ) \varpropto \sqrt { \det \mathcal { I } ( \theta ) }$$

for its density function, where I ( θ ) is the Fisher infor- mation metric defined as

$$\underset { \text{vant} { i o r } } { s } \text{ of } \quad \mathcal { I } ( \theta ) _ { i j } = \mathbb { E } \left [ \left ( \frac { \partial } { \partial \theta _ { i } } \log f ( X | \theta ) \right ) \left ( \frac { \partial } { \partial \theta _ { j } } \log f ( X | \theta ) \right ) \left | \theta \right ] \\ \text{ntial } \quad \text{ for a union libalibroad function of non-annotore } f / \text{$Y$} / \text{$\lambda$} \right \rangle$$

for a given likelihood function of parameters f ( X θ | ), which is the probability density of a set of observed data X given parameters θ . We see that the Jeffreys prior is the volume measure with respect to the Fisher information metric. It is therefore diffeomorphism invariant and has the significant advantage that it does not depend on the choice of coordinates for model parameters.

We will now proceed to derive an appropriate Jeffreys prior for Γ . ℓ The time until decay τ for a particle given a decay rate Γ is governed by the PDF

$$f ( \tau | \Gamma ) = \Gamma \, e ^ { - \Gamma \tau } \,. \quad \quad \quad ( 3 0 )$$

Under conditions in which a likelihood function is twice differentiable and has a vanishing expectation value of the gradient of the log-likelihood, a simpler definition of the Fisher information can be used. By inspection, we see that f ( τ | Γ) is twice differentiable with respect to Γ. It only remains to check that

$$\mathbb { E } \left [ \frac { \partial } { \partial \Gamma } \log f ( \tau | \Gamma ) \right ] & = \mathbb { E } \left [ \frac { 1 } { \Gamma } - \tau \right ] \quad \quad ( 3 1 ) \\ & = \frac { 1 } { \Gamma } - \frac { 1 } { \Gamma } \\ & = 0 \,,$$

FIG. 5. ( Left ) Stasis configurations for N = 100 species initialized with power-law prior draws across 100 realizations for Ω (0) ℓ and Γ , which correspond to the model of stasis introduced in [1]. ℓ Ω (0) ℓ and Γ ℓ are additionally sort-correlated before entering the simulation. The (non-differentiable) sliding-window stasis finder was used with a 10%-tolerance. We see that the maximum number of e -folds for this prior distribution is N ∼ 8 4 . e -folds for a more matter-dominated cosmology. ( Right ) Stasis configurations for N = 100 species initialized with exponential prior draws across 10 realizations, which correspond to sorted samples from a log-uniform distribution. The axes values correspond to Ω (0) ℓ ∼ Log-U(10 α , 10 ) and Γ 0 ℓ ∼ Log-U(10 γ , 10 ), with 0 the values chosen to illustrate the transition in N and Ω M . The maximum number of e -folds for this prior is N ∼ 55 1 . e -folds, a noticeably longer stasis duration than the power-law distribution and completely matter-dominated. Both distributions also feature a disallowed region in Ω (0) ℓ , in which the abundance spectrum becomes sufficiently stressed that stasis is not possible.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-b69d9308-97b1-481a-a48f-78618f1a7c13.jpg)

where we have invoked that E [ τ ] = 1 Γ, as / τ follows an exponential distribution with respect to Γ. With these conditions satisfied, we can use a simpler form of the Fisher information

$$\mathcal { I } ( \theta ) _ { i j } = - \mathbb { E } \left [ \frac { \partial ^ { 2 } } { \partial \theta _ { i } \partial \theta _ { j } } \log f ( X ; \theta ) \Big | \theta \right ] \quad ( 3 2 )$$

in deriving the Jeffreys prior for Γ . ℓ

The likelihood of observing a set of particles for a period of time τ = { τ 1 , τ 2 , ..., τ N } before they decay with i.i.d decay rates Γ = Γ { 1 , Γ 2 , ..., Γ N } is given by

$$p ( \tau | \mathbf \Gamma ) = \prod _ { i = 1 } ^ { N } f ( \tau _ { i } | \Gamma _ { i } ) = \prod _ { i = 1 } ^ { N } \Gamma _ { i } \, e ^ { - \Gamma _ { i } \tau _ { i } } \,. \quad ( 3 3 )$$

Exploiting properties of logarithms, we arrive at

$$\log p ( \tau | \mathbf P ) = \sum _ { i = 1 } ^ { N } \left ( \log \Gamma _ { i } - \Gamma _ { i } \tau _ { i } \right ) \, \quad \ \ ( 3 4 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

which allows one to compute the Fisher information matrix I ( Γ ) ij as

$$I ( \mathbf P ) _ { i j } = - \delta _ { i j } \mathbb { E } _ { \tau } \left [ \frac { \partial ^ { 2 } \log p ( \tau | \mathbf P ) } { \partial \Gamma _ { i } \partial \Gamma _ { j } } \right ] \Rightarrow \quad ( 3 5 ) \quad \text{$to$ a} \\ \text{$mat}$$

$$I ( \mathbf P ) _ { i i } = - \mathbb { E } _ { \tau } \left [ - \frac { 1 } { \Gamma _ { i } ^ { 2 } } \right ] = \frac { 1 } { \Gamma _ { i } ^ { 2 } } \,, \quad \quad ( 3 6 ) \quad \underset { \dots } { \text{ure} }$$

where δ ij is the Kronecker delta, enforcing that the Γ 's i are i.i.d. distributed. Therefore, we see that a Jeffreys prior for Γ ℓ must obey

$$p ( \mathbf R ) \, \in \, \sqrt { \det I ( \mathbf \mathbf P ) } = \prod _ { i = 1 } ^ { N } \frac { 1 } { \Gamma _ { i } }. \quad \quad ( 3 7 ) \, \quad \, \quad \, \text{$\stackrel{ \dots } { \text{the} \choose \text{pos} }$}.$$

Such a property is exactly true for a log-uniform prior, as defined by its probability density function (PDF)

$$f ( x ) = \frac { 1 } { x \log \left ( \frac { b } { a } \right ) }, \quad \text{for} \ a \leq x \leq b ; \ a > 0 \quad \text{(38)} \quad \text{specific} \quad \text{there} \quad \text{is see}$$

where the interval [ a, b ] is known as the support and the term log( b/a ) appears after normalizing (37).

It is remarkable that the distribution on Γ ℓ motivated by flow on Γ ℓ using the differentiable simulation aligns exactly with the Jeffreys prior for decay rates! Wehave thus established that such a distribution for decay rates is worthy of study, on both numerical and statistical grounds, with physics motivation in Section V.

## B. Random Stasis from Physical Priors

It is now natural to study in detail how much stasis can emerge for Γ ℓ and Ω (0) ℓ spectra simply drawn from powerlaw and log-uniform priors. We study the mean e -folds of stasis that emerge across several configurations of the log-uniform and power-law prior, which are related to a new exponential model of stasis and the original model, respectively.

We initialize the power-law priors according to draws Γ ℓ ∼ ℓ γ and Ω (0) ℓ ∼ ℓ α , where γ ∈ [ -1 7] , and α = 2 k for k ∈ [ -4 4]. , The procedure for generating power-law draws is given in II B. Similarly, the exponential model prior is initialized according to Γ ℓ ∼ Log-Uniform(10 γ , 10 ) 0 and Ω (0) ℓ ∼ Log-Uniform(10 α , 10 ) where 0 γ, α = 5 k for k ∈ [ -9 0]. , These samples are sort-correlated before entering the simulation. Experiments are run across 100 realizations for N = 100 and with Γ N -1 /H (0) = 0 1 while simultane-. ously enforcing max (Γ ) = 1 ℓ M p .

We see in the top panel of Figure 5 that the powerlaw prior achieves a maximum mean of N = 7 8 . e -folds of stasis for α = 2 -3 and γ = 7. Indeed, we see that stretched spectra in Γ ℓ and compressed spectra in Ω (0) ℓ result in more matter-dominated cosmologies with longer epochs of stasis. Similarly, the exponential model prior achieves a maximum mean of N = 55 e -folds for γ = -45 and α = 0. It is important to note that after normalization, as Ω (0) ℓ ( t (0) ) = 1, this configuration corresponds to a constant mass spectrum where Ω (0) ℓ = 1 /N . This is indeed also a configuration in which the cosmology is matter dominated, as seen in the bottom panel of Figure 5. All configurations in Figure 5 were computed using the sliding-window algorithm for computing ϵ -Stasis, with ϵ = 0 1. .

The lower bounds of α and γ in Figure 5 are chosen to illustrate an effective cutoff in Ω (0) ℓ in which the species completely decouple from the stasis phenomena. That is, their abundance spectra are so stretched that the individual species achieve their peak and decay such that any possible stasis epoch would occur with little to no matter in the Universe. In other words, in order for to achieve stasis with a non-trivial mixture of components, the prior on abundances should allow for an effective number of species to be abundant enough at their peaks such that there is a possibility of a stasis epoch with Ω M &gt; 0. This is seen clearly in both choices of priors, with the cutoff occurring for α ≤ -30 and α ≥ 4.

It is remarkable that stasis can emerge in some robust configurations with random (albeit sort-correlated) abundances and decay rates! We emphasize that there is no optimization or gradient information being used here; this is simply the result of sampling from distributions and sort-correlating the Γ ℓ and Ω (0) ℓ spectra before entering the simulation. Further, it is interesting to see how persistent in e -folds an exponential model of stasis can be compared to a power-law; yielding close to 60 e -folds in certain cases with just N = 100 species.

## IV. STASIS-CONDITIONED BAYESIAN POSTERIORS

We have so far operated in a data-driven setting without any explicit statistical modeling. In this section, we wish to perform a Bayesian analysis on the complete 2 N -dimensional space of Ω (0) ℓ and Γ . Probabilistic inference ℓ on such a high-dimensional space using traditional techniques such as Markov Chain Monte Carlo would be computationally infeasible. Nonetheless, in the current age of machine learning, we can leverage techniques that exploit gradient information to make such tasks more accessible. To this end, we employ stochastic variational inference (SVI), a gradient-based inference technique that approximates complex probability distributions which are often intractable.

We would like to utilize SVI on the task of searching the full input parameter space of θ = Γ { ℓ , Ω (0) ℓ } (keeping in mind that we fix the ratio Γ N -1 /H (0) ) by modeling

the posterior distribution p θ ( |N ) conditioned on optimizing N . At first glance, this seems like a tall order. The parameter space in question is 2 N -dimensional; it must be searched, optimized on, and result in an easilysampled posterior. The crux of what makes this feasible is in the expressivity of neural networks, combined with the information of gradients. In a modern setting, SVI can transform Bayesian inference into a neural-network based optimization problem, thereby accelerating statistical modeling in high-dimensional spaces.

We begin with a review of the essentials of Bayesian and stochastic variational inference, including the evidence lower bound, and then apply these techniques in the context of stasis.

## A. The Evidence Lower Bound

We proceed to define the fundamental object in SVI, which defines the objective function that is regressed on. All Bayesian inference problems begin with Bayes' Theorem, from which we have

$$p ( \theta | \mathcal { N } ) = \frac { p ( \mathcal { N } | \theta ) p ( \theta ) } { p ( \mathcal { N } ) } \, \quad \quad ( 3 9 ) \quad \text{$seen$} \\ \text{$\text{not$}$} \\ \text{qui}$$

where p ( N| θ ) is the likelihood and p θ ( ) is the prior distribution on θ . In a traditional setting, where there are observations to guide the posterior towards, the likelihood encodes the probability of observed data given input parameters. It is sought to be maximized. For our purposes, the likelihood is computed via the simulation S ( θ ) which outputs N e -folds of stasis. Since solving Boltzmann equations is deterministic, the simulation likelihood is exactly

$$p ( \mathcal { N } | \Gamma _ { \ell }, \Omega _ { \ell } ^ { ( 0 ) }, H ^ { ( 0 ) } ) = \delta ( \mathcal { N } - \mathcal { S } ( \Gamma _ { \ell }, \Omega _ { \ell } ^ { ( 0 ) }, H ^ { ( 0 ) } ) ) \, \ \ ( 4 0 ) \quad \text{when}$$

where δ is the Dirac-delta function. We will see that the simulation likelihood is disadvantageous as far as introducing information that aids in searching the full parameter space. For this reason, we adopt a simple optimization likelihood with the definition

$$p ( \mathcal { N } | \Gamma _ { \ell }, \Omega _ { \ell } ^ { ( 0 ) }, H ^ { ( 0 ) } ) \, \in \, e ^ { \kappa \cdot \mathcal { N } } \,, \quad \ \ ( 4 1 ) \, \quad \ \ s u r.$$

which can be potentially subject to a matter abundance constraint similar to that used in equation 23. We emphasize the use of ∝ as this optimization likelihood is, of course, not a normalized probability density. Further, when implementing SVI, we find that the use of a numerical pre-factor κ to increase the 'strength' of the likelihood term is beneficial. This surrogate (unnormalized) likelihood serves the utility of encouraging the posterior distribution to explore parameter space that yields stasis; the exact motivation behind this optimization likelihood will be expanded upon in section IV C. The denominator of equation 39, p ( N ) = ∫ p ( N| θ p θ dθ ) ( ) is known as the Bayesian evidence or marginal likelihood and is in general computationally intractable, as it requires an integral over the entire parameter space.

The crux of SVI is bypassing the direct application of Bayes' theorem by introducing a variational family q ϕ ( θ |N ) with variational parameters ϕ , which should be distinguished from the parameters θ of the statistical model. The variational family can range from anything as simple as a standard Gaussian distribution with a trainable mean and variance parameter, to a complex neural network modeling a density on θ with millions of trainable parameters ϕ . In this work, we will choose the latter.

We would like to maximize the similarity between the variational family and the true posterior by minimizing the Kullback-Leibler (KL) divergence between the two distributions, generally defined as

$$D _ { K L } ( P | | Q ) = \int P ( x ) \log \left ( \frac { P ( x ) } { Q ( x ) } \right ) d x \quad ( 4 2 )$$

̸

for two continuous distributions P and Q . The KL divergence is a type of statistical 'distance', measuring how different the two distributions are. It is always positive semi-definite, satisfying D KL ( P || Q ) ≥ 0; however, it is not a metric in a formal sense as it is not a symmetric quantity (i.e. D KL ( P || Q ) = D KL ( Q P || )).

Proceeding with our derivation, after some algebra and substitutions of Bayes' Theorem we arrive at

$$\underset { \text{$i$leli-} } { \text{$D_{K}L}(q_{\phi}(\theta)||p(\theta|\mathcal{N})}) } = \, - \, \mathbb { E } _ { q _ { \phi } ( \theta ) } \left [ \log p ( \mathcal { N } | \theta ) - \log q _ { \phi } ( \theta ) \\ \text{For} \quad & \quad + \log p ( \theta ) \right ] + \log p ( \mathcal { N } ) \, \quad \quad ( 4 3 )$$

where the first term, known as the Evidence Lower Bound (ELBO), is the only term on that RHS that depends on the variational parameters. This is a fundamental object in SVI. It is easy to interpret the purpose of the ELBO when written as

$$\text{at the} \quad \text{ELBO} ( \phi, \theta, \mathcal { N } ) = \mathbb { E } _ { q _ { \phi } ( \theta ) } \left [ \log p ( \mathcal { N } | \theta ) - D _ { K L } \left ( q _ { \phi } ( \theta ) \| p ( \theta ) \right ) \right ] \,, \\ \text{at the} \quad \text{(44)}$$

where the first term is the expected log-likelihood, introducing information from the likelihood (simulation) to the posterior, and the second term is the KL-divergence between the variational distribution and prior, which ensures the posterior distribution has knowledge of the prior distribution, as required by Bayes' theorem. We are then free to rewrite equation (43) as

$$\underset { \widetilde { \text{her} } } { l i s }, \quad D _ { K L } ( q _ { \phi } ( \theta ) | | p ( \theta | \mathcal { N } ) ) = - \text{ELBO} ( \phi, \theta, \mathcal { N } ) + \log p ( \mathcal { N } ) \.$$

Thus, it is clear to see that minimizing the KL divergence is equivalent to maximizing the ELBO. Moreover, since the KL divergence is positive-definite, the ELBO is the lower bound of the evidence. If q ϕ ( θ ) approximates p θ ( |N ) well, computing the intractable evidence is no longer needed.

We see that the ELBO is a function of both variational family parameters ϕ and simulation parameters θ . It is a

loss function that the variational family is regressing on, so its gradients must also be accessible. What must be computed to optimize q ϕ ( θ ) is the following:

$$& \bar { \nabla } _ { \phi } \text{ELBO} ( \phi, \theta, \mathcal { N } ) = \left \{ \underbrace { - \mathbb { E } _ { q _ { \phi } ( \theta ) } \left [ \nabla _ { \phi } \log q _ { \phi } ( \theta ) \right ] } _ { \text{Variational Distribution} } \quad & \text{Despi} \\ & \underbrace { \mathbb { E } _ { q _ { \phi } ( \theta ) } \left [ \nabla _ { \phi } \log p ( \mathcal { N } | \theta ) \right ] } _ { \text{Likelihood Term} } \right \}, \\ & \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \wedge \quad \cdot \quad \cdot$$

where the expectation value is taken over q ϕ ( θ ). To backpropagate through the likelihood term in updating the variational family parameters, we use the chain rule

$$\nabla _ { \phi } \log p ( \mathcal { N } | \theta ) = \nabla _ { \theta } \log p ( \mathcal { N } | \theta ) \cdot \nabla _ { \phi } \theta \quad ( 4 7 ) \quad \text{can} \,.$$

where we see the first term requires the gradients with respect to θ :

$$\nabla _ { \theta } \log p ( \mathcal { N } | \theta ) = \left \{ \underbrace { \frac { \partial \log p ( \mathcal { N } | \theta ) } { \partial \Gamma _ { \ell } } } _ { \text{Decay Rates} }, \underbrace { \frac { \partial \log p ( \mathcal { N } | \theta ) } { \partial \Omega _ { \ell } ^ { ( 0 ) } } } _ { \text{Abundances} } \right \}. \quad \begin{array} { c } \text{NN a} \\ \text{block} \\ \text{at a ti} \\ \text{sions.} \end{array} \\ \text{One can recognize from the optimization likelihood equa-} \quad \begin{array} { c } \text{NN a} \\ \text{putat} \end{array}.$$

One can recognize from the optimization likelihood equation 41 that log p ( N| θ ) = N , which is precisely the output of the differentiable simulation S θ ( ). It is thus clear to see mathematically why gradients are necessary when implementing SVI - in order to backpropagate through the likelihood to update ϕ , the gradients with respect to θ must be known.

## B. Normalizing Flows

There are several options one can consider for constructing a variational family in SVI. One of the most expressive options is the normalizing flow (NF) [15]. This method operates by transforming samples z from a simple probability density (e.g., Gaussian or Uniform) q 0 ( z ) into a complex posterior density q ϕ ( x ). Here, x represents the transformed variables in the target distribution, achieved through a series of invertible transformations. The invertibility is crucial because it allows for both forward and inverse mappings between the base distribution and the complex target distribution, generally enabling the calculation of probability densities.

NFs are a class of bijective neural networks with trainable parameters ϕ , where q ϕ ( x ) : R d → R d . To learn arbitrarily complex invertible functions, NFs are constructed as the composition of a series of N invertible and bijective functions f :

$$q _ { \phi } ( x ) = f _ { N } \circ f _ { N - 1 } \circ \dots \circ f _ { 1 } ( z ) \,, \quad \ \ ( 4 9 )$$

leveraging the fact that the composition of a set of invertible functions is itself still invertible.

The mapping between two density functions q ϕ ( x ) and q 0 ( z ) is related via the absolute value of Jacobian of the transformation

$$q _ { \phi } ( x ) = q _ { 0 } ( z ) \cdot \left | \det \left ( \mathbf J _ { q _ { \phi } ( z ) } \right ) \right | \,. \quad \quad ( 5 0 )$$

Despite its simple functional form, the Jacobian computation can be prohibitively expensive for highdimensional data or more complex architectures, where it is computed as

$$\left | \det \left ( \mathbf J _ { q _ { \phi } ( z ) } \right ) \right | = \prod _ { i = 1 } ^ { N } \left | \det \left ( \frac { \partial f _ { i } ( z _ { i - 1 } ) } { \partial z _ { i } } \right ) \right | \,. \quad ( 5 1 )$$

This matrix can in general be dense, yielding this computation O ( N · d 3 ) expensive.

There are a variety of methods and architectures one can choose when constructing a NF. To bypass the Jacobian computation expense, we implement the Block Neural Autoregressive Flow (BNAF) [16], a variation of the neural autoregressive flow (NAF) [17]. BNAFs employ a NN architecture where the transformations are applied in blocks, each one generating one dimension of the output at a time, conditioned on the previously generated dimensions. The autoregressive property and block structure of BNAFs leads to more stable training and mitigated computational expense, which are essential for implementing SVI for high-dimensional problems.

The crux of NAFs that makes them so efficient is that one can construct q ϕ ( x ) such that its Jacobian is lowertriangular, and thus its determinant that must be computed in equation 50 is a simple product:

$$\left | \det \left ( \mathbf J _ { q _ { \phi } ( z ) } \right ) \right | = \prod _ { i = 1 } ^ { d } \frac { \partial f _ { i } ( z _ { i - 1 } ) } { \partial z _ { i } } \. \quad \ \ \ ( 5 2 )$$

The Jacobian is then computed with backpropagation which is O ( N · d ) expensive, a clear advantage over the non-autoregressive NF. The typical NAF architecture is a set of functions f ( ) i , where each f ( ) i can be decomposed into 'conditioners' c ( ) i and invertible 'transformers' t ( ) i

$$f _ { \phi } ( x _ { < i } ) = t _ { \phi } ^ { ( i ) } ( x _ { i }, c _ { \phi } ^ { ( i ) } ( x _ { < i } ) ) \. \quad \quad ( 5 3 )$$

Despite this structure satisfying all the basic needs of a flow, we see that the number of parameters scales quadratically. To bypass this expense, the BNAF structure models each t ( ) i ϕ directly as an NN with no accompanying conditioner. These are all necessary ingredients for a high-dimensional problem like stasis; however, this comes at the expense of some functionality. Despite being theoretically invertible, accessing the inverse of the BNAF, which can be used to compute posterior samples' probability densities, is not currently computationally feasible.

## C. Searching for Stasis Theories with SVI

We have now gathered all the essential components to conduct SVI on the stasis simulation, and will use it to

FIG. 6. Stochastic variational inference pipeline. For a given experiment, a Bayesian prior in parameters is chosen which acts as a form of regularization in the ELBO. During training, parameters are sampled from the variational family q ϕ ( θ ), in this work chosen to be a Block Neural Autoregressive Flow (BNAF). Samples are differentiably sorted before entering the stasis simulation, which solves the set of N +1 coupled Boltzmann equations using diffrax to preserve the flow of gradients. The following Ω M ( ) curve is passed into the differentiable stasis finder to isolate the stasis t e -folds N and the asymptotic matter abundance Ω M . The stasis value is used in the likelihood calculation which is factored into the ELBO loss, which is used to iteratively optimize q ϕ ( θ ).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-65aaa43b-86c5-48e4-b650-c1fabb82244b.jpg)

sample configurations of rates and abundances that yield stasis, from which we will be able to understand trends of stasis configurations. For this, we employ numpyro [18], a probabilistic programming library that leverages jax for automatic differentiation in the implementation of SVI. Together with diffrax , these packages can utilize GPU computation, significantly speeding up the inference process.

We now revisit our choice of the surrogate optimization likelihood in equation 41 instead of the simulator likelihood in equation 40. The likelihood function is theoretically a normalized PDF that encodes the probability of the observed data given the model parameters. However, in this theoretical particle cosmology setting, there are no direct observations; instead, we aim to understand the input parameter space that produces epochs of stasis. This raises the question: how can SVI be made appropriate here?

The answer lies in the definition of the optimization likelihood given in equation 41. In a typical Bayesian inference setting, the likelihood function captures the probability of the observed data given the model parameters. This can often be interpreted in terms of residuals between simulation outputs and observed data. Minimizing this residual corresponds to maximizing the loglikelihood which enters the ELBO, and therefore results in a concrete optimization objective. However, instead of following this traditional approach, we adopt a utilityoriented perspective of the likelihood, treating it purely as an optimization objective. The simulation likelihood given in equation 40 is unsuitable for optimization as it is singular and lacks a well-defined gradient. For this

## Algorithm 2: Searching Stasis with SVI

Require: ϕ variational parameters, θ = Γ { ℓ , Ω (0) ℓ , H (0) } simulation parameters, p (Γ ) prior over decay rates, ℓ p (Ω (0) ℓ ) prior over initial abundances, Γ N -1 /H (0) fixed early time decays relative to Hubble constant, ξ early-stopping threshold

1: while not converged do

2: z 0 ∼ q 0 ( z )

▷

Sample initial latent variables

3: (Γ ℓ , Ω (0) ℓ ) ← q ϕ ( z 0 )

▷

Normalizing Flow

4: Ω (0) ≤ Ω (0) ≤ . . . ≤ Ω

1 2 (0) N ▷ Differentiably sort Ω (0) ℓ

5: Γ 1 ≤ Γ 2 ≤ . . . ≤ Γ N

▷

Differentiably sort Γ

ℓ

6: N ← S (Γ ℓ , Ω (0) ℓ , H (0) )

▷

Run simulation

7: p ( N| Γ ℓ , Ω (0) ℓ , H (0) ) ▷ Compute surrogate likelihood

8: L ← ELBO( ϕ, Γ ℓ , Ω (0) ℓ , N )

▷

Compute ELBO

9:

∇ ∝ {∇

ϕ

Ω

(0)

ℓ

L ∇

,

Γ

ℓ

L}

▷ Backpropagate through

stasis simulation

10: Γ ℓ ← clip(Γ ℓ , 0 1) ,

▷ Ensure Γ ℓ physical

11: Γ N - 1 /H (0) ← 0 1 . ▷ Enforce initial decays time scale

12: ∆ ϕ ∝ -∇ L ϕ (Γ ℓ , Ω (0) ℓ )

▷

Update

ϕ

13:

if

L

does not improve over

ξ

epochs or

NaN

loss

then

14:

break

▷

Early-stopping criterion

15: end if

16: end while

reason, we adopt the optimization likelihood.

A schematic of the SVI pipeline is provided in Figure 6, highlighting the four major components of the SVI pipeline: the prior distribution, the variational family, the simulation S ( θ ) and stasis-finder used in the likelihood calculation, and gradient update. A more detailed description of how we employ SVI in this setting is given

FIG. 7. Prior and posterior comparison for choices of power-law prior with Γ ℓ ∼ ℓ 7 and Ω (0) ℓ ∼ ℓ and log-uniform prior with choices of γ = -62 and α = -15 for N = 50. Each individual heat map is a depiction of 1000 samples and their stasis configuration in ( N , Ω M ) space. For both choices of priors, there is a higher degree of stasis in the posterior in both mean and maximum value. The power-law posterior features a mean stasis value of 13.59 e -folds and maximum of 25.24 e -folds, while the log-uniform posterior has a mean stasis value of 31.06 e -folds and maximum of 96.5 e -folds. Over 1000 samples, both priors had &lt; 1% of samples achieve a stasis epoch of more than 10 e -folds, while the power-law and exponential posteriors have 76% and 99%, respectively. While SVI is able to find non-trivial distributions of Γ ℓ and Ω (0) ℓ that results in epochs of stasis, it is clear to see the effect of the prior regularization in optimization in the large discrepancy between posterior configurations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-776fe3da-bf95-4c8b-a4b5-c170961b6728.jpg)

in Algorithm 2, where we are again reminded the importance of preserving differentiability. A methodology such as SVI requires that one be differentiable at every step , from differentiably sorting Ω (0) ℓ and Γ , ℓ to solving the Boltzmann equations while preserving their gradients, and calculating the stasis duration and matter abundance in a differentiable way. This leads to an uninterrupted flow of gradients through the pipeline.

Bayesian machine learning paradigm offers both the benefits of expressivity of neural networks with the statistical rigor of Bayesian analysis. We have introduced an alternative formulation of SVI that allows one to operate without direct observations, yet still exploit the properties and benefits of SVI. This approach generally enables the study of otherwise prohibitively difficult highdimensional parameter spaces, and can be adapted to any physical system that can be simulated.

Conducting these numerical analyses under the

FIG. 8. Optimizing stasis with SVI and the ELBO loss (equation 44) for log-uniform ( α = -15, γ = -62) and power-law ( α = 1, γ = 7) initializations for N = 50 and Γ N -1 /H (0) = 0 1. . A BNAF with two hidden layers and a hidden layer width of 8 neurons was trained for 2000 epochs with Adam optimizer and early-stopping. A batch size of 10 was used in training. ( Left ) A comparison of model fit for 1000 power-law prior and posterior samples. A flow toward a strictly more power-law posterior is seen for Γ , with Ω ℓ (0) ℓ becoming more exponential in the posterior, with a comparable mean R 2 score for both power-law and exponential fits in the posterior. ( Right ) A comparison of model fit for 1000 exponential prior and posterior samples. Posteriors for both Γ ℓ and Ω (0) ℓ become strictly more exponential. This posterior additionally has a much larger mean stasis e -folds and maximum stasis e -folds than the power-law prior.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-00e0901e-506e-4287-99a0-0a0642dde83f.jpg)

## D. Stasis Results with SVI

The analyses of the previous sections exhibited a flow in the space of parameters that preferred an exponential model of stasis. Additionally, it was shown that an exponential model had manifestly more robust stasis epochs than the power-law prior corresponding to the model in [1]. Using both sets of priors, it is then natural to consider whether such flows persist when using neural networks to do inference on the full parameter space.

We run SVI with both priors (power-law and loguniform) for N = 50 and Γ N -1 /H (0) = 0 1. . We denote the power law prior samples as being drawn from Ω (0) ℓ ∼ ℓ α and Γ ℓ ∼ ℓ γ . Similarly, we denote the sort-correlated log-uniform samples as being drawn from Ω (0) ℓ ∼ Log-U(10 α , 10 ) 0 and Γ ℓ ∼ Log-U(10 γ , 10 ). 0 A numerical prefactor of κ = 10 was used in the surrogate likelihood. We utilize a BNAF with two hidden layers and a hidden layer width of 8 neurons. A batch size of 10 was used for all experiments, corresponding to 10 loglikelihood evaluations per training step. We use Adam [19] optimizer with an initial learning rate of η = 01 for . 2000 epochs of training with an early-stopping threshold of ξ = 200 epochs. All training was done on a single NVIDIA A100-80GB GPU.

and Ω (0) ℓ that result in epochs of stasis for both choices of priors, as shown in the N Ω M space heat maps for 1000 posterior samples in Figure 7. Samples from the power-law prior exhibit a mean stasis value of 4.61 e -folds and maximum of 16.51 e -folds, which after optimization is 13.59 e -folds and 25.24 e -folds in the posterior. The posterior/prior discrepancy of the exponential prior experiment is much more drastic, with prior samples exhibiting a mean of 7.46 e -folds and maximum of 19.38 e -folds, which is 31.06 e -folds and 96.49 e -folds in the posterior. While the SVI optimization has worked in both scenarios, it is clear by the discrepancies in posterior configurations the effect that the choice of prior has in the optimization, as shown in equation 44. This was indeed expected from the numerical study of prior configurations shown in Figure 5; it is therefore interesting to study the model dependence of posterior samples and their flow from model dependence in their priors.

We find that the power-law prior was susceptible to a mode collapse of either completely matter dominated or radiation dominated stasis configurations using the likelihood defined in equation 41. To address this, an additional constraint on the likelihood was imposed, similar to the regularization condition when just using gradients in the differentiable simulation, shown in equation 23. The experiments with the log-uniform prior did not exhibit this mode collapse.

We find that SVI is able to learn distributions of Γ ℓ

The flow in the space of parameters of Γ ℓ and Ω (0) ℓ for both experiments is shown in Figure 8. For the powerlaw prior, it is seen that Γ ℓ and Ω (0) ℓ are more confidently fit with a power-law model, as expected. After training, there is a clear drift towards a more power-law model in Γ . ℓ Interestingly, it is only Ω (0) ℓ that displays a drift towards an exponential model with the Ω (0) ℓ posterior featuring a mean fit value of R 2 = 0 95 for both power-. law and exponential, with the mean power-law R 2 value decreasing from R 2 = 0 99 in the prior and the mean ex-. ponential fit R 2 value increasing from R 2 = 0 87 in the . prior. This flow of parameters is contrary to the result in which just the differentiable simulation was used (Figure 3), where both Γ ℓ and Ω (0) ℓ drift towards an exponential. The posterior in this case is a qualitative hybrid of both models of stasis, and even in such conditions the

emergence of stasis epochs is possible. Stasis epochs are not physically exclusive to single models, whether they be power-law or exponential in parameters.

For the log-uniform prior, both Γ ℓ and Ω (0) ℓ posteriors display clear preference towards an exponential model with a mean score of R 2 = 1 for an exponential fit. The mean power-law fit R 2 score is invariant after optimization, but it is seen in Figure 8 that the scatter has significantly decreased.

These results highlight that in a properly Bayesian setting, both a power-law and exponential model for Γ ℓ and Ω (0) ℓ are valid solutions when conditioning on producing epochs of stasis. The drift towards an exponential model of Ω (0) ℓ even when having a power-law prior can indicate that completely power-law distributions of Γ ℓ and Ω (0) ℓ are more sparse in the full parameter space. Additionally, the discrepancy in the statistics of stasis in the two posteriors illuminate a sharp difference in the persistence of stasis epochs for the two types of models - a mean difference of 18 e -folds.

It could be criticized that the flow of the posterior towards an exponential model is purely due to the choice of prior. We see that this is not strictly the case as with the Ω (0) ℓ parameter flow for the power-law prior experiment, in which the parameter flow goes distinctly against the prior regularization in the ELBO. Thus, the prior plays an important role in the ELBO, but SVI can lead to results that are qualitatively at odds with the prior, due to the condition intrinsic to the posterior.

## V. MODELS OF STASIS

Thus far, we have remained as model-agnostic about stasis as possible, choosing only the prior distribution on rates and abundances for both our gradient ascent and Bayesian inference analyses. In both analyses, optimization was chosen to optimize stasis, with a clear preference for exponential models over power-law models.

In this section we wish to study an exact exponential model both analytically and numerically, compare it to exact power-law models, and provide a preliminary analysis of stasis in the String axiverse; more thorough analysis will be reserved for future work. By 'exact' in both contexts, we mean a model with fixed formulae for Γ ℓ and Ω (0) ℓ and no noise, as opposed to the intrinsically noisy draws of previous sections.

## A. An Exact Exponential Model of Stasis

Motivated by our numerical results, we introduce an exponential model of stasis that yields parametrically more stasis than power-law models, as a function of the number of species. We will show that such a model obeys the dynamical equations of stasis and will result in robust, finite epochs of stasis. We additionally derive al- gebraic relations with exact predictions for N and Ω M from our model.i

We begin my parameterizing this model as

$$\Gamma _ { \ell } = \Gamma _ { N } e ^ { \gamma ( \ell - N ) }, \quad \Omega _ { \ell } ^ { ( 0 ) } = \Omega _ { N } ^ { ( 0 ) } e ^ { \alpha ( \ell - N ) } \quad ( 5 4 )$$

for a spectrum of N species indexed by ℓ and general scaling factors α , γ , and Γ N . As in the original model of stasis, α and γ control the hierarchies in abundances and rates, this time with exponential dependence. We are also reminded that Ω (0) ℓ ( t (0) ) = 1, which requires the overall normalization factor Ω (0) N = [ ∑ N ℓ =1 e α ℓ ( -N ) ] -1 . However, in contrast with the model in [1], where the fundamental object was the mass spectrum m ℓ , we choose to parameterize this model of stasis directly with the decay rates and abundances. We are thus ready to define the five parameter model inspired by an exponential:

$$\{ \alpha, \gamma, \Gamma _ { N }, \Omega _ { N } ^ { ( 0 ) }, t ^ { ( 0 ) } \} \,,$$

with which shall proceed to show obeys the dynamical equations required of a stasis epoch.

We begin by reminding ourselves that a model of stasis must satisfy the constraint equations given in equations 10 and 15. Focusing on the time evolution of the total Ω (0) ℓ , we get

$$\sum _ { \ell } \Omega _ { \ell } ( t ) & = \Omega _ { \ell } ^ { ( 0 ) } h ( t ^ { ( 0 ) }, t ) \sum _ { \ell = 1 } ^ { N } e ^ { \alpha ( \ell - N ) } e ^ { - \Gamma _ { N } ( t - t ^ { ( 0 ) } ) e ^ { \gamma ( \ell - N ) } } \,, \\ \text{where} \, \text{$two$ how need the result from annotation $1M$ but}$$

where we have used the result from equation 14, but replaced the time dependence with a more general net gravitational redshift factor h t ( (0) , t ) such that there is no assumption that we are in a period of stasis. h t ( (0) , t ) invariant across all species (i.e. is independent of ℓ ); a more in-depth discussion of this h -factor can be found in [1]. Moving to the N → ∞ limit, we simplify this calculation by transforming the sum into an integral

$$\begin{smallmatrix} i \\ \mathbf y \\ \mathbf n \\ \mathbf n \end{smallmatrix} & \sum _ { \ell = 1 } ^ { N } e ^ { \alpha ( \ell - N ) } e ^ { - \Gamma _ { N } ( t - t ^ { ( 0 ) } ) e ^ { \gamma ( \ell - N ) } } \\ \mathbf r & \to \int _ { 0 } ^ { \infty } e ^ { \alpha ( \ell - N + 1 ) } e ^ { - \Gamma _ { N } ( t - t ^ { ( 0 ) } ) e ^ { \gamma ( \ell - N + 1 ) } } d \ell \,, \quad ( 5 7 )$$

which can be computed by invoking the Gamma function identity Γ( k /b ) k = ∫ ∞ 0 z k -1 e -bz dz, yielding

$$\sum _ { \ell } \Omega _ { \ell } ( t ) = \frac { \Omega _ { N } ^ { ( 0 ) } } { \gamma } & \Gamma \left ( \frac { \alpha } { \gamma } \right ) h ( t ^ { ( 0 ) }, t ) \\ & \times \left [ \Gamma _ { N } ( t - t ^ { ( 0 ) } ) \right ] ^ { - \alpha / \gamma } \,. \quad ( 5 8 )$$

We can solve for the quantity ∑ Γ Ω ( ) similarly, resultℓ ℓ t

ing in

$$\sum _ { \ell } \Gamma _ { \ell } \Omega _ { \ell } ( t ) = \frac { \Omega _ { N } ^ { ( 0 ) } \Gamma _ { N } } { \gamma } & \Gamma \left ( \frac { \alpha } { \gamma } + 1 \right ) h ( t ^ { ( 0 ) }, t ) \\ & \times \left [ \Gamma _ { N } ( t - t ^ { ( 0 ) } ) \right ] ^ { - ( \alpha / \gamma + 1 ) } \,. \ ( 5 9 )$$

It is lastly straightforward to compute the ratio of these conditions, where we invoke the Gamma function identity Γ( z +1) Γ( ) = / z z ,

$$\frac { \sum _ { \ell } \Gamma _ { \ell } \Omega _ { \ell } ( t ) } { \sum _ { \ell } \Omega _ { \ell } ( t ) } = \frac { \alpha } { \gamma } \left ( \frac { 1 } { ( t - t ^ { ( 0 ) } ) } \right ) \,. \quad \quad ( 6 0 ) \quad \,.$$

We see that this result is power-law in the difference ( t -t (0) ), dissimilar to the t -1 dependence from the stasis dynamical equations in equations 10 and 15. Indeed, this discrepancy also arises in the power-law model of stasis. Qualitatively, this means that the stasis epoch must emerge some time after the initial species production time, when t ≫ t (0) and hence ( t -t (0) ) -1 ≈ t -1 . This is when the edge effects, most apparent when Γ N /H (0) ≫ 1, have died away. This further indicates that our model produces finite epochs of stasis, with a natural beginning of the stasis epochs emerging due to the presence of edge effects, and an end when all species have decayed.

With this result in hand, we can proceed to compare coefficients in 15 and 60, in which we find

$$\frac { \alpha } { \gamma } = \frac { 2 ( 1 - \overline { \Omega } _ { M } ) } { 4 - \overline { \Omega } _ { M } } \. \quad \quad \quad ( 6 1 )$$

Wecan further invert this to obtain the prediction for Ω M during stasis in terms of our model parameters, resulting in

$$\overline { \Omega } _ { M } = \frac { 2 ( \gamma - 2 \alpha ) } { 2 \gamma - \alpha } \, \quad \ \ \ ( 6 2 ) \quad \stackrel { \dots } { \text{ex} }$$

from which we can also find a parameter restriction for α and γ by enforcing 0 &lt; Ω M &lt; 1:

$$0 < \alpha < \frac { \gamma } { 2 } \. \quad \quad \ ( 6 3 ) \quad m _ { \ t w }$$

With this result, we can identify the model parameters that result in matter-radiation equality (MRE). Solving equation 61 for Ω M = 1 2, we get the condition /

$$\frac { \alpha } { \gamma } = \frac { 2 } { 7 } \, \quad \quad \quad \left ( 6 4 \right ) \quad \, \quad \, \ t i$$

for MRE. We will study this result numerically and further compare it to the power-law model of MRE.

It remains now to check that the constraints for ∑ ℓ Ω ( ) (equation 10) and ℓ t ∑ ℓ Γ Ω ( ) (equation 15) are ℓ ℓ t individually satisfied. To that end, without loss of generality we can simply show that the condition in equation 10 is satisfied, as we have already checked the condition of their ratio. We begin by operating under t ≫ t (0) , for which we have

$$h ( t ^ { ( 0 ) }, t ) & = h ( t ^ { ( 0 ) }, t _ { * } ) h ( t _ { * }, t ) \\ & = h ( t ^ { ( 0 ) }, t _ { * } ) \left ( \frac { t } { t _ { * } } \right ) ^ { 2 - 6 / ( 4 - \overline { \Omega } _ { M } ) } \,, \quad ( 6 5 )$$

where we have inserted the factor for h t ( ∗ , t ) from equation 14, and t ∗ ≫ t (0) is some fiducial time much later than the initial species production. Inserting this into equation 60, we arrive at

$$\sum _ { \ell } \Omega _ { \ell } ( t ) = \frac { \Omega _ { N } ^ { ( 0 ) } } { \gamma } & \Gamma \left ( \frac { \alpha } { \gamma } \right ) h ( t ^ { ( 0 ) }, t _ { * } ) \left ( \frac { t } { t _ { * } } \right ) ^ { 2 - 6 / ( 4 - \overline { \Omega } _ { M } ) } \\ \intertext { \i c e } _ { \i d, } & \times \left [ \Gamma _ { N } ( t - t ^ { ( 0 ) } ) \right ] ^ { - \alpha / \gamma } \,,$$

from which we see using the result in equation 61 and considering that t ≫ t (0) , the time dependence falls out, ensuring that we are in a period of stasis in which Ω M ≡ ∑ ℓ Ω ( ). ℓ t Our exponential model of stasis is thus consistent with the conditions for stasis as defined from their corresponding dynamical equations.

We are now in a position to estimate exactly the number of e -folds of stasis expected from a configuration of our model, given by:

$$\mathcal { N } \equiv \log \left [ \frac { a ( t = \tau _ { 1 } ) } { a ( t = \tau _ { N } ) } \right ] & = \frac { 2 } { 4 - \overline { \Omega } _ { M } } \log \left ( \frac { \Gamma _ { N } } { \Gamma _ { 1 } } \right ) \\ & = \frac { 2 \gamma } { 4 - \overline { \Omega } _ { M } } \left ( N - 1 \right ) \,, \quad \left ( 6 7 \right ) \\. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

where we have used the result for the time evolution of a t ( ) in equation 12.

We see that the exponential model of stasis yields parametrically more e-folds than the power-law model, which exhibited N ∼ log( N ) scaling. A numerical comparison of e -fold data and the theoretical prediction for γ = 1, α = 2 7, and Γ / N = 0 01 yielding a stasis . epoch with matter-radiation equality (MRE), is shown in Figure 9. Further, see Figure 10 for a comparison of stasis plots with N = 300 species in the exponential and power-law model. As expected, there is excellent agreement between data and theory in the limit N → ∞ . The figure also illustrates that more than 60 e -folds of MRE can be achieved with just 100 species in the exponential model. This result shows that a relatively small number of species is needed to achieve inflation-level e -folding, which may have important phenomenological implications.

We additionally compare e -fold scaling with the powerlaw model of stasis in Figure 9, where α = δ = 1 and γ is being varied. The figure demonstrates the qualitative difference in e -folds for stasis epochs between the two models, even with strong power-law scaling of Γ ℓ ∼ ℓ 7 corresponding to MRE. Specifically, the exponential model achieves the same e -folds of stasis with just N ∼ 50 species, compared to the 400 species required for Γ ℓ ∼ ℓ 7 .

FIG. 9. Comparison between exponential model of stasis with γ = 1, α = 2 7, and Γ / N = 0 01 yielding matter-radiation . equality, and the power-law model of stasis with α = δ = 1. The theoretical prediction of e -fold scaling for the exponential model from equation 67 is shown in the red dashed line, in which we see that in the N → ∞ limit there is exact agreement between numerical data and the theoretical prediction. MRE for the power-law model with γ = 7 is shown in green, with the power-law model exhibiting logarithmic scaling with N and the exponential model exhibiting linear scaling with N as N →∞ . It is also seen that more than 60 e -folds of MRE (black dashed line) is achieved with just N ∼ 100 species for the exponential model.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-ac39f5ab-e257-4951-9f57-4e78a735d0bd.jpg)

The logarithmic scaling in the power-law model and the linear scaling at high N for the exponential model are both evident. Therefore, while both models produce epochs of stasis, it is clear that the exponential model results in longer epochs of stasis.

We lastly recall that in [1], it was shown that a stasis state is a global attractor of the dynamical system. This was demonstrated to be true for the power-law model, and we can proceed to show minimally that it holds for the exponential model as well. In doing so, it is sufficient to see that for the exponential model

$$\sum _ { \ell } \Gamma _ { \ell } \Omega _ { \ell } & = \left ( \frac { \alpha } { \gamma } \right ) \frac { \Omega _ { M } } { t - t ^ { ( 0 ) } } \\ & = \left [ \frac { 2 ( 1 - \overline { \Omega } _ { M } ) } { 4 - \overline { \Omega } _ { M } } \right ] \frac { \Omega _ { M } } { t - t ^ { ( 0 ) } } \, \quad \ \ ( 6 8 ) \quad \ t \text{tion} \quad \text{a to}$$

where we see that this is exactly the result of equation 4.1 in [1] up to a constant factor. It then follows that the subsequent analyses showing that the stasis state is an attractor also apply to our exponential model of stasis.

## B. Stasis and the String Axiverse

Axions arise readily in string theory via the dimensional reduction of higher-form gauge fields; see, e.g., [20, 21]. Notably, they need not be the QCD axion or axions that couple to any dark gauge sector. Accordingly, string theory axions can be relevant in a number of phenomenological or cosmological roles, including as the inflaton [22-25], as the reheaton [26], as the QCD axion [20, 27-29], as particles that couple to photons [30, 31], and more. The largest number of axions in a known Type IIB / F-theory compactification is 181,820 [32] in a single geometry, whereas [33, 34] and [35, 36] give rise to large ensembles of geometries that typically have thousands or hundreds of axions, respectively.

Could axions give rise to stasis? Features that enable robust periods of stasis include a spectrum of N particles that may decay to radiation and non-trivial hierarchies in the decay rates. The former is clearly satisfied in many corners of the string landscape, but in lieu of a detailed analysis we would like to comment on an essential feature that affects rates: in many string constructions, axion masses are generated non-perturbatively and the ℓ th axion masses appears schematically as

$$m _ { \ell } \in e ^ { - c T _ { \ell } }$$

where T ℓ is the volume of the cycle wrapped by the string instanton that generates the mass. Decays to radiation often occur via dimension five operators, in which case we have

$$\Gamma _ { \ell } \in e ^ { - 3 c T _ { \ell } } \quad \Omega _ { \ell } ^ { ( 0 ) } \in e ^ { - \alpha c T _ { \ell } }, \quad \ \ ( 7 0 )$$

where α depends on the production mechanism as discussed in [1]. Though the rates and abundances have power-law dependence on masses, the exponential dependence of the non-perturbative mass causes it to depend exponentially on the cycle volumes, which are themselves sensitive to the details of moduli stabilization. However, due to the non-perturbative effects there is a clear opportunity for exponential hierarchies, and in fact in string theoretic constructions of the models similar to the Standard Model, the exponential dependence on T QCD sets Λ QCD exponentially below the Planck scale. Furthermore, if moduli stabilization distributes T ℓ uniformly, then the rates and abundances are log-uniform and can achieve many e -folds of stasis, as we have shown.

## C. Stasis and the Emergent String Conjecture

The Emergent String Conjecture [37] is a generalization of the Swampland Distance Conjecture [38] in which a tower of either Kaluza-Klein or string states

$$m _ { \ell } ( \phi ) = m _ { \ell } ( \phi _ { 0 } ) e ^ { - \lambda d ( \phi, \phi _ { 0 } ) } \text{ \quad \ \ } ( 7 1 )$$

becomes exponentially light upon going a distance d ϕ, ϕ ( 0 ) from ϕ 0 toward ϕ near the boundary of moduli space. The towers still have power-law masses, and

therefore the number of e -folds of stasis satisfies

$$\mathcal { N } \infty \log ( N ), \quad \ \ ( 7 2 ) \quad \ t h$$

but the exponential suppression of the lowest mass state yields an exponentially large number of states below a fixed cutoff. We therefore have

$$\mathcal { N } \varpropto \frac { \lambda d ( \phi, \phi _ { 0 } ) } { \delta }, \quad \quad \quad ( 7 3 ) \quad \text{$\text{$\text{$fol$}$}$}$$

where m ℓ ∝ ℓ δ . We therefore expect stasis to become increasingly important as the boundary of moduli space is approached, potentially leading to experimental constraints or observational consequences.

## VI. DISCUSSION

In this paper, we have studied M → γ cosmological stasis. In this scenario, a tower of matter states initially dominates the energy density of the universe, but subsequent decays to radiation balance against Hubble expansion, leading to a cosmological epoch with approximately constant matter and radiation abundances.

A central result of this paper is a new exponential model of stasis that was motivated by our numerical approach, searching the full 2 N -dimensional parameter space of decay rates Γ ℓ and initial abundances Ω (0) ℓ using gradients from a differentiable Boltzmann solver, as well as stochastic variational inference with neural networks. This model leads to a longer duration of stasis in e -folds relative to a power-law model and could be motivated by non-perturbative effects in string theory. We elaborate on and review our main results in what follows.

A stasis epoch could have significant cosmological implications [1], potentially affecting dark matter production, large-scale structure formation, and the overall understanding of cosmological evolution and the age of the Universe. Of particular interest for M → γ stasis is the possibility of stasis epochs occurring after a matterdominated epoch at the end of inflation or after a later matter-dominated epoch prior to nucleosynthesis. In the first case, the decay of the inflaton can instead source ϕ ℓ states. The hierarchical decays of those states during stasis would then be the source of inflationary reheating, with the universe entering radiation domination with the conclusion of the decays. In the second case, the tower of matter fields can lead to an early matter dominated epoch prior to BBN, but well after inflation, the decays of which give rise to stasis. There are a number of potential phenomenological implications of such an epoch, but a stasis state is an attractor regardless of the arena in which ϕ ℓ are produced.

In this work, we have taken a model-agnostic datadriven approach to studying theories of stasis, employing analytic methods only in the final stages. Such an analysis is complementary to a model-driven approach, as e.g. in the initial model of stasis that was inspired by Kaluza-Klein excitations, motivating a mass spectrum that follows a power-law in the species index. Together, the model-driven and model-agnostic approaches demonstrate that stasis is a very general phenomenon.

Our methodology began with constructing a differentiable Boltzmann solver, which is the crux of our approach. We maximized the number of e -folds of stasis by following gradients in our differentiable simulation with model-agnostic (aside from the prior) random initializations of Γ ℓ and Ω (0) ℓ . Such stasis-maximizing trajectories through the space of rates and abundances motivated the study of log-uniform distributions, which we compared to power-law distributions by taking random samples and solving the Boltzmann equations. This comparison provided an initial indication of the significant discrepancy in the duration of stasis epochs between the two models. With distributional priors established, we employed SVI with normalizing flows for a Bayesian analysis of the full, high-dimensional parameter space; there are some essential differences from standard SVI that we discuss in the text. From the posteriors, we again saw that stasis maximization generally prefers log-uniform distributed Ω (0) ℓ and Γ , corresponding to a model that is exponential in ℓ the species index ℓ . This motivates a hybrid power-law / exponential model beyond the strict power-law or exponential models we have studied, further indicating the generality of stasis.

The culmination of our numerical analyses is the analytic exponential model of stasis studied in Section V. Although our study of stasis does not assume specific physical mechanisms to source the ϕ ℓ , we find in Section VB that such a model of stasis could be naturally realized within the string axiverse or near the boundaries of moduli space in accord with the emergent string conjecture. This new model of stasis is also an attractor.

In addition to the analytic model it inspired, our numerical methodology allows us to conclude with the following observations:

- · Deviations from strict stasis (e.g., ϵ -stasis that could have small oscillations), which can still be phenomenologically relevant, indicate that Ω (0) ℓ and Γ ℓ drawn from certain families of distributions can result in epochs of stasis.
- · Stasis epochs can arise without strictly correlated rates and abundances, with the species participating in stasis having abundances increasing with decay rates.
- · There is a numerical preference for an exponential model when maximizing stasis; this was observed in both gradient ascent and SVI experiments.
- · An exponential model of rates and abundances yields more e -folds of stasis than the power-law model, with both analytics and numerics demonstrating that N ∝ N in the exponential model whereas N ∝ log ( N ) in the power-law model.

FIG. 10. Comparison of matter-radiation equality (MRE) between the power-law (left) and exponential (right) models of stasis, both for N = 300 species. MRE in the power law model from [1] corresponds to α = δ = 1 and γ = 7 and with Γ N -1 /H (0) = 0 01. . MRE in the exponential model corresponds to γ = 1, α = 2 7, and Γ / N /H (0) = 0 01. . We see that the power-law model achieves a stasis configuration of ∼ 17 e -folds with 300 species, while the exponential model achieves ∼ 165 e -folds, demonstrating the qualitative difference in e -folds between the two models, attributable to different scaling of e -folds with N .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00835v3-figure-e15cca08-3abc-468e-9440-e90dccb98248.jpg)

In our opinion, these model-agnostic numerical analyses together with model-driven analytic considerations together make stasis a very compelling phenomenon from a theoretical perspective. The generality of the alternative cosmological histories implies that in complex cosmologies, such as those in string theory, one should not make assumptions about the history and instead one must simply solve the Boltzmann equations.

In future work, it would be interesting to study other flavors of stasis, notably Λ → γ , Λ → M , and triple stasis with Λ → M → γ , using the numerical techniques presented in this paper. To that end, extending the differentiable methodology to include other energy pumps, such as the overdamped/underdamped transition described in [2] used to model Λ → M , would be essential. From a numerical perspective, these transitions can be challenging to implement due to their instantaneous (i.e., nondifferentiable) nature; however, one can consider a family of numerical approximations (e.g., tanh , sigmoid ) to model such transitions in a differentiable way.

It would also be interesting to change the optimization objective in future work. All of our numerical analyses were aimed at ensuring robust periods of stasis by optimizing the number of stasis e -folds N . Of course, this optimization is different from being physically optimal, which depends crucially on theory priors, as well as phenomenological and cosmological viability. The latter considerations motivate different optimization objectives that could easily be implemented by adapting our open-source code /github , provided that the objective is implemented in a differentiable way. Doing so would open up new avenues for different physical studies of stasis.

We have seen that studies of theories with highdimensional parameter spaces using differentiable sim- ulators, despite being arduous work, can be essential to understanding the physics. Specifically, we have demonstrated the power of gradients from these simulators to direct flows in parameter space, both with and without neural networks. The methodology presented here is flexible and easily adaptable. Indeed, these techniques, in conjunction with neural networks, can provide profound physical insights, motivate new physical models, and have far-reaching implications in various scientific fields, including those beyond particle cosmology.

## ACKNOWLEDGMENTS

The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University. We thank Siddharth Mishra-Sharma and Carolina Cuesta-Lazaro for initial guidance on SVI, Keith Dienes for discussions regarding stasis, and the authors of [1] for comments on a draft of this manuscript. This work is supported by the National Science Foundation under CAREER grant PHY-1848089 and Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions).

- [1] K. R. Dienes, L. Heurtier, F. Huang, D. Kim, T. M. Tait, and B. Thomas, Stasis in an expanding universe: A recipe for stable mixed-component cosmological eras, Physical Review D 105 , 10.1103/physrevd.105.023530 (2022).
- [2] K. R. Dienes, L. Heurtier, F. Huang, T. M. P. Tait, and B. Thomas, Stasis, stasis, triple stasis (2023), arXiv:2309.10345 [astro-ph.CO].
- [3] K. R. Dienes, L. Heurtier, F. Huang, D. Kim, T. M. P. Tait, and B. Thomas, Primordial black holes place the universe in stasis (2023), arXiv:2212.01369 [astroph.CO].
- [4] K. R. Dienes, L. Heurtier, F. Huang, T. M. P. Tait, and B. Thomas, Cosmological stasis from dynamical scalars: Tracking solutions and the possibility of a stasis-induced inflation (2024), arXiv:2406.06830 [astro-ph.CO].
- [5] S. S. Schoenholz and E. D. Cubuk, Jax m.d. a framework for differentiable physics, in Advances in Neural Information Processing Systems , Vol. 33 (Curran Associates, Inc., 2020).
- [6] J.-E. Campagne, F. Lanusse, J. Zuntz, A. Boucaud, S. Casas, M. Karamanis, D. Kirkby, D. Lanzieri, A. Peel, and Y. Li, Jax-cosmo: An end-to-end differentiable and gpu accelerated cosmology library, The Open Journal of Astrophysics 6 , 10.21105/astro.2302.05163 (2023).
- [7] M. Hoffman, D. M. Blei, C. Wang, and J. Paisley, Stochastic variational inference (2013), arXiv:1206.7051 [stat.ML].
- [8] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, Hybrid Monte Carlo, Phys. Lett. B 195 , 216 (1987).
- [9] P. Kidger, On Neural Differential Equations , Ph.D. thesis, University of Oxford (2021).
- [10] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
- [11] A. Kværnø, Singly diagonally implicit runge-kutta methods with an explicit first stage, BIT Numerical Mathematics 44 , 489 (2004).
- [12] P. Stumm and A. Walther, New algorithms for optimal online checkpointing, SIAM Journal on Scientific Computing 32 , 836 (2010).
- [13] Q. Wang, P. Moin, and G. Iaccarino, Minimal repetition dynamic checkpointing algorithm for unsteady adjoint calculation, SIAM Journal on Scientific Computing 31 , 2549 (2009).
- [14] F. Petersen, C. Borgelt, H. Kuehne, and O. Deussen, Differentiable sorting networks for scalable sorting and ranking supervision (2021), arXiv:2105.04019 [cs.LG].
- [15] D. J. Rezende and S. Mohamed, Variational inference with normalizing flows (2016), arXiv:1505.05770 [stat.ML].
- [16] N. D. Cao, I. Titov, and W. Aziz, Block neural autoregressive flow (2019), arXiv:1904.04676 [stat.ML].
- [17] C.-W. Huang, D. Krueger, A. Lacoste, and A. Courville, Neural autoregressive flows (2018), arXiv:1804.00779 [cs.LG].
- [18] D. Phan, N. Pradhan, and M. Jankowiak, Composable effects for flexible and accelerated probabilistic pro-
- gramming in numpyro, arXiv preprint arXiv:1912.11554 (2019).
- [19] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization (2017), arXiv:1412.6980 [cs.LG].
- [20] P. Svrcek and E. Witten, Axions In String Theory, JHEP 06 , 051, arXiv:hep-th/0605206.
- [21] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D 81 , 123530 (2010), arXiv:0905.4720 [hep-th].
- [22] S. Dimopoulos, S. Kachru, J. McGreevy, and J. G. Wacker, N-flation, JCAP 08 , 003, arXiv:hep-th/0507205.
- [23] J. E. Kim, H. P. Nilles, and M. Peloso, Completing natural inflation, JCAP 01 , 005, arXiv:hep-ph/0409138.
- [24] E. Silverstein and A. Westphal, Monodromy in the CMB: Gravity Waves and String Inflation, Phys. Rev. D 78 , 106003 (2008), arXiv:0803.3085 [hep-th].
- [25] L. McAllister, E. Silverstein, and A. Westphal, Gravity Waves and Linear Inflation from Axion Monodromy, Phys. Rev. D 82 , 046003 (2010), arXiv:0808.0706 [hepth].
- [26] J. Halverson, C. Long, B. Nelson, and G. Salinas, Axion reheating in the string landscape, Phys. Rev. D 99 , 086014 (2019), arXiv:1903.04495 [hep-th].
- [27] J. P. Conlon, The QCD axion and moduli stabilisation, JHEP 05 , 078, arXiv:hep-th/0602233.
- [28] M. Demirtas, N. Gendler, C. Long, L. McAllister, and J. Moritz, PQ axiverse, JHEP 06 , 092, arXiv:2112.04503 [hep-th].
- [29] N. Gendler and D. J. E. Marsh, Qcd axion dark matter in string theory: Haloscopes and helioscopes as probes of the landscape (2024), arXiv:2407.07143 [hep-th].
- [30] J. Halverson, C. Long, B. Nelson, and G. Salinas, Towards string theory expectations for photon couplings to axionlike particles, Phys. Rev. D 100 , 106010 (2019), arXiv:1909.05257 [hep-th].
- [31] N. Gendler, D. J. E. Marsh, L. McAllister, and J. Moritz, Glimmers from the axiverse (2023), arXiv:2309.13145 [hep-th].
- [32] Y.-N. Wang, On the elliptic calabi-yau fourfold with maximal h1, 1, Journal of High Energy Physics 2020 , 1 (2020).
- [33] J. Halverson, C. Long, and B. Sung, Algorithmic universality in F-theory compactifications, Phys. Rev. D 96 , 126006 (2017), arXiv:1706.02299 [hep-th].
- [34] W. Taylor and Y.-N. Wang, Scanning the skeleton of the 4D F-theory landscape, JHEP 01 , 111, arXiv:1710.11235 [hep-th].
- [35] M. Kreuzer and H. Skarke, Complete classification of reflexive polyhedra in four-dimensions, Adv. Theor. Math. Phys. 4 , 1209 (2000), arXiv:hep-th/0002240.
- [36] M. Demirtas, C. Long, L. McAllister, and M. Stillman, The Kreuzer-Skarke Axiverse, JHEP 04 , 138, arXiv:1808.01282 [hep-th].
- [37] S.-J. Lee, W. Lerche, and T. Weigand, Emergent strings from infinite distance limits, Journal of High Energy Physics 2022 , 1 (2022).
- [38] H. Ooguri and C. Vafa, On the geometry of the string landscape and the swampland, Nuclear physics B 766 , 21 (2007).
James Halverson Sneh Pandya
On the Generality and Persistence of Cosmological Stasis
Hierarchical decays of $N$ matter species to radiation may balance against
Hubble expansion to yield stasis, a new phase of cosmological evolution with
constant matter and radiation abundances. We analyze stasis with various
machine learning techniques on the full $2N$-dimensional space of decay rates
and abundances, which serve as inputs to the system of Boltzmann equations that
governs the dynamics. We construct a differentiable Boltzmann solver to
maximize the number of stasis $e$-folds $\mathcal{N}$. High-stasis
configurations obtained by gradient ascent motivate log-uniform distributions
on rates and abundances to accompany power-law distributions of previous works.
We demonstrate that random configurations drawn from these families of
distributions regularly exhibit many $e$-folds of stasis. We additionally use
them as priors in a Bayesian analysis conditioned on stasis, using stochastic
variational inference with normalizing flows to model the posterior. All three
numerical analyses demonstrate the generality of stasis and point to a new
model in which the rates and abundances are exponential in the species index.
We show that the exponential model solves the exact stasis equations, is an
attractor, and satisfies $\mathcal{N}\propto N$, exhibiting inflation-level
$e$-folding with a relatively low number of species. This is contrasted with
the $\mathcal{N}\propto \log(N)$ scaling of power-law models. Finally, we
discuss implications for the emergent string conjecture and string axiverse.