mochi_class: Modelling Optimisation to Compute Horndeski In class

We introduce mochi_class, an extension of the Einstein-Boltzmann solver hi_class, designed to unlock the full phenomenological potential of Horndeski gravity. This extension allows for general input functions of time without the need for hard-coded parametrisations or covariant Lagrangians. By replacing the traditional $\alpha$-parametrisation with a set of stable basis functions, mochi_class ensures that the resulting effective theories are inherently free from gradient and ghost instabilities. Additionally, mochi_class features a quasi-static approximation implemented at the level of modified metric potentials, enhancing prediction accuracy, especially for models transitioning between a super- and sub-Compton regime. mochi_class can robustly handle a wide range of models without fine-tuning, and introduces a new approximation scheme that activates modifications to the standard cosmology deep in the matter-dominated era. Furthermore, it incorporates viability conditions on the equation of motion for the scalar field fluctuations, aiding in the identification of numerical instabilities. Through comprehensive validation against other Einstein-Boltzmann solvers, mochi_class demonstrates excellent performance and accuracy, broadening the scope of hi_class by facilitating the study of specific modified gravity models and enabling exploration of previously inaccessible regions of the Horndeski landscape. The code is publicly available at https://github.com/mcataneo/mochi_class_public


INTRODUCTION
The current cosmological paradigm (ΛCDM) is capable of describing the dynamics of the universe, the Cosmic Microwave Background (CMB) and the Large-Scale Structure (LSS) of the universe with only six parameters (Aghanim et al. 2018).On top of standard model particles, it assumes two exotic components whose nature is unknown: (i) a cosmological constant (Λ) and some form of Cold Dark Matter (CDM).The gravitational interactions are regulated by General Relativity (GR).This modelling explains the observed accelerated expansion of the universe (Riess et al. 1998;Perlmutter et al. 1998), and it provides a good fit to the data (e.g., Heymans et al. 2020;eBOSS Collaboration et al. 2020;DES Collaboration et al. 2021;Miyatake et al. 2023).
However, the standard cosmological model is unsatisfactory as it suffers from at least three drawbacks: (i) the value of the observed cosmological constant is too small to be explained by fundamental physics (Martin 2012), (ii) GR is tested very precisely on local scales (Will 2014) but it is extrapolated over 15 orders of magnitude in length scale to be applied to cosmological data (Peebles 2004), and (iii) comparing different datasets produces tensions for the measured parameters (Valentino et al. 2021;Abdalla et al. 2022).In particular, the most significant is the Hubble tension, arising when combining * mcataneo@uni-bonn.de"early" and "late" times measurements of the expansion rate of the universe (Verde et al. 2019).
To alleviate these problems several possible extensions have been proposed, either promoting Λ to a dynamical field (Dark Energy, DE) (see, e.g., Copeland et al. 2006, for a review) or modifying the laws of gravity at cosmological scales (MG) (see, e.g., Clifton et al. 2011;Joyce et al. 2014Joyce et al. , 2016, for reviews), for reviews).A unifying framework is realised by the Horndeski theory, that encompasses many of the scalar-tensor theories proposed in the literature, e.g.Quintessence, f (R) gravity, Brans-Dicke, Galileons.Horndeski gravity provides one of the most general scalar-tensor theories with at most second-order derivatives in the equations of motion on any background (Horndeski 1974;Deffayet et al. 2009a;Kobayashi et al. 2011).It gained significant popularity due to its rich phenomenology while remaining analytically tractable.
To study the Horndeski landscape, there are two potential strategies (e.g., Ezquiaga & Zumalacárregui 2018).It is possible to: (i) select a specific model based on its covariant Lagrangian formulation and analyse the effects of the additional degree of freedom across different regimes, from cosmological scales to compact objects and black holes; or (ii) adopt an effective-theory approach (EFT), focusing solely on the energy scales relevant to cosmology.
While strategy (i) allows for a detailed investigation of particular scalar-tensor theories, it demands signifi-cant resources and time investment.In the absence of a compelling alternative to GR, strategy (ii) provides a more general and efficient method to detect or constrain deviations from standard gravity (e.g., Ezquiaga & Zumalacárregui 2017;Creminelli & Vernizzi 2017;Baker et al. 2017;Sakstein & Jain 2017).
Given the large classes of models that the Horndeski Lagrangian encapsulates, there have been countless attempts to get cosmological constraints on specific parametrisations.Following each one of the two distinct approaches described above, in the literature it is possible to find both constraints based on sub-classes of Horndeski (see, e.g., Joudaki et al. 2022;Zumalacárregui 2020;Ye & Silvestri 2024), and on different parametrisations of the EFT functions (see, e.g., Bellini et al. 2015;Kreisch & Komatsu 2017;Noller & Nicola 2018;Mancini et al. 2019;Seraille et al. 2024;Castello et al. 2024;Chudaykin & Kunz 2024).A promising but less explored direction is to get model independent constraints.The idea is to try to "reconstruct" gravity without imposing an a priori time dependence for the EFT functions (Zhao et al. 2009(Zhao et al. , 2010;;Hojjati et al. 2012;Raveri 2019;Pogosian et al. 2021;Raveri et al. 2021;Park et al. 2021).
The public version of hi class (Zumalacarregui et al. 2016;Bellini et al. 2019), an extension of the Einstein-Botzmann solver class (Blas et al. 2011;Lesgourgues 2011a), implements both approaches, allowing users to easily implement new models belonging to Horndeski scalar-tensor theories or the EFT framework with the basis proposed in Bellini & Sawicki (2014).In a second release of the code, Bellini et al. (2019) included the quasi-static approximation scheme for improved computational efficiency and to extend the code's applicability to models where following the full dynamics of the scalar field can be extremely challenging.
Even with these significant advances, further steps are required to reliably extend hi class across a broader range of the Horndeski landscape.First, we need to simplify the implementation of specific models and parametrisations, which currently requires modifications to the source code.Additionally, efficient selection of stable models is crucial to avoid extensive searches for parameter combinations free from ghost and gradient instabilities.The accuracy of the QSA implementation in hi class can be improved.Furthermore, addressing numerical noise from limited precision in the calculation of the scalar field speed of sound is essential to prevent incorrectly labeling healthy models as unstable.
To address these issues, in this paper we introduce mochi class, a numerical tool to facilitate the study and statistical analysis of Horndeski models expressed in the effective-theory language.To improve flexibility, mochi class externalises the models definition, allowing to pass to the code precomputed arrays for the time evolution of the background and the EFT functions.In addition, it implements a new stable parametrisation, to satisfy the stability conditions by construction.Furthemore, we implement a new integrated approximation scheme to switch on gravity modifications only after a certain (user specified) redshift.This fixes the early time evolution to ΛCDM, avoiding numerical instabilities when deviations from standard gravity are irrelevant.Finally, mochi class implements an additional QSA approximation implemented at the level of modified metric poten-tials in Newtonian gauge.This offers a cross-check with the standard hi class implementation and improves the agreement with the exact solution.
Throughout the paper we assume the speed of light c = 1, use the class normalisation for the components of the stress-energy tensor (e.g., 3ρ class = 8πGρ standard ), and adopt the following values for the cosmological parameters: for the fractional energy-density of the cold dark matter and baryons, we set Ω c h 2 = 0.120108 and Ω b h 2 = 0.022383; for the Hubble constant, h = 0.6781; for the amplitude and slope of the primordial power spectrum, 10 9 A s = 2.10055 and n s = 0.96605; for the optical depth at reionisation, τ = 0.054308.In what follows we also assume that the background is described by a flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric with signature (−+++).
The paper is organised as follow.In Section 2 we introduce the Horndeski models and their EFT description.In Section 3 we describe the code structure and implementation, and in Section 4 we validate it against other codes.In Section 5 we present a new parametrisation, and leave Section 6 for the summary and outlook.
2. HORNDESKI'S THEORY AND PARAMETRISATIONS In our analysis we focus on sub-classes of the Horndeski theory, one of the most general scalar-tensor theories with at most second-order equations of motion on any background (Horndeski 1974;Deffayet et al. 2009a;Kobayashi et al. 2011).Its action can be written as L 2 = G 2 (ϕ, X) , (1b) Here, g µν is the metric tensor, L m is the matter Lagrangian and G 2, 3, 4, 5 are arbitrary functions of the scalar field ϕ and its canonical kinetic term X = −ϕ ;µ ϕ ;µ /2.In a cosmological setup, there are two complementary approaches that are commonly used.Each one has its own advantages and motivations.For completeness we quickly review both, as it will be useful for the rest of the paper.
A first approach consists in directly specifying the field dependence of the G 2, 3, 4, 5 functions.This is the case of e.g.quintessence or f (R) gravity.Assuming a FLRW background one can get the evolution of the matter and the scalar field at the background level.Then, it is possible to get the evolution of the perturbations on top of the background obtained.The key advantage of this approach is self-consistency.For a given covariant La-grangian formulation one can derive the necessary equations in every regime (background vs perturbations, weak vs strong gravity) and apply the result of one regime to the others.This is extremely important, since the laws of gravity and DE should be universal.However, it is rather cumbersome to jump from one model to another, because the background evolution is non trivial and the effort to setup the machinery needed to solve the physical system under investigation has to be repeated for each model.Then, the main limitation of this approach is that without a favorite model in mind it is difficult to scan a significant part of the parameter space of gravity.
Alternatively it is possible to use an Effective Field Theory (EFT) approach.This framework compresses the information content of the full Equation (1a) into few functions of time.Assuming a FLRW background and up to linear order in perturbation theory a common choice is to describe the linear evolution of cosmological perturbations in Horndeski with (Bellini & Sawicki 2014) Here, w ϕ = p ϕ /ρ ϕ is the equation of state of the scalar field (or we can use any other function describing the background evolution of DE), α K has been dubbed kineticity, α B is called braiding, α T is the tensor-speed excess and M 2 * an effective Planck mass.M 2 * can also be replaced with its run-rate together with the initial condition M * ,ini .Here and throughout, a represents the scale factor of the FLRW metric.These functions describe all the possible dynamics of the Horndeski Lagrangian and they are independent of each other.w ϕ is the only function affecting the background evolution (and the perturbations through a different expansion history), while the other functions in Equation (2) modify only the evolution of the perturbations.As a consequence, a commonly used strategy is to fix the background evolution to the one predicted by some fiducial model (typically ΛCDM) and focus on the phenomenology of the perturbations.In ΛCDM the EFT functions in Equation 2 take the values (4) One of the key advantages of the EFT approach is that it is very powerful to test deviations from the standard cosmology.Indeed, the usual (but not unique) approach is to choose a time dependence of the EFT functions such that ΛCDM is recovered at early times.The EFT framework can be seen as a bridge between the underlying Horndeski theory and observations (see Figure 1 of Bellini et al. 2019), thus facilitating the inference of cosmological constraints.Then, for a given time evolution of the EFT functions it is possible to get a particular Horndeski covariant model following the strategy introduced in (Kennedy et al. 2017(Kennedy et al. , 2018)).
In particular, we choose to work within a sub-class of Horndeski models, specifically those that have α T = 0.The nearly-simultaneous detection of the Gravitational Wave (GW) signal (GW170817) and the corresponding gamma-ray burst (GRB 170817A) emitted from the merger of a binary neutron star system (LIGO Scientific Collaboration et al. 2017) suggests that the speed of GW (c 2 T ≡ 1 + α T ) has to be equal to the speed of light for any cosmological purpose (Baker et al. 2017;Creminelli & Vernizzi 2017;Ezquiaga & Zumalacárregui 2017;Sakstein & Jain 2017).While we make this assumption for simplicity, it is important to stress that there are a number of caveats that may invalidate it.In particular the speed of tensors can be frequency dependent, while the LIGO/VIRGO detection only probes > O(1)Hz frequencies.Then, it is possible to accommodate all observations by assuming a time varying cosmological α T and α T = 0 on LIGO/VIRGO frequencies.In addition, one should discuss about the validity of the EFT framework itself at the LIGO scales (de Rham & Melville 2018).However, given that this has to be intended as a first exploration of model optimization with mochi class, we prefer to keep the parameter space under investigation as little as possible.
As explained above, the background of this class of models can be fully specified by one function of time.It can be the scalar field energy-density, the scalar field equation of state, or directly the Hubble rate.This property is manifest by construction when working within the EFT framework.If working with the Horndeski G i functions one can project them into the background function chosen (see, e.g., Equations (A.6) and (A.7) in Zumalacarregui et al. 2016).
The evolution of the perturbations depends on both the chosen expansion history and on the other EFT functions in a non-trivial way.An in depth discussion that adopts the same notation used here can be found in, e.g., Bellini & Sawicki (2014).Here it is important for the rest of the discussion to report the definition of the de-mixed kinetic term for the scalar field perturbation, and of its effective sound speed where H is the Hubble parameter, ρ m and p m are the energy-density and pressure, respectively, of the matter species (excluding the scalar field), and primes denote conformal time derivatives.A definition of the α i 's as a function of the Horndeski G i 's can be found in, e.g., Bellini & Sawicki (2014).
The models we are going to focus on are: Cubic Covariant Galileon.The covariant Galileon (Deffayet et al. 2009b) model corresponds to the sub-class of Horndeski theories, Equation (1a), that possesses the Galilean symmetry in a flat spacetime, i.e. ∂ µ ϕ → ∂ µ ϕ + b µ , with b µ a constant 4-vector (Nicolis et al. 2008).While the most general version of the covariant Galileon does not respect the condition α T = 0, here we focus on the cubic Galileon case, for which the G i functions read Here, as usual, we have set Λ 3 3 ≡ H 2 0 M Pl .On top of the initial conditions (IC) for the scalar field ϕ and its time derivative φ ≡ dϕ/dt, c 2 and c 3 are extra parameters.In particular, c 2 can be normalized to ±1 through a convenient redefinition of the scalar field.The sign of c 2 is crucial for the expansion history of the universe.If c 2 = +1 one needs a cosmological constant to drive the accelerated expansion of the universe.On the contrary, if c 2 = −1 the model self-accelerates.Since the model is shift symmetric, the IC chosen for the scalar field itself is irrelevant, while the IC for its derivative can be fixed by following the attractor solution H φ = const.(Felice & Tsujikawa 2010a).Additionally, c 3 is fixed by requiring that the sum of the fractional energy-densities today is equal to one in a flat universe.Finally, we choose the self-accelerating branch (c 2 = −1), which gives where Ω ϕ,0 is the fractional energy-density of the scalar field today, and the Hubble rate is the only timedependent quantity.Kinetic Gravity Braiding.The Kinetic Gravity Braiding (KGB) model in its original formulation (Deffayet et al. 2010) is a cubic Horndeski theory.Here we focus on a model (nKGB) defined as (10) where Λ (4n−1) ≡ M (2n−1) Pl

H 2n
0 .The sign of the standard kinetic term has been chosen negative, to have a self-accelerating universe.g and n are the two free parameters of the theory, and the exponent of g has been chosen such that gΩ ϕ (a = 1) ∼ O(1) for every choice of n.It is easy to see that Equation ( 10) is a generalisation of Equation ( 8), obtained by fixing n = 1 and g = 4c 2 3 .In nKGB the only non-zero α i functions are and here overdots represent derivatives with respect to cosmic time.In our analysis we fix n = 5, while g is obtained from requiring the sum of the fractional energy densities of all species to be one today.f (R) gravity.This scalar-tensor theory extends the Einstein-Hilbert action by adding a non-linear function of the Ricci scalar, f (R).After defining ϕ ≡ 1 + f R , with f R ≡ df /dR, it can be reformulated in terms of the Horndeski G i functions as follows (see, e.g., Felice & Tsujikawa 2010b): From the definition of the α i functions (Bellini & Sawicki 2014), we can immediately see that: (14) Note that here and in the remainder of this paper overdots denote derivatives with respect to the natural logarithm of the scale factor.
The evolution of the scalar field and the background expansion are related by the modified Friedmann equation (see, e.g., Hu & Sawicki 2007) coupled with the equation for the Ricci scalar In Equation ( 15), f RR = df R /dR, and ρ m includes all matter species, that is, radiation, non-relativistic matter, neutrinos etc.Thus, we can either specify the expansion history and solve for f (R), or provide an explicit functional form for the extension to the Einstein-Hilbert action and derive the associated background evolution.In the first case, we adopt the so-called designer approach (Song et al. 2006;Pogosian & Silvestri 2007), and in this work, we fix the background to that of a ΛCDM cosmology.For this family of models, deviations from standard gravity are generally parametrised in terms of the dimensionless background quantity and we will consider a model with B 0 ≡ B(z = 0) = 0.01.Alternatively, we can adopt an explicit Lagrangian formulation, with the high-curvature limit of the Hu & Sawicki (2007) form being one of the most extensively studied models.Here c 1 , c 2 and n are free dimensionless parameters, and the mass scale m 2 ≡ ρ c + ρ b , where ρ c and ρ b are the background energy-density of cold dark matter and baryons, respectively.We can eliminate one parameter by requiring convergence to a cosmological constant in the early universe (i.e R ≫ m 2 ), which gives c 1 /c 2 = 6 ΩΛ /Ω m , with ΩΛ = λΩ Λ , and λ is a real positive number that accounts for deviations from the late-time value Ω Λ = 1 − Ω m − Ω r .We can also express the combination c 1 /c 2 2 in terms of the initial scalar field value, f R,ini , allowing Equation ( 18) to be rewritten as where R ini ≡ R(a ini ).
Fixed-form parametrisation.Using the EFT approach we use a simple, yet effective, parametrisation.The background is fixed to ΛCDM, i.e. w = −1, and the alpha functions defined as Here Ω Λ (a) is the fractional energy density of the cosmological constant, and the αi 's are free parameters that fix the amplitude of the EFT functions.The idea is to have a ΛCDM universe before the onset of DE, and allow for modifications at late times.Equation ( 20) does not pretend to emulate realistic Horndeski models (with the exception of simple Quintessence), but rather to provide a minimal extension capable of capturing the phenomenology of all the alpha functions.
Given that the background evolution is fixed to be ΛCDM, on top of the standard cosmological parameters the only extra ones are to which we assign αK = 1, αB = 0.2 and αM = 0.1 for our test cosmology.
2.1.Stable parametrisation One challenge when considering the entire class of Horndeski models is their potential for instability in response to perturbations.These instabilities result from unsuitable background solutions, and it is crucial to identify and exclude any parameter combinations leading to such unstable evolution.Physical instabilities come in three kinds (e.g.Hu et al. 2013;Bellini & Sawicki 2014): (i) gradient instabilities occur when the square of the speed of sound for perturbations turns negative as the background evolves.This negative value can cause perturbations to grow exponentially at small scales, destabilizing them over timescales that are comparable to the cutoff of the theory; (ii) ghost instabilities arise when the sign of the kinetic term for background perturbations is incorrect.This issue is often considered in discussions of quantum stability.If ghost modes are present and dynamically interact, they can destabilize the vacuum, resulting in the production of both ghost and normal modes; (iii) tachyonic instabilities manifest when the effective mass squared of scalar field perturbations is negative, leading to power-law instabilities at large scales growing faster than the Hubble scale.While the absence of ghost and gradient instabilities is an integral requirement for a healthy theory, tachyonic instabilities are not necessarily harmful.In fact, imaginary effective masses can still produce observationally viable models, and considering them as pathological a priori would amount to an overly conservative approach (see Zumalacarregui et al. 2016;Frusciante et al. 2018;Gsponer & Noller 2021, for different treatments and interpretations of this instability).Here, we will only impose that a theory must be free from ghost and gradient instabilities, which is expressed by the conditions Another class of instabilities that can plague Horndeski's models are the mathematical (or classical) instabilities (Hu et al. 2014), which manifest as exponentially growing modes in the perturbations.The scalar field fluctuations, V X ≡ aδϕ/ϕ ′ , follow the equation of motion ) where the expressions for the time-dependent coefficients, {A, . . ., E}, are detailed in Appendix D. To prevent significant growth of unstable modes over cosmological timescales, the following conditions must be satisfied: (24) Here, ξ is a parameter controlling the level of instability in units of the Hubble constant, H 0 , with larger values of ξ allowing models to exhibit faster growing rates of these pathological modes.
For simple fixed-form parametrisations of the αfunctions (e.g., Equation 20), or even for theories based on a covariant formulation (e.g., Equations 8-13), finding stable models is just a matter of solving for the background evolution and check that the physical (and mathematical) conditions are satisfied.In these scenarios, we must only vary a handful of parameters and the search process is rather efficient.However, for more complex parametrisations producing a wider variety of functional forms, this trial-and-error approach can rapidly become sub-optimal.A mathematically equivalent, yet computationally more advantageous strategy consists in describing the Horndeski's Lagrangian at the level of background and linear perturbations using a basis that guarantees no-ghost and no-gradient instabilities from the start (Kennedy et al. 2018;Lombriser et al. 2018).In practice, this can be achieved by replacing the α i 's with together with a function describing the evolution of the background, such as the energy-density of the scalar field, ρ ϕ , its equation of state, w ϕ , or the Hubble parameter.The initial condition, α B0 ≡ α B (z = 0), is used to integrate the non-linear differential equation for the braiding s , plotted against the scale factor, a, for a designer f (R) gravity with w = −1 and B 0 = 0.01.The blue solid line illustrates the result derived from the Horndeski's α i functions, the Hubble parameter, and the background densities and pressures (Equation 7).The black dashed line indicates the expected theoretical value of c 2 s = 1.The early-time oscillations in the derived quantity are the consequence of cancellation errors further amplified by the small magnitude of D kin in this model.These rapid variations and, more importantly, the violation of the stability condition c 2 s > 0 erroneously flag this viable model as unstable, motivating the need for an alternative parametrisation.Right: Effect of the QSA on the z = 0 matter power spectrum for our test cosmologies.Each line represents the ratio of the matter power spectrum computed by selectively activating the QSA (see discussion in Section 3.1) to that obtained by following the full dynamics (FD) of the scalar field fluctuations.Solid lines indicate the QSA implemented in hi class, while dashed lines represent the same approach as implemented in MGCAMB.The two QSA strategies are generally consistent, although the methodology based on the metric potentials performs significantly better for the f (R) gravity cosmologies on scales k ≳ 0.1 h/Mpc.
where the source term is with the scalar field pressure, p ϕ , obtained either from inverting the continuity equation (for more details see Section 3.1) or, if w ϕ is provided, we must first integrate the continuity equation (see Appendix D) and then compute , and α M is given by Equation (3).The Hubble parameter is obtained from the first Friedmann equation, H = ρ tot , where the total energy-density includes the contribution from the scalar field.Similarly, the second Friedmann equation gives Ḣ = −3(ρ tot +p tot )/2H.Note that Equation (27) follows from Equation ( 7) after expressing H ′ in terms of the density and pressure of the matter and scalar fields, and upon transforming conformal time derivatives into derivatives with respect to ln a. Once a solution for α B is found, the kineticity, α K , can be derived from Equation ( 6).
It has been argued that stable basis functions offer no distinct advantages over the standard α-functions in terms of computational efficiency or modeling capabilities (e.g., Denissenya & Linder 2018).Contrary to this, our work suggests that reformulating Horndeski's gravity using the stable basis, Equation ( 26), enhances our ability to efficiently sample viable theories and improves numerical stability.The left panel of Figure 1 illustrates how parametrisation with the α-functions might misclassify models as unstable.Specifically, designer f (R) gravity on a ΛCDM background, which is inherently stable (c 2 s = 1), shows rapid variations in the speed of sound due to inexact numerical cancellations when using Equation ( 7), violating the no-gradient condition multiple times.In contrast, the stable basis directly incorporates the speed of sound as an input function, which, in this example, remains constant at unity throughout.Furthermore, using physically motivated parametrisations of the stable basis functions alongside highly optimised root finders and fast integrators (such as those implemented in modern Einstein-Boltzmann solvers) typically leads to rapid converge towards a solution for the braiding with the expected early-time evolution (see Section 5 for details).Depending on the complexity of the parameter space, iteratively solving a differential equation is often more efficient than searching for a stable model starting from the α-functions.

Quasi-static approximation
The quasi-static approximation (QSA) is key to the study of cosmological perturbations within modified gravity and dark energy models.This approximation rests on the premise that the primary timescale of influence is the Hubble rate, and on scales smaller than the sound horizon (c 2 s k 2 ≫ a 2 H 2 ) the time derivatives of the scalar field perturbations can be neglected owing to the rapid response of this degree of freedom to changes in the system.This simplifies the equations of motion for the scalar field into algebraic constraints, thus offering substantial computational benefits.In particular, for large values of the effective mass (i.e. the term proportional to V X in Equation 23) the scalar field perturbations undergo high-frequency oscillations that necessitate small integration steps.These can slow down computations significantly or even cause the solver to stall when the step size falls below machine precision.By employing the QSA, we can focus on the secular evolution of these perturbations, effectively excluding the high-frequency components.
However, the QSA must be applied with care depending on the scale and the specific dynamics of the model.Einstein-Boltzmann solvers like MGCAMB (Hojjati et al. 2011;Zucca et al. 2019b;Wang et al. 2023) and mgclass (Sakr & Martinelli 2021) employ this approximation across all scales and throughout the entire evolution of the perturbations -an approach that can potentially affect the accuracy of the predictions on large scales (Sawicki & Bellini 2015).At the other end of the spectrum are solvers like EFTCAMB (Hu et al. 2013) that avoid using the QSA entirely, a choice that can significantly degrade performance, in particular for certain subclasses of models 2 .The hi class solver (Zumalacarregui et al. 2016;Bellini et al. 2019) adopts a hybrid strategy that leverages the computational efficiency of the QSA along with the accuracy provided by the full scalar field dynamics.Within this scheme, the QSA is selectively activated or deactivated for each time and scale individually, ensuring its use is optimized.The scalar field fluctuations and their time derivatives can be read off Equation ( 23) once the terms proportional to V ′′ X and V ′ X are neglected, It should be noted that, even within the quasi-static regime, V X remains time-dependent, implying that its time derivative does not vanish.Consequently, terms proportional to V ′ X (or time derivatives of other metric perturbations) in the Einstein equations, which are typically not reduced by small coefficients, must be accounted for in our calculations.
Although the QSA implemented through Equation (29) provides a good description of the scalar field fluctuations and their derivatives, it can lead to smallscale inaccuracies in the matter perturbations for some models.This is illustrated in the right panel of Figure 1, where the matter power spectrum for two of our test f (R) gravity models (pink and yellow solid lines), computed by selectively activating the QSA in hi class, deviates by up to 1% for k ≳ 0.1 h/Mpc.For analyses restricted to the linear regime, these deviations are negligible.However, the linear power spectrum is also used as input for methods that predict non-linear structure formation (e.g., Bose & Koyama 2016;Aviles & Cervantes-Cota 2017;Cusin et al. 2017;Cataneo et al. 2018;Giblin et al. 2019;Valogiannis et al. 2019), and such small-scale inaccuracies can propagate to observables probing the non-linear regime, potentially exceeding the accuracy requirements established for Stage IV surveys (e.g., Taylor et al. 2018).
A different QSA strategy consists in working directly with the modified metric potentials (Bloomfield et al. 2012;Pace et al. 2020).In the Newtonian gauge, the linearised modified Einstein field equation for the Newtonian potential, Ψ, and the intrinsic spatial curvature, Φ, read (Zhao et al. 2008;Hojjati et al. 2011) 2 See Hu et al. (2014) for the strategy adopted in EFTCAMB to mitigate the impact of this choice.
where the total comoving-gauge density perturbations are summed over all matter species I ∈ {c, b, r, ν, . . .}3 , ρ m ∆ m ≡ I ρ I ∆ I , with ∆ I being the comoving density contrast for the species I, and (ρ m + p m )σ m ≡ I (ρ I + p I )σ I with σ I denoting the anisotropic stress of the individual matter fluids.Changes to the gravitational coupling experienced by non-relativistic particles are encoded in µ, while the gravitational slip, γ, sources differences between the two metric potentials also in the absence of anisotropic stress.General Relativity is recovered for µ = γ = 1.
For the Horndeski's models described by the Lagrangian in Equation ( 5) the two phenomenological functions take the simple forms (Pogosian & Silvestri 2016) consistent with the notation of Pace et al. (2020).Definitions for the quantities {µ p , µ ∞ , µ Z,∞ } are detailed in Appendix D. The extension of hi class introduced in this work, mochi class, bypasses the QSA formulation from Equation ( 29) in favor of Equations ( 30) and ( 31) in synchronous gauge, similar to MGCAMB4 .However, it activates the QSA selectively based on scale and time, when the scalar field's effective mass and the discrepancy between full and quasi-static solutions meet specific user-defined criteria (see Bellini et al. 2019, and Section 3 for details).The dashed lines in the right panel of Figure 1 demonstrate that the QSA approach implemented through the metric potentials is less affected by the neglected dynamics of the scalar field fluctuations compared to the implementation using Equation (29).

CODE STRUCTURE AND IMPLEMENTATION
The code developed for this study, mochi class, builds on hi class, which itself extends the capabilities of the Einstein-Boltzmann solver class (Blas et al. 2011;Lesgourgues & Tram 2011) to evolve the background and linear perturbations of scalar-tensor theories within the Horndeski's framework.This section explores the new features introduced in mochi class.For a comprehensive review of the functionalities in hi class we refer the reader to Zumalacarregui et al. (2016); Bellini et al. (2019).
mochi class enhances hi class by modifying the input, background, and perturbations modules to ensure greater numerical stability and facilitate a more comprehensive investigation of Horndeski's gravity.In the following, we provide a detailed description of these modifications.

Input and precision parameters
By setting the gravity model variable to stable params, users can specify a path to a file containing the functions ∆M 2 * , D kin , and c 2 s , as well as the ln a sampling times, using smg file name.Alternatively, these quantities can be directly assigned to the corresponding code variables: Delta M2, D kin, cs2, and lna smg.This direct assignment method is particularly recommended when running mochi class via the classy Python wrapper.To mitigate the impact of interpolation errors within the computational pipeline, we recommend using a time-sampling step size no larger than ∆ ln a = 7.5 × 10 −4 for the range ln a ∈ [−2, 0].In mochi class, the hi class variable parameters smg serves a container for the initial condition, α B0 .
Users must also select an expansion model to define the background evolution of the scalar field.They can choose from the two hi class analytic expansion models, {lcdm, w0wa}, or opt for the new non-parametric models {wext, rho de}.In the wext model, the scalar field background density is integrated using the continuity equation (Equation D1), with the background pressure defined as p ϕ = w ϕ ρ ϕ .For the rho de model, the approach starts with the normalised background density ρϕ ≡ ρ ϕ /ρ ϕ (z = 0), which is transformed to ρ ϕ = Ω ϕ H 2 0 ρϕ .Here, the Hubble constant, H 0 , is an input parameter, and Ω ϕ = 1 − I Ω I represents the present-day normalised background energy-density in a flat universe.The scalar field background pressure in this configuration is then computed as It is important to note that hi class employs a rootfinding algorithm to determine the value of a user specified parameter necessary to achieve the present-day energy density, Ω ϕ H 2 0 (see Section 3.3 in Bellini et al. 2019).It can be the initial conditions of the scalar field, ϕ(τ ini ), its time derivative, ϕ ′ (τ ini ), or any other parameter that contributes to the energy-density of the scalar field.This step is essential for models based on a covariant Lagrangian where ρ ϕ is not known in advance, requiring the cosmological background to be solved self-consistently with the scalar field and the α-functions.However, this becomes redundant in cases where the scalar field energydensity evolution is predefined, as often happens in the effective-theory approach.To optimize computational efficiency, mochi class has been modified to bypass the root-finding process.Instead, it directly uses Ω ϕ from the closure relation, Ω tot = 1, to rescale the normalized energy density, ρϕ .
Similarly to the stable basis functions, wext and rho de can read the necessary data either from a text file containing two columns {ln a, w ϕ or ρϕ } located at a path specified in expansion file name, or directly from input arrays lna de and de evo.To minimise interpolation errors, it is recommended to use a time-sampling step-size ∆ ln a ≤ 7.5 × 10 −4 for the range ln a ∈ [−2, 0].In addition, when utilizing mochi class, users must activate the GR approximation scheme (see Section 3.3 below) by setting method gr smg = on.This feature allows the code to solve the standard system of equations down to a user-specified redshift, z gr smg.In practice, for redshifts greater than z gr smg, the scalar field behaves like a cosmological constant with ρ Λ = ρ ϕ (z gr smg), and vanishing α i 's (see Equation 4).Given our focus on late-time modifications to GR, the default value for this transition epoch is set deep within the matter-dominated era, specifically at z gr smg = 99.
While the stable basis in mochi class ensures the positivity of the de-mixed kinetic term, the Planck mass, and the speed of sound, it does not automatically guarantee that the model is free from the α B = 2 discontinuity (see Equation A.14 in Zumalacarregui et al. (2016)), as this requires a solution for the braiding first.Therefore, we can safely set skip stability tests smg = yes, as mochi class always tests that the braiding never crosses 2, which was not implemented in hi class, and terminates returning a computation error5 .Furthermore, setting skip math stability smg = no allows the code to check for exponentially growing modes by enforcing the mathematical stability conditions outlined in Equations ( 24) and ( 25).Users can control the allowed instability growth rate via the exp rate smg variable, which defaults to 1.
Similar to hi class, mochi class has the option to treat the scalar field as a dynamical degree of freedom throughout the perturbations' evolution, or to apply the QSA whenever is safe to do so6 .These configurations can be selected by setting the method qs smg input variable to fully dynamic for the former, or automatic for the latter.In automatic mode, the QSA is activated for specific time intervals and wavenumbers based on whether the effective mass of the scalar field exceeds the threshold trigger mass qs smg, and the amplitude of scalar field fluctuations has reduced by a factor defined in eps s qs smg (Bellini et al. 2019) 7 .Default values are set at trigger mass qs smg = 100 and eps s qs smg = 0.01, which, even for large departures from the standard cosmology, affect the computed power spectra by ≲ 0.1%.mochi class can accurately solve the quasi-static equations down to z = 0 (see Section 2.2), eliminating the need to enforce the full dynamic of the scalar field at low redshifts.Consequently, the corresponding control variable z fd qs smg can be safely set to 0.
mochi class also includes new precision parameters.Below is a list describing their scope and default values: • background Nloga smg and loga split: the background time-sampling strategy adopted for a ≳ 0.05 is crucial for achieving accurate interpola-tion of the quasi-static quantities {µ p , µ ∞ , µ Z,∞ }.
To this end, mochi class employs a denser sampling for ln a > loga split, with the number of points controlled by background Nloga smg.
For ln a < loga split the number of sampling points reverts to the standard set by the class precision variable background Nloga.The default settings are background Nloga = 3000, background Nloga smg = 3000, and loga split = -3.
• eps bw integration rho smg: when using the expansion model = wext, where users define the equation of state for the scalar field, its background density is calculated from the continuity equation.This equation is integrated from a = 1 down to the scale factor (1 − eps bw integration rho smg)/(1 + z gr smg)8 .To ensure smooth derivatives of the scalar field pressure and related quantities at z gr smg, the default setting for eps bw integration rho smg is 0.1.
• eps bw integration braid: the equation for the braiding, Equation ( 27), is integrated from a = 1 backward to the scale factor (1 − eps bw integration braid)/(1 + z gr smg).To ensure smoothness at z gr smg for all functions derived from α B and its derivatives, we recommend setting eps bw integration braid to at least 0.1.
• tol background bw integration: this variable governs the integration accuracy of the continuity equation, Equation (D1), and of the braiding equation, Equation ( 27).Based on our analysis, setting tol background bw integration = 1e-10 provides a good balance between computational speed and solution accuracy across all models examined in this study.
• braid activation threshold: Solving Equation (27) into the matter-dominated epoch poses numerical challenges, particularly for models that rapidly converge to GR at early times.While inaccuracies at high redshifts generally do not impact the late-time modified gravity phenomenology, they can introduce spurious instabilities that may incorrectly render a model non-viable.To minimize these effects, modifications to the Einstein equations are triggered only when the braiding solution reaches the braid activation threshold.Practically, this means adjusting the transition redshift, z gr smg, to match the redshift at which |α B | first equals braid activation threshold.Note that this threshold is not enforced for models with |α B0 | ≤ braid activation threshold, e.g.quintessence and k-essence9 .The default value of braid activation threshold = 1e-12 is chosen to ensure precise predictions across the broad range of deviations from GR examined in this study.
• dtau start qs: at the onset of modified gravity, marked by z gr smg, mochi class introduces the scalar field as a new degree of freedom by incorporating the equation for the field fluctuations (Equation 23) and its couplings to the metric perturbations into the Einstein-Boltzmann system.Initial conditions are defined by the QSA outlined in Equation ( 29).The code then waits a conformal time specified by dtau start qs before assessing whether to activate the QSA, effectively turning the scalar field fluctuations into algebraic constraints.We found that dtau start qs = 1e-7 works well for all models tested.
3.2.Background Once the input arrays are prepared, the initial step involves calculating the background density of the scalar field.Depending on the chosen expansion model, the code performs one of the following actions: (i) integrates the continuity equation and stores both ρ ϕ and p ϕ for later use; (ii) rescales the input normalized background energy-density of the scalar field to compute ρ ϕ , and calculates p ϕ using Equation (32), saving both values for future reference; (iii) analytically determines the scalar field density and pressure, along with other background quantities, as already implemented in hi class.
Subsequently, with access to the stable basis functions, the scalar field density and pressure, as well as H and Ḣ, the code solves for α B by backward integrating Equation ( 27) from its present-day value, α B0 .It calculates the kineticity, α K , using Equation (6), and if necessary updates z gr smg based on the braid activation threshold.However, it is important to acknowledge a slight inconsistency in our approach: ρ ϕ , p ϕ , and associated variables are not adjusted to emulate a cosmological constant over the interval between the updated and the original z gr smg values.This methodological choice is made to save computational resources while accounting for the very small effect of the modified background on the perturbations.Within the Horndeski's framework we expect changes to the growth of structure to be linked to changes to the background expansion.Therefore, viable models exhibiting significant departures from standard gravity only in recent epochs, presumably also closely match a ΛCDM expansion history until recently.
Finally, with all relevant background functions now available, mochi class computes the quasi-static quantities {µ p , µ ∞ , µ Z,∞ } and their derivatives.These are also stored in the background table and forwarded to the perturbations module for further processing.

Perturbations
The evolution of the perturbations closely follows that implemented in hi class, except for the new GR approximation scheme and the effective field equations used for introducing unwanted inaccuracies in the evolution of perturbations in these modified gravity models, we recommend setting braid activation threshold to 0, ensuring that z gr smg remains fixed at the original value specified by the user.
the QSA.While hi class follows the dynamics of the scalar field and its coupling to the metric tensor already prior to recombination, mochi class activates the additional degree of freedom only after the time z gr smg.This is achieved by enabling method gr smg, which directs the solver to apply the standard class equations up to the transition redshift.
Once the GR approximation is deactivated, the initial conditions for the scalar field are derived from the quasistatic expressions given in Equation ( 29), transitioning from the standard class equations to those of hi class.This approximation scheme is crucial to manage the discontinuity at z gr smg, which otherwise leads to computational instability by forcing the integrator to attempt reductions in the step size to below double precision levels.Consistent with other approximation schemes used, the initial conditions for the metric and stress-energy tensor perturbations are set from the integration time step immediately preceding the deactivation.
With the initial conditions now defined, mochi class can use the QSA scheme of hi class for the evolution of the perturbations.Should the scalar field effective mass and fluctuations amplitude satisfy specific criteria (Bellini et al. 2019), the QSA may be selectively activated or deactivated up to six times.However, within mochi class, this scheme is only applicable to each wavenumber, k, following a designated time, τ (z gr smg) + dtau start qs.Furthermore, rather than applying the QSA directly to the scalar field fluctuations, V X , mochi class targets the metric potentials, thereby enhancing numerical accuracy.Accordingly, the quasistatic part of hi class has been restructured to align with the methodology developed for MGCAMB.

CODE VALIDATION AND PERFORMANCE
Before exploring the new capabilities of mochi class, we first compare its predictions with those from the established codes hi class (Zumalacarregui et al. 2016;Bellini et al. 2019), EFTCAMB (Hu et al. 2013(Hu et al. , 2014)), and MGCAMB (Hojjati et al. 2011;Zucca et al. 2019b;Wang et al. 2023).This comparative analysis allows for a thorough assessment of all test models and ensures the robustness of our results across different parametrisations, approximations, and integration strategies.Although not discussed here, similar agreements between mochi class and the other Einstein-Boltzmann solvers are achieved when including massive neutrinos.This is especially true for the small, non-factorisable contributions to the growth of structure caused by their coupling to modified gravity.

CMB spectra
Figure 2 displays the CMB multipoles for temperature (upper left), polarization (upper right), their crosscorrelation (lower right), and the lensing potential power spectrum (lower left).The lower sub-panels contain the relative deviations of the mochi class predictions (coloured lines) compared to the results from other Einstein-Boltzmann solvers (black lines).As reference spectra, we use hi class without the QSA for cubic Galileon, nKGB, and propto omega models; EFTCAMB for designer f (R) gravity; and MGCAMB for Hu-Sawicki f (R) gravity.These quantities were generated using the preci-sion settings detailed in Appendix A for both CAMB and class.
First, it is important to note that the differences between mochi class and hi class are barely visible at this scale.This suggests that the newly implemented GR approximation at early times, the reparametrisation using stable basis functions, and the application of the quasi-static approximation until late in the perturbation evolution, all introduce minimal numerical errors.
The agreement between mochi class and the CAMBbased codes is very good, with differences remaining well within ∼ 0.25%.Further analysis reveals that mismatches exceeding 0.05% come from variations between class and CAMB, similarly affecting ΛCDM predictions (not shown).This validates the QSA implementation in mochi class, which closely aligns with MGCAMB, and confirms the reliability of the QSA/full dynamics switching feature in hi class.Note that despite EFTCAMB and MGCAMB both deriving from CAMB, discrepancies around 0.2% are observed between these codes for the low-ℓ temperature multipoles.Initially, these could be attributed to the activation of the QSA in mochi class, potentially affecting the accuracy of the integrated Sachs-Wolfe (ISW) effect calculations, particularly in the designer f (R) model where deviations from GR are more pronounced.However, further investigation indicated that these differences are due to the use of different CAMB versions, with MGCAMB employing a more recent release, and also persist for the standard cosmology (not shown).

Matter power spectrum
Figure 3 presents the evolution of the matter power spectrum ratio, P MG m (k)/P GR m (k), for the three cosmologies characterised by scale-independent growth on subhorizon scales.The lower sub-panels show that for these models mochi class (coloured lines) deviates by an average of approximately 0.01% from hi class (black lines).Differences are negligible for scales k ≲ 0.01 h/Mpc, where the QSA is never activated.
To isolate the impact of the new mochi class features, Figure 4 also shows the matter power spectrum of the f (R) gravity predictions relative to same quantity in ΛCDM.For the Hu-Sawicki models the agreement between mochi class (coloured lines) and MGCAMB (dotted lines) is excellent across all redshifts.Differences larger that 0.2% are visible only for k ≳ 1 h/Mpc, regardless of model and time.Using an independent Mathematica notebook employing the QSA, we computed the linear growth functions both in HS f (R) gravity and in ΛCDM (see, e.g., Eq. 81 in Brax & Valageas 2013), and by taking their ratio we verified that these small differences originate from MGCAMB.
Moreover, Figure 4 illustrates the (albeit small) inaccuracies introduced by the QSA in the designer model.The blue lines, generated using a hybrid method that activates the QSA under certain conditions (detailed in Section 3), should be compared with the EFTCAMB predictions (dashed lines), which do not employ the QSA.As expected, the small scales are most affected by this approximation, due to the conditions for its activation being met for more extended time intervals.To verify that deviations from EFTCAMB are caused by the QSA, we forced mochi class to follow the full dynamics of the scalar field (green line).The lower sub-panels confirm Figure 2. Comparison of the CMB auto-and cross-power spectra as computed by mochi class (colored lines) against predictions from well-validated Einstein-Boltzmann solvers (black lines).The dashed lines for cubic Galileon and nKGB models, as well as propto omega, are generated by hi class following the full scalar field dynamics.For designer f (R) gravity, the dot-dashed lines are obtained with EFTCAMB, while the dotted lines for Hu-Sawicki (HS) f (R) models are produced by MGCAMB using the quasi-static approximation.The fractional differences, ϵ, shown in the lower sub-panels quantify the agreement between mochi class in automatic mode and the reference spectra obtained with the other Einstein-Boltzmann solvers.Notably, all discrepancies are below 0.25% except in the TE cross-spectra, where the reference spectra's zero crossings amplify the deviations.All discrepancies larger than 0.05% for the f (R) models stem primarily from differences between class and CAMB, and are also present in ΛCDM (not shown).
that mochi class predictions now align with EFTCAMB at a level of 10 −4 .

Performance
Here we assess the impact of the mochi class patch on execution time, examining both individual modules and the overall runtime.Table 1 summarises the relative changes in execution time compared to hi class for the input, background, and perturbations modules, as well as their combined effect and the total runtime variation, including the other unmodified modules.To ensure a fair comparison, we evaluate two models that activate different features in hi class: the cubic Galileon model, which requires preparatory steps to iteratively determine suitable initial conditions and map the Horndeski G i 's to the α-basis, and the propto omega model, which involves direct analytic expressions for both the α's and the scalar field background evolution.
Contrary to intuition, mochi class shows reduced runtime for both models at the input level, despite propto omega only needing to fill the relevant arrays.This is because hi class unnecessarily applies the root finder for the initial conditions to this model-a step not performed by mochi class.Unsurprisingly, the execution of the background module takes twice as long as in hi class, due to the integration of the differential equation for the braiding function.However, this slowdown is less severe than suggested by Denissenya & Linder (2018).Assuming a uniform scanning strategy of their two-dimensional parameter space and a stable region amounting to 22% of the prior volume, the number of unstable models is 3.5 times that of the viable models.Therefore, starting from the stable basis ensures an overall 56% speedup for the background module compared to using the α-functions in this particular scenario.This efficiency gap becomes even wider for complex parametrisations, such as those discussed in Section 5, where stable models can be particularly challenging to find.
The third column of Table 1 highlights the advantage of enabling the QSA for z < 10: it can reduce the com- Figure 3. mochi class accuracy quantification for the matter power spectrum through comparison with hi class across various redshifts (z = 0, 0.5, 1, 2).Each main panel displays the matter power spectrum ratio in modified gravity models relative to ΛCDM as calculated by mochi class (colored lines) and hi class (dashed lines).The sub-panels show the fractional deviation of mochi class from hi class, indicating an agreement within 0.05% for the considered redshifts.This level of accuracy demonstrates the robustness of the computational strategy adopted in mochi class, which includes an alternative parametrisation of the Horndeski's functions, a QSA formulated in terms of modifications to the equations for the metric potentials, and the activation of modified gravity deep into the matter-dominated era.
puting time for the evolution of perturbations by up to 50%.Since this part of the code contributes most significantly to the cumulative execution time of the three modified modules (fourth column), any performance improvement here directly impacts the total runtime of the Einstein-Boltzmann solver (last column).The difference between the two codes arises from enabling the activation of the QSA in mochi class throughout the perturbation evolution, while in hi class we apply it only up to some high redshift (i.e., z = 10), as recommended by its developers.
5. EXTENDED PARAMETRISATION FOR HORNDESKI GRAVITY After establishing mochi class accuracy and efficiency using extensively studied modified gravity cosmologies, we can now leverage its new features to input general functions of time and explore the phenomenology of Horndeski gravity in greater detail.For this purpose, a parametrisation that generalises propto omega and includes all covariant theories analysed in Section 3 could be extremely valuable (Linder 2016).
A common feature of all our test models is that, deep in the matter-dominated era, their stable basis functions either approach the GR limit as power laws, A J a ζ J , or tend to constant values, C J , to a very good approximation.For instance, Figure 5 illustrates that, regardless of the specifics governing the dynamics of the scalar field, the de-mixed kinetic term is well described by a power law for a ≲ 0.05 (shaded area).Therefore, any late-time evolution of the stable functions can be interpreted as deviations, ∆ J , from this early-time power-law behavior.Specifically, we can express the stable functions as where x ≡ ln a, b J ≡ ln |A J |, sgn denotes the sign func-   The reference power spectra for f (R) gravity are computed with EFTCAMB (dashed) and MGCAMB (dotted).The lower panels focus on the fractional deviation of mochi class from the CAMB-based solvers, providing a direct comparison for the matter power spectrum ratios shown above.The Hu-Sawicki (HS) models demonstrate a close match, with deviations from MGCAMB remaining under 0.025% across all scales.In contrast, for the designer f (R) model, mochi class shows a consistent discrepancy that increases at smaller scales (blue line), a result predominantly due to the quasi-static approximation.This difference, in fact, significantly reduces when mochi class always integrates the full scalar field dynamics (green line).
tion, and for x ≪ 0 we have the limiting values The constant C ρ is not a free parameter; it is determined by the constraint ρϕ (x = 0) = 1.For negative A M and A cs , we must ensure that ∆M 2 * > −1 and c 2 s > 0 throughout the evolution.Given that the power-law term for the speed of sound generally represents a small correction to the constant C cs (as discussed below), it can be neglected when generating new models.Therefore, the gradient-free condition implies simply that C cs > 0.
5.1.Parametrising departures from constants and power laws Next, we need to find a sufficiently general yet simple parametrisation for the functions ∆ J .A commonly employed strategy to parametrise the EFT functions or the α i 's involves expanding them as ratios of polynomials, i.e., Padé approximants (see, e.g., Gleyzes 2017;Peirone et al. 2017).Here, we opt for an alternative methodology grounded in Gaussian Processes (GPs) (see, e.g., Rasmussen & Williams 2006;Aigrain & Foreman-Mackey 2022, for detailed discussions) and Principal Component Analysis (PCA) applied to the deviations ∆ J .In broad terms, we proceed as follows: 1. We generate training data such that the GP is consistent with Equation (34).Practically, this requires an array of zeroes over the interval x ∈ [−5, −3].
2. After choosing a GP kernel, K(x i , x j ), for a particular basis function, and values for its hyperparameters, we let the GP learn from the training data.
After the training, the GP posterior collapses to a delta function centered around zero over the range x ∈ [−5, −3], while the prior increasingly dominates the more we move into the untrained region.Note that all our GPs have zero mean.
3. We draw thousands of samples from the trained

Table 1
Variations in execution time for the input (I), background (B), and perturbations (P) modules in mochi class, relative to hi class.The values in parentheses are the mochi class's execution times in seconds averaged over 10 single-thread runs on an Apple M1 Pro chip.Both codes used the automatic mode, with hi class forcing full dynamical evolution of the scalar field from z = 10, as recommended by its developers.The models considered-Cubic Galileon, based on the covariant Lagrangian, and propto omega, based on effective theory-are representative of the two approaches implemented in hi class.The column labeled (I)+(B)+(P) lists the cumulative variations across the three modules, emphasising that the evolution of the perturbations is the primary bottleneck.The final column indicates that, overall, enabling the QSA at low redshifts can lead to substantial performance improvements, with minimal impact on accuracy (see Figure 1).For a detailed comparison between the execution time with and without QSA see Fig. 11 in Bellini et al. (2019)., across a variety of modified gravity theories.At sufficiently early time (shaded area) D kin closely follows a power-law behavior for all models.A similar behavior (not shown) characterises departures from the cosmological strength of gravity, ∆M 2 * , and the background energy density of the scalar field, ρ ϕ .
GP, calculate their mean, and subtract it from each sample so that the transformed samples have zero mean.
4. We perform PCA on the transformed samples and retain the first N J principal components.
5. We project the specific de-trended basis function (i.e.we first remove the constant and/or power-law contributions in Equation 33) of each test model (cubic Galileon, nKGB, etc.) onto the selected principal components and quantify the accuracy of the reconstructed approximation.
6.For each basis function, we independently iterate over steps 2-5 to identify the kernel, hyperparameter values, and number of principal components that achieve percent-level reconstruction accuracy for x > −3 across all test models simultaneously.
It is reasonable to assume that a covariant theory of gravity in the Horndeski class should be described by basis functions that vary smoothly with time.This assumption restricts the choice of the GP kernel to either the squared exponential form, or the rational quadratic form, Here, the variance σ 2 represents the average distance of the generated samples away from the mean, the lengthscale ℓ determines the oscillatory features of the samples, and the parameter ϑ in the rational quadratic kernel controls the relative weighting of large-scale and small-scale variations10 .Intuitively, we expect that the most relevant hyperparameters controlling the functional forms generated by the GP are ℓ and ϑ, as σ 2 acts as a simple rescaling factor11 .Therefore, we fix σ 2 = 10, and only vary ℓ and ϑ together with the number of principal components to achieve our target reconstruction accuracy of 1%.As an example, the left panel of Figure 6 displays random functions, ∆ ρ , sampled from a trained GP using the rational quadratic kernel, as described by Equation ( 36), with ℓ = 1 and ϑ = 5.The right panel of the same figure shows the first six principal components derived from a total set of 10,000 samples, collectively accounting for 90% of the cumulative explained variance, which meets the desired accuracy threshold12 .Similar considerations apply to the other basis functions; Table 2 details the kernel, hyperparameter values, and number of principal components used for each.To keep interpolation errors in mochi class under control, our time-sampling strategy for the random curves (and consequently for the principal components) closely follows that described in Section 3.1.Specifically, we use 3,000 sampling points for the range [−3, 0] and 200 sampling points for the range [−5, −3]. Figure 6.Left: Samples from a Gaussian Process (GP) representing the correction factor, ∆ρ, which quantifies deviations in the scalar field background energy-density from a baseline power-law behavior as a function of time, x ≡ ln a.The GP's hyper-parameters are detailed in the main text, and its mean is set to 1.The shaded region marks the training range, over which the GP is constrained to unity, coinciding with the interval where the background energy-density of the scalar field closely follows a power-law behaviour.The various trajectories of the 20 curves outside the training range represent possible extrapolations of ∆ρ beyond the power-law regime.Right: The first six principal components (PCs) extracted from 10,000 samples generated by the Gaussian Process described in the left panel.The PCs illustrated here are sufficient to reconstruct the scalar field background energy-density of all test models to better than 1%, thus demonstrating the efficiency of the dimensionality reduction achieved by the principal component analysis (PCA).Similar performance apply to PCAs for the other input functions within the stable basis.), and the right panel showing the normalized absolute deviations from the standard cosmological strength of gravity, |∆M 2 * /∆M 2 * ,max |.Coloured lines indicate the functions as reconstructed from the principal components and power-law approximation, while the dashed lines are the exact quantities.The lower sub-panels quantify the reconstruction accuracy, which is consistently better than 1% for the entire range where modified gravity significantly affect the growth of structure a > 0.1.

Reconstruction of test models
Table 2 Kernel form (rational quadratic (RQ) or squared exponential (SE)), corresponding hyperparameters, and number of retained principal components (# PCs) used to model each of the deviation functions, ∆ J , in Equation ( 33).
Figures 7 to 9 illustrate the reconstruction accuracy of the principal components combined with the power law and/or constant approximation for the four basis functions in all of our test models.To project the exact functions onto the principal components, we first compute the natural logarithm of their absolute values13 .We then remove the power-law and/or constant trends from this quantity and finally apply the projection.For instance, to reconstruct the de-mixed kinetic term for any particular model, we do the following: , where the w i 's are the N D projection weights and ψ i is the i-th principal component.
5. Define the approximated function using Equation (33) by replacing ∆ D with ∆D .
From the lower panels of Figures 7-9 we observe that our PCA reconstruction matches the exact stable functions to within 1% for a ≳ 0.05.In a few cases, however, the relative difference, ϵ, exhibits a peak exceeding the target accuracy around a = 0.1.Due to negligible departures from ΛCDM at those redshifts, this is not concerning and has no measurable effect on the cosmological observables.For comparison, Figure 14 in Appendix B demonstrates that when using the same number of free parameters, alternative parametrisations of the ∆ J functions rooted in polynomial expansion underperform compared to the GP+PCA strategy adopted here.
To test how well the reconstructed stable functions describe the original MG cosmologies, we can use them as input to mochi class to predict the background expansion and linear perturbations.To assess the performance of the background reconstruction, we examine the Hubble parameter and conclude that the reconstructed expansion history resembles the exact one to better than 0.2% for all models.However, having the full set of stable functions is not sufficient to compute the perturbations, as we also need the initial condition, α B0 , for the integration of the braiding parameter.We determine this by minimizing the mean squared error where N k is the total number of wavenumbers, k, in the interval [10 −4 , 10] h/Mpc, and P rec.m (k|α B0 ) and P exact m (k) are the present-day matter power spectra of the reconstructed and exact model, respectively.The derived bestfit initial conditions are typically within 0.1-1% of their exact counterparts.Due to the high sensitivity of f (R) gravity to small inaccuracies in the reconstructed functions (mainly ∆M 2 * and ρϕ ), we could not find a satisfactory α B0 that produced stable perturbations for the Hu-Sawicki n = 4 model.We do not reconstruct trivial stable functions: the speed of sound in f (R) gravity is fixed at 1, ∆M 2 * = 0 for the cubic Galileon and nKGB models, and ρϕ = 1 for propto omega.
Figure 10 confirms that by combining GPs and PCA we can predict the secondary CMB anisotropies (i.e.lensing and ISW effect) sourced by the late-time modified growth of the test models to precision levels well within the expected uncertainties of CMB-S4 experiments (see, e.g., Ade et al. 2018).Matter power spectra for cubic Galileon, nKGB and propto omega models are also reconstructed to ≲ 0.05% precision across all redshifts, as shown in Figure 11.On the other hand, even though the reconstruction accuracy of the stable functions for f (R) gravity never exceeds 1%, Figure 12 highlights that the derived matter power spectra deviate from the truth by ≲ 0.2% only up to k ≈ 1 h/Mpc, with the agreement degrading to 1% at k ≈ 8 h/Mpc.Despite these relatively large differences on small scales, uncertainties in the nonlinear evolution (Winther et al. 2015) and baryonic feedback (Chisari et al. 2019) have a much greater impact in this regime, effectively making the reconstructed models observationally indistiguishable from the true f (R) gravity theories.

Exploring Horndeski's landscape
By combining mochi class with the newly introduced parametrisation, Equation (33), we can venture beyond the phenomenological boundaries defined by the ΛCDM extensions discussed in Section 2. For this first exploratory task, we focus on two families of models: 1. Scalar Horndeski (SH): here, we exploit the full functional freedom provided by Equation ( 33) and use the GP+PCA basis presented in Section 5.1.
For simplicity, we fix the number of principal components to six for all ∆ J 's, and given the small contribution of the power-law term to the speed of sound (see Figure 8) we set A cs = 0. Therefore, we have a total of 31 parameters accounting for the power laws and PC weights.
2. f (R)-like (FR): this represents a particular subspace of the SH cosmologies, where sampling can be extremely challenging without imposing specific relations for the stable functions.The freedom is limited to ∆M 2 * and ρϕ , with the remaining functions being c 2 s = 1 and D kin = 3α 2 M /2.Models with α K ≪ α 2 B closely resemble f (R) gravity.To describe this sub-class we use six PC weights and two power-law parameters for each stable basis function.
In both cases we sample the parameter space using a Latin hypercube strategy with 5,000 sampling points for Figure 9.Comparison between the reconstructed background energy-density of the scalar field (coloured lines) and the exact calculations (dashed lines).The left panel shows the evolution of the energy-density for two models minimally coupled to gravity, whereas the right panel displays ρ ϕ for two variants of the Hu-Sawicki (HS) f (R) gravity model which feature a conformal coupling to gravity.For both types of models, the reconstructed quantities closely align with the exact references, maintaining the accuracy within 0.5% for scale factors a > 0.02.
SH and 1,500 for FR.The parameters ranges defining the hypercube can be found in Appendix C.
As with the reconstruction of the test models, we need to choose a value for the initial condition, α B0 .While in the former case we could use the exact calculations as reference and minimise Equation ( 37), here we are generating new models, and α B0 should be interpreted as an additional free parameter.However, if we want the early-time evolution of the braiding to be driven by the stable functions, we can impose an additional constraint to find suitable initial conditions.To see this, we start by linearising Equation ( 27) around α B = 0.This approximation is valid for times x ≤ x 0 , with x 0 sufficiently deep in the matter-dominated era, when we expect α B → 0. By applying the integrating factor technique, we can express the solution to the linearised non-homogeneous differential equation as where represents the time-varying part of the homogeneous solution.The function V contributes to the particular solution that depends on the source term, S. It is calculated by integrating the following differential equation with the associated initial condition: Note that this equation contains no information about α B0 .In contrast, α B (x 0 ) does depend on the value of the braiding at x = 0.In the regime of interest here, beside applying the approximations in Equation ( 34), we can also neglect radiation by setting Ω r = 0.This allows us to simplify the stable functions as power laws plus constants, write ( H /H0) 2 ≈ Ω m e −3x , and approximate M 2 * (x) ≈ M 2 * (x 0 ) ≈ 1. Equation ( 40) can then be integrated analytically for all times x ≤ x 0 , and its solution consists of two parts: a constant term, v, that depends on x 0 , and a sum of power laws.In the limit x ≪ x 0 , the power laws vanish, leaving v as the only remaining contribution.
Now, let us assume a value α B0 such that α B (x 0 ) is either larger or smaller than −v.The early-time evolution of the approximate solution, Equation (38), will be controlled by the homogeneous solution, Equation ( 39), which ultimately scales as ∼ ± √ e x .This scaling is completely independent of the cosmological parameters and late-time dynamics defined by the stable functions.To avoid this spurious behaviour, we must adjust α B0 such that α B (x 0 ) = −v, which yields the following early-time solution: As a case in point, we examine a designer f (R) gravity cosmology with a ΛCDM background and B 0 ≪ 1.For this model, Equation (41) simplifies significantly, as we have c 2 s = 1, D kin = 3α 2 M /2 ≪ α M , and ρϕ = 1.From the definition of the running, Equation (3), and by imposing the condition α B = −α M , we derive a quadratic equation for the slope of ∆M 2 * , with the positive solution being ζ M = (5 + √ 73)/4.Reassuringly, this result matches the value derived by Pogosian & Silvestri (2007) in the limit Ω r → 0 (see also Lombriser et al. 2018).
Finally, for each set of principal component weights, constants, power-law slopes, and amplitudes, we compute the stable functions and use Equation ( 41) as an early-time target for the integrated braiding solution to obtain α B0 .In practice, we find the zero of the function defined as the definite integral over the range x ∈ [−5, −4] of the difference α Figure 10.Comparison of the CMB auto-and spectra obtained with mochi class using the reconstructed stable basis (dashed lines) against the exact predictions (coloured lines).For the cubic Galileon and nKGB models, as well as propto omega, the reference spectra are generated by hi class following the full scalar field dynamics (FD).The exact calculations for the two f (R) gravity models are performed using mochi class with the automatic method and precisely determined stable basis functions.For all reconstructed models, the initial conditions for the braiding parameter, α B0 , are optimised to minimise the matter power spectrum's relative deviation at z = 0 compared to the exact models (refer to Equation 37).The small fractional differences, ϵ, shown in the lower sub-panels indicate that the reconstructed models are observationally indistinguishable from their exact counterparts.
is the solution to Equation ( 27).For this step, we run only the background module of mochi class and search for α B0 in the interval [0, 2]14 .The root finder typically takes 6-7 iterations to reach convergence.Consequently, this process slightly degrades the performance presented in Section 4.3 since more time must be spent on the background calculations.
The left panel of Figure 13 shows the present-day matter power spectrum ratio for 12 models selected from approximately 200 stable cosmologies generated using the SH Latin hypercube.A viability rate of 4% might seem surprising given that all 5,000 sample models pass the stability criteria, Equation ( 22).However, only a small fraction of these models feature stable perturbations not affected by exponentially growing modes, high-lighting the importance of starting from the stable basis functions-using the α i 's to sample this parameter space would be even more challenging.Similar arguments apply to the right panel, where we find a higher stability rate of 20%, owing to the considerably smaller parameter volume.
Returning to the left panel, the variety of shapes for the matter power spectrum ratios reflects the rich phenomenology of the SH landscape.Besides the familiar behavior displayed by models such as cubic Galileon, nKGB, and propto omega, the growth of structure can be slower than in ΛCDM, and can also exhibit non-trivial scale dependence arising from the interplay between the braiding scale, the Compton scale, and the scalar field sound horizon (see, e.g., Sawicki et al. 2012;Bellini & Sawicki 2014, for detailed discussions).
As expected, the power spectrum ratios shown in the right panel of Figure 13 generally follow the trend of our test f (R) cosmologies (see, e.g., Figure 4).However, the particular scale dependence can depart significantly from the commonly adopted designer and Hu-Sawicki models.Figure 11.Evolution of the modified gravity-to-ΛCDM power spectrum ratio for three test models with scale-independent linear growth on sub-horizon scales.The colored lines represent the output from hi class solving for the full dynamic (FD) of the scalar degree of freedom.The dashed lines are the mochi class predictions for the reconstructed models, where the stable parameters, {∆M 2 * , D kin , c 2 s , ρϕ }, and the initial conditions, α B0 , are computed as detailed in the main text.The excellent agreement between the exact and reconstructed power spectrum ratios, shown in the lower sub-panels, demonstrates that the combination of principal components with a power-law approximation successfully captures the linear growth of structure in these models.
This behaviour can be primarily attributed to (i) a different, possibly time-dependent, gravitational coupling in the sub-Compton limit (i.e., µ ∞ ̸ = µ f (R) ∞ = 4/3), and (ii) a more complex evolution of the Compton wavelength.

SUMMARY AND OUTLOOK
Horndeski gravity provides a comprehensive theoretical framework to explore scalar-tensor extensions of GR, offering rich phenomenology through diverse background expansion and structure formation scenarios.To study this extensive Horndeski landscape, we can restrict ourselves to the energy scales relevant to cosmology by adopting an effective-theory approach.In this work, we have introduced a numerical tool to facilitate the study and statistical analysis of Horndeski models expressed in the effective-theory language.
At the level of the background and linear perturbations, models within the scalar Horndeski sub-class are entirely described by four functions of time.The common α-parametrisation, proposed by Bellini & Sawicki (2014), neatly separates the phenomenology of the background from that of the growth of structure and is im-plemented in Einstein-Boltzmann solver hi class.This code also includes the quasi-static approximation scheme for improved computational efficiency and to extend applicability to models where following the full dynamics of the scalar field is challenging.
Despite these key computational advances, to reliably extend the use of hi class to a broader portion of the Horndeski landscape, we must first address the following critical points: • Simplify the implementation of specific models, which currently still requires modifications to the source code.
• Generalise the input functions to more flexible parametrisations, going beyond {w 0 , w a } for ρ ϕ , and ∝ Ω ϕ (a) or a s for the α-functions.
• Efficiently select stable models, eliminating the need to search for parameter combinations free from ghost and gradient instabilities.
• Improve on the accuracy of the QSA implementation in hi class.
Figure 12.Evolution of the modified gravity-to-ΛCDM power spectrum ratio for two f (R) gravity models-designer and Hu-Sawicki (HS).The colored lines represent the output from mochi class using the exact input stable parameters, {∆M 2 * , D kin , c 2 s , ρϕ }, and initial conditions, α B0 .The dashed lines are predictions from mochi class using the reconstructed input functions, with initial conditions optimised to compensate for reconstruction inaccuracies.The lower sub-panels show that at all examined redshifts the reconstructed power spectra remain within 0.2% of the exact models for k ≲ 1 h/Mpc, and within 1% for scales up to k = 8 h/Mpc.However, reconstructing the HS model with n = 4 (not shown) proved extremely difficult; despite sub-percent errors in the reconstructed stable functions (refer to Figs. 7 and 9), we could not determine initial conditions leading to stable perturbations.
• Correct the erroneous labeling of healthy models as unstable due to limited numerical precision in the calculation of c 2 s .The hi class extension released with this work, mochi class, effectively resolves all the points above by: • Replacing the α-functions directly with the parameters used in the stability conditions, {D kin , M 2 * , c 2 s }.To ensure the absence of ghosts and gradient instabilities, users need only to provide positive input functions.
• Accepting general functions of time as input, with the constraints c 2 sN → 0, M 2 * → 1, and ρ ϕ ≪ ρ m deep in the matter-dominated era.In other words, the model must resemble the ΛCDM cosmology at early times.
• Including mathematical (or classical) stability conditions to catch exponentially growing modes, as done in EFTCAMB.
• Embedding a QSA implementation based on modified Poisson-like equations for the metric potentials.
We validated mochi class's accuracy and performance by comparing its output power spectra and computational efficiency to those of the publicly available Einstein-Boltzmann solvers hi class, EFTCAMB, and MGCAMB.For all tested cases, we found that the predicted changes to the CMB spectra due to new physics matched the output from the other solvers within 0.05%, while for the matter power spectrum, the agreement reached 0.01%.Thanks to its alternative QSA implementation and to the manifestly stable EFT basis, mochi class generally improves on the accuracy and numerical stability attained by hi class, with f (R) gravity being an extreme case that hi class cannot even solve.
In addition, we proposed a new parametrisation for the stable functions and the background energy-density of the scalar field.This approach uses a simple early-time Figure 13.Ratios of the matter power spectrum at z = 0 for randomly generated modified gravity models to that of ΛCDM.Left:: 12 stable models selected from a 31-dimensional parameter space, sampled via a Latin hypercube strategy.We use 6 principal components for each of the input functions, 2 power-law parameters to describe the early-time evolution of ∆M 2 * , D kin and ρϕ , as well as 1 scaling parameter setting the amplitude of c 2 s deep in the matter-dominated era.The initial conditions, α B0 , are determined to ensure the solution to the braiding non-linear ODE (Equation 27) closely follows the expected behaviour at early times (Equation 41).The observed variety of power spectrum shapes demonstrates that the parametrisation of the stable functions proposed here can capture a richer phenomenology than the fixed-form parametrisations, like propto omega, or traditional models, such as Galileons and nKGB.Right: 8 realisations drawn from a pool of stable models obtained by sampling a 16-dimensional parameter space with a Latin hypercube strategy.In this case we only parametrise ∆M 2 * and ρϕ using 6 principal components and 2 power-law parameters.We set c 2 s = 1, and D kin is given by 3α 2 M /2.Initial conditions for the braiding follow the same methodology as the left panel.In this sub-space of the full Horndeski's landscape, models with α K ≪ α 2 B are closely related to f (R) gravity theories.Indeed, the morphology of power spectrum ratios is similar yet more varied than the f (R) gravity test models used in this work.
power-law (or constant) evolution and models late-time smooth deviations with Gaussian Processes, combined with principal component analysis for dimensionality reduction.We showed that this parametrisation can accurately reconstruct our test cosmologies, both in terms of input functions and output power spectra.As an initial application, we used it to generate Horndeski models with either all input functions free or by fixing relations and values to closely resemble f (R) gravity.The linear growth of structure in these models can vary significantly, deviating from the well-known behavior of our test models.This diversity emphasises the necessity of exploring beyond the commonly employed parametrisations (see also Raveri et al. 2021, for a discussion on biased constraints and information loss).
The high versatility of class's Python wrapper, classy, ensures mochi class seamless integration into the pipeline of cosmological analyses using data probing both the background evolution and the linear regime of structure formation.This will provide an alternative avenue of investigation to the reconstruction strategies employed in, e.g., Raveri (2019); Zucca et al. (2019a) and Pogosian et al. (2021); Raveri et al. (2021), and will extend the analysis of studies based on fixed-form parametrisations of the αor EFT-functions (e.g., Salvatelli et al. 2016;Noller & Nicola 2018;Aghanim et al. 2018;Seraille et al. 2024).In a future release, we will also include a special case of the stable parametrisation enforcing the relation α B = bα M , which can be more effective for analyses restricted to models in the generalised Brans-Dicke sub-class (Bergmann 1968), where α B = −α M , or to investigate no-slip gravity modifications (Linder 2018), defined by α B = −2α M .
The larger volume of Horndeski theory space covered by the new parametrisation will likely require scalable sampling methods capable of efficiently handling the high dimensionality of the problem (Mootoovaloo et al. 2024).One such method is Hamiltonian Monte Carlo (e.g., Hajian 2006), a gradient-based sampling technique that leverages the derivatives of the likelihood function to significantly enhance the acceptance rate of samples in highdimensional spaces.To implement this, we could rewrite mochi class using the JAX framework, taking advantage of its automatic differentiation capabilities, as demonstrated by Hahn et al. (2023).Alternatively, a simpler approach might involve using mochi class to generate training data for neural network emulators targeted to specific summary statistics (Mancini et al. 2021;Bonici et al. 2023).These emulators can then be integrated into cosmology libraries that feature automatic differentiation (Piras & Mancini 2023;Campagne et al. 2023).
Finally, for general Horndeski models, mochi class's calculations represent a fundamental first step towards accessing the cosmological information contained on small scales (Heymans & Zhao 2018).For instance, in the halo model reaction framework (Cataneo et al. 2018;Giblin et al. 2019), the linear matter power spectrum and background evolution are key ingredients for obtaining the mean halo statistics and non-linear matter power spectrum of non-standard cosmologies.This method, implemented in the ReACT code (Bose et al. 2020), can be applied to Lagrangian-based theories of gravity (Cataneo et al. 2018;Bose et al. 2021;Atayde et al. 2024;Bose et al. 2024), to phenomenological modifications to GR (Srinivasan et al. 2021(Srinivasan et al. , 2023)), to an interacting dark sector (Carrilho et al. 2021), and to a broader class of parametrised ΛCDM extensions within the Horndeski realm (Bose et al. 2022).ing a private version of the code.We are grateful to Noemi Frusciante and Alessandra Silvestri for their help with EFTCAMB.We would also like to thank Lucas Lombriser for insightful discussions on the inherently stable parametrisation, and Lucas Porth for brainstorming strategies to optimise initial conditions.Finally, a special thank to Alexander Eggemeier for suggesting a memorable name for the code released with this work.We acknowledge the use of the software packages numpy, scipy, matplotlib, and scikit-learn.), while the right panel displays the normalized absolute deviations from the standard cosmological strength of gravity (|∆M 2 /∆M 2 * ,max |).The coloured lines represent functions reconstructed by fitting Padé approximants to their exact counterparts (dashed lines), once the early-time power-law trend has been subtracted.We maintain an identical number of free parameters as in the PCA-based approach to ensure a fair comparison.The lower sub-panels show the reconstruction accuracy of the fittings, indicating that the Padé approximants result in larger deviations from the exact models compared to the PCA-based reconstructions.
C. LATIN HYPERCUBE BOUNDARIES Table 3 and 4 list the parameter ranges defining the Latin hypercubes used to generate random models belonging to the Scalar Horndeski and f (R)-like classes, respectively (see Section 5.3).Here we remind the reader that overdots denote derivatives with respect to ln a, and primes (used below) represent derivatives with respect to conformal time.The choice between ln a or τ as the time variable is dictated by the original class structure, with the former used in the background module and the latter in the perturbations module.
D.2.Scalar field fluctuations The coefficients and the source term in Equation ( 23) can be found in Appendix A.2 of Zumalacarregui et al. (2016).For convenience, we provide them here: where D.3.Quasi-static approximation In the effective-field formulation of the quasi-static approximation, the gravitational coupling for the non-relativistic matter and the gravitational slip can be expressed as in Equation ( 31).The small-scale limits can be expressed in terms of the α's as Here, we also present a typo-corrected equation for the metric perturbations in synchronous gauge, η, used in our MGCAMB patch to hi class: a 2 9ρ m a 2 µγ(1 + w m )/2 + k 2 3ρ m (1 + w m )µγθ 1 + 3 where H = aH is the conformal Hubble parameter.The corrections compared to Equation (26) in Zucca et al. (2019b) include removing an extra power of two from H ′ in the first line, and changing the minus sign to a plus sign in front of the γ ′ µ term in the last line.

Figure 4 .
Figure 4. Extending Figure3, here we consider modified gravity models characterised by a scale-dependent growth on sub-horizon scales.The reference power spectra for f (R) gravity are computed with EFTCAMB (dashed) and MGCAMB (dotted).The lower panels focus on the fractional deviation of mochi class from the CAMB-based solvers, providing a direct comparison for the matter power spectrum ratios shown above.The Hu-Sawicki (HS) models demonstrate a close match, with deviations from MGCAMB remaining under 0.025% across all scales.In contrast, for the designer f (R) model, mochi class shows a consistent discrepancy that increases at smaller scales (blue line), a result predominantly due to the quasi-static approximation.This difference, in fact, significantly reduces when mochi class always integrates the full scalar field dynamics (green line).

Figure 5 .
Figure 5. Evolution of the de-mixed kinetic term, D kin , normalised to its maximum value, D (max) kin

Figure 7 .
Figure 7. Reconstruction of the stable basis functions, with the left panel displaying the normalized de-mixed kinetic term (D kin /D (max) kin

Figure 8 .
Figure 8. Reconstructed (coloured) and exact (dashed) speed of sound for the test models exhibiting non-trivial evolution (i.e., c 2 s ̸ = 1).The lower panel demonstrates that the PCA-based reconstruction can achieve better than 1% accuracy at all times.1. Compute the natural logarithm of the exact D kin .2. Fit a power-law over the range [−5, −3] to find the best fit parameters, ζ D and b D ; 3. Subtract the power law from the exact ln(D kin ) to obtain ∆ D ; 4. Project ∆ D onto the six principal components to get ∆D = Additionally, Mathematica was utilized for comparison purposes and to compute the stable basis functions of some test cosmologies.MC acknowledges the sponsorship provided by the Alexander von Humboldt Foundation in the form of a Humboldt Research Fellowship for part of this work.He also acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research.EB has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement No 754496.